基于SBAS-InSAR技术的采空区残余变形规律分析

梁思语，胡海峰

(太原理工大学矿业工程学院，山西 太原 030024)

2020年，由于去产能政策调整和煤炭行业结构优化，全国关闭煤矿达428家，其中山西省关闭了32家[1]。而关闭矿山的环境恢复、改造和发展已成为矿业城市社会经济发展的重要问题[2-3]。当矿井排水设备停止运行时，地下水位上升将改变采空区破碎岩石的应力和承载力。这导致地表发生二次变形，威胁到在封闭矿井顶部建造的建筑物的安全[3-4]。矿山关闭前后发生的形变相互关联，随着时间的推移，地表变形呈现出连续性。因此，获取矿区关闭或者停止开采后发生的地表变形，对于确定地表变形规律、完善开采沉陷理论、土地利用规划和建筑物稳定性评价具有重要意义。传统的监测手段(如水准测量、全球卫导航卫星系统等)成本高、耗时长、地表变形的空间覆盖度低，只能获得点位的变形，很难获得大范围的地表变形的时空发展规律，而且监测点很难长期保存，受外界环境影响较大。与传统方法相比，合成孔径雷达干涉测量(InSAR)具有效率高、成本低、范围大、所有天气条件下可全天候使用等优点[5]。因此，SAR图像对于探明矿区关闭前后地表变形的位置和范围具有很好的应用价值。基于DInSAR技术的小基线集，InSAR技术实现了长时间、大范围区域、连续微小形变的监测，且很大程度提高了地表形变监测的精度，得到了广泛的应用[6]。张艳梅等[7]将多期哨兵数据采用SBAS-InSAR技术对西安市城区地表进行沉降监测，并利用已有研究资料，证明了方法的有效性；李达等[8]通过采用SBAS-InSAR技术对陕西一矿区进行了残余变形的监测，并用DInSAR的结果验证了SBAS-InSAR监测结果的精度；余昊等[9]用SBAS-InSAR技术对徐州西部地区关闭矿井的地表变形规律进行了研究，发现矿井关闭的地区地面有抬升的趋势；何荣等[10]对大采深条带开采的矿区进行地表沉降的监测研究，发现条带工作面开采会诱发相邻老采空区地表再次发生变形。

现有研究多针对矿区开采过程，对工作面开采结束后地表变形的研究较少，且对残余变形的研究大多只针对某一个工作面，时间跨度小，无法反映矿井关闭后长期的变化情况。本文以山西省某矿区为研究对象，工作面的停采时间跨度为2005—2018年，选取24景Sentinel-1A影像，时间跨度为2018年1月—2018年6月、2018年10月—2019年4月，应用SBAS-InSAR技术对该矿区进行了地表形变监测，获取了自2018年1月—2018年6月和2018年10月—2019年4月该地区地表形变信息。通过对2018年1月—6月的InSAR监测结果与实测数据进行对比分析并计算差值，发现SBAS-InSAR反演出的下沉变化情况与实际情况相符，结果表明采用SBAS-InSAR技术反演矿区残余变形具有较好的可靠性。对2018年10月—2019年4月的反演结果进行分析，得到了其形变规律，为预测和评价采空区残余形变提供了基础。

1 研究区概况

研究区位于山西省长治市某矿区内，该研究区域的地表由第四系覆盖，煤层上依次是二叠系山西组、下石盒子组上段、下石盒子组下段、上石盒子组下段。区内地形以低山-丘陵为主，地势西高东低，为典型的北方山地丘陵地貌，属东亚季风区半干旱大陆性气候，四季分明，夏季多雨，春秋季多风少雨，冬季寒冷。

研究区域选取的是该矿区内已开采的9101工作面、90105工作面、90207工作面～90209工作面，上述工作面均为已经停采工作面，开采时间为2005—2018年，位置如图1所示。本次实验选取90208工作面、90209工作面、90105工作面以及9101工作面进行研究分析，开采煤层为9号煤层和10号煤层，开采深度约为150 m，开采平均厚度为2.4 m，煤层倾角为2°～5°。本次研究工作面的范围跨度大，停采时间涵盖了10 a以上，可以获取比较完整的残余形变规律。

