利用 Sentinel-1 InSAR 推断天津-廊坊地区的减速地面沉降和地下水储存动态

原文：Inferring decelerated land subsidence and groundwater storage dynamics in Tianjin–Langfang using Sentinel-1 InSAR

作者：Xuguo Shi, Tongtong Zhu, Wei Tang, Mi Jiang, Houjun Jiang, Chen Yang, Wei Zhan, Zutao Ming & Shaocheng Zhang

To cite this article: Xuguo Shi, Tongtong Zhu, Wei Tang, Mi Jiang, Houjun Jiang, Chen Yang, Wei Zhan, Zutao Ming & Shaocheng Zhang (2022) Inferring decelerated land subsidence and groundwater storage dynamics in Tianjin–Langfang using Sentinel-1 InSAR, International Journal of Digital Earth, 15:1, 1526-1546,

DOI: 10.1080/17538947.2022.2122610

To link to this article: https://doi.org/10.1080/17538947.2022.2122610

摘要：为满足日益增长的社会经济发展需求，全世界都在从承压含水层中抽取大量地下水。华北平原在过去六十年间经历了相当严重的地下水枯竭和沉降。在本研究中，我们利用 2015 年至 2020 年的 Sentinel-1A/B SAR 图像绘制了天津-廊坊地区的地面沉降图。确定了以广阳、武清-霸州和静海为中心的三个地面沉降带，其最大沉降速率分别为 98.1、121.8 和 104.7 mm/年。利用连续小波变换从时间序列沉降和水力测量中分离出季节信号和长期信号，以检索含水层参数。长期沉降与指数衰减模型非常吻合，但在我们的研究区域明显减缓。弹性骨架储存系数介于 0.52×10-3 和 9.66×10-3 之间。然后，我们检索了地下水总储量、可恢复地下水储量和不可逆地下水储量的时空变化。地下水储量枯竭率明显下降，这得益于南水北调工程的运行和当地的地下水管理措施。

 一、介绍

        近年来，随着社会经济的发展，由于不可持续地开采地下水而导致的土地沉降现象已在世界各地的特大城市广泛出现，例如北京（Du 等，2018 年）、天津（Tang 等，2021 年）、雅加达（Chaussard 等，2013 年）、科尼亚（Orhan，2021 年）、马什哈德（Khorrami 等，2020 年）和河内（Bui 等，2021 年）。作为中国的政治、经济和文化中心，华北平原对工业、农业和生活用水的需求量很大。

        地下水占华北平原总供水量的 70% 以上（Qin，2019 年）。自 20 世纪 70 年代以来，地下水被大量抽取，在华北平原形成了大规模的锥形洼地（guo等人，2020 年）。1984-2014 年期间，在国家重点工程中东部城市观测到水头下降 40-110 米（Yang 等，2021 年）。地下水的快速下降降低了孔隙压力，增加了含水层的有效应力，导致含水层骨架变形。因此，国家地质公园已成为最大的沉降影响区，沉降速度很快（Guo 等，2017 年；Zhang 等，2014 年）。1995-2007 年间，北京的最大累计沉降量为 1.1 米（Zhang 等，2014 年）。自 1958 年以来，天津平原已发现累计沉降大于 1 米的地区面积超过 4000 平方公里（Wang 等，2010 年）。当地政府出台了各种限制和禁止开采地下水的政策（Tang 等，2021 年；Guo 等，2017 年），南水北调中线和东线工程也已投入运行（Dong 等，2021 年；Ha 等，2020 年）。许多地方的地下水位已趋于稳（Long et al.）。

        

        干涉合成孔径雷达（InSAR）是以毫米级精度鉴定区域形变时空变化的有效工具（Xiao 等，2022 年；Jiang、Zhao 和 Shi，2022 年；Shi 等，2022 年）。利用多波段合成孔径雷达数据集（如 C 波段 ERS1/2 和 C 波段 ERS1/2），多时相合成孔径雷达分析方法已被用于描述北太平洋或平原主要城市的沉降特征。C 波段 ERS1/2（Zhang 等人，2016 年）、Envisat ASAR（Jiang、Lin 和 Zhao，2010 年；Perissin 和 Wang，2011 年）、RADARSAT-2（Zhou 等人，2019 年；Zhou 等人，2016 年）、Sentinel-1（Dong 等人，2021 年；Shi 等人，2020 年；Tang 等人，2021 年）。2020；Tang 等 2021）、X 波段 TerraSAR-X（Luo 等 2014；Qin 等 2017）和 L 波段 ALOS PALSAR（Luo 等 2014；Liu 等 2016）数据集。地下水开采导致的地面沉降测量值与地下水储量（GWS）变化和含水层参数有关（Chaussard 等，2014 年；Hoffmann、Galloway 和 Zebker，2003 年）。因此，在沧州和北京平原，InSAR 测量到的沉降也与水头测量一起用于检索 NCP 中精细尺度的 GWS 变化（Jiang 等，2018 年；Zhang 等，2022 年；Bai 等，2022 年）。此前在天津-廊坊地区的研究主要集中在地面变形监测（如 Luo 等 2014 年；Liu 等 2016 年；Zhang 等 2016 年；Zhou 等 2020 年；Tang 等 2022 年）和含水层参数估计（如 Zhang 等 2019 年；Liu 等 2016 年）。与此同时，在 NCP 中观测到了区域性 GWS 枯竭的减少（Long 等，2014 年；Liu 等，2016 年；Zhang 等，2016 年；Zhou 等，2020 年；Tang 等，2022 年）。

        在本研究中，我们首先通过 Sentinel-1 InSAR 分析绘制了 2015 年 2 月至 2020 年 9 月期间天津-廊坊地区的时空沉降图。利用连续小波变换 (CWT) 分离出沉降的季节和长期趋势。然后，利用长期趋势估算含水层衰减系数，并利用季节性沉降和水力测量结果估算骨架贮存系数。根据 InSAR 推算的沉降图和水力测量结果，检索并分析了研究区域的细尺度 GWS 变化。最后，我们讨论了 SNWTP 对地下水变化和地面变形的影响。

