基于InSAR相位梯度疊加的毛爾蓋庫岸滑坡隱患快速識別

唐光民, 戴可人,,3, 卓冠晨,, 沈 月,4, 陳 晨, 許 強(qiáng)

(1. 成都理工大學(xué) 地球科學(xué)學(xué)院, 四川 成都 610059;2. 地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室, 四川 成都 610059;3. 長安大學(xué) 地質(zhì)工程與測繪學(xué)院, 陜西 西安 710054;4. 中國地質(zhì)大學(xué)地球物理與空間信息學(xué)院, 湖北 武漢 430074)

0 引言

庫岸滑坡失穩(wěn)可能會(huì)沖毀沿岸建筑物,或阻塞河道,對庫岸地質(zhì)環(huán)境及沿岸居民生命財(cái)產(chǎn)安全造成巨大威脅[1]。國內(nèi)外由庫岸滑坡引發(fā)的重大災(zāi)害時(shí)有發(fā)生,例如:意大利瓦依昂滑坡不僅造成整個(gè)水庫堵塞,山體滑坡引發(fā)的洪水更是造成2 000余人死亡,多個(gè)村莊和城鎮(zhèn)夷為平地[2];法國馬爾帕塞拱壩事故造成死亡和失蹤的人數(shù)超過了500人[3];國內(nèi)湖南柘溪水庫塘巖光庫岸滑坡造成64人淹溺死亡,24人受傷[4];湖北的黃龍灘水庫滑坡造成2 316戶居民搬遷,給庫區(qū)經(jīng)濟(jì)帶來巨大損失[5]。水庫水位短期的劇烈變動(dòng),勢必改變庫區(qū)原有地質(zhì)環(huán)境,從而引發(fā)滑坡等地質(zhì)災(zāi)害的發(fā)生。庫岸滑坡災(zāi)害分布廣泛,多發(fā)頻發(fā),傳統(tǒng)測量手段和群測群防方式耗時(shí)耗力,難以進(jìn)行有效排查[6]。因此,如何快速識別庫岸滑坡從而進(jìn)行有效防災(zāi)減災(zāi)是亟待解決的重要問題。

近些年,合成孔徑雷達(dá)干涉測量(Interferometric Synthetic Aperture Radar,InSAR)因其具有覆蓋范圍廣、監(jiān)測精度高、全天時(shí)、全天侯、穿云透霧等特點(diǎn)已成為滑坡探測領(lǐng)域的重要技術(shù)手段[6-10]。常用的時(shí)間序列InSAR方法,例如小基線集干涉測量 (Small Baselines Subset InSAR,SBAS-InSAR)、干涉影像堆疊(Stacking)技術(shù)等,被廣泛應(yīng)用于滑坡隱患識別中,大量成功的案例證明了時(shí)序InSAR技術(shù)具備識別滑坡隱患的能力[11-15]。但是由于時(shí)序InSAR技術(shù)涉及相位解纏[16]、時(shí)序建模解算[17]等復(fù)雜步驟,操作相對繁瑣,對操作者經(jīng)驗(yàn)要求較高,需要消耗大量時(shí)間,尤其面對海量數(shù)據(jù),識別效率不足。相比之下,基于InSAR相位梯度疊加方法受到參數(shù)設(shè)置的影響較小,同時(shí)對微小形變具有更高的敏感性,能夠揭示更小尺度的形變變化。此外,由于該方法直接從差分干涉對中計(jì)算梯度,避免了相位解纏和時(shí)序建模解算等步驟,能夠有效縮短數(shù)據(jù)處理時(shí)間,提高數(shù)據(jù)處理效率。例如,Fu等[18]利用梯度疊加算法獲取西南山區(qū)InSAR干涉圖距離向和方位向的相位梯度,通過訓(xùn)練YOLOv3神經(jīng)網(wǎng)絡(luò)進(jìn)行快速識別,成功探測西南地區(qū)3 979處潛在滑坡。Shen等[19]利用改進(jìn)的Sobel算子求取D-InSAR處理后的9張差分干涉圖相位梯度,去除各項(xiàng)形變誤差后,采用閾值分割的方法對滑坡隱患進(jìn)行自動(dòng)識別,整體精度與人工解譯結(jié)果相比達(dá)到81%,然而其僅針對D-InSAR處理后兩景影像進(jìn)行梯度處理,無法探測長時(shí)間序列的滑坡隱患。目前通過InSAR相位梯度快速識別長時(shí)間序列滑坡隱患值得進(jìn)一步研究。

