基于时序InSAR技术的甘肃金川铜镍矿区地表沉降研究

李强，王卫红，杨宁，熊高翔 ,雷欢欢，李凯

(1.四川省冶金地质勘查院，成都 610051；2.国家遥感中心绵阳科技城分部，四川 绵阳 621010；3.西南科技大学环境与资源学院，四川 绵阳 621010)

0 引言

甘肃省金川铜镍矿区位于阿拉善台块的南部边缘隆起区，南邻北祁连山加里东地槽边缘拗陷带。矿区生产采用无底柱分段崩落式采矿方法，在采矿爆破时会加速地表的沉陷，形成地裂缝。同时地下开采会引起岩层及地表产生移动和变形，诱发各种地质灾害[1]。因此，掌握地表变形规律，能够有效保证人民群众的生命财产安全、国家经济的可持续发展以及社会的稳定[2]，对矿区灾害预警、开采沉陷控制和治理有重要指导意义。

矿区地表沉陷是一个长期过程，具有地形地质条件复杂、沉降速度快等特点，所以利用GPS测量、大地水准测量和全站仪测量等传统的地表形变监测方法监测矿区沉陷区面临诸多难题。第一，工业法采矿多采用放炮方式，导致水准点稳定性存在问题，加之采矿是一个长期过程，传统测量周期长，工作量大，耗时耗力[3]；第二，传统监测方式只是在局部得到沉降信息，不利于体现整个矿区的沉降规律[4]；第三，传统测量具有地域局限性，需要预先知道沉降位置才能布置观测点。合成孔径雷达干涉(Interferometric Synthetic Aperture Radar)技术，即InSAR测量技术是近30年来迅速发展起来的一种空间对地观测技术，基于合成孔径雷达复数影像的相位信息获取地表变化信息，可用于监测厘米级或更微小的地球表面形变[5]。相比可见光和红外遥感等传统测量技术，合成孔径雷达干涉技术具有覆盖范围广、精度高、不受天气状况影响、空间采样密度高等优势。目前，InSAR技术已广泛应用于区域地表形变监测[6]、隐患点探测[7]、冰川冻融监测[8]、洪水淹没区变化监测[9]等多个领域。

对于InSAR技术在矿山监测方面的应用，国内外学者展开了大量研究，并取得高精度形变的观测结果。1993年，Massonnet等[10]利用ERS-1卫星的SAR图像，对加利福尼亚州兰德斯地震产生的运动进行监测研究，监测精度高于以前的空间成像技术，精度优于3 cm；2005年，Tomás等[11]在用D-InSAR技术监测1993—1995年期间塞古拉河周边地面沉降时发现，某些区域的沉降达8 cm，对预测和治理该区域地面沉降有一定指导作用；2010年，Herrera等[12]利用相干像素技术(CPT)绘制和监测La Union矿区的地面运动，证明此方法可用于废弃矿区大部分地面运动；2014年，Saygin等[13]分析了D-InSAR的缺陷，利用PS-InSAR技术对土耳其西北部成交地区采煤沉陷区域进行监测，充分说明了PS-InSAR技术在矿区沉降监测具有指导性意义；2017年，朱建军等[14]着重对InSAR监测方法进行了分类，分析了PSInSAR、SBAS-InSAR及D-InSAR，提出改进的干涉测量技术可以用于矿区的沉降监测；2018年，李达等[15]利用SBAS-InSAR技术对陕西省榆林市某矿区西北部进行地表沉降监测，证明了SBAS-InSAR技术在矿山沉降监测的可靠性，为矿区地表沉降监测提供新手段；2019年，黄长军等[16]采用持续散射干涉测量法(PSI)，避免了传统干涉测量的限制，有效揭示了山东葛亭煤矿开采沉陷的演变过程，证明PSI法是测量采矿区沉陷的有效方法，可为矿山灾害提供预警。

1 研究区概况与数据来源