二、研究区和数据

        国家冲积平原主要由黄河、淮河、海河和滦河冲积而成。根据沉积的不同，可将其划分为山前平原、冲积扇和洪积平原以及滨海平原（Chen 等，1996 年）。第四纪冲积沉积厚度为 200-600 米，主要由卵石、砾石、砂、粉砂和粘土组成（Guo 等，2015 年）。平原中东部第四纪地层中的含水层系统包括四个含水层组，由粉砂或粘土层分隔（Guo 等，2015 年）。第一和第二组为非承压或半承压含水层，第三和第四组为承压含水层。第一、第二、第三和第四含水层组的下界深度分别为 40-60、120-170、250-310 和 350-550 米（Zhang 等，2009 年；Yang 等，2021 年）。

        我们的研究区域（图 1）位于天津市和廊坊市交界处。该地区深层地下水开采和沉降历史悠久。1973 年和 1997 年，第三含水层组和第四含水层组的地下水凹陷锥体被发现，此后，凹陷锥体从城市地区进一步扩大到农村地区（Ha 等人，2020 年）。第三含水层组和第四含水层组的最大地下水位降幅在 1980-2015 年间分别为 60 米和 1997-2015 年间为 30 米。天津的地面沉降可追溯到 20 世纪 20 年代（Guo 等，2015 年），由于地下水超采，天津的地面沉降逐渐严重，地面沉降速率高达 15-170 毫米/年（Tang 等，2021 年；Zhang 等，2019 年；Zhao 等，2019 年；Guo 和 Zhan，2020 年）。最近，南水北调中线工程为天津市和廊坊市的居民供水（Tang et al 2021）。天津市区的水头以每年 1-2 米的速度恢复（Yang et al 2021）。

       

