考慮尾礦材料參數空間變異性的尾礦壩穩定可靠度分析*

蔣水華，朱明明，黃勁松

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

0 引言

隨著國民經濟和工程技術的快速發展，我國尾礦庫數量越來越多，規模越來越大。目前我國尾礦庫超過了12 000座，其中95%的尾礦庫采用上游法筑壩，尾礦庫同時也是礦業風險的主要來源之一[1]。近年來，地震、洪水等各種原因導致我國尾礦庫潰壩事故頻繁發生，危害極其嚴重。另外，尾礦材料是性質極其復雜的地質介質，在施工過程中人為與非人為因素的作用下，其物質組成和內部結構會發生一定的變化，尾礦材料性質在空間上相差較大，即呈現一定的空間變異性。尾礦材料參數空間變異性的客觀存在給尾礦壩穩定性分析帶來了一定的困難，常用的確定性分析方法因不能反映這種空間變異性的影響，會導致評價結果存在一定的偏差。

目前針對尾礦壩穩定可靠度問題開展了大量有益的研究工作，如胡平安等[1]采用JC法探究了不同尾礦材料力學參數的變異性對湖北省漁潭尾礦壩可靠度的影響；劉迪等[2]利用貝葉斯網絡模型建立尾礦壩穩定性因果關聯分析模型，分析了各因素變化對尾礦壩穩定的影響；曹志松等[3]采用二次多項式構建響應面模型，計算了黏聚力、內摩擦角和重度服從正態分布情況下的尾礦壩穩定可靠度；張揚等[4]采用蒙特卡羅模擬(MCS)方法探討了尾礦材料黏聚力、內摩擦角、密度等和孔隙水壓力的不確定性對獨木垅尾礦壩抗滑穩定的影響。雖然這些研究工作推動了可靠度理論在尾礦庫安全性評價中的應用，但是大多采用隨機變量模型模擬尾礦材料參數的不確定性，不能有效考慮空間不同點局部和整體尾礦材料物理力學性質的差異，無法客觀地評價尾礦材料參數空間變異性對尾礦壩穩定可靠度的影響。相比之下，隨機場理論將土層不同位置上的土壤參數模擬為一個服從某一概率分布的隨機變量，并與相鄰位置上的隨機變量相關，非常適合表征尾礦材料的空間變異性，目前在邊坡工程中得到了廣泛應用[5-7]，但是較少應用到尾礦材料參數空間變異性表征中。此外，對于需考慮尾礦材料參數空間變異性的復雜尾礦壩穩定問題，MCS等傳統可靠度分析方法的計算精度和效率欠佳，因此需要發展一種高效的考慮參數空間變異性的尾礦壩穩定可靠度分析方法。

本文結合隨機場理論、代理模型和結構可靠度理論，采用考慮尾礦材料參數空間變異性尾礦壩可靠度分析的非侵入式隨機有限元法，通過Karhunen-Loève(K-L)展開方法模擬尾礦材料物理力學參數空間變異性，借助Hermite隨機多項式展開建立尾礦壩安全系數代理模型，在此基礎上計算尾礦壩邊坡穩定可靠度，并以實際尾礦壩工程為例驗證了提出方法的有效性[5-7]。

1 尾礦材料參數空間變異性模擬

采用隨機場模擬尾礦材料參數空間變異性時，重要的一步是進行隨機場離散，本文選用計算精度和效率較高的K-L級數展開方法離散尾礦材料參數隨機場。

以離散尾礦滲透系數ks與內摩擦角φ各向異性對數正態參數隨機場為例，簡要介紹基于K-L級數展開方法的尾礦材料參數空間變異性模擬過程。K-L展開方法離散隨機場是以譜分解有界，對稱且正定的相關函數為基礎，將土體參數隨機場的離散轉化為求解Fredholm積分方程的特征值問題[6]，即

(1)

式中：(x1,y1)和(x2,y2)表示計算區域Ω中的任意2個坐標點；ρ[(x1,y1),(x2,y2)]為任意2點(x1,y1)，(x2,y2)處隨機場特性值之間的自相關系數；fi(.)和λi分別為自相關函數對應的特征函數和特征值。

采用描述參數空間自相關性的高斯型自相關函數所模擬的參數隨機場分布平滑度和連續性較好[8]，故本文選取高斯型自相關函數，表達式為：

(2)

式中：lh和lv分別表示水平和垂直方向上的相關距離，用于表征參數空間變異性的自相關程度，相關距離越大對應的尾礦材料參數的空間自相關性越強。

根據式(2)可計算得到原始空間不同點處的隨機場特性值之間的自相關系數。一般水平相關距離是垂直相關距離的10倍左右，水平相關距離在10～40 m范圍內，垂直相關距離變化范圍為1～3 m[9]。

當采用式(2)高斯型自相關函數時，式(1)沒有解析解，為此本文采用wavelet-Galerkin技術[10]進行數值求解?；谧韵嚓P函數的特征解可將隨機場H(x,y)展開為：

(3)

