白龙江流域大型高位滑坡成灾动力过程模拟研究*

冯 振，游 杨，陈 亮，王立朝

(1.中国地质环境监测院，北京 100081；2.湖北省地质局第二地质大队，恩施 445000)

白龙江流域位于青藏高原东北缘、西秦岭构造带，长期地质构造发展过程中沿北西构造线形成与白龙江大致互相平行的挤压带，构造运动强烈，活动断裂发育，曾发生过震级较大的地震(图1)。

图1 立节北山滑坡构造环境背景(断裂带数据来源于参考文献[1])

受秦岭纬向构造体系的影响，历经燕山运动和喜马拉雅山造山运动，白龙江复式背斜在长期挤压、扩张、褶皱和不断复合过程中，形成多期性断裂构造，次级褶皱非常发育，岩体极其破碎。区域性深大断裂-光盖山-迭山南麓断裂带由多条分支走滑逆冲断裂组成，分支断裂断错在泥盆系与志留系之间。光盖山-迭山南麓断裂带总体走向300°，倾向NE，倾角60°～80°，其左旋走滑速率约1.4 mm/a、挤压逆冲速率约3.7 mm/a。

地质构造发育，地层岩性复杂，水系切割强烈，地形具有山势陡峻、相对高差大等特点，使白龙江流域发育大量高位滑坡灾害，并且具有大型滑坡成群集中的特点，是我国高位地质灾害风险极高的地区之一[2-4](图2)。大型高位滑坡发育于高陡斜坡中上部，蕴含巨大的重力势能，滑坡失稳发生重力势能与动能转化，极易形成高速远程滑坡-碎屑流灾害，具有超远距离的滑动位移、强大的破坏与致灾能力[5-10]。调查发现，白龙江流域舟曲段分布由45处大型-特大型高位滑坡，包括咀疙瘩滑坡、锁儿头滑坡、泄流坡滑坡、南峪乡江顶崖滑坡等，曾发生过多起高位高速远程滑坡堵江事件，造成重大的经济损失[11-15]。因此，有必要在该区域开展大型高位滑坡风险调查评价工作，为保障山区城镇化与重大工程建设安全提供技术支撑。

图2 白龙江流域舟曲段滑坡分布及立节北山滑坡位置图(底图来源于Bigemap软件，灾害点数据来源于参考文献[16])

立节北山古滑坡位于光盖山-迭山南麓断裂带分支断裂上盘、白龙江左岸的立节镇北侧山体上部，体积约240万 m3，高程2 210～2 430 m，剪出口与坡脚的立节场镇、白龙江江面高差约700 m，属于典型的大型高位滑坡(图3)。本文以白龙江流域舟曲立节北山滑坡为例，通过资料搜集、遥感调查与解译、野外踏勘等手段，查明了滑坡所处的地质环境条件和变形破坏特征，采用基于光滑质点流体动力学与等效流体模型的DAN3D动力分析软件，开展滑坡后破坏运动过程模拟，获得运动学参数与动力致灾特征。

图3 立节北山滑坡全貌(2021年9月无人机倾斜摄影DOM)

1 地质环境条件

1.1 地形地貌

舟曲县立节镇北山滑坡区域构造上位于青藏高原东缘、西秦岭西翼与岷山山脉交会地带。古滑坡所在山体地形北高南低，坡顶高程约2 780 m，白龙江江面高程约1 510 m，最大高差近1 270 m，属于强烈上升的侵蚀构造地貌(图3)。依据坡度、高程及变形破坏特征，可将山体划分为5段(图4)：高程大于2 440 m为稳定缓坡区，坡度14°～19°，主要为耕地；高程2 210 m～2 440 m为古滑坡变形区，坡度14°～60°，呈整体坡度32°，上陡下缓的地形，为梯级旱地；高程1 870 m～2 210 m为浅表层滑坡区，坡度34°～41°，基岩裸露，分布零星灌木杂草；高程1 600 m～1 870 m为坡面侵蚀沟谷区，坡度41°～52°，沟道呈“V”字形；高程小于1 600 m为白龙江河谷阶地，坡度0～17°，为建筑用地及农田(图2)。

