Stacking,InSAR和IPTA技术应用于探测宁夏隆德县滑坡隐患

王 川，涂 宽，谌 华，耿 丹，王文龙，李樵民

（1.二十一世纪空间技术应用股份有限公司，北京 100096；
2.宁夏回族自治区遥感调查院，宁夏银川 750021；
3.高分辨率对地观测系统宁夏数据与应用中心，宁夏银川 750021）

滑坡是在重力、降雨、地震等因素的作用下，坡体作剪切运动的自然现象，通常可划分为推移式滑坡、土流型滑坡、土质滑坡和古崩滑坡堆积体等［1-2］。其产生条件包含内因和外因两类，内因为地质构造和地貌特征，外因则为降雨、地震和人类工程活动等。近年来由于全球气候波动，异常天气条件增多，导致更多滑坡灾害在各地发生［2-4］。

由于滑坡灾害的发生严重威胁了人类的生命和财产安全，对人类社会环境、经济环境、生态环境造成了巨大破坏，因此需要对滑坡隐患进行及时的探测和监测。通常，对滑坡进行探测和监测的手段有：常规三角网测量、伸缩计深部位移监测、应力应变测量、GPS遥测监测网、地震监测、近景摄影测量以及多种技术的联合等［2，5］。上述探测和监测手段的准确性高，测量结果可靠，但测量主要集中在小区域，不利于开展大规模滑坡隐患探测。

近年来，随着科技的不断进步，出现了以测量可见光光谱、热辐射、微波辐射等为目标的多种遥感技术［6-7］。特别是合成孔径雷达（Synthetic Aperture Radar，SAR）技术的出现。微波遥感干涉合成孔径雷达（Interferometry Synthetic Aperture Radar，InSAR）测量地表形变技术测量结果准确度高，覆盖面大，在滑坡隐患探测中得到了普遍应用［8］。最先发展起来的应用电磁波相位值对地表起伏进行量测的是差分干涉技术（Differential Interferometry Synthetic Aperture Radar，DInSAR），但是单独使用两景影像进行干涉易产生各种去相干（时间、空间、体）等，且易受大气、系统热噪声、处理误差等影响，因此目前主流技术是多个干涉对共同测 量 的 时 序InSAR技 术（Multi-Temporal Interferometry Synthetic Aperture Radar，MT-InSAR）［2，9］。在MT-InSAR中，典型的是：（1）永久散射体技术（Persistent Scatter InSAR，PS-InSAR），该方法只选择回波信号稳定的点（永久散射体）进行分析，可以摆脱时空基线的限制，相干性保持好；
（2）小基线集技术（Small Baseline Subsets InSAR，SBAS-InSAR），该技术选择分布式散射体（如沙地、草地等）进行分析，更适用于自然环境［10-12］。以此两种思想为核心，产生了相干点目标分析（Interferometry Point IPTA）、SqueeSAR等技术，同时也出现了以计算形变速率为目标的干涉叠加技术（Stacking）和不依赖时序相位特征测量地表形变的多孔径雷达技术（Multiple Aperture Interferometry，MAI）［13-19］。InSAR技术不可避免的会受到植被、地形起伏和大气等的影响，使用单一技术往往只能克服部分缺陷而导致不能更全面的识别滑坡隐患，同时基于光学影像对滑坡隐患进行专家判读缺少坡体滑移的数据支持，部分老（古）滑坡已趋于稳定。因此需要联合多种InSAR方法并结合光学影像进行滑坡隐患探测。

IPTA技术针对点目标进行回归分析以分离相位（大气、地形、非线性形变等），计算效率和精度显著提高，而Stacking技术利用加权平均以抑制噪声使其能够快速获取大范围的形变信息，在现有各类InSAR技术中，此2种技术可操作性最强，同时结果可互相补充。隆德县境内以黄土丘陵区为主，较少出现大型滑坡隐患，小型滑坡隐患分布广泛，利用Stacking技术能够实现对大部分滑坡隐患的探测，基于IPTA技术能够去除各种误差，探测滑坡隐患的可靠性高。

