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

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

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

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

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

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

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

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

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

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

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)為當前水位，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ù)?？傻?/p>

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

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

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

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

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

式中，f(x)T=[f1(x)，f2(x)，…，fn(x)]，為回歸模型；β=[β1，β2，...，βn]T，為回歸系數(shù)；z(x)為隨機誤差。Kriging模型與響應(yīng)面法最主要的不同點在于z(x)是一個隨機過程，服從正態(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)為任何兩個樣本點w和x的空間相關(guān)函數(shù)，它對模擬的精確程度起決定性作用。

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

假定隨機抽樣m個坐標點組成的矩陣為S=[s1，s2，…，sm]T；Y=[y1，y2，…，ym]T，為相應(yīng)坐標點的響應(yīng)值矩陣。定義R=(Rij)m×m為S中坐標點的相關(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ù)測待測點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ù)進行擬合，實現(xiàn)對土石壩壩坡穩(wěn)定可靠度的分析，其計算流程圖如圖1所示。

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

為確定建立Kriging模型的最佳樣本點數(shù)，分別建立不同樣本點數(shù)模型，記錄Morgenstern-Price法(M-P法)并與2萬次蒙塔卡洛法(MC法)計算結(jié)果相對比，見表1。易知當樣本點數(shù)為40時，相對誤差最小，故本文選取的樣本點數(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ù)量遠低于下游滑坡[20]，但在實際工程中必須防止因水位驟降引起的滑坡事故發(fā)生，本文選取土石壩上游水位驟降情況建立有限元模型。

4.1 模型的建立

某均質(zhì)土石壩，壩高50 m，壩基標高50 m，壩基影響深度與壩高相等，壩基寬486 m，壩頂寬10 m。大壩運行期正常蓄水位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模型估算法進行估算，壩體土體積含水量函數(shù)及滲透系數(shù)函數(shù)關(guān)系曲線見圖2、3。

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

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

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

4.3 計算結(jié)果分析

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

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

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

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

為驗證Kriging模型計算結(jié)果的準確性，分別采用已建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ù)，在水位下降過程中，土體含水量隨時間逐漸變化，可靠度指標隨時間逐漸變化。

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

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