Cesium+GPU计算：实现实时泥石流地形侵蚀模拟的实践

做地质灾害数字孪生的同行，大概都有过这种体验：Cesium 场景本身行云流水，加影像、叠倾斜摄影、切 3D Tiles 都是很成熟的用法，可一旦想让地形“活”起来——不是播放一段预烘焙好的动画，而是实时演算泥石流如何掏蚀坡脚、冲刷沟道、在沟口堆出冲积扇——CPU 端的模拟立刻成了那块短板。分辨率上到 512×512，单帧迭代动辄几十毫秒，Cesium 的主线程又被渲染和交互占得满满当当，整个项目瞬间从流畅的六十帧掉到个位数。

这篇文章要聊的，就是我在一个流域级地质灾害模拟项目里被这块短板逼出来的解法：在 Cesium 中实现 GPU 计算的泥石流地形侵蚀。核心思路一句话概括——把过去 CPU 上逐格点更新的流体侵蚀方程，全部塞进 WebGL 的片段着色器里，用地形高度图和流动参数纹理做输入，渲染到纹理，逐帧迭代，再把计算结果叠加回 Cesium 的三维地形或模型表面。这样做的直接收益是：模拟过程从“卡到没法交互”变成“能跟着镜头走”，而且侵蚀结果天然就是一张纹理，做热力图叠加、剖面分析、动画推演都特别顺手。

1. 为什么把侵蚀模拟搬进 GPU：CPU 算力逼出来的技术选型

1.1 “逐帧演化”要求的算力，CPU 确实给不起

泥石流地形侵蚀和普通地形可视化的本质区别在于“演化”二字。可视化只需要把静态数据画出来，而侵蚀模拟要在地形上持续做时间推进：每个时刻的水流方向、水深、剪切应力、泥沙浓度、地形高程都在互相影响，下一帧的结果依赖上一帧的所有场量。

我算过一笔账：一个 256×256 的规则网格，总共 65,536 个单元格。每个单元一次迭代要做水流积累、坡度计算、流向判断、输沙能力计算、侵蚀/沉积更新，保守估计 150~300 次浮点运算。65,536 乘以 200，单次迭代就是约 1,300 万次浮点操作。你可能会觉得这对现代 CPU 不算什么，但问题是 JavaScript 里跑这个循环不光是浮点运算，还有数组访问、对象属性读取、分支判断、临时变量分配，一套下来单次迭代实际耗时经常在 20~50 毫秒。

而泥石流模拟不是迭代一次就完事。模拟 10 分钟的泥石流过程，时间步长受流速限制通常只有 0.05~0.2 秒，也就是说要跑 3,000 到 12,000 次迭代。即便把每次迭代优化到 20 毫秒，完整的模拟也要 60 秒以上，这还是在 CPU 不算其他任务的情况下。

到了 GPU 这边情况完全不同。同样的 256×256 网格在 GPU 上就是一个 256×256 的纹理，片段着色器天然按像素并行执行，一次迭代的所有单元格同时计算。主流独立显卡有上千个流处理器，这种规则数据流的负载几乎是 GPU 最擅长的场景，单次迭代通常能压到 0.5~2 毫秒。差距是两个数量级起步，这也是我把方案确定为 GPU 计算的直接原因。

1.2 为什么是 Cesium：它提供场景，却不该承担计算

Cesium 在全球坐标、地形瓦片、3D Tiles、倾斜摄影这些数据底座上的积累确实是同类方案里最完整的。数字孪生项目要的是一个能承载“真实地理位置 + 真实地形 + 真实业务数据”的容器，Cesium 在这个层面的优势无可替代。

但 Cesium 骨子里是个渲染引擎，不是仿真引擎。它的地形是基于 Quantized-Mesh 的瓦片化网格，服务端生成好之后客户端再加载，你没法在每帧里直接改某个顶点的高程并且指望它高效重算。所以我的做法是：Cesium 负责“世界表达”，GPU 计算负责“物理演化”，两者通过纹理和自定义 Primitive 建立连接，各干各的活。

这里有一个容易被忽略的前提：在 Cesium 里做自定义 GPU 计算，必须跟 Cesium 共享同一个 WebGL 上下文，所以不能随便搞一套独立的 WebGL2 上下文再叠上去。我采用的模式是在 Cesium 的帧回调里挂一个自定义渲染函数，等 Cesium 这一帧的地形和模型都画完之后，再执行模拟 Pass，把结果渲染到离屏 Framebuffer。下一帧再用这个结果做可视化叠加，或者继续参与模拟迭代。

