木鱼包滑坡形变特征的InSAR监测分析

2022-04-15王尚晓李士垚牛瑞卿

长江科学院院报 2022年4期

王尚晓, 李士垚, 牛瑞卿

(1.中国地质调查局南京地质调查中心，南京 210016； 2.中国地质大学地球物理与空间信息学院，武汉 430074)

1 研究背景

三峡地区是我国地质灾害重大隐患区，而滑坡又是三峡最常见、分布最广的地质灾害。由于区域构造运动活跃，地形陡峻，河谷深切，大量黄土和松散堆积物在山坡和沟谷内堆积。近年来受到城镇化进度加快和不合理的人类活动的影响，生态环境破坏严重。特别是受三峡工程建设的影响，库区水位周期性波动不仅复活了许多古滑坡，而且使一系列新的滑坡灾害出现，种种因素造成了三峡库区存在大量地质灾害隐患[1]。地质灾害的存在不仅对三峡库区当地人民的生命与财产安全和重要设施等构成严重威胁，而且还影响到长江上下游乃至更大范围的地区。由于滑坡是复杂地质现象，由多种因素共同作用，目前我们对滑坡产生和发展不确定性的认知还存在许多不足，无法对此做出确定性决策[2]。常用的边坡变形监测技术主要分为2种：接触式和非接触式变形监测[3-4]。传统的监测传感器包括张力计、测斜仪、全球卫星导航系统(GNSS)监测设备和压力计等，可以测量滑坡体的形变值。然而监测网的建设往往费时费力，维护费用也很高。传统的监测方法是点测量，测量结果不能反映整个滑坡体的变形特征。此外研究区气候类型为季风气候，光学传感器在获得有效的时间序列测量时往往受到限制，难以获得效果较好的连续影像。总的来说，传统监测手段工作效率低，很难在大范围内实现对滑坡目标的高精度监测，尤其是一些缓慢变形的滑坡，其监测面临更大的挑战。

近些年来，随着技术水平的发展，合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)逐渐成为解决上述问题的有效方案之一，可以克服以往手段的不足，现已成为滑坡监测与识别的热点研究方法[5-8]。InSAR技术能够在大区域监测厘米级或毫米级变形，且不受气候条件影响，提供监测区域高分辨率时空地表形变信息[9-10]。时序InSAR分析可以实现对地表位移的远程探测和特征表示，精度可达亚厘米，分辨率可达数米或者数十米级，特别适用于活动性慢速滑坡的研究[11-13]。2001年Ferretti等[14]以长时间范围内SAR数据中相位和幅度都保持高相干的稳定点目标为研究对象，提出了永久散射体干涉测量(PS InSAR)技术，相比于DInSAR (Differential Interferometry Synthetic Aperture Radar)技术，此方法克服了时间和空间失相干等不利因素的影响，能够减弱大气扰动的影响，最终获取长时序的缓慢地表形变场，可以达到米级的高程精度，在理论上可提高雷达视线方向的变形精度到毫米级，甚至亚毫米级。2002年Berardino等[15]改进了时间序列InSAR模型，提出了一种名为小基线集干涉测量(SBAS InSAR)的新算法，SBAS InSAR算法通过设置较小的空间基线、时间基线阈值，将SAR影像进行组合由此组成InSAR干涉对，保证各干涉图的高相干性，有效地减弱了长时间基线对失相干的影响，最后基于奇异值分解(SVD)方法，多个干涉图子集使用最小二乘法求解得到，提高了变形监测的时间分辨率，并保留了传统DInSAR方法大面积覆盖、高时间分辨率的优点。2008年，Wang等[16]在三峡巴东地区利用QPS-InSAR技术对滑坡进行了初步分析，成功发现了两处活动滑坡。Bianchini等[17]于2012年基于多期TerraSAR-X数据，在意大利的Calabria地区使用PS InSAR技术进行了滑坡识别、状态更新和滑坡编目等研究，结果表明在研究区内己知的42个滑坡中，有10个滑坡的边界发生了变化，并提取出了5个新的滑坡。2013年Handwerger等[18]在旧金山东湾山(EBH)地区，对1992—2006年不同卫星InSAR数据集的分析定义已知活动滑坡的空间范围，求解一致的平均滑坡位移率(25～39 mm /a)，揭示了冬季最大降水与滑坡峰值速度之间的平均时间间隔(1～3个月)。

