基于相干性的InSAR時間序列方法追溯四川雅安地區漢源滑坡災前形變

徐曉雪,季靈運,3*,張文婷,劉傳金

(1. 中國地震局第二監測中心，陜西 西安 710054; 2. 國家遙感中心 地質災害部，北京 100036;3. 防災科技學院 地球科學學院，河北 三河 065201)

0 引 言

合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar，InSAR)技術由于全天時、全天候、覆蓋范圍廣、監測精度高和空間分辨率高的優點，在山區監測活躍滑坡方面顯示出巨大的潛力。InSAR技術可以提供滑坡位置、邊界和運動特征的關鍵信息，為滑坡的災前時空形變追溯提供了新的技術途徑。自Achache等首次使用差分合成孔徑雷達干涉測量(Different InSAR,D-InSAR)技術對法國阿爾卑斯地區La Clapiere 滑坡進行監測以來，為了提高形變測量的精度和可靠性，在隨后的研究中諸如永久散射體InSAR(PSI)方法、InSAR小基線子集(Small Baseline Subset,SBAS)方法和混合時序InSAR方法等多種時間序列InSAR方法被提出來。其中，SBAS方法具有獲取微小形變信息和長時間序列緩慢地表變形的優勢，能夠在一定程度上克服時間失相干和空間失相干的不利影響，在復雜山區滑坡形變監測方面具有較為廣闊的應用前景。很多學者采用SBAS方法進行滑坡的早期識別與形變監測，如Manunta等使用SBAS方法對意大利中部翁布里亞地區進行大范圍滑坡識別與監測；Zhao等使用SBAS方法對美國西部森林覆蓋區域的滑坡進行調查，獲得了滑坡的時間序列形變特征；Liu等采用SBAS方法對美國加利福尼亞州三只熊滑坡的滑前形變進行監測，清晰地厘定了森林山區緩慢移動滑坡的運動特征。隨著SAR衛星的相繼發射和大量SAR數據的存儲，基于雷達遙感數據的SBAS方法在未來將會發揮更加重要的作用。

滑坡是中國西南山區最主要的地質災害之一，由于地形復雜多變且滑坡破壞區域較大，往往難以防范，滑坡的高精度監測對于了解其運動學特征和減少滑坡災害具有重要意義。2020年8月21日，四川省雅安市漢源縣富泉鎮中海村6組發生山體滑坡(簡稱“漢源滑坡”)，滑坡體總方量達到500×10m，屬特大型滑坡。此次滑坡造成群眾房屋受損及部分人員傷亡，省道435交通中斷，給人民生命財產安全帶來威脅。追溯滑坡體災前形變狀態，對后續開展針對性的工程地質調查及滑坡誘因分析具有重要作用。因此，本文基于升軌、降軌哨兵一號(Sentinel-1)衛星影像，求解滑坡體的InSAR時間序列，研究四川雅安地區漢源滑坡的災前演變情況。由于漢源滑坡發生在大渡河流域高山河谷地區，坡體表面作物覆蓋，極大地影響了InSAR數據處理中干涉圖的相干性，所以采用基于相干性的小基線集時間序列算法，以適應低相干分析，更好地獲得形變場，從而厘定其時空變化規律，為災后穩定性評估和未來滑坡災害的監測預警提供思路。

1 研究區概況

漢源縣位于四川省西南部，隸屬雅安市管轄(圖1中矩形框代表Sentinel-1影像覆蓋范圍)，為四川盆地與青藏高原之間的過渡地帶，省道306經縣境西北順著大渡河穿過，與成昆線相接。從地形來看，漢源縣位于大渡河流域高山河谷地區，河道下切，山高坡陡，是滑坡易發區。此次滑坡所在的區域屬典型中低山斜坡地貌，為臺階狀地形，斜坡前緣較陡。坡面上既有陡坡又有緩坡，滑坡體整體表面坡度約25°，最陡地方的坡度達40°，為利于滑坡失穩的地形地貌條件。滑坡體整體呈長條形舌狀，寬約280 m，斜長約880 m。

紅色五角星及紅色圓圈為滑坡位置；藍色方框為SAR數據覆蓋范圍圖1 四川雅安地區漢源滑坡位置和SAR數據覆蓋Fig.1 Location of Hanyuan Landslide in Ya’an Area of Sichuan and Coverage of the SAR Datasets

