CN114359500A

CN114359500A - 一种面向滑坡危险范围预测的三维建模与可视化方法

Info

Publication number: CN114359500A
Application number: CN202210229728.6A
Authority: CN
Inventors: 叶健; 吕奕杰
Original assignee: Southwest Jiaotong University
Current assignee: Southwest Jiaotong University
Priority date: 2022-03-10
Filing date: 2022-03-10
Publication date: 2022-04-15
Anticipated expiration: 2042-03-10
Also published as: CN114359500B

Abstract

本发明公开了一种面向滑坡危险范围预测的三维建模与可视化方法,将滑坡概化为流化模型，采用SPH法构建滑坡危险范围数值模型；提出顾及滑坡运动堆积特征的滑坡可视化优化算法，该算法采用Metaball方法渲染SPH数值模型，LOD方法优化滑坡地形建模，以解决大型滑坡计算资源不足的问题。本发明方法提出的三维SPH与Metaball融合建模方法在取得较好的数值模拟结果的同时，可有效保证滑坡可视化模拟的效果，能够提供具有较高准确性的滑坡危险范围预测结果，提出的四叉树多尺度地形渲染算法能有效简化地形渲染，提升程序运行速度。

Description

一种面向滑坡危险范围预测的三维建模与可视化方法

技术领域

本发明涉及滑坡预测模拟技术领域，具体涉及一种面向滑坡危险范围预测的三维建模与可视化方法。

背景技术

高速远程滑坡冲淤距离及危险范围预测对于滑坡灾害防治至关重要，具有很高的研究价值。高速远程滑坡具有滑动时间短、速度快、滑动距离长等特点，对下游社区和财产造成重大损害。SPH模型可应用于高速远程滑坡冲淤距离及危险范围的预测。然而，现有的SPH模型在实际仿真应用中存在诸多不足。

从高速远程滑坡的运动堆积特征来看，滑体从启动到静止既有崩滑后破碎的过程，也有受到雨水冲刷后崩滑体粘合的作用。尽管SPH在数值模拟方面具有一定的优势，但SPH滑坡数值模型是由数万个SPH离散粒子组成的滑体，因此完全基于SPH方法的高速远程滑坡堆积过程三维可视化，无法真实地模拟滑坡运动堆积特征。

利用滑坡模拟预测滑坡危险范围是滑坡灾害防治、减少人员和财产损失的重要手段。尽管滑坡危险范围数值模拟预测的准确性需要重点考虑，但能够反映滑坡危险性关键参数及其运动堆积过程的可视化模拟也不容忽视。此外，由于当前的滑坡模拟通常将大规模地形数据和对应的影像数据一次性调入内存并全部渲染，这是导致滑坡模拟程序运行缓慢且仅能模拟小尺度滑坡区域的重要原因。

通过以上分析不难看出，当前缺乏一种既能保证滑坡危险范围预测准确性及可靠的运行速度，又能通过三维可视化的方法真实再现滑坡堆积过程。

发明内容

为解决现有技术中存在的问题，本发明提供了一种面向滑坡危险范围预测的三维建模与可视化方法，将滑坡概化为流化模型，采用SPH法构建滑坡危险范围数值模型；提出顾及滑坡运动堆积特征的滑坡可视化优化算法，该算法采用Metaball方法渲染SPH数值模型，LOD方法优化滑坡地形建模，以解决大型滑坡计算资源不足的问题，解决了上述背景技术中提到的问题。

为实现上述目的，本发明提供如下技术方案：一种面向滑坡危险范围预测的三维建模与可视化方法，包括如下步骤：

S1、采用SPH法构建滑坡危险范围数值模型模拟滑坡运动；

S2、采用Metaball与SPH融合可视化建模，为SPH粒子生成多边形包围体，再现滑坡运动堆积特征；

S3、使用四叉树多尺度地形渲染算法进行地形渲染，采用不同尺度的地形模型表示地形，利用塔形四叉树结构管理地形数据，通过动态调度地形数据块以及不同尺度地形数据块的组合渲染。

