1. 引言
矿山开采是人类社会主要的经济活动之一，可以为社会发展带来显著的效益，但开采过程可能会造成矿区及其周边地表出现不均匀形变，进而导致边坡失稳、构筑物沉陷、泥石流、洪涝等灾害隐患 [1] [2] 。因此，针对矿山及其所在地区进行地表形变调查和监测，是维护和指导矿山安全生产必不可少的手段。
矿山开采地表形变监测传统上主要依赖经纬仪坐标测量、水准仪测高等技术，但这些传统手段存在监测成本高、空间覆盖度低等劣势。近年来，全球导航卫星系统(Global Navigation Satellite System, GNSS)、合成孔径雷达差分干涉(Differential InSAR, DInSAR)等现代大地测量技术在矿山地表形变监测中逐渐得到重视。特别地，相较于存在基准点稳定性不足、观测周期长、垂直向精度相对较低等问题的GNSS技术，DInSAR技术具有测量效率高、覆盖空间范围大、结果分辨率高等优势 [3] 。目前，InSAR技术经过数十年的发展，基于多景SAR影像的时序InSAR (Time Series DInSAR, TS-InSAR)技术已基本成熟，能够有效克服常规DInSAR技术普遍存在的时空失相干、大气延迟和相位噪声等问题，地表形变监测的精度可达到毫米级。其中，小基线集InSAR (Small Baseline Subsets InSAR, SBAS-InSAR)技术具有很强的代表性，进而成为了面向地表沉降监测的主要手段之一 [4] 。
目前，已有大量学者将 SBAS-InSAR技术应用于大范围地表形变监测研究中，并在矿山及其周边地区取得了可靠的监测结果。梁思语等 [5] 以山西省某矿区为研究对象，应用SBAS-InSAR技术对该地区进行了地表形变监测，结果与实测数据对比表明SBAS-InSAR技术在矿区残余变形反演方面具有良好的可靠性。陈志达等 [6] 基于SBAS-InSAR技术，使用ALOS-PALSAR和RADARSAT-2数据处理得到了山西省浑源县形变速率与形变时间序列，并分析得出该地区内矿区地表形变与煤矿开采范围在空间上具有较高的一致性。张香凝等 [7] 以辽宁省沈阳市蒲河煤矿为研究对象，使用SBAS-InSAR技术得到了该地区因采煤引起地表形变的结果。此外，研究还结合数值模拟技术对观测结果进行了模拟分析，结果表明SBAS-InSAR技术处理结果与模拟分析结果具有一致的分布规律，反演结果与实际情况接近。
鞍山市地处辽宁省中部，是我国主要的钢铁工业基地，辖区内拥有东鞍山铁矿、齐大山铁矿及大孤山铁矿等多个矿区，矿产资源丰富。鞍山露天矿区长期处于铁矿开采作业过程中，矿区开采在推动当地经济与资源发展的同时，也带来了地质灾害风险。因此，针对鞍山地区，一些学者对其地表形变进行了基于Sentinel-1A影像的监测。曹翠华等 [1] 使用SBAS-InSAR技术对14景Sentinel-1A数据进行了处理，获得了鞍山市及周边的2017~2018年地表形变速率，并结合相关资料分析了引起地表形变的主要原因包括铁矿、地下水开采和地质条件的影响。李鸣庚等 [2] 利用2020~2021年获取的50景Sentinel-1A数据，使用优化分布式目标InSAR (Distributed Scatterers InSAR, DS-InSAR)技术对鞍山地区齐大山、大孤山和鞍千矿矿区进行了形变监测。此外，2015年苏畅 [8] 利用1996~2000年共17景ERS-2数据，通过SBAS-InSAR等三种TS-InSAR技术对鞍山露天矿边坡形变进行了监测，并基于重力轨道对配准精度的影响和DEM数据的应用效果，揭示了TS-InSAR在露天矿区边坡形变监测中的有效性和可行性。田辉 [9] 在2015年使用2006~2010年22景ALOS-PALSAR数据。基于SBAS-InSAR技术对辽阳-鞍山地区典型区域进行了形变监测和沉降机理分析。
然而，以上研究都是基于中等分辨率SAR影像进行InSAR监测(如表1所示)，并且在2017年之前由于SAR卫星在研究区采样相对较少，对于2011年至2017年间鞍山市及周边地表形变还没有文献进行报道，相关研究内容也没有对矿区地表的精细形变进行揭示(如表2所示)。因此，选择相应时段内的高分辨率SAR影像对鞍山地区进行地表形变监测能够一定程度上填补该时段内研究区形变监测的空白，同时也可以探究高分辨率SAR影像在露天矿区进行地表精细形变监测的可行性和可靠性。
Table 1. Comparison of resolution of SAR images used in existing studies in study area
表1. 研究区已有研究使用SAR影像分辨率对比(改自文献 [10] )
Table 2. Comparison of using SAR image periods and the main study in existing studies in study area
表2. 研究区已有研究使用SAR影像时段与主要研究内容对比
基于以上研究意义，本文利用SBAS-InSAR技术，选取2015年5月至2016年8月共37景COSMO-SkyMed高分辨率SAR影像，获取了鞍山地区的年均形变速率及形变时间序列；此外，结合Google Earth卫星影像对东鞍山、齐大山矿区及市区铁西区地表覆盖信息的解释，本文还分析了研究区地表形变可能的原因。
2. SBAS-InSAR技术原理
SBAS-InSAR技术由Berardino等 [11] 在2002年提出，是一种基于最小二乘(LS)的多主影像TS-InSAR算法 [12] ，其利用对相干目标(coherenttarget, CT)的识别提取、相位解缠及分析得到时序形变，通过时空基线阈值的最优选择选取差分干涉对，并对CT点使用线性相位变化模型进行建模与解算，在消除常规DInSAR失相干影响的基础上获得地面的形变时间序列 [11] [13] [14] [15] [16] 。
SBAS-InSAR首先计算研究区域不同时段多景SAR影像的时空基线，选取最佳的时空基线阈值得到干涉对；之后将得到的干涉对进行差分干涉及相位解缠；最终根据干涉图自由组合所形成子集的情况，将全部干涉图组成相位方程，通过最小二乘或奇异值分解估计形变参数。在数据处理中，通常会利用时空滤波对大气延迟进行去除，从而分离出非线性形变。研究区域的形变信息则为该非线性形变与估计低频形变的总和 [15] 。
SBAS-InSAR主要技术流程如图1：
Figure 1. Flow chart of SBAS-InSAR technology
图1. SBAS-InSAR技术流程图
3. 研究区概况与数据源
3.1. 研究区概况
研究区域(图2)位于鞍山市境内，地理坐标为E122˚84'~E123˚16'，N41˚01'~N41˚26'。区域内主要包括立山区、铁西区、铁东区以及千山区部分区域，以及东鞍山铁矿和齐大山铁矿等多个矿区。研究区整体地势平坦，东南高西北低，自东南向西北呈倾斜态势，整体海拔为11~528米。
Figure 2. Topographic map of the study area
图2. 研究区地形图
3.2. COSMO-SkyMed数据
COSMO-SkyMed卫星系统由意大利航空局与国防部合作研发，共包括4颗在轨卫星，每颗卫星均配备工作于X波段的多模式高分辨率合成孔径雷达 [17] [18] 。COSMO-SkyMed可提供分辨率1米、3米、15米、30米和300米的雷达数据，具备全天时、全天候和全球范围的对地观测能力，在时序InSAR领域被广泛应用 [19] [20] 。
基于其分辨率高、覆盖范围广和重访能力高等优势，本研究选取2015年5月14日至2016年8月12日获取的37景COSMO-SkyMed影像对研究区的地表形变进行监测，共生成128个干涉对。如图3所示为SBAS-InSAR处理过程中的时空基线图，其中最大垂直基线为433米，最大时间基线为160天。
Figure 3. SBAS-InSAR spatio-temporal baseline map
图3. SBAS-InSAR时空基线图
本文基于GAMMA软件进行SAR数据处理，其中DEM采用由高级陆地观测卫星ALOS (Advanced Land Observing Satellite)获取的AW3D 30 (ALOS World 3D-30 m)数据，其水平分辨率约为30米，高程精度约为5米，能够覆盖全球所有的土地尺度 [21] [22] [23] 。AW3D 30数据在SBAS-InSAR数据处理过程中，本研究使用费舍尔信息矩阵(Fisher Information Matrix)作为最小二乘反演权重，空间相干性掩膜阈值为0.4，时间相干性掩膜阈值为0.7。
4. 实验结果与分析
4.1. 研究区地表形变场
利用SBAS-InSAR技术在研究区域内获取的形变速率图如图4所示。在2015年5月14日至2016年8月12日期间，研究区域内除矿区外整体形变较为稳定，而矿区及其周边地区则沿卫星视线向的沉降现象显著研究区形变速率范围在−6.0 cm/yr~5.6 cm/yr之间。
为进一步分析形变结果，选取图4标注的A、B、C三个典型形变区作为重点区域进行讨论。
Figure 4. Surface deformation rate of Anshan City acquired based on SBAS-InSAR
图4. 基于SBAS-InSAR获取的鞍山市地表形变速率
4.2. 重点形变区结果分析
重点形变区A (以下简称A区)主要包括了齐大山铁矿及周边区域，其Google Earth卫星影像与形变速率结果如图5所示，形变速率范围在−6.0 cm/yr~2.9 cm/yr之间。图6则为A区内典型形变点a1、a2在2015年5月14日至2016年8月12日期间的形变时间序列图。
A区内，齐大山铁矿及其周边存在以沉降为主的形变，且在铁矿内部和铁矿边界沿线较为显著。由于区域内铁矿面积较大，因而沉降区域分布也较为广泛。此外，由图6可知，a1、a2在2015年5月14日至2016年8月12日期间在视线向存在明显沉降，平均形变速率分别为−4.6 cm/yr和−4.7 cm/yr，沉降程度较为接近，累积沉降量约分别达到6.1 cm和6.5 cm。通过观察时序曲线可知，a1在2015年8月至2016年2月处于较大程度的持续沉降。而a2在2015年8月至2015年12月存在大程度持续沉降，但在2016年1月内曲线呈短暂上升。结合卫星影像可推测在2015年12月至2016年1月间，a2及其周边地表可能受附近采矿工期调整或其它因素等原因影响，从而导致了地表持续沉降的中断。
重点形变区B(以下简称B区)主要覆盖东鞍山铁矿及周边区域，该区域形变速率范围在−3.8 cm/yr~4.7 cm/yr之间，Google Earth卫星影像及形变速率结果如图7所示。与A区类似，B区内，东鞍山铁矿及其周边也存在以沉降为主的不同程度形变。在铁矿内部，北侧和西南侧形变较为明显，且越靠近铁矿内部形变速率越快。此外，在铁矿外东南侧也存在显著的LOS向沉降信号。结合卫星影像可知，该区域沉降沿远离铁矿方向沉降速率更快。
Figure 5. Google Earth satellite image of Qidashan Iron Mine (area A) (a) and surface deformation rate of Qidashan Iron Mine (area A) from SBAS-InSAR analysis (b)
图5. 齐大山铁矿(A区) Google Earth卫星影像(a)和基于SBAS-InSAR获取的齐大山铁矿(A区)地表形变速率(b)
Figure 6. Time series accumulative subsidence of points a1 and a2 in area A from SBAS-InSAR analysis
图6. SBAS-InSAR获取的2个A区特征点a1、a2时序沉降量
图8为B区内典型形变点b1、b2在2015年5月14日至2016年8月12日期间的时序累积沉降量，b1、b2在研究时段内呈明显视线向沉降趋势，平均形变速率分别为−3.8 cm/yr和−2.2 cm/yr。由图8可知，b1、b2在2015年8月至2015年11月及2016年1月至2016年3月存在较大程度的持续沉降。至2016年8月，b1和b2累积沉降量约分别达到了5.5 cm和4.5 cm。
Figure 7. Google Earth satellite image of East Anshan Iron Mine (area B) (a) and surface deformation rate of East Anshan Iron Mine (area B) from SBAS-InSAR analysis (b)
图7. 东鞍山铁矿(B区) Google Earth卫星影像(a)和基于SBAS-InSAR获取的东鞍山铁矿(B区)地表形变速率(b)
Figure 8. Time series accumulative subsidence of points b1 and b2 in area B from SBAS-InSAR analysis
图8. SBAS-InSAR获取的2个B区特征点b1、b2时序沉降量
重点形变区C(以下简称C区)主要覆盖鞍山市铁西区，范围内以工业区为主，并包括鞍山西站所在区域，其Google Earth卫星影像与形变速率结果如图9所示，形变速率范围在−3.8 cm/yr~3.6 cm/yr之间。C区内，东北侧工业区整体呈明显沉降，且越靠东侧形变越为明显。结合卫星影像分析，C区形变情况可能与工业区和鞍山西站及其周边工程建设密切相关。
在2015年5月14日至2016年8月12日期间，C区沉降漏斗中心位置c累积沉降量为4.46 cm (图10)，平均形变速率约为−3.8 cm/yr。
Figure 9. Google Earth satellite image of Tiexi District (area C) (a) and surface deformation rate of Tiexi District (area C) from SBAS-InSAR analysis (b)
图9. 铁西区(C区) Google Earth卫星影像(a)和基于SBAS-InSAR获取的铁西区(C区)地表形变速率(b)
Figure 10. Time series accumulative subsidence of the center of the subsidence funnel c in area C
图10. C区沉降漏斗中心位置c时序沉降量
通过对A、B、C区域的分析以及相关研究可知，矿区开采会导致矿区及其周边地区地表发生显著沉降。此外，过度开采会对矿区地质结构造成破坏，进而加速矿区的地表沉降，并可能导致泥石流、滑坡和崩塌等地质灾害的发生。随着矿区地表沉降面积的扩大，农业及建设用地将会显著减少，从而进一步制约矿区的经济发展，并可能对人民生命财产安全带来威胁。因此，使用时序InSAR技术对矿区及其周边地区的地表沉降情况进行长时段持续监测是极其必要的。
此外，本文与覆盖鞍山地区的已有形变监测结果[1-2, 8-9]相比，整体形变分布在地理范围上基本一致，但形变程度存在一定差异。从2000年至2010年和2016年至2018年，鞍山地区地表沉降呈加速趋势；而从2010年至2016年和2018年至2021年，地表沉降则呈减缓趋势。由此可推测，在鞍山地区地表沉降加速阶段，铁矿开采作业较为频繁，且开采程度较大；而在地表沉降减缓阶段，铁矿开采频率较弱，开采程度也较小。与此同时，以上对比也证明了本文使用高分辨率COSMO-SkyMed影像对露天矿区进行地表精细形变监测的方法具备可行性与可靠性。
5. 结论
本文选取2015年5月至2016年8月覆盖鞍山市的37景COSMO-SkyMed影像，基于SBAS-InSAR技术对该区域的地表形变进行了监测与分析。监测结果表明，研究区的地表形变主要集中在露天开采矿区，鞍山市区局部地方有沉降信号，研究区形变速率范围在−6.0 cm/yr~5.6 cm/yr之间。
此外，本文选取了齐大山铁矿(A区)、东鞍山铁矿(B区)和铁西区(C区)及其周边区域作为重点形变区进行了讨论。A区形变整体呈沉降趋势，且在铁矿内部及铁矿边界沿线较为显著，速率范围在−6.0 cm/yr~2.9 cm/yr之间。且通过推测可知，在2015年12月至2016年1月，A区南侧特征点a2周边地表可能受采矿工期调整或其它地质因素等原因影响，导致了地表持续沉降态势中断。B区形变与A区类似，整体呈沉降趋势，速率范围在−3.8 cm/yr~4.7 cm/yr之间，在铁矿内部及周边均存在以沉降为主的形变，且铁矿外东南侧存在显著的LOS向沉降信号。C区东北侧工业区整体呈沉降趋势，且有明显沉降漏斗，最大沉降速率约为−3.8 cm/yr。
结合相关研究可知，本文结果符合矿区开采会造成所在地区地表沉降的规律，证明了利用高分辨率COSMO-SkyMed影像对露天矿区进行地表精细形变监测的可行性与可靠性。同时也进一步验证了利用SBAS-InSAR等时序InSAR技术对矿区及其周边地区地表进行长时段持续监测的必要性。
NOTES
*通讯作者。