图1 研究区位置Fig.1 Location of the study area

2 数据与方法

2.1 数据源

实验数据为24景分辨率为5 m×20 m的Sentinel-1A卫星影像数据，影像时间跨度为2018年1月—2018年6月，2018年10月—2019年4月，数据为C波段，波长为5.6 cm，具体的影像信息见表1。同时，为了减少地形相位的影响，实验选取了分辨率为30 m×30 m的DEM数据。

表1 Sentinel-1A卫星影像数据Table 1 Satellite image data of Sentinel-1A

2.2 方法

SBAS-InSAR技术能够有效地削弱时空失相干对时序沉降结果的影响，从而使得到的形变图在时间和空间上更为连续[11]。 由于哨兵数据的波长是C波段(5.6 cm)， 矿区夏季植被茂密， 为避免出现严重的失相干现象，本次研究将数据处理分为两次重复性实验，避开失相干严重的时期。其中处于采空区下沉时期的影像做一次数据处理，获得与水准测量时间重合的形变情况；将沉降衰退期结束后的影像做第二次数据处理，获取到工作面的残余形变加以分析。

首先，对影像进行预处理，将ESA提供的原始数据转换成SLC格式的数据，提取覆盖研究区的公共burst，并对研究区域进行裁剪。利用轨道参数、SAR影像的强度信息来配准图像，使配准的精度达到千分之一个像元。将主影像、从影像分别去斜，再利用精密轨道数据对轨道信息进行修正。为了提高监测的精度，将时间基线阈值设置为90 d，空间基线阈值设置为150 m，通过对干涉图滤波、相位解缠、基线精细化后，剔除掉相干性差的干涉对，分别组成30对和33对干涉对进行差分干涉处理，干涉对情况如图2所示。为了避免低相干点引起的误差，在差分干涉相位图中，利用影像的相位稳定性、振幅离散指数和空间相干性等作为影像参考因子，寻找高相干的目标点。在此基础上，对相干点进行差分、滤波、解缠[12]，通过逐次迭代，去除高程相位、大气误差、噪声误差，提取真实的形变相位，最后利用奇异值分解获得高相干点的时间序列的沉降速率。将解得的各时段相位速率在时间域上积分，即可得到整个观测时段的形变时间序列。

图2 时空基线组合Fig.2 Spatiotemporal baseline combination

3 可靠性分析

为了验证SBAS-INSAR反演结果的可靠性，本次实验共选取工作面开采过程中布设的观测线上的九个点用于验证Sentinel-1A的监测结果，水准点的位置如图3所示。由于从InSAR获得的形变是LOS向的形变，精度分析时通过将LOS形变转换为垂直方向的形变而不考虑水平运动来进行的。

图3 沉降衰退期SBAS累计沉降值Fig.3 Cumulative settlement value of SBAS during settlement decline period

提取InSAR结果中对应水准点位置的沉降值[13]，InSAR处理结果与水准测量结果对比情况如图4所示。由图4可知，1号点、2号点、4号点～8号点的Sentinel-1A图像产生的时间序列累计沉降趋势和大小与水准数据均一致。虽然3号点和8号点的沉降趋势与水准数据一致，但在大小上存在显著差异。3号点的变形最大值达到62 mm，9号点的变形最大值达到81 mm，而InSAR监测得到的变形值只达到35.7 mm和46.9 mm。

图4 SBAS-InSAR结果与水准测量结果对比Fig.4 Comparison between results of SBAS-InSAR and leveling

为了定量评价监测结果的准确性，将矿区的InSAR时间序列形变值线性拟合后进行插值，得到与水准数据时间重合的沉降值。采用最大偏差(MaxD)、最小偏差(MinD)、均方根误差(RMSE)、标准差(STD)和相关系数(R2)对InSAR结果的准确性进行评价，评价结果见表2。

