D-InSAR 技术与SBAS-InSAR 技术在矿区地面沉降监测中的应用研究

王 鹏，程海强，王 翔

（1.山东省三河口矿业有限责任公司地质测量科，山东 济宁 277600；2.枣庄矿业（集团）有限责任公司地质测量部，山东 枣庄 277000；3.山东建筑大学测绘地理信息学院，山东 济南 250101）

0 引言

煤炭资源作为我国储量丰富的矿产资源，在推动我国经济快速发展的同时，也导致矿区范围内的地面发生不同程度的沉降，对其上覆房屋、桥梁以及水利设施等建（构）筑物造成不同程度的损毁，因此，对矿区沉降情况进行监测具有重要意义[1-2]。传统地面沉降监测是利用水准测量、GNSS 测量等方法进行，不仅工作量大、成本高，而且无法获取大范围、全沉降盆地的沉降信息[3]。InSAR 技术是利用主动微波成像技术获取覆盖研究区重复轨道的SAR 影像，结合差分干涉测量技术获取地面沉降信息，具有全天时、全天候、大范围、高精度等优点，克服了传统监测方法的不足[4]。但D-InSAR 技术监测精度受时空去相干和大气延迟等因素影响，进行长时间序列地面沉降监测时并不能提供可靠的沉降信息。为克服相关问题，InSAR 技术通过对长时间跨度内多景SAR 影像联合处理，在提高地面沉降监测精度的同时，提供长时间序列的地面沉降监测结果，使得该技术在矿区沉降监测中得到广泛应用[5-6]。目前，已有不少学者基于不同的InSAR技术对矿区沉降情况进行了监测与分析[7]。冷红伟等[8]采用一种融合多视和全分辨率的D-InSAR 方法对济宁矿区下沉盆地进行了沉降监测，结果表明2017 年2 月17 日—3 月11 日期间，研究区内最大沉降出现在济宁市新驿煤矿，最大沉降量达到-74 mm，最大沉降速率达到-3.36 mm/d。张树衡等[9]集成D-InSAR 技术和SBAS-InSAR 技术对淮南市丁集煤矿进行了地面沉降监测，有效解决了SBAS-InSAR 技术在矿区沉降监测中沉降漏斗边缘精度低和部分区域像元空洞的问题，获得了连续高精度的矿区地面沉降结果。

为验证D-InSAR 技术和SBAS-InSAR 技术在矿区地面沉降中的监测能力，本文以济宁市某矿区为研究区域，收集了2022 年11 月2 日—2023 年8 月29 日期间共26 景Sentinel-1A 影像，分别采用D-InSAR 技术和SBAS-InSAR 技术进行数据处理，获取了该时段内研究区地面沉降监测结果，并结合研究区开采工作面内的实测水准数据对其精度进行验证，对比分析了两种技术在矿区地面沉降监测中的应用能力。

1 原理与方法

1.1 D-InSAR 技术基本原理

D-InSAR 技术通过两次重复轨道观测获取同一地区不同时刻的两幅SAR 影像，结合外部DEM 数据去除干涉图中的地形相位而得到差分干涉图，进而获取研究区范围内地面沉降监测结果情况[10]，主要原理见式（1）。

式中：φflat为平地相位；φtop为地形引起的相位；φdef为视线方向的形变引起的相位；φatm为大气延迟相位；φnoise为噪声相位；k为整周模糊度。其中，φflat可根据基线进行估算[11]，φtop可根据已有的DEM 估算，若忽略大气相位和噪声相位的影响，通过将估算出的φflat和 φtop去除，则可估算出地面沉降后的干涉相位 φdef，见式（2）。

式中：R1为主影像斜距；R2为辅影像斜距；Δr为沉降导致的地面点位移；等式右侧第一部分-（R1-R2）为形变前的干涉相位；第二部分Δr为由地面沉降引起的形变相位。

1.2 SBAS-InSAR 技术基本原理

SBAS-InSAR 技术通过设定时间-空间基线阈值，将满足阈值要求的时间序列SAR 影像数据分为多个具有较短时空基线的相互独立的干涉子集，以保证干涉对的相干性，进而得到更准确的地面沉降监测结果[12]。

假设SAR 卫星在研究区内获得了不同时段的N+1 幅SAR 图像，利用设定的时间和空间阈值可以获得k幅差分干涉图，则第k幅影像在方位-距离像元坐标系（x,r）中的差分干涉相位 δφk（x，y）见式（3）。

