基于相干系数的DInSAR地表形变信息提取方法

乐 颖，夏元平，刘媛媛，胡 昊

(东华理工大学测绘工程学院，江西 南昌 330013)

0 引言

20世纪70年代左右，合成孔径雷达干涉(interferometric synthetic aperture radar，InSAR)技术出现，Rogers A E等于1969年首次利用 InSAR 技术研究了金星表面的地形特征[1]。随后，Gabriel等于1989年首次提出合成孔径雷达差分干涉测量(Differential InSAR，DInSAR)技术[2]。该技术是基于InSAR技术提取地表形变信息的一门极具潜力的新兴技术，由于在一定程度上克服了光照和天气的影响，可全天时全天候监测，具有高精度、大区域、快速等优势，因此受到了学者的广泛关注[3]。目前，DInSAR技术精度甚至可达毫米级，己经成为地震、火山、滑坡、地表形变信息等应用领域的主要监测手段[4]。

高分辨率SAR卫星的快速发展，使DInSAR技术已经成为目前获取地表形变信息的研究热点[5]。文献[6]比较了两种配准方法的特点，并将DInSAR技术应用于神木矿区的沉陷监测，获取最新的矿区沉陷信息及对应的时间序列。文献[7]利用DInSAR技术监测吕梁山区某矿区地表，通过两个剖面数据分析地表沉陷规律，表明能够得到对应的形变信息。文献[8]将DInSAR技术运用于陕西彬长矿区，并与实际情况相结合，有效获取了矿区的形变信息。文献[9]基于SAR 影像的振幅信息提出了一种改进的PS探测方法，通过对上海市局部区域进行研究，准确地提取了研究区域的形变信息。文献[10]基于DInSAR技术，结合GPS水平形变结果，可快速获得研究区地表三维形变，为后期的研究提供有效的参考。文献[11]针对不均匀地表形变问题，利用TS-DInSAR(time series DInSAR)技术对Sentinel-1A数据进行研究，结果表明该方法对大范围地表形变监测具有可靠的精度。然而，以上文献为了校正卫星轨道和相位偏移，去除处理过程中残余相位，在运用DInSAR技术处理影像时凭借人为主观经验选取控制点进行轨道精炼，这将影响相位信息转化为形变信息的准确性。若研究区域缺少地面控制点信息时，随意选取控制点将影响变形信息结果提取的精度。本文针对此问题，提出基于相干系数的DInSAR地表形变信息提取方法。

1 DInSAR技术

1.1 基本原理

DInSAR形变监测技术是经由InSAR技术衍生发展而来，主要包括二轨、三轨与四轨三种。其中，“二轨法”差分干涉技术运用领域最广，其表达式如下[12]：

φ=φflat+φtopo+φdefo+φatm+φnoise

(1)

式(1)中，φflat为平地相位；φtopo为地形起伏引起的地形相位，可以将地形信息恢复；φdefo为地表发生形变造成的相位；φatm为大气发生延迟的大气相位；φnoise为观测时产生的噪声相位[12]。简而言之，就是将式(1)中φflat、φtopo、φatm与φnoise等相位去除，最后只保留φdefo形变相位的技术[13]。

1.2 数据处理流程图

由于雷达数据与光学影像不同，属于主动遥感，影像包含的是振幅、相位和极化等多种信息。目前，常用的技术方法主要有合成孔径雷达干涉、极化分析、幅度追踪、层析建模和立体量测等五种方法，而其中合成孔径雷达干涉方法是最常用也是运用范围最广的技术。合成孔径雷达差分干涉通过处理覆盖同一地区不同时刻获取的两幅SAR影像的相位信息，能够测量地表厘米级的形变。在DInSAR数据处理过程中，主要包括：影像配准、干涉图生成、干涉图自适应滤波和相干性计算、最小费用流相位解缠、控制点选取、轨道精炼、重去平、相位转形变及地理编码等，每个数据处理步骤过程中产生的误差对干涉结果的精度和可靠性都具有显著影响，具体的处理流程见图1。DInSAR地表形变监测的数据准备主要包括SAR数据选取，外部DEM和精密轨道数据准备。其中，外部DEM是为了利用其模拟地形起伏相位分量，精密轨道数据则可以有效削弱SAR卫星轨道误差导致的残余参考相位。

图1 DInSAR数据处理流程图Fig.1 DInSAR data processing flow chart

2 基于相干性系数的 DInSAR 地表形变信息提取方法

2.1 干涉相干性测度

相干图主要用于体现两幅影像的相关程度，可以用于评价SAR影像精确配准后所得到干涉条纹图的质量，主要是通过计算干涉条纹图的相干系数实现。相干系数最初由普拉蒂于1993年提出，相干系数的取值范围为[0,1]，值越接近1说明相干性越高，其理论模型表达式如下[14]：