本文在不涉及相位解纏和時(shí)序建模解算的情況下,對InSAR差分干涉對進(jìn)行梯度運(yùn)算,識別出相位信息不連續(xù)的區(qū)域。通過去除植被等低相干性區(qū)域、幾何畸變、水體、積雪等形變誤差后,采用相位梯度疊加方式獲取相位不連續(xù)區(qū)域。采用閾值分割生成二值圖,并進(jìn)行邊緣提取,以獲取發(fā)生形變的區(qū)域,將獲取的形變區(qū)域與Stacking-InSAR和SBAS-InSAR處理的結(jié)果進(jìn)行對比驗(yàn)證分析,探究本文方法在滑坡隱患識別方面的特點(diǎn)。最后,統(tǒng)計(jì)三種方法從數(shù)據(jù)配準(zhǔn)到識別出滑坡隱患所耗時(shí)間,總結(jié)本文方法的優(yōu)勢。

1 研究區(qū)域與數(shù)據(jù)

1.1 研究區(qū)域

毛爾蓋水庫隸屬四川省阿壩藏族自治州黑水縣,地處青藏高原東南緣,屬于該區(qū)橫斷山脈中段地區(qū)[圖1(a)]。南東面與茂縣相鄰,南西面則與理縣相接,北東毗鄰松潘[圖1(b)],地理坐標(biāo)為:102°35′～103°30′E、 31°35′～32°38′N[20],屬于“北亞熱帶”氣候、川西高原氣候區(qū),每年有旱季和雨季,降雨分布不均,主要集中在夏季[21]。

圖1 研究區(qū)概況Fig.1 Overview of the study area

庫區(qū)內(nèi)修建的毛爾蓋水電站位于四川省阿壩州黑水縣麻窩鄉(xiāng)境內(nèi)[圖1(c)],系岷江上游最大的一級支流黑水河流域梯級開發(fā)的中游“龍頭”水庫電站[22],山勢陡峭、溝谷縱橫,谷底寬廣,呈不對稱“U”型,壩址河段開闊,山勢雄偉,地形完整,坡度為40°至50°,呈“V”型谷。地勢西北高東南低,屬于典型高山峽谷地貌[23]。死水位為2 063 m,水庫正常蓄水位2 133 m,計(jì)劃2008年3月開工,于2011年11開始投入使用。庫區(qū)出露的基巖受“5.12”汶川特大地震[24]以及近年來強(qiáng)降雨等影響,加之庫區(qū)水位線短期快速的變動(dòng),觸發(fā)或加劇了庫岸的滑坡災(zāi)害[19]。

1.2 研究數(shù)據(jù)

本次研究采用數(shù)據(jù)為Sentinel-1衛(wèi)星數(shù)據(jù),Sentinel-1衛(wèi)星是歐洲航天局(European Space Agency,ESA)發(fā)射運(yùn)營的兩顆對地觀測衛(wèi)星組成。搭載波長為5.6 cm的SAR,重訪周期12 d,覆蓋范圍達(dá)到250 km×170 km[25],可實(shí)現(xiàn)全天時(shí)全天候監(jiān)測,在地質(zhì)災(zāi)害監(jiān)測方面被廣泛應(yīng)用[26-28]。

