三峽庫區典型堆積層滑坡變形滯后時間效應研究*

高晨曦 劉藝梁 薛 欣 易慶林 陳健翔

(①三峽大學防災減災湖北省重點實驗室， 宜昌 443002， 中國)

(②三峽大學三峽庫區地質災害教育部重點實驗室， 宜昌 443002， 中國)

(③三峽大學湖北長江三峽滑坡國家野外科學觀測研究站， 宜昌 443002， 中國)

0 引 言

我國是個深受滑坡災害困擾的國家，每年由滑坡造成的經濟損失異常慘重(許強等， 2004)。其中：堆積層滑坡是三峽水庫運行過程中的重要地質災害，其變形演化往往滯后于庫水位的變化，表現出時間滯后效應，這給滑坡災害精準預測和災害警情準確發布造成極大困擾。堆積層滑坡變形滯后時間的研究對實現這類滑坡災害準確預報有重要作用，同時對降低滑坡災害造成的經濟損失及人員傷亡具有重要意義(亓星等， 2020)。繆海波等(2014)以三峽庫區滑坡為研究對象，通過滑體變形與滑體內部孔隙水壓力變化間相關性的分析，揭示了滑體內部超孔隙水壓力變化的外界因素觸發機制，探討了降雨和庫水位變化作為外界因素對水壓力變化的影響，該機制的揭示為解決滑坡位移預測模型考慮滯后效應提供了有益參考。楊巧佳等(2015)運用統計學方法，建立了滑坡變形滯后與庫水位及降雨的回歸模型，認為庫水位降速對變形滯后期的影響占主導作用。石愛紅等(2013)選取月降雨量與庫水響應滯后影響下的庫水位月平均變化量為誘發因子，采用時間序列分析方法進行分析，研究表明誘發因子對滑坡的位移變化具有重要的影響。張建等(2018)對雨水入滲規律進行分析，研究土體滲透參數與土體強度參數及遲滯現象的關聯性。陳亮青等(2018)以樹坪滑坡為例，研究了誘發因素對于滑坡位移變形的滯后影響。前人對研究庫水位與滑坡變形滯后響應之間關系做了大量研究，對于定性分析該方面問題的述評不夠全面，不能很好地闡述清楚當前研究中的不足。因此，本文在前人的基礎上，以三峽庫區典型堆積層滑坡——樹坪滑坡為例，以實時監測數據為依據，采用集對分析法并結合層次分析法，構建滑坡加權位移向量計算模型，在滑坡加權位移演化與庫水位波動相互關系定性分析的基礎上，尋找滑坡加權位移與庫水位下降速率相關性達到最大時的平移步數，計算出滑坡變形滯后具體時間，實現對三峽庫區樹坪滑坡變形滯后時間效應定量研究。研究成果可以為樹坪滑坡的監測預警防治工作提供參考，對重大水利工程涉水滑坡監測預警具有一定借鑒意義。

1 基于集對分析的時間滯后模型

1.1 集對分析法基本原理

集對分析法是趙克勤于1989年提出的一種解決確定-不確定問題的新方法，能較好地實現不確定信息的整體處理(趙克勤， 2000)。劉曉等(2009)將集對分析法引入巖土變形監測分析領域，開展了滑坡變形對水庫蓄水過程滯后響應的定量研究。

集對分析在一定的問題背景下，討論的兩個集合所具有的特征，進行同、異、反分析并加以定量刻畫，得出這兩個集合的同、異、反聯系度表達式(劉曉等， 2009； 桂蕾等， 2012)：

μ=a+bi+cj

(1)

式中：μ為聯系度；a為同一度；b為差異度；c為對立度；i為差異度系數；j為對立度系數。

1.2 構建滑坡加權位移向量計算模型

設有Z1、Z2、Z3、…、Zm共m個監測點，第k個監測點t時刻的位移記為xt, k(t=1, 2,…，n；k=1, 2…，m)，構成多點位移時序矩陣X：

(2)

考慮到每個滑坡位移監測點所處空間位置不同，在滑坡整體位移變形中所占權重大小不同，因此設監測點的權重向量W為：

(3)

則滑坡加權位移向量可表示為：

S=XWT

(4)

1.3 考慮庫水位下降速率的滑坡變形滯后時間計算模型

