郑顾平，张 萌，邵绪强，伍朝辉

真实感水力侵蚀地形建模与绘制

郑顾平1，张 萌1，邵绪强1，伍朝辉2

(1. 华北电力大学控制与计算机工程学院，河北 保定 071003；2. 交通运输部科学研究院基础研究创新中心，北京 100019)

水力侵蚀地表是现实生活及虚拟场景中常见的一种自然现象。为研究水力对地表的侵蚀现象，提出一种采用光滑粒子流体动力学(SPH)模拟流体对地形表面的湿润以及侵蚀过程的真实感仿真模型。通过对自然界中真实溅蚀现象进行合理简化描述，建立一种包含3个阶段的地形水力溅蚀模型：首先，利用弹丸冲击平板靶材变形模型模拟水流冲击地表变形效果；其次，构建基于达西定律的水流渗透量计算方法以确保地表吸水与水分扩散的真实性；最后，提出将地表湿润度因素对地表侵蚀的影响考虑在内的新型计算模型，实现与冲蚀相融合模拟积水继续运动对地表冲蚀作用。实验结果表明，该模型能够真实地模拟出不同方向的水流对地表的侵蚀效果以及侵蚀过程中地表吸水与湿润扩散现象。

水力侵蚀；光滑粒子流体动力学；溅蚀；达西定律；湿润扩散

水力侵蚀地形现象时常发生在大自然中，大到泥石流、山洪，小到滴水穿石，都是流体和固体发生交互产生的结果。水力侵蚀地形真实感绘制在防治类似泥石流等自然灾害、游戏与影视特效等多个方面有着广泛的应用。随着计算机性能的飞速提升，基于物理的侵蚀现象动画模拟技术取得很大进展。但是如何在保证高效计算速度的前提下更为真实的模拟出侵蚀现象一直是这些年来的研究热点。

水力侵蚀在地质学与物理学也被广泛研究，多数研究者将重点集中于水力对地形土壤层结构与土壤流失量的影响或者是对生态环境的影响的评估与预测[1-3]。水力侵蚀现象的研究对象可以分为流体和固体。在计算机图形学领域，基于物理的流体模拟方法可以分为欧拉网格法[4-5]和拉格朗日粒子法[6-8]。由于欧拉网格法是在规则的网格上进行数值求解，所以遇到较大的形变时计算效率很低并且在模拟效果上容易失真。相比而言，拉格朗日粒子法从研究流体的各个微团物理属性(密度、压强、速度等)着手，更容易模拟出如水花四溅等细节，其中光滑粒子流体动力学(smoothed particle hydrodynamics，SPH)方法[9]在流体模拟中应用最广泛。

水力侵蚀过程中固体和流体交互伴随着碰撞、湿润扩散以及侵蚀现象的产生。在固液交互的研究中，柳有权等[10]利用GPU加速模拟出气流在迷宫和城市中扩散的场景，但由于内存有限并不能进行大规模问题求解；NARAIN等[11]实现了固液之间的双向交互，其仅将固体作为流体的运动边界，并未考虑固液物理性质的变化；LENAERTS和PHILIP[12]提出了一种多孔介质和流体的耦合方法，实现了干旱的土壤被雨水浸润成泥浆或沙结构被海浪侵蚀的现象；常元章等[13]将Navier-Stokes (N-S)方程添加弹性应力项实现了牛顿流体与粘弹性流体的统一模拟；SHAO等[14]通过将PCISPH和LSM方法进行耦合逼真的模拟出速度差异较大的流体与可变形固体之间的交互；TAMPUBOLON等[15]扩展了基于SPH方法的多相流体框架，包括可变形固体、颗粒材料和多种流体。

