高位滑坡动力学特征分析

徐杨青，乔龙腾，杨龙伟，江强强，王孝臣

(1.煤炭科学研究总院，北京 100013；2.中煤科工集团武汉设计研究院有限公司，湖北 武汉 430064)

近年来，受到全球极端气候变化和人类工程活动等影响，高位滑坡地质灾害呈现频发状态，给人民生命财产造成巨大威胁[1]，如2017年6月24日发生的四川茂县新磨村滑坡造成83人死亡[2]，2019年7月23日发生的贵州省水城县鸡场镇坪地村滑坡造成43人死亡[3]，2021年2月发生的印度北阿坎德邦查莫利高位冰川崩裂灾害造成超过200多人失踪[4]。由于高位滑坡地质具有急剧启动，破坏性大和运动范围广等特点，其防灾减灾难度较大，一直是工程地质灾害领域研究的热点和痛点[5]。

对于高位滑坡灾害定义目前尚未统一，大量研究成果表明，其主要是指滑坡剪出口与滑坡前缘的垂直高差大于50m，且滑源区多分布在高山峡谷等地区，植被茂密，这导致其滑坡演化过程十分隐蔽[6]。目前围绕高位滑坡灾害的研究主要聚焦于成灾机制、滑坡动力学和远程运动堆积等方面，其中高位滑坡灾害动力学是分析高位滑坡动力灾害效应的重中之重，这是因为滑坡动力学计算分析得到的滑体速度、堆积体厚度和最远运动距离等都可为滑坡防灾减灾等提供重要参考[7，8]。广大专家学者围绕滑坡动力灾害学开展了卓有成效的工作，殷跃平等基于滑块-弹簧模型的基础上提出了高位远程滑坡动力侵蚀犁切模型，并计算了滑体运动过程中受到的犁切阻力[9]。胡伟等通过现场调查和室内试验发现鸡尾山高位滑坡的底部摩擦温度高达790℃，提出了高位远程滑坡滑带高温热分解机理[10]。高杨等基于Hertz接触模型和等效流体理论计算分析了动力冲击铲刮对高位滑坡体积效应等影响，并提出了四种动力冲击铲刮模式[11]。严炎和崔一飞等通过对高位滑坡灾害震动信号分析，反演了高位滑坡运动过程[12]。罗刚和程谦恭等总结了高位高能岩崩灾害的动力冲击力及工程防护措施等[13]。邢爱国基于离散元模型方法计算了贵州纳雍高位崩塌运动过程及其震动信号特征[14]。Zhang Ming探讨了重庆鸡尾山滑坡的岩体强度降低可能与矿物质分解有关，进而诱发了滑坡急剧启动[15]。Watkins通过遥感技术分析了火星上的远程滑坡运动特征，发现水化硅酸盐润滑基底滑动带这一机制起到重要作用[16]。上述研究为高位滑坡动力学灾害研究提供了重要指导方向，也为开展类似研究提供有力支撑。

基于此，本文探讨了高位滑坡灾害运动学特征，并以菲律宾 Guinsaugon 滑坡为例，通过基底流变模型动力学方法计算分析滑坡运动速度、堆积体厚度和最远运动距离等，并与雪橇模型计算结果进行对比，最后总结该滑坡的动力学过程，为类似的高位滑坡地质灾害研究提供借鉴。

1 菲律宾 Guinsaugon 滑坡概况

2006年2月17日，菲律宾莱特岛Saint Bemard镇Guinsaugon村发生高位滑坡地质灾害，运动过程中掩埋一所小学，此次灾害共造成1800人死亡，震惊世界[17-19]。事后调查发现，该研究区地处Philippine 断层，其活动速率达到2.5cm/a，且该地区地震活跃，这些都加剧了岩体破碎程度。该地区降雨丰沛，月均降雨量在600～800mm，滑坡发生前该地区出现连续16d降雨，累计达到700mm。由此可见，地质营力活跃、岩体结构破碎和降雨是该滑坡发生的主要诱发因素。滑源区体积约为4×106m3，滑坡后缘高程为800m，剪出口高程约为260m，前缘高程为30m，最远距离为4100m，呈现出现典型的高位远程滑坡灾害特征，详情如图1所示[17-19]。根据Guinsaugon滑坡的运动堆积和空间分布特征，将其划分为滑源区、碎屑流区和堆积区。

