李如仁 孫加瑤

(沈陽建筑大學(xué)交通與測繪工程學(xué)院,遼寧 沈陽 110168)

尾礦庫作為一種具有人造高勢能泥石流危險(xiǎn)源的特殊工業(yè)建筑,一旦發(fā)生潰壩或泄露,會(huì)對社會(huì)安全造成巨大威脅,同時(shí)因尾礦中含有重金屬物質(zhì),也將嚴(yán)重污染環(huán)境[1]。據(jù)統(tǒng)計(jì),在2001—2015年全國尾礦庫共發(fā)生99起潰壩和泄漏事故[2]。2008年9月8日山西新塔礦業(yè)有限公司尾礦庫潰壩直接導(dǎo)致泥石流發(fā)生,造成下游居民277人死亡、4人失蹤,直接經(jīng)濟(jì)損失9 619.21萬元[3]。2015年11月5日,巴西淡水河谷公司尾礦庫發(fā)生潰壩,造成19人死亡,數(shù)千人失蹤[4]。因此對尾礦庫進(jìn)行長期的形變監(jiān)測及預(yù)警具有十分重要的意義。

目前,針對尾礦庫的監(jiān)測手段主要依賴于傳統(tǒng)的水準(zhǔn)測量、全球定位系統(tǒng)(Global Positioning System,GPS)、測量機(jī)器人等。通過該類技術(shù)手段只能實(shí)現(xiàn)對監(jiān)測站點(diǎn)繪制出沉降剖面曲線圖并進(jìn)行相關(guān)分析,如需對監(jiān)測區(qū)域進(jìn)行大范圍、全方位的分析則必須大量布設(shè)沉降監(jiān)測點(diǎn),并通過擬合、插值等技術(shù)手段繪制沉降圖,這種大范圍沉降監(jiān)測方式不僅存在工作量大、費(fèi)時(shí)、費(fèi)力等不足,還會(huì)由于測點(diǎn)難以保存產(chǎn)生更大的誤差[5]。近年來,合成孔徑雷達(dá)干涉測量(Interferometric Synthetic Apeture Rader,InSAR)技術(shù)作為一種新型的對地觀測技術(shù),可以根據(jù)衛(wèi)星的軌道信息與觀測點(diǎn)的幾何關(guān)系,測算出觀測點(diǎn)的空間信息和形變,理論上可以獲取高精度的數(shù)字高程模型(Digital Elevation Model,DEM)和毫米級(jí)的地表形變信息[6]。為獲取連續(xù)形變值,提高測量精度,有學(xué)者提出了時(shí)間序列InSAR技術(shù),很好地解決了時(shí)空失相關(guān)和大氣延遲等問題,例如永久散射體合成孔徑雷達(dá)干涉測量(Permanent Scatterer Synthetic Aperture Radar,PSInSAR)技術(shù)[7],小基線集合成孔徑雷達(dá)干涉測量(Small Baseline Subsets Interferometric Synthetic Aperture Radar,SBAS-InSAR)技術(shù)[8]等,目前這些方法被廣泛用于尾礦庫形變監(jiān)測。例如,王麗娟[9]利用DInSAR技術(shù)對萬年溝尾礦庫進(jìn)行監(jiān)測,通過實(shí)測數(shù)據(jù)證明了該方法的可靠性;陳婭男等[10]利用SBAS-In-SAR技術(shù)獲取了卡房尾礦庫的形變演化規(guī)律,研究發(fā)現(xiàn)卡房尾礦庫形變受人為工程和雨季影響較大;吳昊等[11]利用DS-InSAR技術(shù)對巴西Brumadinho尾礦庫進(jìn)行監(jiān)測,分析了其演化過程及潰壩原因。雖然InSAR技術(shù)在尾礦庫形變監(jiān)測方面的應(yīng)用研究已有較多案例,但是針對尾礦庫形變預(yù)測模型的研究涉及較少,在很大程度上影響了對尾礦庫災(zāi)害的實(shí)時(shí)預(yù)警。