由于三峽庫區庫水位的不斷變化，使得庫區滑坡演化受控于庫水位的漲落過程。

為考察相鄰監測時段滑坡加權位移演化，引入某時間段內滑坡加權位移平均變化速率閾值q(單位為mm·d-1)，將滑坡加權位移變化按照“波動q1～q2內”、“大于q2”和“小于q1”，劃分為“平穩增長”、“加速增長”和“減速增長”3種狀態。考察庫水位下降期間的波動情況，引入某時間段內庫水位平均下降速率閾值p(單位為m·d-1)，將相鄰時段內庫水位的下降速率按照“波動p1～p2內”、“大于p2”和“小于p1”，劃分為“平穩下降”、“加速下降”和“減速下降”3種狀態。確定閾值q和p的波動范圍是進行相關性分析的關鍵。

通過定性分析滑坡加權位移變化與庫水位下降速率相互關系，建立同異反關系判別表，如表1 所示。

表1 滑坡加權位移變化與庫水位下降速率同異反關系判別

選取2012年、2013年、2014年汛雨期地表位移與庫水位下降速率的監測數據，構建滑坡加權位移S={St}(t=1, 2, 3,…,n){Xt}與庫水位波動L={Lt}(t=1, 2,…,n)時序，從t=2開始，可根據表1 進行n-1次判別，設其中取同、異、反的個數分別為na、nb、nc，其按式(1)聯系度μ的系數為：

(5)

聯系度μ刻畫了某段時間滑坡加權位移和庫水位下降速率的相關性，其a值越大，說明所構成的集對相關性越好。設平移步數為step，對滑坡加權位移進行等步長Δt平移，即考察t→t+1區間滑坡加權位移和庫水位下降速率的相關性。取step=0， 1， 2…，搜索不同step下a的值，當a達到最大時，其對應的step與Δt的乘積即為滑坡變形滯后時間T：

T=step·Δt

(6)

2 典型滑坡分析

2.1 樹坪滑坡概況

樹坪滑坡地處三峽庫區長江南岸，距三峽大壩約47km。樹坪滑坡屬于古崩滑堆積體，呈南北向展布，向北傾斜，發育于沙鎮溪背斜南翼，其巖層產狀為120°～173°∠9°～38°，滑坡體平均坡度為22°，屬逆向坡。巖性主要為泥巖、粉砂巖夾泥灰巖。樹坪滑坡總體呈圈椅狀，后緣高程380～400m，前緣高程約60m，南北長約800m，東西寬約700m，面積約55×104m2，厚約30～70m，平均厚約50m，總體積約2750×104m3(盧書強等， 2014)。

2012年6月，建設了實時相對位移自動專業監測點。根據樹坪滑坡多年變形特征，將其劃分為主滑區和影響區，主滑區滑體上現有6個位移監測點分別為SP-2、SP-6、ZG85、ZG86、ZG87和ZG88，變形影響區有2個位移監測點為ZG89和ZG90，監測點布置如圖1 所示。因監測點所處滑坡體的位置、高程和監測數據變化量等因素不同，本文選取ZG85、ZG86、ZG88和SP-2這4個監測點的滑坡變形數據作為數據來源，選取監測點的變形量如表2 所示。

表2 樹坪滑坡專業監測點變形分析

圖1 樹坪滑坡監測點平面布置圖

2.2 庫水位變化對樹坪滑坡變形的影響

由于三峽庫區庫水位不斷變化，使得庫區滑坡演化受控于庫水位漲落過程。從圖2 可知，自三峽庫區蓄水以來，當庫水位下降或低水位運行期(每年4月至9月)，滑坡累計位移曲線會出現一個明顯的變形增長階坎，此時各監測點位移變化幅度增大。而在庫水位上升或高水位運行階期(每年10月至次年3月)，滑坡累計位移曲線相對趨于平穩，即各監測點位移變化幅度減小。因此，庫水位下降與樹坪滑坡變形之間具有較強的相關性。

圖2 樹坪滑坡地表累積位移-庫水位-降雨量-時間關系

