牛眠沟地震滑坡碎屑化全过程离散元模拟

毕钰璋，付跃升，何思明，吴 永

(1.福州大学紫金矿业学院，福建福州 350116;2.中国科学院山地灾害与地表过程重点实验室，四川成都 610041;3.中国科学院水利部成都山地灾害与环境研究所，四川成都 610041)

0 引言

地震是引起地表大面积移动以及诱发大规模滑坡的主要因素，其引起滑坡的范围甚至可以达到100000 km2［1-2］，给国家经济以及人民财产安全造成重大损失，滑坡灾害在水利水电工程、矿业工程、交通工程等方面也存在很大隐患［3-6］，因而研究此类灾害有着深远的意义。滑坡碎屑化是指高速远程滑坡或崩塌在运动过程中解体破碎形成的碎屑流［7-10］。在过去十年中，国内外学者对引发滑坡的主要因素(诸如材料强度、外力、孔隙流体压力等)作了详尽的研究［11-13］。

2008年的汶川地震引发了大量的滑坡灾害，截止到2008年6月18号，四川省确定的汶川地震诱发的山体滑坡灾害就有4504起［14］。其中有至少23起滑坡引发的碎屑流灾害造成了20000人次的人员伤亡，其中最多一次致死30余人［15］。

针对地震引发的滑坡灾害，国内外学者从各个方面对此作出了详细的研究。国外学者Sato根据滑坡的总体积，将地震诱发的山体滑坡分为岩崩和岩滑［16］。Hungr根据离散元数值模拟，采用DAN对滑坡机理及碎屑流成因进行了详尽的阐述［17］。He和Li使用PFC2D研究了东河口堰塞湖的形成机制，得出了滑程与摩擦系数的关系［18］。Tang［19］在研究台湾集集地震引发的草岭山体碎屑流灾害中，提出了大型远程滑坡的摩擦系数决定于滑坡的垂直位移和水平位移之比，并认为此摩擦系数远远小于材料试验所得的摩擦系数。东京大学的Wang团队在调查了大型地震诱发的山体滑坡之后总结了以下三点结论［20］:(1)所有滑坡都是在很短的时间内发生，其滑移速度相当迅速;(2)大型远程滑坡主要发生在板岩、泥岩和页岩的地质构造条件下，岩崩型滑坡则主要形成于白云石的地质条件下;(3)水文地质条件是引发远程滑坡的重要因素，尤其是在滑坡地质周围有河流或沟谷。Huang根据地震诱发的滑坡灾害研究提出了以下两点［21］:(1)90%的滑坡是发生在地震断层带30 km的范围内;(2)虽然随着滑坡体积的增大，等效摩擦系数也趋向于减少滑坡动能，但是大部分的大型滑坡仍然具有相当大的动能。

滑坡碎屑流的特点有:规模大、运动快、滑程远，其运动过程受到地质、地形以及流化作用等多方面原因的影响。基于各方面原因，室内试验根本无法满足且还原其运动过程，因此数值模拟就成为研究滑坡碎屑流的主要手段。由于碎屑流的破坏能力往往和它的滑程息息相关，为了更好地研究其破坏机理以及分析其滑程，就需要建立适当的连续模型或不连续模型。许多学者致力于有限单元法来模拟大规模的滑坡灾害［19-21］，但是滑坡碎屑流中含有大量成分复杂的颗粒物质，例如岩块、土块、碎石等等。正是基于有限元本身的局限性，因而不能作为研究滑坡碎屑流的有效手段。然而，近年来基于离散元方法的数值模拟的迅猛发展给解决此问题带来了可能。特别是对于地震诱发的滑坡灾害的研究过程中，已被证明是有效的途径。大量学者运用此方法就岩石破坏机理、边坡稳定性分析以及滑坡灾害动力演化等多方面问题进行了详尽的研究［19-20，25-27］。

1 牛眠沟地震滑坡灾害

牛眠沟滑坡位于成都以北73 km的北川县映秀镇(图1、图2)，是汶川地震引发的三大滑坡碎屑流灾害之一［32］，巨型的滑体从海拔约1.8 km的何家山村西山头以雷霆万钧之势瞬间下滑，并且造成了沟内房屋、良田悉数被毁，以及近50余人伤亡的巨大损失。

