一種基于靜態(tài)貝葉斯的降雨誘發(fā)尾礦壩失穩(wěn)評估模型

胡靖雯 聶 聞 蘆 松 王震豪 Pooya Saffari

(1.中北大學(xué)理學(xué)院,山西 太原 030051;2.中國科學(xué)院福建物質(zhì)結(jié)構(gòu)研究所,福建 福州 350000;3.金屬礦山安全與健康國家重點(diǎn)實(shí)驗(yàn)室,安徽 馬鞍山 243000;4.青島城市學(xué)院土木工程學(xué)院,山東 青島 266071)

尾礦壩是一種用于堆放礦山廢棄物的形似堤壩的構(gòu)筑物,其中可能含有高濃度的可溶性重金屬化合物。 作為一個(gè)高勢能的危險(xiǎn)源,其對礦山安全生產(chǎn)、生態(tài)環(huán)境和人民生命財(cái)產(chǎn)安全造成了巨大的威脅[1]。 我國現(xiàn)存尾礦壩12 655 座,具有數(shù)量大、新增多、高危壩分布多且廣的特點(diǎn)[2]。 在所有尾礦壩失穩(wěn)的致災(zāi)因素中,降雨是導(dǎo)致其失穩(wěn)的主要因素之一,占據(jù)了尾礦壩失穩(wěn)事故總量的24%[3]。 例如2010 年9 月21 日,強(qiáng)降雨導(dǎo)致了廣東省銀巖錫礦尾礦壩倒塌,造成22 人死亡,523 間房屋倒塌,直接經(jīng)濟(jì)損失達(dá)1 900 萬元[4]。 2010 年10 月4 日,由于異常降雨,匈牙利Ajka 赤泥大壩發(fā)生災(zāi)難性坍塌,波及下游1 017 ha 的土地,致使367 處不動(dòng)產(chǎn)受損,多人受傷[5]。 鑒于國內(nèi)外尾礦壩潰壩的慘痛教訓(xùn),我國對其安全運(yùn)行愈發(fā)重視,對尾礦壩失穩(wěn)的機(jī)理研究和預(yù)測愈發(fā)重要。

近年來,不少學(xué)者對尾礦壩的潰壩機(jī)制展開研究,為尾礦壩的安全管理提供了理論支撐。 宋志飛等[6]明確了尾礦壩子壩高度是尾礦沉積規(guī)律的主要外部影響因素。 YIN 等[7]在剖面模型試驗(yàn)的基礎(chǔ)上,對尾礦壩不同高度、不同工況下的穩(wěn)定性進(jìn)行了分析,并考慮了尾礦壩堆積材料多層次的復(fù)雜特征。 部分學(xué)者運(yùn)用數(shù)值仿真模擬方法預(yù)測潰壩,繆海波等[8]、侯永莉等[9]通過GEO-Studio 軟件進(jìn)行滲流模擬,李海港等[10]通過數(shù)值模擬試驗(yàn),分析了尾礦壩的潰壩演化規(guī)律,然而數(shù)值仿真的邊界條件理想化易造成結(jié)果失真。 隨著機(jī)器學(xué)習(xí)算法的發(fā)展,學(xué)者們將其應(yīng)用到尾礦壩領(lǐng)域:NIE 等[11]通過建立改良的水箱模型預(yù)測孔隙水壓變化;DU 等[12]運(yùn)用InSAR 時(shí)間序列對尾礦壩進(jìn)行風(fēng)險(xiǎn)評估;王飛躍等[13]擬合浸潤線觀測孔水位與庫水位的函數(shù)曲線,求得尾礦壩浸潤線矩陣用于進(jìn)行穩(wěn)定性分析;李輝等[14]通過層次分析法和信息熵理論計(jì)算主客觀權(quán)重,得到組合因子的辨識(shí)向量,根據(jù)置信度準(zhǔn)則對尾礦壩的穩(wěn)定性進(jìn)行了評價(jià)。 由于尾礦壩失穩(wěn)時(shí)其內(nèi)部屬性相互作用的復(fù)雜性與不確定性,傳統(tǒng)的分析方法如故障樹[15]、層次分析法[16]、事故樹法[17]在解決問題的不確定性方面略顯不足。 上述方法均無法有效解決各級事件之間的關(guān)聯(lián)性,而分析各級事件之間的相關(guān)性在實(shí)際工程中十分必要。

