近年来,全国发生多起重大地质灾害,主要类型包括滑坡、泥石流和崩塌等。受地理环境、气候变化和人类活动等多方面因素影响,地质灾害在全国分布具有明显的地区性,南部多雨地区因其地质结构复杂、季节性降雨量大、山区丘陵地貌占比大等,属于地质灾害高发区。中华人民共和国自然资源部印发的《2024年全国地质灾害防治工作要点》中明确指出,需重点关注浙江等地的地质灾害防范,优化完善并持续开展全国地质灾害高、中易发区重大隐患综合遥感识别,逐步从发现隐患向重要区域和重大隐患动态监测转变^[1]。

传统的滑坡等地质灾害识别与监测多采用人工巡查与水准测量、全球导航卫星系统(Global Navigation Satellite System, GNSS)等相结合的方式进行,该方法一般适用于已知滑坡且人力可达的局部区域,对于大范围区域尺度的定期监测而言,其人力物力成本高^[2]。以近景摄影测量(包括无人机航空遥感)、高分辨率光学影像变化监测、三维激光雷达扫描为代表的遥感技术,能获得滑坡体地表覆盖和三维信息,在滑坡等地质灾害监测中的应用日益普遍,但是容易受到山区复杂地形环境、天气以及植被覆盖等因素影响^[3-5]。此外,自动化倾斜仪、变形监测仪等自动化监测设备能够有效监测滑坡体的动态变化^[6],但其仅适用于监测高风险滑坡单体,且在地形陡峭山区仪器布设难度大。

近年来,美国地质调查局、我国自然资源部、欧洲地质调查机构等多国机构积极推动结合多种先进技术手段的综合遥感地质灾害监测方法,通过将干涉雷达(interferometric synthetic aperture radar, InSAR),无人机激光雷达(light detection and ranging, LiDAR),GNSS、地震监测和自动化设备等多种监测手段结合,形成全方位、多层次的地质灾害监测网络^{[1,6 -8]}。我国早在2017年已提出并倡导了面向重大地灾隐患早期识别,由高精度遥感结合InSAR的“普查”、机载LiDAR结合无人机航拍的“详查”、地面调查核实的“核查”共同组成的“空-天-地”一体化的“三查”体系,以此来突破传统灾害隐患排查的极限和瓶颈^[9]。目前,“三查”体系与综合遥感的应用已在个别地区的灾害调查和灾后应急中发挥了显著作用,但面对大范围、规模化的监测需求,回答从形变到致灾的问题,仍需进一步深入研究^[10]。

以浙江省宁波地区为代表的华东地区,其地质灾害点呈散落状分布,从近年发生的滑坡分析发现,新生地质灾害具有高位、隐蔽性、人不能至等特点,仅靠传统监测方法无法全面识别和发现地灾滑坡隐患。综合遥感技术以其宏观、快速等优势可弥补上述方法的不足,特别是InSAR技术,可以全天时全天候地对地表进行大范围长时间的观测,获取厘米级乃至毫米级精度地表形变信息,可作为常规滑坡监测方法的重要补充^[11]。本研究以浙江省宁波地区作为试验区,结合InSAR、高分光学影像与无人机LiDAR的综合遥感技术在解决地质灾害形态、形变和形势的调查与判断方面的优势,对成灾状况、当前变形状况、潜在发展趋势与致灾形势进行评估,实现宁波全区域地质灾害的高精度、多层次、周期性的科学监测管理。

浙江省宁波市地处华东沿海,呈低山丘陵地貌(图1),属亚热带季风气候,降水具有明显的双峰特征,主要集中在3—7月的“春雨连梅雨”季节及8—9月的“台汛”期,具有雨量集中、降水强度大、持续时间长等特点。境内地貌类型多样,受不同时期地质构造运动和海侵、海退影响,地质环境条件总体较为复杂。山地约占陆域面积的54.7%,为火山岩大面积覆盖,沟谷纵横、岩石较为破碎,或上覆松散堆积层较厚,梅汛期中热带气旋、强对流等产生的暴雨容易引发滑坡、崩塌、泥石流等地质灾害。滑坡大多发育于低山丘陵、山间盆地边缘,滑坡体主要由松散残坡积层和基岩风化层组成,大中型滑坡基本分布于区内山势较大的地段,以土质滑坡为主,多发生在土层厚度大(一般大于3 m)的陡坡地,单纯岩质滑坡少见,人为削坡和汛期强降雨是引发滑坡的主要因素。

