矿山开采后地表沉降边界的D-InSAR 监测与分析*

赵 恒，李亚星，杨可明*，窦 智

（1.淮北双龙矿业有限责任公司地测科，安徽 淮北235000；2.中国矿业大学（北京）地球科学与测绘工程学院，北京100083）

地下煤层开采会对其上覆的岩层结构产生破坏，进而引起地面塌陷与沉降变形，煤矿区的采动严重影响威胁了采区地表的土地、水体、建（构）筑物等的安全。为了保障煤矿区生产和生活安全的同时尽可能多回收煤炭资源，需要科学地确定开采沉降范围，以便明确资源开采位置和界定相关责任，降低煤矿企业的开采风险。

矿山开采引发地表沉降影响范围的常规分析主要是依据开采工作面的实际地质采矿条件、地表岩移观测数据以及岩移规律等资料，运用实测数据分析、岩层移动角几何图解、概率积分法（Probability Integral Method，PIM）预计和数值模拟计算等方法[1-4]，但这些常规方法存在较多且明显的不足，如：岩移观测站实测数据分析与几何图解主要体现在开采工作面的走向与倾向方向上，观测点数量与图解范围相当有限；PIM 与数值模拟方法虽然能大范围确定地表沉降边界，但不能反映不同的开采工作面存在不同的上覆岩层岩移与物理参数而引发实际沉降边界的不同变化，导致所确定的边界与实际位置存在一定甚至较大的偏差。而微波遥感的合成孔径雷达干涉测量（Interferometric Synthetic Aperture Radar，InSAR）技术并不需要矿山开采的地质采矿资料，而是基于侧视雷达电磁波的地面反射回波相位干涉处理技术来获取地表的沉降信息，常用技术有1989 年Gabriel 等首次提出的差分干涉测量（Differential InSAR，D-InSAR）[5]、2000 年Ferretti正式提出的永久散射体干涉测量（Permanent Scatterer InSAR，PS-InSAR）[6-7]和2002 年Berardino 提出的小基线集干涉测量（Small Baseline Subset InSAR，SBAS-InSAR）[8]等技术。InSAR 技术由于可以克服因地下开采条件复杂、地表移动变形，以及易受地下开采、覆岩岩性及结构、地下水等多种因素影响而难以准确界定地下开采所形成影响边界的局限，所以现今已成为矿山开采形变监测[9-10]、地质灾害分析与环境影响评价[11-13]等开采沉陷研究的前沿技术。

D-InSAR 是目前较为先进的一种地表形变监测方法，与常规方法相比在监测地表微小形变方面有着显著的优势，它不仅能以厘米量级甚至毫米量级监测地表形变，且能对大面积测区进行快速监测，从而得到采区主要形变区域与范围[14-15]。本文将结合安徽省淮北矿区某煤矿井下311 工作面煤层开采下保护地面重要建（构）筑物运营安全等实际需求，利用D-InSAR 干涉测量技术监测煤层工作面采动影响的地表沉降边界，以下沉10mm 的位置为确定边界的依据，并基于现场实测数据、SBAS-InSAR 技术分析结果和概率积分法变形预测等对所监测到的地表沉降边界进行验证，进而能为煤矿安全开采和最大程度回收地下煤炭资源提供科学的理论依据与技术保障。

1 理论与方法

1.1 D-InSAR 技术原理

D-InSAR 技术以其大范围、低成本、速度快、高精度等优势，已被广泛应用于各种地表形变监测中。如图1 所示，D-InSAR 的基本思想从两幅SAR 影像的干涉相位中分离出地表形变所引起的相位分量，从而对影像区域的地表沿雷达视线向地表形变进行制图。干涉相位是两幅配准过后的SAR 影像对应像元的相位差，其组成可以表示如下[16]：

图1 平地-高程相位函数模型

