基于InSAR的黄土高原延河流域滑坡降雨阈值模型构建

01 摘要

由图1可知，延河流域（108^°39^'30^″—110^°28^'46^″E，36^°20^'39^″—37^°19^'18^″N）位于我国西北内陆半湿润向半干旱过渡地区，属于我国黄土丘陵沟壑区第二副区，流域面积7 725 km²。延河是黄河一级支流，全长286.9 km。延河流域地势西北高东南低，多年平均气温为9.2 ℃，降水量为500 mm，土壤类型主要为黄绵土，流域水土流失严重，河流含沙量大，多年平均输沙量为7.16×10⁷ t，其黄土结构疏松，厚度约为几十米至上百米，最厚处达150 m左右。延河流域地形错综复杂，地势起伏较大，河流和沟壑分布密集，岩石和土壤受到强烈的侵蚀作用，长期受到季节性热浪和降雨的影响，导致该区域滑坡、崩塌、不稳定斜坡等重力地质灾害频发。

本文主要采用的数据包括哨兵一号（Sentinel-1A） 影像、DEM和降雨数据。其中Sentinel-1A数据基于时序InSAR计算延河流域的地表形变速率和地表形 变时间序列；DEM数据用于模拟地形相位和生成坡度等地形因子；降雨数据用于分析降雨与滑坡形变的关系从而计算降雨阈值模型。

（1）SAR数据。Sentinel-1是由欧洲空间局（European Space Agency，ESA）发射运营的两颗对地观测卫星组成，工作波段为C波段，波长约为5.6 cm，每颗卫星重访周期为12天，渐进式扫描地形观测（terrain observation by progressive scans，TOPS）模式下单景覆盖范围达到250 km×170 km，可提供全天时、全天候的高分辨雷达卫星数据。由于覆盖本研究区最早的一景Sentinel-1A影像时间为2017年3月12日，因此，本研究选用时间覆盖范围为2017年3月12日至2021年12月4日的Sentinel-1A升轨影像数据，共115景。由于原始Sentinel-1A影像覆盖范围较大，为提高数据处理效率，裁剪出与研究区相关的区域进行处理。其中在InSAR数据处理中，为减小基线误差，本文还使用Sentinel-1A影像对应的POD精密轨道卫星星历（Precise Orbit Ephemerides，POD）数据，数据详细参数见网址https：//slqc.asf.alaska.edu/aux_poeorb/。

（2）DEM数据。本文采用美国航空航天局（NASA）和国防部国家测绘局（NIMA）联合发布的30 m分辨率SRTM DEM数据，数据下载网址为https：//earthexplorer.usgs.gov/。

（3）降雨数据。本文采用延河流域2017—2021年5年的GPM降雨数据。全球卫星降水计划（GPM）是由NASA和JAXA合作开展的国际卫星任务，利用多传感器多卫星多算法结合卫星网络和雨量计反演更高精度的降水数据，时间分辨率为1天。

干涉图堆叠技术（Stacking）最初由SANDWELL等于1998年提出，其基本思想是：假设地表形变为线性形变，而大气延迟在时间域上为随机信号，当对N幅干涉图进行相位平均时，线性的地表形变信号增强了N倍，而随机分布的大气延迟只增加了N倍，从而提高形变信号的信噪比。其数学模型为：

由图2可知，基于InSAR的滑坡识别主要分为地理编码、形变速率计算和滑坡识别3部分。具体流程为：

（1）选取主影像，将外部DEM与主影像配准，获取高程相位；

（2）首先将主影像与副影像配准，为了抑制噪声影响、提高干涉图的相干性，在距离向和方位向采用10×2的多视比，其次设置空间基线75 m、时间基线85天组成干涉对，基于流程（1）获取的高程相位，对干涉对进行差分干涉、滤波、相位解缠及大气改正处理。最后对447个干涉对进行筛选，筛选出相干性较好、解缠误差较少的305个高质量解缠图进行Stacking解算，获取沿雷达视线（LOS）向的形变速率；

（3）首先结合研究区地形和地质灾害的活动特 性，选择人工建筑物作为稳定区，并计算标准差和最 大形变速率，最终选定形变速率阈值为10 mm/a、斜坡坡度阈值为5°进行滑坡的初步识别；其次将初步识别结果导入谷歌地球中，与多时相光学遥感影像进行叠加，分析疑似滑坡对应的地形地貌特征进行滑坡最终识别，包括滑坡位置、范围等。

