基于干涉點目標分析方法的大光包滑坡形變監測

于 冰 王 楊 馬德英 蔣榮乾 張 過 周志偉

1 西南石油大學土木工程與測繪學院，成都市新都大道8號，6105002 西南石油大學油氣藏地質及開發工程國家重點實驗室，成都市新都大道8號，6105003 西南石油大學油氣空間信息工程研究所，成都市新都大道8號，6105004 武漢大學測繪遙感信息工程國家重點實驗室，武漢市珞喻路129號，4300795 中國科學院精密測量科學與技術創新研究院大地測量與地球動力學國家重點實驗室，武漢市徐東大街340號，430077 6 中鐵二院工程集團有限責任公司，成都市通錦路3號，610031

干涉點目標分析(interferometric point target analysis, IPTA)方法[1-2]是對PS-InSAR技術的發展，其通過建立相位模型，計算并分離PS點的各相位分量，再根據結果不斷精化模型，迭代求解研究區的形變信息。IPTA方法在地表沉降、基礎設施形變等方面應用較為廣泛[3-5]，但在滑坡形變監測方面應用較少。

本文利用C波段Sentinel-1A影像數據，采用IPTA方法對大光包滑坡進行形變監測，提取該區域2016～2018年的形變信息，并結合全球降水測量(global precipitation measurement, GPM)計劃日均降水量數據，對特征點進行時序形變分析，將該方法獲得的結果與StaMPS-SBAS方法的處理結果進行對比，分析IPTA方法應用于大型滑坡形變監測的可行性與潛力，為大光包滑坡的監測與預警提供技術和信息支撐。

1 IPTA方法基本原理

IPTA方法的核心思想是對相干點目標進行處理，通過構建和更新二維線性相位回歸模型，迭代求解相鄰相干點的線性形變速率增量和高程改正值，直至逼近線性形變速率和高程誤差的最優解[1-2]。在IPTA處理中，根據影像幅度和相位信息，計算每幅影像的光譜相干系數和時序振幅穩定性，綜合選取高質量的相干點目標[6]。光譜相干系數可表征光譜多樣性的大小，聚焦良好的相干點目標具有與分布式目標不同的光譜特性，相似的相干點目標具有相似的光譜特性，表現為低光譜多樣性，光譜相干系數較高。振幅穩定性可表征觀測目標對雷達波后向散射的穩定程度，相干點目標具有較高的散射穩定性。本研究中時序振幅穩定性閾值設為1.4，光譜相干系數閾值設為0.5，分別獲得2個候選相干點目標點集，然后對二者進行合并，剔除同名點，共獲取32 758個初始相干點。

相干點的差分干涉相位主要包含形變、殘余地形、軌道誤差、大氣延遲和噪聲等相位。IPTA以二維線性相位回歸模型為基礎，即在空間域對相干點與局部參考點的差分干涉相位進行二次差分，以削弱大氣延遲和非線性形變的影響；在時間域對差分干涉相位進行回歸處理來分析相位變化過程[7]。相干點目標的二次差分相位可表示為：

(1)

式中，φdiff,i-j為相干點目標間的差分相位差，Δv、Δh分別為建立連接關系的點目標的線性形變速率增量和高程改正值增量，BT、B⊥分別為差分干涉圖的時間基線和空間垂直基線，λ為雷達中心波長，R為雷達至地面目標的斜距，θ為雷達入射角，φres,i-j為殘余相位，由軌道誤差相位φorb_res,i-j、大氣延遲相位φatm,i-j和噪聲相位φnoi,i-j組成。

在IPTA中引入殘余相位標準差來判斷回歸分析的最優解，并采用逐步回歸、迭代改進的方式不斷計算修正值，更新線性形變相位、高程相位和殘余相位，計算所有相干點目標間的線性形變速率增量和高程誤差改正值增量，選定一個相對穩定點作為參考點，從而獲取每個相干點相對于參考點的線性形變速率和高程改正值。根據IPTA官方說明書，殘余相位標準差小于1.2的相干點可滿足精度要求。由于本文研究區地形起伏較大，且無先驗信息，因此將第一次回歸分析的高程改正值閾值設為±2 000 m，形變速率閾值設為±150 mm/a。較大的搜索范圍會解算出較大的形變速率和地形誤差，后續逐步改進閾值，縮小范圍。在最終的回歸分析中，高程改正值閾值設為±500 m，形變速率閾值設為±80 mm/a。最終輸出的相干點目標相位標準差均低于1.2，符合實驗精度要求。

從差分干涉相位中扣除線性形變和高程誤差相位可得到殘余相位。其中，殘余軌道誤差在InSAR中表現為基線誤差，可通過基線精化改正[8]；噪聲在空間維為高頻信號，在空間域對殘余相位進行低通濾波可去除；大氣延遲為時間維的高頻信息，非線性形變為時間維的低通成分，在時間域進行低通濾波可得到非線性形變。最后將線性形變與非線性形變進行組合，計算得到每個相干點的形變時間序列。IPTA數據處理流程如圖1所示。

