利用Sentinel-1A的时序干涉测量探测蒲县滑坡形变

高二涛， 李 豪， 雍 琦， 付波霖， 范冬林

(1.桂林理工大学测绘地理信息学院，桂林 541006； 2.广西空间信息与测绘重点实验室，桂林 541006；3.山西省工业设备安装集团有限公司，太原 030000)

中国地质环境复杂、气候条件时空差异大，地震、滑坡、泥石流等地质灾害类型多且分布广，难以有效防范[1-2]。在上述自然灾害中，滑坡是一种危害程度仅次于地震，排名第二的环境地质灾害。滑坡广泛发育于河流、湖泊、海洋和沟渠的岸坡，地形高差大的山区、峡谷区，道路的斜坡地段及强降雨地区等，具有隐蔽性、突发性、不确定性等特点，且破坏力巨大[3]。蒲县属于晋西黄土高原地质灾害区，地貌复杂多样，此地区应重点加强对居民居住区和采煤区的黄土滑坡、崩塌等地质灾害的防范。滑坡在突发前一般伴随长期缓慢的地表变形，在此背景下，基于变形监测为主的滑坡灾害早期识别与持续监测可以为预报预警与防灾减灾工作提供科学依据，对于主动防范地质灾害意义重大[4]。

水准测量、全球导航卫星系统(global navigation satellite system，GNSS)观测、摄影测量法和位移计等传统滑坡监测方法以其点观测、高精度等特点，目前已被广泛应用于滑坡监测中[5-7]。近年来，空间对地观测技术迅速发展，以永久散射体雷达差分干涉测量技术(persistent scatterer interferometric synthetic aperture radar，PS-InSAR)与短基线集雷达干涉技术(small baseline subset，SBAS)为代表的时间序列合成孔径雷达干涉(time series interferometric synthetic aperture radar，TS-InSAR)技术，可以远程监测地球表面以获取全球的高精度地形和地表微小形变信息，具有全天时、全天候、时空分辨率高、覆盖范围广等特点，弥补了传统滑坡形变监测方法计算量大且监测成本高的不足，在滑坡、地震、城市地面沉降等地质灾害的普查与监测预警中展现出巨大的发展潜力[8-12]。现采用PS-InSAR与SBAS技术分别处理覆盖山西蒲县地区24景哨兵-1A(Sentinel-1A)数据，获取研究区域2019年年平均形变速率和累计形变位移量，通过结合已有的资料和谷歌地球(Google Earth)等数据，进一步分析是否为滑坡隐患点及其可能性成因。

1 研究区概况与数据简介

1.1 研究区概况

蒲县(36°11′32″N～36°38′13″N，110°51′09″E～111°23′36″E)位于山西省西南部，临汾盆地中心，吕梁山西南侧，东部属于吕梁山和太行山断裂带，西部属于鄂尔多斯断裂带。蒲县东接洪洞，西靠大宁，南依吉县，北与隰县、汾西毗邻。境内由山地、丘陵、河谷等地形组成，地势剖面为东高西低，北、东、南三面环山，植被茂盛。西部与中部为典型的黄土高原地貌，第四纪黄土发育良好，地表大多被不同厚度的黄土覆盖，地形破碎、垂直节理发育、沟谷深切超过200 m，易引发崩塌、滑坡、泥石流等多种次生地质灾害。蒲县是山西煤炭开采的重要基地之一，含煤面积占全县总面积的九成以上[13]。蒲县位于北半球中纬度大陆内部，是典型的温带大陆性气候。春季多风沙，降雨量少；夏季多暴雨，雨热同期；秋季短暂、温和凉爽；冬季漫长、寒冷干燥。蒲县地理位置如图1所示。

图1 蒲县地理位置

1.2 实验数据

哨兵1号(Sentinel-1)卫星是欧洲航天局(ESA)发射的载有C波段合成孔径雷达的地球监测卫星[14]，可全天时、全天候提供连续图像。哨兵1号是哥白尼计划的一个重要组成部分，具有1A和1B双系统。2014年4月发射的Sentinel-1A卫星是太阳同步极地轨道卫星，卫星的重访周期为12 d，轨道高度是693 km。它的主要用途是监测地球环境变化，存在四种成像模式、两种极化方式，具有生产能力强、重访周期短的优势。现采用Sentinel-1A卫星获取的24景干涉宽幅模式(interferometric wide swath，IW)的单视复数图像数据(single look complex，SLC)。获取影像的时间范围是2019-01—2019-12，数据的详细参数如表1所示。

表1 Sentinel-1A TOPS SAR影像数据参数

本文研究采用日本宇宙航空研究开发机构获取的30 m分辨率的AW3D30(ALOS World 3D-30 m)作为外部DEM(digital elevation model)数据。DEM的范围在35°N～38°N、108°E～112°E，运用其生成地质背景构造图和去地形相位处理，以消除地形起伏的影响[15]。

