黄河上游曲哇加萨滑坡“9·20”动力学过程模拟与分析

周 保，马 涛，魏正发，王文沛

（1.青海省地质环境监测总站,青海 西宁 810008；2.中国地质环境监测院（自然资源部地质灾害防治技术指导中心）,北京 100081）

0 引言

曲哇加萨（军功）滑坡位于青海省果洛藏族自治州玛沁县拉加镇，是发育于黄河上游新近系岩层中的巨型老滑坡。滑体前缘长期受黄河强烈侵蚀切割，以及修建公路、坡脚削坡建房等人类工程活动影响，滑体变形迹象明显，历史上曾多次发生局部滑动（图1）。

图1 玛沁县拉加镇曲哇加萨滑坡位置图Fig.1 Location map of Quwajiasa landslide in Laga Town,Maqin County

1985年7月21日老滑坡体中段发生滑坡。2011年8月12日当地居民在老滑坡体前缘削坡建房引发局部滑动。2013年6月1日、2018年5月2日受降雨影响，老滑坡体局部变形加剧。2019年9月20日老滑坡体上H1、H2、H4 滑坡发生滑动，造成国道G227上山路段严重变形，损毁公路400 m，前缘18 户居民受灾，51 间房屋倒塌。同时，造成老滑坡体中前部，即拉加镇军功路南侧坡体出现变形迹象。2020年3月1日，老滑坡体中前部拉加镇军功路南侧坡体前缘局部变形加剧，威胁住户120 户、商铺65 户，其中19 户139 间房屋出现不同程度倾斜、地基下沉、墙体开裂等情况[1−4]。

针对曲哇加萨滑坡复杂的变形破坏特征，程强等[5]提出：这种特殊的新近系以来沉积的高原红层，成岩时间短，力学强度相对较弱，易发生红层与上覆土体界面的大型滑坡灾害。张永康等[6]通过现场勘查提出青海高原红层滑坡的具有多区、多级、多层的变形破坏特征，且各滑体滑动速度不同。吴永刚等[7]通过物理模型试验认为青海高原红层滑坡河谷下切引起的坡体卸荷回弹变形使得浅层滑面位移大于深层滑面位移，雨水浸润引起材料强度软化进一步加剧了变形。殷志强等[8]指出：黄河上游巨型滑坡主要发育于气候的温暖湿润期和气候变化的快速转型期，具有多期次滑动过程特征。Wang 等[9−13]利用安全系数云图分析法模拟出了茂县梯子槽高位滑坡多层滑面安全系数分布，并据此进行防治结构设计，提出了针对性的小口径组合桩群治理方案。李滨等[14−19]提出此类特大高位地质灾害易形成泥石流、堵江等灾害链，需加强调查与识别能力、监测预警与风险防范能力以及防灾减灾能力建设。

文章以2019年9月20日19 时发生于曲哇加萨老滑坡东北段中前缘 H1 滑坡为例，在野外调查、形变数据以及稳定性分析的基础上，研究该滑坡的变形破坏特征，并进一步通过动力学模拟分析滑坡的成灾范围，可为后续的综合防治方案提供科学参考。

1 研究区地质环境条件

曲哇加萨老滑坡所处位置为黄河强烈下切侵蚀区，河谷形态为峡谷。拉加镇段因黄河呈急剧“S”拐弯形成相对开阔的一小型河谷型盆地，其凸岸发育有黄河Ⅱ、Ⅶ级阶地，阶面地形较为平缓。凹岸坡体和河床直接接触，长期遭受强烈的侵蚀切割，形成高几十米至数百米不等的陡坡、陡崖。滑坡区地形地貌主要为黄河Ⅱ、Ⅶ级阶地组成的河谷平原地貌。黄河Ⅱ级阶地分布于滑坡前缘，现为居民区，黄河Ⅶ级阶地分布于老滑坡后缘。老滑坡西侧坡脚长期受黄河冲刷、侵蚀下切（图1）。

滑坡区北侧发育有拉家压扭性逆断层（F1），该断层位于拉加北山，NE60°方向延展，为逆冲断层，断层破碎带宽度10～30 m。挽近时期的隆拗运动在滑坡区的表现也较明显，其隆拗的长轴方向继承了老构造断裂带的走向，多呈东西向展布。该地区地震活动频繁，地震基本烈度Ⅶ度。

老滑坡整体呈宽簸箕形，坡体发育有3～6 级滑移台坎。东、南侧滑坡后壁明显，呈陡坡状，局部为陡崖，高度40～50 m，西侧滑坡边界以深切的塔尔隆沟为界；滑体宽约1 900 m，长700～900 m，厚30～100 m，方量约1.67×108m3，整体坡度约25°，主滑方向307°。滑体后缘高程3 225～3 340 m，前缘高程3 040 m，相对高差185～300 m。老滑坡后缘出露地层上部为“二元结构”黄土状土和卵石，下部为新近系泥岩，泥岩产状NE5°∠25°，组成顺向缓倾坡（图2）。