目前,常用的預(yù)測方法主要有BP神經(jīng)網(wǎng)絡(luò)[12]、支持向量機(jī)回歸(Support Vector Regression,SVR)[13]、循環(huán)神經(jīng)網(wǎng)絡(luò)(Recurrent Neural Network,RNN)[14]等。然而上述預(yù)測模型存在過擬合、梯度消失、梯度爆炸[15]等問題,近年來通過RNN改進(jìn)的長短期記憶(Long Short-Term Memory,LSTM)神經(jīng)網(wǎng)絡(luò)在各方面時(shí)間序列預(yù)測中被廣泛應(yīng)用[16],并表現(xiàn)出良好的長期及短期預(yù)測能力。本研究在現(xiàn)有成果的基礎(chǔ)上融合SBAS-InSAR技術(shù)、LSTM神經(jīng)網(wǎng)絡(luò)和網(wǎng)格搜索(Grid Search,GS)算法,對鞍山市西果園尾礦庫進(jìn)行形變監(jiān)測及分析,總結(jié)形變規(guī)律,將所獲沉降結(jié)果作為學(xué)習(xí)訓(xùn)練樣本,建立預(yù)測動(dòng)態(tài)模型,實(shí)現(xiàn)監(jiān)測與預(yù)測的一體化,為尾礦庫災(zāi)害防治提供參考。

1 形變監(jiān)測與預(yù)測原理

1.1 SBAS-InSAR基本原理

SBAS-InSAR作為一種在InSAR基礎(chǔ)上發(fā)展的新型時(shí)序分析手段,極大地改善了常規(guī)InSAR技術(shù)的空間、時(shí)間失相關(guān)等問題,可以有效降低相位噪聲和誤差[17]。其基本思路是對空間基線與時(shí)間基線進(jìn)行限制,設(shè)置閾值,將SAR影像按規(guī)定組合生成若干個(gè)小基線集合從而獲取干涉對,使同一子集內(nèi)的基線相對較短,集合之間的基線相對較長。在各短基線子集內(nèi)部利用最小二乘法(Least Squares,LS)計(jì)算時(shí)間序列形變,各集合之間使用奇異值分解法(Singular Value Decomposition,SVD)聯(lián)合求解具有高精度的完整時(shí)間序列形變值。

假設(shè)在時(shí)間t0至?xí)r間tn內(nèi)獲取同一研究區(qū)的N幅SAR影像,選擇一幅影像作為主影像,然后根據(jù)所設(shè)組合閾值,形成M幅干涉圖,且滿足下式:

對于ta與tb(ta＜tb)時(shí)刻生成的第j幅干涉圖像,在去除平地效應(yīng)和地形相位的影響后,若方位向坐標(biāo)以x表示,距離向坐標(biāo)以y表示,則干涉圖j中任意一點(diǎn)的相位可表示為

式中,φtb(x,y)、φta(x,y)分別表示坐標(biāo)為(x,y)的像元在其對應(yīng)時(shí)刻的相位;λ為雷達(dá)波長;dtb、dta分別代表在tb和ta時(shí)刻該像元在雷達(dá)視線方向上的形變值;φtopo、φatm、φnoise分別為干涉圖中殘余的地形相位、大氣延遲相位、噪聲相位。

暫不考慮大氣延遲相位與噪聲相位,地形相位可通過高精度DEM消除。因而式(2)可以簡化為

因此,解纏后的所有差分干涉圖可以表示為

式中,φ為所有影像數(shù)據(jù)相位值所組成的向量;δφ為所有差分干涉相位值所組成的向量;A為M×N階矩陣。當(dāng)所有數(shù)據(jù)均屬于一個(gè)基線集時(shí),M≥N,矩陣A的秩為N,此時(shí)可以通過最小二乘法求出φ的估計(jì)值(AT·A)-1·AT·δφ。

當(dāng)數(shù)據(jù)集分散在幾個(gè)子集中時(shí),矩陣A的秩小于N,為秩虧矩陣,方程組有無窮多組解。為求得唯一解,SBAS-InSAR技術(shù)采用奇異值分解法,將多個(gè)基線集聯(lián)合求解,求出最小范數(shù)意義上的最小二乘解,從而獲取累計(jì)形變量。

