面向 InSAR 时间序列的基于奈奎斯特的多步校正与平滑框架：黄河三角洲案例研究

《Estuarine, Coastal and Shelf Science》：A Nyquist-based Multi-step Correction and Smoothing Framework for InSAR Time Series: Case Study of the Yellow River Delta

　　

•多级阈值/波峰-波谷校正抑制InSAR时序异常
•聚类和PCA揭示8年土地沉降异质性，沉降速率减缓且范围缩小
•人类活动变化是时空沉降变异性变化的主要驱动力

1. 引言

土地沉降目前已成为全球广泛存在的环境与地质灾害。在中国，每年用于防治沉降的支出高达15亿美元（Herrera-García et al., 2021）。河口及沿海地区通常人口稠密，经济价值高。然而，由于地质条件及人类活动的影响，这些地区经历了严重的显著沉降。除了对人工结构的破坏性影响外，沿海地区的土地沉降还会通过与极端天气和海平面上升的相互作用，显著增加沿海地区自然灾害的风险（Nicholls et al., 2021; Wu et al., 2022）。一方面，沿海地区的土地沉降会导致沿海堤坝和排水站等工程设施的沉陷和开裂，从而导致沿海地区防洪和防潮能力的下降（Gracia et al., 2019; Wang et al., 2020）。另一方面，它会显著增加局部地区的海平面高度（Nicholls et al., 2021; Tay et al., 2022; Li et al., 2023），加剧沿海洪水、土壤盐渍化、海岸侵蚀和湿地退化的风险（Syvitski et al., 2009; Shirzaei et al., 2021）。这些效应削弱了沿海地区应对海洋灾害的能力，对这些地区的人类安全和环境保护构成了更为严峻的挑战（Nicholls and Cazenave, 2010; Cazenave and Cozannet, 2014; Syvitski et al., 2022; Ao et al., 2024）。

目前已知，土地沉降可由一系列自然因素或人类活动引起。常见的自然因素包括火山活动（Doke et al., 2020）、地震（Ali et al., 2021; Hanebuth et al., 2022; Orellana et al., 2022）、地下溶洞坍塌、干旱（Welch et al., 2024）以及软土的压缩或固结（Steckler et al., 2022）。然而，在人类活动强烈的地区，土地沉降的主要驱动因素通常是人为的，包括抽取地下水（Xu et al., 2012; Yu and Michael, 2019; Wang et al., 2022）、石油（Hu et al., 2020）、煤炭（He et al., 2020）、天然气（Chaussard et al., 2013）和地热资源所导致的地面塌陷和沉积物固结。此外，城市建设及大型工程造成的过度荷载也会引起土地沉降（Herrera-García et al., 2021）。由此可见，土地沉降的驱动因素复杂，受自然环境和人类活动的双重影响。沉降过程可能随着周围条件的改变而发生变化，包括沉降的加速和减速、影响范围的扩大或收缩，以及回弹等趋势。因此，识别这些驱动因素并理解其变化规律，对于分析驱动因素和制定防治措施至关重要。

随着2014年"哨兵-1"（Sentinel-1）卫星的发射以及数据源的持续更新，利用多时相干涉合成孔径雷达（MT-InSAR）研究黄河三角洲（YRD）土地沉降取得了质的飞跃。MT-InSAR的监测结果在覆盖范围和空间采样率方面显著增加，精度也得到了极大提高。

基于Sentinel-1数据和MT-InSAR技术，Peng et al. (2020)和Zhang et al. (2019)详细描述了黄河三角洲近期的土地沉降情况，识别出东北和东南地区的若干显著沉降漏斗。Wang et al. (2022)利用MT-InSAR、主成分分析（PCA）和多源数据定量分析了黄河三角洲土地沉降的驱动因素。Zhang et al. (2024)指出黄河三角洲土地沉降存在显著的时空变化。Huang et al. (2022)也指出，黄河三角洲目前正面临严重的土地沉降威胁，其沉降速率在所有已知沉降三角洲中是最快的，这为理解黄河三角洲土地沉降提供了新的认识。

利用MT-InSAR分析土地沉降的驱动因素已成为重要的研究热点。目前的研究表明，沉积物固结和压缩是黄河三角洲地区广泛但沉降幅度相对较低的土地沉降的主要原因（Zhang et al., 2014; Liu et al., 2015; Zhang et al., 2018; Liu et al., 2021），而更为严重且具有局地性的沉降事件通常归因于卤水开采和石油开采活动（Higgins et al., 2013; Zhang et al., 2019; Peng et al., 2020）。这些发现凸显了沉降模式与其驱动因素之间的强相关性（Raspini et al., 2014; Bui et al., 2021）。此外，一些研究人员尝试利用整合InSAR反演变形数据与地质先验知识的地球物理模型，定量研究沉降源、沉降量级与空间分布之间的关系（Liu et al., 2015; Liu et al., 2016）。总体而言，这些研究为深入了解黄河三角洲土地沉降过程奠定了坚实基础。

尽管小基线 subsets 干涉合成孔径雷达（SBAS-InSAR）已被证明能够提供高精度的变形测量，但其结果通常以长期平均速度的形式呈现。这种表示方式往往掩盖了短期的动态变化。例如，一年内发生的显著季节性沉降可能在多年平均图中被平滑，导致单个平均沉降速率无法反映时间变异性（Martina et al., 2023; Ziyan et al., 2024）。因此，为了更全面地了解土地沉降动态，有必要进一步分析SBAS-InSAR生成的完整变形时序数据。然而，大多数现有研究尚未从年沉降速率的角度充分探讨土地沉降的时空动态。此外，沉降的时空演化模式也未被强调为考察其驱动机制的关键视角。缺乏这种时间分辨率和模式导向的重点关注，重要的变化信号可能未被探测到。

此外，一些研究依赖于有限数量的代表性点来反映区域沉降趋势，从而未能充分利用完整沉降时序中蕴含的丰富时空信息。这种简化可能导致遗漏关键的局地特征。Xue et al. (2020)对时序InSAR技术的综述强调了长期累积变形数据在理解地面运动规律性模式中的价值，这与充分利用时序数据潜力的需求相一致。最后，缺乏连续的原始观测通常会在InSAR时序中引入时间不连续性，使得变形趋势的解释变得复杂。由大气延迟、处理误差或实际地面事件引起的突跳和离群波动等问题，已被广泛认为是时序变形分析中的关键挑战（Doin et al., 2009; Zhiwei et al., 2019）。

