基于TSS-DInSAR方法的注浆采区地表动态沉降分析

王剑，董祥林，杨可明，姚树一，石晓宇

(1.中国矿业大学(北京)地球科学与测绘工程学院，北京，100083；2.淮北矿业股份有限公司通防地测部，安徽淮北，235000)

地下采矿活动引起的地表大梯度快速沉降会导致上覆建筑物、管线以及道路等人工设施发生损坏[1]，因此，常采取相应科学有效的采矿工艺或开采措施以减弱地表下沉带来的危害，确保矿区开采的可持续性和安全性。同时，也必须提供精确可靠的地表沉降监测方法来获取矿区地表下沉盆地范围以及地表沉降量。差分干涉合成孔径雷达测量(differential interferometry synthetic aperture radar,DInSAR)作为新兴的空间对地观测技术，可以敏感地获取地表微小变化，其监测精度可达厘米级，甚至毫米级，并且具有全天候、全天时、大范围便捷监测的优势[2]。DInSAR 技术及由其发展而来的时序InSAR 技术已成为监测矿区地表沉降的有效热门手段之一[3-5]。目前，国内外学者利用InSAR 技术监测矿区地表沉降主要有3 种：PS(permanent scatterers，永久散射体)、SBAS(small baseline subsets，小基线集)等时序InSAR 技术[6-7]，POT(Pixel offset-tracking，像素偏移量追踪)技术[8-9]以及DInSAR 技术[10-11]。虽然时序InSAR 技术对形变梯度较小的地表形变具有很高的监测精度[12]，但得到的结果不连续，只能获取部分高相干点目标的形变信息[13]，并且在对矿区地表大梯度沉降进行监测的过程中存在结果明显偏低的缺陷[14-15]。POT技术虽然可以监测到地表较大沉降量[16]，但数据处理时间较长，且结果精度低，仅为像元分辨率的1/10～1/30[14,17]。而DInSAR 技术处理时间短，所需影像少，得到的地表形变结果连续，可监测的地表沉降最大梯度适中，得到的结果精度较高，更适用于监测矿区地表沉降[18]。此外，利用上述InSAR技术获取煤矿区地表沉降信息研究大多数是针对开采的煤层上覆岩层在自然垮落情况下的矿区地表沉降监测，对开采工作面注浆充填影响下的矿区地表动态沉降监测研究很少[6,8-11]；且大多采用的是商业L、X波段数据[3,17-18]，若要得到地表长期的动态沉降信息则所需成本较高。因此，本文作者利用哨兵1 号A 星(Sentinel-1A)的C 波段影像数据，以淮北矿业集团临涣煤矿Ⅱ2 采区1021 和Ⅱ1021 工作面采动影响范围为试验区，采用时序叠加DInSAR (time series stacking DInSAR, TSSDInSAR)方法并结合现场实测水准数据来监测分析工作面开采时在注浆充填影响下的矿区地表动态沉降，以探寻注浆充填的地表减沉效果以及注浆充填影响下地表动态沉降规律，同时为矿区提供大范围、低成本的地表长期动态沉降监控分析新手段。

1 原理与方法

TSS-DInSAR 方法是建立在二轨法DInSAR 基础上的。二轨法DInSAR 最早由MASSONNET等[19]提出，其主要思想为利用形变前后各一景SAR影像以及参考DEM数据，通过差分处理生成包含不同组分相位信息的差分干涉图，然后将除地表真实形变相位之外的其他相位信息剔除，即：

式中：λ为合成孔径雷达波长。由于垂直沉降是矿区地表变形的主要贡献分量[20]，故忽略水平移动分量对雷达视线方向形变的贡献，由三角函数关系得地表垂直方向的下沉量W为[6,20-21]

式中：θ为雷达成像几何入射角。

对按时序获取的同一位置的n 景影像A=[a1，a2，…，an](各影像对应的获取时间分别为t1，t2，…，tn，且满足t1＜t2＜…＜tn)，每相邻两景重复进行二轨法DInSAR处理，得到每相邻2个影像时段内的地表沉降，即a1与a2得到w1～2，a2与a3得到w2～3，…，an-1与an得到wn-1～n，从而有n-1期地表下沉量：

