三峡库区三门洞滑坡潜在涌浪风险研究*

霍志涛，黄波林，王鲁琦，卢书强，范意民

(1.中国地质调查局武汉地质调查中心，湖北 武汉433002；2.防灾减灾湖北省重点实验室，三峡大学，湖北 宜昌443002；3.中国地质大学(武汉)，湖北 武汉433000)

对于涉水滑坡，其失稳所引发的涌浪次生灾害将进一步扩大其影响范围，对航道过往的船只以及上下游沿岸居民的生命安全造成极大的威胁[1-2]。例如，1963年意大利瓦伊昂水库滑坡[3]和1985年长江新滩滑坡[4]。瓦伊昂水库产生的涌浪高出坝顶140 m多，3 000多人死亡，新滩滑坡产生涌浪高达54 m，涌浪范围涉及上下游42 km。2003年三峡库区千将坪滑坡滑移之后产生了39 m高的涌浪，造成岸坡上13人死亡[5]；2008年龚家方滑坡产生的涌浪在传播至4.2 km外的巫山新城时的涌浪高度仍然超过1 m[6-7]；2015年红岩子滑坡产生的涌浪导致多艘船只翻沉及2人死亡[8]。因而，预测滑坡涌浪的演化过程并确定其威胁范围，对于涌浪次生灾害的防控极为关键[9-10]。

三峡库区三门洞滑坡总体积为548万m3，属大型滑坡，一旦滑移破坏，所产生的涌浪将威胁青干河航运安全，确定其滑坡涌浪威胁范围及风险控制具有重要的意义。本文在分析三门洞滑坡破坏机理的基础上确定了滑坡的运动特征，进而采用库区崩塌滑坡涌浪灾害的快速评价系统软件(FAST，Fast assessing system for Tsunami)进行涌浪计算，该软件是基于四阶Boussinesq水波动力学模型二次开发得到的，在滑坡涌浪分析方面取得了很好的应用效果[11-12]。在完成滑坡涌浪的模拟之后，根据涌浪波幅的分布对航道及沿江斜坡的风险预警区进行了划分，为三峡库区类似滑坡涌浪灾害的防控提供依据。

1 滑坡概况

三门洞滑坡位于三峡库区秭归县青干河右岸边坡上(110°34′52.4″E，30°58′17.9″N)，距河口8.50 km，距三峡坝址50.2 km(图1)。

图1 三门洞滑坡位置及概貌

滑坡体平面形态为舌形状，顺向坡，滑坡前缘较缓，中部和后缘稍陡，西高东低(图1)。滑体后缘高程350 m，前缘高程125 m，后缘呈圆弧状，以基岩为界，右侧以陡坎临空面为界，左侧以基岩山脊为界，平均坡度15°，均宽300 m，滑体长830 m，滑体面积24.9万m2，滑坡体积548万m3，主滑方向60°。滑坡一旦成灾，将中断通往乡镇的公路，滑坡产生的涌浪将威胁青干河航运安全。

滑体由崩坡积物构成，物质成分为碎块石土，物质结构松散，厚度约20 m。滑带为土体与基岩接触带。滑床为侏罗系中下统紫红或浅灰色泥岩和泥质砂岩，岩层产状110°∠32°，为斜顺向坡，泥岩和砂岩强度较低，为易滑地层，属微透水岩体，利于上覆坡积物滑动。长江三峡库区支流青干河从滑坡前缘流过，水位在145～175 m间变化，淹没滑坡体前缘20～50 m，河水对滑坡坡脚进行冲蚀及掏空，库水位的升降导致库水与滑坡地下水互为补给。

2 滑坡变形特征

为了实时监控三门洞滑坡的变形趋势，在滑坡体上布设有1条监测纵剖面(剖面1-1′)，该剖面方向与滑坡主滑方向一致。滑坡于2006年10月开始实施专业监测，滑坡体上共布设7个GPS监测点，监测点布置详见图2，历年的监测数据如图3所示。

图2 三门洞滑坡专业监测布置图

图3 现场监测数据

