基于SBAS-InSAR的南水北调中线工程禹州采空区沉降监测与分析

Settlement Monitoring and Analysis of Yuzhou Goaf in the Middle Route of South-to-North Water Diversion Project Based on SBAS-InSAR

摘要: 南水北调中线一期工程总干渠禹州段采空区稳定性是该工程面临的重要技术挑战。本研究利用2016年2月至2017年4月获取的19景C波段Sentinel-1A升轨雷达影像数据，采用短基线集干涉测量(SBAS-InSAR)技术对新峰煤矿、郭村煤矿、工贸煤矿和福利煤矿等4个采空区的地表形变进行监测分析。监测结果显示：郭村煤矿和工贸煤矿采空区地表沉降最为显著，最大沉降速率分别达到−22 mm/a和−20 mm/a，其中郭村煤矿采空区东侧累计沉降量最大为−22 mm；4个采空区地表变形参数分析表明，最大倾斜变形(−0.782 mm/m)和最大曲率变形(−0.041 mm/m²)均出现在郭村煤矿采空区。通过与二等水准测量数据对比验证，SBAS-InSAR监测结果与水准数据的均方根误差(RMSE)小于2 mm，相关系数高达0.96，两者在形变趋势和量级上均表现出良好的一致性，充分验证了SBAS-InSAR技术在南水北调渠段采空区安全监测中的可靠性和适用性。

Abstract: The stability of goafs in the Yuzhou section of the main canal of the first phase of the South-to-North Water Diversion Middle Route Project is an important technical challenge faced by the project. In this study, 19 scenes of C-band Sentinel-1A ascending orbit radar image data obtained from February 2016 to April 2017 were used, and the short-baseline set interferometry (SBAS-InSAR) technology was adopted to monitor and analyze the surface deformation of four goafs, including Xinfeng Coal Mine, Guocun Coal Mine, Gongmao Coal Mine, and Fuli Coal Mine. The monitoring results show that the surface subsidence of the goafs in Guocun Coal Mine and Gongmao Coal Mine is the most significant, with the maximum subsidence rates reaching −22 mm/a and −20 mm/a respectively, among which the maximum cumulative subsidence on the east side of the goaf in Guocun Coal Mine is −22 mm. The analysis of surface deformation parameters of the four goafs shows that the maximum inclined deformation (−0.782 mm/m) and the maximum curvature deformation (−0.041 mm/m²) both occur in the goaf of Guocun Coal Mine. By comparing and verifying with the second-order leveling measurement data, the root mean square error (RMSE) between the SBAS-InSAR monitoring results and the leveling data is less than 2 mm, and the correlation coefficient is as high as 0.96. Both show good consistency in the deformation trend and magnitude, fully verifying the reliability and applicability of the SBAS-InSAR technology in the safety monitoring of goafs in the South-to-North Water Diversion canal section.

1. 引言

南水北调中线工程作为全球规模最大的跨流域调水工程，在缓解我国北方地区水资源短缺、优化国家水资源配置格局、支撑经济社会可持续发展等方面发挥着不可替代的战略作用[1]。然而，工程沿线部分渠段面临着复杂的地质环境挑战：地下水超采、矿区开采活动、煤矿采空区沉降以及区域地质构造活动等多重因素相互叠加，导致局部区域地质稳定性显著降低[2]。特别是在煤矿采空区等特殊地质段，潜在的地面沉降风险可能直接影响渠道结构的完整性和工程输水安全。为此，采用高精度形变监测技术对渠道穿越的煤矿采空区实施全方位、立体化监测及预测，不仅是保障工程长期稳定运行的技术需求，更是维护国家重大水利基础设施安全的战略举措。

