基于MatDEM的烟家沟滑坡演化过程数值模拟分析

栗晓松，范 文，2，曹琰波，全倬梁

1.长安大学 地质工程与测绘学院，陕西 西安710054；2.工业和信息化部电子综合勘察设计院，陕西 西安710054

0 引言

2015年8月12日，陕西省山阳县烟家沟村发生一起突发性滑坡灾害，坡体的总体积达168×104m3，经过烟家沟支沟左、右两侧山体的先后阻挡，滑体的运动方向先后两次发生偏转，最终沿沟谷下游运动，给当地居民的生命财产造成了严重损失［1］.

正确认识滑坡的运动、堆积特征是对滑坡灾害进行合理的评价，提出合理的治理措施的前提.滑坡体滑动速度、滑动距离和堆积特征是对滑坡灾害进行有效评估的重要参考指标.滑坡的发生一般较突然且危害程度大［2-3］，许多学者针对滑坡运动速度、堆积特征等方面开展了广泛的研究.鲁晓兵等从理论分析出发，主要研究滑坡坡角、土体的内摩擦角对斜坡失稳下滑的运动特征和最终堆积分布的影响［4］；樊晓一等从滑槽模型试验入手，对滑坡体碎屑流前缘运动速度研究，重点考虑坡度、颗粒级配及坡角等因素［5］；郝明辉等在现场试验的基础上，结合调查结果，具体地研究了滑坡体在下滑过程中颗粒分选性、反序列堆积的问题［6］；Hungr等人采用数值模拟的方法，基于DAN模拟了滑坡发生过程，得出沟谷形态是影响滑坡运动特征和堆积分布的重要影响因素［7］.目前，虽大多数学者采用多种研究方法，如现场调查研究、理论分析计算、模型试验和数值模拟等对滑坡运动过程进行研究，但是总体来说针对滑坡运动发生过程的直观性和整体性呈现还相对较弱.

基于此，本文以山阳县烟家沟滑坡作为研究对象，利用高性能离散元软件MatDEM［8-9］，依据无人机现场实拍，获得烟家沟滑坡高程信息，经ArcGIS初步栅格化处理后，结合滑坡实际特征建立初始三维地质模型，并对模型赋予真实的力学性质参数，模拟真实情况下滑坡启动发生的全过程.

1 滑坡区域地质条件

1.1 地形地貌

滑坡现场遥感照片及斜坡剖面如图1、2.滑坡区位于商洛市山阳县中村镇，为大西沟和烟家沟的交界处.滑坡区总体上属基岩中山陡坡区，地质构造发育，处于倒转褶皱的一侧［10］.滑坡区原始斜坡地形陡峻，坡度约为35°，高程1010～1300 m，高差290 m.前缘大西沟沟谷深切，坡向与地层产状一致，呈突出的山梁，构成陡倾顺层边坡.

1.2 地层岩性

滑坡区岩性复杂，地层出露较多，主要地层为震旦系灯影组、寒武系下统水沟口组、寒武系中统岳家坪组和在沟谷地带大量分布的第四系［11］.详见表1.

图1 山阳烟家沟滑坡遥感图像Fig.1 Remote sensing image of Yanjiagou landslide in Shanyang County

图2 滑前斜坡地质结构剖面图Fig.2 Geological profile of the slope before sliding

1.3 地质构造

滑坡区地处秦岭南部，属薄皮逆冲推覆构造带，是舒家坝－刘岭逆冲推覆构造带和凤县－镇安逆冲推覆构造带的接合部位［12］，断层和褶皱较发育，且整个区域地层层序倒转，使得水沟口组地层处于震旦系灯影组的地层之下，呈平行不整合接触.坡体结构上部坚硬，下部软弱.在区域构造背景下，断层和褶皱呈东西向分布，地应力较为集中，岩石风化程度高，斜坡破坏现象经常发生.

1.4 滑坡运动过程特征

据野外调查发现，滑源区岩土体在重力作用下失稳破坏以66°方向开始滑动，在受到烟家沟左侧山体碰撞后朝105°方向继续滑动，进而受到烟家沟右侧山体阻挡后转向60°方向运动直至停止（图3）.滑坡运动过程呈现典型的右旋特征.

表1 滑坡区地层岩性Table 1 Stratigraphic lithology of landslide area

图3 滑坡运动过程分析Fig.3 Movement process of landslide

2 数值模型建立

2.1 颗粒流离散元程序

离散元法（DEM）是由Cundall等［13］首次提出，主要是用以研究颗粒状物质的运动及相互作用.随着近年来计算机技术进步和离散元理论的发展，离散元法在多种领域得到广泛应用［14］，包括固体力学领域中的颗粒堆积体力学性质研究［15］，地质领域中各类构造的形成和演化研究等.近年来很多研究人员［16-17］采用离散元数值模拟方法对滑坡运动过程进行反演分析，取得了较好的成果.因此，本文选择矩阵离散元法进行烟家沟滑坡运动数值模拟.

