雷达干涉测量对金川西二矿地表二维形变监测

王志岗, 周文韬

(1.四川科技职工大学 应急管理系，四川 成都 610000；2.西南科技大学 环境与资源学院，四川 绵阳 621010；3.国家遥感中心绵阳科技城分部，四川 绵阳 621010)

矿产资源的开采一方面可以促进我国经济建设和社会发展，同时也会破坏矿区地下岩层原有的应力结构，从而造成一系列地质灾害[1-4]。2016年3月，甘肃省金昌市金川西二矿1 610 m水平面5～7行范围的胶结充填体发生大面积垮塌并贯通至地表，从而导致停产[5]。从2019年起，该区域开始使用无底柱分段崩落法进行采矿，但该方法对地表影响较大[6]。截至2021年7月，该区域地表已出现3个大小不一的塌陷坑，对附近矿区人员的生命财产、建筑设施等具有一定的安全隐患。因此，开展地表形变监测工作，对后续工作开展具有重要的意义。

不同于传统大地测量方法，合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar, InSAR)技术自1989年首次提出以来，就因其具有空间分辨率高、覆盖范围广等优点在矿区地表形变监测中受到广泛关注[7-10]。经过国内外学者的研究，一些新的InSAR技术逐渐出现，例如合成孔径雷达差分干涉测量(Differential InSAR, D-InSAR)技术、永久散射体合成孔径雷达差分干涉测量(Persistent Scatterer InSAR, PS-InSAR)技术、短基线集合成孔径雷达差分干涉测量(Small Baseline Subset InSAR, SBAS-InSAR)技术和人工角反射器差分干涉测量(Corner Reflector InSAR, CR-InSAR)技术等[11-15]。然而，上述技术仅能监测研究区地表沿SAR卫星视线向(Line of Sight, LOS)形变，而矿区复杂的地形信息难以靠单一LOS向形变获得[16-18]。

因此，本文将利用SBAS-InSAR技术和InSAR几何分解原理，基于升降轨Sentinel-1A数据解算甘肃省金昌市金川西二矿地表二维形变场，对该区域地表形变情况进行分析，为矿区开采沉陷监测和预警工作提供思路和方法。

1 研究区概况及数据基础

1.1 研究区概况

金昌市(101°4′35″～102°43′40″E，37°47′10″～ 39°00′30″N)位于我国西北部，是甘肃省下辖地级市，也是我国重要的硫化镍铜矿资源地之一，被誉为“祖国的镍都”[19]。本文研究区选为金昌市金川西二矿，位于金昌市区西南部，如图1所示。该区域降雨量少，日照充足，平均海拔达1 700 m左右，地下水较少，生态环境总体偏弱。矿区主要以超基性岩型矿体以及围岩中接触交代的矿体构成。由于地下长期采矿和F8断层不断活动，研究区地表形成较大塌陷区，对地下开采产生极大的影响。

1.2 数据基础

本文采用由欧洲航天局(European Space Agency, ESA)针对哥白尼全球观测计划研制的Sentinel-1卫星星座的Sentinel-1A数据作为实验数据。自2014年4月发射以来，该卫星可以为用户提供免费的雷达影像，并能够对全球范围的地表形变、海洋、冰川等进行监测[20-22]。

实验分别选取Sentinel-1A升轨25景和降轨27景影像，统一时间跨度为2020.08.19～2021.06.27，实验数据参数如表1所示。实验处理过程统一时间基线阈值120 d，临界基线阈值45%。其中升轨数据以2020.12.17期影像为主影像，生成187个干涉对；降轨数据以2020.08.19期影像为主影像，生成215个干涉对。同时，在数据处理过程中引入外部数字高程模型(Digital Elevation Model, DEM)数据和精密轨道数据来消除地形起伏和轨道误差带来的影响。

2 技术与方法

2.1 SBAS-InSAR技术原理

在过去几十年里，InSAR技术因具有监测范围广、精度高等优点而被广泛应用。但随着技术的进步，人们发现传统InSAR同样伴随大气效应、时空失相干等误差的缺点[23-24]。SBAS-InSAR技术由Berardino和Lanari于二十世纪初首次提出[25]，该方法可以在一定程度上克服传统InSAR遇到的难题，有效改善时空失相干和大气延迟等问题。

SBAS-InSAR技术的主要原理如下：假设有t0,t1，…,tn景按时间序列排序的N+1幅SLC (Single Look Complex)影像，以任意一景影像为超级主影像对所有影像进行配准和重采样，通过设定时间基线阈值和临界基线阈值可得到M个干涉对如下：

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

(1)

然后利用精密轨道数据和外部DEM数据去除平地效应和地形效应，得到M个差分干涉图。若有时间tA>tB，以t0为参考时间，则第j幅干涉图在x处的干涉相位值为：

(2)

