基于時序InSAR的滑坡早期識別
——以貴州省水城縣為例

廖 軍，吳彩燕，王立娟，譚秋焰，朱新婷

（1.西南科技大學環境與資源學院，四川綿陽 621010；2.四川省安全科學技術研究院，四川成都 610045）

引言

滑坡是指巖體、碎屑、土壤沿斜坡自上而下滑動的現象，也指斜坡在重力作用下沿著坡體內部一個或者多個面做剪切運動的現象［1～3］，在我國的山區和峽谷地帶廣泛分布，尤其是西南高山區和三峽庫區，發生頻率最高，造成的損失最嚴重［4～5］。常規的滑坡識別手段主要有野外實地勘查和遙感影像識別，野外實地勘查能準確圈定滑坡要素，但其工作既費時又低效，并且勘測區域受到地形的極大限制［6］；遙感影像識別滑坡大多是利用專家經驗知識通過光學遙感影像的光譜、空間、紋理等信息進行的［7～8］，盡管已經取得了相當不錯的成就，但光學影像容易受到天氣因素的限制如大霧、降雨等，并且在監測和識別滑動不明顯的潛在滑坡能力不足［9～10］。而新出現的合成孔徑雷達干涉測量（InSAR）技術不僅能精確獲取地表微小形變信息，且具有全天候、全天時、時效性高等優點，成為了國內外學者研究的熱點［11～14］。

近年來，InSAR技術逐漸應用于滑坡災害識別，并取得了一些成功的案例。廖明生等以三峽庫區滑坡為對象，使用高分辨SAR數據對其進行監測，成功獲得了滑坡發生的時間、位置及形變情況，證明了高分辨率雷達數據有效應用在滑坡監測中［15］；Tofani V等［16］采用永久散射體干涉測量（Persistent Scatterer，PS-InSAR）技術對意大利北部的Santo Stefano d′Aveto滑坡進行分類與監測，結果表明PS技術結合地面數據能夠準確判別滑坡的邊界和運動狀態；Bianchini等［17］以意大利的Calabria地區為研究區，利用PS-InSAR技術對該地區的滑坡進行了識別及滑坡編目的更新；趙超英等［18］利用多源SAR數據來識別黑方臺地區的潛在滑坡，結合實地調查和光學影像驗證了識別結果的可靠性，并分析典型滑坡的失穩模式；Dong等［19］使用相干散射體（Coherent Scatterer，CS-InSAR）技術更新了丹巴縣甲居滑坡的活動邊界，繪制了更詳細的滑坡位移圖。現有的研究雖然在一定程度上可以對潛在滑坡進行識別，但是這些方法大多是從雷達視線方向（LOS）角度出發進行識別研究，雖然可以識別滑坡發生的時間和地理位置，但LOS方向與地表和山坡之間存在一定的夾角，不能直接反映滑坡的變形程度，也不能得到滑坡變形的真實時間序列，難以識別滑坡前后的具體變形趨勢。

本研究以滑坡災害多發的水城縣為研究區，通過改進的SBAS-InSAR技術，獲得水城縣坡度方向形變信息，結合水城縣雷達可視分區及谷歌影像剔除無效觀測值，利用Anselin Local Moran′s I對獲取的監測值進行異常值分析和聚類處理，獲取低值聚集區域，最后結合谷歌影像的地表特征，實現水城縣潛在滑坡的識別，用于指導該區域滑坡災害的預測預警。

1 研究區概況及數據源

1.1 研究區概況

水城縣位于貴州省六盤水市，鄰接云南省，在水系和構造運動作用下，西北地勢較東南高。該區屬于亞熱帶濕潤季風氣候地區，降雨集中在5～10月份，其中6～8月降雨量占全年53.03%［20］，縣內河網密閉，水量集中，是造成境內地質災害的主要誘因之一。水城縣地質構造復雜，存在多個構造變形區，地層巖性以二疊系、三疊系以及第四系為主。

水城縣境內滑坡廣泛分布在飛仙關組、峨眉山組及龍潭組等地層中，該地層中節理發育，巖層破碎，極易在人類活動或降水等因素下發生滑坡等地質災害［20～21］。

圖1 研究區地理位置Fig.1 Geographical location of the study area

