汪鈺晴，唐伶俐，王新鴻?，周增光，李子揚，李傳榮

(1 中國科學院空天信息創(chuàng)新研究院 中國科學院定量遙感信息技術重點實驗室， 北京 100094； 2 中國科學院大學電子電氣與通信工程學院， 北京 100049) (2020年6月12日收稿； 2020年7月6日收修改稿)

中巴經濟走廊(簡稱“中巴走廊”)北起中國新疆喀什，南至巴基斯坦瓜達爾港，全長約3 000 km，連接著北部的“絲綢之路經濟帶”和南部的“21世紀海上絲綢之路”，是貫通南北絲路的一條重要貿易走廊。中巴走廊沿線滑坡、泥石流、崩塌、溜石坡等地質災害頻發(fā)，不僅給沿線人民生活造成極大不便，也給走廊的工程建設和運營帶來了很大的挑戰(zhàn)。然而，由于該地區(qū)地形地質條件復雜、自然環(huán)境惡劣、安全狀況差、交通不便，不適合開展長期的野外調查活動，導致該區(qū)域相關研究工作進展緩慢。

20世紀70年代發(fā)展起來的衛(wèi)星合成孔徑雷達干涉測量技術(interferometric synthetic aperture radar，InSAR)是一種極具應用潛力的微波遙感技術，具有監(jiān)測范圍大、空間分辨率高、幾乎不受云雨天氣制約等優(yōu)點，已在地震形變[1-2]、火山運動[3]、冰川漂移[4]、城市沉降[5-6]、山體滑坡[7-8]以及線性工程形變[9]等研究領域得到廣泛應用。得益于星載SAR的發(fā)展，在相同研究區(qū)可積累大量的重復軌道SAR數據，這為開展長時間序列InSAR地表形變反演研究提供了可能。自本世紀初意大利米蘭理工大學Ferretti提出永久散射體方法后，一系列基于點目標的時間序列分析方法先后被提出，包括永久散射體(permanent scatterers，PS)方法[10]，小基線集(small baseline subsets，SBAS)方法[11]，IPTA(interferometric point target analysis)方法[12]，StaMPS(stanford method for permanent scatterers)方法[13]等，這些方法通過分析時間序列SAR影像，提取出相位和幅度信息在長時間范圍內保持穩(wěn)定的離散點，即PS點或相干目標候選點(例如城市建筑、道路、裸露的巖石等地物)，利用這些點目標的相位特征，進行長期緩慢地表形變的反演。時序InSAR技術不僅可以克服傳統差分合成孔徑雷達干涉測量技術(differential InSAR，D-InSAR)易受時間、空間去相干影響的問題，還能在一定程度上估計大氣影響和優(yōu)化初始數字高程模型(digital elevation model，DEM)精度，達到抑制大氣延遲相位和地形誤差相位的效果。大量研究表明，時序InSAR技術能提取大范圍、高精度的長時間緩慢地表形變信息，可廣泛應用于城市地面沉降[14]、潛在滑坡形變監(jiān)測[15-16]、地質災害識別[17]等領域。

本文選擇中巴走廊沿線地質災害最為發(fā)育的蓋孜河谷段為研究區(qū)，利用2015年4月5日至2018年12月9日共30景Sentinel-1A干涉寬幅模式下的雷達影像，分別采用PS-InSAR技術和SBAS-InSAR技術對研究區(qū)進行地表形變信息提取，比較兩種時序InSAR技術所得形變序列結果，對地質災害形變序列特征進行初步分析。

1 研究區(qū)概況與數據

1.1 研究區(qū)概況

中巴走廊穿過帕米爾高原腹地，地質構造上屬于印度-巴基斯坦板塊與歐亞板塊強烈碰撞和擠壓的區(qū)域，是地殼隆起迅速、地震活動強烈、溝谷侵蝕劇烈的典型地區(qū)之一[18]。因構造運動劇烈、地勢陡峭、冰川廣泛發(fā)育，中巴走廊沿線地質災害頻發(fā)，總體上具有災種多、規(guī)模大、頻度高、分布廣的特點[19-20]。

