南水北调中线区域地面沉降SBAS-InSAR监测研究

2022-11-30张双成宋明鑫张彬玲樊茜佑司锦钊张雅斐

大地测量与地球动力学 2022年12期

张双成余静宋明鑫张彬玲樊茜佑司锦钊张雅斐

1 长安大学地质工程与测绘学院，西安市雁塔路126号，710054 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需满足：

(1)

设初始时刻为t0，任意时刻ti(i=1,…,N)相对于初始时刻的相位差是未知参数φ(ti)，干涉处理获取的差分干涉相位δφ(tk)(k=1,…,M)为观测量。若不考虑其他相位误差，则第k景差分干涉图的像元(x，y)地表形变相位δφdef可表示为：

δφdef=φ(tB,x,y)-φ(tA,x,y)=

(2)

式中，tA、tB(tA

将所有差分干涉图进行自由组合，φ为N幅SAR影像上的干涉相位值组成的矩阵，δφ为M幅差分干涉图上相位组成的矩阵。则有：

(3)

(4)

将上式写成矩阵形式：

Aφ=δφ

(5)

式中，A为M×N维矩阵。

(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

2 研究区域及实验数据来源

2.1 研究区域概况

南水北调中线工程输水干渠源自丹江口水库，流经河南、河北，最后到达北京、天津，总长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数据去除地形相位。

3 中线区域地面沉降时空演变特征分析

由于研究区域较大，需进行数据分块处理。将每景影像经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

4 重点区域时序特征分析

4.1 北京市沉降结果及验证

为探究区域地面沉降在南水北调中线工程通水后的时空演变特征，对北京市的地面沉降进行重点分析。北京市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监测结果精度较高。

4.2 地面沉降时空特征分析

为揭示北京市地面沉降形变特征和演化规律，获取北京市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个特征点呈不同速度的非线性沉降趋势，经历了从快速到逐渐减缓的沉降形变过程，由此可知，北京市的地面沉降正在减缓。

表2 InSAR结果与GPS观测值比较

图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。

为进一步分析地面沉降速率随时间的演化过程，选取地面沉降分布区域内的典型特征点。特征点的选取基于局部沉降区域中较为严重的沉降中心，即图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个点的沉降速率均减小，一定程度上反映了北京市的地面沉降程度正在减缓。

4.3 地面沉降与地下水演化相关性分析

地下水的过量开采会使含水层孔隙度逐渐变小，导致有效应力增加，从而引发地面沉降。北京市由于水量短缺及水体污染导致可用的地表水所剩无几，不得不长期超采地下水，从而引发大范围的地面沉降。众多学者结合北京市的地面沉降现状，分析了其与地下水演化的相关性[8-9]。

图10 北京市2016～2020年LOS向形变速率Fig.10 Mean deformation velocity in LOS of the Beijing from 2016 to 2020

图11 主要沉降区域特征点形变速率对比Fig.11 Comparison of deformation velocity of feature points in main subsidence areas

北京市密云水库长期担负着城市生活及农业生产用水的重要任务，仅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，彻底扭转了供水量入不敷出的局面。因此，由北京地下水水位与地面沉降的对应关系可知，南水北调工程与该区域的地面沉降缓解密切相关。