EN中文

Journal of Geo-information Science >

2023 , Vol. 25 >Issue 12: 2402 - 2417

DOI: https://doi.org/10.12082/dqxxkx.2023.220779

High Precision Extraction of Surface Deformation Information Based on Principal Component Spatiotemporal Analysis and Time-series InSAR: Taking Xuzhou as an Example

  • CHEN Yu , 1, 2, * ,
  • CHEN Si 2 ,
  • LI Jie 2 ,
  • LI Huaizhan 1, 2 ,
  • GAO Yandong 1, 2 ,
  • WANG Yong 3 ,
  • DU Peijun 4
Expand
  • 1. Key Laboratory for Land Environment and Disaster Monitoring (Ministry of Natural Resources), China University of Mining and Technology, Xuzhou 221116, China
  • 2. School of Environment and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China
  • 3. Navigation and Positioning Center, Jiangsu Provincial Surveying and Mapping Engineering Institute, Nanjing 210013, China
  • 4. School of Geography and Ocean Science, Nanjing University, Nanjing 210023, China
*CHEN Yu, E-mail:

Received date: 2022-09-30

  Revised date: 2023-01-18

  Online published: 2023-12-05

Supported by

National Natural Science Foundation of China(42171312)

National Natural Science Foundation of China(42001409)

Jiangsu Geology & Mineral Exploration Bureau Science and Technology Plan Project(2021KY08)

Fold

Abstract

Urban areas often suffer from varying degrees of land surface deformation due to infrastructure construction and resources exploitation, which threatens the safety of residents' lives and property. So regular monitoring of urban surface deformation is of great significance for preventing related geological disasters. However, urban surface deformation has the characteristics of small-scale and continuous-slow change, it is necessary to process the error carefully in order to improve the monitoring accuracy. This paper proposes a high-precision surface deformation extraction method combining the principal component spatiotemporal analysis and time-series Interferometric Synthetic Aperture Radar (InSAR). Through the mining and analysis of time-series InSAR signals, a surface deformation model combined with polynomial functions is constructed to realize the hierarchical estimation of error and noise signals. Then the high-precision, small-scale surface deformation information is extracted. Taking Xuzhou, a typical city prone to geological disasters, as the research area, the results show that the proposed method can accurately separate the surface deformation information and error in the time-series InSAR signal, and the deformation monitoring accuracy is 10%~57% higher than other existing methods. The deformation rate from 2018 to 2022 is about -17~35 mm/a in Xuzhou, which is mainly distributed in the urban area, along the subway and in the old goaf. In recent 8 years, urban construction has continuously triggered local subsidence areas, the secondary deformation of the old goaf can last for more than 6 years, and the surface of several mining areas is still in an unstable state. The results can provide important technical support and decision support for high-precision monitoring of urban surface deformation and prevention of potential geological disasters.

Key words： principal components spatiotemporal analysis; time-series interferometric synthetic aperture radar technology; land surface deformation; high precision; Xuzhou; old goaf; urban construction; secondary deformation of closed mine

Cite this article

CHEN Yu , CHEN Si , LI Jie , LI Huaizhan , GAO Yandong , WANG Yong , DU Peijun . High Precision Extraction of Surface Deformation Information Based on Principal Component Spatiotemporal Analysis and Time-series InSAR: Taking Xuzhou as an Example[J]. Journal of Geo-information Science, 2023 , 25(12) : 2402 -2417 . DOI: 10.12082/dqxxkx.2023.220779

1 引言

地表形变是在自然或人为活动的作用下，地表因受力不均匀形成的升降、倾斜和错动等现象[1]。城市地区常因基础设施建设、地下水和矿产资源开采等产生不同程度的地表形变情况，严重时会引发地裂缝、塌陷等地质灾害，威胁着城市居民的生命与财产安全。截至2021年，全球约有1 596个城市存在地表沉降现象[2]。我国华北平原2006—2010年形成了以城市为中心向外扩展的大范围沉降状态，以天津为沉降中心的平均沉降速率高达 -64.2 mm/a[3]，长江三角洲城市群2007—2010年形变速率在-27.4 ~4.9 mm/a之间，其中重大市政建设中存在明显漏斗状沉降现象[4]，墨西哥的墨西哥城2007—2011年由地下水开采导致的形变速率达-30 mm/a[5]、土耳其的伊斯坦布尔1992—2017年岩层因硬度较低和人为因素造成形变速率为-17.3 ~-3.8 mm/a[6]等。因此，对城市地区，特别是有地质灾害隐患的城市进行定期的地表形变监测十分必要。
传统测绘手段，如三角高程测量、水准测量、全球卫星导航系统（Global Navigation Satellite System, GNSS）技术等[7-8]，因其高精度的特点被广泛应用于地表形变监测领域，但其离散点监测的方式有时不能满足广域地表形变监测的需求。合成孔径雷达干涉测量[9]（Interferometric Synthetic Aperture Radar, InSAR）技术能全天时、全天候地获取高精度、高空间分辨率的地表形变信息，特别是近年来发展的时序InSAR技术，如永久散射体InSAR技术[10]（Permanent Scatterer InSAR，PS-InSAR）、短基线集干涉测量技术[11]（Small Baseline Subsets InSAR，SBAS-InSAR）和多时相InSAR技术[12]（Multitemporal InSAR, MT-InSAR）等，在保证空间分辨率的基础上可获取长时间序列地表形变信息，有助于刻画和揭示地表形变时空演变规律[13⇓-15]。不同于地震、火山活动引发的大梯度、高强度地表形变，城市地区地表形变具有量级较小和连续缓慢的特点，因此，对大气延迟、轨道残余等误差和噪声的精细估计与控制尤为重要。现有城市地表形变时序InSAR监测应用研究中多采用时空滤波方法去除大气延迟误差和噪声信号，但存在引入新噪声或形变信息被误判为噪声被去除等问题[16-17]。InSAR通用大气校正在线服务系统（Generic Atmospheric Correction Online Service for InSAR，GACOS）[18-19]也是当前大气延迟改正的主流方法，但对外部数据的依赖度较高。主成分分析（Principal Component Analysis, PCA）是一种常用的数据统计分析方法，在海量空间数据处理与信息挖掘中得到了广泛应用。Kositsky和Avouac[20]基于PCA提出一种GPS时序数据分析反演方法（Principal Component Analysis-based Inversion Method, PCAIM），并用于断层滑动模型反演。Lin等[21]采用PCAIM方法联合反演电子测距仪和InSAR数据，初步验证了该方法对InSAR数据的处理能力。Remy等[22]采用PCAIM联合处理火山区域InSAR和GPS数据，并将第一主成分用于数值模拟反演火山形变源。Chen等[23]基于PCAIM挖掘出火山侧翼的指数形变特征。与时空滤波方法及GACOS方法相比，PCA方法自身具有信号分析和分离能力，能够在不借助外部数据的情况下，基于SAR自身信号深度挖掘干涉信息时空分布规律，构建精确的时序形变模型，实现形变从干涉信号的预分离，并同时消除轨道、高程等误差，避免因形变与大气延迟信号耦合造成的过矫正。综合前人研究可知，当前PCA在InSAR领域的应用多面向大梯度形变监测及模型反演（如火山活动、地震等），较少聚焦于微小形变（如城市地区）的高精度提取和误差分析。因此本研究基于PCA方法挖掘时序InSAR信号，通过对主成分的时空分析，提取时序地表形变模型，进而构建城市地区高精度时序InSAR地表形变信息提取策略。
本研究选择江苏省徐州市作为研究区。徐州曾是我国重要的煤炭产地，煤炭资源的大规模开发产生了一系列由于矿区地表塌陷、采空区生态环境破坏等引发的生态环境问题。 《江苏省地质灾害防治“十四五”规划》明确指出徐州是地面沉降、地面及采空区塌陷、滑坡等地质灾害的易发地区，提出了“落实专业化监测，提升预警能力”的针对性要求。此外，依据现有研究， 2015—2018年，徐州城市轨道交通的持续建设导致沿线多个站点周边均产生了局部沉降[24]。因此定期监测徐州市地表形变并揭示其驱动机制对预防相关地质灾害、促进城市安全可持续发展具有重要意义。然而，目前对徐州地区2018年以来的地表形变监测，特别是高精度时空动态监测的相关研究工作较少。
因此，以徐州地区Sentinel-1A卫星SAR影像为主要数据源，利用GAMMA和StaMPS（Stanford Method for Persistant Scatterers Algorithm）进行时序InSAR处理分析。在此基础上，融合PCAIM软件包[20]和PCA时序InSAR分析方法[23-24]对时序信号进行分析，构建时序地表形变模型和误差估计策略，提出融合主成分时空分析与时序InSAR的高精度地表形变提取方法，获取徐州地区高精度地表形变信息，并将本方法与现有误差分析与改正方法进行对比分析。基于研究结果综合分析近年徐州市地表形变时空演变规律及其驱动因素，讨论研究方法的有效性和适用性，以期为区域地质灾害防治及城市地表形变高精度监测提供参考与理论支持。

