基于SBAS-InSAR技術(shù)的金沙江流域沃達(dá)村巨型老滑坡形變分析*

馮文凱 頓佳偉 易小宇 張國強(qiáng)

(地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室(成都理工大學(xué))，成都 610059，中國)

0 引 言

隨著人類活動(dòng)深度和廣度的日益加深和拓展，我國西部地區(qū)大型突發(fā)地質(zhì)災(zāi)害日益頻發(fā)。這些災(zāi)害表現(xiàn)出隱蔽性強(qiáng)、危害性高，具有顯著的突發(fā)性，嚴(yán)重影響了社會經(jīng)濟(jì)發(fā)展和人民安居樂業(yè)，給地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)帶來了巨大挑戰(zhàn)。

2017年6月24日清晨6時(shí)許，四川省茂縣疊溪鎮(zhèn)新磨村新村組后山約450×104￣￣m3山體突發(fā)順層高位滑動(dòng)，造成80余人失蹤或遇難，引起國內(nèi)外廣泛關(guān)注(許強(qiáng)等， 2017)。災(zāi)害發(fā)生后，意大利Tre Altamira (2017)通過InSAR(Interferometric Synthetic Aperture Radar，合成孔徑雷達(dá)干涉)技術(shù)對2014年以來45景Sentinel-1雷達(dá)衛(wèi)星影像進(jìn)行了分析，證實(shí)新磨村滑坡的顯著形變主要發(fā)生在滑前的數(shù)月之內(nèi)。2018年10月11日凌晨4時(shí)許，金沙江流域西藏昌都市江達(dá)縣波羅鄉(xiāng)白格村發(fā)生特大庫岸巖質(zhì)滑坡，總方量約3165×104￣￣m3，滑坡體迅速阻斷金沙江，回水造成多條道路中斷，川藏滇等地緊急轉(zhuǎn)移安置2.1萬余人(馮文凱等， 2019)。許強(qiáng)等(2018)通過分析白格滑坡災(zāi)前雷達(dá)衛(wèi)星數(shù)據(jù)，指出該滑坡災(zāi)害發(fā)生前一年，滑源區(qū)斜坡最大位移量約25 m。這些案例充分說明InSAR技術(shù)對處于高位、交通困難、人跡罕至的隱蔽性地質(zhì)災(zāi)害的識別和監(jiān)測研究具有顯著優(yōu)勢。

InSAR技術(shù)是近20 a來興起的一種新型對地觀測技術(shù)，具有全天候工作、覆蓋范圍廣、探測精度高、非接觸和綜合成本低等優(yōu)點(diǎn)，已被廣泛地應(yīng)用于滑坡(李凌婧等，2017)、地面沉降(張靜等，2018)、地裂縫(趙超英等，2011)等地質(zhì)災(zāi)害的形變監(jiān)測。InSAR技術(shù)中的小基線(Small Baseline Subsets， SBAS)技術(shù)具有獲取微小形變信息和長時(shí)間序列緩慢地表變形的優(yōu)勢(莫玉娟等， 2018)，在復(fù)雜山區(qū)地質(zhì)災(zāi)害預(yù)警和監(jiān)測領(lǐng)域有較為廣闊的應(yīng)用前景。在基于SBAS-InSAR技術(shù)對滑坡早期識別與形變監(jiān)測的研究領(lǐng)域中，很多學(xué)者取得了相應(yīng)的實(shí)質(zhì)性成果。余睿(2014)采用SBAS-InSAR技術(shù)對甘肅省西和縣慢速型滑坡進(jìn)行全局和局部的監(jiān)測，揭示了影響西和縣滑坡體移動(dòng)的主導(dǎo)因素； 劉筱怡等(2017)以鮮河水?dāng)嗔褞?nèi)蠕滑型滑坡為研究對象，采用SBAS-InSAR技術(shù)獲取蠕滑型滑坡時(shí)間序列形變特征，分析其發(fā)展變化規(guī)律，對鮮水河斷裂帶以及類似構(gòu)造活動(dòng)地區(qū)的地質(zhì)災(zāi)害防治起到了重要作用； 聶兵其等(2018)基于SBAS-InSAR技術(shù)對丹巴縣城進(jìn)行滑坡形變探測及隱患識別，驗(yàn)證了InSAR監(jiān)測結(jié)果與現(xiàn)場調(diào)查的一致性?？梢哉f，SBAS-InSAR技術(shù)為地質(zhì)災(zāi)害隱患點(diǎn)的普查和分析研究提供了重要支撐。在充分利用星載雷達(dá)對地觀測技術(shù)的前提下，開展針對性的工程地質(zhì)調(diào)查和復(fù)核工作，一方面有利于正確認(rèn)識InSAR解譯數(shù)據(jù)的規(guī)律性，另一方面也有利于災(zāi)害預(yù)警研判的準(zhǔn)確性，從而大大提升我國科學(xué)防災(zāi)減災(zāi)的技術(shù)實(shí)力。

