考慮滲流時滯及參數(shù)不確定性的土壩壩坡失穩(wěn)概率分析

徐鎮(zhèn)凱，蔡 磊，魏博文，王 鋒

(南昌大學(xué)建筑工程學(xué)院，江西南昌330031)

壩坡穩(wěn)定問題是土石壩設(shè)計及運(yùn)行管理關(guān)注的重點(diǎn)內(nèi)容之一。由于滲流問題具有復(fù)雜性，其影響因素繁多，水位、土體參數(shù)等都是影響其安全系數(shù)大小的因素。國內(nèi)外學(xué)者針對壩前水位的升降引起壩體內(nèi)部滲流場發(fā)生變化，進(jìn)而誘發(fā)的壩坡穩(wěn)定問題進(jìn)行了大量研究，鄭穎人等[1]、時衛(wèi)民等[2]研究了由水位下降導(dǎo)致的浸潤線變化為基礎(chǔ)引起邊坡失穩(wěn)的影響；唐輝明等[3]基于非飽和土強(qiáng)度理論，研究了庫水位在不同下降速率下的土石壩壩坡穩(wěn)定性，得出穩(wěn)定性系數(shù)極小值的出現(xiàn)與水位下降速度密切相關(guān)；Terzaghi[4]就庫水位變化時對滑坡體穩(wěn)定性的影響進(jìn)行了專門分析，并指出庫水位驟降坡腳會出現(xiàn)滲透變形的情況；文獻(xiàn)[5- 6]通過建立耦合模型分析庫水位快速下降引起的滑坡失穩(wěn)，同時分析對比了數(shù)值模擬結(jié)果與現(xiàn)場檢測數(shù)據(jù)。趙瑞欣等[7]通過分析三峽庫區(qū)涼水井滑坡，對水位變化速率分析，得出在同一滲透系數(shù)的前提下，水位驟降的風(fēng)險最大。

綜上所述，庫水位變化對浸潤線變化、土體強(qiáng)度的影響是一個復(fù)雜的過程，其導(dǎo)致土石壩壩坡土體經(jīng)歷飽和-非飽和的循環(huán)轉(zhuǎn)化，最終影響壩坡。現(xiàn)有研究中，蒙特卡洛法、一階矩和二階矩等方法被廣泛應(yīng)用于結(jié)構(gòu)可靠性分析中[8- 9]。其中，蒙特卡洛法雖模擬計算精度高，但其計算量巨大；而一階矩法和二階矩法計算精度亟待提升。壩坡穩(wěn)定分析多基于極限平衡法和極限分析法[10]進(jìn)行數(shù)值模擬，再結(jié)合蒙特卡洛模擬法計算結(jié)果，分析過程繁瑣，限制了其在工程中的應(yīng)用。

本文在GEOSTUDIO軟件建立有限元分析模型的基礎(chǔ)上，建立了考慮土石壩上游水位驟降的壩坡穩(wěn)定可靠度模型。通過對壩體進(jìn)行應(yīng)力分析，結(jié)合Kriging方法，推導(dǎo)了土石壩壩坡穩(wěn)定安全系數(shù)解析式與相應(yīng)可靠度功能函數(shù)，分析了土石壩壩坡穩(wěn)定性。工程實(shí)例分析表明，本文所建模型計算精度高、收斂速度快、計算分析簡便，可為土石壩健康服役提供理論依據(jù)與決策支持。

1 土石壩壩坡失穩(wěn)概率的分析方法

1.1 安全系數(shù)解析式

實(shí)際工程中，結(jié)構(gòu)正常使用極限狀態(tài)中抗力效應(yīng)R和荷載效應(yīng)S的功能函數(shù)關(guān)系式為

Z=g(R，S)=R-S

當(dāng)壩坡處于極限狀態(tài)時，R-S=0。國內(nèi)學(xué)者常結(jié)合安全系數(shù)和可靠度理論建立可靠度分析功能函數(shù)對壩坡進(jìn)行穩(wěn)定可靠性分析[11]，基于穩(wěn)定安全系數(shù)的可靠度分析功能函數(shù)表達(dá)式為

Z=Fs-1

式中，Z為功能函數(shù)；Fs為壩坡穩(wěn)定安全系數(shù)。