漢源縣區域巖土體可分為松散巖、碎屑巖、碳酸鹽巖、變質巖和巖漿巖5種巖類，易滑巖組的存在在降雨誘發下易于發生區域滑坡事件。漢源滑坡發生后，多部門開展了現場成因調查。從地質結構來看，滑坡體表層的土壤為耕作土和砂性土，結構十分松散；土層下面的巖石是風化很強烈的砂泥巖，容易失穩變形。漢源縣由于所處的獨特地理位置，具有年降雨豐沛、降雨集中和連續多日強降雨的特點。從2020年8月10日起，漢源縣經歷了兩輪強降雨天氣，巖土體含水率接近飽和，斜坡整體穩定性降低，容易發生地質災害。降雨通過表層土壤下滲到巖石，慢慢浸潤巖土體，超過穩定臨界值之后就發生了滑坡。強降雨可能是此次滑坡發生的誘發因素。

2 數據來源與分析方法

2.1 研究區SAR數據與干涉處理

本文所采用的SAR數據為Sentinel-1衛星影像。Sentinel-1衛星影像不僅可以實時免費下載，而且覆蓋范圍廣，重訪周期短，軌道精度高，可以進行全天候、全天時的高分辨率監測。由于SAR采用側視成像方式，而且西南山區地形復雜多變，致使覆蓋漢源滑坡區域的降軌Sentinel-1 SAR影像在滑坡區域形成陰影，雷達接收不到后向散射信息，嚴重影響滑坡監測的準確度。因此，本文只針對覆蓋漢源滑坡區域的升軌Sentinal-1 SAR影像進行討論，影像覆蓋范圍如圖1中矩形框所示。本研究所采用的數據分別為59景、56景兩個升軌的Sentinel-1影像數據，不同軌道SAR數據的時間跨度統一為滑坡發生之前的2018年10月至2020年8月，時間間隔為12 d。使用美國航空航天局(NASA)發布的30 m空間分辨率SRTM1 DEM地形數據模擬并消除地形對干涉相位的貢獻，與SRTM3數據相比具有更高的空間分辨率。

數據處理使用開源InSAR處理軟件GMTSAR進行。由于此次滑坡的空間范圍較小，在雷達坐標距離向和方位向不進行降采樣，多視比設置為1∶1；經過地理編碼后，單個像素的空間分辨率約為25 m。采用較小的多視比得到的差分干涉圖相干性較低，存在較多斑點噪聲，因此，對原始差分干涉圖進行濾波，以減少噪聲，提高相干性。對于軌道誤差改正，通過建立二次多項式模型來消除。由于滑坡位于復雜地貌山區，茂密的植被使SAR影像難以保持較好的相干性，也很難在植被茂密的滑坡上提供有用信號，干涉對在很大程度上受到時間失相干的影響，特別是臨滑前夏季影像生成的干涉對失相干現象明顯，導致了長時間跨度干涉圖的時間失相干。在去除低相干干涉對后，最終選擇P128軌道的54個干涉對和P26軌道的41個干涉對進行后續InSAR時間序列分析。這些干涉對的基線通常小于100 m，SAR數據的具體獲取日期與時空基線關系如圖2所示。

20181216是指2018年12月16日，其他依此類推圖2 SAR數據時空基線圖Fig.2 Spatial and Temporal Baselines of SAR Datasets

2.2 InSAR時間序列分析

SBAS方法采用較小的時空基線，在提取滑坡沿雷達視線方向(Line of Sight,LOS)的時間序列形變方面較為常用。SBAS方法的主要思路是基于一定的時空基線閾值，選取多個主影像組成自由組合的干涉對，若干涉構網形成一個唯一子集并充分連接，可采用最小二乘法求解時間序列信息。當存在多個子集時，依據各相干像元相位與觀測時間的關系，利用奇異值分解(Singular Value Decomposition,SVD)方法將多個小基線集合數據聯合起來求解，有效地解決了各數據集之間的時間不連續問題，從而求得時間序列形變、平均形變速率結果。然而，由于研究區植被覆蓋茂密，傳統SBAS方法容易受到失相干的影響，低相干像元通常使用相干性閾值進行掩膜，降低了時間序列結果的點位密度。