貝葉斯網(wǎng)絡(luò)(Bayesian Network,BN)是由美國加州大學(xué)PEARL[18]首次完整提出的一種不確定知識(shí)模型,該方法可以很好地彌補(bǔ)以上方法的不足。 BN作為一種概率模型,提供了一種直觀方式表達(dá)大量相互關(guān)聯(lián)變量的聯(lián)合分布,變量代表了可能發(fā)生故障的原因、不同的故障模式及可能的后果[19]。 貝葉斯網(wǎng)絡(luò)通過更新給定的觀測值分布,以一致的概率方式表示不確定的結(jié)果,在電力系統(tǒng)評估[20]、房地產(chǎn)市場[21]、煤與瓦斯突出[22]等領(lǐng)域得到應(yīng)用。 目前,貝葉斯分析先驗(yàn)概率的確定仍是憑借主觀猜測[23],尋找先驗(yàn)概率科學(xué)的求解方法對提高結(jié)果精度至關(guān)重要。 本研究通過GeNIe 軟件進(jìn)行貝葉斯分析,預(yù)測尾礦壩在不同降雨量下各空間維度的失穩(wěn)情況,分析不同降雨事件下尾礦內(nèi)部屬性的變動(dòng)及相互作用機(jī)制,研究其失穩(wěn)機(jī)理,并對其穩(wěn)定性進(jìn)行評價(jià)分析。

1 模型構(gòu)建

1.1 基本原理

貝葉斯網(wǎng)絡(luò)是一種非循環(huán)有向圖及相關(guān)參數(shù)的集合,可通過邏輯推理解決不確定性問題[24]。 BN 由模型結(jié)構(gòu)和相關(guān)參數(shù)組成。 圖1 中每個(gè)節(jié)點(diǎn)對應(yīng)一個(gè)系統(tǒng)變量,箭頭代表直接的定性概率依賴,節(jié)點(diǎn)X1和X2是節(jié)點(diǎn)X3的父節(jié)點(diǎn),節(jié)點(diǎn)X3是X1和X2的子節(jié)點(diǎn)。

圖1 BN 網(wǎng)絡(luò)示意Fig.1 Schematic of BN network

貝葉斯公式為

式中,P(X1)為標(biāo)準(zhǔn)化常量;P(X2)為先驗(yàn)概率;P(X1|X2)為似然率;P(X2|X)1為后驗(yàn)概率。

先驗(yàn)概率P(X2)反映的是在觀測數(shù)據(jù)之前對其的認(rèn)知[25]。 似然率PX1|X2()反映了在固定參數(shù)的情況下得到的某觀測值的可信程度。 后驗(yàn)概率P(X2|X)1是貝葉斯網(wǎng)絡(luò)預(yù)測結(jié)果,反映了在固定參數(shù)和模型的情況下對事件的知識(shí)。 標(biāo)注化常量P(X1)可視為歸一化系數(shù)。

1.2 概率擬定

由于BN 節(jié)點(diǎn)是布爾屬性(Bool)的特點(diǎn),即只有0、1 兩種狀態(tài)[26],故需將數(shù)據(jù)轉(zhuǎn)譯成概率代入BN 網(wǎng)絡(luò)。 通過信息增益[27]選取特征在決策樹分類算法中使用廣泛,技術(shù)較為成熟,其通過計(jì)算所有特征值化分?jǐn)?shù)據(jù)集得到,信息增益最高的特征值為最優(yōu)值。

熵就是信息的期望值,即數(shù)據(jù)集的無序度,其值越小,變量的不確定性就越小,可進(jìn)行如下計(jì)算:

式中,Hs為數(shù)據(jù)集的熵值;n為分類數(shù)目;pi為選擇分類i的概率。

信息增益可由下式計(jì)算:

式中,D為數(shù)據(jù)集;A為特征;g(D,A) 為特征A對訓(xùn)練數(shù)據(jù)集D的信息增益。

