李 仁 江,魏 明 垒,赵 小 铭,董 星 辰,黄 波 林

(1.中国三峡集团公司,四川 成都 610000; 2.中国水利水电第七工程局有限公司,四川 成都 610081; 3.防灾减灾湖北省重点实验室,湖北 宜昌 443002; 4.三峡大学 土木与建筑学院,湖北 宜昌 443002)

0 引 言

大型水库库区中往往存在着大型甚至巨型滑坡体[1]。水电站运行后,库水位将在较大的幅度范围内周期性升降,在此过程中,滑坡会随着库水位的变化而逐渐变形,甚至发生破坏[2-4],如三峡库区中的千将坪滑坡[2]、青石滑坡[4]以及水布垭库区的新塘滑坡[5]等。

滑坡变形在地下水渗流作用下会表现出一定的滞后效应,这对滑坡的变形预测造成了较大的影响。例如,石爱红等[6]对滑坡0.9 a内的监测数据进行了分析,认为在单个蓄水周期内,滑坡变形滞后于库水位5 d;陈亮青等[7]对滑坡2.4 a间的监测数据进行了分析,认为滑坡变形滞后于库水位升降1～2 d;高晨曦等[8]则得到滑坡变形滞后于库水位升降0.2～1.4 d的结论。值得注意的是,杨何等[10]在使用不同时段的监测数据对同一个滑坡进行滞后性分析时发现,滑坡变形滞后于库水位的滞后天数在不同的时间段内是并不相同的,这表明滑坡变形的滞后性并不是固定不变,而是会随着时间而改变,这会导致在预测时出现一定的误差。

在对滑坡的位移进行预测时,目前的预测方法虽然考虑了滑坡变形滞后性的影响,但均未考虑到滑坡变形的滞后性随时间发生了改变。为解决滑坡变形滞后性随时间改变而导致的预测困难问题,本文以溪洛渡库区雨林二组滑坡为例,结合长达9 a的库水位和滑坡位移监测数据,对滑坡的变形机制、滑坡变形滞后性随时间的变化规律以及考虑滞后性的滑坡位移预测方法进行了研究。基于研究成果,提出了一种考虑滞后性变化的滑坡变形预测方法,并以金沙江白鹤滩库区中的王家山滑坡为例[9],对预测方法的普适性进行了验证,可为库区其他涉水滑坡的变形预测提供借鉴。

1 溪洛渡库区雨林二组滑坡概况

如图1所示,溪洛渡库区雨林二组滑坡位于金沙江溪洛渡库区右岸,永善县务基镇捏池村,距坝38.7～39.1 km,且处于坝下游水库影响区域内。滑坡平面上呈矩形,纵向上呈陡缓陡的折线形态,上部和下部坡度在30°～35°之间,中部为约15°、长近200 m的缓坡平台。滑坡顺河长近700 m,后缘高程 790 m,前缘处于库水之下。540 m高程附近滑坡宽约450 m,最大厚度约100 m,主滑方向近于垂直河道。

图1 雨林二组滑坡地理位置Fig.1 Geographical location of the Yulinerzu landslide

雨林二组滑坡地质剖面图如图2所示。滑坡区出露奥陶系(O3t)和志留系(S1l)地层,岩性为灰色、深灰色、灰黄色、黄绿色砂岩、粉砂岩、泥岩、页岩、砂质页岩夹泥灰岩、泥质灰岩;滑坡下游为寒武系上统二道水组(∈3e)的灰-深灰色、粉-细晶白云岩、白云质灰岩夹少量砂岩、粉砂岩,基岩出露较多,多形成陡坎地形,地层产状为N15°～20°W/SW∠60°～70°,构成顺斜向岸坡。

图2 雨林二组滑坡地质剖面图Fig.2 Geological section of the landslides

