基于SBAS-InSAR技術的大型滑坡時序微小形變分析及其災變時間預報

馬闖, 繆海波, 楊犇, 朱隆奇, 安廣強, 楊冰穎

(安徽理工大學土木建筑學院, 淮南 232001)

在高海拔山區及峽谷地帶,因強降雨和/或地震誘發的大型滑坡往往具有隱蔽性強、突發性高、危害性大等特點,給當地人民生命安全和國民經濟發展造成了嚴重威脅,同時給該區域的地質災害防范預警、防治和風險管控帶來了極大的挑戰。例如,2017年四川茂縣新磨村滑坡導致80人遇難[1],2018年金沙江流域發生的白格特大巖質滑坡造成了金沙江堵塞,影響了2.1萬余人[2]。因此,針對此類滑坡的早期識別和災變時間預報,是避免或減少人員傷亡和財產損失,提高地質災害防災減災水平的重要途徑。

由于大型災難性滑坡的演化過程通常跨越數年、數十年甚至上百年,且在災變前均會產生明顯的地表形變跡象[3]。因此,針對大型滑坡的地表時序形變解析結果不僅可以描述整個孕育期間的演化規律,而且對其運動學機制解析和災變前兆信息檢測(如滑坡破壞前的突然加速)等至關重要。目前,合成孔徑雷達干涉測量技術(interferometric synthetic aperture radar,InSAR)因其非接觸、超遠距離、高探測精度和全天候等特點,可快捷獲取地表(mm級別)變形跡象,在發生大規模滑移(m級別)前給予災變信號;該技術已在地表沉降、地裂縫發展、滑坡等地質災害的監測方面得到了廣泛運用[4],尤其是在交通困難的高海拔山區和深切峽谷地帶的隱蔽性古滑坡解譯與復活識別、新生型滑坡災變前兆信息檢測中已展現出獨特的優勢。

InSAR技術中的小基線算法(small baseline subsets,SBAS)能獲得更為精準的數字高程信息與細微的地表形變信息,近年來在復雜地貌形態的山區地表侵蝕、形變檢測及地質災害早期識別和監測預警中得到了深入發展[5]。許東麗等[6]以西藏龐村滑坡為例,獲取75景Sentinel-1A降軌數據,利用SBAS-InSAR技術對該滑坡的變形區域、形變速率進行劃分,并與地面監測結果進行比對,驗證了SBAS-InSAR技術在特大型滑坡監測中的可行性。高二濤等[7]利用SBAS-InSAR技術獲取蒲縣地區疑似滑坡體位置、范圍及滑坡 速率,在此基礎上并對兩種監測結果進行對比,深入分析滑坡隱患的可能性成因,為該地區滑坡誘因 提供科學的依據,主動規避自然地質災害風險,對 未來滑坡災害的監測與預警提供思路與方法。周呂等[8]基于SBAS-InSAR 技術研究了武漢中心城區的地表沉降特征,分析并初步建立了該地區地表沉降與城市建設、降雨量和長江水位變化等因素的關系。樂穎等[9]通過 SBAS-InSAR 技術完善研究礦區形變中出現的空缺值,從而更為準確地估算礦區地表形變結果,全面、有效建立起礦區的監測預警體系。上述研究成果表明,SBAS-InSAR技術憑借其在獲取持續的地表微小變形方面的優勢,能滿足隱患區域識別、城市地表形變監測及礦區災害預警等需求,因此具有廣闊的應用前景。

目前,準確地進行滑坡災變時間預報仍然是懸而未決的世界性難題。滑坡時間預報的起步是以齋藤蠕變模型的提出為標志[10],并由Hayashi等[11]進一步發展。Fukuzono[12]通過大量滑坡破壞過程的物理模型試驗,推導出在大多數情況下,臨近破壞階段的速度倒數與滑坡發生時間呈線性關系,Voight[13]也提出了相同的滑坡時間預報模型,此后Fukuzono-Voight速度倒數模型在滑坡、地震、火山爆發等預報中得到了廣泛的運用。然而,Fukuzono-Voight速度倒數模型預報結果的精確性需要持續位移時序數據的支撐[14-15]。對于隱蔽性強、突發性高的大型災難性滑坡,往往由于缺乏持續的形變時序數據使其預警預報變得異常困難,而SBAS-InSAR作為獲取持續且細微的地表形變信息的有效手段,其在大型突發性滑坡災變時間預報中的應用目前尚不廣泛。