滑源区：滑源区海拔高程范围为380～800m。母岩岩性主要为角砾岩、碎屑岩和沉积岩等为主。岩体发育多组反倾节理，这些与岩层面一起将岩体切割成不对称“楔形体”，高度约为420m。在滑坡发生前，岩体后缘裂缝已贯通至基岩处，为滑坡发生提供了优势入渗通道。

碎屑流区：碎屑流区位于海拔50～380m，高差约为330m，斜长约为2027m。根据现场目击者描述，滑体撞击坡面并散落堆积，顺着原冲沟进行运动，呈现管道型碎屑流。

堆积区：堆积区坡度较缓，其水平长度约为1000m，其前缘已接近河流，平面形状呈现“扇形”。由于其运动路径上有稻田覆盖，降低了滑体运动过程中受到的摩擦阻力，加速了滑体的远程运动。

图1 Guinsaugon滑坡

2 高位滑坡运动学特征分析

高位滑坡地质灾害由于其滑源区多位于脊状山顶，在降雨和地震等诱发因素的作用下，岩体结构进一步破碎，加速滑坡失稳。高位滑坡急剧启动后，由于滑坡前后缘间的垂直落差大，经过重力势能转化成动能，导致滑体的运动速度加快，平均速度大于20m/s，也使得滑体运动的最远距离一般超过1km，形成典型的高位远程滑坡灾害，极易造成群死群伤事件。

在运动过程中，根据滑坡运动过程的空间演化特征及速度变化趋势，将滑坡运动路径主要分为滑源区，碎屑流区和堆积区。滑源区由于植被覆盖茂密，其岩体裂化演化过程十分隐蔽，这使得滑坡发生具有突发性。此外，滑体在剪切启动过程中释放的能量也加剧了滑坡运动动能，使得滑体在启动后，短时间内速度增长至2～3m/s。碎屑流区内地表多覆盖松散堆积物，这使得滑体在运动至该区域时会铲刮并裹挟地表松散堆积物，实现了滑坡“体积效应”。由于碎屑流区内坡度逐渐由陡变缓，且碎屑流区内地表微地貌起伏变化，这导致滑体在该区域内的速度波动起伏，平均速度在20m/s以上。在堆积区内，由于滑体在运动过程中受到河流、对岸山体等阻碍，其速度逐渐变缓直至停止。从碎屑流区到堆积区的过程中，由于地形逐渐开阔，滑体从碎屑流运动状态渐变至散落堆积，由于滑体之间的相互挤压碰撞，使得滑体速度呈现差异化，进而在散落堆积过程中形成定向排列、横向脊和纵向脊等堆积地貌特征。

综上所述，高位滑坡孕灾过程隐蔽，急剧发生、动力冲击和远程成灾等特点。加强高位滑坡动力学过程分析，对于高位滑坡风险防控显得尤为重要。

3 高位滑坡的动力学模型

3.1 模型原理

高位滑坡在快速运动过程中，与地表产生巨大碰撞并粉碎成块状，进行形成滑坡碎屑流，铲刮地表松散堆积物，使得滑体体积增大。本研究主要利用等效流体计算方法，主要是基于拉格朗日解析解来计算组成滑体的一系列块体在运动过程中受到的基底剪切阻力(图2和图3)[20]。在刻画高位滑坡动力学过程中，Friction模型(F)和Voellmy模型(V)被广泛运用，其中Friction模型可以较好地阐释块体受到的基底阻力，Voellmy模型适用于湿型滑坡泥石流的流变特性。本文的研究重点是高位滑坡碎屑流的动力学特征，且受到连续降雨影响，Guinsaugon滑坡呈现出一定的滑坡泥石流的状态，故采用Friction模型和Voellmy模型来进行分析(表1)，其主要是在曲线坐标系中滑体划分为多个块体，对各个块体建立物理平衡方程进行求解，其中各区域块体受到的摩擦阻力根据其Friction模型和Voellmy模型进行选取，进而求解滑体的运动距离和运动速度，具体公式见式(1)—(5)[21]，上述公式已被很多专家学者编写成相应的应用程序，为高位远程滑坡运动计算提供了支撑。

