五龙沟矿区时序InSAR地表形变监测

杨 利，薛东剑，王海方，付 林，张 婷

(1.成都理工大学地球科学学院，四川 成都 610059；2.西南交通大学地球科学与环境工程学院，四川 成都 611756)

矿产资源的开发利用是我国经济的重要组成部分。然而矿产资源开发引发的环境地质问题日益突出，特别是矿山开采易造成崩塌、滑坡、矿震、地面沉陷等地质灾害，严重威胁到矿区周边人民的生产生命安全以及国家的财产安全[1]。对矿区地表进行动态监测可以有效获取地表移动变化规律，最大限度地降低地表沉陷带来的损失[2]。因此对矿区地表沉降进行有效监测，分析沉降特征，对开采沉陷预计与治理具有重要意义[3]。

合成孔径雷达(synthetic aperture radar，SAR)具有全天时、全天候、多波段、多极化、高分辨率、不易受气象条件制约、对地表具有一定的穿透力等优势[4]。 干涉合成孔径雷达(interferometric synthetic aperture radar，InSAR)是利用目标与两天线的位置几何关系，通过计算两次过境时SAR影像的相位差来获得数字高程模型的技术[5]。GABRIEL等[6]提出的差分干涉合成孔径雷达(differential InSAR，D-InSAR)技术则是引入外部DEM或三轨差分、四轨差分实现对地表厘米级的微小变形的监测，高精度、高时效、低成本、大范围是该技术的监测优势。然而时空去相干及大气延迟等因素使得干涉条纹无法生成或不连续，极大制约着该技术监测的精度和运用，较难完成高精度、长时间间隔地表监测[7]。为克服这些因素的影响，BERARDINO等[8]提出了小基线集方法(small baseline subset，SBAS)。 SBAS-InSAR技术通过影像的自由组合得到较多干涉，选择空间和时间基线较小、具有高相干性的干涉图，以减小地形对差分的影响，提高变形监测的准确性，能够更好地消除D-InSAR技术中的时空失相关、大气延迟相位以及地形相位误差等影响，获取研究区域毫米级的形变信息[9-10]。

本文针对矿区的地表沉降问题，结合地形考虑Sentinel-1A影像的可视情况，采用SBAS-InSAR技术获取研究区的地表形变信息。利用地学统计方法，分析沉降区在空间上的分布，结合形变速率和累计形变量，对工作区范围内的典型沉降区进行时序分析。

1 研究区概况及数据准备

1.1 研究区概况

研究区处于萤石沟-红旗沟脆韧性剪切带及其所形成的断裂构造集中带之中东段，是五龙沟矿区的一个重要构造发育区，隶属五龙沟矿区三大主要控成矿构造区带之一。研究区内出露地层主要以下元古界金水口群、上元古界青白口系丘吉东沟组和下古生界奥陶系祁曼塔格群变火山岩组为主，中元古界长城纪小庙组次之，沟谷和山前有大面积的第四系冲洪积物分布。 研究区中心坐标：95°56′4.14″E，36°12′24.2″N，面积约5 144 377 m2，总体构造线呈北西-南东向展布。矿石采矿均采用地下平硐开拓方式，随着地下矿石开采量的增加，研究区水闸东沟-黄龙沟方向(北西-南东向)形成了大范围的地下采空区，导致研究区及其附近存在崩塌、泥石流、地面塌陷等地质灾害。

图1 研究区地质构造简图

1.2 数据源

本文采用SRTMDEM 90M分辨率原始高程数据作为参考DEM为SBAS-InSAR提供地形信息和地理位置参考，利用从欧洲航天局质量控制中心下载的精密轨道数据对轨道信息进行修正，去除因轨道误差引起的系统性误差，以保证结果的准确性。由于观测几何的限制，不同方向形变的贡献值不相同。结合地形考虑影像的可视性，采用的Sentinel-1A卫星以IW模式获取的30期降轨单视复数产品进行时序反演，该数据辐射精度为1 dB，相位误差为5°，地面分辨率为5 m×20 m(距离向×方位向)，时间跨度为2017年11月18日—2019年9月21日。

