考慮滲流時滯及參數不確定性的土壩壩坡失穩概率分析

徐鎮凱，蔡 磊，魏博文，王 鋒

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

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

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

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

1 土石壩壩坡失穩概率的分析方法

1.1 安全系數解析式

實際工程中，結構正常使用極限狀態中抗力效應R和荷載效應S的功能函數關系式為

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

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

Z=Fs-1

式中，Z為功能函數；Fs為壩坡穩定安全系數。

采用簡化畢肖普方法計算壩坡的抗滑穩定安全系數，即

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

1.2 壩體內部滲流的時滯特性

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

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

w1，w2，…wi，…wn，)

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

式中，x1為庫水位滯后時間；x2為庫水位影響時間；A為調整參數。在一定時間范圍內，x1，x2為定值，A常視為常數?？傻?/p>

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

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

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

2 基于Kriging模型的可靠性功能函數

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

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

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

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

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

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

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

3 可靠度分析流程

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

圖1 基于Kriging模型的壩坡穩定可靠度分析流程

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

表1 M-P法與MC法對比

4 案例考證

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

4.1 模型的建立

某均質土石壩，壩高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，下游設排水棱體。

圖4 土石壩有限元模型

4.2 材料參數的選取

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

表2 材料參數統計特性

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

圖2 壩體土體積含水量函數

圖3 壩體土滲透系數函數

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

4.3 計算結果分析

圖5給出了水位驟降開始前和結束后一段時間后的孔隙水壓力分布曲線。由圖5可以看出，在水位下降過程中，由于壩體水的排出，涂料陷落壓實，孔隙水壓有一個短暫的上升，隨后壩體水排出，孔隙水壓迅速下降，土體由飽和狀態逐漸變為非飽和狀態，上游壩面負壓區范圍變大，正壓變小。

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

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

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

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

表3 壩坡穩定可靠度分析結果

5 結 論

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

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

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