本研究旨在利用InSAR技术对木鱼包滑坡进行地表变形监测,并对其进行时空变形模式分析，以期能够替代传统的监测手段，为滑坡的监测与治理提供了新的思路。

2 InSAR技术处理方法

永久散射体是一种高相干点，具有稳定的散射特性，不受时空失相干的影响。永久散射体方法利用时间序列图像中这些永久散射体的相位信息提取变形信息和数字高程模型(DEM)误差。PS InSAR不受基线去相关的影响，即使基线比临界基线长，也可以形成一个完整的干涉对，可以实现分米DEM精度和毫米级的地表变形精度[19]。

首先从覆盖研究区的N幅SAR影像中选取一景作为主影像，并配准主影像与剩余的所有影像，配准精度不低于1/8个像元可以满足配准精度。然后与InSAR的原理一致，对每组影像对进行差分干涉处理，得到满足条件的M幅差分干涉对。

此时在干涉对上进行PS点的识别，主要基于时间序列SAR图像中相同像元强度值的离散特性来实现PS点的识别，该特性可用幅度离散指数DΔA表示，其计算公式为

(1)

式中σΔA与μA分别表示干涉幅度的标准差和均值。振幅离差指数值越小则表示在整个监测周期内像素点越稳定。根据此原理设置相应的阈值，则可以选择出满足条件的PS点。

干涉相位包含地形、地面运动、大气和噪声等相位的影响，而地表形变信息获取只关注地面变形引起的相位，因此需要去除其他相位。参考面分量可以通过相关网站发布的卫星精密轨道数据进行估计和去除，DEM数据可以模拟出地形相位，根据模拟的结果可以实现地形相位的去除。残余相位包括线性和非线性形变相位、噪声相位、高程误差和大气延迟相位。

一般认为，在同一幅图像中，小研究区(1 km邻区)大气延迟相位的空间相关性很强。因此，将选取PS点连接形成Delaunay不规则三角网，沿网格边界差分可以减弱甚至消除大气延迟相位的影响。

而对于并无统一规律的地表形变相位φdefo，通常分为线性变形和非线性变形两类进行求解,即

式中：φlinear表示线性形变相位；φnonlinear表示非线性形变相位；λ表示雷达波长；T为干涉对的时间基线；Δv表示相邻PS点之间的线性形变速率差。通过沿Delaunay不规则三角形网格边缘差分，可以得到相邻两点之间的相位差Δφ为

(3)

式中：B表示干涉对的空间基线；R为天线与地面目标点间的斜距；θ表示雷达波的入射角；Δε表示相邻PS点间的高程误差增量；Δφres为残余相位。

对于M幅干涉图，每条PS点组成的网格边都有由M个方程组成的差分相位方程组，根据二维周期图解算方法，当|Δφres|<π时，便能得到Δv和Δε的解。剩余的相位为非线性形变相位、大气延迟相位和噪声相位，在解算过程中引入时相相干因子，利用解空间搜索法，根据大气延迟相位和非线性变形相位的不同时空特性，采用时域和空域滤波方法估计大气延迟相位并且有效降低噪声相位，便可得到非线性形变信息。最后将非线性形变与线性形变相加便可得到总的形变量。

3 研究区概况及数据源

3.1 研究区概况