本文采用基于相干性的SBAS方法進行InSAR時間序列分析，將干涉圖的相位相干性引入到SBAS時間序列處理算法中，不丟棄干涉圖中存在的部分噪聲數據，盡可能多地保留干涉圖中的像素，通過相干性對觀測相位進行加權以減弱失相干的影響，提高了時間序列產品的空間分辨率。在時間序列處理過程中，通過SNAPHU軟件進行相位解纏，選擇較低的相干性閾值0.16，以獲得盡可能多的相位信息，這可能會帶來相位解纏誤差，但是在時間序列分析中失相干像元點會被降權處理。這種基于相干性的SBAS方法使用加權矩陣根據每個差分干涉圖中每個像素的相干性對觀測相位進行加權，對差分干涉圖中的噪聲不太敏感，降低了失相干相位的權重，從而可以得到相干信號的密集空間覆蓋。每個像素的加權最小二乘反演可以表示為

(1)

=diag{,,,,}

(2)

2.3 坡向形變提取

采用InSAR技術只能得到三維地表變形沿雷達視線方向的投影，然而大多數滑坡運動通常發生在平行斜坡方向。在一定的假設條件下，雷達視線方向變形可以轉移到由特定的坡度角和方位角所確定的斜坡方向。本文將雷達視線方向測量結果投影到滑坡的斜坡方向，獲得滑坡沿坡向的運動情況，平行斜坡運動可以更好地識別滑動體，提高對滑坡實際運動的量化分析。

經過投影轉換后，坡向矢量大小總是不小于雷達視線方向。從雷達視線方向到坡向的轉換可以很容易地通過放大比例因子來實現。放大比例因子()表達式為

(3)

式中：=[sinsin-sincoscos]，為雷達視線方向的形變單位向量；=[-coscos-cos·sinsin]，為坡向的單位向量；和分別表示雷達飛行方位角和入射角；和分別為坡角和坡向。

3 災前形變結果分析

3.1 雷達視線方向形變速率與時間序列

采用基于相干性的InSAR時間序列方法成功監測到四川雅安地區漢源滑坡的災前運動特征，各獨立軌道的雷達視線方向平均形變速率結果如圖3所示。其中，正值表示滑坡靠近衛星運動，負值表示滑坡遠離衛星運動。升軌P128和P26軌道在2018～2020年變形分布相似，但變形幅度略有不同(圖3)。這主要是因為滑坡位于P128軌道的近距離范圍，而位于P26軌道的遠距離范圍，相對于同一目標而言，其雷達視線方向存在9.9°的偏差。兩個升軌的平均形變速率結果顯示，滑坡存在明顯的滑前形變信號。在滑坡發生前的近兩年時間，P128軌道的最大平均形變速率可達到-42 mm·年，P26軌道達到-42.1 mm·年，非形變區大多數監測點的平均形變速率位于-10～10 mm·年，兩個升軌的形變值為負值表明滑坡遠離衛星運動。從兩個軌道的平均形變速率結果可以看出，位移變化較大處主要發生在河流沿岸，這與現場識別的活動滑坡位置一致。

為了更好地分析此次滑坡的滑前時空運動特征，對2018～2020年兩個升軌的形變時間序列進行研究。圖4為P128和P26兩個軌道在雷達視線方向上的時間序列變形。從圖4可以看出，不同軌道滑坡運動的時空變化規律基本一致，滑坡在整個監測期間存在持續變形。總的來說，滑坡在2018年底至2019年4月移動迅速，而在之后移動緩慢，基本處于穩定狀態。滑坡運動在空間上并不均勻，隨著時間的積累，滑坡形變區域逐漸形成形變中心，滑坡中心的累積形變可達到97 mm，而位于周邊區域的累積形變相對較小，最大累積形變僅為70 mm。

圖3 雷達視線方向平均形變速率結果Fig.3 Images of Average LOS Deformation Rate

兩個軌道的時間序列形變結果表明，漢源滑坡的災前形變過程經歷了形變快速積累和穩定兩個階段。由于缺少可用的夏季干涉對，時間序列變形分析未獲取2019年7月至9月及臨滑前的時間序列信息，災前形變分析并沒有捕捉到臨滑前的形變加速信號。

20181221是指2018年12月21日，其他依此類推圖4 雷達視線方向累積形變結果Fig.4 Images of LOS Cumulation Deformation

圖5 雷達視線方向單點累積形變統計結果Fig.5 Statistical Results of LOS Cumulation Deformation for Single Point