現(xiàn)成的數(shù)據(jù)集可以較大程度上減少參數(shù)化BN所需的工作。 但是當(dāng)數(shù)據(jù)集很小時(shí),許多條件案例由于沒有數(shù)據(jù)記錄,沒有足夠的數(shù)據(jù)基礎(chǔ)支撐條件概率的學(xué)習(xí)。 針對這一問題,ONISKO 等[28]提出了Noisy-OR 來減少條件概率數(shù)據(jù)需求的方法,并給出了Leaky Noisy-or gate 模型。 節(jié)點(diǎn)X的條件概率可表示為

式中,Xa為事件X的父節(jié)點(diǎn)之一;Pa為事件X與Xa的連接概率;Xc為所有未知因素的綜合;Pc為事件X與Xc的連接概率;Xp為事件X的父節(jié)點(diǎn)的集合。

2 實(shí)例分析

2.1 工程概況

落木坑尾礦庫地理坐標(biāo)為東經(jīng)114° 21′16″,北緯25°23′08″,根據(jù)《選礦廠尾礦設(shè)施設(shè)計(jì)規(guī)范》(ZB J1-90),該尾礦庫為二等庫,其基本參數(shù)取值見表1。

表1 落木坑尾礦庫概況Table 1 General situation of Luomukeng tailings pond

該尾礦庫所處區(qū)域氣候溫暖濕潤,屬于中亞熱帶季風(fēng)濕潤氣候,年降水強(qiáng)度1 563 mm,全年最大降雨強(qiáng)度為61.7 mm/h,日降雨量最大為121 mm,年蒸發(fā)量為1 445.47～1 846.8 mm,降雨量及溫度如圖2 所示。 尾礦壩所處區(qū)域降雨量較大,雨水對其穩(wěn)定性的影響不容忽視。

圖2 大余縣2020 年降雨量及氣溫Fig.2 Rainfall intensity and temperature in Dayu County in 2020

2.2 模型試驗(yàn)

本研究采用物理堆壩模型試驗(yàn)重建尾礦壩的堆積過程[29],落木坑尾礦庫的鳥瞰圖如圖 3(a)所示。根據(jù)尾礦庫現(xiàn)場的情況和實(shí)際研究的可行性,確定模型比尺為1 ∶100。 模型長約13 m,寬7 m,高1.5 m,傾角約35°。 尾礦壩由1 個(gè)初級壩和4 個(gè)子壩組成。建壩初期,采用自下而上遞增的方法,修建了4 級子壩。 初級壩采用紅黏土建造,4 個(gè)子壩的材料是從落木坑尾礦中收集的尾礦。 試驗(yàn)?zāi)Ｐ腿鐖D3(b)所示。

圖3 落木坑尾礦庫Fig.3 Luomukeng tailings dam

在堆壩模型試驗(yàn)的基礎(chǔ)上,采用專用的降雨模擬裝置,再現(xiàn)降雨過程對尾礦壩的影響。 本研究尾礦壩模型試驗(yàn)系統(tǒng)主要由4 部分組成,即山地模型、放礦系統(tǒng)、供水系統(tǒng)、降雨系統(tǒng),如圖4 所示。

圖4 試驗(yàn)?zāi)Ｐ驼掌現(xiàn)ig.4 Photos of the experimental model

為了模擬大壩尾礦的堆積過程,將尾礦浸泡在少量的水中,并放置在一個(gè)大塑料桶中。 尾礦沿塑料桶底部直徑8 cm 的PVC 管注入尾礦池。 降雨設(shè)備由6組相同直徑的噴嘴組成(每組6 個(gè)噴嘴),降雨強(qiáng)度可以在0 到100 mm/h 之間變化。 參考中國氣象局頒布的降雨強(qiáng)度等級,結(jié)合試驗(yàn)實(shí)際情況劃分標(biāo)準(zhǔn),試驗(yàn)的降雨強(qiáng)度分為15、25、35 mm/h 3 個(gè)等級。 降雨設(shè)定持續(xù)時(shí)間為10 min,關(guān)閉降雨設(shè)施后,每級壩20 cm 深處隨機(jī)取樣測量實(shí)時(shí)含水率。 取樣方法為每級壩隨機(jī)取3 個(gè)樣品,共12 個(gè)樣品,樣品質(zhì)量為3 kg,取樣位置如圖3(b)所示,將樣品裝入塑料袋中進(jìn)行后續(xù)的試驗(yàn)測試。 取樣后,取樣產(chǎn)生的孔洞用相同的材料填充。 等待模型排水完畢,進(jìn)入下一個(gè)降雨模式,重復(fù)前面的步驟,直到完成3 組降雨測試試驗(yàn)。分別在每級壩體內(nèi)部埋設(shè)4 個(gè)孔隙水壓計(jì),位置沿壩體中軸線剖面布設(shè),以初級壩為基準(zhǔn)點(diǎn),計(jì)算獲得不同降雨量下不同區(qū)域的尾礦初始浸潤線高度。

