基于多源遙感數據的山西省交口縣某尾礦庫堆載過程回溯與失穩模式分析

周勝森, 李為樂, 陸會燕, 單云鋒, 郁文龍, 許善淼, 張 攀, 李維敏

(成都理工大學地質災害與地質環境保護國家重點實驗室, 四川 成都 610059)

0 引言

尾礦庫是礦山三大基礎工程之一,是筑壩攔截或圍地構成的、用以儲存金屬或非金屬礦山尾礦或其他工業廢渣的場所[1]。同時,尾礦庫也是金屬、非金屬礦山的重大危險源,一旦發生潰壩,將對下游居民生命財產安全和周邊環境造成嚴重威脅,因此尾礦庫在世界93種事故、公害中位列第18位[2],而造成百人以上死亡的尾礦庫事故已不鮮見[3]。例如,2019年1月25日,巴西東南部布魯馬迪紐市發生鐵礦石尾礦庫潰壩事故,造成165人死亡、160人失蹤,經濟損失達數千億[4]。據統計,我國目前仍有尾礦庫1.2萬座,約占世界尾礦庫數量的60%,幾乎遍布全國各地,且總體安全狀況不容樂觀[5],表1列出了近年來我國發生的重大尾礦庫潰壩事故。2008年山西襄汾新塔礦業公司“9·8”特別重大尾礦庫潰壩事故造成281人死亡,直接經濟損失達9 619.2萬元,造成了較為惡劣的社會影響,也將尾礦庫監測與安全評估工作提上戰略高度[6]。

傳統的尾礦庫監測主要是利用GNSS(Global Navigation Satellite System)獲取壩體位移數據,并結合浸潤線、干灘、降雨量及水位等數據實現災害預警[7]。GNSS能夠對尾礦庫壩體位移實現高精度、實時采集,但費時費力,監測范圍有限。近年來,得益于光學遙感衛星技術不斷發展,遙感影像時空分辨率不斷提高,為尾礦庫的遙感監測提供了有效的技術手段。有學者利用高分辨率國產光學衛星完成了大范圍尾礦庫遙感調查和環境影響分析[8-10]。2014年之后,歐空局免費開放了C波段Sentinel-1A雷達影像,重訪周期12天,且雷達影像對天氣的不敏感性使得周期性獲得雷達影像成為可能。相較于GNSS,合成孔徑雷達(Interferometric Synthetic Aperture Radar,InSAR)技術具有成本低、能夠大范圍監測等優勢,逐漸成為尾礦庫形變監測的有效手段。陳婭男等[11]利用SBAS-InSAR(Small Baseline Subsets InSAR)技術獲取了卡房尾礦壩的形變,分析表明該尾礦庫受人為工程和季節性降雨的影響較為嚴重。吳昊等[12]利用DS-InSAR(Distributed Scatterers InSAR)獲取了巴西 Brumadinho 尾礦庫時間序列地表形變信息,認為降雨對尾礦庫的變形影響顯著。

國內外對于尾礦庫的穩定性研究較為深入,表現為利用理論計算、物理實驗和數值模擬方法分析降雨[13]、地震[14]及凍融條件[15]下尾礦庫壩體力學性質和穩定性,以及潰壩災害發生后尾礦砂的運動特征[16],其中CFD-DEM、Fluent、MatDEM等數值模擬軟件在尾礦庫穩定性評價中的應用已日漸成熟[17-19]。在眾多數值模擬軟件中,GeoStudio具有操作簡便、界面交互式可視化等優點,Slope/W 模塊可廣泛用于巖土邊坡的穩定性分析。梁遠驥[20]應用GeoStudio軟件計算發現,廣東云浮某尾礦庫強降雨和連續降雨2種工況條件下尾礦壩體最小安全系數均大于尾礦庫安全規程明確的最小安全系數,表明該尾礦壩體穩定性滿足生產要求,不存在重大生產安全事故隱患。王峰等[21]利用GeoStudio軟件模擬分析得到遼寧某尾礦壩加高前后安全系數高于允許最小安全系數。

