0.   引 言

地震滑坡是对人类生活影响最大的一类滑坡，也是自然界中最强的滑坡^[1]。例如在2008年汶川地震期间，北川县城城西发生的滑坡造成了1600人遇难，数百间房屋遭到破坏，这是汶川地震引发的最为严重的滑坡灾害，极为罕见^[2]。而高位远程滑坡是指运动高差大、剪出口高、滑程远，并在特定情况下可引发堰塞坝和涌浪等次生效应的地质灾害^{[3 − 4]}。在不同诱导因素下，会形成不同种类的高位远程地质灾害^[5]，这就导致其通常包含多种基本地质灾害。历史上，在地震因素的诱导下发生的高位远程地质灾害破坏性极大。例如1950年，在察隅地震的影响下，南迦巴瓦峰西侧的则隆弄冰川发生了大规模且迅速的移动，导致下方的村庄被掩埋，雅鲁藏布江也因此被阻断^[6]。

近年来，国内外大量学者投身于地震滑坡动力响应的研究。地震诱发滑坡的失稳特征的研究方法有多种，包括动力时程分析法、动力强度折减法^{[7 − 8]}、拟静力法^{[9 − 10]}等。周兴涛等^[11]使用有限差分法研究了地形地貌对地震波的放大效应规律。祁生文等^{[12 − 13]}使用FLAC^3D软件研究了坡高、坡角对单面均质边坡动力响应的影响。王文沛等^[14]以高陡斜坡为研究对象，运用FLAC软件探究了高陡岩质斜坡地震失稳机制。刘铮等^[15]运用FLAC软件对双面斜坡的动力响应规律进行了研究，得到了坡体的卓越频率并不是唯一的，在不同部位其卓越频率不尽相同。言志信等^[16]通过有限元法研究了边坡在双向地震作用下的共振规律以及固有频率。杨国香等^[17]利用大型振动台进行实验，研究了大型岩质边坡，分析了顺层和均匀结构岩质边坡的动力加速度响应特性，以及动力输入参数如何影响边坡的特性。Wu等^[18]根据滑坡发生密度的分析，发现滑坡发生与地震和地貌因素之间存在强相关性。Mitani等^[19]通过Abaqus有限元软件分析和研究了在不同坡高、坡角的均质斜坡的地震放大效应。Hovius等^[20] 和Marc等^[21]发现瞬时地震动后，地震影响逐渐衰弱，但地震的累积效应会削弱区域结构稳定性。Fan等^[22]在Hovius等^[20]和Marc等^[21]的基础上研究了在地震滑坡发生后，地区地质灾害的易发性。蔡国军等^[23]使用离散元法，对强震作用下的岩质边坡表面放大效应进行研究，结果表明边坡拐点放大系数呈较大值。虽然国内外已经在地震滑坡领域做了非常多的工作，但是在地质环境复杂、沟谷深切、高山常年积雪、断裂带遍布的西藏高山峡谷地区，对启动区达到数千米的超高位地震滑坡的探究还有待完善。

综上所述，地震诱发的超高位远程地质灾害规模巨大，严重威胁西藏人民人身和财产安全。本文以西藏则隆弄沟为例，采用FLAC^3D的动力模拟方法，探究超高位边坡在强震动力情况下的失稳特征，对于与则隆弄地区类似的山区地震动响应研究起一定的启示意义。

1.   1950年震前剖面绘制

由于数值模拟需要用到1950年察隅地震前的地质剖面，所以需要根据前人的记载对则隆弄沟剖面进行恢复。接下来将从多次自然灾害的记录中对则隆弄冰川冰舌的厚度和海拔进行反推和绘制。

1.1   则隆弄地质特征

则隆弄沟位于西藏自治区林芝市米林县派镇直白村旁，该沟由两条沟组成，分别为主沟和支沟，物源形成区呈现“Y”字形。沟域形状类似沙漏形，面积约为53.61 km²，堆积区附近有一个村庄（图1）。

图  1  则隆弄沟平面卫星图（引自omap 2023）

Figure  1.  The aerial map of Zhelongnong gou (from omap 2023)

下载: 全尺寸图片 幻灯片

