邢台市2019~2021年地表形变SBAS-InSAR监测及分析

SBAS-InSAR Monitoring and Analysis of Surface Deformation in Xingtai City from 2019 to 2021

摘要: 长期超采地下水导致邢台市地表发生了大范围的沉降，这严重威胁了邢台市社会经济发展。本文以邢台市为研究区，采用2019年1月3日至2021年1月28日年的64景Sentinel-1A影像数据，利用短基线集(Small Baseline Subset Interferometric Synthetic Aperture Radar, SBAS-InSAR)技术对邢台市地面沉降进行监测，以获取该区地面沉降速率和累积沉降量，探究邢台市地面沉降状况，为邢台市地面沉降灾害的防治及资源开发利用规划提供基础依据。研究发现，邢台市东部存在一个范围较大的沉降区域，沉降面积约为1689 km²，沉降速率约为80 mm/year，沉降范围没有扩大，但呈现出沉降加剧的趋势；此外，邢台市西部有多个面积较小的沉降漏斗，主要分布在邢台县(Q1)和内丘县(Q2)，最大沉降速率达到150 mm/year。在此基础上，对邢台市东部地下水埋深变化及西部矿区分布与沉降区分布进行了对比分析，结果表明东部地下水埋深高的区域地面沉降呈加剧趋势，西部矿区开采对邢台市的地面沉降产生一定的影响。

Abstract: The long-term over-exploitation of groundwater has led to a large-scale subsidence of the surface of Xingtai City, which seriously threatens the economic and social development of Xingtai City. This paper takes Xingtai City as the research area, and uses 64-scene Sentinel-1A image data from January 3, 2019 to January 28, 2021. Small Baseline Subset Interferometric Synthetic Aperture Radar (SBAS-INSAR) technology was used to monitor the land subsidence in Xingtai City, so as to obtain the land subsidence rate and cumulative subsidence volume in this area, explore the land subsidence situation in Xingtai City, and provide a basis for the prevention and control of land subsidence disasters and resource development and utilization planning in Xingtai City. The study found that there was a large subsidence area in the east of Xingtai City, with a subsidence area of about 1689 km² and a subsidence rate of about 80 mm/year. The subsidence range did not expand, but showed a trend of intensified subsidence; in addition, there are several smaller sedimentation funnels in the western part of Xingtai City, which are mainly distributed in Xingtai County (Q1) and Neiqiu County (Q2), with a maximum settlement rate of 150 mm/year. On this basis, this paper compares and analyzes the changes of groundwater burial depth in eastern Xingtai City and the distribution of mining areas in western China and subsidence areas. The results show that the land subsidence in the region with high groundwater depth in the east is increasing, and the mining in the west has a certain impact on the land subsidence in Xingtai City.

1. 引言

地面沉降是地表在一定时期内不断降低的环境地质现象，是一种不可逆的地质灾害。地面沉降的危害主要表现为毁坏建筑物和基础设施、造成海水倒灌、不利用建设事业和资源开发等。地下水位下降是引起华北地区地面沉降的主要原因。区域性地下水过度开采容易诱发地面沉降、地裂缝、岩溶塌陷等地质灾害，破坏建筑基础设施，而且还增加了发生洪水的风险 [1]。地下水过度开采后导致水位急剧下降，使得储存地下水的砂土或黏土等介质受到挤压，地面发生沉降。地下水过度开采所导致的地面沉降是一个缓慢的过程，一般在较长时间尺度上发生并且在较大区域范围引起了地面高程的逐渐降低 [2]。长期过度开采地下水，会导致土层固相结构承压超限，增大了地质灾害风险 [3]，由地下水开采造成的地面沉降的现象在我国北方地区尤为显著。

邢台市位于华北平原西部边缘，是京津冀城市群节点城市，是冀中南重要的物流枢纽和制造业基地。邢台市受到温带大陆性季风气候控制，年降水量约为490 mm [4]，地表水资源缺乏，社会经济发展受到地下水的严重制约 [5]。根据河北省水务局公布的地下水公报显示邢台市存在长期的地下水超采，这容易导致地面沉降和地面裂缝等地质灾害，制约着邢台市的基础设施建设和社会经济发展。因此对邢台市进行地表时序地表形变监测是十分有必要的。

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

因此本文基于SBAS-InSAR技术，利用2019年1月3日至2021年1月28日年的64景Sentinel-1A影像对邢台市地面沉降进行监测，获取该区地面沉降速率和累积沉降量，分析地面沉降地下水之间的关系，进而揭示邢台市地面沉降状况。

2. SBAS-InSAR原理

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

$\frac{N+1}{2}\le M\le N\left(\frac{N+1}{2}\right)$ (1)

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

$\delta {\varphi }_j={\varphi }_B\left(x,r\right)-{\varphi }_A\left(x,r\right)\approx \frac{4\pi }{\lambda }\left[d\left(t_B,x,r\right)-d\left(t_A,x,r\right)\right]$ (2)