前人從光學遙感識別、InSAR 監測和數值模擬等角度對尾礦庫進行了研究,為尾礦庫的科學管理提供了思路,但研究手段較為單一,很少將遙感分析和穩定性計算相聯系。2022年3月17日17時,位于山西省呂梁市交口縣的某鋁礦尾礦庫發生潰壩事故,尾礦砂掩埋下方部分廠區并阻斷公路,造成較大經濟損失。尾礦庫潰壩災害發生后,在國家應急管理部的號召下,筆者所在團隊收集了多期高分辨率光學遙感衛星影像、立體像對數據和雷達衛星影像等多源遙感數據。首先,利用多源光學衛星影像和立體像的數據對溫泉鄉尾礦庫堆載過程進行回溯分析;然后,采用SBAS-InSAR技術持續監測尾礦庫潰壩前壩體的累積形變;最后,通過GeoStudio中Slope/W模塊對尾礦庫潰壩前不同時期的壩體邊坡進行穩定性評價。本文從遙感分析和數值模擬兩個角度分析該尾礦庫潰壩的主要影響因素,以期為類似尾礦庫工程的監測與安全評估提供支撐。

1 區域地質背景與尾礦庫概況

該尾礦庫位于山西省西部、呂梁市南端的交口縣溫泉鄉柴家莊村(111°25′49″ E,37°6′20″ N),地處黃土高原東緣、呂梁山脈中段,黃土覆蓋深厚,為典型的梁狀黃土丘陵地貌。研究區域屬于大陸性季風氣候,夏季降雨十分集中,占全年的70%以上,交口河流經尾礦庫東南方,為汾河的季節性支流[22-23]。研究區地處鄂爾多斯和華北地塊交界處,雖無大型斷裂發育,但地震活動較為頻繁,最近一次較大地震為2016年4月7日原平ML4.7地震[24-25]。尾礦庫周圍主要出露中奧陶統峰峰組(O2f)、上石炭統本溪組(C2b)和中更新統(Q2)等地層,其中上石炭統本溪組地層為頁巖和鐵鋁巖系,是重要的產鋁巖層[26](圖1)。

該尾礦庫于2017年12月左右開工建設,2018年投入使用,前緣高程約1 160 m,屬于山谷型尾礦庫,設計壩高50 m,庫容約91.2萬m3,為三級尾礦庫,用于堆存品質較低的鋁料。因下游有公路和廠房等設施,建設之初即形成“頭頂庫”。從潰壩后的光學影像上看(圖2),該尾礦庫礦砂顏色呈明顯暗紅色,潰壩后整體呈不規則形狀,潰壩范圍約7.2×105m2。壩體潰口位于尾礦壩中部,長約85 m,尾礦潰壩后掩埋了下方的公路和部分廠房。

2 尾礦庫潰壩遙感回溯分析

尾礦庫潰壩災害發生之后,為查明尾礦庫堆載過程、庫容變化和潰壩特征、壩體形變歷史,以及潰壩原因和機理,收集了潰壩前后各平臺多時相遙感資料和研究成果,對尾礦庫潰壩進行了多源遙感回溯分析。本文技術路線如圖3所示。

圖3 本研究采用的技術路線圖Fig.3 The technical flowchart adopted in this study

2.1 數據與方法

2.1.1 衛星影像數據

(1) 光學衛星影像

基于高分辨率光學衛星影像,從庫面紋理、幾何形狀可以解譯尾礦庫的堆載過程、潰壩形態等。本次收集到2017年5月6日—2021年5月27日期間共5期尾礦庫潰壩前的Google Earth平臺光學衛星影像,以及尾礦庫潰壩前2021年12月22日和潰壩后2022年4月6日高分一號(GF-1)衛星影像。具體參數如表2所列。

表2 光學衛星影像參數

(2) 衛星立體像對數據

衛星立體像對是兩個攝影基站對同一地物拍攝的兩張具有不少于60%重疊度的影像,通過對立體像對進行相對定向、絕對定向、核線影像生成、空三測量可構建測區空間坐標和數字高程模型(Digital Elevation Model,DEM)[27]。本次收集了尾礦庫堆載前的ALOS DEM (12.5 m分辨率),以及開始堆載后的2019年12月29日美國WorldView-1立體像對(0.5 m分辨率)、2021年12月22日和2022年4月6日國產高分7號立體像對(0.65 m分辨率)。

(3) 雷達衛星影像

