馬闖, 繆海波, 楊犇, 朱隆奇, 安廣強(qiáng), 楊冰穎

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

在高海拔山區(qū)及峽谷地帶,因強(qiáng)降雨和/或地震誘發(fā)的大型滑坡往往具有隱蔽性強(qiáng)、突發(fā)性高、危害性大等特點(diǎn),給當(dāng)?shù)厝嗣裆踩蛧?guó)民經(jīng)濟(jì)發(fā)展造成了嚴(yán)重威脅,同時(shí)給該區(qū)域的地質(zhì)災(zāi)害防范預(yù)警、防治和風(fēng)險(xiǎn)管控帶來(lái)了極大的挑戰(zhàn)。例如,2017年四川茂縣新磨村滑坡導(dǎo)致80人遇難[1],2018年金沙江流域發(fā)生的白格特大巖質(zhì)滑坡造成了金沙江堵塞,影響了2.1萬(wàn)余人[2]。因此,針對(duì)此類滑坡的早期識(shí)別和災(zāi)變時(shí)間預(yù)報(bào),是避免或減少人員傷亡和財(cái)產(chǎn)損失,提高地質(zhì)災(zāi)害防災(zāi)減災(zāi)水平的重要途徑。

由于大型災(zāi)難性滑坡的演化過(guò)程通?？缭綌?shù)年、數(shù)十年甚至上百年,且在災(zāi)變前均會(huì)產(chǎn)生明顯的地表形變跡象[3]。因此,針對(duì)大型滑坡的地表時(shí)序形變解析結(jié)果不僅可以描述整個(gè)孕育期間的演化規(guī)律,而且對(duì)其運(yùn)動(dòng)學(xué)機(jī)制解析和災(zāi)變前兆信息檢測(cè)(如滑坡破壞前的突然加速)等至關(guān)重要。目前,合成孔徑雷達(dá)干涉測(cè)量技術(shù)(interferometric synthetic aperture radar,InSAR)因其非接觸、超遠(yuǎn)距離、高探測(cè)精度和全天候等特點(diǎn),可快捷獲取地表(mm級(jí)別)變形跡象,在發(fā)生大規(guī)?；?m級(jí)別)前給予災(zāi)變信號(hào);該技術(shù)已在地表沉降、地裂縫發(fā)展、滑坡等地質(zhì)災(zāi)害的監(jiān)測(cè)方面得到了廣泛運(yùn)用[4],尤其是在交通困難的高海拔山區(qū)和深切峽谷地帶的隱蔽性古滑坡解譯與復(fù)活識(shí)別、新生型滑坡災(zāi)變前兆信息檢測(cè)中已展現(xiàn)出獨(dú)特的優(yōu)勢(shì)。

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

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

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

1 滑坡工程地質(zhì)概況

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

圖2 滑坡全貌Fig.2 Landslide panorama

區(qū)域地質(zhì)背景顯示,白格滑坡所在區(qū)域處于金沙江構(gòu)造帶上,受多期構(gòu)造運(yùn)動(dòng),滑坡巖體結(jié)構(gòu)破碎明顯,且經(jīng)多期變質(zhì)作用與強(qiáng)烈風(fēng)化后,碎裂巖體遇水易崩解軟化。根據(jù)文獻(xiàn)[2],滑坡區(qū)出露的地層為元古界熊松群(Ptxn2)深灰色斜長(zhǎng)片麻巖和華力西期蛇紋巖(φw4),其中蛇紋巖帶主要位于滑源區(qū),綠泥石化嚴(yán)重?；露逊e物主要為碎塊石及其破碎物,碎屑成分較為復(fù)雜[16-17]。

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

2.1 數(shù)據(jù)來(lái)源

使用寬幅干涉(interferometric wide swath,IW)模式、單視復(fù)數(shù)(single look complex,SLC)的SAR影像,數(shù)據(jù)來(lái)自歐洲航天局雷達(dá)觀測(cè)站點(diǎn)的哨兵1號(hào)(Sentinel-1)衛(wèi)星。選取23幅C波段Sentinel-1A升軌影像。裁剪后的覆蓋區(qū)域如圖3所示。

