基于SBAS_InSAR的新建广湛铁路采空区形变监测及分析

基于SBAS_InSAR的新建广湛铁路采空区形变监测及分析

杨修强   

中铁二院工程集团有限责任公司 四川成都610031

摘要：

长期的矿产资源开采在佛山市高明区留下了许多采空区，采空区的存在严重威胁广湛线建设的开展。本文以广湛线采空区以及周边区域作为研究区，采用2015年至2019年的118景Sentinel-1A影像数据，利用短基线集(Small Baseline Subset Interferometric Synthetic Aperture Radar, SBAS-InSAR)技术对采空区及周边区域地面沉降进行监测，以获取该区地面沉降速率和累积沉降量，探究该区域地面沉降状况，为广湛线建设的顺利开展提供支持。研究发现采空区年形变速率最大为-6cm/年。此外在周边区域存在多个沉降漏斗，其中位于广湛线南侧的E、F、G、H点三年形变累积量都达到了10cm，位于广湛线北侧B点和C点近三年累积形变量也接近8 cm和6cm，且累积形变量逐年增大。结果表明为了保证广湛线的顺利开展应加强对该范围的深入勘察，并加大对D点周边山体稳定性的巡视力度，及时对不稳定区域其展开治理工作。

0引言

近年来，中国高速铁路工程建设高速发展，高速铁路里程和铁路网密度持续增加，因此不可避免的会遇到各种复杂的地质情况。其中，采空区作为一种在中国广泛分布的不良地质，容易诱发地裂缝、岩溶塌陷、地面沉降等地质灾害，对铁路工程建设及运营安全具有严重的威胁[1, 2]。因此对铁路建设采空区及周边区域开展实时的地表形变监测具有重要意义。

新建广湛铁路采空区位于珠江三角洲西翼，佛山市高明区。高明区矿产资源丰富，已发行矿产高达14种，是珠江和西江交汇的重要节点。前期在该区域的矿产资源开发留下了许多采空区，这些采空区的存在为该区域广湛线的建设带来了难度。因此有必要对该区域新建铁路附近的采空区地表形变进行时序监测。

传统的大地测量方法主要利用水准测量、全球导航卫星系统(Global Navigation Satellite System, GNSS)等测量手段，这些常规的方法存在观测范围小、周期长、 效率低、时空分辨率不高、耗费大量人力物力等缺点，很难适应大范围地表形变监测[3]。新发展起来的InSAR技术可以获取毫米级大范围地表形变，具有时空分辨率高，全天时、全天候、受云雾雨等恶劣天气的影响较小等优点，已经广泛的被用于大范围地表形变监测工作中[4-6]。小基线子集(Small-Baseline Subset, SBAS)技术，通过使用短时空基线的差分干涉对(DInSAR结果)实现对地表形变的时序解算[7, 8]。这种方法通过提高目标的干涉相干性实现对更多观测点的获取，进而得到更多、更全面的形变信息。SBAS-InSAR技术已成功应用于滑坡、火山、冻土等时序形变监测中[9-11]。

因此本文基于SBAS-InSAR技术，利用2015年6月-2019年12月共118景Sentinel 1A卫星数据对佛山市高明区新建广湛铁路采空区地面沉降进行监测，获取该区地面沉降速率和累积沉降量，进而揭示新建广湛铁路采空区地面沉降状况。

1 SBAS-InSAR原理

短基线干涉测量技术(Small Baseline Subset Interferometric Synthetic Aperture Radar, SBAS-InSAR)是由Berardino等[7]于2002年提出的一种用于监测地面微小形变的时间序列分析方法。SBAS-InSAR方法利用差分干涉图的适当组合，其特征在于设置恰当的时空基线阈值，利用小基线限制空间去相关现象。在对时序干涉对解缠相位进行处理时，SBAS-InSAR在假定时序形变是缓慢线性变化的前提下，通过对时序形变相 位速率和 DEM 高程误差构建线性解算模型，并采用最小二乘的方法对时序形变和 DEM 高程误差进行估计。若DInSAR时序差分干涉对存在时域上不连续的子集(即存在单幅 SAR 影像仅参与一次干涉对的构建)，求解方程则会出现秩亏，此时可采取奇异值分解 (SVD)的方法代替最小二乘进行求解。其基本流程图具体如下：