因此,現以西藏昌都市江達縣波羅鄉白格滑坡為例,首先搜集其災變前的升軌哨兵一號(Sentinel-1A)2017年3月18 日至2018年10月3日的23景雷達影像,然后利用短基線雷達干涉測量技術(SBAS-InSAR)獲取該滑坡在此時間段內的地表垂直方向時序形變數據,進而分析該滑坡的形變速率、形變分布規律以及典型觀測點的累積位移變化趨勢,在此基礎上利用Fukuzono-Voight速度倒數模型對其災變時間進行預報;最后討論SBAS-InSAR技術在大型災難性滑坡的早期形變檢測和災變時間預報方面的應用前景,以期為今后類似滑坡的預警預報提供有價值的參考。

1 滑坡工程地質概況

白格滑坡(98°42′17.98′′E,31°4′56.41′′N)位于西藏自治區江達縣波羅鄉白格村金沙江右岸,分別于2018年10月11日和11月3日先后發生兩次高位滑動(其中第2次滑動為第1次滑動后形成的后緣強變形區),堵塞金沙江,形成堰塞湖。滑坡區屬于典型的構造侵蝕高山峽谷地貌。滑坡體所在斜坡的坡向為80°～110°,平均坡度為50°～65°。根據前人調查結果[2,16],整體滑動的主滑動方向為80°～105°,滑坡體長約1 413 m,平均寬度約560 m,平均厚度約40 m,最厚處達57 m,面積約為7.9×105m2(其中滑源區面積約為5.0×105m2),體積約為3 165×104m3。滑坡后緣為一山脊,高程3 720 m,前緣為金沙江,平水位高程2 880 m,垂直落差約840 m,前緣剪出口高程約為3 000 m。滑坡破壞后形成長約1 100 m,寬約500 m,面積約33×104m2的堰塞壩,阻斷金沙江。白格滑坡的地理位置、滑坡破壞后的照片如圖1和圖2所示。

圖2 滑坡全貌Fig.2 Landslide panorama

區域地質背景顯示,白格滑坡所在區域處于金沙江構造帶上,受多期構造運動,滑坡巖體結構破碎明顯,且經多期變質作用與強烈風化后,碎裂巖體遇水易崩解軟化。根據文獻[2],滑坡區出露的地層為元古界熊松群(Ptxn2)深灰色斜長片麻巖和華力西期蛇紋巖(φw4),其中蛇紋巖帶主要位于滑源區,綠泥石化嚴重。滑坡堆積物主要為碎塊石及其破碎物,碎屑成分較為復雜[16-17]。

2 SAR數據獲取及處理方法

2.1 數據來源

使用寬幅干涉(interferometric wide swath,IW)模式、單視復數(single look complex,SLC)的SAR影像,數據來自歐洲航天局雷達觀測站點的哨兵1號(Sentinel-1)衛星。選取23幅C波段Sentinel-1A升軌影像。裁剪后的覆蓋區域如圖3所示。

圖3 哨兵一號升軌數據裁剪覆蓋范圍Fig.3 Data clipping coverage map of Sentinel-1A

遙感數據獲取的時間跨度為2017年3月18日—2018年10月3日,時間間隔為24 d或者36 d,累積時間基線及空間基線長度如表1所示。同時為了提高SBAS-InSAR的干涉精度并減少干涉相位圖不連續產生的誤差,本文選取每一景影像采集時間前后對應的23條Sentinel-1A精密軌道數據校正影像。由于研究區內地形起伏較大,為消除地形相位對形變相位的影響,本文研究選取全新的全球30 m分辨率的DEM數據,其來自地理空間數據云上的美國國家航空航天局(NASA)。并借助ArcGIS中鑲嵌工具進行鑲嵌處理,最后在將ENVI格式下對DEM進行轉換,作為計算地理相位分量的輔助數據。