由表2可知，除Z2点外，所有点的InSAR结果与水准测量结果的相关系数都达到了0.9，Z2点相关性较低是由于Z2点所处位置坡度变化较大，沉降值受到了地形影响，而SBAS-InSAR监测结果由于分辨率较低，并不能准确地反映Z2点这一点位的形变值。 除Z9点外，所有点的精度均小于20 mm。所有点的平均均方根误差和标准差分别为11.05 mm和20.22 mm，整体来看，影响InSAR技术影响监测结果精度的因素有以下几个方面。①监测期间只采集了少量的Sentinel-1A图像，且由于处于工作面下沉的衰退期，进行水准测量的次数较少。虽然C波段具有较大的可检测变形梯度，且对植被和时间去相关相对不敏感，但Sentinel-1A图像之间较长的时间间隔导致了较大的形变梯度。②进行对比分析的数据处于衰退期，沉降速率在逐渐减小，使用线性插值来比较SBAS-InSAR的结果与水准数据将不可避免地导致较低的精度水平。③每个水准点的InSAR形变值因为分辨率的原因，并不能准确的反映水准点也就是观测点处的形变值，而是通过对水准点附近的相干点的变形进行平均得到的，这影响了对比结果。尽管仍存在偏差，但Sentinel-1A图像监测结果的精度可以达到或优于厘米级，足以监测开采引起的地表残余变形。

表2 SBAS-InSAR结果精度分析Table 2 Accuracy analysis of SBAS-InSAR results

4 试验结果及分析

4.1 沉降速率

通过对影像处理，获得了该矿区残余形变时期2018年10月—2019年4月期间沿雷达视线方向的年平均沉降速率，将雷达视线方向的沉降速率换算到垂直方向上，得到形变速率图如图5所示。矿区的右侧紧邻处于开采过程中的工作面，根据开采沉陷中沉陷盆地影响半径计算公式见式(1)。

(1)

式中：H为开采深度，m；tanβ为主要影响角正切值。

图5 研究区域垂直沉降速率Fig.5 Vertical settlement rate in the study area

90207工作面和90208工作面的停采线方向受到其影响，无法准确获得其沉降规律，90209工作面与其距离大于其开采沉陷半径，未受到其影响，可以对整个工作面进行规律分析。结果显示，2018年2月停采的90209工作面沉降速率较大，下沉速度为50～88 mm/a，2017年和2016年停采的90208工作面、90207工作面沉降速率较小，为30～59 mm/a。90105工作面开采时间为2013年，沉降速率约为10 mm/a。早期开采的9101工作面沉降速率接近0，说明地面已经基本稳定，附近没有正在进行开采的工作面，在9号煤层和10号煤层下方的煤层也未进行开采，整体处于稳定状态。

4.2 时序累计沉降值分析

图6是以2018年10月28日为起始时间，其他时间相对于起始时间的形变结果。右侧沉降区域的范围在向南扩大，90209工作面未受正在开采的工作面的影响。90209工作面于2018年1月开采结束，下沉值连续6个月小于50 mm，但衰退期结束后仍在不断下沉。由图6可以看出，自2018年10月28日起，90209工作面出现明显沉降，沉降区域在逐渐变大，逐渐出现一个下沉盆地，最终盆地范围与采空区的范围相近，沉降值最大为32 mm，从下沉量来看，90209工作面的残余下沉仍会持续一段时间。90208工作面和90207工作面忽略邻采的影响，靠近开切眼方向的沉降值为10 mm左右，最大的沉降值为13 mm，90105工作面沉降值较小，仅8 mm，9101工作面基本没有沉降。通过分析可知该地质采矿条件下工作面在开采结束1 a内沉降较为明显，在开采沉降结束1 a以上，5 a以内，仍然会有较小的沉降，开采结束5 a以上，沉降值接近于0。

图6 研究区域时序累计沉降Fig.6 Time series cumulative settlement in the study area

图7展示了观测线位置。如图7所示，选取工作面在开采过程中布设的观测站作为时序累计沉降分析的观测点，并延伸至下沉盆地边缘生成东西向和南北向的剖面线，提取剖面线的时序沉降值。忽略受植被和噪声影响的观测点，根据90209工作面的走向剖面线和倾向剖面线上的点，绘制下沉曲线如图8所示。

