2022, Vol. 42
2. 地理信息工程国家重点实验室，西安市雁塔路中段1号，710054;
3. 西安市勘察测绘院，西安市南二环东段29号，710054
地面沉降是指在自然因素和人为因素的作用下，由于地壳表层未固结土体压缩变形而引起的区域性地面高程降低的环境地质现象[1]，具有持续时间长、影响范围大和多种因素共同作用等特点[2]。南水北调中线工程是世界上最大的跨流域调水工程，意在解决华北平原区域性缺水的问题，保障干渠的安全运营是推进受水区人民生活与经济发展水平的重要措施[3]。而工程沿线部分渠段受城市地下水开采、矿区开采及区域地质构造活动等多种威胁，区域地质结构稳定性较差[4]。为避免沿线城市地面沉降影响到中线工程正常通水运行，对沿线区域进行长时间、高精度的地表形变监测具有重要意义。
与传统监测技术相比，合成孔径雷达干涉测量InSAR技术具有获取数据速度快、可覆盖范围大、全天时全天候监测及精度高等优势，被广泛用于各个领域[5]。由于受到时空失相干和大气延迟等的影响，常规差分干涉测量D-InSAR技术的监测精度较低，为突破D-InSAR技术的限制，小基线集技术SBAS-InSAR[6]应运而生。该技术大幅提高了形变监测精度，并且在探测长期累积的缓慢形变方面表现出极大的潜力。
本文利用SBAS-InSAR技术对南水北调中线工程沿线2015-04~2020-11的地表形变进行研究，得到沿线区域形变特征分布及时空演变规律，并结合有关资料着重探讨了北京市地面沉降在南水进京后的时空演变特征，其结果可为相关单位在防灾减灾决策方面提供参考。
1 SBAS-InSAR数据处理方法 1.1 SBAS-InSAR基本原理2002年，Berardino等[6]首次提出小基线集SBAS技术，该方法广泛应用于长时间序列的地表形变监测中，基本原理如下[7]：若有同一研究区域的N+1景SAR影像，则可能得到M个干涉对，M需满足：
| $ \frac{N+1}{2} \leqslant M \leqslant \frac{N(N+1)}{2} $ | (1) |
设初始时刻为t0，任意时刻ti(i=1, ⋯, N)相对于初始时刻的相位差是未知参数φ(ti)，干涉处理获取的差分干涉相位δφ(tk)(k=1, ⋯, M)为观测量。若不考虑其他相位误差，则第k景差分干涉图的像元(x，y)地表形变相位δφdef可表示为：
| $ \begin{array}{c} \delta \varphi_{\mathrm{def}}=\varphi\left(t_B, x, y\right)-\varphi\left(t_A, x, y\right)= \\ \frac{4 \pi}{\lambda}\left(r_{t_B}-r_{t_A}\right)=\frac{4 \pi}{\lambda} \delta_r \end{array} $ | (2) |
式中，tA、tB(tA＜tB)为干涉图k(k=1, ⋯, M)对应的SAR影像获取时刻，δφdef为tA到tB时刻之间视线向的形变相位，λ为雷达波长，rtB、rtA分别为tA和tB时刻干涉图的任意像元对应视线向的地表形变，δr为该时段的地表形变量。
将所有差分干涉图进行自由组合，φ为N幅SAR影像上的干涉相位值组成的矩阵，δφ为M幅差分干涉图上相位组成的矩阵。则有：
| $ \boldsymbol{\varphi}=\left[\begin{array}{lll} \varphi\left(t_1\right) & \cdots & \varphi\left(t_N\right) \end{array}\right]^{\mathrm{T}} $ | (3) |
| $ \delta \boldsymbol{\varphi}=\left[\begin{array}{lll} \delta \varphi\left(t_1\right) & \cdots & \delta \varphi\left(t_{M}\right) \end{array}\right]^{\mathrm{T}} $ | (4) |
将上式写成矩阵形式：
| $ \boldsymbol{A} \boldsymbol{\varphi}=\delta \boldsymbol{\varphi} $ | (5) |
式中，A为M×N维矩阵。
当小基线集合数L=1时，矩阵A的秩为N，采用最小二乘法即可求解出累积形变相位的估值
| $ \hat{\varphi}=\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}\right)^{-1} \boldsymbol{A}^{\mathrm{T}} \delta \boldsymbol{\varphi} $ | (6) |
当基线组L>1时，式(5)秩亏，秩亏数为N－L+1，为此对系数矩阵A进行奇异值分解，求出累积形变相位φ最小范数意义上的最小二乘解。
1.2 数据处理设置一定的时空基线阈值，满足条件的SAR影像自由组合生成干涉图集。将DEM数据与主影像进行配准，模拟相位对去平后的干涉图进行地形相位去除，得到差分干涉图。对差分干涉图进行滤波来降低由于失相干引起的噪声等，并对滤波后得到的干涉图进行相位解缠。构建矩阵方程，基于奇异值分解(SVD)方法估算解缠相位的高程误差和平均形变速率，通过时间域和空间域滤波对残余相位中的失相干噪声相位和大气延迟相位进行去除，得到非线性形变速率。结合线性形变速率和非线性形变速得到不同影像时间序列间的地表形变速率结果。SBAS方法的数据处理技术路线如图 1所示。
|
图 1 SBAS数据处理技术路线 Fig. 1 Data processing technique flowchart by SBAS |
南水北调中线工程输水干渠源自丹江口水库，流经河南、河北，最后到达北京、天津，总长1 432.8 km。自2014年底通水后，中线工程已全面进入运行阶段，截至2021-07-19，南水北调中线一期工程自陶岔渠首累积调水入渠水量达400亿m3，分别向河南省、河北省、天津市、北京市供水135亿m3、116亿m3、65亿m3及68亿m3，沿线直接受益人口增加至7 900万人，比2015年通水1周年时的3 800万受益人口增加1倍多。
2.2 实验数据来源Sentinel-1A卫星是欧空局继ERS和Envisat系列卫星后发射的环境监测卫星，基于C波段的成像系统采用条带成像、干涉宽幅、超宽幅及波浪模式4种成像模式，其中干涉宽幅模式(IW)采用中等分辨率5 m×20 m获取幅宽250 km的影像，成为陆地覆盖的默认模式。本文收集IW模式下Sentinel-1A数据，其覆盖情况如图 2所示，详细数据参数如表 1所示。除此之外，利用欧空局提供的POD精密轨道数据提高SAR卫星影像轨道精度，利用美国宇航局提供的分辨率为30 m的SRTM DEM数据去除地形相位。
|
图 2 研究区概况和数据覆盖情况 Fig. 2 Overview and data coverage of the study area |
|
表 1 Sentinel-1A数据参数 Tab. 1 Data parameter of Sentinel-1A |
由于研究区域较大，需进行数据分块处理。将每景影像经SBAS-InSAR技术处理后的结果转换到地理坐标系下，在ArcGIS软件中依据一定的规则进行栅格镶嵌拼接，得到中线工程沿线区域形变结果。
图 3为中线工程沿线区域2015-04~2020-11 LOS向的平均形变速率，由图可见，整个沿线区域存在多处地面沉降，主要分布在京津冀地区，包括北京市、保定市、天津市、廊坊市、邢台市及邯郸市等。大部分沉降分布在河北省东南部，几乎接连成片，最大地面沉降速率达139 mm/a，位于邢台市南宫市沉降中心。北京市最大地面沉降速率达133 mm/a，位于朝阳区沉降中心。天津市西南部最大地面沉降速率为81 mm/a，相对较小。位于河北省东南部的地面沉降由于距离中线干渠有一定的距离，对输水影响较小，但天津市支线经过2个小范围沉降区，应引起重视。
|
图 3 中线沿线区域2015~2020年LOS向平均形变速率 Fig. 3 Mean deformation velocity in LOS of the middle route project from 2015 to 2020 |
由于覆盖各个区域的影像日期不统一，在计算累积形变时以2016-12为时间基准。图 4为以2016-12为时间起点，至2017-12、2018-12、2019-12、2020-11获取的中线沿线区域LOS向累积形变分布。在此期间，整个中线沿线的地表累积形变与形变速率具有相似的空间分布格局，且随着时间的推进呈持续发展趋势。最大累积沉降量位于北京市朝阳区沉降中心，2016-12~2017-12的最大沉降量为158 mm，2016-12~2018-12的最大累积沉降量为287 mm，相比于2016-12~2017-12，2016-12~2018-12的累积沉降量增量为129 mm；2016-12~2019-12的最大累积沉降量为402 mm，比2016-12~2018-12增加了115 mm；2016-12~2020-11的最大累积沉降量为476 mm，比2016-12~2019-12增加了74 mm。
|
图 4 中线沿线区域2016-12~2020-11 LOS向累积形变量 Fig. 4 Cumulative deformation in LOS of the middle route from December 2016 to November 2020 |
以天津支线主干渠两侧5 km缓冲区范围为研究对象，利用SBAS-InSAR技术获取覆盖天津支线的2015-09~2020-11年平均形变速率(图 5)。沿线区域的最大形变速率为-94 mm/a，整条线贯穿地面沉降区域，主要沉降区有2个：沉降区Ⅰ位于保定市与廊坊市交界处，区内最大平均形变速率为-94 mm/a；沉降区Ⅱ位于廊坊市与天津市交界处，区内最大平均形变速率为-81 mm/a。地面沉降的不均匀分布容易对地下箱涵造成破坏，产生裂缝，从而导致渗水。相比于干渠，沉降区Ⅰ、Ⅱ内沉降的不均匀分布较为明显，应对上述典型地段进行实地调查，检查箱涵的健康状况，防患于未然。
|
图 5 天津支线2015-07~2020-11 LOS向平均形变速率 Fig. 5 Mean deformation velocity in LOS of the Tianjin branch line from July 2015 to November 2020 |
为探究区域地面沉降在南水北调中线工程通水后的时空演变特征，对北京市的地面沉降进行重点分析。北京市2015-11~2020-11 LOS向平均形变速率如图 6所示，由图可见，北京市地面沉降空间分布差异性很大，沉降区域主要分布在昌平区、海淀区、顺义区、朝阳区、通州区及大兴区。各个行政区内形成多个沉降漏斗，分别为HD(海淀)、CP(昌平)、SY(顺义)、CY(朝阳)、TZ(通州)、DX(大兴)，其中CY的年均沉降速率最大，达到133 mm/a，HD、CP、SY、TZ、DX的年均沉降速率分别为89.72 mm/a、59.07 mm/a、57.01 mm/a、95.67 mm/a和69.15 mm/a。
|
图 6 北京市2015-11~2020-11 LOS向平均形变速率 Fig. 6 Mean deformation velocity in LOS of the Beijing from November 2015 to November 2020 |
为评定InSAR监测结果精度，获取2017~2019年北京市10个GPS站点的监测结果，依据GPS站点位置选择InSAR结果的相应点进行对比。为保证时间的一致性，获取该时段内LOS向InSAR形变速率分布，将GPS监测结果转换到LOS向与其对应的InSAR结果进行对比分析，如表 2(单位mm/a)所示。由表可见，两者之间具有较高的一致性，互差为0~7 mm，GPS观测值与InSAR结果之间的RMSE为2.47 mm/a，并表现出明显的线性相关，最大线性相关系数R2达0.974，说明本次InSAR监测结果精度较高。
|
|
表 2 InSAR结果与GPS观测值比较 Tab. 2 Comparison between time series InSAR measurements and GPS measurements |
为揭示北京市地面沉降形变特征和演化规律，获取北京市2015-11-03~2020-11-24 LOS向时序累积形变量。图 7为2015-11-03~2020-11-24北京市地面沉降空间特征演化过程，图 8为该时段的累积形变分布。由图 7和8可见，北京市的地面沉降处于持续加重状态，最大累积沉降量为697 mm，位于朝阳区。为对比沉降区域的形变特征是否具有相似趋势，选取4个点(点位见图 8)并提取其线性形变时序。由图 9可见，最大累积沉降点位于F4，累积沉降量为507.14 mm，最小累积沉降点位于F2，累积沉降量为342.11 mm。4个特征点呈不同速度的非线性沉降趋势，经历了从快速到逐渐减缓的沉降形变过程，由此可知，北京市的地面沉降正在减缓。
|
图 7 北京市2015-11-03~2020-11-24地面沉降空间特征演化过程 Fig. 7 The evolution of spatial characteristics of ground subsidence in Beijing from November 3, 2015 to November 24, 2020 |
|
图 8 北京市2015-11-03~2020-11-24累积形变 Fig. 8 Cumulative deformation in Beijing from November 3, 2015 to November 24, 2020 |
|
图 9 4个点的形变时序 Fig. 9 Time series deformation of 4 points |
图 10为2016~2020年北京市LOS向年均形变速率分布。由图可见，不同时段内的地面沉降分布相似，研究时段内沉降速率最大的区域始终位于朝阳区，2016年最大沉降速率为167 mm/a，2017年为175 mm/a，2018年为145 mm/a，2019年最大沉降速率为136 mm/a，2020年最大沉降速率持续减缓，变为98 mm/a。
|
图 10 北京市2016~2020年LOS向形变速率 Fig. 10 Mean deformation velocity in LOS of the Beijing from 2016 to 2020 |
为进一步分析地面沉降速率随时间的演化过程，选取地面沉降分布区域内的典型特征点。特征点的选取基于局部沉降区域中较为严重的沉降中心，即图 6中的HD、CP、SY、CY、TZ、DX，分别提取各点年均速率进行对比分析。由图 11可见，CP的沉降速率变化最为明显，由134.69 mm/a变为19.16 mm/a；HD的沉降速率从132.95 mm/a变为89.72 mm/a；SY的沉降速率从91.9 mm/a变为48.29 mm/a；CY的沉降速率从165.94 mm/a变为67.76 mm/a；TZ的沉降速率从135.08 mm/a变为54 mm/a；DX的沉降速率从74.86 mm/a变为63.82 mm/a。虽然部分特征点在中间时段的沉降速率有所增加，但相比于2016年，2020年6个点的沉降速率均减小，一定程度上反映了北京市的地面沉降程度正在减缓。
|
图 11 主要沉降区域特征点形变速率对比 Fig. 11 Comparison of deformation velocity of feature points in main subsidence areas |
地下水的过量开采会使含水层孔隙度逐渐变小，导致有效应力增加，从而引发地面沉降。北京市由于水量短缺及水体污染导致可用的地表水所剩无几，不得不长期超采地下水，从而引发大范围的地面沉降。众多学者结合北京市的地面沉降现状，分析了其与地下水演化的相关性[8-9]。
北京市密云水库长期担负着城市生活及农业生产用水的重要任务，仅1999~2003年密云水库水量就萎缩了3/4，全市超过70%的用水量只能靠抽取地下水维持，因此北京平原地区的地下水位以每年1 m的速度持续下降。尽管2003年后的10 a里，北京通过各项节水措施使用水量下降近七成，22%的用水也已被再生水替代，但地表水稀缺的现实、用水量增长的趋势难以改变，地下水位仍在逐年下降。2014年底南水北调中线工程进入全面输水阶段，2015年北京地下水水位自1999年以来首次上升，截至2018年，地下水水位已上升了1.2 m[10]。截至2020-08底，进入密云水库的南水累积超过5亿m3，蓄水量达到23亿m3，水库水位超过147 m，彻底扭转了供水量入不敷出的局面。因此，由北京地下水水位与地面沉降的对应关系可知，南水北调工程与该区域的地面沉降缓解密切相关。
5 结语1) 利用SBAS-InSAR技术可快速获取研究区地面沉降信息，通过与GPS观测值对比可知，InSAR监测结果精度较高，可为研究区域地面沉降现状及时空演化特征提供参考依据。
2) 北京市地面沉降虽处于持续发展阶段，但沉降速率呈逐渐减缓趋势，南水北调工程在一定程度上解决了华北平原的缺水问题，地下水的开采量也随之减少，从而缓解了地面沉降。
3) 工程沿线的地面沉降分布在空间上极不均匀，河北省东南部的地面沉降虽几乎连接成片，但不会影响正常输水。而天津市支线经过2个小沉降区，输水存在安全隐患，为保障南水北调中线工程的正常运行，需对其进行深入研究。
| [1] |
许强, 蒲川豪, 赵宽耀, 等. 延安新区地面沉降时空演化特征时序InSAR监测与分析[J]. 武汉大学学报: 信息科学版, 2021, 46(7): 957-969 (Xu Qiang, Pu Chuanhao, Zhao Kuanyao, et al. Time Series InSAR Monitoring and Analysis of Spatiotemporal Evolution Characteristics of Land Subsidence in Yan'an New District[J]. Geomatics and Information Science of Wuhan University, 2021, 46(7): 957-969)
(  0) |
| [2] |
邱志伟, 岳建平, 汪学琴, 等. 基于短基线集技术的城市地表沉降监测研究[J]. 测绘通报, 2016(7): 25-29 (Qiu Zhiwei, Yue Jianping, Wang Xueqin, et al. Research on Urban Surface Subsidence Monitoring Based on SBAS[J]. Bulletin of Surveying and Mapping, 2016(7): 25-29)
( 0) |
| [3] |
张永光, 田凡, 李迎春, 等. 时序InSAR在南水北调中线形变监测中的应用[J]. 长江科学院院报, 2021, 38(8): 72-77 (Zhang Yongguang, Tian Fan, Li Yingchun, et al. Application of Time Series InSAR to Deformation Monitoring in Central Line Project of South-to-North Water Transfer[J]. Journal of Yangtze River Scientific Research Institute, 2021, 38(8): 72-77)
( 0) |
| [4] |
马超, 屈春燕, 孟秀军. 南水北调总干渠中线工程豫北段基础稳定性的InSAR时序分析[J]. 地震地质, 2014, 36(3): 749-762 (Ma Chao, Qu Chunyan, Meng Xiujun. Embankment Stability of the North Henan Section of Middle Route Project of South-to-North Water Diversion Based on InSAR Time Series Analysis[J]. Seismology and Geology, 2014, 36(3): 749-762 DOI:10.3969/j.issn.0253-4967.2014.03.016)
( 0) |
| [5] |
Fryksten J, Nilfouroushan F. Analysis of Clay-Induced Land Subsidence in Uppsala City Using Sentinel-1 SAR Data and Precise Leveling[J]. Remote Sensing, 2019, 11(23): 2 764 DOI:10.3390/rs11232764
( 0) |
| [6] |
Berardino P, Fornaro G, Lanari R, et al. ANew Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(11): 2 375-2 383 DOI:10.1109/TGRS.2002.803792
( 0) |
| [7] |
朱建军, 李志伟, 胡俊. InSAR变形监测方法与研究进展[J]. 测绘学报, 2017, 46(10): 1 717-1 733 (Zhu Jianjun, Li Zhiwei, Hu Jun. Research Progress and Methods of InSAR for Deformation Monitoring[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(10): 1 717-1 733)
( 0) |
| [8] |
周朝栋, 宮辉力, 张有全, 等. 北京平原区地面沉降PS-InSAR监测[J]. 遥感信息, 2017, 32(1): 17-22 (Zhou Chaodong, Gong Huili, Zhang Youquan, et al. Land Subsidence Monitoring in Beijing Plain with PS-InSAR Techniques[J]. Remote Sensing Information, 2017, 32(1): 17-22)
( 0) |
| [9] |
田芳, 罗勇, 周毅, 等. 北京地面沉降与地下水开采时空演变对比[J]. 南水北调与水利科技, 2017, 15(2): 163-169 (Tian Fang, Luo Yong, Zhou Yi, et al. Contrastive Analysis of Spatial-Temporal Evolution between Land Subsidence and Groundwater Exploitation in Beijing[J]. South-to-North Water Transfers and Water Science and Technology, 2017, 15(2): 163-169)
( 0) |
| [10] |
Lü M Y, Ke Y H, Lin G, et al. Change in Regional Land Subsidence in Beijing after South-to-North Water Diversion Project Observed using Satellite Radar Interferometry[J]. Mapping Science and Remote Sensing, 2020, 57(1): 140-156
( 0) |
2. State Key Laboratory of Geo-Information Engineering, 1 Mid-Yanta Road, Xi'an 710054, China;
3. Xi'an Surveying and Mapping Institute, 29 East-Nan'erhuan, Xi'an 710054, China