Stochastic geomorphological transport for terrain erosion simulation

一句话总结

提出一种并行、随机、基于粒子的”地貌搬运”（geomorphological transport）算法，把水、泥沙、碎屑乃至速度本身都当作沿流线搬运的守恒量来求解，从而在地质时间尺度上做出带动量守恒的大尺度地形侵蚀模拟，首次自然涌现出蜿蜒河、辫状河与三角洲等形态。

研究背景

• 领域现状：大尺度山地地形的物理侵蚀模拟主流依赖 Stream Power Law（SPL），把河流侵蚀写成流量与坡度的函数。它靠”水严格沿最陡坡下流”这一假设，只需对地形做拓扑遍历即可算出流量，因而能用很大的时间步长、模拟百万年尺度的地貌演化，且有 GPU 加速方案（如 Fastflow）。
• 核心痛点：SPL 这类方法为了效率牺牲了物理完整性——它丢掉了动量守恒，速度被强行绑定到地形梯度上。结果是无法表现蜿蜒河、辫状河、三角洲、碎屑流冲沟这些恰恰源于动量相互作用的关键形态；同时基于欧拉网格的离散会带来明显的方向性伪影和强分辨率依赖。
• 本文 idea：作者观察到”动量守恒才是这些特征地貌的驱动力”。于是把侵蚀重新表述为对若干搬运量（流厚、泥沙厚、碎屑厚、速度）的统一搬运问题，并让速度场本身也被搬运以实现动量守恒；再用一套源自基本守恒定律推导的随机粒子格式来高效求解。

方法

整体框架：全流程建立在”地貌搬运”这一核心近似上——因为水/泥沙穿越地形只需几天，而地形演化要上千年，两个时间尺度相差数个量级，所以在侵蚀的时间步内可把搬运视为准静态（局部 \(\partial/\partial t \approx 0\)）且长程（物质一步内即从源头走到域边界）。在此近似下，质量搬运、泥沙/碎屑搬运、动量守恒都能写成同一种线性（或可线性化）守恒律，于是可以用一个统一的求解器处理，最后落到 GPU 上的蒙特卡洛粒子积分。

flowchart LR
  A["初始平地 + 边界条件<br/>(抬升/降雨/反照率)"] --> B["计算侵蚀/沉积速率"]
  B --> C["解动量守恒 → 更新速度场 v"]
  C --> D["解质量/泥沙/碎屑搬运 → 更新各搬运量"]
  D --> E["更新高程 z"]
  E --> F["随机粒子蒙特卡洛积分 Alg.1<br/>沿流线累加上游贡献"]
  F --> B

关键设计：

1. 把侵蚀重写成”依赖速度”的搬运模型：作者回到 Stream Power Law 的剪应力原始推导，让河流侵蚀显式依赖壁面剪应力 \(\tau_f = \tfrac{1}{8} f_D \rho_f \lVert \boldsymbol{v}_f \rVert^2\)，从而侵蚀 \(E_f = k_f \tau_f^{a}\) 直接由速度决定，而不是由坡度经验幂律给出。碎屑流则采用地球物理的力平衡模型，把侵蚀/沉积表达成按流速归一化的形式，并用混合了 Mohr-Coulomb（干颗粒）与 Bingham（黏性）的壁面剪应力刻画从干碎屑到饱和泥石流的连续过渡。滑坡（图形学里常称热侵蚀）作为碎屑流的触发项，在坡度超过休止角 \(\theta\) 时启动。

2. 统一的准静态线性守恒律：所有搬运量都满足一般形式 \(\nabla \cdot (\phi \boldsymbol{v}) = S - R\phi\)，其中 \(S\) 是源、\(R\) 是衰减。水搬运里源是降雨、衰减是蒸散；动量守恒里源是坡度驱动力 \(-g\nabla z\) 加黏性项、衰减来自摩擦，而且速度以 \(\boldsymbol{v}^2\) 形式非线性出现，通过 \(\boldsymbol{v}(t{+}\Delta t)\cdot\boldsymbol{v}(t)\) 线性化处理。关键是速度自己也作为一个被搬运量，这正是动量守恒得以贯穿全局的机制。

3. 上游贡献的叠加积分：求解上把 \((\phi\boldsymbol{v})\) 表达成沿流线的”衰减后上游贡献”之叠加。引入沿流线的累积衰减函数 \(M\) 与衰减比 \(\alpha(y \to x) = M(x)/M(y)\)，把源与衰减解耦；再借助散度定理，把穿过某曲线 \(C\) 的通量写成对上游流域 \(\Omega_C\) 内衰减源项的面积分。该衰减比实际上是这条守恒律的格林函数。