金川铜镍矿床由4个矿区组成，龙首矿是主力矿山之一，20世纪60年代建成投产，在生产镍铜的同时回收铂族金属，是我国铂族金属的主要生产基地。金川铜镍硫化物矿床的金属矿体产于超基性岩体中，部分矿体产于含矿岩体底板与围岩接触带附近。在矿体下方存在断裂构造，由于开采强度较大，故产生的地质灾害分布广。龙首矿西二采区位于Ⅰ矿区以西，矿内断裂构造发育，其中F8断裂规模最大，角砾岩层厚20～40 m，在4行勘探线直接与矿体接触，地质环境复杂。受断裂影响，在2016年出现明显的地面开裂现象，地面裂隙主要分布在5行—7行这一区段，呈环形闭合环，直接威胁矿山的生产安全，地下的采矿作业因此停采。在停采后，为监测其地面沉陷情况，矿区管理人员对沉陷区进行维护，在裂缝附近安装控制点进行重点监测，直到2018年年底沉陷区趋于稳定后，选用无底柱分段崩落法进行试验性复采。在此之前，金川矿区的地表沉降监测研究通常采用传统变形测量方法，故采用InSAR技术对金川矿区进行地表监测具有必要性。

本次研究工作利用2018年1月—2018年11月的Sentinel-1A降轨影像和同期的升轨影像对西二采区进行沉降监测。所有影像均为宽幅单视复数影像数据(slc)，影像覆盖范围如图1所示，斜距向和方位向分辨率为3.7 m和22.7 m，单景影像覆盖范围约为251 km×251 km。本文采用了C波段升轨影像共26景，并采用地面分辨率为30 m的SRTM DEM数据来消除干涉相位中的地形相位。

图1 金昌市区域Sentinel-1A雷达卫星覆盖示意及区位图

2 数据处理

2.1 SBAS-InSAR和时序D-InSAR基本原理

SBAS-InSAR(小基线集)技术是对于长时间序列上的一组SAR复数图像根据一定的基线约束条件进行分组，通过控制空间基线的长度来提高干涉图的相干性，对差分干涉图进行多视处理降低噪声，提取高相干性单元，使用奇异值分解法求得影像序列间地表形变速率的最小范数最小二乘解[17]，SBAS详细技术原理见文献[18]。

假设按时间顺序(t1，…，tn)在同一地点获取了N+1幅SAR影像，同时每一短基线子集至少有2幅影像，可生成M个干涉图，M满足以下不等式

(N+1)/2≤M≤[N(N+1)]/2

(1)

假设j(j=0，1，…，N)干涉图是由tA、tB时间获取的影像生成(tB>tA)，去除平地效应与地形相位影响后, 干涉图j中任意一个像素的干涉相位可表示为

Δφj=φ(tB)-φ(tA)=φdef，j+φatm，j+φflat，j+φnoise，j

(2)

式中，φ(tB)、φ(tA)分别表示其时刻相位值；φdef，j表示视线向(LOS)时刻的形变相位；φfiat，j表示残余地形相位误差；φatm，j表示大气相位误差；φnoise，j表示噪声误差。因此，相对基准时间t0，每幅SAR影像对应的差分相位时间序列φ与每个差分干涉图相位时间序列Δφ存在以下关系

Δφ=Aφ

(3)

式中，A为M*N系数矩阵，每个差分干涉图相位时间序列Δφ=[Δφ(t0)，Δφ(t1)，…，Δφ(tM)]T，每幅SAR影像对应的差分相位时间序列φ=[φ(t0)，φ(t1)，…，φ(tN-1)]T。

再通过矩阵的奇异值分解法(SVD)求出最小范数意义上的最小二乘解，所以观测时刻形变量：

(4)

式中，A+=(ATA)-1AT。

合成孔径雷达差分干涉测量技术(D-InSAR)是干涉测量技术和合成孔径雷达技术的结合，其技术的核心思想是利用重复轨道观测得到干涉相位，通过差分处理去除2次观测相位中的共有量[19](地形相位、平地相位、大气相位以及噪声相位)，传感器视线方向上的沉降值Δr：

Δr=(λ/4π)φdef

(5)

式中，λ为雷达传感器工作波长；φdef为有2次观测期间沿着雷达实现方向地表形变引起的形变位移。

将公式(5)求出的Δr转化成垂直方向上的沉降值Δh：

Δh=Δr/cosθ

(6)

式中，θ为传感器入射角。

2.2 数据处理