2 研究区概况与数据来源

2.1 研究区概况

徐州市位于我国华北平原东南部、江苏省西北部、苏皖豫鲁四省交界处。研究区覆盖徐州市区范围（33.96°N—34.46°N，116.97°E—117.56°E）（图1）。地貌主要以冲积平原为主，海拔约20~50 m；局部有少量丘陵分布，海拔约100~300 m。徐州受中生代燕山运动影响强烈，形成弧形构造及纵向逆断层结构，有利于岩溶发育，导致岩溶洞穴与地下河流形成广泛，形成丰富的地下水资源[25]。同时，该地区的煤层也十分发育，煤炭开采活动与地下水开采是造成当地地质灾害的主要因素。近年来，为探索资源枯竭型城市的转型之路，徐州陆续关闭各处矿 井[26]，开展采煤塌陷地区的生态修复和综合利 用[27-28]。徐州西北矿区（图1区域1）和东北矿区（图1区域2）的生态修复获得显著成效并成为全国采煤塌陷地生态修复示范项目，其中潘安湖生态修复工程于2017年12月得到了习近平总书记亲临实地考察时的高度肯定[29]。2020年徐州市政府针对两区域生态修复现状，以潘安湖湿地公园为中心规划建设潘安湖新城，以九里湖为中心设立徐州淮海国际港务区，承担徐州老工业城市转型和发展的职责[30]。虽然随着区域治理与生态修复的推进，地下水开采受到严格限制，相关地质灾害发生显著减少，但是徐州市的地表稳定性问题仍然不容忽视。特别是近年徐州城市迅速扩张，建筑物体量日益增加，且自2014年以来开始建设城市地铁网络，包括目前已经建成并运行的1-3号线一期工程及正在建设的3号线二期和6号线一期工程。因此，影响徐州地表稳定性的因素多而复杂，定期开展徐州地区地表形变时空动态监测十分必要。
图1 徐州市特征地区地理位置及空间分布

Fig. 1 Geographical location and spatial distribution of the characteristic areas in Xuzhou

2.2 研究数据

影像数据为欧洲航天局提供的于2018年8月19日至2022年3月31日期间获取的66景C波段升轨哨兵1号A星（Sentinel-1A）SLC SAR影像[31]和对应的精轨数据，Sentinel-1A数据方位向和距离向分辨率分别为14 m和12 m，重访周期为12 d，轨道号142，极化方式VV，影像中心入射角为40.12°，SAR数据集时间序列信息见表1。利用美国地质调查局提供的30 m航天飞机雷达地形任务数字高程模型[32]（SRTM DEM）模拟与去除地形相位，并对干涉图地理编码。利用研究区内4个连续运行（卫星定位服务）参考站（Continuously Operating Reference Stations, CORS）站点数据（2018年4月2日—2021年12月31日）和5个水准站数据（2019年6月25日—2021年2月21日）验证InSAR地表形变监测结果的精度。
表1 徐州地区2018年8月19日—2022年3月31日SAR影像时间序列

Tab. 1 Time series of SAR images in Xuzhou from August 19, 2018 to March 31, 2022

