基于時序InSAR技術的滑坡隱患識別對比研究

董繼紅， 馬志剛， 梁京濤， 劉 彬， 趙 聰， 曾 帥， 鄢圣武， 馬曉波

(1.四川省地質調查院，稀有稀土戰(zhàn)略資源評價與利用四川省重點實驗室，成都 610081； 2.四川省智慧地質大數據有限公司，成都 610081； 3.四川省國土空間生態(tài)修復與地質災害防治研究院，成都 610081)

0 引言

滑坡作為頻繁發(fā)生、破壞強烈的全球自然災害之一，對于人民生命財產安全造成了巨大的威脅，同時也制約著當地經濟的發(fā)展[1-2]。四川省位于中國西南山區(qū)，地質環(huán)境條件復雜，構造發(fā)育，是地質災害高易發(fā)區(qū)，同時也是我國地質災害防治工作的重點區(qū)域[3]。川西片區(qū)不僅地質構造復雜，活動斷裂和強震也極為發(fā)育，被多條斷裂帶穿越，同時該區(qū)域海拔變化劇烈，降雨充沛，使得該區(qū)域地質災害尤為頻發(fā)。2017年6月24日，四川省茂縣疊溪鎮(zhèn)新磨村發(fā)生高位隱伏滑坡，造成13人死亡，73人失蹤[4-6]； 2020年8月21日雅安市漢源縣富泉鎮(zhèn)發(fā)生山體滑坡，造成7人死亡2人失蹤。

面對頻發(fā)的地質災害，合成孔徑雷達干涉測量技術(interferometric synthetic aperture Radar，InSAR)作為一種大范圍地表形變監(jiān)測的新技術，因其非接觸、大范圍、空間覆蓋范圍廣、監(jiān)測精度高等優(yōu)勢，已經廣泛用于地震、地面沉降、冰川、地裂縫、滑坡等地質災害調查與監(jiān)測研究中，并取得了很好的效果[7-10]。趙超英等[11]基于InSAR技術對甘肅省黑方臺區(qū)域的黃土滑坡利用多源合成孔徑雷達(synthetic aperture Radar,SAR)數據進行滑坡編目，長時間序列監(jiān)測，并以新塬2號滑坡體為例進行失穩(wěn)模式識別研究，將識別與監(jiān)測采用的方式方法進行區(qū)分與介紹; 劉星洪等[12]利用光學遙感解譯技術與InSAR技術對巴塘—芒康段的地質災害進行了調查，并通過野外地質調研、地理信息系統(tǒng)(geographic information system,GIS)空間分析和工程類比等工作，證明綜合遙感技術方法在青藏高原高山峽谷區(qū)的應用，文中通過對光學與InSAR結果對比進行分析，討論了InSAR技術對隱患的識別效果; 韓冬建等[13]利用小基線集(small baseline subset,SBAS)技術和干涉圖堆疊技術基于不同波段的SAR數據對樟木口岸在尼泊爾地震引發(fā)的山體滑坡進行識別和形變特征研究，成功識別出了13處滑坡隱患，文中只是對不同波段的SAR數據識別能力進行了對比，其中SBAS技術只是用來監(jiān)測典型滑坡隱患，并沒有探討SBAS技術對滑坡隱患的識別能力差異性。

Stacking技術利用對解纏圖進行相位堆疊，可以快速獲取地表形變速率，廣泛用于滑坡隱患的早期識別，SBAS技術通過時間域和空間域的解算，可以有效消除對地觀測過程中各種誤差對結果的影響，在獲取地表形變速率的同時可以獲取地表形變點長時間序列監(jiān)測信息。因此本文選擇在植被密集區(qū)域，開展不同時序InSAR技術在滑坡隱患識別數目和準確率方面對比研究，探討這2種技術用于滑坡隱患識別效果的優(yōu)劣，并選擇典型滑坡隱患變形特征進行對比分析。

1 研究區(qū)概況及數據源

1.1 研究區(qū)概況