式中：(x,y)為隨機場區域中的任意坐標點，(x,y)∈Ω；ξXi,j為參數隨機場離散的獨立標準正態隨機變量；μlnXi和σlnXi分別為相應參數隨機場lnXi的均值和標準差；n為截斷項數，與隨機場離散的隨機變量數目對應。截斷項數n取決于精度要求和自相關函數的形式，文獻[11]～[12]建議采用隨機場期望能比率因子ε≥95%作為確定截斷項數n取值的依據。

2 尾礦壩可靠度非侵入式隨機分析方法

通過K-L展開方法得到參數隨機場特性值后,需要進行含大量隨機變量的尾礦壩邊坡可靠度分析。尾礦壩滲流分析和邊坡穩定一般借助有限元軟件實現，因此尾礦壩安全系數與尾礦材料參數之間是非線性隱式函數關系。傳統的MCS、一次二階矩等可靠度方法不能有效求解含大量隨機變量的極限狀態函數是隱式的尾礦壩可靠度問題。為此，本文將近幾年發展起來的非侵入式隨機有限元法[8]拓展到尾礦庫工程領域，采用Hermite隨機多項式展開(Hermite Polynomial Chaos Expansion，PCE)代替尾礦壩邊坡安全系數FS與尾礦材料參數ξ之間的非線性隱式函數關系：

(4)

式中：n為隨機變量的數目；ξ=(ξ1,ξ2,…,ξn)T為獨立標準隨機向量，與式(1)參數隨機場離散的隨機變量一一對應；Γp(.)表示為p階Hermite隨機多項式展開；a0,ai1,ai1i2,ai1i2i3,…為待定系數，其數目為(n+p)!/(n!p!)。

接著，采用拉丁超立方抽樣技術(Latin Hypercube Sampling，LHS)[13]產生輸入參數隨機樣本點，代入有限元程序或巖土商業軟件計算安全系數，再根據式(4)建立線性代數方程組求解多項式展開系數。一旦確定了尾礦壩安全系數和輸入參數ξ的顯式函數關系，即安全系數代理模型，便可建立新的顯式函數表達的尾礦壩穩定可靠度分析極限狀態函數：

G(ξ)=FS(ξ)-1.0

(5)

基于采用Hermite隨機多項式展開建立的代理模型，利用100 000次MCS方法計算尾礦壩失穩破壞概率。顯然，該方法不僅實現了可靠度分析和確定性有限元分析的不耦合，而且直接使用代理模型計算安全系數FS，無需進行大量的有限元計算，提高了計算效率。

3 工程應用

以西木尾礦庫工程[14]為例，探討尾礦壩材料滲透系數和內摩擦角空間變異性對尾礦壩穩定可靠度的影響。西木礦區位于加拿大魁北克省魯恩-諾蘭達以東約40 km的鮑斯凱鎮，于2004年被發現。西木尾礦壩壩高15 m，大壩在建設結束時的幾何形狀，初始模型中各層材料分布、層厚和層寬如圖1所示。各個土層的物理力學參數見表1。土層下限是基巖和過密的冰磧層，在大壩建設之后保持平衡，被認為是穩定地層。

圖1 尾礦壩材料剖面及有限元網格Fig.1 Material profile and finite element mesh of the tailing dam

表1 尾礦材料物理力學參數Table 1 Physical and mechanical parameters of tailings materials

注：γ，φ，c和ks分別為尾礦材料的重度、內摩擦角、黏聚力和滲透系數；a,n和m為土水特征曲線擬合參數。

3.1 有限元模型及計算參數

西木尾礦壩有限元模型如圖1所示，一共劃分為10 685個節點和10 407個網格大小為0.5 m的四邊形和三角形混合單元。邊界條件為：上游水位(左側)13 m，下游水位(右側)8 m，并在下游邊坡(右側)設置了潛在滲流面?；赟EEP/W模塊對西木尾礦壩進行飽和滲流分析，再將得到浸潤線和孔隙水壓分布等滲流結果導入到SLOPE/W模塊，采用摩根斯坦-普萊斯法計算的安全系數FS=1.931，最危險滑動面如圖2所示。所計算的安全系數和臨界滑移面與Coulibaly等[14]計算的安全系數1.967和臨界滑移面基本一致，說明了該確定性有限元分析的有效性。

因受基質吸力、含水量、滲透率等因素的影響，因此對西木尾礦壩進行飽和-非飽和滲流分析十分必要。本文將填石和尾細砂視作非飽和材料，采用目前應用廣泛的Fredlund-Xing[15]模型表征的土水特征曲線(SWCC)模擬體積含水量與基質吸力之間的函數關系。模型擬合參數a為進氣值相關的土性參數；n為與SWCC斜率相關的土性參數；m為與殘余含水率相關的土性參數[16]，參數取值見表1。圖3(a)和3(b)分別給出了尾細砂土水特征曲線與滲透系數函數。將飽和-非飽和滲流分析結果導入到SLOPE/W模塊，采用摩根斯坦-普萊斯法計算得到西木尾礦壩安全系數FS=1.944，大于在所有尾礦材料為飽和狀態下計算的安全系數。

