CN109992830A - 基于物质点方法的山体滑坡灾害场景模拟方法 - Google Patents

Abstract

本发明公开了一种基于物质点方法的山体滑坡灾害场景模拟方法。其步骤为：1)在计算机图形学领域对山体滑坡这一灾害现象进行了建模与绘制。2)在传统MPM方法的基础上改进了屈服面准则，引入修改后的剑桥模型。3)提出了基于MPM的颗粒状材质与刚体作用的双向耦合模型，使得山体滑坡过程中的土石相互作用更加真实，提高了大规模复杂灾害场景模拟的效率。本发明解决了在计算机图形学中山体滑坡灾害场景模拟缺失的问题，能够模拟山体滑坡发生过程中的土颗粒之间的运动、石块刚体之间的碰撞以及土颗粒与石块刚体的耦合现象。

Description

基于物质点方法的山体滑坡灾害场景模拟方法

技术领域

本发明涉及基于粒子的颗粒流模拟和流固双向耦合模拟领域，尤其涉及一种基于土力学与以及图形学的基于物理的山体滑坡模拟框架。

背景技术

基于物理的自然场景仿真一直是计算机图形学领域研究的一大热点。其内容主要可分为用于游戏娱乐的实时仿真以及用于影视特效的离线仿真两大分支。对于仿真效率以及仿真效果的评判角度而言，前者更注重计算的效率,即实时性；而后者更为注重视觉效果，即真实感。但是即便如此，在影视界的真实感特效制作中，若能够提高计算绘制效率、减少艺术创作者的等待时间也是十分关键的，这有助于缩短制作周期、大大降低人工制作成本，从而提高市场竞争力。

山体滑坡是一种常见的地质灾害，是指斜坡上的岩土在重力或者静水压力作用下，发生整体移动而产生的灾难性地质现象。山体滑坡包含广泛的地面运动，包括岩崩、深层滑坡、泥石流、土流等。山体滑坡能够发生在各种各样的环境中，以陡峭的坡度为特征，从山脉到海岸悬崖，甚至能够发生在水下，形成海底滑坡。重力是滑坡发生的主要驱动力，但影响边坡稳定性的其他因素也会导致土样破坏从而发生滑坡。在许多现象情况下，滑坡是由某一特定事件引起的，如暴雨、地震、为修路而开挖的斜坡等，但是确定滑坡的原因往往是十分困难的。

MPM(Material Point Method,即：物质点方法)是一种基于拉格朗日粒子法以及欧拉网格法的混合方法，作为一种数值计算方法，可以用于计算固体、可变性物体以及其他连续介质等。其也被认为是FLIP方法的广义化的方法，而FLIP(Fluid Implicit Particle，即：流体隐式粒子法)作为一种流行的网格方法，已经被广泛地应用于流体仿真中。在MPM中，一个连续介质能被许多拉格朗日粒子来描述，这些粒子被称为物质点。这些物质点被一个背景网格包围，而该网格仅用于计算梯度项，如变形梯度。与其他基于网格的方法比如有限元法FEM、有限体积法FVM或者有限差分法FDM不同，MPM并不是一种基于网格的方法，而是一种无网格或者说是连续粒子法，因为与粒子法的结合，尽管存在背景网格，MPM并没有遇到基于网格方法(如平流项误差)的缺点，是计算力学中的一个很有前途以及强大的工具。

MPM最初是由新墨西哥大学的Deborah L.Sulsky等人在1990年提出的，是将FLIP方法进一步扩展到计算固体动力学的一种方法，而近些年来，MPM已经被应用于多种图形学模拟中，如2013年，Stomakhin等人使用MPM方法模拟雪，首次将MPM方法带入了大众的视野，其模拟的雪的动力学效果逼真，并应用于《冰雪奇缘》(Frozen)电影之中，产生了巨大的经济效益。基于上述的工作，Daviet等人于2016年用MPM方法引入Drucker-Prager屈服面来模拟沙粒的流动，Ram等人用MPM方法来模拟泡沫以及非牛顿流体。Stomakhin使用MPM来模拟固态以及液体物体之间相变以及Pradhana与Tampuboloneta等人用MPM方法来模拟多孔沙-水混合物的交互运动。

