基于SBAS/PS-InSAR技術的滑坡遙感監測對比研究

付波霖 解淑毓 李 濤 李 豪 左萍萍 高二濤

1 桂林理工大學測繪地理信息學院，桂林市雁山街319號，541006 2 廣西空間信息與測繪重點實驗室，桂林市雁山街319號，541006

合成孔徑雷達干涉(interferometric synthetic aperture radar, InSAR)測量技術因具有能夠全天時、全天候工作、以cm量級甚至更小尺度精確獲得地表形變等優勢，而被用于滑坡隱患的早期識別與形變監測中[1-2]。在此基礎上發展的永久散射體(permanent scatters InSAR, PS-InSAR)測量技術和小基線集(small baseline subset InSAR, SBAS-InSAR)測量技術能有效克服傳統D-InSAR技術時空去相關的影響，已被廣泛應用于大面積、連續性的滑坡監測[3-5]。

本文以西藏自治區江達縣波羅鄉白格滑坡點為研究區，以2015～2019年12景ALOS-2 PALSAR2影像和2018～2019年38景Sentinel-1A影像為主要數據源，結合90 m分辨率的SRTM3 DEM，經PS-InSAR和SBAS-InSAR技術處理后獲取白格滑坡的形變信息，分析并驗證這2種監測技術在滑坡形變監測應用中的有效性。

1 研究區及數據簡介

1.1 研究區概況

研究區白格滑坡位于西藏江達縣白格村區域的金沙江岸，是一個地質活動頻繁、山體結構相對破碎的地帶，其地理位置為98°42′15.21″E、31°4′56.26″N(圖1中黃色點處)。白格滑坡地處青藏高原東部，橫斷山脈北部，位于芒康山與沙魯里山之間的金沙江河谷地帶，坡體最高處海拔達到3 719 m。區域地質調查資料顯示，白格滑坡位于金沙江構造結合帶，主要由構造混雜體組成。由于地質活動較多，容易受到降雨的影響，而該地區多年年平均降雨量達660 mm，時間集中在6～9月[6]。

1.2 數據源

本文圍繞白格滑坡點建立研究區，選取12景L波段ALOS-2 PALSAR2影像和38景C波段Sentinel-1A影像作為實驗數據。表1為2種數據的具體參數，2種影像數據范圍如圖1所示。

表1 SAR影像參數

圖1 SAR數據成像范圍Fig.1 SAR data imaging range

本研究下載的90 m分辨率DEM是由美國太空總署(NASA)和國防部國家測繪局(NIMA)聯合測量的SRTM數據，SRTM數據包括2種類型：一種是SRTM1，分辨率為30 m；另一種則是本文所使用的SRTM3，分辨率為90 m。

2 數據處理流程及結果

雷達干涉的模式包括重復軌道干涉法、交叉軌道法和沿軌道干涉法3種，其中最常用的是重復軌道干涉法，其干涉相位表達式為[7]：

φint=φfla+φtop+φdef+φnoi+φatm

(1)

式中，φfla為平地相位，φtop為地形造成的相位誤差，φdef為地表形變引起的相位，φnoi為觀測噪聲引起的相位，φatm為大氣延遲誤差。由于噪聲相位的隨機性，采用低通相位濾波方法可以去除噪聲相位。同時，大氣延遲相位具有較高的地形空間相關性，在相位模型中呈現低頻信息，因此，該相位的高通濾波可以減弱大氣相位效應。圖2為本文基于SBAS/PS-InSAR的技術方案。

圖2 技術方案Fig.2 Technical solution

2.1 PS-InSAR技術處理流程及結果

PS-InSAR技術由Ferretti等[8]首次提出，該技術利用長期保持穩定反射特性的永久散射體降低數據時空相干性，糾正大氣效應。其主要用于分析點目標、結果與點目標線性形變的相關性，因此需要干涉條件和輻射條件比較穩定的區域的影像數據才能保證結果質量。數據處理過程主要為：先選取1景影像為主影像，其余為輔影像，設置500%的基線閾值，生成11個數據對(ALOS-2)和37個數據對(Sentinel-1A)，得到的空間基線如圖3所示。

圖3 PS-InSAR時空基線分布Fig.3 Distribution of PS-InSAR spatial-temporalbaselines

然后根據數據對的連接關系對所有的影像進行干涉生成。在候選PS點的時候利用振幅離差法設置的閾值0.75進行穩定目標點的選擇。利用第1次模型的反演獲得殘余地形和位移速率，該模型為線性模型，可簡潔地表達為：

Disp=V(t-t0)

(2)

式中，Disp為t時間的形變，V為形變速率。經高低通濾波去除大氣相位后，通過時序分析便得到PS點在時間序列上的形變信息。

2.2 SBAS-InSAR技術處理流程及結果