然后，将这n-1期地表沉降结果与实地地物进行严格配准，并以第1 期的地表沉降结果w1～2为基准按时序依次叠加，即w1～3=w1～2+w2～3，w1～4=w1～2+w2～3+w3～4，…，w1～n=w1～2+w2～3+…+wn-1～n，从而得到从t1到tn时段内的地表时序累积下沉量：

基于式(4)和式(5)结果分析在注浆充填影响下的矿区地表动态沉降规律。

2 研究区概况与研究数据

2.1 研究区概况

选取中国安徽省北部淮北市境内的淮北矿业集团临涣煤矿Ⅱ2 采区作为研究区，其地表覆盖较多的植被、农作物，地理位置如图1所示，研究区内布置有Ⅱ923 里、1021 和Ⅱ1021 共3 个采煤工作面，均使用综合机械化长壁垮落法结合覆岩关键层注浆充填开采，详细参数见表1。图1 中，黑色多边形区域表示工作面，在研究时段(2016-11-21—2017-02-25)内，有注浆活动的注浆孔共有10个，用同心环状符号表示。

2.2 研究数据

图1 研究区内的工作面与注浆孔位置Fig.1 Working-face and grouting-hole positions in study area

经大量试验发现，受研究区地表浓密植被的影响，春、夏、秋季得到的C波段雷达影像差分干涉对相干性很低，很大程度上影响了DInSAR正确提取地表形变相位，导致结果精度较差[22]，而且矿区地表大梯度沉降会导致相位缺失，最短重复周期的SAR 影像对可以使时间失相关减至最小，因此，本文在尽量满足时间基线最短的前提下获取研究区2016-11-21—2017-02-25 的8 景Sentinel-1A影像，影像数据主要参数见表2，组成的7个干涉对及其相关信息如表3 所示。可见：除了2017-01-08—2017-02-01 的影像对时间基线为24 d 外，其余影像对的时间基线均为Sentinel-1A 最短重复周期12 d。用于与研究结果作对比、验证的现场实测水准数据由研究区煤矿提供，其时间范围为2016-11-16—2017-02-26，比研究时段多6 d，地面水准测量点数为40个。

3 地表动态沉降规律分析

3.1 下沉盆地时序沉降分析

利用式(1)～(4)对获取的8 景Sentinel-1A 影像进行DInSAR处理，得到7期地表下沉速率，结果如图2所示，对其统计得到单期的最大下沉量Wmax与最大下沉速率Vmax见表4。由图2 可知：在研究时段内，研究区地表形成了明显的下沉盆地，与工作面位置相吻合；随着工作面的推进，每期地表下沉速率最大点也向工作面推进方向移动，并且符合在近水平煤层开采条件下，地表沉降具有超前影响以及最大下沉速率点滞后于开采活动中心的规律[1]。结合表4 可知：1)第5 期的地表最大下沉量Wmax最大，达47 mm，但其最大下沉速率Vmax在监测时段内最小，为1.96 mm/d，这是由于第5 期的DInSAR 时间基线为24 d，在此期间内发生的地表沉降应该远大于此期监测结果，但其实际下沉量超出了哨兵1号数据的最大监测能力，因而只监测到了47 mm的下沉量；2)第1期的地表最大下沉量最小，仅为30 mm；3)第3期的地表最大下沉速率最大，为3.50 mm/d。综合各期数据分析可得，平均地表最大下沉量为38 mm，平均最大下沉速率为2.89 mm/d。

表1 研究区工作面参数Table 1 Parameters of working faces in study area

表2 影像数据主要参数Table 2 Main parameters of data

表3 干涉对组成及基线信息Table 3 Composition of interferometry pairs and baseline information

图2 地表时序下沉速率Fig.2 Time series subsidence velocity of surface

结合实际注浆资料，第4期与第5期之间经历了较剧烈的注浆活动，因而相较于前几期，第5期地表沉降中心区下沉速率明显降低，并且各期最大沉降速率围绕平均最大沉降速率波动较为剧烈，在第3期到达形变速率峰值后开始减小，至第5期减小为形变速率低谷，而第6期又到达形变速率次高峰，表明受注浆充填活动影响的采区地表沉降与无注浆充填也即岩层自然垮落引起的采区地表沉降规律有较大差异(正常情况下单点下沉速率呈现类似二次函数抛物线形状)。另外，因受到注浆充填的影响，即使有些建筑物处于开采中心区，相对于未受到注浆充填影响的地表，其中心及其附近地表下沉速率也要小很多，因而也使得DInSAR监测的每期地表下沉速率图以开采工作面走向线为中轴线表现出明显的非对称性(见图2)，这些结果均表明注浆减沉效果显著。

