四川省汉源县中海村滑坡动力学特征数值分析

徐奕梓， 樊晓一,2， 张友谊， 田述军， 温 翔， 刘浩南， 郑榆枫， 廖洪阳

(1.西南石油大学土木工程与测绘学院，四川 成都 610500； 2.西南石油大学工程安全评估与防护研究院，四川 成都 610500； 3.西南科技大学土木工程与建筑学院，四川 绵阳 621010)

0 引言

四川省汉源县地处横断山脉北段东缘，为川西高原与四川盆地之间的过渡地带，区内地形起伏大，地质构造复杂且地震活动强烈，是滑坡灾害的频发区[1-2]。近年来，发生了多起地质灾害，导致严重的人员伤亡和财产损失[3-5]。

灾难性滑坡常具有隐蔽性、突发性特征，现场监测和实时获取灾害参数存在较大的难度，模型试验和数值模拟成为了理解和分析这类灾害的有效方法。滑坡运动距离和堆积特征是滑坡防灾减灾和工程防护的主要技术参数，前人利用室内实验模型探究了滑坡运动距离与滑坡体积之间的关系以及滑坡运动距离的影响因素[6-7]，总结了体积、垂直高度、基底摩擦角、地形条件和不同偏转角度等因素对堆积特性的影响[8-9]。计算机的发展为滑坡运动过程的模拟带来更多的可能性，DAN-W模型、离散元模型、DAN3D模型和Massflow模型是4种基于不同理论方法和本构模型建立起来的模型，均能对滑坡失稳破坏和运动堆积过程进行演化分析。国内学者分别采用这4种模型对舟曲泄流滑坡、三溪村滑坡、水城滑坡、甘肃黑方台滑坡的运动特征、堆积特征、能量转化特征和冲击铲刮效应展开数值模拟，取得了较好的研究成果[10-13]。

降雨、地震、水库蓄水、人类工程活动、斜坡结构、特殊地层结构、地质构造等是滑坡发生的主要因素[3,13-15]，地形条件、岩土体组成则影响了滑坡的运动过程和致灾区域。2020年8月21日发生的中海村滑坡在地形条件、岩土体组成上不同于该区域的其他滑坡，本文通过现场调查并结合滑坡区域地形地貌、地层岩性和气候水文条件，分析滑坡灾后特征，在此基础上建立Massflow数值模型，反演整个滑坡动态过程并对滑坡动力学参数和运动堆积形态进行深入探讨，旨在揭示影响滑坡启动的关键因素和动力学特征，为该类滑坡运动机理的研究提供参考和借鉴。

1 滑坡区地质环境条件

1.1 地理位置及地形地貌

滑坡位于距汉源县城约5 km的富泉镇中海村，地处大渡河左岸，地貌类型属于侵蚀构造中低山地貌，海拔高程820～1 060 m，斜坡坡向N15°E，滑坡区两侧各发育1条冲沟，冲沟在坡体下部省道S306处汇聚。斜坡平均坡度17°，滑坡区分别在海拔高程940 m、1 000 m处分布2级缓坡台地(图1)。

图1 研究区地理位置及地形地貌

1.2 地层岩性

滑坡区地层岩性单一，主要以第四系昔格达组半成岩为主(图2)。半成岩是介于软岩和松散土体之间的特殊岩体，含泥量较高，具有明显的水平层理构造。根据现场钻孔结果显示，中海村滑坡山体岩性在垂直方向自上而下依次为强风化层—泥岩半成岩—泥岩(图3)，强风化层和泥岩半成岩主要呈黄褐色，泥岩呈灰白色。出露在滑坡区域周围的地层由老至新依次有寒武系川盆西南筇竹寺组—娄山关组、中上二叠统喀大崩组+奔子栏组、上二叠统扬子西南缘峨眉山玄武岩组+宣威组/黑泥哨组。

图2 昔格达组分布示意图[16]

图3 滑坡现场岩石特征

1.3 降雨特征

汉源县属于亚热带季风湿润气候，冬暖夏凉，四季分明。2010—2020年期间，研究区年平均气温16 ℃，年降雨量600 mm，其中2020年降雨量为858.9 mm (图4a)。2020年全年降雨量不均匀分布，7—8月为雨季，降雨量分别为167.5 mm和351.5 mm，约占全年降雨量的60%(图4b)。8月强降雨主要集中在滑坡发生的前7天，16日出现最大降雨，日降雨量达56.5 mm (图4c)。经实地调查，滑坡灾后现场表面出现了流水特征，认为滑坡发生前受连续降雨的影响，雨水沿裂缝结构面下渗使得动静水压力增加，饱和岩土体强度降低并引起应力的变化，导致天然岩体结构遭到破坏产生变形，而半成岩因成岩度不高且具遇水易软化膨胀的特性，对于强风化层来说渗透性较大，因此，强降雨和半成岩特性是导致滑坡灾害发生的关键因素。

(a) 年平均气温及降雨量 (b) 月降雨量 (c) 日降雨量

2 滑坡基本特征

2.1 滑坡规模及运动演变特征