在天然状态下,滑坡堆积体基本稳定,未见变形迹象。溪洛渡水库水位抬升后,堆积体前缘部分土体及滑带处于库水位以下,受水的浸泡影响,土体强度指标降低,加之堆积体前缘坡度较陡,在库水作用下,堆积体前缘变形塌滑,进而诱发了滑坡堆积体的整体复活变形。如图3所示,雨林二组滑坡在蓄水后的2013年6月19日前后出现变形迹象。至同年6月26日,后缘拉裂宽30～60 cm,最大约1.2 m,深度大于5 m(见图3(a),位于图4中蓝色虚框位置),错台高20～50 cm,最大超过1 m(见图3(b),位于图4中紫色虚框位置)。后缘和前缘次级裂缝较多,后缘次级裂缝长度20～50 m,宽度5～20 cm,局部有5～10 cm的错台;前缘次级裂缝相对较小,长10～20 m,宽0.5～3.0 cm,其中临江部位已出现数处垮塌。堆积体组成物质主要为志留系的泥页岩,后缘高程790 m,前缘高程约为500 m。水库蓄水前,河流水位约为410 m,在天然状态下堆积体基本稳定,未有变形迹象。

图3 滑坡后缘裂缝及错台Fig.3 Cracks and dislocations in trailing edge of the landslide

图4 位移监测点布置Fig.4 Position of displacement monitoring points

2 滑坡变形情况分析

为保障水库的安全运行,在图4所示位置设立了9个位移监测点,分别于滑坡前缘、中部和后缘,其中TP02、TP05、TP06和TP09位于滑坡前缘,TP01、TP04以及TP08位于滑坡中部,而TP03、TP07位于滑坡后缘。从2013年9月2日起,对滑坡的位移情况进行了长达9 a(库水位共升降变化9个周期,1,2,3,…,9)的监测。在9 a里监测记录了各监测点位在水平面上的位移L和重力方向的位移H,以及库水位变动情况和降雨情况,以监测点TP09(滑坡前缘)、TP01、TP04(滑坡中部)和TP03(滑坡后缘)为例,绘制出如图5所示的位移曲线。

图5 库水位、降雨及监测点位移曲线Fig.5 Reservoir water level,rainfall and monitoring point displacement curves

为方便对位移情况进行描述,对各监测点在水平面和重力方向上的位移进行编号,字母代表位移方向,而数字代表监测点编号,如H03表示监测点TP03在重力方向上的位移。

从位移情况来看,滑坡在水平面上产生最大位移的点为位于滑坡前缘的TP09,达到14 847.6 mm;在重力方向上的最大位移点为位于滑坡后缘的TP03,达到11 987.4 mm。滑坡前缘在水平面上的位移明显大于滑坡中部和滑坡后缘,而滑坡后缘在重力方向上的位移远大于滑坡前缘和滑坡中部。即从滑坡前缘到滑坡后缘,水平面上的位移逐渐减小,而重力方向上的位移逐渐增加。

滑坡前缘于2016年6月4日最先进入到相对稳定的阶段(历经1 003 d),随后滑坡中部和后缘几乎同时进入到相对稳定的阶段(滞后于前缘约301 d);滑坡整体最先再次出现小变形迹象(稳定变形约1 430 d)。

从库水位升降情况来看,除第1蓄水周期外,其余蓄水周期的最高水位均为600 m,每个周期持续时间基本相同,均值为372 d。每个周期库水位上升速度和库水位下降速度也基本相等,均值分别为0.41 m/d和-0.21 m/d。

从降雨情况来看,2015～2018年间,每年的降雨情况基本相同,而从2019年开始,降雨量稍有增加。最大降雨量出现在2020年的9月7日,达到407 mm/d。

3 雨林二组滑坡变形机制分析

由地质剖面图(见图2)可知,雨林二组滑坡剖面中部平缓,而上下部较陡,是典型的“靠椅状”滑坡。前缘消落带处于滑坡剖面线的陡缓交界,且地势较陡,给滑坡提供了较好的变形条件。“靠椅状”土质滑坡中前部的平缓段对滑坡体起到阻滑作用,而滑坡体物质结构稍密,透水性较差。