图形学领域，流固耦合方法的使用往往伴随着流体模拟。从计算对象上可以分为固体对流体单向耦合(one-way solid-to-fluid coupling)，流体对固体单向耦合(one-wayfluid-to-solid coupling)以及双向耦合(two-way coupling)。在固体对流体单向耦合方法中，固体按照设定的状态运动，流体不会对固体的运动路径产生影响。相对应的，流体对固体单向耦合方法中，通过对接触面上流体压力和扭转力的积分确定流体对固体的运动影响。然而这两种方法只能简单地作为固体或流体的边界条件，无法真实地模拟固体与流体交互。而双向耦合方法由于交互过程在物理上十分复杂，不存在解析解，不同的方法之间存在各种差异，适用范围小。

下面先介绍已有的流固双边耦合模拟方法：

1)任意拉格朗日-欧拉(Arbitrary Lagrangian-Eulerian)方法

任意拉格朗日-欧拉方法结合了传统的拉格朗日和欧拉方法，允许边界网格或面片顶点的任意运动，从而有效跟踪物质边界的运动，在内部网格上，根据定义参数求解网格，在独立于物质实体同时使得网格不发生严重的畸变。

2)内嵌边界(Immersed Boundary)方法

IB方法，通过在固体点所在的位置计算质量力来模拟流体域固体在边界的交互。质量力产生的运动约束被用来维持接触面上速度和应力的连续性。

发明内容

本发明的目的在于探索一种基于物理的山体滑坡灾害场景的模拟框架，其由物质点方法作为框架，并采用GPU加速的方法进行效率提升，并且综合了多学科优点，生成真实的山体滑坡模拟场景

基于物质点方法的山体滑坡灾害场景模拟包括以下步骤：

1)在计算机图形学领域对山体滑坡这一灾害现象进行建模，在建模过程的每一个时间步进行如下计算:将土粒子数据转移到网格，计算网格受力以及更新速度，处理网格碰撞，更新变形梯度，处理粒子碰撞，处理塑性流动以及硬化；其中在第一时间步还需计算粒子体积与密度；

2)在MPM方法的基础上改进屈服面准则，引入修改后的剑桥模型，通过计算平均主应力以及偏应力，将土粒子的变形梯度分为弹性变形梯度部分以及塑性变形梯度部分；再由软化硬化参数更新土粒子的屈服面；

3)使用水平集方法计算土粒子、网格与石块刚体的碰撞，计算得到土粒子的速度，从而计算出更新后的土粒子、网格位置以及刚体的角速度以及位置，得到MPM的颗粒状材质与刚体作用的双向耦合模型。

步骤1)中的1.3)计算网格受力以及更新速度中计算网格受力的离散化处理具体如下：

离散化处理，通过计算粒子上势能的和来离散地计算出势能：

上式中为粒子的体积，ψ为势能密度，由算法在第一个时间步的时候进行初始化。其中仅有变形梯度的弹性部分F_Ep会产生势能。

如果系统的状态能够被有限数量的位置来描述，可以类比为一系列质点通过弹簧连接在一起。若移动这些质点中的其中一个点会引起能量的变化(拉伸或者压缩弹簧需要能量)。此时弹簧将对这些质点产生力来释放积累的能量。通过这种方法，粒子j受到的力f_j可以描述为:

其中x_j为粒子j所处于的网格位置。基于应力的MPM空间离散化可以等价于该能量相对于欧拉网格结点材料位置的离散近似的微分。然而，实际上并没有变形欧拉网格，因此可以认为网格节点位置的变化是实际上由网格节点速度决定的。也就是说，如果当前网格结点的速度是v_i，结点变形后的位置可以被近似地表示为考虑到不同的结束位置意味着到达那里的速度不同。这些速度计算粒子的速度梯度：

同样，本发明把对总弹性势能的MPM近似写为:

上式中，是与点x^*相关的变形梯度F的弹性部分，其计算公式为:

综上，可以推导出基于MPM的计算网格结点力所需要的计算公式为:

还要注意，所有的变形被假设为弹性的，当计算该力时，塑性流动的影响被忽略。

重力产生的力以相同的方法被离散为:

步骤1)中1.3)计算网格受力以及更新速度中的速度更新详细如下

理论上，计算网格的弹塑性变化应该在欧拉网格结点进行。然而，正如之前所述，实际上网格并没有进行变形，因此，将x^*＝x^*(v)定义为与速度v值相关的函数。通过这样的设置，可以得到；

f_i ⁿ＝f_i(x^*(0))

f_i ⁿ⁺¹＝f_i(x^*(vⁿ⁺¹))

通过上面式子定义的求导方程，可以定义隐式更新为:

因此，通过上式可以推导出质量堆成的

其中上式方程的右侧项为:

在上面两式中，参数β用于控制速度更新方式，其中当β取值为0时，表示速度以显式方式进行更新，当β取值为1/2时，表示速度以梯形半隐式的方法进行更新，当β取值为1时，表示速度以后向欧拉的方式隐式更新。