3.2 下沉盆地时序累积沉降分析

为研究矿区地表长期动态沉降，将7幅地表沉降图所示的阶段沉降结果按时间先后顺序依次进行叠加，得到地表时序累积沉降图，结果如图3所示，对其统计得到的时序累积最大下沉量与最大下沉速率见表4。从图3 可以得出：在监测期间，

地表下沉盆地的范围随时间推移有逐步扩大的迹象，主要表现在1021 与Ⅱ1021 工作面新开采的区域，并且在这2个工作面开采活动中心地表沉降异常明显，形成了2个沉降中心，下沉量最大的沉降中心位于Ⅱ1021工作面采空区上覆地表，累积最大下沉量达到18 5 mm，另一个沉降中心位于1021工作面西北侧地表，累积最大下沉量达到169 mm，以下沉量10 mm 为下沉边界起算点，最终下沉盆地面积约为1.35 km2。由于Ⅱ923 里工作面已于2016-04-28 收作，距离开始时间(2016-11-21)已达7月，故II923里工作面上覆地表仅在靠近2个开采工作面的区域有较小沉降，最大下沉量为42 mm。研究区地表最终累积最大下沉速率为1.93 mm/d，平均累积最大下沉速率为2.18 mm/d，均大于地表移动持续时间三阶段划分阈值1.67 mm/d，表明该研究时段内地表下沉盆地处于下沉活跃期(开始期为从下沉10 mm 开始至下沉速率达到1.67 mm/d；活跃期为下沉速率大于1.67 mm/d；衰退期为下沉速率小于1.67 mm/d 且6月内下沉量小于30 mm)。

表4 地表时序沉降统计Table 4 Statistics of time series subsidence of surface

图3 地表时序累积沉降Fig.3 Time series cumulative subsidence of surface

基于图3并结合实际注浆资料分析，由于1021工作面及其东南侧储灰场在研究时段内持续实施注浆充填，故位于1021 工作面采空区上方的地面按无充填开采沉陷规律来说，应该处于沉降中心区，但事实上其下沉量比未受到注浆减沉影响的1021 工作面西北侧区域下沉量小；另外，研究区地表下沉也无明显向东南侧扩张的迹象，下沉盆地无论是下沉范围还是下沉量都表现出显著的非对称性，这些均有力表明注浆充填减缓地面沉降的效果显著。

3.3 下沉盆地主断面时序累积沉降分析

为进一步量化了解研究区地表动态沉降规律，对下沉盆地主断面上的沉降特征进行分析。在研究区1021和Ⅱ1021工作面上选取了经过2个工作面累积沉降最大值中心的线段AA1近似作为走向线，经过Ⅱ1021 工作面累积下沉量最大值中心的线段BB1作为倾向线1，经过1021 工作面累积下沉量最大值中心的线段CC1作为倾向线2，如图4 所示。基于这3条线段提取出3组位于下沉盆地主断面上的时序累积沉降曲线，结果如图5 所示(其中坐标横轴表示以A，B，C 为起点，主断面上各点至起点的距离)。

图4 所选主断面位置示意图Fig.4 Location map of selected main sections

1)在AA1的走向主断面，以下沉量10 mm为下沉边界起算点，则下沉盆地沉降范围为距离A 点330～1 970 m 之间，记为A-330～1 970 m，长约1 640 m。Ⅱ1021 工作面的最大沉降中心和1021 工作面西北侧的最大沉降中心分别位于850 m附近和1 350 m 附近，对应的累积最大下沉量分别为185 mm和169 mm。另外，在800 m处的地表本应处于最大下沉中心区，却存在凸起现象，结合实地注浆资料以及水准测量数据分析，这是由于该处紧邻Ⅱ1021 工作面2 号注浆孔，而该注浆孔在DInSAR 第4 期与第5 期之间进行过剧烈的注浆活动，因此，地表出现了一定程度的抬升，进而使5～7 期时序累积沉降曲线在此处均存在局部峰值，与紧邻此注浆孔的实地水准测点数据吻合。由于950～1 150 m 范围内的盆地距离2 个开采活动中心较远，因而，其下沉量比最大沉降中心的下沉量小。1 500～1 900 m 范围内的地表升降波动也由第4～5期之间的地表注浆活动引起。

