徐州矿区2016—2018年地面沉降监测与分析

李梦梦， 范雪婷， 陈 超， 李倩楠， 杨 锦

(江苏省基础地理信息中心，南京 210013)

0 引言

随着科学的进步和经济的发展，地面沉降问题日益突出，已成为21世纪的主要地质灾害之一。地面沉降可破坏基础设施，改变土壤环境和水文环境，引起沉降附近社区的居民恐慌，影响其未来居住规划等[1-2]。造成地面沉降的原因错综复杂，其中固体矿产的开采是影响地面沉降的原因之一[3-4]。被誉为“华东地区的煤炭海洋”的徐州，煤炭资源丰富，但在煤炭开采的同时也造成了地面沉降、地面塌陷等地质环境问题，严重制约了该地区的可持续发展[5-6]。传统的地面沉降监测主要以水准测量和全球卫星导航定位系统(global navigation satellite system，GNSS)测量为主，这些方法虽然单点精度高，但是外业周期长、连续性差、成本较高，且不能获得整个监测区域的整体沉降趋势[7]。

矿区本身的地质条件比较特殊，且现场测量安全隐患大，水准作业难度较高，传统地面沉降监测方法对不定期发生的沉降和塌陷难以做到实时监测[8]，但时间序列合成孔径雷达干涉测量(interferometric synthetic aperture Radar，InSAR)技术可以安全高效地获取大范围、高精度地面沉降数据[9-10]，因此基于InSAR技术来监测矿区地面沉降是较佳选择。

Liu等[5]基于Sentinel-1影像，监测徐州沛县部分矿区地面沉降状况，监测结果表明： ZSL煤矿和SHJ煤矿沉降值与时间呈线性关系，并基于Landsat8影像的分类结果，研究沛县煤矿的沉降与土地利用分类的关系，发现农村地区比城市地区沉降严重，但影像的分类结果还需进一步提高，以便更精确地研究沉降与土地利用分类的关系； 李达等[11]基于小基线集(small baseline subset，SBAS)技术对13景TerraSAR-X影像进行时间序列处理，并基于合成孔径雷达差分干涉测量(differential InSAR，D-InSAR)技术，验证SBAS技术得到的监测结果，发现二者监测结果相近，但SBAS技术得到的监测结果更精确，由于缺乏研究区域的实测数据，监测结果的纵向正确性还有待验证； Chatterjee等[12]基于6对多频C和L波段的SAR影像，通过D-InSAR技术研究矿区的地面沉降，监测结果表明，C波段可以监测到缓慢或中度沉降的矿区，L波段可以监测到快速沉降的矿区，但影像数量较少，基线长度相对较大； Grzovic等[13]基于时间序列永久散射体干涉测量技术(permanent scatter interferometry，PSI)和SBAS技术处理ERS-1/2和ALOS PALSAR数据，得到Springfield地区煤矿开采区的地面沉降监测结果，由于ALOS PALSAR影像数量较少，其监测结果有失真现象； Bateson等[14]基于间断相干小基线集(intermittent small baseline subset，ISBAS)技术处理ERS-1/2数据，提高了监测点密度，监测结果表明，矿区存在地面抬升现象，抬升速率达1 cm/a，其原因可能是相较于开采期，停采期矿区水位恢复，但缺乏矿区开采数据，无法准确确定开采初期的矿井水位，其结果分析只能根据假定的开采数据来进行； Gupta等[15]利用InSAR技术获取研究区域数字高程模型(digital elevation model,DEM)数据来分析矿区地面高程值变化，高程精度偏差为23.07 m，均方根误差为±2.31 m； 陆燕燕等[16]利用ALOS PALSAR影像对徐州张双楼煤矿进行地表形变监测，得到2011-01-16—2011-03-03的地表形变分布图，监测结果表明张双楼煤矿出现3处沉降漏斗，最大沉降量达420 mm，且沉降漏斗与矿区分布一致，由于实验数据比较少，只能采用D-InSAR方法来监测矿区沉降变化，并未考虑时空去相关以及大气非均匀性影响，使得地表形变的监测精度受到限制。

