邸 明 慧，张 丽 云，贺 月 娟，徐 宁，乔 良*，李 晓 婧

(1.河北省科学院地理科学研究所/河北省地理信息开发应用工程技术研究中心，河北 石家庄 050011；2.河北交通职业技术学院，河北 石家庄 050035)

0 引言

地下矿产资源开采易导致地面塌陷，不仅严重破坏土地(特别是耕地)和地面建筑物，还会影响交通运输和水资源等[1],因此，矿区地表形变监测对于矿区安全生产、地质环境治理和灾害防控等具有重要作用[2]。InSAR(Interferometric Synthetic Apeture Radar)技术具有不受气候制约和大范围、高精度的监测优势，已被广泛应用于各种地质灾害的监测研究中[3，4]，尤其是SBAS-InSAR(Small Basedline Subsets InSAR)技术克服了传统测量方法存在的时(空)间失相干和大气效应限制，适用于长时序缓慢非线性形变监测[5]，并且在采空塌陷区应用较好。

太行山地区是河北省重要的矿产资源分布区之一，其北段为Fe-Mo、Cu-Pb、Zn-Au成矿区，南段为Fe-煤-铝土矿-石膏成矿区,矿产资源的大量开采导致该区出现严重的地质灾害和生态环境问题。武安地区是河北太行山南段矿产资源集中分布区，因采矿活动造成多处不同程度的地面塌陷[6]，尤其在中部丘陵和盆地区域地面塌陷更普遍[7]。张振生等[8]基于D-InSAR技术发现1993-2004年武安矿山存在11处以煤矿为主的塌陷区，沉降范围由小变大再趋于平稳；张景发等[6，9]利用InSAR技术提取武安市慧兰村矿区塌陷区，发现塌陷区范围和沉降量均随开挖时间推移而变化；周洪月等[10，11]基于时序InSAR处理发现武安市存在一个沉降中心。但以上研究均未针对武安市域内的塌陷区开展区域性监测和分析。1996-2014年武安市大量小矿被关停整合，2009年启动沉陷区综合治理工程，2017年推出矿区生态修复政策，许多矿区得到整治。在此背景下，停采和修复矿区的地面塌陷是否存在继续扩大以及是否有新的塌陷区形成，均直接关系到地方政府相关政策的制定，亟须从全区出发开展地面塌陷区的形变监测研究。因此，本文以武安市为研究区，利用2016年10月-2017年3月和2017年10月-2018年3月两时段各15景降轨Sentinel-1A卫星影像，通过SBAS-InSAR技术反演得到研究区雷达视线方向(Line Of Sight，LOS)的时序累积形变量，结合塌陷区和地裂缝灾害点数据对反演结果进行验证，分析塌陷区的空间分布和形变特征，为地方政府防灾减灾提供科学依据。

1 数据与方法

1.1 研究区与数据

武安市位于河北省邯郸市西北(113°45′～114°22′E，36°28′～37°1′N)，处于太行山隆起向华北平原沉降带的过渡区，地势西北高(最高约为1 900 m)、东南低，连续下降地形被紫山—鼓山打断，形成宽缓的武安盆地(图1)。武安市是全国58个重点产煤县(市)和全国四大富铁矿基地之一，“邯邢式铁矿”因铁矿资源在邯郸—邢台一带集中分布而得名。邯邢式铁矿集中区可以划分为“三带两盆”，其中“三带”自西向东依次为符山成矿带、固镇—武安—矿山—綦村成矿带和鼓山—洪山—新城成矿带，“两盆”指“三带”间的阳邑盆地和武安盆地[12]。在《河北省地质灾害防治“十三五”规划》中，武安—沙河铁、煤矿区属于地面塌陷高易发区。本文将研究区划分为M1、M2和M3共3个区域(图1)，以便进行更细尺度的分析。

图1 研究区概况Fig.1 Survey of the study area

