基于时序InSAR技术的焦作地区地表形变监测研究

Study on Surface Deformation Monitoring in Jiaozuo Area Based on Time Series InSAR Technology

摘要: 焦作地区地面沉降显著，为了准确掌握地面沉降动态，利用时序InSAR技术进行研究区的地表形变监测。文中使用短基线集干涉测量(SBAS-InSAR)技术对2020年9月~2021年11月的20景Sentinel-1A SAR影像进行干涉测量处理，从而生成形变时间序列；然后用相同时期的水准数据对地面沉降速率结果数据进行验证，通过进行参考基准的统一和整体偏差补偿，得到了经补偿以后的地面沉降速率数据，再对地面沉降速率数据与水准测量数据两者之间进行线性回归，比较的结果是二者趋势一致，线性关系明显，相关系数为0.703762，量值吻合。

Abstract: The land subsidence in Jiaozuo area is remarkable. In order to accurately grasp the land subsidence dynamics, the time series InSAR technology is used to monitor the land deformation in the study area. In this paper, 20 Sentinel-1A SAR images from September 2020 to November 2021 are processed by short baseline interferometry (SBAS-InSAR) technology, so as to generate deformation time series. Then, the leveling data of the same period are used to verify the rate of land subsidence result data, and the compensated rate of land subsidence data are obtained by unifying the reference datum and compensating the overall deviation. Then, the rate of land subsidence data and leveling data are linearly regressed, and the comparison results show that the two trends are consistent, the linear relationship is obvious, the correlation coefficient is 0.703762, and the values are consistent.

1. 引言

传统的监测方法如水准测量、全球卫星导航系统(GNSS)测量，这些方法精度很高，但其针对小范围离散点观测，不仅观测周期长，费时费力，且难以在短时间内快速获取大范围面状地表形变信息[1] [2]。时序InSAR技术的提出和发展有效地解决了这些影响，可以获得毫米级的地表形变监测精度[3]，在对长期累积的缓慢形变监测方面表现出了极大的潜力[4]-[6]。目前，以永久散射体干涉测量(PS-InSAR)技术和短基线集干涉测量(SBAS-InSAR)为代表的时序InSAR技术得到了广泛应用[7]。不少学者在此方面开展了不少的研究和应用。许兵等[8]以山西大同矿区为例，获取了2022年12月23日至2023年5月20日共25景LT-1条带模式影像数据，分别进行SBAS-InSAR技术和PS-InSAR技术数据处理；范军等[9]使用PS-InSAR和SBAS-INSAR两种技术相结合进行了地面沉降监测的对比分析；鲁魏等[10]利用SBAS-InSAR和PS-InSAR对52景升轨Sentinel-1A雷达影像，分别获取了研究区2017年6月至2020年11月的地表形变速率及时序形变量；董继红等[11]利用时序InSAR技术(Stacking技术、小基线集技术)开展滑坡隐患识别对比分析，通过对比基于Sentinel-1数据利用不同时序InSAR技术获取地表形变速率图，发现SBAS结果受各种误差影响较小，监测效果较Stacking结果更好；侯靖钥等[12]基于小基线集(SBAS)技术对24景Sentinel-1A影像进行时间序列处理，得到安源区2018~2019年间地面沉降监测结果，并讨论分析了影响地面沉降的主要因素；刘辉等[13]采用2016~2018年的53景Sentinel-1数据，基于永久散射雷达干涉测量(PS-InSAR)和小基线集(SBAS)技术获取秦皇岛段明长城地表形变信息。刘媛媛等[14]利用时序InSAR技术实现大范围地形变监测，有效揭示了区域性地表形变的空间分布特征，对于宏观揭示盆地尺度的地表形变的空间变化特征具有重要意义。

文章使用时序InSAR技术，对研究区的20景Sentinel-1A SAR影像，进行地表形变信息提取，对提取的结果通过同一时期的水准数据进行对比验证，即进行参考基准的统一与整体偏差补偿，得到了经补偿以后的地面沉降速率数据，最后再分析地面沉降速率数据与水准数据之间的线性关系。

2. 研究区概况和数据来源

2.1. 研究区概况

