0 引言

煤炭资源的开发对我国经济持续健康发展具有重要意义。但煤炭资源的大规模开发为各行各业带来巨大经济效益的同时，也给矿区的生态环境和土地利用造成了巨大影响，制约了矿业城市的可持续发展，其中矿区开采引起的区域地面沉降地质灾害尤为突出。针对这一问题，对煤矿开采沉陷进行实时动态监测，研究煤矿区的地表形变特征和时空演化规律，进而为地面沉降防控提供依据资料是一项重要工作。传统的水准测量、基岩标和分层标测量优点是方法简单、精度高，不足之处是观测数据时空分辨率低，测量成本高。全球导航卫星系统(Global Navigation Satellite System，GNSS)测量的水平方向形变精度较高，不足之处是对垂直向敏感度较差，且目前受限于空间分辨率严重不足。合成孔径雷达干涉测量(Synthetic Aperture Radar Interferometry，InSAR)技术具有全天时全天候成像、高时空分辨率、大范围连续监测、厘米甚至毫米的监测精度等特点，这些特点使得InSAR技术在地震同震、震间和震后形变监测、地质灾害隐患早期识别、矿区地表沉陷监测、城市基础设施监测等方面开展了诸多成功实践(敖萌等，2020；张景发等，2008；许强等，2019；廖明生等，2012；单新建等，2017；伍吉仓等，2020)，为煤矿区的开采沉陷监测提供了有力手段(张学东等，2012；徐小波等，2020)。

鹤岗是黑龙江四大煤城之一，为国家发展提供了重要的能源保障。随着鹤岗煤层开采不断延伸，地质环境问题日益突出。截至2016年，鹤岗采空区面积75.7km²，占煤田开发区域的46.5%，引起煤田区出现大面积非均匀地面沉降，威胁采空区居民的生命财产安全(龙精华，2017)。相关学者对鹤岗地区的煤矿开采沉陷进行了监测研究。金丽华等(2021)使用26景3m分辨率的COSMO-SkyMed数据，利用小基线集合成孔径雷达干涉测量(Small Baseline Subset Interferometric Synthetic Aperture Radar，SBAS-InSAR)算法提取2013年1月—2016年12月鹤岗煤矿区160余万个监测点的形变信息，并对沉降结果精度和时空分布进行了统计分析。武丽梅等(2020)利用2015年6月—2018年7月的34景Sentinel-1数据，提取了鹤岗市采煤沉陷区的地面沉降信息，监测到最大沉降速率达-150mm/a，沉降区范围呈逐年减小趋势。赵喜江等(2019)利用2017年3月—2018年10月的26景Sentinel-1B数据对鹤岗矿区进行地面沉降时序监测，发现鹤岗矿区地面沉降较为严重，最大累计沉降量达到了-190mm，4个明显的沉降漏斗均位于主采区上方。前人研究均对鹤岗地区2019年前的采煤沉陷情况进行了时序监测，使用的SAR数据较少(30景左右)，导致干涉相对时间基线较长，且仅分析了地面沉降在空间上的分布情况与变化趋势。本文利用覆盖鹤岗全部煤矿区的2019年8月31日—2021年8月20日59景Sentinel-1BSAR数据，开展鹤岗煤矿区地面沉降时间序列InSAR监测研究，发现在研究期内，鹤岗煤矿区的形变范围与前人结果基本相似，前人研究给定的最大形变速率位于南山煤矿，而本文提取最大形变速率位于新陆煤矿，形变速率也进一步增大。我们发现在兴安矿和峻德矿局部开采面形变出现了正向抬升，呈现SN向条带状分布。本文使用了更多SAR数据进行小基线集解算，提取的形变速率更加连续完整，未出现明显的空洞现象，在此基础上，进一步分析了鹤岗煤矿区最新的地面沉降特征，尤其是鹤岗南部煤矿区，如富力煤矿、兴安煤矿、峻德煤矿的时空演化规律。

1 研究区概况及数据准备 1.1 鹤岗煤矿区概况