SBAS-InSAR技術是由Berardino等[9]提出的一種時間序列分析方法，主要用于分析分布式目標，得到一個空間上連續的面結果。首先通過選取1景主影像(ALOS-2為2017-11-27,Sentinel-1A為2018-08-04)生成43對ALOS-2和153對Sentinel-1A干涉像對；然后通過手動選取GCP點對所有數據進行去平處理。為了估算和去除干涉處理中解纏后的相位圖中依然存在的恒定相位，點對應盡量選擇在地形平坦、沒有形變條紋的區域，因為控制點是被認為形變為0的點；最后設置相關系數閾值和第2次解纏相干系數閾值為0.2，采用矩陣奇異值分解方法(SVD)進行第1次估算形變速率和殘余地形。圖4為1對像對在第2次解纏后的結果，其中紅圈為滑坡形變區域。

圖4 Sentinel-1A像對(20180313/20180205)第2次相位解纏Fig.4 Second phase unwrapping of sentinel-1A image pair(20180313/20180205)

設置相關系數閾值為0.2、時間域高通濾波窗口為365 d、空間域低通濾波窗口為1 200 m，進行大氣相位去除，得到時間序列上的形變信息。

3 白格滑坡形變監測精度對比分析

3.1 基于PS-InSAR技術的白格滑坡形變監測分析

在采用PS-InSAR技術進行數據處理時，從2種數據結果中的滑坡發生區域選取194個地理位置相同的點進行對比分析(圖5)，其中正值表示向衛星方向運動，運動趨勢是沿雷達視線向(line of sight，LOS)上升；負值表示遠離衛星方向運動，運動趨勢是沿LOS向下降。

圖5 Sentinel-1A和ALOS-2數據194個共同點的PS-InSAR技術監測結果對比Fig.5 Comparison of PS-InSAR monitoring results at 194 common points of Sentinel-1A and ALOS-2 data

由圖5可知，2種數據的PS-InSAR處理結果吻合度很高，其中，沿LOS向上升最大速率為37.9 mm/a(ALOS-2)和24.2 mm/a(Sentinel-1A)，沿LOS向下降最大速率為-68.9 mm/a(ALOS-2)和-64.5 mm/a(Sentinel-1A)。2種數據結果中，差值最大為28.9 mm/a，最小為0.1 mm/a，差值的平均值為14.5 mm/a。計算這194組數據的標準差，得到其范圍為0.06～14.4 mm/a，均值為7.2 mm/a。對于同一種PS-InSAR技術手段，Sentinel-1A和ALOS-2數據的結果存在一定的差異，但對于如此大形變量的滑坡來說，這個監測結果已經很理想，同時也驗證了這2種數據用于滑坡監測的可行性。

3.2 基于SBAS-InSAR技術的白格滑坡形變監測分析

從2種數據的SBAS-InSAR結果中提取197個地理位置相同的點進行對比分析(圖6)。由圖6可見，ALOS-2數據的平均形變速率范圍為-84.2～-40.0 mm/a，Sentinel-1A數據的平均形變速率范圍為-84.0～-13.0 mm/a。2種數據結果中差值最大為29.9 mm/a，最小為0.03 mm/a，差值的平均值為17.1 mm/a。計算這197組數據的偏離程度，得到的標準差范圍為0.02～15.0 mm/a，絕大部分點位的標準差均在10 mm/a以下，均值為8.5 mm/a。

圖6 Sentinel-1A和ALOS-2數據197個共同點的SBAS-InSAR技術監測結果對比Fig.6 Comparison of SBAS-InSAR monitoring results at 197 common points of Sentinel-1A and ALOS-2 data

由圖5、6發現，2種影像數據的形變趨勢大致相吻合，對于2種技術，Sentinel-1A和ALOS-2數據結果的差異較為穩定。在所統計的共同點中，ALOS-2數據的大部分監測結果的絕對值略微大于Sentinel-1A的結果，可能是因為ALOS-2影像L波具有波長優勢，使得在一些植被覆蓋率高的區域監測能力稍好。

3.3 同一數據源白格滑坡PS-InSAR與SBAS-InSAR形變監測結果對比分析

為了檢驗同一數據源在不同技術手段下結果的精度情況，在ALOS-2數據結果中選取4個特征點進行時序分析。圖7、8為4個特征點的詳細信息。

圖7 ALOS-2結果中4個特征點的空間分布Fig.7 Spatial distribution of 4 feature points in ALOS-2 results

圖8 ALOS-2結果中4個特征點的累積形變量Fig.8 Cumulative deformation of 4 feature points in ALOS-2 results

觀察發現，4個特征點上PS-InSAR和SBAS-InSAR結果的總體形變趨勢一致。在2015-07～2018-10期間，A1特征點上形變速率呈加速趨勢，累積形變量最大達到-288.1 mm(SBAS-InSAR)和-287.6 mm(PS-InSAR),隨后在2018-11、12出現沿LOS向小幅上升，上升量為51.6 mm(SBAS-InSAR)和45.6 mm(PS-InSAR)。此后，A1點位整體保持沿LOS向下降的趨勢。