根據(jù)《試驗(yàn)室玻璃器皿密度計(jì)》(GB/T 21785—2008)進(jìn)行了比重計(jì)試驗(yàn),試驗(yàn)儀器為TM85 尾礦密度儀。 根據(jù)《煤炭篩分試驗(yàn)方法》(GB/T 477—2008)進(jìn)行篩分試驗(yàn),裝置為GZS-1 型高頻振動(dòng)篩。 大雨時(shí)部分尾礦的級配曲線如圖5 所示。 由圖5 可知:不同高度尾礦的材料非均質(zhì)特性明顯,各級壩的級配曲線在尾礦粒徑為0.002～0.600 mm 時(shí)差異較大。 尾礦位置越高,顆粒級配曲線越陡,尾礦的粒徑比較均勻。

圖5 大雨時(shí)尾礦的粒徑級配累積曲線Fig.5 Cumulative curve of particle size gradation of tailings during heavy rain

根據(jù)《土壤檢測第23 部分:土粒密度的測定》(NY/T 1121.23—2010),利用環(huán)刀法測量尾礦的密度。 直剪試驗(yàn)依據(jù)的規(guī)范為《土工試驗(yàn)方法標(biāo)準(zhǔn)》(GB/T 50123—2019),采用的設(shè)備為ZJ 應(yīng)變控制直剪儀,通過直剪試驗(yàn)獲得樣品的黏聚力c和摩擦角φ。 限于篇幅,部分參數(shù)取值見表2。 由表可知:降雨量越大,浸潤線位置越高;高處的尾礦黏聚力較小,摩擦角略小,含水率較大,密度較低。

表2 部分參數(shù)取值Table 2 Values of part parameters

3 基于靜態(tài)貝葉斯的尾礦壩穩(wěn)定性分析

基于尾礦庫各參數(shù)自身的性質(zhì)將其劃分為兩類:物理參數(shù),包括浸潤線、密度、含水率;力學(xué)參數(shù),包括黏聚力、摩擦角、抗滑力和滑動(dòng)力。 通過查閱相關(guān)資料[30],基于各參數(shù)的邏輯關(guān)系,建立了如圖6 所示的模型。 為了豐富樣本數(shù)據(jù),提高信息增益的準(zhǔn)確性,針對每個(gè)節(jié)點(diǎn)進(jìn)行拉格朗日多項(xiàng)式插值[31]。 通過選取信息增益最高值劃分?jǐn)?shù)據(jù)集,并將其作為預(yù)警值判斷每個(gè)節(jié)點(diǎn)的狀態(tài),據(jù)此計(jì)算先驗(yàn)概率。 由于貝葉斯模型以安全因子為標(biāo)準(zhǔn)評判壩體是否失事,故選擇該值作為評價(jià)指標(biāo)劃分各個(gè)參數(shù)區(qū)間。 以密度為例說明求取過程,當(dāng)安全因子為1 時(shí)密度值就是該參數(shù)的預(yù)警值。 按照式(4)計(jì)算信息增益,密度最小值為2.32,為了保證結(jié)果精度,從2. 32 開始,每次增加10-6,直至計(jì)算到最大值2.99。 當(dāng)信息增益值最大達(dá)到2.750 001 時(shí),該值即為密度的區(qū)間劃分值。 計(jì)算出密度大于該值的概率占總體的31.42%,該值即為尾礦失穩(wěn)時(shí)密度的先驗(yàn)概率。 各個(gè)節(jié)點(diǎn)狀態(tài)如圖6所示,箱子下方左右兩側(cè)數(shù)據(jù)分別表示數(shù)據(jù)的最大和最小值,中間數(shù)據(jù)表示該節(jié)點(diǎn)的預(yù)警值,箱中百分比為先驗(yàn)概率。