雅安位于四川盆地與青藏高原的過渡地帶，其西側邊緣地帶的山地海拔高達5 000余m，東側低山丘陵地帶海拔僅1 000 余m，區(qū)內地形高差極大。龍門山斷裂帶、鮮水河斷裂帶、安寧河斷裂帶三大活動活動斷裂帶均穿過雅安地區(qū)，導致區(qū)內構造活動強烈，地震頻發(fā)。同時，雅安是我國內陸降雨最為充沛的地區(qū)之一，素有“雨城”之稱。極為復雜的地質條件以及充沛的降雨構成了雅安地區(qū)脆弱的地質環(huán)境，加之日益增多的人類工程活動，使得雅安成為全國地質災害最易發(fā)的區(qū)域之一[14-15]。

1.2 數據源

本文收集了覆蓋研究區(qū)的Sentinel-1A衛(wèi)星SAR影像數據，其覆蓋范圍如圖1所示。由于研究區(qū)域面積較大，升軌數據需要4個軌道數據(圖1綠色框)，降軌數據需要2個軌道數據(圖1藍色框)可以完全覆蓋研究區(qū)域，時間范圍為2018年10月—2020年10月，單個軌道是91景數據。表1列出了本文研究中所使用的SAR數據集的基本參數。

圖1 研究區(qū)域示意圖Fig.1 Schematic diagram of the study area

表1 Sentinel-1A衛(wèi)星參數與數據選用時間范圍Tab.1 Sentinel-1A satellite parametersand data selection time range

使用AW3D30數字表面模型來消除InSAR干涉處理過程中的地形相位，同時該數據也被用來輔助SAR影像進行地理編碼及計算SAR數據的疊掩與陰影區(qū)域。Sentinel-1衛(wèi)星的精密軌道星歷(precise orbit ephemerides,POD)數據被用來輔助Sentinel-1數據的預處理和基線誤差改正[9]。

同時在數據處理中對于與高程不相干的隨機大氣擾動誤差，本文利用李振洪教授團隊研制的通用型衛(wèi)星雷達在線大氣改正系統(tǒng)(generic atmospheric correction online service for InSAR, GACOS)來去除[16]。

2 時序InSAR技術數據處理流程及結果分析

φint=φfla+φtop+φdef+φnois+φatm+φorb，

(1)

式中：φfla為平地相位；φtop為DEM誤差引入的殘余地形相位；φdef為地表形變相位；φnois為噪聲相位；φatm為大氣延遲相位；φorb為軌道誤差引起的相位。因此最終想從φint中獲取到形變相位φdef，需要經過相位解纏及一系列的誤差去除處理。

基于國內外研究[17-21]并結合本文的研究內容，本文采用SBAS技術與相位堆疊技術(Stacking)進行滑坡隱患識別對比分析研究。SBAS技術是由Berardino等[22]于2002年提出的，是針對覆蓋同一區(qū)域的多幅SLC(single look complex)數據，按照規(guī)范設置一定的垂直基線閾值和時間基線閾值進行組合，通過這種形式可以有效降低因垂直基線和時間基線過大引起的失相干現象，同時利用空間濾波提高相干性，通過犧牲分辨率提高信噪比的方式進行多視處理，進一步減弱失相干對地表監(jiān)測的影響； 利用差分干涉計算獲得差分干涉圖和相干圖，對差分干涉圖進行相位解纏，并對解纏相位進行時間域和空間域的形變量估算，減去趨勢誤差； 通過時空域濾波來去除大氣誤差的影響，從而更高精度地獲取相對于第一幅SAR影像的累積形變時間序列及平均形變速率。

Stacking技術是由Sandwell等[23]在1998年提出的，該技術原理是對多幅干涉圖進行加權平均解算，以達到削弱空間上不相關的噪聲的影響，在本文的數據處理中對基線精化之后的解纏圖進行相位堆疊解算，通過地理編碼獲取其在地理坐標系下的平均形變速率圖[24]。圖2為本文數據處理流程圖。通過不同時序InSAR技術獲取的地表形變速率圖，依據形變閾值，結合DEM、解纏圖、差分干涉圖及光學影像等信息進行滑坡隱患的早期識別解譯。

圖2 時序InSAR技術數據處理流程Fig.2 Time-series InSAR technology data processing flow

2.1 基于Stacking InSAR技術地表形變監(jiān)測結果分析

