基于时间序列Insar的黄土高原城镇滑坡识别研究

殷宗敏，韩守富，杨 环，马金辉

(兰州大学资源环境学院，甘肃 兰州 730000)

绥德县位于陕北黄土高原的东北部，区内总的地势为：北高南低，中东部无定河与黄河的分水岭高于两侧，无定河两侧梁峁地势高于河谷谷地，黄河谷地最低。县域位于鄂尔多斯地台东部，构造活动比较弱，区内仅发育少量平缓开阔的褶皱。随着人类改造自然能力的增强，人类活动对自然环境的改变越来越大[1]。绥德县与地质灾害关系密切的人类工程活动主要有交通建设、水利设施建设、群众修建房屋时削坡过陡等，这些建设生产过程导致滑坡情况时有发生，造成人员伤亡、财产损失。

图1 绥德县城周边发育滑坡

滑坡是区内最为严重的地质灾害之一。境内滑坡主要发育在河谷两侧边坡上和黄土梁峁区的陡坡地段，稳定性差，大多为黄土滑坡，如图1所示。其发生一般预兆不明显，运动速度快，发生过程短促，一旦发生，往往造成人员财产损失或者影响交通，因此对该区域进行地表形变监测尤为重要。利用时序InSAR技术进行地表形变监测，提取研究区长时间序列地表变形数据[2]，基于地表变形速率圈定形变速率较大的不稳定区域，该区域可能为滑坡区。通过实地调查等方式，对这些区域进行甄别，以便进行相应的防灾减灾工作。由于城区人口聚集且存在发生滑坡的条件，因此对于该区域的监测尤为重要，故选择的研究范围如图2所示。

图2 研究区范围

1 数据来源及研究范围

1.1 数据基本信息

本文选择Sentinel-1A IW SLC数据，轨道为11，下载网址为：https://scihub.copernicus.eu/。影像截止日期为2017年9月20日，覆盖研究区的影像共38景，数据基本特征如表1所示。

Sentinel-1A数据大的空间基线、多普勒中心和数据的时间分布情况如图3所示，其中空间基线以0为中心，正态分布于两侧，值在150 m以内，这为图像高精度配准提供了保障，多普勒中心分布于0～-0.1之间，影像的获取时间以2016年12月30日为中心，且中心时间点以后的影像更为密集。

1.2 数据时间分布

Sentinel-1A的11轨道为升轨道数据，以2016年12月30日影像为主影像，影像的时间序列分布情况相对均匀，其时间序列分布以及相应的空间基线、时间基线、多普勒中心的分布情况如表2所示。

表1 所用数据基本信息

表2 所用Sentinel-1A IW SLC数据时间分布表

2 数据处理方法

2.1 时序Insar参数优化和处理

研究区地表覆盖变化大，失相干现象严重，难以捕获到永久散射体[3]。滑坡变形既包括缓慢蠕动，又包括快速变形，不仅仅存在线性变形，更多的是非线性变形。针对InSAR处理过程中的主要环节和目前所面临的问题，为提高滑坡区InSAR测量精度，选择合适的图像连接拓扑、滤波、解缠方案及参数；进行大气延迟效应估算，获得适合研究区地形和大气特征的大气延迟效应参数[4]。

2.2 图像连接拓扑、滤波

通过不同方案和参数的设定，提高滑坡区InSAR测量精度[5]，扩大InSAR测量范围，针对特定的SAR数据集，建立适当的图像拓扑，如：基于单主图像的星型结构、基于空间基线最佳的多主图像拓扑结构、基于时间基线的多主图像拓扑结构，分析其弃除大气效应的不同效果从而选择最佳图像拓扑，本文最终选择星型结构；选择适合研究区地表覆盖、变形规律的滤波方式，如：Boxcar Goodstein和Adaptive Filter等，并进行参数调整，分析其对干涉效果和最终变形结果的影响，本文最终选择Adaptive Filter滤波方法。

2.3 剔除或减轻地形影响

选择了不同来源和分辨率的DEM数据[6]，测试其分离地形效应的效果，获得能有效分离地形效应的DEM精度大小和分辨率等指标，这里选择90米分辨率dem；采用多角度观测来减缓因地形产生的SAR图像几何畸变(透视收缩、阴影、叠掩)，增加SAR测量点密度。

