SBAS-InSAR和高分辨率光学遥感相结合的矿区地表形变时空演化特征研究

耿 丹，涂 宽，谌 华，郑 健

（二十一世纪空间技术应用股份有限公司，北京 100096）

引言

矿产开采沉陷所引发的地表大面积沉降会给矿区及周边环境带来严重的负面影响。对矿区地表形变进行长时间序列监测，分析其时空演化特征，可为矿区地面沉降、矿洞坍塌、泥石流、滑坡等地质灾害防治提供有效支撑，减少对矿区周边人民生命财产及工程设施的威胁，提高矿区土地的使用价值，促进矿区的可持续发展［1-2］。

高分辨率光学及雷达遥感技术已广泛应用于地表形变监测、地质灾害监测、信息提取等方面［3-6］，其中合成孔径雷达差分干涉测量（Differentialinterfer Interferometric Synthetic Aperture Radar，D⁃InSAR）技术能够快速获取大面积的整体地表形变信息，且监测成本低，能克服各种天气进行全天时监测，在地表形变监测方面已开展广泛应用［7-9］，但是该方法受时间、空间失相干和大气延时相位的影响，难以在长时间序列的地表缓慢变形监测中得到理想结果［10-14］。小基线集差分干涉（Small Baseline Subsets Interferometric Synthetic Ap⁃erture Radar，SBAS⁃InSAR）同时具有D⁃InSAR技术的优点，而且受垂直基线的影响较小，能有效减小空间失相干和大气相位影响，监测长时间大范围地表缓慢形变，监测有效性及精度较高，同时能够获得研究区的形变时间演化特征，在地表形变监测方面应用广泛［15-16］。孙晓鹏等［17］利用SBAS⁃InSAR方法和ENVISAT ASAR数据，监测成都平原2008年到2010年的地表形变；彭米米等［18］采用SBAS⁃InSAR方法和哨兵1A数据监测大西安地区2015年到2017年的地表形变并对时序特征进行研究；华怡颖等［12］以大宁矿区为例，利用小基线集技术对研究区9景PALSAR数据进行时序处理，得到研究区沉降速率图、时序累计沉降图。上述学者都利用SBAS⁃InSAR技术对地表沉降进行监测和时序分析，但却不能有效的对地表沉降的空间演化特征进行分析。

标准差椭圆（Standard Deviational Ellipse，SDE）是分析地理要素空间分布特征的经典方法之一，已广泛应用于经济、人口、农业、城市发展等各行各业。近年来，已有学者将标准差椭圆方法应用于自然灾害的空间分布特征研究。如：熊俊楠等［19-20］利用标准差椭圆等方法对重庆市1950⁃2015年、西藏自治区1983⁃2015年山洪灾害进行时空分布特征分析；周超凡等［21］利用标准差椭圆方法对2003⁃2010年北京市地面沉降进行空间演化特征研究。

综上所示，为研究矿区地表沉降时空演化特征，本文选取阜新某矿区为试验区，采用SBAS⁃InSAR技术提取矿区长时间序列地表形变信息，基于形变剖面提取典型特征点分析矿区沉降时间演化特征，利用标准差椭圆方法分析矿区形变空间分布演化特征，结合高分辨率光学遥感影像对监测结果进行验证，为矿区的地质灾害防治提供参考依据。

1 试验区及数据源介绍

1.1 试验区介绍

本文选取阜新市某矿区作为试验区（图1），该试验区地势由东南向西北倾斜，海拔100～621 m，中部地势较低。矿区地层由三部分构成：上部主要为砂岩、砂砾岩；中部是含有少量砂岩、砾岩和砂质页岩等的煤层；下部基本是细砂岩［22］。根据已有调查成果，矿区有多处矸石山和采坑，矸石山与采坑相互连接，整体沿北东-南西向呈条带状展布［23］。

图1 试验区范围Fig.1 Scope of test area

1.2 数据源介绍

本文采用由欧空局免费提供的Sentinel-1B星Level-1级的单视复数图像（Single Look Complex，SLC）作为数据源［24-26］，共获取2019年10月14日-2020年3月30日期间的15景影像数据。同时采用2020年12月获取的北京二号0.8 m高分辨率光学遥感数据产品进行专家判识。雷达及光学卫星相关参数见表1。为提高地面形变结果精度，选用SRTM 30m分辨率DEM数据来减少地形因素对SBAS⁃InSAR监测结果的影响。

表1 影像相关参数Table 1 Relevant parameters of image

2 研究方法