优选的，所述步骤S1中的SPH法具体是指：通过连续性密度法对粒子密度进行更新，公式如下：

；

粒子的运动方程如下式所示：

；

将人工粘度

引入到计算中，保证流体模拟的数值稳定性，公式如下：

；

其中，

的定义如下：

；

由此，得到粒子的加速度，加速度改变粒子速度，速度改变粒子位置，通过使用蛙跳法进行时间积分，即可完成对滑坡运动模拟。

优选的，所述步骤S2中Metaball与SPH融合可视化建模具体包括：设整个三维空间中存在一个能量场，能量场中任意一个空间点具有的能量值等于所有SPH粒子对该点造成的影响的总和，通过使用与SPH法中的密度求和法相同的式子进行计算，能量值计算公式如下：

；

设定一个能量值

，则可以在SPH粒子外部生成一个由具有

能量值的空间点构建的等值面，当SPH粒子互相靠近时，空间点的能量值随之发生变化，进而影响生成的等值面的形状，形成融合的效果。

优选的，所述的等值面的构建采用MC算法，所述MC算法包括如下步骤：

S21、为SPH粒子群生成一个矩形包围盒，并按设置的单元宽度

将该包围盒分割成一个个规则的正方体单元；

S22、遍历所有正方体单元，对于每一个正方体单元，计算其八个顶点各自的能量值，并将这些能量值分别与预先设定的能量值

相比较，若大于或等于

，则说明该顶点位于等值面的内部，标记为0；若小于

，则说明该顶点位于等值面的外部，标记为1；

将八个顶点的标记情况按顶点编号顺序存入一个八位二进制数中，该二进制数的每一位分别对应一个顶点的标记情况，

建立边搜索表EdgeTable和三角形搜索表TriTable，EdgeTable保存正方体单元的边与等值面相交的情况，TriTable保存在正方体单元内部生成三角形面片的情况的三角形顶点连接顺序；

S23、将步骤S22中得到的八位二进制数转为十进制数并代入EdgeTable中得到一个十二位二进制数，新的十二位二进制数的每一位分别对应正方体单元的一条边，表示该正方体单元的哪些边与等值面相交，然后进行线性插值，计算出这些边各自与等值面的交点坐标；

S24、查找TriTable，依据TriTable提供的连接顺序，将步骤S23中所得的交点连接成一个个三角形面片，用这些三角形面片组成的多边形面来逼近等值面。

优选的，所述线性插值的公式如下：

；

其中，

为正方体单元的边与等值面的交点三维坐标，

和

分别为边的两个端点的三维坐标，

和

分别为边的两个端点的能量值。

优选的，所述步骤S3中的使用四叉树多尺度地形渲染算法进行地形渲染具体包括如下步骤：

S31、创建队列1和队列2, 队列1用于存放正在处理的当前层次的地形节点（每一个地形数据块可看作一个四叉树节点，即地形节点），队列2用于存放待处理的下一个层次的地形节点；

S32、采用广度优先遍历四叉树的原则从根节点开始对四叉树进行遍历；在遍历四叉树时，一个四叉树节点即为一个地形节点，该节点所对应的高程点集直接从滑坡区域DEM数据中选取；

S33、将遍历到的节点送入队列1，判断该节点在视景体内是否可见，若是不可见的，则直接出队列，不予处理；若是可见的，则使用节点细分评价体系判断该节点是否需要细分，若不需要细分，则直接渲染该节点；若需要细分，则将该节点出列，进行细分，并将该节点细分产生的下一个层次的四个子节点送入队列2，等待下一轮的处理；

S34、当程序运行到队列1为空时，则交换队列1和队列2中的节点，并重复步骤S31至步骤S33，直到所有需要处理的节点均已处理完毕为止；对于那些在视景体内可见的叶子节点，则直接进行渲染。

优选的，所述步骤S32中选取的具体步骤如下：

将DEM数据存入一个二维数组

中，每个地形数据块可由9个高程特征点构成，用一个

大小的二维数组

进行存储，通过行索引值

和列索引值

进行访问，则任意地形数据块点集

