基于SBAS-InSAR技術的白格滑坡形變監測研究

(1.山東理工大學 建筑工程學院，山東 淄博 255049； 2.北京科技大學 土木與環境工程學院，北京 100083；3.燕山大學 建筑工程與力學學院，河北 秦皇島 066064)

1 研究背景

我國是地質災害多發國家，作為典型地質災害，滑坡每年都給國家和人民造成了嚴重危害。2018年10月11日07：10，西藏昌都市江達縣波羅鄉白格村發生山體滑坡，造成金沙江斷流并形成堰塞湖，所幸未造成人員傷亡。11月3日再次發生滑坡，此次滑坡點與原“10·11”山體滑坡點位置相同，滑坡后金沙江干流被阻斷，導致上游水位持續上漲，形成堰塞湖，堰塞湖水位累計上漲57.44 m，推算堰塞湖蓄水量約為4.69億m3，轉移群眾2.5萬人?？梢妼聻暮M行有效的監測預警至關重要，而滑坡的早期識別及趨勢監測是滑坡災害監測的重要內容之一。傳統的滑坡地表形變監測方法需先知道滑坡在哪，其主要方法有大地測量方法、GPS法、自動伸縮計法和分布式光纖法等。這些監測方法都是“點”式數據采集方法，無法全面完整地反映區域內滑坡變形的整體面狀狀況，不能應用于人員難以進入的區域、處于高風險和尚未形成的滑坡地段，且往往容易造成“測者未災，災者未測”的現象。

合成孔徑雷達差分干涉測量技術(D-InSAR) 是最近幾十年來發展起來的可以監測地表形變的一種遙感技術。1989年A.K.Grabriel等首次驗證了D-InSAR技術可監測厘米級精度地表形變。早期研究發現，該技術主要應用于形變量比較明顯的地震、火山活動等監測，后期隨著理論方法的成熟和技術研究的深入，研究重點逐漸轉移至地面沉降、山體滑坡等細微的地表位移監測。該技術具有全天時、全天候、遠距離、面狀、厘米甚至毫米級精度的優勢，目前已應用于各個領域。劉云華等采用Sentinel-1A及ALOS-2為數據源，獲得了阿克陶地震震源參數反演及發震構造[1]。李達等基于SBAS-InSAR技術提取出各觀測點的時序沉降值，并對工作面時序沉降進行量化分析,證明了SBAS-InSAR技術在礦區地表沉降監測與分析方面具有良好的應用前景[2]。李珊珊等采用InSAR技術獲得了青藏高原季節性凍土形變趨勢[3]。張金芝等將SBAS時序分析技術應用于現代黃河三角洲地面沉降監測，提取了1992～2000年現代黃河三角洲地面沉降速率，結果與水準觀測數據基本保持一致[4]。我國也已開展了利用D-InSAR技術監測滑坡等地質災害形變的研究，但起步較晚，大多集中于對D-InSAR技術的方法、原理及其在滑坡地表形變監測應用潛力的分析研究?？祦喌炔捎眯』€集InSAR技術對金坪子滑坡進行監測，不僅獲得該滑坡的空間分區特征，也獲取重點滑坡區的時間序列結果，并且與地面監測結果一致[5]。Zhao C等以ALOS/PALSAR為數據源對加利福利亞北部區域進行廣域滑坡監測[6]。王立偉基于D-InSAR技術分析了高山峽谷區域滑坡位移識別技術路線與方法研究[7]。上述研究證明了利用D-InSAR技術監測滑坡形變的可行性。

然而傳統D-InSAR技術更側重于研究拍攝時間相隔較短的形變，且容易受到空間、時間失相干以及大氣等誤差因素的影響，無法獲取時間上連續的地面沉降場。王桂杰等通過ALOS衛星PALSAR傳感器獲得了三景合成孔徑雷達數據,再利用D-InSAR技術獲得研究區域內地表高精度形變位移值，并對不同方法的適用性進行了探討[8]。王平等利用D-InSAR測量和高密度電阻率剖面揭示了焦作市王封煤礦老采空區地面沉降機制[9]。劉曉菲等采用D-InSAR技術監測老采空區殘余變形[10]。但都因受到傳統D-InSAR技術缺點的限制，精度有所影響。為了克服傳統D-InSAR技術的缺點，一些高級D-InSAR技術逐漸出現，如小基線子集法(SBAS)和永久散射體干涉測量技術(PS-InSAR)。