本研究使用了從2019年1月至2022年3月覆蓋研究區(qū)的97景升軌Sentinel-1A衛(wèi)星數(shù)據(jù),主要參數(shù)如表1所列。為了保證數(shù)據(jù)有良好的相干性,設(shè)置的時(shí)間基線為36天,空間基線為200 m,生成干涉對組合。結(jié)合歐空局的精密軌道數(shù)據(jù)減少參考橢球相位的影響。此外,數(shù)字高程模型(Digital Elevation Model,DEM)使用SRTM(Shuttle Radar Topography Mission),可以有效減弱地形相位的影響[29]。

表1 SAR數(shù)據(jù)主要參數(shù)

2 研究方法

本研究利用相位梯度對局部變形非常敏感的特點(diǎn),采用經(jīng)過改進(jìn)后的Sobel算子,對InSAR差分干涉對進(jìn)行梯度處理,獲取InSAR差分干涉對的相位信息。去除植被等低相干性區(qū)域、幾何畸變區(qū)域、水體、積雪等形變誤差的影響,通過梯度疊加的方式減弱噪聲和大氣的影響,獲取相位不連續(xù)區(qū)域。這些區(qū)域經(jīng)過閾值分割生成二值圖,最后通過邊緣提取快速識別滑坡隱患。具體步驟如圖2所示。

圖2 技術(shù)路線Fig.2 Technical route in this study

首先,是數(shù)據(jù)預(yù)處理階段。利用收集的Sentinel-1影像結(jié)合精密軌道數(shù)據(jù)生成SLC數(shù)據(jù),從中裁剪出研究區(qū)域并進(jìn)行數(shù)據(jù)配準(zhǔn)。結(jié)合SRTM-DEM去除干涉對中參考橢球相位和地形相位的影響,選用合適的濾波對干涉對進(jìn)行濾波處理。為了減少誤差的引入,需選取干涉質(zhì)量較好的干涉對,因此對濾波后的干涉對進(jìn)行精煉,剔除效果不好的干涉對。

其次,是梯度處理階段。利用改進(jìn)的Sobel算子對濾波、精煉后的干涉對進(jìn)行梯度運(yùn)算。普通的Sobel算子只采用水平和垂直兩個(gè)方向的梯度模板進(jìn)行檢測,雖然具有普適性和計(jì)算速度快等優(yōu)勢,但對紋理復(fù)雜的圖像邊緣特征識別不足,且無法連續(xù)檢測邊緣。因此,本文利用改進(jìn)后的Sobel算子即加入了45°和135°兩個(gè)對角方向的模板識別,使滑坡邊界檢測更為連續(xù)和精細(xì)[30-33]。由于在進(jìn)行梯度處理可能受到水體[34]、積雪[35]、衛(wèi)星的成像幾何帶來的幾何畸變[11]、植被相干性較低[36]的區(qū)域等影響,導(dǎo)致相位梯度不連續(xù),形成形變誤差,處理后呈現(xiàn)為相位梯度的突變,類似于滑坡隱患。因此需要對這些形變誤差進(jìn)行有效去除。本文結(jié)合目標(biāo)點(diǎn)相干性,利用相干系數(shù)選取合適的閾值,對受植被、噪聲等影響嚴(yán)重區(qū)域進(jìn)行掩膜處理。通過多次實(shí)驗(yàn),確定相關(guān)性閾值為0.7,可以有效去除植被、噪聲等影響。利用衛(wèi)星參數(shù)等信息生成研究區(qū)幾何畸變數(shù)據(jù),通過掩膜幾何畸變區(qū)域,去除幾何畸變的影響。收集毛爾蓋庫區(qū)水位情況與積雪線,通過DEM去除水體和積雪的影響。再利用梯度疊加50對干涉對識別出相位不連續(xù)區(qū)域。

