0 引 言
煤炭资源是我国能源结构的重要组成部分，同时也是国家经济发展的重要物质基础[1]。随着煤炭资源的开采，矿区土地质量下降、水土盐渍化、水土流失等现象逐步加剧，同时地表沉陷对矿区附近的建筑物、农用耕地以及交通道路和铁路等产生了不可逆转的损坏[2]。因此，监测矿区长时序地表形变有助于宏观掌握矿区地表及构筑物损坏情况，为政府决策提供依据，保障人民的生命和财产安全。
由于煤炭资源开采导致的地表形变具有覆盖范围广、长时序性、沉降量级大和瞬时性等特点，因此传统的全站仪、GNSS和水准测量方法在监测矿区地表形变时具有工作量大、时效性低、难以保障测量员的生命安全等难点[3]。
自20个世纪90年代以来，合成孔径雷达技术凭借着全天时、全天候等优势逐步被国内外学者广泛地应用于对地观测。进入21世纪以来，单主影像时序InSAR技术(称PSInSAR)[4]和多主影像时序InSAR技术(称SBAS)[5]相继被提出，基于PSInSAR和SBAS技术的核心思想，国内外学者相继提出IPTA[6]、StaMPS[7]、TCPInSAR[8]等时序InSAR技术，这些时序InSAR技术被广泛应用于城市地表沉降监测[9]、冰川移动[10]、地下水开采[11]、矿区形变监测[12]等领域。由于矿区多分布在非城镇区域，建筑物等高反射率物体较少，上述传统的时序InSAR技术在监测矿区的地表形变时获取到的点目标较少，不能获取完整的形变场。针对这一问题，国外学者提出了融合分布式目标进行处理的时序InSAR方法，如SqueeSAR[13]、CAEInSAR[14]等，将分布式目标纳入PSInSAR处理流程，显著提高了非城镇区域的点目标密度。目前，融合分布式目标的时序InSAR方法在地表形变监测领域的应用较少。
基于矿区地质条件的差异，矿区地表塌陷监测方法主要分为2种：当煤炭资源埋深较小时，工作面开采引发的地表沉陷下沉速率大，短时间内沉降达到较大量级，因此学者大多基于开采沉陷规律和地质采矿参数预计地表沉降[15-17]；当煤炭资源埋深较大时，工作面开采后地表沉陷速率相对较为缓慢，时序InSAR方法无需获取矿区工作面信息，即可监测矿区井田边界内的地表沉陷范围。
江苏省沛县矿区煤炭资源丰富，煤炭资源的埋深为1 000 ～1 100 m，沛县矿区工作面开采引发的地表沉降面积逐年增加，因此监测矿区地表沉降意义重大。笔者采用融合分布式目标和“多主影像策略”的SBAS方法，基于14景Sentinel-1A影像，设置时间和空间基线阈值，将时序SAR影像组成若干干涉对子集，获取了张双楼矿区的时序形变，并将试验结果与普通SBAS方法的监测结果进行对比。
1 沛县矿区区域概况与数据
1.1 矿区域概况
江苏省徐州市沛县是我国沿海地区的主要煤炭生产基地，地处山东和江苏两省交界处，东邻微山湖，铁路公路等交通运输较为发达。张双楼煤矿位于江苏省徐州市沛县北部，经度116°42′—116°55′，纬度34°45′—34°51′，矿区东西走向长约13.8 km，南北走向长约4.7 km，矿区边界内分布着大量的村庄和农田，一条矿区自建铁路自西向东横穿矿区。张双楼矿区自20个世纪60年代开始开采煤炭资源，如今煤炭资源逐渐枯竭，目前矿区布设的工作面的采深约为1 100m[18]，张双楼煤矿的地理位置如图1所示。
图1 矿区区域地理位置示意
Fig.1 Geographical map of mining area
1.2 试验数据
选用了2016年10月到2017年4月时间段内的14景Sentinel-1A卫星影像。Sentinel-1A卫星是欧空局于2014年4月发射的对地观测卫星，卫星采用了极地太阳同步轨道，轨道高度693 km，重访周期为12 d，工作波段为C波段。Sentinel-1A卫星提供了多种极化方式和4种成像模式：条带模式、超宽幅模式、宽幅干涉模式和波模式[19]。试验采用宽幅干涉模式下获取的SAR影像，影像的距离向分辨率约为5 m，方位向分辨率约为20 m，影像幅宽约250 km。试验中去除地形相位的外部DEM数据来自于美国宇航局航天飞机测图任务获取的90 m分辨率的DEM数据。
2 试验方法
2.1 分布式目标法
针对非城镇区域点目标密度较低的现象，众多学者将具有分布面积较广、中等反射率等特点的裸地、荒漠等视为分布式目标，将分布式目标与人工建筑、道路等高反射率物体进行联合处理，达到增加点目标密度，获取完整形变场的目的。
分布式目标的处理主要包括2个方面：同质点目标识别和相位优化。使用快速同质点识别法(FaSHPS)[20]对SAR影像的像元进行同质点识别。假设有N景SAR影像，目标像元j的平均强度
可以表示为
(1)
式中：Ai(j)为目标像元在第i景SAR影像的强度。
以目标像元为中心的一定窗口内且平均强度落入如下区间范围内的像元被视为目标像元的同质点为
(2)
式中：P{·}为概率；L为视数；M为参考像元的平均强度；Z1-α/2为标准正态分位点；α为拒绝概率。
结合公式同质点数目的计算结果，设定同质点的数目阈值，针对同质点数目大于一定阈值的像元，结合像元的同质点，采用主成分分析法对像元进行相位优化为
(3)
式中：C(j)为SAR影像中目标像元的相位序列；λk为特征值；μk为特征向量；H为共轭转置。
当特征值λk越大，其对应的特征向量越占据目标像元回波信号的主导优势。由于试验中对Sentinel-1A卫星数据进行了5∶1多视，影像的分辨率约20 m，分辨率单元内的地物目标并不均匀，因此将最大特征值对应的特征向量作为像元优化后的相位。
2.2 SBAS方法原理
针对干涉图受空间时间失相干影响的现象，Berardino学者于2002年提出多主影像时序InSAR技术—SBAS，SBAS技术通过设置时间空间基线阈值，将SAR影像数据集分若干个干涉对子集，利用最小二乘法将各个数据集进行连接，最终利用奇异值分解法解算出影像时间跨度内的形变信息。
假设分别在(t1,t2,…，tN)时刻获取覆盖研究区域的N景SAR影像，设定基线阈值后组成M组干涉对。假设干涉对i由tX和tY两景SAR影像组成，则干涉对中任意一像元j的干涉相位Δφ可用以下方程表示为
(4)
式中：λ为雷达波长；d(tY,i,j)和d(tX,i,j)分别为相应时间相对初始SAR影像沿视线方向的累计形变量。
差分干涉相位主要由以下几部分组成，即
Δφ=Δφdef+Δφatm+Δφorb+Δφtopo+Δφnoise
(5)
式中：Δφdef为雷达视线向相位；Δφatm为大气延迟相位；Δφorb为轨道误差相位；Δφtopo为地形相位；Δφnoise为噪声相位。
N景SAR影像组成M组干涉对，因此可以列出M组方程，可用矩阵方程式表示为
Δφ=AM×NΦ≈(4π/λ)BM×NV
(6)
式中：AM×N为M×N维度的矩阵，矩阵中每1行对应的主影像元素为1，副影像元素为-1，其余元素均为0；Φ是N×1维度的矩阵，对应元素分别为SAR影像的缠绕相位；BM×N为M×N维度的矩阵，矩阵中对角线元素为干涉对的时间基线；V是N×1维度的矩阵，对应元素分别为干涉对时间内的平均沉降速率。
若矩阵B的秩不小于N时，可用最小二乘法求解形变速率；若矩阵B的秩小于N时，利用奇异值分解法求解形变速率。通过对干涉对时间内的形变速率进行积分获取研究区域视线向的时序形变，同时忽略水平移动的影响，最终将视线向的时序形变转化成垂直向的时序形变。
2.3 试验处理流程
完整技术路线如图2所示。
图2 技术路线流程
Fig.2 Flow chart of technical routes
1)获取SAR数据集对应的强度影像，利用FaSHPS方法判别每一个像元的同质点。
2)结合像元的同质点，利用主成分分析方法进行相位优化。
3)设定时间空间基线，将影像组成若干小基线集，结合外部DEM数据，获取差分干涉图。
4)分离出大气相位和DEM误差相位，获取研究区域的形变信息。
3 试验结果与分析
3.1 试验结果
针对矿区地表植被覆盖较广、干涉图受时间空间失相干影响严重等问题，试验利用融合分布式目标的SBAS方法对2016年10月到2017年4月的14景Sentinel-1A影像，将相邻的两景SAR影像组成干涉对进行时序处理，获取了研究区域的形变信息(表1)，并利用SBAS方法获取的形变信息作为参考，验证试验结果的准确性。
试验如流程图2所示，试验分别利用传统SBAS方法和融合分布式目标的SBAS方法获取的研究区域的形变信息，获取的沉降速率和累计沉降量分别如图3和图4所示。
由图3、图4可知，由于工作面的开采，张双楼煤矿在研究时间段内形成了2个较为明显的圆形下沉盆地，分别位于矿区的中部和右上部位，矿区其余地方无明显沉降区域，形变较为稳定。其中，融合分布式目标的SBAS方法显著地增加了矿区内的点密度，增加的点目标主要集中在农田、道路和建筑等区域，可以明显地辨别出矿区下沉盆地的边界，但由于矿区下沉盆地中心的沉降量级较大、沉降速度快，因此下沉盆地中心内点数量增加的不多。
图4 张双楼煤矿累积沉降
Fig.4 Accumulative subsidence of Zhangshuanglou Mine
图3 张双楼煤矿沉降速率
Fig.3 Subsidence velocity of Zhangshuanglou Mine
研究时间段内，SBAS方法共获取了19 936个点目标，融合分布式目标的SBAS方法共获取了85 632个点目标，融合分布式目标的SBAS方法的点目标数目约是SBAS方法的4.3倍。SBAS方法获取的地表最大沉降速率和最大累积沉降量分别为-301.0 mm/a和-202.4 mm;融合分布式目标的SBAS方法获取的地表最大沉降速率和最大累积沉降量分别为-308.7 mm/a和-209.3 mm。
形变结果表明，融合分布式目标的SBAS方法和SBAS方法的监测结果较为接近。试验将进一步计算2种方法间形变结果的误差，验证融合分布式目标的SBAS方法监测结果的可靠性，并分析下沉盆地的时序变化情况。
3.2 结果验证与分析
试验分别将矿区左侧的下沉盆地和右侧下沉盆地标记为A和B，为了直观地分析2种方法监测结果间的差异，试验分别在A和B区域选取一个共有的点目标，分析了2种方法获取的点目标时序形变之间的差异，同时沿下沉盆地构造剖面线，分析下沉盆地的时序形变情况。
黑色的点为选取的共同点，紫色线条为构造的剖面线如图5所示。试验分别在下沉盆地A和下沉盆地B选取一点，通过分析该点的时序形变情况，证明融合分布式目标的SBAS方法监测结果的可靠性，相同点处，2种方法获取的时序形变如图6所示。
图5 张双楼煤矿下沉盆地累积沉降
Fig.5 Accumulative subsidence of subsidence basins in Zhangshuanglou Mine
图6中紫色线条代表融合分布式目标的SBAS方法监测到的时序形变，棕色线条代表SBAS方法监测到的时序形变。在选取的共同点处2种方法监测到的时序变化趋势较为一致，地表沉降量均随时间的推移在逐渐增大。其中，下沉盆地A处共同点的两种监测结果沉降量最大误差为11.7 mm，平均误差为6.5 mm；下沉盆地B处共同点的2种监测结果沉降量最大误差为7.4 mm，平均误差为5.0 mm。试验结果证明，融合分布式目标的SBAS方法监测结果与SBAS方法监测结果间的误差较小，因此融合分布式目标的SBAS方法的监测结果是可靠的。
图6 张双楼煤矿下沉盆地共同点的时序形变
Fig.6 Time-series subsidence of common point in Zhangshuanglou Mine
为分析矿区地表的时序演变规律，试验分别沿两处下沉盆地构造剖面线，如图5所示。试验选取了距离剖面线5 m范围内的点目标对下沉盆地的形变进行时序分析。
如图7所示，在下沉盆地A和B中，随着到剖面线端点距离的增加，点目标的下沉呈先逐步增加后逐步减小的趋势，与矿区盆地下沉趋势较为符合。由于下沉盆地B区域地表分布着部分水域，因此剖面线后半段点目标的密度较低。
图7 张双楼煤矿下沉盆地的时序形变
Fig.7 Time-series subsidence of subsidence basins in Zhangshuanglou Mine
矿区地表的下沉与工作面的开采进度和矿区的地质条件有较直接的关系，因此矿区地表的沉降并非呈线性规律。根据图7所示，其中2016-10-16—2016-11-09，2016-12-03—2017-01-08和2017-02-01—2017-03-09三个时间段内地表的沉降量较大。研究时间段内，下沉盆地A剖面线上点目标的最大累积沉降量达到了-118.1 mm；下沉盆地B剖面线上点目标的最大累积沉降量达到了-205.1 mm，两处下沉区域的沉降量均有逐步增加的趋势，因此需要对该地区进行长时序监测，从而进行灾害预警。
4 结 论
1)融合分布式目标的SBAS方法监测的点目标密度是SBAS方法的4.3倍，而且2种方法监测结果较为接近，因此融合分布式目标的SBAS方法比SBAS方法更具有优势，为矿区地表形变监测提供了新的手段。
2)监测结果表明，2016年10月4日至2017年4月2日时间段内张双楼煤矿内的最大地表沉降量约为-209.3 mm，且沉降区域的下沉有着逐步增加的趋势，因此需对该矿区进行长时序的监测，提供灾害预警，最大程度减少损失。
3)虽然DSInSAR技术可以获取高密度的地表监测点位，但仍难以提取开采沉陷积水区的地表形变，有必要结合水下地形测量、概率法预测模型等技术方法，实现开采沉陷全盆地信息的监测。
参考文献(References):
[1] 赵 敏,杨伟红,王国平,等.我国煤炭资源的战略储备研究[J].中国矿业,2017,26(10):90-92.
ZHAO Min,YANG Weihong,WANG Guoping,et al.Study on the strategic reserve of coal resources in China[J].China Mining Magazine,2017,26(10):90-92.
[2] 邓喀中,谭志祥,姜 岩,等.变形监测及沉陷工程学[M].徐州:中国矿业大学出版社,2014:10-13.
[3] 杨建图,姜衍祥,周 俊,等.GPS测量地面沉降的可靠性及精度分析[J].大地测量与地球动力学,2006(1):70-75.
YANG Jiantu,JIANG Yanxiang,ZHOU Jun,et al.Analysis reliability and accuracy of subsidence measurement with GPS technique[J].Journal of Geodesy and Geodynamics,2006(1):70-75.
[4] FERRETTI A.Nonlinear subsidence rate estimation using permanent scatters in differential SAR interferometry[J].IEEE Transactions on Geoscience &Remote Sensing,2000,38(5):2202-2212.
[5] 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 and Remote Sensing,2002,40(11):2375-2383.
[6] WERNER C,WEGMULLER U,STROZZI T,et al.Interferometric point target analysis for deformation mapping[C]// IEEE International Geoscience &Remote Sensing Symposium.IEEE,2003.
[7] HOOPER A,ZEBKER H,SEGALL P,et al.A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers[J].Geophysical Research Letters,2004,31(23):1-5.
[8] ZHANG L,DING X,LU Z.Ground settlement monitoring based on temporarily coherent points between two SAR acquisitions[J].Journal of Photogrammetry and Remote Sensing,2011,66(1):146-152.
[9] 王 艳,廖明生,李德仁,等.利用长时间序列相干目标获取地面沉降场[J].地球物理学报,2007,50(2):598-604.
WANG Yan,LIAO Mingsheng,LI Deren,et al.Subsidence velocity retrieval from long-term coherent targets in radar interferometric stacks[J].Chinese Journal of Geophysics,2007,50(2):598-604.
[10] HU J,LI Z,LI J,et al.3-D movement mapping of the alpine glacier in Qinghai-Tibetan Plateau by integrating D-InSAR,MAI and Offset-Tracking:a case study of the Dongkemadi Glacier[J].Global &Planetary Change,2014,118(4):62-68.
[11] QU F,ZHANG Q,LU Z,et al.Land subsidence and ground fissures in Xi\"an,China 2005-2012 revealed by multi-band InSAR time-series analysis[J].Remote Sensing of Environment,2014,155(9):366-376.
[12] FAN HD,XU Q,HU ZB,et al.Using temporarily coherent point interferometric synthetic aperture radar for land subsidence monitoring in a mining region of western China[J].Journal of Applied Remote Sensing,2017,11(2):026003.
[13] FERRETTI A,FUMAGALLI A,NOVALI F,et al.A new algorithm for processing interferometric data-stacks:Squee SAR[J].IEEE Transactions on Geoscience &Remote Sensing,2011,49(9):3460-3470.
[14] FORNARO G,VERDE S,REALE D,et al.Caesar:an approach based on covariance matrix decomposition to improve multi-baseline multi-temporal interferometric SAR Processing[J].IEEE Transactions on Geoscience &Remote Sensing,2014,53(4):2050-2065.
[15] 李 强,王继仁,杨庆贺.浅埋深煤层开采沉陷预测方法应用及研究[J].煤炭科学技术,2019,47(5):175-181.
LI Qiang,WANG Jiren,YANG Qinghe.Application and research on prediction method for mining subsidence in shallow buried deep coal seam[J].Coal Science and Technology,2019,47(5):175-181.
[16] 张 凯,李全生,戴华阳,等.矿区地表移动“空天地”一体化监测技术研究[J].煤炭科学技术,2020,48(2):207-213.
ZHANG Kai,LI Quansheng,DAI Huayang,et al.Research on integrated monitoring technology and practice of “Space-Sky-Ground” on surface movement in mining area [J].Coal Science and Technology,2020,48(2):207-213.
[17] 贾新果.采煤沉陷区地表残余沉陷时间函数模型研究 [J].煤炭科学技术,2018,46(11):157-162.
JIA Xinguo.Study on time function model of residual subsidence at surface coal mining subsidence area [J].Coal Science and Technology,2018,46(11):157-162.
[18] 鹿 璐,范洪冬,郑美楠,等.基于TCP-InSAR的采动区铁路动态沉降监测与时序分析[J].大地测量与地球动力学,2019,39(2):164-168.
LU Lu,FAN Hongdong,ZHENG Meinan,et al.Dynamic settlement monitoring and time series analysis of Railway in mining area based on TCP-InSAR[J].Journal of Geodesy and Geodynamics,2019,39(2):164-168.
[19] 杨 魁,杨建兵,江冰茹.Sentinel-1卫星综述[J].城市勘测,2015,(2):24-27.
YANG Kui,YANG Jianbin,JIANG Bingru.Sentinel-1 satellite overview[J].Urban Geotechnical Investigation &Surveying,2015(2):24-27.
[20] 蒋 弥,丁晓利,何秀凤,等.基于快速分布式目标探测的时序雷达干涉测量方法:以Lost Hills油藏区为例[J].地球物理学报,2016,59(10):3592-3603.
JIANG Mi,DING Xiaoli,HE Xiufeng,et al .FaHPS-InSAR technique for distributed scatterers:a case study over the Lost Hills oil field,California[J].Chinese Journal of Geophysics,2016,59(10):3592-3603.