对于侵蚀现象的研究，ŠT’AVA等[16]提出了以水力为主要影响地表变化因素的侵蚀交互模拟方法，且将地表划分成多层土壤层，在重力作用下模拟了滑移侵蚀的效果；KRIŠTOF等[17]通过定义SPH粒子携带固体粒子的体积分数表示固液两相流的状态，并且提出供体-受体模型以求解地表土壤的溶解量，同时模拟溶解的物质由于水的流动而在不同的地方形成沉积的现象；MEI等[18]采用GPU加速计算腐蚀的方法模拟了大规模的地形侵蚀现象，在保证模拟逼真度的前提下提升了求解的效率；CORDONNIER等[19]从生态环境因素考虑，引入一种支持侵蚀和植被双向交互的框架，能够模拟出植被、砂石等对侵蚀现象的影响。

溅蚀是水力侵蚀的一种，与其他侵蚀不同的是溅蚀现象由于具有较大动能的流体与地表发生碰撞造成地表变化。当水滴从一定的高度落地，动能作用于地表土壤，导致土壤结构破坏，表层土壤孔隙减少或者堵塞，土壤团粒被分散、剥离，并使部分土粒沿抛物线的运动轨迹向四周飞溅。为了方便在计算机动画中模拟此类现象，本文忽略侵蚀对土壤层结构的影响，以SPH方法为基础，提出了包含3个阶段的近似简化模型：首先，水滴溅落到未湿润的地表上，由于地表还未来得及吸收液体，细小土颗粒被水滴的动能冲撞开[20]；其次，随着降水量的增大，地表持续吸收水分，固液混成的泥浆受水滴冲击飞溅[21]；最后，表面粒子吸水到达饱和，产生积水，积水的流动对地表造成冲蚀的效果。与此同时利用一种流体在多孔固体颗粒渗透力学的方法来模拟地表吸水湿润以及湿润的扩散过程，并且在侵蚀的过程中考虑地表湿润程度对地表侵蚀的影响。通过分析水流冲撞地面各个阶段的成因，将溅蚀与冲蚀相融合来模拟不同阶段地形湿润、变形以及湿润的扩散现象。

1 基于物理的固液建模

1.1 SPH方法流体模型

SPH方法是拉格朗日数值求解方法的一种，其将模拟的流体离散成粒子，每个粒子携带质量、位置、速度、密度等属性，通过研究各个粒子团的运动来研究流体整体的运动。流体的运动状态主要由N-S方程求解得到，N-S方程包括质量守恒方程与动量守恒方程，即

由式(1)和式(2)简化可得

其中，m，分别为粒子的质量与密度；(–,)为径向对称的核函数；为光滑核半径。属性()的梯度算子可由核函数(–,)的梯度算子求得

同理可求得属性()的拉普拉斯算子

1.2 多孔固体颗粒模型

由于地面的土壤缝隙间存在空气，水滴浸润时先填满土壤中空气所占的空间，达到一定饱和度后产生积水。为了模拟这种吸水与积水的过程，本文利用简化多孔固体颗粒模型代替土壤颗粒，每个粒子属性有质量、饱和度等，即

对于地形的模拟采用了高度图方法，即一张二维灰度图，每个像素的灰度值表示该点的高度，然后将其映射到三维空间中，即可得到地形。为了保证模拟效果并提高计算速度，需对地形表面进行多孔固体颗粒采样。如图1所示，将地形投影到平面的一个二维网格，在每个小正方形网格内均匀采样5个多孔固体颗粒。

图1 地表多孔固体颗粒采样图

流体刚开始接触多孔固体颗粒时产生浸润，同时颗粒中的水分向周围扩散，这种现象从微观上看是流体小分子在多孔固体颗粒之间传输的过程。DEAR[22]从大量的实验统计与总结中得出描述这种渗流力学现象的达西定律

其中，为渗透系数；为距离对渗水量的影响权重。多孔固体颗粒的毛细势P与饱和度呈正相关，即

其中，c为大于0的毛细系数。对于任一多孔固体颗粒而言，当其与周围颗粒毛细势不相同时，高毛细势和低毛细势产生压力差，迫使毛细势高的多孔固体颗粒中的液体向周围毛细势低的多孔固体颗粒中渗透。因此利用这种模型可以宏观的模拟出整个湿润的扩散现象，水分子在多孔固体颗粒间扩散如图2(b)所示。

