陈俊伊, 李为乐, 陆会燕, 李鹏飞, 铁永波
(1.成都理工大学地质灾害防治与地质环境保护国家重点实验室,四川 成都 610059;2.中国电建集团贵阳勘测设计研究院有限公司,贵州 贵阳 550081;3.中国地质调查局成都地质调查中心,四川 成都 610081)
水力发电是我国能源供给的重要组成部分,是我国水电工程建设的主战场,主要集中在西部山区,特别是青藏高原东南部横断山脉的三江地区。该区域具有地势高差大、河谷切割深、地质构造复杂、地震频发等特点,导致该地区斜坡稳定性差,地质灾害分布广泛,严重影响了水电工程的施工建设和运营安全[1-2]。例如: 2018年10月11日和11月3日金沙江上游白格滑坡发生两次失稳,形成了滑坡-堰塞湖-溃决洪水灾害链,不仅威胁了下游居民的生命财产安全,还影响了川藏铁路及下游水电站的规划建设[3]。澜沧江丰富的水力资源对于我国的水电工程有着非常重要的意义,但复杂的地质环境也为水电工程的建设增添了难度,为了避免出现滑坡等灾害威胁工程建设,地质灾害的早期识别尤为重要。
InSAR(interferometric synthetic aperture rader)技术通过分析干涉信息来获取地表形变信息,可以在地质灾害识别工作中常规地面调查手段难以快速大面积实施的区域进行崩滑隐患识别[4-5]。InSAR技术对地表形变的探测方法主要有雷达差分干涉测量(D-InSAR)技术及时序InSAR技术等。时序InSAR技术基于雷达影像序列探测地表形变的时空演化规律,该技术使InSAR技术对地表形变探测的精度达到亚厘米级[6]。时序InSAR技术包括永久性散射体合成孔径干涉测量(PS-InSAR)[7]、小基线集干涉测量(SBAS-InSAR)[8]以及干涉堆叠法(Stacking-InSAR)[9]等。PS-InSAR技术需要通过选取永久散射体(permanent scatters,PS点)提高相干性,减少空间失相干对结果的影响,要求PS点在长时间序列上保持稳定散射[10],在山区环境PS点目标主要为裸露的岩土体,其分布并不广泛,因而利用PS-InSAR技术进行广域崩滑隐患识别效果并不理想。相对来说,Stacking-InSAR和SBAS-InSAR技术主要利用在山区环境广泛分布的分布式目标(distributed scatters,DS点)进行干涉,在广域崩滑隐患识别和监测应用中效果较好[1,11-13]。因此本文选用SBAS-InSAR技术和Stacking-InSAR技术对覆盖研究区的SAR影像数据进行处理,对澜沧江卡贡乡—如美镇段进行崩滑隐患识别和重点区域时序形变监测,为下游水电工程的建设提供保障,也为研究具有较强隐蔽性的高位、高危崩滑灾害提供参考。
本文主要以澜沧江干流卡贡乡—如美镇段及干流两侧5 km的范围内为研究区。澜沧江发源于青藏高原唐古拉山脉,位于我国西南部[14]。研究区卡贡乡—如美镇段为澜沧江中游,位于青藏高原东南部西藏自治区昌都市境内,属于横断山系,此处因金沙江、澜沧江、怒江自北向南并行流淌,也称为三江并流区[15](图1)。
研究区属于高山峡谷区,河谷深切于横断山脉,地形起伏大且岸坡高陡,局部甚至为陡壁[16]。区域内的地形地貌由地质构造运动等控制,而构造运动则会在区域内形成活动断裂。根据地震活动断层探察数据中心提供的数据可知,研究区内的活动断层主要为发育于全新世的巴青—类乌齐断裂和澜沧江断裂,其中澜沧江断裂贯通了研究区,且该断裂活动较为强烈,曾多次发生4.7级以上地震,对研究区山体的稳定性造成较大影响[17]。研究区按照区域地层划分,主要属于华南地层大区,羌北—昌都—思茅地层区,唐古拉—昌都地层分区。研究区沉积岩分布较广,出露地层主要有泥盆系、石炭系、二叠系、三叠系、侏罗系、白垩系、古近系及新近系,第四系出露相对较少[18]。区域内属于典型的高原温带半湿润季风性气候,干湿季分明,6—9月集中降雨,其余时间降水量极少,降水峰值月的多年平均降水量比降水低值月的多年平均降水量多152.8 mm,降水量分配极为不均[19]。
图1 研究区位置和雷达影像覆盖范围
本研究主要收集了覆盖研究区的Sentinel-1雷达卫星影像数据。利用GAMMA软件通过Stacking-InSAR以及SBAS-InSAR两种方法对收集的雷达影像数据进行处理,探测出研究区地表形变异常区,进而将探测出的地表形变区与2009年—2022年5月的Google Earth三维卫星影像平台(分辨率0.5~1.0 m)和天地图二维卫星影像平台(分辨率1.0 m)多时相高分辨率光学遥感影像进行对比,利用崩滑隐患的形态特征和形变迹象,圈定研究区崩滑隐患点。
本次研究采用欧空局发布的Sentinel-1雷达卫星影像数据[20]。本文选用的雷达卫星数据为IW工作模式,通过渐进式扫描地形观测获取空间分辨率为5 m×20 m的3个子条带,使其单景覆盖范围达250 km×170 km[21]。由于研究区为高山峡谷地区,为防止产生的几何畸变加大形变探测的难度,甚至形成监测盲区,采用联合升轨、降轨探测结果对研究区进行解译,提高隐患探测的准确性[22]。本次研究共收集了覆盖研究区2018年1月5日—2021年6月30日的103景升轨影像数据,以及覆盖研究区2018年1月12日—2021年6月25日的104景降轨影像数据,影像覆盖范围如图1所示,同时收集了30 m分辨率的SRTM-DEM数据以提供地形地貌特征参数。
为避免大气效应等因素对解译结果造成影响[13],同时相对快速地获取地表形变信息,首先通过Stacking-InSAR方法进行大范围形变探测,并获取隐患点的空间形态特征以及年均形变速率,再通过SBAS-InSAR抑制大气效应以及DEM误差的干扰对典型隐患点进行详细形变探测,获取隐患点在时间序列上的形变特征。
2.2.1 基于Stacking-InSAR的形变探测方法
由于地表形变相位与时间成线性变化,而大气效应的影响与时间无关,因此通过对多幅差分干涉图所获取的解缠相位进行加权平均,可以得到更高精度的地表形变相位。具体步骤为: ①通过空间基线与时间基线设置,将配准后的SAR影像进行差分干涉,获得N个干涉对; ②对干涉对进行滤波、解缠; ③选取解缠相位图较为理想、解缠结果在工作区无明显跳变的组合,加权求取平均值,估算线性相位速率,计算线性相位速率的公式为[23]
。
(1)
式中:Vφ为线性形变对应的相位变化速率;N为参与Stacking计算的干涉对数;ΔTi为第i个干涉对的成像时间间隔,d;pi为第i个差分干涉图的解缠相位。
2.2.2 基于SBAS-InSAR的形变探测方法
Stacking-InSAR技术的计算过程较为简单,能够快速筛查出地灾隐患点,但不能反映隐患点在时间序列上的形变特征。因此,为了更好地了解其形变过程,需要对典型形变点进行SBAS-InSAR处理。
。
(2)
利用外部DEM模拟地形相位,去除地形相位和参考椭球面引起的平地相位,得到m幅差分干涉图,对干涉图进行空间滤波后选取相干性较高的目标点进行相位解缠。根据时间基线、空间基线与干涉相位估算线性形变速率以及残余相位,其中残余相位中主要包括大气相位、非线性形变相位以及噪声; 对残余相位以空域低通、时域高通的方式进行滤波,分离出大气相位并将其去除,获取非线性形变相位及噪声; 通过多视及滤波等方法去除噪声,得到非线性形变相位,将其与线性形变相加得到完整的形变相位,获取累积形变量[11]。
Stacking-InSAR技术通过设置时空基线阈值形成小基线对,提高探测结果的准确性。本文设置升轨影像数据的最大时间间隔阈值为72 d,最大空间基线阈值为186.2 m,共生成500对干涉对,通过对解缠结果进行筛选,筛除141对有明显相位跳跃及大范围大气噪声的干涉对。设置降轨影像数据的最大时间间隔阈值为84 d,最大空间基线阈值为199.5 m,共生成505对干涉对,通过对解缠结果进行筛选,筛除118对有明显相位跳跃及大气噪声的干涉对。分别对得到的升轨、降轨Stacking-InSAR结果进行解译,并将解译结果与Google Earth卫星影像上获取的高精度多时序光学影像图进行对比,提高解译结果的精确性并判别隐患点类型。滑坡隐患在光学影像中的主要解译标志为: 后缘出现高亮裂缝,侧缘边界出现剪切裂缝,且有扩大贯通趋势,部分坡体上出现新鲜下错台坎,局部有垮塌现象,整体形态多呈圈椅状、梨状、舌状,有时呈不规则形状[24-25](图2)。崩滑隐患一般发育于地形陡峭的山崖或斜坡上部,岩体节理裂隙较为发育[26](图3)。
(a) Stacking-InSAR年均形变相位 (b) 光学影像
(a) Stacking-InSAR年均形变相位 (b) 光学影像
研究区在观测期内升轨、降轨影像共探测出明显隐患点共149处(图4、图5),其中升轨影像共探测出隐患点68处,降轨影像共探测出隐患点95处,升轨、降轨探测出的相同隐患点14处,且隐患点多分布于研究区东岸。澜沧江卡贡乡—如美镇段沿江共探测出隐患点65处,这些隐患点皆位于澜沧江干流的陡壁上,范围广,形变量较大,具有堵江风险,严重威胁了沿线水电站的施工建设和运营安全。
图4 Sentinel-1升轨影像重点区Stacking-InSAR年均形变相位
图5 Sentinel-1降轨影像重点区Stacking-InSAR年均形变相位
为了进一步分析研究区隐患点的形变特征,本文对研究区内的隐患点进行筛选,选取其中3处较为典型的隐患点进行SBAS-InSAR处理,获取隐患点在时间序列上的形变信息。
YH1滑坡隐患位于澜沧江干流左岸,整体呈长舌状,为古滑坡堆积体复活。隐患点长约3 337 m,宽约1 827 m,面积约2.9 km2,隐患点底部海拔3 197 m,顶部海拔4 302 m。隐患点坡体后缘有明显下错迹象,坡体上无植被覆盖,冲沟遍布,坡体上多处可见崩塌,受风化侵蚀较严重,岩体质量较低,为泥石流发育提供了充足的物源(图6(a))。根据SBAS-InSAR探测结果(图6(b)),YH1滑坡隐患点视线向年均形变速率最大达到30.2 mm/a,如图6(c)所示。在隐患点内选取3处特征点P1、P2及P3绘制时间序列曲线(图7),3处特征点的形变趋势相同。在监测期内该隐患处于匀速形变阶段,最大累积形变量已达到102.6 mm,表明该滑坡隐患形变较为强烈,整体不稳定,且该隐患点范围大,具有极大堵江风险,影响了坡体上的房屋安全,同时对澜沧江下游水电站工程建设具有较大威胁。
(a) 光学影像(b) Stacking-InSAR年均形变相位 (c) InSAR年均形变速率
图7 YH1滑坡隐患特征点时间序列曲线
YH2滑坡隐患点位于澜沧江右岸,整体呈圈椅状,坡体岩层属于反倾岩层,其变形破坏模式为倾倒变形,隐患点长约1 174 m,宽约1 285 m,底部海拔3 113 m,顶部海拔3 744 m,面积约1.3 km2。光学影像上可见隐患点坡体后缘有下错迹象,坡体上冲沟遍布,无植被覆盖,坡体受风化侵蚀较严重,前缘多处发生崩塌(图8(a))。选取隐患点内的3处特征点P1、P2及P3绘制时间序列曲线(图9),时间序列曲线显示3处特征点均呈匀速变形阶段,在曲线上某些观测值出现小幅度跳动,可能由大气效应等因素形成的噪声未彻底清除引起。根据SBAS-InSAR结果可知(图8(b)),该隐患点形变速率较大,视线向最大年均形变速率达54.8 mm/a(图8(c)),在观测时间内特征点的最大累积形变量达127.1 mm,可见该隐患点形变较为强烈,且形变范围较大,有堵江风险,对道路、对岸居民及下游水电站建设具有较大威胁。
(a) 光学影像(b) Stacking-InSAR年均形变相位(c) InSAR年均形变速率
图9 YH2滑坡隐患特征点时间序列曲线
YH3滑坡隐患点位于澜沧江左岸,宽度约813 m,长度约1 114 m,面积约0.8 km2,隐患点底部海拔2 828 m,顶部海拔3 699 m,地形起伏大。隐患点整体呈圈椅状,坡体上有冲沟,坡体局部有浅层崩塌,坡脚受江水掏蚀较剧烈,致使坡脚较陡(图10(a))。根据SBAS-InSAR技术对隐患点进行形变探测结果(图10(b)),坡体最大视线向年均形变量达41.7 mm/a(图10(c)),图中红色点代表远离卫星视线方向运动。利用SBAS-InSAR技术在坡体上探测出一处强烈形变区,与Stacking-InSAR探测结果一致(图10(b)),强烈形变区最宽处约448 m,长度约为532 m,面积约为0.2 km2。选取强烈形变区中的3处特征点P1、P2及P3绘制时间序列曲线(图11),时间序列结果显示,3处特征点的形变趋势相同,整体呈匀速变形,在监测时间内最大累积形变量达135.0 mm。对比所收集的光学影像,强烈形变区后缘出现裂缝,区域内有局部浅层崩塌,证明该区域处于不稳定状态。结合光学及InSAR解译结果,可以认定该隐患点具有较大的堵江风险。
(a) 光学影像 (b) Stacking-InSAR年均形变相位 (c) InSAR年均形变速率
图11 YH3滑坡隐患特征点时间序列曲线
本文利用Stacking-InSAR技术对覆盖澜沧江卡贡乡至如美镇段的Sentinel-1升轨、降轨影像进行处理,结合光学影像及Stacking-InSAR结果解译出研究区范围内的崩滑隐患点。通过SBAS-InSAR技术对研究区范围内的3处典型滑坡隐患点进行处理,获取其时间序列曲线,进一步分析隐患点的形变特征。综合Stacking-InSAR技术、SBAS-InSAR技术及光学影像的结果,主要获得以下几点认识:
(1)研究区内升轨、降轨影像共探测出隐患点149处,其中升轨影像共探测出隐患点68处,降轨影像共探测出隐患点95处。
(2)研究区遥感影像探测出的隐患类型主要为崩塌和滑坡,崩塌隐患共54处,滑坡隐患共95处。
(3)研究区滑坡隐患点多处于匀速变形阶段,坡体多有明显变形迹象,坡体岩体较为破碎,多有局部崩塌现象,视线向年均形变速率与累积形变量较大,坡体较不稳定,且面积较大,有堵江风险,对澜沧江及下游水电站建设威胁较大。
通过Stacking-InSAR技术及SBAS-InSAR技术可以较为快速准确地探测出研究区崩滑隐患点的数量及形态,但由于卫星飞行方向固定及侧视成像的特点,会不可避免地产生几何畸变,尤其是在高山峡谷地区。研究区处于高山峡谷地貌,地形起伏较大,本文所解译的隐患点结果可能由于几何畸变等因素的影响出现遗漏,同时大气效应等噪声未能彻底清除也会对探测结果产生影响。因此,为了能确保探测结果的准确性,还需要后期通过机载LiDAR探测以及野外勘察等手段进行验证。