截至2020年底,全市地质灾害隐患点共34处,其中滑坡12处、崩塌10处,泥石流12处,共威胁人口735户1 789人; 地质灾害风险防范区825处,共影响人口3 779户9 595人; 中心城区累计沉降量大于100 mm的沉降面积169.4 km²,2020年平均沉降速率9.9 mm/a^[12]。针对宁波地质灾害的成因机理、影响因素、沉降特征以及沉降监测,以往已经开展了较多工作,地质灾害防控成效显著,但考虑到未来城市地质安全风险防控需求,仍需要持续开展山区等地带的地质灾害早期识别监测,掌握地质灾害沉降时空分布特征,预测未来的发展趋势。

本研究采用覆盖宁波全区的31景Radarsat-2数据(5 m分辨率,C波段)进行CMT-InSAR形变反演,数据详细信息如表1所示。

本研究所使用高分辨率光学影像为2018—2021年采集获取的宁波地区0.5 m分辨率高分航空遥感影像,以及2021年原始高分-2卫星遥感影像(多光谱); 所使用DEM数据为分别基于2020与2021年获取的高分辨率航空遥感影像制作生成的分辨率为2 m的DEM数据。由2 m的DEM渲染得到坡体单元影像,并制作生成山体坡度影像。

本研究基于CMT-InSAR获取的形变场,结合坡度数据,按照“年形变率≤-5 mm/a”和“山区坡度≥15°”且靠近居民地房屋等建筑物的原则进行疑似灾害点筛选,兼顾形变点空间分布选取适量疑似灾害点开展野外验证核查。结合以上数据分析及野外核查情况,选取了2处重点关注坡体,开展了无人机LiDAR数据采集、处理与综合分析研究。

本研究采用的综合遥感滑坡灾害隐患点识别作业模式为: 高分影像获取承灾体地形→InSAR形变反演提取潜在风险区域→数字高程模型(digital elevation model,DEM)获取坡体单元、坡度图解译潜在风险区域→机载LiDAR辅助精细化识别等高精度地理信息内外业核查。实现了由面到点的聚焦,由表征到机理的探索,将有助于提高面域范围上地质灾害早期识别的效率,减少人工外业核查的工作量和安全风险,进一步大大减少物力、财力投入。整体流程见图2。

本研究采用基于永久散射体(persistent scatterer,PS)与分布式散射体(distributed scatterer, DS)联合构网的CMT-InSAR方法,进行地表形变反演。此方法集成了主要针对城区人工强散射体PS和主要针对自然地表的分布式散射体DS这2种相干目标点的提取方法^[13-15],可有效增加山区丘陵地貌条件下的高相干点密度和空间覆盖度,提高形变监测的精准度。此外,CMT-InSAR 方法通过联合PS和DS散射体的相干信息,在构建干涉网时增加了冗余观测,以提高时序相位解算的稳定性; 同时,该方法利用DS点的统计分布特性补充了传统PS-InSAR方法中因相干目标不足导致的空间盲区问题,从理论上有效降低了误差累积的风险。利用这些高相干散射体上的相位信息,来反演地表形变场并获取形变时空过程信息。相对于传统MT-InSAR方法,CMT-InSAR方法对植被覆盖区域的丘陵地表具有更好的适用性。此外,CMT-InSAR方法通过对每个DS对应的相干矩阵进行处理得到一个最优相位值向量,将PS和DS联合处理而无需对经典PS-InSAR流程做重大改变。

将时序SAR影像进行数据解压和标准格式转换后,依据轨道与成像信息生成时空基线分布图。依据数据采集时天气晴好、影像尽可能位于时间和空间基线分布中心等原则选取主影像,然后将所有影像进行配准和裁剪。在数据预处理完成后,通过干涉图生成、多视处理、去平地效应、去地形相位等步骤,生成差分干涉图。