牛眠沟滑坡区地貌类型以高山为主并且地形高差大，其后缘高程约1.79 km，前缘高程约0.89 km(图3)，相对高差0.87 km。此次滑坡碎屑流的横向移动距离约3.0 km。最终在距离岷江0.25 km处停积下来。沟谷发育，其区内的侵蚀强烈。斜坡的平均坡度范围为30°～40°，局部坡度甚至可达70°。

图1 (a)显示了牛眠沟滑坡所在的具体地理位置(b)反应了研究区域的的地形地貌，红线所围区域为牛眠沟滑坡碎屑流灾害区，红色虚线为断裂带［15］Fig．1 (a)Location map of Niumiangou Creek rock avalanche(b)Topography of the study area，red line is the boundary of the rock avalanche;red dashed line is Main Fault［15］

图2 牛眠沟滑坡—碎屑流主滑方向的剖面图，下部碎屑流堆积是在残余摩擦系数为0.1的情形Fig．2 Longitudinal cross-section of the Niumiangou Creek rock avalanche，the accumulation mass is depicted in the situation of the residual friction coefficient equals 0.1

滑坡位于龙门山华夏系构造体系中的九顶山华夏系构造带之内，滑坡区内主要断层是映秀——北川断层(图1)。断层两侧的地层岩性明显不同:映秀——北川断层的上盘为晋宁期花岗岩，且滑坡启动区多为灰白色斜长花岗岩，其表部节理极度发育;映秀——北川断层的下盘主要为中生界三叠系上统须家河组(T3x2)，出露岩层多为泥质板岩和厚层砂岩，岩体内的裂隙发育，完整性差。

Yin对牛眠沟滑坡碎屑流运动过程进行了详细的阶段划分，分别是:滑坡启动与粉碎阶段，凌空飞行阶段、碎屑流阶段，此三个阶段发生的区域分别分布在滑源区、莲花心沟和牛眠沟(图2、图3)。并认为滑坡的失稳是因为水平方向力和垂直方向力耦合的结果。对牛眠沟滑坡碎屑流运动的远程效应，Yin［15］给出了如下的解释:(1)由于在凌空飞行的过程中，碎屑物质遭受了强大的气垫效应，从而使滑体飞行距离加大;(2)碎屑物质流化作用是导致其能在水平坡面运动千米的最主要原因。

对于汶川地震诱发的牛眠沟滑坡碎屑流的运动机理，已经有学者对其进行相关研究［20］。对于牛眠沟滑坡碎屑流的运动过程，也有学者对其进行了部分的探讨，但其研究范围只是注重于滑坡碎屑流的堆积阶段，并未完整地反映其整体的运动过程。基于此，本文的主要目的就是运用离散元数值模拟的手段来研究牛眠沟滑坡碎屑流的运动机理并力求还原其运动过程。尤其是对于滑坡模型力学参数之间关系的研究，例如滑动面的摩擦、滑坡体内聚力和滑坡位移的关系。

图3 研究区域的牛眠沟滑坡碎屑流航拍影像［33］Fig．3 Remote sensing image of the Niumiangou Creek rock avalanche

2 牛眠沟地震滑坡碎屑化过程的PFC2D模型

离散元方法是由Cundall和Strack［28］在1979年创立的一种基于二维圆盘和三维圆球排列建立数值模型的方法。该方法基于粒子间接触的位置关系以及相互作用的力学关系来求得粒子的运动状态。离散元最初被用于模拟颗粒材料，但是随着其本构模型的发展，粘结模型和平行接触模型的应用使得离散元方法能用于研究固体材料的力学特性中去。本文数值模型的建立，就是基于上述本构模型的二维颗粒流方法(PFC)。PFC作为离散元方法中比较成熟的一种数值模拟手段，可以通过颗粒的固结模型较好地用于模拟材料的静力特性和动力特性。