MatDEM采用矩阵离散元计算法和三维接触算法，在MATLAB语言运行环境下自主建立模型并完成离散元生态下的能量守恒计算.该算法运行速度更快，模拟效果更直观，对模拟滑坡、岩爆、撞击破坏、桩土作用等具有突出的优点.MatDEM软件可将一定数量的具有特定力学性质的单元按照相应组合方式胶结在一起，建构出多种符合实际结构性的三维岩土体（图4）.本文采用线弹性接触方式建立滑坡模型，单元间胶结力可用法向弹簧力和切向弹簧力来表示（图5）.弹簧接触和切向弹簧接触.

图4 颗粒胶结三维球体模型Fig.4 3D sphere model of particle cementing

通常情况下，颗粒单元之间的法向力和法向变形、剪切力和剪切变形都维持在稳定的限度内.满足：

图5 线弹性接触方式示意图（据文献［18］）Fig.5 Schematic diagram of linear elastic contact mode（From Reference［18］）

其中，Fn为法向力，Kn为法向刚度，Xn为法向变形；Fs为剪切力，Ks为切向刚度，Xs为剪切变形.颗粒单元间的最大法向力满足：

Xd为断裂变形.当单元在法向上外力逐渐增加至超过Fnmax时，单元间胶结断开，拉力消失.单元胶结所能承受最大切向力满足：

其中，Fs0为初始抗剪力，μp为摩擦系数.在实际模拟过程中，切向力大于Fsmax时，单元胶结亦发生断裂，单元产生错动，此时单元间只剩滑动摩擦－μpFn.

MatDEM数值模拟计算符合力－位移定律和牛顿第二定律，模型时间步（dT）很小，此间颗粒的加速度及速度保持不变，通过时间步迭代来计算颗粒的受力、加速度、速度和位移，通过反复迭代最终实现离散元法动态模拟［18］.

2.2 滑坡模型建立

建立准确的三维模型是模拟真实滑坡的基础，通过无人机航拍获取高精度高程坐标，经ArcGIS软件栅格化后，导入MatDEM软件进行后续处理分析，最终建立起山阳烟家沟滑坡的三维地形模型（图6）.通过滑前、滑后数字高程模型对比并结合现场实际勘察发现，滑坡体堆积分布分区较明显，无明显刮铲区，且岩体单元的半径小，坡体前缘破碎严重，故依据坡体形态和物质组成特点把滑坡区域概括为滑源区和堆积区.

3 模拟参数选取

颗粒单元的力学性质通常会在其重新堆积和胶结后发生较大改变.在本研究中，单元的力学参数主要通过宏微观参数转换公式计算以及数值模拟测试来确定.通过采用这样的参数建立的模型，能够使得岩土体在最大程度上符合天然状态下的结构性和复杂性.

图6 滑坡前、后数字高程模型及分区Fig.6 Digital elevation model and zoning before and after landsliding

宏微观参数转换公式计算：线弹性接触模型的5个微观力学参数，分别是摩擦系数（μp）、初始抗剪力（Fs0）、正向劲度系数（Kn）、切向劲度系数（Ks）和断裂位移（Xb），这些力学参数经由所选材料的5个宏观力学性质（泊松比v、杨氏模量E、抗压强度Cu、抗拉强度Tu和内摩擦系数μi）结合宏微观参数转换公式计算得到.转换公式如下所示［19］：

其中，d为颗粒直径，I为比例系数，

获得参数具体步骤：①经过室内物理实验测得岩体力学参数值（表2），结合转换公式求取初始数值模型的物理力学参数，将参数赋值给颗粒单元并建立滑坡堆积模型；②运行MatDEM内置的抗拉强度、单轴压缩和抗压强度测试程序，求取模型测试的杨氏模量及抗拉强度等性质；③将数值测试的模型力学性质与实际测量的岩体力学性质对比得到一个尺度值；④最后把实际测量的岩体力学参数与尺度值对比，再次放到上述转换公式中计算，把此次的单元力学参数赋值给颗粒单元.

表2 实测滑坡体岩石力学性质Table 2 Mechanical properties of measured landslide dolomite

经过至少3次②③④过程可获取误差小于2%的单元力学参数.多次迭代计算后得出用于滑坡模型单元的宏微观力学参数值（如表3）.

表3 宏、微观力学参数值Table 3 Macro and micro mechanical parameters

根据前人的试验研究及其相关数值模拟结果［20］得知，岩土体的抗压强度、抗拉强度及弹性模量值有“尺度效应”的现象.在试验室当中实测的小尺寸岩体的物理力学性质通常情况下要比现场的大尺度岩体的物理力学性质偏大.本文数值模拟计算过程实际使用的强度参数值约等于室内实测试验值的40%.

4 模拟结果及分析

4.1 滑坡滑移距离变化特征

结合表3中所得模拟参数进行迭代计算，对烟家沟滑坡运动情况进行反演分析，模拟过程如图7所示.图7a、b、c、d分别截取了第6、10、16、28 s时的位移距离.经读取模型数据得，滑坡堆积体长度约为510 m，最大滑移距离约为570 m，与现场调查得到的最远滑动距离600 m较接近［11］，进一步验证了选取参数的合理性，有效地再现了烟家沟滑坡滑动和堆积过程.