在雅安市數據處理中，在距離向和方位向采用10×2的多視比，垂直基線不大于150 m，時間基線不大于48 d進行干涉對組合，然后對所有干涉對開展差分干涉、自適應濾波、相位解纏處理及大氣改正。解纏方法采用最小費用流法，解纏參考點選擇在相干性高、且穩(wěn)定的區(qū)域，利用GACOS在線大氣改正系統(tǒng)對解纏圖進行殘余大氣延遲誤差改正[25]，最后對解纏圖進行篩選，篩選標準依據相干性較好、解纏誤差較少，對篩選之后的解纏圖進行Stacking解算，獲取了覆蓋研究區(qū)域沿雷達視線向的形變速率圖，如圖3所示。從圖3中可以看出升軌數據結果在石棉縣、滎經縣出現較多空白區(qū)域，在整個實驗區(qū)的左半部分出現較多的解纏誤差，主要是為了讓植被更加密集的區(qū)域(滎經縣、漢源縣、石棉縣等)可以獲取更多的觀測信息，通過對覆蓋該區(qū)域的平均相干圖的相干性進行查看，選擇大于相干性閾值面積占比80%以上的值作為相干性閾值，對低于相干性的區(qū)域進行掩模處理。然而在這些區(qū)域相干性過低，使得存在較多的解纏誤差； 因為研究區(qū)域不能被單景Sentinel-1數據所覆蓋，這也就使得兩景之間影像邊界處存在顏色的差異性，后期在圖幅拼接的時候未做校正處理，但是對滑坡隱患的識別、解譯并不影響。將石棉縣區(qū)域結果放大，可以很明顯發(fā)現因升降軌數據衛(wèi)星飛行方向不同，石棉縣等地區(qū)域的結果出現明顯的差異性，正是由于升降軌數據收集處理分析，避免了單個軌道數據對隱患點的漏判。

2.2 基于SBAS InSAR技術地表形變監(jiān)測結果分析

在SBAS數據處理過程中，參數設置與Stacking計算過程設置基本一致。通過時間域和空間域的解算獲得了沿衛(wèi)星雷達視線向方向的地表形變速率結果(圖4)，其中負值表示遠離衛(wèi)星飛行方向，正值表示接近衛(wèi)星飛行方向。

通過解譯圖4發(fā)現升軌數據SBAS結果中存在較大的空白區(qū)域，主要集中在寶興縣、天全縣和滎經縣，結合該區(qū)域植被覆蓋情況，分析認為這與地表植被覆蓋情況相關，造成這些區(qū)域在時間域上相干性不連續(xù)造成空缺。同時一些區(qū)域存在部分解纏誤差，分析原因是在相位解纏過程中為了顧及更多的信息，基于相干性設置的掩模閾值過低所導致。

在對高速公路的橋梁進行養(yǎng)護的過程中，其實不僅僅只針對于橋梁和高速公路本身，還要兼顧高速公路周圍的生態(tài)環(huán)境以及生活服務，因為它們是一個整體，不可割裂開來，只有做到全面養(yǎng)護才能保證高速公路橋梁的持久利用。

通過對比Stacking結果與SBAS結果發(fā)現，在SBAS結果中山脊處存在部分空白區(qū)域，以寶興縣降軌結果表現最為明顯，主要原因是在SBAS處理過程中剔除了疊掩、陰影區(qū)域。同時升軌數據SBAS結果較升軌數據Stacking結果有較大范圍空白，主要是由于空白區(qū)域在時間域上相干性不連續(xù)所造成的。SBAS結果明顯好于Stacking數據結果，分析原因主要是Stacking技術只是簡單的進行相位加權平均，對一些殘余大氣誤差、DEM誤差等未進行去除，而SBAS技術通過時間域和空間域的解算，能有效地消除或削弱解纏粗差、大氣誤差以及DEM誤差等因素的影響，因此后者結果優(yōu)于前者。

通過以上對比分析，從地表形變監(jiān)測結果和處理效果來分析，在雅安市SBAS技術獲取結果明顯優(yōu)于Stacking技術獲取結果。

3 時序InSAR結果對比分析

3.1 漢源縣區(qū)域地表形變監(jiān)測結果分析