中某一高程点的高程值

与原始地形数据块（未经分层分块的地形数据块）的高程值

相对应，因此可通过映射关系计算出某一地形数据块行索引值

和列索引值

在原始地形数据块中对应的行索引值row和列索引值col，计算公式如下：

；

优选的，所述步骤S33的节点细分评价体系具体公式如下：

；

式中，

是节点细分评价值，

为视点到地形节点中心的距离，

为调节因子，

、

、

和

分别对应地形节点四条边的中点处高程值与边的两个端点处高程值的平均值之差。

本发明的有益效果是：本发明方法提出的三维SPH与Metaball融合建模方法，有机结合metaball技术、SPH建模以及LOD地形优化算法，充分发挥三项技术的各自优势，使得滑坡模拟不会因为SPH与Metaball融合建模技术由于模拟计算量过大，加之地形渲染耗费大量的计算资源，导致滑坡模拟无法实现，更为重要的是本发明在取得较好的数值模拟结果的同时，有效保证了滑坡可视化模拟的效果。

附图说明

图1为本发明方法步骤流程图；

图2为实施例中SPH粒子互相靠近时外部等值面的变化情况图；

图3为实施例中SPH滑坡粒子群外部生成矩形包围盒并分割包围盒示意图；

图4为实施例中滑坡运动堆积特征的可视化模拟图，图4（a）为滑坡碎屑化过程，图4（b）为滑坡碎屑物质重新粘合过程；

图5为实施例中地形节点点集示意图；

图6为实施例中滑坡地形简化渲染效果图；

图7（a）为0s时运动堆积状态厚度截图，图7（b）为10s时运动堆积状态厚度截图，图7（c）为20s时运动堆积状态厚度截图，图7（d）为40s时运动堆积状态厚度截图，图7（e）为60s时运动堆积状态厚度截图，图7（f）为100s时运动堆积状态厚度截图；

图8 为实施例中模拟滑体运动足迹与实际滑坡堆积范围对比图；

图9为实施例中茂县滑坡灾害场景模拟与实际情况对比图，图9（a）是当t=0s时对茂县滑坡灾害场景的可视化模拟效果，图9（b）是当t=20s时对茂县滑坡灾害场景的可视化模拟效果，图9（c）是当t=30s时对茂县滑坡灾害场景的可视化模拟效果，图9（d）是当t=40s时对茂县滑坡灾害场景的可视化模拟效果，图9（e）是当t=100s时对茂县滑坡灾害场景的可视化模拟效果，图9（f）滑坡实际情况；

图10为实施例中滑坡模拟程序渲染帧率随运行时间变化图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

与三维数值模拟相比，二维数值模拟具有便捷和易于实现的优势，因此目前滑坡数值模拟仍以二维数值模拟为主，可以更有效地研究滑坡灾害机理，但在滑坡危险范围的预测和可视化表达方面仍然存在局限性，特别是二维数值模拟展现的二维可视化信息难以展现滑坡至灾的三维细节信息。三维SPH数值模拟结合三维可视化模拟是滑坡模拟和预测向着多维度、多源信息表达方向发展可行有效的方案，具有很好的发展前景。

在滑坡模拟三维可视化方面，Metaball技术是一种适用于建立可变形表面的技术，与SPH融合建模时常被用于构建三维流体表面，本发明人尝试将metaball技术与SPH融合建模，从三维可视化模拟效果上看，metaball技术与SPH融合建模可为SPH粒子生成多边形包围体，当SPH粒子相互靠近时，相邻SPH粒子的包围体会互相融合生成新的包围体，由于包围体在可视化时能够模拟滑坡碎屑物质的质感且便于进行纹理映射，因此可弥补SPH粒子在滑坡可视化表达方面的不足，是对SPH滑坡数值模拟的重要补充，但是该项融合建模技术由于模拟计算量过大，且地形渲染将耗费大量的计算资源，因此简单的融合metaball技术与SPH建模方法无法在大规模滑坡模拟中得以实际应用。