1.2 LSTM算法原理

傳統(tǒng)的RNN模型存在梯度消失、梯度爆炸等問題,會(huì)遺忘距離較遠(yuǎn)的信息,而LSTM模型作為RNN的一種特例[18],在細(xì)胞狀態(tài)也稱作記憶單元中增添遺忘門、輸入門和輸出門[19],可以有效控制信息的保留、遺忘和傳遞更新,使其具有長期記憶的能力,大幅提升了預(yù)測精度。圖1所示為一個(gè)LSTM單元結(jié)構(gòu)[20]。

圖1 LSTM單元結(jié)構(gòu)Fig.1 Structure of LSTM cell

當(dāng)上一時(shí)刻隱藏層輸出的ht-1和當(dāng)前時(shí)刻輸入的xt流入記憶單元后,遺忘門可以控制信息的取舍,通過Sigmoid函數(shù)得到當(dāng)前時(shí)刻遺忘門的輸出值ft:

式中,ft取值為0～1,表示對上一細(xì)胞狀態(tài)輸出信息的遺忘程度;Wf為這一步驟的權(quán)重系數(shù)矩陣;ht-1為前一時(shí)刻輸出值;xt為當(dāng)前時(shí)刻輸入值;bf為偏置項(xiàng);σ·()為Sigmoid函數(shù);t表示時(shí)刻。

輸入門接受到當(dāng)前時(shí)刻輸入的信息后,通過Sigmoid層和Tanh層共同篩選將新信息儲(chǔ)存到細(xì)胞狀態(tài)中。Sigmoid函數(shù)決定將要更新的輸出值it作為長期記憶。Tanh函數(shù)創(chuàng)建一個(gè)候選向量作為短期記憶。公式為

式中,Wi、Wc為當(dāng)前步驟Sigmoid層與Tanh層的權(quán)重系數(shù)矩陣;bi、bc為Sigmoid層與Tanh層的偏置項(xiàng)。

為了忘記部分之前的信息和記住部分新信息,當(dāng)前細(xì)胞狀態(tài)Ct可用上一狀態(tài)Ct-1與ft、it、共同表示:

輸出門確定細(xì)胞輸出的有效信息Ot,而后當(dāng)前細(xì)胞狀態(tài)Ct通過Tanh函數(shù)進(jìn)行處理生成新的隱含層輸出結(jié)果:

式中,Wo為當(dāng)前步驟輸出門權(quán)重系數(shù)矩陣;bo為偏置項(xiàng)。

2 形變監(jiān)測結(jié)果與分析

2.1 研究區(qū)概況與試驗(yàn)數(shù)據(jù)

西果園尾礦庫位于鞍鋼集團(tuán)鞍山礦業(yè)公司東鞍山燒結(jié)廠東南11 km處的西果園村[21]。庫區(qū)中心地理坐標(biāo)為東經(jīng)122°59′52″、北緯40°58′42″,東側(cè)距離鞍隆公路1.9 km,西側(cè)緊鄰排土場,庫南北長2.3 km,東西長1.6 km。該尾礦庫由3座庫組成,分別為西溝、南峪溝和井峪溝尾礦庫,西溝和南峪溝尾礦庫始建于20世紀(jì)70年代,于1979年投入使用,原設(shè)計(jì)標(biāo)高為260.0 m。1990年在南峪溝東側(cè)新建了井峪溝尾礦庫,該尾礦庫原設(shè)計(jì)堆積標(biāo)高260.0 m。尾礦庫采用上游法筑壩,堆積壩平均坡比為1∶5,庫中尾礦砂主要成分為石英、磁鐵礦、赤鐵礦、鈦鐵礦、石榴子石、白云母、方解石,還有少量穩(wěn)定和非穩(wěn)定礦物。在標(biāo)高225.0 m之前3座庫為獨(dú)立的尾礦庫,225.0 m之后三者堆積壩連成一體。后續(xù)為滿足生產(chǎn)需要加高擴(kuò)容后,設(shè)計(jì)最終標(biāo)高為275.0 m,總庫容為7 349.5×104m3,總壩體最高處為158 m,按照相關(guān)標(biāo)準(zhǔn),該尾礦庫為二等庫,研究區(qū)域如圖2所示。

