Cesium 城市风场 CFD：北斗体素网格、LBM 与 GPU 粒子

城市风场案例不是给建筑贴一层流动纹理，而是先把城市表面转成碰撞体素，再在有限求解域内运行 LBM，最后用粒子把速度场画出来。

实现思路导图

flowchart TD
  ROOT[城市风场 CFD]
  ROOT --> INPUT[空间输入]
  INPUT --> AREA[绘制求解区域]
  INPUT --> SAMPLE[采样 3D Tiles 高度]
  ROOT --> GRID[空间离散]
  GRID --> MASK[建筑固体掩膜]
  GRID --> LATTICE[北斗三维点阵]
  ROOT --> SOLVER[LBM 求解]
  SOLVER --> COLLIDE[碰撞与迁移]
  SOLVER --> BOUNDARY[边界与阵风]
  ROOT --> OUTPUT[速度图集与粒子可视化]

从用户区域建立求解域

用户左键绘制多边形，双击结束。程序先得到经纬度包围盒，再等待 3D Tiles 加载并采样表面高度：

const sample = await sampleBBoxHeights(viewer, bbox, {
  sampleRes: params.sampleRes,
  tileset: osmTileset,
  onProgress: (text) => {
    status.value = text;
  },
});

采样值同时包含地形与建筑表面。structureMin 用来判断某个高度突起是否应当视为建筑结构，groundWindow 则控制局部地面估计窗口。

北斗三维网格承担空间离散

案例按指定层级建立三维点阵。区域过大或纹理预算不足时会自动降低层级：

lattice = buildLattice(bbox, params.level, hMin, hMax, {
  maxXY: 200,
  maxZ: 48,
  maxAtlas,
});

随后生成两个掩膜：domainMask 表示多边形内的有效求解列，solid 表示地形和建筑占据的体素。碰撞体素可以用半透明盒或线框显示，便于检查采样是否把建筑识别出来。

LBM 求解与可视化分离

BeidouWindLBM 接收点阵、固体掩膜和入流参数。案例在每个高度层运行 D2Q9，并通过连续性关系得到垂直速度。求解结果进入纹理，粒子再从纹理中采样速度：

solver = new BeidouWindLBM(viewer, {
  lattice,
  solid: mask.solid,
  domain: domainMask,
  inflow: params.inflow,
  omega: params.omega,
  windDeg: params.windDeg,
  particleTexN: params.particleN,
});

粒子数量和实际绘制箭头数分别控制。前者影响速度场采样密度，后者影响屏幕上的几何提交量。

参数修改为什么有不同成本

透明度、步长和切片高度可以直接更新求解器。粒子纹理尺寸、播种网格或风向变化会改变资源布局，需要重建求解器。风场高度变化只重建点阵和掩膜，可以复用已经完成的表面采样。

结果的适用边界

这是浏览器内的交互式近似求解。网格分辨率受 GPU 最大纹理尺寸和实时帧率限制，不能替代工程级 CFD 网格、边界条件和收敛验证。它更适合方案预览、交互演示和局部趋势判断。

截图中的结果应该怎么看

新的截图已经完成一次真实求解，不是初始化画面。青色四边形是用户绘制的有效区域；区域内部密集的绿色、青色短线是粒子在当前速度场中的运动轨迹。左侧统计显示这次计算使用第 24 级北斗网格，实际点阵是 158 × 140 × 32，总计 707840 个格点。

截图关闭了红色碰撞体素，否则建筑占用网格会遮住粒子。关闭的只是体素可视化，solid 掩膜仍然参与 LBM 求解。这个区别很重要：面板里的“显示碰撞网格”控制画面，不控制障碍物是否进入计算。

三套离散空间不能混在一起

实现中同时存在三套网格：

表面采样格是输入，不是最终 CFD 网格。程序先得到 72×72 的表面高度，再用双线性方式上采样到北斗点阵的 nx × ny。如果把两者当成同一分辨率，调整北斗层级时就必须重新请求 3D Tiles，交互成本会明显增加。

三维点阵也不能直接作为一张三维纹理使用，因为 WebGL 设备的三维纹理尺寸和兼容性限制更严格。案例把每个高度层排进二维图集：

const maxAtlas = Math.min(
  4096,
  viewer.scene.context.maximumTextureSize || 4096,
);