式中，φint表示干涉相位；φtop、φflat、φorb、φatm和φnoise分别表示高程相位、平地相位、轨道误差相位、大气相位和相位噪声；φdef表示形变相位。轨道误差可以通过精密定轨和模型拟合去除，噪声一般通过滤波来抑制。矿区沉降空间范围较小，大气相位常常可以忽略不计。高程相位和平地相位可通过函数模型来去除。

图1 中S1、S2分别表示主辅影像成像时的天线位置。B 为空间基线，B||和B⊥分别表示平行和垂直于雷达入射方向的基线分量，称为平行基线和垂直基线。α 为空间基线与水平方向的夹角，θ 为雷达入射角，Δθ 为高程引起的雷达入射角变化。R1、R2表示天线与目标点P 之间的斜距。R 表示P 点等斜距投影到参考椭球面的R1距离（R=R1）。h 为目标点P 的高程。根据图中几何关系和泰勒公式，有：

式中，λ 为波长。因此可以得出平地相位和高程相位的表达式：

形变相位φdef与地表视线向形变dlos之间的关系如下：

1.2 地表沉降边界分析方法

D-InSAR 是在InSAR 的基础上发展起来的监测地表微小形变的技术，在监测工矿区地表形变方面具有连续空间覆盖，高度自动化和高精度的优势。一般常用二轨法D-InSAR 处理流程，即采用外部的数字高程模型（Digital Elevation Model，DEM）数据模拟并去除地形相位。二轨法D-InSAR 的基本处理流程如图2 所示。

图2 二轨法D-InSAR 基本处理流程

1.3 概率积分法地表移动变形预计思想

地表移动与变形计算的概率积分法（PIM）是指把岩层移动看作服从统计规律的随机过程，从而将开采引起的地表下沉剖面表示成概率密度函数积分公式的预计方法。已有大量理论推导及实例验证表明，当地下开采工作面的顶板管理方法采用自然垮落时，概率积分法得到的地表沉降盆地稳定后的预计结果与实测数据具有较高的拟合优度和精度，是用于矿区开采沉陷预计的一种行之有效且广泛使用的方法。

2 结果与分析

2.1 研究区工作面概况与遥感数据源

2.1.1 311 工作面概述

图3 为研究区煤矿的311 工作面及其地表走向观测线的沉降监测点位示意图。311 工作面位于该矿井Ⅱ1 采区上山西翼，工作面可采走向长550m、面长240m。该工作面整体呈一单斜构造，煤层倾角8°～18°，平均倾角为10°，煤层厚度为2.6m～4.2m，平均采厚为3.0m。该工作面标高-160m～-280m，地面标高为+31m～+32.5m，2018 年2 月9 日开始回采，于2019 年7 月27 日开采结束。

图3 井下开采的311 工作面及地表走向观测线点位图

2.1.2 遥感数据源

所使用的SAR 数据为Sentinel-1A（哨兵-1A）IW 成像模式SLC 数据，此SAR 数据为C 波段数据，空间分辨率为5m×20m（距离向×方位向），其他详细参数见表1 所示。根据311 工作面开采的时间节点，获取了2018 年2月8 日至2019 年9 月7 日间的46 景哨兵1A 影像数据。

表1 所用SAR 数据主要参数

2.2 基于D-InSAR 技术的开采沉降边界分析