对于粘结模型而言，颗粒的力学参数包括了正应力、剪应力、粘结的法向刚度、粘结的切向刚度以及库伦摩擦系数等等。这些颗粒的力学参数很大程度上决定了颗粒间的力学作用机理。因此颗粒的力学参数的合理选择，就显得十分重要。

2.1 牛眠沟地震滑坡碎屑化模型

本文所建立的牛眠沟滑坡模型的尺寸为:长3 500 m，高1 090 m。其中滑坡体下部的岩石数值模型采用 PFC2D中的“墙单元”来表示，滑坡体用PFC2D中的众多随机排列的颗粒模型所组成。通过对滑坡的原位观测，我们设定了滑坡模型的颗粒半径范围为1.28到2.56 m，颗粒的总数为3647个。通过设定颗粒间摩擦系数以及和墙体的摩擦系数变化范围为0.1～0.5，我们确定了岩块的残余摩擦系数变化范围。

Zhou［29］通过能量守恒演算出了牛眠沟滑坡的运动摩擦系数0.24以及滑坡启动的初始速度为67 m/s，He［14］在使用PFC2D研究东河口滑坡的过程中，通过反复数值模拟试验确定了20m/s的滑坡初始速度符合实际滑坡的滑程。Yin［15］在研究牛眠沟滑坡碎屑流的动力过程中得出了碎屑流在距岷江约250m处停止，并给出了碎屑流凌空飞行阶段所在的区域范围(莲花心区)。在综合考虑这些学者的研究成果后，为了更好地还原牛眠沟滑坡碎屑流的运动过程，将滑坡分为四段并渲染不同颜色(图4)，这样便于观测滑坡的启动粉碎阶段。研究残余摩擦系数在0.1～0.5范围内的碎屑流滑移过程。

2.2 颗粒的力学参数

在离散元方法中，材料的宏观运动行为取决于颗粒间接触的力学参数。然而对于参数的选取却没

图4 牛眠沟滑坡碎屑流启动及粉碎化阶段数值模拟，不同时步滑坡体的断裂粉碎情况Fig．4 The different stages about the simulation of the failing and disintegrating of the Niumiangou Creek rock avalanche

有一个完全有效的方式。较多的研究方式是通过物理模型和数值模型的破坏模式相匹配，从而确定其力学参数。由于没有完备的理论依据来确定宏观条件下的运动过程和颗粒力学参数之间的选取关系，因此我们在前人对牛眠沟滑坡研究时参数选取的基础上［30］，通过反复的双轴压缩数值模型实验(图5右)来观测破坏模式，以期使其与实际的物理实验相匹配。在经过多次的数值模拟试验之后，确定了PFC2D的力学参数(表1)。

图5 左:根据实际地形剖面图建立的牛眠沟离散元模型，长为3500m，宽为1090m;右:双轴压缩实验的数值模拟离散元模型Fig．5 Left:DEM model of Niumiangou landslide．Length and height of model are 3500 and 1090 m，respectively，which is the actual field scale;Right:principal stresses exerted，in cross-sectioned cylindrical test sample．

2.3 地震荷载的选取和处理

关于PFC2D中地震动载荷的选取方法，也有学者给出了不同的思路［30-31］。考虑到模型的选取、地质剖面的走向，本文选取的地震加速度时程曲线为靠近牛眠沟的汶川卧龙台站当天的EW方向分量(图6(a))。图中反映了强震加速度的几个特点:频率快、峰值大、持续时间长，并且其地震加速度有两个波峰区。图6(b)是通过积分手段获得的地震速度随时间的变化曲线，在给模型赋予地震加速度时通过混合加载的方式，即在边界墙体和边界颗粒单元上施加图6(b)的速度时程曲线。

表1 通过数值模型的双轴压缩试验确定的颗粒力学参数Table 1 Microproperties macroproperties of models，as suggested by the results of a series of biaxial tests

图6(a):汶川地震加速度时程曲线;(b):通过积分得到的汶川地震速度时程曲线Fig．6 (a)Acceleration diagram during the Wenchuan earthquake．(b)Velocity diagram during the Wenchuan earthquake

3 牛眠沟地震滑坡碎屑化全过程模拟