圖2 尾礦壩穩定性分析結果Fig.2 Stability analysis results of the tailings dam

圖3 尾細砂土水特征與滲透系數函數曲線Fig.3 Soil-water characteristic curve and permeability coefficient function curve of tailings fine sand

另外，尾礦壩上的局部坍塌和滑移大多是由于地震導致尾部或地基土液化引起的[14]。為此，本文考慮Ⅶ級地震作用，進一步對該尾礦壩進行穩定可靠度分析。根據建設部頒布的《關于統一抗震設計規范地面運動加速度設計取值的通知》規定，Ⅶ級地震設計基本地震加速度取0.1 g，豎向地震影響系數的最大值取水平地震影響系數最大值的65%[17]。采用擬靜力法模擬Ⅶ級地震作用，同樣通過非飽和滲流穩定分析得到安全系數FS=1.418，圖4給出了最危險滑面位置。根據尾礦設施設計規范[18]可知，該尾礦壩在Ⅶ級地震作用下仍處于穩定狀態。

圖4 地震作用下尾礦壩穩定性分析結果Fig.4 Stability analysis results of tailings dam under effect of earthquake

3.2 尾礦壩穩定可靠度分析

如前所述，受到尾礦壩長期多循環水力充填以及固結沉降的作用，尾礦材料強度參數存在較大的離散性和空間變異性，與尾礦材料抗剪強度參數空間變異性相比，尾礦材料重度空間變異性較小，一般小于10%，可視作常量。故本例僅考慮尾礦材料滲透系數ks和內摩擦角φ的空間變異性，并假設尾細砂、粉土、粉質黏土和填石的ks，φ均服從對數正態分布，尾礦壩材料參數統計特征見表2。

表2 尾礦材料參數統計特征Table 2 Statistical characteristics of tailings materials parameters

圖5 尾細砂滲透系數隨機場的典型實現Fig.5 Typical realization of random field for permeability coefficient of tailings fine sand

在此基礎上，利用LHS法對輸入的變量進行抽樣，根據參數統計特征生成隨機變量和隨機場數據，分別賦到尾礦壩模型，然后通過滲流有限元和穩定分析得到FS，再采用Hermite展開多項式建立FS代理模型(即FS與隨機變量和隨機場離散的隨機變量之間的顯式函數關系)。圖6比較了由不同方法計算的極限狀態函數累積分布函數(CDF)曲線。由圖6可見，分別進行2 000，4 000和10 000次滲流有限元及穩定性分析構建FS代理模型，對應的二階PCE+MCS方法計算的極限狀態函數的CDF曲線非常吻合，并且與10 000次直接LHS方法的計算結果僅僅在低概率區域存在微小差別。其原因是10 000次直接LHS方法對于低失穩破壞概率(<10-3)水平尾礦壩可靠度問題，計算精度不夠。

圖6 極限狀態函數累積分布曲線的比較Fig.6 Comparison on cumulative distribution curves of limit state functions

由圖6極限狀態函數的CDF曲線可以直接計算尾礦壩失穩破壞概率。表3比較了不同方法計算的尾礦壩失穩破壞概率。從表3可以看出，西木尾礦壩失穩破壞概率水平很低，二階PCE+MCS方法(Np=2 000)，二階PCE+MCS方法(Np=4 000)和二階PCE+MCS(Np=10 000)方法計算的破壞概率分別為8.6×10-4,8.1×10-4，7.6×10-4，可見3種方法計算的破壞概率非常吻合，和10 000次直接LHS方法計算的9.9×10-4也基本一致，說明了本文非侵入式隨機有限元法的有效性。此外，與10 000次直接LHS方法相比，二階PCE+MCS方法只需要進行2 000次尾礦壩滲流有限元和穩定性分析就可以滿足計算精度要求，計算量是直接LHS方法的1/4，可見本文提出的方法具有較高的計算效率，為解決考慮多個尾礦材料參數空間變異性的低概率水平尾礦壩穩定可靠度問題提供了一條有效路徑。

表3 尾礦壩可靠度分析結果Table 3 Reliability analysis results of the tailings dam

4 結論

1)采用隨機場理論模擬西木尾礦材料參數空間變異性，采用Hermite隨機多項式建立了尾礦壩安全系數與不確定性輸入參數之間的代理模型，將尾礦壩滲流有限元及穩定性分析視作黑箱子直接調用，實現了尾礦壩可靠度分析與確定性穩定分析的不耦合。

2)與直接LHS方法相比，提出方法滿足計算精度要求的同時，計算量減少了80%，為解決考慮多個尾礦材料參數空間變異性的低概率水平尾礦壩穩定可靠度問題提供了一條有效的路徑。

3)盡管西木尾礦壩在正常和地震工況下的安全系數很高，失穩破壞概率低至10-4量級，與尾礦壩設計規范給出的允許安全系數和破壞概率相比，該尾礦壩在Ⅶ級地震作用下發生失穩破壞的可能性較小，但是仍然需要設計一些簡單的尾礦壩安全防護抗震措施。