4.2 滑坡速度变化特征

通过提取滑坡体活动单元的速度、滑坡表面跳跃单元的速度及对应时间，生成滑坡岩土体速度－时间变化曲线（图8）.滑坡体内部在应力条件下发生失稳破坏，在5～12 s时间内完成了启动、加速过程，最大速度达到44 m/s；滑坡启动约9 s时，滑体平均速度达到最大值，约23.5 m/s，滑坡启动28 s后滑体平均速度几乎收敛于0，滑体运动趋于停止，仅有零星的表面单元向下跳跃运动.

根据滑坡运动速度、滑移距离及堆积位置，将滑坡堆积过程分为加速碰撞堆积阶段、整体滑动堆积阶段和减速堆积阶段.

图7 滑坡运动过程模拟Fig.7 Simulation of landsliding process

图8 滑坡模型单元速度－时间变化曲线Fig.8 Speed vs.time curves of landslide model unit

加速碰撞堆积阶段（图7b）：滑源区前缘颗粒单元碰撞解体滑落，在内部应力释放及重力作用下沿着坡向运动，到达坡脚下“V”型谷时平均速度达到约23 m/s.由于受到沟谷左侧山体的阻挡，大部分滑体单元在此处堆积；其余颗粒单元在惯性力影响下以较高速度向前冲出，并顺着沟谷方向继续运动.

整体滑动堆积阶段（图7c）：剩余滑体单元整体以较高速度下滑并到达第二次碰撞点，在此阶段沟谷宽度变大，部分滑坡单元能量得以充分释放并沿途堆积，其余部分滑体单元会沿沟谷发生偏转继续运动.

减速堆积阶段（图7d）：在16～28 s时滑体单元速度逐渐降低，直至到达坡底，主体仅发生局部滑动和调整，堆积基本完成.

4.3 滑坡堆积厚度分布特征

根据滑坡前后高程数据对比，获取滑坡堆积厚度分布（图9）.烟家沟滑坡发生后，滑源区单元碰撞解体滑落，故滑源区堆积厚度为负值，从图中可以看出滑源区最大深度约为43 m.随着滑体下滑，在运动方向会受到沟谷左侧山体的阻挡，颗粒之间相互挤压碰撞，导致大部分的颗粒单元能量损耗［21］.颗粒大量在沟谷处堆积，使得平均堆积厚度大幅增加，最大堆积厚度为48 m.另外部分单元在惯性作用力下保持原有运动趋势向下游方向继续运动，此阶段沟谷宽度变大，颗粒堆积厚度减小，且有小范围“爬高现象”.剩余颗粒单元会沿下游方向继续运动并不断沿途堆积，堆积厚度逐渐减小.堆积区堆积厚度呈现出先增加后减少的总体变化趋势.

图9 沿滑动方向的堆积厚度分布Fig.9 Distribution of accumulation thickness along the sliding direction

5 讨论

本文基于MatDEM将巨量的颗粒单元用于数值建模，借助无人机航拍等手段获得准确高程信息，于三维模型中呈现出复杂的滑坡要素特征，改进了传统二维模型在滑坡运动发生过程中直观性和整体性呈现不足的缺陷，很好地再现了烟家沟滑坡的运动、堆积过程.尽管精确的滑坡高程信息和MatDEM软件高效的矩阵离散元计算法在滑坡演化过程模拟中具有无可比拟的优势，但在地质构造、坡度、颗粒级配等因素如何具体影响滑坡的启动发生问题上还有待进一步研究.实际的滑坡运动及堆积过程非常复杂，影响因素众多，后续还需要进一步结合理论推导、模型实验等手段来进一步分析.本文旨在通过以上的数值模拟研究，为以后类似滑坡的运动和堆积特征分析、成灾范围划定及灾害防治提供一定的技术参考.

6 结论

通过上述研究取得了以下结论：

1）山阳县烟家沟滑坡从启动到最终静止整个过程持续时间约28 s，滑坡体单元最大速度可达44 m/s，平均速度最高可达23.5 m/s，故将烟家沟滑坡归为高速滑坡；滑坡堆积体长度约为510 m，最大滑移距离约为570 m，最大堆积厚度约为48 m.该类滑坡对人员及建筑物具有较强的破坏性.

2）根据滑坡运动速度和距离将滑坡堆积过程分为加速碰撞堆积阶段、整体滑动堆积阶段和减速堆积阶段.顺着滑坡运动路径的方向，堆积厚度呈现先增加后减小的总体变化趋势，堆积厚度在堆积体的后部最大.

3）经过MatDEM数值模拟后所呈现出的堆积特征和现场调查结果完全符合，证明矩阵离散元法模拟滑坡演化过程结果具有可靠性.

致谢：本文使用的高性能离散元软件MatDEM由南京大学刘春团队自主研发，模拟过程中使用的代码为在其源代码基础上进行的改编.