lattice = buildLattice(bbox, params.level, hMin, hMax, {
  maxXY: 200,
  maxZ: 48,
  maxAtlas,
});

buildLattice 会检查图集预算。当前层级放不下时自动降级，而不是创建一张超过 GPU 上限的纹理后等待 WebGL 报错。

如何从城市表面得到固体掩膜

sampleBBoxHeights 对包围盒内的经纬度格点执行场景高度采样。Cesium 的 sampleHeightMostDetailed 会同时考虑当前场景里的 3D Tiles 和地形，因此返回值更接近“从地心向外遇到的最高表面”。

采样完成后，程序先估计局部地面，再判断建筑结构：

mask = buildSolidMask(lattice, surfaceHeights, {
  structureMin: params.structureMin,
  groundWindow: params.groundWindow,
});

groundWindow 用邻域低值估计地面基准。表面高程比基准至少高出 structureMin，这一列才按建筑处理。随后从求解域底部填充到表面高度，得到 Uint8Array 固体掩膜。

这个方法能把明显建筑转成障碍物，但它不是建筑语义识别。高架桥、树冠、悬挑结构和多层空洞都会被压成按列占用的体素。文章必须明确这一点，否则读者会误以为模型保留了完整建筑拓扑。

D2Q9 在每个高度层算什么

LBM 不直接求解连续的速度方程，而是在每个格点维护九个离散方向的分布函数。一次迭代分成碰撞和迁移：

碰撞：fᵢ* = fᵢ - ω(fᵢ - fᵢᵉᵠ)
迁移：fᵢ(x + eᵢ, t + 1) = fᵢ*(x, t)

其中 ω 是面板里的松弛频率。它决定分布函数回到平衡态的速度，也间接影响数值黏性。ω 太接近稳定边界时，速度场容易振荡；数值更小则流场更黏、更平滑。

案例没有运行完整 D3Q19，而是在每个高度层执行 D2Q9。水平速度来自各层 LBM，垂直速度再通过连续性近似耦合。这样显著降低纹理通道和计算 pass 数，代价是垂直涡结构只能近似表达。

固体格点使用反弹边界。粒子遇到固体时不会继续穿过建筑，而会按速度纹理和播种规则重置或绕流。

入流、风廓线和阵风怎样进入求解

windDeg 决定入流面及水平方向。inflow 是格子单位下的数值速度，refWind 只负责把无量纲速度标定成图例中的米每秒，两者不是同一个参数。

风廓线使用高度幂律控制近地面和高处的速度差：

u(z) = uref × (z / zref)^α

shear 对应指数 α。值越大，近地面减速越明显。gust 则让入流随时间发生周期或噪声扰动，它改变的是边界条件，不是简单让粒子动画加速。

速度场如何变成粒子轨迹

粒子位置也保存在 ping-pong 纹理中。每帧读取旧位置，在速度图集中做三线性近似采样，再写入新位置。下一帧交换两张纹理：

positionA --读取--> 粒子更新 pass --写入--> positionB
positionB --读取--> 下一帧更新   --写入--> positionA

particleN = 256 表示状态纹理边长为 256，因此总粒子数是 65536。visibleArrows 只限制最终绘制的箭头数。计算粒子和显示箭头分开后，可以保留较密的粒子状态，同时控制屏幕几何量。

life 是粒子寿命，dt 是推进步长，tailLen 是按格子尺度计算的箭头长度。增大 dt 不等于提高真实风速；它只是让粒子每帧跨过更远距离，过大时会穿过薄障碍或产生断裂轨迹。

一次完整执行的状态流

1. startDraw 清理旧求解器和旧体素，进入绘制状态。
2. finishDraw 计算经纬度包围盒并调用 runPipeline。
3. 等待 3D Tiles，就绪后完成场景表面采样。
4. 构建北斗点阵、区域掩膜和固体掩膜。
5. Worker 生成碰撞体素几何，主线程创建 Primitive。
6. BeidouWindLBM 创建速度纹理、求解 pass 和粒子资源。
7. 渲染循环持续迭代速度场并更新粒子。

busy 只防止用户重复启动流程。surfaceSample 则是可复用的采样快照。修改求解域高度时可以复用它；重新绘制区域或修改表面采样分辨率才需要重新读取城市表面。

销毁为什么要分顺序