A2點在2019-02-18和2018-10-29達到最大累積形變量-319.4 mm(SBAS-InSAR)和-315.4 mm(PS-InSAR)，一直到2018-11-26，A2點位的SBAS-InSAR和PS-InSAR結果均出現沿LOS向小幅上升，形變量為11.1 mm(SBAS-InSAR)和39.4 mm(PS-InSAR)。

A3點在2015-07-27～2016-07-25期間呈現加速形變趨勢，形變量約-150 mm。在2018-07-23～11-26期間存在沿LOS向小幅上升，至2019-04-15，該點位累積形變量達到最大值-328.3 mm(SBAS-InSAR)和-315.9 mm(PS-InSAR)。

A4點到2018-05-28為止累積形變量已達到-295 mm(SBAS-InSAR)和-331 mm(PS-InSAR)。2018-07-23～2019-03-18其SBAS-InSAR結果達到最大累積形變量為-357.7 mm。2018-07-23～2019-02-18其PS-InSAR結果最大累積形變量為-342.6 mm。

為了能夠更加準確地分析同數據源不同數據結果的精度情況，利用圖8中的差值數據，使用均方根誤差(RMSE)進行分析：

(3)

利用式(3)計算A1～A4點2種技術結果的RMSE。經計算，ALOS-2數據的A1～A4點的RMSE范圍分別為0～16.6 mm、0～21.8 mm、0～17.2 mm、0～15.3 mm，均值分別為8.6 mm、8.3 mm、8.4 mm、8 mm。其中，A1點上的最大誤差為16.6 mm，A2點上的最大誤差為21.8 mm，A3點上的最大誤差為17.2 mm，A4點上的最大誤差為15.3 mm。觀察得出，2種InSAR技術結果變化趨勢總體較為一致。

同樣從Sentinel-1A數據結果中選取4個特征點進行時序分析，圖9、10為4個點的詳細信息。

圖10 Sentinel-1A結果中4個特征點的累積形變量Fig.10 Cumulative deformation of 4 feature points in Sentinel-1A results

對比發現，在2018-01-12～07-23期間，4個點的SBAS-InSAR和PS-InSAR監測結果較為吻合，呈現近乎線性趨勢。在2018-07-23之后，兩者均出現波動，但不影響總體形變趨勢。到監測結束(2019-04-13)為止，B1點上的最終累積形變量為-76.1 mm(SBAS-InSAR)和-76.7 mm(PS-InSAR)，二者差值為0.6 mm；B2點上的最終累積形變量為-79.8 mm(SBAS-InSAR)和-79.2 mm(PS-InSAR)，二者差值為0.6 mm；B3點上的最終累積形變量為-98.2 mm(SBAS-InSAR)和-96.2 mm(PS-InSAR)，二者差值為-2 mm；B4點上的最終累積形變量為-102.4 mm(SBAS-InSAR)和-93.1 mm(PS-InSAR)，二者差值為-9.3 mm。

同樣，為了能更精確地了解同為ALOS-2數據源情況下，基于SBAS-InSAR和PS-InSAR技術監測結果的精度差異，利用式(3)計算出B1～B4點的RMSE進行比較分析。經計算，Sentinel-1A數據的B1～B4點的RMSE范圍分別為0～7.4 mm、0～8.2 mm、0～9.5 mm、0～5.2 mm。其中，B1點上的最大誤差為7.4 mm，日期為2018-08-28；B2點上的最大誤差為8.2 mm，日期為2018-11-08；B3點上的最大誤差為9.5 mm，日期為2019-04-01；B4點上的最大誤差為5.2 mm，日期為2018-10-03。4個特征點的RMSE均值分別僅為3 mm、2.9 mm、2.9 mm、2 mm。觀察圖10可發現，Sentinel-1A數據結果的累積形變量值在-110～-80 mm之間，稍小于ALOS-2影像結果。可能是因為C波段的波長短，植被穿透能力相對L波段要弱，收到的回波信息量要少，使其結果受到一定程度的影響。

綜上，Sentinel-1A數據的4個特征點均保持線性形變趨勢，且PS-InSAR和SBAS-InSAR結果趨于一致，驗證了2種時序InSAR技術應用于滑坡監測中的可靠性。

4 結 語

本文以西藏江達縣波羅鄉白格滑坡點為研究對象，以12景ALOS-2 PALSAR2影像以及38景Sentinel-1A影像為主要數據源，通過SBAS/PS-InSAR技術反演研究區域的地表形變速率及時間序列形變信息。得到結果如下：基于ALOS-2數據的PS-InSAR監測到的最大平均形變速率為-68.9 mm/a，SBAS-InSAR監測到的最大平均形變速率為-84.2 mm/a；基于Sentinel-1A數據的PS-InSAR監測到的最大平均形變速率為-64.5 mm/a，SBAS-InSAR監測到的最大平均形變速率為-84.0 mm/a。