图1 SBAS-InSAR数据处理流程概略图

根据上述有关SBAS-InSAR形变建模与解算的描述，在进行DInSAR处理时，N+1景SAR影像可组合的差分干涉对数量M为：

(1)

任意差分干涉对的差分相位可表示为：

(2)

式中和分别表示雷达坐标系为(x，r)的像素点在tB时刻与tA时刻形变相位，为雷达波长，和分别表示雷达坐标系为(x，r)的像素点在tB时刻与tA时刻形变量。由此我们可以对所有时序差分干涉对的形变相位构建观测方程为：

(3)

表示M个差分干涉对的形变相位数据集，表示其它N幅SAR影像相对于参考影像的累积形变相位，表示示每个差分干涉对残余误差相位，A表示维系数矩阵，具体形式如下：

(4)

式中，-1代表主影像，1代表副影像，0表示SAR影像之间不存在连接组合关系。如果在理想条件下，所有差分干涉对在时间上是连续的，矩阵B的秩为N，那么形变相位利用最小二乘法可表示为：

(5)

当差分干涉对不连续时，观测方程存在秩亏的情况，利用奇异值分解方法(Singular Value Decomposition, SVD)求解式(5)的最小范数最小二乘解，矩阵A的SVD分解可以表示为：

(6)

其中U和V分别表示M阶和N阶酉矩阵，S是一个M×N阶对角矩阵，其对角线元素为矩阵A的奇异值，于是形变相位的最小范数最小二乘解可以表示为：

(7)

最后，可利用最小二乘解算的方法求解出各未知参数，进而得到扣除DEM高程残差的时序形变相位信息。

2研究区域概况和数据

研究区位于佛山市高明区，高明区西、北部以低山、丘陵为主，东部和东北部为广阔的冲积平原。形成西、北面环山，西南向东北走向由高至低的狭长地形。该区域矿产资源丰富，矿产资源丰富，已发现的矿种达14种，富含花岗岩、石灰石、铁、锰、金、银等。长期的矿产开采造成该区域存在多个踩空区。如图2所示，新建广湛线经过区域存在大范围采空区，给广湛线的建设带来了诸多不确定因素。为了保证广湛线在该区域的顺利施工，有必要对该区域开展时序监测。

图2 研究区域概况图

本文的数据采用欧空局网站发布的，2015年6月-2019年12月共118景C波段升轨影像。Sentinel-1A影像的空间分辨为5×20米，时间采样间隔为12天，Sentinel-1A影像的基本参数如表1所示。采用AWSD30产品作为外部DEM以去除地形效应，其高程精度为5米，空间分辨率为30米[12]。

表1 Sentinel-1A IW模式影像基本参数

3结果和分析

本文采用的是2015年6月至2019年12月共118景Sentinel-1 A影像数据，由于该卫星在2017年3月份进行了轨道调整，为了影像数据的高精度配准及减少变轨带来的误差影响，本文以2017年3月12日为分界点将118景影像分成两个阶段，采用SBAS-DInSAR时序形变解算方法，设定垂直基线阈值为186米，时间基线阈值为48天，分别获得85个和227个差分干涉组合，并分别得到了各干涉组合的差分结果。本文所采用的SBAS-DInSAR方法影像集垂直基线和时间基线图如图3-1、3-2所示。

图3 干涉网络时空基线图

3.1基于SBAS-DInSAR的广湛线采空区形变