沃達(dá)村巨型老滑坡是區(qū)域地質(zhì)災(zāi)害調(diào)查隱患數(shù)據(jù)庫的隱患點(diǎn)，位于江達(dá)縣巖比鄉(xiāng)沃達(dá)村、金沙江右岸，距離白格滑坡直線距離42 km。本文在搜集沃達(dá)村老滑坡區(qū)域升軌哨兵一號(Sentinel-1A)2017年3月30日至2019年9月28日的39景數(shù)據(jù)的基礎(chǔ)上，采用SBAS-InSAR技術(shù)分析了老滑坡堆積體的形變特征，并進(jìn)一步針對性地開展了現(xiàn)場工程地質(zhì)調(diào)查和復(fù)核工作，在此基礎(chǔ)上分析了SBAS-InSRAR技術(shù)在復(fù)雜山區(qū)地質(zhì)災(zāi)害預(yù)警和監(jiān)測領(lǐng)域的適用性，為類似老滑坡監(jiān)測預(yù)警提供了新的思路與借鑒。

1 沃達(dá)村老滑坡概況

沃達(dá)村老滑坡位于西藏自治區(qū)昌都市江達(dá)縣巖比鄉(xiāng)沃達(dá)村，坐標(biāo)東經(jīng)E93°26′06″、北緯N31°26′06″(圖 1)。前期區(qū)域地質(zhì)災(zāi)害調(diào)查表明，該地段地質(zhì)歷史時(shí)期曾發(fā)生滑坡災(zāi)害，其滑坡堆積體在降水、地震和水位波動(dòng)情況下易發(fā)生復(fù)活，應(yīng)對其重點(diǎn)監(jiān)測研究。

圖 1 滑坡地理位置圖Fig. 1 Geographic location map of the landslide

圖 2 綜合光學(xué)遙感解譯圖Fig. 2 Integrated optical remote sensing interpretation map

從綜合光學(xué)遙感解譯圖(圖 2)分析，老滑坡體縱長約1.9 km，橫寬約1.7 km，主滑方向?yàn)镹E30°； 老滑坡堆積體局部復(fù)活，復(fù)活部分縱向長約0.88 km，橫寬約1.08 km，面積約0.925 km2?；麦w平面形態(tài)為半圓狀，前緣呈由緩變陡的陡崖形態(tài)，緩部與江面高差最大約240 m，前端臨空面陡傾，坡腳即為金沙江河谷，高程約2950 m，中后部相對平緩。據(jù)前期調(diào)查顯示，該老滑坡上部覆蓋層為第四系滑坡堆積物，含碎石、礫石黏性土，厚約20～30 m，下伏基巖為變質(zhì)頁巖，巖層產(chǎn)狀約為240°∠80°。

圖 3 數(shù)據(jù)覆蓋范圍圖Fig. 3 Data coverage chart

2 SAR數(shù)據(jù)獲取及處理方法

2.1 數(shù)據(jù)選取及覆蓋概況