基于Sentinel-1A数据，对在采矿区与停产矿区长时间序列沉降状况监测与对比分析的研究较少，结合地理国情监测数据研究高相干点分布的研究也比较少，本文利用Sentinel-1A影像数据基于多主影像相干目标小基线InSAR方法(multiple master-image coherent target small-baseline interferometric SAR，MCTSB-InSAR)，结合江苏省2016年地理国情监测数据，开展徐州部分在采矿区与停采矿区2016—2018年的地面沉降监测研究与对比分析。

1 研究区概况与实验数据

1.1 研究区概况

徐州市地处苏鲁豫皖四省接壤地区，属温带季风气候，四季分明，夏无酷暑，冬无严寒，是资源富集且组合条件优越的地区之一。其中煤炭储量约69亿t，是江苏省唯一的煤炭产地，也是江苏省重要的能源基地。本文根据矿区停采时间，共选择6个矿区作为研究对象，矿区基本信息见表1，影像、水准点和矿区分布见图1，图中白色框为矿区放大影像。

表1 矿区基本信息Tab.1 Basic information of mining area

图1 影像、水准点和矿区分布Fig.1 Images, bench mark and mining map

1.2 实验数据

本文采用的实验数据是Sentinel-1A雷达卫星影像。Sentinel-1是由欧洲委员会和欧洲航天局针对哥白尼全球对地观测项目研制的一个载有C波段合成孔径雷达的双星星座。Sentinel-1的A星于2014年4月3日发射成功，B星于2016年4月25日发射成功，单星重访周期为12 d，两星同时工作其重访周期为6 d。Sentinel-1A卫星有4种工作模式： 条带模式(strip map，SM)、干涉宽幅模式(interferometric wide swath，IW)、波模式(wave，WV)和超宽幅模式(extra wide swath，EW)，本文使用IW模式的数据。IW模式幅宽为250 km，空间分辨率为5×20 m，采用了TOPS(terrain observation with progressive scans SAR)成像技术。TOPS技术不仅能在保证空间分辨率的情况下增大地面覆盖范围，还解决了ScanSAR的图像不均匀问题，但是该成像模式使影像在方位向的多个脉冲序列(burst)之间易产生多普勒中心不一致问题[17]，引入的干涉相位偏差可表示为：

Δφ=2πfηcΔη

(1)

式中： Δφ为干涉相位偏差；fηc为burst多普勒中心频率变化值； Δη为方位向主辅影像配准偏移量。干涉测量通常要求干涉相位偏差Δφ不超过3°[18]，则影像配准精度需达到0.001个像素[18-19]。

本研究共使用58景研究区范围内的Sentinel-1A升轨数据，影像获取时间跨度为2015年11月—2018年12月。

2 数据处理

基于GDEMSI5.0软件采用MCTSB-InSAR技术实现对徐州矿区的时间序列监测。MCTSB-InSAR融合了永久散射体InSAR技术(permanent scatterer InSAR，PS-InSAR)、SBAS和相干目标等方法的优势，针对Sentinel-1A数据的TOPS成像模式特性，数据处理主要包括： 数据预处理、线性形变反演和非线性形变反演3个部分。

2.1 数据预处理

SAR数据预处理主要是SAR影像差分干涉处理的过程。针对Sentinel-1A影像特点，预处理主要包括轨道加密、精确配准、干涉相对选择、生成差分干涉图等，技术路线如图2所示。

图2 技术路线图Fig.2 Technique flow chart