因此，本研究利用Sentinel-1卫星影像，采用IPTA和Stacking这2种InSAR技术对隆德县全境滑坡隐患进行探测，分析和评价2种InSAR技术的探测结果。该项工作对隆德县滑坡隐患普查及监测具有十分重要的指导意义。

研究区隆德县位于宁夏回族自治区的南部，六盘山西侧，全县总面积992.39 km2。在地形地貌特征上，隆德县处于昆仑秦岭地槽褶皱区东段，曾受到燕山期、喜马拉雅期构造运动的影响，其主要地貌类型有：红层岩丘陵地貌、中低山地貌、河谷平原地貌、黄土丘陵地貌。在气候特征上，隆德县属于温带大陆性季风气候，受东高西低的地势差异影响，致使东部湿润寒冷而西部干燥温暖，冬季干燥而夏季多雨。每年7～9月降雨最为集中，约占全年降水量的60%。降雨和降雪是隆德县地质灾害发生的重要诱因。由于人类工程活动的影响，如毁林开荒、削坡修路、兴修水利和过度放牧等，破坏了生态系统平衡，加剧了地质环境恶化，更进一步的促使了滑坡灾害的发生。

隆德县滑坡包含黄土滑坡和岩质滑坡2种，其中黄土滑坡包含黄土层内滑坡和黄土—基岩滑坡，岩质滑坡主要是新近系泥岩残坡积物滑坡。滑坡规模通常为小型滑坡，中、大型滑坡发育较少。滑坡和滑坡隐患多分布于渝河、南川河的冲蚀岸边及其支流或支沟两侧，呈条带状或树枝状分布，在地貌类型上多集中在丘陵地貌区［20］。滑坡以浅层和中层为主，厚层滑坡较少。

Sentinel-1卫星影像幅宽广（IW成像模式可达250 km），重访周期相对较短且固定（单星重访周期12 d），因此文中研究选择2019年1月至2020年5月在重访周期过境的全部卫星影像用于测量地表形变，表1详细列出了所使用的卫星影像参数信息。图1为研究区概况图和Sentinel-1卫星影像覆盖范围。此外，为了纠正地形影响和提取地形信息，文中研究还同时获取了SRTM 30 m的DEM作为辅助数据。

表1 Sentinel-1卫星影像参数表Table 1 Information of Sentinel-1

注：在Sentinel-1假彩色合成图中，绿色区域反映出地物主要由体散射比较强烈的散射体如植被构成，深粉色区域反映出地物主要由粗糙的表面如裸地构成，浅粉色区域表示地物由去极化能力非常弱且后向散射强烈的建筑物组成，水体由于镜面散射而呈现黑色。

本次实验所获取的SAR数据能够完整覆盖2019年全年，时间和空间基线满足形变测量需求，配准精度达到InSAR干涉测量要求（1/8像元精度），DEM数据空间分辨率与SAR数据接近且已对空值进行插值，能够利用全部卫星影像和DEM数据开展后续研究工作。

2.1 InSAR成像原理

根据图2几何关系，P点的基线分解为：

同样，参考平面的点P0的基线分解为：

假设2次观测时没有地表起伏变化，同时忽略地物散射和辐射传输中产生的随机噪声，则P点的干涉相位可以视作2次传输路径之差所对应的相位：

对路径差ΔR求导，可知相位随路径的变化关系为：

对于星载SAR卫星，当2次观测视角的差异较小时，基于微波辐射传输的远场平行近似［21］，可以将ΔR等效替换为水平基线B∥，当观测距离不变且视角由θ0开始变化产生的干涉相位变化为：

若定义参考面点P0的干涉相位为参考相位，其与P点到传感器S1的距离相同，因此可以视其为P点的平地相位（排除地形起伏的影响，仅由观测视角θ0造成）：