这种“先展示后计算”的顺序有个额外好处：模拟计算不会阻塞 Cesium 场景渲染，用户操作相机、加载瓦片、拾取对象都保持流畅。代价是计算结果会延后一帧显示，但对于泥石流这种时效性要求没那么苛刻的过程模拟，一帧延迟根本感知不到。

1.3 RenderTarget Ping-Pong：GPU 模拟的地基

在 GPU 上做逐帧迭代模拟，最经典的架构就是双缓冲渲染到纹理，行话叫 Ping-Pong。原理和 CPU 端的双缓冲一样：准备两个 Framebuffer 对象（FBO），每个 FBO 绑定一张纹理作为颜色附着点。奇数帧从 A 纹理读数据，计算结果写入 B 纹理；偶数帧反过来从 B 读、往 A 写。读写对象永远分离，避免同一张纹理同时作为输入和输出导致的未定义行为。

Cesium 里实现这一步要注意的是：不要自己创建一个全新的 WebGL2 渲染上下文，而是拿到 Cesium 的上下文，在其上创建 Framebuffer 和纹理。Cesium 官方 API 不直接暴露 context，但可以通过 scene.context 拿到底层对象，或者直接用 scene.postRender 事件回调里执行自己的 GL 调用。这要求你对 WebGL 状态管理有一定功底，因为 Cesium 每帧会设置一堆 scissor、viewport、depth test、blend 状态，你做完模拟 Pass 之后必须恢复，否则下一帧 Cesium 的地形会出现莫名裁剪或者半透明叠加异常。

2. 塞给 GPU 的“地形”长什么样：从 DEM 到 RGBA 高度场

2.1 数据来源与范围裁剪：Cesium 地形不等于模拟地形

很多人的第一反应是直接用 Cesium 自带的地形数据，比如 Cesium World Terrain，省去自己找数据的麻烦。但这里有个硬伤：Cesium 的 Quantized-Mesh 地形瓦片是不规则三角网，每个瓦片内部、瓦片之间的顶点密度都不一样，没法直接转成规则网格高度场供 GPU 计算使用。

正确做法是从数据源头拿规则网格 DEM（数字高程模型），来源可以是本地 GeoTIFF 文件，也可以是地形服务商提供的栅格高程接口，或者从 Cesium 地形瓦片里先做一次离线栅格化导出。关键点是：模拟区域的 DEM 必须裁剪成矩形，并且空间分辨率一致。

我处理的时候会把 DEM 先裁剪到流域边界的外接矩形，再在矩形范围内重采样。重采样方法上，简单项目用双线性内插就够，但如果地形起伏大，为了保证坡度和流向计算不出现伪边界，建议用带抗锯齿的高阶重采样。这一步在 CPU 上离线做就行，做一次后面全部复用。

2.2 高程编码：单通道浮点纹理与 RGBA 打包之争

准备好规则 DEM 之后，要把它写成 GPU 能采样的纹理。这里有一个关键选择：用单通道浮点纹理（R32F），还是把高程编码到 RGBA 8 位通道里。

WebGL2 环境下的 Cesium 项目，理论上可以直接用浮点纹理。但实际踩坑下来，不是所有设备的浮点纹理都支持线性过滤，有些移动端 GPU 对 R32F 的 Framebuffer 完整性检查会直接报错，导致模拟 Pass 无法运行。如果你的目标平台是桌面端 WebGL2，R32F 可以优先考虑，写法干净，精度也足够。但如果要考虑移动端或者跨浏览器兼容，RGBA 打包是更稳的方案。

我最终在项目里采用了 RGBA 打包方案，核心解码函数是这样：

vec4 encodeFloat(float v) { vec4 enc = vec4(1.0, 255.0, 65025.0, 16581375.0) * v; enc = fract(enc); enc -= enc.yzww * vec4(1.0 / 255.0, 1.0 / 255.0, 1.0 / 255.0, 0.0); return enc; } float decodeFloat(vec4 rgba) { return dot(rgba, vec4(1.0, 1.0 / 255.0, 1.0 / 65025.0, 1.0 / 16581375.0)); }