F=ΥHiBidSsinα+P-T

(1)

Friction模型：T=σ(1-γu)tanφ

(4)

式中，F是高位滑坡下滑力；Υ是滑体重度；H是滑块的高度；B是滑块的宽度；P是相邻块体间的切向压力；T是滑体运动过程中受到的阻力；v是滑体速度；M是滑块的动量通量；S是曲线位移；i是滑块的编号；j是滑块间的边界次序；σ是滑体运动过程中块体受到的竖向总应力；γu是孔隙水压力系数；φ是滑块的动摩擦角。ξ是湍流项系数。

图2 菲律宾Guinsaugon滑坡数值模型

图3 块体之间的受力分析

表1 运动路径模型选择

3.2 模型建立及参数选取

根据高位滑坡动力学特征，建立高位滑坡动力学基底流变计算模型(图2和图3)，结合滑体的运动轨迹，将滑体的运动路径分为滑源区、碎屑流区和堆积区等。由于很难获取滑体的岩土体强度参数，根据已有的大量岩质滑坡计算案例和试错法可知[22-24]，对于高位滑坡远程动力学基底流变剪切计算，多采用Friction-Voellmy-Friction计算模型组合。由式(4)可知，Frictional模型主要受到动摩擦角φ和孔隙水压力γu的影响，根据已有大量的岩质滑坡案例的数值试验，确定了模拟参数见表2。

表2 Guinsaugon滑坡数值模拟建议参数

4 计算结果分析

4.1 滑坡前后缘运动特征

计算结果显示(图4和图5)，Guinsaugon滑坡运动总时长约为90s，滑体在水平距离X=1896m处时的运动速度为最大，约为52.95m/s，这是由于该处地处坡脚，滑体经过重力势能转化成动能，到达该处时速度最大。

图4 滑坡前后缘随时间变化

图5 滑体前后缘随滑程变化

由图4和图5可见，滑体前缘和后缘的速度显示滑体运动主要经历三个阶段，以滑坡前缘为例：①滑坡急剧加速阶段，滑体启动后，在水平X=960～1000m范围内，滑体速度急剧增加，从0增长到33m/s，这主要是由于滑体急剧启动，能量处于急剧释放的过程；②滑体稳定加速阶段，在X=1000～3000m范围内，滑体的速度保持稳定，基本维持在45～53m/s，这主要是由于受地形限制，滑体呈现管道状碎屑流状态，一直保持高速运动状态；③滑体减速阶段，在X=3000～4100m范围内，滑体速度逐渐降低，直至停止，这是由于滑体运动至平原地区，坡度急剧变缓，滑体动能逐渐消耗。

4.2 堆积体形态分析

滑体的堆积体平均厚度约为20～40m，如图6所示。其中滑源区只有少量残余碎屑流堆积，其厚度不足20m；碎屑流区域由于坡度相对较缓，运动路径上停积了一定厚度的滑体，其厚度约为30m；在堆积区内，堆积体的形态呈现“动车”头型，平均厚度约为15m。此外，滑坡前缘运动至4100m处停止，与实际情况相近。

图6 堆积体形态

4.3 典型特征点分析

根据对于滑坡灾害分析和已有的现场调查资料分析，本文选取了水平距离X=1100m为典型特征点，因为该点地处滑坡剪出口附近，能显示滑坡急剧启动后速度变化过程和滑源区堆积体的残余变化情况。由图7和图8可知，该点的速度最大为45.5m/s，最终的堆积体厚度约为27m。

图7 典型点速度随时间变化

图8 典型点堆积体厚度随时间变化

4.4 模型对比和成灾模式

4.4.1 模型对比

高位滑坡动力学方法研究的主要手段包含数值模拟法和统计学方法。Scheidegger统计了大量的高位滑坡案例，提出了“等效视摩擦角”概念(图9)[25]，即滑坡前缘最远点与滑坡后缘点之间的连线的斜率，其本质上就是滑体势能转化成动能(该方法也称为“雪橇模型”)，见式(6)。

式中，ν0为滑体的速度；g为重力加速度；H为滑体运动点与滑坡后缘之间的垂直高差；f为等效视摩擦系数；L为滑坡计算点与滑坡后缘点之间的水平距离。