2 时序InSAR数据处理

2.1 基于PS-InSAR技术的蒲县滑坡形变监测

首先，利用PS-InSAR技术监测蒲县地区滑坡形变过程，主要处理步骤如图2所示。实验程序自动选取2019-06-30的影像数据为公共主影像，其余23景从影像逐个与之配准并建立连接关系，共得到23个干涉像对。SAR数据时空基线分布如图3所示，具体时空基线如表2所示。可知，时间基线范围为-180～168 d，空间基线范围为-135～110 m。

表2 Sentinel-1A影像数据获取时间及基线参数信息

图2 PS-InSAR处理流程

图3 PS-InSAR时空基线分布

使用AW3D30 DEM对干涉相位进行去平处理，实验采用的DEM精度和分辨率越高，去除地形相位的效果越理想。打开Goldstein滤波器，将距离向多视数设置为5，方位向多视数设置为1，输出差分干涉图。选择永久散射体(PS点)时，应选取有稳定散射特性的点，如桥墩、大坝、房顶等人造地物或者裸岩等自然地物。通过密集分布的PS点消除由于信号在大气中传播延迟造成的波动，从而有效改善实验精度。为提高稳定散射体点的准确性，采用振幅离差指数法确定PS点，将相干系数阈值设置为0.75，对研究区域内的稳定散射目标进行有效识别[16]，如图4所示，在图4中用绿点标注，PS点目标主要分布在裸露的岩石、公路上，而植被覆盖率高的地区几乎没有PS点。

图4 PS-InSAR技术永久散射体(PS点)分布

利用重去平的干涉图建立线性模型，估算变形率、残余高度等。大气效应的估算和清除通过两种滤波程序进行，分别是大气高通滤波和大气低通滤波。对大气校正后的结果进行地理编码，即可得到蒲县地区2019-01—2019-12年平均沉降速率(图5)。对PS-InSAR技术所得结果图进行观察分析可知，研究区域内形变速率范围为-53～48 mm/a。研究区域绝大部分地区处于稳定状态，形变速率范围为-15～10 mm/a。

沿雷达视线方向LOS (line of sight)地表抬升的地区用红色表示，发生地面沉降的区域显示为蓝紫色

2.2 基于SBAS技术的蒲县滑坡形变监测

为验证PS-InSAR方法监测的准确度与可靠性，采用SBAS方法在同一数据源的基础上对山西蒲县地区进行滑坡监测，SBAS处理流程如图6所示。

图6 SBAS处理流程

采用蒲县地区2019-01—2019-12的24景Sentinel-1A数据，设置临界基线最小、最大百分比分别为0和2，时间基线阈值设定为60 d，系统自主选取2019-06-18的影像为超级主影像。根据设定得到阈值，共生成170个干涉像对，图7为SBAS技术时空基线分布图。

图7 SBAS时空基线分布

为使实验结果在距离向和方位向的分辨率尽可能保持相同的效果，实验将方位向和距离向多视数分别设置为1和4，与超级主影像保持一致。选用常被应用于植被密集或环境潮湿地区的最小费用流法(minimum cost flow， MCF)进行相位解缠，将解缠相关系数阈值设置为0.23，利用Goldstein滤波方法以提高研究区中低相干位置的滤波强度。

打开存放差分干涉成果的文件夹，仔细查看每一个像对的相干性图片以及解缠后的图片，用连接图编辑工具剔除相干性低和含有明显大气相位的像对。轨道精炼与重去平是为了估算差分后依旧残存的恒定相位和解缠后还留有的相位坡道，利用GCP点对所有的数据对执行重去平以消除地形相位。在实验过程中选取了远离形变区域、孤立相位、残余地形条纹等位置，位于相对稳定地点的70余个GCP点。在第一次反演中，对相干点建模并构建方程组矩阵，由于线性模型的稳定性最佳，因此，实验选用线性模型求解所有影像对的变形率和新的DEM等。基于第一次反演所得形变速率，第二次反演是通过设置大气高通和大气低通来估计和清除实验中的大气影响，获得更为准确纯净的时间序列中的位移结果并进行地理编码，将其转换到地理坐标系，获得形变信息和年平均位移速率，如图8所示。研究区域形变速率范围为-70～52 mm/a，有4处地表沉降较为明显的疑似滑坡隐患地点，大部分地表形变在-10～10 mm/a，趋于稳定状态。通过SBAS方法采集的蒲县地区2019年清晰明了的时序地表变化过程如图9所示。

红色表示雷达视线方向存在抬升，蓝色表示下沉

图9 基于SBAS技术的蒲县地区2019年地表累计形变

3 时序形变结果分析

3.1 结果对比