2)在BB1的倾向主断面，以下沉量10 mm为下沉边界起算点，则其下沉盆地范围为B-470～1 570 m，长约1 100 m。在观测期间，Ⅱ1021 工作面地表最大沉降中心从第1期的16 mm下沉量，以平均下沉速率2.01 mm/d 下沉到最后一期的185 mm，1 100 m 附近的5～7 期地表升降波动也由地表注浆充填引起。

图5 下沉盆地主断面时序累积沉降曲线Fig.5 Time series cumulative subsidence curves of main sections of subsidence basin

3)在CC1的倾向主断面，以下沉量10 mm为下沉边界起算点，则其下沉盆地范围为C-500～1 620 m，长约1 120 m，与BB1主断面的沉降长度接近。在观测期间，1021 工作面西北侧地表最大沉降中心从第1期的20 mm下沉量，以平均下沉速率1.76 mm/d 下沉到最后一期的169 mm。在5～7期，出现在1 300 m 附近的曲线峰值与1 500 m 附近的曲线低谷同样是由于注浆活动引起的，并且，离最大沉降中心更近的1 300 m附近的地表下沉量小于1 500 m附近的地表下沉量，进一步证明了注浆减沉效果显著。

由BB1倾向主断面与CC1倾向主断面对比可得，虽然Ⅱ1021工作面属于重复采动并且倾向长度要比1021工作面的大，但因其为新开采的工作面，地表沉降有一定的延迟时间，即存在起动距，因而在前4期监测结果中，Ⅱ1021工作面最大沉降中心的下沉量一直比1021 工作面西北侧最大沉降中心的小，但随着采掘工作面继续向前推进，从第5期起，Ⅱ1021 工作面采空区走向长度达到155 m，Ⅱ1021 工作面最大沉降中心的下沉量大于1021 工作面西北侧最大沉降中心的下沉量。

总体上看，这3组时序累积沉降曲线具有双工作面同时开采的地表下沉特征，并且地表最大下沉量均显示出增加的趋势，表明研究区未达到充分采动状态。另外，上述结果也间接反映了若研究区地表未进行注浆减沉，则下沉盆地范围与下沉最大值皆会增大。

4 精度验证与误差分析

4.1 精度验证

为了验证利用TSS-DInSAR方法来监测分析矿区地表动态沉降规律的可靠性，采用了地面实测水准数据进行对比验证，水准测点位置见图6。因本研究时段为2016-11-21—2017-02-25，而水准数据监测时段为2016-11-16—2017-02-26，二者时间间隔并不完全相同，为了尽可能减少因时间基准不一致带来的误差，根据各位置时间段内地表下沉速率等相关信息将水准数据起止日期内插至与研究时段相同的日期，即水准数据与TSSDInSAR监测结果起止时间一样。

图6 水准验证数据点位Fig.6 Point positions on verification data of leveling

TSS-DInSAR与水准测量的地表累积下沉量对比如图7 所示。由图7 可知：除了5～8 点差值较大(差值分别为-109，-96，-121 和-75 mm)以及38号，39 号点曲线趋势不一致外，其余点位差值都比较小，拟合度较好，均方根误差RMSE 为35.8 mm。这是由于5～8 号点位于下沉盆地中下沉梯度最大的位置，以4号和5号点为例，其下沉梯度约为2.0 mm/m，这远超出C 波段Sentinel-1A 数据可监测到的最大下沉梯度1.4 mm/m[21]，导致这些点的DInSAR 监测结果与水准测量结果差值较大。此外，对于39号点，其位置比38号点更靠近最大沉降中心，但水准测量数据却出现抬升异常，综合研究区的多次DInSAR处理试验结果，该异常可能是由水准测量中人为误差传递导致。于是，以2.5倍RMSE为限差，剔除观测相位失相干点5～7 后，得到新的RMSE 为19.7 mm，精度提高了45%。

为进一步验证本研究精度，将去除失相干点5～7 后的TSS-DInSAR 累积下沉量与水准数据累积下沉量进行回归分析，并进行了F检验，结果如图8所示。从图8可知：决定系数R2达到0.910 4，二者之间的标准化残差值95%分布于-2～2，且服从正态分布；F′=355.53，F=F0.05(1,35)=4.13，F′＞F，且p＜0.001。这些结果表明：去除失相干点后的最终DInSAR累积下沉量与水准测量累积下沉量之间具有很强的线性关系，回归模型具有高拟合优度与置信度，利用TSS-DInSAR方法来监测分析矿区地表动态沉降规律是切实可靠的，满足实际应用需求。