表1 白格滑坡Sentinel-1A數據時間序列列表Table 1 Time series of Sentinel-1A data of the Baige landslide

2.2 技術原理及形變速率獲取方法

SBAS-InSAR技術的基本原理是:首先獲取N+1景 SAR影像(成像時間為t0,t1, …,tN),通過影像匹配,獲得M幅差分干涉像對,且滿足(N+1)/2 ≤M≤N(N+1)/2,然后利用 DEM數據去除地形相位差產生差分干涉圖,再根據多次不同閾值實驗選取合理的初始閾值對其進行相位解纏,最后選取研究區內穩定區域或者已知形變量的地面控制點(ground control point,GCP)對解纏之后的差分干涉圖進行校正。若某一差分干涉圖為j(j= 1, 2, …,M),則像元(x,r)的差分干涉相位δφ(x,r)在不考慮殘余的地形相位、大氣延遲相位、噪聲相位的情況下,表達式[18]為

δφj(x,r)=φB(x,r)-φA(x,r)

(1)

式(1)中:x為方位向坐標;r為距離向坐標;λ為雷達波長;tA為從影像時間;tB為主影像時間;φA(x,r)為從影像相位;φB(x,r)為主影像相位;d(tB,x,r)和d(tA,x,r)分別為相對于參考時間t0的雷達視線方向(LOS方向)累積形變量。

若主、從影像按時間序列排列,則差分干涉相位可寫成矩陣形式[18],即

Aφ=δφ

(2)

式(2)中:A為M×N階系數矩陣,當基線集中含有多個子集時,矩陣A秩虧;δφ為像素點(x,r)處的干涉相位。

鑒于SBAS技術是在配準過程中是利用多幅主影像進行配準,故采用奇異值分解(singular value decomposition,SVD)的方法可以解決矩陣A的秩虧問題,且基于最小范數準則獲得式(2)中φ的解。 若將差分干涉相位的求解轉化為相位變化速率v的求解問題,則式(2)寫成矩陣形式為

Bv=δφ

(3)

式(3)中:B為M×N階矩陣。對B進行奇異值分解,得到各時間段內的相位變化速率v,然后計算相位時間序列并進一步得到形變時間序列。

SBAS-InSAR技術的具體處理方法如下。

步驟1在 ENVI的 SARscape模塊中,利用 SBAS方法對裁剪后的Sentinel-1A進行了干涉濾波處理,共生成了73個干涉像對,其中時間基線與空間基線閾值分別設置為120 d和2%。像對的空間基線與時間基線之間的關系如圖4和圖5所示。

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

圖5 影像時間基線Fig.5 Image time baseline

步驟2對干涉相對進行干涉處理,過程包括:相干性生成、干涉去平、濾波(選取Goldstein濾波方法,增強濾波效果)和相位解纏(選取最小費用流方法)。

步驟3軌道精煉和重去平。考慮到研究區地形復雜,山脈較多,故借助Google Earth 選取穩定點(如穩定交通樞紐等構筑物)作為地面控制點(ground control point,GCP),以此來降低因為選點不當而帶來的誤差。

步驟4SBAS第一次的反演,其中包括根據第一次反演閾值進行形變速率及殘余地形估算,以及對第二步生成的解纏圖進行優化,以便后續的處理。

步驟5地理編碼以及柵格轉矢量,生成地表雷達視線(line of sight,LOS)方向的年平均形變速率。考慮到LOS方向的形變信息不能代表滑坡沿滑動方向的實際變形特征,且白格滑坡的實際滑動方向存在轉折[3],故本文研究參考文獻[19],利用垂直方向的形變信息來分析白格滑坡在災變之前的演化特征。研究區域各點不同時間的累積形變量最終由柵格轉矢量獲得。

3 白格滑坡災變前形變結果分析

3.1 垂直方向形變速率分析