序号 成像日期 序号 成像日期 序号 成像日期
1 2018.08.19 23 2019.12.24 45 2021.04.29
2 2018.08.31 24 2020.01.29 46 2021.05.11
3 2018.09.12 25 2020.02.22 47 2021.05.23
4 2018.10.06 26 2020.03.17 48 2021.06.28
5 2018.11.11 27 2020.04.22 49 2021.07.10
6 2018.12.05 28 2020.05.16 50 2021.08.03
7 2018.12.29 29 2020.06.09 51 2021.08.15
8 2019.01.10 30 2020.06.21 52 2021.09.08
9 2019.02.03 31 2020.07.15 53 2021.10.02
10 2019.02.27 32 2020.07.27 54 2021.10.14
11 2019.03.23 33 2020.08.20 55 2021.10.26
12 2019.04.16 34 2020.09.01 56 2021.11.07
13 2019.05.10 35 2020.09.25 57 2021.11.19
14 2019.06.15 36 2020.10.07 58 202112.01
15 2019.07.09 37 2020.10.31 59 2022.01.06
16 2019.07.21 38 2020.11.12 60 2022.01.18
17 2019.08.14 39 2020.12.06 61 2022.01.30
18 2019.08.26 40 2021.01.11 62 2022.02.11
19 2019.09.19 41 2021.01.23 63 2022.02.23
20 2019.10.01 42 2021.02.16 64 2022.03.07
21 2019.10.25 43 2021.03.12 65 2022.03.19
22 2019.11.06 44 2021.04.05 66 2022.03.31

3 研究方法

本研究针对城市地区形变量级较小，形变信号易受误差和噪声干扰，影响地表形变监测精度的 问题，构建融合主成分时空分析与时序InSAR的高精度地表形变提取策略。主要包括3个部分：时序InSAR处理、InSAR信号主成分时空分析和高精度地表形变信息提取。具体的研究技术路线如图2所示。
图2 融合主成分时空分析与时序InSAR的高精度地表形变信息提取技术路线

Fig. 2 Flow data processing chart of high precision extraction of surface deformation information based on principal component spatiotemporal analysis and time-series InSAR

3.1 时序InSAR处理

鉴于城市地区后向散射特性较稳定，且形变缓慢连续，采用的影像数据集时间间隔密集，能够保证干涉影像具有较高的相干性，因此选取StaMPS的SBAS方法[11]进行时间序列分析。首先，在GAMMA软件中使用精确轨道文件校正Sentinel-1A影像数据，并进行多视处理。利用SRTM DEM数据进行地理编码，经方位向和距离向频谱滤波后，以2020年1月29日采集的影像为超级主影像进行配准。其次，设定时间基线阈值<60 d，空间基线阈值<200 m，生成150对差分干涉对，时空基线情况如图3所示。然后，利用StaMPS软件，基于振幅离差与相位稳定性原则，筛选出地面共计216 739个高相干性点，通过3D解缠法[33]进行相位解缠，并估计和去除了由DEM、时空基线引起的相位误差，获得时序InSAR地表形变初始结果。
图3 短基线干涉对时空基线

Fig. 3 Perpendicular and temporal baselines of small baseline interferograms

3.2 InSAR信号主成分分析

为了进一步挖掘InSAR信号组分，运用PCAIM工具包[20]对时序InSAR初始结果进行PCA分解。首先中心化时序InSAR初始结果 X m × n = x 1 , x 2 , x 3 , ⋯ , x n，其中 m为像素数， n为时间序列影像幅数， x j为第 j幅影像相对第一幅影像高相干性点形变信息， j为时序编号， x j均值为 μ j j = 1,2 , ⋯ , n， x j ' = x j - μ j，得到中心化后的矩阵为 X m × n '。再通过SVD算法对 X m × n '初步分解如下：
X m × n ' = S m × m ∙ Σ m × n ∙ T n × n T
式中： S m × m和 T n × n为通过SVD解算后得到的标准正交矩阵； T n × n T为 T n × n的转置； Σ m × n中每一项都是该矩阵的奇异值。此时，时序信号被分解成若干个主成分的线性组合。
时序InSAR中心化后的结果数据降维至 k，生成主成分分解方案，自变量 k为主成分空间中主成分的数量。生成的主成分空间由空间矩阵（ S）、权重矩阵（ W）和时间矩阵（ T）组成。空间矩阵表示各主成分的空间分布情况（与干涉图的空间分布相对应），权重矩阵表示主成分分量的重要性指数，而时间矩阵描述了主成分分量的时间演变趋势。
X m × n ' ≅ S m × k ∙ W k × k ∙ T k × n T
式中： S m × k和 T k × n分别为空间矩阵和时间矩阵， W k × k为权重矩阵，其值表示每个分量的重要性， T k × n T为 T k × n的转置。因奇异值在矩阵中从大到小排列，因此可用前 k个奇异值近似描述矩阵，即式（1）和（2）中：
S m × m ≅ S m × k T n × n ≅ T k × n
时序InSAR信号的主要信息在主成分分解后被集中在若干主成分中。将初始主成分分量数目参数设定为3，依次分析分解后的主成分。本研究中，当分量数为3~5时主成分分量存在信号杂糅现象，导致主成分的时空规律难以识别；当分量数为7或更多时出现较多低重要性的冗余成分，且存在不同低重要性成分表现出相似时空特征的情况，主要表现为多主成分分量同时表现季节性非线性形变和噪声相位特征；当主成分分量个数为6时，各成分时空规律具有较明显的识别特征。如图4所示，第一主成分（PC1）与第二主成分（PC2）均呈明显的时间相关性，即整个时间段内， T 1和 T 2随时间呈线性变化趋势，观察其空间分布图（图4（a）—（b））可以看出，在研究区西北部和东北部以及市区均存在区别于背景的信号，这些信号与时序InSAR初始结果中的形变信号具有相似的空间格局，并且，此信号也与研究区地表形变空间分布的先验知识相吻合[34⇓⇓-37]，因此判定PC1和PC2为线性地表形变相关分量。第三主成分（PC3）未表现出明显的时间相关性，整个研究时段的 T 3值在零值上下波动，其 S 3 × W 3值具有一定的空间相关性，但与高程并无 明显相关性，可推测其主要为残余轨道或长波长误差。第四主成分（PC4）与第六主成分（PC6）随时间呈现出一定的季节性变化，其中PC4的 S 4 × W 4值与高程略成相关性，但在空间上未表现出明显的形变特征，因此可能由季节性大气扰动引起；PC6的空间分布图中存在形变特征，可能与季节性非线性形变相关。第五主成分（PC5）的时间函数在零值附近随机波动，空间分布未表现出明显的形变特征，可能为噪声相关分量。
图4 时序InSAR初始信号的PCA分解结果