蓋孜河谷橫切西昆侖山脈，北靠蓋昆山(最高海拔6 670 m)，南依公格爾山(最高海拔7 649 m)，河谷大致呈北東東(NEE)走向，平均海拔2 500 m，最寬處不到2 km，最窄處不足50 m，是中巴走廊的咽喉地段。蓋孜河谷位于歐亞大陸腹地，降水相對稀少而蒸發(fā)較為強烈，屬于典型的大陸性氣候。由于氣候寒冷干旱、冰川廣泛發(fā)育、風化作用強烈，該區(qū)域植被稀少、山體裸露。滑坡、崩塌、泥石流等地質災害的頻繁發(fā)生，給走廊建設和人民生活帶來很大威脅[21]。

使用星載SAR數據進行時序InSAR處理，考慮到整景SAR影像覆蓋范圍較大，為提高數據處理效率，在進行時序InSAR處理前先對SAR影像做空間位置裁剪，以提取感興趣的蓋孜河谷段區(qū)域，裁剪范圍為74°58′～75°30′ E，38°40′～38°53′ N(約45 km×25 km)，所得研究區(qū)地理位置如圖1所示。

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

1.2 SAR數據

選擇Sentinel-1A衛(wèi)星SAR數據作為數據源。Sentinel-1A衛(wèi)星于2014年4月3號發(fā)射，采用太陽同步軌道(軌道高度693 km)，具有條帶、干涉寬幅、超寬幅及波浪4種成像模式[22]，主要參數如表1所示。

以中巴走廊蓋孜河谷段部分區(qū)域為研究區(qū)，共選取30景Sentinel-1A衛(wèi)星SAR影像(干涉寬幅(IW)模式、升軌、VV極化)。影像獲取時間如表2所示，時間范圍為2015年4月5日至2018年12月9日。

表1 Sentinel-1A衛(wèi)星主要參數Table 1 Main parameters of Sentinel-1A

表2 SAR數據時間分布表Table 2 Date distribution of the SAR data

此外，Sentinel-1A衛(wèi)星發(fā)布了精密定軌星歷參數(precise orbit ephemerides，POD)用于修正軌道信息，定位精度優(yōu)于5 cm，將它作為輔助數據可以有效去除軌道誤差對形變測量的影響。

1.3 DEM數據

采用空間分辨率為30 m的ASTER GDEM V2數據(下載網址https:∥www.gscloud.cn)作為InSAR處理的輔助數據。該數據集根據NASA的EOS-Terra衛(wèi)星的全球觀測結果制作而成，其數據覆蓋南北緯83°之間的所有陸地區(qū)域。該高程數據集的垂直精度為20 m，水平精度為30 m，可在InSAR處理過程中較好地模擬地形相位。

2 地表形變提取方法

分別采用PS-InSAR技術和SBAS-InSAR技提取研究區(qū)地表形變信息，術總體流程如圖2所示。

PS-InSAR技術提取地表形變的主要步驟包括：主影像選擇、差分干涉處理、PS候選點選擇、估算候選點的平均形變速率及DEM校正系數、大氣相位的估計和去除、形變時間序列估計等。本文選取2017年3月13日的SAR影像作為公共主影像，其他29幅影像作為輔影像。圖3(a)為干涉對的時間基線和空間基線連接圖，從中可見主輔影像干涉像對空間基線分布于-100～105 m，時間基線范圍為-708～636 d。選取公共主影像之后，將其余29幅輔影像配準重采樣到主影像的像素空間上并進行主輔影像干涉。利用Sentinel-1A衛(wèi)星精密定軌星歷參數與輔助DEM數據模擬并去除參考橢球面相位及地形相位后，得到29幅差分干涉圖。

采用振幅離差指數方法，將振幅離差指數閾值設置為0.4，對研究區(qū)內PS點進行識別。PS點的干涉相位包含線性形變、高程誤差、非線性形變、大氣延遲和噪聲等相位分量。首先用線性模型從所有差分干涉圖中估算線性形變速率和殘余高程信息，再在時域和空域上做相應的濾波處理將非線性形變相位分離出來，將它與線性形變分量相加，得到每個PS點的總地表形變信息。最后，將視線向形變進行地理編碼，獲得地理坐標系下的地表時序形變信息。

圖2 研究區(qū)地表形變提取流程圖Fig.2 Flow chart of surface deformation extraction in the study area