因此，在沉降监测和分析中推进InSAR结果的利用仍存在巨大的潜力。针对当前黄河三角洲土地沉降研究中存在的问题，如年际变化趋势不明确、发展模式和空间转移不清晰等，本研究引入了一种专门设计的误差校正算法，用于缓解InSAR变形时序中不合理的跳变。为进一步分析土地沉降的年际变化并考察黄河三角洲的整体时间趋势，本研究采用了聚类和PCA等方法。通过这种方法，本研究旨在为黄河三角洲土地沉降的时空模式提供更全面的理解，并识别这些变化背后的驱动因素。

2. 研究区

黄河三角洲在渤海经济圈中发挥着至关重要的作用，是中国油气和卤水生产的重要基地。除了工业重要性之外，黄河三角洲还分布着广袤的湿地，为候鸟提供必要的栖息地，凸显了其高度的生态重要性。然而，该地区的生态系统尤为脆弱，其地质地层中普遍分布着盐碱土，淡水资源匮乏。由于这些脆弱性，黄河三角洲及其周边地区对土地沉降和海平面上升（SLR）尤为敏感，使减轻土地沉降成为该地区面临的最紧迫的环境挑战之一。

黄河三角洲的地形、土地覆盖分类和沉积物厚度如图1所示。作为一个典型的沿海地区，黄河三角洲以密集的人类活动为特征，包括石油开采、卤水开采、化工行业和农业。从地质学角度来看，该地区底层由黄河搬运的厚层沉积物构成，具有显著的固结潜力——为土地沉降的发展提供了有利条件。对地形图、土地覆盖类型、沉积物厚度和卫星影像的分析表明，自然因素构成沉降的基本前提条件，而人类活动则充当主要驱动力。加强研究方法论、评估土地沉降现状、识别时空模式并研究其驱动机制，对于有效预防和管理黄河三角洲土地沉降具有重要意义。

下载：下载高分辨率图像（1MB）下载：下载全尺寸图像

图1. 黄河三角洲的地形图、地貌图、土地分类图和沉积层厚度图，以及哨兵数据覆盖区域。(a) 黄河三角洲的SRTM数字高程模型。(b) Google地图和Sentinel-1A数据覆盖区域。(c) 土地利用与土地分类（Wang et al., 2024），(d) 沉积物厚度（Zhang et al., 2018）。

3. 数据集与方法

本研究采用Sentinel-1A数据集和MT-InSAR方法推导变形时序和沉降速率。时序中的缺失值使用移动平均（MA）方法进行插补，随后进行异常校正以提高数据可靠性。随后，对校正后的时序应用聚类分析和PCA，以探索时空变形模式。最后，我们考察了黄河三角洲的变形过程及其时空特征。以下详细介绍了数据来源和分析方法。

3.1. 数据集

本研究使用的数据集列于表1中，包括Sentinel-1A数据（图1b）、基于Landsat-8影像的黄河三角洲土地覆盖分类以及SRTM数字高程模型。Sentinel-1数据集源自欧洲航天局（ESA）。尽管Sentinel-1A仍在工作且数据持续更新，但Sentinel-1B已于2021年底退役。因此，本研究使用Sentinel-1A数据，时间跨度从2016年1月9日至2023年11月16日，共包含194景影像。本研究使用的土地覆盖分类结果基于先前发表的研究成果（Wang et al., 2024）。此外，SRTM数字高程模型用于在MT-InSAR处理过程中去除地形相位分量。

表1. 数据集介绍。

| 数据集 | 应用 |
|--------|------|
| Sentinel-1A | 计算土地沉降（2016年1月9日-2023年11月16日，194景） |
| SRTM-DEM | 去除地形相位 |
| 土地覆盖分类 | 辅助识别人类活动和自然因素 |

3.2. MT-InSAR参数与主要操作

合成孔径雷达（SAR）数据具有全天候、全天时获取、空间覆盖广以及高时空分辨率等优势。MT-InSAR作为一种先进的SAR数据处理技术，用于变形监测时可达亚厘米级精度，已成为土地变形分析领域不可或缺的工具。在过去十年中，众多研究应用该技术调查黄河三角洲的土地沉降。

黄河三角洲主要由河流-冲积沉积物组成，农业为主要土地利用方式。除城市化区域外，该地区通常缺乏天然永久散射体。考虑到该区域的相干性特征，本研究采用SBAS-InSAR技术来评估黄河三角洲土地沉降的时空动态。MT-InSAR参数如表2所示。

表2. MT-InSAR参数。

| 参数 | 阈值 |
|------|------|
| 解缠阈值 | 0.4 |
| 参考点坐标 | 37.547983, 118.376349 |
| 解缠算法 | Goldstein |
| 多视 | 20个距离向视 × 5个方位向视 |
| 时间基线 | 12至120天 |
| 最大空间基线 | 170.8 m |

如图S1所示，共通过干涉合成孔径雷达科学计算环境（ISCE）生成了960幅原始干涉图，空间分辨率约为100米（20个距离向视 × 5个方位向视）（https://github.com/isce-framework/isce2）。其中，950幅干涉图被保留用于后续SBAS-InSAR处理，10幅因低相干性等质量问题被排除。2016年至2023年的平均变形速率通过MintPy软件（Yunjun et al., 2019）（https://github.com/insarlab/MintPy）推导得出。

需要注意的是，MintPy揭示的形变被重采样至90 m的空间分辨率，且所得形变速率为视线（LOS）方向的形变速率。由于研究区域水平运动不显著，我们假设观测到的LOS位移主要归因于垂直形变。因此，可根据公式（1）从LOS速率中分解出垂直形变速率。

$$V_{ver} = \frac{V_{LOS}}{\cos\theta} \tag{1}$$

式中，$V_{ver}$为垂直形变速率，$V_{LOS}$为LOS方向形变速率，$\theta$为每个像素的入射角。

### 3.3 时间序列空隙填充

基于12天的获取间隔，2016年1月9日至2023年11月16日之间预期共有234景Sentinel-1影像。然而，由于不可避免的数据缺失，最终仅获得194景影像用于分析。如图S2所示，特定时段存在大量数据缺口。为保持分析的连续性并适应滤波算法的时间依赖性，本研究采用窗口大小为3的移动平均（MA）算法对缺失数据点进行插值。需要指出的是，由于连续性验证的局限性，2016年的插值结果可能可靠性较低。

### 3.4 基于奈奎斯特定律的时间维度多步校正与平滑（NMCS）

区分时域校正与空间相位解缠至关重要。空间解缠算法（如MCF、SNAPHU）可以利用空间连续性绕过相位不连续区域，从而有效恢复空间域中的大形变梯度。然而，时间采样稀疏且离散（Sentinel-1重访间隔为12天）。当两个相邻时间点之间的形变超过λ/4时，单幅干涉图的相位缠绕会导致绝对模糊信息的丢失。尽管如此，SBAS等多时相算法可以利用冗余观测网络在时间维度上部分约束这种模糊性，从而仍能估计形变量。然而，由于相位解缠误差的引入和反演问题的病态性，估计结果的精度和可靠性显著降低，表现为形变误估和时间序列波动。因此，在实际应用中，λ/4通常被视为一个"安全阈值"；超过该值的结果应谨慎对待。

