基于SAR 技术的开采沉陷全盆地分区形变提取方法

邹英杰，范洪冬，孙 叶，许耀宗

（中国矿业大学 自然资源部国土环境与灾害监测重点实验室，江苏 徐州 221116）

煤炭工业作为我国重要的基础工业，在带来巨大经济效益的同时，地下煤层开采活动也会引发地表人工建构筑物损坏和一系列的生态环境问题[1-2]。因此，对矿区地表进行精细化形变监测十分必要。常规地表形变监测手段，如水准测量和GNSS 等方法观测精度虽高，但存在工作量大、观测点位密度低且易被破坏等缺点。

合成孔径雷达差分干涉测量（DInSAR）技术具有全天时、全天候、空间分辨率高、覆盖区域广等优点[3]，其利用雷达视线向（LOS）相位信息能够监测到地表厘米级甚至亚厘米级的形变[4]；但是该技术受失相干影响较为严重，且受限于相位解缠的理论形变梯度阈值[5]，只能监测地表小量级形变。偏移跟踪（OT）方法基于SAR 影像的强度而非相位信息，因此其不受时空失相干和形变梯度的限制，多用于地震、冰川、滑坡、矿区等大量级形变监测领域[6-10]；OT方法的监测精度会受到匹配窗口大小的影响[11]，实现匹配窗口的自适应选取可有效提升OT 方法的监测精度。

近年来，DInSAR 和OT 相结合的方法在矿区地表形变监测方面已经取得了相应成果，通过设定相干系数阈值或掩膜大形变区是将地表形变分为大小2 个量级的主流方法[12-16]。在此基础上，相关学者根据先验知识和目视解译等方法划分出中等量级形变区域[11,17-19]，并对该区域进行不同的融合处理。上述划分形变区域方法的适用性和范围会受到现场实测数据和人为主观因素的影响，所采用的OT 方法均未考虑到匹配窗口的自适应选取。因此，通过相干系数和地表形变梯度划分出中等量级形变区域，联合该区域累积DInSAR 方法和自适应OT 方法的监测结果进行加权融合，实现对矿区沉陷盆地大、中、小不同量级形变的分级与精细监测。

1 开采沉陷盆地分区形变提取方法

1.1 累积DInSAR 方法原理

SAR 卫星通过向地面发射电磁波信号，接收机接收到地面的反射信号后，经处理得到SAR 单视复数影像，其中包含地物的辐射强度信息和相位信息。DInSAR 技术利用形变前后观测的相位进行干涉测量，其沿LOS 向路径差△r 可以用干涉相位的相位差表示，即：

式中：Φ 为干涉相位的相位差；Φflat、Φtopo分别为平地相位和地形相位，可以利用SAR 卫星成像参数和外部高精度DEM 进行去除；Φdef为地表形变引起的形变相位；Φnoise为其他因素引起的噪声相位。

SAR 卫星可以以一定的周期重复获取同一地区的时间序列影像，随着SAR 数据的积累，一系列时序InSAR 方法得到快速发展，如SBAS-InSAR、PS-InSAR、StaMPS-InSAR 等。累积DInSAR 方法作为1 种最简单的时序InSAR 技术，对一系列SAR 数据进行差分干涉处理，将差分处理后的解缠相位进行累积，表达式如下：

式中：Φsum为累积后的形变相位；Φiumw为第i 个差分干涉处理后的解缠相位；n 为SAR 序列影像的数量。

最后，可根据式（2）计算出成像期间地面点沿雷达LOS 向的累积形变量。

1.2 自适应OT 方法原理

传统OT 方法利用光学影像匹配中的互相关系匹配算法（NCC），对主影像和从影像按照一定窗口进行精确匹配，通过式（3）计算互相关系数和寻找互相关系数矩阵峰值的位置来确定2 幅影像在方位向和距离向上的像素级偏移量。为达到亚像素级的配准精度，可通过插值、曲面拟合等方法对相关系数进行拟合，将极值点的位置作为匹配的最佳位置。