本次研究以白格滑坡為研究區域，6景災前ALOS-2衛星PALSAR影像為數據源，采用SBAS-InSAR技術對數據進行處理，最終獲得監測時間段內的年平均沉降速率和累計沉降值。從結果可知，白格滑坡在發生之前已有較長時間的緩慢移動過程，再次驗證了以SAR數據為數據源，采用SBAS-InSAR技術可以有效地對高山峽谷區滑坡進行監測，實現發現變動和監視趨勢的兩個目的。

2 數據介紹

ALOS-2衛星是ALOS-1的后續星，于2014年5月24日發射升空，是唯一一個利用L波段頻率的高分辨率機載合成孔徑雷達的衛星，波段越長抗干擾性越強。它能很好地用于監測地殼運動和地球環境，且不受氣候條件和時間的影響。相比ALOS-1，JAXA在傳感器的設計和制造上有了更進一步的改善，更短的重訪周期(14 d)，雙向側擺的觀測能力，可以提供3種工作模式：Spotlight模式，Stripmap模式，ScanSAR模式。以6期災前ALOS-2衛星SAR影像為數據源(2015年7月27日至2018年7月23日)，拍攝模式為Stripmap模式，數據類型為FBD雙極化類型(HH+VV)，因自然地物對HH極化電磁波會產生較強的回波信號，地質災害形變監測優選選用HH極化模式。拍攝范圍如圖1所示，圖中紅色圈出區域為影像覆蓋范圍，綠色所圈區域為滑坡范圍。

圖1 研究區域Fig.1 Map of study area

3 數據處理流程

常規D-InSAR技術受時間和空間失相干的影響，其運用受到較大程度限制。PS-InSAR[11-12]方法雖然可有效解決失相干和大氣延遲等問題，但對SAR 影像數目要求較多，一般需要15景以上。而SBAS-InSAR方法[13-14]對SAR 影像數目要求相對來說不是很嚴格，在數據量相對有限情況下，該方法可優先選擇。但SBAS-InSAR方法獲得的形變監測結果精度較PS-InSAR技術要低，數據量至少需要6景影像。SBAS-InSAR方法通過連接具有長基線距且相互獨立的SAR影像，形成短基線 SAR 影像集合，從而增加數據獲取的采樣率，可在已有的SAR影像數據集中形成若干小集合，每個小集合內SAR影像間的基線較小，集合間 SAR影像的基線較大。SBAS-InSAR方法流程如圖2所示，具體處理流程如下所述。

3.1 生成連接圖

因數據范圍很大，而滑坡區域只是一小部分，因此對SLC數據進行了裁剪處理，從而減少數據，節省處理時間。裁剪后對所有圖像建立對應關系，每對主從影像空間基線距設為臨界基線距的45%，時間基線距設為1 000 d，共得到15對像對。圖3是生成的數據配對圖表，圖中像對用線連接，輸入的數據用點表示。每個像對的空間基線距和時間基線距如圖4所示。

圖2 SBAS-InSAR處理流程Fig.2 Process flow chart of SBAS-InSAR

圖3 數據配對Fig.3 Paired data

圖4 數據對時間和空間基線距Fig.4 Time baseline and spatial baseline of paired data

3.2 干涉工作流生成

根據像對的連接關系，對每一對像對進行干涉工作流處理(干涉圖生成、去平、濾波、相干性計算、相位解纏)，生成一系列解纏之后的相位圖。以2018年5月28日和7月23日兩景影像形成的像對為例，經干涉圖生成、去平、濾波、相干性計算、相位解纏，生成的相干系數如圖5所示，相干系數越大表示擾動性越小。相位解纏后干涉條紋圖如圖6所示。所有的干涉條紋圖最終都與超級主影像進行了配準，為下一步軌道精煉、重去平以及SBAS-InSAR的反演做準備。