鹤岗市位于黑龙江省东北部，东临松花江，西邻伊春市，南接佳木斯市，北部与俄罗斯隔江相望，地处“两江一岭”(黑龙江、松花江、小兴安岭)围成的“金三角”区域，地势西北高东南低，是三江平原向小兴安岭山地过渡的上升地段。鹤岗市辖8个县级行政区，包括6个市辖区(向阳区、工农区、南山区、兴安区、东山区、兴山区)和2个县(萝北县、绥滨县)。鹤岗市煤矿资源丰富，煤矿区主要分布于市区东部，从北到南依次由兴山矿、益新矿、振兴矿、鸟山矿、南山矿、新陆矿、富力矿、兴安矿和峻德矿九大煤矿组成。为分析鹤岗煤矿区地面沉降在空间上的分布情况及时序演化特征，利用覆盖鹤岗煤矿区的Sentinel-1数据(图 1)开展鹤岗南部煤矿区地面沉降时序监测研究。

1.2 数据准备

选取2019年8月31日—2021年8月20日59景134轨道的Sentinel-1B数据进行雷达干涉处理，分析鹤岗煤矿区近年来的地面沉降分布情况及时空演化特征。Sentinel-1B的空间分辨率为5m×20m，C波段降轨VV极化，主要参数如表 1所示(赵超英等，2019)。去除地形相位所用的外部DEM数据为美国航空航天局提供的SRTM3 DEM数据，其余数据由N47E129和N47E130两幅影像拼接而成，空间分辨率30m。

2 SBAS InSAR监测方法

小基线集合成孔径雷达干涉测量技术由Berardino等(2002)提出。SBAS方法将SAR影像序列按照一定规则组合成若干子集，同一子集内干涉像对的时间和空间基线较小，而各子集之间的时间基线和空间基线均较大。利用最小二乘法获取每个小基线集的地表形变时间序列，进而使用奇异值分解方法对多个子集进行联合求解，最终获得地表形变速率结果。SBAS-InSAR基本原理是(杨成生等，2014；何秀凤，2012；廖明生等，2014)：假设有按照时间序列t_i (i=1，2，3，…，N)顺序获取的N+1幅覆盖同一地区的单视复数SAR影像，从中选出一幅影像作为主影像，将其他影像进行配准，设定合适的时间基线阈值和空间基线阈值，将时间与空间基线都小于该阈值的影像组成若干集合，其中每幅影像至少可与其他N幅影像中的一幅形成干涉像对。对每个集合内的影像对分别进行差分干涉处理，得到M幅差分干涉图，可表示为

$ \frac{N+1}{2} \leqslant M \leqslant \frac{N(N+1)}{2} $

(1)

以t₀时刻作为初始时刻，对于任意差分干涉图j(j=1，2，3，…，M)，任意时刻t_i相对于t₀时刻的差分相位为φ(t_i)，在去除平地相位和地形相位后，SAR坐标系下任一像元(x，r)对应的解缠相位可表示为(假设t_B>t_A)

$ \delta \varphi_j(x, r)=\varphi\left(t_B, x, r\right)-\varphi\left(t_A, x, r\right) \approx \frac{4 \pi}{\lambda}\left[d\left(t_B, x, r\right)-d\left(t_A, x, r\right)\right] $

(2)

式中，λ为雷达信号波长；d(t_B，x，r)和d(t_A，x，r)为在t_B和t_A时刻相对于t₀时刻的雷达视线向地表形变，即d(t₀，x，r)≡0；用d(t_i，x，r)(i=1，2，3，…，N)来表示需要提取的形变序列，对应的相位为φ(t_i，x，r)，则有

$ \varphi\left(t_i, x, r\right) \approx \frac{4 \pi}{\lambda} d\left(t_i, x, r\right) $

(3)

假设IE=[IE₁，…，IE_M]，IS=[IS₁，…，IS_M]分别为干涉处理时按时间序列排序主影像序列和辅影像序列，且满足IE_j>IS_j(j=1，2，3，…，M)，则所有的差分干涉图相位可以组成观测方程，即

$ \delta \varphi_j=\varphi\left(t_{E_j}\right)-\varphi\left(t_{S_j}\right), j=1, 2, 3, \cdots, M $

(4)

式(4)定义了一个由M个等式构成的方程组，含有N个未知数，可将差分干涉相位的方程组改写为矩阵形式表示，即

$ \boldsymbol{A}{\varphi}=\delta \varphi $

(5)

其中，A为M×N的系数矩阵。若所有数据同属一个小基线子集中，有M≥N，矩阵 A的秩为N，式(5)的解在最小二乘约束下可用矩阵形式表示，即

$ \hat{\varphi}=\left(A^T A\right)^{-1} A^T \delta \varphi $