然后,是快速識別階段。識別到的相位不連續(xù)區(qū)域仍然存在由于衛(wèi)星自身成像等因素引起的形變誤差,采用均值濾波對結(jié)果進(jìn)行平滑處理,如形變識別階段所示。本文通過設(shè)定一定閾值生成二值圖,進(jìn)行內(nèi)部填充后提取相位不連續(xù)區(qū)域邊界,考慮到發(fā)生滑動(dòng)的隱患點(diǎn)區(qū)域面積在一定的范圍內(nèi),而較大面積的區(qū)域可能由于大氣效應(yīng)引起的干涉相位的不連續(xù),因此再次通過閾值分割去除掉面積較大區(qū)域,最終得到真實(shí)形變區(qū)域。

最后,是對比驗(yàn)證階段。將本文識別結(jié)果與相應(yīng)時(shí)間段的傳統(tǒng)時(shí)序InSAR結(jié)果進(jìn)行對比驗(yàn)證,以確保識別結(jié)果的準(zhǔn)確性。同時(shí)總結(jié)本文方法對滑坡隱患識別的特點(diǎn)和能力。

3 毛爾蓋庫岸滑坡隱患識別結(jié)果

利用InSAR相位梯度疊加的滑坡隱患快速識別方法針對毛爾蓋庫區(qū)開展滑坡隱患識別的結(jié)果如圖3(a)所示,共識別出23處滑坡隱患點(diǎn)。綜合谷歌光學(xué)影像可以明顯看出,基于相位梯度疊加的滑坡隱患快速識別方法識別出的滑坡隱患,能夠看出較為明顯形變跡象,其滑坡邊界清晰,后緣可見明顯的下錯(cuò),說明該處已產(chǎn)生較大形變[見圖(c)、(e)、(a)、(i)];在空間分布上,這些識別的災(zāi)害隱患點(diǎn)主要分布于庫岸兩側(cè)[圖3(i)],這一現(xiàn)象意味著庫區(qū)水位的快速變動(dòng)導(dǎo)致土體物理特性下降,對庫岸坡體的穩(wěn)定性產(chǎn)生顯著影響;從威脅對象來看,庫區(qū)兩岸坡體可能出現(xiàn)失穩(wěn)情況,發(fā)生較大規(guī)模滑坡,從而出現(xiàn)堵江風(fēng)險(xiǎn),影響水電站正常運(yùn)行。此外,有5處隱患點(diǎn)附近分布有居民點(diǎn)及道路等基礎(chǔ)設(shè)施,如果發(fā)生垮塌,將對當(dāng)?shù)鼐用竦纳?cái)產(chǎn)安全造成巨大威脅。

(a為InSAR相位梯度疊加快速識別總體結(jié)果;b、d、f、h為典型區(qū)域放大結(jié)果;c、e、g、i為典型區(qū)域Google Earth光學(xué)影像)圖3 InSAR相位梯度疊加快速識別結(jié)果Fig.3 Rapid identification results of the InSAR phase gradient stacking

為了驗(yàn)證本文方法的精度,將快速識別的結(jié)果與2019年1月至2022年3月的SBAS-InSAR結(jié)果進(jìn)行對比驗(yàn)證。如圖4所示,獲取了升軌軌道視線方向上的形變,紅色區(qū)域(負(fù)值)表示遠(yuǎn)離雷達(dá)視線方向移動(dòng),藍(lán)色區(qū)域(正值)表示靠近雷達(dá)視線方向移動(dòng),綠色區(qū)域則表示相對穩(wěn)定。研究區(qū)幾何畸變區(qū)域較少,植被覆蓋不是特別茂密,且獲取數(shù)據(jù)時(shí)間范圍較長(2019—2022年),獲取的SBAS-InSAR結(jié)果精度較高。SBAS-InSAR技術(shù)共識別到24處滑坡隱患,從整體識別結(jié)果來看,識別出的形變區(qū)域主要分布在庫岸兩側(cè),其中區(qū)域C兩處滑坡形變量級較大,達(dá)到109 mm/年。基于InSAR相位梯度疊加的滑坡隱患快速識別方法共識別到23處滑坡隱患[圖4(a)中白色區(qū)域],與SBAS-InSAR方法所識別結(jié)果具有高度一致性。從精度驗(yàn)證結(jié)果可知本文提出的方法進(jìn)行的滑坡隱患識別準(zhǔn)確率較高。在對廣域突發(fā)性高的滑坡隱患進(jìn)行快速動(dòng)態(tài)識別時(shí)具有一定的潛力。