本研究试验区所在的宁波山区地表大都有植被覆盖,PS点在此类研究区域普遍偏少,通过同质像元滤波和相位优化可提取部分DS点,从而尽可能增加测量点空间分布密度。PS候选点的提取采用经典的幅度阈值法进行。对于DS候选点,为准确估计分布式目标的相位与相干性,本研究采用自适应的非规则窗口识别在地表均质分布的分布式散射目标,进行稳健估计相干矩阵,减小窗口内系统相位带来的信号非平稳性。该方法同时可减小多视处理降低图像几何分辨率的影响,以及减少矩阵估计不准的影响。

CMT-InSAR形变参数反演核心流程见图3,本研究采用迭代算法进行参数估计和优化。

1)算法初始化。去除地形相位时间域解缠(参数求解)残余相位残余干涉图。

2)迭代过程。残余干涉图滤波空间域解缠非缠绕残余相位滤波低频相位(大尺度大气、基线误差)去除地形残余相位与低频相位时间域解缠(参数再求解)残余相位残余干涉图。

3)终止迭代。残余地形估计趋于0,模型拟合残差小,残余干涉图平滑无突变。

该方法的主要思路是采用高相干边增强相位解缠的可靠性,以及模型拟合与参数估计的精度; 通过迭代去除残余地形相位,降低残余相位的缠绕度,增强空间域解缠的可靠性; 去除低频相位,把模型不能拟合的噪声相位去除,提高模型拟合度。

采用CMT-InSAR获取形变场,并提取出中高风险的地表异常形变区作为疑似灾害点。在此基础上,通过对比不同时期的高分辨率光学影像,解译地表覆被变化来辅助分析引起地表异常形变的影响因素。此外,在高分影像上提取试验区居民地、道路等承灾体数据,从而分析疑似灾害点与居民地、道路等重要设施之间的相对位置关系,评估潜在风险对社区和基础设施的影响。

基于DEM数据生成坡度、坡向、高程等要素,结合高分辨率遥感影像及地形地质条件,综合评估易致灾区域。地质灾害的发生与坡体形态、山体坡度和地表植被覆盖情况均具有一定的相关性,尤其斜坡单元,是滑坡、崩塌、泥石流等地质灾害发育的基本单元。如图4所示,在已有的地质灾害收集中,滑坡多发生在圈椅状地形(a)、大肚子斜坡(b)、扇面状(c)、双沟同源(d)等坡体形态,泥石流多发生在漏斗状(e)、长条形(f)、树枝形(g)等坡体形态; 崩塌多发生在坡度大于30°(h)以上的坡体形态。因此,选择坡体的形态作为地质灾害判识的一个重要的因素。

坡体单元的提取可通过对DEM产品图进行人工解译获取。首先对 DEM 进行重采样,利用 ArcGIS工具生成山体阴影产品,根据山体坡度的形态解译出单元坡体。

此外,在降雨期,坡面上的积水易汇聚于冲沟区域,从而为地质灾害的发生提供了诱发条件,因此判断坡面上是否存在冲沟也是判识坡体在强降雨过程中是否稳定的一个重要因素。本研究通过遥感影像纹理与山体阴影图相结合来提取冲沟区域。

基于CMT-InSAR形变监测获取的地表异常形变区,结合高分光学影像对比分析地表覆被,以及DEM数据坡体、坡度分析结果,经野外实地核查疑似灾害点,选取2处典型坡体开展实地无人机LiDAR数据采集与处理,将所采集的点云数据进行滤波处理分离出地面点,生成林下高精度DEM,进而开展坡体单元分析,辅助分析滑坡风险。

通过CMT-InSAR处理得到的宁波地表形变速率图(图5)可见,宁波地区整体地表较为平稳,局部沿海区域因填海建设等活动导致地表形变较大,沉降速率超过-20 mm/a。在平原地区,江北区、海曙区东北部、奉化区东北部、镇海区东南部等区域地表形变较为显著,部分区域地表形变年均速率超过-20 mm/a。在山区,宁海县、奉化区中西部、余姚市、北仑区北部等区域地表形变较为显著,部分区域地表年均形变速率达到-20~-5 mm/a。

