基于InSAR技术对汕头潮阳区地质灾害体识别分析

方 刚，叶 波，纪振付，张明亮，李建颖，李晓淳

(广东省有色金属地质局九三一队，广东 汕头 515041)

合成孔径雷达干涉测量(InSAR)技术是利用多次重复观测的合成孔径雷达获取地表微小变形，利用具有一定视角差的两部天线(或一部天线重复经过)来获取同一地区两幅具有相干性的SAR单视复数影像，并根据其干涉相位数据提取地表的高程信息。InSAR技术可以快速地提供被监测区域内大范围、准确的地形与地表形变信息[1]，便于揭示地学现象的时空变化规律，在地质灾害监测与识别中发挥了重要作用。

1 数据源获取

本次InSAR地质灾害识别范围为广东省汕头市潮阳区全区[2]，总面积666.63km2，重点调查区面积103.47km2。选用SAR影像数据源为Sentinel-1A雷达卫星数据[3]，该卫星由欧空局于2014年4月发射的一颗雷达卫星，与2016年4月发射的Sentinel-1B卫星组成双星伴飞模式，重访周期从12天缩短到6天。IW模式采用了最新的TOPs成像技术，能够解决宽幅成像时出现的scalloping效应并增强成像辐射性。Sentinel-1卫星配备了新一代C波段合成孔径雷达系统(波长5.6cm)，相较于L波段，波长更短，形变探测能力更强。

在欧空局官方网站上下载覆盖研究区C波段Sentinel-1A卫星IW模式下level-1A级影像，分辨率均为5m×20m，覆盖研究区Sentinel-1升轨影像，时间跨度分别为2017年1月1日至2022年8月15日共168景。采用DEM数据主要是30m分辨率的Aster GDEM，其标称高程精度10m～25m，处理Sentinel-1A数据前，下载与之对应的精密轨道数据并对轨道信息进行修正，有效去除因轨道误差造成的系统性误差。

2 数据可视性分析

SAR地形可视性取决于SAR采集数据的几何参数、视线(LOS：light of sight)方向的角度和地形坡度和坡向，根据这些参数对研究区叠掩和阴影区分和提取。建立SAR地形可视性分析流程(图1)，对Sentinel-1影像数据进行地形可视性分析(图2)，结果表明所选用的影像阴影和叠掩分别较少，仅在山区少量分布，阴影区域占0.2%，叠掩区域占0.01%。

图1 SAR 地形可视性分析流程图

3 InSAR数据处理流程

3.1 D-InSAR处理流程

该技术是InSAR形变监测的基本处理方式，利用两景高分辨率SAR影像可有效进行坡体短期厘米级变形特征识别，包括滑移面积、形变量信息。处理流程见图3。判别一定范围内潜在的滑移坡体，要求地质灾害体表面保持良好的相干性[4-6]。两景高分辨率SAR影像时间间隔一般1～2个月，适合短期内厘米级滑坡变形监测。其缺点是对于长期毫米级缓慢形变的滑坡体，差分干涉SAR厘米级形变监测精度难以识别，单一监测角度也不足以克服SAR影像叠掩和阴影影响，难以准确判定长期毫米级缓慢形变滑坡的整体范围、边界、强度。D-InSAR干涉过程中生成干涉图(图4、图5)是识别形变区域的重要工具，对所有干涉图浏览，筛选出多处具有明显干涉条纹的点位。

图4 部分全域干涉图

图5 部分局部干涉图

3.2 SBAS-InSAR处理流程

SBAS-InSAR技术是基于D-InSAR技术基础上的发展，其核心是时空基线阈值设置，对于试验区而言，由于高相干点较少，因此除了需要确定时空基线阈值外，适宜的空间范围也尤为重要。通过对调查区SBAS-InSAR处理[7-10]，反演地表形变速率。具体流程(图6)。

图6 SBAS-InSAR数据处理工作流程图

3.2.1 数据预处理

