黃會寶,江德軍,劉恒,羅忠平

(1.四川大學(xué)水利水電學(xué)院,成都 610065; 2.國能大渡河流域水電開發(fā)有限公司,成都 610095)

作為江河蓄水防旱、抗洪減澇的重要樞紐工程,水電站通過調(diào)控水資源時空分布來實(shí)現(xiàn)水資源合理調(diào)控與配置優(yōu)化[1-2],給中國社會經(jīng)濟(jì)帶來了極大便利。但是水電站的建設(shè)與運(yùn)營會引起周邊環(huán)境與地質(zhì)構(gòu)造變化,影響大壩、庫區(qū)及壩下整個水利樞紐系統(tǒng)的正常運(yùn)維[3],特別是持續(xù)的地面變形可能造成壩體損壞、水系基礎(chǔ)設(shè)施破壞等危害,嚴(yán)重威脅人類生命安全和公私財產(chǎn),造成經(jīng)濟(jì)損失[4]。因此,對大壩及其周邊地區(qū)地表形變的時空變化進(jìn)行監(jiān)測、分析和預(yù)測是十分必要的。

傳統(tǒng)的大地測量技術(shù)如全球?qū)Ш叫l(wèi)星系統(tǒng)(global positioning system,GPS)[5]、水準(zhǔn)測量[6]、三維激光掃描[3]由于受局部單點(diǎn)測量、空間分辨率低、成本高的限制,在水電站大壩區(qū)域變形監(jiān)測存在分辨率低、成本高等劣勢。近10年來,干涉合成孔徑雷達(dá)(interferometric synthetic aperture radar,InSAR)技術(shù)已被廣泛應(yīng)用于城市地面沉降[7]、地鐵及高速公路沿線[8-9]和滑坡災(zāi)害[10]等。多時間干涉合成孔徑雷達(dá)(multi-temporal interferometric synthetic aperture radar,MT-InSAR)技術(shù)[11-12]可實(shí)現(xiàn)大規(guī)模、高精度、高密度和低成本的地表形變監(jiān)測,并在水電站及周邊地區(qū)形變監(jiān)測方面得到了廣泛應(yīng)用[13-17]。典型的MT-InSAR技術(shù)之一是永久散射體技術(shù)(persistent scatterer interferometric synthetic aperture radar,PS-InSAR)[18],該技術(shù)選擇時序穩(wěn)定的點(diǎn)目標(biāo)并基于SAR(synthetic aperture radar)圖像中的相位信息進(jìn)行時間序列建模和分析,從而大大降低了大氣延遲引起的誤差,使變形監(jiān)測精度達(dá)到毫米級,這些穩(wěn)定的點(diǎn)一般為建筑物、橋梁、大壩和暴露的巖石等。另一種有效的監(jiān)測方法是小基線子集干涉合成孔徑雷達(dá)(small BAseline subset InSAR,SBAS-InSAR)技術(shù)[19],該技術(shù)是選擇時空基線短的干涉圖并形成干涉網(wǎng)絡(luò)并進(jìn)行形變解算,可有效減少大氣延遲和相關(guān)噪聲誤差。這兩種方法在大壩形變監(jiān)測中取得了很大的成功。然而,時間序列InSAR測量的形變結(jié)果很少用于地面沉降的預(yù)測。