图1：我们研究区域的位置；插图为国家核电厂以及国家水利水电工程局中线和东线的位置。灰线为不同类型平原的边界。

        

        1978 年，廊坊市城区发现了一个面积为 18 平方公里的地下水凹陷锥体，随后由于地下水的快速抽取，凹陷锥体面积大幅扩大（Yi 等，2005 年）。2001-2015 年期间，深承压含水层的地下水年开采率为 0.33 米/年（Hou 和 Li，2020 年）。直到 2015 年，地下水凹陷锥面积达 300 平方公里（Hou 和 Li，2020 年）。沉降速率从 1979-1983 年间的 5 毫米/年（Yi 等，2005 年）增加到 2001-2015 年间的 ∼ 80 毫米/年（Hou 和 Li，2020 年；He 等，2012 年）。

        我们的研究区域属于四季分明的半干旱和半湿润大陆性季风气候。年降雨量约为 500-600 毫米，其中 75% 集中在夏季（Guo 等，2015 年；Yang 等，2021 年）。然而，年蒸发量约为年降雨量的两倍（Yang et al 2021）。永定河、南运河、北运河和子牙河（图 1）是海河的支流，通过降雨进行季节性补给。

        

        我们共收集了124幅Sentinel-1A/B合成孔径雷达图像，这些图像是在路径47的下降轨道上获取的，时间跨度为2015年2月至2020年9月。这些图像以逐行扫描地形观测（TOPS）模式获取，空间分辨率为 20 × 5 m。选取 2017 年 12 月 21 日获取的合成孔径雷达图像作为参考图像。每幅图像最多允许有三个时间基线小于 120 天的连接来生成干涉图（图 S1）。站点（117.17°，39.1°）的每日降水记录用于分析影响因素。我们在 2019 年 4 月和 2020 年 9 月期间在承压含水层的 7 口观测井（图 1 中的 Q1-Q7）获得水头测量数据，以分析含水层参数。在 2016-2019 年期间，我们在三个连续 GPS 站（KC02、CH01 和 JHAI）获取了地面变形测量数据，以验证我们的 InSAR 变形结果。这三个站点的原始 GPS 测量数据使用 GAMIT/GLOBK 10.40 软件（Herring、King 和 Mcclusky，2010 年）进行处理，并以 2008 年 ITRF 为参照基准。我们使用 ALOS World 3D-30 m (AW3D30) 数字高程模型生成差分干涉图和地理编码。

 三、方法

        在本节中，我们首先通过时间序列 InSAR 分析得出时间序列变形。然后，我们从得出的形变中分离出季节形变和长期形变，以检索含水层参数，这对 GWS 变化估算非常有用。工作流程见图 2。

        3.1 Insar数据处理

        多普勒中心频率在 TOPS SAR 图像的每个突发中都有很大差异；因此，我们采用了一种增强型光谱分集方法来避免方位角配准误差（Jiang，2020 年；Jiang 等，2017 年）。我们使用 GAMMA 软件对所有突发图像进行配准并拼接成无缝图像（Wegmüller 等，2015 年）。使用 AW3D30 DEM 去除地形相位。为了提高低相位稳定性分布式散射体的信噪比，采用了同质像素自适应滤波器（Jiang 等，2015 年；Zhao 和 Jiang，2022 年）。图 S2 中给出了两个以 SAR 图像坐标叠加振幅的干涉图作为示例。

        我们最初使用宽松的振幅离散阈值（<0.6）来包含候选像素。然后使用时间一致性阈值（>0.3）来排除噪声候选像素，以便进行时间序列分析（Jiang 和 Guarnieri，2020 年）。对所选像素进行空间二维相位解包。我们使用双线性模型去除每个干涉图中的相位斜坡。通过线性拟合相位和垂直基线，估算并去除每个像素中与高程残差相关的相位。使用 InSAR 产品的通用大气校正在线服务（Yu、Li 和 Penna，2018）对每个干涉图上的大气信号进行校正。