圖3 PS及SBAS處理中SAR數據對的時空基線連接圖Fig.3 Temporal-spatial baseline connection diagram of SAR data pairs in PS and SBAS processing

SBAS-InSAR技術提取地表形變的主要步驟包括：連接圖生成、差分干涉處理及相位解纏、軌道精煉與重去平、大氣相位估計與去除、形變時間序列估計等。本文選取2015年8月27日的SAR影像作為公共主影像。經過主輔影像配準得到像空間坐標系一致的30幅SLC影像后，設置時空基線閾值。本文使用的SAR數據總時間跨度為1 344 d，為保證干涉質量，將空間最大臨界基線百分比設置為45%，時間基線閾值設置為180 d，一共生成94個干涉對，最長的空間基線約為210 m。圖3(b)為生成的時空基線連接圖，這里每幅影像平均有6.27個連接，連接得比較均勻，且每個時相都與其他時相建立了連接，表明數據配對良好。

對滿足時空基線閾值的所有干涉像對進行差分干涉處理，基于衛(wèi)星精密定軌星歷參數與輔助DEM數據，去除參考橢球面相位及地形相位，然后對差分干涉圖進行濾波處理，改善干涉條紋的清晰度，之后通過相位解纏得到2次觀測期間同一目標的真實相位差。最后經過軌道精煉與重去平、形變速率與DEM校正速率估計、大氣效應相位估計與去除等處理，得到所有PS點的年平均形變速率及時間序列形變信息。本文采用ENVI軟件的SARscape模塊完成上述時序InSAR處理過程。

3 結果與分析

3.1 PS-InSAR與SBAS-InSAR形變結果比較分析

圖4(a)和4(b)分別為利用PS-InSAR及SBAS-InSAR技術獲得的研究區(qū)PS點分布及其年平均形變速率圖，結合光學影像(這里的背景底圖)可以看出，兩種方法提取PS點的分布范圍較為一致，主要分布在海拔較低的中巴公路沿線、西部盆地及坡度較緩山區(qū)，在冰雪覆蓋區(qū)基本無PS點，受Sentinel-1A衛(wèi)星雷達波入射角和入射方向影響，由于陰影效應，雷達波束照射不到的北東(NE)向陡峭坡面上也沒有形成PS點。

PS-InSAR處理時全區(qū)共提取了305 853個PS點，所測形變?yōu)閭鞲衅饔^測方向的形變，即雷達視線向(line of sight，LOS)形變量(地面向衛(wèi)星運動方向為正，遠離衛(wèi)星運動方向為負)，PS點的LOS形變速率主要分布于-98～98 mm/a區(qū)間范圍內。PS點的形變速率直方圖大致為正態(tài)分布，如圖5(a)所示，均值約為-0.32 mm/a，標準差約為5.25 mm/a。從圖4(a)可以看出，研究區(qū)西部盆地及中巴公路沿線大部分PS點年平均形變速率都較小，表示這些區(qū)域在監(jiān)測時間段內(2015年4月—2018年12月)未發(fā)生明顯形變；而在中巴公路沿線部分路段及公路兩側冰川前緣，分布有年平均形變量較大的PS點(圖4(a)中矩形區(qū)域)，表示這些區(qū)域在監(jiān)測時間段內可能發(fā)生了由于地質災害而引發(fā)的形變。此外，提取出的PS點中還存在一些孤立的年平均形變速率絕對值異常大的點簇(圖4(a)中橢圓形區(qū)域)，由于研究區(qū)地處高寒干旱山區(qū)，高程落差大且氣候變化明顯，這些LOS形變異常高值PS點可能是在PS-InSAR數據處理過程中由于大氣相位或DEM殘差相位分離不當造成的。

基于相同的30幅雷達影像，在同一區(qū)域采用SBAS-InSAR技術共提取了670 757個PS點，測得PS點的LOS年平均形變量主要位于-46～38 mm/a區(qū)間范圍內。PS點的形變速率直方圖接近正態(tài)分布，如圖5(b)所示，均值約為0.15 mm/a，標準差約為2.86 mm/a。從圖4(b)可以看出，研究區(qū)西部盆地及中巴公路沿線大部分PS點在監(jiān)測時間段內保持穩(wěn)定，僅在公路沿線部分斜坡及兩側冰川前緣分布有形變量較大的PS點，這與PS-InSAR所得結果較為一致，兩種時序InSAR技術可以相互印證。