截止2019年12月，课题组选用了从2015年6月-2019年12月共118景Sentinel 1A卫星数据，采用SBAS-DInSAR时序形变解算方法，通过设定时间和空间阈值以保证相干散射点散射信号的稳定性和形变结算结果的高可靠性，进而获得了广湛线在近三年的形变速率和时序形变曲线结果。图4-1为广湛线及周边区域2015-2017年形变速率图，图4-2为广湛线及周边区域2017-2019年形变速率图，其中正值代表沿LOS方向靠近卫星运动，负值代表沿LOS方向远离卫星运动。由年形变速率图可以看出，由于植被覆盖密集，以及检测区域内存在水体等，使得用SBAS-DInSAR时序形变解算方法解算出的效果不是特别好，很大一部分区域都没有得到有效点，从而无法检测出其形变。

图4-1 广湛线及周边区域2015-2017年形变速率图

图4-2 广湛线及周边区域2017-2019年形变速率图

3.2基于Stacking的广湛线采空区形变

为了更好地得到该区率的形变情况，所以又采用了Stacking的进行解算。图5-1和图5-2为用Stacking方法进行解算的广湛线在2015-2017和2017-2019年的形变速率和时序形变速率图。通过图5-1和图5-2所示，可以看到在监测区域内有几个明显的沉降区域。

图5-1 广湛线及周边区域2015-2017年形变速率图

图5-2 广湛线及周边区域2017-2019年形变速率图

为了更进一步的对广湛线形变进行详细深入的分析，项目组在线路周边分别选取了2015-2017时间段和2017-2019时间段各8个具有代表性的特征点（空间分辨率约30米，即代表30米半径范围内的平均观测值），并对这些特征点的时序形变曲线进行提取，以更好的描述山体的形变特征。这些点的位置如表2和3所示。山体及周边区域2017-2019年形变特征点选取位置如图6所示，其中图(a)表示2015-2017年时序形变选点图，图(b)表示2017-2019年时序形变选点图。各特征点2015-2017时序形变曲线如图7所示，其中图(a)至图(e)表示1号点至5号点的时序形变曲线图。各特征点2017-2019年时序形变曲线如图8所示，其中图(a)至图(h)表示A点至H点的时序形变曲线图。

表2 2015-2017年时序形变特征点经纬度坐标

表3 2017-2019年时序形变特征点经纬度坐标

图6广湛线及周边区域时序形变选点图。

图7 2015-2017年1号点至5号时序形变曲线图

从各特征点时序形变曲线图可以看出，2015年6月至2017年3月期间，距离广湛线D1K59+800左侧约530米的1号点形变量较大，达到了9 cm。其次是位于广湛线D1K59+000左侧560米的2号点和D1K58+800左侧约320m的3号点都发生了4 cm的累积形变，4、5号点的形变量较小，处于稳定状态。

图8 2017-2019年A点至H点时序形变曲线图

2017年3月至2019年11月期间， A点位于广湛线北侧，距离广湛线约5.4公里，其近三年累积形变量达到了16 cm。B点位于广湛线北侧，距离广湛线4.5公里，近三年累积形变量接近8 cm。C点位于广湛线北侧，距离广湛线0.5公里，近三年累积形变量接近6 cm。D点位于广湛线D1K59+800左侧400米附近，近三年累积形变量达到了14 cm，且呈现出形变加速的趋势，根据谷歌光学影像判断是采掘区域，所以其形变对线路的施工会带来隐患。E点位于线路南侧约2公里，近三年累积形变量超过了10 cm。F点位于线路南侧约1公里处，近三年累积形变量超过了10 cm。G点位于线路南侧约4公里处，近三年累积形变量超过10 cm。H点位于线路南侧约5.7公里处，近三年累积形变量超过12cm。

4总结

本文重点针对广湛线的形变情况进行监测分析，目采用雷达差分干涉测量技术，获取山体地表高时空分辨率、高精度、高可靠性形变信息，反映山体边坡的地表形变时空特征规律，建立该区域地表形变监测的基础数据，以便为长期性开展监测工作打下基础。