采用簡化畢肖普方法計算壩坡的抗滑穩(wěn)定安全系數(shù)，即

式中，c′i和φ′i分別為土的有效黏聚力和內(nèi)摩擦角；Wi和bi分別為土條的質(zhì)量和寬度；αi為土條滑面的傾角。

1.2 壩體內(nèi)部滲流的時滯特性

庫水位對土石壩滲流性態(tài)的影響是一個逐變過程，在水位變化情況下，引入水位滯后影響函數(shù)來反映庫水位時滯作用。假設(shè)Hi(i=1，2，3…，n)為當(dāng)前水位，Hd為能產(chǎn)生與上述水位等效滲流性態(tài)的等效水位

Hd=y(H1，H2，…Hi，…Hn，

w1，w2，…wi，…wn，)

由于滲流效應(yīng)是由庫水位引起的時變連續(xù)函數(shù)，故等效水位也為時變連續(xù)函數(shù)[12]。設(shè)相應(yīng)的影響程度分布密度函數(shù)為w(t)。根據(jù)大量的統(tǒng)計分析發(fā)現(xiàn)w(t)的一般規(guī)律接近正態(tài)分布(μ，σ2)，因此w(t)可表示為

式中，x1為庫水位滯后時間；x2為庫水位影響時間；A為調(diào)整參數(shù)。在一定時間范圍內(nèi)，x1，x2為定值，A常視為常數(shù)。可得

設(shè)在t=t0時刻其對應(yīng)水位的影響程度分布密度函數(shù)為w(t)，對應(yīng)等效水位為Hd，則

式中，H(t)為t時刻的水位。上式結(jié)合統(tǒng)計學(xué)理論即可求出x1，x2。

由于水庫長期蓄水，庫水位驟降時壩體內(nèi)的孔隙水不能及時從土體內(nèi)排出，壩前水位的突然下降，導(dǎo)致了壩前水壓突然減小。隨著時間變化，壩體水逐漸流失，致使壩體浮托力減小，同時孔隙水壓的增大，導(dǎo)致土料間的有效應(yīng)力的降低，抗剪強(qiáng)度隨之降低。上述兩方面影響因素協(xié)同驅(qū)動引起土體邊坡的穩(wěn)定性發(fā)生變異，導(dǎo)致壩坡失穩(wěn)破壞。

2 基于Kriging模型的可靠性功能函數(shù)

Kriging方法[13- 15]最早由南非學(xué)者Krige提出并應(yīng)用于地質(zhì)統(tǒng)計學(xué)研究。Romero等[16]分析對比結(jié)合了拉丁超立方抽樣的Kriging方法的結(jié)構(gòu)可靠度計算與其他插值技術(shù)方法的可靠度計算。我國學(xué)者張崎等[17]提出了一種用于結(jié)構(gòu)可靠度計算且具有較高計算效率的基于kriging模型的重要抽樣方法。佟操等[18]通過結(jié)合蒙特卡洛法和Kriging方法提出了一種有效提高可靠度計算效率及其精度的可靠度主動學(xué)習(xí)算法。Kriging方法是一種半?yún)?shù)化模擬方法，其表達(dá)式為

式中，f(x)T=[f1(x)，f2(x)，…，fn(x)]，為回歸模型；β=[β1，β2，...，βn]T，為回歸系數(shù)；z(x)為隨機(jī)誤差。Kriging模型與響應(yīng)面法最主要的不同點(diǎn)在于z(x)是一個隨機(jī)過程，服從正態(tài)分布N(0，σ2)，z(x)的協(xié)方差矩陣為

cov [z(w)，z(x)]=σ2R(θ，w，x)

式中，R(θ，w，x)是以θ為參數(shù)的相關(guān)模型，其中w，x∈Rn。R(θ，w，x)為任何兩個樣本點(diǎn)w和x的空間相關(guān)函數(shù)，它對模擬的精確程度起決定性作用。

線性函數(shù)、高斯函數(shù)、指數(shù)函數(shù)以及廣義指函數(shù)等都應(yīng)用于Kriging方法的建模過程，其中使用高斯函數(shù)相關(guān)函數(shù)的計算效果最佳[19]。本文采用高斯函數(shù)，則有

