杜志强，李 静

(1.武汉大学 测绘遥感信息工程国家重点实验室，湖北 武汉 430079；2. 地球空间信息技术协同创新中心，湖北 武汉 430079；3.香港中文大学 太空与地球信息科学研究所，香港 999077)

0 引 言

极端天气频发，自然灾害多有发生，其危害面广，破坏性大，是对人民生命财产安全的最大威胁和对社会经济发展的重大制约因素[1-2]。滑坡属于突发性自然灾害，其危害不言而喻[3]。由于突发性自然灾害不可预测，爆发突然，其拍摄往往比较困难[4]，所以，突发性自然灾害模拟仿真尤为重要，其有助于增强人们对自然灾害的认知，并帮助决策者采取正确的防灾救灾措施。

当前滑坡模拟的方法种类很多，其中将滑坡运动假想为单点运动的集中质量模型[5]和基于连续介质的力学模型[6]是由Hunger[7]于1995年归纳的滑坡两大类模拟仿真模型；Heim提出的“雪橇模型”是最著名的集中质量模型，它是描述滑坡运动的第一个模型，其将滑坡假想为一个单一的点，模型的优点在于它的简单明晰，虽然可计算出滑坡运动的重力，但无法预测滑坡的滑动方向；基于连续介质力学模型在滑坡滑动距离及相关参数预测方面有重要作用，Sassa[6]利用三维模型来模拟滑坡体的分布特征，在研究过程中将滑坡体视为像泥石流一样的流体。理论上来讲，虽然可以应用模型定量地计算出滑坡的滑动体积、影响范围，但由于模型自身的精确度和基本假定等因素影响，滑坡的实际运动不可能完全符合模型的计算结果，另一方面，滑坡本身的诱发因素很多，在滑动过程中也有较强的随机性和模糊性。而粒子系统是一种可以生成模糊及不确定对象的计算模型的方法[8]，到目前为止，其被认为是用来模拟不规则物体运动的最成功的图形生成算法。

本文进行滑坡模拟的主要思路是利用PhysX物理引擎中的粒子系统中的粒子来代表滑坡体中的泥石等成分，通过粒子在模拟灾害环境中的运动，来模拟滑坡体在真实环境中的分布及运动情况。

1 粒子系统基本原理

粒子系统采用大量具有生命周期和特征属性的微小粒子作为基本单元来描述模糊对象的运动特征[9]。在粒子系统中，每个粒子单元都具有形状、大小、颜色、位置、速度大小和方向、生命周期等属性，其中生命周期用于描述粒子在系统中的生存时间[10]，图1为粒子在系统中随时间变化的流程图。

图1 粒子生存周期Fig.1 Life span of particle

其将滑坡假想为一个单一的点，模型的优点在于它的简单明晰，虽然可计算出滑坡运动的重力，但无法预测滑坡的滑动方向；基于连续介质力学模型在滑坡滑动距离及相关参数预测方面有重要作用，Sassa[6]利用三维模型来模拟滑坡体的分布特征，在研究过程粒子在系统中的运动以系统中的参数为基准，并添加一些随机变量进行控制，使粒子的运动变化带有很强的随机性，因此，可以很逼真地模拟出模糊对象的运动特征。

粒子系统是由很多细小的粒子的集合共同表示一个模糊对象所构成的[11]。用粒子系统来模拟滑坡动态过程的具体方法如下：

1）产生新的粒子加入系统中；

2）赋予每一新粒子以一定属性；

3）删除那些已经超过其生命周期的粒子；

4）根据粒子的动态属性对粒子进行移动和变换；

5）绘制并显示由有生命的粒子组成的图像。

1.1 粒子属性

粒子的运动状态源于其属性的设置，以下为本文中所用到的粒子属性[11]。

1）粒子初始位置

粒子的位置信息用三维坐标（X,Y,Z）进行描述。本文中粒子的初始位置设定为灾害发生的源点。

2）粒子初始速度

粒子系统生成时的形状决定了粒子的运动方向。如果为球形，粒子将从粒子系统的原点向外运动，如果是圆形或者矩形，粒子将向XY平面的上方运动，同时，可以调节粒子的发射角度来控制粒子的运动方向。本文中粒子的运动方向由滑坡该点处的坡向决定。

3）粒子生命周期

生命周期是用时间的方式描述了粒子的生命历程，本文中通过控制时间间隔来控制粒子生命周期。

1.2 粒子动力学特征

粒子产生后，在三维空间中开始运动，此时，粒子的属性在运动过程中会不断改变。粒子的动力学特征是指每个粒子的属性值在每一帧都要更新变化。

粒子在n+1帧的属性由n帧的属性决定，在第n+1帧时粒子的运动位置[12]为：

式中，Xn+1表示在第f+1帧时粒子的X坐标，VXn表示在第f帧时粒子的速度，△t表示第f帧到第f+1帧的时间差，其他同理。

第n+1帧时粒子的速度[12]为：

式中，Vn+1表示在第f+1帧时粒子的速度，a为粒子的加速度。

2 滑坡模拟模型设计

本文进行滑坡模拟的主要思路是利用PhysX物理引擎中的粒子系统中的粒子来代表滑坡体中的泥石等成分，通过粒子在模拟灾害环境中的运动，来模拟滑坡体在真实环境中的分布及运动情况。