另一方面，大规模滑坡滑程远、堆积范围广，而地形是滑坡模拟的边界，因此滑坡模拟存在的另一个重要问题是滑坡地形数据的调度和可视化问题。当前,复杂地形边界条件下的滑坡模拟常常将大规模数字高程模型数据一次性加载到计算机内存中并全部渲染，再加上滑坡模型复杂度高、科学计算复杂，因此常造成滑坡模拟程序运行速度缓慢甚至无法运行的后果。基于以上问题，发明人通过长时间的实验与论证，克服了上面的问题，具体方案如下：

本发明提供一种技术方案：一种面向滑坡危险范围预测的三维建模与可视化方法，如图1所示，包括如下步骤：

采用SPH法构建滑坡危险范围数值模型模拟滑坡运动；

采用Metaball与SPH融合可视化建模，为SPH粒子生成多边形包围体，再现滑坡运动堆积特征；

使用四叉树多尺度地形渲染算法进行地形渲染，采用不同尺度的地形模型表示地形，利用塔形四叉树结构管理地形数据，通过动态调度地形数据块以及不同尺度地形数据块的组合渲染。

SPH法滑坡建模

SPH方法是一种无网格拉格朗日粒子法，其基本思想是将连续的流体离散为一组粒子，每个粒子都有其各自的密度、速度等属性，空间中任意粒子

的属性值

可由其周围粒子的属性值内插近似求得：

（1）

式中，

，

和

分别代表粒子

和粒子

的位置坐标，

为粒子

的属性值，

为粒子

的属性值，

为粒子

的质量，

为粒子

的密度，

为搜索半径为

的光滑核函数。

在SPH方法中有两种密度近似法：一种是对密度直接应用SPH近似法，称为密度求和法；另一种方法是连续性密度法。由于连续性密度法对于具有强间断问题的模拟（如爆炸、高速冲击等）可以得到较好的结果。

因此，本发明选择连续性密度法（Monaghan 1994）对粒子密度进行更新：

（2）

式中，

和

分别代表粒子

和粒子

的速度。

粒子的运动方程（Becker and Teschner 2007）如下式所示：

（3）

式中，

代表粒子

的速度，

和

分别代表粒子

和粒子

的压强，

代表外力作用，通常为重力。

为保证流体模拟的数值稳定性，通常将Monaghan型人工粘度

（Monaghan 2005）引入到数值计算中：

（4）

其中，

的定义如下：

（5）

式中，

，

为一可调粘度参数，

为流体中的声速，

的引入是为了防止分母等于0。

由式（3）和式（4）可得到粒子的加速度，加速度改变粒子速度，速度改变粒子位置，如此，通过使用简单易用的蛙跳法（Leap-Frog，LF）（Desbrun and Cani1996）进行时间积分，即可模拟滑坡运动。

可视化渲染

从高速远程滑坡的运动堆积特征来看，滑体从启动到静止既有崩滑后破碎的过程，也有受到雨水冲刷后崩滑体粘合的作用。尽管SPH在数值模拟方面具有一定的优势，但SPH滑坡数值模型是由数万个SPH离散粒子组成的滑体，因此完全基于SPH方法的高速远程滑坡堆积过程三维可视化，无法真实地模拟滑坡运动堆积特征。尽管Metaball方法无法提升滑坡危险范围数值模拟的正确性，但通过使用Metaball技术能够从可视化角度真实再现滑坡运动堆积特征。

Metaball与SPH融合可视化建模方法如下：

设整个三维空间中存在一个能量场，能量场中任意一个空间点具有的能量值等于所有SPH粒子对该点造成的影响的总和，通过使用与SPH法中的密度求和法相同的式子进行计算，能量值计算公式如下：

（6）

式中，

为该点具有的能量值，

为该点的位置坐标，

和

分别为某一 SPH粒子的质量和位置坐标，

为搜索半径为

的光滑核函数。

如图2所示，通过预先设定一个能量值

，则可以在SPH粒子外部生成一个由具有

能量值的空间点构成的等值面。当SPH粒子互相靠近时，空间点的能量值随之发生变化，进而影响生成的等值面的形状，形成融合的效果。