圖2 研究區(qū)及尾礦庫范圍Fig.2 Scope of the study area and tailings pond

本次試驗(yàn)采用60景Sentinel-1A雷達(dá)衛(wèi)星升軌數(shù)據(jù),入射波為C波段,微波波長為5.6 cm,入射角為39°,重訪周期為12 d,時(shí)間跨度為2018年7月18日—2020年6月25日,極化方式為VV同極化,同時(shí)結(jié)合SRTM精度為30 m的DEM以及對應(yīng)時(shí)期的精密軌道文件進(jìn)行地形相位和軌道誤差的去除。

2.2 時(shí)序InSAR監(jiān)測結(jié)果

利用SARscape軟件,采用SBAS-InSAR技術(shù)對所獲取的影像進(jìn)行處理,選取日期為2019年4月8日的影像作為超級(jí)主影像,生成了318個(gè)干涉對;利用精密軌道數(shù)據(jù)進(jìn)行去平地效應(yīng),利用DEM數(shù)據(jù)消除地形相位,采用最小費(fèi)用流法(Minimum Cost Flow)進(jìn)行相位解纏;在進(jìn)行軌道精煉與重去平后,進(jìn)行形變速率反演與地理編碼,以獲取雷達(dá)視線方向的年平均形變速率值。

2018年7月—2020年6月研究區(qū)年平均形變速率分布如圖3所示。由圖3可知:尾礦庫堆積壩形變量較小,初期壩及壩底較為平穩(wěn),主要沉降區(qū)為尾礦庫西側(cè)排土場與尾礦庫頂部干灘區(qū)域,最大沉降速率為95 mm/a,年平均形變速率為-95～30 mm/a。圖中,A、B點(diǎn)為研究區(qū)內(nèi)2個(gè)GPS觀測站,G1、G2、G3點(diǎn)為后期建立預(yù)測模型時(shí)選取的點(diǎn)位,這3個(gè)點(diǎn)分別分布在3座尾礦庫頂部位置,且十分靠近最上層堆積壩,一旦發(fā)生劇烈形變,其監(jiān)測數(shù)據(jù)最具參考意義。

圖3 年平均形變速率Fig.3 Annual average deformation rate

由于本研究使用的影像較多,雖然已經(jīng)獲取了研究區(qū)整體形變速率,但為了更加直觀地展示出整個(gè)監(jiān)測時(shí)間段內(nèi)研究區(qū)的形變演化過程,繪制了等間隔的12幅時(shí)間序列形變圖,間隔時(shí)間為60 d,所有圖像均顯示以2018年7月18日為參考時(shí)間的時(shí)序累計(jì)形變值,具體如圖4所示。由圖4可知:隨著時(shí)間推移,尾礦庫區(qū)域的壩體頂部沉降量不斷增大,沉降范圍開始從中心向外擴(kuò)散,最終覆蓋整座尾礦庫,3座尾礦庫沉降趨勢及變化規(guī)律大致相同;在2 a內(nèi),研究區(qū)最大的地面累計(jì)沉降量達(dá)到-181 mm,最大累計(jì)抬升量為60 mm。

圖4 時(shí)序累計(jì)形變值Fig.4 Time series cumulative deformation values

為了檢驗(yàn)SBAS-InSAR技術(shù)的監(jiān)測精度,本研究在A、B兩觀測站點(diǎn)中獲取了與SAR影像時(shí)間相互對應(yīng)的10期GPS觀測數(shù)據(jù),將其與SBAS-InSAR監(jiān)測結(jié)果進(jìn)行對比驗(yàn)證。但由于SBAS-InSAR技術(shù)得到的形變?yōu)槔走_(dá)視線方向的形變,屬于真實(shí)形變在入射角上的投影,可以根據(jù)入射角的余弦值對視線方向形變進(jìn)行解算,得到垂直方向的形變值。GPS觀測值和SBAS-InSAR垂直方向沉降值結(jié)果對比如圖5所示。由圖5可知:A點(diǎn)最大絕對誤差、均方根誤差分別為2.8 mm、1.6 mm;B點(diǎn)最大絕對誤差、均方根誤差分別為2.3 mm、1.5 mm。兩種監(jiān)測方法所呈現(xiàn)的形變趨勢總體上較為一致,擬合度較高,說明SBAS-InSAR技術(shù)在該區(qū)域監(jiān)測精度較理想,適用于尾礦庫沉降監(jiān)測。