在更新完粒子位置后，本发明对之前计算的弹性变形梯度F_E进行处理，进行塑性流动部分处理以及应用材质的硬化以及软化性质，使得每一个粒子都拥有属于自己的独一无二的屈服面。塑性流动是通过将变形梯度F投影到屈服面上(定义此投影为函数PJ())来计算弹性变形部分，然后剩余的部分就是塑性部分。并且由于塑性部分不应该改变总的变形梯度F,因此可以通过下面两式来更新塑性部分:

然后使用更新变形梯度

最后用于改变硬化参数的方程为式中，为时间步长为n+1时的参数p_c，为时间步长为n时的参数p_c，其中dε为塑性应变，而k为材料的硬化参数，其控制在塑性变形发生时，屈服面变化程度的大小；为了使硬化函数能够同时适用于软化以及硬化现象，定义dε为：

本发明的优点在于：

流固耦合问题涵盖了计算流体力学、数学、计算机等多门学科，往往包含大量的物理性质和数学运算，需要进行合理准确的抽象和建模，并通过计算机高效准确地模拟仿真。

基于土力学、材料学与计算机图形学等多学科的知识，首次在计算机图形学领域对山体滑坡这一灾害现象进行了基于物理的建模与绘制，这在灾难防治与救援、影视特效等领域都有较大的应用前景；

在传统MPM方法的基础上改进了屈服面准则，引入修改后的剑桥模型，从而能更为真实地模拟山体滑坡中土颗粒流的不同硬化特性，能模拟土石流的流动及堆积等特性效果；

提出了基于MPM的颗粒状材质与刚体作用的Level Set双向耦合模型，使得山体滑坡过程中的土石相互作用的效果更加真实，并采用GPU对碰撞检测过程进行了加速，提高了模拟的效率。

最后，本发明将上述方法整合，形成了一个山体滑坡灾难场景的模拟方法，成功解决了计算机图形学领域对此涉及较少的问题。说明书最后对实际山体滑坡场景进行了模拟，展示出了本方法的场景真实感。

附图说明

图1本发明生成的山体滑坡场景不同阶段展示；

图2本发明的步骤展示；

图3本发明流固耦合场景模拟结果图；

图4本发明(下行)与SPH方法(上行)以及传统颗粒流方法(中行)对比图；

图5本发明(右)与传统颗粒流方法(左)细节对比；

图6本发明实际应用：山体滑坡灾难场景模拟。

具体实施方式

基于物质点方法的山体滑坡灾害场景模拟，其具体实施方法在于包括以下步骤

1)在计算机图形学领域对山体滑坡这一灾害现象进行建模，在建模过程的每一个时间步进行如下计算:将土粒子数据转移到网格，计算网格受力以及更新速度，处理网格碰撞，更新变形梯度，处理粒子碰撞，处理塑性流动以及硬化；其中在第一时间步还需计算粒子体积与密度；

2)在MPM方法的基础上改进屈服面准则，引入修改后的剑桥模型，通过计算平均主应力以及偏应力，将土粒子的变形梯度分为弹性变形梯度部分以及塑性变形梯度部分；再由软化硬化参数更新土粒子的屈服面；

3)使用水平集方法计算土粒子、网格与石块刚体的碰撞，计算得到土粒子的速度，从而计算出更新后的土粒子、网格位置以及刚体的角速度以及位置，得到MPM的颗粒状材质与刚体作用的双向耦合模型。

所述的步骤1)具体为：

初始化土粒子以及场景信息，在每一时间步，迭代计算土粒子的位置与速度，具体如下：

1.1)将土粒子数据转移到网格

将土粒子的质量和速度从土粒子转移到空间网格,质量数据通过权重函数来进行转移其中为第n时间步的网格质量、m_p为粒子质量、为权重函数；速度数据使用正则化后的权重函数来进行转移:其中为第n时间步的网格体积、为第n时间步的的粒子体积；判断当前时间步是否为第一个时间步，若为第一个时间步，执行步骤1.2)，否则直接执行步骤1.3)；

1.2)在第一个时间步，计算粒子体积与密度

假设一个网格的一个细胞的密度为因此可以计算出粒子的密度为其中为粒子密度,h为粒子与网格的距离,从而估计一个粒子的体积为其中为粒子的体积，m_p为粒子密度；

1.3)计算网格受力以及更新速度