D-InSAR 处理的具体技术路线如图4 所示，所涉及的主要技术内容包括：（1）影像预处理。包括数据格式转换和目标区裁剪。主要是将原始影像数据处理为标准的SLC 数据，然后根据研究区范围的大小，对影像进行目标区裁剪。（2）基线估算。估算结果包括时间基线、空间基线、多普勒偏移、一个相位变化周期代表的高程变化等信息，主要目的是评估像对干涉质量。（3）生成干涉图。首先根据影像的精密星历计算影像中心偏移值，并设置合适的窗口大小以进行自动配准，精度可以达到1 个像元的水平；然后将主辅影像与已有的外部基准DEM 进行配准，再经过采样处理之后使配准精度达到0.001 个像元的水平；在此基础上根据配准好的主辅影像的相位信息进行干涉，生成干涉条纹图。（4）滤波和相干性计算。滤波处理主要目的是去除大气延迟和平地效应所引起的误差相位，同时得到相干系数图。（5）相位解缠。使影像相位信息与地形信息相对应。（6）控制点选择。控制点主要用于轨道精炼和重去平，常在去平后的干涉图上进行选择，选择应遵循以下原则：避免在有地形相位未去除和变化的区域选择控制点；选择相干性高的区域；控制点应分布于整个范围内；避免解缠错误的区域，如相位孤岛等。（7）轨道精炼和重去平。目的是估算和去除残余的恒定相位和解缠后还存在的相位坡道。（8）相位转形变和地理编码。解缠后的相位经过转换就可得到相对于解缠起始点的视线向形变信息。地理编码的目的是将所有的结果编码到所需的坐标系下。

图4 D-InSAR 技术路线图

基于获取的46 景哨兵1A 影像组成45 个干涉对，利用叠加二轨法D-InSAR 技术叠加了2018 年2 月8 日至2019 年9 月7 日的累积沉降，将每对D-InSAR 得到的沉降监测结果按时序叠加即可获得研究区内的时序沉降，以此可确定311 工作面开采下沉影响边界，如图5 所示。311 工作面中心区域出现下沉梯度过大导致的地表形变失相关，因此监测到的最大下沉为-302mm，其原因是因为淹水或下沉梯度过大区域，InSAR 监测失相干，无法监测实际最大下沉量。

2.3 验证分析与下沉边界的确定

以下沉值-10mm 作为开采沉陷盆地边界，联合D-InSAR、SBAS-InSAR 技术以及概率积分预计法，以多手段、多角度来综合权衡确定311 工作面南侧实际下沉边界。

2.3.1 基于D-InSAR 技术确定的沉降边界位置

如图5 所示，利用D-InSAR 技术监测到的下沉边界大约在II311 工作面走向观测线上的D013 点南侧、距II311 工作面采空区边缘的切眼位置平均距离为165m。

图5 D-InSAR 技术确定的311 工作面实际下沉边界

2.3.2 基于SBAS-InSAR 技术确定的沉降边界位置

基于2018 年2 月8 日至2020 年9 月7 日的46 景哨兵1A 影像，再利用SBAS-InSAR 技术对311 工作面开采影响的下沉边界进行监测，监测到的下沉边界大约在II311 工作面走向观测线上的D013 点、距II311 工作面采空区边缘的切眼位置平均距离为160m。

2.3.3 基于概率积分法确定的沉降边界位置

如图6 所示，可看出基于概率积分预计模型所获得的下沉边界在II311 工作面走向观测线上的D013 点、距II311 工作面采空区边缘的切眼位置平均距离为160m。

2.3.4 基于InSAR 及概率积分法预测的下沉边界综合确定

通过综合上述D-InSAR、SBAS-InSAR 技术与概率积分预计三种方法分析311 工作面实际下沉边界，得到如图6 所示多方法确定的311 工作面实际下沉边界位置，与实际观测线上的点位相符合。

图6 多方法综合的311 工作面南侧实际下沉边界对比图

3 结束语

工矿区地表形变监测过程复杂，其破坏及危害性不可避免。传统的变形监测技术费时费力，难以得到大范围内的整体形变结果，随着微波遥感的快速发展与应用，D-InSAR 和SBAS-InSAR 等技术已广泛应用于地表沉降动态监测。本研究主要是利用D-InSAR 技术对淮北矿区某煤矿311 工作面开采的采动影响边界进行确定研究，并通过概率积分法预计与SBAS-InSAR 干涉测量技术对D-InSAR 监测结果进行验证，提取到311 工作面开采后的沉降边界，同时根据311 工作面走向观测线的观测点位，在走向观测线上指定了沉降边界的具体位置，即在D013 测点处，该测点至311 工作面采空区切眼巷道的平面距离为162m。