基于InSAR技术的大同市云冈矿区地面沉降监测

黄德华

(福建省地质测绘院,福建 福州 350011)

开采地下矿产资源容易破坏岩层结构,可能引发地面沉降等一系列地质问题,从而对沉陷区的地表建筑、公共设施、耕地、水环境等人类生存环境、经济财产安全造成极大的影响[1-4]。为了区域的可持续发展,科学合理地进行开发、生态治理,对矿区长时间序列的动态监测十分必要。

由于矿区形变周期性长、范围广,而常规的监测手段具有监测点离散、周期长、成本高的特点,难以进行时空连续性监测。合成孔径雷达差分干涉测量技术是通过短周期内的影像对干涉相位的差分计算,可大范围观测地表形变,其监测成本低,精度可达厘米至毫米级[5]。2000年CARNEC C等[6]首次将D-InSAR技术用于煤矿区地表形变场监测;2005年吴立新等[7]对中国东部典型工矿区进行了矿区地表沉陷D-InSAR监测研究。在D-InSAR技术的基础上,2001年FERRETTI A等[8]提出了PS-InSAR技术;2002年BERARDINO P等[9]随后提出了SBAS-InSAR技术,此后更是发展了众多时序监测算法。近年来,多国学者利用InSAR技术对矿区地面开展了形变监测,例如:2015年GRZOVIC M等[10]使用PS-InSAR与SBAS-InSAR方法对美国伊利诺伊州斯普林菲尔德地下煤矿采区进行了地表形变监测,监测结果显示地下矿井坍塌导致了该区域的地表形变;2023年虎小强等[11]使用SBAS-InSAR方法对新疆阿希矿区地面沉降时空变化特征进行监测,分析了沉降主要驱动因素,等等。现有研究表明,SBAS-InSAR与PS-InSAR是目前比较有代表性的时序InSAR算法,可有效监测矿区地面沉陷。

本文阐述了InSAR技术原理,运用D-InSAR、时序InSAR技术对云冈矿区开展地面沉降监测,其中在2019年1月—2020年12月时段采用24景哨兵-1A卫星影像数据进行SBAS-InSAR处理获取形变信息,在2020年1月—2020年9月时段采用25景Sentinel影像数据进行PS-InSAR处理获取形变信息,并选取2020年1月—2020年9月两种技术方法不同形变信息进行对比分析,验证两种时序InSAR技术在矿区地表沉降监测应用中的可靠性。

1 研究区及数据概况

1.1 研究区概况

研究区位于山西省北部、大同市西部云冈区内,面积约137.16 km2,西北部为煤矿开采区,79%以上为地下开采区,其余为露天开采区。东南部为大同市城市区域(图1)。

(a).大同市范围;(b).研究区范围图1 研究区概况图Fig. 1 The overview of the study area

1.2 实验数据

本文采用2019年1月—2020年12月的哨兵-1A卫星数据及SRTM 30 m分辨率数字高程数据作为数据支撑,其中SBAS-InSAR和PS-InSAR采用数据情况详见表1。

表1 哨兵-1A卫星影像主要信息表Table 1 Main information of Sentinel-1A satellite image

2 数据处理原理及流程

InSAR技术中相位信息由平地相位、轨道误差、地形相位、形变相位、大气相位以及噪声相位共6部分相位信息组成[12],公式如下:

φdiff=φtopo+φdef+φorb+φatm+φnoise,

式中:φtopo为参考数字高程模型(Digital Elevationg Model, DEM)残余地形相位;φdef为形变相位,由线性形变及非线性形变组成;φorb为去平引入的轨道误差相位,空间上表征为一定的线性特性,可包含于φatm中;φatm为大气延迟相位;φnoise为噪声相位。

对形变前后影像进行干涉处理,引入参考DEM模拟地形相位,去平、降噪、减弱大气延迟后,剔除模拟地形相位,即可获取形变相位从而生成差分干涉图,反演一维地表视向形变信息[13]。

相比D-InSAR技术的短期形变场获取,多时相合成孔径雷达干涉测量技术(Multi-temporal Synthetic Aperture Radar Interferometry, MT-InSAR)增大影像数量,通过空间域、时间域滤波,缓解大气延迟效应等误差,能对长时间序列下的影像集进行综合分析,得到地表长期演变趋势[14]。永久散射体技术(PS-InSAR)及小基线集技术(SBAS-InSAR)是2个目前比较有代表性的MT-InSAR算法[15]。

