基于SBAS-InSAR 的采煤区地表形变监测与分析

宋旭斌，路 鑫

（潞安化工集团有限公司 古城煤矿，山西 长治 046000）

0 引 言

煤炭作为我国的主体能源，长期担负着国家能源安全和经济高速发展的重任。但大规模、高强度的地下开采往往会破坏地层的原始应力平衡，而且随着采空区范围扩大，直接顶在自身重力和上覆岩层荷载作用下破断，最终引起地表的移动变形[1]。因开采导致的地表沉陷易造成建筑物墙体开裂、交通道路损坏、农田积水等地质灾害问题，严重制约了矿区城市的可持续发展。因此，对采空区地表进行实时精密形变监测，获取动态沉降信息势在必行。

常规开采监测手段主要有精密水准测量和全球导航系统（Global Navigation Satellite System,GNSS），但在实际工程中往往面临观测点易丢失、观测精度易受干扰、成本投入巨大、时空分辨率低等问题，无法准确获取矿区整体的形变特征。合成孔径雷达干涉测量(Interferometry Synthetic Aperture Radar,InSAR)是一种新型的对地观测技术，它仅对同一区域的多景SAR 影像进行复共轭相乘即可提取大范围、毫米级的地面形变信息[2]，具有全天时、高精度、自动化等优势。而小基线集（Small Baseline Subsets InSAR,SBAS-InSAR）技术在差分InSAR（Differential InSAR,D-InSAR）基础上提出，有效克服了后者易受大气效应、时空失相干影响的缺陷[3]，它通过对SAR 影像上的高相干点进行相位模型建立和相位变化分析，即使在永久散射体PS（Persistent Scatterer）密度较低的矿区也能获取更为精确的长时序形变结果。

许多学者都曾将SBAS 技术运用于矿区的形变监测，并取得了一系列丰硕的成果。李达等[4]利用SBAS 技术对13 景TerraSAR-X 进行处理，获取了榆林某矿的地表沉降信息，发现地表下沉值与时间呈线性关系；栾元重等[5]借助19 景Sentinel-1A 影像和SBAS 技术对沉陷漏斗区域内的工作面进行时序监测，并利用水准数据验证了SBAS 技术的可靠性；柴华彬等[6]提出利用反距离加权插值方法（Inverse Distance Weighted,IDW）将实测水准数据与SBAS 解译结果相融合，以解决部分区域因失相干导致InSAR 形变值缺失的问题；张香凝等[7]以蒲河矿区为例，结合SBAS 技术和数值模拟共同探究了煤矿区地面沉降的时空变形规律和机制。

山西长治地区煤炭资源丰富，但因地下开采引起的地质灾害频发，亟需开展相关的地表沉陷监测工作。本文在前人研究的基础上，借助SBAS-InSAR 技术对29 景Sentinel-1A 影像进行解译处理，以期获取古城煤矿S1306 工作面在2021年10 月12 日至2022 年10 月7 日期间的时序沉降结果和开采影响范围。研究成果对矿区开采规划和地表塌陷的预警防治均具有重要的现实意义。

1 研究区与数据源

1.1 研究区概况

本文研究区位于长治市古城煤矿东南方向一采区，整体地势呈西高东低，地貌类型以河谷平原为主。研究区分布着大面积农田、村庄房屋以及电力线等重要的基础设施，且覆盖有较厚的第四系黄土沉积；地下山西组3 号煤层目前有多个工作面正在回采，其中S1306 工作面走向长约1 982 m，倾向宽约363 m，平均采厚6.4 m，煤层倾角约为0～6°，底板标高-352—-397 m。开采方法以倾斜长壁式综采法为主，全部垮落法管理顶板。目前，该采区周边已出现大量明显的地裂缝及路面损坏情况，亟需对该区域进行地表沉陷监测工作。研究区的地理位置如图1 所示。

图1 研究区地理位置Fig.1 Geographical location of the study area

1.2 数据源