本文以Sentinel-1B SAR影像和北京二号高分辨率光学遥感影像为数据源，基于SBAS⁃InSAR技术和标准差椭圆方法开展矿区地表形变时空演化特征研究。首先提取矿区形变时间序列和矿区年均形变速率，并基于GIS的空间分析功能对矿区地表形变的时间和空间演化特征进行分析，最后利用北京二号光学影像对监测结果进行验证。研究技术路线见图2。

图2 矿区地表形变时空演化特征研究技术路线Fig.2 Technical route of research on temporal and spatial evolution characteristics of land subsidence in mining area

2.1 SBAS-InSAR技术

SBAS⁃InSAR小基线集差分干涉技术将所有雷达影像数据组合成多个集合，集合内数据之间的时间和空间基线距尽量小，再解算每个小集合的地表形变时间序列［27⁃38］。主要步骤包括：（1）按照时间顺序排列研究区的雷达单复数影像数据，设定时空基线阈值，进行差分干涉；（2）选取一定数量的高相干点进行轨道精炼和重去平；（3）在稳定区域选择GCP控制点，对差分干涉影像进行相位解缠；（4）利用SVD方法解算解缠相位，求解形变速率和高程系数，去除大气相位和地形残余相位；（5）地理编码得到WGS84坐标系下的矿区年平均形变速率、形变时间序列信息。

2.2 标准差椭圆方法

1926年Lefever提出了标准差椭圆的方法，其基本思想是以椭圆的中心、面积、长轴、短轴、方位角来分析矿区形变的中心、范围、形状、方向等演化过程和趋势［29⁃31］。以矿区地表形变量为权重，求解标准差椭圆的中心、长轴、短轴和方位角的计算公式如下［32］：

（1）标准差椭圆的中心

式中，（SX,S Y）即为标准差椭圆的中心，（x i，y i）是时序点i的坐标是所有时序点平均中心，n是时序点的数量。

（2）标准差椭圆的方位角

式中，θ即为标准差椭圆的方位角和是时序点坐标（x i，y i）与所有时序点平均中心的偏差。

（3）标准差椭圆长轴和短轴的标准差

式中，σx和σy分别为沿长轴和短轴的标准差。

3 结果分析与验证

3.1 SBAS-InSAR处理结果

本文利用SBAS⁃InSAR技术对15景Sentinel-1B数据处理，对影像进行多视、滤波、相位解缠、干涉对挑选等，最终获取86个干涉对，像对的时空基线见图3，部分差分干涉像对见图4。

图3 干涉像对时空基线连接图Fig.3 Time space baseline connection diagram of interference image pair

图4 差分干涉像对Fig.4 Differential interference image pair

矿区2019年10月-2020年3月的地表形变年平均速率见图5。整个监测期内，矿区整体存在不同程度的形变，地表不均匀形变严重，形成了多个形变中心。矿区整体形变沿北东-南西方向呈条带状分布，与开采坑口空间展布方位一致。矿区内形变严重区域形变速率达到50～60 mm/a，主要分布在形变的中心位置。

3.2 矿区形变时序演化特征分析

矿区整体存在不同程度的形变，部分区域形变量较大，造成在形变中心及周边区域失相干现象严重，在图5中表现为无值区，无法分析及确定形变中心位置。为分析矿区整体形变时序演化，利用Kriging插值方法对2019年10月-2020年3月的累计形变量进行空间插值分析。Kriging插值是一种能够准确产生预测表面的地统计方法，在土壤科学和地质学中尤为适用。

图5 矿区年平均速率分布图Fig.5 Distribution of annual average rate in mining area

基于矿区地表累计形变分布图及插值图，根据形变中心位置（图5中无值区）沿矿区的北东-南西方向作一条纵向贯穿整个矿区的剖面线AA’，尽可能多的通过分布在该方向上的形变中心；依据以上原则，再沿北西-南东方向作四条垂直于剖面线AA’的剖面线BB’、CC’、DD’、EE’，如图6黑色线所示（CC’和DD’剖面线之间区域由于形变中心聚集在AA’剖面线附近，因此不再绘制剖面线进行分析），并得到五条剖面线的形变剖面图（图6），进而对整个矿区的地表形变时序形变特征进行分析。

图6 矿区地表累计形变量插值图Fig.6 Interpolation map of accumulated surface settlement in mining area