Fig. 4 PCA decomposition results of the initial time-series InSAR signals

3.3 高精度时序地表形变提取

基于上述分析，采取如下误差估计策略： ① 使用线性函数拟合形变相关主成分（PC1和PC2）的时间演变规律，得到最优时间函数 T d '； ② 将得到的 T '与对应成分的空间分量（ S d）和重要性指标（ W d）联合，基于式（4）的重构时序地表形变模型（ D），并将时序地表形变模型从时序InSAR信号中预去除；③ 同理与上一步，对误差项（PC3和PC5）重构，并将其从时序InSAR信号中去除；④ 利用二次多项式函数（式（5））估计时间序列上每一干涉对残余信号中可能存在的残余大气、轨道等误差并将其去除；⑤ 将预先去除的时序形变模型加回到包含非线性形变的残余时序InSAR信号中，得到最终的时间序列地表形变结果。
D = S d ∙ W d ∙ T d ' T
式中： T d ' T为最优时间函数 T d '的转置； d为相应的主成分集合。
ϕ i , j = a x i + b y i + c x i y i + d x i 2 + e y i 2 + f z i + g z i 2 + h
式中： i为干涉图中的像素编号； j为时序编号； ϕ i , j为InSAR初始干涉相位； x i和 y i代表像素点坐标； z i为像素点的高程； a, b,…, h为拟合参数。

4 结果与分析

4.1 精度验证

为验证所提出方法获得的时序InSAR地表形变监测结果的精度，将其与研究区内4个CORS 站和5个水准站相近时期的观测数据进行对比， 获得地表形变散点图（图5）。由图5（a）可知，时序InSAR估算的地表形变结果与地面观测数据具有较强的一致性，相关系数达0.89，均方根误差（RMSE）为3.85 mm；图5（b）显示了时序InSAR初步结果和本研究方法得出的地表形变速率标准差分布情况，可见，改正后的形变速率标准差分布更为集中，整体向低值方向移动，均值由2.2 mm/a减少到0.9 mm/a，说明整体上降低了误差水平。综上，本研究方法所获取的时序地表形变结果具有较高的精度和可靠性。
图5 时序InSAR地表形变结果误差分析及精度验证

Fig. 5 Error analysis and accuracy verification diagram of time-series InSAR surface deformation results

分别采用目前应用较广泛的时空滤波方法和GACOS方法估计和改正3.1节获取的原始时序InSAR数据的误差。选择100 d与80 m×80 m的时空滤波窗口，利用大气相位空间低频和时间高频的特性，对干涉图直接估计大气相位并滤除。利用GACOS系统解算出相应时间和区域的高空间分辨率天顶对流层延迟改正图像，对所有差分干涉图进行大气延迟改正。将2种方法获取的改正结果与本研究结果进行对比分析。图6展示了研究区中3个子区域的原始差分干涉图及利用3种方法改正后的结果。由图可知，时空滤波方法和GACOS均能在一定程度上削弱误差和噪声影响，但在局部区域仍存在较明显的残余误差，或引入额外的噪声。相较之下，本研究方法能够更精准地去除残余误差及和噪声相位。
图6 不同误差改正方法效果比较

Fig. 6 Comparison of the results derived from different correction methods

表2统计了无形变区域的RMS及在4个CORS站点位置采用不同方法得到的形变结果与GNSS数据差值的RMS，及本方法相对其余方法精度的提高比率。由表2可知，时空滤波方法和GACOS均能使InSAR结果更逼近于GNSS形变数据，且两者效果相当，除0001号站对应像元外，其他GNSS站对应像元的RMS均被降低到10 mm以内。相较之下本文的方法进一步降低了研究区地表形变的误差，将地表形变计算精度提高了10%~57%。
表2 不同误差改正方法形变结果的误差分析

Tab. 2 Error analysis of deformation results from different correction methods

形变量/mm 相对时空滤波提高比率/% 相对GACOS法提高比率/%
原始数据 时空滤波 GACOS 本方法
无形变区域 2.9 2.7 2.8 1.4 48.2 50.0
DAWU 21.4 9.4 9.5 7.9 16.0 16.8
BNTG 9.9 7.1 7.4 4.8 32.4 35.1
0001 32.7 18.5 15.0 7.9 57.3 47.3
0003 16.4 6.4 6.8 5.8 9.4 14.7

4.2 徐州市地表形变速率综合分析

图7为徐州市时序InSAR地表形变平均速率图。结果表明，研究时段内，研究区地表形变具有覆盖范围大且空间分布不均匀的特点，形变速率范围为-17~35 mm/a。主要形变区域可划分为城区形变区、形变区域1和形变区域2，且主要集中在研究区北部。其中，城区形变区的地表形变主要分布在地铁沿线及城市建筑区，以地表沉降为主。最显著的沉降区域位于北部的祥和小区、孟家沟村和南部新城区地铁6号线的惠民家园站附近，沉降速率最高达15.76 mm/a。形变区域1和区域2为徐州市的2个老采空区，形变主要以地表抬升为主，伴随局部地面沉降。区域1内最显著的抬升分布在淮海国际港务区及其西部的夹河煤矿内，抬升速率达26.89 mm/a；区域2的抬升主要分布在潘安湖及其东部的权台煤矿范围，抬升速率达31.32 mm/a。
图7 研究区2018—2022年时序InSAR地表形变速率

Fig. 7 Displacement rate derived from InSAR time series from 2018 to 2022

4.3 城市建设区地表形变时空分析

城市扩张期间的建筑施工和地铁施工期间的隧道开挖是城市地面沉降的主要原因。结果表明，徐州城区地表形变幅度整体偏低，约在-17~5 mm/a区间。为详细分析地表沉降与地铁建设之间的关系，提取了地铁沿线1.5 km缓冲区范围的地表形变信息，如图8所示。除孟家沟村区域外（形变速率约为-15.76 mm/a），其他城区地表形变重点区均落在地铁沿线缓冲区范围内。地铁沿线缓冲区内地表沉降主要分布在地铁1号线人民广场站至徐医附院站区间、地铁2号线九里山站至庆云桥站区间、地铁3号线下淀站至后蟠桃村站区间和6号线紫金路站至奥体中心站区间。其中沉降最为显著的区域为2号线奔腾大道站附近的祥和小区（图8，P1）和6号线惠民家园站附近（图8，P2），形变速率分别约为 -13.67 mm/a和-14.54 mm/a。
图8 地铁沿线1.5 km缓冲区范围内2018—2022年地表形变时空变化