经过数值模型的建立和参数的合理选定，对牛眠沟滑坡碎屑流的动力过程进行了数值模拟。并分别对启动阶段、凌空飞行阶段、碎屑流运动堆积阶段的数值模拟结果进行了讨论和分析。特别是对于碎屑流运动堆积阶段，我们不仅讨论了不同摩擦系数条件下碎屑流的运动距离，而且对牛眠沟地震滑坡的堆积特征进行了详细的研究。

3.1 启动及碎屑化阶段数值模拟

碎屑流的启动区位于莲花心沟的西南方向，地质成分主要是花岗岩。牛眠沟滑坡的滑坡源区厚度约为100 m，长300～400 m，宽同样约为300～400 m，总体积约为7.5×106m3。由于牛眠沟滑坡的破坏是受汶川地震水平作用力和纵向作用力共同作用的结果，并迅速拉开一组结构面使滑坡解体并高速滑出［15］。

模拟结果(图4)显示，在第5 s时滑坡体开始出现断裂。当进行到15 s时，滑坡体呈现明显断裂面，并被分为上部滑坡体和下部滑坡体，部分颗粒有喷射而出的现象。当进行到30 s时，滑坡体被断裂面解体为上、中、下三个部分，且下部滑坡体朝自由面喷射而出。

3.2 凌空飞行阶段数值模拟

许多学者证明了牛眠沟滑坡碎屑流在莲花心区是呈凌空飞行状态，并有学者经过实地考察给出了以下两点证明［15］:(1)莲花心沟的物质堆积不具有“反序”现象，相反的是，此处底部堆积物为大块堆积，上部为细粒物质堆积;(2)现场调查过程中，发现山体被撞出很大的一个豁口，证明碰撞体应具有巨大能量，而一般滑坡不具备如此之大的动能。有学者根据现场观测推测了牛眠沟滑坡碎屑流水平速度近于100m/s，垂直速度近于12m/s。

图7所示:本次数值模型试验模拟了碎屑流在莲花心区凌空飞行阶段的全过程，并对各个时步的凌空飞行范围予以划分。120s时下部滑坡体呈粉碎状，上部滑坡体和中部滑坡体部分成粉碎状。至150s时，碎屑状滑坡体高速运动并席卷至莲花心区的1/3处，此时仍有颗粒固结成的大块存在。至185s时，碎屑状滑坡体运动至莲花心区2/3处，伴随运动过程中颗粒碰撞碎屑程度加重。至225s时，碎屑状滑坡体运动至莲花心区和牛眠沟区交界处，即莲花心区末端。此时的数值模拟结果显示:莲花心区的坡面已有碎屑状颗粒堆积，仍保持凌空飞行的碎屑状颗粒已无明显固结大块存在。至350s时，凌空飞行阶段完全结束，碎屑流因保持了一定速度，将在重力作用下向牛眠沟区滑动。整个滑坡体约1/2被高速抛出，剩下的滑坡体在初速度以及重力作用下仍向下滑动。

3.3 摩擦系数对地震滑坡碎屑流运动距离的影响

图7 牛眠沟滑坡碎屑流凌空飞行阶段数值模拟，不同时步的碎屑流凌空飞行情况Fig．7 The different stages about the simulation of the airborne movement of the Niumiangou Creek rock avalanche

现场实地的考察证明了牛眠沟滑坡碎屑流最终在距离岷江250m处停止，最终停止处和岷江之间隔着都汶公路。此次滑坡灾害并未像东河口滑坡那样形成堰塞湖，但是我们通过数值模拟试验预测不同地质构造对滑坡碎屑流滑程影响对于预防灾害仍有一定的指导意义。根据现有的滑坡地质构造数据，我们可以通过数值模拟试验来研究坡面平均摩擦力以及坡体固结力对于滑程的影响。在相同峰值加速度条件下，设定了不同残余摩擦系数下的滑坡条件(图8)，其残余摩擦系数分别为:0.5，0.3，0.2，0.1。