图9 等效视摩擦系数概念示意

通过图9和式(6)，可以计算Guinsaugon滑坡的滑程、垂直高差与等效视摩擦系数之间的关系。由此，计算Guinsaugon滑坡在滑坡剪出口(水平距离约为1000m)处的运动速度：

根据目击人描述，Guinsaugon滑坡的平均速度约为27～38m/s[26]，其最大速度预估为40.5～57m/s。目击证人的预估速度、雪橇模型和基底流变模型三者进行比较，雪橇模型计算结果也显示滑坡运动过程中速度急剧增加，波动性高速运动和减速三个阶段。此外，雪橇模型计算结果偏大，这是由于雪橇模型在计算过程中，忽略了滑体运动过程中受到的基底剪切阻力，且对滑动路径微地貌进行了简化。传统的雪橇模型适用于应急抢险过程中的应急估算，而基于基底剪切阻力的计算模型可以精细化刻画滑坡动力学运动过程，与实际速度相近。

表3 Guinsaugon滑坡的雪橇模型计算

4.4.2 成灾模式探讨

Guinsaugon滑坡的母岩以角砾岩、碎屑岩等为主，岩体相对破碎，滑源区坡度较陡，约为48°，滑体呈现楔形体，其高度约为450m，并发育有多组节理，使得岩体呈现碎裂结构。在滑坡发生前16d连续降雨的影响下，岩土体强度进一步降低，使得楔形体失稳滑动，撞击坡面并呈现碎块状，以碎屑流的形式继续向坡角运动。随着坡度变缓，滑体向四周散落，最终直至停止。基于此，建议将“楔形体”块体作为野外高位滑坡识别的重要标志之一，这是因为楔形体失稳后，滑体撞击地面呈现出巨大的动力冲击力，并扰动(动力侵蚀)滑动路径上松散的堆积物，进而加剧滑体体积效应。这种楔形体岩质滑坡多发育在我国西南山区线型工程两侧，严重威胁交通安全运行。对于楔形体岩质滑坡的风险防控，多采用清方+支护措施，并加强监测。

与Guinsaugon滑坡类似，我国西南山区的高位远程滑坡碎屑流灾害频发，诸如2009年发生的重庆鸡尾山滑坡、2017年发生的贵州纳雍滑坡、2022年9月5日四川泸定发生6.8级地震诱发的大量高位远程滑坡碎屑流灾害。这种灾害大多数在降雨、地震或者人类工程活动诱发因素的影响，滑源区岩体强度降低，进而滑坡急剧启动、动力撞击坡面后破碎解体，滑坡-碎屑流在向前运动过程中动力侵蚀地表松散堆积体，最后遇到沟谷和河流等阻挡后形成滑坡坝，坝体在堰塞湖的作用下溃决形成洪水灾害[27]。对于这种高位远程滑坡链式灾害应加强地质灾害风险调查评价，进而减少灾害损失[28]。

本次高位滑坡灾害动力学分析主要以计算滑体运动速度和堆积体厚度等为主，建模过程主要采用条块法，以等效流体的方式来进行计算，缺乏在建模过程中考虑岩体节理、且岩土体参数主要以试错法来取得，这些与真实情况有一定差距，这也是今后研究的重要方向。

5 结 论

1)Guinsaugon滑坡属于典型的高位远程滑坡灾害，将其运动路径划分为滑源区、碎屑流区和堆积区。滑源区岩体形状为楔形体，碎屑流区呈现管道状和堆积区呈扇形。

2)Guinsaugon滑坡的主要诱发因素为地质构造运动活跃、岩体结构破碎和降雨，其运动阶段主要划分为急剧启动-撞击碎屑化-远程堆积的模式，其速度运动阶段主要为急剧加速-波动性维持-减缓停积的过程。

3)模型计算结果显示，滑坡运动时长约为90s，在水平距离X=1896m处时的运动速度为最大，约为52.95m/s，堆积体的平均厚度约为15～20m，最远运动距离约为4100m，与实际情况相近。与传统的雪橇模型计算结果相比，基于基底剪切阻力计算模型可以更好的精细化刻画滑体动力学过程。