基于InSAR相位梯度叠加的毛尔盖库岸滑坡隐患快速识别

2023-10-18唐光民戴可人卓冠晨

地震工程学报 2023年5期

唐光民, 戴可人,,3, 卓冠晨,, 沈月,4, 陈晨, 许强

(1. 成都理工大学地球科学学院, 四川成都 610059;2. 地质灾害防治与地质环境保护国家重点实验室, 四川成都 610059;3. 长安大学地质工程与测绘学院, 陕西西安 710054;4. 中国地质大学地球物理与空间信息学院, 湖北武汉 430074)

0 引言

库岸滑坡失稳可能会冲毁沿岸建筑物,或阻塞河道,对库岸地质环境及沿岸居民生命财产安全造成巨大威胁[1]。国内外由库岸滑坡引发的重大灾害时有发生,例如:意大利瓦依昂滑坡不仅造成整个水库堵塞,山体滑坡引发的洪水更是造成2 000余人死亡,多个村庄和城镇夷为平地[2];法国马尔帕塞拱坝事故造成死亡和失踪的人数超过了500人[3];国内湖南柘溪水库塘岩光库岸滑坡造成64人淹溺死亡,24人受伤[4];湖北的黄龙滩水库滑坡造成2 316户居民搬迁,给库区经济带来巨大损失[5]。水库水位短期的剧烈变动,势必改变库区原有地质环境,从而引发滑坡等地质灾害的发生。库岸滑坡灾害分布广泛,多发频发,传统测量手段和群测群防方式耗时耗力,难以进行有效排查[6]。因此,如何快速识别库岸滑坡从而进行有效防灾减灾是亟待解决的重要问题。

近些年,合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar,InSAR)因其具有覆盖范围广、监测精度高、全天时、全天侯、穿云透雾等特点已成为滑坡探测领域的重要技术手段[6-10]。常用的时间序列InSAR方法,例如小基线集干涉测量 (Small Baselines Subset InSAR,SBAS-InSAR)、干涉影像堆叠(Stacking)技术等,被广泛应用于滑坡隐患识别中,大量成功的案例证明了时序InSAR技术具备识别滑坡隐患的能力[11-15]。但是由于时序InSAR技术涉及相位解缠[16]、时序建模解算[17]等复杂步骤,操作相对繁琐,对操作者经验要求较高,需要消耗大量时间,尤其面对海量数据,识别效率不足。相比之下,基于InSAR相位梯度叠加方法受到参数设置的影响较小,同时对微小形变具有更高的敏感性,能够揭示更小尺度的形变变化。此外,由于该方法直接从差分干涉对中计算梯度,避免了相位解缠和时序建模解算等步骤,能够有效缩短数据处理时间,提高数据处理效率。例如,Fu等[18]利用梯度叠加算法获取西南山区InSAR干涉图距离向和方位向的相位梯度,通过训练YOLOv3神经网络进行快速识别,成功探测西南地区3 979处潜在滑坡。Shen等[19]利用改进的Sobel算子求取D-InSAR处理后的9张差分干涉图相位梯度,去除各项形变误差后,采用阈值分割的方法对滑坡隐患进行自动识别,整体精度与人工解译结果相比达到81%,然而其仅针对D-InSAR处理后两景影像进行梯度处理,无法探测长时间序列的滑坡隐患。目前通过InSAR相位梯度快速识别长时间序列滑坡隐患值得进一步研究。

本文在不涉及相位解缠和时序建模解算的情况下,对InSAR差分干涉对进行梯度运算,识别出相位信息不连续的区域。通过去除植被等低相干性区域、几何畸变、水体、积雪等形变误差后,采用相位梯度叠加方式获取相位不连续区域。采用阈值分割生成二值图,并进行边缘提取,以获取发生形变的区域,将获取的形变区域与Stacking-InSAR和SBAS-InSAR处理的结果进行对比验证分析,探究本文方法在滑坡隐患识别方面的特点。最后,统计三种方法从数据配准到识别出滑坡隐患所耗时间,总结本文方法的优势。