(2)

式(2)中，E[·]表示数学期望，u*表示共轭复数，u1与u2为主副图像的信号。在实际的InSAR处理过程中，在同一像素中无法得到计算所需数量的采样值。因此，根据式(2)的理论模型，基于SAR影像复数数据计算获得相干系数的标准表达式为[15]：

(3)

式(3)中，u1(n,m)，u2(n,m)分别表示主、副影像数据块内某个坐标(n,m)处的复数值；|•|2表示对应的二阶范数；M与N表示计算相干性的数据块尺寸大小；m与n表示数据块内对应的行列号。

2.2 轨道精炼控制点的选取

利用DInSAR技术获取地表形变信息时，引入稳定的地面控制点可以有效估计残留相位并去除平地效应，从而提高结果的可靠性。一般而言，选取控制点要远离形变区，认为形变为0，陡峭的地形区域和有残余地形相位区域，地形起伏大的山区，最好是选择山谷底部的平地区域[16]。为了更加快速准确地选择轨道精炼的控制点，本文基于相干系数阈值选取所需的控制点，相干性系数图的像元分布情况见图2。

从图2可看出，2015—2020年干涉对像元个数在相干性系数为0.5时达到了峰值，2015、2017和2020年的整体相干系数数值比其他年份高。由于相干系数图上的像元个数数量众多，选取轨道精炼时需要的像元数量很少。相干性阈值设置过低会导致相关性高像元数目剧增，大大增加数据计算的时间，从而增加后续筛选工作的难度。因此，在保证高相干性又考虑数据量大小的基础上，通过反复对比像元个数进行设置合适的阈值，最终将相干系数阈值设定为0.98，每年干涉相干性系数的像元分布情况见表1。

图2 2015—2020年相干性系数图像元分布情况Fig. 2 Distribution of image elements of coherence coefficient from 2015 to 2020

表1 每年干涉相干性系数的像元分布情况

2.3 控制点选取验证

为了更系统地验证本文方法选取控制点的可靠性，首先统计2015—2020年期间每年干涉相干性系数图符合要求的像元数量，然后根据相干系数阈值分别选取25～35个轨道精炼控制点，最后计算各点绝对残余误差，其计算公式为：

(4)

式(4)中，γ是干涉相干性测度，H是高程模糊度。

以2019年为例，对比分析DInSAR方法与本文方法用于选取轨道精炼控制点的绝对残余误差，结果见图3。本文方法用于选取轨道精炼控制点整体的绝对残余误差较小，只有个别点的绝对残余误差偏大，误差符合要求；而DInSAR方法选取控制点的绝对残余误差整体偏大，后续需要筛选绝对残余误差较大的点。此外，文中为了对比分析DInSAR方法与本文方法用于选取轨道精炼控制点的准确性，引入了标准差指标，并计算得到两种方法绝对残余误差的标准差分别为1.33和0.36，而指标值越小说明准确性越高。因此，在缺少地面控制点信息或研究区形变先验认识的情况下，若需要选择可靠的地面控制点，可借助相干性系数图进行剔除低相干性点。

图3 2019年选取轨道精炼控制点绝对残余误差Fig.3 Absolute residual error at the control point of orbit refining in 2019

3 实例验证

3.1 研究区域概况

赣州市的稀土储量占全国最多，而定南县是赣州市离子型稀土主产县之一，足以表明定南县稀土矿矿区扮演着重要的作用。基于此，本文以定南县作为研究区域，其地理坐标为东经114°47′49″E～115°22′48″E，北纬24°33′37″N～ 25°03′21″N，总面积为1 321 km2。县区内分布众多河流，主要隶属于定南水和桃江两水系[17]。同时，定南县地处南岭与武夷成矿带的中间部位，故而拥有良好的成矿条件，矿产资源非常丰富。在矿山开采过程中，露天矿坑开挖会直接破坏地表植被，严重影响着当地的可持续发展和生态功能恢复。因此，研究并提取定南县的地表形变信息尤为重要[18]。

3.2 实验数据

由于雷达数据与光学影像不同，属于主动遥感，影像包含的是振幅、相位和极化等多种信息。目前，合成孔径雷达干涉方法是最常用也是运用范围最广的技术。本文使用的雷达影像数据为C波段的Sentinel-1A卫星的SLC数据，选取了时间跨度为2015—2020年期间每年10月份共12景的影像数据，具体雷达影像数据的基本参数见表2。

表2 雷达影像数据的基本参数

3.3 结果分析