由图3可知，自专业监测以来，滑坡持续变形，且具有阶跃型增长的变形演化特征。累计位移曲线在每年的5-9月出现阶跃上扬，之后趋于平稳，这种现象在2007年和2012年的变形中尤为明显。截止2015年7月，监测点的累计位移在600～2 300 mm之间，其中，位于滑坡前缘的两个监测点ZG365和ZG362变形最大，中部变形次之，后缘最小，整体呈现前缘快、后缘慢、指向青干河的牵引式运动特征。

由于滑坡坡面为松散状碎石土，降雨极易入渗至滑坡浅表土体，造成浅表土体的增重及局部变形坍塌；另一方面，由于地下水位上升及水的润滑、软化作用，对滑体物质和滑带土产生软化效应，直接降低了滑带土的抗剪强度；再者地下水的动水压力效应将增加沿渗流方向的滑动力，不利于滑坡体的局部稳定。此外，三峡水库的高水位蓄水(最高水175m)以及每年30m的大幅度水位升降成为常态，滑坡的涉水部分在江水的浸泡和冲刷作用下强度降低，会诱发滑坡前缘产生局部塌岸，对滑坡的稳定性造成一定影响。

3 滑坡失稳运动特征

通过滑坡宏观地质巡查和地表GPS监测可知，三门洞滑坡的变形特征为：滑坡由前缘到后部变形逐渐减小，具牵引式滑坡特征。水库蓄水时滑坡趋稳，退水和强降雨时，滑坡变形加剧，判断滑坡目前处于等速变形阶段，属于欠稳定状态。在降雨或水库快速降蓄水等诱发因素的影响下，该滑坡存在大规模滑移的可能性。由于三门洞滑坡体积较大且是顺向斜坡上覆的堆积层滑坡，可能出现较快速的滑动。

为了确定滑坡涌浪的规模和范围，需确定滑坡的失稳运动特征。目前有3种方法可以用于预测滑坡的运动速度，分别为工程地质类比法、重心质点的能量法估算以及美国土木工程师协会推荐法。定量估算的基础是假定滑动后滑体入水，采用等体积法(亦即滑动后的体积略大于滑动前的体积)来预估滑动后的位置，得到滑动前和滑动后重心位置，估算重心滑移距离和高差。需要注意的是，本文所述的滑坡速度为最差工况下的最大滑坡速度。

(1)工程地质类比法

从岩性和地形坡度上来看，三门洞滑坡与秭归马家坝滑坡岩性相同，两个滑坡同样发育于侏罗系中统聂家山组，且均为粉砂岩夹泥岩(图4)。但是，由于马家坝滑坡的坡度(25°)稍陡于三门洞滑坡(15°)，且马家坝滑坡体积(3000万m3)大于三门洞滑坡(548万m3)，三门洞滑坡的滑移速度应与马家坝滑坡相近或稍小。马家坝滑坡各区的运动速度不一，最大滑动速度约为2.5 m/s[13]，则三门洞滑坡最大运动速度也应约为2.5 m/s。

图4 对比剖面图

千将坪滑坡距三门洞滑坡约2 km，该滑坡是一个基岩滑坡，最大滑动速度为16 m/s。基岩滑坡与堆积层滑坡的类比性要差一些。因此，三门洞滑坡在发生大规模滑动时的最大速度应在2.5～16 m/s范围内。

(2)质点的能量法估算

能量法是根据牛顿运动定律提出的一种简单的物理计算方法，它概化了崩滑体的运动方式，忽略了岩土体运动过程的复杂性。该方法的来源清晰，物理意义明显，操作性强。对确定的摩擦系数而言，直线性滑动面的运动速度vs为[14]：

(1)

式中：g为重力加速度；ΔZsc为崩滑体的重心与计算点间的高度，即重心下降高度重心下落高度(m)；α为斜坡坡角(°)；f为运动摩擦系数。

根据该公式(1)，马家坝滑坡的vs为2 m/s，α为平均坡角25°，ΔZsc落差为50 m左右，求得马家坝滑坡的摩擦系数为0.466。由于三门洞平均坡度仅15°，如果将马家坝滑坡的摩擦力系数引入三门洞滑坡，滑坡无法启动。按照坡度对摩擦系数进行线性折减，可以得到三门洞滑坡滑动时的摩擦力系数为0.264。取三门洞主剖面为计算剖面，主剖面滑体的重心位置(高程230 m)为计算起点，ΔZsc落差为30 m，滑坡坡角约为15°。根据公式(1)计算可得三门洞滑坡整体运动速度约为2.9 m/s，该计算值与工程地质类比法所得结论基本相符，略大于2.5 m/s。