这个编码的核心思想是把一个 float 的二进制表示拆成四段，分别存进 RGBA 四个 8 位通道。解码时做加权求和。要注意的是，这种编码方式对取值范围很敏感：如果高程是 1,200 米，直接编码没问题，但如果模拟过程中高程可能出现负值或者超过编码上界（比如 16,777,215 米的极限值），就要先做“归一化+偏移”，把数据映射到 0 到 1 的安全区间。我在项目里把高程统一减去区域最低点，再除以区域最大高差，保证编码过程不越界。

2.3 分辨率与真实尺寸的换算：每个像素代表多少米

这一步是很多人忽略但特别关键的。上传到 GPU 的纹理虽然只是一个二维数组，但它每个像素都对应着真实世界里的一块地面面积。模拟计算里的流速、剪切应力、侵蚀量全部依赖真实的物理尺寸，纹理像素和实际距离的换算关系必须清清楚楚。

计算方法很简单：模拟区域东西向长度除以纹理宽度，就是 X 方向的像素分辨率；南北向长度除以纹理高度，就是 Y 方向的像素分辨率。比如模拟区域是 2 公里 × 2 公里，纹理用 256×256，那每个像素对应约 7.8 米。Cesium 世界里还要考虑经纬度在不同纬度上的距离差异，所以我通常先把经纬度坐标投影到本地 ENU（东北天）坐标系，再建立 ENU 坐标与纹理像素坐标的映射。

着色器里计算坡度的时候，采样步长要用真实的物理距离而不是像素距离：

float dx = u_pixelSizeX; // 每个像素对应的真实米数 vec2 slope = vec2( (getHeight(uv + vec2(dx, 0.0)) - getHeight(uv - vec2(dx, 0.0))) / (2.0 * dx), (getHeight(uv + vec2(0.0, dx)) - getHeight(uv - vec2(0.0, dx))) / (2.0 * dx) );

如果这一步不做真实距离换算，直接用像素坐标算坡度，模拟出来的泥石流流速和侵蚀分布会严重失真，看起来像模像样，实际上完全不对。

3. 泥石流模拟的核心：GPU 里到底在算什么

3.1 模型选型：从完整流体力学到工程可用

泥石流严格来说是固液两相流，里面有水、有泥沙、有大石块，还要考虑屈服应力、粘塑性本构关系，完整求解 Navier-Stokes 方程在这个场景下完全不现实。即使简化成浅水方程，也要面对地形剧烈起伏带来的激波问题。

工程项目的现实是：我不需要一个严格符合流体力学教科书的理论模型，我需要一个在 GPU 上跑得动、趋势正确、参数可调、结果能说服业务方的简化模型。所以我在项目里采用的是一个分级近似方案：

• 用深度积分的浅水方程描述泥浆流动，忽略垂直方向的速度分量
• 引入 Bingham 塑性模型描述泥石流的屈服特性：当底部剪切应力小于屈服应力时，泥浆不流动；超过屈服应力后，流速与剪切应力超出的部分近似线性相关
• 用经验侵蚀公式估算泥浆对沟床和沟岸的冲刷能力，结合输沙能力判断侵蚀还是沉积

这个方案的好处是每一步都有明确的物理或经验公式支撑，同时又足够简单，可以全部塞进片段着色器。

3.2 核心方程与着色器实现思路

先说水流场部分。在 GPU 里，每个像素代表一个地形单元格，相邻像素代表空间相邻的地形单元。我维护三张核心纹理：高程纹理、水深纹理、泥沙浓度纹理。

水流速度采用曼宁公式近似：

v = (1 / n) * R^(2/3) * S^(1/2)

其中 n 是曼宁糙率系数，R 是水力半径（浅水情况下可以用水深 h 近似），S 是坡度的模长。这个公式在 GPU 里的优势是全部是局部量，不涉及跨单元格的复杂差分，非常适合片段着色器。

泥石流和普通水流的关键区别是屈服应力。我在速度计算里加了一个阈值判断：

float tau_b = rho * g * h * slopeMag; // 底部剪切应力 float velocity = 0.0; if (tau_b > u_yieldStress) { float tau_excess = tau_b - u_yieldStress; velocity = (1.0 / n) * pow(h, 2.0 / 3.0) * sqrt(slopeMag); velocity *= smoothstep(0.0, u_yieldStress * 0.5, tau_excess); // 用平滑过渡替代硬分支 }