1 研究区域与数据

1.1 研究区域

毛尔盖水库隶属四川省阿坝藏族自治州黑水县,地处青藏高原东南缘,属于该区横断山脉中段地区[图1(a)]。南东面与茂县相邻,南西面则与理县相接,北东毗邻松潘[图1(b)],地理坐标为:102°35′～103°30′E、 31°35′～32°38′N[20],属于“北亚热带”气候、川西高原气候区,每年有旱季和雨季,降雨分布不均,主要集中在夏季[21]。

图1 研究区概况Fig.1 Overview of the study area

库区内修建的毛尔盖水电站位于四川省阿坝州黑水县麻窝乡境内[图1(c)],系岷江上游最大的一级支流黑水河流域梯级开发的中游“龙头”水库电站[22],山势陡峭、沟谷纵横,谷底宽广,呈不对称“U”型,坝址河段开阔,山势雄伟,地形完整,坡度为40°至50°,呈“V”型谷。地势西北高东南低,属于典型高山峡谷地貌[23]。死水位为2 063 m,水库正常蓄水位2 133 m,计划2008年3月开工,于2011年11开始投入使用。库区出露的基岩受“5.12”汶川特大地震[24]以及近年来强降雨等影响,加之库区水位线短期快速的变动,触发或加剧了库岸的滑坡灾害[19]。

1.2 研究数据

本次研究采用数据为Sentinel-1卫星数据,Sentinel-1卫星是欧洲航天局(European Space Agency,ESA)发射运营的两颗对地观测卫星组成。搭载波长为5.6 cm的SAR,重访周期12 d,覆盖范围达到250 km×170 km[25],可实现全天时全天候监测,在地质灾害监测方面被广泛应用[26-28]。

本研究使用了从2019年1月至2022年3月覆盖研究区的97景升轨Sentinel-1A卫星数据,主要参数如表1所列。为了保证数据有良好的相干性,设置的时间基线为36天,空间基线为200 m,生成干涉对组合。结合欧空局的精密轨道数据减少参考椭球相位的影响。此外,数字高程模型(Digital Elevation Model,DEM)使用SRTM(Shuttle Radar Topography Mission),可以有效减弱地形相位的影响[29]。

表1 SAR数据主要参数

2 研究方法

本研究利用相位梯度对局部变形非常敏感的特点,采用经过改进后的Sobel算子,对InSAR差分干涉对进行梯度处理,获取InSAR差分干涉对的相位信息。去除植被等低相干性区域、几何畸变区域、水体、积雪等形变误差的影响,通过梯度叠加的方式减弱噪声和大气的影响,获取相位不连续区域。这些区域经过阈值分割生成二值图,最后通过边缘提取快速识别滑坡隐患。具体步骤如图2所示。

图2 技术路线Fig.2 Technical route in this study

首先,是数据预处理阶段。利用收集的Sentinel-1影像结合精密轨道数据生成SLC数据,从中裁剪出研究区域并进行数据配准。结合SRTM-DEM去除干涉对中参考椭球相位和地形相位的影响,选用合适的滤波对干涉对进行滤波处理。为了减少误差的引入,需选取干涉质量较好的干涉对,因此对滤波后的干涉对进行精炼,剔除效果不好的干涉对。

其次,是梯度处理阶段。利用改进的Sobel算子对滤波、精炼后的干涉对进行梯度运算。普通的Sobel算子只采用水平和垂直两个方向的梯度模板进行检测,虽然具有普适性和计算速度快等优势,但对纹理复杂的图像边缘特征识别不足,且无法连续检测边缘。因此,本文利用改进后的Sobel算子即加入了45°和135°两个对角方向的模板识别,使滑坡边界检测更为连续和精细[30-33]。由于在进行梯度处理可能受到水体[34]、积雪[35]、卫星的成像几何带来的几何畸变[11]、植被相干性较低[36]的区域等影响,导致相位梯度不连续,形成形变误差,处理后呈现为相位梯度的突变,类似于滑坡隐患。因此需要对这些形变误差进行有效去除。本文结合目标点相干性,利用相干系数选取合适的阈值,对受植被、噪声等影响严重区域进行掩膜处理。通过多次实验,确定相关性阈值为0.7,可以有效去除植被、噪声等影响。利用卫星参数等信息生成研究区几何畸变数据,通过掩膜几何畸变区域,去除几何畸变的影响。收集毛尔盖库区水位情况与积雪线,通过DEM去除水体和积雪的影响。再利用梯度叠加50对干涉对识别出相位不连续区域。