图7 TSS-DInSAR与水准测量的累积下沉量对比Fig.7 Comparison of TSS-DInSAR cumulative subsidence and leveling cumulative subsidence

4.2 误差分析

沉陷中心区域TSS-DInSAR监测结果与水准测量结果差值较大并且研究中存在误差，其产生原因主要有：

1) 由于Sentinel-1A 数据波长较短、分辨率较低，能够监测到的最大下沉梯度远小于研究区地表实际的最大下沉梯度，因此在下沉梯度较大的地区出现相位丢失现象，难以探测出准确的下沉量。

图8 TSS-DInSAR与水准累积下沉量回归分析Fig.8 Regression analysis on TSS-DInSAR cumulative subsidence and leveling cumulative subsidence

2)地面实测水准数据起止时段与DInSAR监测起止时段有些差异，在一定程度上会影响数据拟合验证分析。

3)单元范围不一致。DInSAR 技术获取的是1个面的结果，而水准测量针对的则是单点。具体而言，DInSAR技术测量的结果是影像上一个分辨单元的沉降量，受影像分辨率的限制，地表上某点的沉降情况会受到分辨单元范围内所有点的影响。因此，用1个“点”来评价1个“面”本身就存在着误差[2]。

4)紧邻研究区西边的沉陷积水区，并且所在区域覆盖大量植被，这降低了DInSAR干涉对的相干性，进而影响解缠结果与最终提取的地表沉降结果。

5)未考虑水平移动分量对雷达视线方向形变的贡献，直接通过投影关系将视线向形变转换至垂直向形变，这必然会带来误差。

6)其他因素。包括所采用的DInSAR处理算法本身存在误差、大气延迟效应等因素。

综合考虑以上误差来源，在仍使用高性价比的Sentinel-1 卫星影像作为数据源的前提下，减小上述误差的方式主要有：

1) 引入Sentinel-1B 卫星影像数据，将时间分辨率提高1倍，从而可以显著减少因矿区地表形变梯度过大而导致的时空失相干，并且能更精细的提取矿区地表形变。

2)优化DInSAR处理算法，比如增大图像配准窗口大小以及配准窗口数量来提高配准精度，对得到的干涉图进行更为精细的滤波以尽量保证下沉盆地中的干涉相位连续且下沉盆地外的噪声较少，相位解缠可使用基于不规则格网的方法来减少MCF解缠结果存在较多的相位不连续现象。

3)在干涉时尽量避免雨雪天气所获取的影像数据，以提高得到的干涉图的质量并减少大气延迟相位的组分，另外，可以借助外部GPS 数据来减弱大气延迟误差。

4)反演出地表三维形变，针对某个方向上的真实形变值进行对比。

5 结论

1)利用TSS-DInSAR方法监测到的研究区地表下沉盆地在研究时段内处于下沉活跃期，其平均最大下沉量为38 mm，平均最大下沉速率为2.89 mm/d，最终累积最大下沉速率为1.93 mm/d，平均累积最大下沉速率为2.18 mm/d，最终下沉盆地面积约为1.35 km2，且在Ⅱ1021工作面和1021工作面西北侧形成了2个沉降中心，二者的最终累积最大下沉量分别为185 mm和169 mm。

2)受注浆充填影响的采区建筑物及其附近地表的下沉速率、下沉量以及下沉范围均显著减小。下沉盆地以工作面走向线为中轴线表现出明显的非对称性，表明结合注浆充填的采煤工艺可以有效减缓控制区域内的地表沉降，有注浆充填活动的采区地表沉降与岩层自然垮落引起的采区地表沉降规律存在较大差异。

3)TSS-DInSAR监测结果与地面实测水准数据就最终累积下沉量进行验证精度，剔除3个相位失相干点后，RMSE 为19.7 mm，回归分析R2为0.910 4，F 检验p＜0.001，表明基于TSS-DInSAR方法的矿区地表沉降监测结果可靠性较高，具有实际应用可行性。同时也反映出Sentinel-1A 数据能监测到地表微小沉降且测量精度很高，而对下沉梯度大于1.4 mm/m 的地表沉降，其监测值会偏低。