在重力作用下，滑坡体通常具有3阶段演化特 征：初始变形阶段、等速变形阶段和加速变形阶段。当滑坡变形发展到一定阶段后，变形速率会呈现出不断加速增长的趋势，直至坡体整体失稳（滑坡）之前，变形曲线近于陡立，这一阶段被称为加速变形阶段。为了确定滑坡加速变形的时间，对SAR数据采用小基线集（Small base line subset，SBAS）InSAR技术进行处理。其基本原理为：使用时空基线较短的影像对任意组成小基线集，并使用最小二乘法求解集合内地表形变，最后使用奇异值分解法将多个独立的小集合联合求解，以获取地表形变时间序列。

由图3可知，地表形变时间序列中，时间节点对应每一景影像获取时间，由于相邻2个景影像的获取时间间隔存在差异，因此本文计算了相邻2个景影像间滑坡点的单位时间内形变量。其单位形变量最大值可设定为滑坡加速变形时刻，该值对应影像的开始时间即为滑坡加速变形时间。图中，比较曲线AB、BC 和CD，BC单位时间形变量最大，故点B对应的时间即为滑坡加速变形时间。由于数据局限性，本文将滑坡加速变形时间近似估计为滑坡发生时间用于后续构建滑坡降雨阈值模型。

降雨对诱发滑坡的影响存在一定的滞后性，随着降雨时间的前移，前期降雨的影响逐渐减弱。本文采用CROZIER有效降雨量模型分析滑坡与前置降雨量的关系，其计算公式为：

本模型中有2个关键因子，即前期有效降雨天数n和有效降雨衰减系数K。本文引入累积降雨量的概念，累积降雨量指滑坡发生当日及前一段时间内的降雨总量，其中将滑坡发生当日的降雨量记为R₀，滑坡发生当日及前n天内的降雨总量记为R_0n。采用皮尔逊相关系数r计算累积降雨量与滑坡形变量间的相关性，确定前期有效降雨天数n。其中，相关系数最大的累积降雨量对应的前期降雨天数设定为前期有效降雨天数n。皮尔逊相关系数r计算公式为：

对于前期有效降雨衰减系数K，本文分别设置不同的前期有效降雨衰减系数K值（从1到0.5依次减小），根据前期有效降雨天数n，计算出n天有效累积降雨量，然后分别计算不同衰减系数下，滑坡形变量与n天有效累积降雨量的相关系数r，最终确定合理的前期有效降雨衰减系数。

由表1、图4可知，本文主要通过前期有效降雨量E、降雨强度I和降雨历时D，建立现阶段主流的3种经验性降雨关系阈值模型E-D、I-D和E-I，其中，前期有效降雨量E根据公式（2）确定；降雨历时D基于李长江等模型，即将降雨事件定义为从日降雨量大于等于2 mm作为降雨事件的开始时间，直到连续2天日降雨量小于2 mm结束的时间段确定为降雨历时D；降雨强度I，短历时取峰值为降雨强度，长历时取平均值为降雨强度。

03 结果与分析

（1）基于时序InSAR技术，对115景Sentinel-1A影像进行数据处理，获取2017—2021年覆盖研究区LOS向的形变速率；选定形变速率阈值为10 mm/a、斜坡坡度阈值为5°进行滑坡的初步识别，并与多时相光学遥感影像进行叠加，分析疑似滑坡对应的地形地貌特征进行滑坡最终识别；采用SBAS方法计算滑坡时间序列形变，确定滑坡的加速变形时间。

（2）计算滑坡形变量与累积降雨量的相关性，结果表明：随着前期降雨天数的增加，累积降雨量与滑坡形变量的相关性呈现先升后降的趋势；当前期有效降雨衰减系数K=1时，9天有效累积降雨量（即滑 坡发生当日及前8日降雨量）与滑坡形变量的相关性最强。

（3）计算延河流域滑坡点的前期有效降雨量E、降雨强度I和降雨历时D，并结合幂函数建立E-D、I-D和E-I降雨阈值模型，对应的阈值模型表达式分别为E=12.19D^0.661、I=10.66D^-0.255和E=3.84I^1.141。