蔣水華，朱明明，黃勁松

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

0 引言

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

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

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

1 尾礦材料參數(shù)空間變異性模擬

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

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

(1)

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

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

(2)

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

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

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

(3)

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

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

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

(4)

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

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

G(ξ)=FS(ξ)-1.0

(5)

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

3 工程應(yīng)用

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

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

表1 尾礦材料物理力學(xué)參數(shù)Table 1 Physical and mechanical parameters of tailings materials

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

3.1 有限元模型及計算參數(shù)

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

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

圖2 尾礦壩穩(wěn)定性分析結(jié)果Fig.2 Stability analysis results of the tailings dam

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

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

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

3.2 尾礦壩穩(wěn)定可靠度分析

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

表2 尾礦材料參數(shù)統(tǒng)計特征Table 2 Statistical characteristics of tailings materials parameters

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

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

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

由圖6極限狀態(tài)函數(shù)的CDF曲線可以直接計算尾礦壩失穩(wěn)破壞概率。表3比較了不同方法計算的尾礦壩失穩(wěn)破壞概率。從表3可以看出，西木尾礦壩失穩(wěn)破壞概率水平很低，二階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次尾礦壩滲流有限元和穩(wěn)定性分析就可以滿足計算精度要求，計算量是直接LHS方法的1/4，可見本文提出的方法具有較高的計算效率，為解決考慮多個尾礦材料參數(shù)空間變異性的低概率水平尾礦壩穩(wěn)定可靠度問題提供了一條有效路徑。

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

4 結(jié)論

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

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

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