增强谱分集配准方法虽然可以达到很高的配准精度，但是有2点缺陷： ①若初始影像之间的配准精度低于0.05个像元，则该方法容易出现迭代死循环的现象，无法满足配准要求； ②若2幅影像时间基线较长，容易造成时间失相干。又因为基于精密轨道和DEM数据的配准不受影像相干性影响，可与基于相关性配准互补，使粗配准精度达0.05个像元。因此本文采用基于精密轨道数据和DEM数据、相关性配准和递进式增强谱分集的多级配准方法。其中递进式增强谱分集就是选择与其相近的，能够与第一主影像配准的影像作为第二主影像，配准后再与第一主影像配准，避免时间失相干[20]。

生成差分干涉图分3步，首先基于典型小区域(一般大小为500像素×500像素)快速生成全部(两两组合)小干涉图； 然后根据小干涉图质量确定初选干涉像对； 最后根据初选像对生成大区域干涉图并进行筛选。针对整景影像，先快速生成1 653个小区域差分干涉图，然后对所有小干涉图的相干性进行定量评估，按平均相干系数由高到低，筛选出影像数3～5倍的干涉相对，最终保留315幅差分干涉图，基线对组合如图3所示。

图3 小基线对基线时空分布Fig.3 Spatial and temporal distribution of small baseline pair

2.2 线性形变反演

线性形变反演主要包括： 点目标提取、线性形变速率和高程误差求解等处理。

本文使用幅度离差、平均幅度、平均相干系数“三阈值”法来提取高相干点目标，分别设为0.35，1.1，0.54，本实验整景影像共提取高相干点目标3 046 739个。

对高相干点目标进行三角网连接，三角网边上两顶点之间的二次差分相位差包括线性形变、非线性形变、高程误差、大气影响、噪声5个部分贡献，具体公式为：

(2)

Δω=ΔRnon-liner+ΔφAPS+Δφnoise

(3)

式中：B⊥和T分别为干涉像对的空间基线和时间基线；λ为波长；R为传感器到目标的距离；θ为雷达入射角； Δv为线性形变速率； Δh为高程误差； ΔRnon-liner为非线性形变相位； ΔφAPS为大气延迟相位； Δφnoise为噪声相位； Δω为残余相位。通常|Δω|<π，可通过整体相干性模型来估算Δv和Δh，即

(4)

式中：γ为相干因子值；M为干涉图数量；i为第i幅干涉图；j表示单位复数。利用二维周期图[21]或者空间搜索[22-23]等优化方法，通过模型解算出Δh和Δv。

通常选择|γ|≥0.7的边作为可靠的连接，剔除不连接的高相干点目标。最后从某一参考点开始对连接边上解算出的速率和高程误差增量进行积分，将得到各相干点上形变速率和高程误差的绝对量[24]。

2.3 非线性形变反演

从差分干涉相位中去除线性形变相位和高程误差相位，得到残余相位，包括大气相位、非线性形变相位和噪声相位。根据残余相位3个分量的不同时空频谱特征，利用时间和空间滤波的方法将三者分离出来[21]。最后叠加线性形变信息和非线性形变信息得到累计地面沉降量。

垂直方向的沉降结果可根据雷达视线方向沉降结果得到，公式为：

(5)

式中： Δh⊥为垂直方向沉降结果； Δr为雷达方向沉降结果。基于GDEMSI5.0软件可直接得到垂直方向的形变结果。

3 精度评定

为验证MCTSB-InSAR方法得到的监测结果，收集江苏省“十二五”和“十三五”2套水准数据，比较实测水准点得到的年均沉降速率与邻近InSAR点目标监测值的差异来评估InSAR沉降监测的精度。本文数据处理基于江苏省域数据处理，监测结果由覆盖江苏全省的多个Frame数据拼接得到，数据验证基于最后的全省拼接结果进行。根据80 m邻近原则，全省共筛选出35个实测水准点用来对监测结果进行精度评定，即： 2期水准数据的高程差值比上日期间隔，得到基于水准数据的点位沉降速率，然后与该水准点80 m范围内最近的InSAR点目标的沉降速率对比分析，如图4所示，年均沉降速率二者差值最大为8.2 mm/a，最小为0.1 mm/a，均方根误差为4.0 mm/a，差值绝对值小于1 mm/a的水准点占比为40%，满足规范[25]的要求，保证监测结果的可靠性。