本文使用的單視復(fù)數(shù)SAR影像數(shù)據(jù)來自歐空局(European Space Agency)2014年4月3日發(fā)射的對地觀測衛(wèi)星Sentinel 1A，數(shù)據(jù)覆蓋區(qū)域如圖 3所示。波段為C波段，所用影像為合成孔徑雷達(dá)99軌道數(shù)據(jù)，時(shí)間跨度從2017年3月30日至2019年9月28日， 24 d一景，共39景，軌道方向?yàn)樯墸上裰行娜肷浣菫?6.2°，成像模式為IW寬幅模式，地面分辨率為5 m×20 m，幅寬為250 km×250 km，極化方式為VV，方位向和距離向采樣間距分別為2.329 m、13.956 m(表 1)。本研究采用了日本JAXA提供的30 m分辨率ALOS World 3D(AW3D30)DEM作為外部數(shù)據(jù)來消除地形相位的影響。

表 1 沃達(dá)村滑坡Sentinel 1A數(shù)據(jù)列表Table 1 Data list of sentinel 1 in Dawo Village landslide

2.2 技術(shù)原理與形變速率獲取

2.2.1 SBAS-InSAR技術(shù)原理

為減少時(shí)間和空間去相干，Berardino et al.(2002)在Casu(2000)研究的基礎(chǔ)上，提出了短基線差分干涉法。該方法的基本原理是利用具有較短時(shí)-空基線的影像對產(chǎn)生干涉圖提高相干性，再進(jìn)行解纏干涉相位。依據(jù)相干像元相位和觀測時(shí)間的關(guān)系，采用奇異值分解方法將多個(gè)小基線集合數(shù)據(jù)聯(lián)合起來求解(許軍強(qiáng)等， 2019)，有效地解決了各數(shù)據(jù)集之間空間基線過長造成的時(shí)間不連續(xù)問題，提高了監(jiān)測的時(shí)間分辨率，從而求得影像序列間地表形變速率的最小范數(shù)最小二乘解(何秀鳳等， 2012)。

其基本原理：假設(shè)有N+1景單視復(fù)數(shù)影像，其成像時(shí)間為(t0，t1，….，tn)，設(shè)某一差分解纏圖為j，某像元(x，r)對應(yīng)的解纏相位可表示為式(1)：

δφj(x，r)=φ(tB，x，r)-φ(tA，x，r)

(1)

(2)

由于在處理過程中產(chǎn)生M幅干涉圖，根據(jù)式(2)可得M個(gè)方程，以矩陣形式表示為式(3)：

δφ(x，r)=Aφ(x，r)

(3)

其中，A為M×N系數(shù)矩陣；φ(x，r)為(x，r)點(diǎn)在N個(gè)時(shí)刻對應(yīng)的未知形變相位構(gòu)成的矩陣，當(dāng)M≥N時(shí)，通過最小二乘法得出式(4)：

φ(x，r)=(ATA)-1ATδφ(x，r)

(4)

當(dāng)M

圖 4 影像時(shí)間基線Fig. 4 Image time baseline

圖 5 影像空間基線Fig. 5 Image spatial baseline

2.2.2 視向形變速率獲取

本文使用ENVI中SArscape模塊，用SBAS算法對Sentinel 1A數(shù)據(jù)進(jìn)行干涉處理，共生成109個(gè)干涉對，時(shí)間基線為75 d，臨界基線百分比為2%，各像對時(shí)間基線和空間基線的連接方式如圖 4和圖 5所示，圖中綠色的點(diǎn)代表 39景影像的成像時(shí)間，黃色點(diǎn)代表作為參考主影像的成像時(shí)間。對合成的干涉圖做第1次反演、第2次反演、地理編碼以及柵矢轉(zhuǎn)換，生成地表沿雷達(dá)視向的平均形變速率(mm·a-1)，最終得到時(shí)間序列上的形變結(jié)果。具體處理流程如圖 6所示。

圖 6 SBAS處理流程Fig. 6 SBAS-InSAR processing flow