主沟的走向为北西—南东，沟口汇入雅鲁藏布江。河床地面高程为2847 m，坡顶高程约为7782 m，高差为4935 m。主沟的长度可达12.48 km，平均坡降为386‰，最大坡降为582‰^[24]。斜坡主要由南迦巴瓦岩群上亚群的派乡花岗片麻岩组成，沟内有大量的冰碛物堆积体，随着海拔的上升冰碛物中的冰川冰不断增加。沟底发育有走滑断裂带、接地-当噶断裂。

则隆弄沟所在区域处于欧亚板块和印度板块强烈碰撞带上，地质构造活动非常强烈，导致地震十分频繁^[25]。多期次的造山运动过后，该区域地势高差不断增加，形成了非常典型的高山峡谷地貌。所处区域发育有青藏高原喜马拉雅东构造结，其中最活跃的断裂为东久—米林断裂，断裂方向为西南—东北。

1.2   历史灾害记载

则隆弄沟有记录以来共发生过4次较大规模的自然灾害，发生的年份分别是1950年、1968年、1984年、2020年^[26]。1950年8月15日，墨脱8.5级大地震，则隆弄冰川受到影响，冰川破碎并加速滑下。次日，人们发现，在雅鲁藏布江大峡谷入口，由于冰体快速移动形成了一个数十米高的冰坝，这导致江水暂时断流。在冰川跃进期，支沟的冰体未能及时流入主沟，反而直接越过海拔4000 m的新冰期和小冰期的侧碛堤，沿北侧边沟迅速下滑，最终摧毁了直白村，将之夷为平地。1968年9月2日午后，由于气温上升使得冰雪融化，冰川迅速前进，导致雅鲁藏布江再度受阻，形成了一个数十米高的冰坝。沟口的一座木桥因此被水流冲毁并埋没，次日上午冰堤决堤。这一次的冰川活动与1950年的首次跃动在空间上的变化程度相似。1984年4月13日，自从1968年冰川活动以后，主要的冰流由于差异性融化被划分为6个部分。在这些部分中第6段发生滑动，结果再次阻塞了雅鲁藏布江^[26]。2020年9月10日，则隆弄沟又暴发泥石流，致使直白大桥损毁严重，但未造成堵江。

1.3   冰舌厚度及海拔推断

根据现有记录可以直接得到1950年冰舌所处海拔为3650 m^[26]，但是缺少冰川冰碛物的大致厚度信息。通过不同时期冰川冰碛物厚度记载，可以进行合理范围内的厚度推断。

如表1所示，在不同时期冰舌厚度均有所不同，其中最薄的时期为距今400～1500年前，厚度为150 m。最厚的时期为末次冰期，厚度为300 m，而现代厚度为200～300 m。运用数理统计规律及冰川的变化规律，结合过去时期和现在的厚度，能大致估计得出一个比较合理的厚度，约为200～300 m。于是可以得到以下剖面图（图2—3）。

表  1  不同时期冰舌海拔位置及冰川冰碛物厚度

Table  1.  The thickness and elevation position of the ice tongue at different times

时期	时间	冰舌海拔/m	冰碛物厚度/m
末次冰期	1万年前	2850	300
新冰期	1500年前	2 920～4 000	>150
小冰期	400年前	2 950	>150
21世纪	2020年	3 750	200～300
21世纪	2022年	3 750	200～300

下载: 导出CSV 

| 显示表格

图  2  不同时期冰川冰碛物厚度及冰舌位置剖面图

Figure  2.  Cross-sectional diagram of glacial moraine thickness and ice tongue position at different times

下载: 全尺寸图片 幻灯片

图  3  1950年则隆弄沟震前地形剖面图

Figure  3.  Topographic profile of the Zelongnong gou ditch before the earthquake in 1950

下载: 全尺寸图片 幻灯片

2.   力学模型及动力输入

2.1   力学模型及参数

根据现场资料可以将力学模型简化为两种材质的黏弹性模型（图4），岩样1（花岗片麻岩）：重度$\rho =27.0{\mathrm{k}\mathrm{N}/\mathrm{m}}^3$，弹性模量$E=80\mathrm{G}\mathrm{P}\mathrm{a}$，泊松比（λ）为0.3，黏聚力$c=20\mathrm{M}\mathrm{P}\mathrm{a}$，内摩擦角$\phi =36\text{°}$；岩样2（冰碛物）：重度$\rho =19.0{\mathrm{k}\mathrm{N}/\mathrm{m}}^3$，弹性模量$E=9\mathrm{M}\mathrm{P}\mathrm{a}$，泊松比（λ）为0.3，黏聚力$c=50\mathrm{k}\mathrm{P}\mathrm{a}$，内摩擦角$\phi =27\text{°}$。

