0 引 言
煤矿开采引发大区域的地面沉降，严重的还会引起地面坍塌、地物严重形变等现象，对采矿区域的居民和环境造成一系列的影响，甚至威胁着居民的人身安全问题[1-2]。InSAR技术是近年来较受追捧的地面沉降监测手段，相比GPS测量、常规大地测量法和近景摄影测量等传统方法，具有监测精度高(可达毫米级)、覆盖范围广、作业效率高、数据处理流程化等优势。目前常用的InSAR技术有D-InSAR、PS-InSAR、SBAS-InSAR、CR-InSAR、CT-InSAR等[3-7]。InSAR技术经过十几年的发展，已经在城市地面沉降监测、山体滑坡监测及矿区沉降监测等方面有了广泛的应用研究[3-5],与D-InSAR技术相比，小基线子集(SBAS)InSAR技术可以有效地克服差分干涉合成孔径雷达(D-InSAR)的空间和时间去相关，并可以获得时间序列的沉降[6-7]。与永久性散射体(PS)InSAR技术相比，SBAS-InSAR技术减少了所需的SAR图像数量，并且精度更高、计算速度更快，因此在地面沉降监测应用中SBAS-InSAR技术具有良好的监测精度和有效性。
国内外相关学者陆续开展了关于InSAR技术的应用研究，文献[8]用D-InSAR技术对法国Vauvert地区的一个盐矿进行了地面沉陷监测，证明了D-InSAR技术对地面沉降监测的可行性，但不足的是时空基线、大气效应等误差影响了精度。而SBAS-InSAR技术可以减弱时空基线、大气效应等误差，通过有限数量的影像可以得到毫米级的监测结果。刘志敏等[9]结合D-InSAR技术和SBAS-InSAR技术对长治矿区进行了试验分析，发现矿区因开采时间、开采方式、采储量以及地形等因素的不同而呈现不同的沉降结果。张艳梅等[10]利用SBAS-InSAR技术的时序地表沉降监测方法对西安市城区及周边进行了监测，对西安城区的沉降范围变化进行了分析。冯文凯等[11]利用SBAS-InSAR技术对金沙江流域的老滑坡进行了分析，取得了较好的试验结果。已有大量学者利用SBAS-InSAR技术对城市地面、长三角洲和滑坡形变等进行了监测实验，且都获得了较好的监测结果，但对矿区工作面的沉降监测相对较少。
笔者以欧空局(ESA)Sentinel-1 A卫星获取的19幅高分辨率SAR影像和精度为30 m分辨率的SRTM DEM作为数据支撑，采用SBAS-InSAR技术对郓城某矿区域进行时序变形监测，通过时序监测数据对沉降漏斗区域内各工作面进行了沉降趋势分析，并与水准测量数据进行验证分析，更加清楚地掌握沉降漏斗区域内工作面的演化规律。
1 SBAS-InSAR基本原理
获取同一研究区域的n幅影像，按时间顺序T0、T1、T2、…、Tn排列，选择适合的时间基线和空间基线，生成m幅差分干涉图，则m应满足
(1)
假设TA时刻获取的影像与TB 时刻获取的影像(TA<TB) 生成第k幅差分干涉图，图上某个像元方位向坐标为x，距离向坐标为 r，那么任意像元(x，r)上TB时刻与TA时刻的差分干涉相位为φA(x,r)、φB(x,r)为未知数，m幅差分干涉图的干涉相位δφk(x,r)(k=1，…，m)则为观测量[12],那么像元干涉相位可以表示为
δφk(x,r)=φB(x,r)-φA(x,r)≈
(2)
式中:λ为雷达信号波长；d(TB,x,r)、d(TA,x,r)为像元(x,r)在TB时刻与TA时刻相对于初始时刻T0在雷达LOS方向的累计形变量。
为便于理解，式(2)并未考虑大气变化引起的相位变化、去相关现象，也未精确去除地形相位。
差分干涉图有m幅，根据式(2)就可以得到m个方程，用矩阵表示方程组为
δφ(x,r)=Aφ(x,r)
(3)
式(3)中：A为m×n矩阵，m对应于干涉图数量，n对应于SAR影像数量，φ(x,r)=[φ1(x,r),…,φn(x,r)]为每一景SAR影像中高相干点对应的相位值所组成的向量，δφ(x,r)=[δφ1(x,r),…,δφm(x,r)]为各差分干涉图对应的解缠相位值组成的向量；当A的秩为n时，利用最小二乘法[13]，解得φ(x,r)为
φ(x,r)=(ATA)-1ATδφ(x,r)
(4)
式(4)中，当A的秩小于n时，利用奇异值分解法(SVD)对方程组进行求解，但是求解得到的形变量一般表现出不连续、跳跃等现象[14]，因此，可以考虑求解2个相邻时刻相位变化平均速率为
(5)
VT=[v1,v2,…,vn]
(6)
根据式(5)得到相位平均速率，最后对各时间段速率进行积分，得到该时间段内的总形变量。
2 矿区概况和数据概述
2.1 矿区概况
选择选用菏泽市郓城县某煤矿作为试验区域。该煤矿位于主要位于田地区域，附近有车楼村、文庄和邵集村等诸多村庄，其地理坐标位于东经115°50′、北纬35°27′左右；该试验区域包含5个开采工作面，如图1所示。
图1 工作面地理分布
Fig.1 Geographical distribution of working faces
工作面编号分别为：2302、1307、1305、1303和4302，工作面的开采厚度范围为2.8～3.4 m,2302工作面走向长837 m，倾向长177 m，开采时间为2016年1月至2016年11月；1307工作面走向长1 456 m，倾向长246 m，开采时间为2016年1月至2017年3月；1305工作面走向长1 629 m，倾向长233 m，开采时间为2015年6月至2016年7月；1303工作面走向长1 602 m，倾向长209 m，开采时间为2012年6月至2015年5月；4302工作面走向长417 m，倾向长399 m，开采时间为2016年11月至2017年12月[15]。
2.2 试验数据
本试验利用19幅干涉宽幅(Interferometric wide swath ,IW)模式、升轨、VV极化方式的Sentinel-1 A SAR卫星影像数据，时间跨度为2016年9月到2017年9月[16]，其空间分辨率为5 m×20 m，其入射角都为38.9°左右[17]，详细影像信息见表1。由于每幅影像数据的幅宽为250 km，为了方便试验处理，所以对原始影像数据进行了裁剪，裁剪范围为:北纬35.57°—35.25°、东经115.53°—116°。差分干涉所用的DEM数据为SRTM1(空间分辨率为30 m),覆盖范围为北纬36°到南纬35°、西经115°到东经116°，覆盖整个煤矿区域[18]。
表1 影像参数
Table 1 Image parameters
序号成像时间时间基线/d垂直基线/m12016-09-27-8425.129822016-10-09-72-23.32333456789101112131415161718192016-10-212016-11-022016-11-142016-11-262016-12-082016-12-202017-01-012017-01-132017-01-252017-02-182017-03-022017-03-262017-04-192017-06-062017-08-052017-08-292017-09-10-60-48-36-24-120+12+24+36+60+72+96+120+168+228+252+26412.135275.105594.743670.98068.2250032.486484.113072.165337.687599.847371.1516-61.2478 14.462841.4459-22.5118 25.3492
3 数据处理和结果分析
3.1 数据处理
试验基于Sentinel-1 A卫星数据、DEM数据和ENVI软件的SARscape，利用SBAS-InSAR技术对郓城某矿进行了沉降反演，主要步骤如下：
1)基线估算及小基线集生成:根据理论计算可得到n(n-1)/2=171个干涉对的空间基线和时间基线。为了提高精度，设置空间基线阈值为临界基线(参考临界基线为5 554 m)的2%，设置时间基线阈值为100 d，最终生成了99个干涉相对。
2)SAR影像配准和小基线集干涉流及相位解缠：选择2016年12月20日的影像为超级主影像，再与其他影像进行配准；小基线集干涉流包括干涉图生成、去平地效应、干涉图滤波(采用Goldstein法)、相干系数计算[19]；相位解缠采用的是最小费用流法(Minimum Cost Flow)，相位解缠相干系数阈值设为0.45。其中干涉处理采用精密轨道星历参数 AUX_POEORB，精度在5 cm以内。
3)轨道精炼和重去平：轨道参数选择AUX_POE-ORB精密星厉，采用3次轨道精炼多项式估算轨道精炼和相位偏移量，消除斜坡相位，通过选择GCP点所有数据进行重去平。
4)估算平均位移速率和DEM纠正：根据式(3)、式(4)建立方程组并求解平均位移速率和DEM纠正，主要过程包括第1次估算平均位移和DEM纠正(通过大气相位消减和估算)、第2次估算平均位移和DEM纠正(时间序列和平均位移速率估算)。
5)地理编码：将上述生成的LOS方向的形变量都投影编码到GCS-WGS-84坐标系下的垂直方向。
3.2 结果分析
3.2.1 沉降速率分析
根据上述处理，获得了该煤矿的年平均沉降速率图(图2)，其中负值代表地面沉降，正值代表地面抬升，其中最大沉降速率为-389 mm/a。
从图2中可以看出一个明显的沉降漏斗，该沉降面积大约为3.13 km2，在监测时间段内,平均沉降速率在-30～-152 mm/a(即绿色区域)的面积为2.51 km2，约占沉降漏斗的80.25%，平均沉降速率在-152～-307 mm/a(即黄色区域)的面积为0.33 km2，约占沉降漏斗的10.42%，平均沉降速率在-307～-389 mm/a(即红色区域)的面积为0.29 km2，约占沉降漏斗的9.33%，可见沉降漏斗区域内，沉降速率大部分发生在绿色、黄色区域，该区域为沉降影响区域；而红色区域为活跃区域，该区域正处于开采或者刚开采完阶段。
图2 年平均沉降速率
Fig.2 Annual average settlement rate
5个工作面都在这个沉降漏斗区域，其中最大沉降速率的区域在1307工作面，其沉降速率为-178～-389 mm/a；2302工作面的沉降速率为-239～-353 mm/a；1305工作面的沉降速率为-89～-345 mm/a；4302工作面的沉降速率为-251～-352 mm/a；1303工作面的沉降速率相对于前4个工作面是最小的，主要原因是其开采结束时间在本试验研究时段的前半年左右，并且该工作面的开采时间相对比较长，它的沉降速率为-90～-273 mm/a，其南部沉降速率大的主要原因是受附近工作面开采的影响。
3.2.2 时序累计沉降分析
以2016年9月27日为初始值的10幅时序累计沉降量如图3所示。从2016年9月到11月期间，2302工作面、1307工作面和1305工作面有明显的沉降，而1303和4302工作面没有明显的沉降现象，与工作面开采的时间吻合；从2016年11月之后，沉降区域逐渐增大，并且2302、1307和4302工作面的沉降量明显增加，1305、1303工作面的南部也有持续的沉降情况，主要原因有2点：①该工作面紧挨1307和4302工作面，这2个工作面在11月一直在开采，导致一定的沉降影响；②该工作面虽然已经开采结束，根据开采沉陷规律，开采深度在300～600 m，地表点移动时间为5年左右，所以有沉降情况是正常的。
图3 时序累计沉降量
Fig.3 Time series cumulative subsidence
每个月的整体沉降量都在增加，其中工作面2302、1307和4302在2017年3月到4月的沉降量有明显变化；直到2017年9月，2302工作面的沉降量为-110～-327 mm，4302工作面的沉降量为-165～-368 mm，这2个工作面几乎都处于沉降漏斗区域，从时序沉降量上来看，沉降情况还会持续一段时间；1307工作面的沉降量为-21～-381 mm，1305工作面的沉降量为5～-297 mm，1303工作面的沉降量为-17～-186 mm，预计1307和1305工作面的沉降趋势为北部的沉降量逐渐缩小，南部沉降量还会逐渐增大；1303工作面并没有达到稳定，沉降还在继续。
从沉降面积来看，发现沉降漏斗的面积一直在扩大，从2016年10月到2017年1月的沉降面积变化为0.4～1.21 km2，2017年1月到2017年3月沉降面积变化为1.21～1.90 km2，2017年3月到2017年9月的沉降面积变化为1.90～3.13 km2。
3.2.3 水准数据对比分析
为了了解工作面的沉降情况，收集到5期2302工作面走向观测线的水准数据，观测日期分别为2016-10-20、2016-12-10 、2017-02-09、2017-05-21、2017-08-20，以2016-10-20为起始时间计算出水准点的累计沉降量，具体水准点分布如图4所示。
图4 观测线布设
Fig.4 Observation line layout
为了与水准数据对比分析，通过最邻近法选取对应的像素点，并且以折线形式绘出，如图5所示。
由图5可看出2302工作面南部沉降相对更活跃，由北到南其沉降量逐渐增加，而南部的沉降量差别不大，其中最大下沉点为A17，其最大下沉量为-372 mm，其最大沉降速率为-352 mm/a，与开采方向是由北向南开采一致，所以形成了勺状沉降漏斗。发现像素点的沉降趋势与水准数据的沉降趋势一致。
图5 2302工作面走向沉降对比
Fig.5 Comparison of dip subsidence in No.2302 working face
为了对比同时期的单点精度，表2选择水准观测线的4期数据(2016-10-20、2016-12-10、2017-02-09、2017-08-20)与像素点数据进行误差分析。由表2可知：单点最大误差为24 mm，最大中误差为9.25 mm，符合DD2014-11《地面沉降干涉雷达数据处理技术规程》中的精度(10 mm) 要求，证明SBAS-InSAR技术沉降监测精度可靠[20]。
表2 观测线时序沉降误差
Table 2 Observation line timing subsidence error
观测日期单点最大误差/mm单点最小误差/mm平均误差/mm中误差/mm2016-10-201315.397.072016-12-10926.135.742017-02-092429.479.252017-08-20815.484.66
5 结 论
1)采用SBAS-InSAR技术得到某煤矿内沉降漏斗的沉降速率，发现2302、1307、1305、1303和4302工作面最大沉降速率分别为353、389、345、273和352 mm/a。
2)通过对各工作面的沉降速率和时序沉降量进行分析，发现1303和1305工作面北部沉降量在减小，南部沉降量还会增加；而4302、2302和1307工作面还处于活跃期，其沉降还会不断增加。
3)通过对比分析2302工作面走向观测线的时序沉降量和水准监测数据，发现最大中误差为9.25 mm，符合DD2014-11《地面沉降干涉雷达数据处理技术规程》中的精度(10 mm) 要求，证明SBAS-InSAR技术沉降监测精度可靠。
参考文献(References):
[1] 梁 涛.利用短基线集InSAR技术监测矿区地表形变[J].测绘通报,2014(S1):82-84．
LIANG Tao.Monitoring the surface deformation of mining area by using short baseset InSAR technique[J].Bulletin of Surveying and Mapping,2014(S1):82-84.
[2] 汤伏全,董龙凯,王宗良,等.基于沉陷对称特征的近水平煤层开采InSAR三维位移反演模型[J].煤炭学报,2019,44(1):210-220.
TANG Fuquan,DONG Longkai,WANG Zongliang,et al.InSAR three-dimensional displacement inversion model of near-horizontal coal seam mining based on subsidence symmetry characteristics [J].Journal of China Coal Society,2019,44(1):210-220.
[3] 鲍金杰,汪云甲.InSAR测量技术在矿区地表沉降监测中的应用[J].煤炭科学技术,2012,40(5):96-99，48.
BAO Jinjie,WANG Yunjia.Application of InSAR measurement technology in surface subsidence monitoring of mining area [J] .Coal Science and Technology,2012,40(5):96-99，48.
[4] 刘晓菲,邓喀中,范洪冬,等.D-InSAR监测老采空区残余变形的试验[J].煤炭学报,2014,39(3):467-472.
LIU Xiaofei,DENG Kazhong,FAN Hongdong,et al.Test of D-InSAR monitoring residual deformation in old goaf [J] .Journal of China Coal Society,2014,39(3):467-472.
[5] 李 强,王继仁,杨庆贺.浅埋深煤层开采沉陷预测方法应用及研究[J].煤炭科学技术,2019,47(5):175-181.
LI Qiang,WANG Jiren,YANG Qinghe.Application and research of prediction methods for mining subsidence in shallow buried deep coal seam [J] .Coal Science and Technology,2019,47(5):175-181.
[6] 刘沂轩,耿智海,杨俊凯,等.基于SBAS技术的概率积分法矿区沉降量提取模型[J].煤炭科学技术,2017,45(2):156-161.
LIU Yixuan,GENG Zhihai,YANG Junkai,et al.Extraction model of mining area settlement based on SBAS technology [J] .Coal Science and Technology,2017,45(2):156-161.
[7] 周 吕,郭际明,李 昕,等.基于SBAS-InSAR的北京地区地表沉降监测与分析[J].大地测量与地球动力学，2016(9):793-797．
ZHOU Lyu,GUO Jiming,LI Wei,et al.Monitoring and analysis of surface subsidence in Beijing based on SBAS-InSAR[J].Journal of Geodesy and Geodynamics,2016(9):793-797.
[8] CARNEC C,DELACOURT C.Three years of mining subsidence monitored by SAR interferometry near Gardanne,France[J].Journal of Applied Geophysics,2000,43(1) :43-54．
[9] 刘志敏,李永生,张景发,等.基于SBAS-InSAR的长治矿区地表形变监测[J].国土资源遥感,2014,26(3):37-42.
LIU Zhimin,LI Yongsheng,ZHANG Jingfa,et al.Surface deformation monitoring in Changzhi mining area based on SBAS-InSAR[J].Remote Sensing for Land and Resources,2014,26(3):37-42.
[10] 张艳梅,王 萍,罗 想,等.利用Sentinel-1数据和SBAS-InSAR技术监测西安地表沉降[J].测绘通报，2017(4):93-97．
ZHANG Yanmei,WANG Ping,LUO Xiang,et al.Monitoring surface subsidence in Xi’an using Sentinel-1 data and SBAS-InSAR technology[J].Bulletin of Surveying and Mapping,2017(4):93-97.
[11] 冯文凯,顿佳伟,张国强,等.基于SBAS-InSAR技术的金沙江流域沃达村巨型老滑坡形变分析[J/OL].工程地质学报:1-9[2020-04-19].https://doi.org/10.13544/j.cnki.jeg.2019-411.
FENG Wenkai,DUN Jiawei,ZHANG Guoqiang,et al.Deformation analysis of giant old landslide in Wuda Village of Jinsha River Basin based on SBAS-InSAR technology [J/OL].Journal of Engineering Geology:1-9 [2020-04-19] .https://doi.org/10.13544/j.cnki.jeg.2019-411.
[12] 朱 猛,董少春,尹宏伟,等.基于SBAS-InSAR方法的苏州地区2007—2010年地表形变时空变化研究[J].地球信息科学学报，2016,18(10):1418-1427．
ZHU Meng,DONG Shaochun,YIN Hongwei,et al.Spatial and temporal variations of surface deformation in Suzhou area from 2007 to 2010 based on SBAS-InSAR method[J].Journal of Geo-information Science,2016,18(10):1418-1427.
[13] 郭乐萍,岳建平,岳 顺.SBAS技术在南京河西地表沉降监测中的应用[J].测绘通报，2017(3):26-28，41.
GUO Leping,YUE Jianping,YUE Shun.Application of SBAS technology in surface settlement monitoring in Hexi,Nanjing[J].Bulletin of Surveying and Mapping,2017(3):26-28,41.
[14] 汪 磊．基于SAR技术的矿区大梯度形变时序监测[D]．徐 州:中国矿业大学，2017．
[15] 李 达,邓喀中,高晓雄,等.基于SBAS-InSAR的矿区地表沉降监测与分析[J].武汉大学学报(信息科学版),2018,43(10):1531-1537.
LI Da,DENG Kzhong,GAO Xiaoxiong,et al.Monitoring and analysis of surface settlement in mining area based on SBAS-InSAR[J].Journal of Wuhan University(Information Science Edition),2018,43(10):1531-1537.
[16] 张晓博,赵学胜,葛大庆.基于Sentinel-TOPS模式Stacking技术监测淮南矿区沉降[J].国土资源遥感，2018,30(4):200-205．
ZHANG Xiaobo,ZHAO Xuesheng,GE Daqing.Monitoring the settlement of Huainan mining area based on Sentinel-TOPS model stacking technology[J].Remote Sensing for Land and Resources,2018,30(4):200-205.
[17] PAOLOB，GIANFRANCOF，RICCARDOL，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):23752383．
[18] 朱建军,李志伟,胡 俊.InSAR变形监测方法与研究进展[J].测绘学报，2017,46(10),1717-1733.
ZHUJianjun,LI Zhiwei,HU Jun.Investigation and research progress of InSAR deformation monitoring[J].Acta Geochimica Sinica,2017,46(10),1717-1733.
[19] 华怡颖,胡晋山,康建荣,等.基于SBAS-DInSAR的大宁矿区多工作面开采地表沉陷规律研究[J].金属矿山,2020(3):177-183.
HUA Yiying,HU Jinshan,KANG Jianrong,et al.Research on ground subsidence law of multi-face mining in Daning mining area based on SBAS-DInSAR [J] .Metal Mine,2020(3):177-183.
[20] 赵喜江,王 斌,张在岩.基于SBAS-InSAR技术的鹤岗矿区地表沉降监测与分析[J].黑龙江科技大学学报，2019,29(1):66-70.
ZHAO Xijiang,WANG Bin,ZHANG Zaiyan.Monitoring and analysis of surface subsidence in Hegang mining area based on SBAS-InSAR technology[J].Journal of Heilongjiang University of Science and Technology,2019，29(1):66-70.