假定隨機(jī)抽樣m個坐標(biāo)點(diǎn)組成的矩陣為S=[s1，s2，…，sm]T；Y=[y1，y2，…，ym]T，為相應(yīng)坐標(biāo)點(diǎn)的響應(yīng)值矩陣。定義R=(Rij)m×m為S中坐標(biāo)點(diǎn)的相關(guān)矩陣，其中Rij=R(θ，si，sj)，(i，j=1，2，…，m)。構(gòu)造m維權(quán)系數(shù)向量，c=[c1，c2，…，cm]T，根據(jù)響應(yīng)值yi(i=1，2，…，m)的線性加權(quán)疊加來預(yù)測待測點(diǎn)x的響應(yīng)值，可得

要建立最優(yōu)Kriging模型即將最優(yōu)Kriging建模過程轉(zhuǎn)換為某一非線性無約束優(yōu)化問題的求解過程，通過求解上述優(yōu)化問題，求得參數(shù)θk后，即可建立最優(yōu)Kriging模型。

3 可靠度分析流程

土石壩有限元模型的建立涉及有限元應(yīng)力分布，土石壩壩坡安全系數(shù)FS同壩體材料物理參數(shù)之間的關(guān)系一般為隱式函數(shù)關(guān)系。利用Kriging方法并結(jié)合軟件GEOSTUDIO對土石壩壩坡隱式功能函數(shù)進(jìn)行擬合，實(shí)現(xiàn)對土石壩壩坡穩(wěn)定可靠度的分析，其計算流程圖如圖1所示。

圖1 基于Kriging模型的壩坡穩(wěn)定可靠度分析流程

為確定建立Kriging模型的最佳樣本點(diǎn)數(shù)，分別建立不同樣本點(diǎn)數(shù)模型，記錄Morgenstern-Price法(M-P法)并與2萬次蒙塔卡洛法(MC法)計算結(jié)果相對比，見表1。易知當(dāng)樣本點(diǎn)數(shù)為40時，相對誤差最小，故本文選取的樣本點(diǎn)數(shù)為40。

表1 M-P法與MC法對比

4 案例考證

土石壩土體的主要參數(shù)包括粘聚力c、容重γ、內(nèi)摩擦角φ。其中，c、φ對土石壩邊坡穩(wěn)定分析的影響程度較大，γ影響程度較小，故本文只考慮c、φ的變異影響。庫區(qū)水位驟降是邊坡穩(wěn)定性最不利的情況，極易導(dǎo)致災(zāi)害的發(fā)生。據(jù)統(tǒng)計，在大壩事故中有1/4是由于土石壩滑坡導(dǎo)致的潰壩事故，就滑坡位置而言，上游滑坡導(dǎo)致的潰壩數(shù)量遠(yuǎn)低于下游滑坡[20]，但在實(shí)際工程中必須防止因水位驟降引起的滑坡事故發(fā)生，本文選取土石壩上游水位驟降情況建立有限元模型。

4.1 模型的建立

某均質(zhì)土石壩，壩高50 m，壩基標(biāo)高50 m，壩基影響深度與壩高相等，壩基寬486 m，壩頂寬10 m。大壩運(yùn)行期正常蓄水位h1=98 m，死水位h2=86 m，大壩下游水位h3=52 m。馬道位于下游壩坡高程70 m處，寬1 m，上游壩坡坡比1∶m=1∶3，下游壩坡坡比均為1∶n=1∶2.5，下游設(shè)排水棱體。

圖4 土石壩有限元模型

4.2 材料參數(shù)的選取

土石壩有限元模型建立的過程中涉及壩體土、排水棱體材料及壩基土三種材料，其中，壩體土滲透系數(shù)k為6.50×10-7m/s、容重γ為20.20 kg/m3；排水棱體滲透系數(shù)k為16.40×10-6m/s、容重γ為19.00 kg/m3；壩基土滲透系數(shù)k為3.04×10-7m/s、容重γ為22.10 kg/m3。參數(shù)統(tǒng)計特性見表2，分布類型均為對數(shù)正態(tài)分布。

表2 材料參數(shù)統(tǒng)計特性