通过上一时间步由屈服面准则计算得到的弹性变形梯度进行受力计算，应用 进行速度计算更新；其中ρ为密度，t为时间，v为速度，σ为柯西应力项，g为重力加速度，Ψ为弹塑性势能密度,F_E为变形梯度F的弹性部分以及J＝det(F)；

1.4)处理网格碰撞

通过处理与地形边界、刚体之间的碰撞，更新速度；先计算网格相对于碰撞刚体的相对速度:v_irel＝v_iOri-v_icoll，式中，v_irel为网格相对于碰撞刚体的速度，v_iOri为网格的速度，v_icoll为与网格碰撞刚体的速度；如果刚体与网格是分离的，对应于法向速度v_in大于等于0的情形，此时没有碰撞发生，其中法向速度v_in的计算公式为:v_in＝v_irel·n，其中n为碰撞体的法向方向，可以由水平集φ的梯度方向来确定:n＝▽φ，若是v_in小于0的情形，此时网格与刚体发生碰撞，此时网格的切向速度的计算公式为:v_ti＝v_irel-nv_in，根据库伦-摩尔第一定律，若是||v_it||≤-μ_fv_in时，网格的作用力不能超过静摩擦力，因此网格不发生位移，其中μ_f为刚体表面的摩擦系数；否则，网格能够克服摩擦力运动，此时摩擦力表现为动摩擦，更新后的速度为：

1.5)更新变形梯度

通过式来更新变形梯度F；其中为前一时间步长计算的变形梯度，x为网格位置，I为单位矩阵；

1.6)将速度从网格转移到粒子；

采用FLIP与PIC混合的方法，粒子更新后的速度为其中PIC部分为FLIP部分为a为典型值0.95；上标n表示第n时间步，上标n+1表示第n+1时间步；

1.7)处理粒子碰撞

先计算粒子相对于碰撞刚体的相对速度:v_prel＝v_pOri-v_pcoll，式中，v_prel为粒子相对于碰撞刚体的速度，v_pOri为粒子速度，v_pcoll为碰撞刚体的速度；如果刚体与粒子是分离的，对应于法向速度v_pn大于等于0的情形，此时没有碰撞发生，其中法向速度v_pn的计算公式为:v_pn＝v_prel·n，若是v_pn小于0的情形，此时粒子与刚体发生碰撞，此时粒子的切向速度的计算公式为:v_pt＝v_prel-nv_pn，根据库伦-摩尔第一定律，若是||v_pt||≤-μ_fv_pn时，粒子的作用力不能超过静摩擦力，因此粒子不发生位移；否则，粒子能够克服摩擦力运动，此时摩擦力表现为动摩擦，更新后的速度为：

1.8)处理塑性流动以及硬化

将1.5)中得到的变形梯度F通过屈服面准则进行投影，更新其弹性部分以及塑性部分；并根据硬化软化参数来更新每一个粒子的屈服面，其中式中，为时间步长为n+1时的硬化软化参数p_c，为时间步长为n时的硬化软化参数p_c，其中dε为塑性应变，而k为材料的硬化参数，其控制在塑性变形发生时，屈服面变化程度的大小。

所述的步骤2)具体为：

2.1)对于场景中的粒子，其受应力为σ，则存在平均主应力p，以及偏应力q,其计算公式为：

其中，σ₁、σ₂、σ₃为三个主方向上的所受应力；

对土粒子构成的土样的持续剪切最终会导致在不改变应力或体积的情况下发生进一步的剪切，土壤在恒定应力状态下变形，体积不变；这种状态称为临界状态，以临界状态线CSL为特征，在p-q空间上CSL是通过原点的斜率等于M的直线来表现；剑桥模型的屈服函数为：式中F_Cy为剑桥模型的屈服函数，M为临界状态线CSL的斜率，p，q分别为平均主应力以及偏应力,而修正剑桥模型的屈服函数为：F_My＝q²+M²p(p+p_c)＝0，式中F_My为修正剑桥模型的屈服函数，p_c为硬化软化参数；在p-q空间中，剑桥模型的屈服面是对数曲线而修正剑桥模型表现为一个椭球，其中参数p_c控制屈服面的大小，参数M是CSL在p-q空间的斜率，并且CSL线与修正剑桥模型的屈服面相交于椭圆曲线q值最大处；

2.2)材料的硬化是由于塑性体积应变或者材料的压实作用造成的，根据屈服面理论，材料在超过屈服面的过程中，会发生屈服面的移动，具体表现为，硬化会导致材料可能在屈服之后需要受到更大的应力作用才有可能在新的点屈服；对于修正剑桥模型，简化其硬化函数，将其定义为：式中，为时间步长为n+1时的参数p_c，为时间步长为n时的参数p_c，其中dε为塑性应变，而k为材料的硬化参数，其控制在塑性变形发生时，屈服面变化程度的大小；为了使硬化函数能够同时适用于软化以及硬化现象，定义dε为：