滑坡災變前的形變速率是分析滑坡演化過程的直觀信息。白格滑坡區因其較低的植被覆蓋密度而表現出良好的干涉效果,因此能夠看到明顯的形變速率變化。本文研究繪制了2017年3月18日—2018年10月3日白格滑坡沿垂直方向的地表微小形變的年平均形變速率,如圖6所示,并約定負值表示滑坡運動遠離雷達(即滑坡向下滑動),正值表示滑坡運動與雷達相向(即滑坡體隆起)。

圖6 白格滑坡垂直方向地表形變速率圖Fig.6 Deformation rate of the Baige landslide in vertical direction

由圖6看出,白格滑坡災變前具有明顯的形變分區。強烈形變區位于滑坡的后緣(紫色區域),垂直方向形變速率約為-39～-67 mm/a。次強烈形變區(紅色區域)主要分布在滑坡后緣、右側中部和前緣非臨江段,最大形變速率可達-39 mm/a,其中前緣非臨江段的形變跡象指示了滑坡剪出口的位置,這與白格滑坡災變時表現出的高位剪出一致[2,16]。一般形變區(粉色和藍色區域)垂直方向的年平均形變速率約在-0.7～-22 mm/a,覆蓋了大部分滑坡體,表明白格滑坡已經表現出整體滑動的趨勢。基本穩定區(黃色區域,主體形變速率>0 mm/a)主要分布在滑坡體后緣邊界、左側中部及右前緣臨江段。根據白格滑坡災變后的現場調查及滑坡演化過程解析[16],該滑坡首先阻滑區脆性剪段,導致主滑區啟動,隨后牽引區失去支撐而整體高位高速剪出。基于前人研究成果,本文研究認為右前緣臨江段因其阻滑作用表現出基本穩定,滑坡左側中部主要是后緣向下滑移引起的擠壓隆起,后緣邊界處則因整體尚未滑動而未失去支撐導致其處于基本穩定狀態。

3.2 垂直方向累積形變分析

圖7為2017年3月18日—2018年10月3日白格滑坡沿垂直方向的累積時序形變圖。該滑坡在在災變之前的變形特征表現出較大的滑移量(m級),而SAR數據僅僅從其微小形變(mm級)維度去捕獲,但二者的變化的共性相同,都是滑坡變形破壞的特征。該時序圖以2017年3月18日影像為基準,即假設該時間點地表垂直方向的形變量為0 mm。根據圖7所示的不同時刻白格滑坡垂方向累積形變量將該滑坡在研究期間內的變形演化過程劃分為3個階段。

第1階段：2017年3月18日—6月22日,白格滑坡的地表垂直方向形變特征表現為滑坡前緣—中前部的向下滑移,且前緣累積滑移量持續增大,垂直方向的時序累積形變量最大達-75～-115 mm。這一結果表明,在此時間段內白格滑坡前緣有向金沙江滑移的趨勢。

第2階段：2017年8月9日—12月31日,白格滑坡的地表垂直方向形變特征表現為滑坡后緣滑移區的出現與持續發展,后緣垂直方向的累積形變量達-32～-50 mm,但前緣在垂直方向上的最大累積形變量則為-14～-50 mm。這一結果表明,在此時間段內,滑坡后緣開始出現持續的牽引滑移,前緣因后緣對中部的擠壓而表現出一定的隆起,此階段前緣起到明顯的阻滑作用。

第3階段：2018年1月24日—10月3日,白格滑坡的地表垂直方向形變特征表現為滑坡后緣的滑移急劇發展,且滑移牽引區不斷擴大,垂直方向的最大累積形變量已達-115 mm。與此同時,滑坡前緣也表現出持續向下滑移且有向中部擴展的趨勢,垂直方向的最大累積形變量達-50 mm左右。這一結果表明,從2018年開始,白格滑坡的后緣牽引滑移區不斷擴大,前緣阻滑段也有持續向下滑動的趨勢,后緣及中部因前緣阻滑段的滑移變形已形成整體滑動的趨勢。

根據3個階段的變形演化特征,在2017年3月18日—2018年10月3日期間內,可將白格滑坡的整體形變模式描述為:后緣滑移牽引(牽引區)—前緣阻滑段漸進滑移(阻滑區)—中后部整體變形(主滑區)。這種形變模式與鄧建輝等[16]對白格滑坡形變區的劃分一致。