这里用 smoothstep 替代 if 判断，是为了减少 GPU 着色器里的动态分支开销。GPU 的并行执行模型对 warp 内分支很敏感，同一个分支路径下的像素越多效率越高，smoothstep 让不同单元格的流速连续变化，既避免了分支惩罚，也让泥石流前锋的推进看起来更自然。

侵蚀和沉积部分我采用输沙能力模型。每个单元格有一个输沙能力 C_eq，它跟流速和水深正相关：

C_eq = K_erosion * pow(velocity, alpha) * pow(h, beta)

当单元格当前的泥沙浓度高于输沙能力，超过部分就沉积下来，地形高程增加；低于输沙能力，就从下方地形侵蚀补充泥沙，地形高程降低。着色器里对应的伪代码是：

float capacity = u_erosionK * pow(velocity, u_alpha) * pow(h, u_beta); float deposit = max(0.0, sediment - capacity) * u_dt; float erode = max(0.0, capacity - sediment) * min(1.0, u_erodeRate * u_dt); height -= erode; height += deposit; sediment += erode - deposit;

这里有一个关键参数 u_erodeRate，它反映了基岩和松散堆积物的抗侵蚀能力差异。泥石流沟道里堆积物厚的地方侵蚀速率高，基岩出露的地方侵蚀速率要压低一个数量级。我习惯在着色器外面准备一张可侵蚀性纹理，不同区域涂上不同的侵蚀系数，效果比全局单一参数真实得多。

3.3 时间步长、数值稳定性与边界条件

GPU 模拟和 CPU 模拟同样受数值稳定性约束，只是它跑得快，所以可以承受更小的时间步长。我用的是显式欧拉时间推进，稳定性条件遵循 CFL 条件：

dt <= 0.2 * cellSize / v_max

举个例子：像素分辨率是 8 米，泥石流最大流速 10 米/秒，那么时间步长不能超过 0.16 秒。实际我会在这个基础上打对折，取 0.08 秒左右，宁可多迭代几步也不冒发散的风险。

边界条件也容易踩坑。我的模拟区域是流域边界裁剪出来的矩形，矩形边上的单元格不能当成普通内部单元格处理。我在边缘加了一圈“不侵蚀缓冲带”：边缘两圈的高程不参与侵蚀更新，水流可以流出区域边界（做吸收边界），但边界地形保持不变。这样既模拟了泥石流冲出流域范围的效果，又不会因为边界条件处理不当产生人工坑洞。

4. 把模拟结果变成 Cesium 里看得见的地形变化

4.1 动态修改 Cesium 地形几乎走不通，我换了思路

刚开始做这个项目时，我一度想直接让 Cesium 的地形跟着模拟变化。试了之后发现这条路很难走通：Cesium 的地形瓦片是服务端预生成的 Quantized-Mesh 数据，客户端的 GlobeSurface 有自己的一整套 Tile 加载、LOD 切换、顶点更新机制，你想在每帧改它某个区域的高程，就得深入 Cesium 的 surface 更新流程，改完还得应对它内部的地形缓存和法线重算。这条路工程复杂度太高，几乎没有团队愿意为这个改动长期维护 fork 版本。

所以我换成了“动态网格叠加”的方案：不用动 Cesium 自带地形，而是自己生成一个覆盖模拟区域的规则网格 Primitive，这个网格的顶点高度每一帧从 GPU 模拟出来的高度纹理里读取，然后叠在 Cesium 地形上方。网格和 Cesium 地形之间用深度偏移错开，视觉上就像地形真的被泥石流改变了。

这个思路的精髓在于：展示层和计算层解耦。Cesium 地形永远是最新加载的原始地形，而模拟带来的地形变化全部体现在叠加网格上。当模拟结束，你甚至可以把这个叠加网格的高度场导出来，作为下一次模拟的初始地形用，实现“多次模拟累积演化”的效果。

4.2 用自定义 Primitive + 顶点着色器实现动态地形叠加

在 Cesium 里实现这个叠加层，我选择自定义 Primitive 而不是简单的 Entity。自定义 Primitive 的好处是几何体和材质都完全可控，而且可以直接在顶点着色器里采样 GPU 模拟纹理。

几何体部分比较简单：根据模拟区域的像素分辨率，生成一个 (N-1)×(N-1) 的规则网格，每个顶点存下对应的纹理采样坐标。地形变化缓存的插值需要高精度，所以 vertex attribute 里用 floats 保存。