将λ/4作为"安全阈值"也很大程度上源于实际考虑。如果大跳跃是真实信号，则将与黄河三角洲已知的渐进沉降特征相矛盾。事实上，我们观察到的跳跃模式多种多样：相邻时间步之间单次大幅振荡形变，以及以连续单调递减后连续单调递增为特征的峰谷序列。虽然这些跳跃在时间序列中出现频繁，但它们均缺乏物理上的渐进性和持续性。真实的土地沉降应表现为缓慢、连贯的趋势变化，而非离散的急剧跳跃或振荡性反转。因此，这类模式更可能表明各种误差的影响，而非真实形变。在缓慢变化的形变假设下，奈奎斯特定律可作为时间一致性约束的合理依据。

在此背景下，本研究提出了一种针对InSAR形变时间序列的后处理策略，旨在在缓慢变化形变的假设下减少误差。该方法可在抑制虚假异常值的同时保留真实形变信号，从而提高缓慢变形区域时间序列分析的可靠性和可解释性。

此外，已有文献充分记录，在InSAR时间序列数据的"三维（空间-时间）相位解缠"或"时间演化"分析中，存在一个假设：任何维度（包括时间维度）上相邻样本间的相位变化不应超过半个周期。例如，Hooper（2010）在其2009年的论文中指出："完整的三维相位解缠算法通常基于一个核心假设，即任何维度上相邻样本点之间的相位差通常小于半个相位周期（奈奎斯特定律）"（Hooper, 2010）。因此，本研究旨在将奈奎斯特原理应用于有效约束时间序列数据中连续获取之间的形变增量。

时间域上的奈奎斯特定律可近似表示为：

$$\Delta\phi \leq \pi \tag{2}$$

其中，$\Delta\phi$为两次观测之间的相位变化，相位与形变之间的关系为：

$$\Delta\phi = \frac{4\pi}{\lambda} d \tag{3}$$

因此可以推断，形变速率$d$必须满足：

$$d \leq \frac{\lambda}{4} \tag{4}$$

其中，$\lambda$为雷达波长，$T$为卫星重访周期。对于工作在C波段、波长λ ≈ 5.6 cm（5.546 cm）的Sentinel-1，奈奎斯特约束下在一个重访间隔（12天）内的形变增量为λ/4 ≈ 14 mm。这一物理基准构成了以下多级校正策略中增量阈值的基础。

所提出的方法始于一种多级阈值校正策略，其中基于相邻时间epoch之间的差分形变定义了一系列增量阈值（例如14 mm、28 mm、42 mm、56 mm等）。如果两次连续获取之间的形变变化超过指定阈值，则该值将被调整至与假设的物理约束相对应的适当范围内。该步骤旨在抑制通常由相位模糊或数据伪影引起的急剧跳跃，从而保留形变时间序列的物理连续性和合理性。

在InSAR形变时间序列中，某些像元可能在连续时段内表现出异常波动模式，如急剧下降后紧接着急剧上升，或反之。这些模式通常缺乏明确的物理基础，往往归因于大气扰动、随机噪声或处理伪影。即使在修正了由相位模糊引起的急剧跳跃后，识别和处理此类异常值（特征为先连续增加后连续减少，或反之）仍然至关重要，以增强形变时间序列的可靠性和物理一致性。

为进一步平滑形变时间序列并消除残余异常，本研究应用了局部均值约束方法。该技术通过参考相邻观测的平均值来调整离群值，从而稳定形变信号的时间演化并最大限度地减少数据不确定性的影响。重要的是，该方法旨在保留时间序列的固有趋势，在不引入可能掩盖有意义形变模式的过度平滑的情况下修正急剧偏差。这一后处理策略增强了形变信号的时间连续性和物理合理性，为后续时空分析提供了可靠基础。整体数据处理流程如图2所示。

图2. NMCS流程图

在校正过程中，所有操作均基于LOS（视线）数据进行，以确保波长与形变之间的对应关系。然而，所有与形变速率相关的后续展示均为转换到垂直方向后的结果。

**步骤1：形变时间序列重组**

InSAR形变时间序列被组织为二维矩阵，其中每一行对应一个独立的空间位置（即研究区域内的单个像元），每一列代表特定时间epoch获取的形变值。理想情况下，形变值表现出逐渐且连续的时间演化，反映了土地位移的缓慢特性。然而，某些时间点可能会出现突发的偏差。这些突变通常由大气扰动、相干性丢失或处理伪影引起，会显著降低形变分析的准确性。

**步骤2：两类异常的识别**

**1）相邻时间步之间急剧跳跃的识别**

异常通常表现为相邻时间步之间的急剧形变变化。因此，需要计算连续观测之间的形变差值，包括当前时间步与前一时间步的差值以及下一时间步与当前时间步的差值。如果两个相邻epoch之间的形变变化超过一定阈值，则当前时间步的值被视为潜在异常值。为更准确地校正此类急剧跳跃，定义了一组多级阈值：14 mm、28 mm、42 mm、56 mm和70 mm。校正策略如下：

- 若形变差值落在(14, 28] mm范围内，对正跳跃施加?14 mm校正，对负跳跃施加+14 mm校正。
- 若差值落在(28, 42] mm范围内，施加?28 mm或+28 mm校正。
- 若差值在(42, 56] mm范围内，校正值为?42 mm或+42 mm。
- 若形变变化超过70 mm，可相应施加更大的校正值。

该步骤确保由相位模糊或测量误差可能引起的急剧异常得到适当缓解，同时保持时间序列的一致性。

**2）时间序列中大振幅反转异常的检测**

定义一系列连续变化长度（例如3、4或5个连续增加或减少）以检测不同尺度的异常段，旨在识别波动模式，如急剧下降后紧接着急剧上升，或反之。如果某段前半部分表现为单调递减，后半部分表现为单调递增，则将其分类为"谷型异常"；反之，如果前半部分表现为单调递增，后半部分表现为单调递减，则识别为"峰型异常"。

在多阈值校正步骤的基础上，进一步根据每种模式的长度和12天间隔内预期的最大形变速率，为每种模式分配总形变阈值。例如，3步、4步和5步段分别使用13 mm、17 mm和21 mm作为阈值。这确保只有振幅显著的波动才被标记为异常。