南水北调工程的安全监测目前主要依赖于水准测量和GNSS技术，然而这些传统方法受限于地形复杂性、气象条件以及时间成本等因素，仅能获取离散点位沉降数据，难以全面反映工程沿线的连续形变特征及整体变形趋势[3]。为解决这一技术瓶颈，20世纪70年代发展起来的合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar, InSAR)技术凭借其全天时、全天候的观测能力，高精度的测量优势以及大范围覆盖的特点，有效弥补了传统监测手段的不足[4]，该技术已在矿区沉陷[5]-[7]、灾害预警[8]-[10]、地面沉降[11]-[12]和地震监测[13] [14]等各个领域得到广泛应用。然而常规差分干涉测量(D-InSAR)技术受时空失相干和大气延迟效应的影响，其监测精度存在明显局限[15]。为突破这一技术限制，Berardino [16]等人在2002年首次提出了短基线集干涉测量(Small Baseline Subset, SBAS)方法，该方法通过优选空间基线和时间基线较短的SAR影像对进行联合处理，显著降低了失相干和大气误差的干扰，实现了大范围、长时序、毫米级微小形变的连续监测。李达[17]等通过采用SBAS-InSAR技术对陕西一矿区进行了残余变形的监测，并用D-InSAR的结果验证了SBAS-InSAR监测结果的精度；余昊[18]等用SBAS-InSAR技术对徐州西部地区关闭矿井的地表变形规律进行了研究，发现矿井关闭的地区地面有抬升的趋势；张金芝[19]等通过SBAS-InSAR技术分析黄河三角洲地面沉降主要与石油开采等因素有关。张永光[20]等利用SBAS-InSAR技术对南水北调膨胀土区域进行精细化分析，结果表明，渠段的沉降与降雨量、土壤湿度等因素有高度相关性。

现有研究大多集中在矿区开采过程中产生的地表形变，对于矿井关闭后形成的采空区研究仍值得关注，其残余变形对上方设施建筑的影响仍是一个重要问题。本文以南水北调中线禹州渠段采空区为研究对象，利用SBAS-InSAR技术监测并分析了研究区域的时空演化规律，通过与二等水准测量结果验证其精度和可靠性，为渠道的安全运行、地下资源的合理开采提供技术依据。

2. 研究区概况及数据来源

2.1. 研究区概况

研究区位于河南省禹州市境内，所在的区域地形地貌以低山丘陵为主，地形平缓。主要包括原新峰矿务局二矿、梁北镇郭村煤矿、梁北镇工贸公司煤矿、梁北镇福利煤矿4个煤矿采空区。通过采空区渠段长3.11 km，采空区煤层多为1层，单层煤层厚度约为1米，埋深多为100~290 m之间。回采率为64.7%~78.3%。覆岩为软质岩夹硬质岩层，软硬岩比例为1:3~1:5，岩层走向近东西，倾向南，倾角12˚~19˚。采空区上覆岩土体的变形特征，总体来看是自下而上逐渐变形呈“漏斗状”沉降，“三带型”是最常见的上覆岩土体变形破坏形式。目前四个煤矿均已关停或废弃。

Figure 1. General situation map of the study area

图1. 研究区概况图

2.2. 研究数据

Sentinel-1A卫星由欧空局于2014年4月在法属圭亚那库鲁发射，轨道高度693 km，搭载的是C波段合成孔径雷达。它具有12 d的重访周期，并且有着多种极化方式和工作模式。本文选用2016年2月至2017年4月覆盖南水北调中线禹州段采空区的19景Sentinel-1A干涉宽幅(IW)模式VV极化的影像进行实验，影像空间分辨率为5 m × 20 m (方位向 × 距离向)。Sentinel-1A升轨影像数据基本参数信息如表1所示。由于IW模式的影像覆盖面积(250 km × 250 km)远大于本文研究区域，对合成孔径雷达(synthetic aperture radar, SAR)影像进行裁剪得到本文研究区域，研究区范围如图1所示。采用Sentinel-1A卫星的精密轨道数据(precise orbit ephemerides, POD)修正轨道误差，采用由美国太空总署(NASA)和国防部国家测绘局(NIMA)联合测量的分辨率为30 m的SRTM DEM数据以去除地形相位。