圖 7 雷達(dá)視線方向和坡度方向幾何示意圖 (Cascini et al., 2010)Fig. 7 Geometric sketches of radar line- of -sight direction and slope direction(Cascini et al.， 2010)

2.2.3 二維形變轉(zhuǎn)換

由于滑坡多沿斜坡面進(jìn)行滑動(dòng)，雷達(dá)視線方向的形變信息無法準(zhǔn)確反映斜坡面的真實(shí)形變情況(戴可人等， 2019)，考慮雷達(dá)視線方向、坡度方向和垂直沉降方向的幾何關(guān)系(圖 7)，假設(shè)運(yùn)動(dòng)沿單位矢量指定的方向發(fā)生，采用式(5)、式(6)(Cascini et al.， 2010)和式(8)(Colesanti et al.， 2006)將視線方向的形變速率轉(zhuǎn)化為坡度方向和垂直方向的形變速率。

(5)

VSlope=VLos/cosβ

(6)

式(6)中，cosβ可用式(7)(Cascini et al.， 2010)計(jì)算；

(7)

(8)

3 形變結(jié)果分析

結(jié)合遙感影像分析，該滑坡中部和后緣部位表面植被覆蓋密度不大，相干點(diǎn)較多，干涉效果較好，有較明顯的速率積累，滑坡部分地方植被覆蓋茂密，造成影像的失相干，速率點(diǎn)較為稀少。

形變速率為負(fù)值表示形變點(diǎn)向遠(yuǎn)離衛(wèi)星傳感器方向運(yùn)動(dòng)，正值表示形變點(diǎn)向靠近衛(wèi)星傳感器方向運(yùn)動(dòng)。由InSAR視向形變分析結(jié)果可見，圖 8中形變值均為負(fù)值，表明滑坡整體向著金沙江方向向下滑動(dòng)，符合滑坡體的運(yùn)動(dòng)規(guī)律。從圖 8中可以看出該滑坡形變較強(qiáng)部位位于滑坡中部和中前緣(圖 8紫色虛線圈內(nèi))，其中紫色虛線圈內(nèi)的紅色區(qū)域形變量最大，形變速率≥80 mm·a-1。

圖 8 2017年3月～2019年9月形變速率圖(LOS方向)Fig. 8 Deformation rate map from March 2017 to September 2019(LOS Direction)

沃達(dá)村滑坡沿坡向年平均速率(VSlope)和垂直方向年平均速率(Vu)分別如圖 9、圖 10所示，解譯出的強(qiáng)變形區(qū)與LOS方向大致一致，但在形變速率量值上存在較大差異，沿坡向形變速率最大。

圖 9 2017年3月～2019年9月形變速率圖(Slope方向)Fig. 9 Deformation rate map from March 2017 to September 2019(Slope Direction)

圖 10 2017年3月～2019年9月形變速率圖(Vertical方向)Fig. 10 Deformation rate map from March 2017 to September 2019(Vertical Direction)

4 實(shí)地工程地質(zhì)調(diào)查與復(fù)核

為復(fù)核解譯數(shù)據(jù)，開展了實(shí)地工程地質(zhì)調(diào)查與復(fù)核工作。經(jīng)實(shí)地踏勘，滑坡體前緣地形陡峭； 中前部坡度較大，樹木生長茂盛； 中部為緩坡平臺，主要為耕地，樹木零星生長； 后部為古滑坡滑源區(qū)(圖 11a)。受沖溝下切侵蝕及前緣陡傾臨空面綜合作用，滑坡前緣發(fā)育多處向沖溝、臨空面方向的崩滑變形體(B01～B04)，崩滑變形體表部裂縫橫張、錯(cuò)坎發(fā)育，裂縫張開2～10 cm，最大錯(cuò)距1～2 m，前端陡崖可見明顯基覆界面(圖 11b)。滑坡復(fù)活區(qū)中部發(fā)育多處小型滑動(dòng)變形體(H01～H03)，坡體發(fā)育橫切梯級陡坎，高度幾十厘米到幾米不等，延伸幾十米到數(shù)百米不等，陡坎上部多耕作，有灌木叢生，表明陡坎形成時(shí)間較長(圖 11c)； 可見明顯馬刀樹，結(jié)合樹齡推測穩(wěn)定時(shí)間至少10 a以上(圖 11d)。