图2：我们的研究流程图

        然后通过奇异值分解对时间序列变形进行反演。先前的研究表明，该区域的水平变形极小（Tang 等，2021 年）。因此，我们通过除以入射角的余弦值，直接将视线变形投影到垂直方向。

        3.2 含水层参数反演

        为进一步检索含水层参数，可将 InSAR 时间序列变形和水头测量值分解为季节和长期趋势信号。小波变换是处理自然非稳态信号的有用工具（Torrence 和 Compo，1998 年）。小波变换可以从时间序列信号中分离出不同时间和频率的频率特性（Miller 和 Shirzaei，2015 年）。在本研究中，我们利用 CWT 从时间序列变形和水头中分离出季节成分。我们首先通过假设一个恒定的速率来去除线性趋势（Jiang 等人，2018 年）。然后，将时间序列残差分解为一系列具有不同时间和频率尺度的信号。通过反滤波器，使用时间尺度在 0.5 至 1 年之间的信号重建短期信号。振幅是年内最大值与最小值之差的一半。然后通过从原始时间序列信号中减去短期季节信号来获取长期趋势。

         长期变形与地下水长期下降造成的承压含水层非弹性压实有关。低导水率的承压含水层的压实通常需要几十年甚至上百年的时间才能达到与地下水位的延迟平衡。根据水动力固结理论，长期变形可以用时间的指数函数来表征（Hu、Lu 和 Wang，2018 年）。

        其中，dl(t) 为 t 时的长期变形量，M 为含水层响应的幅度系数，K 为[-1，0]之间的衰减系数。M 和 K 采用最小二乘法估算。无量纲弹性骨架储水系数反映了含水层或含水层对水头短期变化的响应，在地下水管理中起着至关重要的作用。它描述了含水层在不引起非弹性变形的情况下，因水头变化而储存（或排泄）的水量（Miller 和 Shirzaei，2015 年；Zhang 等人，2019). 含水层系统的弹性骨架储水系数（Ske）可通过下式求得：

        其中，Δds 和 Δhs 分别为季节性弹性垂直变形和水头变化。一般来说，季节性垂直变形与水头变化之间存在时滞，可通过相关分析方法确定（Jiang 等，2018 年）。利用估算出的恒定时滞来校核季节变形。同时，对季节变形进行线性插值，以获得与水头观测结果同步的变形。弹性骨架储水系数采用最小二乘法求得。

        3.3 GWS 估算

        地下水总储量（TGWS）的变化等于承压含水层系统中含水层的压实/膨胀（Jiang 等，2018 年），可以通过 InSAR 获取的垂直变形进行估算。这一假设在我们的研究中是合理的；天津的沉降主要是由 300 米以下含水层系统的压实造成的（Guo 等，2017 年）。

其中，ΔV 为 TGWS 变化量，Δdt 为总垂直变形量，A 为含水层面积。TGWS 可分为可恢复（ΔVr）和不可逆（ΔVi）GWS。可恢复地下水储量（RGWS）可在含水层得到补给时恢复。可利用弹性储水系数和水头变化对其进行估算。

其中，Δht 为总水头变化。不可逆地下水储量（IGWS）可通过去除 RGWS 从 TGWS 中估算出来。

4 结果与分析

4.1 下沉与季节性变形

        图 3（a,b）分别给出了哨兵 1 号数据集显示的 2015-2020 年研究区域的下沉率和季节变形幅度。图 S3 给出了视线方向的原始位移率。共探测到 6,226,271 个测量点，密度为 ∼798 个点/平方公里。在图 3（a）中，我们确定了以广阳、武清-霸州和静海为中心的三个大型沉降区。下沉速率超过 30 毫米/年的区域用多边形勾勒，面积分别为 357.1 平方公里、1,470.8 平方公里和 317.9 平方公里。沉降面积和季节变形的分布并不完全相同。在农业区发现了明显的季节性振幅（10-15 毫米）（图 3(b)）。灌溉和降水导致的地下水季节性变化引起了这些季节性变形（Zhou 等，2020 年；Tang 等，2022 年）。

        廊坊市广阳区的最大下沉率在 2003-2007 年间为 65 毫米/年（Ge 等，2008 年），在 2015-2020 年间增加到 98.1 毫米/年（图 3(a)）。虽然这两个时期的下沉速率接近，但下沉范围从北部扩大到几乎整个廊坊市区。图 3（c）给出了 AA'剖面的累积沉降。廊坊市城区的地下水开采量自 2009 年开始趋于稳定（Hou 和 Li，2020 年），并在国家污水处理厂运行后有所下降，沉降速率约为 -60 毫米/年。目前，沉降主要集中在农业区和新开发的居民区。2015 年和 2020 年的累计位移分别达到 468.1 毫米和 502.3 毫米。