式中：ρ（x，y）为像元点（x，y）的互相关系数；I1、I2分别为模板窗口和搜索窗口内像素点的像元强度；I¯1、I¯2分别为模板窗口和搜索窗口内所有像素点的像元强度平均值；（u，v）为像素点在从影像中相对于主影像的距离向和方位向偏移量。

自适应OT 算法旨在实现匹配窗口的自适应选取，提高OT 方法的监测精度。因此，需要寻找1 个最合适的匹配窗口，为了衡量不同窗口间的匹配结果的可信度，采用信噪比SNR 作为衡量标准，其定义如式（4）：

式中：ρmax为相关系数峰值；ρ¯为3×3 窗口内相关系数平均值。

SNR 值越大，代表相干系数峰值所在位置越突出，越有利于进行图像的精确匹配，将SNR 最大值的窗口大小作为OT 方法匹配的最佳窗口。自适应OT 方法通过寻找以初始步长为b（通常为8）的一定窗口区间［a，d］内的SNR 最大值的窗口大小a+bi，并继续在［a+b（i-1）a+b（i+1）］窗口区间内以步长为b/2 寻找新的SNR 最大值的窗口大小，以此类推，直到找到最终的最佳匹配窗口；在窗口大小自适应选取的基础上还需完成OT 方法后续相关系数插值以及距离向与方位向亚像素级偏移量解算等步骤。自适应OT 方法的实现过程如图1，

图1 自适应OT 方法流程图Fig.1 Processing flow of adaptive OT

1.3 累积DInSAR 和自适应OT 联合加权方法原理

地下煤层开采过程中，上覆地表形变量级通常是从下沉中心到边界逐渐减小的，从前面的分析可以看出，累积DInSAR 方法可以测量下沉盆地边缘小尺度地表形变，自适应OT 方法可以测量下沉盆地中心部分大尺度地表形变。对累积DInSAR 和自适应OT 方法监测精度较低的中等量级形变可以利用二者得到的结果进行加权融合，实现对矿区多量级地表形变精细监测。

对于小量级形变区域，采用累积DInSAR 方法，由于相干系数越高代表地表形变的解算结果越可靠，多数文献及GAMMA 软件说明书在进行差分干涉处理时，一般取相干系数0.3 作为执行相位解缠的阈值[13，20]。对于不同分辨率SAR 影像的处理阈值并无明确说明，考虑到大柳塔矿52303 工作面地表形变大而导致利用相位解算地表形变误差大，在进行试验时，曾选取多个相干性阈值进行试验，当阈值为0.45 时，所划分的小量级形变区域效果最佳。

由于相邻像元形变梯度大于1/4 个波长的区域超出了单次相位解缠的极限要求，而矿区因开采导致的地表沉降超过2 m 的情况非常普遍，若要保证本研究累积DInSAR 方法能够参与地表形变解算，只能在矿区地表形变梯度小于（n-1）/4 个波长的区域。概率积分法作为矿区开采沉陷适用性最广的方法[21-22]，本研究根据地质采矿参数和概率积分法预计并预计出相邻像元梯度大于（n-1）/4 个波长的区域作为大量级形变区域，该区域采用自适应OT 方法进行处理。

对于其它区域则作为中等量级形变采用累积DInSAR 与自适应OT 方法按照一定的权值进行融合处理，表达式如下：

式中：DLOS为本方法最后得到的地表LOS 向形变信息；P1、P2分别为DInSAR 方法和自适应OT 方法所占权值，一般依据先验定权的方式或经验确定；C 为相干系数；Div 为概率积分法预计LOS 向形变梯度；ΣDDInSAR为累积DInSAR 监测结果；DOT为自适应OT 方法监测结果；n 为SAR 影像数量；λ 为雷达波长。

根据现场实测资料先验定权的方式，其权值确定的具体公式如下：

式中：σΣDInSAR、σOT分别为累积DInSAR 方法和自适应OT 方法相对于实测水准数据或GPS 数据的均方根误差，并将其作为监测精度。

2 试验数据概况

2.1 矿区概况