為對尾礦庫壩體歷史形變做進一步分析,本次研究選擇歐空局的哨兵一號(Sentinel-1A)雷達衛星影像作為數據源,利用時序InSAR技術對尾礦庫潰壩前壩體的變形進行監測。經查詢,2018年6月1日—2022年3月17日共有109景升軌影像覆蓋尾礦庫區域,其影像距離向分辨率為5 m,方位向分辨率為20 m。

2.1.2 形變監測方法

為了獲取尾礦庫潰壩前的時序形變,利用SBAS-InSAR技術對壩體變形情況進行監測。SBAS-InSAR技術通過選擇合適的空間基線和時間基線閾值組成差分干涉對,并且選取相干目標點利用線性相位變化模型進行建模和解算,并通過時空濾波去除大氣延遲,獲取地表的形變時間序列[28-29]。其主要技術流程包括:

(1) 獲取N+1幅覆蓋相同區域的雷達影像,時間序列為[t0,t1,…,tn]T,任選其中一幅為主影像進行配準,并設置合適的時間和空間基線閾值,在所有差分干涉對中選取符合時空閾值的干涉對,假設得到M幅差分干涉圖,則滿足

(1)

式中:N為雷達影像數量;M干涉圖數量。

(2) 以t0時刻為參考時刻,則任意時刻ti(i=1,2,…,N)相對于t0時刻的差分相位φ(ti)為未知數,數據處理過程中所獲取的差分干涉相位δφ(tk)為觀測量。去除失相關、高程誤差以及大氣延遲等因素,則第i幅差分干涉圖中像元的相位值為

(2)

式中:λ為雷達波長;d(tA)和d(tB)分別為像元在時間tA和tB沿雷達視線方向的形變。

(3) 利用奇異值分解法估算出形變速率后,將平均形變速率在時間域上積分,即可得到時序形變。

2.2 堆載過程回溯