图2 曲哇加萨滑坡平面图Fig.2 The engineering geological plane of Quwajiasa landslide

2 2019年9月20日滑坡特征

2019年9月20日，曲哇加萨老滑坡东北段中前缘H1、H2、H4 滑坡发生滑动，造成路面隆起、损毁公路400 m，18 户居民受灾，51 间房屋倒塌。H4 滑坡滑动后造成原来公路的20 根抗滑桩裸露，悬臂5～10 m，出现桩间土坍塌变形，桩后 H3 滑坡内也形成多条纵向裂缝。灾害发生后，采取了拆除房屋和回填压脚措施[2]（图3—6）。

图3 曲哇加萨滑坡发生前后对比Fig.3 Comparison before and after the occurrence of Quwajiasa landslide

图5 前缘鼓胀引起房屋倒塌Fig.5 The swelling of the front edge caused the houses to collapse

由图4—5 可以看出，临近公路位置H1、H2 滑坡前缘发育多处土体解体，并挤压公路产生多处放射状鼓胀张拉裂缝，裂缝宽度30～60 cm，深度2 m，裂缝走向与滑坡方向平行或呈小角度相交。公路外侧挡墙产生严重的鼓胀变形，裂缝宽度达5 cm。公路下方由于滑坡滑动鼓胀导致地表隆起，造成51 间房屋倒塌。滑坡前缘影响范围至公路挡墙和军功路之间。

图4 曲哇加萨“9·20”滑坡前缘Fig.4 Frontal edge of “9·20”landslide in Quwajiasa

由图6 可以看出，H1 滑坡后缘陡坡下错约2.0 m，侧界清晰，滑坡呈现蠕滑特征，表面裂缝遍布。现场后缘可见水体入渗迹象，土体含水量较高。

图6 曲哇加萨“9·20”H1 滑坡后缘Fig.6 Back edge of “9·20”H1 landslide in Quwajiasa

为避免滑坡进一步变形致灾，灾害发生后，采取了拆除房屋、回填压脚及截排水措施，9 天后，该滑坡逐渐趋于稳定状态。

实际上H1、H2、H4 历史上曾出现多次变形。2011年8月12日，由于省道S101 线（现为国道G227）修建时开挖该区西侧滑体前缘和当地居民削坡建房等工程活动，引发老滑坡前缘部分滑动，使得H1、H2 滑坡后缘和右侧缘形成连续的圈椅状陡壁，高2～6 m。滑坡后缘拉张裂缝密集发育。

根据2013年8月—2014年10月地表变形监测数据监测数据分析，监测点的位移量在2014年6月、2014年8月出现两次明显阶跃，最大水平变形累计达到900 mm[1−4]。根据2014年7月—2019年6月时间序列InSAR 监测数据，获得滑坡的年平均形变速率超过70 mm，说明曲哇加萨滑坡一直在变形。

从图7 可以看出，8月27日—9月22日近1 个月累计降雨量为91.5 mm，其中18—20日3日连续降雨量为23.5 mm，占比25.7%。降雨沿着密集发育的裂缝下渗，加速了地下水的渗流作用，进一步降低了岩土体的强度，最终导致北侧滑坡发生滑动。

图7 曲哇加萨滑坡2019年8月27日—9月22日降雨量Fig.7 Rainfall from August 27 to September 22,2019

为更好说明东北段中前缘滑坡发生的相关机理，选择H1 滑坡进行具体分析。结合物探、钻探及探井资料，H1 滑坡由上至下共发育4 层滑面（图8），从上至下分别为滑面1：位于碎石土与粉质黏土交界层前部，深度约5～10 m，饱水，呈泥团状；滑面2：位于碎石土与粉质黏土交界层附近，深度约20～30 m，含水率较高，呈软塑状，局部可见擦痕；滑面3：位于黏土层与粉质黏土交界层处，深度约35～50 m，含水率较高，呈可塑-软塑状，可见擦痕；滑面4：位于基岩与黏土层的交界处，深度约60～80 m，含水率较高，呈硬塑-可塑状。

图8 曲哇加萨滑坡典型工程地质剖面（A—A′）Fig.8 Typical engineering geological profile of Quwajiasa landslide (A—A′)

选取图8 典型剖面进行滑坡稳定性计算，确定降雨对滑坡稳定性的影响程度，根据程柯力等[2]计算结果表明浅层滑面1、滑面2 的稳定系数为0.94、1.02，处于不稳定状态和欠稳定状态，易于继续发生变形破坏（图9）。模拟结果很好验证了此次东北段中前缘滑坡主要是由于浅层滑面1 蠕滑形成，同时滑面2 的滑动可能性也较高，需做相关的动力学预测分析。而深层滑面3、滑面4 的稳定系数均为1.35，处于稳定状态，故不需做相关动力学预测分析。

图9 滑坡典型剖面安全系数分区Fig.9 Safety factor zoning of typical landslide profile