圖6 節(jié)點(diǎn)狀態(tài)Fig.6 Node status

當(dāng)節(jié)點(diǎn)值大于預(yù)警值時(shí),即Xi=1,節(jié)點(diǎn)處于危險(xiǎn)狀態(tài);反之,Xi=0,節(jié)點(diǎn)處于安全狀態(tài)。 將壩體穩(wěn)定即“壩體失穩(wěn)”節(jié)點(diǎn)的狀態(tài)為“0”作為標(biāo)準(zhǔn)學(xué)習(xí)各個(gè)節(jié)點(diǎn)的狀態(tài)。 以抗滑力為例說明建模過程,設(shè)含水率為事件X,設(shè)其父節(jié)點(diǎn)降雨,高度分別為X1,X2。 基于模型試驗(yàn),將降雨分為小中大3 個(gè)等級,分別為X1=1,X1=2,X1=3,假設(shè)發(fā)生各降雨量的概率相等,每個(gè)等級的降雨量發(fā)生概率為33. 33%。 尾礦壩為四級壩,分別為X2=1,X2=2,X2=3,X2=4,將每級壩的概率設(shè)為25%。 統(tǒng)計(jì)X=1 時(shí)各個(gè)節(jié)點(diǎn)的狀態(tài),計(jì)算的概率參數(shù)見表3。

表3 含水率的概率參數(shù)Table 3 Probability parameters of the water content

將含水率的概率參數(shù)代入式(5)中,計(jì)算含水率的條件概率。 含水率在一二級壩時(shí)的條件概率見表4。

表4 含水率的條件概率Table 4 Conditional probability of water content

根據(jù)《尾礦庫安全規(guī)程》(GB 39496—2020),將尾礦的安全系數(shù)劃分為4 個(gè)等級(表5)。

表5 尾礦安全系數(shù)及其狀態(tài)Table 5 Safety factor and status of tailings

將所有節(jié)點(diǎn)的條件概率輸入建立的網(wǎng)絡(luò)中,各節(jié)點(diǎn)的右下角會(huì)出現(xiàn)“√”腳標(biāo),得到如圖7 所示的模型,State 0 表示該節(jié)點(diǎn)仍處于安全狀態(tài),State 1 表示該節(jié)點(diǎn)處于危險(xiǎn)狀態(tài)。 本研究尾礦壩表層發(fā)生失穩(wěn),滑動(dòng)力、抗滑力及安全系數(shù)的計(jì)算采用無限邊坡模型[32]:

圖7 貝葉斯條形圖Fig.7 Bar graph of Bayesian

式中,β為潛在的滑動(dòng)面傾角,(°);ρ為土體平均密度,g/cm3;H為土面至破壞面的深度,m;c′為潛在滑動(dòng)面的有效土壤黏聚力,N/m2;φ′為潛在滑動(dòng)面的有效摩擦角,(°);PWP為潛在滑動(dòng)面的孔隙水壓力,Pa。

為了模擬潛在失穩(wěn)效果,在試驗(yàn)初期設(shè)置了較高的尾礦庫庫區(qū)水位線,因此計(jì)算的失穩(wěn)概率會(huì)普遍偏高。 展示各個(gè)降雨事件下不同空間位置的尾礦參數(shù)的貝葉斯結(jié)果如圖8 所示。 小雨時(shí),各級壩幾乎處于穩(wěn)定運(yùn)行狀態(tài);中雨時(shí),尾礦壩雖然穩(wěn)定概率較高,但與小雨相比,其穩(wěn)定運(yùn)行的概率有所降低,尾礦壩有局部失穩(wěn)風(fēng)險(xiǎn)且失穩(wěn)的概率略有增加;與小雨和中雨條件下的尾礦壩情況相比,大雨時(shí)其失穩(wěn)的概率增高,一級壩穩(wěn)定運(yùn)行的概率較大,其余子壩極有可能處于臨界失穩(wěn)狀態(tài)。 綜合分析得到如下結(jié)論:

圖8 尾礦壩失穩(wěn)分析結(jié)果Fig.8 Results of instability analysis of tailings dam