木鱼包滑坡地处长江右岸，距三峡坝址56 km，位于湖北省秭归县沙镇溪镇范家坪村(110°29′52.45″E，31°01′55.9″N)。滑坡体平均长、宽分别为1 500 m和1 200 m，面积为180万m2，厚度平均50 m，体积在9 000万m3左右，主滑方向为20°，边坡为顺向坡，平均坡度20°，但坡体上地形起伏较大(图1)。木鱼包滑坡是一个古老的大型岩质滑坡，滑带土龄相当于晚更新世早期。据推测，水平滑动距离>500 m。松散堆积层构成滑坡表层，下层为扰动层状石英砂岩。表层主要由冲洪积壤土和泥质碎石层、残坡积壤土和崩坡积块石组成；下层滑坡主体为原香溪组中部石英砂岩、含砾石英砂岩。第二层为弱强透水层，香溪组下段石英砂岩、泥岩、粉砂岩等地层，地表堆积物渗透性差，岩石破碎，连通性差，裂隙发育，透水性好，排水通畅[20]。

图1 木鱼包滑坡全貌Fig.1 Overall view of Muyubao landslide

3.2 数据源

日本宇航局(JAXA)于2014年3月24日将ALOS-2卫星成功发射。与其他合成孔径雷达卫星不同，ALOS-2卫星采用L波段，并且由于具有高分辨率，可以获得通过植被覆盖面的反射数据，因此山区林区InSAR的相干性不会因此降低，更有利于提取坡面运动[21]。ALOS-2卫星重访周期为14 d，具有多种成像模式，可进行灵活部署与观测，最高空间分辨率可达1 m×3 m。

结合区域特征和存档数据情况选取Ultra-Fine模式3 m分辨率的ALOS-2数据，时间周期为2016年8月至2017年10月之间共22景，获取时间见表1。

表1 ALOS-2数据列表Table 1 List of ALOS-2 data

研究区高分辨率遥感影像基于获取到的多时相GF-2和GF-1卫星数据。气象和库水位数据分别从中国气象数据网(http:∥data.cma.cn/)和中国长江三峡集团有限公司网站(https:∥www.ctg.com.cn/)获取并进行预处理。

DEM数据基于美国航空航天局(NASA)发布的30 m分辨率的SRTM产品，能够精确提供研究区域的高程数据，提高去平地结果的可靠性。

地质灾害数据从三峡库区地质灾害监测预警指导中心灾害地质数据库中获取，地层岩性和地质构造数据以区域地质图为底图进行矢量化得到。

4 试验结果及分析

4.1 滑坡形变特征

根据SAR影像成像几何可知，InSAR技术处理得到的雷达视线方向(LOS向)形变结果可以分解为LOS向在水平面的投影(地距向)和垂直向两部分。根据公式VLOS=VHcosθ(VLOS为LOS向形变，θ为雷达入射角，VH为垂直方向形变)，可将LOS向的形变换算到垂直方向[22]。

图2是基于ALOS-2数据处理得到的木鱼包滑坡垂直方向形变速率。其中负值表示地面沉降速率，正值表示地面抬升的速率，绝对值表示地面形变速率的大小。由形变速率结果可知，在整个监测周期内，其变形速率范围为-113～135 mm/a，最大形变速率为135 mm/a，平均形变速率值为10.4 mm/a；滑坡的后缘形变速率最大，其次是滑坡的中部，滑坡前缘的形变速率较小。

图2 木鱼包滑坡形变速率Fig.2 Deformation rate of Muyubao landslide

将InSAR形变结果投影至与GPS相同的方向，图3为InSAR结果与典型GPS监测点形变值的对比结果。通过两者相比，可以发现InSAR的结果与GPS的结果吻合度相对较高，每个监测点与InSAR结果间的形变量平均差值最大为5.3 mm，最小为0.12 mm，其中点ZGX297间标准差最大，为10.11 mm，ZGX295点间的标准差最小，为4.1 mm。由于卫星为近南北飞行，飞行方向与滑坡的主滑方向近似，造成InSAR结果对南北方向的形变不够敏感，点间的误差主要由于起算基准和形变方向不完全相同而造成的，2种数据的结果并不能完全一致。但通过2种数据的对比，可以发现点间的误差在可以接受的范围内，也能证实InSAR结果的可靠性，2种结果都表明该滑坡呈现线性缓慢变形的趋势。