現(xiàn)有的大多數(shù)形變預(yù)測方法主要有基于數(shù)學(xué)統(tǒng)計(jì)模型[20]、經(jīng)驗(yàn)?zāi)Ｐ蚚21]和傳統(tǒng)的人工神經(jīng)網(wǎng)絡(luò)模型[22]。其中數(shù)學(xué)模型包括回歸分析模型和時間序列方法。但是回歸分析方法的準(zhǔn)確性受樣本數(shù)量的影響較大[23],時間序列的方法容易受到歷史數(shù)據(jù)的質(zhì)量可能會導(dǎo)致預(yù)測性能下降[24]。而經(jīng)驗(yàn)?zāi)Ｐ头椒ㄝ斎雲(yún)?shù)較多且計(jì)算復(fù)雜[25]。人工神經(jīng)網(wǎng)絡(luò)模型雖然具有較高的自學(xué)習(xí)能力,但是其輸出結(jié)果受初始權(quán)值和閾值的影響較大。此外,該方法無法準(zhǔn)確捕捉時間序列形變的波動,往往無法獲得令人滿意的預(yù)測結(jié)果。近10年來,深度學(xué)習(xí)在遙感圖像處理與分析領(lǐng)域的快速發(fā)展顯示出廣闊的應(yīng)用潛力。Anantrasirichai等[26]和Valade等[27]卷積神經(jīng)網(wǎng)絡(luò)(convolutional neural network,CNN)被用于探測火山的緩慢變形和預(yù)測短期InSAR變形圖。Nukala等[28]提出了一種基于Sentinel-1圖像遞推神經(jīng)網(wǎng)絡(luò)(recurrent neural network,RNN)預(yù)測時間序列形變圖的新方法,并取得了良好的預(yù)測性能。長短時記憶網(wǎng)絡(luò)(long short-term memory,LSTM)網(wǎng)絡(luò)解決了循環(huán)神經(jīng)網(wǎng)絡(luò)(RNN)變體在學(xué)習(xí)數(shù)據(jù)的長期依賴性時可能遭受的梯度爆炸和消失的限制,Chen等[29]利用LSTM網(wǎng)絡(luò)建立了時間序列InSAR變形預(yù)測模型,并指出該網(wǎng)絡(luò)具有較好的預(yù)測性能。陳毅等[30]基于時序 InSAR數(shù)據(jù)采用堆疊式LSTM模型來預(yù)測香港國際機(jī)場的地面沉降的短期預(yù)測。然而,這些深度學(xué)習(xí)模型在水電站大壩地區(qū)的地表形變預(yù)測沒有得到應(yīng)用,且在提取大壩地區(qū)InSAR時間序列的形變特征能力較差。此外,InSAR時間序列形變數(shù)據(jù)集的獲取也非常困難。

現(xiàn)采用LSTM神經(jīng)網(wǎng)絡(luò)對瀑布溝水電站的地表時間序列InSAR形變進(jìn)行預(yù)測。首先利用MT-InSAR技術(shù)對2018年1月至 2020年4月大渡河瀑布溝水電站的上行軌道哨兵一號(Sentinel-1)圖像進(jìn)行時間序列形變監(jiān)測;然后基于時序InSAR形變數(shù)據(jù)建立LSTM神經(jīng)網(wǎng)絡(luò)并對利用LSTM神經(jīng)網(wǎng)絡(luò)對研究區(qū)的時序形變進(jìn)行預(yù)測,最后預(yù)測出研究區(qū)在2020年5—9月累積沉降量。為瀑布溝水電站的長時序形變監(jiān)測提供定量評估的數(shù)據(jù),也為該地區(qū)的地面沉降災(zāi)害預(yù)警和輔助決策提供理論依據(jù)。

1 研究區(qū)和數(shù)據(jù)集

1.1 研究區(qū)介紹

瀑布溝水電站位于四川省雅安市漢源縣和涼山彝族自治州甘洛縣交界處(圖1),是大渡河中游控制性水電工程[31]。電站為發(fā)電、防洪、攔沙等多項(xiàng)利用結(jié)合的大型工程,同時它是國家“十五”重點(diǎn)建設(shè)項(xiàng)目,也是西部大開發(fā)的標(biāo)志性工程。瀑布溝水電站總裝機(jī)容量360萬kW,年均發(fā)電量147.9億kW·h,總庫容53.9億m3。它的大壩是深厚覆蓋層上建成的礫石土心墻堆石壩,最大壩高186 m。

圖2(a)所示為大渡河位于青藏高原與四川盆地的過渡地帶,地勢陡峭,高程范圍為564～5 947 m,如圖2(b)所示。研究區(qū)由于地處揚(yáng)子地臺西部古老的康滇地軸北段東側(cè),且水庫周圍分布了多條大型活動斷裂帶,如水庫以北的龍門山斷裂帶和,水庫西北的鮮水河斷裂帶以及水庫以南的安寧河-則木河斷裂和大涼山斷裂帶,這些斷裂帶自有記錄以來均發(fā)生過中等以上地震[31-33],因此研究區(qū)地質(zhì)構(gòu)造活動頻繁。

