Stacking InSAR和IPTA技術應用于探測寧夏隆德縣滑坡隱患

王 川，涂 寬，諶 華，耿 丹，王文龍，李樵民

（1.二十一世紀空間技術應用股份有限公司，北京 100096；2.寧夏回族自治區遙感調查院，寧夏銀川 750021；3.高分辨率對地觀測系統寧夏數據與應用中心，寧夏銀川 750021）

引言

滑坡是在重力、降雨、地震等因素的作用下，坡體作剪切運動的自然現象，通常可劃分為推移式滑坡、土流型滑坡、土質滑坡和古崩滑坡堆積體等［1-2］。其產生條件包含內因和外因兩類，內因為地質構造和地貌特征，外因則為降雨、地震和人類工程活動等。近年來由于全球氣候波動，異常天氣條件增多，導致更多滑坡災害在各地發生［2-4］。

由于滑坡災害的發生嚴重威脅了人類的生命和財產安全，對人類社會環境、經濟環境、生態環境造成了巨大破壞，因此需要對滑坡隱患進行及時的探測和監測。通常，對滑坡進行探測和監測的手段有：常規三角網測量、伸縮計深部位移監測、應力應變測量、GPS遙測監測網、地震監測、近景攝影測量以及多種技術的聯合等［2，5］。上述探測和監測手段的準確性高，測量結果可靠，但測量主要集中在小區域，不利于開展大規?；码[患探測。

近年來，隨著科技的不斷進步，出現了以測量可見光光譜、熱輻射、微波輻射等為目標的多種遙感技術［6-7］。特別是合成孔徑雷達（Synthetic Aperture Radar，SAR）技術的出現。微波遙感干涉合成孔徑雷達（Interferometry Synthetic Aperture Radar，InSAR）測量地表形變技術測量結果準確度高，覆蓋面大，在滑坡隱患探測中得到了普遍應用［8］。最先發展起來的應用電磁波相位值對地表起伏進行量測的是差分干涉技術（Differential Interferometry Synthetic Aperture Radar，DInSAR），但是單獨使用兩景影像進行干涉易產生各種去相干（時間、空間、體）等，且易受大氣、系統熱噪聲、處理誤差等影響，因此目前主流技術是多個干涉對共同測 量 的 時 序InSAR技 術（Multi-Temporal Interferometry Synthetic Aperture Radar，MT-InSAR）［2，9］。在MT-InSAR中，典型的是：（1）永久散射體技術（Persistent Scatter InSAR，PS-InSAR），該方法只選擇回波信號穩定的點（永久散射體）進行分析，可以擺脫時空基線的限制，相干性保持好；（2）小基線集技術（Small Baseline Subsets InSAR，SBAS-InSAR），該技術選擇分布式散射體（如沙地、草地等）進行分析，更適用于自然環境［10-12］。以此兩種思想為核心，產生了相干點目標分析（Interferometry Point IPTA）、SqueeSAR等技術，同時也出現了以計算形變速率為目標的干涉疊加技術（Stacking）和不依賴時序相位特征測量地表形變的多孔徑雷達技術（Multiple Aperture Interferometry，MAI）［13-19］。InSAR技術不可避免的會受到植被、地形起伏和大氣等的影響，使用單一技術往往只能克服部分缺陷而導致不能更全面的識別滑坡隱患，同時基于光學影像對滑坡隱患進行專家判讀缺少坡體滑移的數據支持，部分老（古）滑坡已趨于穩定。因此需要聯合多種InSAR方法并結合光學影像進行滑坡隱患探測。

IPTA技術針對點目標進行回歸分析以分離相位（大氣、地形、非線性形變等），計算效率和精度顯著提高，而Stacking技術利用加權平均以抑制噪聲使其能夠快速獲取大范圍的形變信息，在現有各類InSAR技術中，此2種技術可操作性最強，同時結果可互相補充。隆德縣境內以黃土丘陵區為主，較少出現大型滑坡隱患，小型滑坡隱患分布廣泛，利用Stacking技術能夠實現對大部分滑坡隱患的探測，基于IPTA技術能夠去除各種誤差，探測滑坡隱患的可靠性高。

因此，本研究利用Sentinel-1衛星影像，采用IPTA和Stacking這2種InSAR技術對隆德縣全境滑坡隱患進行探測，分析和評價2種InSAR技術的探測結果。該項工作對隆德縣滑坡隱患普查及監測具有十分重要的指導意義。

1 研究區概況及研究數據

研究區隆德縣位于寧夏回族自治區的南部，六盤山西側，全縣總面積992.39 km2。在地形地貌特征上，隆德縣處于昆侖秦嶺地槽褶皺區東段，曾受到燕山期、喜馬拉雅期構造運動的影響，其主要地貌類型有：紅層巖丘陵地貌、中低山地貌、河谷平原地貌、黃土丘陵地貌。在氣候特征上，隆德縣屬于溫帶大陸性季風氣候，受東高西低的地勢差異影響，致使東部濕潤寒冷而西部干燥溫暖，冬季干燥而夏季多雨。每年7～9月降雨最為集中，約占全年降水量的60%。降雨和降雪是隆德縣地質災害發生的重要誘因。由于人類工程活動的影響，如毀林開荒、削坡修路、興修水利和過度放牧等，破壞了生態系統平衡，加劇了地質環境惡化，更進一步的促使了滑坡災害的發生。

隆德縣滑坡包含黃土滑坡和巖質滑坡2種，其中黃土滑坡包含黃土層內滑坡和黃土—基巖滑坡，巖質滑坡主要是新近系泥巖殘坡積物滑坡。滑坡規模通常為小型滑坡，中、大型滑坡發育較少?；潞突码[患多分布于渝河、南川河的沖蝕岸邊及其支流或支溝兩側，呈條帶狀或樹枝狀分布，在地貌類型上多集中在丘陵地貌區［20］?；乱詼\層和中層為主，厚層滑坡較少。

Sentinel-1衛星影像幅寬廣（IW成像模式可達250 km），重訪周期相對較短且固定（單星重訪周期12 d），因此文中研究選擇2019年1月至2020年5月在重訪周期過境的全部衛星影像用于測量地表形變，表1詳細列出了所使用的衛星影像參數信息。圖1為研究區概況圖和Sentinel-1衛星影像覆蓋范圍。此外，為了糾正地形影響和提取地形信息，文中研究還同時獲取了SRTM 30 m的DEM作為輔助數據。

表1 Sentinel-1衛星影像參數表Table 1 Information of Sentinel-1

注：在Sentinel-1假彩色合成圖中，綠色區域反映出地物主要由體散射比較強烈的散射體如植被構成，深粉色區域反映出地物主要由粗糙的表面如裸地構成，淺粉色區域表示地物由去極化能力非常弱且后向散射強烈的建筑物組成，水體由于鏡面散射而呈現黑色。

本次實驗所獲取的SAR數據能夠完整覆蓋2019年全年，時間和空間基線滿足形變測量需求，配準精度達到InSAR干涉測量要求（1/8像元精度），DEM數據空間分辨率與SAR數據接近且已對空值進行插值，能夠利用全部衛星影像和DEM數據開展后續研究工作。

2 基于Stacking和IPTA技術的滑坡體探測方法

2.1 InSAR成像原理

根據圖2幾何關系，P點的基線分解為：

同樣，參考平面的點P0的基線分解為：

假設2次觀測時沒有地表起伏變化，同時忽略地物散射和輻射傳輸中產生的隨機噪聲，則P點的干涉相位可以視作2次傳輸路徑之差所對應的相位：

對路徑差ΔR求導，可知相位隨路徑的變化關系為：

對于星載SAR衛星，當2次觀測視角的差異較小時，基于微波輻射傳輸的遠場平行近似［21］，可以將ΔR等效替換為水平基線B∥，當觀測距離不變且視角由θ0開始變化產生的干涉相位變化為：

若定義參考面點P0的干涉相位為參考相位，其與P點到傳感器S1的距離相同，因此可以視其為P點的平地相位（排除地形起伏的影響，僅由觀測視角θ0造成）：

根據圖像的幾何關系：

則視角由θ0開始變化產生的高程變化：

結合式（7）和式（10）可得：

圖2 InSAR原理圖Fig.2 InSAR schematic diagram

此處??if的變化僅與高程相關（P0為參考面起始觀察點），即由P和P0點的干涉相位差分得到，所以稱其為高程相位?h（忽略大氣、形變等噪聲相位）。

若P點在觀測期間產生了位移（地表形變），則干涉相位可表示為：

式中：?d表示地表形變產生的相位，與之對應的雷達視線向的形變量為ΔRd，由地表形變相位與視線向形變量的關系可以看出，對于C波段SAR衛星，其地表形變的測量精度通常能達到毫米級。

2.2 研究方法

本研究選擇Stacking和IPTA這2種較為成熟的InSAR方法用于測量隆德縣境內的地表形變值。Stacking技術采用對干涉相位加權平均的方式計算相位的形變速率繼而得到地表延視線向的形變速率，其處理方法相對簡單且能夠在一定程度上對大氣相位等噪聲進行抑制，因此能夠快速的對大區域在長時間序列上得到地表形變速率的測量結果［22］。IPTA方法是一種基于點目標的干涉測量方法，通過迭代回歸分析將大氣相位、高程相位殘余等進行分離，能夠同時獲取地表形變速率和地表非線性形變，與SBAS和PS等方法相比具有計算效率高、節省存儲空間的特點［23］。

在進行InSAR處理前，需要首先對Sentinel-1衛星影像進行預處理，主要包括：讀取原始數據并基于外部DEM數據對SLC圖像進行配準，將配準后的圖像裁剪至研究區范圍，對DEM進行重投影至SAR坐標系。兩種InSAR方法的主要處理流程如下：

2.2.1 Stacking InSAR技術流程

（1）計算各個像對間的垂直基線，設置閾值為-200 m～200 m，計算時間基線，設置閾值50 d，生成基線連接圖和干涉對列表（圖3）。

（2）根據干涉對列表對SLC圖像進行差分干涉，對差分干涉圖進行濾波以提高相干性，對干涉圖進行解纏，對解纏圖進行一定的空間濾波部分消除大氣、軌道誤差、高程誤差等噪聲相位。解纏方法：最小費用流法（MCF）［24］。

（3）基于加權平均的方法計算相位形變速率：

式中：φrate為相位形變速率；Δtj表示第j個干涉相對的時間間隔；N表示干涉像對數，然后按照λ/4π將相位形變速率換算為延視線向的地表形變速率（單位：mm/年）。由于單軌干涉測量存在著多種去相干（如時間、幾何、體等），在解纏過程中需要對相干性較低的像元進行掩膜（文中研究相干性閾值為0.2），使得在進行Stacking時，各像素點所使用的樣本數（干涉像對數）均不同，為了使結果更具可靠性，文中研究設定最小可接受樣本數為108，即當像素點干涉像對數低于108則對結果進行掩膜。圖4為Stacking InSAR總體處理技術流程圖。

圖3 Sentinel-1衛星影像基線連接圖Fig.3 Sentinel-1 perpendicular base connection diagram

圖4 Stacking InSAR處理技術流程圖Fig.4 The flow chart of Stacking InSAR

2.2.2 IPTA技術流程

（1）分別基于像素點在時序上的穩定性和相干性挑選像元并生成點目標文件。計算兩兩像對間的垂直基線，設置閾值為-200～200 m，計算時間基線，設置閾值50 d，生成基線連接圖和干涉對列表（圖3）。根據干涉對列表對點目標進行差分干涉。

（2）在第j幅差分干涉圖，點目標i其相位可以表示為：

圖5 IPTA處理技術流程圖Fig.5 The flow chart of IPTA

式中：Δ?表示差分干涉相位；vi表示第i個點目標的地表形變速率；Δtj表示第j個干涉對的時間間隔；δi,j表示第i個點目標在Δtj時間段內的非線性地表形變值；B⊥j表示第j個干涉像對的垂直基線；Δzi表示第i個點目標的高程誤差；λ表示波長；ρi表示微波從雷達到達第i個點所對應的地物再返回的路徑（斜距）；θi表示入射角；?atm表示大氣相位；n表示噪聲相位。通過非線性回歸分析分別求解出其中的差分大氣相位和高程誤差，對地表高程進行補償并對差分大氣相位進行去除后再次進行回歸分析，進行多次迭代后可將差分大氣相位逐步分離，并對地表高程進行充分補償。

（3）通過回歸分析獲得解纏的差分干涉相位后，使用奇異值分解（SVD）的方法對多主影像的相位進行求解計算出延時間序列的相位形變［25］，對其在空間上和時間上進行濾波，并換算成地表沉降值（單位：mm）。在得到時序相位形變后，再次進行回歸分析計算出相位的形變速率，同樣換算為地表形變速率（單位：mm/年）。圖5為IPTA總體處理技術流程圖。

2.2.3 有效探測滑坡隱患和未有效探測滑坡隱患的定義及有效探測率

在開展文中研究工作之前，已先期對隆德縣滑坡隱患進行了調查，此次工作聯合多源遙感技術對隆德縣全境滑坡隱患進行了普查，在此基礎上對所識別出的滑坡隱患已完成了實地的調查和驗證工作。該項工作由“寧夏科技廳重點研發計劃項目”支持開展［26］。

根據InSAR地表形變測量結果識別滑坡隱患，若圈定的滑坡隱患與隆德縣滑坡隱患遙感監測工作成果中的滑坡隱患相對應，則稱其為該InSAR技術（Stacking和IPTA）的“有效探測滑坡隱患”，對于使用該InSAR測量方法未解譯出但在隆德縣滑坡隱患遙感監測工作成果中出現的其他滑坡隱患，則稱其為“未有效探測滑坡隱患”。對于一種InSAR技術（IPTA或Stacking），其有效探測滑坡隱患占所有滑坡隱患的百分比文中研究稱其為“有效探測率”。

3 Stacking和IPTA滑坡隱患探測結果

通過與隆德縣滑坡隱患調查結果進行核實，2種InSAR方法對47處滑坡隱患進行了有效探測［26］（圖6，圖7（a）、（b）），在形變測量中正值表示滑坡隱患在延視線向靠近衛星滑動（滑坡隱患和衛星間的距離縮短），負值表示滑坡隱患在延視線向遠離衛星滑動（滑坡隱患和衛星間的距離增大）。表2詳細列出了滑坡隱患的基本信息，其中滑坡隱患面積達到2.0×106m2的共6處，最大滑坡隱患面積6.2×106m2。

Stacking技術獲得的形變結果中有效探測滑坡隱患44處，有效探測率93.6%，IPTA有效探測滑坡隱患共計37處，有效探測率78.7%，Stacking技術略高。2種InSAR技術共有效探測滑坡隱患34處，在相同數據源情況下監測重合度高，對最大形變速率進行相關性分析（圖8），線性相關性強（R20.56），測量結果一致性強（例如圖7（c）和（d）所示滑坡隱患），可相互印證2種InSAR技術測量結果可靠。

圖6 隆德縣境內滑坡隱患分布圖Fig.6 Distribution map of potential landslide in Longde County

圖7 隆德縣Stacking和IPTA測量形變結果圖及典型形變滑坡隱患Fig.7 The displacement map measured by Stacking and IPTA and the typical potential landslide in Longde County

表2 隆德縣解譯滑坡隱患基本信息表Table 2 Information of interpreting potential landslide in Longde County

續表

圖8 不同InSAR方法測量滑坡隱患最大形變量線性相關性分析Fig.8 Correlation analysis of the maximum displacement for different InSAR methods

圖9 有效探測滑坡隱患最大形變速率分布特征Fig.9 Characteristics of the maximum displacement of effectively detected potential landslides

3.1 Stacking和IPTA有效探測滑坡隱患結果分析

3.1.1 有效探測滑坡隱患的最大形變速率分布特征

圖9描述了滑坡隱患最大形變速率分布特征。箱形圖中間箱體上邊沿表示上四分位數（3/4分位數），下邊沿表示下四分位數（1/4分位數），上下四分位數之差為四分位距（IQR），箱體中間橫線表示中位數，箱體中圓圈表示均值，由箱體延伸出上下邊緣線（箱體外橫線），上邊緣線為上四分位數加1.5倍四分位距，下邊緣線為下四分位數減1.5倍四分位距，超過上下邊緣線的視為異常值。2種InSAR方法對不同形變速率滑坡隱患探測能力基本一致。其中一半的滑坡隱患最大形變速率集中在-20～20 mm/年之間（箱的上邊沿和下邊沿分別表示3/4和1/4分位數）。Stacking所識別滑坡隱患形變速率中位數（箱內橫線）高于IPTA。

3.1.2 有效探測滑坡隱患的地形分布特征

研究主要針對地形中的坡度和坡向進行分析，其中坡度劃分為平坡、緩坡、斜坡和陡坡4個等級（表3）［27］。坡向劃分為北、東北、東、東南、南、西南、西和西北8個方向（表4）［27］。IPTA和Stacking有效探測滑坡隱患均集中在緩坡和斜坡上，對于坡度較大的平坡（1.7°～5°）和坡度較小的陡坡（25°～31.6°）亦有少量分布（圖10（a）），這與InSAR的識別原理相符合，當坡度較小時滑移速率較慢無法滿足InSAR測量精度，坡度過大則滑移速率過快易產生失相干導致無法測量。北坡方向未出現有效探測滑坡隱患，與北坡所受光照輻射少其微氣候環境不易出現滑坡相關，在其他方向上未出現明顯集中分布特征（圖10（b）），因此綜合2種InSAR方法能夠較全面的對各個方向的滑坡隱患進行有效探測。

表4 地形-坡向統計結果Table 4 Statistics of topographic aspects

表3 地形-坡度統計結果Table 3 Statistics of topographic slopes

坡度級劃分來源于“基于3PGS-MTCLIM模型模擬根河林區火后植被凈初級生產力恢復及其影響因子”。

圖10 滑坡隱患地形分布特征Fig.10 Characters of topography of potential landslides

3.2 滑坡隱患滑動特征分析

圖11 IPTA有效探測滑坡隱患的時序形變量Fig.11 The displacement of potential landslides detected by IPTA method

IPTA技術不僅能夠獲取滑坡隱患的年平均形變速率，同時能夠獲取時序累計形變量，因此可以通過IPTA技術了解滑坡隱患的滑動狀態。文中研究對IPTA獲取的37處滑坡隱患的滑動過程做一基本描述（圖11）。

LDhp1、LDhp3、LDhp4、LDhp5、LDhp9、LDhp13、LDhp26、LDhp28、LDhp37、LDhp44均以一定的滑動速率發生線性形變，滑動過程在整個監測期內較穩定，沒有表現出受其他因素干擾而滑動過程受到擾動的情況。

LDhp2、LDhp6、LDhp8、LDhp11、LDhp12、LDhp18、LDhp22、LDhp30、LDhp32、LDhp43在前期滑動不明顯或者滑動速率保持穩定，在經過一段時間后滑動速度突然加劇并且會維持至監測結束，這一加速過程發生在2019年6月～9月時間段內，這一時間段是一年中降水集中時期，因此推測是由于降雨沖刷坡體導致滑動過程加劇。

LDhp7、LDhp10、LDhp14、LDhp16、LDhp24、LDhp25、LDhp27、LDhp29、LDhp39同樣在2019年6月～9月出先了較明顯的加速過程，但隨后在2019年10月～12月份加速過程突然減緩進入新的相對慢速滑動的過程，這一特征應該與該地區溫度降低導致土層的穩定性發生改變相關，其滑動過程是否隨季候變化呈現周期性還需要在更長時間周期內監測進行驗證。

LDhp15、LDhp17、LDhp19、LDhp21、LDhp23、LDhp38、LDhp42、LDhp47其最主要的特征是在2019年6月～12月相對于其他時間段始終保持緩慢滑動的狀態，而在2019年1月～5月份或是2020年1月～5月份則出現了明顯加速滑動的過程，這種在非生長季滑動過程加快應于坡體上植被的季候變化相關。

4 結論

文中研究基于Sentinel-1衛星影像進行Stacking和IPTA處理，對隆德縣境內的滑坡隱患進行探測，核實探測到的滑坡隱患共計47處，其中，Stacking技術有效探測滑坡隱患44處，IPTA有效探測滑坡隱患37處，2種InSAR技術共同識別34處?；?種InSAR技術結果：

（1）Stacking和IPTA技術所探測到的滑坡隱患最大形變速率線性相關性強，形變測量結果可相互印證，可以認為基于InSAR技術進行滑坡隱患的探測具有較高的可靠性。

（2）由于IPTA需要選擇高質量的點目標從而導致部分滑坡隱患無點目標覆蓋，而使用Stacking方法能夠對這些滑坡隱患進行有效識別。點目標的缺失與野外環境復雜有關，如植被等地表覆蓋物造成相干性降低，地形起伏在局部地區出現大量地形相關的大氣相位。后續應聯合多種InSAR技術，如采用DS InSAR方法增加同質點的選?。?8］，采用BM3D進行更好的濾波降噪，提高IPTA分析的點目標覆蓋范圍和測量精度，從而提高滑坡隱患探測的可靠性。

（3）在形變速率上，隆德縣以黃土丘陵區為主要地貌特征，較少存在劇烈的坡體滑動過程，因此Stacking和IPTA技術能夠對該地區不同形變速率的滑坡隱患進行有效探測。同樣由于IPTA的選點機制導致滑動較的坡體候選點減少，其所能識別的滑坡隱患形變速率中位數略低于Stacking技術。

（4）在地形特征上，Stacking和IPTA所探測到的滑坡隱患主要集中在緩坡和斜坡上，同時除北坡不易出現滑坡外，在各個坡向的滑坡隱患均可探測。

（5）基于IPTA獲取的時序形變量表明，隆德縣境內的滑坡隱患存在4種典型滑動類型：（a）以穩定的速率持續滑動；（b）在2019年6月～9月份出現明顯加速滑動過程并持續到監測結束；（c）在2019年6月～9月份出現明顯加速滑動過程，10月～12月份滑動過程逐漸減緩；（d）在2019年或2020年1月～5月滑動過程明顯，其它時間坡體較穩定。

總之，綜合Stacking和IPTA這2種InSAR技術能夠全面的識別隆德縣境內的滑坡隱患，在后續滑坡隱患探測中應聯合多種InSAR方法，以獲取更可靠的結果。下一步將首先基于GNSS等技術對Stacking InSAR和IPTA技術得到的平均形變速率和時序形變結果進行精度評價，并聯合其他InSAR技術提高IPTA的探測能力。