融合SBAS-InSAR與GS-LSTM的尾礦庫沉降監(jiān)測與預測

李如仁 孫加瑤

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

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

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

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

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

1.1 SBAS-InSAR基本原理

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

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

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

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

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

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

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

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

1.2 LSTM算法原理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 時序累計形變值Fig.4 Time series cumulative deformation values

為了檢驗SBAS-InSAR技術(shù)的監(jiān)測精度,本研究在A、B兩觀測站點中獲取了與SAR影像時間相互對應的10期GPS觀測數(shù)據(jù),將其與SBAS-InSAR監(jiān)測結(jié)果進行對比驗證。但由于SBAS-InSAR技術(shù)得到的形變?yōu)槔走_視線方向的形變,屬于真實形變在入射角上的投影,可以根據(jù)入射角的余弦值對視線方向形變進行解算,得到垂直方向的形變值。GPS觀測值和SBAS-InSAR垂直方向沉降值結(jié)果對比如圖5所示。由圖5可知:A點最大絕對誤差、均方根誤差分別為2.8 mm、1.6 mm;B點最大絕對誤差、均方根誤差分別為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 時序InSAR形變分析

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

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

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

圖7 尾礦庫形變量波動與降水的關系Fig.7 Relationship between deformation fluctuation and precipitation of tailings pond

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

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

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

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

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

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

表3 G3點各預測模型結(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種預測方法與沉降值均較為一致,為更直觀地分析及評價各模型的精度,分別計算各預測模型的平均絕對誤差(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。在兩種評價指標上GS-LSTM模型預測結(jié)果誤差均為最小,優(yōu)于傳統(tǒng)的SVR模型及GM(1,1)模型。

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

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

4 結(jié) 論

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

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