图  4  动力模型及监测点

Figure  4.  Dynamic models and monitoring points

下载: 全尺寸图片 幻灯片

模型的边界条件采用静止边界和自由边界，其中自由边界为软件自动生成，参数特性为主要网格边界的材料参数。模型阻尼使用局部阻尼作用于整个地质体，阻尼比设置为0.1571。网格尺寸按照Kuhlemeyer等提出的公式进行计算，公式如下：

式中：$c_{\mathrm{s}}$——剪切波速；

G——剪切模量/（N·m⁻²）；

$f_{max}$——地震波最高频率。

按照规定，网格尺寸应当小于等于输入地震波的最高频率对应波长的1/10～1/8，通过计算本文的网格尺寸的长度取20 m，在满足计算精度的同时效率更高。

2.2   动力输入

在选择动力的输入源时，由于缺少1950年察隅地震波的详细记录，特选择离研究区位置较近且震级也比较大的M6.9米林地震波，该地震波由拉萨地震监测站测得。滤波和时程转换前，将速度记录数据进行傅里叶变换，得到功率幅值谱。水平向及竖直向高能频率范围均为0.1～0.5 Hz，此为水平向及竖直向主频。由于地震波频率分布范围处于较低水平，在过滤高频杂音的同时，也防治漏掉可能存在的高频震波，特截取0～5 Hz内频率。

滤波后进行时程转换，得到水平向的峰值加速度为7.5×10⁻⁵ m/s²，竖向加速度峰值为8.0×10⁻⁵ m/s²。为还原强震工况，特将源加速度记录进行放大处理，将水平向和竖直向峰值加速度均放大到0.4 g，见图5（a）（b）。水平向和竖直向峰值加速度均在震动开始后10 s左右。

图  5  地震动力加速度记录

Figure  5.  Seismic dynamic acceleration record

下载: 全尺寸图片 幻灯片

3.   则隆弄沟斜坡动力响应模拟

3.1   高陡斜坡动力稳定性

为模拟1950年则隆弄沟强震作用下的坡体稳定性，在FLAC^3D中经过静力平衡计算之后，清除静力阶段的位移和速度，继续进行动力荷载的计算。得到地震作用下该力学模型存在一条潜在滑面（剪切口在坡高4000 m位置处），该滑面位置与静力计算结果所得相似，滑面上覆冰川冰碛物，下部为花岗片麻岩（图6）。其塑性破坏区域呈弧形，滑坡后缘为拉张破坏区，前缘为挤压破坏区，底部为剪切破坏区域，纵观整个塑性区分布图，得出坡体高度区间4000～5800 m有下滑的趋势。由此可知，在发生如此地震的时候该研究区危险系数较高，滑坡可能性较大。

图  6  动力条件下斜坡塑性区分布图

Figure  6.  Distribution map of plastic zones in slopes under dynamic conditions

下载: 全尺寸图片 幻灯片

为进一步探究该坡体在地震荷载下是否失稳，特选取1个内部监测点（P124）和5个坡面监测点（S114、S126、S137、S147、S157），图7所示5个坡面监测点位移均为与坡内监测点（P124）位移做差所得。用坡面监测点的水平、竖直位移减去内部监测点的位移，即为残余变形量（图7）。通过计算，可得坡面监测点的残余位移均较大，而坡内监测点残余位移极小，说明滑动主要出现在坡体表面。斜坡表面的残余变形一直处于增加的状态，说明坡体表面一直处于滑动状态。其中水平向最大残余位移为55 m，竖直向最大残余位移为25 m。

图  7  坡面监测点残余位移时程曲线

Figure  7.  Residual displacement time history curve at slope surface monitoring points

下载: 全尺寸图片 幻灯片

为查明该滑坡的失稳类型，同时进一步佐证滑面所处位置。现分析输入动载荷之后的最大切应变增量云图（图8），从图8中可以得到：①当输入米林地震波后，斜坡的最大剪应变增量为0.49，位置位于岩样1和岩样2之间，也就是花岗片麻岩和冰川冰碛物之间，坡高约5100 m处；②结合前文塑性区分布图以及抗拉强度为0的区域，可以判断米林地震波输入之后该潜在滑面发生了拉-剪破坏，从而失稳。