等值面的构建采用Marching Cubes（MC）算法，该算法最早由Lorensen和Cline于1987年提出（Lorensen and Cline 1987），是从三维离散数据场中提取等值面的经典算法，具有生成等值面质量好、实现简单等优点。

该算法分为四步：第一步，为SPH滑坡粒子群生成一个矩形包围盒，并按预先设置的单元宽度

将该包围盒分割成一个个规则的正方体单元，如图3所示。

第二步，遍历所有正方体单元，对于每一个正方体单元，计算其八个顶点各自的能量值，并将这些能量值分别与预先设定的能量值

相比较，若大于或等于

，则说明该顶点位于等值面的内部，标记为0；若小于

，则说明该顶点位于等值面的外部，标记为1。接着将八个顶点的标记情况按一定的顶点编号顺序存入一个八位二进制数中，该二进制数的每一位分别对应一个顶点的标记情况。

由于每个顶点存在两种标记情况，因此八个顶点共存在有256种情况，基于这256种情况建立一个边搜索表EdgeTable和一个三角形搜索表TriTable，EdgeTable保存正方体单元的边与等值面相交的256种情况，TriTable保存在正方体单元内部生成三角形面片的256种情况的三角形顶点连接顺序。

第三步，将第二步中得到的八位二进制数转为十进制数代入EdgeTable中得到一个十二位二进制数，新的二进制数的每一位分别对应正方体单元的一条边，表示该正方体单元的哪些边与等值面相交，然后通过式（7）进行线性插值，计算出这些边各自与等值面的交点坐标：

（7）

式中，

为正方体单元的边与等值面的交点三维坐标，

和

分别为边的两个端点的三维坐标，

和

分别为边的两个端点的能量值。

最后一步，通过查找TriTable，依据TriTable提供的连接顺序，将第三步中所得的交点连接成一个个三角形面片，用由这些三角形面片组成的多边形面来逼近等值面。理论上来说，单元宽度

设置得越小，则分割出来的正方体单元越多，最后生成的多边形面越接近等值面，但也意味着计算量的急剧上升，因此选择一个合适的

也很重要。

由于等值面完全包裹SPH粒子，因此在三维空间中可称之为包围体。对于单个SPH 粒子，Marching Cubes算法中的单元宽度

设置得越小，则生成的多边形包围体越近似于球体。在模拟水流时，通常会取较小的

以模拟光滑的水珠。但对于滑坡来说，其碎屑物质通常是不规则的且棱角分明的。因此在模拟滑坡时，

可以取较大的值，如此既可生成粗糙的包围体，模拟碎屑物质的质感，又可减少Marching Cubes算法的计算量，提升程序运行速度。

由于Metaball的特性，当SPH粒子相互远离时其对应的多边形包围体会互相分离，基于这个特性，可模拟高速远程滑坡崩滑后由于滑体内部粘聚力减弱以及重力作用下不断破碎、碎屑化的可视化过程，如图4中图4（a）所示。与之相反的，当SPH粒子相互靠近时其对应的多边形包围体会互相融合，依据这个特性，可模拟滑体受到雨水冲刷作用后滑体碎屑化以后重新粘合的可视化过程，如图4（b）所示。

相较于粒子，用多边形包围体对滑体进行渲染更符合实际情况，既可利用相邻包围体互相融合的特性表达滑坡整体运动趋势，又可利用相互远离的包围体互相分离的特性模拟部分滑坡滑动后碎屑化的过程，从可视化角度反映滑坡运动堆积特征。

滑坡地形建模

地形是滑坡模拟的边界，是滑坡场景的重要组成部分。本发明在渲染场景时，利用滑坡区域数字高程模型DEM数据构建地形，将滑坡区域数字正射影像图（DigitalOrthophoto Map，DOM）作为纹理映射至地形上，生成滑坡地形。