由上式可知，如果屈服发生的横坐标值小于临界状态线CSL与屈服面交点的横坐标值时，材料发生软化，此时材料在发生塑性流动的同时，屈服面发生变化，p_c值减小，材料在更小的应力条件下就会发生屈服；如果屈服发生在临界状态线CSL与屈服面交点的右侧时，材料发生硬化，此时材料在发生塑性流动的同时，屈服面发生变化，p_c值增大，材料将在更大的应力条件下才会发生屈服。

所述的步骤3)具体为：

3.1)土粒子与石块刚体双向耦合以及摩擦力的作用：对于每一时间步，对碰撞刚体处理两次碰撞；第一次处理碰撞是在当力作用于网格，网格速度刚更新完的时候，第二次处理碰撞是在粒子更新位置之前，再次应用碰撞来消除粒子和网格速度的微小差异；在两次碰撞中，碰撞操作采用相同的方法来处理，并且所有的碰撞体都将其视为刚体，将碰撞刚体以水平集φ的方式进行处理，判断粒子是否与刚体碰撞；

3.2)由上步可知，当水平集φ的值为负时，说明粒子与刚体发生了碰撞，首先计算粒子或者网格相对于碰撞刚体的相对速度:v_iprel＝v_ipOri-v_ipcoll，式中，v_iprel为粒子相对于碰撞刚体的速度，v_ipOri为粒子或者网格的速度，v_ipcoll为碰撞刚体的速度；如果刚体与粒子或者网格是分离的，对应于法向速度v_ipn大于等于0的情形，此时没有碰撞发生，其中法向速度v_ipn的计算公式为:v_ipn＝v_iprel·n，，若是v_ipn小于0的情形，此时粒子或者网格与刚体发生碰撞，此时粒子或者网格的切向速度的计算公式为:v_ipt＝v_iprel-nv_ipn，根据库伦-摩尔第一定律，若是||v_ipt||≤-μ_fv_ipn时，粒子或者网格的作用力不能超过静摩擦力，因此粒子或者网格不发生位移，其中μ_f为刚体表面的摩擦系数；否则，粒子或者网格能够克服摩擦力运动，此时摩擦力表现为动摩擦，更新后的速度为：

3.3)根据动量守恒，同时计算出粒子或者网格对于刚体的平移速度以及角速度：v_refRIG＝v_lineRIG+ω_angleRIG×(x-x_mass)，其中v_lineRIG为粒子或者网格对刚体影响的平移速度，ω_angleRIG是粒子或者网格对刚体作用的角速度，x_mass为粒子或者网格的质量，x为刚体的质量，v_refRIG为粒子或者网格对刚体影响的相对速度。

所述的步骤3.2)中，计算粒子或者网格与刚体碰撞的的具体实现为：

a)由于SDFs的预计算是在刚体空间进行的，需要先对SDFs做刚体坐标到世界坐标的变换，使其处于同一坐标系，再进行接下来的操作，把离散后的SDFs以f_SDF来表示，通过纹理内存的性质，可以通过硬件加速的GPU插值办法，将SDFs进行插值，可以生成更为圆滑的碰撞效果，由于SDFs存在与等值面相互垂直的梯度，因此可以把梯度等同于曲面的法线，通过取值需要采样点位置上下、前后、左右三轴上的6个点来计算出曲面的法线方向，并且得益于纹理内存的局部性原则，需取值位置的附近纹理也同时被读入缓存中，因此访问这些纹理并不需要多少计算量，其计算公式为下:

f_x＝f_SDF(point+(step,0,0))-f_SDF(point-(step,0,0))

f_y＝f_SDF(point+(0,step,0))-f_SDF(point-(0,step,0))

f_z＝f_SDF(point+(0,0,step))-f_SDF(point-(0,0,step))

normal＝normalize(f_x,f_y,f_z)

上面式子中，normal为表面法向,f_x,f_y,f_z为采样点相对于碰撞点在x,y,z三轴上的偏移量，f_SDF为SDF函数，point为土粒子与刚体碰撞点的位置，计算完法线后，便可以计算粒子与刚体之间的相互耦合；