2 矿区形变监测

SBAS技术是将输入影像自由组合，以多主影像的干涉对为基础，基于高相干点恢复研究区域的时间序列形变信息[12]。原理是根据每个图像中基线的时间和位置，将研究区中的SAR图像组合成小的基线子集，进行差分干涉处理，得到较多的干涉像对，保证了较高的空间密度和时间采样[13]。然后利用奇异值分解算法联合多个小基线子集求解地表形变相位最小范数意义下的最小二乘解[14]。该算法有效解决了总体法方程秩亏和单个集合时间采样率较低的问题[15]。本文利用SBAS-InSAR技术进行数据处理的流程如图2所示。

图2 SBAS-InSAR技术流程图

基于青海省的地形、气候和研究区的植被茂密程度分析，将时间基线阈值设为550 d，以避免出现无相干数据对；为避免完全空间失相关，将空间基线阈值设为45%，并进行了3D解缠，以便估算相干性低的区域、建立新的连接。利用参考DEM模拟并去除干涉图中的地形相位，采用自适应滤波方法平滑相干斑噪声，去除大气延迟相位，提高干涉条纹的清晰度。选择Delaunay MCF算法对干涉图进行相位解缠，并剔除解缠结果不理想的像对。选取相位稳定、未发生形变的68个地面控制点对干涉结果进行轨道精炼和重去平，消除了斜坡效应，对相位偏移进行修正，估算和去除残余的恒定相位。去除地形残余相位、大气延迟相位、噪声相位是SBAS技术的核心，常分成两步进行。第一次形变反演是建立线性模型计算出所有像对的残余地形及地表形变速率，利用奇异值分解法解算相关参数，得出干涉像对平均形变速率和高程误差；第二次形变反演是通过大气滤波估算和去除大气相位和地形残余相位，获取更加准确的时间序列最终位移图。经过SBAS两次反演得到的地表形变信息是在SAR坐标系中，将反演结果通过地理编码转化为地理坐标，最终获得研究区的地表形变图。

3 实验结果分析与讨论

在雷达图像干涉处理的基础上，利用GIS工具，结合研究区SBAS反演的年平均沉降速率、累计形变等信息，进而分析工作区的沉降分布、沉降量及典型区域沉降中心的形变时间序列，并对沉降区的剖面进行时序形变分析。

3.1 矿区整体形变分析

图3 研究区形变速率图

在监测时间范围内，获得的五龙沟矿区年平均速率如图3所示，该范围内地表移动速率相对来说比较稳定，但沉降区表现明显。 工作区的平均沉降速率为-1.40 mm/a，最大沉降速率为-11.74 mm/a。如图4(a)所示，-2～0 mm/a为主要的形变速率区间，占工作区面积的72.35%，而形变速率为-11.8～-10 mm/a的强烈形变区域仅占工作区面积的0.54%。

五龙沟矿区以2017年11月18日为基准，在2018-05-05期开始出现较明显的沉降区域，随着时间的推移，采矿活动的继续开展及矿区自身重力的影响，沉降区域沉降量增加，沉降范围逐渐扩大，到2019年9月21日，五龙沟矿区范围内最大沉降量可达-21.68 mm。如图4(b)所示，工作区范围内的整体型变量不大，主要集中在-3～0 mm，占该工作区面积的61.57%，其中，沉降量超过-6 mm的面积约为364 344 m2，占工作区面积的7.09%。

图4 形变面积占比图

3.2 典型区形变分析

目前工作区有原矿堆存场、选矿厂、尾矿库、地下采矿区(水闸东沟采区和黄龙沟采区)、矿区道路、硐口工业场地等设施。如图5所示，工作区范围内主要存在两个明显沉降区，且分布在工作区的西部，这与实际情况相符。在此基础上圈定两个典型沉降区域，分别为沉降区A和沉降区B。