3.3 典型監測點垂直方向地表累積形變分析

鑒于InSAR技術獲取較大形變量(m級別)較為困難,及超出(mm級別)規定的界限的形變量很難捕獲。則本文只能從該技術微小形變的角度去深入分析白格滑坡的地表微小形變的特征,繪制了白格滑坡垂直方向地表形變速率矢量圖(圖8),并在滑坡體上選取典型監測點獲取其地表累積形變量。

圖8 白格滑坡垂直方向形變速率矢量圖Fig.8 Deformation rate vector graph of the Baige landslide in vertical direction

其中,1～3號監測點位于滑坡邊界,4～6號監測點位于滑坡牽引區,7～9號監測點位于滑坡主滑區,10～11號點位于滑坡阻滑區。各監測點在2017年3月18日—2018年10月3日這一時間段內的累積形變量如圖7所示。需要指出的是,圖7中的起始形變為0 mm,未考慮2017年3月18日之前的歷史累積形變量。

由圖8和圖9可知,在2017年3月18日—2018年10月3日這一時間段內有以下結論。

圖9 白格滑坡典型監測點垂直方向累積位移Fig.9 Vertical cumulative displacements of the typical monitoring points of the Baige landslide

(1)位于滑坡邊界的1～2號監測點其垂直方向形變速率為-5.1～6.8 mm/a,累積形變量為-5～-10 mm,而3號監測點的垂直方向形變速率為-5.1～-17.0 mm/a,累積形變量約為-30 mm。值得注意的是,3號監測點在2017年和2018年的雨季具有較為明顯的下滑跡象。

(2)位于滑坡牽引區的4～6號監測點,其垂直方向形變速率為-17.0～-67.0 mm/a,結合氣象降雨數據,后緣監測點受降雨影響較小,累積形變量與時間呈線性關系。最終累積形變量持續增大,最終為-60～-90 mm,表明該區滑坡體呈持續下滑狀態。這一結果與楊成業等[20]對白格滑坡地表垂直方向的形變量監測結果接近。

(3)位于滑坡主滑區的7～9號監測點,其垂直方向形變速率約為-5.1～ 26.6 mm/a,累積形變量呈現先負后正,最終為0～8 mm,結合氣象降雨數據,監測點在雨季累積形變量波動較大,受降雨因素影響明顯,整體地表形變以抬升為主。表明該區滑坡體在滑坡災變前出現地表隆起。究其原因,牽引區滑坡體的下滑擠壓以及滑坡體表面松散巖土體滑落堆積導致了主滑區部分區域地表隆起。

(4)位于滑坡阻滑區的10～11號監測點,其垂直方向形變速率為-17.0～ 6.8 mm/a,累積形變量呈現先增大(向下滑移)后減小(坡體隆起)再大幅增加(向下滑移)的特征,最終約為-50 mm。結合氣象降雨數據,前緣監測點受降雨影響敏感,尤其在2018年5月到災變發生前累積形變量都處在持續增大的趨勢。

這一形變特征表明,阻滑區的坡體演化過程較為復雜。首先,阻滑區滑坡體在2017年4—5月出現明顯的下滑跡象,此后隨著牽引區的持續下滑,受主滑區向下擠壓表現出地表隆起,發揮了阻滑功能。其后,在2018年4月后,阻滑區出現加劇下滑跡象,且這一現象與雨季恰好對應,表明降雨對阻滑區的穩定性有較強烈的影響。一旦阻滑區脆性剪段,則主滑區整體高速剪出,緊接著牽引區啟動,白格滑坡發生災變。

4 白格滑坡災變時間預報

Fukuzono[12]基于大量降雨誘發滑坡的大型場地試驗,提出了用速度倒數法預報滑坡破壞時間,即通過線性擬合加速蠕變階段的速度倒數-時間曲線,得到擬合線與時間軸的交點,把交點所對應的時間作為預報的滑坡破壞時間。Voight[13]用經驗公式描述滑坡加速蠕變階段速度的變化趨勢,并依據速度和加速度的關系推導滑坡失穩時間。Fukuzono-Voight速度倒數模型預報滑坡破壞時間的表達式為