圖2 研究區(qū)Sentinel-1數(shù)據(jù)覆蓋范圍與SRTM DEM圖

1.2 數(shù)據(jù)集

為加強(qiáng)對大渡河流域瀑布溝水電站的形變監(jiān)測與預(yù)測,獲取2018年1月2日至2020年4月21日VV極化的31景Sentinel-1 SAR數(shù)據(jù)以反演研究區(qū)的歷史形變信息。研究區(qū)Sentinel-1數(shù)據(jù)覆蓋范圍圖如圖2(a)所示。Sentinel-1衛(wèi)星的重復(fù)周期為12 d(2014年4月3日發(fā)射),波長約為5.6 cm[34],這些圖像是通過漸進(jìn)掃描地形觀測(terrain observation by progressive scans,TOPS)的干涉寬幅(interferometric wide swath,IW)模式在4年內(nèi)沿上行軌道T26采集的。Sentinel-1圖像的幅寬為250 km,視線(line of sight,LOS)的入射角在30°～47°變化。同時獲取地面分辨率為30 m的1弧秒分辨率的航天飛機(jī)雷達(dá)地形測量數(shù)字高程模型(shuttle radar topography mission digital elevation model,SRTM DEM) 用來模擬InSAR技術(shù)形變解算中的地形相位,并利用它來進(jìn)行形變結(jié)果的地理編碼。

2 研究方法

采用基于LSTM的模型[35]對MT-InSAR技術(shù)反演的InSAR時序形變信息來預(yù)測瀑布溝水電站未來時期SAR圖像采集所對應(yīng)的時序位移。圖3顯示了本文提出方法的概述和詳細(xì)步驟。

圖3 基于多時相InSAR與LSTM方法瀑布溝水電站形變監(jiān)測與預(yù)測流程圖

2.1 Sentinel-1 TOPS模式預(yù)處理

將Sentinel-1 TOPS模式數(shù)據(jù)導(dǎo)入生成SLC格式作為輸入數(shù)據(jù),并由歐空局提供精密軌道數(shù)據(jù)用于精煉軌道參數(shù)信息;接著進(jìn)行Burst分割,將條帶的單視復(fù)數(shù)據(jù)(single-look complex,SLC)分割成9個Burst進(jìn)行后續(xù)處理。

采用30 m分辨率的SRTM數(shù)據(jù)進(jìn)行Sentinel-1數(shù)據(jù)的DEM配準(zhǔn)。由于TOPS模式SAR圖像存在著多普勒質(zhì)心的變化,圖像配準(zhǔn)精度需要在方位向上優(yōu)于千分之一像素,因此采用增強(qiáng)譜分集(ESD,enhanced spectral diversity)方法進(jìn)一步細(xì)化方位向的偏移量。在從圖像重采樣到主圖像的坐標(biāo)系框架之前,對從圖像進(jìn)行去斜和解調(diào)。同時,在每個Burst級別上執(zhí)行去斜、插值和反去斜等流程。然后采用ESD方法對整體方位向偏移量進(jìn)行改正,方位向的配準(zhǔn)精度可小于千分之一個像素。最后,將每個條帶的所有Burst拼接成有一個子條帶,并對該條帶進(jìn)行感興趣區(qū)(region of interest,ROI)的提取,獲取配準(zhǔn)后時間序列SAR數(shù)據(jù)集。

2.2 多時相InSAR處理