1.2 數據源介紹

文中采用的數據源包括歐空局新一代C波段的中高分辨率衛星Sentinel－1A SAR數據，選取了2018年1月～2018年12月共20景降軌數據，入射角為39.5°，訪問周期為12 d，拍攝模式為IW，空間分辨率為5×20m（Rg×Az），極化方式為VV，并使用POD精密定軌星歷數據（Precise Orbit Ephemerides，POD）進行軌道矯正，定位精度優于5 cm；數字高程模型（DEM）為30 m分辨率的SRTM DEM（Shuttle Radar Topography Mission，SRTM）；從Google Earth中獲取的研究區光學遙感影像數據；野外調查資料。

2 研究方法

本研究選取2018/1～2018/12共20景Sentinel－1A SAR降軌數據，利用永久散射體干涉測量（PS-InSAR）技術提取獲得相干性高且相對穩定的散射體用做地面控制點（Ground Control Point，GCP），引入到小基線集干涉測量（SBAS-In-SAR）技術處理流程中，獲取長時間序列的地表形變信息（LOS向），并轉換為坡度方向形變信息；結合水城縣地形可視性分區及谷歌影像，剔除不可視區（陰影、透視收縮、疊掩）、高密度林地及灌木林、河流等不可靠的形變信息，保留其他區域可靠形變信息；對獲取的可靠形變信息采用Anselin Local Moran′s I指數進行空間異常值分析和聚類處理，獲取低值聚集區域，并結合谷歌影像的地表特征，從而對潛在滑坡進行識別。滑坡識別路線如圖2。

圖2 潛在滑坡識別流程Fig.2 Landslide identification process

2.1 結合PS-InSAR技術的GCP點選擇

在SBAS-InSAR處理流程中，選取穩定的地面控制點（GCP）進行軌道精煉和相位偏移的計算，可以有效地估算和去除殘余相位及可能存在的斜坡相位，進而消除平地效應，對相位轉形變的可靠性起著重要作用［16－18］。

由于人工選點會對形變速率引入誤差，文中利用PS-InSAR技術的優勢［22］，依次設置振幅離差閾值（≤0.25）、相干性閾值（≥0.8）以及空間域的低通濾波和時域的高通濾波，減小干涉相位中的大氣延遲相位誤差，得到可靠的穩定散射點的形變信息，并對該結果設置速率閾值選取在－0.1～0.1 mm的穩定點作為GCP點，通過地理坐標到SAR斜距坐標系的轉換，即可應用于在SBAS處理流程中。

振幅離差指數（D）是依據SAR圖像上同一位置像元的強度值在時間序列上的離散特性進行目標的識別，其計算公式如下［23］：

式中：μ是時序上的某個點的振幅平均值；σ是時序上的振幅的標準差。

2.2 SBAS-InSAR技術原理及方法

小基線集干涉測量（SBAS-InSAR）是在差分干涉測量（Differential Interferometry Synthetic Aperture Radar，D-InSAR）和永久散射體干涉測量（PS-InSAR）發展起來的一種時序監測方法，該方法能有效消弱常規D-InSAR中大氣效應，時空失相干等因素對監測結果的影響，同時彌補了PS-InSAR技術不適用于相干點稀疏地區的不足［10，16～17］。在2002年，Berardino等［24］首次提出了小基線集時序監測方法，SBAS-InSAR能夠獲取非線性的形變結果，在滑坡形變監測中具有良好的適應性。

小基線集的原理是通過設置時間、空間基線的閾值，將多期SAR影像進行自由配對組合生成多個子集，每個子集由超主圖像連接，從而形成一個大集合，保證在每個小集合中SAR影像間的基線小，集合間的SAR影像基線大，從而提高數據的采樣率，克服空間基線過長導致的空間失相干影響。配合外部高精度的DEM消除地形相位，生成一系列差分干涉圖，對差分干涉圖進行篩選得到高相干性像元，然后利用地面控制點進行軌道精煉，剔除低相干性和解纏錯誤的干涉像對，采用奇異值分解法和最小二乘法對多個小基線集聯合求解，從而得到研究區的時間形變序列及地表平均形變速率［25－26］。