圖4 PS-InSAR及SBAS-InSAR技術提取的蓋孜河谷LOS形變速率圖Fig.4 LOS deformation rate map of the Gaizi valley by PS-InSAR and SBAS-InSAR techniques

圖5 兩種時序InSAR技術所提取PS點的形變速率分布圖Fig.5 Distribution of LOS velocity obtained by the two time series InSAR techniques

值得注意的是，PS-InSAR技術是通過計算每個像素點的振幅離差指數來選擇PS候選點，進而通過設定多時間相干系數閾值得到所有PS點的，這里多時間相干系數表示有多少形變趨勢符合線性模型，該閾值設置越大，所得PS點越少(本文將此閾值設定為0.8)。與PS-InSAR技術利用強度信息確定PS點不同，SBAS-InSAR技術是通過設定相關系數閾值來獲取高相干點的，本文處理過程中將相關系數閾值設置為0.35，低于該閾值的像素以空值輸出。由于PS點的識別方法不同，這兩種時序InSAR技術所獲得的PS點在數量及分布上有較大差距。

相較而言，SBAS-InSAR技術更充分地利用了小基線數據子集內的良好相干性，提取出的時序形變結果在空間上更為連續(xù)，減少了出現形變異常結果的可能性；此外，在測量形變類型上，PS-InSAR主要適用于線性形變，而SBAS-InSAR對線性及非線性形變均適用。因此，下面將基于SBAS-InSAR提取結果對研究區(qū)部分潛在地質災害進行形變特征分析。

3.2 研究區(qū)地質災害形變特征分析

結合PS點正態(tài)分布圖、PS點誤差等綜合分析，可將形變速率劃分為7個區(qū)間，其中極小值-46～-20 mm/a與極大值20～38 mm/a用來識別LOS形變異常高值點，中間區(qū)間-5～5 mm/a可以看作背景值，表示該區(qū)域地表較為穩(wěn)定，其余4個區(qū)間：-20～-10、-10～-5、5～10、10～20 mm/a可以用來標記產生一定程度形變的區(qū)域。

從圖6形變速率劃分結果可發(fā)現，研究區(qū)西部盆地與中巴公路沿線大部分區(qū)域年平均形變速率均不超過±5 mm/a，表示研究區(qū)大部分區(qū)域在監(jiān)測時間段內保持穩(wěn)定，未發(fā)生明顯地表形變；而在中巴公路沿線兩側斜坡上，識別出較多年平均形變量大于±5 mm/a的PS點，由于公路沿線兩側山體植被稀少、風化作用強烈，易產生崩塌、溜石坡等斜坡災害；在布倫口南北兩側冰川前緣區(qū)域，分布有大量年平均形變速率較大PS點，部分PS點形變速率超過±20 mm/a，這是由冰川運動或冰川融水引發(fā)泥石流所造成。

結合研究區(qū)光學影像及前人相關研究成果和經驗，根據地表LOS形變速率提取出研究區(qū)內多處不穩(wěn)定斜坡及冰川運動，見圖6。其中形變速率較大的不穩(wěn)定斜坡主要分布在中巴公路兩側的坡面上(圖6中橢圓區(qū)域)，存在大量不穩(wěn)定PS點，個別不穩(wěn)定斜坡與公路距離較近，可能發(fā)生損毀道路的安全隱患。公格爾山是研究區(qū)內海拔最高的山峰，冰川發(fā)育較為顯著，由于冰川運動，在冰川前緣提取到較多形變速率相對高值PS點(圖6中矩形區(qū)域)，此外，在公路北側木吉縣方向也有比較明顯的冰川運動。

圖6 蓋孜河谷區(qū)域SBAS-InSAR處理所得PS點分布及形變速率Fig.6 Distribution and velocities of PS points obtained by SBAS-InSAR processing in the Gaizi valley

3.2.1 不穩(wěn)定斜坡形變特征