为了保证滑坡模拟的顺利进行，需要将DEM数据一次性全部读入内存中，如果在模拟程序运行时将DEM数据全部渲染将会大大降低程序运行速度。使用LOD技术是解决上述问题的有效途径。常用的LOD算法有：基于四叉树的LOD算法、实时优化自适应网格算法、递进网格算法以及状态无关单通过自适应细化算法。其中，由于四叉树结构方便对地形数据进行分层分块和管理，因此基于四叉树的LOD算法便成了首选。

为此，本发明以基于四叉树的LOD算法思想为框架，提出了一种四叉树多尺度地形渲染算法，采用不同尺度的地形模型表示地形，依据人类视觉原理，在靠近观察者的区域用大量小尺度地形模型渲染高分辨率地形，在远离观察者的区域用少量大尺度地形模型渲染低分辨率地形，运用塔形四叉树结构管理地形数据，通过动态调度地形数据块以及不同尺度地形数据块的组合渲染，达到简化地形渲染的目的，在保证可视化效果的前提下，提高地形渲染效率，提升程序运行速度。

滑坡地形数据的组织与调度

本发明采用塔形四叉树结构对DEM和DOM数据进行组织与管理。将整个滑坡地形区域按四叉树结构划分成一个个小的地形数据块，每一个地形数据块可看成一个四叉树节点，称为地形节点。从四叉树根节点开始，除了叶子节点外，每个节点均可细分产生下一个层次的四个子节点。

每个地形节点的索引值都依据其所属层次的地形数据块数

，从0开始至n−1依次编写。由当前地形节点细分产生的下一层次子节点的索引值

的计算公式（Lyu et al. 2020）为：

（8）

式中，

为当前地形节点的索引值，

为当前地形节点的层次值。

每个地形节点所对应的高程点集直接从滑坡区域DEM数据中选取，如图5所示，DOM数据同理。

选取方法如下：将DEM数据存入一个二维数组

中，数组中的每个元素用于存储其对应点位的高程值，数组的行数和列数与DEM数据的行数和列数相同。本发明采用三角形扇的形式绘制每个地形节点，每个地形节点均由9个高程点构成，可用一个

大小的二维数组

进行存储，通过行索引值

和列索引值

进行访问。则任意地形数据块点集

中某一高程点的高程值

与原始地形数据块（未经分层分块的地形数据块）的高程值

和列索引值

在原始地形数据块中对应的行索引值row和列索引值col，计算公式如下：

（9）

式中，变量

和

的指代与式（8）相同。

地形节点的调度主要由两个方面的内容来决定：①地形节点的细分程度；②地形节点的取舍。

针对第一点，本发明构建了一个地形节点细分评价体系（Liu et al. 2010），用于对地形节点细分与否做出判断，具体公式如下：

（10）

式中，

是节点细分评价值，

为视点到地形节点中心的距离，

为调节因子，

、

、

和

分别对应地形节点四条边的中点处高程值与边的两个端点处高程值的平均值之差，

、

同理，对应的边为对角线。

针对第二点，本发明在每次绘制地形节点前，都对节点进行了视景体裁剪，仅保留处于视景体内部的节点，而舍弃位于视景体外部的节点，并由此判断需要调度的节点对象。另外，由于地形节点本身的不规则性，本发明使用节点包围盒（Decoret et al.1999）来替代原有节点进行视景体裁剪判断。

四叉树多尺度地形渲染算法

为了满足滑坡数值模拟过程中边界处理的需要，在程序运行过程中，需要将滑坡区域DEM数据完整读入内存中，但在渲染时没有必要渲染全部数据，仅在保证可视化效果的前提下渲染部分数据即可。

为此，本发明提出了一种四叉树多尺度地形渲染算法控制地形渲染：

1）创建两条队列（队列1用于存放正在处理的当前层次的地形节点，队列2用于存放待处理的下一个层次的地形节点）。

2）采用广度优先遍历四叉树的原则从根节点开始对四叉树进行遍历。在遍历四叉树时，一个四叉树节点即为一个地形节点，该节点所对应的高程点集直接从滑坡区域DEM数据中选取，DOM数据同理，选取方法如上所述。