(1)SAR主影像选择和像对组合。a)计算所有影像像对间时间和空间基线，生成时间和空间基线分布图；b)采用时间和空间基线均满足给定阈值的像对组合生成差分干涉图集。(2)所有SAR影像对一景影像进行配准、裁剪，并组合生成时间序列干涉图集。a)选择非夏季、时空基线尽量居中的影像作为配准参考影像，所有影像对其进行配准；b)将所有数据裁剪成一致的区域；c)对所有配准好的干涉像对，按时间和空间基线限制条件，选择像对组合。逐像元计算干涉相位，生成时间序涉图集；(3)将DEM与配准参考影像进行配准，将DEM范围裁剪成与配准参考影像一致区域；(4)将所有主辅影像前置滤波，计算干涉相位，生成干涉图；(5)相干系数计算；(6)相干点目标选取，对时间序列干涉图集的像元进行相干点目标筛选。a)相干点目标选取；b)相干点目标干涉相位序列生成。将满足上述条件要求的辅影像与主影像进行相位干涉处理，提取相干点目标干涉相位序列图。

3.2.2 差分干涉计算

(1)平地和地形相位去除，对相干目标点组成的干涉图，进行平地和地形相位去除；(2)差分干涉图滤波；(3)相位解缠。

3.2.3 时间/空间域变形估算

对干涉图的差分干涉相位应进行时间和空间域线性变形相位估计，如有要求，还应进行非线性变形相位估计，去除大气、噪声等残余相位，得到点目标时间序列变形相位。

SBAS-InSAR计算步骤：a)相邻点间参数估计；b)线性变形相位和残余高程计算；c)残余相位低通滤波。从差分干涉相位中减去步骤：a)中两项相位分量后得到残余相位，并进行空间域低通滤波得到滤波后的残余相位；奇异值分解处理是根据短基线像对组合关系，对步骤b)得到的滤波后残余相位进行奇异值分解(SVD)处理，求解每个影像对应时刻的大气相位和非线性变形相位；大气相位和非线性变形相位计算。对奇异值分解得到的大气相位和非线性变形相位进行空间域高通滤波，得到大气相位，并对滤波后的相位序列进行时域低通滤波，得到非线性变形相位；时间序列变形相位计算是将b)步骤中线性变形相位和大气相位和非线性变形相位计算中非线性变形相位相加，结合时间基线参数，得到每个相干点目标的时间序列变形相位。

3.3 关键参数设置

本次试验采用SBAS-InSAR(小基线集)技术数据处理，其核心是时空基线阈值的设置，对于试验区而言，由于高相干点较少，因此除了需要确定时空基线阈值外，适宜的空间范围也尤为重要。依据多次试验结果，最终选择时空基线分别为24d和2m，超级主影像为2017/1/25影像，共634组干涉对。连接图如图7所示，图中方形点代表超级主影像，圆形点为其他影像。小基线干涉对组合图能清晰地反映超级主影像与其它影像的相互位置以及图像参与运算的情况。由于本研究采用的SAR影像时序较多，因此干涉对组合情况较好，每个时相都与其它时相建立了连接，且连接较为均匀。

图7 影像时间连接图

3.4 干涉工作流

InSAR受时间过长而失相干、几何位置由于DEM影像或卫星影像错位失相干、多普勒效应造成的中心频率失相干和由于温度产生热噪声等失相干影响。干涉配准后获得的干涉图表面含有随机噪声，以灰白色圆点或者直接失相干而失去影像，严重时会干扰条纹条带分布或失去干涉条纹。需要对干涉图滤波处理，主要目的是减少相位解缠中粗差。

本次研究数据量较多，选用Goldstein滤波方法处理。干涉工作流是对所有配对的干涉像对干涉处理，包括相干性生成、去平、滤波和相位解缠等过程，最终所有数据对都会被配准到超级主影像上，以便为轨道精炼与重去平、第一次反演和第二次反演做准备。干涉工作流处理过程中需设置多视的视数，可在SAR影像文件中查看相应参数。估算多视的视数，若输入的数据分辨率不同，程序会对每个像对多视处理，得到与超级主影像一致的分辨率。多视可以增加干涉图的信噪比，提供可靠的相干性，也可提高运算速度。此处设置为4：1。制图分辨率为20m，解缠方法为Minimum Cost Flow，解缠相关系数阈值为0.2，滤波方法为Goldstein。高程精度阈值为20，速率精度阈值为32。