式中：φ（tb,x,r）、φ（ta,x,r）为tb时刻、ta时刻的相位；d（ta,x,r）、d（tb,x,r）为tb时刻、ta时刻视线向累计沉降量。去除地形残余相位[13]、大气延迟相位[14]以及各种噪声相位[15]后，地面沉降的平均速率vt可以表示为式（4）。

则相位为式（5）。

式中：EJ为主影像获取时间；SJ为从影像获取时间；vk为k时刻像元形变速率。

由此，定义A（j,k）=tk-tk-1，且A是一个m×n的秩亏矩阵，得到矩阵方程，见式（6）。

通过奇异值分解法和最小二乘法[16]可求出平均沉降速率相位值，进而可以获取地面沉降量。本文数据处理流程如图1 所示。

图1 数据处理流程图Fig.1 Flow chart of data processing

2 研究区概况及数据源

济宁市位于鲁西南地区，处在黄淮海平原与鲁中南山地交接地带，地理坐标介于东经115°54′～117°06′，北纬34°25′～35°55′，地势东高西低，地形以平原洼地为主。济宁市矿产资源极其丰富，已发现和探明的矿产储量达70 多种，以煤炭为主，含煤面积占全市面积的45%，经探测，济宁市煤炭储量达260 亿t，占山东省煤炭总储量的50%，主要分布在兖州区、曲阜市和微山县等地区，是全国重点建设的14 个亿吨级大型煤炭基地之一[17]。本文研究区地理位置如图2 所示，该采矿工作面所处区域周边断层、褶曲一般发育，岩层产状变化不大，地质构造条件中等。工作面赋存标高-244.8～-308.3 m，平均采厚5.1 m；煤层倾角4°～12°，平均倾角7°左右；走向长1 070～1 120 m，倾向长102～226 m，面积22.69 万m2。工作面采煤方法为倾向长臂后退式采煤法，顶板管理方法为全部冒落法。

图2 研究区地理位置图Fig.2 Geographic location of the study area

本文选取覆盖研究区的26 景VV 极化Sentinel-1A数据开展实验分析，SAR 影像时间跨度为2022 年11 月2 日到2023 年8 月29 日，具体参数见表1。为去除地形效应的影响，本文收集了覆盖研究区范围的SRTM DEM 数据，该产品为美国航天航空局“航天飞机雷达地形测绘任务”获取的产品数据，空间分辨率为30 m，具有高精度、高覆盖率等优点。此外，为了准确获取SAR 影像轨道信息，本文收集了可以提供高精度的卫星位置、卫星速度信息及卫星姿态数据等信息的POD 轨道数据，该数据是目前常用的SAR 影像定轨数据，可用来处理数据中矫正卫星运动误差和大气延迟误差，精密轨道数据日期与Sentinel-1A 数据日期一一对应。

表1 Sentinel-1A 数据IW 成像模式基本参数Table 1 Basic parameters of IW imaging mode of Sentinel-1A data

3 数据处理

3.1 D-InSAR 技术

在D-InSAR 技术数据处理中，每间隔约36 d 选择一幅影像，在2022 年11 月2 日—2023 年8 月29日期间内共有9 幅影像。将相邻的两幅影像组合为干涉对，依次对所获得的8 个干涉对进行干涉处理。SAR 影像配准直接影响干涉图的质量，本文先基于轨道信息对影像进行粗配准，再基于频谱差异法进行精配准，然后通过不断迭代，最终达到0.001 像元的配准精度[18]。研究区地表主要被鱼塘、农田和村庄覆盖，当雷达波束照射到地面时，在表面光滑的水体区域会产生类似镜面反射的现象，而在植被覆盖区域会受到植被结构的影响发生多次的散射和反射，导致这些地物特征区域的目标后向散射系数较低，回波信号弱，相干性差且干涉相位噪声严重。为了抑制这些噪声对干涉相位的影响，以10∶2 的比例对SAR 影像进行多视处理[19]，并采用自适应滤波方法[20]对差分干涉图进行滤波处理。进一步通过设定相干性阈值，掩膜掉部分相干性较差的区域（主要为水体区域）。在这一过程中，虽然矿区沉降中心的大幅度地表形变也导致该区域的相干性降低，但是为了获取全面的矿区形变结果，没有对这一主要研究区域进行掩膜。采用最小费用流法[21]进行相位解缠，获取解缠后的差分干涉相位图并将雷达视线方向的解缠相位转换为形变量，通过地理编码将获取的形变结果转换到地理坐标系下，将8 个干涉对获取的形变结果导入ArcGIS 软件中，用栅格计算器将其进行叠加，获得相应时间段内的累积沉降量。D-InSAR技术获取的截至2023 年8 月29 日的累积沉降结果如图3 所示。