選取雅安市漢源縣區(qū)域單獨展示對比分析。使用升軌數據獲取了其年平均形變速率圖(圖5)。發(fā)現Stacking結果與SBAS結果具有較大的差異性，從Stacking結果中并不能看到明顯形變信息，結合SBAS解譯的結果在Stacking結果中識別出2處形變區(qū)域，并且主要集中在城鎮(zhèn)區(qū)域，且形變信息不明顯，受大氣影響嚴重。在山區(qū)出現許多嘈雜點和顏色拉伸現象，如圖5(a)中黃色方框所示，分析原因一方面是這些區(qū)域相干性差，解纏閾值過低，存在解纏誤差； 另一方面是疊掩現象的存在，經地理編碼之后，疊掩區(qū)域出現拉伸現象。圖5(b)為利用SBAS技術獲取的結果，可以看到整個結果干凈，但是前期處理中設置的相干性掩模閾值過低，因為剔除了一些誤差，所以結果中沒有出現左圖對應位置的嘈雜點，但是在山區(qū)還是表現出一些效果差的區(qū)域； 同時對疊掩陰影區(qū)域的剔除，使得SBAS結果中在一些地形起伏區(qū)域出現空白，但是結果更加精確。基于SBAS結果初步可以解譯出5處明顯的疑似隱患形變點(紅色圓圈所示)，中間圓圈所圈的3處隱患點，一處形變信息與另外兩處相反，主要是該結果為雷達視線向形變信息。

選擇圖5(b)左邊第一個隱患點，即紅色箭頭所示，利用SBAS InSAR技術進行時序分析，獲取了該隱患點在時間域上的累積形變曲線圖，如圖6所示。

圖6 疑似隱患點累積形變Fig.6 Cumulative deformation of suspectedhidden danger points

由圖6可知，該點從2018年10月—2020年10月一直處于不斷變形狀態(tài)，并且目前處于加速變形期，累積形變量超過了0.07 m，佐證了此點為隱患點。結合圖3—6對比分析發(fā)現在雅安區(qū)域SBAS技術更加適合于滑坡隱患的早期識別及監(jiān)測工作。

3.2 典型滑坡隱患識別對比分析

選取一經野外驗證的滑坡隱患點進行不同時序InSAR技術識別對比分析，該點位于寶興縣隴東鎮(zhèn)先鋒村。初步通過利用升軌數據基于Stacking技術及SBAS技術處理獲取了形變速率圖，結合光學影像進行解譯，獲取了該滑坡的邊界信息，如圖7所示。從SBAS結果中可以發(fā)現明顯形變信息，形變區(qū)與滑坡范圍基本一致，最大形變速率超過了-80 mm/a，最大形變區(qū)域位于滑坡體左側及前部區(qū)域，同時將解譯的邊界信息疊加至Stacking結果上，發(fā)現在Stacking結果中變形信息并不明顯，與SBAS結果變形速率最大區(qū)域所對應位置一致。對比可知該隱患點的SBAS技術識別效果明顯優(yōu)于Stacking技術的識別效果，SBAS結果中形變信息空缺區(qū)域是因為在時間域上相干性不連續(xù)所造成的。因為該隱患點在降軌數據中處于疊掩區(qū)域，故沒有展示降軌結果。

圖7 典型滑坡隱患解譯信息

從光學影像上可知解譯區(qū)植被分布茂密，可見多處明顯小規(guī)模溜滑跡象，土地裸露。滑坡呈舌狀，兩側以沖溝為界，前緣臨溝，外凸較明顯，無新近堆積，后緣可見早期滑坡壁，坡體上分布有較多農戶。

通過野外調查發(fā)現該點為已知點滑坡隱患： 唐包滑坡，地理位置為E102°42′53.1″，N30°26′52.7″。滑坡變形區(qū)域主要集中在滑坡前緣右側和中部，變形強烈區(qū)域與InSAR識別結果吻合性較好(圖8)。圖8(a)為滑坡上半部分照片輪廓，從圖中可以看到植被發(fā)育茂密，植被以灌木為主，喬木次之，可以看到多處滑塌區(qū)域，土地裸露。P1位于滑坡前緣右側，附近居民點院壩開裂長約10 m，寬3～4 cm，房屋也見多條裂縫存在，寬約5～6 cm不等； 公路錯開25 cm，裂縫寬約30 cm，長約5 m，局部下陷25 cm。P2位于滑坡前緣中部，地面可見網狀拉裂，寬3～8 cm，公路外側見明顯裂縫，路面破壞長度約30 m。

