基于SBAS/PS-InSAR技術(shù)的滑坡遙感監(jiān)測對比研究

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

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

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

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

1 研究區(qū)及數(shù)據(jù)簡介

1.1 研究區(qū)概況

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

1.2 數(shù)據(jù)源

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

表1 SAR影像參數(shù)

圖1 SAR數(shù)據(jù)成像范圍Fig.1 SAR data imaging range

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

2 數(shù)據(jù)處理流程及結(jié)果

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

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

(1)

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

圖2 技術(shù)方案Fig.2 Technical solution

2.1 PS-InSAR技術(shù)處理流程及結(jié)果

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

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

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

Disp=V(t-t0)

(2)

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

2.2 SBAS-InSAR技術(shù)處理流程及結(jié)果

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

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

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

3 白格滑坡形變監(jiān)測精度對比分析

3.1 基于PS-InSAR技術(shù)的白格滑坡形變監(jiān)測分析

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

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

由圖5可知，2種數(shù)據(jù)的PS-InSAR處理結(jié)果吻合度很高，其中，沿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種數(shù)據(jù)結(jié)果中，差值最大為28.9 mm/a，最小為0.1 mm/a，差值的平均值為14.5 mm/a。計算這194組數(shù)據(jù)的標(biāo)準(zhǔn)差，得到其范圍為0.06～14.4 mm/a，均值為7.2 mm/a。對于同一種PS-InSAR技術(shù)手段，Sentinel-1A和ALOS-2數(shù)據(jù)的結(jié)果存在一定的差異，但對于如此大形變量的滑坡來說，這個監(jiān)測結(jié)果已經(jīng)很理想，同時也驗證了這2種數(shù)據(jù)用于滑坡監(jiān)測的可行性。

3.2 基于SBAS-InSAR技術(shù)的白格滑坡形變監(jiān)測分析

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

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

由圖5、6發(fā)現(xiàn)，2種影像數(shù)據(jù)的形變趨勢大致相吻合，對于2種技術(shù)，Sentinel-1A和ALOS-2數(shù)據(jù)結(jié)果的差異較為穩(wěn)定。在所統(tǒng)計的共同點中，ALOS-2數(shù)據(jù)的大部分監(jiān)測結(jié)果的絕對值略微大于Sentinel-1A的結(jié)果，可能是因為ALOS-2影像L波具有波長優(yōu)勢，使得在一些植被覆蓋率高的區(qū)域監(jiān)測能力稍好。

3.3 同一數(shù)據(jù)源白格滑坡PS-InSAR與SBAS-InSAR形變監(jiān)測結(jié)果對比分析

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

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

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

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

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

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

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

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

(3)

利用式(3)計算A1～A4點2種技術(shù)結(jié)果的RMSE。經(jīng)計算，ALOS-2數(shù)據(jù)的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技術(shù)結(jié)果變化趨勢總體較為一致。

同樣從Sentinel-1A數(shù)據(jù)結(jié)果中選取4個特征點進(jìn)行時序分析，圖9、10為4個點的詳細(xì)信息。

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

對比發(fā)現(xiàn)，在2018-01-12～07-23期間，4個點的SBAS-InSAR和PS-InSAR監(jiān)測結(jié)果較為吻合，呈現(xiàn)近乎線性趨勢。在2018-07-23之后，兩者均出現(xiàn)波動，但不影響總體形變趨勢。到監(jiān)測結(jié)束(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數(shù)據(jù)源情況下，基于SBAS-InSAR和PS-InSAR技術(shù)監(jiān)測結(jié)果的精度差異，利用式(3)計算出B1～B4點的RMSE進(jìn)行比較分析。經(jīng)計算，Sentinel-1A數(shù)據(jù)的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可發(fā)現(xiàn)，Sentinel-1A數(shù)據(jù)結(jié)果的累積形變量值在-110～-80 mm之間，稍小于ALOS-2影像結(jié)果?？赡苁且驗镃波段的波長短，植被穿透能力相對L波段要弱，收到的回波信息量要少，使其結(jié)果受到一定程度的影響。

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

4 結(jié) 語

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