用于识别峰谷异常的13、17和21 mm阈值（分别对应3步、4步和5步连续变化段）是根据黄河三角洲的内在形变特征确定的。这一推导遵循保守原则：对于绝大多数稳定变形的像元（约99%），在整个8年观测期间的累积沉降不超过1000 mm。鉴于Sentinel-1A的12天重访周期，总获取间隔数为233个，可以估算短时间窗口内的最大合理形变增量。据此，3步、4步和5步的阈值计算为(1000/233) × n，其中n = 3、4和5，分别约为13、17和21 mm。这些值因此代表了渐进、稳定沉降预期增量的上界。在此类窗口内超过这些限制的波动被视为物理上不合理，并标记为异常。

重要的是，这些阈值不是普遍可移植的。它们明确针对黄河三角洲的特定形变模式进行校准。当将NMCS框架应用于其他三角洲区域或不同研究区域时，用户应根据当地关于最大累积沉降、形变速率和时间采样间隔的先验知识重新校准这些参数。例如，沉降较快的区域（例如研究期间累积超过1000 mm）需要按比例更高的阈值，而形变较小的区域则受益于更严格的约束。一般指导原则是将阈值设为(D_max / N) × n，其中D_max为场地特定的最大合理累积形变，N为时间间隔数。这确保了NMCS方法在不同地质和人为环境下的适应性和物理合理性，同时保留了其缓慢变化形变的核心假设。

对于每个检测到的模式，计算下降段和上升段（或反之）的形变量。如果两段的绝对值均超过各自的振幅阈值，则确认该波动为异常。为校正异常，该段内的值被替换为该段首尾值之间的线性插值。这在不改变时间序列整体趋势的情况下消除了不真实的峰值或谷值。

此外，需要澄清NMCS方法如何区分真实季节性信号与异常跳跃。季节性形变——如冷季的冻胀融沉，或由补给周期驱动的地下水位波动——通常表现出周期性复发（如年周期）和在数月时间尺度上的渐进演化。相比之下，NMCS针对的峰谷异常特征是在短时间窗口内（3–5个连续12天间隔，即36–60天）的急剧反转，且没有预期会重复出现。

此外，季节信号展开的时间尺度远超检测窗口（通常>90天），而被标记的异常值则具有瞬时性，即缺乏时间持续性或周期规律性。振幅阈值（13、17和21毫米）进一步确保仅对急剧的大振幅反转——即远超长三角地区典型季节变动的部分——进行修正。因此，正常的季节模式得以保留，只有物理上不合理的短持续时间波动才会被修改。

步骤3：局部均值约束平滑
在对异常值进行初步修正后，时间序列中可能仍残存一些不规则的波动。为进一步增强时间平滑性并减少噪声，应用了局部均值约束平滑方法。

如公式5)、6)和7)所示，计算时间步t?1和t+1处的变形均值。若时间t处的变形值与此局部均值偏离显著——超过预设阈值——则用均值替换该值。否则，保留时间t处的原始值。

设原始时间序列为D[d1, d2, ..., dn]，其中dt表示时间t处的变形值。此外，设定阈值θ。

计算局部均值：
(5)
计算当前值dt与局部均值mt之间的偏差：
(6)
根据偏差与阈值θ之间的关系确定平滑值：
(7)

该方法有助于抑制孤立异常值，同时保留时间序列的潜在趋势。

3.5. 地面沉降时间序列的K-means聚类分析
由于高维性、时间依赖性、固有噪声以及大数据量，InSAR沉降时间序列的解读面临重大挑战。为有效分析变形演化模式、揭示地面沉降的时空异质性，并增强时间序列数据的可视化，本研究将K-means聚类算法应用于InSAR时间序列。

K-means聚类特别适用于大规模数据处理。它将具有相似变形模式的时间序列归为聚类，便于数据压缩、隐藏结构识别、异常检测以及对复杂时空系统的深入理解。此外，对于大规模数据集，直接可视化可能不切实际或过于杂乱。通过将数据组织为不同聚类，K-means使得能够使用高级可视化技术更清晰地呈现数据中嵌入的结构模式，从而提高可解释性和分析洞察力。

影响聚类结果的一个关键因素是聚类数量的选择。本研究使用肘部法则确定最优聚类数量，结果表明五类最为合适。关于K-means聚类应用于InSAR时间序列数据的详细介绍，请参阅作者的先前出版物（Wang等，2024）。

3.6. 时间序列的主成分分析
尽管平均趋势线提供了简单、直观且易于实施的显示时间序列趋势的方法，但它们主要反映整体平均水平，可能无法捕捉数据中的重要波动或结构特征。此外，MT-InSAR沉降时间序列本质上十分复杂，由实际变形信号与大气伪影和噪声交织而成。

为更好地从该复杂数据集中提取潜在模式，本研究对聚类后的时间序列应用了PCA。PCA作为一种有效的降维工具，能够揭示变形过程的内在结构和主导特征，增强对时空变异性的解释。

本研究采用PCA以深入了解每个聚类中的变形特征，并验证误差校正程序的有效性。PCA的理论基础及其在InSAR沉降时间序列中的应用已在作者的先前研究中进行了全面讨论（Wang等，2022）。此外，PCA还可以验证时间序列误差校正的合理性和有效性。

4. 结果
本研究对变形时间序列在误差校正前后的聚类分析和PCA进行了比较。在展示误差校正后结果显著改善的基础上，我们从校正后的时间序列中提取了年沉降速率。以下各节呈现与这些发现相对应的结果和分析。

4.1. 时间序列校正前后的聚类结果比较
通过K-means聚类，我们将长三角地区的沉降时间序列分为五个聚类，如图3所示。该分析揭示了近年来长三角地区地面沉降存在明显的时空变异性。图3a展示了误差校正前的聚类结果，图3b展示了误差校正后的结果。聚类1至聚类5代表沉降程度递增，聚类1相对稳定，聚类5沉降最为严重。从图中可以明显看出，时间序列校正影响了聚类结果的空间分布。校正后，稳定区域（聚类1）显著扩大，变得更加集中和连续。

图3. 校正前后的K-means结果。

通过统计各聚类的像素数量，如图4所示，可以看出误差校正后，聚类1代表的稳定区域显著增加，而聚类2代表的轻微沉降明显减少。聚类3和4代表的大规模变形也有所减少，而聚类5代表的严重沉降略有增加。

图4. 误差校正前后聚类像素数比较。

4.2. 8年平均沉降速率
2016年1月至2023年11月的平均沉降速率如图5所示。多年平均沉降速率图显示长三角地区存在严重地面沉降，幅度大且具有多个显著集中的沉降区（沉降漏斗）。为验证我们InSAR导出的速度场，我们将结果与两项近期独立研究进行了比较。我们的变形模式与Elliott等（2026）的跨区域速度场以及Xu等（2026）的长三角地区多传感器InSAR时间序列总体一致，特别是在主要沉降漏斗的空间分布方面。这种与独立数据集和处理策略的交叉验证支持了我们速度估计的稳健性。