式中 ${\varphi }_B\left(x,r\right)$ 和 ${\varphi }_A\left(x,r\right)$ 分别表示雷达坐标系中，坐标为(x，r)的像素点在t_B时刻与t_A时刻形变相位， $\lambda$ 为波长， $d\left(t_B,x,r\right)$ 和 $d\left(t_A,x,r\right)$ 分别表示雷达坐标系中，坐标为(x，r)的像素点在t_B时刻与t_A时刻形变量。由此所有时序差分干涉对的形变相位构建观测方程为：

$\Delta \varphi \left(x,r\right)=A\cdot \varphi \left(x,r\right)+\Delta {\varphi }_{\epsilon }\left(x,r\right)$ (3)

$\Delta \varphi \left(x,r\right)$ 表示M个差分干涉对的形变相位数据集， $\varphi \left(x,r\right)$ 表示其它N幅SAR影像相对于参考影像的累积形变相位， $\Delta {\varphi }_{\epsilon }\left(x,r\right)$ 表示每个差分干涉对残余误差相位，A表示 $M\times N$ 维系数矩阵，具体形式如下：

$A=\left[\begin{array}{cccccccc}0& -1& 0& 0& 0& \cdots & 0& 0\\ & 0& +1& 0& & \cdots & 0& 0\\ ⋮& ⋮& 0& & & \cdots & & ⋮\\ & & ⋮& ⋮& ⋮& \cdots & +1& 0\\ 0& 0& 0& \cdots & 0& -1& 0& +1\end{array}\right]$ (4)

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

$\hat{\varphi }\left(x,r\right)={\left(A^{\text{T}}A\right)}^{-1}{A}^{\text{T}}\Delta \varphi \left(x,r\right)$ (5)

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

$A=US{V}^{\text{T}}$ (6)

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

$\hat{\varphi }=V{S}^{-1}{U}^{\text{T}}\Delta \varphi \left(x,r\right)$ (7)

最后，利用最小二乘方法求解各未知参数，进而获取不含DEM高程残差的时序形变相位信息。

3. 研究区域概况和数据

3.1. 研究区域概况

邢台市位于河北省南部，毗邻石家庄市、衡水市和邯郸市，面积约为12456 km² [10]，邢台市拥有耕地67.56 hm²，其中50.83 hm²的耕地为水浇耕地 [11]。邢台市河道众多，但是除留垒河外均发源于西部山区，且多为季节性河道，基流少，汛期河水暴涨，历时却很短 [10]。除此之外，资料显示邢台市多年平均年降水量约为490 mm，蒸发量高达1160 mm [12]，年内降水集中于6~9月，约占全年降水的78%~80%，因此地下水成为邢台市用水的主要来源。

长期的地下水过度开采导致了邢台市发生大范围的地表沉降，为了及时掌握该地区地表沉降情况，探究邢台市地下水开采和地表沉降的相关性。本文基于SBAS-InSAR方法，利用64景sentinel-1A升影像获取邢台市201~2020年的地表时序形变。在进行SBAS处理前对SAR影像进行裁剪获取本文的研究区域如图1所示，其中黄色实线表示SAR影像覆盖范围，红色实线表示裁剪后的影像范围。

Figure 1. Overview of the study area

图1. 研究区域概况图

3.2. 数据来源

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

Table 1. Basic parameters of Sentinel-1A IW mode image

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

4. 结果与分析

本文采用2019年1月3日至2021年1月28日的64景C波段Sentinel-1A升轨数据对邢台市进行沉降监测。选择2020年2月3日的影像作为公共主影像，使其他的63景影像与其配准。使用SBAS-InSAR时序形变解算方法，以最大程度地减少大气误差、提高形变监测精度。设定垂直基线阈值为200米，时间基线阈值为48天，分别获得245个差分干涉组合，并分别得到了各干涉组合的差分结果，其时空基线分布如图2所示。图中蓝色实线表示平均空间相干性大于0.7的干涉对，红色虚线表示不参与形变解算的平均空间相干性低于0.7的干涉对。

Figure 2. Interferometric network space-time baseline map

图2. 干涉网络时空基线图

4.1. 邢台市地面沉降现状

本文采用SBAS-InSAR技术获取邢台市在2019年1月至2021年1月的时序形变，获取的形变速率图如图3所示，从图中可以看出邢台市东部存在一个范围较大的沉降区域，沉降面积约为1689 km²沉降速率约为80 mm/year。此外，邢台市西部存在多个面积较小的沉降漏斗，这些沉降漏斗主要分布在邢台县(Q1)和内丘县(Q2)，最大沉降速率达到150 mm/year。将2019年12月5日的影像作为参考影像，即地表形变量为0，其他影像以该影像为基准，获取与参考影像之间的累积形变图。如图4所示，在2019年1月至2021年1月期间，邢台市东部(红色框)存在持续的沉降，沉降范围没有扩大，但是却呈现出沉降加剧的趋势。

Figure 3. Figure of annual average deformation rate of Xingtai City

图3. 邢台市年平均形变速率图

Figure 4. Cumulative deformation sequence of Xingtai City from January 2019 to January 2021

图4. 邢台市2019年1月至2021年1月累积形变序列