(6)

当矩阵 A的秩小于N时，最小二乘解不唯一，此时需要对矩阵 A进行奇异值分解，求解最小范数意义上的最小二乘解，将各时间段上的形变速率在时间域上积分，获得整个监测时间段的形变时序(闫怡秋等，2021)。SBAS-InSAR处理技术流程如图 2所示。

3 鹤岗煤矿区SBAS-InSAR数据处理 3.1 SAR影像配准

在进行SBAS-InSAR数据处理过程中，首先找到鹤岗煤矿区所在的Sentinel-1B SAR图像中的burst单元，将覆盖研究区的burst单元生成TOPS影像的单视复数产品(Single Look Complex，SLC)。Sentinel-1卫星轨道控制在100m的空间管道内，保证了足够小的空间基线，故在本文中不考虑SAR影像的空间基线的影响，选择SAR成像时间段内时间相对居中的2020年8月13日观测数据作为参考主影像，将其余58景辅影像和外部DEM数据全部配准至参考主影像上。利用GAMMA软件中的S1_coreg_TOPS脚本进行TOPS SLC数据配准，根据精密轨道信息和地形信息计算配准查找表，再利用强度匹配和频谱差异方法精化配准查找表。迭代强度匹配直到方位向校正小于0.01个像素，然后迭代频谱差异方法直至方位向偏移值小于0.0005个像素(屈春燕等，2014；吴文豪等，2017；祝杰等，2022)。

3.2 干涉相对组合选取

鉴于鹤岗煤矿区为采空区，地表形变量较大导致图像整体相干性不高，因此设定比较严格的时空基线阈值(空间基线150m、时间基线36天)进行短基线干涉影像对组合，生成的123个干涉像对连接图如图 3所示。从干涉连接图可以看出123个干涉像对连续未中断，分布比较均匀，空间基线最小的干涉对为2019年8月31日—2019年9月24日的1.2m，最大干涉对为2020年1月10日—2020年2月15日的146.9m，时间基线最小12天，最大36天。

3.3 差分干涉处理

对123个干涉相对进行雷达干涉处理，距离向方位向视数比设置为5:1，保证数据解算中输出的相位图保持高分辨率。利用欧洲航天局提供的POD精密定轨星历数据进行SLC数据轨道校正，去除轨道误差。利用SRTM DEM去除地形相位，得到123个差分干涉图。采用Goldstein自适应滤波方法对差分干涉图进行滤波，减弱大气和各类噪声误差。采用最小费用流算法进行相位解缠，相位解缠中相干性阈值设置为0.4，对于相干性阈值以下的区域不进行解缠计算，采用高、低通滤波方法进一步消除大气延迟相位影响，得到差分解缠相位。

3.4 时序形变计算

检查差分干涉产出的123个相位解缠图的质量，剔除相位解缠错误和质量较差的干涉相对。然后基于解缠的相位建立观测方程，采用奇异值分解方法求解形变速率的最小二乘范数解，对求取的形变速率在时域上进行积分，得到形变时间序列结果。

4 地面沉降InSAR时序监测结果及时空演化特征分析 4.1 InSAR时间序列结果及空间演化特征分析

运用SBAS-InSAR时间序列技术提取鹤岗煤矿区2019年8月31日—2021年8月20日的地表形变信息(杨梦诗等，2017)。鹤岗工农区欣虹小学周边环境较稳定，在所有干涉相对中相干性高，故将其设置为相位解缠起点。59景Sentinel-1 SAR数据提取的地表沿雷达视线方向(Line of Sight，LOS)的平均形变速率如图 4所示，由图可以明显看出SBAS-InSAR提取的高相干点位空间分布密度高，典型形变区域的高相干点位呈现清晰、连续的分布状态，未发生明显空洞。