图3 InSAR结果与GPS数据对比Fig.3 Comparison between InSAR results and GPS data

图4表示转换为垂直向的时序累积形变。从图4可以看出监测结果能够覆盖整个滑坡范围，区别于点监测，InSAR能够得到面监测的效果。其中在2016年10月之前，木鱼包滑坡整体没有明显的地表形变迹象。2016年10月后，滑坡中东部部分地区开始出现变形迹象，累计形变逐渐增加。随着时间的推移，形变逐渐增大，最大累积形变达到100 mm左右。而滑坡靠近前缘的区域一直保持稳定，前缘可能由于临空作用，在整个周期内部分区域也出现变形，但形变一直较小。从累积形变图上可以看出，木鱼包滑坡的形变具有较为明显的空间差异特性，滑坡中部和右侧变形最强，后缘滑体次之，左侧和前缘滑体形变最小。

图4 木鱼包滑坡累积形变Fig.4 Cumulative deformation of Muyubao landslide

为了进一步了解木鱼包滑坡空间形变特征，以木鱼包滑坡的形变速率图为基础，做出4条纵剖面，将4条纵剖面的监测点形变速率值与距离后缘的距离作堆积折线图，可将不同剖面的各监测点形变速率值大小进行比较(图5)。剖面A整体形变速率较小且一致，而近前缘位置的形变速率稍小。剖面B、C和D的位移速率特征基本一致，滑坡中部距后缘200～300 m处的形变速率最大，形变随距后缘的距离增大而减小，靠近前缘的位置形变稍有增大，前缘位置的形变速率明显比后缘和中部要小。这与该滑坡的坡体结构和地质力学模式有直接的关系，该滑坡为典型的滑移-弯曲型滑坡溃决后形成的古岩质滑坡，滑坡体中后部岩层倾角为25°～27°，前部岩层变为平缓甚至反倾状，因而该滑坡在前部形成阻滑段，后缘驱动块体在下滑力作用下发生位移，越靠近锁固段位移量越小；中后部沿软弱滑带顺层蠕滑，位移量明显较大，为坡体主要下滑段，形成驱动力。滑坡前缘具有良好的临空条件，不断被库水侵蚀，导致位移量增大。

图5 木鱼包滑坡剖线位置及累计形变速率Fig.5 Section of Muyubao landslide and cumulativedeformation rate

选取滑坡体不同部位的特征点目标对其空间变形特征进行进一步分析。图6为InSAR特征监测点目标的空间位置分布及其时间序列累积位移曲线。由时序形变曲线可知，监测点P6、 P10和P11累积位移最大， 2016年8月—2017年10月的累积位移超过100 mm，这3个点形变规律一致且线性形变特征明显；位于滑坡中前部的监测点P2、 P3、 P5、 P13、 P14和P15累积位移中等，在监测周期内位移达到60～70 mm。特征点的时序运动特征与剖线揭示的规律一致，都能够证实木鱼包滑坡的形变具有明显的分区特征，滑坡变形与其自身结构密切相关，并受外界因素影响。

图6 木鱼包滑坡特征点分布及其时序形变Fig.6 Time-series deformation of characteristic pointsof Muyubao landslide

4.2 变形机理分析

首先做出木鱼包滑坡特征点位移与库水位和降雨量曲线(图7)，根据选取的特征点时序位移特征，可将滑坡运动变形分为Ⅰ、Ⅱ、Ⅲ、Ⅳ和Ⅴ共5个阶段。

图7 木鱼包滑坡特征点位移-库水位-降雨量关系Fig.7 Displacement-water level-rainfall relationshipof characteristic points of Muyubao landslide