3）将遍历到的节点送入队列1，判断该节点在视景体内是否可见。若是不可见的，则直接出队列，不予处理；若是可见的，则使用节点细分评价体系（见上面的内容）判断该节点是否需要细分。若不需要细分，则直接渲染该节点；若需要细分，则将该节点出列，进行细分，并将该节点细分产生的下一个层次的四个子节点送入队列2，等待下一轮的处理。

4）当程序运行到队列1为空时，说明某一个层次的所有节点均已处理完毕，则交换队列1和队列2中的节点，并重复如上处理步骤，直到所有需要处理的节点均已处理完毕为止。对于那些在视景体内可见的叶子节点，则直接进行渲染。

通过利用以上算法控制地形渲染，可有效减少不必要的渲染资源开销，提升程序运行速度，保证滑坡场景的渲染效率。对于不同层次地形节点拼接处产生的裂缝（Zhang etal. 2005），采用跳点法（Decoret et al.1999）进行消除。如图6所示，图6用格网模式展示了滑坡地形简化渲染效果。

实验验证分析

本发明以2017年四川茂县滑坡为例，搭建了滑坡灾害场景模拟系统，对滑坡危险范围进行预测。实验平台硬件参数如下：处理器：英特尔 Core i5-9400F，显卡：NvidiaGeForce GTX 1050 Ti。编程环境为C++和OpenGL。

1）滑坡危险范围数值模拟与分析

茂县滑坡从小规模垮塌到主滑体启动、失稳运动、到最终停积整个过程用时约 100s，其中主滑坡仅持续60s，滑坡启动体积约为

，停止时平均堆积厚度大于 10m，最大堆积厚度大于32m，图7中黑色轮廓线圈出了滑坡堆积范围（Ouyang et al. 2017b; Xu et al. 2017）。

本发明对茂县滑坡的模拟结果如图7（a）-（f）所示：图中展示了模拟滑体在0s、10s、20s、40s、60s、100s这些不同时刻的运动堆积状态(时间序列按照从左到右从上到下的次序)，滑体的不同颜色代表不同的堆积厚度。0s至10s期间失稳的山体在重力作用下开始逐渐垮塌，10s至40s期间滑体沿着斜坡坡面高速运动并不断解体破碎，40s左右滑体前部到达坡脚平坦的地面，由于受到阻挡，滑体运动速度开始减缓，40s至60s期间滑体在平坦的地面上不断堆积，60s后整个滑坡趋于停止，60s至100s期间是滑坡的最后调整阶段，仍有一些碎石掉落，100s后滑坡几乎停止。模拟滑体在t=100s时的平均堆积厚度为7.3m。模拟滑体的运动堆积过程与滑坡实际情况吻合较好。

为了准确度量模拟滑坡堆积范围与实际滑坡堆积范围的吻合程度，本发明在滑坡模拟过程中每间隔1s便对模拟滑体的运动状态进行截图（参照图7（a）-（f），在模拟结束后将这些截图叠加合成一张模拟滑体运动足迹图，该图中记录有每个SPH粒子的运动足迹。接着，这里引入一个模拟滑坡拟合度的概念，用以描述模拟滑坡与实际滑坡的相似程度，如图8所示，其具体数值通过对比模拟滑体运动足迹图与实际滑坡堆积范围图的相似程度来确定，计算公式如下：

（11）

式中，

为实际滑坡在图像中所占的像素数量（对应图8中的黑色区域加灰色区域），

为模拟滑坡在图像中重合于实际滑坡的像素数量（对应图8中的灰色区域），

为模拟滑坡在图像中不重合于实际滑坡的像素数量（对应图8中的白色区域），

为模拟滑坡拟合度值，该值越接近1说明模拟滑坡堆积范围与实际滑坡堆积范围吻合程度越高，模拟滑坡与实际滑坡相似程度越高，滑坡模拟效果越好，反之则越差。

通过使用式（11）进行计算可得本发明模拟滑坡拟合度值为0.8542，模拟滑坡堆积范围与实际滑坡堆积范围的吻合程度较高，滑坡模拟效果较好。

2）可视化渲染与分析