图3：变形图和剖面图。(a) 长期沉降率和 (b) 2015-2020 年期间的季节性变形幅度，(b-d) AA'、BB'和 CC'剖面的时间序列沉降。a) 中的三角形为参考点。

        天津市武清区和廊坊市霸州市是两个重要的第二产业基地。大量的深层地下水开采导致了严重的地面沉降，在 2005-2006 年期间，最大沉降速率超过了 200 毫米/年（Zhang、Li 和 Xia，2008 年）。两个明显的地面沉降中心分别位于武清王庆坨镇和霸州市胜芳镇，前者的地面沉降速率为 90 至 122 毫米/年，后者的地面沉降速率为 80 至 100 毫米/年。如图 3（d）所示，在 2015-2020 年期间，剖面 BB'的王庆坨和胜芳沉降中心的累计沉降量分别为 609.5 毫米和 651.7 毫米。

        静海大部分地区的沉降速率小于 60 毫米/年，而新建居民区的最大沉降速率为 104.7 毫米/年。图 3（e）中的 “剖面 CC ”显示，2015-2020 年间最大累计下沉量为 551.1 毫米。以永清县曹家坞镇和固安县牛坨镇为中心的两个小规模沉降区的沉降速率约为 60 毫米/年。造成这种情况的原因是与灌溉有关的地下水抽取。

4.2 InSAR 与 GPS 测量结果的比较

        我们将 GPS 的三维位移投影到哨兵-1 数据集的 LOS 方向。对 GPS 观测站 50 米范围内的 InSAR 测量值取平均值，并与 GPS 测量值进行比较（图 4）。KC02、CH01 和 JHAI 的累积 LOS 变形量分别为 -12.2、-23.5 和 -165.1 毫米。两种技术的测量结果非常吻合，从而验证了 InSAR 变形测量结果的可靠性。在 KC02、CH01 和 JHAI，InSAR 和 GPS 在相同日期的测量值差异的平均值和标准偏差均小于 14 毫米（图 4）。InSAR 和全球导航卫星系统形变时间序列之间的差异可能是由于形变表示方法造成的，InSAR 反映的是 50 米范围内的形变，而全球导航卫星系统反映的是全球导航卫星系统支柱的局部形变（Bui 等，2021 年）。

4.3 长期衰减系数

        图 5 给出了衰减系数 K 和振幅系数 M。指数衰减模型的均方根误差（RMSE）见图 S4。衰减系数描述了地面沉降对水头变化的长期延迟响应。M 的负值和正值分别表示地面的隆起（含水层恢复）和下沉（含水层压实）（图 5(b)）。所有像素的均方根误差值均小于 30 毫米。在本研究中，为确保稳健性，对振幅系数小于 10,000 毫米的像素（占 4b 中给出的像素总数的 88.4%）进行了分析。K 值越小，意味着对水头变化的响应速度越快，含水层系统的水力动态达到平衡的时间越短。以 D1-D5 为中心的衰减系数为-0.5 至-0.15，明显低于广阳和武清-霸州其他地区，表明对水头变化的响应速度较快。静海的平均衰减系数为-0.1，没有明显的群集低值。天津市宜兴埠镇（以 D6 为标志）的衰减系数小于-0.5，表明该地区含水层系统的导水性较高。宜兴埠镇在 20 世纪 70 年代至 90 年代期间是一个沉陷盆地（Ha 等人，2020 年），2015-2020 年期间的沉陷速率为 5.4 毫米/年（图 3(a)）。

        图 6 给出了图 5 中标注的六个典型地点的长期下沉情况及其使用指数衰减模型的建模参数。长期沉降可以完全模拟。这些点的沉降幅度从 85.9 毫米到 698.9 毫米不等，衰减系数为-0.6 到-0.2。所有选定点均处于减速过程中。廊坊和天津的用水量正逐渐被从南水北调中线调水所取代。地下水开采或回采量减少导致地面沉降减速。对长期水头变化的延迟响应的时间尺度大致可以从衰减系数的倒数的绝对值来估算（胡、卢和王 2018 年）。

        对长期水头变化的延迟响应的时间尺度大致可从衰减系数的倒数的绝对值来估算（胡、卢和王，2018 年）。D1 至 D6 的估算时间尺度分别为 3.3 年、3.2 年、5.3 年、2.0 年、3.8 年和 1.6 年。图 4 中较大的均方根误差可能与异常加速度或线性位移趋势有关（胡、卢和王，2018 年）。D3 及其周边区域的均方根误差（>20 毫米）与王庆坨水库（WOTR）建设导致的变形趋势的突然变化有关，将在第 5.1 节进一步讨论。

