郓城矿区地表形变时空演变DS-InSAR监测

王新玲，杜玉玲，蒋金雄，赵 峰，闫世勇

（1.青岛滨海学院 机电工程学院，山东 青岛266555；2.中国矿业大学 自然资源部国土环境与灾害监测重点实验室，江苏 徐州221116）

我国是世界上煤炭资源最丰富的国家之一，随着煤炭资源的持续开采，矿区地面沉降日益严重，导致沉降区内地面损毁、庄稼及道路被淹、桥体开裂等危害发生，给当地人民的生产、生活带来了极大的安全隐患[1-2]。因此，及时掌握矿区地表形变及其动态变化尤为重要。

传统的监测技术如精密水准测量和GNSS等方法耗时耗力、空间采样低，难以满足矿区大范围长时序地表沉降监测及其稳定性评估的需求[3-4]。随着新型空间对地观测技术的不断发展建设，合成孔径差分干涉测量（D-InSAR）技术因其测量精度高、观测范围大、不受天气限制等技术优势，在火山[5]、地震[6]、冰川[7]、滑坡[8]和矿区沉降[9]等领域得到了广泛的应用，但D-InSAR技术易受各类失相干及大气相位延迟等因素影响，不适用于形变速度缓慢、量级较小的地表形变监测研究。时序InSAR分析技术[10-11]则在一定程度上克服了上述不足，为广域小量级形变的监测提供了可能的方法。近年来，为进一步提高时序形变估计质量，融合DS和PS目标的时序InSAR技术受到广泛的关注[12-14]，围绕同质像元的识别[12,15]和相位优化[14]等重要方面开展了相关研究，并逐渐应用于多个领域的地表形变监测[11-19]，但目前针对大采深煤矿引起地表形变的应用较少，针对性研究有所不足。为此在已有研究的基础上，采用基于相干矩阵特征值分解的DS-InSAR技术，开展郓城矿区地面沉降及其时空演变的高精度监测分析，实现对矿区地面采空塌陷区域的稳定性评估，为塌陷区域的恢复治理提供必要的技术与数据支撑。

1 研究区概况与数据源

1.1 矿区概况

郓城矿区位于山东省西南部、菏泽市东北部，面积约624.6 km2，属于温带半湿润季风区海洋-大陆性气候，四季分明，春旱多风，夏热多雨。该地区大地构造属中朝陆台山东台背斜的西南边缘凹陷带，处于“鲁西断块”中的“鲁西南块陷”上，地表全为新生界覆盖，没有基岩出露。郓城地区煤炭矿产资源丰富，主要有李楼煤矿、彭庄煤矿、郭屯煤矿、赵楼煤矿等4个生产矿井，地表覆盖以耕地、城镇工矿和水域水利设施为主，矿区开采沉陷对耕地及交通、水利等基础设施造成了一定的破坏[20]，严重威胁了当地人民的生命财产安全，因此对该区域地表形变进行快速高效监测显得至关重要。郓城矿区位置示意图如图1。

1.2 数据源

选取了2018年1月3日至2019年11月30日期间Sentinel-1卫星获取的55景TOPS模式SAR影像，并以2019年3月11日的SAR数据为主影像进行数据的准确配准。欧洲航空局提供的精密轨道数据用于干涉处理中的轨道精估计和平地效应去除。此外，采用30 m分辨率的SRTM DEM数据进行地理编码和地形相位补偿。

图1 郓城矿区位置示意图Fig.1 Location distribution map of Yuncheng Mining Area

2 研究方法

不同于传统基于高相干点目标的时序InSAR技术，融合分布式目标的DS-InSAR技术能够提取到DS点和PS点目标进行形变解算，从而提升时序InSAR测量中相干目标的空间分布密度，其中同质像元（SHPs）识别[21]和相位优化估计[22]是DS-InSAR的2个核心内容。融合分布目标的时序InSAR技术主要流程如图2。

图2 DS-InSAR数据处理流程图Fig.2 DS-InSAR data processing flowchart

2.1 同质像元识别

同质像元的概念由Ferretti首次提出[23]，同质像元的识别是在假设检验的基础上，识别与样本数据具有相同/相似散射特性的分布式目标像元作为同质像元的集合。目前基于后向散射信息的参数假设检验是同质像元的识别常用方法，其在已知样本分布函数的前提下，判断样本数据的总体分布统计量是否存在显著性差异来进行同质像元的识别，选择双总体T检验进行同质像元的识别。

2.2 相位优化估计

针对分布式目标多视干涉相位，需要采取一定的相位优化算法在满足相位一致性的条件下构建出1组单一主影像优化相位值，获取经自适应多视降噪处理后的1组最佳拟合相位，以尽可能地减弱分布式目标失相干现象的影响。所选用的特征值分解方法是通过从相干矩阵中有效地分离、估计出多元散射机制并得到所对应的相位分量，实现对单一主影像干涉相位的优化估计。为评估DS优化相位的质量，通过拟合比较原相干矩阵中自适应多视干涉相位和优化相位得到的干涉相位，即用时间相干性衡量γDS分布式目标受时间去相干的影响程度，同时也可用于估计相位优化的质量。

式中：γDS为分布式目标的时间相干性；φ^om、φ^on为分布式目标在单一主影像下的优化相位；φmn为自适应多视干涉相位；m、n、o为影像序列号。