图4 立节北山滑坡剖面A-A’(剖面位置见图2)

1.2 地层岩性

1.3 岩体结构

图5 立节北山滑坡工程地质条件图(底图为2021年1月无人机倾斜摄影DOM)

2 滑坡体基本特征

立节北山滑坡具有典型的凹陷负地形地貌，其平面形态为圈椅状，剖面呈阶梯状，边界清晰(见图5)。调勘察结果显示，立节北山滑坡前缘高程为2 223 m，后缘高程为2 432 m，相对高差209 m，主滑方向203°，滑体最大长度391 m，平均宽度362 m，面积约14.15×104m2，平均厚度17 m，总体积约240×104m3，为中层大型堆积层滑坡。

滑坡后缘与山体后部稳定缓坡区发生地形转折，坡度陡峭，坡度约40°～45°，局部出露厚层状灰岩，小型崩塌时有发生(图5、图6)。滑坡中下部坡度约20°～28°，为条带状梯田，各梯田高差3～5 m，陡坎坡度68°～84°。东西两侧为与后壁相连的近南北向山梁，前缘均以为冲沟界，滑体变形导致两侧道路断陷，东侧前缘黄土覆盖区揭露明显的陡立滑动面(图7)。滑坡前缘地形陡立，坡度32°～45°，发生大面积的浅层滑动，地表凌乱，分布大量的裂缝与错坎(图5、图8)。

图6 小型崩塌-碎屑流(镜向345°，位置见图5点a)

图7 立节北山滑坡东侧边界(左图镜向330°，右图镜向40°，位置见图5点b)

图8 立节北山滑坡前缘变形(镜向280°，位置见图5点c)

3 变形破坏特征

立节北山滑坡长期处于蠕滑变形，滑体上分布规模各异的裂缝，局部浅层滑动、崩塌不断发生，强降雨与地震的变形加剧效应尤为显著[17]。1979年滑坡首次发生复活变形，前缘东侧发生局部滑动，经长期降雨侵蚀形成深切冲沟。1982年滑坡前缘出现裂缝，中部民房变形开裂。1984年、1998年强降雨导致前缘不断发生小面积滑动。2008年汶川地震使滑坡前缘浅层滑动加剧，东侧冲沟出现小型黄土崩塌。2013年受岷漳地震及降雨影响，滑坡西侧前缘发生浅层滑动，东侧冲沟坡面不断崩塌滑落。2014-2020年，滑坡各个区域变形不断加剧，局部崩滑时常发生。直至2020年8月17日暴雨，滑坡前部发生大面积滑动，滑体沿坡面冲沟转化为泥石流冲入白龙江，造成白龙江拥堵过半。2021年1月14日，立节北山滑坡变形加剧，坡体多处出现裂缝合错落，滑体前部村道局部变形破坏(见图5)。

根据无人机倾斜摄影与现场调查，截至2021年底，立节北山滑坡范围内形成8处崩塌、39条裂缝、8个次级滑动(图5)崩塌源自滑坡局部出露的厚层状碎裂灰岩，岩体风化程度高，节理裂隙发育，垂向裂隙切割岩体发生错落式崩塌，沿陡坡形成碎屑流，最终堆积于梯田、陡坎与道路上(图5、图6)。裂缝倾向临空方向，倾角近直立，均为张拉性质，根据发育分布特征可分为三组(图5)。一组沿滑坡边界展布，走向与滑坡边界基本平行，包括：后缘4条弧形裂缝，平均宽度15 cm，平均深度50 cm；西翼3条、东翼13条羽状裂缝，平均宽度15 cm，平均深度70 cm，局部有错动陡坎，高差5～20 cm，其中滑坡下部东侧侧壁处裂缝落距约2.1 m，在表层黄土中揭露平直的滑动面。一组位于滑坡中部，共4条，长25～70 m，整体走向北西，其中3条构成西侧次级滑动的东侧边界。一组位于滑坡前缘，共16条，长20～120 m，整体走向与滑动方向近垂直，是前缘3处局部滑动形成的张拉裂缝。