纵向贯穿矿区的剖面线AA’形成了7个沉降漏斗，见图7（a），分别对应图6中的特征形变点A1-A7点。分析这7个特征形变点的地表形变时间序列变化（图7（b））发现，特征点的地表形变量在监测期内呈现出逐渐增大的趋势，且形变量基本都到达150 mm以上，最大累计形变量为217 mm，位于A5号特征点。随着各特征点形变范围的不断扩大，相邻沉降漏斗很可能会合并，如A3和A4特征点。贯穿矿区的剖面线BB’形成了2个沉降漏斗，见图7（c），分别对应图5中的特征形变点B1、B2。由这2个特征形变点的形变时间序列变化图（图7（d））发现，监测期内B1和B2都呈现出逐渐增大的趋势，且B2点从2020年1月开始形变速度加快，超过B1点，最终形变量达到197 mm。剖面线CC’形成了3个沉降漏斗，见图7（e），分别对应图5中的特征形变点C1、C2、C3。由这3个特征形变点的形变时间序列变化图（图7（f））发现，监测期内3个特征点都呈现出逐渐增大的趋势，且C2点从2019年12月开始形变速度明显加快，最终形变量达到212 mm。剖面线DD’形成了5个沉降漏斗，见图7（g），分别对应图5中的特征形变点D1-D5，由这5个特征形变点的形变时间序列变化图（图7（h））发现，监测期内特征点都呈现出逐渐增大的趋势，且D1点从2020年3月18日形变量超过D5特征点，达到205 mm。剖面线EE’形成了2个沉降漏斗见图7（i），分别对应图5中的特征形变点E1、E2，由图7（j）可以看出监测期内特征点都呈现出逐渐增大的趋势，且形变速率基本一致，形变量均超过200 mm。

图7 矿区剖面图及剖面上特征点累计形变量Fig.7 Mining area profile and cumulative settlement of characteristic points on the profile

3.3 矿区形变空间演化特征分析

采用标准差椭圆方法进一步分析矿区形变的空间分布及形变演化趋势，以地表形变量为权重，计算标准差椭圆来分析地表形变随时间推移的空间演变趋势。

以2019年10月-2020年3月各月份的地表累计形变图制作矿区各月份地表形变的标准差椭圆，见图8。从2019年10月到2020年3月，形变量在持续增大，且在矿区不均匀分布。本文从标准差椭圆的圆心变化、面积、长短轴变化以及方位角变化等几个角度来定量分析2019年10月-2020年3月期间矿区地表形变的空间分布特征和发展趋势。2019年10月-2020年3月标准差椭圆参数数值见表2，各参数变化趋势见图9。

图8 2019年10月-2020年3月矿区地表形变空间演化图Fig.8 Spatial evolution of land subsidence in mining area from October 2019 to March 2020

表2 2019年10月-2020年3月标准差椭圆参数列表Table 2 Standard deviation ellipse parameter list from October 2019 to March 2020

图9 标准差椭圆参数变化图Fig.9 Variation of standard deviation ellipse parameters

（1）矿区形变中心变化分析

从矿区2019年10月-2020年3月年各月份的累计形变标准差椭圆的中心经纬度以及迁移轨迹（图10），2019年10月至2020年3月，标准差椭圆的中心逐渐向南西方向移动，矿区形变中心在这6个月内整体向南西方向移动了505.62 m。

图10 标准差椭圆中心位置变化图Fig.10 Variation of center position of standard deviation ellipse

（2）矿区形变空间范围变化分析

标准差椭圆的面积及长轴变化反应矿区形变空间范围的变化。矿区2019年10月-2020年3月累计形变量的标准差椭圆面积统计结果及长轴长度见表2，椭圆面积变化图及长轴变化见图9（a）（b）。2019年10月-11月，标准差椭圆的面积减少，对应的长轴长度也相应减少，表明在此期间矿区地表形变的范围在减少；但从2019年11月-2020年3月，标准差椭圆的面积由8.06 km2增大到9.32 km2，长轴长度由2.21 km增大到2.64 km，表明在此期间矿区形变面积在逐渐扩大，且长轴方向的形变也在逐渐增大。

（3）矿区形变空间形状变化分析

标准差椭圆的短轴比长轴值的变化可以表示矿区形变空间分布形状的变化。比值越接近于1，表明矿区地表形变在各个方向演化较为均匀；而比值越接近于0，表明矿区地表形变的方向性分布特征越明显。矿区2019年10月-2020年3月的地表累计形变标准差椭圆短轴与长轴的比值见图9（c）。2019年10月到同年11月，短轴比长轴值变大，说明矿区地表形变在北东-南西方向形变减缓，而在北西-南东方向形变加剧；而在2019年11月-2020年3月，短轴与长轴比值逐步下降，矿区地表形变空间分布形状表现出扁化趋势，表明地表形变在长轴方向（北东-南西）上逐步发展，而在短轴方向上（北西-南东）发展相对减缓。