试验选取榆林大柳塔煤矿52303 工作面为研究区域。大柳塔煤矿地处陕西省榆林市神木市与内蒙古鄂尔多斯市伊金霍洛旗交界处，研究区地理位置及GPS 点位信息如图2，红色的圆点为地表移动观测GPS 点。该地平均海拔在1 250 m，属于北温带大陆性半干旱气候，冬季寒冷，夏季干热，春季多风，秋季凉爽。52303 工作面于2012 年12 月开始试采，工作面长度约301.5 m，推进长度约为4 443 m，平均采深250 m，平均采高6.6 m，煤层倾角1°～3°，采用走向综合机械化长壁采煤法，顶板管理方法为全部垮落法。作为1 座年产量超2 000 万t 的超大型现代化煤矿，长久的地下采煤活动造成地表不同程度的地表变形，危及矿区生产安全和周围人民生命财产安全。因此，对该矿区进行大范围、多量级的地表形变精细监测具有重要的意义。

图2 研究区地理位置及GPS 点位信息Fig.2 Location of the study area and the information of GPS points

2.2 SAR 数据概况

研究采用德国航天中心发射的TerraSAR-X 卫星高空间分辨率雷达影像，该卫星波长3.1 cm，卫星重访周期为11 d，其距离向和方位向空间分辨率分别为0.91 m 和1.82 m。为了便于试验精度评定，选取TerraSAR-X 卫星成像日期为2014-04-05—2014-05-08 日的4 景影像，现场实测GPS 数据日期为2014-04-06—2014-05-08 日，可认为SAR 卫星和GPS 二者为同期观测。所选取的SAR 影像对数据具体参数见表1。

表1 TerraSAR 影像对参数Table 1 TerraSAR image pairs parameters

3 试验结果及分析

3.1 数据处理结果

累积DInSAR 试验使用4 景时间序列SAR 影像做DInSAR 处理，在DInSAR 处理过程中引入外部STRM 分辨率为30 m 的DEM 数据去除地形相位。为了保证在进行差分干涉处理过程中起始解缠相位不产生相对变化，在每次DInSAR 处理过程中选择同1 个远离形变区域、保持高相干性的相位解缠起始点。将解缠后的形变相位进行时序叠加后，根据式（2）转换为LOS 向形变，累积DInSAR 方法得到的矿区沉陷盆地边缘小形变区域监测结果如图3。

图3 累积DInSAR 方法形变监测结果Fig.3 Deformation monitoring results of cumulative DInSAR

利用2014-04-05 和2014-05-08 日2 景影像进行自适应OT 处理，本研究自适应OT 方法设置初始窗口区间为［73-121］，初始步长b 为8。同时，根据现场实测资料以及矿区采煤地质条件，判断出该工作面地表LOS 向形变值应处于-0.5～4.5 m 的范围内，因此需要对自适应OT 方法监测结果内大于4.5 m 和小于-0.5 m 的形变值进行滤除处理。研究使用均值滤波对滤除后空洞进行填充，得到自适应OT方法在矿区沉陷盆地中心大量级形变亚像素级监测结果如图4。

图4 自适应OT 方法形变监测结果Fig.4 Deformation monitoring results of adaptive OT

研究提出的加权融合方法以相干系数0.45 和预计形变梯度0.023 为阈值，选择出中等量级形变区域。结合该区域现场实测GPS 数据分别对自适应OT 方法和累积DInSAR 方法进行精度验证，累积DInSAR 方法的监测精度为8.5 cm，自适应OT 方法的监测精度为8.2 cm，根据式（6）进行定权，确定中等量级形变区域累积DInSAR 方法和自适应OT 方法所占权值分别为0.48 和0.52。根据式（5）对累积DInSAR 方法和自适应OT 方法的监测结果进行加权融合，得到中等量级形变监测结果。

对3 种量级形变进行线性叠加后，得到矿区完整多量级地表形变场，其监测结果如图5。

图5 本研究形变监测结果Fig.5 Deformation monitoring results of the method in this paper

3.2 试验结果分析