Table 1. Satellite image data

表1. 卫星影像数据

参数	数值
获取卫星	Sentinel-1A
轨道	升轨
分辨率/(m × m)	5 × 20
雷达波长/cm	5.6
极化方式	VV
重放周期/d	12
入射角/(˚)	41.5
影像时间	2016-02~2017-04
影像数量/景	19

3. 研究方法与数据处理

3.1. SBAS-InSAR基本原理

短基线干涉测量技术是Berardino [16]等于2002年提出的一种用于监测地面微小形变的新型时间序列分析方法。SBAS-InSAR技术能有效克服合成孔径雷达差分干涉测量技术在数据处理过程中的限制因素(时空失相干、大气效应等)，可连续监测地面微小形变，监测精度可以达到毫米级。基本原理如下：若有同一研究区域的N+1景SAR影像，则可能得到M个干涉对，M需满足：

$\frac{N+1}{2}\le M\le \frac{N\left(N+1\right)}{2}$ (1)

设初始时刻为$t_0$ ，任意时刻($t_i=1,\cdots ,N$ )相对于初始时刻的相位差是未知参数$\phi \left(t_i\right)$ ，干涉处理获取的差分干涉相位$\delta \phi \left(t_k\right)$ ($k=1,\cdots ,M$ )为观测量。若不考虑其他相位误差，则第$k$ 景差分干涉图的像元$\left(x,y\right)$ 地表形变相位$\delta {\phi }_{\text{def}}$ 可表示为：

$\delta {\phi }_{\text{def}}=\phi \left(t_{\text{B}},x,y\right)-\phi \left(t_{\text{A}},x,y\right)=\frac{4\pi }{\lambda }\left(r_{t_{\text{B}}}-r_{t_{\text{A}}}\right)=\frac{4\pi }{\lambda }\delta r$ (2)

式中，$t_{\text{A}}$ 、$t_{\text{B}}$ ($t_{\mathrm{A}}<t_{\text{B}}$ )为干涉图k($k=1,\cdots ,M$ )对应的SAR影像获取时刻，$\delta {\phi }_{\text{def}}$ 为$t_{\text{A}}$ 到$t_{\text{B}}$ 时刻之间视线向的形变相位$\lambda$ 为雷达波长，$r_{t_{\text{B}}}$ 和$r_{t_{\text{A}}}$ 分别和为$t_{\text{B}}$ 和$t_{\text{A}}$ 时刻干涉图的任意像元对应视线向的地表形变，$\delta r$ 为该时段的地表形变量。

将所有差分干涉图进行自由组合，$\phi$ 为N幅SAR影像上的干涉相位值组成的矩阵，$\delta \phi$ 为M幅差分干涉图上相位组成的矩阵。则有：

$\phi ={\left[\phi \left(t_1\right)\cdot \cdot \cdot \phi \left(t_N\right)\right]}^{{\rm T}}$ (3)

$\delta \phi ={\left[\delta \phi \left(t_1\right)\cdot \cdot \cdot \delta \phi \left(t_M\right)\right]}^{{\rm T}}$ (4)

将上式写成矩阵形式：

${\rm A}\phi =\delta \phi$ (5)

式中，A为M × N维矩阵。

当小基线集合数L = 1时，矩阵A的秩为N，采用最小二乘法即可求解出累积形变相位的估值$\widehat{\phi }$ ：

$\widehat{\phi }={\left({{\rm A}}^{{\rm T}}{\rm A}\right)}^{-1}{{\rm A}}^{{\rm T}}\delta \phi$ (6)

当基线组L > 1时，式(5)秩亏，秩亏数为$N-L+1$ ，为此对系数矩阵A进行奇异值分解，求出累积形变相位$\phi$ 最小范数意义上的最小二乘解。

