于 冰 王 楊 馬德英 蔣榮乾 張 過 周志偉

1 西南石油大學(xué)土木工程與測繪學(xué)院，成都市新都大道8號，6105002 西南石油大學(xué)油氣藏地質(zhì)及開發(fā)工程國家重點(diǎn)實(shí)驗(yàn)室，成都市新都大道8號，6105003 西南石油大學(xué)油氣空間信息工程研究所，成都市新都大道8號，6105004 武漢大學(xué)測繪遙感信息工程國家重點(diǎn)實(shí)驗(yàn)室，武漢市珞喻路129號，4300795 中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院大地測量與地球動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢市徐東大街340號，430077 6 中鐵二院工程集團(tuán)有限責(zé)任公司，成都市通錦路3號，610031

干涉點(diǎn)目標(biāo)分析(interferometric point target analysis, IPTA)方法[1-2]是對PS-InSAR技術(shù)的發(fā)展，其通過建立相位模型，計算并分離PS點(diǎn)的各相位分量，再根據(jù)結(jié)果不斷精化模型，迭代求解研究區(qū)的形變信息。IPTA方法在地表沉降、基礎(chǔ)設(shè)施形變等方面應(yīng)用較為廣泛[3-5]，但在滑坡形變監(jiān)測方面應(yīng)用較少。

本文利用C波段Sentinel-1A影像數(shù)據(jù)，采用IPTA方法對大光包滑坡進(jìn)行形變監(jiān)測，提取該區(qū)域2016～2018年的形變信息，并結(jié)合全球降水測量(global precipitation measurement, GPM)計劃日均降水量數(shù)據(jù)，對特征點(diǎn)進(jìn)行時序形變分析，將該方法獲得的結(jié)果與StaMPS-SBAS方法的處理結(jié)果進(jìn)行對比，分析IPTA方法應(yīng)用于大型滑坡形變監(jiān)測的可行性與潛力，為大光包滑坡的監(jiān)測與預(yù)警提供技術(shù)和信息支撐。

1 IPTA方法基本原理

IPTA方法的核心思想是對相干點(diǎn)目標(biāo)進(jìn)行處理，通過構(gòu)建和更新二維線性相位回歸模型，迭代求解相鄰相干點(diǎn)的線性形變速率增量和高程改正值，直至逼近線性形變速率和高程誤差的最優(yōu)解[1-2]。在IPTA處理中，根據(jù)影像幅度和相位信息，計算每幅影像的光譜相干系數(shù)和時序振幅穩(wěn)定性，綜合選取高質(zhì)量的相干點(diǎn)目標(biāo)[6]。光譜相干系數(shù)可表征光譜多樣性的大小，聚焦良好的相干點(diǎn)目標(biāo)具有與分布式目標(biāo)不同的光譜特性，相似的相干點(diǎn)目標(biāo)具有相似的光譜特性，表現(xiàn)為低光譜多樣性，光譜相干系數(shù)較高。振幅穩(wěn)定性可表征觀測目標(biāo)對雷達(dá)波后向散射的穩(wěn)定程度，相干點(diǎn)目標(biāo)具有較高的散射穩(wěn)定性。本研究中時序振幅穩(wěn)定性閾值設(shè)為1.4，光譜相干系數(shù)閾值設(shè)為0.5，分別獲得2個候選相干點(diǎn)目標(biāo)點(diǎn)集，然后對二者進(jìn)行合并，剔除同名點(diǎn)，共獲取32 758個初始相干點(diǎn)。

相干點(diǎn)的差分干涉相位主要包含形變、殘余地形、軌道誤差、大氣延遲和噪聲等相位。IPTA以二維線性相位回歸模型為基礎(chǔ)，即在空間域?qū)ο喔牲c(diǎn)與局部參考點(diǎn)的差分干涉相位進(jìn)行二次差分，以削弱大氣延遲和非線性形變的影響；在時間域?qū)Σ罘指缮嫦辔贿M(jìn)行回歸處理來分析相位變化過程[7]。相干點(diǎn)目標(biāo)的二次差分相位可表示為：

(1)

式中，φdiff,i-j為相干點(diǎn)目標(biāo)間的差分相位差，Δv、Δh分別為建立連接關(guān)系的點(diǎn)目標(biāo)的線性形變速率增量和高程改正值增量，BT、B⊥分別為差分干涉圖的時間基線和空間垂直基線，λ為雷達(dá)中心波長，R為雷達(dá)至地面目標(biāo)的斜距，θ為雷達(dá)入射角，φres,i-j為殘余相位，由軌道誤差相位φorb_res,i-j、大氣延遲相位φatm,i-j和噪聲相位φnoi,i-j組成。