研究显示了不同残余摩擦系数对于滑坡碎屑流滑程的影响是巨大的。当残余摩擦系数大于0.3时，没有碎屑流滑动至岷江。当残余摩擦系数为0.2时，岷江和都汶公路之间有明显的碎屑流堆积。因此得出实际摩擦系数应当在0.2～0.3的范围内，有学者通过研究得出其值为0.24［29］，跟其结果相对吻合。当残余摩擦系数为0.1时，都汶公路被碎屑流完全淹没，岷江有明显碎屑流堆积，其堆积体积约10%左右。

图8 不同残余摩擦系数条件下的滑坡碎屑流滑程Fig．8 Different runouts resulting from numerical tests，as a function of the different residual friction parameters between sliding surface and sliding mass

3.4 牛眠沟地震滑坡堆积特性分析

当摩擦系数为0.1的时候，灾害范围大，堆积范围广，并且都汶公路和岷江内均有大量的碎屑流堆积体存在。研究最坏情况下的碎屑流运动堆积情况，对于灾害防治工作的借鉴更有指导意义。因此为了研究滑坡体的滑移堆积规律，我们在残余摩擦系数为0.1的滑坡条件下做了以下两个工作:(1)将滑坡体沿垂直边坡方向分成颜色不同的四个条状(图4 t=0s)，宏观地观测不同条状的运动堆积情况。(2)对滑坡体设置12个监测点(图10)，以期对滑坡体各部分运动规律进行细化的研究。

如图9所示，是碎屑流的最终堆积结果图。我们可以看出:黄色和蓝色圆盘主要堆积在靠近碎屑流的龙头处，绿色圆盘则完全分布在靠近碎屑流的龙尾处。这是因为黄色和蓝色标记部分的滑坡体大部分被抛射出来，并且其破坏时初速度远大于其他部分的滑坡体。由此可以推测滑坡体深部的岩石主要运动方式是滑移为主。通过图9可以得出，浅部滑坡体运动的距离更远，深度滑坡体运动的距离较近。

图9 牛眠沟地震滑坡数值模拟的滑坡体最终堆积结果Fig．9 The deposit result of the Niumiangou Creek rock avalanche

如图10(a)所示，我们在滑坡体内设置了12个监测点，圆盘1到圆盘12之间的水平距离约为400m，堆积后的情况如图10(b)所示，其水平距离约为2000m。由于碎屑流凌空飞行阶段获得了很大的初速度，导致了部分监测点(圆盘 2、5、6、9、12)的运动距离远远大于其他监测点的运动距离。而在这五个监测点中，圆盘2的运动距离最远(约3100m)，这是因为在相同的动能条件下，圆盘2拥有更大的势能。五个观测点中的圆盘6并未遵循这一规律，应该是由于凌空飞行阶段和其他圆盘撞击而使其能量减少的原因。

对于没有凌空飞行阶段的滑坡体，其滑移堆积也满足这样的规律:圆盘3的水平运动距离最远(在圆盘1、3、4、7、8、10、11)，约为 1900m。圆盘 3 因其上部滑坡体高速喷射而出之后而获得自由面，由于没有上部滑坡体的摩擦，并且其拥有较大的势能，故而使其水平运动距离相对其他6个监测点更远。

图10 设置的12个监测点在滑坡体初始时的位置和滑坡体堆积完成后的位置Fig．10 Comparison between original and final positions for twelve monitored disks

4 结论

为了更深入地研究滑坡碎屑流的运动机理和堆积情况，我们根据2008年汶川地震引发的牛眠沟滑坡碎屑流建立了离散元数值模型。本次模型的建立是根据参考资料在实际地质条件的基础上建立的(包括地质结构、地质形态以及岩石力学特性)。由于二维模型的局限性，这种运动过程只能局限在二维空间里。但这并不能妨碍我们通过此例对牛眠沟滑坡碎屑流的碎屑化全过程进行定性地理解。通过研究我们可以发现:

(1)启动及碎屑化阶段发生时间快，且靠近自由面的滑坡体部分碎屑化程度高。

(2)凌空飞行阶段发生在部分滑坡体，且动能大，飞行范围广;伴随着凌空飞行阶段的进行，飞行的滑坡体碎屑化程度加剧。