圖5 SBAS-InSAR與GPS形變監(jiān)測結(jié)果對比Fig.5 Comparison of the deformation monitoring results of SBAS-InSAR and GPS

2.3 時(shí)序InSAR形變分析

為了清晰地觀察西果園尾礦庫沉降趨勢,繪制出G1、G2、G33個(gè)點(diǎn)包含60個(gè)沉降記錄的沉降曲線,可以看出庫內(nèi)處于下沉狀態(tài),且形變量存在周期性波動(dòng),如圖6所示。

圖6 西果園尾礦庫沉降曲線Fig.6 Subsidence curves of Xiguoyuan Tailings Pond

在圖6所示的沉降曲線中可以發(fā)現(xiàn)其劇烈波動(dòng)時(shí)間段恰好與降雨時(shí)間吻合。本研究將沉降曲線進(jìn)行差分處理,得到形變的周期項(xiàng),其波動(dòng)原因主要受外部因素影響[22],所以結(jié)合研究時(shí)段內(nèi)的月降水量進(jìn)行分析,結(jié)果如圖7所示。

圖7 尾礦庫形變量波動(dòng)與降水的關(guān)系Fig.7 Relationship between deformation fluctuation and precipitation of tailings pond

由圖6及圖7可知:形變量的波動(dòng)與降水量呈正相關(guān),在2018年及2019年8月前后降水量最大時(shí),周期項(xiàng)波動(dòng)最為劇烈,波動(dòng)最大值達(dá)到12 mm,在2020年6月雨季開始時(shí),曲線開始出現(xiàn)抬升現(xiàn)象,而在10月至6月非雨季時(shí)間段內(nèi),其周期項(xiàng)在零值上下小幅波動(dòng)。這主要是因?yàn)榕磐翀黾皫靺^(qū)內(nèi)主要物質(zhì)為廢土及尾礦砂,大多含有黏土顆粒[23],在雨季期間降水量大時(shí),地表吸水膨脹,并且西果園尾礦庫為山谷型尾礦庫,大量的降水會(huì)降低壩體的黏結(jié)系數(shù),使得壩體穩(wěn)定性降低,短期內(nèi)出現(xiàn)抬升現(xiàn)象,隨后由于質(zhì)量增加,沉降速率較降水前也會(huì)變大,所以沉降曲線會(huì)出現(xiàn)先抬升、再加速下降的現(xiàn)象,在降水減少時(shí),庫內(nèi)水分由于蒸騰作用逐漸流失,沉降速率及其波動(dòng)曲線也趨于平穩(wěn)。由此可見,降水對庫區(qū)沉降影響較大,在后續(xù)的預(yù)測中應(yīng)考慮這一波動(dòng)現(xiàn)象。

3 GS-LSTM模型構(gòu)建及預(yù)測結(jié)果分析

在構(gòu)建LSTM網(wǎng)絡(luò)模型時(shí),眾多超參數(shù)的選取對模型的預(yù)測精度會(huì)產(chǎn)生較大影響,為了提高預(yù)測精度,本研究采用網(wǎng)格搜索算法搜索全局最優(yōu)解[24]。網(wǎng)格搜索算法是窮舉搜索的一種,在所有候選的超參數(shù)選擇中,通過循環(huán)遍歷,將超參數(shù)互相組合,嘗試每一種可能性,以得到全局最優(yōu)的超參數(shù)組合。