在IPTA中引入殘余相位標(biāo)準(zhǔn)差來判斷回歸分析的最優(yōu)解，并采用逐步回歸、迭代改進(jìn)的方式不斷計算修正值，更新線性形變相位、高程相位和殘余相位，計算所有相干點(diǎn)目標(biāo)間的線性形變速率增量和高程誤差改正值增量，選定一個相對穩(wěn)定點(diǎn)作為參考點(diǎn)，從而獲取每個相干點(diǎn)相對于參考點(diǎn)的線性形變速率和高程改正值。根據(jù)IPTA官方說明書，殘余相位標(biāo)準(zhǔn)差小于1.2的相干點(diǎn)可滿足精度要求。由于本文研究區(qū)地形起伏較大，且無先驗(yàn)信息，因此將第一次回歸分析的高程改正值閾值設(shè)為±2 000 m，形變速率閾值設(shè)為±150 mm/a。較大的搜索范圍會解算出較大的形變速率和地形誤差，后續(xù)逐步改進(jìn)閾值，縮小范圍。在最終的回歸分析中，高程改正值閾值設(shè)為±500 m，形變速率閾值設(shè)為±80 mm/a。最終輸出的相干點(diǎn)目標(biāo)相位標(biāo)準(zhǔn)差均低于1.2，符合實(shí)驗(yàn)精度要求。

從差分干涉相位中扣除線性形變和高程誤差相位可得到殘余相位。其中，殘余軌道誤差在InSAR中表現(xiàn)為基線誤差，可通過基線精化改正[8]；噪聲在空間維為高頻信號，在空間域?qū)堄嘞辔贿M(jìn)行低通濾波可去除；大氣延遲為時間維的高頻信息，非線性形變?yōu)闀r間維的低通成分，在時間域進(jìn)行低通濾波可得到非線性形變。最后將線性形變與非線性形變進(jìn)行組合，計算得到每個相干點(diǎn)的形變時間序列。IPTA數(shù)據(jù)處理流程如圖1所示。

圖1 數(shù)據(jù)處理流程Fig.1 Flow chart of data processing

2 研究區(qū)概況與實(shí)驗(yàn)數(shù)據(jù)

2.1 研究區(qū)概況

大光包滑坡為2008年汶川地震觸發(fā)的最大規(guī)模滑坡，汶川地震后，該地區(qū)海拔由3 047 m降至2 964 m，其長達(dá)1.8 km的剪滑光面和高達(dá)690 m的滑坡堰塞壩引起國內(nèi)外學(xué)者的廣泛關(guān)注[9-11]。大光包滑坡位于四川省綿陽市安州區(qū)高川鄉(xiāng)，中心坐標(biāo)為104°07′0.2″E、31°38′19.8″N，處于汶川地震發(fā)震斷層上盤，距地震中心85 km，與地震地表破裂帶相距4.5 km[9]，其地理位置如圖2所示。

圖2 大光包滑坡地理位置Fig.2 Geographical location of Daguangbao landslide

大光包滑坡地質(zhì)結(jié)構(gòu)復(fù)雜，根據(jù)滑坡要素、滑坡堆積體形態(tài)、長度和延伸方向、巖性分布，可將滑坡體分為12個區(qū)域(圖3)。Ⅰ為側(cè)向前沖式碎屑流區(qū)，表層物質(zhì)易隨白果林溝沖刷流動形成堰塞湖；Ⅱ?yàn)轱w躍前沖式碎屑流區(qū)，表層物質(zhì)向川林溝流動形成堰塞湖；Ⅲ為左側(cè)滑坡擾動體，巖體破碎嚴(yán)重，呈碎塊狀；Ⅳ為主滑坡體堆積區(qū)，由大光包頂拋射降落于該處，主要為泥盆系沙窩組白云巖；Ⅴ為溝道碎屑流區(qū)，多為黑色灰?guī)r碎塊，沿黃洞子溝向下流動堆積；Ⅵ為滑坡后緣二次滑塌體，在滑動時與主滑坡體分開后停留于此；Ⅶ為坡面碎屑流區(qū)，巖性主要為灰黑色灰?guī)r；Ⅷ為雨水沖刷滑坡表層物質(zhì)形成的溝口沖積扇；Ⅸ為白云巖和板巖構(gòu)成的滑坡左側(cè)斷壁；Ⅹ為滑坡后緣斷壁倒石錐，巖體松散破碎；Ⅺ為滑坡右側(cè)主滑動面，滑床為震旦系燈影組；Ⅻ為滑坡后緣斷壁，出露大光包粉砂巖至白云巖的地層剖面[10-11]。

