0 引言
长期以来，我国富煤少油。煤炭作为能源消费的主体，其重要性不言而喻。近年来，随着国家政策调整，大量矿井被关闭，有关报道显示：2006—2008年全国关闭小煤矿11 618个，预计到2030年，关闭矿井可达15 000个[1-2]。目前，很多矿井在关闭以后缺乏持续的监测，处于无人监管的状态[3]。矿井关闭后，煤岩体在应力、地下水、氧气等多种因素的作用下发生风化劣化、强度降低，采动破裂岩体的应力和承载能力发生改变，易导致采空区地表发生二次形变。这种形变具有隐蔽性、复杂性、突发性和长期性等特点，对工程建设造成巨大的潜在威胁[4]。
传统的测量方法如水准测量不仅耗费大量的人力物力，而且只能对矿区内少数离散点进行监测，不能反映出矿区地表的整体形变结果。20世纪70年代发展起来的合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar, InSAR)技术以其全天时、全天候、高精度、大范围覆盖等突出优点很大程度上弥补了传统测量方法的不足[5-6]。常规InSAR技术易受时空去相干及大气效应的扰动，在长时间地表形变监测中难以获得高精度的测量结果。为解决这一问题，各种基于高相干点目标的经典时序InSAR技术被相继提出，主要包括永久散射体方法[7]、最小二乘法[8]、小基线集(Small Baseline Subset, SBAS)方法[9]、相干目标方法[10]等。早期的时序InSAR技术主要用于城市地表沉降监测，随着该技术应用领域的拓展，在矿区形变监测领域也取得了丰富的研究成果[11-14]。D. Walter等[15]利用永久散射体干涉测量对德国某开采中煤矿的地表形变进行了监测；S.Samsonov等[16]基于167景合成孔径雷达(SAR)影像，利用融合的时序InSAR方法监测了位于德国和法国边界因煤炭开采引起的地表沉降；C. Zhao等[17]将InSAR与偏移量追踪相结合的方法应用于矿区形变监测；L. Zheng等[18]将分布式目标方法与永久散射体方法相结合，利用L波段的SAR影像绘制了鄂尔多斯市某矿区危险区域边界。
现有研究多是集中在矿区开采过程中的地表形变监测，针对矿井关闭后的地表形变监测研究较少。为探明矿井关闭后的地表形变规律，本文基于2015-12-21—2019-12-24的62景Sentinel-1A影像，采用短基线集干涉测量(SBAS-InSAR)技术对徐州西部关闭矿井地表进行形变监测，得到徐州西部矿区关闭矿井后4 a内的时间序列地表形变信息和形变规律，研究结果可为安全、有效、合理地利用关闭矿井土地提供决策依据。
1 SBAS-InSAR技术原理
SBAS方法通过设定时空基线阈值将所获得的SAR数据组成若干个集合，利用最小二乘法获取每个集合里的地表形变时间序列。假设N+1幅影像按照时间序列t0,t1,…,tN排列，选取一幅主影像对其他影像进行配准，按照小基线干涉组合原则，生成M幅多视差分干涉图，M的取值满足以下条件：
(1)
对于在tA，tB(A，B=0，1，…，N，A≠B)时刻生成的第K幅干涉图，其任意像元的干涉相位值为
δØK(x,r)=Ø[tB,x,r]-Ø[tA,x,r]≈
(2)
式中：Ø[tA,x,r],Ø[tB,x,r]分别为tA，tB时刻相对于初始时刻t0的形变相位，x，r分别为方位向与距离向坐标；λ为雷达波长；d(tA,x,r)和d(tB,x,r)为相对于初始时刻的视线方向(Line of Sight, LOS)的形变量。
在所有干涉图中，设按照时间顺序排列的主影像序列IE=[IE1,IE2,…,IEM]，辅影像序列IS=[IS1,IS2,…,ISM]，tIEj,tISj(j=1,2,…,M)分别为主影像与辅影像时间序列，且tIEj>tISj，则干涉相位为
δØj=Ø(tIEj)-Ø(tISj)
(3)
式中Ø(tIEj),Ø(tISj)分别为对应时刻主影像与辅影像的形变相位。
由式(1)可知，一共生成M个干涉对，故根据式(3)可得到由M个等式构成的方程组，其矩阵形式可表示为
AØ=δØ
(4)
式中：A为M×N的系数矩阵，其中主影像对应的系数为1，辅影像对应的系数为-1，其余值为0；Ø为待求的相位矩阵；δØ为差分干涉图的相位矩阵。
当所有的影像都被分在一组且M≥N时，矩阵A的秩为N，即可利用最小二乘法求解出Ø的估计值
(5)
若影像被分为多个基线集合，即A是秩亏的,则对应的ATA是奇异矩阵，可利用奇异值分解法求其最小范数的最小二乘解，进而获取地表形变信息。
由于试验区域为关闭矿井，总体地表形变量较小，在不考虑水平移动的情况下，可利用式(6)得到垂直方向的沉降量W:
W=d(tend,x,r)/cos θ
(6)
式中：tend为最终时刻；θ为卫星入射角。
2 研究区域与数据源
徐州西部矿区位于铜山区与泉山区境内，该区域为旧黄河泛滥形成的冲积平原，地势较为平坦，属南温带鲁淮区，具有长江流域和黄河流域夏季多雨、冬季干寒的过渡性特点。徐州西部矿区内，石炭系、二叠系地层是其含煤地层，在井田内均被第四系冲积层覆盖。整个矿区由张小楼矿、夹河矿、庞庄矿组成。张小楼矿井田东西长约5.0 km，南北宽约3.4 km，于2015年12月闭矿；夹河矿井田东西长约5.5 km，南北宽约4.5 km，于2015年12月闭矿；庞庄矿井田东西长约6.1 km，南北宽3 km，于2012年12月闭矿，总面积约为59.99 km2。夹河矿东南侧边界与拾屯矿接壤。研究区概况如图1所示。
图1 研究区位置及闭矿时间
Fig.1 Research area location and mine closure time
共收集了62景Sentinel-1A IW模式的SAR卫星影像，时间跨度为2015-12-21—2019-12-24。SAR影像方位向分辨率为5 m，距离向分辨率为20 m。影像除2017-05-08—2017-06-13、2018-12-27—2019-02-27、2019-11-18—2019-12-24间隔为36 d，2016-01-14—2016-03-02间隔为48 d外，其余获取时间间隔均为24 d。公共参考主影像时间为2016-10-04，部分影像获取时间和时空基线信息见表1。采用数字高程模型(DEM)为美国航天局提供的STRM 1数据进行分析，分辨率为30 m。
表1 部分影像获取时间和时空基线信息
Table 1 Partial image acquisition time and time-space baseline information
编号主影像时间辅影像时间空间基线/m时间基线/d12016-01-149.82422016-03-0299.77232015-12-212016-03-2672.99642016-04-1940.5120︙︙︙︙︙2662019-11-18170.9242672019-10-252019-12-2415.1602682019-11-182019-12-24157.336
3 关闭矿井地表形变规律分析
3.1 沉降速率分析
SAR数据的垂直基线和时间基线阈值分别设置为197 m与120 d，最终生成可用干涉对268对，时空基线分布如图2所示。
图2 干涉对时空基线分布
Fig.2 Time-space baseline distribution of interference pair
首先对生成的268个干涉对进行差分干涉处理，采用最小费用流方法进行相位解缠，根据干涉图的相干性及解缠结果筛选干涉对；然后选取地表稳定的控制点进行轨道精炼和重新去除平地效应；最后去除大气相位误差，采用最小二乘法计算得到平均沉降速率，结果如图3所示。
从图3可知：研究区大部分地表呈下降趋势；庞庄矿东南部E区域由于开采结束不久，地表沉降速率较大，可达-48 mm/a；夹河矿和张小楼矿的最大沉降速率分别为-25.8，-29.6 mm/a；庞庄矿东北部A区域(黑色圆圈范围内)出现抬升，抬升速率达11.3 mm/a。
图3 研究区垂直方向年平均沉降速率
Fig.3 Annual average subsidence rate in the vertical direction of the research area
3.2 时间序列形变分析
为进一步分析关闭矿井后各矿区地表时空形变规律，在获取沉降速率的基础上对其进行时间维度上的积分，获取各矿区时间序列累计沉降量的变化，结果如图4所示。因篇幅所限，将时间间隔调整为4个月。由图4可知：在闭矿1 a后，即2015-12-21—2016-12-27期间，夹河矿已经出现明显的下沉趋势，最大沉降量可达-126 mm；庞庄矿东南侧受拾屯矿开采影响也出现了下沉，最大沉降量可达-110 mm。2016-12-27—2018-12-17期间，庞庄矿、夹河矿、张小楼矿3个矿区大部分区域均出现明显的下沉现象，下沉影响范围进一步扩大，最大沉降量为-154 mm。此外，图4中的A、B、C、D区域出现了抬升现象，庞庄矿A区域最大抬升值达47 mm。至2019-12-24，张小楼矿与庞庄矿中部沉降区域范围继续变大，夹河矿与庞庄矿东南侧下沉范围趋于稳定，但沉降量还在继续变大，说明地表下沉还未停止，将继续下沉。
(m) 2015-12-21—2019-12-24
图4 研究区时序沉降分布
Fig.4 Time series subsidence distribution in the research area
为了进一步研究矿区沉降规律，在每个矿区已停采的工作面周围提取部分沉降点(位置如图1所示)进行时间序列分析，结果如图5所示。由图5可知：夹河矿在2015-12-21—2017-05-08期间，各点均出现明显下沉现象，最大沉降量可达-105 mm；2017-05-08—2017-10-11期间，地表整体沉降表现平缓；在2017-10-11后，所选点均出现缓慢抬升现象，出现这一现象的原因可能是由于地下水上升，造成土地回弹。张小楼矿在2015-12-21—2019-08-14期间下沉明显，最大沉降量可达-97 mm，2019-08-14后下沉趋势减缓。庞庄矿在2015-12-21—2019-07-21期间，所选点皆呈现出下沉趋势，最大沉降量为-88 mm，此后逐渐趋于稳定。
(a) 夹河矿
(b) 张小楼矿
(c) 庞庄矿
图5 沉降点下沉曲线
Fig.5 Sink curves of subsidence points
为进一步评估地表受损情况，在每个矿区沉降量最大的区域及其周边区域分别提取对应的观测线(如图1所示)，并利用式(6)将视线向的沉降转换为垂直向沉降。根据开采沉陷理论公式，即式(7)和式(8)，分别求出观测线上对应各点的倾斜与曲率值，其曲线如图6所示。
(7)
(8)
式中：i(k-1)-k为观测点k-1和k之间的倾斜，mm/m；Wk-1，Wk为观测点k-1和k的沉降量，mm；l(k-1)-k为观测点k-1和k之间的水平距离，m；u(k-1)-k-(k+1)为点k-1，k和k+1之间的曲率，mm/m2。
(a) 夹河矿倾斜分布
(b) 夹河矿曲率分布
(c) 张小楼矿倾斜分布
(d) 张小楼矿曲率分布
(e) 庞庄矿倾斜分布
(f) 庞庄矿曲率分布
图6 各矿区内倾斜和曲率分布曲线
Fig.6 Distribution curves of tilt and curvature of each mine
由图6可知：夹河矿观测线上最大倾斜与曲率出现在2019-08-14，分别为0.51 mm/m和0.012 4 mm/m2；张小楼矿倾斜与曲率最大值分别为0.689 mm/m和-0.029 mm/m2；庞庄矿倾斜与曲率最大值分别为1.70 mm/m和-0.039 mm/m2。以上形变值均在“三下采煤规范”地表允许Ⅰ级破坏变形值范围(倾斜≤3 mm/m，曲率≤0.2 mm/m2)以内，说明目前地表形变暂时不会对地表上方建(构)筑物造成危害。同时从图6可以得知，夹河矿、张小楼矿观测线上的倾斜与曲率均在2019-08-14达到最大值，随后逐渐变小并趋于稳定，而庞庄矿观测线上的倾斜与曲率在2019-12-24达到最大值，因受拾屯矿开采的影响，还有不断增大的趋势，预测未来可能对矿区周围的村庄产生一定影响。
3.3 沉降面积统计分析
为分析研究区域内地表沉降的年变化情况，对获得的沉降值进行叠加，得到以年为间隔的沉降值，然后进行反距离权重插值，得到研究区域内完整的年沉降面积变化情况，如图7所示。对沉降面积进行统计分析，得到各矿区沉降面积变化柱状图，如图8所示。由图8可知：在4 a时间里，夹河矿沉降面积由23.83 km2下降为19.59 km2，抬升面积由1.15 km2增长到5.42 km2，由此推测夹河矿依旧会有抬升的趋势；张小楼矿沉降面积由16.23 km2下降为15.43 km2，抬升面积增长了0.8 km2，整体抬升缓慢，较为稳定；庞庄矿沉降面积由14.67 km2下降为13.22 km2，抬升面积增长了1.45 km2。徐州西部矿区具体下沉与抬升变化情况见表2。
(a) 2016年
(b) 2017年
(c) 2018年
(d) 2019年
图7 各矿区2016—2019年沉降面积变化情况
Fig.7 Changes in the subsidence area of each mine from 2016 to 2019
(a) 夹河矿
(b) 张小楼矿
(c) 庞庄矿
图8 各矿区2016—2019年沉降面积变化柱状图
Fig.8 Histogram of changes in the subsidence area of each mine from 2016 to 2019
表2 徐州西部矿区下沉与抬升变化
Table 2 Changes of subsidence and uplift in the western mine of Xuzhou
形变量/mm面积/km22016年2017年2018年2019年-120～-1000.130.19↓0.050.22-100～-500.720.430.971.05-50～-1023.812.092.74 0.45 -10～030.11↓9.69↓2.79↓2.150～104.745.28↓1.14↓0.4910～500.011.700.270.91
注：↓表示下沉。
4 结论
(1) 基于自2015-12-21—2019-12-24的62景Sentinel-1A影像，利用SBAS InSAR技术获取了徐州西部矿区矿井在关闭以后的地表沉降信息。其中庞庄矿东南部E区域最大沉降速率可达到-48 mm/a，累计下沉178 mm；夹河矿最大沉降速率达到-25.8 mm/a，累计下沉141 mm，张小楼矿最大沉降速率达到-29.6 mm/a，累计下沉133 mm。
(2) 夹河矿观测线最大倾斜与曲率出现在2019-08-14，分别为0.51 mm/m和0.012 4 mm/m2；张小楼矿最大倾斜与曲率分别为0.689 mm/m和-0.029 mm/m2；庞庄矿在2019-12-24最大倾斜达1.70 mm/m，曲率最大值为-0.039 mm/m2，因受拾屯矿开采的影响，仍有不断扩大的趋势。根据“三下采煤规范”，以上形变值均在地表允许变形值范围内，在闭矿后的4 a内未对地表及建筑物造成危害。
(3) 采用反距离权重方法对研究区进行插值统计分析，得出结论：徐州西部矿井关闭后，2016—2019年，整个矿区抬升区域面积由4.95 km2上升到11.29 km2，沉降区域面积由55.04 km2下降到48.7 km2。
[1] 贺佑国，叶旭东，王震.煤炭工业发展形势及“十三五”展望[N].中国能源报,2015-02-02(012).
HE Youguo,YE Xudong,WANG Zhen.The development situation of the coal industry and the prospect of the “13th Five-Year Plan” [N].China Energy News,2015-02-02(012).
[2] 袁亮,姜耀东,王凯,等.我国关闭/废弃矿井资源精准开发利用的科学思考[J].煤炭学报,2018,43(1):14-20.
YUAN Liang,JIANG Yaodong,WANG Kai,et al.Precision exploitation and utilization of closed/abandoned mine resources in China[J].Journal of China Coal Society,2018,43(1):14-20.
[3] 王金华,康红普,刘见中,等.我国绿色煤炭资源开发布局战略研究[J].中国矿业大学学报,2018,47(1):15-20.
WANG Jinhua,KANG Hongpu,LIU Jianzhong,et al.Layout strategic research of green coal resource development in China[J].Journal of China University of Mining & Technology,2018,47(1):15-20.
[4] 冯东梅.闭矿后矿区环境地质灾害分析[J].辽宁工程技术大学学报,2002(4):526-527.
FENG Dongmei.Analysis of environmental geological disasters in mining area after mine closure [J].Journal of Liaoning Technical University,2002(4):526-527.
[5] ROSEN P A,HENSLEY S,JOUGHIN I R,et al.Synthetic aperture radar interferometry[J].Proceedings of the IEEE,2000,88(3):333-382.
[6] SIMONS M,ROSEN P A.Inteferometic symnthetic aperture radar geodesy[J].Treatise on Geophysics,2007(3):391-446.
[7] FERRETTI A,PRATI C,ROCCA F.Permanent scatterers in SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing,2001,39(1):8-20.
[8] USAI S.A least-squares approach for long-term monitoring of deformations with differential SAR interferometry[C]//IEEE International Geoscience & Remote Sensing Symposium,Toronto,Ontario,2002.
[9] BERARDINO P,FORNARO G,LANARI R,et al.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].IEEE Transactions on Geoscience & Remote Sensing,2002,40(11):2375-2383.
[10] MORA O,MALLORQUI J J,BROQUETAS A.Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images[J].IEEE Transactions on Geoence & Remote Sensing,2003,41(10):2243-2253.
[11] CHEN B,LI Z,YU C,et al.Three-dimensional time-varying large surface displacements in coal exploiting areas revealed through integration of SAR pixel offset measurements and mining subsidence model[J].Remote Sensing of Environment,2020,240:1-12.
[12] TANG Fuquan.Mining subsidence monitoring using the method of combining InSAR and GPS technology[J].Journal of Coal Science and Engineering(China),2011,17:133-136.
[13] LI Shengyi,LIU Yingchun,ZHU Rongbo,et al.Use of D-InSAR technique for monitoring ground subsidence in the Yanzhou coal mining area(China)[J].Applied Mechanics and Materials,2010,34-35:756-760.
[14] JIA Xiuming,MA Chao,ZHAO Anyuan.Environmental investigation and evaluation of land subsidence in the Datong coalfield based on InSAR technology[J].Acta Geologica Sinica,2008,82(5):1035-1044.
[15] WALTER D,HOFFMANN J,KAMPES B,et al.Radar interferometric analysis of mining induced surface subsidence using permanent scatterers[C]//Proceedings of the 2004 Envisat & ERS Symposium,Salzburg,2004.
[16] SAMSONOV S,D OREYE N,SMETS B.Ground deformation associated with post-mining activity at the French-German border revealed by novel InSAR time series method[J].International Journal of Applied Earth Observation & Geoinformation,2013,23(8):142-154.
[17] ZHAO C,LU Z,ZHANG Q,et al.Mining collapse monitoring with SAR imagery data:a case study of Datong Mine,China[J].Journal of Applied Remote Sensing,2014,8(1):083574.
[18] ZHENG L,ZHU L,WANG W,et al.Land subsidence related to coal mining in China revealed by L-band InSAR analysis[J].International Journal of Environmental Research and Public Health,2020,17(4):1170-1187.