通过对12景覆盖研究区的SAR影像数据进行DInSAR处理，借助相干系数阈值获取相干性高的像元，并将其作为轨道精炼的控制点引入后续的处理中，最终可得到赣州市定南县乡镇在2015—2020年的地表形变信息，其形变结果图如图4所示。监测结果表明：研究区在监测时段内地表形变变化相对稳定，整体形变量主要集中在-0.07～0.1 m，部分区域出现了一定范围的地表沉降，沉降量集中在-0.11～0.11 m，定南县在2017、2019和2020年内的地表形变较2015、2016和2018年更稳定。

定南县共有7个镇，分别为鹅公镇、黄香镇、岿美山镇、老城镇、历市镇、岭北镇和天九镇。根据雷达SAR影像数据处理得到的结果，得到定南县各镇每年的地表形变变化情况，所得结果如图5所示。从图5的监测结果可以得出：定南县各个乡镇在2017、2018、2019和2020年内变化趋势基本一致，2015年形变量较大的乡镇为岭北镇，形变量达到了-25.26 mm；2016年形变量较大的乡镇为岭北镇、黄香镇和岿美山镇，形变量分别为5.26、5.26和26.29 mm；整个研究区域在监测时段内最大沉降量和最大上升量分别为-84.82和84.816 mm，分别位于岭北镇和鹅公镇。经过查找资料得知，赣州市定南县内具有丰富的离子型吸附性稀土矿矿产资源，且这些矿区位置大多位于岭北镇，其中包含了木子山、甲子背、长坑尾等稀土矿区[17]。由于岭北镇的平均形变量与其他乡镇的平均形变量相差-20.20～-10.10 mm，稀土矿的开采可能会引起少量的地表形变。

图4 2015年-2020年研究区域地表形变结果图Fig.4 Results of surface deformation in the study area from 2015 to 2020

图5 2015—2020年定南县乡镇地表形变信息Fig.5 Surface deformation information of towns in dingnan county from 2015 to 2020

为了探寻稀土矿的开采与地表形变之间的关系，并获取稀土矿开采区域的地表形变细节信息和地表形变发展趋势信息，本文结合稀土矿权在已有矿权范围内寻找了一些相关性高且能代表形变规律的特征点，特征点选取分布见图6，最后对特征点进行时间序列分析。

分别计算选取的30个特征点在2015—2020年间的地表形变量，图7绘制了地表形变细节信息和发展趋势。

图6 30个特征点选取分布图Fig.6 Distribution of 30 feature points

图7 30个特征点的地表形变情况Fig.7 Surface Deformation of 30 feature points

从图7中不难发现，选取的30个特征点在同年份的变化趋势基本一致，在2015年的地表形变发展整体为下降趋势，沉降量范围在-25.5～-15.5 mm，结果表明这些稀土矿区均存在开采现象；2016年和2017年的地表形变在0附近波动，地表变化幅度分别为-5.26～15.77 mm和-12.82～4.27 mm，由于这段时间国家开始加大对非法开采行为的打击力度，使得一些稀土私矿被查封；2018年整体呈下降趋势，部分点为向上抬升，地表形变在-32.984～14.136 mm之间，赣州地区对稀土矿山环境进行了综合治理以解决历史遗留问题，将稀土矿区平整为建设用地投入到当地经济发展中去；2019年和2020年的地表变化幅度分别为-26.03～18.07 mm和-3.81～3.81 mm，说明赣州定南地区开展的废弃稀土矿山治理取得效果显著，对当地的生态环境有了明显的修复。监测结果表明：由于离子型吸附性稀土矿的矿区范围不大，而且分布稀疏，从而导致稀土矿的开采引起的地表形变较缓慢。由于稀土矿的开采严重破坏了生态环境，江西赣州自2013年来不断推进稀土产业整顿，从加大稀土矿非法开采行为打击力度到对生态环境进行综合治理，很大程度上促进了矿山生态环境保护与稀土资源开发利用的协调发展。

4 结论

本文提出了基于相干性系数的 DInSAR地表形变信息提取方法。该方法通过设置一定的相干系数阈值，筛选出适合的高相干性像元作为控制点进行轨道精炼，最后提取出地表形变信息。以赣州市定南县作为研究区域，借助 2015—2020 年间的Sentinel-1A卫星数据，利用相干性系数形变信息提取方法，监测出研究区范围内-25.5～18.07 mm的地表形变信息。实验结果表明，该方法相较于常规的DInSAR方法，用于提取控制点的标准差提高了72.9%，定量表明了基于相干性系数的DInSAR方法用于选取控制点的准确性更高，能够实现对研究区域地表形变信息的动态监测。由于本文采用的数据量不足，研究过程中没有考虑大气因素的影响，导致得到的结果仍表现出一定的偏差。将雷达影像数据与光学遥感影像相结合，并与外业水准测量结果进行对比，将是下一步要进行的工作。