庫水變動下水庫滑坡變形預測研究

程江濤

中化地質(zhì)礦山總局湖北地質(zhì)勘查院，湖北 武漢 430075

水庫滑坡是岸坡內(nèi)地下水動力學條件改變而產(chǎn)生的一種地質(zhì)災害，主要發(fā)生在水庫建設(shè)與運營期間，由水位上升、水位驟降或降雨等導致[1]。水庫滑坡是我國三峽庫區(qū)較為常見的一種地質(zhì)災害，據(jù)不完全統(tǒng)計，三峽庫區(qū)范圍存在體積大小不等的各類滑坡、崩塌體約2490處，其中在175m水位范圍內(nèi)的滑坡體約有1630處[2]，占不穩(wěn)定斜坡體數(shù)量的65.5%。大量滑坡的存在勢必給三峽水庫的正常運營帶來巨大的安全隱患，其變形失穩(wěn)的危害主要體現(xiàn)在兩個方面[3-4]：①滑坡失穩(wěn)形成的大量巖土體涌入水庫內(nèi)，導致水庫不能正常運營；②高速滑入水庫內(nèi)的巖土體將會產(chǎn)生巨大涌浪，當滑坡體位于近壩庫岸時，由涌浪產(chǎn)生的巨大沖擊荷載將可能導致大壩在瞬間垮塌，同時涌浪波及范圍內(nèi)的城鎮(zhèn)居民往往來不及轉(zhuǎn)移，給庫區(qū)人民的生命財產(chǎn)造成嚴重損失。

國內(nèi)外有關(guān)水庫滑坡失穩(wěn)的實例屢見不鮮：1963年10月發(fā)生于意大利瓦依昂（Vajont）水庫左岸巨型滑坡，20s內(nèi)有超過 2億 m3的滑體以28m/s的速度飛越寬約 80m的河谷，造成水庫下游約2600人死亡，同時滑體將壩前 1.8km范圍內(nèi)的庫容全部填滿，導致整個瓦依昂水庫完全報廢[5]；2003年7月發(fā)生于三峽庫區(qū)支流青干河左岸的千將坪特大型巖質(zhì)滑坡，約有1500萬m3滑體在5min內(nèi)滑入水庫內(nèi)，導致14人死亡，10人失蹤，造成的直接經(jīng)濟損失高達8000多萬元[6-7]。以上典型事例表明，水庫滑坡失穩(wěn)造成的后果往往是災難性的。

水庫滑坡一直是水利水電工程與工程地質(zhì)領(lǐng)域研究的熱點與難點問題，并得到了廣大學者與業(yè)內(nèi)人士的高度重視[8-9]：在1992年召開的第六屆國際滑坡專題研討會議上，水庫滑坡的穩(wěn)定性問題被作為一個重點進行了專題討論[10]；2002年ICOID為水庫滑坡的勘察與評價建立了專門的分析指南[11]?；碌淖冃晤A測預報問題的難度就更大，滑坡變形監(jiān)測資料是反映內(nèi)外系統(tǒng)相互作用的有效信息來源，如何通過監(jiān)測資料揭示滑坡內(nèi)外因素相互作用的過程是滑坡變形預測的關(guān)鍵。自 1968年日本學者齋騰迪孝提出基于蠕變破壞三階段的經(jīng)驗預報模型以來，經(jīng)國內(nèi)外學者近40余年的苦心探索，已開拓出以統(tǒng)計分析、非線性科學、系統(tǒng)綜合分析以及實時動態(tài)跟蹤為理論基礎(chǔ)的綜合預測預報模型[12-15]。