立节北山滑坡属于古滑坡堆积体复活变形，具有典型的土质滑坡特征与要素，包括圈椅状地貌、圆弧形滑动面、“上陡下缓”的地形等特点，后缘拉张错断、前缘塌陷溜滑、地表裂缝交错等各类变形迹象，是滑坡蠕滑变形的结果。滑坡下部斜坡由千枚岩、板岩构成，强度低、性质软、风化程度高，坡面冲沟发育，长期溯源侵蚀作用导致沟坡及浅表层滑动不断发生，规模与面积逐渐扩大，形成多处次级滑动(图9)。由于前缘逐级失稳，滑体阻滑力不断减小，滑坡发生渐进破坏、产生蠕滑变形，受局部地形控制，地表发育规模与走向各异的裂缝，随着裂缝贯通形成2个次级滑动。总而言之，立节北山滑坡变形是前缘坡面冲沟溯源侵蚀引起的高位古滑坡堆积体复活，在自然条件下持续蠕滑(图10)，强降雨与强震可能导致整体剧烈失稳，形成高速远程滑坡，对坡脚的立节镇及白龙江造成重大威胁。

图9 立节北山滑坡前缘次级滑动(2021年9月无人机倾斜摄影三维立体影像)

图10 滑坡前缘裂缝扩展特征(位置见图5点d)

4 动力致灾过程模拟

4.1 DAN3D基本原理

为了预测立节北山滑坡整体失稳后的滑动运动过程，选择DAN3D软件开展滑坡动力学数值模拟。DAN3D软件在连续介质模型的基础上，将地质灾害等效为具有流变性质的流体，采用基于光滑质点流体动力学(Smoothed Particle Hydrodynamics，SPH)方法，对结合流体深度的圣维南方程进行求解，运用连续动态方程，在复杂三维地形上模拟地质灾害运动、铲刮及堆积过程，实现反演或预测分析[19]。根据前人研究结果表明，在流变模型的选取上，Frictional、Voellmy两种流变模型对于大型高位远程滑坡的运动过程模拟有较好的结果[20-23]。两种模型的流变关系式如下。

1)Frictional流变模型：该模型适用于计算颗粒材料在基底阻力作用下的情况。当流体基底阻力T仅作用于基岩上有效的正应力时，可以使用该模型。在这种情况下，应力与滑体高度、单位重力和孔隙水压力相关，其阻力关系如下：

T=AγH(cosα+ac/g)(1-ru)tanφ。

(1)

其中：T为流体基底阻力；A为底面积；γ为流体重度；H为流体厚度；ac为滑面的离心加速度；ru为孔隙水压力系数；φ为内摩擦角；g为重力加速度。

2)Voellmy流变模型：该模型由颗粒摩擦和湍流运动两种运动形式组合而成，主要适用于雪崩、泥石流、滑坡碎屑流等。

T=A[γH(cosα+ac/g)tanφ+γv2/ξ]。

(2)

其中：ξ为湍流系数；v为流体运动速度；其他参数与公式(2)中一致。

4.2 模型建立与参数选取

为了建立DAN3D运动分析数值模型，首先需要得到三个地形数据文件，包含滑移路径地形(Path Topography)、滑体地形(Source Topography)和侵蚀地形(Erosion Map)。

1)滑移路径地形：滑体运动的下垫面地形。基于等高线生成的数字高程DEM，导入到Surfer软件，生成DAN3D软件可运行的格式(.grd)。

2)滑体地形：即滑源区滑体厚度。根据调查得到的滑坡边界、钻孔揭露的滑体厚度，利用Surfer软件，将滑前和滑后DEM数据的对比，从而得到滑源区的滑体厚度。此次模拟的滑体不考虑滑坡前缘的浅层次级滑动，主要考虑主变形区的滑体(图11)。

图11 滑坡厚度模型(图例与图5一致)

3)侵蚀地形：用于定义滑移路径的下垫面材料分布。根据野外调查和工程地质图(图3)，确定地层岩性的分布。