图5. 2016年至2023年平均沉降速率。三角形位于多年平均沉降速率最快的点处。(a)为LOS变形速率图，(b)为LOS变形速率对应的标准差。(c)为垂直变形速率图，(d)为垂直变形速率对应的标准差。

长三角地区最大平均沉降速率可达308毫米/年，位于长三角地区的卤水开采区。最大沉降的具体位置由图5中的三角形标出，坐标为37°56′47.16879″，118°47′48.5263″。整个长三角地区的总平均沉降速率为16.1毫米/年，标准差为15.1毫米/年。

4.3. 2016年至2023年年度沉降速率
本研究基于校正后的沉降时间序列和MA方法生成了2016年至2023年的年度沉降速率。年度沉降速率结果（图6）揭示了长三角地区地面沉降的显著空间变化，沉降漏斗有明显位移。特别是2019年之后，位于古隆油田北部的主要沉降中心几乎停止了沉降。

图6. 统一色标的年度沉降速率。

2018年沉降速率超过100毫米/年的面积最大。虽然此后沉降仍在继续，但明显可见总体趋势逐渐减缓，沉降速率最快区域的空间位置发生了转移。

最严重的沉降发生在2017年、2018年和2019年，这由年度沉降速率图的极端值和平均沉降速率统计分析所表明。结合图6和表III的观测，可以清楚地看到，尽管存在一些波动，长三角地区地面沉降的整体趋势呈现明显的减速模式。表III中的最大值和平均值进一步证实了地面沉降的递减趋势，而标准差则表明空间变异性在逐步减小。

表III. 年度沉降速率统计。
年份 最大值（毫米/年） 平均值（毫米/年） 标准差 坐标
2016 -369 -23 3.3 37°57'01"N, 118°45'43"E
2017 -483 -23 34.6 38°00'37"N, 118°42'07"E
2018 -448 -17 37.4 37°56'28"N, 118°47'55"E
2019 -446 -18 29.8 37°56'49"N, 118°47'52"E
2020 -377 -15 28.3 37°57'34"N, 118°47'34"E
2021 -372 -12 23.6 37°57'25"N, 118°47'40"E
2022 -336 -13 22.3 37°57'28"N, 118°48'07"E
2023 -271 -15 19.1 37°57'01"N, 118°48'13"E

图7突出了沉降最快点、无人机数据采集点以及典型沉降区域的位置。

图7. 地面沉降分析关键位置索引。蓝色实心圆点为图8所示时间序列点的位置。年度极端沉降点（Year系列），典型沉降区域（T系列）。红色十字标记代表无人机测量点，在此现场拍摄了航空图像。绿色三角形标示多年沉降最快点的位置。黑色标记代表具有不同驱动因素F的典型沉降区域。

根据图6和图7，我们可以发现：

古隆油田区域严重沉降：如图6所示，2017年，位于古隆油田北部（图7中的T6）的盐滩区域（图7中的T5）在幅度和范围上都经历了显著沉降。这种严重沉降影响了古隆油田及其沿海堤坝。然而，近年来该区域沉降显著减缓，几乎未再观察到进一步沉降。

黄河以南沿海区域沉降停止：黄河以南沿海区域，围绕东营胜利机场（图7中的T8），地面沉降持续至约2021年，之后完全停止。

古隆油田西北部沉降加剧：在古隆油田西北部，靠近北孤水库（图7中的T3）的盐滩区域，沉降显著加剧，沉降漏斗的覆盖面积明显增加。到2021年和2022年，它成为该地区最突出的沉降漏斗。

黄河入海口附近沉积物抬升：除沉降漏斗区域外，黄河入海口海岸线附近（图7中的T7）观察到抬升现象。其他非沿海区域则普遍存在沉积物压实。

广饶油田附近保持稳定：广饶油田附近（图7中的T4，靠近仙河镇和孤岛镇）地面沉降多年来保持相对稳定，波动较小。

稳定农业区域：在没有卤水或石油开采活动的农业地区（如河口以西，图7中T1附近），沉降通常较小，速率主要在一倍半标准差以内，表明该区域相对稳定。

新兴沉降漏斗：图7中T2附近形成一个新的沉降漏斗。该区域此前未被识别，由于缺乏显著沉降且存在严重去相关，导致仅保留少数像素，难以及时检测。

5. 讨论
通过对沉降速率图和时间序列数据的分析，观察到长三角地区地面沉降多年来并未无限扩大。相反，沉降的时空特征发生了显著变化。这些变化有效反映了人类活动的影响，并凸显了随时间实施的灾害缓解措施的有效性。

5.1. NMCS与低通、高斯及卡尔曼滤波方法的比较
InSAR时序分析导出的原始变形时间序列通常具有显著的高频噪声。在低相干区域，这通常表现为突出的随机噪声。此外，时间序列可能表现出非标准随机高频扰动，如相位解缠错误导致的突变阶梯或异常值。虽然经典滤波技术——如时序低通滤波、高斯滤波和卡尔曼滤波/递归估计——通常能有效抑制常规随机高频噪声，但在准确识别和去除非标准随机高频扰动方面往往力不从心。

本研究选取了四个沉降程度不同的代表性点位进行对比分析。

如图8所示，可以直观地看出，低通滤波和高斯滤波方法导致时间序列过度平滑，可能会损害变形趋势。相比之下，卡尔曼滤波和非均值残差聚类（NMCS）滤波方法不仅能够保留关键的变形趋势和细节，还能有效抑制一定程度的高频噪声。特别是在前三行的对比中，NMCS在抑制非常规高频噪声的同时保持适当的平滑度方面展现出卓越的能力，在一定程度上优于卡尔曼滤波。

下载：下载高分辨率图片（875KB）
下载：下载全尺寸图片

图8. 时间序列滤波方法对比：低通、高斯、卡尔曼和NMCS。上述四个点分别位于R1、R2、R3和R4处，如图7所示。

高频抑制能力（HF Supp）衡量滤波器从时间序列中去除高频成分（通常被视为噪声）的能力。其基本原理是，一个理想的滤波器在有效消除噪声的同时，会降低信号的高频能量。因此，其噪声抑制能力可以通过计算滤波后信号高频能量的衰减程度来量化。HF Supp值越接近1，表明对干扰成分的抑制越强，从而增强信号的平滑度和可解释性。

斜率相似度（SlopeSim）评估滤波后时间序列与原始时间序列在长期趋势方面的一致性。它关注滤波器是否有效地保留了代表物理过程的长期、缓慢变形信号，同时去除噪声。SlopeSim值接近1意味着滤波过程有效地保留了潜在的变形趋势，同时减少了噪声。