图  8  斜坡最大剪应变增量云图

Figure  8.  Contour map of maximum shear strain increment in the slope

下载: 全尺寸图片 幻灯片

对于失稳形态，在海拔4000 m处出现挤压破坏区域，失稳区间集中于海拔4000～6000 m之间，而并非在更低处的冰碛物出现挤压破坏，且最大剪应变增量位于海拔5500 m冰川冰碛物与花岗片麻岩的交界处，想要弄明白造成这种失稳现象的原因，将从坡体固有频率展开探究，同时，在这一框架下将进一步探讨固有频率对加速度放大效应的影响以及最高点的加速度放大系数。

3.2   边坡固有频率

地震发生时，在一些震级和烈度均较低的区域却发生了让人难以理解的边坡或者房屋的破坏情况。这种类型的破坏，多和震坏对象的固有频率有关。由此可知边坡的固有频率对于研究地震震害具有非常重要的意义。

依据前人研究的成果可知，边坡的弹性模量与固有频率呈正相关，且弹性模量以1 e¹¹Pa为界，之前增幅较缓，之后则增加较快。岩土体的密度越大，边坡的固有频率越小（密度范围为2000～3500 kg/m³）。几何特性方面，在其他参数不变的情况下，随着坡高的增加，边坡的固有频率将减小^[9]。边坡的固有频率可以用以下公式进行计算：

式中：$f$——边坡固有频率/Hz；

$h$——坡高/m；

$\rho$——材料密度/（kg·m⁻³）；

$E$——弹性模量/Pa；

e——自然指数。

固有频率与坡度和泊松比呈正相关关系，但其对固有频率变化影响的幅度比较低。通过式（3）进行计算，可以得到该边坡的固有频率约为0.09～1.89 Hz。

3.3   边坡共振

坡体是一个自由度比较高的系统，其本身的固有频率非常多，当边坡受到与自身固有频率相近的输入波频率激励时，坡体很可能会发生较强的振动，这种振动将会使坡体出现破坏^[9]。为进一步确定和缩小边坡可能产生共振的频率范围，通过在力学模型中最可能发生滑动的区域安插坡面监测点S114、S126、S137、S147和S157，共5个监测点，用于监测在不同频率震波作用下的位移情况。

通常情况下地震波包含了多种频率，比如本次模拟所用的米林地震波，其本身频率范围为0.01～5.0 Hz。为进一步确定处于高位的冰川冰碛物固有频率范围，特取0.01～5.0 Hz频率简谐波，从模型底部输入，对模型进行激励。结果表明：在坡体高度4400～5500 m，该范围内的固有频率处于0.2～0.6 Hz，与输入地震波主频存在一定的重合，说明该范围内将出现一定程度的同频共振（图9），图9内黄框为坡体地震动响应位移峰值，代表坡体对在黄框范围内的震动频率更敏感。

图  9  多频率激励下水平位移频响曲线

Figure  9.  Horizontal displacement frequency response curve under multi-frequency excitation

下载: 全尺寸图片 幻灯片

3.4   斜坡动力放大效应

引入概念Rg，边坡地震峰值响应加速度（peak resp-ond acceleration，PRA）与输入地震动峰值加速度（peak ground acceleration，PGA）的绝对值之比为Rg，见式（4）：

考虑到在实际工程环境中纯花岗岩地层存在的可能性，在研究带冰川冰碛物模型的同时，特添加一组均质的纯花岗片麻岩力学模型作为参照，用来验证坡体固有频率与材质的关系。一般来讲坡体放大效应与坡体材料存在相关性，这种相关性既体现在材料的强度上，同时也应体现在材料的固有频率上。在坡体材料发生共振时，材料对放大效应的影响程度如何。于是，在纯均质花岗岩力学模型中增加竖向监测点P124—P126，在带冰川冰碛物的底部为花岗片麻岩的力学模型中同样使用监测点P124—P126，这三个监测点用于监测在动力作用下不同材料对动力响应效果。在体现坡高对固有频率影响，从而导致边坡同频共振的时候，再在纯均质花岗片麻岩力学模型和带冰川冰碛物的花岗片麻岩力学模型中新增监测点P114—P123，这10个监测点同样用于监测水平加速度。为了解在强震工况下南迦巴瓦峰顶部的动力响应情况，再选取一组监测点P54—P71，用于监测水平加速度响应。