库水位在540～600 m之间周期性涨落,库水位升降直接引起了滑坡地下水渗透压力的变化。雨林二组滑坡的前缘坡脚较陡段大部分处于600 m以下,且透水性较差,导致地下水位升降速率小于库水位升降速率。在库水位下降过程中,地下水来不及排出坡体,滞后于库水位下降,从而形成指向坡体外部的水头差,产生渗透水压力效应,形成向坡体外的动水压力作用,使滑坡产生了明显变形。在库水位上升期间,因雨林二组滑坡的前缘坡体结构较密,渗透性较差,地下水上升滞后于库水位上升,形成向坡体内的水头差,产生的渗透水压力反推坡体,有利于滑坡稳定。在蓄水初期,每年水位上升时,滑坡产生的变形小于水位下降时滑坡产生的变形。由此可见,雨林二组滑坡是较为典型的动水压力型滑坡。

据历年降雨和变形数据可知,在约9 a时间内,降雨对滑坡变形的影响较小,具体表现为:降雨前后,滑坡未表现出明显的变形迹象,如在2020年9月7日,降雨量一度达到407 mm/d,各监测点位的位移曲线均未出现明显变化。因此,认为该滑坡体变形主要与库水位升降有关,而降雨对其影响较小。

4 考虑滞后性影响的滑坡变形预测方法

4.1 滑坡变形的滞后性分析

众多研究表明,滑坡变形往往会滞后于库水位升降[6-10]。而滑坡变形的滞后性往往难以进行估量,原因在于滞后性可能随着时间推移而发生改变[11-12]。根据长时间的监测数据,将位移所对应的时间后移,再与库水位进行相关性分析,找到相关性最高的后移时间,即找到了库水位升降导致滑坡主要变形的时间段[13]。如式(1)所示,通过皮尔逊相关法可以对多个变量间的相关密切程度进行计算[14],相关程度通过相关系数r表现。

(1)

式中:σx表示变量X的标准差,σy表示变量Y的标准差,而σxy表示X和Y的协方差。

计算得到的r值介于-1和1之间,|r|越接近于1意味着两者的相关性越强,r的正负分别代表正相关和负相关。当|r|≥0.8时认为变量间高度相关;0.5≤|r|<0.8时认为两变量中度相关;当0.3≤|r|<0.5 时则认为两变量低相关;而当|r|<0.3时,则认为两者相关程度极低,基本不相关[13]。

结合监测曲线可知,位移曲线总是随着时间增长而单调递增的,而库水位则呈周期性变化有升降起伏的。水位上升时,库水位与滑坡变形将呈现正相关;水位下降时,库水位与滑坡变形呈负相关。由此,以TP03在重力方向的监测数据(H03)为例,9个周期内水位升、降时对应位移的时间后移ty=1～15 d,再对每个周期内的位移和库水位变化进行相关性分析,便可得到滑坡变形与库水位变化的相关性随滞后天数的变化规律。如图6所示,图中红框标识了H03与库水位变化相关性最高的滞后天数ty及相关系数r。

图6 不同蓄水周期内H03与库水位变化的相关性随滞后天数的变化规律Fig.6 Variation of correlation between H03 and reservoir water level with lag days in different periods

从相关性来看,H03与库水位变化始终保持高度相关关系,在整个变形过程中|r|几乎都为0.8以上。且随着时间的推移,库水位升高与滑坡变形的相关性会进一步提升(见图6(a)),|r|由T=2时的0.85上升至T=5时的0.96;而库水位下降与滑坡变形的相关性却恰好相反,随着时间的推移,相关性逐渐降低(见图6(b)),|r|由T=2时的0.94下降至T=5时的0.74。

从滞后性来看,滑坡变形的滞后性在水位上升和下降时的规律不同。在水位上升时,变形的滞后性会随着蓄水时间逐渐显现;而在水位下降时,几乎没有表现出滞后性。TP03在第2,3,4和第5个周期水位上升时,ty=0,4,9 d和11 d时的|r|最大。即在蓄水的第2周期,滑坡的主要变形和库水位上升几乎在当天发生,而在第4和第5周期,滑坡的主要变形分别滞后于库水位上升9 d和11 d。而在水位下降时,除第3周期以外,ty均为1 d时的|r|最大,即滑坡主要变形几乎在降低水位的当天内发生。