利用光學衛星影像進行目視解譯發現,該尾礦庫始建于2017年5月6日—2017年12月4日期間。2017年5月6日,該區域僅為一小型沖溝,植被覆蓋較高,還未開始堆放尾礦[圖4(a)]。2017年12月4日,沖溝內已堆放了少量尾礦[圖4(b)]。2018年4月16日,溝口第一道壩(1#壩)已有一定高度,壩體后方傾倒大量尾礦,并開始在庫區中部修筑第二道壩(2#壩)[圖4(c)]。2019年9月4日,尾礦庫堆積范圍接近失穩前的范圍,壩體接近最大高度[圖4(d)]。2021年5月27日,1#壩與2#壩之間的子庫已停止排放尾礦,并進行了場地平整,2#壩后方子庫較2019年尾礦堆載體積增加不明顯[圖4(e)]。2021年12月22日,2#壩后方子庫較2021年5月又有大量新堆載尾礦[圖4(f)]。2022年4月6日,2#壩中間形成了長約85 m的潰口,泄流尾礦將第一道壩沖毀后,翻越溝口山梁,掩埋了溝口道路和廠房,形成災害(圖2)。

2.3 庫容變化分析

利用ArcGIS軟件對多期DEM進行差分計算(圖5,圖6),發現2019年尾礦庫堆積方量約83萬m3,平均堆積厚度約13 m,最大堆積厚度位于2#壩后方,約33 m。2021年,尾礦庫堆積方量增加至96萬m3(大于設計庫容91.2萬m3),平均堆積厚度約13.8 m,最大堆積厚度約36 m。2019—2021年期間,新增堆積區主要位于2#壩以上子庫區,新增加方量約13萬m3,平均新增加堆積厚度約3.7 m,新增最大堆積厚度約11 m。

圖6 尾礦庫1-1′剖面圖Fig.6 1-1′ profile of the tailings pond

尾礦庫失穩物質損失區主要位于2#壩以上庫區,平均失穩深度約10 m,起動方量約13萬m3,與2019年以后新增尾礦體積相當。潰決起動后,碎屑流在2#壩與1#壩之間以及1#壩到溝口均主要以堆積為主,平均堆積厚度約4 m,總堆積方量約14萬m3。

2.4 壩體形變監測

通過SBAS-InSAR技術對雷達影像的處理,得到2018年6月至潰壩前尾礦庫區域的時序變形場,由于受不定期尾礦傾倒影響,尾礦堆積區域內InSAR監測可靠性相對較差,不易準確判識監測信號的成因,因此該區域數據做了掩膜處理,得到壩體、邊坡和周圍環境的形變場,如圖7所示。可以看出,庫區與周圍坡體存在明顯形變量差異,主要的形變區域位于中部的區域。分別取位于1#壩體和2#壩體上的3個特征點得到累計形變曲線,其中A點位于2#壩左側頂部、B點位于2#壩左側低部、C點位于2#壩右側底部,D、E、F點均位于1#壩。可以看出各點形變曲線總體呈現周期性下降趨勢,其中C點形變量最大,達到110 mm。

圖7 尾礦庫壩形變監測結果圖Fig.7 Deformation monitoring results of the tailing dam

圖8 2018年6月—2019年9月特征點形變曲線Fig.8 Deformation curves of characteristic points from June 2018 to September 2019

圖9 2019年9月—2022年3月特征點形變曲線Fig.9 Deformation curves of characteristic points from September 2019 to March 2022

2019年9月,尾礦庫堆積范圍已經接近失穩前的范圍,因此將累計形變曲線分為2018年6月—2019年9月、2019年9月—2022年3月兩個階段,分別對這兩個階段進行分析:(1)2018年6月—2019年9月時間段:A～F號點的形變信號不明顯,累積形變量約為30 mm,形變速率約18 mm/a;

(2)2019年9月—2022年3月時間段:A～F號點均發生較大形變,呈較強的增長趨勢。A、B、C三點位形變量最大,在2019-09—2022-03時間中達到80 mm左右的形變量。

3 2019年與2021年尾礦庫穩定性對比分析

目前,尾礦壩穩定性分析計算還是將其當作邊坡來處理,一般還是沿用土力學的傳統理論進行分析。本次穩定性分析采用二維有限元計算軟件GeoStudio中的Slope/W模塊,利用畢肖普極限平衡方法(Bishop)分別探測出2019年和2021年尾礦庫邊坡的最危險滑面,并求得安全系數。簡化的畢肖普安全系數Fs公式為:

(3)

3.1 幾何模型

根據2019年12月29日(以下簡寫為2019年)和2021年12月22日(以下簡寫為2021年)DEM尾礦庫堆載過程中的1-1′剖面數據,建立用于穩定性計算的幾何模型。該模型寬約610 m,高約95 m,2#壩體高約40 m。2019年,2#壩后方子庫最大堆積厚度約33 m[圖10(a)]。2021年,尾礦庫在2#壩后部進行了新的堆載,因此幾何模型坡面高程相比于2019年的幾何模型有所提升,最多處增加約11 m,最大堆積厚度約36 m[圖10(b)]。壩體力學參數與庫內沉積尾礦有顯著差異,故在模型中將子壩與庫內尾礦進行了區分。且為了簡化計算模型,未對尾礦砂進行分級。確定尾礦庫材料分為基巖、尾礦砂和壩體3種類型7個區域,如圖10所示。

圖10 計算模型材料區域劃分Fig.10 Material area division of the calculation models

3.2 材料參數

該尾礦庫下伏地層為奧陶統峰峰組(O2f)灰巖,尾礦庫為壓實的廢棄礦石和礦砂,庫內基本為自然條件下沉積的較為干燥的尾礦砂。為保證基巖強度足夠高,將其設置為“高強度”材料。尾礦砂和壩體物理力學性質變化不大,采用Mohr-Coulomb材料模型。通過文獻[30]和類似材料的參數取值,確定了本次數值模擬中材料的容重、黏聚力和內摩擦角等靜力學參數(表3)。

表3 材料靜力學參數

3.3 穩定性計算結果

2019年和2021年尾礦庫穩定性計算結果如圖11所示。2019年,2#壩體的安全系數為1.125,處于基本穩定狀態。2021年,尾礦庫2#壩后方子庫堆載量增加,2#壩體安全系數降至0.991,壩體已處于不穩定狀態。可以看出,在不考慮孔隙水壓力的條件下,礦砂堆載會導致壩體穩定性明顯降低,從而造成潰壩失穩。

圖11 最危險滑面和靜力安全系數Fig.11 The most dangerous sliding surface and static safety factor

4 討論

交口縣某尾礦庫潰壩原因主要為降雨和超載引起的尾礦壩力學強度降低,其中超載是誘發因素,降雨為影響因素。

4.1 子庫超載

尾礦庫不斷堆載,礦砂荷載持續增大,作用在壩體上的水平附加應力也不斷增大,導致壩體發生水平向變形。同時,快速堆載導致尾礦庫內水體不能排出,形成超孔隙水壓力,減小了壩體潛在滑移面顆粒間的有效應力,從而導致尾礦壩抗剪強度降低。當礦砂荷載超過壩體最大抗剪強度時,不均勻的形變將引發壩體產生剪切破壞。從光學衛星歷史影像來看,2019年5月尾礦庫已達到其最大庫容量,此時尾礦壩安全系數為1.125,處于基本穩定狀態。2021年12月2#壩后部子庫繼續堆載,尾礦壩安全系數降至0.991,子庫超載是該尾礦庫潰壩的誘發因素。

4.2 降雨作用

降雨入滲增大了壩體自重應力,壩體下滑力增大。同時,壩體內部孔隙水壓力增加,材料力學強度參數減小,壩體抗滑力降低。不能滲入壩體的雨水在坡表形成徑流,侵蝕壩體。如果發生短時強降雨,易形成洪水,導致漫頂發生。通過InSAR時序形變曲線與月降雨量數據(呂梁站:111.1°E,37.5°N)對比,尾礦庫壩體形變在雨季出現明顯加速現象,可以看出降雨對尾礦庫穩定性有較大影響(圖12)。但2022年3月交口縣累計降雨量僅有7.11 mm,由此看出,降雨是本次尾礦庫潰壩的影響因素而非誘發因素。

圖12 特征點形變與月降雨量關系Fig.12 Relationship between the deformation of characteristic point and monthly rainfall

4.3 尾礦庫失穩模式

尾礦庫失穩的誘發因素不同,失穩過程和模式也往往有較大差異。例如持續降雨和尾礦加載造成浸潤線埋深減小,尾礦壩體內發生管涌,造成尾礦庫內部裂縫和局部坍塌,同時壩坡表面發生流土,導致壩坡侵蝕,模式為滲流場變化直接誘發尾礦壩失穩。壩基巖體破壞造成壩基內形成滑移面,整個壩體發生深層滑移,模式為壩基失穩導致尾礦庫潰壩。短時強降雨和洪水造成庫區水位短時上升,尾礦砂被水流大量帶走,模式為漫頂導致尾礦壩坡失穩。地震荷載作用下,壩體發生局部液化;同時壩體出現較多裂縫,形成滲漏通道,降低壩體穩定性,模式為地震作用下導致尾礦庫潰壩[2]。根據尾礦庫潰壩發生條件和失穩過程對比,交口縣某尾礦庫潰壩屬于典型的滲流場變化誘發尾礦壩失穩,誘發因素為尾礦庫超載。

5 結論

本文利用多時相光學衛星立體像對數據對山西交口縣某尾礦庫堆載過程及庫容變化進行了回溯分析,通過SBAS-InSAR技術監測分析了尾礦庫潰壩前的壩體和庫區邊坡的微小形變,最后采用Geo-Studio軟件Slope/W模塊評價了潰壩前不同時間邊坡的穩定性。文中通過遙感分析和數值模擬兩個角度分析了該尾礦庫潰壩的主要影響因素,主要得到以下結論:

(1) 該尾礦庫壩體形變量隨著堆載量的增加逐漸增大,庫容于2019年9月前已堆積至最大庫容,2021年12月前在2#壩后方子庫堆載至96萬m3,超出設計庫容4.8萬m3,從而降低了壩體的穩定性;

(2) 潰壩物源主要來自2#壩以上庫區,在2#壩與1#壩之間以及1#壩到溝口區域均主要以堆積為主;

(3) 該尾礦庫潰壩屬于典型的滲流場誘發尾礦壩失穩,尾礦加載導致其壩體邊坡穩定性安全系數由1.125降至0.991,是該尾礦庫潰壩的直接誘發因素。

本文仍存在以下不足之處:首先,通過SBAS-InSAR技術得到的形變是雷達衛星視線方向的投影,能判斷形變趨勢,但與真實形變有一定差異。其次,在建立尾礦庫幾何模型時,為了簡化計算模型,對尾礦砂未進行分級,計算結果可能存在一定偏差。