由于试验区域为关闭矿井，总体地表形变量较小，在不考虑水平移动的情况下，可利用式(7)得到垂直方向的沉降量W：

$W=\frac{d\left(t_{\text{end}},x,r\right)}{\text{cos}\theta }$ (7)

式中，$t_{\text{end}}$ 为最终时刻，$x$ ，$r$ 分别为方位向与距离向坐标，$d\left(t_{\text{end}},x,r\right)$ 为相对于初始时刻的视线方向(Line of sight, LOS)的形变量，$\theta$ 为卫星入射角。

3.2. 数据处理

(1) 数据预处理：由于原始数据较大，为了提高数据处理速度及效率，对所有原始Sentinel-1A影像进行预处理，并裁剪得到本研究区域范围。

(2) 连接图生成：综合考虑空间基线、时间基线和多普勒质心频率基线等因素，选取(2016-10-02)获取的影像为超级主影像，其他影像均为辅影像，将时间基线阈值设置为180 d，空间基线阈值设置为临界基线的2%，避免由单一主影像引起的时空失相干误差。共生成100个小基线干涉像对，用于后续差分干涉处理，像对组合连接情况如图2时空基线图所示。

(3) 干涉处理：采用多级配准方式完成主辅影像的精确配准，通过复共轭相乘生成干涉图。利用精确轨道数据模拟并去除平地相位后，采用自适应滤波方法降噪并生成相干系数图。相位解缠过程使用最小费用流(Minimum Cost Flow, MCF)奇异值分解法(singular value decomposition, SVD)对所有的干涉对进行解缠，解缠相干系数阈值设为0.2。

(4) 轨道精炼与重去平：相位解缠后仍存在残余的恒定相位和相位坡道，在没有残余地形条纹、相位跃变及远离形变条纹的区域选择合适数量的GCP点，通过GCP点将上述残余相位去除。

(5) 形变反演：采用分步优化策略确保结果精度。首先基于线性模型开展核心反演计算，通过最小费用流方法进行二次相位解缠，精确估算形变速率和残余地形相位。在此基础上，进一步实施空间域低通和时间域高通滤波处理，有效分离并去除大气延迟相位干扰，最终获得高精度的地表形变反演结果。

(6) 地理编码：将SAR坐标系下的累积沉降量和形变速率转换到地理坐标系，得到地理坐标系下的累积沉降量和形变速率。

SBAS-InSAR方法的数据处理技术路线如图3所示。

Figure 2. Time-space baseline map

图2. 时空基线图

Figure 3. Technical route for SBAS-InSAR data processing

图3. SBAS-InSAR数据处理技术路线

4. 可靠性分析

为了验证SBAS-InSAR技术监测南水北调中线工程禹州段采空区形变结果的可靠性，选取图1所示的研究区内8个水准点进行对比验证。提取水准点位置的形变时间序列值，将监测点视线向的形变投影到垂直方向上，同时根据时间段将二等水准观测值归零化处理，为保证时间的一致性，通过插值拟合将水准监测数据转换到与其对应的InSAR结果进行对比，采用折线图进行分析。SBAS-InSAR监测结果和二等水准测量结果的时间序列位移量变化如图4所示，形变趋势吻合，充分证明了SBAS-InSAR技术的可靠性。

为了进一步定量分析SBAS-InSAR形变监测结果的可靠性，利用均方根误差(Root Mean Square Error, RMSE)评价SBAS-InSAR监测结果与水准测量结果的吻合程度。由表2可知，两者之间具有较高的一致性，互差为0~2 mm，二等水准观测值与SBAS-InSAR结果之间的RMSE值均小于2 mm。

$\text{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n\text{Δ}{h}_i^2}}$ (8)

式中：h为第i个水准数据与SBAS-InSAR监测数据结果沉降量之差；n为选取点的个数。

Figure 4. Comparison between SBAS-InSAR monitoring results and second-order leveling measurement results

