基于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级的具体变形的特征,后续尚需进一步研究。