图7 观测线位置示意图Fig.7 Position diagram of observation line

由图8可以看出，90209工作面最大下沉值是32 mm，逐渐趋于稳定。从下沉曲线整体上看，在工作面内有下沉，在工作面边界以外基本无下沉或者下沉值较小。工作面的边界位于27号点和28号点以及58号点和59号点之间，通过走向剖面线提取的线上可以看出，在31号点和55号点，即工作面边界附近的点出现了较大的下沉值，这是由于在开采结束后采用全部垮落法处理顶板，在工作面的边缘会由于上方岩石的抗拉性等出现空洞，可能会由于上方覆岩体应力作用下，下方的空洞处的岩石失稳，空洞被岩石压填，覆岩移动，造成该区域下沉较大。

图8 下沉盆地剖面时序变化Fig.8 Time series change of subsidence basin profile

由图8可以看出，下沉曲线保持了较好的连续性。采用二维高斯曲线拟合沉降值，并根据曲线形态分析下沉盆地的沉降特征，拟合结果如图9所示。根据拟合结果可知，沉降值与高斯拟合曲线的相关系数为0.911 8，高斯拟合曲线与SBAS结果达到较高拟合度。该矿地表变形特征在空间上类似于高斯曲线，表明矿井地表运动特征在一定范围内符合概率积分法模型的特征。

图9 下沉盆地高斯拟合结果Fig.9 Gauss fitting results of subsidence basin

4.3 研究区单点沉降特征分析

为了分析研究区域的沉降性质，选取下沉盆地中心的6个点作为研究对象，其中1号点～3号点位于90209工作面内，4号点～6号点位于90208工作面内，考虑相邻工作面开采的影响，4号点～6号点选取的位置靠近开切眼。提取出这6个点的时间序列沉降值，并用线性拟合模型拟合下沉值与下沉时间的关系。拟合后发现，90209工作面的下沉值与下沉时间呈明显的二次相关关系，但是二次项的系数非常小，接近于线性关系。90208工作面的下沉值与下沉时间为一次线性关系，得到的线性关系式如图10所示。

图10 研究区域单点沉降特征Fig.10 Characteristics of single point settlement in the study area

根据多个工作面下沉值的提取发现，停采时间在1 a以内的工作面的下沉值呈二次相关关系减小，停采时间在1 a以上的工作面下沉值趋于线性形变，停采时间越长，短期内的下沉值越符合线性关系，最终下沉值都趋近于0。

5 结 论

本文采用SBAS-InSAR方法对采空区的残余变形进行研究，通过采用SBAS-InSAR方法对24景Sentinel-1A影像分别进行两次数据处理，分别得到了工作面开采结束后衰退期的累计沉降值，并与水准观测得到的沉降值进行了精度对比；同时得到了地表衰退期结束以后的年平均沉降速率、时序累计沉降值等高空间分辨率的时间序列形变信息。通过分析，得到结论如下所述。

1) 将SBAS-InSAR得到的监测结果与水准测量得到的结果进行对比分析，发现二者得到的结果大小、趋势一致，所有点的平均均方根误差和标准差分别为11.05 mm和20.22 mm，并且与开采沉陷规律一致，说明SBAS-InSAR在进行矿区残余变形监测方面具有较高的可靠性。

2) 研究表明，煤层采深和采厚均较小的工作面，在无外界因素干扰(如重复开采、地震等)条件下，开采结束后一年内沉降值与沉降时间符合二次相关关系，并逐渐趋于线性关系；开采结束1 a以上，沉降值与时间趋于线性关系；开采结束时间达到2 a以上，5 a以内的工作面仍然会有微小的沉降；开采结束时间在5 a以上的工作面，年沉降值接近于0。

3) 工作面的残余变形经过线性拟合，符合二维高斯曲线拟合的结果，说明工作面的残余形变仍然符合概率积分法。

综上所述，SBAS-InSAR技术在研究采空区残余变形时具有很好的应用前景，可以有效监测到矿区在开采结束后的地形的变化规律，为今后进行重复开采等工作提供参考。