圖3 哨兵一號(hào)升軌數(shù)據(jù)裁剪覆蓋范圍Fig.3 Data clipping coverage map of Sentinel-1A

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

表1 白格滑坡Sentinel-1A數(shù)據(jù)時(shí)間序列列表Table 1 Time series of Sentinel-1A data of the Baige landslide

2.2 技術(shù)原理及形變速率獲取方法

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

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

(1)

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

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

Aφ=δφ

(2)

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

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

Bv=δφ

(3)

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

SBAS-InSAR技術(shù)的具體處理方法如下。

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

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

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

步驟2對(duì)干涉相對(duì)進(jìn)行干涉處理,過(guò)程包括:相干性生成、干涉去平、濾波(選取Goldstein濾波方法,增強(qiáng)濾波效果)和相位解纏(選取最小費(fèi)用流方法)。

步驟3軌道精煉和重去平?？紤]到研究區(qū)地形復(fù)雜,山脈較多,故借助Google Earth 選取穩(wěn)定點(diǎn)(如穩(wěn)定交通樞紐等構(gòu)筑物)作為地面控制點(diǎn)(ground control point,GCP),以此來(lái)降低因?yàn)檫x點(diǎn)不當(dāng)而帶來(lái)的誤差。

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

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

3 白格滑坡災(zāi)變前形變結(jié)果分析

3.1 垂直方向形變速率分析

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

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

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

3.2 垂直方向累積形變分析

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

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

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

第3階段：2018年1月24日—10月3日,白格滑坡的地表垂直方向形變特征表現(xiàn)為滑坡后緣的滑移急劇發(fā)展,且滑移牽引區(qū)不斷擴(kuò)大,垂直方向的最大累積形變量已達(dá)-115 mm。與此同時(shí),滑坡前緣也表現(xiàn)出持續(xù)向下滑移且有向中部擴(kuò)展的趨勢(shì),垂直方向的最大累積形變量達(dá)-50 mm左右。這一結(jié)果表明,從2018年開(kāi)始,白格滑坡的后緣牽引滑移區(qū)不斷擴(kuò)大,前緣阻滑段也有持續(xù)向下滑動(dòng)的趨勢(shì),后緣及中部因前緣阻滑段的滑移變形已形成整體滑動(dòng)的趨勢(shì)。

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

3.3 典型監(jiān)測(cè)點(diǎn)垂直方向地表累積形變分析

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

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

其中,1～3號(hào)監(jiān)測(cè)點(diǎn)位于滑坡邊界,4～6號(hào)監(jiān)測(cè)點(diǎn)位于滑坡?tīng)恳齾^(qū),7～9號(hào)監(jiān)測(cè)點(diǎn)位于滑坡主滑區(qū),10～11號(hào)點(diǎn)位于滑坡阻滑區(qū)。各監(jiān)測(cè)點(diǎn)在2017年3月18日—2018年10月3日這一時(shí)間段內(nèi)的累積形變量如圖7所示。需要指出的是,圖7中的起始形變?yōu)? mm,未考慮2017年3月18日之前的歷史累積形變量。

由圖8和圖9可知,在2017年3月18日—2018年10月3日這一時(shí)間段內(nèi)有以下結(jié)論。

圖9 白格滑坡典型監(jiān)測(cè)點(diǎn)垂直方向累積位移Fig.9 Vertical cumulative displacements of the typical monitoring points of the Baige landslide

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