2.1 SBAS-InSAR

SBAS-InSAR技术将所有覆盖同一地区的SAR影像组成若干个子集,子集内的影像基线距(包括时间基线距和空间基线距)小,子集间的影像基线距大,运用奇异值分解方法获取最小形变速率,基于最小形变速率标准获得小基线干涉图。SBAS-InSAR方法限制了长基线导致的几何去相干,而且使更多的SAR图像参与到形变计算,增加了时间上的采样。提出SBAS方法是为解决影像数量不足、基线失相干无法进行相位分离的问题,基于一定的空间时间基线的若干干涉对集进行后续相位解缠等操作,从而产生解缠相位图用以分析地表形变情况,以减少时间和空间去相干因素的影响,可以得到监测区域的平均形变速率[16]。

SBAS-InSAR主要处理流程(图2)包括:估算影像对基线,选择配准的主影像,完成时序影像配准;按照短空间基线组合原则,生成短基线干涉对数据集,进行常规差分干涉步骤,去除平地、地形相位,获取差分干涉图集;利用相干系数图,选取高相干点,进行相位解缠;通过相干性设定阈值进行分布式目标分析提取;采用奇异值分解法模型求解形变参数和高程误差值即残余地形的最小范数解;进行高通、低通滤波进行第二次反演,计算非线性形变和大气相位并进行分离,从而反演形变总和,将结果进行地理编码生成时序形变信息[17]。

图2 SBAS-InSAR方法基本流程图Fig. 2 Basic flowchart of SBAS-InSAR

2.2 PS-InSAR

PS-InSAR技术选取影像集采用长时间序列下各个分辨单元内部保持稳定、散射特性强烈的永久散射体,利用相位信息可靠的特征,对这些离散点进行时间序列分析,从而提取目标形变量,探测地表形变信息[18]。

PS-InSAR技术(图3)利用K+1幅SAR单视复数影像集,结合外部DEM,计算得到K幅差分干涉图,并通过幅度离差系数法等算法提取稳定的N个PS候选点[19]。

图3 PS-InSAR方法基本流程图Fig. 3 Basic flowchart of PS-InSAR

根据大气相位时序随机、空间相关的特点,构建差分网络,对相邻的两PS点干涉相位进行3次差分,以降低大气噪声相位影响。利用二维周期图估计高程以及形变参数,根据求解结果进行相位解缠,得到的相位残差中包括残余大气信号、非线性形变信号。由于大气相位空间域上低频信号、时间域上高频信号,进行空间域低通滤波、时间域高通滤波,将大气相位与非线性形变相位分离,将线性形变相位与非线性形变相位相加即为视向形变相位,最终反演形变场[20]。

3 结果分析

3.1 D-InSAR沉降结果分析

对2020年6月2日—2020年7月8日的2景哨兵-1A影像进行D-InSAR处理,生成差分干涉图(图4)。结果表明,研究区西北部荣华皂村、栗庄村、永定庄村、里南沟村等地可见明显沉降漏斗,最大视向沉降量达4.6 cm/36 d,沉降区域普遍形变量为2 ～3 cm/36 d;东南部城市区域地表稳定,无明显形变。

图4 D-InSAR监测研究区沉降结果图Fig. 4 Subsident results of D-InSAR monitoring area

3.2 SBAS-InSAR沉降结果分析

2019年1月—2020年12月的24景哨兵-1A影像进行SBAS-InSAR处理,获得平均形变速率图(图5)。研究区西北部存在12处明显沉降,分布范围基本与D-InSAR结果重合,且由外缘至中心形变量逐渐增大,沉降区最大平均沉降速率达180 mm/y,其余区域平均沉降速率为5～0 cm/y,2019—2020年最大累积沉降量为333 mm,存在沉降中心超出形变测量梯度的情况;东南部城市区域地表稳定,无明显形变。

图5 SBAS-InSAR技术监测研究区平均形变速率图Fig. 5 Average deformation rate of SBAS-InSAR monitoring area

根据所属地理位置将沉降漏斗划分为荣华皂村沉降区、晋华宫街道沉降区、忻州窑村沉降区、曹家窑村沉降区、里南沟村沉降区以及永定庄村沉降区6个沉降区域,其中里南沟村沉降区沉降现象最为严重,沉降区域面积达1.27 km2。