当降水量达到一定程度，从地表脱落的沙土与流体混合形成泥浆。泥浆是一种典型的固液两相流[26-27]，对于泥浆的密度mix、粘性系数mix等物理性质由泥浆的浓度所决定

其中，为泥浆浓度(与上述饱和度不同的是，饱和度是多孔粒子吸水量，而浓度是液体粒子携带固体量)；，分别为流体粒子的密度和粘性系数；，分别为多孔固体颗粒的密度和粘性系数。

2 侵蚀模型

与已有对侵蚀的研究相比，本文扩展了侵蚀的种类。在流体对地形侵蚀阶段，基于KRIŠTOF等[17]的冲蚀模型，提出了一种包含溅蚀过程的侵蚀模型。同时引入达西定律保证水分在地表扩散物理过程的真实性，并且提出的侵蚀模型计算方法能够描述湿润因素对地表变形影响。

2.1 溅蚀

地表变形阶段可采用弹丸冲击平面靶材变形模型进行模拟。如图3所示，当一定速度的弹丸撞击平面靶材时，由于冲击力的作用平面靶材发生形变，形变的宽度a和深度d与速度之间关系的计算方法为[28]

其中，为弹丸的半径；为平面靶材的屈服强度。对于溅蚀的某个区域地表变形与其周围多孔固体颗粒的湿润程度有关，地表越湿润越容易发生形变，所以地表的屈服强度会随着湿润值的增大而减小，基于此性质提出了描述屈服强度与地表多孔固体颗粒饱和度之间相互关系的计算模型

溅起或者卷入水中的沙土与液体形成泥浆，其密度与粘滞力由式(10)和式(11)求得，同时对于泥浆粒子，若其与周围粒子有浓度差，会产生浓度扩散现象，该现象可通过以下模型进行求解

其中，为扩散系数。由于泥浆粒子中的固体微颗粒所受的力主要是重力，所以不能忽略重力对沉淀物沉降作用的影响，文中采用供体-受体体系[17]来描述重力影响下的沉降现象。如图4所示，对于每个浆体粒子而言，其与周围粒子之间的供受体关系由当前泥浆中固体微颗粒沉降速度与粒子间距离夹角决定。若·≥0，即粒子的浓度向粒子方向沉降，反之则粒子向粒子方向沉降。重力影响下SPH形式的浓度扩散式为

图3 弹丸冲击平面靶材变形图

泥浆与融入积水中的固体颗粒会在沉降过程中附着在地表形成堆积。对于堆积量的求解可由供体-受体体系分析得到，在此仅仅将地表处多孔固体颗粒作为受体，泥浆作为供体。堆积是由供体浓度向受体传递产生的，即计算·<0的所有项。堆积量变化率与浓度变化率之间的关系式为

地表多孔固体颗粒吸水达到饱和时产生积水，积水沿着地表流动的动能相对于从空中降落时减少了很多，这一阶段的侵蚀主要是冲蚀。

2.2 冲蚀

冲蚀效果是由流体对地表的剪切力造成的。当流体与地形表面发生接触，地表处多孔固体颗粒受到半径为范围内的流体粒子剪切力的作用，该力将固体颗粒卷入流体中形成泥浆。剪切力的计算方法为

其中，为剪切力的大小；K为剪切力常数；为剪切率，即

其中，K为侵蚀常数；τ为标准剪切力大小。故式(22)可以求得边界处的多孔固体颗粒质量的减少

其中，L为多孔固体粒子之间的距离。

其中，A为小网格的面积。

3 实验结果及分析

