利用基于游戏引擎的SPH-DEM方法显著提升泥石流快速模拟效果

《Computers & Geosciences》：Enhanced rapid simulation of debris flow using game engine-based SPH-DEM methodology

　　袁-Hung Damiel Chiu | Tzu-Yin Kasha Chen | Hervé Capart
**摘要**
在山区，碎屑流和河流搬运作用常常在陡峭的支流与主干河交汇处形成圆锥形沉积体。这些碎屑和冲积扇通常表现出明显的坡度或海拔与距扇顶距离之间的关系。然而，当扇形受到山谷地形限制时，这种关系就变得更加难以描述，因为最陡峭的坡度路径需要绕过障碍物弯曲。在之前的研究中，我们开发了一种正向建模方法，根据假设的海拔-距离剖面来重建这些复杂的扇形表面。在本文中，我们补充了一种逆向方法：一种系统且快速的方法，可以从观测到的扇形地形或其轮廓曲线推断出海拔-距离剖面。该方法利用可见性多边形来计算从扇顶到所有观测点的最短路径距离，并考虑了障碍物的影响。然后对处理后的海拔-距离数据对进行二次函数拟合。通过对由正向算法生成的合成扇形表面进行测试，我们验证了逆向方法能够精确恢复假设的海拔-距离曲线，误差在可接受范围内。接着，我们从意大利和台湾的实际扇形表面中提取出海拔与最短路径距离之间的经验关系，并证明它们比以往研究中假设的海拔与直线距离关系更能准确描述扇形地形。

**引言**
在支流交汇处和山区前沿，碎屑流和河流搬运作用常形成圆锥形沉积体。这些碎屑和冲积扇不仅具有重要的沉积学价值（Harvey, 2012; Nicovich et al., 2023），同时也对人类和基础设施构成威胁（Davies and McSaveney, 2008; Calvo and Savi, 2009; Khan et al., 2013）。扇形地貌受到主导形成过程的影响显著：以沉积重力流为主的碎屑流扇形通常更陡峭且崎岖，而由流体重力流形成的冲积扇形则坡度较缓且表面更平滑（Bertrand et al., 2013）。无论形成机制如何，由于偶发的碎屑流和洪水脉冲，扇形可能迅速堆积（Blair and McPherson, 2009; Hsieh and Capart, 2013; de Haas et al., 2018），因此预测其堆积高度和受影响范围至关重要。

**划定危险区域的方法**
从简单的行程距离估算（Rickenmann, 2005）到利用数值流体模拟的详细分析（例如Hsu et al., 2010; Wu et al., 2013; Yanagisaki et al., 2016; Shiu et al., 2023）都有应用。这些模型能够处理复杂地形，并详细计算流动过程中的速度和深度。然而，这类计算资源消耗大，且难以仅凭有限的地形数据进行校准。我们提出了一种基于几何的简化方法，可以直接从扇形地形数据进行校准。该方法假设流动沿最陡路径进行，并产生具有稳定距离-海拔关系的沉积体。沉积体的几何形态取决于其是否受地形限制（Stanistreet and McCarthy, 1993; Al-Farraj and Harvey, 2005）：不受限制的扇形会在开放地形上自由扩展，而受限制的扇形则会受到周围地形的约束。

**方法细节**
1. **不受限制的扇形**：流动沿从扇顶辐射出的无障碍直线路径进行，其海拔与径向距离的关系唯一决定了沉积体的几何形态。这种关系通常是向上凹的，即坡度随距离增加而减小（Bull, 1977; Bowman, 2019）。这种关系可以通过将海拔与距扇顶的直线距离绘制出来很容易确定（Bull, 1964; Nicovich et al., 2023; Iacobucci et al., 2024），随后可以用指数或二次函数来近似实际剖面（Troeh, 1965; Rice and Church, 2001）。基于水力或扩散原理的理论模型 también 提供了类似的剖面形状（Parker et al., 1998; De Chant et al., 1999; Densmore et al., 2007; De Chant et al., 2021）。
2. **受限制的扇形**：流动仍沿最陡路径进行，但路径必须绕过地形障碍物弯曲，因此沉积体不再呈辐射对称。尽管如此，只要沿弯曲路径测量距离，仍然可以定义海拔-距离剖面。在Chen和Capart（2022）的研究中，我们解决了正向问题，即使用假设的海拔-距离剖面在不规则地形上重建扇形表面，结果与实际测量结果吻合良好。然而，从测量数据中反推海拔-距离剖面的问题尚未解决。本文探讨了逆向问题，即利用观测到的扇形地形或其轮廓来推导其背后的海拔-距离关系。