本文提取了漢源滑坡形變中心區域P點[圖3(a)]的雷達視線方向累積形變時間序列，結果如圖5所示。單點形變時間序列結果顯示，兩個升軌表現出明顯的滑前形變信號。滑坡運動經歷了兩個時段的形變：Ⅰ階段為2018年10月至2019年4月，形變加速積累，其中P128軌道累積形變量達97 mm，P26軌道累積形變量達86 mm；Ⅱ階段為2019年4月以后，滑坡運動速率幾乎為0，并在2020年8月21日發生滑坡。由于失相干現象，缺少夏季的干涉對，在2019年7月至2019年9月和臨滑前(即2020年6月至8月)并沒有獲得有效形變信號。2019年4月以后滑坡運動緩慢的具體原因還需結合后續工程地質調查展開。

3.2 坡向平均形變提取與結果驗證

為了分析漢源滑坡滑前沿坡向的運動特征，首先通過研究區DEM數據提取滑坡的坡度和坡向。雖然不同軌道數據集測得的雷達視線方向的形變量級不同，但是平行坡向的分量是相同的。因此，在提取坡度和坡向之后，結合式(3)分別計算了P128軌道和P26軌道的放大比例因子。為了控制坡向的偏差，剔除了絕對放大值大于一定閾值(本文取值為2)的孤立點，便于更好地分析滑坡沿坡向的運動規律。

圖6展示了P128軌道和P26軌道經放大比例因子逐像素放大校正后得到的平行坡向的平均形變速率結果。從圖6可以看出，滑坡發生前在斜坡方向上存在明顯的災前形變信號，滑坡方向的最大平均形變速率可以達到43.9 mm·年，顯著大于雷達視線方向。滑坡體中心為滑前形變最劇烈的區域，整體呈現出明顯的漏斗狀特征。滑坡沿坡向的形變結果顯示，形變中心區域長約250 m，寬約250 m，形變中心海拔高度1 017 m，該區域是在滑坡發生前發生顯著形變的范圍。

由于沒有獲取研究區地面監測數據，本文使用相同時間跨度的獨立SAR數據集來驗證InSAR監測結果的精度。通過統計2個軌道獨立觀測的SAR數據結果重疊點的差值，驗證InSAR結果之間的一致性。統計結果顯示，2組相鄰軌道滑坡區域重疊點的坡向速率差值符合正態分布，速率差值的標準差為2.9 mm·年，表明InSAR監測結果是可靠的。

綜合此次滑坡發生前InSAR時間序列結果，漢源滑坡存在明顯的災前形變積累。在滑坡發生地為高山河谷的地形和砂泥巖地質結構背景下，經過強降雨的加載作用，可能打破邊坡的穩定階段。足夠的降水浸透滑坡體基底面，降低有效抗剪強度；當邊坡剪切應力超過剪切強度時，就可能觸發滑坡。結合滑坡發生前該區域經歷的強降雨過程，推測在災前形變積累情況下，強降雨可能是此次滑坡發生的誘發因素。

4 結 語

本文利用多軌道Sentinel-1衛星影像，采用基于相干性的InSAR時間序列方法，對四川雅安地區漢源滑坡的災前形變情況(2018年10月至2020年8月)進行追溯。首先，獲取了滑坡沿雷達視線方向的平均速率與時間序列，結果顯示滑坡雷達視線方向的最大平均形變速率為42.1 mm·年；漢源滑坡從2018年底便發生了變形，至滑坡發生前最大累積形變量達97 mm，共經歷了形變加速和穩定兩個階段，形變在空間分布不均勻。然后，提取了滑坡沿坡向的形變，結果顯示在坡向上同樣存在明顯形變信號，最大坡向平均形變速率達43.9 mm·年。本文研究結果為此次漢源滑坡提供了第一手、相對完整的形變測量資料，所揭示的滑坡時空運動特征對于評估滑坡危害及進一步實地調查至關重要。上述基于InSAR技術的災前形變追溯方法也可以用來識別其他滑坡可能發生點的前兆運動。

圖6 漢源滑坡平行坡向的平均形變速率結果Fig.6 Images of Average Deformation Rate of Hanyuan Landslide Parallel to Aspect

值得注意的是，由于研究使用的Sentinel-1 SAR影像數據植被穿透能力較弱，在研究區植被茂密的情況下相干性會受影響。本文采用的基于相干性的InSAR時間序列方法雖然在一定程度上可以提高形變結果的點位密度，但由于夏季植被非常茂密，臨滑前干涉對相干性很差，所以本文形變追溯結果并未捕捉到臨滑前明顯的加速變形信號。這也說明InSAR技術在滑坡監測預警中容易受到觀測環境影響，可能存在時間采樣率不足等缺點。

歐洲航天局提供了Sentinel-1衛星影像，在此表示感謝!