(4)

(5)

(6)

(7)

在確定滑坡的加速階段后,根據式(7)利用速度倒數擬合直線與時間軸的交點即可確定滑坡的破壞時間tf。

由InSAR形變結果(圖8和圖9)可知,白格滑坡的后緣牽引區形變最為顯著,故選取位于該區6號監測點的形變數據進行滑坡災變時間的預報。考慮到白格滑坡形變演化的漸進性,累積時間長度按照表1取0～564 d(即2017年3月18日—2018年10月3日)。6號監測點垂直方向的累積形變量和形變速率倒數隨時間變化的散點圖如圖10和圖11所示。Zhou等[15]指出滑坡災變時間預報的精度取決于加速蠕變階段起點的確定合理與否,因此,根據圖11選取2018年5月24日作為白格滑坡進入加速變形的起始時間點,進而利用該時間點及之后的速度倒數進行線性擬合,擬合方差R2= 0.93,擬合效果良好。根據該擬合曲線與時間軸的交點,確定白格滑坡的預報災變時間為2018年10月19日,與其實際災變時間2018年10月11日僅相差8 d,預報精度較高。這一預報結果表明,利用SBAS-InSAR技術獲取的滑坡地表時序形變數據在開展大型滑坡的災變時間預報和防災減災中具有良好的應用前景,可為決策管理者判斷類似大型滑坡處于何種演化階段以及制定防災減災預案提供有價值的參考。

圖10 白格滑坡6號監測點累積形變量Fig.10 Cumulative deformation of No.6 monitoring point of Baige landslide

圖11 白格滑坡6號監測點垂直方向速度倒數-時間擬合曲線Fig.11 Fitting curve of the inverse velocity for the monitoring point No.6 of the Baige landslide

5 結論

以西藏昌都市江達縣波羅鄉白格滑坡為例,通過采集滑坡災變前2017年3月18日—2018年10月3日時間段內的Sentinel-1A升軌數據,利用短基線雷達干涉測量技術(SBAS-InSAR)獲取了該滑坡沿垂直方向的地表形變時序數據,從年平均形變速率、累積形變量和典型監測點的時序累積位移3個方面,分析了該滑坡災變前的形變演化特征,在此基礎上利用速度倒數模型對該滑坡的災變時間進行了預報。得出如下結論。

(1)白格滑坡災變前表現出明顯的形變分區。強烈形變區的垂直方向年平均形變速率為-39～-67 mm/a,次強烈形變區的垂直方向最大形變速率達-39 mm/a,一般形變區的垂直方向年平均形變速率為-0.7～-22 mm/a,基本穩定區主要位于右前緣臨江段,該部位在災變前起阻滑作用。

(2)白格滑坡災變前的形變可分為3個階段,其中在2018年開始后強烈形變區的垂直方向最大累積形變量達-115 mm,已形成整體滑動的趨勢。形變模式可概括為:后緣滑移牽引—前緣阻滑段漸進滑移—中后部整體變形,與此模式對應形成牽引區、阻滑區和主滑區,與現場調查結論一致。降雨對阻滑區的形變有較強烈的影響。

(3)基于Fukuzono-Voight速度倒數模型預報白格滑坡的災變時間為2018年10月19日,與其實際災變時間相比滯后8 d,表現出良好的預報效果。目前國內外通過速度倒數法基于孔徑雷達的邊監測對災變進行預報案例還比較少,如何進一步改善預報模型和方法,驗證其實際工程應用仍需要更加深入考究。

(4)需要說明的是,白格滑坡災變前的微小形變演化特征及災變時間預報與實際情況吻合度高,得益于良好的影像干涉效果。由于白格滑坡在近20年中滑動量過大,超出了SAR影像監測的范圍(mm級別),其僅可以獲取微形變,而無法獲得滑坡體m級的具體變形的特征,后續尚需進一步研究。