图3 轨道数据集基本参数图

2.4 大气延迟效应估算和剔除

Processing Parameters分区是SARProZ软件反演大气相位的主要功能模块，其通过干涉计算获得各像素点的干涉相位，并依据研究区域的地表形变情况以及干涉相位中各分量的物理及统计特性建立其与高程误差、线性形变、热形变等因子参数的缠绕相位函数模型，根据大气相位分析像素点网邻近永久散射体目标(像素点)间的相位差分关系，采用相关算法考察相邻像素点差分干涉相位关系对此函数模型进行求解，解算出各像素点的高程误差、线性形变等贡献因子相位，并对残余相位中的非线性形变与失相干噪声相位进行分解，最终计算获得精确的大气延迟分量缠绕相位值。大气延迟相位是InSAR测量的主要干扰误差源[7]，需要对该误差进行剔除，其处理步骤为：通过调整参数，有效的选取PS候选点；对选取的PS候选点构建合理的分析像素点网；结合研究区域的现状及数据集的情况，给出合理的参数，分析大气延迟相位[8]；选择合理的参考点并检验各参数相位分量的合理性。

图4 研究区形变速度图

3 结果与分析

3.1 实验结果

将覆盖研究区的不同时段的多期SAR影像按照成像时间先后顺序，以选取的公共主影像为基准，将其余所有影像都配准并取样到主影像空间[9]，总共可形成K个干涉对。通过相位差分处理，可得到K幅干涉相位图。然后借助提取影像中永久散射体(PS)目标，通过分离地形数据误差、大气延迟误差等一系列的误差项，从而提高PS的地表形变估计精度，最终得到可靠的地表形变结果。将获取的可靠的时间序列形变结果导出为csv和kml文件，将kml结果在谷歌地球中显示，其结果如图4所示，其中形变速度值范围为-15mm/year至15mm/year,颜色越深表示形变速度越大，接近-15mm/year为红色，接近15mm/year为蓝色。

结合地形指标、ps点的形变速度值和聚集情况，勾画出可能的滑坡区。发现其多分布于山坡部分且一侧为道路或者河流谷底(见图4)，符合自然条件下滑坡的发生区域范畴，说明基于时间序列的PS点能够识别滑坡。

3.2 结果分析

视线方向的形变结果在空间呈现聚集程度不一的点要素，其中空间分布呈现聚集状态的点要素，其聚集区域往往代表地表变化大的区域，在人工活动区域多表现为建筑开挖，而非人工影响区域则多表现为滑坡，在谷歌地图中对点聚集区域沿着山坡绘制可能滑坡区，并进行滑坡的野外验证。对滑坡前缘进行调察时，对比滑动面、滑坡前缘边界，需要注意的是，滑坡前缘一般存在剪出口(剪切破坏明显凸起甚至反翘的部位)，如果滑坡是由于坡脚开挖引起，前缘边界应在坡脚填土前缘。

野外验证发现：先前圈定的滑坡区域，在现场基本都找到了滑动痕迹(见图5)。其中，图片1现场为施工区域，坡面呈现滑动迹象，且前方道路塌陷，表示该区域地表发生了较大形变，与Insar计算的结果吻合；图片2为一个滑坡体前的窑洞内部情况，由于山体滑动，向前产生推力，导致墙体夯土向窑洞内侧脱落；图片3为无定河一侧山谷中山体滑动后形成的阶梯状痕迹，其下窑洞存在较大的安全风险；图片4为公路一侧的山体，其存在明显的滑坡后壁，滑动痕迹明显。通过实地的考察对比分析，证明利用时序Insar技术进行地表形变监测，对形变速度较大且形变点聚集的区域进行圈定，能够找出滑坡区域，指导防灾减灾工作。

图5 野外滑坡验证对比图

在众多滑坡体中，选定3个明显的单体滑坡(单体滑坡轮廓见图6)，获取该区域的形变曲线，进行滑坡单体变形研究。定位滑坡体上的形变点，并在时序Insar处理过程中获得的csv文件里面找出该点对应的各时间段形变值，将其绘制成折线图(见图7)。