三種材料的體積含水量函數(shù)及滲透系數(shù)函數(shù)均采用Van Genuchten模型估算法進(jìn)行估算，壩體土體積含水量函數(shù)及滲透系數(shù)函數(shù)關(guān)系曲線見圖2、3。

圖2 壩體土體積含水量函數(shù)

圖3 壩體土滲透系數(shù)函數(shù)

土石壩幾何模型參數(shù)所建模型見圖4，土石壩進(jìn)行有限元網(wǎng)格的劃分時主要采用三角形及四邊形網(wǎng)格，橫斷面共劃分單元數(shù)1 288個，節(jié)點(diǎn)數(shù)1 385個。

4.3 計算結(jié)果分析

圖5給出了水位驟降開始前和結(jié)束后一段時間后的孔隙水壓力分布曲線。由圖5可以看出，在水位下降過程中，由于壩體水的排出，涂料陷落壓實(shí)，孔隙水壓有一個短暫的上升，隨后壩體水排出，孔隙水壓迅速下降，土體由飽和狀態(tài)逐漸變?yōu)榉秋柡蜖顟B(tài)，上游壩面負(fù)壓區(qū)范圍變大，正壓變小。

圖5 水位驟降前后孔隙水壓分布

利用MATLAB程序計算壩坡可靠度隨時間變化趨勢如圖6所示。由圖6分析可知：在水位驟降過程中，壩坡可靠度是時刻變化的，由于水位的變化，上游壩坡失效概率逐漸增大，可靠度逐漸降低，由于計算時間步長呈指數(shù)變化，導(dǎo)致計算結(jié)果形成突變，總體而言，水位驟降導(dǎo)致壩坡失穩(wěn)，可靠度逐漸降低，隨著時間變化，孔隙水的排出，可靠度有所增大，整體呈現(xiàn)出先減小后稍增大的趨勢，但變化不大，趨于穩(wěn)定。在采用Kriging法對上游游壩坡穩(wěn)定可靠度進(jìn)行分析時，上游水位的下降，壩坡失效概率逐漸增大，相應(yīng)的可靠度指標(biāo)逐漸減小，在水位下降過程中，存在最危險水位，但該水位難以確定，當(dāng)水位下降至該水位時，失效概率最大，相應(yīng)的可靠度指標(biāo)最小。

圖6 可靠度隨時間變化曲線

為驗證Kriging模型計算結(jié)果的準(zhǔn)確性，分別采用已建Kriging模型及蒙特卡洛(MC)法分析在穩(wěn)定滲流條件下，上游正常蓄水位、下游無水情況下的飽和—非飽和土石壩下游壩坡穩(wěn)定可靠度問題，其計算結(jié)果見表2。其中，MC法抽樣次數(shù)為10萬次。由表2中的計算結(jié)果可知，計算結(jié)果相對誤差為0.12%，該計算工況下的相對誤差較小，驗證了本文提出的利用Kriging方法計算可靠度結(jié)果精度滿足工程要求。

表3 壩坡穩(wěn)定可靠度分析結(jié)果

5 結(jié) 論

(1)研究了土石壩在水位驟降情況下，上游水位時滯效應(yīng)下的土壩滲流變化規(guī)律，得出了在水位驟降情況下等效水位的計算公式。分析得出體積含水量是關(guān)于孔隙水壓力的函數(shù)，在水位下降過程中，土體含水量隨時間逐漸變化，可靠度指標(biāo)隨時間逐漸變化。

(2)本文利用Kriging方法對土石壩壩坡穩(wěn)定可靠度問題進(jìn)行分析求解，Kriging方法擁有較為成熟的迭代規(guī)則，只需進(jìn)行較少運(yùn)算便可獲得較為精確的結(jié)果。與蒙特卡洛模擬法相比，有效地降低了有限元模型計算次數(shù)，提高了計算效率，與極限平衡法相比，提高了計算精度。

(3)綜合考慮了土石壩土體中基質(zhì)吸力與滲透系數(shù)之間的非線性關(guān)系，較為真實(shí)地反應(yīng)了土石壩應(yīng)力狀態(tài)，建立了土石壩可靠度功能函數(shù)，不僅適用于本文在水位驟降情況下的壩坡可靠度分析，同時適用于大部分水力邊界條件下復(fù)雜壩坡可靠度分析。