(2)位于滑坡?tīng)恳齾^(qū)的4～6號(hào)監(jiān)測(cè)點(diǎn),其垂直方向形變速率為-17.0～-67.0 mm/a,結(jié)合氣象降雨數(shù)據(jù),后緣監(jiān)測(cè)點(diǎn)受降雨影響較小,累積形變量與時(shí)間呈線性關(guān)系。最終累積形變量持續(xù)增大,最終為-60～-90 mm,表明該區(qū)滑坡體呈持續(xù)下滑狀態(tài)。這一結(jié)果與楊成業(yè)等[20]對(duì)白格滑坡地表垂直方向的形變量監(jiān)測(cè)結(jié)果接近。

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

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

這一形變特征表明,阻滑區(qū)的坡體演化過(guò)程較為復(fù)雜。首先,阻滑區(qū)滑坡體在2017年4—5月出現(xiàn)明顯的下滑跡象,此后隨著牽引區(qū)的持續(xù)下滑,受主滑區(qū)向下擠壓表現(xiàn)出地表隆起,發(fā)揮了阻滑功能。其后,在2018年4月后,阻滑區(qū)出現(xiàn)加劇下滑跡象,且這一現(xiàn)象與雨季恰好對(duì)應(yīng),表明降雨對(duì)阻滑區(qū)的穩(wěn)定性有較強(qiáng)烈的影響。一旦阻滑區(qū)脆性剪段,則主滑區(qū)整體高速剪出,緊接著牽引區(qū)啟動(dòng),白格滑坡發(fā)生災(zāi)變。

4 白格滑坡災(zāi)變時(shí)間預(yù)報(bào)

Fukuzono[12]基于大量降雨誘發(fā)滑坡的大型場(chǎng)地試驗(yàn),提出了用速度倒數(shù)法預(yù)報(bào)滑坡破壞時(shí)間,即通過(guò)線性擬合加速蠕變階段的速度倒數(shù)-時(shí)間曲線,得到擬合線與時(shí)間軸的交點(diǎn),把交點(diǎn)所對(duì)應(yīng)的時(shí)間作為預(yù)報(bào)的滑坡破壞時(shí)間。Voight[13]用經(jīng)驗(yàn)公式描述滑坡加速蠕變階段速度的變化趨勢(shì),并依據(jù)速度和加速度的關(guān)系推導(dǎo)滑坡失穩(wěn)時(shí)間。Fukuzono-Voight速度倒數(shù)模型預(yù)報(bào)滑坡破壞時(shí)間的表達(dá)式為

(4)

(5)

(6)

(7)

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

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

圖10 白格滑坡6號(hào)監(jiān)測(cè)點(diǎn)累積形變量Fig.10 Cumulative deformation of No.6 monitoring point of Baige landslide

圖11 白格滑坡6號(hào)監(jiān)測(cè)點(diǎn)垂直方向速度倒數(shù)-時(shí)間擬合曲線Fig.11 Fitting curve of the inverse velocity for the monitoring point No.6 of the Baige landslide

5 結(jié)論

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

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

(2)白格滑坡災(zāi)變前的形變可分為3個(gè)階段,其中在2018年開(kāi)始后強(qiáng)烈形變區(qū)的垂直方向最大累積形變量達(dá)-115 mm,已形成整體滑動(dòng)的趨勢(shì)。形變模式可概括為:后緣滑移牽引—前緣阻滑段漸進(jìn)滑移—中后部整體變形,與此模式對(duì)應(yīng)形成牽引區(qū)、阻滑區(qū)和主滑區(qū),與現(xiàn)場(chǎng)調(diào)查結(jié)論一致。降雨對(duì)阻滑區(qū)的形變有較強(qiáng)烈的影響。

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

(4)需要說(shuō)明的是,白格滑坡災(zāi)變前的微小形變演化特征及災(zāi)變時(shí)間預(yù)報(bào)與實(shí)際情況吻合度高,得益于良好的影像干涉效果。由于白格滑坡在近20年中滑動(dòng)量過(guò)大,超出了SAR影像監(jiān)測(cè)的范圍(mm級(jí)別),其僅可以獲取微形變,而無(wú)法獲得滑坡體m級(jí)的具體變形的特征,后續(xù)尚需進(jìn)一步研究。