首先根據(jù)配準(zhǔn)后的SAR數(shù)據(jù)集生成相干系數(shù)圖和差分干涉圖,并采用平均相干系數(shù)法和振幅離差法進(jìn)行高相干點(diǎn)提取,其次根據(jù)空間和時間閾值選擇符合條件的小基線干涉圖,將這些高相干點(diǎn)進(jìn)行Delaunay三角網(wǎng)的構(gòu)建,建立相鄰相干點(diǎn)之間的差分相位模型,然后利用最小二乘法求解線性形變率差和DEM改正值差。接著選取參考相干點(diǎn)計(jì)算各相干點(diǎn)的DEM改正值和線性形變率。最后根據(jù)殘余相位的時空特性,采用時空濾波的方法提取大氣相位。該步驟可分為殘余相位最小費(fèi)用流(minimum cost flow,MCF)相位解纏、殘余相位分離和大氣相位估計(jì)等子步驟。線性速率和高程改正能有效地反映形變場的分布,還需要對原始相位進(jìn)行修正,主要考慮到線性速率和高程誤差初始值對迭代的影響。在獲得這兩項(xiàng)以后,對殘余相位進(jìn)行濾波,濾波完成后利用殘余相位與濾波后相位之差,可以認(rèn)為其為大氣和噪聲的影響,剩余的則為非線性形變量。將此加入線性形變量中,就可以得到真實(shí)的形變量。在選取的相干目標(biāo)上,建立觀測方程。采用奇異值分解(singular value decomposition,SVD)的方法進(jìn)行求解,即可得到各時間區(qū)間的沉降速率。各時段沉降速率在時間域上進(jìn)行積分即可得到各個時間點(diǎn)的累積時序形變量。

2.3 多時相InSAR時序形變預(yù)測

對上述獲取的瀑布溝水電站InSAR時序形變數(shù)據(jù)來構(gòu)建LSTM網(wǎng)絡(luò),然后對LSTM構(gòu)建訓(xùn)練集、測試集、驗(yàn)證集來進(jìn)行模型訓(xùn)練和測試,最終根據(jù)LSTM模型輸出結(jié)果獲取預(yù)測的InSAR時序形變。

LSTM神經(jīng)網(wǎng)絡(luò)是傳統(tǒng)RNN的變體,并被廣泛應(yīng)用于時間序列建模任務(wù)。LSTM神經(jīng)網(wǎng)絡(luò)的兩個主要模塊幫助它從數(shù)據(jù)中學(xué)習(xí)時間序列特征。門模塊,包括輸入門、輸出門和遺忘門,可有效地訓(xùn)練全連接層控制單元的狀態(tài),以響應(yīng)來自數(shù)據(jù)的新輸入和建模過去的輸出[36]。相比于RNN網(wǎng)絡(luò),LSTM的門模塊有效地克服了傳統(tǒng)RNN的梯度消失或爆炸、無法學(xué)習(xí)長期依賴的問題。

InSAR監(jiān)測結(jié)果可表示為x={x1,x2,…,xn}。該預(yù)測模型可以有效地學(xué)習(xí)時間序列變形數(shù)據(jù)的變化特征,通過LSTM預(yù)測模型可以得到地表變形的預(yù)測結(jié)果y,即y={y1,y2,…,yn}。圖3顯示了構(gòu)建的LSTM模型,模型由2層LSTM層、2層Dense層和3層Dropout層組成。圖3顯示了LSTM神經(jīng)網(wǎng)絡(luò)的基本結(jié)構(gòu)。LSTM神經(jīng)網(wǎng)絡(luò)的第一階段決定細(xì)胞狀態(tài)下的信息是遺忘還是記憶,如果當(dāng)前數(shù)據(jù)重要,則遺忘門的sigmoid層允許遺忘之前的數(shù)據(jù)。如果當(dāng)前數(shù)據(jù)不需要,它可記住以前的數(shù)據(jù),這解決了RNN算法存在長期依賴問題的缺點(diǎn)。遺忘門的計(jì)算公式為

ft=σ[Wf(ht-1,xt)+bf]

(1)

式(1)中:ft為遺忘門;σ為一個sigmoid激活函數(shù);ht-1為t-1時刻的隱藏狀態(tài);xt為t時刻的輸入向量;Wf和bf分別為遺忘門的權(quán)值向量和偏差向量。如果ft的輸出值接近于0,則表示之前的數(shù)據(jù)已經(jīng)被遺忘,但是如果它接近1,并不意味著之前的數(shù)據(jù)已被記住。

it=σ[Wi(ht-1,xt)+bi]

(2)

(3)

式中:it為當(dāng)前時間步t的輸入門的輸出;σ為一個sigmoid激活函數(shù),它是輸入門;Wi和bi分別為輸入門的權(quán)值向量和偏差向量;Wc為更新的權(quán)重;bc為偏差。