根据图像的几何关系：

则视角由θ0开始变化产生的高程变化：

结合式（7）和式（10）可得：

图2 InSAR原理图Fig.2 InSAR schematic diagram

此处∂ϕif的变化仅与高程相关（P0为参考面起始观察点），即由P和P0点的干涉相位差分得到，所以称其为高程相位ϕh（忽略大气、形变等噪声相位）。

若P点在观测期间产生了位移（地表形变），则干涉相位可表示为：

式中：ϕd表示地表形变产生的相位，与之对应的雷达视线向的形变量为ΔRd，由地表形变相位与视线向形变量的关系可以看出，对于C波段SAR卫星，其地表形变的测量精度通常能达到毫米级。

2.2 研究方法

本研究选择Stacking和IPTA这2种较为成熟的InSAR方法用于测量隆德县境内的地表形变值。Stacking技术采用对干涉相位加权平均的方式计算相位的形变速率继而得到地表延视线向的形变速率，其处理方法相对简单且能够在一定程度上对大气相位等噪声进行抑制，因此能够快速的对大区域在长时间序列上得到地表形变速率的测量结果［22］。IPTA方法是一种基于点目标的干涉测量方法，通过迭代回归分析将大气相位、高程相位残余等进行分离，能够同时获取地表形变速率和地表非线性形变，与SBAS和PS等方法相比具有计算效率高、节省存储空间的特点［23］。

在进行InSAR处理前，需要首先对Sentinel-1卫星影像进行预处理，主要包括：读取原始数据并基于外部DEM数据对SLC图像进行配准，将配准后的图像裁剪至研究区范围，对DEM进行重投影至SAR坐标系。两种InSAR方法的主要处理流程如下：

2.2.1 Stacking InSAR技术流程

（1）计算各个像对间的垂直基线，设置阈值为-200 m～200 m，计算时间基线，设置阈值50 d，生成基线连接图和干涉对列表（图3）。

（2）根据干涉对列表对SLC图像进行差分干涉，对差分干涉图进行滤波以提高相干性，对干涉图进行解缠，对解缠图进行一定的空间滤波部分消除大气、轨道误差、高程误差等噪声相位。解缠方法：最小费用流法（MCF）［24］。

（3）基于加权平均的方法计算相位形变速率：

式中：φrate为相位形变速率；
Δtj表示第j个干涉相对的时间间隔；
N表示干涉像对数，然后按照λ/4π将相位形变速率换算为延视线向的地表形变速率（单位：mm/年）。由于单轨干涉测量存在着多种去相干（如时间、几何、体等），在解缠过程中需要对相干性较低的像元进行掩膜（文中研究相干性阈值为0.2），使得在进行Stacking时，各像素点所使用的样本数（干涉像对数）均不同，为了使结果更具可靠性，文中研究设定最小可接受样本数为108，即当像素点干涉像对数低于108则对结果进行掩膜。图4为Stacking InSAR总体处理技术流程图。

图3 Sentinel-1卫星影像基线连接图Fig.3 Sentinel-1 perpendicular base connection diagram

图4 Stacking InSAR处理技术流程图Fig.4 The flow chart of Stacking InSAR

2.2.2 IPTA技术流程

（1）分别基于像素点在时序上的稳定性和相干性挑选像元并生成点目标文件。计算两两像对间的垂直基线，设置阈值为-200～200 m，计算时间基线，设置阈值50 d，生成基线连接图和干涉对列表（图3）。根据干涉对列表对点目标进行差分干涉。

（2）在第j幅差分干涉图，点目标i其相位可以表示为：

图5 IPTA处理技术流程图Fig.5 The flow chart of IPTA