图4：InSAR 与 GPS 测量结果的比较

图5：使用指数衰减模型计算的长期沉降的 (a) 衰减系数 K 和 (b) 幅值系数 M

 图6：中标注的 D1-D6 的长期沉降以及用指数衰减模型拟合的沉降趋势

 4.4 水头和弹性骨骼储存系数

        图 7 给出了 Q1 至 Q7 点（图 1 中标出）的地面沉降时间序列及其相应的水头变化。Q2 和 Q7 点的水头在下降，这与这两个地点的地面沉降一致。当水头急剧下降时，含水层发生了压实，然后出现了地面沉降。虽然 Q1、Q3、Q4 和 Q6 的水压水头在 2017-2020 年期间显示出明显的恢复趋势，但沉降仍呈减速趋势。这一发现表明，这些地点目前的水头低于固结前的水头（Bai 等，2022 年）。图 S5 给出了变形和水头的季节性变化。在所有水头记录中，季节性变化都很明显，振幅在 0.2-1 米之间，而与长期沉降相比，季节性地面变形更为微小（小于 10 毫米，图 7 和图 S5）。通过相关分析（表 1）计算出了季节运动与水力变化之间的时滞。时滞范围为 12-204 天，相关性大于 0.43。不同的时滞与含水层系统垂直导水性的变化有关（Hoffmann、Galloway 和 Zebker，2003 年）。Q6 的时滞仅为 12 天，这表明含水层的渗透和排水速度快于其他地点。

         表 1 列出了 Q1-Q7 的弹性蓄水系数。Q1-Q6 位于冲积扇和冲积平原，而 Q7 位于沿海平原（图 1）。Q1-Q6 处的弹性骨架储量系数介于 0.52 × 10-3 和 9.66 × 10-3 之间，与之前抽水试验得出的数值一致，例如 1.0 × 10-3 -8.0 × 10-3 （Zhang 等，2009 年）和 2.0 × 10-6 -4.0 × 10-3 （Yang、Meng 和 Li，2015 年）。沿海平原 Q7 的系数为 4.85 × 10-3，与天津市滨海新区的系数（2.0 × 10-4 -5.0 × 10-3）一致（杨、孟、李，2015 年）。

        我们对表 1 中的系数进行了反距离加权插值（图 8），分辨率为 ∼ 380 米。广阳区的弹性储水系数变化显著。广阳沉陷中心南部和东北部的含水层分属永定-九马承压含水层系统的两个不同子系统（Chen，1999 年）。子系统中不同的水文地质条件，如产水量特性和含水层厚度，可能导致系数的变化。南部地区的系数最大，这表明在我们的研究区域，给定水头变化所产生的地下水交换量最大。