滑坡形态整体呈长条形状，最大水平运动距离约600 m，最大宽度约为260 m，前后高差约180 m(后缘高程1 035 m，前缘高程855 m)，滑坡运动等效摩擦系数为0.3。滑坡启动时，约15×104m3的山体在前期连续强降雨的影响下失稳并快速滑出，随即破碎解体，在上部缓坡平台的下部，滑坡体向下加速运动并铲刮坡体上表面的松散土体，使体积逐渐增加至20×104m3； 运动过程中受地形影响，运动方向多次发生改变。滑坡体冲击掩埋了山坡平台上和山脚的民房，部分滑坡边缘区房屋严重受损，最终导致S360省道道路中断，造成7人死亡2人失联，8栋民房被毁(图5)。

(a) 滑坡前 (b) 滑坡后

2.2 滑坡分区及运动堆积特征

该滑坡发生在古滑体的残坡积层[13]。对灾后现场的勘察显示，滑坡区域呈阶梯状分布，通过当地居民获悉，1989年滑坡源区后缘曾发生滑坡。在滑坡现场，利用无人机遥感技术生成滑后DEM数据，对滑坡前后高程变化进行了差分计算(图6)。

图6 滑坡前后的高程变化

结合实地调查结果和高程变化特征，将致灾范围沿滑坡运动方向分为滑源区、堆积区Ⅰ、铲刮区和堆积区Ⅱ共4个部分。

2.2.1 滑源区及堆积区Ⅰ

滑源区位于S360省道北侧，主滑方向为N8°E，据现场测量，滑源区宽度约为165 m，滑坡后壁的最大下错距离为14 m，且壁面完整平滑，形成了典型的“圈椅状”滑坡壁，具有明显的剪切摩擦痕迹(图7)。野外滑坡痕迹及特点显示，滑源区土体失稳向下滑动，滑源区后部瞬间缺少支撑，受到岩土体向下运动时的张拉和剪切作用，使得滑源区后缘部出现2条明显的裂缝。第一条为直线形裂缝，裂缝高程为1 051 m，长度为25 m，裂缝发展方向为E40°S和W40°N，属于张拉裂缝。第二条裂缝为弧形裂缝，裂缝高程为1 043 m，长度为11 m，裂缝发展方向为东南E30°S和W40°N。由于第二条裂缝的位置临近滑源区，受到张拉力的同时也受到了滑源区坡体向下运动时的剪切力，使得第二条裂缝出现下错，最大下错距离约为1 m，属于剪切裂缝。

堆积区Ⅰ位于斜坡3级平台上，堆积方向保持不变，区域内最大堆积厚度约8 m。在滑坡发生后平台上的一栋居民房屋完全被掩埋，导致2位老人失联。

图7 滑坡裂缝

2.2.2 铲刮区

铲刮区位于2级平台上。部分滑坡体在堆积区Ⅰ开始堆积后，另一部分滑坡体沿着下垫面向下运动，运动方向发生改变，主滑方向为N24°E。该段滑坡下垫面场地为梯田，表层土体为强风化层，由于常年受到人类活动的干预，土体较为松软，且残积坡积物较厚，滑坡体向下运动时不断铲刮表面的松散土体使滑坡体体积增加，铲刮表层土体的同时携卷滑坡区两侧土体。滑坡体向下铲刮一段距离后由于地形的阻挡作用运动方向再次偏转，发生第二次堆积。

2.2.3 堆积区Ⅱ

第二次运动方向偏转后，即滑坡体进入最终堆积阶段(1级平台)。主堆积方向为N22°W，沿堆积方向的最大水平距离约为120 m，最大堆积厚度为5 m。在堆积区Ⅱ存在坡度较大的陡坎(图8)，滑坡体抛射加速后冲毁并掩埋了公路内外侧的房屋。根据滑坡的原始地形和堆积特征，中海村滑坡运动地形为阶梯型地形。阶梯型滑坡相对于凹面型滑坡和坡脚型滑坡而言，最大水平运动距离最小，但在坡脚位置具有最大冲击力[17]。因此，虽然滑坡在堆积区Ⅱ运动距离短，但导致了严重的房屋损坏和掩埋。

图8 滑坡现场

3 滑坡数值模拟

3.1 Massflow介绍

动态反演分析是在获取滑坡基本参数和地形数据条件下，借助数值分析工具进行的动态特征分析。Massflow从流体力学的N-S方程出发，基于深度积分的连续介质力学理论，提出具有二阶精度的高效能数值仿真模型，其高效且扩展自如的特点能让软件使用者根据自身需求使用Fortran语言对源代码进行添加和修改[18]。该软件可用来模拟如山体滑坡、碎屑流、泥石流、堰塞湖、溃坝、雪崩等多种山区地质灾害[19]，不仅可有效反演灾害动态过程，模拟结果还能为工程实践提供参考。近年来，该方法在地质灾害研究中取得了丰硕的成果，已应用在如香港滑坡[20]、北京房山滑坡[21]、国外的诺拉泥石流[20]、国内的红椿沟泥石流[22]以及冷渍沟泥石流[23]的工程治理效果评价中。