基于PS-InSAR技术和SBAS方法，利用哨兵升轨SAR数据分别获取蒲县地区年平均形变速率，通过目视解译选取出两种技术所得结果中共同存在的4处较为明显和典型的形变位置进行对比分析，所选位置分布如图10所示，包括太林乡半沟村罗克线沿线、水泉窊015乡道两侧、郭家山周边、井子洼村赵克路附近，具体地理位置如表3所示。图11为4处形变速率放大图，从放大图中可以清楚地看到4个疑似隐患点的沉降信息。除图11所示4处较明显的疑似滑坡隐患点，研究区域仍存在部分升降不一的形变点，为使结果更加精准可靠，只对上述4处共同存在的疑似滑坡点做进一步对比分析。

图11 4处疑似滑坡点SBAS技术形变速率放大图

表3 4处疑似滑坡点详细地理位置

图10 4处疑似滑坡点位置

为了更加精准地分析PS-InSAR技术与SBAS技术对于同一数据源的结果精度对比情况，基于4处疑似滑坡点累计形变量的详细信息，计算每一期数据的差值，并用两种数据所得结果的平均值代替真值计算其均方根误差(RMSE)，RMSE代表每个数据与真值间的偏离程度，均方根误差越大，表示两种技术获取的结果差异越大，可靠性越低；反之，则表示测量的精度越高，进一步说明时序InSAR技术识别并监测滑坡的有效性，如表4所示。

由表4可得，一号疑似滑坡点PS-InSAR技术与SBAS技术所得最终累计沉降量分别为-51.23 mm 和-52.32 mm，差值为1.09 mm，RMSE范围在0～2.25 mm，RMSE均值为0.75 mm；二号疑似滑坡点PS-InSAR技术与SBAS技术所得最终累计沉降量分别为-48.03 mm和-49.51 mm，差值为1.48 mm，RMSE范围在0～2.06 mm，RMSE均值为0.77 mm；三号疑似滑坡点PS-InSAR技术与SBAS技术所得最终累计沉降量分别为-40.87 mm和-42.42 mm，差值为1.55 mm，RMSE范围在0～3.25 mm，RMSE均值为1.22 mm；四号疑似滑坡点PS-InSAR技术与SBAS技术所得最终累计沉降量分别为-46.91 mm和-39.80 mm，差值为7.11 mm，RMSE范围在0～3.55 mm，RMSE均值为1.05 mm。由上述数据可得PS-InSAR技术与SBAS技术监测结果的均方根误差在0～4 mm，均值在0.5～1.5 mm，充分说明了PS-InSAR技术与SBAS技术的一致性较高。

表4 4处疑似滑坡点累计形变量的详细信息

为了使两种时序InSAR技术所得监测结果对比的更为直观简明，基于上述4处疑似滑坡点的累计形变量数据，对选取出的4处点绘制累计形变折线图进行对比分析。图12为识别出的4处疑似滑坡点2019-01-01—2019-12-15的LOS向累积形变位移量对比。总体看来，两种技术所得结果的沉降趋势较为一致，说明2019年蒲县地区部分地区发生地表沉降，有疑似滑坡的可能，更进一步说明了时序InSAR识别并监测滑坡的可行性与准确性。1号疑似滑坡点在1—5月地表相对稳定，6月出现小幅沉降至20 mm，7月地面出现略微抬升，8—10月缓慢沉降，并在11、12月急剧下沉至-50 mm左右。2号疑似滑坡点1—6月缓慢沉降至-21 mm左右，并在7、8月出现小幅抬升，随后逐渐沉降至-49 mm左右。3号疑似滑坡点沉降成线性变化，4号疑似滑坡点1—4月地表处于稳定状态，随后呈线性下沉模式。

图12 4处疑似滑坡点累计形变位移量PS-InSAR与SBAS监测结果对比

3.2 可能性诱因分析

蒲县是典型的温带大陆性气候，雨热同期。由上述分析可知，疑似滑坡点7、8月地面出现略微抬升，可能是由于雨季降雨量骤增的缘故，地表水渗入地下，导致地下含水层增厚。11、12月地表下沉明显，可能是由于冬季降水量过少，当地居民大量使用地下水的缘故导致地下水位下降。3号疑似滑坡隐患点位于郭家山周边，郭家山地区存在采煤区，该点的沉降疑似受到煤炭开采等人类活动的影响。

4 结论

利用时序InSAR技术探测研究蒲县地区的滑坡形变，获取研究区域内2019年年平均形变速率和累计形变位移量，确定太林乡半沟村罗克线沿线、水泉窊015乡道两侧、郭家山周边、井子洼村赵克路附近等4处位置为疑似滑坡点，并结合煤炭开采、人类活动、降雨量信息等，对该地区分布的疑似滑坡隐患进行早期识别与监测预警，为该地区滑坡诱因提供科学的依据，主动规避自然地质灾害风险，对未来滑坡灾害的监测与预警提供思路与方法。