b)考虑到自然灾害场景中存在山体地形与土石的碰撞，为了模拟大规模的地形，采用高度图加上法向贴图的方法进行处理；高度图是一张二维单通道灰度图，其像素值表示的是该点的高程数据，其映射为0.0-1.0；将法向贴图以及高度图存储成一个三通道的纹理，其中r，g通道存储法向贴图的前两个通道，b通道存储高度图中的灰度信息，因此，此处仅需要取一次纹值，便可以得到碰撞信息以及法向信息。

通过以上步骤，图1-6给出了本发明实际应用模拟及对比效果图。图1为本发明在设定初始土粒子位置、刚体石块位置以及设置场景的高度地形的情况下由本发明方法模拟的山体滑坡场景不同阶段的展示；图2为本发明的具体实现步骤展示；图3为本发明方法模拟的刚体石球落入土粒子中的流固耦合场景模拟结果图；图4为本发明与SPH方法以及传统颗粒流方法模拟效果的对比图，场景为相同粒子数目的土堆从45度角的坡滑落的情形；图5为本发明(右)与传统颗粒流方法(左)表面堆积的细节对比，本场景为土流从上而下流动而产生的堆积效果；本图6为本发明在设定初始土粒子位置、刚体石块位置以及设置场景的高度地形的情况下由本发明方法模拟的山体滑坡运动过程各阶段的图象：其中：(a)为滑坡前的景象；(b)为滑坡刚开始阶段的落石和裂隙出现；(c)为山体开始崩落；(d)为山体崩塌继续；(e)为滑坡流开始吞没建筑物；(f)为滑坡流冲向公路上的汽车；(g)为滑坡流吞没汽车；(h)为；滑坡流最终堆积成型。

以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有许多变形。本领域的普通技术人员能从本发明公开的内容直接导出或联想到的所有变形，均应认为是本发明的保护范围。

Claims (5)

1.一种基于物质点方法的山体滑坡灾害场景模拟方法，其特征在于包括以下步骤：

1)在计算机图形学领域对山体滑坡这一灾害现象进行建模，在建模过程的每一个时间步进行如下计算:将土粒子数据转移到网格，计算网格受力以及更新速度，处理网格碰撞，更新变形梯度，处理粒子碰撞，处理塑性流动以及硬化；其中在第一时间步还需计算粒子体积与密度；

2)在MPM方法的基础上改进屈服面准则，引入修改后的剑桥模型，通过计算平均主应力以及偏应力，将土粒子的变形梯度分为弹性变形梯度部分以及塑性变形梯度部分；再由软化硬化参数更新土粒子的屈服面；

3)使用水平集方法计算土粒子、网格与石块刚体的碰撞，计算得到土粒子的速度，从而计算出更新后的土粒子、网格位置以及刚体的角速度以及位置，得到MPM的颗粒状材质与刚体作用的双向耦合模型。

2.根据权利要求1所述的一种基于物质点方法的山体滑坡灾害场景模拟方法，其特征在于所述的步骤1)为：

初始化土粒子以及场景信息，在每一时间步，迭代计算土粒子的位置与速度，具体如下：

1.1)将土粒子数据转移到网格

将土粒子的质量和速度从土粒子转移到空间网格,质量数据通过权重函数来进行转移其中为第n时间步的网格质量、m_p为粒子质量、为权重函数；速度数据使用正则化后的权重函数来进行转移:其中为第n时间步的网格体积、为第n时间步的的粒子体积；判断当前时间步是否为第一个时间步，若为第一个时间步，执行步骤1.2)，否则直接执行步骤1.3)；

1.2)在第一个时间步，计算粒子体积与密度

假设一个网格的一个细胞的密度为因此可以计算出粒子的密度为其中为粒子密度,h为粒子与网格的距离,从而估计一个粒子的体积为其中为粒子的体积，m_p为粒子密度；

1.3)计算网格受力以及更新速度

通过上一时间步由屈服面准则计算得到的弹性变形梯度进行受力计算，应用 进行速度计算更新；其中ρ为密度，t为时间，v为速度，σ为柯西应力项，g为重力加速度，Ψ为弹塑性势能密度,F_E为变形梯度F的弹性部分以及J＝det(F)；

1.4)处理网格碰撞