本文从欧空局（European Space Agency,ESA）哥白尼开放数据访问中心（https://scihub.copernicus.eu/dhus/#/home）收集了29 景Sentinel-1A 升轨单视复数影像，时间跨度为2021 年10 月12 日至2022 年10 月7 日，工作模式均为干涉宽幅（Interferometric Wide-swath，IW）模式，该模式基于递进式地形扫描(Terrain Observation with Progressive Scans SAR，TOPSAR)技术，能够将3 个子条带合并成1 景高质量的SAR 雷达数据[11]，其地面分辨率可达5 m×20 m。实验影像的其他参数见表1。另外，在数据处理过程中，引入了POD精密轨道文件用来去除卫星轨道的系统性误差，并利用30 m 分辨率的STRM1 DEM 数据消除了干涉图中地形相位的影响。

表1 Sentinel-1A 基本参数Table 1 Basic parameters of Sentinel-1A

2 技术原理与数据处理流程

2.1 SBAS-InSAR 原理

小基线集技术最早由国外学者Berardino 等人在2002 年提出，其基本原理如下[12-15]：假设选取（t0,t1,t2,…,tN）时段内覆盖同一区域的N+1 幅SAR影像，基于设定的时空基线阈值，组合成M 幅多视差分干涉对，其中M 的取值范围满足：

若在tA和tB（tB>tA）时刻对2 幅SAR 影像进行干涉处理，并生成了第j 幅已解缠的差分干涉图，在去除地形相位、大气延迟相位和噪声误差之后，任意像元位置（x,r）（x 和r 分别为方位向和距离向坐标）处的差分干涉相位可表示为：

式中：λ 为雷达波长；φ(tB,x,r)和φ(tA,x,r)分别表示2 幅影像上的相位值；d(tB,x,r)和d(tA,x,r)分别为tA和tB时刻相对于初始时刻t0在雷达视线向(Line of Sight,LOS)上的形变，一般假设t0时刻的形变相位值为0；vj为tA至tB时段内的平均形变速率。

将用于干涉处理的M 对影像数据按照时序进行排列，则所有干涉图相位可表示为如下矩阵形式：

式中：B 为1 个M×N 的系数矩阵；V 为速度矢量。

由于矩阵B 的秩小于N，根据最小二乘法无法得到唯一的结果，因此需要采用奇异值分解法（Singular Value Decomposition,SVD）联合多个基线集，来解决法方程的秩亏问题，并得到速度V 的广义逆解。之后将各个时间段内的形变速率在时间域上进行积分，即可获取整个监测时段内地表在LOS 向上的累积时序形变结果。

2.2 SBAS-InSAR 数据处理流程

SBAS-InSAR 数据处理关键步骤和参数设置如下：①选取2021 年11 月29 日的影像作为主影像，其余从影像与之配准，并设置空间基线百分比阈值为45%，时间基线阈值为60 d，多视比为4∶1（距离向∶方位向），最后共生成104 个干涉对；②利用最小费用流(Minimum Cost Flow,MCF)方法进行相位解缠，Goldstein 滤波技术去除干涉图噪声，并逐个检查调整不理想的像对；③参考干涉处理步骤生成的相干图和解缠相位图，选取20～30 个稳定且分布均匀的地面控制点（Ground Control Point,GCP）完成轨道精炼，去除残余轨道误差相位；④选择线性模型（Linear Model）进行二次解缠，并估算平均位移速率和残余地形相位；⑤根据大气延迟相位、形变相位和噪声相位表现出的不同的时空特性，利用时间域的高通滤波（365 d）和空间域的低通滤波（1 200 m）分离并剔除大气延迟误差，得到覆盖整个观测时间的地表形变时间序列；⑥将最终解译结果从斜距坐标系转换至WGS-84 坐标系，即完成地理编码。本文的数据处理流程如图2所示。

图2 SBAS 数据处理流程Fig.2 Data dealing process of SBAS

3 解译结果与分析

3.1 空间演化分析

利用SBAS-InSAR 技术对覆盖S1306 工作面的29 景影像进行解译处理，获取了监测时段内研究区LOS 向的时序累积形变结果。之后忽略水平位移的影响，利用公式dv=dLOS/cosθ 将解译结果从视线向投影至垂直向，其中θ 为卫星入射角。另外，解译结果中的正值表示地面抬升，负值表示地面下沉。

S1306 工作面自2022 年1 月开始回采，目前仅开采到距离始采线约850 m 的位置。但图3 显示，在监测时段内，该工作面周边地表已经发生了明显沉陷。

图3 S1306 工作面时序形变Fig.3 Time series deformation of No.S1306 Face