切换案例时先终止 Worker，再销毁求解器和事件处理器，然后移除体素、3D Tiles，最后销毁 Viewer。求解器持有 framebuffer、纹理和 Primitive；如果先销毁 WebGL 上下文，再让这些对象执行清理，就可能出现已失效资源访问。

实现主线是：3D Tiles 表面采样 → 北斗三维点阵 → 固体与区域掩膜 → LBM 速度场 → 粒子和切片可视化。

从零复现：城市风场 CFD 教程

案例源码位于 src/demos/cesiumdemos/city-wind-cfd/。建议先用小网格跑通，再增加表面采样格和求解层级。截图中的高密度参数不适合作为第一组调试参数。

第 1 步：让用户定义求解域

使用 ScreenSpaceEventHandler 收集多边形顶点，双击结束。计算经纬包围盒后在中心建立 ENU 坐标系，后面的网格尺寸、建筑高度和风速都使用米制局部坐标。

const frame = Cesium.Transforms.eastNorthUpToFixedFrame(center);
const inverseFrame = Cesium.Matrix4.inverseTransformation(
  frame,
  new Cesium.Matrix4(),
);

完成标志：把用户顶点转到 ENU 后，区域中心应接近 (0, 0, 0)，东西和南北长度应与量测距离一致。

第 2 步：等待并采样建筑表面

先等待目标 3D Tiles 可见，再按规则网格生成屏幕或世界射线，用 scene.pickFromRayMostDetailed、pickPosition 等方式取得表面高度。没有命中的点回退到地形高度。

采样结果保存为规则二维高度场，不要只保存散点。后续结构体素判定需要按格点快速索引，而不是每个格子重新拾取一次场景。

第 3 步：建立固体掩膜

对每个水平格点比较体素中心高度与采样表面高度。结构实体模式把表面以下的层标记为固体；表面壳层模式只保留靠近表面的若干层。

solid[index3D(x, y, z)] = cellHeight <= sampledSurface[x + y * nx];

先打开“显示碰撞网格”。若建筑区域没有红色体素，问题在采样或坐标映射；此时不要继续调 LBM。

第 4 步：初始化 D2Q9 分布函数

每个流体格点保存九个方向分布 f0...f8。初始速度设为入口风速，密度设为 1，再用平衡分布函数填充九个通道。二维切片的每一层独立求解，层间通过入口剖面和可视化组合。

feq_i = w_i · ρ · [1 + 3(e_i·u) + 4.5(e_i·u)² - 1.5(u·u)]

松弛时间必须大于 0.5。过于接近 0.5 会降低数值黏性，但更容易在建筑尖角附近发散。

第 5 步：实现碰撞、迁移和边界

一次迭代先由九个分布求密度与速度，再向平衡态松弛；随后按离散方向把结果迁移到相邻格点。固体格点执行反弹边界，把进入墙面的分布反向返回。

for (let iteration = 0; iteration < stepsPerFrame; iteration += 1) {
  collide(readAtlas, collisionAtlas);
  stream(collisionAtlas, writeAtlas, solidMask);
  swapAtlases();
}

入口写入风向和风速，出口使用零梯度近似，上下边界按案例需要设为滑移或无滑移。完成标志：没有建筑时速度场接近均匀；加入单个方盒后，迎风面减速、两侧加速、背风面形成尾流。

第 6 步：用粒子解释速度场

LBM 输出速度图集后，在 GPU 中推进 3D 粒子。粒子先根据高度选择图集层，再采样水平速度；进入固体或离开区域时重生。颜色由速度大小映射到色带。

截图中的青色框是求解域，彩色点线是速度粒子。只有两者同时出现且统计显示已完成迭代，才算得到可解释结果。

第 7 步：管理异步状态

把状态明确分为 idle → sampling → building-mask → solving → ready。每次重新绘制区域时递增 requestId；Worker 或采样回调结束时，只允许当前 id 提交。

取消操作需要终止 Worker、丢弃待提交结果并恢复面板状态。只隐藏 Primitive 会留下仍在运行的求解任务。

验收与排错

• 风场穿过建筑：检查固体掩膜纹理层和 LBM 图集层是否使用同一索引。
• 速度突然变成 NaN：降低松弛强度、时间步或入口风速，检查密度是否为零。
• 结果上下翻转：检查图集行号到高度层的换算。
• 粒子很多但看不出流向：先降低寿命并增加拖尾，再检查速度色标范围。