圖8 典型滑坡野外調查照片

3.3 雅安市滑坡隱患識別對比分析

基于不同時序InSAR技術獲取的結果進行滑坡隱患識別解譯，并通過野外進行驗證結果的準確性，統(tǒng)計了雅安市不同InSAR技術識別滑坡隱患分布情況(圖9)。如圖9所示，圖中藍色“共同識別”為基于Stacking技術和SBAS技術均識別出來，并經過野外調查驗證為正確的隱患點； “SBAS識別正確”為單獨由SBAS技術識別出來并經野外調查驗證為正確的隱患點，“SBAS識別錯誤”為僅由SBAS技術識別并經野外調查驗證為錯誤的隱患點； “Stacking識別正確”為單獨由Stacking技術識別出來并經野外調查驗證為正確的隱患點，“Stacking識別錯誤”為僅由Stacking技術識別并經野外調查驗證為錯誤的隱患點。

圖9 雅安市滑坡隱患識別分布圖Fig.9 Distribution map of landslide hiddendanger identification in Ya’an City

對識別的隱患點進行統(tǒng)計分析(圖10)。 可知基于Stacking技術共識別出隱患點83處，經野外驗證、核查，共計有48處隱患點被識別出來。Stacking技術識別準確率為57.8%。對識別錯誤的隱患點進行分析發(fā)現，造成隱患識別錯誤的原因一方面是植被茂密，相干性較差，解纏過程中為了獲取更多的信息設置較低的相干性掩模閾值，使得在結果中有較多的解纏誤差存在； 另一方面是Stacking技術只是單純對解纏圖進行加權平均，結果中存在較多的殘余誤差。

圖10 雅安市InSAR技術識別地質災害統(tǒng)計圖Fig.10 InSAR technology in Ya’an City toidentify geological hazards statistical map

基于SBAS技術共識別出隱患點54處，經野外驗證、核查，共計有38處隱患點被正確識別出來。SBAS技術識別準確率為70.3%。對識別錯誤的隱患點進行分析發(fā)現，造成隱患識別錯誤的原因一方面植被茂密，相干性較差，造成結果中有一些異常值； 另一方面由于干涉圖在時間域相干性不連續(xù)，對一些隱患點不能提取時序分析判識，造成誤判，以及由其他因素引發(fā)的變形，經野外核查不能歸為滑坡隱患。

基于Stacking技術和SBAS技術共同識別的隱患點為44處，經野外驗證、核查，共計有38處隱患點被兩種技術均被識別出來，識別準確率為86.4%。

經過上述討論可知，在植被密集區(qū)域，地形條件復雜，基于InSAR技術進行滑坡隱患識別，影響和限制因素較多，在顧及識別隱患點數目、準確率下，建議SBAS技術與Stacking技術相結合的形式用于該區(qū)域地質災害隱患識別中。

4 結論

本文選擇在植被密集區(qū)開展利用2018年10月—2020年10月升降軌Sentinel-1數據進行基于Stacking技術和SBAS技術滑坡隱患早期識別對比分析實驗，以雅安市為研究區(qū)。并選取典型滑坡體結合野外驗證資料進行時間序列分析。主要有以下結論：

1)利用Sentinel-1數據基于時序InSAR技術用于雅安市滑坡隱患早期識別具有較強的可行性，通過獲取的結果經對比分析發(fā)現升降軌數據共同監(jiān)測可以有效避免滑坡隱患的漏判，增加觀測信息。

2)通過統(tǒng)計發(fā)現基于Stacking技術識別的滑坡隱患數目最多，結合野外驗證資料發(fā)現基于SBAS技術識別的滑坡隱患準確率最高。

3)通過對時序InSAR技術識別的滑坡統(tǒng)計分析發(fā)現，在顧及識別隱患點數目、準確率等情況下，建議在植被密集區(qū)域采用SBAS技術與Stacking技術相結合的形式開展滑坡隱患識別。