在进行SAR数据处理时,通常影像数量越多、相干性越高、方差越小,可达到的线性形变速率反演精度越高^[16]。本研究采用GAMMA软件结合二次开发的方式实现CMT-InSAR算法,在通过长时序InSAR构建的相位模型解算形变速率参数时,模型会同时估计平均形变速率的反演误差,用于参考判断平均形变速率反演的精度。以整体平均相关系数为0.7考虑(实际植被覆盖山区地貌的整体相干系统更低),影像为30景时,方差约0.9 mm/a^[16]。图6为GAMMA软件计算平均形变速率反演误差统计直方图截图。由图6可见,其标准差为±0.176 mm/a,满足小于2倍的方差(±1.8 mm/a)的精度要求。

本研究重点关注宁波南部山区的地质灾害隐患点分布情况,因此本研究利用最新高分影像,按照1∶2 000比例尺精度要求,获取了2022年宁波市南部山区居民地、道路等承灾体数据,如图7所示。

通过 ArcGIS 软件对数字高程模型处理生成坡度、坡向、高程等要素。结合数字高程模型、遥感影像及地形条件,分析斜坡的高度、坡度、坡向等特征,判断易致灾体。

1)疑似地质灾害隐患点位筛选。本研究重点关注宁波山区地表形变导致的地质灾害隐患问题,因此采用年均形变速率和地形坡度这2个约束条件来提取山区地表的形变风险点。基于CMT-InSAR形变反演得到的8 345 557个形变点,引入2 m分辨率DEM数据,按照“年平均形变率≤-5 mm/a”“山区坡度≥15°”2个阈值,共提取7 749个候选风险点位。在此基础上,结合居民地、道路等承灾体数据,重点聚焦那些易造成人员伤亡及经济损失的风险区,剔除掉远离村落、工厂、基础设施等人类活动区域,得到形变反演疑似地质灾害隐患点1 092个(图8), 为野外开展巡查提供了依据。

2)高风险疑似地质灾害点位提取。本研究采用2 m的DEM分析渲染得到直观的坡体单元形态,以及山体坡度值,并将两者叠加,形成坡度与坡体单元结合图(图9)。参照图5所示的地质灾害坡体形态进行识别,并结合坡度情况(≥30°)和0.5 m高分影像的植被覆盖情况(图10),及邻近居民设施等条件,解译得到灾害隐患区域位置,共识别出高风险疑似点位142个(图11)。

3)野外实地核查。依据以上提取142处高风险地质灾害隐患点,兼顾形变点空间分布,本试验共选取22处点位开展了野外验证核查,经核验所选点位均存在已发生过的滑坡痕迹(InSAR数据时间覆盖范围为2020—2022年),部分点位已经进行了施工加固,如图12至15所示。

根据图13所示的时间序列相对形变量线状图可以看到,2020—2022年期间,该处持续存在沉降形变,且形变趋势呈现线性特征,年平均形变速率超过7 mm/a。从坡度与坡体单元结合图(图15)看,该处坡度陡峭,为57°(大于30°),该处为水坝附近的护坡,野外踏勘发现沿线道路存在明显裂纹,边坡呈现明显危岩体,已采取防护网等加固措施(图14)。

4)与无人机LiDAR数据采集分析。结合以上数据分析及野外核查情况,选取了位于宁海县陈家坑村和象山县王岙村的2处重点关注坡体,开展了无人机LiDAR数据采集、处理。

① 宁海县陈家坑村。根据CMT-InSAR解算结果,发现形变量较大的点位靠近疑似滑坡处(图16)。进一步结合0.5 m高分影像,发现陈家坑村西侧山坡,相比2018年,2019年开始有2处植被明显稀疏,疑似出现滑坡,至2021年植被明显增多(图17)。

对此,采用无人机进行了LiDAR数据采集,获取高精度的DEM和坡度等数据。结合点云渲染图和坡体、坡度图可以看出(图18),在影像上疑似滑坡的位置存在2处明显的冲沟,坡度在35°左右,且形变点位于南北2处冲沟中间的坡体上,属于双沟同源型,易造成滑坡。