通过处理与地形边界、刚体之间的碰撞，更新速度；先计算网格相对于碰撞刚体的相对速度:v_irel＝v_iOri-v_icoll，式中，v_irel为网格相对于碰撞刚体的速度，v_iOri为网格的速度，v_icoll为与网格碰撞刚体的速度；如果刚体与网格是分离的，对应于法向速度v_in大于等于0的情形，此时没有碰撞发生，其中法向速度v_in的计算公式为:v_in＝v_irel·n，其中n为碰撞体的法向方向，可以由水平集φ的梯度方向来确定:n＝▽φ，若是v_in小于0的情形，此时网格与刚体发生碰撞，此时网格的切向速度的计算公式为:v_ti＝v_irel-nv_in，根据库伦-摩尔第一定律，若是||v_it||≤-μ_fv_in时，网格的作用力不能超过静摩擦力，因此网格不发生位移，其中μ_f为刚体表面的摩擦系数；否则，网格能够克服摩擦力运动，此时摩擦力表现为动摩擦，更新后的速度为：

1.5)更新变形梯度

通过式来更新变形梯度F；其中为前一时间步长计算的变形梯度，x为网格位置，I为单位矩阵；

1.6)将速度从网格转移到粒子；

采用FLIP与PIC混合的方法，粒子更新后的速度为其中PIC部分为FLIP部分为a为典型值0.95；上标n表示第n时间步，上标n+1表示第n+1时间步；

1.7)处理粒子碰撞

先计算粒子相对于碰撞刚体的相对速度:v_prel＝v_pOri-v_pcoll，式中，v_prel为粒子相对于碰撞刚体的速度，v_pOri为粒子速度，v_pcoll为碰撞刚体的速度；如果刚体与粒子是分离的，对应于法向速度v_pn大于等于0的情形，此时没有碰撞发生，其中法向速度v_pn的计算公式为:v_pn＝v_prel·n，若是v_pn小于0的情形，此时粒子与刚体发生碰撞，此时粒子的切向速度的计算公式为:v_pt＝v_prel-nv_pn，根据库伦-摩尔第一定律，若是||v_pt||≤-μ_fv_pn时，粒子的作用力不能超过静摩擦力，因此粒子不发生位移；否则，粒子能够克服摩擦力运动，此时摩擦力表现为动摩擦，更新后的速度为：

1.8)处理塑性流动以及硬化

将1.5)中得到的变形梯度F通过屈服面准则进行投影，更新其弹性部分以及塑性部分；并根据硬化软化参数来更新每一个粒子的屈服面，其中式中，为时间步长为n+1时的硬化软化参数p_c，为时间步长为n时的硬化软化参数p_c，其中dε为塑性应变，而k为材料的硬化参数，其控制在塑性变形发生时，屈服面变化程度的大小。

3.根据权利要求1所述的一种基于物质点方法的山体滑坡灾害场景模拟方法，其特征在于所述的步骤2)为：

2.1)对于场景中的粒子，其受应力为σ，则存在平均主应力p，以及偏应力q,其计算公式为：

其中，σ₁、σ₂、σ₃为三个主方向上的所受应力；

对土粒子构成的土样的持续剪切最终会导致在不改变应力或体积的情况下发生进一步的剪切，土壤在恒定应力状态下变形，体积不变；这种状态称为临界状态，以临界状态线CSL为特征，在p-q空间上CSL是通过原点的斜率等于M的直线来表现；剑桥模型的屈服函数为：式中F_Cy为剑桥模型的屈服函数，M为临界状态线CSL的斜率，p，q分别为平均主应力以及偏应力,而修正剑桥模型的屈服函数为：F_My＝q²+M²p(p+p_c)＝0，式中F_My为修正剑桥模型的屈服函数，p_c为硬化软化参数；在p-q空间中，剑桥模型的屈服面是对数曲线而修正剑桥模型表现为一个椭球，其中参数p_c控制屈服面的大小，参数M是CSL在p-q空间的斜率，并且CSL线与修正剑桥模型的屈服面相交于椭圆曲线q值最大处；

2.2)材料的硬化是由于塑性体积应变或者材料的压实作用造成的，根据屈服面理论，材料在超过屈服面的过程中，会发生屈服面的移动，具体表现为，硬化会导致材料可能在屈服之后需要受到更大的应力作用才有可能在新的点屈服；对于修正剑桥模型，简化其硬化函数，将其定义为：式中，为时间步长为n+1时的参数p_c，为时间步长为n时的参数p_c，其中dε为塑性应变，而k为材料的硬化参数，其控制在塑性变形发生时，屈服面变化程度的大小；为了使硬化函数能够同时适用于软化以及硬化现象，定义dε为：

由上式可知，如果屈服发生的横坐标值小于临界状态线CSL与屈服面交点的横坐标值时，材料发生软化，此时材料在发生塑性流动的同时，屈服面发生变化，p_c值减小，材料在更小的应力条件下就会发生屈服；如果屈服发生在临界状态线CSL与屈服面交点的右侧时，材料发生硬化，此时材料在发生塑性流动的同时，屈服面发生变化，p_c值增大，材料将在更大的应力条件下才会发生屈服。