圖5 相干系數圖Fig.5 Coherent coefficient

圖6 相位解纏后的條紋圖Fig.6 Interferogram after phase unwrapping

3.3 軌道精煉和重去平

這一步是估算和去除殘余的恒定相位和解纏后還存在的相位坡道。其中至關重要的是GCP的選擇，選擇的標準為：沒有殘余地形條紋；沒有形變條紋，遠離形變區域，除非已知這個點的形變速率；沒有相位躍變，如果GCP點位于一個孤立相位上，并且解纏的值非常差，這個位置可能是斜坡相位的一部分，那么選擇的這個GCP是不對的。由于在SBAS中很難找到完美的GCP可以全部用在所有的數據對中(因為數據對擁有不同的相干性)，因此建議多選擇一些GCP，至少20～30個點。

3.4 第一步反演

這一步是SBAS-InSAR反演的核心步驟，第一次估計位移速率和殘余地形，用來對合成的干涉圖進行去平，重新作相位解纏和精煉處理，二次解纏，生成更優化的結果，用于下一步計算。

3.5 第二步反演

這一步的核心是計算時間序列上的位移，在第一步得到的形變速率基礎上，進行定制的大氣濾波，從而估算和去除大氣相位，得到更加純凈的時間序列上的最終位移結果。

3.6 地理編碼

將SBAS-InSAR反演的結果進行地理編碼，獲得衛星拍攝視線向平均沉降速率和累計沉降值，如圖7～8所示。

圖7 年平均沉降速率(單位：cm/a)Fig.7 Average annual settlement rate

圖8 累計沉降值(單位：cm)Fig.8 Cumulative settlement value

4 結果分析

從圖7中可看出滑坡體后緣呈現下降趨勢，年平均沉降速率介于13～31 cm/a之間?；麦w中部和前緣部分區域都有抬升，抬升速率介于3～21 cm/a。從圖8中可知在2015年7月27日至2018年7月23日3 a時間內，滑坡后緣累計沉降值為81～155 cm，中部和前緣累計抬升值為6～52 cm。在圖8中選取A處進行時序分析。歷史形變見圖9，可以看出，A點累計沉降值為75 cm，2015年7月27日至2016年7月25日1 a時間沉降值較大，2016年7月25日至2018年5月28日形變較緩，而2018年5月28日至2018年7月23日約兩個月時間形變速率明顯加快。B點歷史形變如圖10所示，累計沉降值為70 cm，整體形變趨勢與A點相似，都為前一年形變較大，中間趨緩，后兩個月形變速率明顯加快。C點歷史形變圖如圖11所示，形變趨勢為先下降后抬升，最終抬升約18 cm，推測原因為滑坡后緣持續滑移導致前部抬升。

圖9 A點歷史形變Fig.9 The history deformation tendency of the Point A

圖10 B點歷史形變Fig.10 The history deformation tendency of the Point B

圖11 C點歷史形變Fig.11 The history deformation tendency of the Point C

5 結 語

本文以白格滑坡為監測對象，采用SBAS-InSAR技術進行沉降監測，獲得白格滑坡災前年平均沉降速率和累計沉降值圖。從監測結果中可知，白格滑坡在發生前已有較大形變，且形變趨勢為先急后緩，在滑坡災害發生前數月形變速率明顯加快。沉降趨勢與實際結果一致，證明了SBAS-InSAR技術用于滑坡早期識別與動態監測的可行性，在廣域地質災害監測具有較廣的發展應用前景。因未對研究區域進行災前持續遙感監測，所以監測結果僅僅是白格滑坡發生后為獲得滑坡發生前的形變趨勢而做得工作，但滑坡的監測預警更至關重要，后期會進一步采用SAR遙感的多期差分干涉持續監測，通過發現災害、災害識別及趨勢分析等關鍵內容的研究，建立發現滑坡變動范圍，監視形變趨勢的理論方法。

受數據量的影響，本次研究未采用精度更高的PS-InSAR技術，將來隨著數據量的累積，擬采用兩種方法相結合進行綜合分析研究，進一步提高監測精度。

致 謝

本文所采用的ALOS-2衛星PALSAR數據由天下圖公司提供。