樹坪滑坡地表位移變形速率與庫水位升降速率間的相互關系圖中(圖3) 顯示，從2012年到2014年，當庫水位平均下降速度達到最大，而在此期間滑坡平均位移速率也達到了最大。這說明庫水位下降對樹坪滑坡變形影響較大，且庫水位下降速率越快，造成的滑坡變形越大(肖捷夫等， 2020)。

圖3 樹坪滑坡地表位移變形速率-庫水位升降速率-時間關系

3 典型滑坡變形滯后時間分析

3.1 樹坪滑坡變形滯后效應整體分析

由于各位移監測點所處于滑坡體空間位置不同，在滑坡整體位移變形中所占權重大小不同。因此，利用層次分析法，對ZG85、ZG86、ZG88和SP-2共4個監測點進行兩兩重要程度分析，并進行賦值，得到判別矩陣A為：

(7)

求解矩陣A的特征值，由此可知判別矩陣A的最大特征值為λmax=4.1256，λmax>n，為保證判定矩陣A準確性，為此需要檢驗其是否達到“滿意一致”的程度(鄧雪等， 2012)。

一致性指標CI：

(8)

一致性率CR：

(9)

式中：RI為隨機一致性指標，可按表3選取(鄧雪等， 2012)。

表3 隨機一致性指標RI

計算得CR=0.047<0.1，于是判別矩陣A通過了一致性檢驗。再采用幾何平均數估計各因素權重。

(10)

經計算可得權重向量為：

(11)

隨著三峽水庫蓄水過程的逐漸完成，樹坪滑坡體的安全儲備逐漸降低，并于2014年8月開始應急治理， 2015年6月全面竣工完成。為了計算結果更加直觀明顯，故選取治理之前的監測數據。因選取的3個時間段內總降雨量分別為10.1mm、141.4mm和64.2mm，降雨量較小，所以本文不考慮降雨這一影響因素。利用集對分析方法分別計算出僅在庫水位下降的情況下滑坡變形滯后時間，具體步驟如下：

(1)對樹坪滑坡2012年監測數據進行觀察分析，在8月1日至8月21日庫水位與滑坡位移相互關系較明顯，庫水位由160.25m下降至146.12m，平均下降速率為0.7m·d-1，即庫水位平均下降速率閾值p=0.7m·d-1，規定波動范圍為±0.1 m·d-1。因此，庫水位下降速率在0.6～0.8m·d-1之間為平穩下降； 庫水位下降速率大于0.8m·d-1為加速下降； 庫水位下降速率小于0.6m·d-1為減速下降。

為考察滑坡位移滯后響應，位移數據選取2012年8月1日至8月31日數據，步長為1d等間距位移時序共計31組數據，該時間段內滑坡位移變形平均速率為3mm·d-1，即平均變化速率閾值q=3mm·d-1，規定滑坡位移變形速率在2～4mm·d-1之間為平穩上漲； 滑坡位移變形速率大于4mm·d-1為加速上漲； 滑坡位移變形速率小于2mm·d-1為減速上漲。

按照表1 判別同異反關系，并搜索不同平移步數下的同一度a。將繪制出的同一度與平移步數關系曲線進行二次多項式擬合，如圖4 所示。以此判斷當step為2步時，同一度a達到最大值，此時滑坡位移與庫水位相關性達到最大。

圖4 2012年同一度與平移步數關系

根據庫水位下降對滑坡變形滯后時間公式T=step·Δt，計算滯后時間為2d。

(2)對樹坪滑坡2013年監測數據進行觀察分析，在5月10日至6月11日庫水位與滑坡位移相互關系較明顯，庫水位由160.17m下降至146.26m，平均下降速率0.43m·d-1，即庫水位平均下降速率閾值p=0.43m·d-1，規定波動范圍為±0.1m·d-1。因此，庫水位下降速率在0.33～0.53m·d-1之間為平穩下降； 庫水位下降速率大于0.53m·d-1為加速下降； 庫水位下降速率小于0.33m·d-1為減速下降。為考察滑坡位移滯后響應，位移數據選取2014年5月10日至6月20日數據，步長為1d的等間距位移時序共計42組數據，該時間段內滑坡位移變形平均速率為4.4mm·d-1，即平均變化速率閾值q=4.4mm·d-1，規定滑坡位移變形速率在3.4～5.4mm·d-1之間為平穩上漲； 滑坡位移變形速率大于5.4mm·d-1為加速上漲； 滑坡位移變形速率小于3.4mm·d-1為減速上漲。