在验证不同材料对固有频率影响，导致坡体共振，从而影响动力响应结果时。对比分析图10（a）（b），可以看到在发生共振的条件下其动力响应效应会和一般情况有一定的不同，当水平地震波到达冰川冰碛物和花岗片麻岩的交界处时Rg值开始降低，从固有频率上看，冰川冰碛物的固有频率在降低，其固有频率降低后，虽然仍处于震波主频中，但是动力响应变差了许多，这说明主频中高能部分对于冰川冰碛物的响应不如花岗片麻岩。

图  10  边坡固有频率与Rg随材质及坡高变化折线图

Figure  10.  Line chart of slope natural frequency and Rg with material and slope height

下载: 全尺寸图片 幻灯片

在分析坡高对于Rg的影响时，图10（a）（b）结果显示Rg总体趋势为逐渐递增。其Rg值随着坡高的增加而不断增加，这与输入地震波主频相关，材质不变的情况下，坡高越大其固有频率越接近输入地震波主频，同频共振现象逐渐显现，坡高与固有频率呈负相关。在分析完坡高和材料对放大效应的影响之后，可以从图10（c）中看到南迦巴瓦峰的峰顶对于米林地震波的动力放大效应是比较强烈的，虽然在坡体内部加速度响应忽高忽低，但是峰顶的最大Rg值已经达到了2.81。

对于以上放大效应归根结底是地形和材质引起的，材质对于固有频率来说起主要作用，同时本次模拟的对象有两个特点，一个是坡度大，另一个是坡高大，高度和坡度作用于坡体相当于拓宽了其固有频率的范围。其次是地震波的主频与该力学模型的固有频率区间吻合程度高，导致产生共振，从而使得整个坡体出现较大的放大效应。

4.   讨 论

（1）局限性主要体现在：第一，计算模型和真实地质体的差异，本文主要采用假三维的地质剖面计算变形破坏特征，且地层均采用均值的岩体参数，本身和野外地质体存在一定差异；第二，使用米林地震波来研究1950年则隆弄动力响应特征，这里理应使用1950年察隅地震波输入分析，但是缺少资料，导致只能折中采用放大后的米林地震波；第三，使用有限差分计算如此庞大的物理模型，以前从未有人做过，因此其动力响应特征也未寻到相应参考，还需要更加深入的研究和探讨。

（2）本次工作的优势主要体现在：第一，在研究方法上选择了有限差分法，这种方法的计算效率和准确性都比较高，在参数赋值和调整上占优势；第二，虽然只做了假三维剖面的计算分析，但为真三维的计算做了铺垫，下一步工作将着手于真三维的地震动响应特征；第三，使用数值模拟能相当程度的节省科研经费，为之后可能进行的物理实验留足选择空间。

5.   结 论

（1）则隆弄沟在地形上属于高山峡谷地形，主要岩性为上覆冰碛物的花岗片麻岩，本文为建立1950年则隆弄沟地质剖面，通过收集历史灾害记录，得到冰舌位置为海拔3650 m^[26]，由多期次冰舌厚度推断得出冰舌厚度为200～300 m。由此绘制了1950年地质剖面图，并搭建了假三维物理计算模型。

（2）该坡体有较强的加速度放大效应，本文通过引入Rg对加速度放大效应进行表征。花岗片麻岩放大效应强于冰川冰碛物，坡顶放大系数为2.81。当输入放大后的米林加速度记录后，发现在海拔4000～6000 m的花岗片麻岩和冰碛物之间形成了一条贯通的滑面，失稳类型为拉-剪。

（3）该滑坡动力响应较复杂，滑面贯通破坏区域高，最大切应变增量在滑体中部介于滑体前缘和后缘之间（海拔约5000 m）达到最大，不排除存在共振可能。根据公式法算出该剖面的固有频率区间为0.09～1.89 Hz。将M6.9米林地震进行傅里叶变换，发现地震波主频与坡体的固有频率范围非常接近，水平向和竖直向主频均为0.1～0.5 Hz。通过多组不同频率的简谐波输入，发现海拔4400～5500 m处固有频率集中在0.2～0.6 Hz，与米林地震波主频存在高度重合关系，此结果验证了坡体确实存在同频共振现象。