因此，综合考虑这两个指标可以对滤波器的去噪和趋势保留性能进行全面评估。这种双视角评估在基于InSAR的地表变形分析中尤为有价值，因为信号保真度和干扰抑制在此均至关重要。

如图9所示，低通滤波和高斯滤波方法在HF Supp和SlopeSim两个指标方面均表现良好。然而，结合图8进行考虑，可以明显看出两种方法都对时间序列施加了过度平滑。通过对图8和图9的联合分析，观察到NMCS和卡尔曼滤波的指标值非常相似，但未引入过度平滑的迹象。值得注意的是，NMCS在抑制高频噪声方面相较于卡尔曼滤波具有比较优势，表明其在保持变形趋势完整性的同时具有更优越的去噪性能。

下载：下载高分辨率图片（1MB）
下载：下载全尺寸图片

图9. 多种滤波方法的残差与评估指标（高频抑制与斜率相似度）。

5.2. 误差校正对聚类结果的影响

根据图10中呈现的K-means聚类结果，我们观察到地面沉降模式并非始终呈线性。相反，它们随季节性变化而变化。对原始时间序列聚类（左列）的分析揭示了以下问题：（1）对原始数据应用K-means聚类后，在同一聚类的尾端观察到显著分歧；（2）季节性波动的幅度过大，变化急剧且突然——相邻区间可表现出超过200毫米的变形差异。

下载：下载高分辨率图片（2MB）
下载：下载全尺寸图片

图10. 以聚类形式展示的沉降时间序列折线图。子图a1至e1和a2至e2呈现了聚类分析的结果，包括带有趋势线、平均线和整体趋势的变形时间序列。

在对校正后的时间序列进行重新聚类后，观察到以下改进：（1）不合理的波动被有效抑制；（2）基于新的聚类结果，主成分分析（PCA）中每个聚类的方差解释比例显著提高。此外，聚类均值曲线的平滑度大大增强，聚类变得更加紧凑，尽管仍然存在一些簇内变异性。

通过分析各个聚类结果，我们观察到以下情况：

聚类1：由图10的子图(a1和a2)中的时间序列代表，该聚类表现出相对较小的沉降，伴有显著的季节性波动。聚类平均趋势线和相应的PCA结果（图11，a1和a2）均表明该区域存在长期趋势。校正后时间序列的聚类结果（图10，a2）显示波动明显减弱，无法解释的突变大幅减少。此外，如图4所示，尽管聚类1中的时间序列数量大幅增加，但PC1解释的方差比例提高了多达12%。这表明尽管数据量更大，但聚类内部的总体一致性和趋势清晰度得到了进一步增强。

下载：下载高分辨率图片（2MB）
下载：下载全尺寸图片

图11. 每个聚类的主成分分析结果，其中每个子图中的三条线分别代表PC1、PC2和PC3。

聚类2：由图10的子图(b1和b2)中的时间序列代表，该聚类在8年间显示出约100毫米的累积平均沉降量，沉降速率适中。如图11 b1所示，误差校正前的季节性波动十分显著。校正后的PCA结果（图11，b2）表明2021年6月左右沉降略微减缓。误差校正后，该聚类中的时间序列数量显著减少，如图4所示。平均曲线和PCA结果均表明时间演化过程更加单调。两种情况下，时间序列在末端均表现出轻微的减速趋势，始于2021年6月左右。此外，PC1解释的方差比例在校正后增加了8.3%，证实了簇内一致性得到了改善。

聚类3：由图10的子图(c1和c2)中的时间序列代表，该聚类在校正前未表现出严格的线性趋势。相反，变形时间序列在不同时间段显示出轻微的加速和减速阶段，并仍受到无法解释的突变影响。误差校正后，观察到更清晰的减速趋势，几乎不受突然异常的影响。聚类3的校正后PCA结果（图11，c2）表明总体呈现减速模式，尽管部分时间序列仍显示出加速沉降趋势。该聚类中的像元数量较之前有所减少，PC1解释的方差比例仅增加了1.9%。

聚类4：该聚类中的未校正变形时间序列表现出总体减速趋势，如图10 d1所示，但季节性波动严重阻碍了对转折点的识别。校正后，不同时间段的加速和减速模式变得更加明显，转换时机的界定也更为清晰，如图10 d2所示。校正后，聚类平均时间序列在变形信号的平滑度方面有了显著改善。校正后，PC1中的减速模式变得更加突出（图11，d2），拐点时机的界定也更加清晰。由于像元数量的增加和极端值的存在，PC1的解释力下降了约0.6%。

聚类5：如图10（e1和e2）所示，聚类5与聚类4类似，其时间序列表现出强烈的分歧。虽然PCA结果（图11，e1和e2）成功捕获了主要时间趋势，但聚类5的特征在于沉降量显著大于聚类4。尽管如此，两个聚类均表现出明显的减速趋势。误差校正后，聚类5的总体变形趋势基本未变。然而，变形时间演化过程变得更加清晰，有助于更准确地识别沉降过程中的转换点。PC1的改善相对较小，方差解释比例增加了1.9%，表明在这一方面增强有限。

从图11左右列的对比可以看出，误差校正显著抑制了季节性异常对时间序列的影响。基于平均变形曲线，聚类1、2和3均表现出接近线性的趋势。从图11 d1和d2所示的平均变形曲线和PCA结果可以看出，聚类4中大多数时间序列显示出减速趋势，尽管部分地区仍存在加速沉降。聚类5的平均曲线揭示了明显的减速趋势，但图11，e1和e2所示的PCA结果表明出先稳定后呈线性趋势的模式，反映了聚类5所代表区域沉降更为复杂的特征。误差校正显著提高了InSAR时间序列的总体一致性和趋势清晰度，特别是增强了数据中主导变形趋势的表达。

我们还观察到，位置越接近沉降中心，PC1解释的方差比例越低。从数学角度来看，PCA主要用于提取数据中的主导趋势。然而，当样本数据表现出复杂多变的模式时（如非线性变化、突然变化或多个沉降机制的叠加），PCA难以识别统一的主要趋势。在这类"无序"的时间序列数据中，PC1的解释力较弱，导致方差解释比例较低。

在沉降中心，时间序列可能表现出非线性沉降、多过程叠加以及季节性和长期趋势混合等模式。此外，人类活动可能导致沉降过程更加复杂和动态。这些自然和人为因素相互耦合，放大了沉降过程的不稳定性，使得PC1难以提取统一的主导变化模式，导致方差解释比例较低。尽管沉降中心变形更为显著，但多样化的趋势和复杂非统一的时间序列波动导致了更大的异质性。