（4）矿区形变空间方向变化分析

标准差椭圆的方位角变化反应矿区地表形变空间发展的主方向变化。矿区2019年10月-2020年3月的地表累计形变标准差椭圆的方位角变化见图9（d）。2019年10月-2020年3月，矿区地表形变空间分布方位角呈逐步增大趋势，方位角由2019年10月的16.7°增大到2020年3月的21.2°。为更直观地表达地表形变空间分布方位角的变化趋势，在图11中以图的形式来示意各期次标准差椭圆角度的增量变化，可以看出各期次标准差椭圆在空间上表现为逆时针旋转，2019年10月到2020年3月方位角增大了4.5°。

图11 准差椭圆方位角变化趋势示意图Fig.11 Variation trend of azimuth angle of standard deviation ellipse

3.4 结果验证

利用2019年12月获取的北京二号0.8 m高分辨率光学遥感影像，叠加矿区地表形变速率见图12，可以看出：矿区整体存在不同程度的形变，部分区域形变量较大，造成在形变中心及周边区域失相干现象严重，在图11中表现为无值区。对图中9个地表形变中心分析（图13），发现形变严重区域在高分辨率光学影像上均表现为已开采矿坑或正在进行开采区域，地表形变中心位置与矿区开采位置一致。

图12 矿区地表形变中心分布图Fig.12 Mine location map

图13 SABS⁃InSAR结果与光学影像对比图Fig.13 Comparison of SABS InSAR results and optical images

另外，以沉降中心位置1为例，利用Google earth上与监测期相近的3期光学影像（2019年8月、2019年9月和2020年6月）进行矿区开采活动验证，如图14。从3期影像的红线范围内可以看出，该位置开采程度随时间推移逐渐加剧。

图14 开采位置1在Google earth中不同时相影像对比图Fig.14 Comparison of different time phase images of mining location 1 in Google Earth

本次试验，利用0.8 m高分辨率的北京二号光学遥感影像和多期次Google earth影像同时验证了利用SBAS-InSAR进行矿区地表形变监测的准确性和技术的可行性。

4 结论

本文利用15景Sentinel-1B干涉宽幅雷达降轨数据及0.8 m高分辨率光学影像，基于SBAS⁃InSAR获取阜新某矿区2019年10月14日-2020年3月30日地表形变信息，并对矿区地表形变的时空演化特征进行分析，利用同期北京二号高分辨率光学影像对SBAS⁃InSAR监测结果进行验证，得到以下结论：

（1）SBAS⁃InSAR监测结果表明，矿区整体存在明显的不均匀地表形变现象，最大地表形变速率达到50-60 mm/a之间，主要分布在矿区形变的中心位置。

（2）对贯穿矿区的5条剖面及剖面线上形变特征点的时序形变分析，共识别出矿区19个沉降漏斗，位置分散，最大形变量达到217 mm。随着特征点形变范围的不断扩大，相邻特征点对应的沉降漏斗在后期很可能会合并。从剖面图也可以看出，随着剖面线方向的推进，形变量在漏斗边缘地区较低，而在距离形变中心较近的区域形变量在不断增大，这与矿区开采的沉降规律一致。

（3）基于标准差椭圆方法对试验区累计形变标准差椭圆的中心、面积、形状和方向4个方面的变化分析发现，2019年10月-2020年3月，试验区地表形变的形变中心由北东向南西方向迁移，形变的空间范围也在逐步扩大；且矿区地面形变在长轴（北东-南西）逐步发展，短轴（北西-南东）发展相对减缓，导致椭圆形状出现扁化趋势。

（4）结合北京二号0.8 m高分辨率光学影像专家判识，发现矿区形变严重区域为已开采或者正在开采区域，验证了利用SBAS⁃InSAR进行矿区地表形变监测的准确性和技术的可行性。

综上所述，本文利用SBAS⁃InSAR技术对试验区15景Sentinel⁃1B降轨数据进行处理，获取由煤矿开采活动引起的地表形变信息和典型形变特征点，定量分析了矿区地表形变的时空演化特征，结合同期高分辨率光学遥感影像对监测结果进行验证，为矿区的地质灾害防治提供了参考依据。

由于本次研究缺乏监测期内的地面实测数据，不能有效检验SBAS⁃InSAR监测结果的精度，后续将重点收集试验区地面实测数据，进一步完善研究成果。