② 象山县王岙村。结合CMT-InSAR解算结果(图19)及外业核查照片(图20)分析,此处属于疑似高位形变点,对此,进一步采用无人机进行了LiDAR数据采集与处理。

结合点云渲染图DEM图与坡体、坡度图(图21)可见,此处坡体冲沟显著,多处坡体属双沟同源型,整体呈扇面状,且形变点位均为冲沟处,坡度为40°,易发生滑坡。

总体看,采用无人机激光点云数据采集与处理,进一步精细的构建了DEM和坡度图,冲沟、坡面更加显著,由此可见采用无人机LiDAR技术,获取林下地形特征,将有助于提升对大面积风险点核查的效率,同时减少人工外业实地核查的安全风险,实现技防替代人工。

本研究将InSAR、高分辨率光学影像、无人机LiDAR等多种遥感手段相结合,对宁波地区进行综合遥感滑坡地质灾害监测,相关研究成果总结如下:

1)形成了一套多层次的综合遥感地质灾害监测方法。利用InSAR监测获取斜坡等地质体地表变形状态,判别灾害体的滑移规律、活动阶段和发展趋势; 利用高分光学影像以及在此基础上制作的DEM数据,分析坡体地表覆被变化与人为活动干扰因素,结合坡体形态与坡度坡向判别坡体单元类型; 利用无人机LiDAR对典型高风险坡体开展高精度地形测绘,提取灾害体林下三维形态细节信息。最后,综合InSAR反演的形变场、地面实地调研核查信息、宁波地质背景资料,获得宁波潜地质灾害隐患点分布信息,并提取典型高风险坡体的灾害隐患点详细位置与坡体形态细节信息。实现了由面到点的聚焦,由表征到机理的探索,将有助于提高面域范围上地质灾害早期识别的效率,减少人工外业核查的工作量和安全风险,进一步大大减少物力、财力投入。

2)通过本研究完成宁波疑似地质灾害隐患点提取。宁波整体地表平稳,从总体区域分布看,所提取的142处高风险点位主要集中在奉化区、宁海县、余姚市和象山县。建议下一步地灾防控防治工作有针对性地开展加固施工,并加强以上区域台风、暴雨天气前后的检查。

3)本研究采用的联合PS与DS的CMT-InSAR方法,相比于传统PS-InSAR方法,提取的形变点数量增加了136%,改善了形变点的空间覆盖度与密度,进而提高了形变监测的精准度。同时,将特征分解作为迭代初值进行组合求解,提高了DS点提取效率,进一步丰富了InSAR用于地质灾害监测反演的技术方法。

在取得以上成果的同时,本研究所反演的InSAR形变场在宁波南部植被密集覆盖的山区地貌,仍存在局部范围的空洞和数据稀少问题,针对此问题,后期研究将着重从2个方面进一步加强:

1)选择L波段的时序SAR数据作为InSAR形变场反演的数据源。从卫星SAR传感器来说,使用波长较长的SAR数据穿透性会更好,可以获取植被层下更加丰富的信息,有利于相干性的准确估计。已有相关研究表明,L波段比C波段具有更强的穿透性,能够更好地抵御失相干的影响,获得更高质量的形变结果。目前,在轨运行的L波段卫星中日本ALOS-2卫星、阿根廷SAOCOM卫星等国外卫星价格昂贵,对宁波大范围开展周期性监测成本较高。我国陆地探测一号(LT-1) 卫星数据的普及与共享将为开展植被区长波InSAR滑坡形变监测提供充足的数据支撑,目前已开启数据收集工作,在后续研究中将考虑使用LT-1卫星数据开展形变监测。

2)宁波山区地表植被覆盖密集,局部山区难以获得足够数量的相干点。针对重点关注的高风险坡体,通过布设人工角反射器可有效增加高相干点数量。目前已在部分重点坡体开展三角反射器布设工作,在后续研究中将引入通过三角反射器获得的高相干点参与形变反演,以提高植被覆盖密集区域的形变点密度和形变反演精度。