在沉降边缘，变形幅度较低，干扰较小。干扰的主要来源可能是区域大气残余或仪器误差，经校正后时间序列更加一致。因此，这些时间序列的变化模式更容易被单个主成分（如线性沉降）所解释，导致PC1解释的方差比例较高。

5.3. 沉降趋势的空间收缩

图12展示了2017年至2023年间沉降速率超过100 mm/yr的区域。从该图中可以得出两个关键观察结果：沉降中心的明显偏移和严重沉降空间范围的显著减少。这些趋势与图10中呈现的时间序列中观察到的时间变形模式一致。特别是古丰油田西北部地区的沉降量级和空间覆盖范围均显著减小。这种减少与禁止卤水开采的监管措施的实施密切相关，突显了政策干预在缓解人为诱发沉降方面的有效性。

下载：下载高分辨率图片（3MB）
下载：下载全尺寸图片

图12. 沉降速率超过100 mm/yr的沉降面积年度变化。

下载：下载高分辨率图片（217KB）
下载：下载全尺寸图片

图13. 沉降速率超过100 mm/yr的像元数量。

针对卤水开采的监管干预并非单一政令，而是一个跨多个行政层级逐步加强的过程。最初阶段始于2019年3月，河口的区人民政府发布了《推动溴产业规范化发展的实施方案》（东港公信发[2019]14号），针对非法卤水开采和未经批准的水资源使用行为。2019年8月，台风"利奇马"为该区域带来强降雨和海水上溯，与长期地面沉降叠加，导致古丰油田北部地区大面积洪涝灾害。这一事件显著提高了公众对沉降相关灾害的认识，并引起了政府机构和油田管理层的广泛关注。随后的调查确定了地下卤水开采为主要致因。2021年11月，全省范围内开展非法取水专项整治行动，明确将天然卤水列为重点调查和取缔对象。

具有决定性意义的法律里程碑出现在2023年5月，当时《山东省黄河三角洲生态保护条例》正式颁布，并于2023年8月1日起施行，为禁采卤水提供了法定依据。随后，2023年11月，东营市政府启动了"地下卤水抽取整改工作方案"，推动了地下卤水井的物理封堵。监管力度的加强与沉降减缓之间在时间上的巧合尤为引人注目。如图12所示，最显著的沉降面积缩减（超过100 mm/yr）出现在2019年之后，与最初的地方整改工作相吻合。2021至2022年，随着全省范围和各区域的整治行动深入推进，沉降降幅进一步扩大。到2023年，随着省级法规的出台和强制封堵措施全面落地，沉降面积已缩减至2017年以来的最小范围。这种时间上的对应关系支持了因果关系的解释：卤水抽取监管力度的逐步收紧直接导致了海岸沉降速度的减缓。

此外，即使在缺乏直接政策干预的情况下，土地沉降预计长期内也会自然缓解。地下卤水在短尺度上是一种不可再生的矿产资源，长期开采导致地下水位和卤水盐度逐渐下降，从而显著增加了开采成本（如电费上涨），迫使许多盐场放弃传统的卤水开采方式。作者在2023年11月进行的实地调查证实了禁采卤水政策的执行。目前，卤水的主要来源已转变为海水，通过连接海洋的人工运河引入内陆。

这一观察强调了两个关键认识：（1）土地沉降的隐蔽性：沉降往往在缺乏专门监测的情况下难以被发现，导致沉降发生与缓解措施实施之间存在滞后。因此，应对措施通常是被动性的，往往在实质性损害发生后才启动。（2）土地沉降的破坏性影响：沉降对生产安全构成严重威胁，迫使决策者放弃重大经济收益以优先进行沉降防控。这是因为土地沉降造成的损害通常不可逆转，无法通过短期经济收益来弥补。

5.4 时间变化及驱动因素

本研究采用聚类分析与主成分分析（PCA）相结合的方法，对InSAR获取的土地沉降时间序列进行了分析，揭示了显著的转折以及沉降整体减缓的趋势。结果总结如下：

聚类1：该区域以极轻微的线性沉降为主，主要受季节波动影响。第一主成分（PC1）与该聚类的平均趋势线高度一致，表明季节物候变化对InSAR信号的影响强于实际地面变形。该区域的人类活动历史悠久、系统性强且强度稳定。结合土地利用分类与卫星影像，该区域的主导活动似乎是农业。

聚类2：位于地下卤水抽取影响范围内，该聚类表现出低幅度但持续沉降的长期趋势。在误差校正前，观测到显著的季节性变化。该模式反映了人为因素（如卤水抽取）与自然因素（如水文或气候）的综合作用。尽管该区域历史上一直经历长期沉降，但PC1显示近期沉降速度有所减缓，可能与地下水抽取管控力度的加强有关。作为紧邻沉降中心的过渡区域，其对政策的响应存在滞后且不够明显。其空间模式与典型地下水抽取引起的沉降相吻合，即变形幅度从中心向外呈径向递减。

聚类3：该聚类显示出持续且日益显著的沉降趋势，是受沉降影响最严重的区域之一（位列第三）。PC1揭示了稳定连续的沉降模式，但2022年6月之后，局部地区沉降速度明显加快，表明某些区域的沉降正在加剧。从地理位置上看，该聚类位于卤水抽取与农业活动主导区域之间的过渡地带。过去卤水抽取的影响依然显著，若不加干预，该区域可能演变为新的沉降漏斗。

聚类4：该聚类的平均时间序列趋势与PCA结果一致。PC1显示从稳定沉降向趋近平衡状态转变。PC2反映了初期相对稳定的阶段，随后沉降开始显现，但其作用相对次要。从空间上看，该聚类毗邻沉降最为严重的聚类5，是其空间延伸部分，反映了从中心向外沉降强度递减的经典梯度。然而，自2021年年中左右起，该区域大部分已进入平台期，沉降活动有所减少，表明主要驱动因素得到了有效控制。尽管如此，某些局部区域仍显示出变形加快的迹象，表明沉降可能存在重新分布的趋势。

聚类5：该聚类对应沉降核心区域，表现出最复杂的时空演化特征。初期，区域部分保持稳定，随后整个区域过渡到线性沉降阶段。平均趋势线揭示了两个主要的分化阶段：第一阶段约在2019年8月，沉降出现明显减缓；第二阶段约在2020年8月，显现出分化行为——部分地区趋于稳定，而另一些地区沉降恢复或加剧。后者的沉降速率回升至2019年8月之前的水平。虽然PCA捕捉了主导性的时间特征，但无法完全反映沉降减缓显著区域的异质性行为。总体模式表明复杂的非同步时空动态，这可能是由于不同区域卤水开采禁令的执行力度不均所致。通过对比可以看出，聚类4在2021年左右进入平台期，表明浅层沉降在卤水抽取停止后约两年内趋于稳定。相比之下，聚类5的持续沉降表明深层含水层的恢复需要超过两年的更长周期。