4.根据权利要求1所述的一种基于物质点方法的山体滑坡灾害场景模拟方法，其特征在于所述的步骤3)为：

3.1)土粒子与石块刚体双向耦合以及摩擦力的作用：对于每一时间步，对碰撞刚体处理两次碰撞；第一次处理碰撞是在当力作用于网格，网格速度刚更新完的时候，第二次处理碰撞是在粒子更新位置之前，再次应用碰撞来消除粒子和网格速度的微小差异；在两次碰撞中，碰撞操作采用相同的方法来处理，并且所有的碰撞体都将其视为刚体，将碰撞刚体以水平集φ的方式进行处理，判断粒子是否与刚体碰撞；

3.2)由上步可知，当水平集φ的值为负时，说明粒子与刚体发生了碰撞，首先计算粒子或者网格相对于碰撞刚体的相对速度:v_iprel＝v_ipOri-v_ipcoll，式中，v_iprel为粒子相对于碰撞刚体的速度，v_ipOri为粒子或者网格的速度，v_ipcoll为碰撞刚体的速度；如果刚体与粒子或者网格是分离的，对应于法向速度v_ipn大于等于0的情形，此时没有碰撞发生，其中法向速度v_ipn的计算公式为:v_ipn＝v_iprel·n，，若是v_ipn小于0的情形，此时粒子或者网格与刚体发生碰撞，此时粒子或者网格的切向速度的计算公式为:v_ipt＝v_iprel-nv_ipn，根据库伦-摩尔第一定律，若是||v_ipt||≤-μ_fv_ipn时，粒子或者网格的作用力不能超过静摩擦力，因此粒子或者网格不发生位移，其中μ_f为刚体表面的摩擦系数；否则，粒子或者网格能够克服摩擦力运动，此时摩擦力表现为动摩擦，更新后的速度为：

3.3)根据动量守恒，同时计算出粒子或者网格对于刚体的平移速度以及角速度：v_refRIG＝v_lineRIG+ω_angleRIG×(x-x_mass)，其中v_lineRIG为粒子或者网格对刚体影响的平移速度，ω_angleRIG是粒子或者网格对刚体作用的角速度，x_mass为粒子或者网格的质量，x为刚体的质量，v_refRIG为粒子或者网格对刚体影响的相对速度。

5.根据权利要求4所述的一种基于物质点方法的山体滑坡灾害场景模拟方法，其特征在于所述的步骤3.2)中，计算粒子或者网格与刚体碰撞的的具体实现为：

a)由于SDFs的预计算是在刚体空间进行的，需要先对SDFs做刚体坐标到世界坐标的变换，使其处于同一坐标系，再进行接下来的操作，把离散后的SDFs以f_SDF来表示，通过纹理内存的性质，可以通过硬件加速的GPU插值办法，将SDFs进行插值，可以生成更为圆滑的碰撞效果，由于SDFs存在与等值面相互垂直的梯度，因此可以把梯度等同于曲面的法线，通过取值需要采样点位置上下、前后、左右三轴上的6个点来计算出曲面的法线方向，并且得益于纹理内存的局部性原则，需取值位置的附近纹理也同时被读入缓存中，因此访问这些纹理并不需要多少计算量，其计算公式为下:

f_x＝f_SDF(point+(step,0,0))-f_SDF(point-(step,0,0))

f_y＝f_SDF(point+(0,step,0))-f_SDF(point-(0,step,0))

f_z＝f_SDF(point+(0,0,step))-f_SDF(point-(0,0,step))

normal＝normalize(f_x,f_y,f_z)

上面式子中，normal为表面法向,f_x,f_y,f_z为采样点相对于碰撞点在x,y,z三轴上的偏移量，f_SDF为SDF函数，point为土粒子与刚体碰撞点的位置，计算完法线后，便可以计算粒子与刚体之间的相互耦合；

b)考虑到自然灾害场景中存在山体地形与土石的碰撞，为了模拟大规模的地形，采用高度图加上法向贴图的方法进行处理；高度图是一张二维单通道灰度图，其像素值表示的是该点的高程数据，其映射为0.0-1.0；将法向贴图以及高度图存储成一个三通道的纹理，其中r，g通道存储法向贴图的前两个通道，b通道存储高度图中的灰度信息，因此，此处仅需要取一次纹值，便可以得到碰撞信息以及法向信息。