基于冲蚀和溅蚀的融合模型，本文实现了一系列的侵蚀现象的真实感绘制。实验硬件配置为：Inter(R) Core(TM) i5-7400 3.00 GHz CPU，8.00 GB内存，NVIDIA GeForce GT 720显卡。编程环境为：Windows 10，Visual Studio 2013，Unity3D 2018。流体浸润及侵蚀地形过程如图5所示，绘制的过程中，先利用Marching Cube算法提取出流体表面的三角形网格，然后利用Unity3D对提取的表面进行渲染。实验中重要参数设置见表1。

图5 浸润及侵蚀算法流程图

表1 重要参数取值表

图6描述了水分浸润材料以及逐渐被吸收的过程，实验中不考虑材料变形，从第100帧开始，每隔100帧提取一张效果图。从序列图中可以看出，湿润从中间向四周扩散，靠近中心的位置颜色深，边缘较浅，颜色的深浅反映不同区域水分的饱和度。当水分完全被材料吸收并扩散开形成第1 200帧的效果，效果图也证实本文计算方法中引入达西定律能够保证湿润及扩散的真实物理效果。

不同屈服强度下地表的变形效果如图7所示。当σ=0时，地表干燥，冲击形成的水坑较浅 (图7(a))；当σ=1时，地表湿润达到饱和，冲击形成的水坑深(图7(c))。本文计算方法效果如图7(b)所示，坑的深浅介于前两者之间，证明该方法可描述出湿润对地表变形的影响。

图8(a)和图9(a)展示了水流刚开始接触地表且湿润地表现象，从图8(b)与图9(b)可以看出水流对地表的侵蚀与地表吸水及湿润扩散的效果，图8(c)和图9(c)是地表完全干涸之后形成的水洼。由此可以看出本文中溅蚀模型具有较好的模拟效果。图8和图9中水流垂直冲撞地表时，其侵蚀现象主要是溅蚀，表2前2行数据也印证了本文溅蚀模型的可靠性。

图6 水分扩散序列图

图7 不同屈服强度形成水坑深度图

图8 单股水流垂直方向侵蚀地表模拟

图9 双股水流垂直侵蚀地表模拟

表2 实验水流速度及溅蚀与冲蚀造成地表变形量

单股水流和双股水流从水平方向流出与地表交互的效果如图10和图11所示。图10(a)和图11(a)为水流湿润流经的地表图；从图10(b)和图11(b)可以直观看出湿润扩散与侵蚀现象；当地表完全干涸后形成的水沟如图10(c)和 图11(c)所示，对比图中的水沟可以发现，水流刚开始接触地表时，由于具有一定的垂直初速度会有溅蚀的现象，水流继续往前运动时，冲蚀效果起主要作用。图11(c)2股流体发生了冲撞，所以侵蚀形成的水沟略短于图10(c)，但由于 图11(c)中冲撞区域水流较多，水沟较图10(c)中相比较深。同时表2中后2行数据也证实了本文中模型较好的将溅蚀与冲蚀现象应用到侵蚀现象的模拟中。

图10 单股水流水平方向侵蚀地表模拟

图11 双股水流水平侵蚀地表模拟

图12展示了字母“JOG”湿润并侵蚀地表，图12(b)和(c)均能看到较为明显的字母侵蚀痕迹。图13(a)水流与书本发生碰撞顺流而下，流经区域呈扇形，在书本的边缘形成3条分支，图13(c)可以看出，水分干涸后呈现的3条浅沟壑。在图14(a)中，水流冲下与水管发生碰撞，由于水流初始位置偏向水管左边，可以看出水管左侧流水量多于右侧，图14(c)能够清晰看出水管左侧侵蚀痕迹比水管右侧明显。

图12 字母“JOG”湿润及侵蚀现象

图13 水流与书本碰撞侵蚀地表模拟

图14 水流与水管碰撞侵蚀地表模拟

4 结束语

本文通过简化自然界中溅蚀的过程提出了一种溅蚀模型来模拟流体侵蚀地表的现象，同时利用多孔固体颗粒来模拟地表湿润的扩散，并考虑地表湿润度对地表变化的影响，通过将提出的溅蚀模型与冲蚀模型相融合实现了水流对地表真实感的侵蚀现象模拟。