在已有研究中存在以下幾方面問題：①基于巖土體微觀變形機制的力學預測模型是建立在室內(nèi)物理力學試驗基礎(chǔ)上的，由于滑坡是一個離散體，位于不同空間部位的土體其物質(zhì)組成與結(jié)構(gòu)是不同的，與之相對應的力學變形機制也是不同的，采用離散部位的室內(nèi)試驗結(jié)果代替復雜的滑坡系統(tǒng)顯然只適用于均質(zhì)土坡；②基于時序分析的統(tǒng)計預測模型是建立在滑坡變形監(jiān)測信息基礎(chǔ)上的，由于變形監(jiān)測信息中不僅包含有滑坡內(nèi)動力物質(zhì)系統(tǒng)產(chǎn)生的變形，而且還包含有滑坡外動力激勵系統(tǒng)產(chǎn)生的變形，而外動力激勵系統(tǒng)是一個隨機的系統(tǒng)，單純依靠時序監(jiān)測數(shù)據(jù)建立的預測模型顯然難以反映滑坡變形演化的內(nèi)在機理；③基于滑坡變形來源機制的動力成因預測模型是建立在滑坡變形影響因素基礎(chǔ)上的，能較好反映滑坡內(nèi)外系統(tǒng)相互作用的過程，但針對水庫滑坡開展這方面的預測模型研究卻較為少見，需進一步進行深入研究。

本文以李家灣滑坡為研究對象，以實測的變形監(jiān)測數(shù)據(jù)為基礎(chǔ)，利用逐步回歸分析的數(shù)學方法，試圖從影響因素與變形之間的相關(guān)性上來構(gòu)建滑坡的變形預測模型。通過研究，不僅能為水庫滑坡的變形預測提供一些新的思路，同時為三峽庫區(qū)水庫滑坡的防災減災及監(jiān)測預警提供一些行之有效的方法。

1 滑坡概況

李家灣滑坡位于長江左岸斜坡地帶，距三峽大壩壩址71.73km，行政區(qū)劃隸屬巴東縣東瀼口鎮(zhèn)舊縣坪村一組。地理坐標為：X 3437798，Y 37439310；E 31°03′34″，N 110°21′51″。研究區(qū)內(nèi)主要的地表水系為自西向東橫穿巴東縣中部的長江水系，西起邊域溪河口，東至?？谌腼鰵w縣境，流經(jīng)縣境長約38km，多年平均流量4487.7 億 m3。

滑坡區(qū)地處長江北岸，地形地貌類型屬構(gòu)造剝蝕中低山峽谷地貌，地勢總體上呈現(xiàn)出北東高、西南低的空間分布特征?；潞缶壐叱?00m，前緣最低高程90m左右，相對高差約210m。坡體中部200～210m高程處修筑有雷太公路，公路上部切坡地段形成4～8m高的陡坡（圖1）。李家灣滑坡區(qū)出露的地層主要為中生界中三疊統(tǒng)巴東組第二段（T2b2）陸相層狀碎屑巖地層和第三段（T2b3）淺海相碳酸鹽巖地層，以及上覆新生界第四系（Q4）松散堆積層。

圖1 滑坡全貌圖Fig.1 Landslide panorama

2 回歸分析模型

滑坡變形監(jiān)測分析的目地主要有兩方面：一方面是了解滑坡變形的主要影響因素；另一方面是找出滑坡變形的演化規(guī)律及其發(fā)展趨勢，在遇到異常情況時能及時發(fā)出安全警報，以確?；轮苓吔ㄖ锛熬用竦陌踩?。目前，大多數(shù)滑坡變形預測模型都是建立在位移數(shù)據(jù)的時序分析和數(shù)學構(gòu)造基礎(chǔ)上的，未能考慮滑坡變形影響因素的作用效果，而工程實踐中大多數(shù)滑坡的位移突變都是由影響因素的作用引起的。為此，有必要從滑坡變形影響因素的角度來建立滑坡變形的預測模型。

按照影響因素作用方式的不同可將滑坡變形分為兩部分，一部分是由滑坡內(nèi)動力物質(zhì)系統(tǒng)控制的以大時間尺度表現(xiàn)的時效變形，該類變形通常具有不可逆性；另一部分是由滑坡外動力激勵系統(tǒng)控制的以小時間尺度表現(xiàn)的隨機變形，該類變形通常具有可逆性，對于水庫滑坡而言，影響滑坡變形的外動力激勵因素主要有庫水變動和大氣降雨兩種。為此，可通過時效[16]、庫水位和降雨量[17]三個影響因子來構(gòu)建滑坡的變形預測模型[18]，其表達形式如下：

