唐光民, 戴可人,,3, 卓冠晨,, 沈 月,4, 陳 晨, 許 強

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

0 引言

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

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

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

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

1.1 研究區(qū)域

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

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

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

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

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

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

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

2 研究方法

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

圖2 技術路線Fig.2 Technical route in this study

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

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

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

最后,是對比驗證階段。將本文識別結果與相應時間段的傳統(tǒng)時序InSAR結果進行對比驗證,以確保識別結果的準確性。同時總結本文方法對滑坡隱患識別的特點和能力。

3 毛爾蓋庫岸滑坡隱患識別結果

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

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

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

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

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

4 討論

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

表2 識別的滑坡隱患分類

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

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

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

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

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

5 結論

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

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

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

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

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