3 滑坡动力学反演模拟及预测评价

3.1 DAN-W 基本原理

为了分析滑坡剪出后动力学特征和评估成灾范围，采用加拿大Hungr 教授开发的DAN-W 二维的动力模拟方法进行正演分析。DAN-W 是一种基于Windows程序，在连续介质模型基础上将滑体等效为具有流变性质的流体，选用不同的流变模型，通过设定滑坡的滑动路径的参数，从而达到模拟滑坡的运动速度、时间、路程以及堆积体特征效果[20−29]。大量的模拟结果表明摩擦准则和Volleymy 准则最能表达滑坡的运动。

摩擦准则是一个单变量的流变准则，其抗剪强度表达式为：

式中：τ——滑坡底部剪应力/Pa；

σ——垂直运动方向的总应力/Pa；

ru——孔隙水压力与总正应力之比；

φ——摩擦角/（°）。

Volleymy 准则的抗剪应力表达式为：

式中：f——摩擦系数；

ρ——滑坡的密度/（kg·m−3）；

g——重力加速度/（m·s−2）；

v——滑坡平均速度/（m·s−2）；

ε——湍流系数/（m·s−2）。

该准则中f和ε 为两个待定的参数。

3.2 模型建立及反演模拟对比

根据现场调查发现，滑坡主要运动模式是蠕滑，并且前部有公路挡墙和房屋阻挡，并未发生远程滑坡，因此选用摩擦流变模型较为合适。当然，为更好的说明摩擦流变模型的有效性，再选取Volleymy 流变模型共同与H1 滑坡1 号滑面浅层滑动现场调查结果进行对比。采用试参法对摩擦流变模型进行模拟[20]，Volleymy流变模型内摩擦角与摩擦流变模型一致，同时，摩擦系数，湍流系数采取工程类比，从而对已发生滑动的曲哇加萨1 号浅层滑面进行运动反演模拟对比分析（表1）。

表1 模型参数选取Table 1 Model parameters

从图10 可以看出，1 号滑面浅层滑坡运动后，若按照摩擦流变模型将运动堆积到623 m 位置，即挡墙至军功路段国道G227 和居民区，暂时不会影响到军功路至黄河段。同时，滑坡体积超过50%基本停留在上部滑面上，其余均匀停留在运动路径上，这与现场调查滑坡堆积特征基本相符。若按照Volleymy 流变模型运动堆积到682 m 位置，即军功路至黄河段居民区，与现场调查结果不符，大于实际运动距离。同时，滑坡体积主要停留在坡脚位置，也与实际情况不符。

图10 1 号滑面浅层滑坡运动前后剖面形态对比图Fig.10 Comparation of longitudinal cross-section of No.1 sliding mass before and after sliding

3.3 欠稳定状态2 号滑面运动预测评价

利用上述反演分析摩擦流变模型的同一套参数（表1）对处于欠稳定状态的2 号浅层滑面进行运动预测模拟分析，得出了堆积厚度、范围等运动特征，并与上述1 号浅层滑面进行对比（图11）。

图11 2 号滑面浅层滑坡运动前后剖面形态图Fig.11 Longitudinal cross-section of No.2 sliding mass before and after sliding

研究发现：2 号滑面浅层滑坡一旦运动后，将运动堆积到673 m 位置，即军功路至黄河段居民区。滑坡体均匀停留在运动路径上。模拟结果说明：目前2 号滑面浅层滑坡一旦发生滑动，威胁军功路至黄河段居民区及商铺房屋，但不至引起堵塞黄河灾害发生。

4 结语

（1）曲哇加萨滑坡坡体呈现多区、多级、多层的变形破坏特征，目前变形破坏强烈，尤其是东北段中前缘滑坡近期频繁出现变形破坏。

（2）东北段中前缘滑坡地表监测曲线、形变速率曲线、钻探数据和现场调查表明，降雨沿着密集发育的裂缝下渗，加速了地下水的渗流作用，进一步降低了岩土体的强度，最终导致滑坡发生变形滑动。

（3）以H1 滑坡为代表，在降雨稳定性计算基础上进行动力学反演，反演分析发现1 号滑面浅层滑坡按照摩擦流变模型运动后，将运动堆积到挡墙至军功路段国道G227 和居民区，这与现场调查滑坡堆积特征基本相符，而Volleymy 流变模型运动堆积到682 m 位置，即军功路至黄河段居民区，与现场调查结果不符，大于实际运动距离。同时，滑坡体积主要停留在坡脚位置，也与实际情况不符。

（4）利用上述反演分析摩擦流变模型的同一套参数对处于欠稳定状态的2 号浅层滑面进行运动预测模拟分析，发现2 号滑面浅层滑坡运动后，将运动堆积到军功路至黄河段居民区，但不至于引起堵塞黄河灾害发生。

（5）本文的研究内容可为黄河上游该类滑坡运动预测评价方法提供一定参考。