图6 滑坡单体轮廓图

图7 坡体变形曲线图

从A坡变形曲线可以看到，坡体A在2016年8月份左右发生了较大规模的滑动，之后呈现较缓的滑动，由于该折线整体呈现下降趋势，说明滑动在持续。

B坡变形曲线中，发生较大规模滑动的时间点为2016年9月份前后，与A滑坡体类似，形变曲线呈现下降趋势，形变曲线两区间段的形变量差值大于A滑坡体，说明其不稳定性比A滑坡体强。

C坡变形曲线与A、B有所不同，其在2016年10月形变量开始趋于稳定，在2017年7月之前，形变速度基本为正值，在7月份以后呈现较大变化，说明较大滑动发生在2017年7月，C滑坡体形变量相对小于A、B，说明C滑坡体较A、B稳定。

综上可以看出，通过时间序列Insar获取的该区域形变点属性值能够反映出滑坡体的变形过程，说明Insar技术能够很好的识别出滑坡体。

4 讨论与结论

4.1 讨论

基于时间序列Insar技术，能够获得缓慢的地表形变过程，其在平原地区具有较好的效果[10]，而在山区由于植被的覆盖以及山体对雷达信号的遮挡，降低了Ps-Insar技术的可靠性，但是可以通过参数的调整，逐步剔除部分不准确相位，提高结果的可靠性。此外，参数的选择具有极强的地域性以及技术性，需要通过不断的试验来获取最佳的参数，而试验又存在主观性，因此只能不断接近最佳的参数集。

由于卫星获取的形变方向为视线方向[11]，卫星的轨道固定，同时山体不同部分的坡度坡向不一致，导致获取的形变结果与坡体形变的真实情况存在一定的误差，而目前还没有很好的方法解决这些偏差，有待进一步考究。

4.2 结论

基于PS-InSAR技术，利用Sentinel影像数据，通过影像裁剪、影像配准、地形相位估计、差分干涉、高相干点选取、相位解缠、变形提取和误差估计等一系列数据处理手段计算的黄土高原重要城镇时间序列地表变形速率具有很高的可靠性，能够很好的识别黄土高原区滑坡区，滑坡体上形变点的形变曲线能较好的反映出坡体的变形情况。

[1]黄润秋，许强. 中国典型灾难性滑坡[M]. 北京：科学出版社.2008.

[2]Michele Crosetto , Oriol Monserrat , María Cuevas-González , Núria Devanthéry, Bruno Crippa ．Persistent Scatterer Interferometry[J]．ISPRS Journal of Photogrammetry and Remote Sensing，2015.

[3]Dong S, Samsonov S, Yin H, et al. Time-series analysis of subsidence associated with rapid urbanization in Shanghai, China measured with SBAS InSAR method[J]. Environmental Earth Sciences, 2014.72(3):677-691.

[4]Schl gel R, Doubre C, Malet J P, et al. Landslide deformation monitoring with ALOS/PALSAR imagery: A D- InSAR geomorphological interpretation method.[J]. Geomorphology.2014.231: 314-330.

[5]陈强．基于永久散射体雷达差分干涉探测区域地表形变的研究[D]．成都：西南交通大学博士论文.2006.

[6]傅文学，田庆久，郭小方，等．技术及其在地表形变监测中的应用现状与发展[J]．地球科学进展.2006.

[7]陈强，李永树，刘国祥．干涉雷达永久散射体识别方法的对比分析[J].遥感信息.2006.

[8]刘国祥，陈强，丁晓利．基于雷达干涉永久散射体网络探测地表形变的算法与实验结果[J]．测绘学报.2007.

[9]汪鲁才．星载合成孔径雷达干涉成像的信息处理方法研究[D]．湖南大学博士论文.2005.11.

[10]李德仁，縻明生，王艳．永久散射体雷达干涉测量技术[J]．武汉大学学报(信息科学版).2004.

[11]彭曙蓉．高分辨率合成孔径雷达干涉测量技术及其应用研究[D]．湖南大学.2009.