式中：Δϕ表示差分干涉相位；
vi表示第i个点目标的地表形变速率；
Δtj表示第j个干涉对的时间间隔；
δi,j表示第i个点目标在Δtj时间段内的非线性地表形变值；
B⊥j表示第j个干涉像对的垂直基线；
Δzi表示第i个点目标的高程误差；
λ表示波长；
ρi表示微波从雷达到达第i个点所对应的地物再返回的路径（斜距）；
θi表示入射角；
ϕatm表示大气相位；
n表示噪声相位。通过非线性回归分析分别求解出其中的差分大气相位和高程误差，对地表高程进行补偿并对差分大气相位进行去除后再次进行回归分析，进行多次迭代后可将差分大气相位逐步分离，并对地表高程进行充分补偿。

（3）通过回归分析获得解缠的差分干涉相位后，使用奇异值分解（SVD）的方法对多主影像的相位进行求解计算出延时间序列的相位形变［25］，对其在空间上和时间上进行滤波，并换算成地表沉降值（单位：mm）。在得到时序相位形变后，再次进行回归分析计算出相位的形变速率，同样换算为地表形变速率（单位：mm/年）。图5为IPTA总体处理技术流程图。

2.2.3 有效探测滑坡隐患和未有效探测滑坡隐患的定义及有效探测率

在开展文中研究工作之前，已先期对隆德县滑坡隐患进行了调查，此次工作联合多源遥感技术对隆德县全境滑坡隐患进行了普查，在此基础上对所识别出的滑坡隐患已完成了实地的调查和验证工作。该项工作由“宁夏科技厅重点研发计划项目”支持开展［26］。

根据InSAR地表形变测量结果识别滑坡隐患，若圈定的滑坡隐患与隆德县滑坡隐患遥感监测工作成果中的滑坡隐患相对应，则称其为该InSAR技术（Stacking和IPTA）的“有效探测滑坡隐患”，对于使用该InSAR测量方法未解译出但在隆德县滑坡隐患遥感监测工作成果中出现的其他滑坡隐患，则称其为“未有效探测滑坡隐患”。对于一种InSAR技术（IPTA或Stacking），其有效探测滑坡隐患占所有滑坡隐患的百分比文中研究称其为“有效探测率”。

通过与隆德县滑坡隐患调查结果进行核实，2种InSAR方法对47处滑坡隐患进行了有效探测［26］（图6，图7（a）、（b）），在形变测量中正值表示滑坡隐患在延视线向靠近卫星滑动（滑坡隐患和卫星间的距离缩短），负值表示滑坡隐患在延视线向远离卫星滑动（滑坡隐患和卫星间的距离增大）。表2详细列出了滑坡隐患的基本信息，其中滑坡隐患面积达到2.0×106m2的共6处，最大滑坡隐患面积6.2×106m2。

Stacking技术获得的形变结果中有效探测滑坡隐患44处，有效探测率93.6%，IPTA有效探测滑坡隐患共计37处，有效探测率78.7%，Stacking技术略高。2种InSAR技术共有效探测滑坡隐患34处，在相同数据源情况下监测重合度高，对最大形变速率进行相关性分析（图8），线性相关性强（R20.56），测量结果一致性强（例如图7（c）和（d）所示滑坡隐患），可相互印证2种InSAR技术测量结果可靠。

图6 隆德县境内滑坡隐患分布图Fig.6 Distribution map of potential landslide in Longde County

图7 隆德县Stacking和IPTA测量形变结果图及典型形变滑坡隐患Fig.7 The displacement map measured by Stacking and IPTA and the typical potential landslide in Longde County

表2 隆德县解译滑坡隐患基本信息表Table 2 Information of interpreting potential landslide in Longde County

续表

图8 不同InSAR方法测量滑坡隐患最大形变量线性相关性分析Fig.8 Correlation analysis of the maximum displacement for different InSAR methods

图9 有效探测滑坡隐患最大形变速率分布特征Fig.9 Characteristics of the maximum displacement of effectively detected potential landslides

3.1 Stacking和IPTA有效探测滑坡隐患结果分析

3.1.1 有效探测滑坡隐患的最大形变速率分布特征