Fig. 8 Surface deformation within a 1.5 km-wide buffer along the subway lines from 2018 to 2022

4.4 老采空区地表形变时空分析

至2016年底，徐州市区内所有煤矿均已关闭。研究表明，矿井关闭后，随着排水系统设备停运、地下水位上升，会改变采动破裂岩体的应力和承载能力，导致覆岩与地表次生移动变形，可能引起严重的环境与地质灾害[38]。本节重点分析徐州市东北和西北矿区老采空区在2018—2022年的地表形变时空演变特征。
图7中的区域1为徐州西北矿区范围，主要包含张集煤矿、张小楼煤矿、庞庄煤矿和夹河煤矿（闭矿时间见表3）。经生态修复整治，该地区已经形成九里湖湿地保护区，并依托该保护区建成徐州淮海国际港务区。图9为港务区在研究时段内的地表形变结果，可见，形变主要发生在该区域南部。其中庞庄煤矿（P3）和夹河煤矿（P4）于2018年7月至2019年中处于平稳状态，之后开始抬升直至研究时间段结束，平均形变速率分别约为7.82 mm/a和26.89 mm/a。两矿区东部的拾西村曾被江苏省国土厅列为徐州市采空地面塌陷隐患点[39]，研究时间段内，该村的总体形变趋势自西北向东南由抬升转为沉降，其中西北部（P5）自研究期初始至2020年1月无明显形变而后开始抬升，平均形变速率为14.04 mm/a，东南部（P6）在研究期间由轻微沉降逐渐趋于平稳，平均形变速率为-2.73 mm/a。拾西村南部中欧尚郡为新建居民区，表现出明显的地面沉降（P7），形变速率约为-28.13 mm/a。除淮海国际港务区外，区域1还包括张集煤矿和张小楼煤矿两处明显形变区（图7）。张集煤矿形变趋势与 庞庄矿、夹河煤矿大致相同，其形变速率约为29.16 mm/a；张小楼煤矿主要形变区位于车村 （图10），主要表现为地面沉降，平均沉降速率约为-15.67 mm/a（图10，P8）。
表3 各矿区位置及关闭时间

Tab. 3 Location of each mine and their closure time

矿区名称 所属区域 关闭时间/年
庞庄煤矿 西北矿区 2012
夹河煤矿 西北矿区 2016
张小楼煤矿 西北矿区 2016
张集煤矿 西北矿区 2016
卧牛山煤矿 西部矿区 2008
新河煤矿 西部矿区 2008
权台煤矿 东北矿区 2014
旗山煤矿 东北矿区 2016
韩桥煤矿 东北矿区 2008
图9 徐州淮海国际港务区2018—2022年地表形变时空变化

Fig. 9 Surface deformation in the Xuzhou Huaihai international port affair area from 2018 to 2022

图10 张小楼煤矿2018—2022年地表形变时空变化

Fig. 10 Surface deformation in Zhangxiaolou Coal Mine from 2018 to 2022

图7中的区域2为徐州东北部矿区范围，主要包括韩桥煤矿、权台煤矿和旗山煤矿（闭矿时间见表3）。整体上该区域地表形变以抬升为主，抬升范围几乎覆盖整个潘安湖新城，影响范围约65 km2，最大形变速率达33 mm/a。权台煤矿自研究初期至2019年11月基本保持平稳，之后开始抬升，平均速率约为31.32 mm/a（图11，P11）；潘安湖新 城东部的蔡庄村为权台煤矿影响区域，同拾西村一起被列为徐州市采空地面塌陷隐患点[39]，在研究时段内持续抬升，其西南部平均形变速率约为24.33 mm/a（图11，P9），东北部抬升速率逐渐减缓，约为6.99 mm/a（图11，P10）。旗山煤矿和韩桥煤矿在研究时段内相对稳定，几乎无明显形变。
图11 潘安湖地区2018—2022年地表形变时空变化

Fig. 11 Surface deformation in Pan'an Lake area from 2018 to 2022

5 讨论

5.1 徐州市地表形变长时序变化趋势及驱动因素分析

Chen等[24]分析了徐州市2015—2018年地表形变时空变化规律，鉴于本研究时段是该研究的后续时段，因此综合2个研究做徐州市近8年长时序地表形变的综合变化趋势分析。如图12所示，南部的形变区在2015—2018年有明显沉降，平均形变速率-10 mm/a，而2018—2022年未表现出明显形变，说明该地区由于城市建设引起的地面沉降已趋于平稳。同时，祥和小区和孟家沟村由前期的平稳转为沉降期。据调查，祥和小区形变地区位于偏离地铁隧道中心轴约1.4 km的位置，其地表为20年以上的老建筑物，地基稳定性随时间推移降低，可能为局部地面沉降的诱因。孟家沟村距地铁隧道较远，但不排除受3号线二期建设影响，另外该区域目前处于棚户区一期工程改造中，此工程也可能是地表沉降的主要原因。近8年来，徐州地铁工程一直处于建设中。2015—2018年间徐州在建三条地铁 线路，地铁建设导致的地表形变主要位于地铁1号线工农路站附近区域（形变速率达到-30 mm/a）和 2号线汉源大道站区域（形变速率为-12 mm/a）；2018—2022年，1—3号线均已通车运营，1号线工农路站附近仍存在一定残余形变。分别选取工农路站和汉源大道站3个随机点位T1—T3、T4—T6做长时序点位分析，如图13（a）、图13（b）所示，沉降速率以封顶期为节点由高到低变化，并逐渐趋于平缓。可见，地铁在修建过程中的地表形变通常要比运营后更显著。另外，本研究期间新建的两条地铁线路中6号线惠民家园站为新增形变区，其他区域暂相对稳定。随着地铁的修建，预计仍会形成新的沉降区。
图12 2015—2018年和2018—2022年徐州地表形变对比

Fig. 12 Comparison of surface deformation in Xuzhou from 2015 to 2018 and from 2018 to 2022

图13 2015—2022年地铁沿线特征点及各矿区形变时间序列