然后,是快速识别阶段。识别到的相位不连续区域仍然存在由于卫星自身成像等因素引起的形变误差,采用均值滤波对结果进行平滑处理,如形变识别阶段所示。本文通过设定一定阈值生成二值图,进行内部填充后提取相位不连续区域边界,考虑到发生滑动的隐患点区域面积在一定的范围内,而较大面积的区域可能由于大气效应引起的干涉相位的不连续,因此再次通过阈值分割去除掉面积较大区域,最终得到真实形变区域。

最后,是对比验证阶段。将本文识别结果与相应时间段的传统时序InSAR结果进行对比验证,以确保识别结果的准确性。同时总结本文方法对滑坡隐患识别的特点和能力。

3 毛尔盖库岸滑坡隐患识别结果

利用InSAR相位梯度叠加的滑坡隐患快速识别方法针对毛尔盖库区开展滑坡隐患识别的结果如图3(a)所示,共识别出23处滑坡隐患点。综合谷歌光学影像可以明显看出,基于相位梯度叠加的滑坡隐患快速识别方法识别出的滑坡隐患,能够看出较为明显形变迹象,其滑坡边界清晰,后缘可见明显的下错,说明该处已产生较大形变[见图(c)、(e)、(a)、(i)];在空间分布上,这些识别的灾害隐患点主要分布于库岸两侧[图3(i)],这一现象意味着库区水位的快速变动导致土体物理特性下降,对库岸坡体的稳定性产生显著影响;从威胁对象来看,库区两岸坡体可能出现失稳情况,发生较大规模滑坡,从而出现堵江风险,影响水电站正常运行。此外,有5处隐患点附近分布有居民点及道路等基础设施,如果发生垮塌,将对当地居民的生命财产安全造成巨大威胁。

(a为InSAR相位梯度叠加快速识别总体结果;b、d、f、h为典型区域放大结果;c、e、g、i为典型区域Google Earth光学影像)图3 InSAR相位梯度叠加快速识别结果Fig.3 Rapid identification results of the InSAR phase gradient stacking

为了验证本文方法的精度,将快速识别的结果与2019年1月至2022年3月的SBAS-InSAR结果进行对比验证。如图4所示,获取了升轨轨道视线方向上的形变,红色区域(负值)表示远离雷达视线方向移动,蓝色区域(正值)表示靠近雷达视线方向移动,绿色区域则表示相对稳定。研究区几何畸变区域较少,植被覆盖不是特别茂密,且获取数据时间范围较长(2019—2022年),获取的SBAS-InSAR结果精度较高。SBAS-InSAR技术共识别到24处滑坡隐患,从整体识别结果来看,识别出的形变区域主要分布在库岸两侧,其中区域C两处滑坡形变量级较大,达到109 mm/年。基于InSAR相位梯度叠加的滑坡隐患快速识别方法共识别到23处滑坡隐患[图4(a)中白色区域],与SBAS-InSAR方法所识别结果具有高度一致性。从精度验证结果可知本文提出的方法进行的滑坡隐患识别准确率较高。在对广域突发性高的滑坡隐患进行快速动态识别时具有一定的潜力。

(a为SBAS-InSAR和InSAR相位梯度叠加快速识别总体结果;b～e 为典型区域放大结果)图4 SBAS-InSAR和InSAR相位梯度叠加快速识别结果Fig.4 Rapid identification results of SBAS-InSAR and InSAR phase gradient stacking