图3 D-InSAR 累积沉降图Fig.3 Cumulative subsidence results from D-InSAR

3.2 SBAS-InSAR 技术

选 取在2022 年11 月2 日—2023 年8 月9 日 期间内的共26 幅SAR 影像，时间间隔12 d。设定时空基线阈值，空间基线阈值最大为310 m，时间基线阈值为36 d，进行干涉对组合，共生成47 个干涉对。在SBAS-InSAR 技术数据处理过程中，与D-InSAR 技术对解缠相位的叠加计算不同，SBAS-InSAR 的关键步骤是对解缠相位进行时间序列分析，从干涉相位中剔除高程残差和大气相位。首先，需要提取每个高相干点对应的解缠相位、图像信息和高程等信息；然后，利用SVD 方法进行小基线集分析，解算出高程改正量和形变序列；最后，采用时空域滤波的方法从形变序列相位中分离出大气延迟相位，即可获得最终的时间序列形变相位。SBAS-InSAR 技术获取的累积沉降结果如图4 所示。

图4 SBAS-InSAR 累积沉降图Fig.4 Cumulative subsidence results from SBAS-InSAR

4 结果分析

4.1 沉降监测结果对比分析

为了对两种技术在矿区沉降监测应用中的可靠性和监测能力进行验证分析。以工作面区域为研究对象（图2 中右图框线所示），以D-InSAR 技术数据处理差分干涉时间为准，获取SBAS-InSAR 技术对应时间的累积沉降量，结果如图5 所示（每组子图中左侧为D-InSAR 结果，右侧为SBAS-InSAR 结果）。由图5 可知，沉降主要集中在该采矿工作面北侧，东经117°04′20″、北纬34°51′50″周围，初期2022 年12月20 日沉降量较小且均匀分散，随着采矿作业的进行，沉降量不断增大，沉降中心逐步显现，至2023 年8 月29 日，从沉降中心向周围扩散且沉降量逐渐减少的沉降漏斗基本形成。在这一过程中，D-InSAR 技术和SBAS-InSAR 技术的监测结果在沉降中心分布及沉降范围方面吻合很好。但是受限于D-InSAR 技术分期独立处理数据，并通过简单累加的方式获取累积沉降量的方法，每个独立干涉对所包含的误差都会体现在最终的累积沉降结果中，致使其累积沉降结果中存在较多误差，在图5 中体现为沉降区域的空间连续性较差。而SBAS-InSAR 技术得益于其时间序列处理能够较好去除大气延迟等误差的特点，得到的累积沉降结果受到的误差影响更小，在空间上的连续性更好，这有益于沉降区域的读取及其对时空演变特征的分析。

进一步对两种技术获取的研究区地面沉降信息进行统计分析（图6 和图7）。由图6 可知，当监测时间较短，小于半年时，SBAS-InSAR 技术监测得到的最大累积沉降量基本大于D-InSAR 技术监测得到的结果；当监测时间更长时，D-InSAR 技术监测得到的最大累积沉降量明显超过了SBAS-InSAR 方法监测到的结果。截至2023 年8 月29 日，D-InSAR 技术监测得到的研究区内最大累积沉降量为69 mm，SBASInSAR 技术监测得到的最大沉降量为59 mm。由图7 可知，在监测初期（2022 年12 月20 日），研究区内没有监测到超过10 mm 的沉降量，在随后的9 个月中，沉降超过10 mm 的区域面积不断增大，截至2023 年8 月29 日，D-InSAR 技术和SBAS-InSAR 技术监测到累积沉降量超过10 mm 的面积分别为0.60 km2和0.87 km2，其中，沉降量大于30 mm 的面积分别为0.10 km2和0.17 km2。在整个监测过程中，SBAS-InSAR 技术识别的沉降量大于10 mm 的区域面积随时间推移而稳定增大，而D-InSAR 技术却存在一定的波动。此外，SBAS-InSAR 技术监测结果中，受植被和水体等地物影响相干性降低而未识别的区域面积为0.13 km2，而D-InSAR 技术则由开始的0.22 km2逐渐增加至0.32 km2。这是由于SBAS-InSAR 技术数据处理中通过小基线集技术获取了更多不同时空基线的干涉对，在后续的时间序列解算中这些在时间维度上有一定重叠的干涉对中的相位信息能够相互补充，弥补了D-InSAR 技术中某一时间间隔内仅依靠单个干涉对获取地表形变相位的不足，得到了更丰富的地面沉降信息。