选取里南沟北西向(C—D)、北东向(A—B)剖面做沉降漏斗分析,得到其空间累积沉降分布图(图6)。A—B剖面(图6(b))相对完整,即累积沉降剖面线由北部边缘-中心-南部边缘呈浅-深-浅漏斗状,该剖面沉降中心随时间变化基本维持不变,各处沉降量随时间变化而增加,总体不断下沉;C—D剖面(图6(c))不完整,边缘至中心沉降逐步加深,由累积沉降剖面可知,2019年沉降中心尚可监测,随着时间推移,该剖面沉降中心逐步向东转移,预计未来沉降中心沉降值将超出探测梯度。

(a).里南沟村沉降区剖面位置图;(b). A—B剖面示意图;(c).C—D剖面示意图图6 2019—2020年研究区累积沉降剖面示意图Fig. 6 Schematic diagram of cumulative subsidence profile from 2019 to 2020 in the study area

3.3 PS-InSAR沉降结果分析

对2020年1月—9月的25景哨兵-1A影像进行PS-InSAR处理,获得平均形变速率图(图7)。研究区西北部存在12处明显沉降,分布范围基本与SBAS-InSAR结果重合,从形变分布上可以看出,两者所获取的形变区域空间分布特征一致,沉降漏斗位置分布吻合,且形变趋势及累积形变量较为一致。

图7 PS-InSAR技术监测研究区平均形变速率图Fig. 7 Average deformation rate of PS-InSAR monitoring area

3.4 SBAS与PS技术结果交叉验证

在3个明显沉降区,选取里南沟村沉降区边缘至沉降中心的特征点进行时序对比分析,采用2020年的SBAS-InSAR监测结果,并提取PS-InSAR中对应时相的监测结果,将2020年1月设为起始月份,绘制特征点2020年1月—2020年9月的累积形变曲线(图8)。

(a). 里南沟村沉降区;(b). A-1点形变量;(c). A-2点形变量;(d). A-3点形变量图8 里南沟村沉降区域3个特征点累积形变图Fig. 8 Cumulative deformation of three characteristic points in area of Linangou Village

由里南沟村沉降区域3个特征点累积形变曲线图(图8)可看出,SBAS-InSAR与PS-InSAR监测结果总体趋势一致,表明2020年以来,该区域持续沉陷,其中A-1特征点SBAS-InSAR技术监测累积形变量最大达35.93 mm,PS-InSAR技术监测累积形变量最大达41.99 mm;A-2特征点SBAS-InSAR技术监测累积形变量最大达49.88 mm,PS-InSAR技术监测累积形变量最大达69.17 mm;A-3特征点SBAS-InSAR技术监测累积形变量最大达57.71 mm,PS-InSAR技术监测累积形变量最大达70.99 mm。用均方根误差(RMSE)对A-1、A-2、A-3点的形变量观测值xsbas,t和xps,t进行比较:

计算得A-1点RMSE值为2.89,A-2点RMSE值为7.48,A-3点RMSE值为5.31。

对两组形变量观测数据做线性拟合(图9),A-1观测点的决定系数R2为0.947 2,A-2观测点的决定系数R2为0.960 1,A-3观测点的决定系数R2为0.962 6。

图9 里南沟村沉降区域3个特征点两组形变量观测结果拟合曲线图Fig. 9 Fitted curve of the two observation groups of shape variables at three characteristic points in area of Linangou Village

由上可知,SBAS-InSAR形变监测与PS-InSAR形变监测结果相近,形变区域吻合,形变趋势基本一致,决定系数均>0.9,验证了SBAS与PS两种时序InSAR技术在区域形变监测具有一致性。

4 结论

大同市云冈区作为全国最大的产煤区之一,其地下煤矿开采区引起地面沉降灾害。本文应用InSAR技术对区域的多景哨兵-1A雷达影像进行处理,获取了地面沉降特征,详细分析了时空演变规律并进行了初步探讨,主要结论如下:

(1) 沉降主要分布于云冈区西北部煤矿开采区,沉降面积大,沉降速率达180 mm/y,最大累积沉降量为333 mm;东南部城市区域没有明显形变迹象。

(2) 采用SBAS-InSAR与PS-InSAR技术对其地表沉降情况进行长时间序列监测,监测结果表明形变空间分布吻合,时序沉降趋势一致,形变速率具有较高的相关性(R2决定系数均>0.9),均能有效监测地表沉降变化。