图9为实施例中茂县滑坡灾害场景模拟与实际情况对比图，图9（a）-（f）中黑色轮廓线圈出了滑坡堆积范围。图9（a）至图9（e）展示了滑坡启动后的0s至100s期间对滑坡灾害场景的可视化模拟效果。图9（a）中用裂缝线圈出了模拟滑体初始化轮廓。图9（a）至图9（e）中使用了SPH与Metaball融合建模方法模拟滑体，相较于图7（a）-（f）中使用纯SPH方法模拟滑体，在进行纹理映射后，前者的可视化效果明显优于后者。图9（e）展示了当t=100s时对茂县滑坡灾害场景的可视化模拟效果，与实际情况（图9（f））相比，二者相似程度较高，滑坡堆积范围吻合较好。滑坡模拟所用SPH粒子数量为45630个，模拟场景中滑坡地形三角面的渲染数量是随着视点移动而动态变化的，但总体控制在10万个以下。滑坡模拟程序渲染帧率变化情况如图10所示。

通过对比未做简化地形渲染操作和使用本发明方法简化地形渲染这两种情况下程序渲染帧率随运行时间变化的情况，可以看出：在使用本发明方法简化地形渲染后，整个程序的运行速度得到明显提升。另外，当程序运行时间超过100s后，SPH数值计算程序停止运行，仅剩可视化渲染程序在持续运行，此时，未做简化地形渲染操作的情况下程序渲染帧率仅维持在18帧/秒左右，而使用本发明方法简化地形渲染的情况下程序渲染帧率能维持在60帧/秒以上，且两者间的差距会随地形数据量的增加而不断扩大，该现象对于程序运行时间未超过100s时亦是如此。

本发明主要由三类建模方法构成，分别为SPH法滑坡建模、Metaball可视化渲染以及LOD滑坡地形建模方法构成。滑坡模拟主要通过数值模拟预测滑坡危险范围，可视化模拟反映灾情信息，数值模拟的准确性是滑坡危险范围预测最为核心的问题，而模拟的准确性和可视化灾情信息都需要模拟发挥模拟计算性能。由于地形（由DEM地形数据和DOM遥感数据构成）是滑坡模拟的边界，因此SPH与地形的交互耦合必不可少，由于大规模滑坡如茂县滑坡模拟时SPH数值模拟的粒子数必须达到一定数量（本次模拟达到了45630个），才能保证模型达到85.42%的危险范围预测的正确率，如果没有LOD地形优化方法，地形渲染将会占用大量的计算资源，导致分配给SPH数值模拟的计算资源大幅减少，势必造成模拟所用的SPH粒子数的大幅减少，滑坡危险范围预测将因为无法发挥模拟计算性能，导致正确率大大降低。同理，SPH/Metaball融合可视化建模时，LOD地形优化方法也是发挥模拟计算性能，有效发挥可视化表达灾情信息最可靠的保证。没有SPH/Metaball与LOD三类建模方法的深度融合，加载的SPH粒子数大幅减少，势必造成可视化模拟滑坡堆积时滑坡岩土的分离和融合的模拟非常粗糙，完全不能模拟出滑坡的灾情信息和运动堆积效果。数值模拟和可视化模拟需要将LOD优化方法、SPH以及Metaball的深度融合，是准确模拟滑坡危险范围、高保真表达滑坡灾情信息最为有效的手段，可创造性的解决了因模拟性能不佳所造成的模拟准确性问题，以及滑坡灾情信息可视化表达力不强的问题。

综上，本发明提出的三维SPH与Metaball融合建模方法在取得较好的数值模拟结果的同时，可有效保证滑坡可视化模拟的效果，能够提供具有较高准确性的滑坡危险范围预测结果，提出的四叉树多尺度地形渲染算法能有效简化地形渲染，提升程序运行速度，本发明方法是行之有效的。通过2017年四川茂县高速远程滑坡危险范围模拟验证了模型在数值模拟和可视化模拟方面的卓越特性，为滑坡模拟及危险范围预测向着多维度、多源信息表达方向发展提供了行之有效的方案。

尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换,凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。