大多數滑坡都是沿坡面滑動的，而小基線集干涉測量法只能監測雷達視線的形變，不足以反映滑坡的真實變形。因此，為了準確地描述滑坡真實的形變情況，將LOS向形變速率轉換為沿坡度方向形變速率［27－28］：

式中：V s l o p e為沿坡度方向的形變速率；V los為LOS向的形變速率；β為雷達視線方向與斜坡坡面的夾角，即視坡夾角；φ為斜坡的坡度/（°）；α為斜坡的坡向/（°）；θ為雷達入射角/（°）；αs為衛星軌道方向與正北方向的夾角，即衛星飛行的方向，升軌數據為負，降軌數據為正。

2.3 聚類和異常值分析原理

采用Anselin Local Moran′s I指數對SBAS-InSAR技術獲取的形變干涉點進行空間域的異常值分析和聚類分析［29～30］，通過計算要素與其他所有要素作為相鄰要素的數值關系，從而獲得數據集中有統計顯著性的高值（熱點）和低值（冷點），自動對形變點中的異常數據進行聚合，確定適當的空間分析范圍，并校正數據的空間相關性，得到形變聚類點置信水平p值和得分，其相關的計算公式如下：

統計量z I i得分的計算方法如下：

式中：I i為形變點i的Anselin Local Moran′s I指數值；x i、x j分別為點i，j的屬性；n為選定的窗口內要素的總數；是對應屬性的平均值；w i,j是要素i和j之間的空間權重；S2和S分別為為選定的窗口內像元的方差、標準差；μ0為總體均值；Z I i為形變聚類點i的得分；p i為形變聚類點i的顯著性置信水平。

3 數據處理

3.1 雷達可視性分析

雷達遙感通過向地面發送脈沖信號，接收地面反射的回波信號，實現對地物的觀測。盡管具有全天時、全天候、大面積等無法比擬的優勢，由于SAR采用側視成像方式，其入射角（θ）與實際地形存在夾角，使得SAR影像會出現陰影、透視收縮、疊掩3類常見的幾何畸變［27］。根據Sentinel－1A衛星降軌條件下的幾何參數及水城縣DEM提取的地形參數（坡度、坡向）可獲得水城縣可視性分區方式，并獲得水城縣可視性區域（見表2）。

表2 SAR影像可見性分區Table 2 SAR image visibility zoning

Sentinel－1A降軌影像在水城縣的可視性區域分布如圖3所示，可見區域約占35.76%，低敏感區域約占24.08%，不可見區域約占40.16%。其中，可見區域以主要分布于東南、西南區域；不可見區域主要分布在中心區域。

圖3 研究區雷達可視性分區Fig.3 Radar visibility zoning of in study area

3.2 基于改進SBAS技術的數據處理

文中使用精密軌道數據（POD）和30m分辨率的SRTM DEM進行影像軌道誤差和地形相位校正。為了提高時間、空間相干性，保證干涉圖質量，將時間基線設置為120 d，空間基線閾值為極限基線的5%（5 162.12×0.05=258.106 m），使用3D解纏。超級主影像為2018年7月3日，共生成了138個干涉像對，如圖4所示。通過查看生成的相干系數圖、去平地干涉圖、相位解纏圖，剔除相干性較差的干涉對，最終使用了115個干涉相對進行后續處理。

圖4 數據集時空基線分布Fig.4 Spatiotemporal baseline distribution of the dataset

圖5（a）是去平地后的相位干涉圖，水城縣由于植被覆蓋密集，地形起伏較大，導致地形相位去除效果一般，為了減小植被等因素的影響，選取了Delaunay MCF方法進行相位解纏；為了獲得較高的地面分辨率，將方位向與距離向的視數比設置為1：4；為了獲得更為清晰干涉條紋，選取了Goldstein濾波進行噪聲去除，濾波后的干涉圖結果如圖5（b），可看出濾波后噪聲誤差減小，干涉條紋的清晰度提高。

圖5 研究區部分干涉處理結果圖Fig.5 Results of partial interference processing in the study area

3.3 坡度方向可靠形變信息獲取