Fig. 13 Time series of characteristic points in subway deformation and each mining area from 2015 to 2022

依据邓喀中等[38]提出的关闭矿井次生沉陷时间规律可知，矿井关闭后的地表形变分为5个阶段：① 初期稳定阶段、② 下沉阶段、③ 中期稳定阶段、④ 抬升阶段和⑤ 最终稳定阶段。地表下沉阶段主要由采动破裂岩体和覆岩结构的变形和失稳造成，而抬升阶段的主要诱因是采动破裂岩体和第四系松散层随着地下水位的升高而吸水回弹[35]。图13（c）为近8年研究区内主要矿区地表形变时间序列图。根据它们的时间变化规律可以看出，夹河、权台、旗山和张集煤矿（表3）均表现出沉降-平稳-抬升的趋势，结合次生沉陷规律可知，3个矿区均处于②③④阶段，且未见明显趋于最终稳定的趋势。虽然权台矿2014年已关闭，但为防止旗山矿突水事故发生，权台矿排水设施直至2016年旗山矿闭井之后才停止运行，因此与其他3个矿区的次生沉陷发生阶段较为同步。相较之下，由于韩桥煤矿（2008年关闭，目前处于⑤阶段）处于区域含煤盆地的盆心，各处煤矿的矿井水及地下径流通过各种渠道向韩桥煤矿汇集，而韩桥煤矿与旗山煤矿间存在F5断层（图11），很大程度上阻止了地下水流入旗山煤矿，加之权台煤矿排水系统帮助排水，导致旗山矿区抬升缓慢。庞庄煤矿关闭时间较早（2012年），虽有轻微的沉降和抬升现象，但从长时间尺度上呈平缓趋势。西部矿区的卧牛山矿和新河煤矿表现出抬升-平稳的趋势，说明其处于④⑤阶段，在经历地下水回升导致的地表抬升后进入最终稳定期。张小楼煤矿（2016年）呈现持续沉降趋势，沉降范围主要集中在车村，该村从2016年起至2020年底持续建设别墅群（图10），因此局部建设工程及建筑物持续载荷作用是车村沉降的重要驱动因素之一。
根据徐州近8年内地表形变趋势可知，城市建设会引起局部的地面沉降，并且会在工程结束后持续一段时间，因此，在当前密集的城市建设背景下，徐州市未来几年仍会出现新的沉降区，且多集中在居民密集区域，因此建设期间及建设完成一段时间内需对施工及其附近区域进行高频率地表形变监测，关注沉降区及附近路基、建筑物地基安全。另外，研究发现关闭矿井的次生形变持续了6年以上，未来几年夹河、权台、旗山和张集煤矿仍有持续抬升的趋势，张小楼煤矿仍有继续沉降的趋势，根据矿井次生形变发展规律，未来需重点关注以上矿区。

5.2 研究方法适用性分析

城市地区地表形变具有量级小、变化连续、进展缓慢的特点，对地表形变误差的精度要求较高。本方法通过对徐州市时序InSAR信号深层分析，精细化分层估计误差相位，提取了高精度时序地表形变信息，具有较高的可靠性。但是，方法的普适性仍需要进一步研究与讨论，这主要依赖于相应研究区地表形变的特点、地形起伏情况、受大气等误差干扰的程度及对地表形变监测精度的需求。如Chen等[23-24]将同类方法（原理与本方法有相似之处，但方法策略不同）应用于火山区域地表形变监测中，研究发现该方法对高程起伏大、形变与大气延迟信号耦合严重的情况具有优于其他多数方法的监测效果，但对于大气湍流引起的局部干扰现象的削弱效果较为有限。另外PCA是基于统计学的信号分析方法，若要明确如何分析和建立各主成分与形变或误差组分之间的关系，则需要在处理过程中反复试验和解译。因此，本方法数据处理的整体思路可以应用于其它领域地表形变监测的研究，但细化的步骤则需根据具体情况进行优化调整。

6 结论

本研究利用C波段Sentinel-1A卫星SAR影像数据，融合主成分时空分析与时序InSAR方法提取了徐州地区2018—2022年的高精度地表形变信息，详细分析了研究区地表形变时空演变规律，并讨论了研究方法的有效性和适用性以及近8年徐州市地表形变的综合变化趋势及其驱动因素，在此基础上得出以下结论：
（1）融合主成分时空分析与时序InSAR的地表形变信息提取方法能对时序InSAR信号进行深度挖掘与分析，通过构建时序地表形变模型，结合多项式函数，提取高精度、小量级的地表形变信息。相较于时空滤波和GACOS方法，本方法将徐州地区长时序地表形变的监测精度提高了10%~57%。本方法数据处理的整体思路可以应用于其它领域地表形变监测的研究，但细化的步骤则需根据相应研究区地表形变的特点、地形起伏情况、受大气等误差干扰的程度及对地表形变监测精度的需求等因素进行优化调整。
（2）2018—2022年，徐州市地表形变主要分布在城区、地铁沿线及西北、东北的煤矿老采空区范围。城区和地铁沿线呈现沉降趋势，最大形变速率达-17 mm/a；老采空区主要呈明显的地表抬升趋势，最大形变速率达35 mm/a。
（3）综合分析近8年（2015—2022）徐州市地表形变时空变化趋势可知，城市建设不断引发新的局部沉降区，随工程结束形变速率逐渐减小直至平缓；西北、东北及西部采空区中， 2012年及之前闭井的采空区当前已处于最终稳定阶段，而2014年后闭井的采空区仍处于持续形变阶段。特别地，断层分布影响了旗山矿的次生形变速率，城市建设减缓了张小楼矿的次生形变进程。因此在未来几年内，徐州地区城市基础设施建设仍可能产生新的沉降区，老采空区的地表仍会存在不稳定现象。为防止地质灾害的发生，仍需定期对相关区域进行地表形变时空监测。
本文研究结果可为城市地表形变高精度监测和潜在地质灾害防治提供重要的技术支撑和决策支持。

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
许才军, 张朝玉. 地壳形变测量与数据处理[M]. 武汉: 武汉大学出版社, 2009.

[ Xu C J, Zhang C Y. Crustal deformation measurement and data processing[M]. Wuhan: Wuhan University Press, 2009. ]