(3)美国土木工程师协会推荐法

美国土木工程师协会推荐滑动速度的计算公式为[14]：

(2)

式中：α为滑面倾角；W为滑体单宽重量；φ、c为滑动时滑面抗剪强度参数；H为滑体重心距离水面的位置；l为滑块与滑面接触面长(沿滑动方向)。饱和滑动土的力学参数φ=16.2°，c=12.6 kPa及滑动平均角度α=15°，重心滑动距离约107 m，当库区水位为175 m和145 m时，滑体运动速度分别为3.62 m/s和4.51 m/s。该计算结果与质点能量法及工程地质类比法分析的结果类似。

从上述计算和工程地质类比来看，速度值分布较为集中。因此，三门洞最大的运动速度约为2.9～4.51m/s。

4 滑坡涌浪数值模拟

(1)数值模型建立

采用FAST中的水库陆地滑坡涌浪源模型进行计算分析，其中，涌浪计算区的地形数据基于NASA开源的矢量化地形图构建(图5)。计算域东西长约18 616 m，南北长约11 856 m，采用26 m×26 m的栅格将计算域划分为716列及456行。计算区域从青干河上游杨家河至青干河河口及谭家河吒溪河的长江主航道段，包括有青干河、锣鼓洞河、泄滩河和吒溪河等支流。长江干流区域为宽谷地貌；支流多为峡谷地貌或不对称V型地貌，不利于涌浪能量消散。

图5 计算区地貌三维示意图

在进行涌浪源计算时，对库水位175 m和库水位145 m两种工况进行分析，不同工况条件下的输入参数如表1所示。

表1 滑坡涌浪源输入参数

计算每时步代表0.239 s，设置计算3 000时步，约717 s(约12 min)。按照波速17.7 m/s的速度计算，涌浪上下游均可以传播约12.7 km，其传播距离未超过计算区域，因此所设置的时步数可以满足河道涌浪的分析。

(2)滑坡涌浪灾害分析

由于库区水位145 m时，其入水速度更大，本文以145 m水位时涌浪计算结果为例，对涌浪灾害进行分析，得到滑坡涌浪传播瞬时河面如(图6)。涌浪模拟表明，在三门洞滑坡完全入江时，滑坡涌浪高度达到最大值，且已经抵达河岸，即滑坡涌浪在入江区域没有形成与爬高的清晰界线，与千将坪滑坡涌浪非常类似。T=11.95 s时，滑坡体完全入水，已有该点的最大涌浪抵达河岸(约21 m)；0.2 s之后，计算区最大的爬高在滑坡对岸的冲沟内形成(约25 m)。此后，波的主要传播方向由最初的滑坡滑移方向转化为沿河道方向，使得河道的纵向流动性大大加强。

图6 滑坡涌浪传播瞬时河面图(145 m水位)

由于涌浪传播区的河道蜿蜒延伸，涌浪在河道中传播时推进方向多变，传播浪衰减强烈。T=119.5 s时，最大涌浪传播至沙镇溪镇后崖邓家湾处，最大浪高约9.0 m，最大爬高约14 m。T=167.3 s时，最大涌浪传播至千将坪滑坡处，最大浪高约5.2 m，最大爬高约7 m。当T=370.5 s时，最大涌浪传播至长江对岸小幺姑沱，最大浪高约0.7 m，最大爬高约1.5 m。当T=561.5 s时，最大涌浪传播至吒溪河河口，最大浪高约0.3 m。当T=705.1 s时，最大涌浪传播至吒溪河内水田坝乡，最大浪高约0.1 m。

可以看出，虽然滑坡发生在青干河内，但长江干流、锣鼓洞河和泄滩河支流都受到了强烈影响。涌浪的分布区与传统认为的远端浪逐渐变小，有些差异。其原因是河流上游一般为峡谷或水很浅，导致了涌浪的雍高。