基于改進的SBAS-InSAR技術，獲取了2018年水城縣LOS方向的地表變形速率，對于西南山區，LOS向的形變速率不足以真實反映斜坡的真正形變情況，因此將LOS向形變速率轉換到坡度方向形變速率。當β接近90°時，cosβ近乎為0，導致V slop e趨于無窮大。為了對產生的異常值進行剔除，Herrera等以cosβ=±0.3為閾值，即V sl ope≤3.33V los。設定0＜C＜0.3時，C=0.3；-0.3＜C＜0時，C=-0.3［27－28，31－32］。V slope表示沿坡度方向的位移速率，塊體移動方向只能沿斜坡向下，故剔除V slope為非負的形變點。

基于水城縣可視性分區結果，先剔除了由透視收縮、疊掩和陰影不可信區域的形變信息，保留了好可視區域和低敏感區域中的形變信息。由于水城縣南部地區常年覆蓋著茂密的植被，盡管處于可見區域內，但C波段地表穿透能力有限，即便SAR數據時間間隔短，也會導致失相干，得到有效監測值較少且監測精度有限。此外，滑坡在河流及平坦的地區發生的可能性較低。因此，在可視性分區結果的基礎上，結合谷歌影像，將密度較高的林地、灌木林、河流以及坡度小于5°不可靠區域的形變信息進行剔除，保留其它區域可靠形變信息，得到水城縣坡度方向年均形變速率（圖6）。

圖6 研究區2018年坡度方向年均形變速率Fig.6 Average annual deformation rate of the slope direction in study area in 2018

4 結果與分析

4.1 水城縣變形監測結果分析

從圖6中可以看出整個水城縣北部的形變速率較為突出，各地形變速率不均，存在多個形變中心，集中分布在陡坡上。區域1形變較為嚴重，年平均速率在－75～0 mm/a，礦區開采處的變速率最明顯，形變速率達到了－142.3 mm/a；區域2整體地形起伏較大，形變程度也較為突出，其年平均形變速率在－35～0 mm/a之間；區域3受到雷達監測制約以及植被覆蓋影響，獲取的監測結果有限，形變速率大多在－10～0mm/a，少部分區域形變速率達到了－20 mm/a左右，整體較為穩定。

由上述分析可知，研究區內的變形區具有分布不均勻的特點，主要是人類活動因素（礦區開采、修筑道路）及差異化的地理環境因素（坡度、坡向、植被覆蓋率等）造成的。根據已有資料分析，區域1分布著較多的礦山，受礦區開采的影響，導致該區域整體形變量較大，易導致滑坡的發生。區域2地形起伏大、山體陡峭、植被覆蓋率低，受地形地貌等地理環境因素的影響，其山體邊坡的穩定性會明顯低于區域3，發生滑坡的可能性相對較高。

對獲取的形變信息進行統計分析，形變速率在－15～0 mm/a的形變點占據總數的40.5%，形變速率在－35～－15 mm/a的形變點占據總數的31.5%，形變速率在－55～－35 mm/a的形變點占據總數的13.5%，形變速率大于－55 mm/a的形變點占總數14.5%。從統計分析結果來看，形變速率大部分集中在－35～0 mm/a，其中形變速率小于－30 mm/a所占比例最大，在一定程度上能夠表明水城縣地表整體較為穩定，形變主要集中在水城縣的北部區域。。

根據上述對水城縣坡度方向形變速率分析，結合國內外眾多學者針對滑坡事故做出的研究成果［13－18，22－23］，文中選取V slop e=15 mm/a為速率閾值，即V slope＞15 mm/a的區域易發生滑坡事故。

4.2 水城縣滑坡識別及驗證

4.2.1 滑坡識別

通過Anselin Local Moran′s I對坡度方向形變點進行空間域異常值分析和聚類處理，會得到HH、HL、LH、LL四種聚類類型，其中HH表示高值聚集，HL表示低值包含高值異常，LH表示高值包含低值異常，LL表示低值聚集。通過獲得低值聚集區域，結合滑坡在谷歌影像上的形狀、分布、微地貌特征，確定滑坡位置及邊界。通過本文的方法最終共識別出了14處潛在滑坡，滑坡詳細信息見表3，滑坡識別結果分布見圖7，局部滑坡識別結果見圖8。

圖8 部分滑坡識別結果圖Fig.8 Some landslide recognition results

表3 研究區滑坡識別結果Table 3 Landslide identification in the study area