圖3 大光包滑坡分區(qū)Fig.3 Zoning of Daguangbao landslide

研究區(qū)地勢西高東低，地形起伏大，地貌以高山為主，具有逆沖-推覆-滑脫-走滑的構(gòu)造特點(diǎn)[9]。由于研究區(qū)常年降雨，地質(zhì)環(huán)境脆弱，該區(qū)在接下來數(shù)十年內(nèi)將會處于不穩(wěn)定狀態(tài)，再次發(fā)生山體滑坡的風(fēng)險高于地震前[12-13]，對其進(jìn)行形變監(jiān)測具有重要意義。

2.2 實(shí)驗(yàn)數(shù)據(jù)

本文采用覆蓋實(shí)驗(yàn)區(qū)的寬幅模式(IW)升軌Sentinel-1A影像，幅寬為250 km，空間分辨率為5 m×20 m，極化方式為垂直同向(VV)極化，入射角為33.83°。研究時間跨度為2016-01～2018-12，影像獲取時間間隔為12 d，共計65幅。干涉組合采用單一主影像模式，空間基線閾值設(shè)定為100 m，圖4為65幅影像的成像時間與空間基線分布。本文選用的DEM數(shù)據(jù)為2000年美國宇航局主導(dǎo)測量的SRTM-3 DEM數(shù)據(jù)，分辨率為90 m，用于地理編碼和差分處理以去除地形相位。

圖4 時空基線分布Fig.4 Distribution of spatiotemporal baselines

3 監(jiān)測結(jié)果與分析

圖5為采用IPTA方法得到的大光包滑坡地區(qū)雷達(dá)視線向年形變速率。從圖中可以看出，監(jiān)測期間大光包滑坡處于形變狀態(tài)，滑坡體左部較為穩(wěn)定，形變速率為-15～2 mm/a。飛躍前沖式碎屑流區(qū)(Ⅱ區(qū))和主滑坡體堆積區(qū)(Ⅳ區(qū))呈現(xiàn)明顯漏斗狀形變，表面物質(zhì)松散不穩(wěn)定，受川林溝、黃洞子溝和雨水沖刷影響，形變量大于其他區(qū)域，雷達(dá)視線向最大形變速率達(dá)-50.590 mm/a，最大累積形變量達(dá)-170.726 mm。

圖5 雷達(dá)視線向形變速率Fig.5 Line-of-sight deformation rate

以72 d為時間間隔輸出2016-01-13～2018-12-28滑坡的形變時間序列，結(jié)果如圖6所示。由圖可知，監(jiān)測時段內(nèi)大光包滑坡右側(cè)形變漏斗逐漸發(fā)育，監(jiān)測時間到684 d時，形變漏斗已初具規(guī)模，雷達(dá)視線向累積形變量為-111.543 mm。隨著時間推移，形變漏斗的范圍逐漸擴(kuò)大，到2018年年底，形變漏斗長度達(dá)767 m，寬度達(dá)687 m，雷達(dá)視線向累積形變量為-170.726 mm。

圖6 大光包滑坡形變時間序列Fig.6 Deformation time series of Daguangbao landslide

在大光包滑坡形變漏斗中心和邊緣區(qū)域選擇6個特征點(diǎn)進(jìn)行時間序列形變分析，特征點(diǎn)位置如圖7所示，形變結(jié)果如圖8所示?？梢钥闯?，6個特征點(diǎn)均呈線性滑動趨勢，但形變速率不同。A、B、C點(diǎn)位于形變量較大的漏斗中心，所在位置的坡度分別為49.73°、53.31°和51.49°，整體形變量和形變速率逐年增大，形變速率約為-48.387～-47.361 mm/a，累積形變量約為-141～-138 mm；D、E、F點(diǎn)位于形變量較小的漏斗邊緣區(qū)域，所在位置的坡度分別為33.67°、25.80°和31.06°，形變發(fā)育較為穩(wěn)定，形變速率約為-17.788～-9.219 mm/a，累積形變量約為-57.432～-31.117 mm。

圖7 特征點(diǎn)位置Fig.7 Location of feature points