4.5 TGWS 变化

        利用哨兵-1 号卫星得出的年度累计下沉量来估算 2015-2020 年期间 TGWS 的变化（图 9）。

 表1：估计时滞、最大相关性和弹性存储系数

 

 图7：承压含水层 Q1-Q7 点的时间序列沉降与水头变化

         等效水厚（EWT）表示单位含水层面积的总含水量变化量。我们将总含水量（TGWS）、可恢复含水量（RGWS）和不可恢复含水量（IGWS）的分辨率统一为弹性骨架储水系数的相同分辨率。含水层面积可通过分辨率单元面积的总和来估算。为了比较年度结果，我们将 2015 年 2 月至 12 月以及 2020 年 1 月至 9 月获取的沉降量转换为年度累计沉降量，采用简单的比例缩放（d/T×365）方法。d 和 T 分别为 2015 年和 2020 年检索到的下沉量和观测天数。2016 年 TGWS 平均耗水率达到 46.6 毫米/年，2020 年降至 15.8 毫米/年。自上海污水处理厂中线和东线开始运行以来，大量的水被输送到我们的研究区域，从而大大减少了深层地下水的消耗。2014-2018 年期间，次区域地下水枯竭有所缓解（Zhang et al 2021）。2016 年，广阳、武清-霸州和静海的 TGWS 平均消耗量分别为 78.7、71.4 和 60.9 毫米。此后，主要消落区的 GWS 消耗有所减缓（图 9）。2020 年，广阳的 EWT 降至 20.0 毫米（图 9（f））。SNWTP-M 经过的武清-霸州地区的耗损自 2019 年起也逐渐缓解（图 9）。其中，2019 年和 2020 年，武清区污水处理厂周围的 TGWS 消耗量分别为 43.7 毫米和 4.8 毫米。2017-2020 年，静海的 TGWS 损耗率为 35.4 毫米/年，主要集中在静海东部新开发区域。

 

 图9：(a-f) 2015-2020 年期间 TGWS 的年消耗量变化

 4.6 遥感预警系统（RGWS）和综合预警系统（IGWS）的时空变化

        利用地下水水头和弹性骨架储水系数（表 1），2018-2020 年期间 RGWS 和 IGWS 的时空变化见图 10（a-f）。2018-2020 年间，RGWS 消耗率分别为 4.4、2.9 和 3.5 毫米/年，空间变化显著。2018 年广阳沉陷区的 RGWS 平均损耗量为 5.9 毫米，2019 年和 2020 年分别逐渐恢复到 8.6 毫米和 3.2 毫米。RGWS的恢复与2019年以来逐步开展的水资源置换有关。武清和霸州的 RGWS 出现了不同的变化。武清的 RGWS 在 2018 年和 2020 年逐渐恢复，恢复速率分别为 12.9、3.3 和 5.3 毫米/年。然而，霸州的 RGWS 变化不大，在 2018-2020 年期间变化范围为 0.5 毫米至 4.1 毫米。造成这种情况的原因可能是廊坊和天津的水资源管理策略不同。此外，霸州地区水头测量数据稀少，也导致 RGWS 估算值存在较大的不确定性。静海沉陷中心的 RGWS 损耗率最小，分别为 3.7、0.5 和 2.9 毫米/年。

        在 2018-2020 年期间，IGWS 的平均损耗率分别为 41.5、22.6 和 19.3 毫米/年。与 TGWS 的变化相似，IGWS 的损耗也有所减轻，这与西南网络图有关。此外，WQTR 附近的 IGWS 损耗（图 10（e,f））在 2019 年为 45.1 毫米，2020 年降至 6.4 毫米，表明 IGWS 损耗显著减少。

5 讨论