[2]
Herrera-García G, Ezquerro P, Tomás R, et al. Mapping the global threat of land subsidence[J]. Science, 2021, 371(6524):34-36. DOI:10.1126/science.abb8549

PMID

[3]
罗三明, 单新建, 朱文武, 等. 多轨PSInSAR监测华北平原地表垂直形变场[J]. 地球物理学报, 2014, 57(10):3129-3139.

[ Luo S M, Shan X J, Zhu W W, et al. Monitoring vertical ground deformation in the North China Plain using the multitrack PSInSAR technique[J]. Chinese Journal of Geophysics, 2014, 57(10):3129-3139. ] DOI:10.6038/cjg20141004

[4]
熊福文, 朱文耀. 长江三角洲地区地形变特征的GPS监测和分析[J]. 地球物理学报, 2007, 50(6):1719-1730.

[ Xiong F W, Zhu W Y. Land deformation monitoring by GPS in the Yangtze Delta and the measurements analysis[J]. Chinese Journal of Geophysics, 2007, 50(6):1719-1730. ] DOI:10.3321/j.issn:0001-5733.2007.06.012

[5]
Chaussard E, Wdowinski S, Cabral-Cano E, et al. Land subsidence in central Mexico detected by ALOS InSAR time-series[J]. Remote Sensing of Environment, 2014, 140:94-106. DOI:10.1016/j.rse.2013.08.038

[6]
Aslan G, Cakır Z, Ergintav S, et al. Analysis of secular ground motions in Istanbul from a long-term InSAR time-series (1992-2017)[J]. Remote Sensing, 2018, 10(3):408. DOI:10.3390/rs10030408

[7]
崔笃信, 王庆良, 李克, 等. 长白山天池火山近期形变场演化过程分析[J]. 地球物理学报, 2007, 50(6):1731-1739.

[ Cui D X, Wang Q L, Li K, et al. Analysis of recent deformation of Changbaishan Tianchi volcano[J]. Chinese Journal of Geophysics, 2007, 50(6):1731-1739. ] DOI:10.3321/j.issn:0001-5733.2007.06.013

[8]
Gao J X, Hu H. Advanced GNSS technology of mining deformation monitoring[J]. Procedia Earth and Planetary Science, 2009, 1(1):1081-1088. DOI:10.1016/j.proeps.2009.09.166

[9]
Madsen S N, Zebker H A, Martin J. Topographic mapping using radar interferometry: Processing techniques[J]. IEEE Transactions on Geoscience and Remote Sensing, 1993, 31(1):246-256. DOI:10.1109/36.210464

[10]
Ferretti A, Prati C, Rocca F. Permanent scatterers in SAR interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(1):8-20. DOI:10.1109/36.898661

[11]
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-2383. DOI:10.1109/TGRS.2002.803792

[12]
Hooper A. A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches[J]. Geophysical Research Letters, 2008, 35(16):L16302. DOI:10.1029/2008gl034654

[13]
朱猛, 董少春, 尹宏伟, 等. 基于SBAS InSAR方法的苏州地区2007—2010年地表形变时空变化研究[J]. 地球信息科学学报, 2016, 18(10):1418-1427.

DOI

[ Zhu M, Dong S C, Yin H W, et al. Spatial-temporal ground deformation study of Suzhou area from 2007 to 2010 based on the SBAS InSAR method[J]. Journal of Geo-Information Science, 2016, 18(10):1418-1427. ] DOI:10.3724/SP.J.1047.2016.01418

[14]
Zeni G, Bonano M, Casu F, et al. Long-term deformation analysis of historical buildings through the advanced SBAS-DInSAR technique: The case study of the city of Rome, Italy[J]. Journal of Geophysics and Engineering, 2011, 8(3):S1-S12. DOI:10.1088/1742-2132/8/3/S01

[15]
Zhou L, Guo J M, Hu J Y, et al. Wuhan surface subsidence analysis in 2015-2016 based on Sentinel-1A data by SBAS-InSAR[J]. Remote Sensing, 2017, 9(10):982. DOI: 10.3390/rs9100982

[16]
Goldstein R M, Werner C L. Radar interferogram filtering for geophysical applications[J]. Geophysical Research Letters, 1998, 25(21):4035-4038. DOI:10.1029/1998gl900033

[17]
Lee J S, Papathanassiou K P, Ainsworth T L, et al. A new technique for noise filtering of SAR interferogram phase images[C]. IGARSS'97.1997 IEEE International Geoscience and Remote Sensing Symposium Proceedings. Remote Sensing - A Scientific Vision for Sustainable Development, 2002:1716-1718. DOI:10.1109/IGARSS.1997.609040

[18]
Bekaert D P S, Walters R J, Wright T J, et al. Statistical comparison of InSAR tropospheric correction techniques[J]. Remote Sensing of Environment, 2015, 170:40-47. DOI:10.1016/j.rse.2015.08.035

[19]
Bekaert D P S, Hooper A, Wright T J. A spatially variable power law tropospheric correction technique for InSAR data[J]. Journal of Geophysical Research: Solid Earth, 2015, 120(2):1345-1356. DOI:10.1002/2014jb011558

[20]
Kositsky A P, Avouac J P. Inverting geodetic time series with a principal component analysis-based inversion method[J]. Journal of Geophysical Research, 2010, 115(B3):B03401. DOI:10.1029/2009jb006535

[21]
Lin Y N N, Kositsky A P, Avouac J P. PCAIM joint inversion of InSAR and ground-based geodetic time series: Application to monitoring magmatic inflation beneath the Long Valley Caldera[J]. Geophysical Research Letters, 2010, 37(23):137-139. DOI:10.1029/2010gl045769

[22]
Remy D, Froger J L, Perfettini H, et al. Persistent uplift of the Lazufre volcanic complex (Central Andes): New insights from PCAIM inversion of InSAR time series and GPS data[J]. Geochemistry, Geophysics, Geosystems, 2014, 15(9):3591-3611. DOI:10.1002/2014gc005370

[23]
Chen Y, Remy D, Froger J L, et al. Long-term ground displacement observations using InSAR and GNSS at Piton de la Fournaise volcano between 2009 and 2014[J]. Remote Sensing of Environment, 2017, 194:230-247. DOI:10.1016/j.rse.2017.03.038