从鹤岗煤矿区2019年8月31日—2021年8月20日近两年的地表平均形变速率可以看出，InSAR技术提取的近垂直向形变呈现明显的不均匀、不规则、漏斗式运动特征(Chen et al，2016)，最大形变位于新陆煤矿开采区(图 4中最大沉降漏斗处)，形变速率约-370mm/a。地表抬升的形变速率最大达75mm/a，位于兴安煤矿区(图 4中形变点4与5处)。InSAR提取的19个地面沉降漏斗呈现离散分布、形态不一、边界清晰的特征，主要分布在鹤岗兴安区和南山区。典型的地面沉降漏斗区域与鹤岗九大煤矿开采区相对应，显示出煤矿资源开采加速了地面下沉运动。兴安矿和峻德矿局部开采面形变出现正向抬升，并呈现SN向条带状，显示出煤矿区开采面调整对缓解开采沉降起到了显著效果。图 5展示的是鹤岗煤矿区2019年8月31日—2021年8月20日地表累计形变量的分布变化(Yang et al，2018)。可以看出，在2019年8月31日—9月12日的12天内地表形变量几乎为零；在2019年8月31日—2020年2月3日的156天时间内，累计形变量分布图上开始出现19个离散且边界比较清晰的形变漏斗区，这些区域对应鹤岗煤矿开采区，漏斗的形变量在-100~-200mm，其中新陆矿地面沉降的面积最大，形变量达-200mm；2019年8月31日—2020年6月26日的300天时间内，鹤岗煤矿区19个形变漏斗的中心位置相对2019年8月31日—2020年2月3日期间基本不变，范围正向外扩散，漏斗面积也随之逐渐增大，漏斗的形变量在-200~-400mm；在2019年8月31日—2020年11月29日的456天时间内，19个形变漏斗区域相对于2019年8月31日—2020年6月26日未进一步向外扩散。在局部区域，新陆矿累计形变量进一步增大，2个大形变漏斗(≤-600mm)雏形显现，漏斗中心形变量达到-600mm的量级。而鹤岗南部的兴安矿和峻德矿局部开采面形变出现正向抬升，抬升形变量约50mm；在2019年8月31日—2021年4月22日和2019年8月31日—2021年8月20日期间，鹤岗煤矿区19个形变漏斗区域相对于2019年8月31日—2020年6月26日同样没有明显的向外扩散。新陆矿2个大形变漏斗进一步向外围扩张，边界逐渐清晰，中心形变量分别达到-800mm和-1000mm的量级。兴安矿和峻德矿局部开采面形变进一步正向抬升，呈现的形态似SN向条带(SN向约4.5km，EW向约0.5km)，南启G201国道与骏德街道办事处交汇区，北至G201国道与兴安路交汇区，最大形变量约150mm，此处出现明显地面抬升，可能与兴安矿和峻德矿开采面调整导致地下水回流有关。综上，鹤岗煤矿区2019年8月31日—2021年8月20日近两年的地面沉降呈现明显的不均匀特征，提取的19个地面沉降漏斗均位于矿井开采区域。在研究期内，自相对于第一景影像的沉降漏斗形成后，后续矿井开采面未在水平方向上延伸。随着煤层开采面进一步延深，漏斗区域的形变量进一步增大，至2021年8月，新陆矿局部区域最大累计形变量增加至-1000mm量级，显著大于鹤岗其他8个煤矿的最大累计形变量，因此有必要对新陆矿的地面沉降趋势持续关注。兴安矿和峻德矿局部开采面出现了约150mm的正向抬升，体现出矿区较好的管控效果。

4.2 鹤岗南部典型煤矿区形变时空特征分析

下面通过沉降漏斗的空间形态特征和特征点的形变时间序列，分别对鹤岗煤田南部3个煤矿区(富力煤矿、兴安煤矿和峻德煤矿)各自的形变时空特征进行分析(许强等，2021；Xu et al，2019；Liao et al，2013)。

(1) 富力煤矿位于图 4中红色标识区。从平均形变速率结果(图 4)可以看出，2019年8—2021年8月，富力煤矿出现了1个明显的NW向长度约2km的大形变漏斗区域，面积约1.69km²，占富力矿区面积的一半，且形变区域有进一步向东扩展的趋势。从1号点位的形变量时间序列(图 6(a))可以看出，沉降漏斗的最大形变量达525mm，形变序列在2019年8月—2021年3月基本呈线性发展趋势，在2021年3月—2021年5月形变速率有所减缓，速率趋近于零，但在2021年5月后形变速率突然快速恢复至2019年8月—2021年3月的水平，且速率有进一步增大的趋势。地面沉降时间演化经历了快速—短期平稳—恢复快速3个沉降变形过程。鉴于富力煤矿沉降漏斗的空间范围正在逐渐向外扩展，因此推测后续仍会有较大的沉降发生，形变趋势达到稳定可能需要更长的时间。有必要对富力煤矿的地表形变趋势重点关注。