5.1 季节性变形的变化

        选取了位于廊坊和天津交界处的一个典型区域来分析季节性变形的变化情况。南水北调西线工程于 2018 年底完工，作为南水北调东线工程的配套基础设施，以保障天津城市和工业用水的稳定供应。

        

 图10：年度（a-c）RGWS 和（d-f）IGWS 的时空变化

        图 S6 给出了 2021 年 10 月 21 日获取的相应光学图像。如第 4 节所述，自 2019 年以来，TGWS 和 IGWS 的耗水量已大幅减少。水质监测站的蓄水为周边地区的含水层提供了充足的水补给。与此同时，王庆坨和东沟岗的地下水消耗也被调水所替代。水头可能会迅速恢复到接近固结前的水头。因此，含水层的压实和地面沉降也会大大减缓（如图 12 中的 W1）。图 12 中的 W1）。天津市王庆坨镇和廊坊市东沽港镇出现异常季节变形（图 12（b），13 毫米）（图 11（f））。该地区的季节信号被高估。当位移无法用指数衰减模型（如图 S4 中 D3 处的较大均方根误差）或线性模型描述时，就会出现较大偏差，尤其是在 2019 年 8 月之后的测量中。不准确的季节信号将通过偏离的残差来构建。不过，由于 Q4 和 Q5 井位于该区域之外，因此异常季节变形不会影响骨架储存系数的可靠性。W1 的季节变形时间序列显示，2019 年之前的振幅小于 5 毫米。外青坨镇的季节变形周期并不恒定。最大下沉在 2 月至 10 月达到最大值（图 12(b)）。

         与王庆坨镇不同，廊坊的盂城镇属于农业区，在 2015-2020 年期间有明显的季节性变形（图 11）。重复周期约为一年（如图 12（c,d）中的 W2）。年振幅在 10-15 毫米之间，这与灌溉和季节性降水有关。冬春灌溉抽取地下水会导致下沉量逐渐增加，并在 7-8 月间达到最大值。然后，含水层和含水层在降雨的自然补给下出现弹性回升。自 2019 年起，灌溉用水逐渐从地下水转变为从新农水厂调水（张，2020 年）。农业区的振幅逐渐减小到 5 毫米以下（图 11 和 12(d)）。

 图11：王庆坨镇和葛渔镇的年季节变形

 

 图12：图 11 标出了 W1 和 W2 的累积沉降和季节变形

 5.2 优势和局限性

        InSAR 的密集测量点可以提供精细的地面沉降图，以探测潜在的风险。Sentinel-1 或未来的 NISAR 任务的重访周期短，使我们能够精确捕捉地面的动态变化。通过密集采样测量，我们可以更好地建立长期和季节变形模型。变形图和水头使我们能够检索衰减系数、骨架储存系数和 GWS 变化，这对地下水管理和降低风险非常重要。

        我们假设地面沉降都与承压含水层系统中的含水层压实/膨胀有关。然而，城市引起的负荷、软土压实（Shi 等，2021 年）和其他原因也会造成地面沉降。因此，可能存在高估 GWS 变化的情况，尤其是在快速发展的地区。同时，空间模式与地质断层（Hu 等，2019 年）和含水层结构有关，我们将在今后的研究中对此进行调查。我们忽略了水平变形，这导致了估算的 GWS 的不确定性。结合升轨和降轨 InSAR 测量可以获得更精确的沉降图。水力测量的密度和分布限制了我们对 RGWS 和 IGWS 的估算，例如我们研究区域内的霸州。可以使用更密集的观测井测量数据或重力恢复和气候实验后续数据来降低不确定性。

6 结论

        在本研究中，我们通过时间序列 Sentinel-1 InSAR 分析绘制了 2015-2020 年期间天津-廊坊地区的时空沉降图。测得的地面沉降与 GPS 测量结果进行了验证。孤立的长期和季节信号使我们能够获得衰减系数和骨架存储系数。长期沉降可以用指数衰减模型很好地描述，这表明我们研究区域的沉降正在减速。水头测量结果和 GWS 变化证实了次区域地下水的恢复。从图 9 和图 10 中的 TGWS、RGWS 和 IGWS 测量结果可以推断，研究区域的地下水开采似乎更具可持续性。这种情况与一系列地下水勘探限制和上海污水处理厂运营政策出台后地下水开采量减少有关。蓄水后，东姑岗镇和王庆坨镇的水头迅速恢复，并改变了其位移模式。我们的研究结果证明，InSAR可以提供有关地面沉降和地下水变化动力学的有价值的信息。