(a為SBAS-InSAR和InSAR相位梯度疊加快速識別總體結(jié)果;b～e 為典型區(qū)域放大結(jié)果)圖4 SBAS-InSAR和InSAR相位梯度疊加快速識別結(jié)果Fig.4 Rapid identification results of SBAS-InSAR and InSAR phase gradient stacking

將本文識別結(jié)果與傳統(tǒng)Stacking-InSAR結(jié)果進(jìn)行比較,對比結(jié)果如圖5所示。基于InSAR相位梯度疊加的滑坡隱患快速識別方法共識別到23處潛在滑坡隱患,而傳統(tǒng)Stacking-InSAR方法共識別到20處潛在滑坡隱患。總體來說兩種方法識別結(jié)果吻合度較高。為了更詳細(xì)的進(jìn)行比較,將兩種方法的能識別到的典型區(qū)域放大,具體如圖5(b)～(e)所示。這些區(qū)域的滑坡隱患表現(xiàn)為形變量級較大,滑動(dòng)跡象明顯,因此Stacking-InSAR方法和基于InSAR相位梯度疊加的滑坡隱患快速識別方法都能夠準(zhǔn)確地識別出來。另外,有3處滑坡隱患,其形變量級較小,滑動(dòng)跡象不明顯,僅由基于InSAR相位梯度疊加的滑坡隱患快速識別方法成功識別出來。這些對比結(jié)果證實(shí)了本文所提方法在滑坡隱患識別方面具有良好的表現(xiàn)。

4 討論

根據(jù)圖4和圖5分析,三種方法共同識別到20處滑坡隱患,這一高度一致的識別結(jié)果充分證明InSAR相位梯度疊加的滑坡隱患快速識別方法的可靠性。然而,不同方法的識別結(jié)果也存在一定差異。下面將該方法和傳統(tǒng)時(shí)序InSAR識別結(jié)果進(jìn)行了分類,以進(jìn)一步探究本文方法的特點(diǎn),我們將識別出的滑坡隱患分為兩大類(表2):

表2 識別的滑坡隱患分類

第一類是所有方法均能識別到的滑坡隱患,三種方法共同識別到20處滑坡隱患,圖6中放大了其中幾處典型區(qū)域,它們反映了第一類滑坡隱患的一般情況。由圖6可知,這幾處滑坡識別的形變速率達(dá)到60 mm/年以上,形變面積均大于400 m2。這類滑坡隱患形變量級大,滑動(dòng)跡象明顯,本文方法和傳統(tǒng)時(shí)序InSAR技術(shù)均能識別。

( a、b、 c為InSAR相位梯度疊加快速識別結(jié)果;d、 e、 f為Stacking-InSAR結(jié)果;g、 h、 i為SBAS-InSAR結(jié)果)圖6 三種技術(shù)均能識別的滑坡隱患典型結(jié)果Fig.6 Typical results of potential landslides that can be identified by the three methods

第二類滑坡隱患是僅由InSAR相位梯度疊加的滑坡隱患快速識別方法探測的滑坡隱患。如圖7所示(白色區(qū)域?yàn)镮nSAR相位梯度疊加的滑坡隱患快速識別結(jié)果),由圖7(a3、b3、c3)可知這幾處滑坡隱患形變速率均在30 mm/年以下,形變速率較小。同時(shí)這幾處滑坡隱患形變范圍均在300 m2以下,形變的范圍較小,滑動(dòng)跡象不明顯。因此傳統(tǒng)Stacking-InSAR難以探測到[14]。這幾處滑坡隱患卻被InSAR相位梯度疊加的滑坡隱患快速識別方法明顯探測到。其原因是InSAR相位梯度對形變具有高度敏感性,即使微小的形變也能引起相位梯度的變化。因此對于形變量級較小的滑坡隱患,基于InSAR相位梯度疊加的滑坡隱患快速識別方法有更好的識別效果。