图5 典型沉降区分布图

沉降区A的总面积约为124 575 m2，占工作区面积的2.42%，该沉降区形变较小，主要集中在-7～-10 mm之间，其所占面积是沉降区A面积的53.67%；沉降区B的总面积为155 311 m2，占工作区面积的3.02%，形变量明显大于沉降区A，其形变量在-22～-19 mm范围内的面积大约是沉降区A的7倍。工作区内两沉降区的形变量面积占比情况见表1。

表1 典型沉降区形变量面积占比

以2017年11月18日为基准，在监测时间范围内，地表监测点的沉降没有出现突变现象，整体沉降位移量最大不超过30 mm。沉降区A的最大形变点位于沉降区中部偏左下方，如图6(a)所示，其最大沉降值为-20.27 mm。该点以近似稳定的速率呈现线性沉降的趋势。沉降区B的最大形变点位于沉降区的左上方，如图6(b)所示，其最大沉降值达到-21.68 mm。该点的沉降呈现非线性特征变化。沉降区A主要采用分段空场法、无底柱分段崩落采矿法回采。在2017-11-18—2019-04-06时间段内，该区域空区已形成，随着矿区地下开采的扩展，采空区的暴露面积、围岩所承受的应力以及围岩蠕变量相应变大，诱发顶板出现冒落，造成地表微小形变。2019-04-06期后较平稳的原因是该区通过自然冒落或强制崩落围岩方式进行处理。沉降区B采矿方式与沉降区A相同。 在监测段内正开采3 505 m中断，平衡拱承受的应力没有达到破坏极限，顶板未出现冒落，因此地表呈现的沉降变化不明显。随着工作的推进，平衡拱受空区倾向和走向跨度扩展的影响，拱的受压逐渐接近极限，超过极限后平衡拱被破坏，且由于矿区自身重力的影响，冒落高度增大，顶板岩层移动变形造成地表微小沉降。相较于沉降区A，沉降区B在局部地段未支护，其稳固性较差，地表沉降趋势较为明显，其中2019-08-04处形变起伏最为显著。

地下空区形成后，采矿区会存在对应的沉陷中心。竖直方向上的变形将导致地面路基、土质的松弛，影响其承受能力，继而引起和加大水平方向的倾斜和拉伸变形。沉降区A以最大形变点位为中心逐渐向四周扩散，其中向东北向扩散的趋势较明显；沉降区B以最大形变点位为中心呈带状沿西北-东南方向逐渐向四周扩散。在沉降区扩散趋势范围内分别绘制沉降区A和沉降区B的剖面线，通过时序分析发现，两沉降区中心沉降(图7)均呈现逐年增强趋势。沉降区A呈漏斗状逐渐向四周扩散，且向四周扩散趋势明显。沉降区B则是呈现盆状逐渐向外围扩散，扩散趋势较为平稳。

图6 沉降区最大形变点形变序列图

4 结 论

1) 利用SBAS-InSAR方法得到了青海省五龙沟矿区2017年11月至2019年9月期间各时间段的地表微小缓慢形变，获取了矿区地面沉降的时间序列图及其空间展布。监测结果显示工作区整体形变不大，但存在两个明显沉降区。

2) 形变干涉结果表明五龙沟矿区总体比较稳定。在监测时间范围内，工作区地表整体平均沉降速率为-1.40 mm/a，其中形变速率为-2～2 mm/a的区域占实验区面积的81.54%；地表形变量为-3～3 mm的区域占工作区面积的79.02%。工作区内最大形变速率为-11.8 mm/a，最大沉降量为-21.68 mm，其中沉降量超过-6 mm的面积为364 344.35 m2，仅占工作区面积的7.09%。

3) 沉降区A呈现线性特征，以最大沉降点为中心向四周扩散，且向东北向扩散的趋势较明显。沉降区B呈现非线性特征，以最大沉降点位为中心，沉降区呈带状沿西北-东南方向逐渐向四周扩散。