基于SAR Scape 软件实现SBAS和时序D-InSAR的处理。时序D-InSAR技术通过影像配准、去平地效应、地理编码、DEM相位差分等一系列步骤得到形变图。图2所示的SBAS技术流程主要包括小基线组合选择生成连接图，生成差分干涉图，Delaunay MCF解缠，选取GCP点，线性估算位移速率和DEM校正系数，大气相位屏的估计和减法，总形变量估计等。

图2 SBAS技术处理流程图

3 结果与分析

研究中使用SBAS技术和二轨法时序D-InSAR技术获取研究区2018年1月—11月的地表沉降速率。为了更直观地展示研究区的地表形变情况，在envi中进行地理编码，将形变结果叠加到谷歌影像得到矿区沉降图(图3)，图3a为SBAS技术沉降结果图，图3b为时序D-InSAR技术沉降结果图，将2018年1月作为形变区域参考，图中1区为西二采区5行—7行，2区为露天采场老矿坑。为了便于观察沉降变化趋势，在1区和2区选择了沉降漏斗中心做时序分析，结果见图4，分别经过D-InSAR与SBAS技术计算，得到其形变速率图(图5，图6)，由于1区相较于2区沉降较小，故分离出1区单独出图(即图5b和图6b)。以上结果表明，1区和2区中均有沉降漏斗出现，1区不明显，2区则很明显。

3.1 西二采区5—7行地表沉降分析

由图3至图6可以看出，2018年1月—11月西二采区整体较为稳定，使用时序D-InSAR方法获得此区域最大沉降量为11 mm，年平均最大沉降速率为-14 mm/a，有明显沉降漏斗的部分，沉降量范围为3 ～11 mm。SBAS方法获取最大沉降量为13 mm，年平均最大沉降速率-14 mm/a，沉降量范围为1～14 mm。在停采之后，2种方法监测结果显示，在2018年此区域仍然存在缓慢的形变，依据我国形变监测基准[20]以及本矿区地质情况，沉降速率<30 mm/a时，处于安全范围，说明5—7行区域呈平稳态势。由于形变较小，2种方法监测结果基本一致。

图3 累积沉降量结果图

图4 累积沉降量曲线

图5 采用SBAS技术获得的沉降速率图

图6 采用D-InSAR技术获得的沉降速率图

3.2 老矿坑地表沉降分析

由图3可以看出，老矿坑东南部有明显沉降，相较于5行—7行，主要是边坡附近受雨水侵蚀等导致沉降并伴随滑坡，基于时序D-InSAR方法获得其最大沉降量在73 mm，最大沉降速率为-83 mm/a；SBAS方法获取其最大沉降量在78 mm，而最大沉降速率为-85 mm/a，两者监测结果均显示老矿坑沉降量大。从图5和图6可以看出，时序D-InSAR方法与SBAS方法得到的沉降速率图明显不同，其大沉降位置在SBAS方法所出得出的结果更为精确。同时试验发现，老矿坑西北部有小部分区域存在轻微抬升的现象，一年抬升量约为20 mm，现场踏勘结果显示边坡有轻微隆起，而西南部有滑坡痕迹，老矿坑北部有轻微沉降。从2区结果可以看出，2种方法均能监测到老矿坑的沉降现象，但由于时空基线以及大气效应等误差[21]，致使时序D-InSAR监测结果差于SBAS监测结果，而SBAS方法能减轻时空失相关、大气延迟等影响[22]能获取更精确的沉降值以及沉降位置。

4 结论

(1)经过2种方法实验结果对比分析，依据当地地质情况以及国家标准，西二采区5行—7行地表总体呈平稳态势，在过去一年最大沉降量为11～13 mm，最大沉降速率为-14 mm/a；老矿坑区域出现了滑坡现象，最大沉降量在75～78 mm，最大沉降速率为-85 mm/a。

(2)停止采矿后，原采矿区域地面沉降明显减缓，基本处于稳定状态，印证了时序InSAR在矿山地表沉降监测的可行性。

(3)虽然2种方法都能得到沉陷区地表沉降规律，但由于SBAS技术能够减轻时空失相关、大气延迟等影响，且从2区域可以看出SBAS监测结果相比时序D-InSAR监测结果更为理想，故SBAS-InSAR方法更能在后期大幅度沉降出现时，精确地反映开采沉陷区沉降速率与沉降位置。