在形变监测的初期阶段，即2021 年10 月12日至2021 年12 月23 日（图3a～图3c），研究区地表总体较为稳定，零星区域发生小幅沉降，量值约-13 mm，产生该现象的原因可归结为工作面所处的农田区域相干性较差，导致SBAS 在反演形变的过程中发生解缠误差，影响了结果的准确性；从2022 年1 月16 日起（图3d），工作面的始采线位置开始出现明显下沉，最大沉降量约-32.6 mm，符合工作面开始回采的时间；随着工作面持续推进，地表沉降范围逐渐朝西南方向（即工作面开采方向）扩张；到了2022 年5 月16 日（图3h），始采线附近地表已经形成了一个较为明显的沉降盆地，盆地中心可探测到的最大下沉值约-178.6 mm；随着采动影响逐渐传递至工作面中部位置（图3i～图3k），始采线区域的下沉速度逐渐趋于平缓，但仍处在下沉状态；直至2022 年10 月7 日（图3m），此时S1306 工作面的开采影响边界稍超过实际开采进度位置，但总体与实际情况相符，说明SBAS 技术在矿区形变监测中具有一定的可靠性和准确性。

3.2 时间演化分析

为了进一步探究研究区地表在采动影响下的时序形变过程，在S1306 工作面的走向和倾向上共选取8 个特征点（编号点A～点H），并以2021 年10 月12 日影像为参考基准，生成各个点的时序累积沉降曲线。具体点位分布如图4 所示。

图4 特征点位及其时序累积沉降值Fig.4 The location of feature points and their time-series cumulative subsidence values

由图4 可知，在2021 年10 月12 日至2022 年1 月16 日期间，走向线上4 个特征点没有发生明显形变，由于点A 和点B 分别位于始采线和沉陷盆地中心，易受开采扰动，这一时段内两点均出现小幅度下沉，最大沉降值约-10 mm；而点C 和点D 分布在工作面中部和停采线的位置，点位表现出反复微弱抬升的迹象，可将这种情况其视为SAR图像噪声带来的解算误差。

2022 年1 月16 日之后，点A 和点B 开始持续加速下沉，且分别于2022 年3 月17 日和2022 年5 月16 日，其形变速率才逐渐放缓，但形变量仍在持续增加，这是由于采动作用对始采线区域地表的影响即将达到“临界点”，该区域的形变经历了“开始下沉—剧烈下沉—缓慢停止下沉”的过程，一定程度上符合常规开采沉陷规律。点C 和点D在2022 年8 月8 日之前，下沉曲线吻合度较高，但之后点C 的沉降量陡增至-87.6 mm，说明工作面中部区域此时开始受采动的强烈影响，与实际工程进度一致。同时，发生迅速沉降的原因还可能是由于点C 处在土质松软的农田，在夏季强降水的渗透下，土层的抗剪强度会大大降低，从而引起地表的湿化崩解和垮落[12]。

倾向线上的4 个特征点均匀分布在工作面的两侧，在2022 年1 月16 日之前，由于点E 和点H远离开采中心，点F 和点G 位于沉陷盆地边缘，因此后者呈现小幅持续下沉的态势，前者沉降与起伏交替，受采动影响较小；之后点F 和点G 开始加速下沉，但点F 的沉降幅度与速率要大于点G。同样，点E 和点H 于2022 年4 月10 日也开始经历一个缓慢下沉的过程，但总体上点E 的沉降量级要高于点H。

因此结合图3 的时序形变结果，可以发现沿开采方向右侧的地表形变值要大于左侧区域，即始采线附近的沉陷盆地实际上有朝西北方向发育的趋势，后期需要对东李高村附近地表进行重点监测。

4 结 论

（1）在采动作用下，S1306 工作面周边地表自2022 年1 月16 日开始出现下沉，之后在始采线附近迅速产生了一个沉陷盆地，且沉陷盆地有朝西北方向发育的趋势。

（2）目前采动影响仅波及至工作面的中部区域，与实际开采到达位置基本吻合，说明SBAS-InSAR 技术在矿区形变监测中具有良好的适用性。

（3）除了利用时序InSAR 技术获取不同时段的地表沉陷情况，仍可以考虑利用UAV 激光雷达扫描技术或建立地面移动观测站来对局部重点区域进行详细监测，真正发挥“空-天-地”协同监测技术在开采沉陷领域的重要作用。