本文引入大光包地區(qū)2016-01-01～2018-12-28的日均降水量數(shù)據(jù)，結(jié)合特征點(diǎn)時序形變進(jìn)行比較，并統(tǒng)計雷達(dá)重訪周期內(nèi)暴雨發(fā)生前后特征點(diǎn)的形變量，結(jié)果如表1所示。從圖8可以看出，大光包地區(qū)的降雨量具有明顯的季節(jié)性特征，每年5～10月雨水充足，8月常有強(qiáng)降雨發(fā)生，滑坡上松散物質(zhì)受雨水沖刷影響，碎石塊在滑坡表面滑動[14]，引起累積形變量增大。2016-08-22降雨量為118.815 mm，2016-08-16～09-09特征點(diǎn)的最大形變增量為-10.829 mm；2017-08-22降雨量為80.400 mm，2017-08-11～23特征點(diǎn)的最大形變增量為-5.877 mm；2018-06-25降雨量為71.849 mm，2018-06-19～07-01特征點(diǎn)的最大形變增量為4.858 mm；2018-07-10降雨量為 60.782 mm，2018-07-01～25特征點(diǎn)的最大形變增量為2.162 mm。2018年強(qiáng)降雨發(fā)生后，滑坡上松散物質(zhì)受雨水影響，碎石塊在滑坡表面形成新堆積體并伴有滑坡體吸水膨脹抬升[14]，引起特征點(diǎn)形變量減小。綜上表明，降水是影響滑坡穩(wěn)定性的主要因素之一。

表1 特征點(diǎn)形變量統(tǒng)計

為進(jìn)一步分析形變漏斗特征，以D點(diǎn)為起點(diǎn)，途經(jīng)A點(diǎn)，F(xiàn)點(diǎn)為終點(diǎn)，作剖面線進(jìn)行研究。從上到下提取該剖面線共74個點(diǎn)的累積形變量，繪制以144 d為間隔的剖面線形變量時序圖(圖9)。

圖9 剖面線形變量時序Fig.9 Time series deformation of the profile

研究區(qū)地形起伏較大，最大高差達(dá)245 m，在監(jiān)測期間內(nèi)該坡面線累積形變量逐年增大，滑坡頂部形變量較小，形變集中在滑坡中部碎屑流區(qū)。該區(qū)域坡向?yàn)闁|南向，坡度為45°～50°，松散堆積體易受到降雨及黃洞子溝沖刷影響而產(chǎn)生較大形變，監(jiān)測時段內(nèi)最大形變量達(dá)-167 mm。

采用StaMPS-SBAS方法對Sentinel-1A影像數(shù)據(jù)進(jìn)行處理，時間基線設(shè)為50 d，空間基線設(shè)為100m，生成212個干涉對，剔除相干性較差的

干涉對后余下183個干涉對，共提取57 225個相干點(diǎn)，所提取的形變速率如圖10所示。對比分析IPTA方法和StaMPS-SBAS方法所得結(jié)果可知，兩者趨勢一致，飛躍前沖式碎屑流區(qū)(Ⅱ區(qū))和主滑坡體堆積區(qū)(Ⅳ區(qū))均呈現(xiàn)出形變漏斗，且形變數(shù)值一致，表明兩者結(jié)果可靠。

圖10 StaMPS-SBAS方法得到的雷達(dá)視線向形變速率Fig.10 Line-of-sight deformation rate obtained by StaMPS-SBAS method

4 結(jié) 語

本文利用IPTA方法對升軌Sentinel-1A影像數(shù)據(jù)進(jìn)行處理，獲取大光包滑坡2016～2018年的形變時間序列，并結(jié)合GPM日均降水量數(shù)據(jù)進(jìn)行特征點(diǎn)累積形變分析，得到以下結(jié)論：

1) 大光包地區(qū)在2016～2018年發(fā)生明顯形變，形變主要集中在飛躍前沖式碎屑流區(qū)和主滑坡體堆積區(qū)，呈線性滑動趨勢，雷達(dá)視線向形變速率在-50.590～19.175 mm/a之間，最大累積形變量達(dá)-170.726 mm。

2) 結(jié)合GPM日均降水量分析發(fā)現(xiàn)，大光包滑坡形變與降雨量變化具有一定相關(guān)性。研究區(qū)每年5～10月雨水充足，該時間段應(yīng)結(jié)合降水?dāng)?shù)據(jù)，加強(qiáng)滑坡監(jiān)測，作好相應(yīng)防范措施，防止滑坡泥石流等災(zāi)害發(fā)生。

3)IPTA方法和StaMPS-SBAS方法均探測出飛躍前沖式碎屑流區(qū)和主滑坡體堆積區(qū)出現(xiàn)形變漏斗，形變量級一致，且形變整體趨勢相同，證明了2種方法的有效性。

4) 作為PS-InSAR方法的發(fā)展，IPTA方法不僅可以減少時空失相關(guān)和大氣延遲的影響，還可以通過矢量化數(shù)據(jù)節(jié)省計算空間，提高數(shù)據(jù)處理速度和效率，迭代求解得到的結(jié)果也更加可靠，在大型滑坡形變監(jiān)測中具有較好的可靠性和廣闊的應(yīng)用前景。