提取出各周期中库水位上升时相关性最高的滞后天数,发现库水位上升时的ty和蓄水周期序数T之间基本满足对数函数规律。通过式(2)进行拟合,得到滑坡变形滞后时间ty随蓄水周期数T的变化规律,如图7所示。

图7 升水位时ty随T的变化规律及拟合曲线Fig.7 The variation law and fitting curve of ty with T at rising water level

ty=alnT+b

(2)

式中:a和b均为拟合参数,与滑坡体条件有关。

4.2 雨林二组滑坡变形预测

通过监测数据对滑坡后续的变形速率进行预测,可以预测滑坡的后续变形[15]。考虑到滑坡变形的滞后性,应将滑坡位移后移相应天数ty后再进行分析:升水位时ty满足式(2),而降水位时ty=0。有研究表明,滑坡变形速率与水位升降速率有关[13]。除此之外,滑坡变形的滞后性还会随着蓄水周期而改变。因此,在进行变形预测时,也应考虑到滞后性的变化。在分析时,则应根据蓄水周期数的不同,采用不同的滞后天数进行计算分析。按照式(3),根据所处的蓄水周期,将滑坡变形后移相应的滞后天数后,再用每周期升、降水位时滑坡的变形速率除以每个周期升、降水位速率,便能得到每个蓄水周期中考虑滞后性影响的滑坡变形速率与水位升降速率关系。

(3)

式中:I为滑坡变形速率与水位升降速率之比;v为滑坡变形速率;vw为水位的升降速率,升高水位时为正,降低水位时为负;ΔH为时间间隔Δt内滑坡的变形量;Hw为同一时间间隔内水位的变化量。

对计算公式进行分析可知,在库水位升高时I值为正(记为Is),在库水位下降时I值为负(记为Ij)。|I| 较大则意味着在相等的水位升降速率时,引起的滑坡变形更快,因此,|I|在一定程度上也能够表示滑坡变形对库水位升降的敏感性。

由监测数据得知,库水位在每个蓄水周期的变化速率基本相同(升高水位速率约为0.41 m/d,降低水位速率约为-0.21 m/d)。仍以TP03的监测数据为例,对每个周期中H03的I值进行计算,发现其随T的变化规律基本服从指数函数规律(式(4)),进行拟合便可得到如图8所示的拟合曲线。

图8 升、降水位条件下I随T的拟合曲线Fig.8 Fitting curve of I with T under the condition of reservoir water level rise and decline

I=ce-fT+I0

(4)

式中:c和f均为拟合参数,与升降水位、滑坡体条件有关;c在升高水位时为正,降低水位时为负;I0表示滑坡发生变形后进入相对稳定阶段的变形速率与库水升降速率的比值。

从图8可以看出,不论水位升高还是降低,|I|均随着蓄水周期的增加而逐渐减小。在库水位升高条件下,|Is|由T=2时的0.0043下降至T=8时的0.00026;在库水位下降时,Ij由T=2时的-0.070下降至T=8时的-0.00083。在前4个蓄水周期中,|I|在降低水位时远大于升高水位时,T=2时,|Is|=0.00043,但|Ij|=0.070;T=4时,|Is|=0.0025,但|Ij|=0.014。

将式(3)和式(4)进行联立,即可根据vw和T求得滑坡变形速率(式(5))。

v=cvwe-fT+I0vw

(5)

再对滑坡变形速率与时间进行积分,便可得到相应的位移计算公式(式(6))和如图9所示的位移预测曲线。

图9 玉林二组滑坡变形预测曲线Fig.9 Deformation prediction curve of Yulinerzu landslide

(6)

式中:H0为上一次变形后滑坡的累计变形值。

从图9可以看出,考虑滞后性影响的滑坡变形预测曲线与监测结果基本吻合,且在第8和第9周期,预测变形与监测数据完全相同(图中蓝色虚线框所标识位置)。随着蓄水周期的增加,每次库水位升、降后,滑坡所产生的变形逐渐减小(a3a1,b2>a2,b3>a3)。

根据监测数据及滑坡位移预测曲线,将库水位滞后性引入到滑坡位移预测模型中是可行、有效的,且预测结果也较为准确。本文认为雨林二组滑坡在后续变形中能继续保持稳定,且不会产生较大变形。

