基于双轨D-InSAR技术的矿区沉降监测

朱文德，王长栓

（1.广西壮族自治区地理信息测绘院，广西 柳州 545006）

D-InSAR 是以合成孔径雷达复数据提取的相位信息为信息源，利用雷达图像的相位差信息来精确测量出图像上每一点的三维位置和提取地面目标点微小地形变化的技术[1]。该技术不同于水准、GPS 测量等传统的地表沉降监测方法，其具有覆盖范围广、观测连续、高效率、高分辨率、低成本、无须布设大量观测点等优点。在地震位移测量、监测火山运动、冰川漂移、地表沉降及山体滑坡等方面展现了其独特的优越性。1989 年，Gabriel[2]等首次利用D-InSAR 技术实现cm 级的地表形变监测。1993 年，Massonnet[3-4]等在Nature 上发表了震惊国际地震界的成果（利用D-InSAR 技术成功获取了1992 年Landers 地区7.3 级大地震的同震形变场，与GPS 测量所获得的结果吻合度非常高）。二十世纪九十年代以来，国内的众多学者也在合成孔径雷达干涉测量领域取得显著的成果。2001 年单新建[5]等基于InSAR 技术提取了西藏玛尼地区的DEM，所获结果优于1∶20 万DEM；2002 年利用D-InSAR 技术通过三通差分干涉处理获取西藏玛尼地震同震形变场[6]。2001 年，刘国祥[7]等利用D-InSAR 技术获取了香港赤腊角机场近乎1 a 内的非均匀沉降场，所获结果与离散水准监测结果相比吻合较好。2004 年，葛大庆[8]等将合成孔径雷达干涉测量技术应用于矿山沉降监测，为煤矿区地表时空演变过程研究和开采沉陷实时动态监测提供了新的技术方法。这些研究在很大程度上推动了我国D-InSAR 技术的发展。

论文阐述了利用双轨D-InSAR 技术进行矿区沉降监测的基本原理和技术方法，选用覆盖山东某矿区2017-02-17、2017-03-11、2017-04-02 三景雷达SAR 影像，获取该时段矿区的沉降信息，并针对沉降情况进行精细化分析。

1 D-InSAR技术原理及方法

1.1 D-InSAR 测量地表形变基本原理

D-InSAR 技术是在单天线模式下，进行重轨干涉测量，即在大致相同的轨道上，于2 个时刻获取同一地区的数据。技术方法包括二轨法（Two-Pass）、三轨法（Three-Pass）和四轨法（Four-Pass）[9]。本文以二轨法为例，简要介绍D-InSAR 测量地表形变的基本原理。

二轨法的基本思想是对同一地区变形前、变形后的两幅SAR 图像进行干涉处理，再引入已知的外部DEM，模拟地形相位，与干涉结果进行差分以去除干涉图中的地形因素，获取地表形变信息。该方法所需的SAR 数据量小，且通过引入外部DEM 消除地形相位，减少数据处理的工作量。但外部DEM 本身存在的误差会传递到最终的形变结果中，此外，利用外部DEM 模拟的SAR 图像需要与主影像进行配准，配准的精度也会影响最终形变监测精度。二轨法监测地表形变的基本原理可以表示为[10]:

式中，φ d表示研究区域发生形变前后两幅SAR 影像生成的干涉图；φsim,t表示采用外部DEM 模拟生成的地形相位；ΔRtow即为去除地形因素后的形变结果。二轨法数据处理流程如图1 所示。

图1 二轨法数据处理流程

1.2 D-InSAR 沉降监测数据处理流程

1）复图像配准。复图像配准的主要目的是使主、辅影像的像素对应同一地面单元。基本操作包括同名点偏移量计算、偏移量拟合及辅影像重采样[11]。主要配准方法有:①相干系数法，选择搜索窗口后，逐一计算每个像元上的相干系数，取其最大值对应的点作为配准点；②平均波动函数法，计算每个像元的干涉相位变化梯度，计算其平均值，值最小的点为配准点。

2）生成干涉图及干涉质量分析。主、辅影像完成数据配准后，进行共轭相乘，生成干涉图，此时的相位是缠绕的，数值范围为[-π,π]。为确保最终结果的精度，需要保证干涉的质量，通常以相干系数作为评定干涉质量的指标:

式中，M×N代表目标窗口大小；u1(m,n)和u2(m,n)表示目标窗口中(m,n)位置上的2 个复型数据；u*2(m,n)表示u2(m,n)的共轭复数。相干系数的值越大，代表干涉的质量越好。

3）去平地效应及相位滤波。基线的存在致使没有高程变化的平地产生干涉相位，为了降低解缠难度，将每一像素初始相位中斜距方向上的平地相位剔除。此外，干涉相位会受到很多噪声的影响，噪声的存在会降低相位解缠的精度，甚至会对最终形变结果产生影响。因此，在相位解缠前必须进行滤波处理。