3.2 模型参数选取

在Massflow已构建的Coulomb-viscous模型、Voellmy模型、Manning模型这3种物理摩擦模型中，选取适用于滑坡、碎屑流和岩崩等地质灾害的Coulomb-viscous模型，该摩擦模型需考虑的物理力学参数有： 密度ρ、黏聚力c、基底摩擦系数μ、孔隙水压力系数ru、基底摩擦角δ和内摩擦角φ。以上模拟参数主要根据文献调研[24-25]和模拟校准确定： 泥岩半成岩的密度范围为1.8～2.3 g/cm3，取2.05 g/cm3，黏聚力c为28 kPa，内摩擦角φ为25°，基底摩擦系数采用等效摩擦系数(f=0.3)，基底摩擦角则利用等效摩擦系数的反正切函数计算(δ约为17°)，孔隙水压力系数控制在0.6～0.9范围。考虑到体积膨胀，将初始体积乘以系数1.2。

在孔隙水压力系数分别为0.6、0.7、0.8和0.9情况下进行试错试验，通过反向分析得到在孔隙水压力系数为0.8时，模拟结果与实际情况更加吻合。为了量化分析数值模拟堆积厚度分布，将模拟结果在滑坡主滑方向、滑源区、堆积区Ⅰ和堆积区Ⅱ进行分析(图6)，使其与航空图像和野外测量数据进行对比验证，对比结果如图9所示。由于模型中没有考虑表面微观地貌的变化和房屋对滑坡体运动的影响，剖面3-3′和4-4′的研究结果存在一定差异，但出现了较为一致的变化特色，因此，数值模拟结果具有可行性。

(a) 剖面1-1′ (b) 剖面2-2′

(c) 剖面3-3′ (d) 剖面4-4′

3.3 滑坡全运动过程模拟

3.3.1 滑坡体运动形态变化特征

采用Massflow模型的仿真模拟结果如图10所示，以10 s间隔绘制不同时间步长的滑坡运动形态，滑坡运动过程共持续约70 s。为了更直观地观察不同时间厚度的变化特征，在主滑方向进行剖切，从整个堆积过程可见，滑坡最大运动距离约600 m，滑坡体主要堆积在3个平台上，最大堆积厚度为8 m(图11)。

图10 中海村滑坡不同时刻运动厚度分布

图11 中海村滑坡不同时刻厚度分布

3.3.2 运动速度特征

滑坡运动速度的大小对整个滑坡致灾强度有着至关重要的影响。图12中显示了不同时间段滑坡体的速度分布特征。10 s时最大运动速度位于滑坡体前部，20 s时最大运动速度位于滑坡体中部，30 s时最大滑坡运动速度位于滑坡体中后部； 在此期间最大运动速度为16.2 m/s，均出现在3级平台上。40～60 s阶段，当滑坡体运动至1级平台时滑坡体前缘速度开始减小，由于惯性作用，后缘颗粒的能量传递到前缘颗粒，使得前缘颗粒继续向前运动，掩埋了对侧房屋，而此时滑源区物质向3级平台上继续堆积，直到70 s时完全停止运动。

图12 中海村滑坡不同时刻运动速度分布

在目前已有的众多滑坡各区段运动速度估算方法中，最为常见且有效的方法是谢德格尔(A. E. Scheidegger)法[26]，其计算公式如下：

。

(1)

式中：V为运动速度，m/s；g为重力加速度,通常取9.8 m/s2；H为研究点距滑坡最高点的垂直距离,m；L为研究点距滑坡后缘的水平距离,m；f为等效摩擦系数。

经计算，整个滑坡运动速度曲线如图13所示。对比2种方法分析的滑坡运动速度，在整个运动速度的演化过程中，滑源区以外的滑坡体运动速度规律基本一致且出现3个峰值，峰值速度均出现在3个缓坡平台的后部，其原因在于坡度由小到大的显著变化致使运动速度增加，从大到小的变化因滑体撞击地面而速度减小。在滑坡整体的速度演化分布上，由于没有考虑地形因素对速度的影响，导致计算值大于数值模拟值。

图13 运动速度对比

4 结论

(1)野外调查和地形高程变化特征显示，中海村滑坡前后缘高差180 m，最大水平运动距离600 m，总堆积体积为20×104m3。滑坡区域总体可分为滑源区、堆积区Ⅰ、铲刮区和堆积区Ⅱ这4个区域。

(2)降雨和半成岩特性是导致滑坡灾害发生的关键因素。滑坡所在区域地层岩性为第四系昔格达组泥岩半成岩，该岩层具有明显的水平层理构造，透水性强，遇水易软化膨胀； 滑坡发生前的连续强降雨增加了地下水的入渗，诱发了滑坡。

(3)数值模拟结果表明，孔隙水压力系数为0.8时模拟结果与实际情况吻合。滑坡运动持续时间约为70 s，最大堆积厚度为8 m，最大运动速度为16.2 m/s。

(4)中海村滑坡速度在运程上的分布演化存在3个峰值，峰值速度和主堆积区均分布在3个缓坡平台上。