阶段Ⅰ时间为2016年8月—2016年10月近2个月期间，位移量平均增加20 mm，对应的库水位由145 m增加到164 m，库水位约增加了19 m。各监测点位移呈小幅度上扬趋势，其变形较明显且时间段刚好与库水位上涨的情况相吻合。

阶段Ⅱ为2016年11月—2016年12月近2个月期间，位移量变化较小，对应的库水位由165 m增加到175 m，库水位约增加了10 m。说明库水位高位状态对木鱼包滑坡的变形影响小，并对稳定性有积极的作用。木鱼包滑坡的坡体主要为类基岩的碎裂块体，节理裂隙发育，渗透性较好，但由于木鱼包前部宽度达1 400 m、厚度达110～140 m，库水位能够形成向坡内的水头差，反而对滑坡体的稳定有利。

阶段Ⅲ为2016年12月—2017年2月近2个月期间，位移量平均增加50 mm，对应的库水位由175 m下降到165 m，库水位约减少了10 m。由此说明三峡水库水位在高水位运行并轻微下降的工况条件下木鱼包滑坡的变形最大，滑坡监测数据形变曲线出现较大的“阶跃”，位移-时间曲线呈阶跃上升，滑坡变形对库水位高水位运行较为敏感，为典型的浸水软化和浮托减重型滑坡。在库水位175 m下降到165 m高位运行期间，木鱼包滑坡前部(抗滑段)浮托力接近最大，滑坡抗滑力下降；加上库水位的下降，形成向外的动水压力，两个因素的叠加，形成了一年中位移最大的变化。

阶段Ⅳ为2017年3月—2017年8月近5个月期间，位移量平均增加40 mm，对应的库水位由165 m下降到145 m，库水位约减少了20 m。期间滑坡的位移曲线处于平缓状态，一般月位移在10 mm以下，木鱼包滑坡受向坡外的水头差和对前部抗滑段的浮托效应减小的双重影响，2个作用相互抵消，导致木鱼包滑坡在此期间位移较小但也未停止变形，一致处于缓慢变形的状态。

阶段Ⅴ为2017年8月—2017年10月近2个月期间，位移量平均增加30 mm，对应的库水位由145 m增加到170 m，库水位约增加了25 m，形变类似于阶段A，说明库水位上涨对滑坡变形的影响较大。

通过分析发现库水位上涨和高水位状态对木鱼包滑坡的影响因素比大气降水更为敏感，这主要是由于岩质滑坡的基本特征控制木鱼包滑坡的变形。当库水上升淹没滑坡前缘时，受地表水入渗影响，滑坡体处于饱和状态，滑坡体抗滑段岩体受浮力作用而减轻。同时，库水位高水位浸泡岩体，岩体会遇水造成软化，特别是滑带土会有软化的情况，滑坡体以及滑带力学性质受高水位的持续浸泡作用影响逐渐变差；并且岩体裂隙发育，充水后产生的静水压力变大，含水量在裂隙的变化直接影响了岩体力学性质，多种因素的叠加导致变形加速。

由木鱼包滑坡的上述特征可以得出，该滑坡为典型的滑移-弯曲型滑坡，是由斜坡的结构和岩性等条件共同控制其变形破坏模式。滑坡为中倾外层状斜坡体，前缘并没有临空，但滑坡受长期重力作用的影响，坡体沿软弱面顺层向下蠕滑变形，造成岩层弯曲隆起出现在坡体的中前部，在滑坡的前部形成很大的抗阻力；从地形地貌角度看，前部的隆起增加了滑坡的抗滑效果，地表冲沟也有利于地表水迅速排出坡外。总之，在自然状况下，这些条件均有利于滑坡体的整体稳定，故坡体的稳定性较高。但在三峡水库水位变化、降雨等条件的叠加影响下，滑坡体存在局部失稳的情况。总的来说，滑坡加剧变形的因素包括大气降雨与库水位的上升，在研究周期内滑坡的变形主要受三峡水库水位上涨作用影响。