外观部分才是核心。Cesium 的 Appearance 允许自定义 vertex shader 和 fragment shader。我在顶点着色器里把高度纹理作为 uniform sampler2D 传入，从纹理里采样出当前高度，在顶点做位移：

attribute vec3 position3DHigh; attribute vec3 position3DLow; attribute vec2 st; uniform sampler2D u_heightTexture; uniform mat4 u_modelViewProjection; void main() { float h = texture2D(u_heightTexture, st).r; vec3 pos = vec3(position3DHigh + position3DLow); pos.z += h * u_heightScale; gl_Position = u_modelViewProjection * vec4(pos, 1.0); }

这里有个关键问题：纹理采样结果在顶点着色器里默认是 GPU 插值过的，对于地形高度这种高频信号，网格密度不够会导致顶点之间出现“拉链状”的锯齿。我的解决办法是让网格密度不低于模拟纹理的分辨率，或者干脆 1:1 对齐，每个顶点对应纹理的一个像素。这样顶点高度就是像素中心的精确值，不依赖插值。

Cesium 的模型坐标系比较复杂，有 position3DHigh 和 position3DLow 两个 attribute，这是为了在 64 位浮点精度下表示全球坐标。做顶点位移时，我是在 ENU 局部坐标里位移完再转回世界坐标。这个细节如果处理不好，叠加网格会在远离场景中心的位置出现明显抖动。

4.3 不只改地形：把模拟结果做成热力叠加和剖面分析

地形位移是最直观的表达，但有些场景下客户更关心的是“哪里侵蚀最严重”“哪里淤积最多”，这时候直接用色带叠加比地形变形更有冲击力。

我把模拟生成的“高程差纹理”在 fragment shader 里做了一次映射：正值（沉积区）用暖色到红色的渐变，负值（侵蚀区）用蓝色到紫色的渐变，零值附近透明。然后把这个半透明材质叠加在 Cesium 地形或倾斜摄影表面上，效果类似于热力图。

这个方案在密级业务里特别实用，因为不需要做任何地形 mesh 修改，也不会有闪烁问题，一张纹理贴上去就完事。侵蚀量、淤积量、水深、流速都可以做成不同的色带，切换观看维度只需要换一张纹理。我还做了一个剖面分析工具：在 Cesium 里拉一条线，沿线采样模拟纹理的值，生成侵蚀前后的地形剖面曲线，业务方看了直点头。

4.4 叠加层与地形的 LOD、闪烁问题

动态网格叠加生效之后，第一个跳出来的问题就是地形闪烁。这是因为叠加网格和 Cesium 原生地形在深度上几乎重叠，GPU 深度缓冲的精度有限，远近交替绘制时就会产生 z-fighting。

我试了三种解决办法。第一种是给叠加网格做 polygonOffset，让它在深度测试时偏近一点，算是见效最快的一种。第二种是修改 Material 的 depthFunc，把默认的 LEQUAL 改成 LESS，强制叠加网格只有严格在地形前面才显示，避免深度相等时的歧义。第三种是把叠加网格抬高一个微小的高度偏移，比如 0.05 米，靠视觉误差掩盖深度冲突。

实际最稳的是 polygonOffset + 微高程偏移组合。原因在于 polygonOffset 只影响多边形光栅化阶段的深度，对自定义 Primitive 的兼容性在不同显卡上表现不一致；而高程偏移则简单可靠，只要偏移量设置得当，看不出地形被“抬起来”的痕迹。

另外一个容易被忽略的问题是 Cesium 地形 LOD 切换。当相机拉近拉远，Cesium 会加载不同级别的地形瓦片，新瓦片刚加载出来的瞬间地形高度和旧瓦片不一致，叠加网格还保持原来的高程，就会突然出现一大片地形“塌陷”或“凸起”的视觉跳变。我的处理是监听 Cesium 的 terrainProvider 的 readyPromise，在瓦片加载完成后强制重建一次叠加网格，同时把新采样出来的高度做一帧线性过渡，避免硬切。

5. 踩坑记录：浮点精度、RTT 复用和令人头秃的“地形抖动”

5.1 RGBA 编码的高程在迭代 100 次后出现了条带断层