从滞后天数随蓄水周期的变化规律可知,滑坡变形的滞后性变化是十分缓慢的:由第1个蓄水周期时的滞后0 d到第8个蓄水周期时的之后15 d,经历了约8 a的时间。

5 方法普适性验证

为验证方法的普适性,选取位于白鹤滩库区王家山滑坡[9]进行变形预测分析。王家山滑坡(见图10)在平面上呈近三角形。滑坡纵长约800 m,剪出口宽约600 m,厚度14.1～87.6 m,总体积约为611万m3。2021年4月起,白鹤滩库区开始阶段性蓄水后,王家山滑坡的变形加剧。

图10 王家山滑坡Fig.10 Wangjiashan landslide

对王家山滑坡在1个蓄水周期内的监测数据进行分析,由于监测时间较短,改用蓄水天数替代蓄水周期进行拟合预测。通过相关性分析得到滑坡变形滞后于库水位变动的滞后天数ty在升水位时为1 d,降水位时为0 d,且在此蓄水周期中基本保持不变。将滞后天数代入并进行计算,便能得到所需的拟合参数。由于白鹤滩库区是首次蓄水,水位的升降速率不恒定,按照式(5),根据水位升降速率是否发生变化,将库水位的升、降过程进行拆分,以获取不同库水位升降速度所对应的Is和Ij随时间的变化曲线。获取各段的拟合参数如表1所列,进而计算得到变形的预测曲线如图11所示。

表1 用于王家山滑坡变形预测的拟合参数Tab.1 Fitting parameters for deformation prediction of Wangjiashan landslide

图11 王家山滑坡变形预测曲线Fig.11 Deformation prediction curve of Wangjiashan

从图11可以看出,在1个蓄水周期中,滑坡变形滞后性变化不明显时,除初次蓄水的A段外,考虑滞后性影响的滑坡变形预测方法仍与监测结果基本吻合(B、C段)。预测结果表明,王家山滑坡在下一个蓄水周期时,仍会保有一定的变形速率继续变形,应加强监测,以防滑坡失稳破坏。

6 结 论

(1) 从监测数据得知,雨林二组滑坡在蓄水初期(1～4周期中)变形较大,且降低水位时滑坡产生的位移远大于升高水位时的滑坡位移,如第2个蓄水周期中,降低水位时产生的变形约为升高水位时的5倍。而后进入到相对稳定的变形阶段(5～9周期),且变形逐渐收敛。在蓄水初期,随着蓄水周期的增加,滑坡每个周期产生的变形量逐渐减小,由第2周期的4 005.7 mm下降至第5周期的78.5 mm。相应的变形速率也逐渐减慢,由第2周期的9.92 mm/d下降至第5周期的0.21 mm/d。

(2) 通过将滑坡变形对应的时间后移,再与库水位升降进行相关性分析,得到了滑坡变形滞后性的变化规律:滑坡变形的滞后性在水位升高的过程中表现明显,而在水位下降的过程中表现不明显;在水位上升时,滑坡变形的滞后性会随着时间逐渐显现,滞后天数由第2蓄水周期的1d增加至第5周期的11d,研究发现滞后天数基本满足对数函数规律。

(3) 在引入滞后性的影响后,对滑坡的变形速率进行拟合分析,得到了变形速率随时间的变化规律,在库水位升降速率相同的情况下,滑坡变形速率基本满足指数函数关系。并进一步预测了滑坡的后续位移。预测结果表明,雨林二组滑坡在后续变形中,变形速率会继续降低,变形会逐渐收敛,与实际监测数据相符。

(4) 通过对王家山滑坡的变形进行预测分析,验证了考虑滞后性变化预测方法的普适性。预测结果表明,王家山滑坡若不加以治理,则仍会以较大的变形速率持续变形。在后续应加强对王家山滑坡的监测,并对滑坡进行治理,以防滑坡失稳破坏。

(5) 本文没有考虑降雨因素对滑坡变形的影响,当降雨与库水位共同作用时,需额外考虑两因素所占权重以及对滑坡滞后性的影响,进而对滑坡变形的滞后性展开更深入的研究。