将本文识别结果与传统Stacking-InSAR结果进行比较,对比结果如图5所示。基于InSAR相位梯度叠加的滑坡隐患快速识别方法共识别到23处潜在滑坡隐患,而传统Stacking-InSAR方法共识别到20处潜在滑坡隐患。总体来说两种方法识别结果吻合度较高。为了更详细的进行比较,将两种方法的能识别到的典型区域放大,具体如图5(b)～(e)所示。这些区域的滑坡隐患表现为形变量级较大,滑动迹象明显,因此Stacking-InSAR方法和基于InSAR相位梯度叠加的滑坡隐患快速识别方法都能够准确地识别出来。另外,有3处滑坡隐患,其形变量级较小,滑动迹象不明显,仅由基于InSAR相位梯度叠加的滑坡隐患快速识别方法成功识别出来。这些对比结果证实了本文所提方法在滑坡隐患识别方面具有良好的表现。

4 讨论

根据图4和图5分析,三种方法共同识别到20处滑坡隐患,这一高度一致的识别结果充分证明InSAR相位梯度叠加的滑坡隐患快速识别方法的可靠性。然而,不同方法的识别结果也存在一定差异。下面将该方法和传统时序InSAR识别结果进行了分类,以进一步探究本文方法的特点,我们将识别出的滑坡隐患分为两大类(表2):

表2 识别的滑坡隐患分类

第一类是所有方法均能识别到的滑坡隐患,三种方法共同识别到20处滑坡隐患,图6中放大了其中几处典型区域,它们反映了第一类滑坡隐患的一般情况。由图6可知,这几处滑坡识别的形变速率达到60 mm/年以上,形变面积均大于400 m2。这类滑坡隐患形变量级大,滑动迹象明显,本文方法和传统时序InSAR技术均能识别。

( a、b、 c为InSAR相位梯度叠加快速识别结果;d、 e、 f为Stacking-InSAR结果;g、 h、 i为SBAS-InSAR结果)图6 三种技术均能识别的滑坡隐患典型结果Fig.6 Typical results of potential landslides that can be identified by the three methods

第二类滑坡隐患是仅由InSAR相位梯度叠加的滑坡隐患快速识别方法探测的滑坡隐患。如图7所示(白色区域为InSAR相位梯度叠加的滑坡隐患快速识别结果),由图7(a3、b3、c3)可知这几处滑坡隐患形变速率均在30 mm/年以下,形变速率较小。同时这几处滑坡隐患形变范围均在300 m2以下,形变的范围较小,滑动迹象不明显。因此传统Stacking-InSAR难以探测到[14]。这几处滑坡隐患却被InSAR相位梯度叠加的滑坡隐患快速识别方法明显探测到。其原因是InSAR相位梯度对形变具有高度敏感性,即使微小的形变也能引起相位梯度的变化。因此对于形变量级较小的滑坡隐患,基于InSAR相位梯度叠加的滑坡隐患快速识别方法有更好的识别效果。

值得注意的是基于InSAR相位梯度叠加的滑坡隐患快速识别方法,可以在不进行相位解缠和时序建模解算的情况下直接求取InSAR差分干涉图的相位梯度。这节省了大量数据处理的时间,提高了滑坡隐患识别效率。如表3所列,统计了三种技术从数据配准到识别出滑坡隐患所需要的时间以及识别的滑坡隐患个数。本次实验统一采用CPU为Inter i7-11700;GPU为RTX2060,32 G内存;4 T机械硬盘,传输速度100 M/s;在常温24°环境下进行。本此研究区面积约为373 km2,共使用了97景升轨哨兵影像。三种方法均能有效识别库区内大多数滑坡隐患。基于InSAR相位梯度叠加的滑坡隐患快速识别方法识别出23处滑坡隐患,用时2h50min,耗时最短,识别速度相较Stacking-InSAR、SBAS-InSAR分别提升了1.4和1.9倍(表3)。与传统时序InSAR相比,基于InSAR相位梯度叠的滑坡隐患快速识别方法具有耗时短、操作简单、精度高等特点。可以有效缩短数据处理的时间,提高滑坡隐患的识别效率,在广域滑坡隐患快速动态识别中具有一定的应用潜力。