圖 11 沃達(dá)村老滑坡照片F(xiàn)ig. 11 Photographs of the old landslide in Woda Village

4.1 滑坡整體變形與形變速率對比分析

圖 12 InSAR時(shí)序形變與降雨響應(yīng)Fig. 12 InSAR time series deformation and rainfall response

復(fù)活區(qū)主要位于老滑坡體中下部，強(qiáng)烈變形區(qū)(Ⅰ現(xiàn)場)后緣位于老滑坡后緣陡緩交界處，兩側(cè)邊界依地形延伸至老滑坡邊界(圖 8)，其內(nèi)裂縫、下錯(cuò)臺坎、局部滑塌現(xiàn)象明顯。現(xiàn)場的變形跡象與SBAS-InSAR方法獲取的形變速率圖(圖 8)特征有著較高的一致性。在滑坡視向形變速率圖(圖 8)形變量較大的紅色區(qū)域取一點(diǎn)A作出該點(diǎn)歷史累積位移形變曲線，與搜集的2017年1月～2019年8月月累計(jì)降雨曲線對比分析(圖 12)：隨著2017-05～2017-06、2018-05～2018-08、2019-06～2019-07 3個(gè)不同時(shí)間段內(nèi)降雨量急劇增加，A點(diǎn)形變量明顯加速； 降雨量較小的時(shí)間段內(nèi)，A點(diǎn)形變量增加變緩，形變趨勢與降雨具有較好的響應(yīng)規(guī)律。據(jù)現(xiàn)場調(diào)查，滑坡所處金沙江水位波動(dòng)最高高程為2930～2950 m，而滑坡局部復(fù)活區(qū)位置位于3220～3580 m，遠(yuǎn)離金沙江水位波動(dòng)影響區(qū)。說明該滑坡中前部復(fù)活區(qū)主要受降雨影響，為老滑坡降雨復(fù)活型滑坡。

圖 13 2017年3～12月形變速率及剖面形變圖Fig. 13 Deformation rate and profile deformation map from March to December 2017

圖 14 2018年形變速率及剖面形變圖Fig. 14 Deformation rate and profile deformation map of 2018

圖 15 2019年1～9月形變速率及剖面形變圖Fig. 15 Deformation rate and profile deformation map from January to September 2019

4.2 滑坡局部變形與形變速率對比分析

現(xiàn)場復(fù)核階段對前緣臨空面發(fā)育的4處崩滑變形體(B01～B04)和滑坡中、下部3處滑動(dòng)變形體(H01～H03)進(jìn)行了細(xì)致調(diào)查。4處崩滑變形體均為表層發(fā)生的溜滑，規(guī)模較小，樹木歪斜的范圍位于沖溝溝口，表明為集中降雨引起的變形。3處滑動(dòng)變形體中， 2處(H02和H03)因?yàn)榍安颗R空面變形導(dǎo)致出現(xiàn)牽引式滑動(dòng)變形， 1處(H01)因降雨和人類工程活動(dòng)因素發(fā)生牽引式滑動(dòng)變形。下面主要分析滑動(dòng)變形體H01的情況。

H01滑動(dòng)變形體位于坡體的中后部，沃達(dá)村公路上側(cè)，滑動(dòng)方向與坡向一致(圖 16a)。據(jù)訪，H01滑動(dòng)變形體前緣為沃達(dá)村村道，村道擋土墻于2013年左右發(fā)生隆起破壞(圖 16b)； 2018年中旬，H01滑動(dòng)變形體發(fā)生了較為強(qiáng)烈的變形，導(dǎo)致?lián)鯄ι戏捷^高處田地出現(xiàn)多處裂縫及下錯(cuò)臺坎(圖 16c)。調(diào)查發(fā)現(xiàn)，其下部緊鄰的村民房屋完好無破壞，表明其僅僅是近年來表層巖土體在降雨等作用下發(fā)生的局部牽引式滑動(dòng)變形。