按照表1 判別同異反關系，并搜索不同平移步數下的同一度a。將繪制出的同一度與平移步數關系曲線進行二次多項式擬合，如圖5 所示。以此判斷當step為5步時，同一度a達到最大值，此時滑坡位移與庫水位的相關性達到最大。

圖5 2013年同一度與平移步數關系

根據庫水位下降對滑坡變形滯后時間公式T=step·Δt，計算滯后時間為5d。

(3)對樹坪滑坡2014年監測數據進行觀察分析，在5月6日至6月4日庫水位與滑坡位移相互關系較明顯，庫水位由160.03m下降至147.34m，平均下降速率0.44m·d-1，即庫水位平均下降速率閾值p=0.44m·d-1，規定波動范圍為±0.1m·d-1。因此，庫水位下降速率在0.34～0.54m·d-1之間為平穩下降； 庫水位下降速率大于0.54m·d-1為加速下降； 庫水位下降速率小于0.34m·d-1為減速下降。為考察滑坡位移滯后響應，位移數據選取2014年5月6日至6月14日數據，步長為1d的等間距位移時序共計42組數據，該時間段內滑坡位移變形平均速率為4.2mm·d-1，即平均變化速率閾值q=4.2mm·d-1，規定滑坡位移變形速率在3.2～5.2mm·d-1之間為平穩上漲； 滑坡位移變形速率大于5.2mm·d-1為加速上漲； 滑坡位移變形速率小于3.2mm·d-1為減速上漲。

按照表1 判別同異反關系，并搜索不同平移步數下的同一度a。將繪制出的同一度與平移步數關系曲線進行二次多項式擬合，如圖6 所示。由圖像可發現，同一度擬合曲線出現兩個峰值，由此判斷，在該時間段內出現了兩次變形滯后響應，故需要分段考慮。又因第一次滯后響應受前面數據的影響，所以本文只考慮第二次滯后響應。因此，當同一度a第二次達到最大值時，相對平移步數為4，此時滑坡位移與庫水位的相關性達到最大。

圖6 2014年同一度與平移步數關系

根據庫水位下降對滑坡變形滯后時間公式T=step·Δt，計算滯后時間為4d。

將2012年、2013年、2014年計算得到的滑坡變形滯后時間與相應時期庫水位平均下降速率繪制成圖表(表4、圖7)。僅對樹坪滑坡，當庫水位下降速率小于等于0.43m·d-1時，滑坡變形滯后時間大于等于5d； 當庫水位下降速率在0.43m·d-1到0.7m·d-1之間時，滑坡變形滯后時間在2d到5d之間； 當庫水位下降速率大于等于0.7m·d-1時，滑坡變形滯后時間小于等于2d。因此，隨著庫水位下降速率不斷增大，滑坡變形滯后時間不斷縮短。造成這種現象的原因是：當庫水位下降，產生向外的滲流，由于坡體滲透性較差，地下水排出緩慢，導致坡體內孔隙水壓力增加，有效應力減小，降低了滑坡的穩定性(尚敏等， 2019； 湯明高等， 2019)。這與向玲等(2014)研究的結論“樹坪滑坡庫水位下降速率越快，地下水位上凸現象越明顯，穩定性系數減小越明顯”相符。因此隨著庫水位下降速率越大，滑坡變形滯后時間越短。

表4 庫水位下降速率與滑坡變形滯后時間對應關系

圖7 庫水位下降速率與滑坡變形滯后時間關系

為進一步驗證，基于集對分析法計算出的滑坡變形滯后時間結果的實用性和準確性，現將理論分析結果與實際監測結果進行對比分析。

選取2013年5月10日至6月20日樹坪滑坡自動監測點(ZG85、ZG86、ZG88和SP-2)的實際累積位移與庫水位數據，繪制滑坡自動監測點累積位移-庫水位-時間關系曲線如圖8 所示。分析圖像可知， 2013年5月13日至6月9日庫水位快速下降，而累積位移在5月17日到6月14日之間急劇增大。因此，滑坡變形滯后于庫水位的下降，本次滯后時間為5d左右，且與前文2013年相對應時間段內計算出的滑坡變形滯后時間基本吻合。