式中，λ为雷达中心波长；φ、μ为干涉相位。则LOS向形变D可表示为：

D(tB,x,r)-D(tA,x,r)=vi(tB-tA)

(3)

式中，vi为从tA至tB的LOS向平均形变速率。由上式可转化为：

Bv=δφ

(4)

式中，B为系数矩阵。

在SBAS-InSAR技术解算过程中，由于多主影像模式容易导致矩阵B秩亏，因此采用奇异值分解法(Singular Value Decomposition, SVD)求解，使得计算有唯一解[26]。该方法的技术流程如图2所示。

2.2 基于升降轨InSAR获取矿区地表二维形变场

根据2.1节的技术方法，我们仅能获取单一轨道沿卫星LOS向的一维形变，LOS向形变的分解关系如图3所示。在复杂的矿区地表形变过程中，LOS向的一维形变往往不能准确反映真实形变特征。当结合多轨道SAR卫星数据联合解算时，则可能获取研究区的多维形变场。

由图3可知，地面某点在东西向、南北向和垂直向的位移是LOS向形变分量。当同时考虑地面点水平和垂直位移时，则有下式为：

DLOS=-sin(θ)cos(α)DE+sin(θ)sin(α)DN+DVcos(θ)

(5)

式中，θ为卫星入射角；α为卫星方位角。

图3 InSAR几何分解原理

令S=[Sx,Sy,Sz]，则Sx=-sin(θ)cos(α),Sy=sin(θ)sin(α)和Sz=cos(θ)分别为三个方向的投影矢量。但通常情况下，SAR卫星采用极轨法飞行(即近南北向飞行)，导致InSAR技术对南北向形变监测并不敏感，因此可忽略南北向形变，利用升轨和降轨监测结果解算矿区地表二维形变[27]，即：

(6)

式中，θ1和θ2分别为升轨和降轨卫星入射角；α1和α2分别为升轨和降轨卫星方位角。

3 结果与分析

基于2.1节的SBAS-InSAR技术原理，本文首先获取了Sentinel-1A数据在升轨和降轨LOS向形变速率图，如图4所示(其中图4(a)为升轨形变速率，图4(b)为降轨形变速率，以抬升为正、下沉为负方向)。可以看出，远离矿体的区域大部分趋于稳定状态，无明显形变特征。由图4(a)可知，升轨监测结果的形变区域主要集中在矿体西南部区域。而在图4(b)的降轨监测结果中，形变区域主要集中在南部以及西南部区域，矿体西部有局部抬升趋势。造成这种差异则是受轨道影响，同时在矿体西部区域有一个明显的山坡，从而导致了这种差异。整体来看，研究区在升轨结果中最大年平均形变速率达到-126.37 mm/y，降轨最大年平均形变速率为-84.03 mm/y，整体呈下沉趋势，且沉降中心汇集于矿体西南部。

图4 升降轨SBAS-InSAR形变速率图

然而，不同轨道的监测结果具有一定差异，不能反映地表的真实形变情况。因此本文采用2.2节的方法，解算了研究区地表二维形变场，如图5所示(其中图5(a)为东西向形变速率，图5(b)为垂直向形变速率，以向东移动为正、向西移动为负方向)。

由图5(a)可知，研究区地表西部区域向东移动，最大年平均形变速率达到103.73 mm/y；东部区域向西移动，且最大年平均形变速率为-32.40 mm/y，最后整体汇集于落矿区域并趋于稳定。由图5(b)可知，研究区整体呈下沉趋势，最大年平均下沉速率达到-135.88 mm/y，最小为-2.47 mm/y，下沉区域集中于矿体西南部。

图5 矿区地表东西向和垂直向形变速率图

根据形变结果来看，矿区形变区域整体位于落矿区域上方及南部区域，不是落矿区域的正上方。经实地考察与分析，总结了造成这种现象的原因，主要有以下两点：

(1)该区域西南部存在F8断层[6]，且地形呈西南高、东北低的特点，因此导致研究区西部和西南部的形变速率更大；

(2)地下开采影响反应至地表需要一定的过程，具有滞后性，因此形变区域与落矿区域具有一定偏差。总体上来说，研究区的地表二维形变结果符合开采沉陷规律。

4 结 论

本文利用Sentinel-1A的升降轨SAR影像，采用SBAS-InSAR技术首先得到了甘肃省金昌市金川西二矿地表LOS向年平均形变速率。通过对升降轨结果差异的分析，结合InSAR几何分解原理获取了研究区地表二维形变场。结果显示，该区域在2020.08.19～2021.06.27期间，地表东西向年平均形变速率达到103.73 mm/y，垂直向年平均形变速率达到-135.88 mm/y，且在地表形成塌陷坑。此外，通过考察与分析，形变结果基本符合开采沉陷规律，验证了方法的可靠性。