结合52303 工作面开采工作进度和现场地表实测数据，此时52303 工作面开采工作已经进入最后阶段直至停采了一段时间，因此沉陷盆地形变中心靠近5-2煤辅运大巷。考虑到矿区开采过程中，地表移动具有一定的滞后性，此期间地表LOS 向形变最大值已超过3 m，所生成的地表沉陷盆地仍可以用本研究方法进行形变监测。通过图3 和图4 可以直观发现，累积DInSAR 所监测到的最大形变小于0.2 m，远小于矿区实际地表形变。自适应OT 方法可以弥补DInSAR 方法无法监测大梯度沉降矿区的缺陷，但是通过图4 可以看出，自适应OT 方法在矿区周边非形变区域出现了错误的监测结果，说明自适应OT 方法在小量级形变区域解算精度并不高。对于矿区大量级沉陷盆地，DInSAR 和OT 方法进行联合监测是非常有必要的，本研究所提出的融合累积DInSAR 和自适应OT 的方法可以有效地提取矿区沉陷全盆地多量级形变场。利用矿区52303 工作面倾向19 个地表GPS 观测点对本研究提出方法所监测的多量级形变结果进行精度分析，将二者结果进行对比，19 个地表观测站形变对比如图6。

图6 地表观测站LOS 向形变对比图Fig.6 Comparison of LOS deformation at the surface observation station

在中等量级形变区域，自适应OT 方法、累积DInSAR 方法和融合方法的监测精度对比见表2。对研究区域内3 种不同量级的形变监测结果进行精度分析，精度对比见表3。

表2 中等量级形变区方法精度分析表Table 2 Accuracy analysis table of medium order deformation region method

表3 研究区不同量级形变精度分析表Table 3 Analysis table of deformation accuracy of different orders in the study area

由表2 可知：采用累积DInSAR 方法，其均方根误差为0.085 m，平均绝对偏差为0.065 m，自适应OT 方法均方根误差为0.082 m，平均绝对误差为0.058 m。根据权值对2 种方法的监测结果进行融合后，将融合加权结果与现场实测结果相比，加权融合方法均方根误差为0.046 m，平均绝对误差为0.036 m，其监测精度明显高于单独使用累积DInSAR 或自适应OT 方法，证明本研究所提出的加权融合方法能够有效提高矿区沉陷盆地边缘中小量级形变监测的精度。

由表3 得到本方法提取的GPS 观测点形变值的均方根误差为0.108 m，最大绝对偏差值点号为Q14，为0.319 m。其中，自适应OT 方法均方根误差为0.154 m，平均绝对误差为0.124 m；累积DInSAR方法均方根误差为0.048 m，平均绝对误差为0.037 m。自适应OT 方法均方根误差最大，这是因为自适应OT 监测形变精度受到像元尺寸的影响，其均方根误差与距离向像素尺寸比值为16.48%。采用融合方法获得的监测精度相比于累积DInSAR 方法并没有显著提高，但是其在矿区地表下沉盆地边缘中小量级形变的监测范围要优于累积DInSAR 方法。

4 结 语

针对DInSAR 方法和OT 方法在中等量级形变监测精度较低的问题，提出1 种融合累积DInSAR和自适应OT 方法，实现对矿区多量级地表形变精细监测。利用4 景TerraSAR 影像分别采用累积DInSAR、自适应OT 以及二者加权融合方法对大柳塔煤矿52303 工作面开采沉陷盆地大量级、小量级和中等量级形变进行监测，得到矿区沉陷全盆地形变信息，与同期现场实测GPS 数据对比评价。

1）提出的融合方法监测结果与19 个GPS 点实测数据对比，证明了融合累积DInSAR 和自适应OT 的方法能够获取矿区完整多量级地表形变的有效性和可靠性，其监测精度优于0.108 m。

2）针对矿区地表形变盆地边缘中小量级形变，累积DInSAR 自适应OT 方法加权融合方法的均方根误差优于0.046 m，其监测精度明显高于单独使用累积DInSAR 或自适应OT 方法，相较于二者均方根误差，分别提高了46%和44%。

3）累积DInSAR 和OT 融合的方法能够得到矿区地表完整的形变场，但是该方法监测精度受到SAR影像分辨率和雷达波长的影响，待未来长波长、高空间分辨率SAR 卫星投入使用后，该方法的监测精度及中小层次形变监测范围将会进一步提高。