本研究將前50期時(shí)序沉降數(shù)據(jù)作為訓(xùn)練樣本,后10期數(shù)據(jù)作為測試樣本進(jìn)行沉降預(yù)測。由于數(shù)據(jù)集不夠龐大,數(shù)據(jù)復(fù)雜度相對較低,為避免模型過擬合,搭建了單層LSTM神經(jīng)網(wǎng)絡(luò),在使用GS算法調(diào)參后,最終設(shè)置批次大小batch_size=2,窗口大小window_size=50,神經(jīng)元個(gè)數(shù)num_units=128,學(xué)習(xí)率為0.001,迭代次數(shù)為4 500次,將均方根誤差作為損失函數(shù),訓(xùn)練過程中采用Adam優(yōu)化算法更新網(wǎng)絡(luò)權(quán)值。為驗(yàn)證該模型預(yù)測結(jié)果的準(zhǔn)確性,將其與傳統(tǒng)的SVR模型和灰色GM(1,1)模型結(jié)果互相比對,逐一分析,結(jié)果見表1至表4。

表1 G1點(diǎn)各預(yù)測模型結(jié)果對比Table 1 Comparison of prediction results of each model at G1 point mm

表2 G2點(diǎn)各預(yù)測模型結(jié)果對比Table 2 Comparison of prediction results of each model at G2 point mm

表3 G3點(diǎn)各預(yù)測模型結(jié)果對比Table 3 Comparison of prediction results of each model at G3 point mm

表4 各模型精度對比Table 4 Comparison of accuracy of each modelmm

由表1至表4可知:3種預(yù)測方法與沉降值均較為一致,為更直觀地分析及評價(jià)各模型的精度,分別計(jì)算各預(yù)測模型的平均絕對誤差(MAE)和均方根誤差(RMSE),其中GS-LSTM模型最大MAE及RMSE為2.51 mm、2.90 mm,SVR模型分別為3.58 mm、4.48 mm,GM(1,1)模型分別為4.01 mm、4.95 mm。在兩種評價(jià)指標(biāo)上GS-LSTM模型預(yù)測結(jié)果誤差均為最小,優(yōu)于傳統(tǒng)的SVR模型及GM(1,1)模型。

3種模型預(yù)測結(jié)果與實(shí)測值對比如圖8所示。由圖8可知:SVR模型和GM(1,1)模型可以較好地預(yù)測沉降趨勢,但是尾礦庫的沉降特點(diǎn)及波動(dòng)性無法得到有效預(yù)測,在G3點(diǎn)表現(xiàn)得尤為突出,沉降波動(dòng)劇烈時(shí)其誤差也隨之增大,而GS-LSTM模型在預(yù)測時(shí)即使波動(dòng)劇烈,仍能保持較高的精度,預(yù)測結(jié)果十分理想,為后續(xù)的預(yù)測模型持續(xù)優(yōu)化與應(yīng)用提供了思路。

圖8 不同監(jiān)測點(diǎn)處3種模型預(yù)測精度對比Fig.8 Comparison of the prediction precision of three models at different monitoring points

4 結(jié) 論

(1)為有效預(yù)防和減少尾礦庫事故的發(fā)生,以鞍山西果園尾礦庫為例,利用SBAS-InSAR技術(shù)對庫區(qū)進(jìn)行沉降監(jiān)測,提取連續(xù)沉降信息,在此基礎(chǔ)上使用GS-LSTM模型對庫區(qū)沉降值進(jìn)行預(yù)測,實(shí)現(xiàn)了尾礦庫沉降監(jiān)測與預(yù)測的一體化。研究表明:本研究所提方法的監(jiān)測值與實(shí)測值的誤差較小,預(yù)測值與傳統(tǒng)模型相比不僅精度更高,并且能記憶尾礦庫的沉降波動(dòng)特點(diǎn),監(jiān)測精度和預(yù)測精度均滿足工程需要,可為尾礦庫沉降規(guī)律研究及沉降預(yù)測提供依據(jù)。

(2)后續(xù)工作中考慮使用層析SAR技術(shù)結(jié)合更高空間和時(shí)間分辨率雷達(dá)數(shù)據(jù)獲取庫區(qū)三維立體形變信息,進(jìn)一步提升尾礦庫監(jiān)測與預(yù)測模型的精度。