图4. SBAS-InSAR监测结果与二等水准测量结果对比

Table 2. Comparison between second-order leveling results and SBAS-InSAR results

表2. 二等水准结果与SBAS-InSAR结果对比

点号	P1	P2	P3	P4	P5	P6	P7	P8
最大差值/mm	0.67	−1.39	4.16	1.61	1.22	−0.52	0.93	1.72
RMSE/mm	0.35	0.63	1.15	0.87	0.45	0.33	0.43	0.61

统计2016年2月至2017年4月SBAS-InSAR与8个水准点二等水准测量整体形变量的相关性，如图5所示，相关系数的计算公式见式(9)。

$R=\frac{{\displaystyle \sum _{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}}{{\displaystyle \sum _{i=1}^n{\left(x_i-\overline{x}\right)}^2}{\displaystyle \sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$ (9)

式中：x为SBAS-InSAR监测值；y为二等水准测量值。

Figure 5. Regression analysis of deformation values monitored by SBAS-InSAR and second-order leveling measurements

图5. SBAS-InSAR监测形变量与二等水准测量值回归分析

经计算，二者相关系数达0.96，二者差异的RMSE为1.147 mm，两者在形变趋势和量值上高度一致，表明了本文SBAS-InSAR结果的可靠性。

5. 监测结果与分析

5.1. 沉降速率分析

由图6可知：新峰煤矿采空区大部分地表呈抬升趋势，最大速率达25 mm/a，位于采空区东侧，郭村煤矿采空区大部分地表呈下沉趋势，最大速率达−22 mm/a，位于采空区东侧，工贸煤矿采空区最大沉降速率达−20 mm/a，距离渠道最近距离275 m，福利煤矿采空区整体沉降速率较小。

5.2. 时序形变分析

为进一步分析关闭矿井后各采空区地表时空形变规律，在获取沉降速率的基础上对其进行时间维度上的积分获取各采空区实现序列累计沉降量的变化，结果如图7所示。由图可知：在2016-02~2017-04期间，工贸煤矿采空区南部、郭村煤矿采空区渠道两侧均有小范围沉降，最大沉降量达−22 mm，位于郭村煤矿采空区最东侧，距离渠道较远；新峰煤矿采空区、福利煤矿采空区沉降不明显；四个煤矿采空区均未形成明显的沉降漏斗。通过对工程相关资料以及监测数据的分析，采空区充填体由于浆体的固结沉降、采空区巷道分布复杂、形状不规则、采煤方法和顶板管理方法落后等原因，充填后的采空区可能存在空洞，渠道在长期运行过程中，随时间推移，基岩内的地下水通过裂隙到达空洞，会在采空区及上覆岩体内形成一定程度的水压。渠道长期运行中，充填体及上覆岩体在水压作用下的蠕变同样会影响地表变形。实际监测数据显示，郭村采空区在渠道运行过程中仍有少量剩余沉降。

为进一步分析研究区沉降性质，在每个矿区沉降量最大的区域及其周边区域分别提取对应的观测线(如图1所示)，并利用式(7)将视线向的沉降转换为垂直向沉降。根据开采沉陷理论公式，即式(10)和式(11)，分别求出观测线上对应各点的倾斜与曲率值，其曲线如图8和图9所示。

$T_{\left(k-1\right)-k}=\frac{W_k-W_{k-1}}{L_{\left(k-1\right)-k}}$ (10)

$C_{\left(k-1\right)-k-\left(k+1\right)}=\frac{T_k-T_{k-1}}{\left(L_{k-\left(k+1\right)}+L_{\left(k-1\right)-k}\right)/2}$ (11)

式中：$T_{\left(k-1\right)-k}$ 为观测点$k-1$ 和$k$ 之间的倾斜，mm/m；$W_k$ 、$W_{k-1}$ 观测点$k-1$ 和$k$ 的沉降量，mm；$L_{\left(k-1\right)-k}$ 为观测点$k-1$ 和$k$ 之间的水平距离，m；$C_{\left(k-1\right)-k-\left(k+1\right)}$ 为点$k-1$ ，$k$ ，和$k+1$ 之间的曲率，mm/m²。

Figure 6. Annual average deformation rate in the LOS direction

图6. LOS向年平均形变速率

Figure 7. Temporal cumulative settlement amount in the study area

图7. 研究区时序累计沉降量

由图8和图9可知：新峰煤矿采空区观测线上最大倾斜和曲率出现在2016-10-26，分别为0.326 mm/m，0.019 mm/m²，郭村煤矿采空区观测线上最大倾斜和曲率在2017-03-07，分别为−0.782 mm/m，−0.041 mm/m²，工贸煤矿采空区观测线上最大倾斜和曲率出现在2016-05-23，2016-10-02分别为−0.381 mm/m，−0.015 mm/m²，福利煤矿采空区观测向上最大倾斜和曲率出现在2017-03-07，分别为−0.363 mm/m，0.017 mm/m²。以上形变值均在“三下采煤规范”地表允许(Ⅰ级破坏变形值范围(倾斜≤3 mm/m，曲率≤0.2 mm/m²)以内[21]，说明目前地表形变暂时不会对地表上方渠道造成影响。同时，可以看出新峰煤矿采空区观测线上的倾斜和曲率在2016-10-26达到最大值后逐渐变小并趋于稳定，郭村煤矿采空区、福利煤矿采空区观测线上的倾斜和曲率均在2017-03-07达到最大值，还有增大的趋势，预测未来可能对采空区周围的村庄产生一定影响。

Figure 8. Inclination distribution curves in each goaf area

图8. 各采空区内倾斜分布曲线

Figure 9. Curvature distribution curves in each goaf area

图9. 各采空区内曲率分布曲线

6. 结论

本文基于2016-02~2017-04期间的19张C波段Sentinel-A升轨的雷达影像数据，采用SBAS-InSAR技术方法对南水北调中线工程禹州段4个煤矿采空区进行监测，分别从研究区内地面沉降的年平均速率、时序累计沉降量、观测线上的倾斜和曲率等进行了分析。结果表明：

(1) 利用SBAS-InSAR技术获取研究区域地面沉降信息，通过与二等水准测量结果对比，8个水准点之间的RMSE值均小于2 mm，相关系数为0.96，表明短时间尺度下SBAS-InSAR技术在该区域的监测结果具有可靠性，为南水北调工程安全监测提供了可行的技术参考。但需指出，本研究验证点数量有限(8个)，且未覆盖所有形变特征区域，未来需增加验证点密度以进一步提升结论的普适性。

(2) 监测期内，郭村煤矿和工贸煤矿采空区表现出较显著沉降特征，最大沉降速率分别为−22 mm/a和−20 mm/a，最大累计沉降量−22 mm出现在郭村采空区东侧。需要说明的是，这些结论仅反映2016-2017年间的形变趋势，由于煤矿采空区沉降具有时间累积效应，更长期的监测数据对于准确评估沉降发展趋势至关重要。

(3) 基于观测点计算的形变参数显示，郭村采空区存在局部形变集中现象(最大倾斜−0.782 mm/m，曲率−0.041 mm/m²)。参照《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》(2017版) [21]，当前形变量值尚未超出地表允许变形阈值。但值得注意的是，本研究时间序列仅14个月，未能涵盖采矿活动的完整周期，后续需结合更长时序数据评估形变的累积效应，尤其需关注曲率变化对渠道结构的潜在影响。

致 谢

在此对欧洲空间局提供的卫星影像数据和美国太空总署、国防部国家测绘局提供的数字高程模型表示感谢。

基金项目

本文获得了华北水利水电大学第十六届研究生创新能力工程(NCWUYC-202416032)的支持。