本项目监测时段为2015年6月15日-2019年12月3日，对该区域近五年间的地表形变演化趋势进行了详细的分析，基本得出以下几项结论：

广湛线及其周边区域存在自东北向西南分布的四处大形变区域，在广湛线点D1K59+800南侧400米附近区域（D点和1号点）存在着较大形变，年形变速率最大达到了-6 cm/y。在离广湛线北侧约5.4公里的A处存在较大的形变，年形变速率达到了-7 cm/y，形变范围沿东北-西南大约1千米，而在广湛线南侧也有明显的形变区域，在广湛线南偏4.0公里、5.7公里的G、H的两处可以看到明显的形变区域，其年形变速率都大约在-6cm/y的范围内，其形变范围大概都在沿东北-西南大约0.5千米。这几个形变区域离广湛线较远，对线路的稳定性造成影响很小。

（2） 针对时序形变监测区域近五年的形变特征进行分析，并通过卫星影像图人工目视解译结果发现，在疑似采矿工厂上方的山体表现为极不稳定，呈逐年下沉的趋势。特别是位于广湛线D1K59+800南侧400-500米附近区域，2015至2017年累积形变量最大超过了9 cm，2015至2019年近三年累积形变量最大超过了12 cm，且呈现出形变加速的趋势。其它位于广湛线南侧的E、F、G、H点三年形变累积量都达到了10cm，位于广湛线北侧B点和C点近三年累积形变量也接近8 cm和6cm，且累积形变量逐年增大。因此，建议加强对该范围的深入勘察，并加大对D点周边山体稳定性的巡视力度，及时对不稳定区域其展开治理工作。

参考文献

[1]马传广, 李翔, 高文峰. InSAR技术在高速铁路采空区地质选线中的应用研究 [J]. 铁道标准设计, 2017, 61(05): 23-7.

[2]闫明. InSAR技术在采空区铁路选线中的应用 [J]. 铁道勘察, 2021, 47(01): 18-22.

[3]侯建国, 张勤, 杨成生. InSAR技术及其在地质灾害中的应用 [J]. 测绘与空间地理信息, 2007, (06): 28-30+5.

[4]胡乐银, 张景发, 商晓青. SBAS-InSAR技术原理及其在地壳形变监测中的应用 [J]. 地壳构造与地壳应力文集, 2010, (00): 82-9.

[5]余景波, 王建波, 王肖露, et al. 基于三轨法D-InSAR技术的矿区地面形变分析 [J]. 华北地震科学, 2012, 30(04): 7-13.

[6]Wang H, Wright T J, Yu Y, et al. InSAR reveals coastal subsidence in the Pearl River Delta, China [J]. Geophysical Journal International, 2012, 191(3).

[7]Berardino P, Fornaro G, Lanari R, et al. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms [J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(11): 2375-83.

[8]Lanari R, Mora O, Manunta M, et al. A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms [J]. Geoscience & Remote Sensing IEEE Transactions on, 2004, 42(7): 1377-86.

[9]Casu F, Elefante S, Imperatore P, et al. SBAS-DInSAR Parallel Processing for Deformation Time-Series Computation [J]. IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing, 2014, 7(8): 3285-96.

[10]Hu B, Li Z. Time-Series InSAR Technology for Ascending and Descending Orbital Images to Monitor Surface Deformation of the Metro Network in Chengdu [J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14: 12583-97.

[11]Barbouchi M, Abdelfattah R, Chokmani K, et al. Soil Salinity Characterization Using Polarimetric InSAR Coherence: Case Studies in Tunisia and Morocco [J]. IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing, 2015, 8(8): 3823-32.

[12]Takaku J, Tadono T, Doutsu M, et al. UPDATES OF 'AW3D30' ALOS GLOBAL DIGITAL SURFACE MODEL IN ANTARCTICA WITH OTHER OPEN ACCESS DATASETS [J]. 2021.