为了揭示邢台市东部持续沉降规律，在邢台市东部依次选取7个特征点，7个特征点大致沿沉降区东北至西南排列，点位编号依次为Point1~Point7，点的位置如图5左列所示。从图5右列可以出7个特征点(Point1~Point7)的累积沉降量分别为165、202、145、160、138、131、122 mm，累积沉降量最小的是Point7；累积沉降量最大的点是Point2。7个特征点虽然累积沉降量不同，却有着相同的沉降趋势，每年3月至7月沉降速率加快，7月至10月形变波动较大，年内其它月份形变速率变缓。

Figure 5. Location and deformation time series of characteristic points in the eastern subsidence area of Xingtai City

图5. 邢台市东部沉降区特征点位置及特征点形变时间序列

研究区沉降面积较小的区域主要集中在Q1及Q2 (如图3所示)。为了进一步探究这些区域的形变规律，本文在这两个区域中提取10个特征点，特征点的位置和形变时间序列如图6所示。在Q1中选取5个特征点(Point8~Point12)的累积沉降量分别为262、208、126、81、181 mm；在Q2中选取5个特征点(Point13~Point17)的累积沉降量分别为222、172、104、141、247 mm。选的特征点均表现出整体下沉的趋势，其中Point8和Point9在2020年8月底形变速率突然增加。

Figure 6. Location and deformation time series of settlement characteristic points in Q1 and Q2 regions

图6. Q1及Q2区域沉降特征点位置及形变时间序列

4.2. 邢台市地面沉降原因分析

邢台市是河北省南部严重缺水的地区，人均水资源占有量仅220 m³，地下水开采一直是邢台市生活用水的主要来源。自上世纪70年代开始邢台市开始大规模开采地下水，长期地下水开采导致地面沉降、咸水下移、环境破坏等一系列问题 [14]。为了防止情况继续恶化，邢台市2014年开始全面开展地下水超采综合治理工作并取得了一定的成果，2014年至2019年地下水供水占比逐年下降 [15]。但是根据河北省地下水超采区地下水位监测情况通报显示2018年至2020年邢台市部分地区地下水埋深依然存在下降的趋势，地下水的持续开采将加剧地面沉降。

邢台市的地表沉降与地下水开采息息相关，为了探究邢台市地面沉降与地下水开采的关系，本文根据河北省水利厅公布的河北省地下水超采区地下水位监测情况通报，获取了位于邢台市东部沉降严重区域的巨鹿县和广宗县的月平均深层地下水埋深变化，如图7所示。从图中可以看出，广宗县和巨鹿县的

Figure 7. Average depth of deep groundwater in Julu County and Guangzong County from 2018 to December 2020

图7. 巨鹿县和广宗县2018年月至2020年12月平均深层地下水埋深

深层地下水埋深具有相同的变化趋势；在3月至6月急剧减小(地下水开采量较大)。深层地下水埋深减小速率变大时，沉降速率也增大，这说明邢台市东部地表沉降主要是由地下水超采造成的。邢台市有50.83 hm²的水浇耕地，3月至5月正值农作物生长期，需要大量灌溉，因此本文猜测农作物灌溉导致地下水的快速消耗。通过和风天气(https://www.qweather.com)获取到邢台市的降雨集中在6至9月份，雨量的增加使地下水得到补充，地下水埋深增加。

此外，图3中Q1与Q2区域存在的小范围地表沉降可能是由矿产资源开采造成的。邢台县和内丘县是邢台市存在大量的矿区 [16]，因此本文推测这两个区域的小范围地面沉降是由矿产资源开采造成的。

5. 总结

SBAS-InSAR可以获取地表毫米级时序形变，已经被广泛的应用于大范围地表时序形变监测。本文采用SBAS-InSAR技术对覆盖邢台市的64景Sentinel-1A数据进行处理，获取了邢台市2019年1月至2021年1月的地表时序形变。研究结果表明：

1) 邢台市地面沉降整体上呈现东西差异，邢台市东部存在大面积的沉降区域，西部存在多个面积较小的沉降漏斗。

2) 邢台市东部的大范围沉降可能是由地下水开采造成的，地下水埋深与地面沉降存在一定相关性。地面沉降速率大的区域地下水埋深高，地面沉降速率小的区域地下水埋深低。

3) 位于邢台市西部的邢台县和内丘县存在多个小范围沉降漏斗，本文猜测这些沉降漏斗是由矿区开采造成的。

本文通过一系列工作取得了一些成果，但是仍然存在着一些问题。首先，仅通过公布公报获取县级行政区域的地下水埋深，并没有找到公开的邢台市地下水监测井站点数据，无法将同一点位的地下水埋深数据与地面沉降数据进行对比。此外，本文只通过资料查阅到邢台县和内丘县存在许多矿区，但并没有获取矿区的具体位置。尽管本文存在这些问题，但是本文的实验结果表明邢台市东部仍然存在范围的地下水超采，造成了严重的地面沉降，带来了严重的安全隐患。本文结果对邢台市地下水综合治理具有一定的价值。