4）相位解缠。差分干涉图中的相位被缠绕在[-π,π] 范围内，必须进行相位解缠来获取反映地面高程的真实相位。常用的相位解缠方法有3 种:①Goldstein 枝切法；②最小费用流法；③加权解缠法。

5）形变图生成及地理编码。通过公式ΔRtow=λ/(4π×φreal)，将解缠后的真实相位φreal转换成雷达视线向上的形变ΔRtow。为了更为直观的分析研究区域的形变，通常将雷达视线向上的形变图进行地理编码，转换到地图投影坐标系下。

2 D-InSAR在矿区沉降监测中的应用实例

2.1 实验区概况

本文选用矿产资源丰富的山东某矿区作为实验区域，位置及范围如图2 所示。该矿区位于山东省西南部，含煤面积约4 826 km2，是全国重点开发的煤炭基地之一。煤炭开采会导致矿区地表出现沉降、房屋坍塌、道路裂缝等问题，为研究矿区的沉降情况，本文选取覆盖研究区域的三景SAR 影像，基本参数如表1 所示。

图2 研究区域示意图

表1 干涉像对参数情况列表

2.2 数据处理关键步骤及成果

引入外部90 m 精度的DEM 数据，采用二轨法对上述3 个相对进行数据配准及差分处理，得到差分干涉图（如图3 所示）。

由图3 可以看出干涉相对1 所得的干涉效果较清晰，沉降区域较明显；干涉相对2、3 存在明显的噪声，影响干涉结果。为过滤噪声，提取精确的相位信息，采用Goldstein 滤波对上述差分干涉图进行滤波处理，滤波后的干涉图如图4 所示。

图3 差分干涉图

滤波后的干涉图，去除了噪声的影响，图像更加清晰，标记位置的干涉条纹发生突变，条纹密集，色度跨越大，可见此区域由于煤矿长期开采，发生了地表沉降，为了获取真实相位值，进一步做相位解缠处理。利用式ΔRtow=λ/(4π×φreal)，将解缠后的真实相位φreal转换成雷达视线向上的形变ΔRtow，再通过地理编码转换为地图投影坐标系下，图5 为所得形变图。

需要说明的是图5 中的形变与图4 表示的区域是一致的，由于所采用的SAR 数据为升轨数据，所以未经地理编码的干涉图的经纬度与实际地理经纬度相反，经过地理编码后的形变图所表示的经纬度与实际地理经纬度一致。

图4 滤波后的干涉图

图5 形变结果

3 矿区沉降结果分析

由图5 可以得知，2017-02-17～2017-04-02 期间，研究区域内出现A、B、C、D 四处明显的沉降盆地，是由于煤矿的长期开采，导致地面出现沉降，且对周边地区造成一定程度的影响。其中D 处的煤矿沉降尤为严重，为更直观的分析矿区沉降情况，以该煤矿为例，作图示剖面线，进行精细化分析，如图6 所示。

图6 剖面线沉降图

通过对比2017-02-17 ～2017-03-11、2017-03-11～2017-04-02、2017-02-17～2017-04-02 三期成像的剖面图可以得到:D 煤矿由于长期开采及工作面间的相互影响，形成2 个沉降漏斗。3 个干涉相对均呈现这一特征，且沉降位置及整体沉降趋势一致。2017-02-17 ～2017-03-11 期间2 个沉降漏斗最大沉降量分别达到3.8 cm、4.7 cm，2017-03-11 ～2017-04-02 期间2 个沉降漏斗最大沉降量分别达到3.0 cm、4.8 cm，2017-02-17 ～2017-04-02 期间2 个沉降漏斗最大沉降量分别为3.0 cm、4.5 cm。干涉相对3 的沉降结果应为干涉相对1、2 之和，但由于时间基线过长，加上各种误差的影响，该时段失相干严重，沉降漏斗周边出现大量不连续的点，造成3 个成像期间沉降量有所差异，说明时间基线过长会降低双轨D-InSAR 的监测精度。

4 结 语

本文论述了基于双轨道D-InSAR 技术进行矿区沉降监测的基本原理和技术方法，采用2017-02-17、2017-03-11、2017-04-02 三景SAR 影像提取该时段矿区的沉降信息，并对沉降情况作进一步的精细化分析。结果表明，在2017-02-17、2017-03-11、2017-04-02成像期间，出现4 个明显的沉降区域，其中D 处煤矿沉降最为严重，形成2 个连续的沉降漏斗，并呈现沉降区域扩大，沉降情况加重的趋势。实验证明D-InSAR 技术可以反映可靠的矿区沉降情况，能够用于矿区沉降的实时监测，但由于外界条件制约，各种系统误差的影响，还需通过限制时间基线，结合水准数据或GPS 等测量结果，以提高成果的精度。