本文所用到的数据包括30 m全国地形数据、高分辨率影像数据（主要为高分数据）以及间接观测数据（如坡度、坡向、水平地震力等）。首先对高分辨率影像数据进行影像纠正、信息提取等处理之后得到滑坡分析的宏观信息。利用间接观测数据对滑坡相关信息进行参数率定，得到有效信息。根据极限平衡分析方法原理，生成滑坡的有效下滑数据。最后，利用有效下滑数据对PhysX物理引擎中的粒子系统进行初始化，粒子发射器在滑坡范围内发射滑坡粒子，滑坡粒子将会受到环境影响而运动。每一时刻，循环所有在生存周期内的粒子，判断滑坡粒子是否触地，若不是，说明滑坡粒子还未落地，则该粒子暂时不处理，若是，则根据粒子位置判断该粒子所在的格网，将对应格网处数组加一，由此可得滑坡粒子数量数组，即落入每个格网中的粒子总数。每个格网的模拟时间即为该格网中的粒子数与时间间隔的乘积，若模拟时间大于0，则说明该格网处于受灾范围区域，由此可得到滑坡灾害范围。

本文进行的滑坡模拟的流程图如图2所示。

图2 滑坡模拟流程图Fig.2 Process of landslide simulation

2.1 多参数率定方法原理

传统的滑坡稳定性分析方法，几乎没有考虑数据获取过程中的人为误差、系统误差对数据的影响。如根据滑坡变形速率的分析方法，没有经过去差方法或只经过简单的线性拟合方法直接利用变形数据特点很容易引起滑坡预警的误报错报[12]。如图3红色变形监测区域，滑坡变形速率突然增大而后又出现减速现象。

图3 传统滑坡监测方法容易引起误报的监测点Fig.3 Monitoring of false positives caused by traditional landslide monitoring method

滑坡稳定性动态分析需要基于多种来源的观测数据，包括地面监测数据、高分辨率SAR数据、光学遥感数据、高时间、高空间分辨率数据、历史灾情数据，在充分分析各种数据的特点和使用目的的基础上利用尺度转换技术、空间数据插值技术和属性数据空间化技术实现多源数据的融合。在利用不同时刻的不能直接测定的土壤参数数据时，可利用线性回归或多参数反演的方法来获取相对时间的参数值，也可采用参数智能率定的方法，确保参数的实时准确性。

2.2 极限平衡分析方法原理

极限平衡分析方法是边坡稳定分析中最常用的方法，它是通过分析在临近破坏状况下，土体外力与内部强度所提供抗力之间的平衡，计算土体在自身和外荷作用下的土坡稳定性程度，通常以边坡稳定系数表示[14]。

基于极限平衡分析方法[15]的计算公式（5），原理如图4所示。

图4 平衡分析方法原理示意图Fig.4 Principle diagram of the limit equilibrium method

参数说明见表1。

表1 参数说明Tab.1 Parameter specification

2.3 滑坡粒子系统模型设计

根据粒子系统的基本原理和滑坡灾害体属性分析，本小节将对滑坡灾害体的粒子系统进行详细设计。滑坡粒子由粒子发射器发射，故粒子发射器的设计是极其重要的，其用于产生新粒子，并初始化粒子属性。由于本文进行的是滑坡的有源模拟，故可选用点发射源。滑坡粒子的特性包括滑坡模拟参数（包括模拟间隔、模拟时间等）、滑坡DEM数据、坡向数据、滑坡静止密度、粘滞系数，粒子点特性（包括大小、颜色、速度等）以及滑坡边坡稳定系数等粒子属性中的参数可直接或间接地从极限平衡分析算法中得到。

3 滑坡模拟仿真实验及结果分析

2015年11月13日22时50分许，浙江省丽水市莲都区雅溪镇里东村发生山体滑坡。山体滑坡塌方量30余万立方米，导致27户房屋被埋，21户房屋进水。截至2015年11月19日12时39分，现场找到第39位失联人员。此次山体滑坡已有38人遇难，1人经抢救后生命体征平稳。

图5为丽水滑坡实地影像图：

图5 丽水滑坡实地影像图Fig.5 The image of landslide in Lishui

利用本软件对从雅溪镇里东村处发生的滑坡灾害进行模拟，滑坡源头设置在北纬28°38′49.2″，东经119°54′36″，模拟开始时间为2015年11月13日22点50分，持续时间0.5 h，模拟时间间隔设置为5 s。图6为模拟结果。

图6 滑坡模拟结果图Fig.6 The result of the landslide simulation

图7 实际灾害与模拟灾害叠加图Fig.7 The overlay chart of the actual disaster with the simulation result

由于模拟采用的DEM数据精度以及时效性等原因导致灾害模拟的范围与实际灾害范围基本吻合，但在细节上存在差异。另外，由于模拟采用的是有源模拟，故源点位置的选择对模拟结果的影响至关重要。本例中选择的泥石流的位置位于里东村。从图7中看出模拟结果与实际灾害范围基本吻合。

4 结束语

在地质及物理学界，已有不少对山体滑坡等自然灾害现象的研究工作，但它们关注的重点往往是这些灾害现象的成因、运动机制及数值模拟等。这些计算模型由于太过复杂并不适合于真实感模拟。本文在传统模拟方法的基础上引入了粒子系统，对滑坡过程进行了动态模拟。模拟结果表明，本文方法能较好地模拟滑坡发生过程，展现滑坡危害区域，提供可靠防灾减灾资料。本文的模拟结果与实际灾情范围具有较强的一致性，但由于DEM精度以及源点位置等因素，导致模拟结果出现部分偏差。