第一次把这个管线跑通，我兴奋地让模拟跑了 200 帧，结果发现高程纹理上出现了明显的条带断层，一条一条的像年轮一样。排查了半天，原因出在 RGBA 编码的精度上。

RGBA 8 位打包方案的理论精度是 1/16777216，看起来很高，但实际上在多次迭代后，高程的微小变化量（比如每一轮侵蚀掉 0.001 米）经过编码、解码、再编码的循环，低位字节的量化误差会被逐步放大。特别是当高程数值比较大的时候，误差更明显。

我的解决方法是加“双通道存储”：一张纹理用 RGBA 打包存高程的整数部分，另一张纹理存小数部分的高精度增量。每次迭代读取两张纹理，叠加后做计算，更新时只把变化量写回增量纹理。缺点是多占一倍显存，但条带问题彻底解决了。后来 WebGL2 环境换用 R32F 纹理后，这个问题自然消失，但移动端的兼容性方案里我仍然保留着双通道编码逻辑。

5.2 Ping-Pong 状态管理翻车：读写了同一张纹理

GPU 模拟最经典的状态管理陷阱就是读写冲突。我一度为了省一次 FBO 切换，尝试用一个 Framebuffer 既作为当前 Pass 的输出，又作为下一个 Pass 的输入，结果出来的地形像被泼了硫酸一样四处扩散。

原因是 WebGL 规范里，帧缓冲的对应纹理如果在同一 Pass 内既绑定为颜色附着点又被采样，属于未定义行为，不同 GPU 驱动处理方式完全不一样。NVIDIA 显卡上可能看起来正常，AMD 或者集成显卡上就会产生不可预期的读写竞争。

正确做法一定是严格双缓冲：两个 Framebuffer，每个绑定独立的颜色纹理，模拟 Pass 之间无条件交换读写角色。纹理创建用gl.DYNAMIC_DRAW或者gl.STREAM_DRAW提示，让驱动优化显存访问路径。

5.3 在 Cesium postRender 里做模拟，viewport 和 scissor 被 Cesium 带偏

这个坑我印象太深了。模拟出来的结果在离屏 Framebuffer 里看完全正常，但一叠加回 Cesium 场景，地形变化区域的边缘就出现奇怪的裁剪，像是被人用剪刀沿着屏幕边界剪了一道。

排查后发现问题出在 Cesium 每帧渲染结束后，viewport 和 scissor 的状态是它最后一次绘制时的状态，不是全屏状态。我在 postRender 回调里直接创建 Framebuffer 做模拟，viewport 和 scissor 还停留在 Cesium 某个局部 UI 窗口的尺寸上，结果模拟只渲染到了纹理的一小块区域。

解决办法是在模拟 Pass 开始前显式设置：

gl.viewport(0, 0, textureWidth, textureHeight); gl.disable(gl.SCISSOR_TEST); gl.disable(gl.DEPTH_TEST); gl.disable(gl.BLEND);

模拟结束后再恢复到 Cesium 期望的状态。如果不恢复，下一帧 Cesium 的 UI 覆盖层会莫名其妙地多出来或者消失，这个 bug 更难查。

5.4 模拟区域跨越多个地形瓦片，网格跟随出现“断崖”

项目里模拟区域覆盖了一个小型流域，横跨了 Cesium 的多个地形瓦片。运行一段时间后，叠加网格的边缘出现了明显的“断崖”，一侧是正常地形，另一侧像被砍了一刀。

原因在于 Cesium 不同 LOD 级别的地形瓦片，在公共边上的顶点高程可能不一致（正常情况误差很小，但接缝处偶尔会有一两米的偏差）。当叠加网格紧贴地形表面渲染时，接缝处的高度差通过叠加网格表现出来了。

我的应对策略是把叠加网格的边界向外扩一圈，让网格边缘离开高精度地形瓦片的接缝，同时对外圈顶点的高程做一个平滑衰减，强制在外圈 10 米范围内把模拟地形过渡回原始地形。这样模拟区域边缘不会出现突然的高程跳变，泥石流沟道变化看起来也更自然——毕竟真实的泥石流也不会在一个整齐的矩形边界上戛然而止。

6. 性能优化：从 512×512 到 2048×2048 的实战调优

6.1 用 GPU 计时器量化，别靠猜

做性能优化第一步不是改代码，而是拿到准确的性能数据。WebGL 里有个扩展EXT_disjoint_timer_query_webgl2，可以用来测量 GPU 上一个 Draw Call 或一组调用花了多少时间。