(2) 兴安煤矿位于图 4中蓝色标识区。2019年8月—2021年8月兴安煤矿出现了3个明显的形变漏斗，呈离散分布，漏斗面积约2.31km²。其中北部2号点位所在的“月牙形”形变区域范围最大，面积约1.38km²，最大形变量达455mm；中部的3号点位形变区涉及的范围较小，面积约0.57km²，最大形变量超375mm。北部2号点位(图 6(b))和中部3号点位(图 6(c))的形变量时间序列特征非常相似，具体表现为：在2019年8月—2020年7月的形变量基本呈线性发展趋势，形变速率相对较大，但在2020年7月形变量时间序列出现明显拐点，后续在2020年7月—2021年8月形变速率持续减小，此期间的形变量仅150mm，2021年1月后形变量持续减小，至接近于零。2号点位和3号点位的地面沉降时间演化经历了快速—趋于平稳2个沉降变形过程。南部4号点位的形变范围在兴安煤矿中最小，面积仅约0.36km²，2019年8月—2021年8月地表形变速率总体保持持续快速状态(图 6(d))。综上，兴安煤矿北部2号点位和中部3号点位所在开采面的地面沉降得到了一定程度的缓解，南部4号点位开采面的沉降趋势仍在持续增强，对于4号点位形变区域也需要进一步关注。

(3) 峻德煤矿位于图 4中黑色标识区。2019年8月—2021年8月，峻德煤矿出现了5个明显的沉降漏斗，呈近NS向离散分布，沉降漏斗面积由北向南分别约0.32km²、0.39km²、0.56km²、0.48km²和0.47km²，总面积约2.22km²。5个沉降漏斗的最大形变量均超过了300mm，其中最大形变量为7号点位的415mm。综合5个点位的形变量时间序列特征发现，在2019年8月—2021年2月，5个点位的形变量基本均呈现线性发展的趋势，形变量比较稳定，期间部分点位的形变速率略有波动但变化量较小，不影响形变量时序曲线的整体趋势。峻德煤矿5个点位的形变量(图 6(e)~6(i))在2021年2月—2021年8月均呈现明显的减小，且逐步趋近于零，其中6号形变点甚至出现地表近垂直向反向向上运动的特征。综上，峻德煤矿地面沉降时间演化经历了快速—减缓—平稳3个变形过程，通过峻德煤矿5个点位的地表形变特征分析，初步认为峻德煤矿现阶段地面沉降在一定程度上得到有效控制，推测后续不会有较大的沉降发生。

5 结语

煤矿资源持续开采造成鹤岗煤矿区出现明显的地面沉降。利用覆盖鹤岗煤矿区的59景Sentinel-1B SAR数据，基于小基线集时间序列InSAR技术，成功提取了鹤岗煤矿区2019年8月—2021年8月的地表形变信息，详细分析了鹤岗煤矿区地面沉降的时空演变规律。在地面沉降分布基础上，进一步详细分析了鹤岗南部富力煤矿、兴安煤矿和峻德煤矿的开采沉陷形势，主要结论如下：

(1) 鹤岗煤矿区在2019年8月—2021年8月出现比较严重的地面沉降现象，地表近垂直向呈现出不均匀、漏斗式的运动形态特征，最大垂直平均形变速率达370mm/a。提取的19个沉降漏斗边界清晰，集中在煤矿矿井开采区，显示了煤矿开采加速矿区开采区地面下沉运动。兴安矿和峻德矿局部开采面形变出现了正向抬升，呈现SN向条带状分布。在非煤矿开采区域基本无地面沉降现象。

(2) 2019年8月至—2021年8月，自相对于第一景影像的沉降漏斗形成后，随着时间推移漏斗的空间形态基本保持不变，但不同矿区沉降漏斗的累计形变量在空间范围上表现略有差异，新陆煤矿矿井区的累计形变量在2020年11月后显著大于鹤岗其他8个煤矿的累计形变量，建议对新陆矿开采沉陷予以管控，预防发生更大的地质灾害。

(3) 2021年5月后，富力煤矿和兴安煤矿4号形变点所在工作面的地表形变趋势仍处于快速下沉状态，预计未来沉降值还会进一步增大，形变趋势达到稳定需要更长的时间，有必要对达到稳定前的地表形变趋势进行持续动态跟踪。