根據表3、圖7分析，識別出的滑坡有一半以上是由于劇烈的人類工程活動誘發導致，85%的滑坡類型為土質滑坡，在強降雨的情況下更容易出現失穩，結合H1、H5、H9的時序形變分析，在6～9月份形變趨勢處于急劇增加的態勢，是因為該期間迎來降雨高峰期，雨水富集在松散土層和相對防水的基巖接觸面之間，增加了土壤的自重，同時軟化了滑帶土，使其內聚力和內摩擦角減小，使得坡體的穩定性降低，為土層滑坡的形成奠定了基礎。

圖7 研究區滑坡識別分布圖Fig.7 Distribution map of landslide identification in study area

4.2.2 驗證分析

根據表3中潛在滑坡的高程值和坡度值進行分析，95%以上的滑坡分布在高程大于1 000 m，坡度大于20°的區域，借助前人對水城縣滑坡分布規律特征的研究成果進行對比驗證［20］，表明識別的潛在滑坡與水城縣滑坡分布規律特征相吻合，在一定程度上說明了識別結果的可靠性。

此外，通過貴州省地質災害排查項目及野外實地調查獲得了水城縣的22處滑坡，用于對本文識別的潛在滑坡進行驗證。驗證結果：識別獲得的14處潛在滑坡，有11處潛在滑坡位置和范圍與野外調查的滑坡較為吻合（占識別結果的78.6%），有3處潛在滑坡范圍與野外調查的結果存在一定偏差（占識別結果的21.4%），部分識別結果見圖8，另有8處滑坡未能識別。

根據圖7，圖8可知，對于成功識別的滑坡點位置，大多是位于坡度和形變程度較大的區域。未能成功識別的8處滑坡點，但也均監測到了部分形變區域，其未能識別的主要原因是：這些區域大多處于雷達不可見和不敏感區域，獲得有效形變點較少；部分地區植被覆蓋嚴重，為避免監測結果受大面積植被因素的影響，剔除了該區域的形變點，減少了該區域的有效時序形變點位提取量，導致無法準確識別。

5 結論與討論

5.1 結論

本文以滑坡多發的水城縣為研究區，探討了一種基于InSAR技術對潛在滑坡識別中的應用。使用歐空局提供的Sentinel－1 A SAR數據，利用改進的SBAS-InSAR方法獲取地表坡度方向形變信息，通過Anselin Local Moran′s I進行聚類分析，獲得高異常值區域，并結合滑坡的形狀、分布、微地貌特征成功識別出了14處潛在滑坡。將識別的潛在滑坡與野外勘探資料相核實，驗證了該時序InSAR用于滑坡隱患識別的有效性及優勢。

本研究受限于單一的降軌數據，對地形可視性差的區域無法進行有效識別，因此無法完全提取研究區內潛在滑坡體。文中是對整個縣城的進行時序監測，在處理數據上所花費時長較大，為了更快速的獲取識別結果，可利用升降軌數據直接對滑坡高易發區進行監測識別，從而有針對性地調查和監測高概率潛在的滑坡范圍，達到防止滑坡災害發生的目的。

5.2 討論

由于受到SAR系統監測方式的制約，在對西南高山區域的滑坡識別應用中，導致可探測的形變區域面積較少，難以全方位對潛在的地質災害隱患進行監測識別，限制了通過InSAR技術對潛在地質災害識別的應用。本研究獲取數據類型單一、有效觀測區覆蓋范圍較小，導致部分滑坡點未能進行監測，因此可以選用聯合升軌和降軌觀測減少觀測盲區，提升滑坡的有效探測率。

現有的民用雷達衛星通常工作于X波段（波長3.1cm）、C波段（波長5.6cm）或L波段（波長2.3.6cm），但是能大量獲取的免費SAR數據為C波段的Sentinel－1 A數據。對于InSAR形變測量來說，雷達波長決定測量靈敏度，即波長越短，監測精度越高。但是在高山峽谷區域，X波段和C波段的波長較短，難以穿透植被冠層，容易導致時間失相干；使用波段更長、穿透力更強的L波段進行觀測時，數據往往具有更好的相干性，在西南山區滑坡識別應用中更具有優勢。