我封装了一个简单的测量工具，在模拟 Pass 前后插入计时查询，把每次迭代的毫秒数打到控制台。真正调优之后发现，性能瓶颈往往和我的直觉不一致：有时候瓶颈在纹理采样数量上，有时候在 FBO 切换频次上，还有一次居然卡在高质量过滤设置上——移动设备对纹理线性过滤的支持比桌面端差很多。

6.2 MRT（多渲染目标）把模拟 Pass 从三个合并成一个

模拟的每一帧迭代，我需要更新高程、水深、泥沙浓度三张纹理。最初我写了三个 Pass：先更新水深，再更新泥沙，最后更新高程。每个 Pass 都要读写一遍纹理，带宽开销很大。

后来我把三个 Pass 合并成一个 Pass，用 GLSL 里的gl_FragData或 WebGL2 的drawBuffers同时输出三张纹理。这种方式叫多渲染目标（MRT），只做一次 FBO 绑定、一次 Draw Call，全部更新在同一趟完成。

实测数据变化很明显：1024×1024 分辨率下，三个 Pass 总耗时约 4.8 毫秒，合并成 MRT 后降到 2.1 毫秒，帧成本直接少了一半多。代价是着色器里逻辑变得复杂一些，要为三个输出分别做计算，但这个复杂度对性能的回报非常划算。

6.3 半浮点纹理 HIghp 到 Mediump 的取舍

WebGL 着色器里的浮点精度也是可以调的。默认情况下 vertex shader 用 highp，fragment shader 用 mediump 或者 highp 看平台实现。泥石流模拟对水深、流速这些物理量的精度要求没那么苛刻，我就把 fragment shader 里大部分中间变量声明成 mediump，高程纹理采样和高程更新仍然用 highp。

这个改动在高分辨率纹理下效果很明显，但也埋了个坑：个别移动平台对 mediump 的实现不标准，某些变量的精度不足以支撑曼宁公式里的指数运算，导致流速计算结果出现异常大的值。后来我的处理是只对泥沙浓度、沉积量这类相对温和的变量降精度，水深、流速、剪切应力保持 highp，平衡精度和性能。

6.4 分辨率选择与帧预算的平衡建议

不同分辨率的实测数据大致如下，给后面做项目的同行一个参考：

这里有个容易误判的地方：堆叠每帧迭代次数时，会显著提高 GPU 占用，导致 Cesium 场景本身的帧率下降。Cesium 场景在镜头旋转、加载倾斜摄影瓦片时本来就吃 GPU，如果模拟把 GPU 算力抢光，用户操作卡顿就是必然的。

所以我的经验是设定每个帧的“模拟预算”上限，比如 4 毫秒。如果当前帧的模拟 Pass 时间快超了，就在 CPU 端通过迭代计数器把这帧的迭代次数减半。这样可以保证模拟在 GPU 空闲时尽量跑满，在 GPU 紧张时自动退让，用户体验始终稳定。

还有一个小技巧：当需要高频模拟时，把模拟分辨率和可视化分辨率分开。模拟用 512×512 保证实时性，可视化叠加层的网格密度用 256×256，靠纹理采样插值出近似效果。人眼对地形变化的敏感度其实没有想象中高，过度追求网格密度只会白白消耗 GPU。

我在实际项目里还试过一次把模拟分辨率降到 128×128 再上采样，用来做超过 30 分钟的长序列泥石流演算。结果发现泥石流沟道的关键形态特征在低分辨率下丢失太多，后续分析阶段被业务方质疑了好几次。所以在分辨率选择上不要一味求低，至少保证沟道宽度能覆盖 3~5 个像素，低于这个阈值就用更高分辨率或缩小模拟区域。

回过头来看，这套“Cesium + GPU 计算 + 动态网格叠加”的架构，核心价值不在于某一个算法有多高级，而在于它把“场景展示”和“物理模拟”这两件天然冲突的事情在 WebGL 的同一个上下文里和谐地放在了一起。只要能控制好纹理编码、双缓冲状态和 GPU 帧预算，泥石流侵蚀模拟在 Cesium 里做到实时交互是完全可行的。后续我还在探索把这个方案扩展到滑坡、洪水演进和多灾种耦合模拟，底层的数据流和渲染架构基本可以复用，改的只是着色器里的物理方程和参数。