[24]
Chen Y, Tan K, Yan S Y, et al. Monitoring land surface displacement over Xuzhou (China) in 2015-2018 through PCA-based correction applied to SAR interferometry[J]. Remote Sensing, 2019, 11(12): 1494. DOI:10.3390/rs11121494

[25]
周东来, 孙亚军, 杨国勇, 等. 徐州市水文地质单元岩溶地下水循环系统分析[J]. 中国矿业大学学报, 2003, 32(1):96-99.

[ Zhou D L, Sun Y J, Yang G Y, et al. Analyses of circulation system of Karst water in Xuzhou hydrogeological unit[J]. Journal of China University of Mining and Technology, 2003, 32(1):96-99. ] DOI:10.3321/j.issn:1000-1964.2003.01.023

[26]
杜培军, 单丹丹, 夏俊士, 等. 北京一号小卫星数据的城市景观格局监测分析——以徐州市城区为例[J]. 地球信息科学学报, 2010, 12(6):855-862.

[ Du P J, Shan D D, Xia J S, et al. Monitoring and analyzing urban landscape pattern change using Beijing-1 small satellite remote sensing data: Taking the urban area of Xuzhou city as an example[J]. Journal of Geo-Information Science, 2010, 12(6):855-862. ] DOI:10.1017/S0004972710001772

[27]
倪庆琳, 侯湖平, 丁忠义, 等. 基于生态安全格局识别的国土空间生态修复分区——以徐州市贾汪区为例[J]. 自然资源学报, 2020, 35(1):204-216.

DOI

[ Ni Q L, Hou H P, Ding Z Y, et al. Ecological remediation zoning of territory based on the ecological security pattern recognition: Taking Jiawang district of Xuzhou city as an example[J]. Journal of Natural Resources, 2020, 35(1):204-216. ]

DOI

[28]
徐州市人民政府. 徐州市国民经济和社会发展第十四个五年规划和2035年愿景目标纲要[EB/OL]. http://www.xz.gov.cn/govxxgk/014051247/2021-05-07/cccf1a4f-d320-47d4-8254-2ed3c562d5c7.html2021-04-19.

[ Xuzhou People's Government. The "Fourteenth Five Year Plan" for Xuzhou's national economic and social development and the outline of vision and objectives for 2035[EB/OL]. http://www.xz.gov.cn/govxxgk/014051247/2021-05-07/cccf1a4f-d320-47d4-8254-2ed3c562d5c7.html2021-04-19. ]

[29]
新华网. 习近平与江苏的故事[EB/OL]. https://baijiahao.baidu.com/s?id=1683256651521471033&wfr=spider&for=pc2020-11-13.

[ Xinhuanet. The story of Xi Jinping and Jiangsu[EB/OL]. https://baijiahao.baidu.com/s?id=1683256651521471033&wfr=spider&for=pc2020-11-13. ]

[30]
江苏省人民政府. 江苏省“十四五”自然资源保护和利用规划[EB/OL]. http://www.jiangsu.gov.cn/art/2021/11/1/art_64797_10093344.html2021-11-10.

[ Jiangsu Province People's Government. The "Fourteenth Five Year Plan" for natural resources protection and utilization in Jiangsu Province[EB/OL]. http://www.jiangsu.gov.cn/art/2021/11/1/art_64797_10093344.html2021-11-10.]

[31]
欧洲航天局[EB/OL]. https://scihub.copernicus.eu/

[European Space Agency[EB/OL]. https://scihub.copernicus.eu/

[32]
Farr T G, Rosen P A, Caro E, et al. The shuttle radar topography mission[J]. Reviews of Geophysics, 2007, 45(2): RG2004. DOI:10.1029/2005rg000183

[33]
Hooper A, Zebker H A. Phase unwrapping in three dimensions with application to InSAR time series[J]. Journal of the Optical Society of America A, Optics, Image Science, and Vision, 2007, 24(9):2737-2747. DOI:10.1364/josaa.24.002737

[34]
Li L P, Zhao C F, Huang J X, et al. Complex surface deformation of the coalfield in the northwest suburbs of Xuzhou from 2015 to 2020 revealed by time series InSAR[J]. Canadian Journal of Remote Sensing, 2021, 47(5):697-718. DOI:10.1080/07038992.2021.1951190

[35]
彭磊. 徐州西部关闭矿井地表次生形变InSAR监测与分析[D]. 徐州: 中国矿业大学, 2021.

[ Peng L. Monitoring and analysis of secondary deformation of closed mine in western Xuzhou based on InSAR[D]. Xuzhou: China University of Mining and Technology, 2021. ]

[36]
郑美楠, 邓喀中, 张宏贞, 等. 基于InSAR的关闭矿井地表形变监测与分析[J]. 中国矿业大学学报, 2020, 49(2):403-410.

[ Zheng M N, Deng K Z, Zhang H Z, et al. Monitoring and analysis of surface deformation in closed mine based on InSAR[J]. Journal of China University of Mining & Technology, 2020, 49(2):403-410. ] DOI:10.13247/j.cnki.jcumt.001133

[37]
卞正富, 张燕平. 徐州煤矿区土地利用格局演变分析[J]. 地理学报, 2006, 61(4):349-358.

[ Bian Z F, Zhang Y P. Land use changes in Xuzhou coal mining area[J]. Acta Geographica Sinica, 2006, 61(4):349-358. ] DOI:10.3321/j.issn:0375-5444.2006.04.002

[38]
邓喀中, 郑美楠, 张宏贞, 等. 关闭矿井次生沉陷研究现状及展望[J]. 煤炭科学技术, 2022, 50(5):10-20.

[ Deng K Z, Zheng M N, Zhang H Z, et al. Research status and prospect of secondary subsidence in closed mine[J]. Coal Science and Technology, 2022, 50(5):10-20. ] DOI:10.13199/j.cnki.cst.2021-1403

[39]
张丽, 黄敬军, 许书刚, 等. 徐州城市规划区煤矿采空区稳定性评价[J]. 水文地质工程地质, 2017, 44(2):124-128.

[ Zhang L, Huang J J, Xu S G, et al. Stability evaluation of goaf in Xuzhou urban planning area[J]. Hydrogeology and Engineering Geology, 2017, 44(2):124-128. ] DOI:10.16030/j.cnki.issn.1000-3665.2017.02.19

Options
Outlines

/

〈 〉