在未来的工作中希望能够从以下几点继续改进：在湿润扩散过程中，不应该仅仅是在地表扩散，水流渗透到地表以下对不同土壤层产生影响，对于湿润的地表还应该考虑温度对水分蒸发的影响；此外对于固液两相流本文并未将固体粒子单独模拟，而是用浓度来表达固液混合状态；最后在碰撞检测及响应方面精度与速度的提高也是未来的研究热点之一。

[1] 唐栋, 李典庆, 周创兵, 等. 考虑前期降雨过程的边坡稳定性分析[J]. 岩土力学, 2013, 34(11): 3239-3248.

[2] 吴宏伟, 陈守义, 庞宇威. 雨水入渗对非饱和土坡稳定性影响的参数研究[J]. 岩土力学, 1999, 20(1): 2-15.

[3] 尹占娥, 许世远, 殷杰, 等. 基于小尺度的城市暴雨内涝灾害情景模拟与风险评估[J]. 地理学报, 2010, 65(5): 553-562.

[4] FOSTER N, METAXAS D. Realistic animation of liquids[J]. Graphical Models and Image Processing, 1996, 58(5): 471-483.

[5] STAM J. Stable fluids[EB/OL]. [2019-07-21]. http://graphics.cs.cmu.edu/courses/15-869-F08/lec/10/ns.pdf.

[6] REEVES W T. Particle systems a technique for modeling a class of fuzzy objects[J]. Computer Graphics, 1983, 17(3): 359-376.

[7] TAKAHASHI T, FUJII H, KUNIMATSU A, et al. Realistic animation of fluid with splash and foam[J]. Computer Graphics Forum, 2003, 22(3): 391-400.

[8] STAM J, FIUME E. Depicting fire and other gaseous phenomena using diffusion processes[C]//Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Techniques-SIGGRAPH’95. New York: ACM Press, 1995: 129-136.

[9] MÜLLER M, CHARYPAR D, GROSS M. Particle-based fluid simulation for interactive applications[C]//Proceedings of the 2003 ACM SIGGRAPH, Eurographics Symposium on Computer Animation. New York: ACM Press, 2003: 154-159.

[10] 柳有权, 刘学慧, 吴恩华. 基于GPU带有复杂边界的三维实时流体模拟[J]. 软件学报, 2006, 17(3): 568-576.

[11] NARAIN R, SEWALL J, CARLSON M, et al. Fast animation of turbulence using energy transport and procedural synthesis[J]. ACM Transactions on Graphics, 2008, 27(5): 1-8.

[12] LENAERTS T, PHILIP D. Mixing fluids and granular materials[J]. Computer Graphics Forum: Journal of the European Association for Computer Graphics, 2009, 28(02): 213-218.

[13] 常元章, 柳有权, 鲍凯, 等. 一种基于粒子的牛顿流体与粘弹性流体统一模拟方法[J]. 计算机学报, 2010, 33(7): 1286-1294.

[14] SHAO X, ZHOU Z, MAGNENAT-THALMANN N, et al. Stable and fast fluid-solid coupling for incompressible SPH[J]. Computer Graphics Forum, 2015, 34(1): 191-204.

[15] TAMPUBOLON A P, GAST T, KLÁR G, et al. Multi- species simulation of porous sand and water mixtures[J]. ACM Transactions on Graphics, 2017, 36(4): 1-11.

[16] ŠT’AVA O, BENEŠ B, BRISBIN M, et al. Interactive terrain modeling using hydraulic erosion[C]// Proceeding of the 2008 ACM SIGGRAPH/Eurographics Symposium on Computer Animation. Aire-la-Ville: Eurographica Association, 2008: 201-210.

[17] KRIŠTOF P, BENEŠ B, KŘIVÁNEK J, et al. Hydraulic erosion using smoothed particle hydrodynamics[J]. Computer Graphics Forum, 2009, 28(2): 219-228.