针对高速远程滑坡滑源区、铲刮流动区、堆积区的分区特点，通过两种流变模型的组合对比发现，F-V-V模型在立节北山滑坡的成灾动力过程模拟中具有较好结果，相关的模型计算参数见表1。

表1 模型计算参数

4.3 模拟结果分析

模拟结果显示，立节北山滑坡整体启动转化为碎屑流，经过高速滑动到堆积完全停止，滑坡运动过程可分为以下三个阶段(图12)。

图12 运动过程滑体厚度变化(底图为无人机摄影DSM数据形成的山体阴影)

1)0～20 s为滑坡启动阶段(图12a)。滑体沿滑面开始向下运动，逐渐脱离滑动面，自2 250 m高程剪出。由于滑动面呈“上陡下缓”的形态，整体倾角约25°，因此滑源区启动速度相对较小(0～15 m/s)，西侧受坡度影响速度较大，在前缘剪出口加速至约30 m/s。

2)20～100 s为高速铲刮阶段(图12b)。滑体脱离滑动面后沿坡面及冲沟向下运动，转化为高速碎屑流。滑体沿途刮铲并挟带浅表层的松散岩土体，使滑体体积快速增大，毁坏沿途植被、耕地等。由于铲刮区斜坡高度大、地形陡，滑体赋存的巨大势能转化为动能，保持加速运动状态，运动速度达24～52 m/s(图13)。由铲刮厚度等值线(图14)可知，碎屑流高速铲刮平均厚度约3 m，最大厚度约5.3 m，体积扩容约77.7万 m3。

图13 最大运动速度分布图

图14 铲刮等值线图

3)100 s后进入减速堆积阶段(图12c)。碎屑流前缘到达坡脚高程约1 600 m处，因地势趋于平缓开阔，碎屑流受到的动摩阻力增大，速度逐渐减小，同时因前缘受阻后，碎屑流向两侧扩散，进入减速堆积状态(图13)。碎屑流前缘最终到达白龙江右岸，约200 s时完全停止运动。据滑体厚度变化图(图12)所示，高速远程滑坡-碎屑流最大运动距离约1 600 m，在坡脚后整体呈扇形堆积，覆盖面积0.22 km2，堆积体积约317.7万 m3，平均厚度7.19 m，最大厚度17.8 m。

国际上一般将大于5 m/s的滑坡归为极速滑坡，将滑坡体重心位置的垂直位移与水平位移的比值H/L(等值摩擦系数)小于0.6的滑坡称为远程滑坡[24]。上述分析可知，立节北山滑坡若失稳启动将形成典型的高速远程滑坡-碎屑流，滑坡等值摩擦系数约0.49，滑体最大运动速度45.7 m/s，具有对立节城镇建设与人民生命财产安全造成严重的灾难的能力。

5 结论和讨论

1)以白龙江流域大型高位滑坡-甘肃舟曲立节北山滑坡为研究对象，在调勘察与遥感解译的基础上，查明了滑坡地质环境条件和变形破坏特征，认为滑坡蠕滑变形是前缘坡面冲沟溯源侵蚀引起的高位古滑坡堆积体复活的结果。

2)立节北山滑坡前后缘高差大，剪出口距坡脚高差极大，拥有巨大的势能，具有转化为高位高速远程滑坡的有利地形和动力条件。近年持续显著的复活变形说明滑坡具有较大的剧滑成灾风险，对坡脚立节镇人民和财产造成严重威胁。

3)基于DAN3D软件，对立节北山滑坡潜在高速远程运动过程进行了模拟预测。滑坡运动过程可分为失稳启动、高速铲刮、减速堆积三个阶段。初始滑体方量240万 m3，滑程约为1 600 m，最大运动速度为45.7 m/s，沿途平均铲刮厚度3 m。滑坡最终堆积体积317.7万 m3，平均厚度7.19 m，最大厚度17.8 m，覆盖面积0.22 km2，约占立节场镇面积的50%。运动过程模拟与分析可为滑坡危险性分区提供理论依据，为类似的高位滑坡风险预测与评价提供参考。