為細(xì)致分析H01滑動(dòng)變形體的滑動(dòng)性質(zhì)與雷達(dá)形變速率信號之間的相關(guān)性，沿H01滑動(dòng)變形體主滑方向布置剖面線，并依次取1#～4#監(jiān)測點(diǎn)，利用圖 9繪制其3年內(nèi)滑坡歷史累積位移形變曲線(圖 16d)。如圖所示， 4#監(jiān)測點(diǎn)歷史累積形變量較大，增加的速率也較大； 1#～3#監(jiān)測點(diǎn)歷史累積形變量增加幅度基本一致。前緣的形變速率遠(yuǎn)大于中部和后部的形變速率，說明H01滑動(dòng)變形體屬于牽引式滑動(dòng)變形。

根據(jù)對滑塌體后緣下錯(cuò)臺坎發(fā)育處沿坡向以及垂向形變速率點(diǎn)的分析：沿坡向的形變速率(圖 9、圖 16a)顏色各異，形變速率差異性很大，形變速率140 mm·a-1、160 mm·a-1、180 mm·a-1均有分布，與裂縫形成特征吻合； 沿垂向形變速率分析(圖 10)，在后緣位置黃色形變速率(-80～-60 mm·a-1)附近出現(xiàn)多級的橙色形變速率(-100～-80 mm·a-1)，說明該區(qū)域發(fā)生多級下錯(cuò)沉降現(xiàn)象。兩種現(xiàn)象印證了滑坡強(qiáng)變形區(qū)歷史出現(xiàn)多級下錯(cuò)臺坎，表現(xiàn)出蠕動(dòng)變形，與地質(zhì)調(diào)查的現(xiàn)象基本相符。

5 結(jié) 論

(1)通過對視向、坡向、垂向方向的形變速率分區(qū)以及形變點(diǎn)之間的差異性分析，初步將沃達(dá)村滑坡分為強(qiáng)烈形變區(qū)(Ⅰ雷達(dá))和均勻形變區(qū)(Ⅱ雷達(dá))，與現(xiàn)場判定的強(qiáng)烈變形區(qū)(Ⅰ現(xiàn)場)和弱變形區(qū)(Ⅱ現(xiàn)場)基本吻合。

圖 16 H01變形體SBAS-InSAR監(jiān)測與野外現(xiàn)象對比Fig. 16 Comparison of SBAS-InSAR monitoring data of H01 landslide with field review

(3)基于坡向形變速率圖、垂向形變速率圖中形變速率值的分布以及量值的差異性特征，初步判定了局部滑塌的性質(zhì)以及多級臺坎發(fā)育的位置，與野外地質(zhì)調(diào)查吻合性較好。

(4)通過SBAS-InSAR技術(shù)分析以及現(xiàn)場調(diào)查復(fù)核，發(fā)現(xiàn)老滑坡復(fù)活變形跡象與SBAS-InSAR技術(shù)解譯成果有著較好的一致性，在復(fù)雜山區(qū)地質(zhì)災(zāi)害預(yù)警和監(jiān)測領(lǐng)域有較為廣闊的應(yīng)用前景，為類似老滑坡的監(jiān)測預(yù)警提供新的思路與借鑒。

致 謝本文在完成過程中得到了中國地震局地殼應(yīng)力研究所、四川省國土空間生態(tài)修復(fù)與地質(zhì)災(zāi)害防治研究院、江達(dá)縣氣象局等單位的幫助。歐洲太空局為本研究提供的Sentinel-1A數(shù)據(jù)，日本JAXA提供的30 m分辨率“ALOS World 3D—30 m” DEM。在此一并表示感謝。