(1)降雨量增加導(dǎo)致尾礦壩失穩(wěn)風(fēng)險(xiǎn)明顯加大。結(jié)合表2 分析原因:降雨量增加使得更多的雨水入滲到尾礦中,導(dǎo)致含水率增長和浸潤線抬升;含水率增長也將進(jìn)一步導(dǎo)致浸潤線上升;雨水入滲壩體,填充了尾礦骨架之間的孔隙,尾礦單元內(nèi)的水量增加,含水率較大。 水的潤滑作用使尾礦的黏聚力降低,抗滑力下降,壩體更易失穩(wěn)。

(2)同等降雨條件下,尾礦位置越高,其失穩(wěn)的概率越大。 結(jié)合圖5 分析可得,不同高度尾礦的材料非均質(zhì)特性明顯,高處尾礦的尾礦粒徑均勻,缺乏小顆粒填充大顆粒之間形成的孔隙,顆粒之間距離比較遠(yuǎn),單位面積上土粒的接觸點(diǎn)少,顆粒之間的咬合作用較弱,黏聚力較小;雨水迸濺和徑流沖刷導(dǎo)致尾礦表面顆粒向下遷徙,加劇了不同空間位置尾礦的差異性。 高處顆粒缺失,密度較低,影響尾礦的摩擦角;因此,尾礦抗滑力降低,失穩(wěn)概率增加。

4 結(jié) 論

本研究以江西落木坑尾礦壩為例,通過1 ∶100 大比尺模型試驗(yàn)獲取不同降雨量下尾礦的浸潤線、黏聚力、摩擦角、含水率和密度的樣本數(shù)據(jù),并利用無限邊坡模型計(jì)算了滑動(dòng)力、抗滑力和安全系數(shù)。 結(jié)合信息增益和Leaky Noisy-or gate 模型,建立了一種新的基于貝葉斯的尾礦壩穩(wěn)定性評估模型,闡明了尾礦壩的失穩(wěn)機(jī)理,預(yù)測了不同降雨事件下不同空間位置的尾礦壩的狀態(tài)。 所得結(jié)論如下:

(1)相較于傳統(tǒng)的貝葉斯分析,通過信息增益計(jì)算先驗(yàn)概率大大降低了模型分析的主觀性。相較于其他不確定性分析方法,使用Leaky Noisy-or gate 模型求各節(jié)點(diǎn)之間的條件概率,降低了因節(jié)點(diǎn)信息缺失或失真造成的偏差,極大提高了模型計(jì)算精度。 通過GeNIe 軟件建模,實(shí)現(xiàn)了結(jié)構(gòu)和參數(shù)的自學(xué)習(xí),發(fā)揮了該軟件在貝葉斯推理方面的顯著優(yōu)勢,并可以通過結(jié)果進(jìn)行反推,追溯壩體破壞的原因。

(2)降雨量對尾礦參數(shù)變動(dòng)有著顯著影響,并且隨著降雨量增加,各參數(shù)之間的相互作用更加明顯。降雨量增加使更多雨水入滲到壩體中,導(dǎo)致尾礦含水率增長和浸潤線抬升,含水率增加又將導(dǎo)致浸潤線進(jìn)一步上升以及黏聚力降低。 由于尾礦壩空間非均質(zhì)特性,壩體高處粒徑均勻,本就容易失穩(wěn),而雨水迸濺和徑流沖刷造成的壩體侵蝕將差異變得更大,高處尾礦的失穩(wěn)風(fēng)險(xiǎn)進(jìn)一步加劇。 降雨通過改變尾礦的物理參數(shù)影響力學(xué)參數(shù),從而改變壩體穩(wěn)定性。

(3)本研究所建模型的分析結(jié)果與試驗(yàn)實(shí)測記錄相符,說明該模型有一定的可靠性,可適用其他工況下的坡體穩(wěn)定性分析。 不足之處是該模型沒有考慮時(shí)態(tài)作用下尾礦壩穩(wěn)定性各影響因素的變動(dòng),鑒于動(dòng)態(tài)貝葉斯在描述隨機(jī)演化不確定關(guān)系方面的優(yōu)勢,下一步考慮將貝葉斯網(wǎng)絡(luò)與時(shí)間信息相結(jié)合,形成一種可處理時(shí)序數(shù)據(jù)的隨機(jī)模型。