通过涌浪最大波幅分布图(图7)和不同流域区段最大浪高图(图8)可知，在主航道上(横剖面方向)，可分为急剧衰减区和平缓衰减区[15]。由于三门洞滑坡处于支流内，河流转弯区域非常多，滑坡产生的涌浪朝向下游，导致上游涌浪衰减慢于下游。对于上游的急剧衰减区，在河道泓深线附近200 m河道范围内涌浪从约15 m下降至8.7 m，平均100 m内涌浪下降高度约为3.15 m。在下游则呈断崖式下降，在河道泓深线附近100 m河道范围内涌浪从约15 m下降至6 m，平均100 m内涌浪下降高度约为9 m。因此，在该滑坡案例中，急剧衰减区仅有300 m长，该区域是涌浪危害航道的重点区域。由于河流蜿蜒，衰减机制不一，平缓衰减区仅近似满足缓斜线形式的下降规律，平均100 m内涌浪下降高度为0.01～0.3 m不等；上下游的下降率相近。沿程河道中的波高呈现复杂的非线性的下降趋势，且传播区支流两岸的浪高也各不相同。因此在平缓衰减区域，浪高一般呈波动状递减，该区域属于滑坡涌浪危害的延伸区域。

图7 流域内最大波幅分布图(145 m水位)

图8 不同流域区段最大浪高图 (145 m 水位)

图9 不同工况涌浪预警分区

表2 三门洞滑坡涌浪预警分区

图8展示了支流河流内涌浪雍高的情况和河流两岸浪高不对称情况。在河道纵深方向，除急剧衰减区外，涌浪传播过程中深水区的浪高低于浅水区，浅水区岸坡波浪壅高和爬高增加。在地形急剧变窄的区域，涌浪高度明显提高，出现放大效应。如支流上游峡谷区或凹槽地形区，涌浪曲线显现出雍高爬升的现象。在河道宽缓区域，涌浪则明显衰减，例如下游河道和长江河道。另外，由图8可知，支流内滑坡涌浪情况复杂多变且浪高较大，譬如锣鼓洞河和泄滩河等涌浪波幅大于0.5 m，锣鼓洞河支流内的涌浪波幅甚至大于3 m。这些支流虽然远离滑坡体，但是涌浪的影响却很大，需要注意相关防护。

(3)滑坡涌浪灾害风险分析

根据对三门洞滑坡145 m水位及175 m水位时的涌浪模拟计算的结果(采用模拟时段内各点的最大波高值)进行航道及沿江斜坡带风险预警区划分(图9、表2)。由于国内尚缺乏内河航道关于涌浪航运预警的方法，因此借鉴国家海洋局发布的《风暴潮、海浪、海啸和海冰灾害应急预案》[16]对本滑坡上下游影响范围进行涌浪风险预警分区。根据《应急预案》将浪高大于3 m时为航道划为红色预警区，当浪高在2～3 m时为航道划为橙色预警区，当浪高在1～2 m时为航道划为黄色预警区，当浪高小于1 m时划为蓝色预警区。

5 结论

根据现场调查及长期监测可知，滑坡目前处于等速变形阶段，属于欠稳定状态，在降雨或水库周期性蓄水等诱发因素的影响下，存在大规模滑移的可能性。针对三门洞滑坡可能产生的涌浪灾害进行了预测分析，并得到了以下结论及建议：

(1)采用工程地质类比法、质点能量法和美国土木工程师协会推荐法对三门洞滑坡入水最大速度进行了估算，得到了各水位条件下滑体入水的最大速度。库水位145 m时，滑坡的最大滑动速度约4.51 m/s；库水位175 m时，最大滑动速度约3.62 m/s。

(2)采用水波动力学数值模拟方法，利用浅水区滑坡涌浪源模型对长约18 km，宽约11 km的河道区域进行了计算预测分析。库水位145 m时产生的最大涌浪高度为21.0 m，对岸最大爬高为25 m；库水位175 m时产生的最大涌浪高度为20.6 m，对岸最大爬高为22 m。

(3)借鉴海啸预警划分方法开展了三门洞滑坡涌浪灾害风险分析。库区145 m水位时，滑坡滑移之后引发的涌浪灾害更为危险，且威胁范围也更广，此时，三门洞滑坡上游1.1 km、下游3.5 km为红色预警区，上游1.6 km、下游及锣鼓洞河部分区域长约1.1 km为橙色预警区，上游3.3 km、下游及锣鼓洞河6.9 km为黄色预警区，这些区域均具有危害风险。