值得注意的是基于InSAR相位梯度疊加的滑坡隱患快速識別方法,可以在不進(jìn)行相位解纏和時(shí)序建模解算的情況下直接求取InSAR差分干涉圖的相位梯度。這節(jié)省了大量數(shù)據(jù)處理的時(shí)間,提高了滑坡隱患識別效率。如表3所列,統(tǒng)計(jì)了三種技術(shù)從數(shù)據(jù)配準(zhǔn)到識別出滑坡隱患所需要的時(shí)間以及識別的滑坡隱患個(gè)數(shù)。本次實(shí)驗(yàn)統(tǒng)一采用CPU為Inter i7-11700;GPU為RTX2060,32 G內(nèi)存;4 T機(jī)械硬盤,傳輸速度100 M/s;在常溫24°環(huán)境下進(jìn)行。本此研究區(qū)面積約為373 km2,共使用了97景升軌哨兵影像。三種方法均能有效識別庫區(qū)內(nèi)大多數(shù)滑坡隱患。基于InSAR相位梯度疊加的滑坡隱患快速識別方法識別出23處滑坡隱患,用時(shí)2h50min,耗時(shí)最短,識別速度相較Stacking-InSAR、SBAS-InSAR分別提升了1.4和1.9倍(表3)。與傳統(tǒng)時(shí)序InSAR相比,基于InSAR相位梯度疊的滑坡隱患快速識別方法具有耗時(shí)短、操作簡單、精度高等特點(diǎn)。可以有效縮短數(shù)據(jù)處理的時(shí)間,提高滑坡隱患的識別效率,在廣域滑坡隱患快速動(dòng)態(tài)識別中具有一定的應(yīng)用潛力。

表3 三種方法識別滑坡隱患所耗時(shí)間

5 結(jié)論

本文提出一種基于InSAR相位梯度疊加的滑坡隱患快速識別方法,結(jié)合Sentinel-1數(shù)據(jù)針對毛爾蓋庫岸滑坡隱患進(jìn)行快速識別。研究結(jié)果表明:

(1) 基于InSAR相位梯度疊加的滑坡隱患快速識別方法識別出23處滑坡隱患,和Stacking-InSAR與SBAS-InSAR方法的識別結(jié)果對比驗(yàn)證,三種不同方法共同識別到了20處滑坡隱患,這一高度一致的識別結(jié)果充分證明了本文方法在滑坡隱患識別中的可靠性。

(2) 本文提出的基于InSAR相位梯度疊加的滑坡隱患快速識別方法成功識別出Stacking-InSAR無法識別的3處發(fā)生微弱形變的滑坡隱患,該方法對變形量級小的滑坡隱患具備更好地探測能力。

(3) 本文提出的基于InSAR相位梯度疊加的滑坡隱患快速識別方法由于直接對差分干涉相位進(jìn)行梯度疊加運(yùn)算,避免了相位解纏和時(shí)序建模解算等步驟,識別速度相較Stacking-InSAR、SBAS-InSAR方法分別提升1.4和1.9倍,極大地縮短數(shù)據(jù)處理時(shí)間,從而有效提升了滑坡隱患識別的效率。

綜上所述,本文提出的基于InSAR相位梯度疊加的滑坡隱患快速識別方法耗時(shí)短、效率高,相較于傳統(tǒng)時(shí)序InSAR,能更清晰探測出微弱變形的滑坡隱患。在廣域滑坡隱患快速動(dòng)態(tài)的識別中具有一定的應(yīng)用潛力,為InSAR技術(shù)在廣域滑坡隱患快速識別中的應(yīng)用提供重要參考。