(3)滑坡碎屑流运动距离和摩擦系数呈负相关，最终推测出实地摩擦系数应在0.2～0.3，和前人学者所得的0.24相符合;并且模拟了出现大规模堵江情况的摩擦系数约为0.1，对指导实际防灾工作有一定的指导意义。

(4)牛眠沟地震滑坡碎屑流的堆积受凌空飞行阶段的影响较大，其堆积情况复杂:受凌空飞行影响的滑坡体部分的堆积情况主要受凌空飞行部分的滑坡体动能影响;而其余部分的滑坡体主要受自重作用下滑，其堆积情况主要受势能及摩擦力的影响。为了得到更精确的结果，我们在后续的工作中应根据具体的地质条件建立三维模型并加以研究。

［1］Keeper D K.Landslides caused by earthquakes［J］.Geological Society of America Bulletin，1984，95(4):406-421.

［2］Keefer D K.Statistical analysis of an earthquakeinduced landslide distribution—the 1989 Loma Prieta，California event［J］.Engineering Geology，2000，58(3):231-249.

［3］李庶林.论我国金属矿山地质灾害与防治对策［J］.中国地质灾害与防治学报，2002，13(4):44-48.LI Shulin.Study on the geological hazard in metal mines and its prevention countermeasures［J］.The Chinese Journal of Geological Hazard and Control，2002，13(4):44-48.

［4］乔建平，吴彩燕.三峡水库区云阳-巫山段斜坡高差因素对滑坡发育的贡献率研究［J］.中国地质灾害与防治学报，2006，16(4):16-19.QIAO Jianping，WU Caiyan.The influence of ralative elevation on landslide occurrence within Yunyang-Wushan section in the Three Gorges Reservoir region［J］.The Chinese Journal of Geological Hazard and Control，2006，16(4):16-19.

［5］刘希林，何思明，乔建平，等.四川泸沽铁矿大顶山矿区泥石流及其防治［J］.中国地质灾害与防治学报，2004，15(3):64-68.LIU Xilin，HE Siming，QIAO Jianping，et al.Debris flow and its control in Dadingshan mining area of Lugu iron mine in Sichuan province［J］.The Chinese Journal of Geological Hazard and Control，2004，15(3):64-68.

［6］李德心，何思明，李军歌，等.降雨作用下山区边坡稳定性分析［J］.中国地质灾害与防治学报，2011，22(1):7-10.LI Dexin，HE Siming，LI Junge，et al.Effect of rainfall on slope stability in mountainous area［J］.The Chinese Journal of Geological Hazard and Control，2011，22(1):7-10.

［7］Crandell D R，Miller C D，Glicken H X，et al.Catastrophic debris avalanche from ancestral Mount Shasta volcano，California［J］.Geology，1984，12(3):143-146.

［8］Cruden D M，Hungr O.The debris of the Frank Slide and theories of rockslide-avalanche mobility［J］.Canadian Journal of Earth Sciences，1986，23(3):425-432.

［9］Nicoletti P G，Sorriso-Valvo M.Geomorphic controls of the shape and mobility of rock avalanches［J］.Geological Society of America Bulletin，1991，103(10):1365-1373.

［10］MudgeM R. Rockfall-avalanche and rockslideavalanche deposits at Sawtooth Ridge，Montana［J］.Geological Society of America Bulletin，1965，76(9):1003-1014.

［11］Crosta G.Regionalization of rainfall thresholds:an aid to landslide hazard evaluation［J］.Environmental Geology，1998，35(2-3):131-145.

［12］Legros F.The mobility of long-runout landslides［J］.Engineering Geology，2002，63(3):301-331.

［13］CrostaG B， ImposimatoS， Roddeman D G.Numerical modelling of large landslides stability and runout［J］. NaturalHazardsand Earth System Science，2003，3(6):523-538.

［14］CUI Peng，CHEN Xiaoqing，ZHU Yingyan，et al.The Wenchuan earthquake(May 12，2008)，Sichuan province， China， and resulting geohazards［J］.Natural Hazards，2011，56(1):19-36.

［15］ZHANG Ming，YIN Yueping，WU Shuren，et al.Dynamics of the Niumiangou Creek rock avalanche trigg-ered by 2008 Ms 8.0 Wenchuan earthquake，Sichuan，China［J］.Landslides，2011，8(3):363-371.