[18] MEI X, DECAUDIN P, HU B G. Fast hydraulic erosion simulation and visualization on GPU[C]//The 15th Pacific Conference on Computer Graphics and Applications (PG’07). New York: IEEE Press, 2007: 47-56.

[19] CORDONNIER G, GALIN E, GAIN J, et al. Authoring landscapes by combining ecosystem and terrain erosion simulation[J]. ACM Transactions on Graphics, 2017, 36(4): 1-12.

[20] 郑顾平, 邢玥, 张荣华. 一种基于GPU的弹坑实时绘制方法[J]. 图学学报, 2016, 37(4): 451-456.

[21] 刘嘉瑞, 杨猛, 吴佳泽, 等. 基于SPH的雨滴打击不规则边界的模拟方法[J]. 图学学报, 2018, 39(3): 411-418.

[22] BEAR J. Dynamics of fluids in porous media[J]. Engineering Geology, 1973, 7(2): 174-175.

[23] 史鑫. 布料湿润过程的模拟[D]. 上海: 上海交通大学, 2015.

[24] FEI Y, MAIA H T, BATTY C, et al. A multi-scale model for simulating liquid-hair interactions[J]. ACM Transactions on Graphics, 2017, 36(4): 1-17.

[25] WANG R, XIAO S J. Motion simulation of wet cloth on contact surface considering air effects and humidity[C]//Proceedings of the 15th ACM SIGGRAPH Conference on Virtual-Reality Continuum and Its Applications in Industry-VRCAI ’16. New York: ACM Press, 2016: 89-95.

[26] 杨志亮, 王章野, 柯晓棣, 等. 多相流灾害场景的真实感建模与绘制[J]. 计算机辅助设计与图形学学报, 2008, 20(8): 1023-1032.

[27] 杨志亮. 泥石流现象的真实感建模与绘制[D]. 杭州: 浙江大学, 2008.

[28] 张远西. 混合水射流冲击强化不同几何特征表面的研究[D]. 郑州: 郑州大学, 2018.

Realistic modeling and rendering of hydraulic erosion terrain

ZHENG Gu-ping1, ZHANG Meng1, SHAO Xu-qiang1, WU Zhao-hui2

(1. School of Control and Computer Engineering, North China Electric Power University, Baoding Hebei 071003, China; 2. Basic Research Innovation Center, Academy of Transportation Sciences, Beijing 100019, China)

Hydraulic erosion topography is a common natural phenomenon in both real life and virtual scenes. In order to study the surface erosion by hydraulic power, a realistic simulation model is proposed by applying the smoothed particle hydrodynamics (SPH) to simulate the wetting and erosion process of the surface affected by fluid. Based on the reasonable simplified description of the real splash erosion in nature, a three stages topographic hydraulic splash erosion model is proposed. Firstly, the effect of water flow on the surface is simulated by the deformation model of projectile impacting flat target. Secondly, the calculation method of water permeability based on Darcy’s law is used to ensure the authenticity of surface water absorption and water diffusion. Finally, a new calculation model is proposed, which takes into account the influence of surface moisture on surface erosion. The experimental results show that the model can simulate the erosion effects of water flow in different directions on surface and the phenomena of surface water absorption and moisture diffusion during the erosion process.

hydraulic erosion; smoothed particle hydrodynamics; splash erosion; Darcy’s law; wet diffusion

TP 391

10.11996/JG.j.2095-302X.2020020169

A

2095-302X(2020)02-0169-09

2019-07-26；

2019-09-13

国家自然科学基金项目(61502168)；中央高校基本科研业务费专项(2018MS068)；北京市自然科学基金项目(4182018)

郑顾平(1960–)，男，河北石家庄人，教授，博士。主要研究方向为计算机仿真、分布式系统和人工智能。E-mail：zhengguping@126.com

邵绪强(1982–)，男，山东泰安人，副教授，博士，硕士生导师。主要研究方向为计算机图形学、虚拟现实等。E-mail：shaoxuqiang@163.com