这段时间在做一个数字孪生方向的灾害可视化项目核心目标是把泥石流对地形的侵蚀过程搬到浏览器里的三维地球上。项目标题里的关键词是“Cesium中实现GPU计算的泥石流地形侵蚀”听起来很硬核但拆开看就三件事地形数据怎么进、侵蚀计算怎么跑、计算结果怎么在Cesium里显示。这篇文章我就把整个思路、关键代码和踩过的坑一起梳理一遍希望能给同样在做WebGIS、数字孪生和灾害模拟的朋友一点参考。先说背景。传统做法是在服务端用CPU跑水文、泥沙模型算完再切成瓦片发布。但地形侵蚀天然是个瞬态过程你想看到一个山坡上的泥浆往下糊、沟谷一点点被掏深、坡脚慢慢堆出扇形堆积体就必须高频迭代。CPU在Python或Node里循环遍历几万个网格一个时间步要几十毫秒操作视角时很容易卡。所以我们把状态全部放到GPU纹理里用片段着色器做并行迭代Cesium只负责把迭代结果转成可视的三维地表。实测下来256x256的地形网格单帧能迭代20到30次场景旋转和缩放基本不会掉帧。这篇内容比较适合已经在用Cesium做数字孪生、WebGIS可视化但对GPU通用计算还不太了解的朋友。我尽量把物理模型简化到能写代码的程度把每个环节为什么这么做讲清楚。同时也会把工程里的坑一个个列出来包括参数调不好导致地形像豆腐渣、readPixels吃性能、坐标精度导致高度纹丝不动之类的细节。跟着走一遍至少能让你的Cesium场景里真正出现一座正在被侵蚀的山区。1. 项目思路泥石流侵蚀模拟到底在算什么1.1 把泥石流简化成“高度场水流沉积物”的物理方程泥石流在数值模拟里有很多建模方式最严格的是SPH、有限体积法求解两相流但在浏览器里、在地形尺度上根本跑不动那么重的物理。实际干活时我用的是“高度场逼近”把研究区域离散成一张二维网格每个格子记录三个量——地表高度H、地表水或泥浆量W、沉积物量S。你想象一个长方形的沙盘沙盘被切成很多小格子每个格子有一个数字表示沙堆高度。降雨从高处落下水带着沙往低处流有的格子被冲刷变低有的格子因为流速下降把沙留下来变高。整个过程的数学核心就是通量守恒水流加速、坡度变陡的地方床面被侵蚀水流减速、坡度变缓的地方沉积物倾向于留下。只要把这个趋势算对视觉上的“泥石流侵蚀”就立住了。我用的简化方程大致是这样先算局部坡度用中心网格和相邻网格的高度差除以格子尺寸然后根据水量的多少和坡度超过休止角的程度计算侵蚀量再根据沉积物量和坡度低于休止角的程度计算沉积量最后把侵蚀和沉积的差累加到高度变化上。水量和沉积物各自做更新形成一个完整的循环。局部坡度: slope max(|H(x,y)-H(neighbor)|) / cellSize 侵蚀量: E k_e * W * max(0, slope - restAngle) 沉积量: D k_d * S * max(0, restAngle - slope) 高度变化: dH -E D 水量更新: W W inflow - outflow 沉积物: S S E - D泥石流和普通水流的差别主要体现在参数上。泥浆粘性大休止角比干沙或纯水大同样坡度下不容易启动一旦启动又容易在坡脚快速堆积。所以代码里的restAngle我一般取0.2到0.3弧度比普通冲沟模拟要更大一点。这样简化之后状态量少、迭代快而且每一步都是逐格独立的局部计算很适配GPU“每个格子一个线程”的工作方式。1.2 为什么选Cesium GPU而不是CPU逐格迭代选Cesium不是因为它能计算而是因为它给了一个现成的三维地球容器。我们需要在不同视点观察泥石流过程叠加影像、地形、倾斜摄影甚至3DTiles模型Cesium负责这些“皮”GPU负责“里”。两者结合才能让物理模拟结果真正长在地球上。CPU模拟慢在哪每个时间步要遍历网格、判断邻域、写入新的数组。Python里一个256x256的网格循环一遍大概几十万次操作看起来不多但我们要在3D地球旋转漫游时保持60帧留给业务的逻辑时间只有十几毫秒CPU很容易被打满。GPU方案把计算写成片段着色器一次draw call就能对整个网格做一次更新本质是数据并行。场景CPU逐格迭代GPU纹理迭代网格规模64x64够用杀鸡用牛刀256x256以上卡顿明显流畅需要每帧迭代多次不推荐合适需要读回数值做分析方便需要readPixels注意性能与Cesium原生结合直接改数组需要处理FBO和纹理我实际测试下来256x256的网格CPU模拟每个时间步大概需要30到50毫秒GPU用WebGL2跑20次迭代加一次读回能控制在10毫秒以内。差距非常明显所以项目一开始就把技术路线定成了WebGL2 Cesium。1.3 整体架构与数据流整个项目的数据流是这样设计的拿到目标区域的经纬度范围通过Cesium的地形采样接口或本地DEM文件得到地表高程点。把高程点插值成规则二维网格写进Float32Array作为初始状态纹理。初始化两张浮点纹理状态A、状态B一个用于读取当前状态一个用于写入下一状态。在渲染循环里调用GPU迭代着色器用乒乓切换的方式不断交换状态纹理。每隔若干帧从GPU纹理中读回高度数据更新Cesium自定义Primitive的顶点高度。根据高度场重新计算法线让地形有明暗变化再叠加流动纹理和粒子表现泥石流的运动方向。这里最关键的设计是“状态纹理是唯一数据源”。Cesium里的Mesh只是GPU状态纹理的一份投影不能反过来回写。这样做的好处是数值迭代完全在GPU里闭环Cesium只管渲染。哪怕后面想换一套更复杂的物理模型只要状态纹理的通道布局不变渲染侧就不用大改。2. 关键工具与前置准备2.1 Cesium地形高程数据的获取与网格化在Cesium里获取地形高程最省事的方式是使用Cesium.sampleTerrainMostDetailed。这个方法可以在指定经纬度列表上返回高程值前提是地形Provider支持量算通常使用Cesium World Terrain或自建的TIN地形。工程里的做法是先选定一块矩形区域比如中心点经纬度和尺寸1000米乘800米按目标分辨率均匀生成采样点网格然后调用sampleTerrainMostDetailed拿高程。把这个过程封装成getHeightGrid函数输入中心点、宽高、行列数输出一个Float32Array。async function fetchTerrainGrid(center, width, height, cols, rows) { const positions []; for (let r 0; r rows; r) { for (let c 0; c cols; c) { // 经纬度偏移换算到米再转经纬度小范围够用 const lon center.lon (c / (cols - 1) - 0.5) * width / 111320; const lat center.lat (r / (rows - 1) - 0.5) * height / 110540; positions.push(Cesium.Cartographic.fromDegrees(lon, lat)); } } const sampled await Cesium.sampleTerrainMostDetailed(terrainProvider, positions); return new Float32Array(sampled.map(p p.height)); }因为地球表面是曲面直接拿经纬度网格做局部模拟会引起格距不均。所以我先把经纬度投影到以中心点为原点的局部ENU坐标系算出每个格子的米制坐标再反过来做双线性插值得到规则网格。后面着色器里要用到的cellSize也在这个阶段计算好单位是米每格子。这个值直接影响坡度计算不能拍脑袋乱填。如果项目里还叠加了MVT格式的矢量切片比如道路、居民地这些矢量要素可以转成掩码栅格。初始化状态纹理时把道路和建筑所在格子的A通道设成0这样泥石流模拟就不会把基础设施“冲塌”了。这个细节对数字孪生场景非常重要。2.2 WebGL环境中做GPU计算的可行方案Cesium底层是WebGL我们做GPU计算可以沿用同一个GL上下文也可以另开一个离屏canvas。我建议使用独立的离屏上下文避免和Cesium的主渲染器抢状态。不然你在计算时突然改了当前绑定帧缓冲Cesium下一帧可能就渲染黑屏了。做GPGPU的核心是渲染到纹理。WebGL1时代需要依赖OES_texture_float和WEBGL_color_buffer_float扩展才可以把高度场写进浮点纹理。WebGL2之后情况好了很多可以直接用EXT_color_buffer_float支持浮点渲染目标。代码里我通常这样初始化const gl canvas.getContext(webgl2, { antialias: false }); const ext gl.getExtension(EXT_color_buffer_float); if (!ext) { console.warn(当前浏览器不支持渲染到浮点纹理会降级到半浮点); }如果目标用户机器太老就自动降级到gl.RGBA16F。半浮点的精度略低视觉上影响不大。但有一个坑半浮点在高度差很小的平地上可能看出台阶所以地形网格要记得减去整体基准高度把高程变化压到合理范围不要让GPU纹理里的数值动辄几万米。2.3 工程搭线与Cesium版本选择Cesium我用的版本是1.117功能上足够没有做源码修改。开发时用Vite 5搭的Cesium通过npm install cesium引入。要注意Cesium默认会加载自己的CSS如果不需要左下角那些控件把baseLayerPicker、animation、timeline等选项全部关掉能省不少资源。工程结构大概是这样erosion-demo/ index.html src/ cesium-main.js // Cesium初始化与Primitive管理 terrain-grid.js // 地形采样与网格生成 gpu-erode.js // WebGL2离屏上下文与着色器 shaders/ erosion.frag // 侵蚀迭代着色器 normal.frag // 法线计算自定义地形更新用Cesium.Geometry加Cesium.Primitive实现。这样不用侵入Cesium源码也方便后期升级版本。网格的顶点数在256x256以内时JS数组更新完全能承担超过512x512再考虑用BufferSubData优化。3. 核心实现把侵蚀方程搬进着色器3.1 状态纹理与双缓冲循环设计GPU计算的“状态”全部放在纹理里。我这里用了RGBA32F纹理四个通道含义如下R通道地表高度HG通道水或泥浆量WB通道沉积物浓度SA通道可侵蚀掩码1表示可侵蚀0表示固定掩码通道在工程里非常有用。道路、房屋基底、护坡这些地方不能因为模拟被冲掉初始化时在对应格子设为0。算法更新高度时乘一下A通道固定区域就不会变化。这一步很关键否则你做的数字孪生场景里泥石流会把道路和建筑物地基都冲出一个坑业务方一看就觉得不专业。双缓冲循环是GPGPU的基础操作。两张尺寸完全一样的纹理交替作为帧缓冲目标let readTex createStateTexture(initialState); let writeTex createStateTexture(initialState); let readFbo createFbo(readTex); let writeFbo createFbo(writeTex);每次迭代绑定writeFbo把readTex作为uniform传入shader绘制一个全屏三角形。绘制完成后交换read和write的引用就完成了一帧迭代。浏览器的一帧里可以连续做多次迭代我通常跑20次。迭代次数越多泥石流侵蚀的传播范围越广但太多次也会让细小噪声被放大所以这个参数要结合地形质量和视觉效果来调。3.2 侵蚀着色器的核心代码侵蚀迭代的着色器不需要太多花哨技巧核心就是“采样邻域、计算通量、更新自身状态”。我给出一个简化版GLSL片段省略了水流方向分配只保留“坡度决定侵蚀或沉积”的主干逻辑。代码里的u_stateTex是当前状态纹理u_stateSize是纹理尺寸u_cellSize是每个格子的米数。#version 300 es precision highp float; precision highp sampler2D; in vec2 v_uv; out vec4 outColor; uniform sampler2D u_stateTex; uniform vec2 u_stateSize; uniform float u_cellSize; uniform float u_dt; uniform float u_restAngle; uniform float u_erosionRate; uniform float u_depositionRate; uniform float u_waterInflow; void main() { vec2 px 1.0 / u_stateSize; vec4 self texture(u_stateTex, v_uv); float h self.r; float water self.g; float sed self.b; float fixedMask self.a; float hL texture(u_stateTex, v_uv - vec2(px.x, 0.0)).r; float hR texture(u_stateTex, v_uv vec2(px.x, 0.0)).r; float hD texture(u_stateTex, v_uv - vec2(0.0, px.y)).r; float hU texture(u_stateTex, v_uv vec2(0.0, px.y)).r; float slopeX (hL - hR) / (2.0 * u_cellSize); float slopeY (hD - hU) / (2.0 * u_cellSize); float slope sqrt(slopeX * slopeX slopeY * slopeY); float erode u_erosionRate * water * max(0.0, slope - u_restAngle); float deposit u_depositionRate * sed * max(0.0, u_restAngle - slope); water u_waterInflow; float dh (-erode deposit) * u_dt; h h dh * fixedMask; water max(0.0, water - 0.001); sed max(0.0, sed erode * u_dt - deposit * u_dt); outColor vec4(h, water, sed, fixedMask); }这段代码的重点是slope直接用中心差分近似坡度。严格来说应该取最大下坡梯度但中心差分会带来更平滑的结果泥石流看起来更“糊”视觉上更像泥浆在漫流。fixedMask把固定区域的修改乘掉保证建筑底下不会塌。水流方向分配在这里被简化了water只做简单的增加和蒸发。如果你需要更真实的沟谷发育可以把相邻格子的高度差转换为流量再分配但那会引入更多迭代次数先跑通再优化。3.3 与Cesium地形Mesh结合的具体实现GPU迭代完了高度数据还躺在纹理里。要把纹理变成Cesium里的三维地形有两种做法。第一种是Mesh更新法。创建一个与网格对应的Cesium.Geometry顶点数目是cols * rows索引连接成(cols-1) * (rows-1) * 2个三角形。每次更新时把从GPU读回来的Float32Array填入顶点position数组再更新Geometry。原型阶段强烈建议先用这种因为开发快、理解成本低还能自由控制更新节奏。我实际代码里第一次创建好Primitive后后续更新会直接替换内部Geometry并重建Primitive。网格规模不超过512x512时重建开销完全可以接受如果到了1024x1024建议用gl.bufferSubData只更新缓冲。第二种是自定义TerrainProvider。Cesium的地形Provider可以返回一张Heightmap图像我们把GPU纹理转换后传给Provider理论上可以让Cesium原生地形块也跟着变化。这个方案更贴合Cesium渲染管线但实现复杂要处理瓦片调度和层级细节适合中后期优化。如果你想上线一个正式的产品最终还是要走向这个方向因为原生地形块可以跨级别调度细节也更丰富。3.4 效果增强动态光照、流动纹理与粒子泥石流不是静态模拟视觉上必须动起来。我叠加了三层效果第一层是流动纹理把渐变噪声图平铺到地形Mesh上UV根据模拟出的水流方向偏移让地表看起来有液体在流。Cesium的Primitive支持自定义Material在fragment shader里采样噪声纹理并按时间旋转UV即可。第二层是动态法线高度场更新后法线要跟着变法线可以用相邻高度差计算体现沟壑和堆积感。这样动态光照打在地形上会有凹凸变化而不是一块平板。第三层是粒子流在泥石流主沟道上方撒一些粒子轨迹沿坡度方向运动模拟泥浆团滚动。这三层效果如果跟高度场更新不同步看起来就会很假。我的做法是把模拟结果的更新时间戳暴露出来流动纹理和粒子的生命周期都绑定这个时间戳。比如粒子从沟道起点启程经过1.5秒到达堆积区刚好对应高度场里“侵蚀量先增后降”的阶段视觉和数值就能对上了。4. 实操记录从零跑到流畅演示的全过程4.1 快速原型256x256网格、纯JS我建议先别接Cesium纯WebGL2把侵蚀迭代跑通再叠加上去。这样定位问题简单。第一步生成一个假地形作为初始高度场用两个高斯叠加模拟山脊和沟谷const size 256; const data new Float32Array(size * size * 4); for (let y 0; y size; y) { for (let x 0; x size; x) { const nx x / size, ny y / size; let h 10 * Math.exp(-((nx - 0.3) ** 2 (ny - 0.4) ** 2) / 0.02) 6 * Math.exp(-((nx - 0.7) ** 2 (ny - 0.6) ** 2) / 0.03); const i (y * size x) * 4; data[i] h; data[i 1] 0.0; data[i 2] 0.0; data[i 3] 1.0; } }第二步创建纹理和帧缓冲。纹理用gl.RGBA32F过滤方式必须用NEAREST否则采样时高度值会被插值破坏数值稳定性。function createDataTexture(gl, data, w, h) { const tex gl.createTexture(); gl.bindTexture(gl.TEXTURE_2D, tex); gl.texImage2D(gl.TEXTURE_2D, 0, gl.RGBA32F, w, h, 0, gl.RGBA, gl.FLOAT, data); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE); return tex; }第三步写一个双缓冲迭代函数。每次迭代绑到写入FBO上绘制一个全屏三角形。全屏三角形的顶点缓冲很简单三个顶点覆盖整个裁剪空间v_uv由顶点位置推导。这个套路在WebGL领域很常见网上也有现成封装。4.2 关键参数调试让侵蚀效果既明显又不失真参数调不好GPU再快也白搭。我整理了一张参数清单适合作为初始值参数含义初始值调整方向u_cellSize每个格子代表米数5到10米越小越精细u_dt时间步长0.02太大会爆炸u_restAngle休止角0.2越小越容易侵蚀u_erosionRate侵蚀系数0.02过大会锯齿u_depositionRate沉积系数0.04过小会让泥沙无限迁移u_waterInflow单次水量补充0.01越大沟谷越明显iterations每帧迭代次数20越多变化越连续调参顺序建议先降u_restAngle让侵蚀发生再看u_erosionRate是否把地形啃穿。如果高度发散先调小u_dt再看。这里有一个很有用的调试方式把迭代时的状态纹理可视化输出到一个2D canvas高度归一化成灰度颜色越亮表示越高。这样你可以实时观察侵蚀形态不用等Cesium刷新。水量可以映射为蓝白沉积物映射为橙色。通过切换显示模式肉眼就能看出状态是否合理比一遍遍跑真3D快得多。4.3 从GPU读回高度数据并更新Mesh要更新Cesium里的地形Mesh最直接的办法是gl.readPixels把浮点纹理拷到CPU。注意必须保证当前绑定的帧缓冲是你要读的那一张不然读出来的是别的纹理数据。const pixels new Float32Array(size * size * 4); gl.bindFramebuffer(gl.FRAMEBUFFER, readFbo); gl.readPixels(0, 0, size, size, gl.RGBA, gl.FLOAT, pixels);readPixels是同步操作会阻塞CPU。我的做法是每5帧读一次读回来只取height通道然后更新Mesh顶点数组。如果场景里还有其他业务逻辑可以把读回频率降到每10帧一次视觉上差异不大。读回后要注意行列方向Cesium默认纹理原点在左下数组通常从下到上遍历不对齐的话地形会上下颠倒。我第一次接Mesh时就在这里花了一整天就是行序没对上。如果要进一步提升性能可以不用readPixels而是把纹理保留在GPU里在顶点着色器里直接采样状态纹理作为顶点高度。这个方案能彻底避免CPU和GPU之间的同步但实现上需要处理Cesium的Primitive渲染打包复杂度会高一些。原型阶段先用readPixels把流程跑通后续再优化。4.4 性能验证与调优验证“GPU比CPU快”不能只看直觉。我在工程里加了一个简易性能计数器在Cesium的scene.preUpdate前后记录侵蚀迭代耗时统计一帧总耗时。实测数据如下256x256网格、20次迭代、每次迭代两次纹理切换整体侵蚀耗时约6到8毫秒readPixels一次约2毫秒但我每5帧才读一次平均下来0.4毫秒Mesh更新约1毫秒。总预算控制在15毫秒以内Cesium自身渲染还有余量。如果发现耗时突然升高优先排查这几个方向。一是FBO切换过于频繁可以尝试同时绑定多个颜色附件并行输出二是纹理过滤误用了LINEAR导致驱动做了额外滤波三是半浮点到浮点转换在CPU端做了耗时翻倍四是Cesium的Primitive每次都重建而不是更新缓冲。后面这几个都是实际踩过的问题列出来帮大家避坑。5. 常见问题与避坑实录5.1 常见问题速查表把我和同事在这个项目里遇到的问题汇总成表基本覆盖从“黑屏”到“地形不更新”再到“浏览器崩溃”的路径。现象可能原因解决办法渲染结果是黑屏未启用浮点渲染扩展FBO不完整检查EXT_color_buffer_float改用RGBA16F地形完全没变化读回时机错误Mesh没绑定对应纹理确认读回帧缓冲和迭代写缓冲是同一个FBO侵蚀越跑越离谱u_dt太大导致数值发散把dt降到0.005重试出现大量条纹或噪点纹理过滤用了LINEAR把MIN和MAG过滤都改成NEAREST地物被冲塌掩码通道没设置初始化时把固定区域A值设为0场景在旋转时崩溃资源泄漏FBO或纹理未释放每次重建前调用gl.deleteTexture和gl.deleteFramebuffer高度差表现不明显直接使用绝对高度单精度尾数不够减去区域最低点高度用相对高度5.2 调试技巧与工具GPU计算没有console.log这是初学的第一道坎。我的经验是把中间结果“画”出来。高度归一化为灰度水量映射为蓝白沉积物映射为橙色。通过切换显示模式肉眼判断状态是否合理。还可以在shader里做“断言”比如检测到高度大于某个阈值就输出红色。if (h 100.0 || h -100.0) { outColor vec4(1.0, 0.0, 0.0, 1.0); }这样数值发散时画面上会立刻出现红点定位是哪一步出了问题。浏览器层面Chrome的WebGL调试器可以看FBO、纹理和着色器编译日志。另外gl.getError()返回的错误码也需要熟悉比如CONTEXT_LOST_WEBGL说明上下文丢失通常是因为页面切换或纹理过大。遇到这种情况要重新创建离屏上下文并恢复所有状态纹理。5.3 项目里最值得提醒的3个坑第一个坑是坐标系统。Cesium使用全球坐标系地面的高度动辄几十万米甚至更多。直接把真实高程塞进RGBA32F单精度浮点的精度不足以表达1米以内的变化。解决办法是先减去一个基准面高度让GPU纹理里的高度值保持在0到100米范围内读回后再加回基准。第二个坑是纹理坐标方向。Cesium的Primitive顶点顺序、纹理UV和GPU纹理行序如果在某一环翻转地形就会出现镜面或错位。排查这类问题的快捷方式是在Cesium场景里放一个带经纬度标注的底图然后看腐蚀后的沟谷是否落在预期位置。第三个坑是不要每帧读回。readPixels虽然方便但会打断GPU流水线。如果每帧同步读回性能会掉一半以上。建议固定节流或者干脆把高度数据留在GPU里通过顶点纹理采样更新Mesh这是更推荐的生产级做法。最后说说我的个人体会。这个项目试下来最大的难点不是把侵蚀方程写对而是让数值循环和三维渲染共享同一份状态。Cesium的地球漫游太成熟了大家都习惯往场景里塞模型、塞瓦片但很少有人把GPU当成一个持续的物理推进器。用纹理做状态表用乒乓缓冲做迭代再用渲染管线同步状态这套思路在数字孪生里会越来越常用。后面的扩展我还会继续做把降雨量数据实时注入加入多岩性层让侵蚀更自然再尝试用3DTiles单体化后的建筑数据做更精细的风险分析。如果你也在Cesium里折腾GPU计算欢迎交流这条路上踩过的坑。