(4)

LSTM神經(jīng)網(wǎng)絡(luò)的最后階段決定輸出什么。第一步是通過sigmoid層輸出Ot,該層用于確定將輸出單元格狀態(tài)的哪一部分。然后,細(xì)胞狀態(tài)經(jīng)過tanh函數(shù);該操作完成后單元格狀態(tài)值在-1～1伸縮。最后,將結(jié)果乘以sigmoid門的輸出來確定輸出,計(jì)算公式為

Ot=σ[Wo(ht-1,xt)+bo]

(5)

ht=Ottanh(Ct)

(6)

式中:Ot為sigmoid層的輸出;ht-1為舊的輸入值;ht為新的輸出值;Ot為輸出門;Wo和bo分別為輸出門的權(quán)值向量和偏移向量。

3 結(jié)果與分析

基于31景Sentinel-1升軌圖像利用MT-InSAR技術(shù)獲得2018年1月2日—2020年4月21日年InSAR LOS方向年平均形變速率結(jié)果如圖4所示。

紅色負(fù)值表示地表沉降

采用MT-InSAR方法提取的Sentinel-1數(shù)據(jù)集中約有4 319個高相干點(diǎn),平均密度為626個相干點(diǎn)/km2。如圖4所示,可看到瀑布溝水電站大壩處的形變模式,年平均形變速率范圍為-34～10 mm/a。如圖5所示,動態(tài)地顯示了瀑布溝水電站形變規(guī)律變化,從2018年1月2日—2020年4月21日瀑布溝水電站大壩附近累積沉降量呈現(xiàn)逐年增大現(xiàn)象。

圖5 瀑布溝水電站4個時期的累積形變量圖

根據(jù)研究區(qū)的形變速率圖,將瀑布溝水電站的沉降區(qū)域劃分為加速沉降區(qū)、線性沉降區(qū)和穩(wěn)定區(qū)。同時,從圖6中提取了3個區(qū)域3個點(diǎn)A、B、C的時間序列InSAR形變量進(jìn)一步分析。

圖6 瀑布溝水電站形變區(qū)域放大圖

首先,獲得3個區(qū)域3個高相干點(diǎn)的歷史形變結(jié)果,利用TensorFlow深度學(xué)習(xí)框架建立研究區(qū)測點(diǎn)地面形變的LSTM預(yù)測模型;利用LSTM預(yù)測模型對InSAR數(shù)據(jù)集進(jìn)行劃分,實(shí)現(xiàn)訓(xùn)練和測試過程。70%的時間序列形變數(shù)據(jù)用于訓(xùn)練,20%的時間序列形變數(shù)據(jù)用于測試,10%的時間序列形變數(shù)據(jù)用于驗(yàn)證;其次,采用最小值和最大值歸一化方法對這些形變點(diǎn)的InSAR時間序列變形數(shù)據(jù)進(jìn)行歸一化處理;然后,對LSTM模型進(jìn)行網(wǎng)絡(luò)初始化步驟,即設(shè)置初始化權(quán)值和偏移向量,并利用Dropout層抑制網(wǎng)絡(luò)訓(xùn)練過程中的過擬合。初始學(xué)習(xí)率為0.000 1,最大訓(xùn)練時間為1 000,優(yōu)化器采用Adam方法。采用網(wǎng)格搜索算法對預(yù)測模型中的超參數(shù)進(jìn)行選擇。

最終利用LSTM模型對2018年1月2日—2019年7月20日22景累積形變量進(jìn)行模型訓(xùn)練,利用2019年8月1日—2020年4月21日進(jìn)行模型測試,利用2020年5月3日—2020年9月12日進(jìn)行模型驗(yàn)證,即預(yù)測出后續(xù)12個時刻的InSAR形變量,形變監(jiān)測和預(yù)測結(jié)果圖如圖7所示。

圖7 點(diǎn)A、B、C的時序形變監(jiān)測與預(yù)測對比圖