式中：S為滑坡總變形量，mm；Sθ為時效變化引起的變形分量，mm；SH為庫水位變化引起的變形分量，mm；SP為降雨量變化引起的變形分量，mm；ε為統(tǒng)計常數(shù)，可通過回歸分析得到。

2.1 時效引起的變形分量

時效變形是由滑坡內(nèi)動力物質(zhì)系統(tǒng)作用所產(chǎn)生的那部分變形，一般由不同時刻滑坡體所受的自重應力和抗滑阻力所決定，該類變形具有不可逆性，其大小與滑坡體所經(jīng)歷的累積時間相關(guān)，通常可用多項式模型來表征，其具體函數(shù)如下：

式中：θ為時效因子，為觀測日期至基準時間的天數(shù)除以 100得到的值；a1、a2、a3為時效變化引起的變形分量回歸系數(shù)。

2.2 庫水位引起的變形分量

庫水位變形是指由庫水變動過程中的地下水水力學效應產(chǎn)生的那部分變形，一般由庫水位漲落引起的水位差所決定，該類變形具有可逆性。已有的研究表明，庫水位變化引起的變形可通過庫水位漲落高差h的一、二、三次函數(shù)來表示，其中庫水位變化引起的撓曲變形與h3有關(guān)，庫水位變化引起的剪切變形則與h2有關(guān)，由水體自重引起的沉降變形與h有關(guān)[19]?？紤]到庫水變動作用對滑坡變形的滯后性，可通過指數(shù)函數(shù)構(gòu)建該變形分量，具體如下：

式中：h為庫水位漲落高差，為觀測當日的庫水位與前期庫水位之差除以100得到的值；b1、b2、b3為庫水位變化引起的變形分量回歸系數(shù)。

2.3 降雨量引起的變形分量

降雨量變形是指由降雨入滲過程中的地下水水力學效應產(chǎn)生的那部分變形，一般由降雨量引起的水頭差所決定，該類變形具有可逆性。已有的研究表明[20]，滑坡發(fā)生前 1～2個月的有效降雨量均會導致滑坡變形的發(fā)展，為此可通過當前降雨量、前1月平均降雨量和前2月平均降雨量作為衡量該變形的指標。考慮到降雨入滲作用對滑坡變形的滯后性，可通過指數(shù)函數(shù)構(gòu)建該變形分量，其具體形式如下：

式中：P0為觀測日期當月平均降雨量除以100得到的值；P1為觀測日期至前1月平均降雨量除以100得到的值；P2為觀測日期至前2月平均降雨量除以100得到的值；c1、c2、c3為降雨量變化引起的變形分量回歸系數(shù)。

通過上述三種變形分量構(gòu)建的滑坡總變形預測模型如下：

根據(jù)上述回歸預測模型，結(jié)合實測的影響因子參數(shù)（自變量）和位移值（因變量），即可通過最小二乘原理獲取上述回歸系數(shù)，進而建立基于影響因子回歸分析的變形預測模型。

3 逐步回歸分析模型的求解

3.1 逐步回歸分析的基本原理

回歸分析是研究因變量y與自變量x=（x1，x2，…，xn）之間相關(guān)關(guān)系的一種數(shù)理統(tǒng)計方法?；伦冃位貧w分析的目地不僅僅是建立因變量與自變量之間的函數(shù)關(guān)系，更重要的是抓住滑坡變形的主要影響因子。在多元回歸分析中，對于不同自變量的組合，通?？梢缘玫讲煌幕貧w方程，為此有必要從眾多的影響因子中遴選出顯著性的自變量，以建立最優(yōu)的回歸方程，同時該回歸方程也能反映出滑坡變形的主要影響因子，逐步回歸分析正是實現(xiàn)上述功能的一種有效統(tǒng)計方法。逐步回歸分析主要通過兩個方面得以實現(xiàn)：一方面，按自變量對因變量的顯著性程度，按照從大到小的順序逐個引入回歸方程；另一方面，當先前引入的自變量由于后面自變量的引入而使得回歸方程變得不顯著時，就將該自變量剔除。由此可見，逐漸回歸分析是分步進行的，在每一步回歸分析過程中有時會引入自變量，有時則會剔除自變量，這取決于引入自變量對回歸方程的顯著性，即逐步回歸分析中的每一步均要作統(tǒng)計檢驗，以確?；貧w方程中只包含顯著因子。由于最終建立的回歸方程中只包含顯著性的自變量，因此其預測效果較好。