(a) 水准点与InSAR点沉降速率差值统计 (b) 不同差值范围内点数统计图4 水准点与InSAR点沉降速率差值及不同差值范围内点数统计Fig.4 Statistics of subsidence rate difference between level point and InSAR point, and point for different difference ranges

4 结果与分析

4.1 矿区地面沉降速率与累计沉降量分析

监测结果显示6个矿区都有沉降发生，各个矿区的最大沉降速率和最大累计沉降量如图5所示，各个矿区的地面沉降速率和累计地面沉降量空间分布如图6所示。其中权台煤矿、张集煤矿和拾屯煤矿的最大沉降速率均介于20～30 mm/a之间，旗山煤矿略大于30 mm/a，三河尖煤矿和张双楼煤矿均大于100 mm/a。权台煤矿、张集煤矿和拾屯煤矿的最大累计沉降量均不超过100 mm； 旗山煤矿略大于100 mm，三河尖煤矿和张双楼煤矿均大于300 mm。

(a) 矿区最大沉降速率 (b) 矿区最大累计沉降量图5 矿区最大沉降速率及最大累计沉降量Fig.5 Maximum subsidence rate and value of mining area

(a) 权台煤矿地面沉降速率 (b) 权台煤矿累计地面沉降量图6-1 矿区地面沉降速率和累计地面沉降量空间分布Fig.6-1 Spatial distribution of land subsidence rate and accumulated land subsidence in mining area

(c) 旗山煤矿地面沉降速率 (d) 旗山煤矿累计地面沉降量

(e) 张集煤矿地面沉降速率 (f) 张集煤矿累计地面沉降量

(g) 拾屯煤矿地面沉降速率 (h) 拾屯煤矿累计地面沉降量

(i) 三河尖煤矿地面沉降速率 (j) 三河尖煤矿累计地面沉降量图6-2 矿区地面沉降速率和累计地面沉降量空间分布Fig.6-2 Spatial distribution of land subsidence rate and accumulated land subsidence in mining area

(k) 张双楼煤矿地面沉降速率(l) 张双楼煤矿累计地面沉降量图6-3 矿区地面沉降速率和累计地面沉降量空间分布Fig.6-3 Spatial distribution of land subsidence rate and accumulated land subsidence in mining area

统计各个矿区内的高相干点，如图7所示，结果表明： 权台煤矿沉降速率大于0 mm/a的点数占比为68.47%。张集煤矿、拾屯煤矿和旗山煤矿沉降速率介于0～10 mm/a的点数占比最大，分别为： 93.90%，69.48%和79.50%。三河尖煤矿沉降速率介于10～30 mm/a的点数占比为92.18%，张双楼煤矿沉降速率介于30～50 mm/a的点数占比为40.15%，介于10～30 mm/a的点数占比为39.68%。结合江苏省2016年地理国情监测数据，把地类分为房屋及道路、水体、植被和裸露地表4类，如图8所示，由于地理国情数据的原因，三河尖煤矿覆盖区域不完全。从图8中可看出，三河尖煤矿内房屋及道路面积明显大于张双楼煤矿，三河尖煤矿有2 844个高相干点落入房屋及道路内，占该矿区总高相干点数的73.66%，张双楼煤矿有672个高相干点落入房屋及道路内，占该矿区总高相干点数的63.33%。由此也说明，高相干点目标提取较为可靠。

图7 矿区内不同沉降速率范围高相干点个数统计Fig.7 Statistic of high-coherence points in different subsidence rate ranges in the mining area

(a) 三河尖煤矿(b) 张双楼煤矿图8 矿区土地利用分类及高相干点分布Fig.8 Land use classification and high-coherence distribution map of mining area

4.2 矿区地表时序沉降空间分析