圖7(a)顯示了截至2020年4月21日瀑布溝水電站位于沉降加速區(qū)的點(diǎn)A最大沉降量為64.62 mm。2020年5月3日—9月12日的預(yù)測結(jié)果顯示了點(diǎn)A有繼續(xù)沉降的趨勢,預(yù)測的最大沉降值為71.29 mm。

圖7(b)顯示了瀑布溝水電站位于線性沉降區(qū)域的點(diǎn)B沉降的最大值為11.70 mm,顯示了該點(diǎn)呈現(xiàn)出線性變化規(guī)律。2020年5月3日—9月12日的預(yù)測結(jié)果顯示點(diǎn)B繼續(xù)出現(xiàn)沉降趨勢,但是沉降幅度較小。圖7(c)顯示了瀑布溝水電站位于穩(wěn)定區(qū)域的點(diǎn)C形變變化穩(wěn)定,累積形變范圍為-8.29～6.67 mm。

為了評價所構(gòu)建的LSTM神經(jīng)網(wǎng)絡(luò)模型的預(yù)測精度,使用均方根誤差(root mean squared error,RMSE)和平均絕對誤差(mean absolute error,MAE)來評價預(yù)測模型的預(yù)測精度。如表1所示點(diǎn)A、B、C的LSTM模型在訓(xùn)練過程中最小的RMSE和MAE值為2.343 mm和2.010 mm,A、B、C的LSTM模型在測試過程中最小的RMSE和MAE值為2.094 mm和1.654 mm,這說明LSTM模型專為處理序列數(shù)據(jù)而設(shè)計(jì),在學(xué)習(xí)時間序列數(shù)據(jù)的特征方面具有獨(dú)特的優(yōu)勢,且對大壩地區(qū)的形變預(yù)測具有一定的應(yīng)用潛力。

4 結(jié)論

采用LSTM模型基于MT-InSAR技術(shù)對瀑布溝水電站的沉降區(qū)域進(jìn)行形變監(jiān)測與預(yù)測研究,得出以下結(jié)論。

(1)MT-InSAR技術(shù)可明顯探測到瀑布溝水電站的沉降趨勢。2018年1月2日—2020年4月21日年的形變結(jié)果顯示LOS方向年平均形變速率范圍為-34～10 mm/a。研究證明了MT-InSAR技術(shù)在地形復(fù)雜的大壩區(qū)域形變監(jiān)測中的應(yīng)用潛力。

(2)結(jié)合MT-InSAR 技術(shù)創(chuàng)新性地提出了LSTM神經(jīng)網(wǎng)絡(luò)時序形變預(yù)測模型,該模型在大壩形變監(jiān)測與預(yù)測領(lǐng)域取得了良好的效果。LSTM模型在訓(xùn)練和測試過程中點(diǎn)尺度的RMSE和MAE最小值分別為2.343 mm和2.010 mm,2.094 mm和1.654 mm。LSTM神經(jīng)網(wǎng)絡(luò)能夠較好地捕獲時間序列InSAR數(shù)據(jù)的形變規(guī)律和波動,該方法適合于大壩地區(qū)的InSAR時間序列沉降數(shù)據(jù)的預(yù)測。

(3)基于2018—2020年InSAR時間序列沉降數(shù)據(jù),采用LSTM神經(jīng)網(wǎng)絡(luò)對2020年5—9月InSAR時間序列形變值進(jìn)行預(yù)測。截至2020年9月,本文提出的方法預(yù)測的最大沉降值將達(dá)到71.29 mm,需要引起相關(guān)部門的重視。LSTM神經(jīng)網(wǎng)絡(luò)是一種有前景的時間序列InSAR地面沉降預(yù)測方法,為大壩地區(qū)的周期性監(jiān)測和預(yù)警提供了重要的數(shù)據(jù)支撐,可服務(wù)于地質(zhì)災(zāi)害的早期預(yù)警和提供防災(zāi)決策信息。隨著Sentinel-1 SAR數(shù)據(jù)指數(shù)級的增長,未來將在LSTM模型中加入更多的數(shù)據(jù),準(zhǔn)確量化大壩地面沉降與外部數(shù)據(jù)的關(guān)系,判斷模型是否適合長期預(yù)測,以提高預(yù)測模型的魯棒性。