圖8 滑坡自動監測點累積位移-庫水位-時間關系

同樣選取2014年5月6日至6月15日樹坪滑坡自動監測點(ZG85、ZG86、ZG88和SP-2)的實際累積位移與庫水位數據，繪制滑坡自動監測點累積位移-庫水位-時間關系曲線如圖9 所示。分析圖像可知， 2014年5月6日至6月7日庫水位快速下降，而累積位移曲線在5月10日之后陡然上升，滑坡變形加劇，到6月11日之后滑坡變形監測曲線趨于平穩。本次滑坡變形滯后于庫水位下降的時間為4d左右，且與前文2014年相對應時間段內計算出的滑坡變形滯后時間基本吻合。

圖9 滑坡自動監測點累積位移-庫水位-時間關系

通過集對分析法計算滯后時間與監測數據圖像分析的滯后時間相互驗證，理論分析結果與實際監測結果吻合較好。若通過監測數據圖像分析滑坡的變形滯后時間，需選取出大量數據繪制曲線圖，才可觀察出整體變化趨勢，從而確定滑坡變形滯后時間。而集對分析法適用于短期滑坡變形滯后時間的計算，相對監測數據圖像分析更方便簡單。因此，運用集對分析法進行滑坡變形滯后時間的定量研究是可行的思路。

3.2 樹坪滑坡不同空間位置的滯后時間分析

為考察在同一滑坡體中，庫水位下降過程對滑坡體空間位置上變形滯后時間長短的影響，根據表5 選取4個監測點，分為主滑區中部監測點SP-2、ZG86和前緣監測點ZG88、ZG85。

表5 樹坪滑坡專業監測點所處位置

利用上述數據運用相同的方法，按照滑坡加權位移-庫水位同異反關系判別表，搜索不同平移步數下的同一度a。繪制同一度與平移步數關系曲線，并根據二次多項式進行擬合得到同一度擬合曲線。

由圖10 可知，庫水位以0.7m·d-1的速率下降，滑坡體前緣的變形滯后時間為2.4d，滑坡體中部的變形滯后時間為3.4d； 庫水位以0.43m·d-1的速率下降，滑坡體前緣的變形滯后時間為5.4d，滑坡體中部的變形滯后時間為5.6d； 庫水位以0.44m·d-1的速率下降，滑坡體前緣的變形滯后時間為4d，滑坡體中部的變形滯后時間為5.4d。

圖10 不同空間位置同一度與平移步數關系曲線

綜上可知，僅對樹坪滑坡而言，如圖11 所示，不同庫水位下降速率對滑坡體空間位置上變形滯后時間影響特點為：越靠近滑坡體前緣變形滯后時間越短，當庫水位下降速率在0.43m·d-1到0.7m·d-1之間時，滑坡前緣和中部的變形滯后時間差在0.2d到1.4d之間。

圖11 不同空間位置上庫水位下降速率與變形滯后時間關系

4 結 論

(1)通過滑坡加權位移的演化與庫水位波動相互關系定性分析的基礎上，尋找滑坡加權位移與庫水位下降速率相關性達到最大時的平移步數，從而計算出滑坡變形滯后時間，且結果與實際情況相符。

(2)當庫水位下降速率小于等于0.43m·d-1時，樹坪滑坡變形滯后時間大于等于5d； 當庫水位下降速率在0.43m·d-1到0.7m·d-1之間時，樹坪滑坡變形滯后時間在2d到5d之間； 當庫水位下降速率大于等于0.7m·d-1時，樹坪滑坡變形滯后時間小于等于2d。

(3)以樹坪滑坡為例，庫水位下降過程對滑坡體空間位置上變形滯后時間影響特點為：越靠近滑坡體前緣變形滯后時間越短，當庫水位下降速率在0.43m·d-1到0.7m·d-1之間時，滑坡前緣和中部的變形滯后時間差在0.2d到1.4d之間。

(4)本文沒有考慮降雨因素對滑坡變形響應的滯后性，當降雨與庫水位共同作用時需考慮降雨影響所占權重的大小，因此如何考慮降雨對滑坡變形滯后響應的影響還需進一步研究。