研究区以焦作市马村区地面沉降带为研究区。研究区东至焦作武陟县龙泉街道，西至焦作孟州市赵和镇，南至焦作孟州市化工镇，北至焦作马村区安阳城街道，总面积1497 km²，地理坐标：东经112˚45'~113˚26'，北纬34˚53'~35˚19'，研究区区内公路、铁路等基础交通设置建设水平良好，不仅实现了对外的高速畅通，同时辖区下的镇村同样具有优良的道路交通环境。见图1地理位置图。

Figure 1. Geographical location map of the study area

图1. 研究区地理位置图

2.2. 数据来源

时序InSAR选取2020年9月11日至2021年11月29日时间区间内20景Sentinel-1A升轨数据，数据产品为SLC，地面分辨率为5 × 20 m。

轨道数据Precise Orbit Ephemerides (POD精密定轨星历数据)在接收到GNSS数据的21天后产生的，为Sentinel-1A数据最精确的精密定轨星历数据文件。

另外，用实地水准观测数据作为参考数据。研究区监测网分布在焦作市的武陟县、孟州市、沁阳市、温县、博爱县；项目组于2021年10月~12月组织完成了2021年度孟州武陟二等水准测量工作，共完成二等水准测量471 km。根据水准测量规范规定，二等水准测量每千米测量的偶然中误差不大于1.0 mm，每千米测量的全中误差不大于2.0 mm。

3. 数据处理方法

3.1. 时序InSAR技术

时序InSAR技术是近年发展较快的技术。该方法不需要提前假定形变模型，而通过三维相位解缠算法获取地表形变信息[15]。时序InSAR技术通过增加干涉对的时间采样率和多普勒频率来增加干涉图相干性，提高干涉图的精度[16]。研究中使用的时序InSAR技术主要为SBAS-InSAR技术。

3.2. 数据处理

SBAS-InSAR数据处理，经过对Sentinel-1A数据分别进行数据导入，数据裁剪，连接图生成，干涉工作流，轨道精炼和重去平及地理编码等步骤[17]，从而获取了研究区的地表形变信息。

Figure 2. Spatiotemporal baseline map of all images

图2. 所有影像的时空基线图

Figure 3. Time baseline map of all images

图3. 所有影像的时间基线图

为了提高结果精度，本次将空间基线阀值设为临界基线的2%，时间基线阀值设为120天，并相应地设置空间基线阀值，进行3D解缠。

得到所有影像的时空基线图，见图2；所有影像的时间基线图，见图3。干涉像对79对，超级主影像为2021年1月9日；4πλ = 226.56。

Figure 4. Coherence coefficient diagram

图4. 相干系数图

干涉处理方面，采用了可以提高干涉条纹清晰度的Goldstein滤波方法，对图像进行滤波处理。干涉流生成的一系列结果，包括相干系数图，见图4；去平后的干涉图，见图5；去平和滤波后的干涉图，见图6；解缠结果图，见图7。

通过地理编码后，获取了时序InSAR雷达视线方向(LOS)的地表形变量图，见图8。

Figure 5. Interference pattern after flattening

图5. 去平后的干涉图

Figure 6. Interferogram after flattening and filtering

图6. 去平和滤波后的干涉图

Figure 7. Results of unwinding

图7. 解缠结果图

Figure 8. SBAS-InSAR rate diagram in the study area

图8. 研究区SBAS-InSAR速率图

4. 分析与结果

4.1. 整体偏差补偿

研究区选择地面沉降速率结果为待验证数据，以2021年的水准观测数据作为参考数据，对地面沉降速率结果进行参考基准的统一和整体偏差补偿。

研究区水准观测点分布较均匀，对研究区进行整体偏差补偿，参考基准统一步骤如下：

1. 坐标转换

利用水准观测值对InSAR形变速率进行参考基准校正，使其统一于同一坐标系下。根据己获取实地观测结果，提取InSAR结果中对应的观测值，比较二者之间的相关关系。

由于InSAR测量的形变速率图进行地理编码时，得到的速率图是WGS_84坐标系；而工作区获取的2021年72个水准观测点为CGCS2000国家坐标。现将InSAR测量的形变速率图投影转换为CGCS2000国家坐标，使其位于同一坐标系下。

2. 回归分析