图6 研究区最大累积沉降量柱状图Fig.6 Histogram of maximum cumulative subsidence in the study area

图7 研究区各阶段沉降面积变化柱状图Fig.7 Histogram of subsidence area change in the study area for different stages

4.2 沉降监测结果精度验证

为了验证两种技术获取的研究区内累积沉降量精度，本文对研究区开采情况进行了调查。该工作面于2022 年11 月开始开采，并在2022 年11 月4 日—2023 年8 月27 日期间进行了8 次水准测量，获取了整个开采过程中主要区域的地面沉降情况。结合精密水准数据，分别对两种技术获取的地面沉降结果进行精度分析，绘制了8 个水准点（图5 中依次为B01、B04、B08、B12、B16、B20、B24、B28），由D-InSAR、SBAS-InSAR 和水准测量三种方法获取的累积形变量折线图，如图8 所示。

图8 水准点累积沉降量折线图Fig.8 Line graph of cumulative subsidence at the level

由图8 可知，B01 号水准点、B04 号水准点、B08号水准点、B12 号水准点和B28 号水准点位于沉降漏斗边缘，监测期间内最大累积沉降量不超过50 mm。两种技术监测到的地面沉降趋势与水准监测结果总体吻合，其中，SBAS-InSAR 技术监测得到的沉降曲线在不同时刻变化相对平缓，但在沉降量级较大区域的变化特征表现不明显；D-InSAR 技术监测得到的沉降曲线变化幅度较大，但其对大量级沉降的探测能力比SBAS-InSAR 技术更敏感。这一现象与两种不同技术形变量获取的方法有关，其中，SBAS-InSAR技术是通过时间序列解算方法去除大气效应等误差影响，提高地面沉降解算结果的精度；而D-InSAR技术的叠加解算方法能够更为明显地反映真实地面沉降的形变特征，但在数据处理过程中往往无法消除如大气效应等误差的影响，进而影响其形变监测精度[22]。B16 号水准点、B20 号水准点和B24 号水准点位于沉降漏斗中心，监测期间内最大累积沉降量均超过300 mm，最大达到1 822 mm。根据InSAR 技术可监测形变梯度理论[23-24]，受雷达数据分辨率和波长等因素的影响，当地面沉降量过大时（超过InSAR技术可监测形变梯度范围），两种技术获取的沉降监测结果不可靠。

5 结论

本文基于2022 年11 月2 日—2023 年8 月29 日期间的26 景Sentinel-1 影像数据，分别通过D-InSAR和SBAS-InSAR 两种技术获取了该阶段内济宁市某矿区的地面沉降信息，并结合煤矿开采工作面的水准测量数据对两种技术得到的监测结果进行分析，得出结论如下所述。

1）在长期的矿区地面沉降监测过程中，D-InSAR技术监测到的累积沉降量大于SBAS-InSAR 技术，但识别到的沉降面积小于SBAS-InSAR 技术。在研究区内，D-InSAR 技术和SBAS-InSAR 技术监测到的最大累积沉降量分别为69 mm 和59 mm，识别出沉降量大于10 mm 的区域面积分别为0.60 km2和0.87 km2，其中，沉降量大于30 mm 的区域面积分别为0.10 km2和0.17 km2。

2）两种技术监测到的沉降区域分布基本一致，D-InSAR 技术容易受到误差影响，得到的沉降图形在空间上较为离散，沉降量波动较大；而SBAS-InSAR技术得到的沉降结果在空间上平滑连续，沉降量相对稳定。

3）由水准测量数据验证可得，两种技术得到的监测结果与真实沉降情况基本吻合，D-InSAR 技术对较大沉降量的监测更为敏感。在沉降漏斗中心地面沉降梯度较大的区域，地面沉降量超出了InSAR 技术可监测范围，两种技术的监测精度都大幅下降。