为减少地表散射机制改变对差分干涉的影响[13]，本研究选取2016年10月-2017年3月和2017年10月-2018年3月各15景降轨Sentinel-1A卫星影像(https://search.asf.alaska.edu/)作为数据源，以30 m分辨率的SRTM1(http://www.gscloud.cn/)地形数据作为外部参考DEM；利用武安市已知灾害点(表1)对反演结果进行验证，灾害点数据来源于《邯郸市2016-2019年地质灾害防治方案的通知》[14]，包括14处塌陷区和3处地裂缝，其中中型规模11处，小型规模6处，均处于欠稳定状态。

表1 研究区塌陷区与地裂缝灾害点和反演结果验证Table 1 Subsidence areas and ground fissures in the study area and verification of inversion results

1.2 SBAS-InSAR技术

相比D-InSAR技术，SBAS-InSAR技术可获取时间序列的干涉结果并削弱大气的影响；相比PS-InSAR技术，其数据量要求少，从而在一定程度上可实现同一季节内地物反射条件基本不变，保证像对干涉质量。该技术实现流程为：假设共有N+1幅SAR影像，根据一定的时(空)间基线选取原则，选择影像组合，生成M((N+1)/2≤M≤N(N+1)/2)景短基线距的干涉图；基于数据基线长度和研究区实际情况，对像对进行短基线干涉处理，再利用外部DEM模拟地形相位去除地形相位的干扰；依据干涉图、相位离散度和相位标准差，选择高相干的点进行相位解缠和残余地形相位估计，从而建立分时段的平均变形速率、高程误差和差分相位的模型方程组；最后利用奇异值分解法(SVD)计算最小二乘解，估计时序非线性形变和DEM误差，利用高斯滤波器进行空间滤波和时间滤波以去除大气相位。

2 结果分析

2.1 整体空间分布特征

本研究应用ENVI软件的SARscape模块，对影像依次进行裁剪、配准、干涉(最大空间基线阈值为5%、最大时间基线阈值为120 d)、滤波、解缠、轨道精炼、形变估计和去大气效应处理，生成两时段的时序形变分布图，结合Google Earth影像去除非塌陷区，得到研究区主要塌陷区LOS向累积形变量空间分布(图2)。结果显示，2016年10月-2017年3月塌陷区最大累积形变量为197.84 mm，平均累积形变量为16.44 mm；2017年10月-2018年3月塌陷区最大累积形变量为172.16 mm，平均累积形变量为14.22 mm，均小于前者，说明研究区塌陷有向好发展趋势。塌陷区的总体空间分布趋势基本一致，主要分布在研究区东部，在西部地区分布比较零散且规模较小，但塌陷区主要形变中心位置和规模都随时间有所变化。

图2 研究区LOS向累积形变量空间分布Fig.2 Spatial distribution of the line of sight accumulative deformation in the study area

利用ArcGIS将17处灾害点(表1)与累积形变量空间分布图叠加(图3)对比，发现12处灾害点分布在反演结果塌陷区内，聚隆煤矿周围地裂缝、丰里村及周围地裂缝和西万年村塌陷区位于反演塌陷区的影响范围内，河西村西山坡塌陷区和龙洞山塌陷区两处未显示形变。由此可见，本文反演结果较准确地反映了区内灾害点的分布。

图3 研究区塌陷区和地裂缝分布Fig.3 Distribution of subsidence areas and ground fissures in the study area

结合研究区主要铁、煤矿分布带可以看出，本研究反演出的塌陷区主要位于铁矿分布带(固镇—武安—矿山—綦村成矿带)和煤矿分布带(鼓山—洪山—新城成矿带)内，与武安地区铁、煤矿资源分布范围高度一致，这也印证了本文反演结果的准确性。固镇—武安—矿山—綦村成矿带下伏基岩地层倾向东南，地表高度从西往东虽然逐渐变小，但盆地内的松散沉积物和岩层厚度自西向东递增[15]，这可能导致西侧矿产资源分布深度比东侧浅，发现和开采时间更早。因此该成矿带西侧塌陷区规模比东侧小，推测西侧矿区多为老矿区或弃采矿区，开采量小或矿区恢复使塌陷区累积形变量逐渐变小；而成矿带东侧为新矿区，正处于开采阶段，塌陷区累积形变量相对较大。

2.2 局部空间特征

由图4和表2可知：在M1分区中，塌陷区1、3、4、9在两时段形变速率变化不大，形变中心转移方向均不相同，塌陷区1向北，塌陷区3向西南，塌陷区4、9向东；塌陷区2、6、8在后一时段形变速率减小，其中塌陷区2和8的规模和累积形变量显著减小，平均形变速率分别由7.62 mm/d、8.61 mm/d减小为2.57 mm/d、1.34 mm/d，塌陷区6规模变化不大，但累积形变量明显变小；塌陷区5、7在后一时段规模变化不大，但形变速率增大，增幅分别为5.67 mm/d和3.19 mm/d，形变中心发生转移。

表2 两时段不同分区塌陷区平均形变速率Table 2 Average deformation rates of different subsidence areas in three sub-regions in two periods mm/d

图4 M1、M2和M3分区LOS向累积形变量Fig.4 Maps of the line of sight accumulative deformation in M1，M2 and M3 sub-regions

在M2分区中，塌陷区6、9在两个时段平均形变速率相近，规模变化不大，但形变中心均向南转移；塌陷区1、4、5、8在后一时段平均形变速率减小，表明开采强度有所减弱，其中塌陷区4形变速率最小，形变中心明显消失，塌陷区1和8范围有所增大，塌陷区5形变速率明显减小；塌陷区2、3、7在后一时段平均形变速率增大，其中塌陷区2范围增大，塌陷区3范围减小且形变中心南移，塌陷区7规模和位置变化不大。

在M3分区中，塌陷区2和5在后一时段形变速率降幅较大。塌陷区7为北洺河铁矿塌陷区，是研究区比较典型的老塌陷区，位于武安市上团城村东北约1 km处，矿区面积20.95万m2，2002年建成投产，2003年2月出现塌陷坑，2004年2月主要采空区上部约1 000 m范围内整体下沉2 m左右，塌陷坑周边地表裂缝增多，最宽超过400 mm，地表倾斜明显，在北洺河(改道)北岸形成1.24 km2的塌陷区[16]。北洺河铁矿塌陷区开挖近20年，但本研究显示该塌陷区目前依然存在形变，两时段的平均形变速率分别为2.60 mm/d和3.47 mm/d。由此可见矿产开挖所形成的塌陷区具有形变时间长、破坏性大和难治理等特点。

由表2可知，M1塌陷区的形变速率为1.15～11.48 mm/d、平均形变速率为5.96 mm/d，M2塌陷区的形变速率为0.28～9.50 mm/d、平均形变速率为5.78 mm/d，M3塌陷区的形变速率为0.76～8.87 mm/d、平均形变速率为3.02 mm/d，可见M1、M2塌陷区的形变量和形变速率明显大于M3塌陷区；3个分区中各塌陷区的形变中心位置变化各异，主要与各矿区实际开挖方向有关，在单时段形变累积图中最大形变中心可视为阶段开挖的起始点，形变中心位置不变、形变加快，这可能是垂向开挖造成的。

3 结论

本文采用小基线集合成孔径雷达干涉测量(SBAS-InSAR)技术，对2016年10月-2017年3月和2017年10月-2018年3月武安市Sentinel-1A影像进行干涉与反演，分别得到两时段的时序地表累积形变量。结果分析表明：反演塌陷区主要分布在固镇—武安—矿山—綦村成矿带和鼓山—洪山—新城成矿带内，且受矿山断裂和鼓山—紫山断裂控制；各塌陷区的范围和形变中心都随时间有所变化，并出现个别塌陷区形变中心的消失和新增；2017-2018年最大形变量和平均形变量均比2016-2017年有所减小，这在一定程度上表明研究区塌陷问题整体有所缓解；M1、M2分区塌陷区范围、累积形变量和形变速率均比M3分区大，建议对M1、M2分区建立长期稳定的InSAR和原位监测体系，加强对塌陷区突发灾害问题的预警。由于缺少塌陷区原位沉降监测数据，本文未对反演速率进行精度验证，未来将结合原位监测手段开展精度验证及数据融合研究。