以2015年11月27日为起始时间，其他时间相对起始时间的研究矿区时序累计沉降量空间分布如图9所示。从图9中可以看出，从起始时间到2016年6月期间，张双楼煤矿、三河尖煤矿、旗山煤矿、拾屯煤矿有10 mm以上的沉降量，张集煤矿、权台煤矿没有发生明显沉降现象，直到2017年6月，权台煤矿开始出现明显沉降现象，张集煤矿稍有沉降现象。随着时间推移，张双楼煤矿和三河尖煤矿的沉降范围迅速扩大，沉降逐渐加深，对周边影响也比较严重。其余煤矿的沉降均发生在矿区范围内，且沉降范围较稳定，没有明显扩张趋势，对其周边地区影响较小。总体来说，矿区在停采一段时间后沉降仍然会继续发生，不同停采矿区的沉降情况不同，但是停采矿区的沉降情况比在采矿区要轻微得多。

(a) 2015年11月27日

(b) 2016年6月6日

(c) 2016年12月27日

(d) 2017年6月25日图9-1 矿区累计沉降量时序变化Fig.9-1 Time-series variation diagram of accumulated subsidence in mining area

(e) 2017年12月22日

(f) 2018年7月2日

(g) 2018年12月5日图9-2 矿区累计沉降量时序变化Fig.9-2 Time-series variation diagram of accumulated subsidence in mining area

4.3 矿区地表沉降时间特征分析

为分析研究区域的沉降时间特征，提取各矿区最大沉降点的时序沉降值，利用简单的线性拟合模型，得到沉降值和时间的关系，如图10所示。从图10中可以看出，矿区的沉降量与时间基本都呈线性关系，R2均大于0.9。线性拟合斜率绝对值从大到小为： 张双楼煤矿(0.365)、三河尖煤矿(0.277)、旗山煤矿(0.097)、拾屯煤矿(0.065)、权台煤矿(0.057)、张集煤矿(0.056)，张双楼煤矿和三河尖煤矿的线性规律比其他煤矿要强，权台煤矿的线性规律较弱。对权台煤矿的时序沉降量做简单的非线性拟合，发现其沉降规律更符合非线性变化，在监测末期，沉降接近平稳状态。虽然旗山煤矿的停采时间与权台煤矿较一致，但旗山煤矿的时序沉降量没有表现出明显的非线性变化，由于影响矿区沉降的因素多种多样，停采矿区的沉降规律还有待进一步研究。

(a) 权台煤矿(b) 旗山煤矿(c) 张集煤矿

(d) 拾屯煤矿(e) 三河尖煤矿(f) 张双楼煤矿图10 矿区累计沉降量变化统计Fig.10 Statistical chart of accumulative subsidence in mining are

5 结论

本文基于Sentinel-1A雷达影像数据利用MCTSB-InSAR方法监测徐州部分矿区2016—2018年间地面沉降状况，经实测水准数据验证，二者差值的均方根误差为4.0 mm/a，符合技术规范要求，监测结果较准确、可靠。

监测结果表明，2016—2018年间张双楼煤矿沉降最为严重，其次为三河尖煤矿，二者的最大年均沉降速率均超过100 mm/a，最大累计沉降量均超过300 mm。结合江苏省2016年的地理国情监测数据发现，本研究提取的高相干点目标较大部分对应房屋、道路等人工地物，揭示了人工地物具有高相干特性。通过分析沉降量的时序特征发现，除权台煤矿外，其他矿区的沉降有明显线性规律，在采矿区的线性规律比停采矿区的线性规律要强一些，权台煤矿呈非线性沉降，具体沉降规律还有待进一步研究。

研究结果表明，针对矿区特殊的地质条件，基于Sentinel-1A数据和MCTSB-InSAR方法，在安全、高效监测矿区沉降方面有较好的应用前景，其监测结果可为生态环境治理、矿区沉降预测、安全隐患区监测等提供数据支撑和技术服务。