图9描述了滑坡隐患最大形变速率分布特征。箱形图中间箱体上边沿表示上四分位数（3/4分位数），下边沿表示下四分位数（1/4分位数），上下四分位数之差为四分位距（IQR），箱体中间横线表示中位数，箱体中圆圈表示均值，由箱体延伸出上下边缘线（箱体外横线），上边缘线为上四分位数加1.5倍四分位距，下边缘线为下四分位数减1.5倍四分位距，超过上下边缘线的视为异常值。2种InSAR方法对不同形变速率滑坡隐患探测能力基本一致。其中一半的滑坡隐患最大形变速率集中在-20～20 mm/年之间（箱的上边沿和下边沿分别表示3/4和1/4分位数）。Stacking所识别滑坡隐患形变速率中位数（箱内横线）高于IPTA。

3.1.2 有效探测滑坡隐患的地形分布特征

研究主要针对地形中的坡度和坡向进行分析，其中坡度划分为平坡、缓坡、斜坡和陡坡4个等级（表3）［27］。坡向划分为北、东北、东、东南、南、西南、西和西北8个方向（表4）［27］。IPTA和Stacking有效探测滑坡隐患均集中在缓坡和斜坡上，对于坡度较大的平坡（1.7°～5°）和坡度较小的陡坡（25°～31.6°）亦有少量分布（图10（a）），这与InSAR的识别原理相符合，当坡度较小时滑移速率较慢无法满足InSAR测量精度，坡度过大则滑移速率过快易产生失相干导致无法测量。北坡方向未出现有效探测滑坡隐患，与北坡所受光照辐射少其微气候环境不易出现滑坡相关，在其他方向上未出现明显集中分布特征（图10（b）），因此综合2种InSAR方法能够较全面的对各个方向的滑坡隐患进行有效探测。

表4 地形-坡向统计结果Table 4 Statistics of topographic aspects

表3 地形-坡度统计结果Table 3 Statistics of topographic slopes

坡度级划分来源于“基于3PGS-MTCLIM模型模拟根河林区火后植被净初级生产力恢复及其影响因子”。

图10 滑坡隐患地形分布特征Fig.10 Characters of topography of potential landslides

3.2 滑坡隐患滑动特征分析

图11 IPTA有效探测滑坡隐患的时序形变量Fig.11 The displacement of potential landslides detected by IPTA method

IPTA技术不仅能够获取滑坡隐患的年平均形变速率，同时能够获取时序累计形变量，因此可以通过IPTA技术了解滑坡隐患的滑动状态。文中研究对IPTA获取的37处滑坡隐患的滑动过程做一基本描述（图11）。

LDhp1、LDhp3、LDhp4、LDhp5、LDhp9、LDhp13、LDhp26、LDhp28、LDhp37、LDhp44均以一定的滑动速率发生线性形变，滑动过程在整个监测期内较稳定，没有表现出受其他因素干扰而滑动过程受到扰动的情况。

LDhp2、LDhp6、LDhp8、LDhp11、LDhp12、LDhp18、LDhp22、LDhp30、LDhp32、LDhp43在前期滑动不明显或者滑动速率保持稳定，在经过一段时间后滑动速度突然加剧并且会维持至监测结束，这一加速过程发生在2019年6月～9月时间段内，这一时间段是一年中降水集中时期，因此推测是由于降雨冲刷坡体导致滑动过程加剧。

LDhp7、LDhp10、LDhp14、LDhp16、LDhp24、LDhp25、LDhp27、LDhp29、LDhp39同样在2019年6月～9月出先了较明显的加速过程，但随后在2019年10月～12月份加速过程突然减缓进入新的相对慢速滑动的过程，这一特征应该与该地区温度降低导致土层的稳定性发生改变相关，其滑动过程是否随季候变化呈现周期性还需要在更长时间周期内监测进行验证。