将InSAR测量值结果与对应的水准点值进行对照比较，见图9，根据两组样本数据的关系，利用线性回归，得到线性关系图，见图10。表1为回归方程式及相对应的相关系数、标准误差、平均误差和中误差。回归方程y = a + bx式中a即为水准与InSAR测量的补偿值，而b的大小反映水准与InSAR测量结果相关性，b值越接近1则两者之间相关性越强。因此，在剔除残差过大的点之后，利用线性关系求解2组数据之间的整体偏差是可行的。

Table 1. Fitting equation of time series InSAR and leveling observation results

表1. 时序InSAR与水准观测结果拟合方程

拟合方程	标准误差(mm)	相关系数	平均误差(mm)	中误差(mm)
y = 0.6665x + 11.529	4.532393	0.536373	±10.25	±11.9245

Figure 9. Comparison of time series InSAR and leveling measurement results

图9. 时序InSAR与水准测量结果比较

Figure 10. Linear regression of time series InSAR and leveling

图10. 时序InSAR与水准测量线性回归

4.2. 补偿后精度评估

比较回归方程的相关系数和方差，可见该组数据线性关系明显，因而采用线性模型进行多项式拟合，即满足条件：

$\{\begin{array}{l}\widehat{y}=a+\widehat{b}x\\ \left|y-\widehat{y}\right|<\epsilon \end{array}$ (1)

式(1)中：x为InSAR测量年均沉降量；$\widehat{\text{y}}$为线性模型拟合值；y为水准测量值。依据该区域地质条件及《DD2014-11地面沉降干涉雷达数据处理技术规程》等相关规范，考虑到图像配准以及InSAR测量值插值处理所引起的误差，取(方差) ε = 20 mm为阀值，比较水准观测值与InSAR拟合结果，剔除超过给定条件及数据处理时距GCP点较远(偏差过大)的异常点位的14个。最终得到线性拟合样本数据，利用该组数据，重新拟合线性方程。

$\sigma =\pm \sqrt{\frac{{{\displaystyle \sum }}_{\text{i}=1}^{\text{N}}{\left(\text{y}-\widehat{\text{y}}\right)}^2}{\text{N}-1}}$ (2)

由式(1)计算得水准观测值与线性回归模型补偿得到的InSAR监测结果值的二组数据的标准差σ = ±3.51382。检验结果表明两者测量结果波动小，即离散程度低，水准观测值与InSAR测量值线性关系明显。图11为经过补偿的样本点比较图，表2为补偿后时序InSAR与水准观测结果的重新拟合方程。

Table 2. Re-fitting equation of time series InSAR and leveling observation results after compensation

表2. 补偿后时序InSAR与水准观测结果重新拟合方程

拟合方程	标准误差(mm)	相关系数	平均误差(mm)	中误差(mm)
Y = 0.8614x + 2.8048	3.462775	0.703762	±2.959495	±3.605558

Figure 11. Comparison between InSAR and leveling results after compensation

图11. 补偿后时序InSAR与水准测量结果比较

4.3. 补偿后地面沉降速率图

利用表1对整景沉降速率图进行补偿，得到新的地面沉降InSAR监测速率图，见图13。以及线性回归图，见图12，可见二者趋势一致，量值吻合。图13即为研究区焦作地区的地面沉降InSAR监测速率图，根据这个结果，结合自然因素(包括构造活动、软弱土层的自重压密固结等)及人类活动因素(包括地下水开采、地热资源开发、石油开采、城市建设等方面)可以分析研究区引起地面沉降的原因。

Figure 12. Linear regression after compensation

图12. 补偿后线性回归

Figure 13. Rate of land subsidence map of the study area after compensation

图13. 补偿后研究区地面沉降速率图

5. 结论

1) 研究以焦作市马村区地面沉降带为研究区，通过时序InSAR方法对20景影像数据进行处理，获取了研究区的2020~2021年地表的时间序列沉降速率形变信息。

2) 研究使用反演结果数据与实测水准数据比较，进行线性回归分析，从而进一步证明时序InSAR技术在监测地表形变的可靠性。

3) 对研究区开展时间序列InSAR技术地面沉降监测研究，有着积极的研究意义。通过研究，分析出形成沉降的人类活动因素和自然因素，将为研究区所在的地区进行建设提供有力的参考价值，对研究区所在的地区发展起到促进作用。

基金项目

河南省地质局局管地灾防治项目(豫地灾防治[2025] 1号)。

NOTES

^*通讯作者。