5.5 空间重新分布的灾害风险

尽管沿海沉降减缓是政策干预的一项显著成效，但这并不意味着灾害风险的消除——相反，风险在空间上发生了重新分布。沉降漏斗向内陆转移对防洪基础设施、土壤条件和风暴潮脆弱性带来了新的威胁，亟需重新关注。

首先，防洪基础设施面临双重压力。尽管孤东油田地区沿海大堤因本地沉降减缓而差异沉降有所降低，但在古北水库（图7中的T3）附近及其他内陆地区新加剧的沉降漏斗，对排水泵站、河堤和城市防洪系统构成了新的威胁。与基于海平面上升设计的沿海防御设施不同，内陆基础设施在建设时往往较少考虑渐进性地面沉降的影响，因此对沉降导致的容量损失可能更为脆弱。

其次，土壤盐碱化可能因沉降向内陆转移而加剧。地表高程降低导致自然排水梯度减小，阻碍咸水地下水和地表水的排出。在一个已具有广泛盐碱土壤和高地下水位的地区，即使沉降速率有所降低，持续的沉降仍会延长渍涝时间并阻碍盐分淋洗，这在排水基础设施有限的农业区域尤为突出。这与以下观察结果一致：即使在没有主动卤水抽取的情况下，厚层沉积物长期固结仍会维持一种背景沉降，使不利的排水条件持续存在。

第三，风暴潮脆弱性并不局限于沿海地区。内陆地形的降低，加上该地区极低的地势（黄河三角洲大部分地区海拔低于4米），意味着风暴潮可通过河道和人工运河向内陆更远处推进。为卤水生产输送海水而修建的运河，可能无意中成为风暴潮传播的通道。随着相对海平面持续上升，新形成沉降漏斗的内陆地区将面临日益增大的复合洪水风险——即风暴潮、河流洪水与降雨洪水交汇叠加。

这些空间转移的风险表明，某一区域的沉降缓解并不等同于整体风险的降低。相反，适应性管理策略必须随着沉降空间模式的变化而调整，将持续监测（如基于InSAR对新出现沉降漏斗的监测）、新受影响内陆区域的基础设施升级以及排水系统加固（同时应对盐碱化和洪水风险）相结合。黄河三角洲的经验表明，三角洲地区的灾害风险降低不仅需要停止沉降的驱动因素，还需要在时空上管理其重新分布的后果。

5.6 NMCS方法的局限性

NMCS校正提高了时间序列的平滑度和信噪比（SNR），同时不会过度修改原始变形信号，适用于地下水抽取引起的沉降等长期缓慢过程。然而，其应用存在以下局限。

首先，该方法假设真实变形增量不超过奈奎斯特（Nyquist）限。这一假设在黄河三角洲地区成立：2017年最大沉降速率约500 mm/yr，对应每12天约16 mm，仅略高于14 mm的阈值。然而，在快速或瞬态变形区域（如地震、滑坡、火山抬升、风暴潮、同震位移和滑坡加速），真实的快速信号可能被错误地校正或被过滤掉。因此，NMCS应严格限制在已知变形缓慢且连续的情况中使用，并且阈值应根据研究区域的具体特征进行调整，以避免过度平滑。

其次，仅使用短基线干涉图（12–120天，约950幅干涉图）适用于低相干的黄河三角洲地区，可最大化相干性，且冗余网络仍可对长期变形进行约束。然而，它对长波长相位偏差（轨道斜坡、残余大气延迟、DEM误差和传播的解缠误差）的约束较弱，这些误差可能污染速度场，尤其是在沿海低相干区域和快速沉降漏斗附近。

6. 结论

本研究提出了一种基于奈奎斯特的多步校正与平滑方法（NMCS），用于解决InSAR时间序列数据中常见的异常波动。NMCS方法被证明是高度有效的，成功消除了显著的异常波动，大幅提高了变形信号的信噪比。

在校正之后，本研究将精化数据与K-means聚类及PCA相结合，分析了黄河三角洲多年及年度的土地沉降速率以及变形时间序列。结果揭示了过去八年中沉降显著的时空变化特征，并为背后的驱动机制提供了新的见解。观测到沉降漏斗明显的内陆迁移趋势：沿海地区的变形显著减缓甚至完全停止，而内陆地区仍持续出现沉降，某些情况下沉降甚至加剧。此外，还检测到沉降时间模式发生了显著变化。

沉降时空动态的这些变化与禁采卤水政策密切相关。实地调查证实，禁采卤水政策的出台不仅是对已观测到的土地沉降的响应，也旨在降低沿海灾害的上升风险和减轻采矿运营日益增加的经济负担。沉降模式的演化反映了人类干预的深远影响。这凸显了人类活动的双重作用：它们既是沉降的主要驱动力，也可能通过监管控制成为解决途径。这些发现强调了可持续土地利用政策和针对性缓解策略的关键重要性。

然而，尽管人为土地沉降在黄河三角洲得到了有效缓解，其残余影响可能仍然存在：1）历史沉降造成的沿海大堤结构性损害可能潜藏于地表之下，且常规遥感方法难以探测。这些潜在隐患对沿海防御基础设施构成长期风险。2）地表高程降低（过去沉降的直接结果）增加了洪水风险，并使盐碱土壤的排水更加困难，可能加速土壤盐碱化及进一步的环境退化。

作者贡献声明：

宋俊：写作——审阅与编辑、验证、监督、资源、概念化。
王国阳：写作——初稿撰写、监督、软件、方法论、调查、正式分析、数据整理、概念化。
李振红：写作——审阅与编辑、验证、监督、概念化。
李广学：写作——审阅与编辑、验证、监督、资源、概念化。
徐志强：监督、概念化。
邢路扬：资源。

董丁：监督、方法论、概念化

利益冲突声明：作者未报告任何潜在利益冲突。

数据可用性声明：所有数据均可通过通讯作者申请获取。

资助信息：本工作受国家自然科学基金（42506149）、山东省博士后科学基金（SDCX-ZG-202400169）；青岛市博士后基金（QDBSH20240102116）；辽宁省科技计划（2024JH2/102400061）；大连市科技创新基金（2024JJ11PT007）；大连市科技创新人才科技计划（2022RJ06）；辽宁省教育厅科研平台建设专项（LJ232410158056）；大连海洋大学基础科学研究基金（2024JBPTZ001）；辽宁省数据中心项目（2025JH27/10100005）以及国家重点研发计划"中国-尼日利亚河三角洲联合实验室"（项目编号2024YFE0116400）资助。

王广阳 | 宋军 | 徐志强 | 李广学 | 董丁 | 邢路阳 | 李振红