LDhp15、LDhp17、LDhp19、LDhp21、LDhp23、LDhp38、LDhp42、LDhp47其最主要的特征是在2019年6月～12月相对于其他时间段始终保持缓慢滑动的状态，而在2019年1月～5月份或是2020年1月～5月份则出现了明显加速滑动的过程，这种在非生长季滑动过程加快应于坡体上植被的季候变化相关。

文中研究基于Sentinel-1卫星影像进行Stacking和IPTA处理，对隆德县境内的滑坡隐患进行探测，核实探测到的滑坡隐患共计47处，其中，Stacking技术有效探测滑坡隐患44处，IPTA有效探测滑坡隐患37处，2种InSAR技术共同识别34处。基于2种InSAR技术结果：

（1）Stacking和IPTA技术所探测到的滑坡隐患最大形变速率线性相关性强，形变测量结果可相互印证，可以认为基于InSAR技术进行滑坡隐患的探测具有较高的可靠性。

（2）由于IPTA需要选择高质量的点目标从而导致部分滑坡隐患无点目标覆盖，而使用Stacking方法能够对这些滑坡隐患进行有效识别。点目标的缺失与野外环境复杂有关，如植被等地表覆盖物造成相干性降低，地形起伏在局部地区出现大量地形相关的大气相位。后续应联合多种InSAR技术，如采用DS InSAR方法增加同质点的选取［28］，采用BM3D进行更好的滤波降噪，提高IPTA分析的点目标覆盖范围和测量精度，从而提高滑坡隐患探测的可靠性。

（3）在形变速率上，隆德县以黄土丘陵区为主要地貌特征，较少存在剧烈的坡体滑动过程，因此Stacking和IPTA技术能够对该地区不同形变速率的滑坡隐患进行有效探测。同样由于IPTA的选点机制导致滑动较的坡体候选点减少，其所能识别的滑坡隐患形变速率中位数略低于Stacking技术。

（4）在地形特征上，Stacking和IPTA所探测到的滑坡隐患主要集中在缓坡和斜坡上，同时除北坡不易出现滑坡外，在各个坡向的滑坡隐患均可探测。

（5）基于IPTA获取的时序形变量表明，隆德县境内的滑坡隐患存在4种典型滑动类型：（a）以稳定的速率持续滑动；
（b）在2019年6月～9月份出现明显加速滑动过程并持续到监测结束；
（c）在2019年6月～9月份出现明显加速滑动过程，10月～12月份滑动过程逐渐减缓；
（d）在2019年或2020年1月～5月滑动过程明显，其它时间坡体较稳定。

总之，综合Stacking和IPTA这2种InSAR技术能够全面的识别隆德县境内的滑坡隐患，在后续滑坡隐患探测中应联合多种InSAR方法，以获取更可靠的结果。下一步将首先基于GNSS等技术对Stacking InSAR和IPTA技术得到的平均形变速率和时序形变结果进行精度评价，并联合其他InSAR技术提高IPTA的探测能力。

猜你喜欢隆德县基线滑动记一次写生速读·上旬(2021年7期)2021-08-05航天技术与甚长基线阵的结合探索科学(2020年5期)2020-11-26一种SINS/超短基线组合定位系统安装误差标定算法中国惯性技术学报(2019年3期)2019-10-15一种新型滑动叉拉花键夹具制造技术与机床(2018年11期)2018-11-23Big Little lies: No One Is Perfect意林(绘英语)(2018年1期)2018-04-28隆德县三通古碑考释西夏研究(2017年4期)2017-08-17隆德县县属国有企业发展现状浅析现代经济信息(2017年3期)2017-04-05一种改进的干涉仪测向基线设计方法舰船电子对抗(2016年5期)2016-12-13滑动供电系统在城市轨道交通中的应用城市轨道交通研究(2015年11期)2015-02-27隆德县农村电影放映创新发展之路西部广播电视(2015年24期)2015-02-26 相关热词搜索：隆德县，宁夏，应用于，