圖1 數據處理流程Fig.1 Flow chart of data processing

2 研究區概況與實驗數據

2.1 研究區概況

大光包滑坡為2008年汶川地震觸發的最大規模滑坡，汶川地震后，該地區海拔由3 047 m降至2 964 m，其長達1.8 km的剪滑光面和高達690 m的滑坡堰塞壩引起國內外學者的廣泛關注[9-11]。大光包滑坡位于四川省綿陽市安州區高川鄉，中心坐標為104°07′0.2″E、31°38′19.8″N，處于汶川地震發震斷層上盤，距地震中心85 km，與地震地表破裂帶相距4.5 km[9]，其地理位置如圖2所示。

圖2 大光包滑坡地理位置Fig.2 Geographical location of Daguangbao landslide

大光包滑坡地質結構復雜，根據滑坡要素、滑坡堆積體形態、長度和延伸方向、巖性分布，可將滑坡體分為12個區域(圖3)。Ⅰ為側向前沖式碎屑流區，表層物質易隨白果林溝沖刷流動形成堰塞湖；Ⅱ為飛躍前沖式碎屑流區，表層物質向川林溝流動形成堰塞湖；Ⅲ為左側滑坡擾動體，巖體破碎嚴重，呈碎塊狀；Ⅳ為主滑坡體堆積區，由大光包頂拋射降落于該處，主要為泥盆系沙窩組白云巖；Ⅴ為溝道碎屑流區，多為黑色灰巖碎塊，沿黃洞子溝向下流動堆積；Ⅵ為滑坡后緣二次滑塌體，在滑動時與主滑坡體分開后停留于此；Ⅶ為坡面碎屑流區，巖性主要為灰黑色灰巖；Ⅷ為雨水沖刷滑坡表層物質形成的溝口沖積扇；Ⅸ為白云巖和板巖構成的滑坡左側斷壁；Ⅹ為滑坡后緣斷壁倒石錐，巖體松散破碎；Ⅺ為滑坡右側主滑動面，滑床為震旦系燈影組；Ⅻ為滑坡后緣斷壁，出露大光包粉砂巖至白云巖的地層剖面[10-11]。

圖3 大光包滑坡分區Fig.3 Zoning of Daguangbao landslide

研究區地勢西高東低，地形起伏大，地貌以高山為主，具有逆沖-推覆-滑脫-走滑的構造特點[9]。由于研究區常年降雨，地質環境脆弱，該區在接下來數十年內將會處于不穩定狀態，再次發生山體滑坡的風險高于地震前[12-13]，對其進行形變監測具有重要意義。

2.2 實驗數據

本文采用覆蓋實驗區的寬幅模式(IW)升軌Sentinel-1A影像，幅寬為250 km，空間分辨率為5 m×20 m，極化方式為垂直同向(VV)極化，入射角為33.83°。研究時間跨度為2016-01～2018-12，影像獲取時間間隔為12 d，共計65幅。干涉組合采用單一主影像模式，空間基線閾值設定為100 m，圖4為65幅影像的成像時間與空間基線分布。本文選用的DEM數據為2000年美國宇航局主導測量的SRTM-3 DEM數據，分辨率為90 m，用于地理編碼和差分處理以去除地形相位。

圖4 時空基線分布Fig.4 Distribution of spatiotemporal baselines

3 監測結果與分析

圖5為采用IPTA方法得到的大光包滑坡地區雷達視線向年形變速率。從圖中可以看出，監測期間大光包滑坡處于形變狀態，滑坡體左部較為穩定，形變速率為-15～2 mm/a。飛躍前沖式碎屑流區(Ⅱ區)和主滑坡體堆積區(Ⅳ區)呈現明顯漏斗狀形變，表面物質松散不穩定，受川林溝、黃洞子溝和雨水沖刷影響，形變量大于其他區域，雷達視線向最大形變速率達-50.590 mm/a，最大累積形變量達-170.726 mm。

圖5 雷達視線向形變速率Fig.5 Line-of-sight deformation rate

以72 d為時間間隔輸出2016-01-13～2018-12-28滑坡的形變時間序列，結果如圖6所示。由圖可知，監測時段內大光包滑坡右側形變漏斗逐漸發育，監測時間到684 d時，形變漏斗已初具規模，雷達視線向累積形變量為-111.543 mm。隨著時間推移，形變漏斗的范圍逐漸擴大，到2018年年底，形變漏斗長度達767 m，寬度達687 m，雷達視線向累積形變量為-170.726 mm。

圖6 大光包滑坡形變時間序列Fig.6 Deformation time series of Daguangbao landslide