［16］Sato H P，Harp E L.Interpretation of earthquakeinduced landslides triggered by the 12 May 2008，M7.9 Wenchuan earthquake in the Beichuan area，Sichuan Province，China using satellite imagery and Google Earth［J］.Landslides，2009，6(2):153-159.

［17］Hungr O，McDougall S.Two numerical models for landslide dynamic analysis［J］. Computers ＆Geosciences，2009，35(5):978-992.

［18］LI Xinpo，HE Siming，LUO Yu，et al.Simulation of the sliding process of Donghekou landslide triggered by the Wenchuan earthquake using a distinct element method［J］.Environmental Earth Sciences，2012，65(4):1049-1054.

［19］TANG Chaolung，HU Jyrching，LIN Minglang，et al.The Tsaoling landslide triggered by the Chi-Chi earthquake，Taiwan:insights from a discrete element simulation［J］.Engineering Geology，2009，106(1):1-19.

［20］WANG Fawu，CHENG Qiangong，Highland L，et al.Preliminary investigation of some large landslides triggered by the 2008 Wenchuan earthquake，Sichuan Province，China［J］.Landslides，2009，6(1):47-54.

［21］HUANG Runqiu，LI weile.Analysis of the geohazards triggered by the 12 May 2008 Wen-chuan Earthquake，China［J］. Bulletin of Engineering Geology and the Environment，2009，68(3):363-371.

［22］Chen H，Lee C F.Numerical simulation of debris flows［J］.Canadian Geotechnical Journal，2000，37(1):146-160.

［23］Chen H，Lee C F.A dynamic model for rainfallinduced landslides on natural slopes ［J］.Geomorphology，2003，51(4):269-288.

［24］ZHENG Hong，LIU Defu，LI ChunGuang.Slope stability analysis based on elasto‐plastic finite elementmethod［J］. InternationalJournal for Numerical Methods in Engineering，2005，64(14):1871-1888.

［25］Potyondy D O，Cundall P A.A bonded-particle model for rock［J］.International journal of rock mechanics and mining sciences，2004，41(8):1329-1364.

［26］Kim W B，Yang H S.Discrete element analysis on failure behavior of jointed rock slope［J］.Geosystem Engineering，2005，8(2):51-56.

［27］WANG Caigen，Tannant D D，Lilly P A.Numerical analysis of the stability of heavily jointed rock slopes using PFC2D［J］.International Journal of Rock Mechanics and Mining Sciences，2003，40(3):415-424.

［28］Cundall P A，Strack O D L.A discrete numerical model for granular assemblies［J］.Geotechnique，1979，29(1):47-65.

［29］ZHOU Jiawen，SU Zhiman，YANG Xingguo.Kinetic Friction Coefficient and Mass Movement Process of Large Rock Avalanches Triggered by the Wenchuan Earthquake［M］//Earthquake-Induced Landslides.Springer Berlin Heidelberg，2013:521-528.

［30］杨冰.边坡动力破坏过程及塌落范围的颗粒力学模拟［D］.北京:清华大学土木工程系，2011.YANG Bing.Granular simulation on the process of slope failure and collapse under seismic load［D］.Beijing:The Civil Engineering Department of Tsinghua University，2011.

［31］Yuan R M，Tang C L，Hu J C，et al.Ejection mechanism of the Donghekou landslide triggered by the 2008 Wenchuan Earthquake revealed by discrete element modeling［J］.Natural Hazards and Earth System Sciences Discussions，2013，1(6):7667-7700.

［32］WU Shuren，WANG Tao，SHI Ling，et al.Study on catastrophic landslides triggered by 2008 great Wenchuan earthquake，Sichuan，China［J］.J Eng Geol，2010，18(2):145-159.

［33］许强，裴向军，黄润秋.汶川地震大型滑坡研究［M］.北京:科学出版社，2009.XU Qiang，PEI Xiangjun，HUANG Runqiu.Largescale landslides induced by the Wenchuan Earthquake［M］.Beijing:Science Press，2009.