选取2021年不同时间段对相同时间基线的干涉对比分析(图8)。结果表明，研究区内冬季相干性较高，干涉条纹鲜亮，可以识别到明显形变区域；而植被较为茂盛区域，不管在任何时间，相干性都较低(0.2～0.3)，这部分区域噪声较多，也是Sentinel-1在植被茂盛区地表形变识别的不足。此外，在解译区东侧与南侧存在水田区域，没有任何相干性。

3.5 轨道精炼和重去平

在研究区选择相对稳定区域，远离形变区的控制点。在轨道精炼和重去平过程中输入GCP点，估算和消除残余的固定相位以及解缠后还存在的相位残差。GCP点既可选择在干涉结果较好像对的解缠结果图上或者干涉图上，刺点生成GCP文件，也可选择外部GCP点。在做地表形变监测时，不能确定刺点位置为稳定区域，尤其是研究区面积大、地形复杂，选择GCP点位置存在极大的不确定性。本次在做地表形变测量时，对比多次遥感影像，将稳定点位置(如稳定建筑物、路桥等构筑物，不存在形变条纹)作为GCP点，以此降低选点不当带来的误差。共选用182处GCP点，保证GCP均匀分布(图9)，控制点标准差小于1。

3.6 反演与地理编码

SBAS反演核心就是第一次反演，主要估算研究区形变速率和残余地形。第一次反演相关系数阈值0.2，解缠相干系数阈值0.2，解缠方法为最小费用流法。第二次反演是获得时间序列上形变量。利用之前轨道精炼和重去平步骤中GCPs点文件移除恒定的相位和斜坡相位，根据第一次反演结果进行大气滤波处理，最后得到时间序列位移结果。第二次反演时相关系数阈值0.3，大于该阈值参与计算。轨道精炼方法选用多项式精炼方法，且精炼残余相位多项式次数设置为3。

4 解算分析及结果

4.1 解算分析

利用时序InSAR结果可得整个研究区域视线向的年均形变速率，Sentinel-1A升轨卫星影像年均形变速率图(图10)，正值表示靠近雷达方向，负值表示远离雷达方向。潮阳区整体形变量较小，大部分区域形变速率处于-5mm/a～20mm/a，处于平稳状态，局部区域存在形变，这和主城区建成区范围广、人类工程活动有显著关系[11]。雷达信号对水体不敏感，干涉效果差，表现为无地表形变，植被和建设用地均可表现出抬升和沉降现象。山区部分区域植被茂盛，造成失相干现象严重。沉降高值区域主要分布于东北部和南西部城镇建设区域以及地表起伏较大区域，沿道路分布。从地质灾害识别的角度出发，居民地及其周边地形起伏区出现明显地表形变异常，在遥感解译时可以重点考虑由于切坡建房引起的斜坡地质灾害。

4.2 解算结果

据2017年1月1日至2022年8月15日168景Sentinel-1影像D-InSAR与SBAS-InSAR处理，按照年形变速率大于20mm/a为强变形区，年形变速率10mm/a～20mm/a为中等强度变形区，年形变速率小于10mm/a为弱变形区。共解译73处疑似地表形变点。大致分为地表开挖(工程施工)、地表侵蚀与地表沉降与斜坡变形三类，主要以地表沉降为主；沉降区多为平原区河流、建成区周围，软土分布广泛，受地表水与地下水影响严重，这类变形往往在地表不出现明显变形，而以区域变形为主。

本次工作识别出4处斜坡变形，其中2处崩塌，光学影像中可见变形破坏迹象；2处滑坡，变形区与地形较为吻合，地表存在一定变形破坏迹象；其余点位推测均为软土受地表水与地下水影响形成的区域地表沉降。其中1处存在持续变形迹象，2017年1月至2021年8月存在持续变形，但存在一定季节性特征，约5.5年累计形变量为116mm；崩塌区季节变形明显，在夏季变形速度增大，曲线波动较大，说明与降雨相关性较高。

5 结 论

本次工作主要采用Sentinel-1进行地表形变探测，由于Sentinel-1为C波段20m分辨率，其分辨率较低，对工程建设开挖导致的地表形变探测效果较好；对于小型滑坡探测能力存在一定限制，而且山区植被茂盛，C波段卫星对植被覆盖区无法有效性探测。建议采用L波段ALOS-2与机载LiDAR相结合对小型地质灾害隐患点进行探测，提高准确性。