在大光包滑坡形變漏斗中心和邊緣區域選擇6個特征點進行時間序列形變分析，特征點位置如圖7所示，形變結果如圖8所示。可以看出，6個特征點均呈線性滑動趨勢，但形變速率不同。A、B、C點位于形變量較大的漏斗中心，所在位置的坡度分別為49.73°、53.31°和51.49°，整體形變量和形變速率逐年增大，形變速率約為-48.387～-47.361 mm/a，累積形變量約為-141～-138 mm；D、E、F點位于形變量較小的漏斗邊緣區域，所在位置的坡度分別為33.67°、25.80°和31.06°，形變發育較為穩定，形變速率約為-17.788～-9.219 mm/a，累積形變量約為-57.432～-31.117 mm。

圖7 特征點位置Fig.7 Location of feature points

本文引入大光包地區2016-01-01～2018-12-28的日均降水量數據，結合特征點時序形變進行比較，并統計雷達重訪周期內暴雨發生前后特征點的形變量，結果如表1所示。從圖8可以看出，大光包地區的降雨量具有明顯的季節性特征，每年5～10月雨水充足，8月常有強降雨發生，滑坡上松散物質受雨水沖刷影響，碎石塊在滑坡表面滑動[14]，引起累積形變量增大。2016-08-22降雨量為118.815 mm，2016-08-16～09-09特征點的最大形變增量為-10.829 mm；2017-08-22降雨量為80.400 mm，2017-08-11～23特征點的最大形變增量為-5.877 mm；2018-06-25降雨量為71.849 mm，2018-06-19～07-01特征點的最大形變增量為4.858 mm；2018-07-10降雨量為 60.782 mm，2018-07-01～25特征點的最大形變增量為2.162 mm。2018年強降雨發生后，滑坡上松散物質受雨水影響，碎石塊在滑坡表面形成新堆積體并伴有滑坡體吸水膨脹抬升[14]，引起特征點形變量減小。綜上表明，降水是影響滑坡穩定性的主要因素之一。

表1 特征點形變量統計

為進一步分析形變漏斗特征，以D點為起點，途經A點，F點為終點，作剖面線進行研究。從上到下提取該剖面線共74個點的累積形變量，繪制以144 d為間隔的剖面線形變量時序圖(圖9)。

圖9 剖面線形變量時序Fig.9 Time series deformation of the profile

研究區地形起伏較大，最大高差達245 m，在監測期間內該坡面線累積形變量逐年增大，滑坡頂部形變量較小，形變集中在滑坡中部碎屑流區。該區域坡向為東南向，坡度為45°～50°，松散堆積體易受到降雨及黃洞子溝沖刷影響而產生較大形變，監測時段內最大形變量達-167 mm。

采用StaMPS-SBAS方法對Sentinel-1A影像數據進行處理，時間基線設為50 d，空間基線設為100m，生成212個干涉對，剔除相干性較差的

干涉對后余下183個干涉對，共提取57 225個相干點，所提取的形變速率如圖10所示。對比分析IPTA方法和StaMPS-SBAS方法所得結果可知，兩者趨勢一致，飛躍前沖式碎屑流區(Ⅱ區)和主滑坡體堆積區(Ⅳ區)均呈現出形變漏斗，且形變數值一致，表明兩者結果可靠。

圖10 StaMPS-SBAS方法得到的雷達視線向形變速率Fig.10 Line-of-sight deformation rate obtained by StaMPS-SBAS method

4 結 語

本文利用IPTA方法對升軌Sentinel-1A影像數據進行處理，獲取大光包滑坡2016～2018年的形變時間序列，并結合GPM日均降水量數據進行特征點累積形變分析，得到以下結論：

1) 大光包地區在2016～2018年發生明顯形變，形變主要集中在飛躍前沖式碎屑流區和主滑坡體堆積區，呈線性滑動趨勢，雷達視線向形變速率在-50.590～19.175 mm/a之間，最大累積形變量達-170.726 mm。

2) 結合GPM日均降水量分析發現，大光包滑坡形變與降雨量變化具有一定相關性。研究區每年5～10月雨水充足，該時間段應結合降水數據，加強滑坡監測，作好相應防范措施，防止滑坡泥石流等災害發生。

3)IPTA方法和StaMPS-SBAS方法均探測出飛躍前沖式碎屑流區和主滑坡體堆積區出現形變漏斗，形變量級一致，且形變整體趨勢相同，證明了2種方法的有效性。

4) 作為PS-InSAR方法的發展，IPTA方法不僅可以減少時空失相關和大氣延遲的影響，還可以通過矢量化數據節省計算空間，提高數據處理速度和效率，迭代求解得到的結果也更加可靠，在大型滑坡形變監測中具有較好的可靠性和廣闊的應用前景。