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

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

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

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

貝葉斯網(wǎng)絡(luò)(Bayesian Network,BN)是由美國加州大學(xué)PEARL[18]首次完整提出的一種不確定知識模型,該方法可以很好地彌補以上方法的不足。 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)用。 目前,貝葉斯分析先驗概率的確定仍是憑借主觀猜測[23],尋找先驗概率科學(xué)的求解方法對提高結(jié)果精度至關(guān)重要。 本研究通過GeNIe 軟件進(jìn)行貝葉斯分析,預(yù)測尾礦壩在不同降雨量下各空間維度的失穩(wěn)情況,分析不同降雨事件下尾礦內(nèi)部屬性的變動及相互作用機制,研究其失穩(wěn)機理,并對其穩(wěn)定性進(jìn)行評價分析。

1 模型構(gòu)建

1.1 基本原理

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

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

貝葉斯公式為

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

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

1.2 概率擬定

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

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

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

信息增益可由下式計算:

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

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

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

2 實例分析

2.1 工程概況

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

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

該尾礦庫所處區(qū)域氣候溫暖濕潤,屬于中亞熱帶季風(fēng)濕潤氣候,年降水強度1 563 mm,全年最大降雨強度為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 模型試驗

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

為了模擬潛在失穩(wěn)效果,在試驗初期設(shè)置了較高的尾礦庫庫區(qū)水位線,因此計算的失穩(wěn)概率會普遍偏高。 展示各個降雨事件下不同空間位置的尾礦參數(shù)的貝葉斯結(jié)果如圖8 所示。 小雨時,各級壩幾乎處于穩(wěn)定運行狀態(tài);中雨時,尾礦壩雖然穩(wěn)定概率較高,但與小雨相比,其穩(wěn)定運行的概率有所降低,尾礦壩有局部失穩(wěn)風(fēng)險且失穩(wěn)的概率略有增加;與小雨和中雨條件下的尾礦壩情況相比,大雨時其失穩(wěn)的概率增高,一級壩穩(wě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)險明顯加大。結(jié)合表2 分析原因:降雨量增加使得更多的雨水入滲到尾礦中,導(dǎo)致含水率增長和浸潤線抬升;含水率增長也將進(jìn)一步導(dǎo)致浸潤線上升;雨水入滲壩體,填充了尾礦骨架之間的孔隙,尾礦單元內(nèi)的水量增加,含水率較大。 水的潤滑作用使尾礦的黏聚力降低,抗滑力下降,壩體更易失穩(wěn)。

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

4 結(jié) 論

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

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

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

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