2021, Vol. 30
PDF
2. 工业和信息化部电子综合勘察设计院, 陕西 西安 710054
2. Electronic Research Institute of Engineering Investigation and Design, MIIT, Xi'an 710054, China
2015年8月12日，陕西省山阳县烟家沟村发生一起突发性滑坡灾害，坡体的总体积达168×104 m3，经过烟家沟支沟左、右两侧山体的先后阻挡，滑体的运动方向先后两次发生偏转，最终沿沟谷下游运动，给当地居民的生命财产造成了严重损失[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 山阳烟家沟滑坡遥感图像 Fig.1 Remote sensing image of Yanjiagou landslide in Shanyang County |
|
图 2 滑前斜坡地质结构剖面图 Fig.2 Geological profile of the slope before sliding |
滑坡区岩性复杂，地层出露较多，主要地层为震旦系灯影组、寒武系下统水沟口组、寒武系中统岳家坪组和在沟谷地带大量分布的第四系[11]. 详见表 1.
|
表 1 滑坡区地层岩性 Table 1 Stratigraphic lithology of landslide area |
2y
1sg3
1sg2
1sg1滑坡区地处秦岭南部，属薄皮逆冲推覆构造带，是舒家坝-刘岭逆冲推覆构造带和凤县-镇安逆冲推覆构造带的接合部位[12]，断层和褶皱较发育，且整个区域地层层序倒转，使得水沟口组地层处于震旦系灯影组的地层之下，呈平行不整合接触. 坡体结构上部坚硬，下部软弱. 在区域构造背景下，断层和褶皱呈东西向分布，地应力较为集中，岩石风化程度高，斜坡破坏现象经常发生.
1.4 滑坡运动过程特征据野外调查发现，滑源区岩土体在重力作用下失稳破坏以66°方向开始滑动，在受到烟家沟左侧山体碰撞后朝105°方向继续滑动，进而受到烟家沟右侧山体阻挡后转向60°方向运动直至停止（图 3）. 滑坡运动过程呈现典型的右旋特征.
|
图 3 滑坡运动过程分析 Fig.3 Movement process of landslide |
离散元法（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]) a-法向弹簧接触（normal spring contact）；b-切向弹簧接触（tangential spring contact） |
通常情况下，颗粒单元之间的法向力和法向变形、剪切力和剪切变形都维持在稳定的限度内. 满足：
| $ {F_{\rm{n}}} = {K_{\rm{n}}}{X_{\rm{n}}} $ | (1) |
| $ {F_{\rm{s}}} = {K_{\rm{s}}}{X_{\rm{s}}} $ | (2) |
其中，Fn为法向力，Kn为法向刚度，Xn为法向变形；Fs为剪切力，Ks为切向刚度，Xs为剪切变形. 颗粒单元间的最大法向力满足：
| $ {F_{{\rm{nmax}}}} = {K_{\rm{n}}}{X_{\rm{d}}} $ | (3) |
Xd为断裂变形. 当单元在法向上外力逐渐增加至超过Fnmax时，单元间胶结断开，拉力消失. 单元胶结所能承受最大切向力满足：
| $ {F_{{\rm{smax}}}} = {F_{{\rm{s0}}}} - {\mu _{\rm{p}}}{F_{\rm{n}}} $ | (4) |
其中，Fs 0为初始抗剪力，μp为摩擦系数. 在实际模拟过程中，切向力大于Fsmax时，单元胶结亦发生断裂，单元产生错动，此时单元间只剩滑动摩擦-μp Fn.
MatDEM数值模拟计算符合力-位移定律和牛顿第二定律，模型时间步（dT）很小，此间颗粒的加速度及速度保持不变，通过时间步迭代来计算颗粒的受力、加速度、速度和位移，通过反复迭代最终实现离散元法动态模拟[18].
2.2 滑坡模型建立建立准确的三维模型是模拟真实滑坡的基础，通过无人机航拍获取高精度高程坐标，经ArcGIS软件栅格化后，导入MatDEM软件进行后续处理分析，最终建立起山阳烟家沟滑坡的三维地形模型（图 6）. 通过滑前、滑后数字高程模型对比并结合现场实际勘察发现，滑坡体堆积分布分区较明显，无明显刮铲区，且岩体单元的半径小，坡体前缘破碎严重，故依据坡体形态和物质组成特点把滑坡区域概括为滑源区和堆积区.
|
图 6 滑坡前、后数字高程模型及分区 Fig.6 Digital elevation model and zoning before and after landsliding |
颗粒单元的力学性质通常会在其重新堆积和胶结后发生较大改变. 在本研究中，单元的力学参数主要通过宏微观参数转换公式计算以及数值模拟测试来确定. 通过采用这样的参数建立的模型，能够使得岩土体在最大程度上符合天然状态下的结构性和复杂性.
宏微观参数转换公式计算：线弹性接触模型的5个微观力学参数，分别是摩擦系数（μp）、初始抗剪力（Fs 0）、正向劲度系数（Kn）、切向劲度系数（Ks）和断裂位移（Xb），这些力学参数经由所选材料的5个宏观力学性质（泊松比v、杨氏模量E、抗压强度Cu、抗拉强度Tu和内摩擦系数μi）结合宏微观参数转换公式计算得到. 转换公式如下所示[19]：
| $ {K_{\rm{n}}} = \frac{{\sqrt 2 Ed}}{{4\left( {1 - 2\upsilon } \right)}} $ | (5) |
| $ {K_{\rm{s}}} = \frac{{\sqrt 2 \left( {1 - 5\upsilon } \right)Ed}}{{4\left( {1 + \upsilon } \right)\left( {1 - 2\upsilon } \right)}} $ | (6) |
| $ {X_{\rm{b}}} = \frac{{3{K_{\rm{n}}} + {K_{\rm{s}}}}}{{6\sqrt 2 \left( {{K_{\rm{n}}} + {K_{\rm{s}}}} \right)}} \cdot {T_{\rm{u}}}{d^2} $ | (7) |
| $ {F_{{\rm{s0}}}} = \frac{{1 - \sqrt 2 {\mu _{\rm{p}}}}}{6} \cdot {C_{\rm{u}}}{d^2} $ | (8) |
| $ {\mu _{\rm{p}}} = \frac{{ - 2\sqrt 2 + \sqrt 2 I}}{{2 + 2I}} $ | (9) |
其中，d为颗粒直径，I为比例系数，
| $ I = {\left( {\sqrt {1 + \mu _{\rm{i}}^2} + {\mu _{\rm{i}}}} \right)^2} $ | (10) |
获得参数具体步骤：①经过室内物理实验测得岩体力学参数值（表 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]，进一步验证了选取参数的合理性，有效地再现了烟家沟滑坡滑动和堆积过程.
|
图 7 滑坡运动过程模拟 Fig.7 Simulation of landsliding process |
通过提取滑坡体活动单元的速度、滑坡表面跳跃单元的速度及对应时间，生成滑坡岩土体速度-时间变化曲线（图 8）. 滑坡体内部在应力条件下发生失稳破坏，在5~12 s时间内完成了启动、加速过程，最大速度达到44 m/s；滑坡启动约9 s时，滑体平均速度达到最大值，约23.5 m/s，滑坡启动28 s后滑体平均速度几乎收敛于0，滑体运动趋于停止，仅有零星的表面单元向下跳跃运动.
|
图 8 滑坡模型单元速度-时间变化曲线 Fig.8 Speed vs. time curves of landslide model unit 1-滑坡表面单元最大速度（maximum velocity of landslide surface unit）；2-滑坡体单元平均速度（average velocity of landslide body 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 1-本文模拟结果（simulation result of the current study）；2-实测结果（measured result） |
本文基于MatDEM将巨量的颗粒单元用于数值建模，借助无人机航拍等手段获得准确高程信息，于三维模型中呈现出复杂的滑坡要素特征，改进了传统二维模型在滑坡运动发生过程中直观性和整体性呈现不足的缺陷，很好地再现了烟家沟滑坡的运动、堆积过程. 尽管精确的滑坡高程信息和MatDEM软件高效的矩阵离散元计算法在滑坡演化过程模拟中具有无可比拟的优势，但在地质构造、坡度、颗粒级配等因素如何具体影响滑坡的启动发生问题上还有待进一步研究. 实际的滑坡运动及堆积过程非常复杂，影响因素众多，后续还需要进一步结合理论推导、模型实验等手段来进一步分析. 本文旨在通过以上的数值模拟研究，为以后类似滑坡的运动和堆积特征分析、成灾范围划定及灾害防治提供一定的技术参考.
6 结论通过上述研究取得了以下结论：
1）山阳县烟家沟滑坡从启动到最终静止整个过程持续时间约28 s，滑坡体单元最大速度可达44 m/s，平均速度最高可达23.5 m/s，故将烟家沟滑坡归为高速滑坡；滑坡堆积体长度约为510 m，最大滑移距离约为570 m，最大堆积厚度约为48 m. 该类滑坡对人员及建筑物具有较强的破坏性.
2）根据滑坡运动速度和距离将滑坡堆积过程分为加速碰撞堆积阶段、整体滑动堆积阶段和减速堆积阶段. 顺着滑坡运动路径的方向，堆积厚度呈现先增加后减小的总体变化趋势，堆积厚度在堆积体的后部最大.
3）经过MatDEM数值模拟后所呈现出的堆积特征和现场调查结果完全符合，证明矩阵离散元法模拟滑坡演化过程结果具有可靠性.
致谢: 本文使用的高性能离散元软件MatDEM由南京大学刘春团队自主研发，模拟过程中使用的代码为在其源代码基础上进行的改编.
| [1] |
杨海龙, 樊晓一, 张友谊, 等. 山阳烟家沟滑坡成因机制与运动特征研究[J]. 路基工程, 2016(6): 30-35. Yang H L, Fan X Y, Zhang Y Y, et al. Study on formation mechanism and motion characteristics of Yanjiagou landslide in Shanyang[J]. Subgrade Engineering, 2016(6): 30-35. |
| [2] |
周明浪. 台风"苏迪罗"期间温州滑坡灾害特征及预警成果分析[J]. 地质与资源, 2017, 26(3): 303-309. Zhou M L. Characteristics of landslide disaster in Wenzhou during typhoon Soudelor and analysis of early warning results[J]. Geology and Resources, 2017, 26(3): 303-309. |
| [3] |
徐晨栋, 苑康泽, 郭子坤, 等. 清子高速某工程滑坡诱发机制及治理模拟[J]. 地质与资源, 2020, 29(2): 196-201. Xu C D, Yuan K Z, Guo Z K, et al. Inducement mechanism and treatment simulation of an engineering landslide on Qingzi expressway[J]. Geology and Resources, 2020, 29(2): 196-201. DOI:10.3969/j.issn.1671-1947.2020.02.012 |
| [4] |
鲁晓兵, 王义华, 王淑云, 等. 碎屑流沿坡面运动的初步分析[J]. 岩土力学, 2004, 25(S2): 598-600. Lu X B, Wang Y H, Wang S Y, et al. The primary analysis on the castic gain fow[J]. Rock and Soil Mechanics, 2004, 25(S2): 598-600. |
| [5] |
Fan X Y, Tian S J, Zhang Y Y. Mass-front velocity of dry granular flows influenced by the angle of the slope to the runout plane and particle size gradation[J]. Journal of Mountain Science, 2016, 13(2): 234-245. DOI:10.1007/s11629-014-3396-3 |
| [6] |
郝明辉, 许强, 杨兴国, 等. 高速滑坡-碎屑流颗粒反序试验及其成因机制探讨[J]. 岩石力学与工程学报, 2015, 34(3): 472-479. Hao M H, Xu Q, Yang X G, et al. Physical modeling tests on inverse grading of particles in high speed landslide debris[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34(3): 472-479. |
| [7] |
Hungr O, McDougall S. Two numerical models for landslide dynamic analysis[J]. Computers & Geosciences, 2009, 35(5): 978-992. |
| [8] |
Liu C, Pollard D D, Gu K, et al. Mechanism of formation of wiggly compaction bands in porous sandstone: 2. Numerical simulation using discrete element method[J]. Journal of Geophysical Research: Solid Earth, 2015, 120(12): 8153-8168. DOI:10.1002/2015JB012374 |
| [9] |
刘春, 范宣梅, 朱晨光, 等. 三维大规模滑坡离散元建模与模拟研究——以茂县新磨村滑坡为例[J]. 工程地质学报, 2019, 27(6): 1362-1370. Liu C, Fan X M, Zhu C G, et al. Discrete element modeling and simulation of 3-dimen-sional large-scale landslide-taking Xinmocun landslide as an example[J]. Journal of Engineering Geology, 2019, 27(6): 1362-1370. |
| [10] |
杨海龙, 樊晓一, 裴向军. 基于离散元法的偏转型滑坡运动堆积特征分析[J]. 长江科学院院报, 2020, 37(2): 106-111, 118. Yang H L, Fan X Y, Pei X J. DEM-based analysis of movement and accumulation characteristics of turning-type landslide[J]. Journal of Yangtze River Scientific Research Institute, 2020, 37(2): 106-111, 118. |
| [11] |
李凯. 陕西山阳县中村钒矿区滑坡形成机理及早期识别研究[D]. 西安: 长安大学, 2017. Li K. Study on formation mechanism and early identification of landslide in Zhongcun vanadium mine of Shanyang County, Shaanxi Province[D]. Xi'an: Chang'an University, 2017. |
| [12] |
苏艳军. 山阳-商南钒矿带矿区斜坡破坏机理分析及地质灾害危险性评价[D]. 西安: 长安大学, 2019. Su Y J. Geological disaster hazard assessments and analysis of slope failure mechanism on Shanyang-Shangnan vanadium mine belt[D]. Xi'an: Chang'an University, 2017. |
| [13] |
Cundall P A, Strack O D L. A discrete numerical model for granular assemblies[J]. Géotechnique, 1979, 29(1): 47-65. DOI:10.1680/geot.1979.29.1.47 |
| [14] |
谢莉, 李渝生, 曹建军, 等. 澜沧江某水电站右坝肩岩体倾倒变形的数值模拟[J]. 中国地质, 2009, 36(4): 907-914. Xie L, Li Y S, Cao J J, et al. Numerical simulation of toppling rock mass deformation in the right dam abutment of a hydropower station on the Lancang River[J]. Geology in China, 2009, 36(4): 907-914. DOI:10.3969/j.issn.1000-3657.2009.04.019 |
| [15] |
Goldenberg C, Goldhirsch I. Friction enhances elasticity in granular solids[J]. Nature, 2005, 435(7039): 188-191. DOI:10.1038/nature03497 |
| [16] |
李祥龙, 唐辉明, 熊承仁, 等. 岩石碎屑流运移堆积过程数值模拟[J]. 工程地质学报, 2011, 19(2): 168-175. Li X L, Tang H M, Xiong C R, et al. Numerical simulation of flow and deposition process of rock avalanche[J]. Journal of Engineering Geology, 2011, 19(2): 168-175. DOI:10.3969/j.issn.1004-9665.2011.02.004 |
| [17] |
施凤根. 基于PFC3D的文家沟滑坡高速远程运动学特征研究[D]. 北京: 中国地质大学, 2014. Shi F G. The study of rapid and long-runout characteristics of Wenjiagou landslide based on PFC3D[D]. Beijing: China University of Geosciences, 2014. |
| [18] |
刘春. 地质与岩土工程矩阵离散元分析[M]. 北京: 科学出版社, 2019: 5-7. Liu C. Matrix discrete element analysis of geology and geotechnical engineering[M]. Beijing: Science Press, 2019: 5-7. |
| [19] |
Liu C, Xu Q, Shi B, et al. Mechanical properties and energy conversion of 3D close-packed lattice model for brittle rocks[J]. Computers & Geosciences, 2017, 103: 12-20. |
| [20] |
周喻, 王莉, 丁剑锋, 等. 多尺度节理岩体力学特性的颗粒流分析[J]. 岩土力学, 2016, 37(7): 2085-2095, 2127. Zhou Y, Wang L, Ding J F, et al. Particle flow code analysis of multi-scale jointed rock mass based upon equivalent rock mass technique[J]. Rock and Soil Mechanics, 2016, 37(7): 2085-2095, 2127. |
| [21] |
胡晓波, 樊晓一, 唐俊杰. 基于离散元的高速远程滑坡运动堆积特征及能量转化研究-以三溪村滑坡为例[J]. 地质力学学报, 2019, 25(4): 527-535. Hu X B, Fan X Y, Tang J J. Accumulation characteristics and energy conversion of high-speed and long-distance landslide on the basis of Dem: a case study of Sanxicun landslide[J]. Journal of Geomechanics, 2019, 25(4): 527-535. |