3.2 逐步回歸分析的求解步驟

首先計算出各自變量之間的相關(guān)系數(shù)矩陣，初步分析各自量之間的相關(guān)性，觀察自變量之間是否有多重共線性現(xiàn)象。假定初始回歸方程為：

利用最小二乘法求得統(tǒng)計回歸方程為：

利用實測值與回歸值計算初始模型的殘差平方和為：

假定從初始回歸方程中剔除x1自變量，則經(jīng)過回歸分析求得的新模型為：

利用實測值與回歸值計算新模型的殘差平方和：

計算兩個模型的殘差平方和之差，該值反映了x1自變量對回歸方程貢獻的大?。?/p>

構(gòu)造統(tǒng)計檢驗量，檢驗x1自變量對回歸方程的顯著性：

式中：k-觀測次數(shù)；n-回歸方程中自變量個數(shù)。

4 預測模型的建立與分析

4.1 建模數(shù)據(jù)的選取與處理

李家灣滑坡的變形監(jiān)測資料顯示，滑坡前部監(jiān)測點的變形主要受庫水變動的影響，滑坡中部監(jiān)測點的變形則受降雨和庫水變動的聯(lián)合作用影響，滑坡后部監(jiān)測點的變形則主要受大氣降雨的影響。為研究不同影響因素作用下的滑坡變形預測模型，選取李家灣滑坡主剖面上的三個地表位移監(jiān)測點（BDT-35、BDT-36、BDT-37）作為研究對象（圖2，編號1、2、3分別代表BDT-37、BDT-36、BDT-35監(jiān)測點），截取 2007-8-8～2009-7-13之間的25組監(jiān)測數(shù)據(jù)作為樣本，其中以2007-8-8～2009-2-13之間的20組監(jiān)測數(shù)據(jù)作為建模樣本，為檢驗模型的預測效果，以2009-3-13～2009-7-13之間的 5組監(jiān)測數(shù)據(jù)作為預測樣本，觀測基準時間為 2007-8-8，各組樣本的數(shù)據(jù)信息如表1所示。利用表中數(shù)據(jù)建立的變形影響回歸因子如表2所示。

圖2 滑坡數(shù)值計算模型Fig.2 Numerical calculation model of landslide

表1 建模樣本數(shù)據(jù)Table 1 Modeling sample data

表2 影響因子參數(shù)Table 2 Influencing factor parameters

續(xù)表2

4.2 預測結(jié)果分析

4.2.1 滑坡后部 BDT-35監(jiān)測點變形預測模型與影響因素分析

通過回歸分析建立的 BDT-35監(jiān)測點變形預測模型如下：

預測結(jié)果見表3。通過該表可以看出，基于影響因子分析建立的回歸預測模型不僅能很好的擬合歷史監(jiān)測數(shù)據(jù)，而且其預測誤差在 1.27%～10.1%之間，具有較好的預測效果，但基于模型的時效性，應不斷利用新的監(jiān)測數(shù)據(jù)建立新的回歸預測模型才能保證預測的可靠性。

表3 BDT-35監(jiān)測點水平位移預測成果表Table 3 Horizontal displacement prediction results of BDT-35 monitoring point

4.2.2 滑坡中部 BDT-36監(jiān)測點變形預測模型與影響因素分析

通過回歸分析建立的 BDT-36監(jiān)測點變形預測模型如下：

預測結(jié)果見表4。通過該表可以看出，基于影響因子分析建立的回歸預測模型不僅能很好的擬合歷史監(jiān)測數(shù)據(jù)，而且其預測誤差在 0.41%～4.69%之間，具有較好的預測效果。