**研究目标**
本研究旨在建立一种从复杂扇形地形中提取海拔-距离剖面的通用方法。为此，我们采用了一种考虑障碍物的最短路径算法，以在地形受限的情况下估算流动路径长度。

**假设与正向模型**
图1展示了我们方法的物理和几何假设。首先考虑一个理想化的无障碍沙堆（图1a）：沙粒从顶点A沿最陡路径流动，形成倾角为θ的圆锥形沉积体。根据这一倾角、初始地形（这里是平面）以及给定的顶点位置和海拔，我们可以推导出沙堆的地形：
$$ z(x,y) = z_A \cdot \tan\theta \cdot s(x,y) $$

**扇形表面重建测试**
本节测试了该方法从部分测量数据重建完整扇形表面的能力。我们选取了一个数据较为完整的沉积事件作为测试对象，从而能够严格验证正向-逆向工作流程的完整性。
所选案例是2021年8月7日台湾Laonong河与Yu-Shui支流交汇处的扇形堆积事件，当时一次碎屑流沉积了约300万立方米的物质。

**多扇形应用**
本节通过将最短路径距离方法应用于多个扇形来评估其有效性。我们分析了意大利Valtellina Valley地区的25个扇形以及台湾50个大面积（>2平方公里）的冲积扇形。之后，我们对比三种剖面提取方法的准确性，并提供实用指南。

**形态测量分析**
本节对完整扇形数据集进行了详细的形态测量分析，使用验证过的最短路径距离方法研究了基本形态参数（如扇形面积、坡度和剖面凹度）之间的定量关系。这一系统分析旨在通过将形态特征与形成过程和终端边界条件等控制因素联系起来，为冲积扇的形态学提供新的见解。

**结论**
我们开发了一种新的逆向方法，通过最短路径距离分析来表征冲积扇的形态。该方法基于扇形边界和顶点生成SPM（Surface Point Model），在直接距离测量无法奏效的情况下仍能有效捕捉地形约束。验证工作包括理论研究和实际案例分析，证明了该方法与正向方法的可逆性（Chen and Capart, 2022），并且适用于不同规模的扇形。

**作者贡献声明**
- 袁-Hung Damiel Chiu：撰写原稿、可视化、验证、软件开发、方法论设计、数据分析、形式化分析、数据整理、概念构建
- Tzu-Yin Kasha Chen：审稿与编辑、验证、监督
- Hervé Capart：审稿与编辑、可视化、验证、监督、方法论设计、资金筹措、概念构建

**计算机代码可用性**
- 库名：fan-inverse-forward-methods
- 许可证：LGPL-3.0
- 首次可用年份：2025年
- 所需硬件：任何标准个人电脑（在Windows系统上测试通过）
- 所需软件：Python 3.11及以上版本，以及GDAL、inpoly-cython、Matplotlib、NumPy、SciPy、Shapely、scikit-image、tqdm和修改版的VisiLibity1
- 程序语言：Python
- 程序大小：共计48 MB（其中GDAL轮包占42 MB）
- 源代码链接：[GitHub](https://github.com/damiel-hub/fan-inverse-forward-methods)

**关于生成式AI和AI辅助技术的使用说明**
在准备本文期间，第一作者使用了Gemini、ChatGPT和Claude工具来辅助初稿的语法和句子可读性，帮助将MATLAB代码转换为Python代码，并进行代码格式化和注释。作者们对内容进行了全面审查和编辑，对最终出版物的内容负全责。

**利益冲突声明**
Herve Capart声明获得了国家科学技术委员会的支持。其他作者声明没有可能影响本文研究的财务利益或个人关系。

**致谢**
该方法的早期开发是在Y.-H. D. Chiu和H. Capart访问鲁汶大学（UCLouvain）的Sandra Soares-Fraz?o教授和Fugro Engineers公司的B. Spinewine博士的研究期间完成的。我们还感谢本研究中使用到的关键算法的开发者。可见性多边形算法使用了VisiLibity1的Python绑定实现（Obermeyer and Contributors, 2008），并进行了修改以记录子扇顶信息。等高线提取功能也是基于VisiLibity1实现的。