研究區(qū)構造運動強烈、河流切割侵蝕、地震發(fā)生頻繁，且氣候高寒、冰川廣泛發(fā)育，公路沿線兩側產生了大量陡坡，易發(fā)生滑坡災害，對公路存在潛在危害。由地震、降雨、風化而引發(fā)的崩塌也是研究區(qū)主要地質災害之一，規(guī)模一般較小，但對公路本身及過往人員車輛有較大威脅。同時威脅公路安全的還有由于巖石風化而形成的溜石坡。由于時序InSAR技術獲取的地表信息有限，較難準確區(qū)分出滑坡、崩塌及溜石坡，這里將這些地質災害統一稱作不穩(wěn)定斜坡。圖7為中巴公路北側西部一處不穩(wěn)定斜坡的形變特征圖，可以看到坡面上分布有大量負向形變量較大的PS點，累積最大負向形變達到-30 mm，且越靠近坡底形變量越小，在2016年6月之前未發(fā)生明顯形變，之后大致呈線性趨勢產生負向形變；坡底分布有少量正向形變PS點，累積正向形變達到20 mm，可能由碎石堆積引起，與不穩(wěn)定斜坡的一般性特點相符。

圖7 不穩(wěn)定斜坡形變特征圖Fig.7 Deformation characteristics diagram of unstable slope

3.2.2 冰川運動形變特征

由于蓋孜河谷獨特的地質地貌和氣候條件，冰川運動及其引發(fā)的冰川融水泥石流是蓋孜河谷區(qū)域的典型地質災害之一，產生的形變相對較大，表現為年均形變量高值區(qū)域，地域上主要分布于各大溝谷中、國道、鄉(xiāng)道公路兩側，時間上多集中于春夏融水強度較大的5—8月。研究區(qū)地勢陡峭、冰雪融水豐富且具有豐富的松散固體源，冰川泥石流是蓋孜河谷危害最為嚴重的地質災害，但由于泥石流運動速度快，形變大，會引起嚴重失相干，在災害點往往無法提取到有效的形變信息，但冰川泥石流會造成坡底碎石堆積，這在平均形變速率圖上表現為正異常值集中區(qū)域。冰川前緣運動一般比較緩慢，可提取到較多有效PS點，故可根據冰川前緣冰磧物的形變特征合理推測冰川運動情況。由圖6底圖光學影像可看出，冰川在蓋孜河谷南側公格爾山及北側蓋昆山均廣泛發(fā)育，圖8為蓋孜河谷北部一處冰川運動形變特征圖，圖8(a)為PS點分布及形變速率，圖8(b)為矩形框和橢圓框內部PS點形變時間序列變化曲線，體現出冰漬物運動引發(fā)的地表形變，其中正負向累積形變分別達到+60 mm和-80 mm。由于該地區(qū)的土質多為季節(jié)性凍土，土質極不穩(wěn)定，冰川運動可能會加劇不穩(wěn)定斜坡的發(fā)生，甚至引發(fā)冰磧物滑坡。

圖8 冰川運動形變特征圖Fig.8 Deformation characteristics diagram of glacial movement

4 結語

本文以中巴經濟走廊蓋孜河谷段為研究區(qū)，利用PS-InSAR及SBAS-InSAR技術分別進行地表形變提取，比較兩種時序InSAR技術所得結果發(fā)現，這兩種方式獲取的形變信息在形變速率的空間變化趨勢上具有較好的一致性，其中SBAS-InSAR方法能更充分地利用小基線數據子集內的良好相干性，所獲取的時序形變速率分布在空間上更為連續(xù)。根據SBAS-InSAR技術提取的結果可知研究區(qū)多處發(fā)生明顯地表形變，主要分布于中巴公路兩側斜坡及布倫口以北的冰川前緣區(qū)域。結合光學影像識別出研究區(qū)內多處不穩(wěn)定斜坡與冰川運動，并對其中典型實例進行了時序形變特征分析。不過，由于研究區(qū)歷史上的實地勘測數據十分稀缺，所獲取的形變結果尚難以進行完全的驗證。本文所做的前期探索性研究初步展現了時序InSAR技術在高寒山區(qū)地質災害監(jiān)測中的應用潛力，有助于今后衛(wèi)星數據與實地勘測同步開展條件下的進一步深入研究，并最終服務于蓋孜河谷區(qū)域地質災害易發(fā)性評價、基礎工程設施建設與防護等工程應用。