表4 BDT-36監(jiān)測點水平位移預測成果表Table 4 Horizontal displacement prediction results of BDT-36 monitoring point

4.2.3 滑坡前部 BDT-37監(jiān)測點變形預測模型與影響因素分析

通過回歸分析建立的 BDT-37監(jiān)測點變形預測模型如下：

預測結(jié)果見表5。通過該表可看出，基于影響因子分析建立的回歸預測模型不僅能很好的擬合歷史監(jiān)測數(shù)據(jù)，而且其預測誤差在0.80%～3.90%之間，具有較好的預測效果，但基于模型的時效性，應不斷利用新的監(jiān)測數(shù)據(jù)建立新的回歸預測模型才能保證預測的可靠性。

表5 BDT-37監(jiān)測點水平位移預測成果表Table 5 Horizontal displacement prediction results of BDT-37 monitoring point

綜上所述，李家灣滑坡后部BDT-35監(jiān)測點的變形主要是由時效變形和降雨量變形引起的，且降雨量對該部位變形的影響較??；滑坡中部BDT-36監(jiān)測點的變形主要是由時效變形、庫水位變形和降雨量變形共同引起的，且降雨量對該部位變形的影響較庫水位大；滑坡前部 BDT-37監(jiān)測點的變形主要是由時效變形和庫水位變形引起的。另外，從時效變形中影響因子的成分上看，BDT-35和BDT-37監(jiān)測點變形回歸方程中均含有θ3因子，而BDT-36監(jiān)測點變形回歸方程中含有θ2因子，這說明李家灣滑坡后部和前部滑體處于加速蠕變階段，而中部滑體則處于勻速蠕變階段，即滑坡中部滑帶尚未貫通。

5 結(jié)論

（1）滑坡變形是滑坡內(nèi)動力物質(zhì)系統(tǒng)與外動力激勵系統(tǒng)共同作用的結(jié)果，按照變形影響因素作用方式的不同可將其分解為時效變形、庫水位變形和降雨量變形；其中時效變形反映內(nèi)動力物質(zhì)系統(tǒng)的作用，該變形具有不可逆性，是時間的函數(shù)，可用多項式函數(shù)來表征；庫水位變形和降雨量變形反映外動力激勵系統(tǒng)的作用，該變形具有可逆性，其中庫水位變形是庫水位漲落高差的函數(shù)，降雨量變形是前期1～2個月累計平均降雨量的函數(shù)，可用指數(shù)函數(shù)來表征，上述三種變形分量的組合即為滑坡的總變形。

（2）BDT-35監(jiān)測點的變形影響因素主要為時效因子和降雨量變化因子，其中時效變形分量隨著時間的推移呈現(xiàn)出三次方增長趨勢；降雨量變形分量與前1個月的累積平均降雨量相關(guān)，其變化幅值僅為1.14mm，與該部位附近地下水位監(jiān)測孔基本為干孔的事實相吻合。

（3）BDT-36監(jiān)測點的變形影響因素主要為時效因子、庫水位變化因子和降雨量因子，其中時效變形分量隨著時間的推移呈現(xiàn)出二次方增長趨勢；庫水位變形分量在水位上升過程中產(chǎn)生趨向坡內(nèi)的變形，在水位下降過程中產(chǎn)生趨向坡外的變形；降雨量變形分量與前2個月的累積平均降雨量相關(guān)，且降雨量變形分量的變化幅值較庫水位變形分量大，與該部位附近地下水位監(jiān)測孔為季節(jié)性孔的事實相吻合。

（4）BDT-37監(jiān)測點的變形影響因素主要為時效因子和庫水位變化因子，其中時效變形分量隨著時間的推移呈現(xiàn)出三次方增長趨勢；庫水位變形分量在水位上升過程中產(chǎn)生趨向坡內(nèi)的變形，在水位下降過程中產(chǎn)生趨向坡外的變形，與該部位附近地下水位監(jiān)測孔與庫水位具有一致性的事實相吻合。