4. 并行随机蒙特卡洛粒子格式：把每个网格单元的出流边作为积分曲线，用蒙特卡洛估计上游积分。核心技巧是样本复用——一个从上游位置发出的粒子沿流线前进，可为它经过的所有下游单元贡献通量，因此采用全域均匀联合采样 \(p(y_i)=1/\lvert\Omega\rvert\) 而非逐单元条件采样。每个粒子边走边按 \(M \leftarrow M\,e^{-\Delta x R/\lVert\boldsymbol{v}\rVert}\) 衰减并累加源项，直到离开域、超时或衰减低于阈值 \(\epsilon\)。该格式无条件稳定，能推广到高维和非规则离散；再配一个时间指数滤波 \(\phi(t{+}\Delta t) \leftarrow \beta\phi(t{+}\Delta t) + (1-\beta)\phi(t)\) 抑制低样本数下的噪声。

实验结果

在锥形汇聚地形上与解析解对比（水搬运，均匀源、无衰减），验证了收敛性：误差随样本数/分辨率比线性下降直到每格 1 个样本，固定采样比时也随分辨率线性收敛。性能上（C++/CUDA，RTX 5070 Ti），运行时在 GPU 线程饱和前基本恒定，之后随样本数线性增长，随分辨率边长线性增长，但与物理尺度无关。下表为不同分辨率与样本数下单次搬运积分的耗时（取自原文表 1）：

样本数 N	512²	1024²	2024²	4096²
2¹⁴	1.1 ms	2.0 ms	3.9 ms	8.6 ms
2¹⁶	1.1 ms	2.1 ms	4.2 ms	11.4 ms
2¹⁷	1.2 ms	3.0 ms	6.2 ms	18.6 ms
2¹⁸	3.0 ms	6.1 ms	12.4 ms	36.1 ms

完整侵蚀模型在 1024² 网格上单次迭代约 51 ms，整场模拟（1000~4000 步）约 1~10 分钟。消融实验表明：去掉动量守恒后蜿蜒河消失、河道退化为无宽度线条；降低沉积率会消除辫状与三角洲；碎屑流去掉动量后沉积只在单一固定坡度堆积。与 SPL 对比显示，本方法消除了方向性伪影、给出更丰富的山体分布，并在多分辨率间保持稳定形态。代价是速度更慢：相比 Fastflow 一次流路由约 2 ms，本方法在同分辨率、2¹⁷ 样本下需约 51 ms（同时求解河流与碎屑流搬运）。

亮点与局限

• 亮点：
  • 从基本质量/动量守恒推导出统一的搬运框架，物理上自洽，且把河流、碎屑流、滑坡乃至外加均匀力（模拟风蚀、类沙丘结构）都纳入同一套方程。
  • 首次在大尺度地质时间侵蚀中让”速度本身被搬运”，从而自然涌现蜿蜒河、辫状河、三角洲、碎屑扇等以往方法做不出的动态形态。
  • 拉格朗日粒子 + 蒙特卡洛格式无条件稳定、易并行、对分辨率稳健，并提供了”样本数换精度”的性能可调旋钮；代码开源。
• 局限：
  • 准静态近似是双刃剑：当侵蚀与搬运时间尺度相当时模型失效，例如单次洪水/滑坡事件（侵蚀太快）或冰川（搬运太慢）。
  • 不建模湖泊/水潭的形成，三角洲只能在高泥沙、低动量的浅水区短暂出现。
  • 动态速度场导致不存在 SPL 意义下的清晰稳态，只有主坡稳定、平缓河道持续摆动的”动态稳态”；且计算成本显著高于 SPL 类方法。

延伸思考

这项工作把图形学地形侵蚀从”坡度驱动”推进到”动量驱动”，本质上是引入了一个针对准静态守恒律的格林函数式蒙特卡洛求解器——这种”沿流线累加衰减上游贡献”的思路，与渲染里沿路径积分的蒙特卡洛估计颇为神似，二者在方差控制、重要性采样、样本复用上的经验或可互相迁移。作者也指出后续方向：改进采样分布、用图归约进一步 GPU 加速、以及扩展到逆向侵蚀（从现有地形反推历史地貌）和沙丘、海岸、落石、地下水等更广的自然现象统一建模。对做程序化地形生成或地学可视化的人来说，”用少量样本快速出稿、加样本再精修”的可调特性，是相对 SPL 很实用的一点。