3 结果分析

3.1 DS-InSAR监测结果

采用融合分布式目标的DS-InSAR时序分析方法能够在研究区内裸地和农田区域选取到更多的观测点，从而更好地达到表征矿区形变空间分布特征的目的，郓城矿区平均沉降速率图如图3。地面沉降速率等值线图如图4。

图3 郓城矿区平均沉降速率图Fig.3 Average subsidence rate chart of Yuncheng Mining Area

2018—2019年郓城矿区内4个矿井均有沉降发生，其中李楼煤矿和郭屯煤矿沉降较为严重，影响范围最广，产生了1个较大的沉降漏斗，形变中心沉降速率达到了100 mm/a；彭庄煤矿和赵楼煤矿沉降量相对较小，最大沉降速率在70 mm/a。年形变速率对比可知，变形均集中在李楼煤矿和赵屯煤矿，相较于2018年1个大的沉降漏斗，2019年的形变在南北方向上表现出3个形变漏斗，形变中心速率分别可达80、80、70 mm/a。整体监测结果表明，郓城矿区煤矿开采导致地表发生了不同程度的沉降且分布不均。

3.2 重点矿区形变分析

为进一步分析煤炭开采引起的地面沉降，结合矿区信息，绘制了地面沉降等值线图，并针对研究区域内郭屯煤矿和彭庄煤矿2个典型矿井，进行了时间和空间维度的精细化解译分析。

图4 地面沉降速率等值线图Fig.4 Contour map of land subsidence rate

郭屯煤矿北临李楼煤矿、东临彭庄煤矿、南临赵楼煤矿，研究时段内该煤矿的开采工作面为一采区1325工作面和四采区的4303和4302工作面，根据2018年和2019年的沉降速率等值线图可知，2018年矿区中部和北部连接形成1个较大范围的沉降区域，其最大沉降量约为100 mm/a，而在2019年，之前的大区域沉降表现为3块小的沉降漏斗，其中矿区北部与李楼煤矿形成1个完整的沉降漏斗，且沉降中心发生了一定的偏移，最大形变量为80 mm/a；中部地区的形变主要集中在一采区和四采区工作面附近，与实际开采情况相一致；还有一处沉降集中在工业广场附近，最大形变量从2018年的90 mm/a变为2019年的70 mm/a。可见，郭屯煤矿的地表形变主要与煤矿开采相关。

彭庄煤矿位于郓城矿区的东侧，在2018-2019年间，该煤矿对东二采区的3303和3304工作面进行开采，期间形成了4个小的沉降漏斗。由平均沉降速率等值线图可知，2018年该煤矿的沉降主要发生在矿区西部，东部边缘仅有小范围沉降，而2019年东西区域内的沉降区域趋于融合形成1个大的沉降漏斗。监测结果表明，彭庄煤矿左侧3处沉降漏斗分布与先前研究发现的采空塌陷区分布基本一致[24]，可见老采空区依旧有持续沉降现象的发生但相对有减弱趋势，而新的采掘活动也逐步形成了新的沉降漏斗。

3.3 累计沉降量分析

为进一步说明各矿区的形变趋势，在矿区显著沉降区内部选取特征点，（图4中点1～点6），绘制的矿区典型沉降中心时序累积沉降图如图5。

图5 矿区典型沉降中心时序累积沉降图Fig.5 Time series cumulative subscidence diagram of typical subsidence center in mining area

由图5可知，李楼煤矿（点1）、彭庄煤矿（点2）和郭屯煤矿（点3）的形变中心均表现为持续下沉特征，其中，点3和点4处的形变量级和形变趋势趋于一致，且形变量均达到了130 mm；点5处的形变相较于点3和点4较为缓和，而矿区东侧（梁宝寺煤矿，点6）的形变中心在2018年虽然有一定波动但总体保持稳定，在2019年则呈现出下沉趋势。资料显示李楼煤矿和郭屯煤矿的生产能力均为240万t/a，彭庄煤矿的生产能力为110万t/a，可见，点3、点4与点5处的不同形变速率可能与开采量有关。而点6处沉降中心位于梁宝寺煤矿，开采塌陷未沉稳[25]，且根据该地的夏季炎热多雨气候特点，沉陷区形变监测结果会在夏季产生较大的形变波动。

4结语

1）基于DS-InSAR技术和Sentinel-1雷达影像获取了郓城矿区2018—2019年的地表沉降分布，并对典型区域的地表形变及其时空演变规律进行了重点分析，结果表明郓城矿区存在多个明显的沉降中心，且最大形变可达145 mm，郓城矿区形变主要体现在李楼煤矿和郭屯煤矿，彭庄煤矿和赵楼煤矿形变量级相对较小。

2）与传统的时序InSAR方法相比，DS-InSAR技术通过同质像元识别、相位优化等步骤能够在裸地和稀疏低矮植被等区域有效识别出分布式目标，提高了测量点数量和空间分布密度，从而准确、全面地反映了矿区开采沉降的影响范围和形变区域的空间分布特征，为矿区稳定性分析和土地合理利用等提供了良好的技术和数据支撑，有助于预防和减少矿山建设及生产活动造成的矿山地质环境及地质灾害。