图1 升轨142轨道数据空间网格化示意图
Fig. 1 Spatial meshing of ascending track 142 data
北京大学学报(自然科学版) 第59卷 第6期 2023年11月
Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 59, No. 6 (Nov. 2023)
doi: 10.13209/j.0479-8023.2023.069
国家自然科学基金(41774008)资助
收稿日期: 2022‒11‒17;
修回日期: 2022‒12‒10
摘要 基于 2015—2019 共 4 年的合成孔径雷达(SAR)卫星影像, 提出一种广域合成孔径雷达干涉测量(In-SAR)时序分析技术, 对华北平原的地表形变进行高精度连续监测。首先对 SAR 影像进行干涉处理, 得到干涉图。在此基础上, 使用经过并行化改进的永久散射体技术斯坦福改进(StaMPS)方法, 提取干涉图中所有永久散射体(PS)像元, 获取研究区域全分辨率的时序形变信息。之后, 使用大气模型校正法与共景叠加法相结合的联合大气校正方法, 估计并去除形变信息中的大气噪声。经过上述处理流程后, 成功地获取华北平原地表大空间尺度、长时间跨度、全空间分辨率和高精度的形变信息, 进而监测到平原内部由长期地下水开采导致的高达 100mm/a 的大范围强烈沉降信号。相较于既有算法, 并行化 StaMPS 方法通过在多个计算节点之间的实时分配, 节约至少 60%的计算时间, 联合大气校正方法则可以去除约 74.3%的大气噪声项, 有效性显著高于两种校正方法的单独使用效果。广域 InSAR 时序分析技术可以有效地实现对大范围地表形变的高精度连续监测。
关键词 广域 InSAR 时序分析技术; 并行化 StaMPS 方法; 联合大气校正方法; 地表形变; 华北平原
华北平原是我国政治与经济的核心区域, 也是重要的粮食和棉花主产区之一, 肩负着重要的生产任务, 快速发展的工农业极大地提高了用水需求。由于降水较少, 且季节波动显著, 地下水成为华北平原农田灌溉、工业生产和生活用水的主要来源, 20 世纪 70 年代以来, 地下水开采量快速增加[1]。地下水的长期过量开采导致华北平原出现严重的大面积地面沉降现象, 产生的数个地下水降落漏斗已影响区域的生态环境, 制约了区域的可持续发展[2]。基于地下水资源在华北平原的重要战略地位, 对华北平原的沉降监测一直备受关注和重视。
借助近年来快速发展的卫星大地测量技术, 如雷达干涉测量技术(interferometric synthetic aperture radar, InSAR)和全球导航卫星系统技术(global na-vigation satellite system, GNSS), 对由抽水导致的大范围沉降信号的监测成为可能, 已经成功地测绘北京[3]、天津[4]和沧州[5]等地的沉降形变场。然而, 这些研究皆受限于测量技术, 因而着眼于局部的沉降信号, 无法形成对覆盖整个华北平原的广域沉降信号特征的整体性认识, 进而制约了地下水保护政策的规划和制定。
差分雷达干涉测量技术(differential InSAR, D-InSAR)是一种通过测量两景 SAR 影像之间的相干相位来获取观测时段内高精度地表形变的大地测量技术[6−7]。传统的 D-InSAR 技术易受时空失相干的影响, 降低相干相位的有效性[8−9], 使用场景受到极大的限制。针对相位失相干问题, 21 世纪初发展出PS-InSAR (persistent scatterer InSAR)技术[10−12]。永久散射体(persistent scatterer, PS)是一类具有较强信号散射能力及相位稳定性的地面标记物, PS-InSAR技术便是通过识别干涉图中所有被永久散射体主导的像元来获得具有高相位稳定性和高相位相干性的点目标干涉图, 进而获取有效的地表形变信息。其中, Hooper 等[13−14]发展的 StaMPS(Stanford method for PS)方法通过在分析振幅特征的基础上增加相位特征分析, 冲破了传统 PS-InSAR 方法在非城镇区域无法提取足量 PS 像元的局限, 极大地扩展了 PS-InSAR 方法的应用范围。
StaMPS 方法在设计之初主要是为了处理方圆数百平方公里范围内的局域形变信号, 并不适宜对广域形变信息进行分析和处理。虽然 Hooper 等[15]进行了一些算法上的相应改进, 但仍然无法满足对广阔的华北平原进行全分辨率广域 InSAR 时序分析的需求。此外, 由于华北平原夏季降雨较多, 大气波动强烈, 相干相位中由对流层效应导致的大气噪声会极大地影响地表形变信号的测量精度[16]。选择一种合适的大气噪声校正方法对大气噪声进行估计和去除, 是广域地表形变监测中的另一项挑战。
本文提出广域 InSAR 时序分析技术的思想, 在已有算法的基础上, 通过使用并行化 StaMPS 方法和联合大气校正方法, 实现对广域地表形变的高精度连续监测。使用欧洲航天局的 Sentinel-1A/1B 卫星从 2015 年 7 月到 2019 年 6 月采集的覆盖华北平原的 92 景升轨 142 轨道 SAR 影像数据, 开展广域InSAR 时序分析, 进而获取华北平原地表形变信息, 验证广域 InSAR 时序分析技术的有效性和计算效率的优越性。
并行化 StaMPS 方法分为两个部分。首先对所使用的 SAR 影像覆盖区域进行空间网格划分, 将整个研究区域划分成若干大小适当的相邻且存在重叠区域的片(patch), 以便降低单次计算任务对计算资源的需求。然后使用任务提交系统, 将不同片的计算任务在单个多核工作站或计算机集群中实时分配到不同的计算节点, 用计算空间换取计算时间, 通过并行化方式降低总的计算时间。
空间网格化划分的原则是确保在每个片内识别出的 PS 像元数量不能过多, 以便降低计算需求的同时, 避免划分出数量过多的片而浪费计算资源。考虑到 SAR 卫星的数据采集模式, 我们将数据沿着卫星飞行方向和卫星观测方向进行等间隔的划分, 每个片与邻近片之间需要保留重叠区, 既满足 PS像元拾取算法的空间平均要求, 也避免后期进行片的拼接时边缘处出现异常的数据跳跃现象。原则上以近似正方形为标准进行片的划分, 可以降低重叠区的总面积, 提高计算效率, 节约计算时间。此外, 如果某些区域识别出的 PS 像元分布密集, 数量远大于所设定的每片 PS 像元数量上限, 这时不需要修改整体的网格化划分方案, 只需在原有片的基础上进行次级片的划分, 这样就可以避免因划分方案调整而引入的额外计算资源消耗。
研究区域由 6 帧影像覆盖, 总共识别出约 3600万个 PS 像元, 升轨 142 轨道的空间网格化划分结果如图 1 所示。本研究采用 30×30 的一级网格划分和部分的次级网格划分, 最终生成 958 个片, 确保每个片中存储不超过 50 万 PS 像元。
如图 2 所示, StaMPS 方法对干涉数据的处理流程大致分为以下 8 个步骤: 1)PS 像元提取; 2)噪声估计; 3)PS 像元筛选; 4)冗余点去除; 5)相位校正; 6)相位解缠; 7)空间相关噪声估计; 8)大气噪声去除。前 5 个步骤只需要利用局域空间信息, 后 3 个步骤则需要用到全域空间信息。第 1 步和第 3 步 PS像元的提取和筛选是 StaMPS 方法的核心步骤, 同时也是对计算资源需求量最大的步骤。因此, 在数据处理过程中, 我们通过任务提交系统, 将不同片的计算任务在单个多核工作站或计算机集群实时分配到不同的计算节点来执行, 并行化地执行前 5 个步骤的计算任务; 确认所有片前 5 个步骤的计算任务均已完成后, 对所有片进行拼接操作, 然后对拼接后的结果进行第 6～8 步的处理, 即可得到最终的时序结果。
图1 升轨142轨道数据空间网格化示意图
Fig. 1 Spatial meshing of ascending track 142 data
图2 并行化StaMPS方法流程
Fig. 2 Flow chart of the parallelized StaMPS method
并行化 StaMPS 方法充分利用 PS 像元提取和筛选过程中依赖局域空间信息的特点, 经空间网格化处理后, 通过多个计算节点间的并行处理实现全分辨率广域 InSAR 时序数据的快速计算。同时, 在片的拼接步骤后, 重新获得大尺度的整体空间信息, 在随后的空间相关噪声估计步骤中, 可以进行大空间尺度的长波误差校正, 通过使用较少的校正参数来避免引入额外的误差。此外, 并行化 StaMPS 方法也可用于局域 InSAR 时序分析中, 以便节省计算时间。然而, 计算效率提升的幅度受空间网格化的方式、所使用的计算节点个数以及保留重叠区的大小等参数影响, 保守地估计, 可以节省约 60%的计算时间。
本研究中使用 shell 脚本, 将并行化 StaMPS 方法分别适配 PBS 集群任务管理系统和 LSF 集群任务管理系统, 并先后分别使用中国地震局地质研究所计算机集群和中国神威·太湖之光超算集群来完成华北平原全分辨率广域 InSAR 时序分析的计算任务, 成功地获得以前无法获得的华北平原全分辨率广域 InSAR 时序信号。对计算机集群来说, 各个计算节点之间的内存和算力是相互独立的, 使用更多的计算节点可以更大程度地节省计算时间, 提高计算效率。考虑到片的拼接操作的有效性, 重叠区的宽度选取至关重要。本研究中各片之间的重叠区宽度取 100 个像元, 故总的数据计算量增加 14.51%。此外, 并行处理步骤中所使用的计算机集群涉及一系列复杂的协同工作, 存储节点负责数据存储, 计算节点负责完成计算任务。在这一过程中, 存储节点与各计算节点之间数据传输的时间消耗是无法避免的。然而, 计算机集群用于节点通信的千兆乃至万兆以太网可以显著地降低所需时间。本研究中共使用 10 个计算节点, 最终节省 80%～85%的计算时间, 效果提升显著。
此外, 本文提出的并行化 StaMPS 方法亦可用于单个多核工作站, 通过将不同片的任务分配到指定的核上进行计算, 提高计算效率。如图 3 所示, 以英特尔至强 E5-2620V4 的 16 核处理器为例, 增加最大同时执行片的数量将大大减少总体计算时间, 但会增大内存读写负担, 降低计算效率的提升。因此, 应根据现有计算资源, 选取适当的同时执行任务数。
图3 计算时间对比
Fig. 3 Comparison of computation time
StaMPS 方法中对大气噪声的估计和去除一直是 InSAR 技术中的主要挑战之一, 目前仍然没有完美无缺的解决方案。本研究尝试使用联合大气校正方法, 即结合大气模型校正法与共景叠加法, 对形变信息中的大气噪声进行最大程度的估计和去除。
连续 GPS 观测台站高时间采样率和高测量精度的时间序列较好地排除了大气噪声等因素, 因此一般将连续 GPS 的观测值视为参照真值(ground truth)。我们使用 20 个 InSAR 数据覆盖范围内的连续 GPS 观测台站 2015—2019 年时间序列观测数据, 视其为地面真实形变值。以连续 GPS 观测台站时间序列为基准, 将 GPS 的三维时间序列数据投影到InSAR 所在的卫星视向线(line of sight, LOS)方向, 通过对比所有 GPS 时间序列与大气校正前后 InSAR时间序列的残差的平均值, 对大气噪声的校正效果进行定量评估。
大气模型校正法是通过计算雷达波传播路径上折射率的变化来计算大气噪声[17]。通常情况下, 空气的折射率 n 相对于真空折射率(n=1)的变化量级特别小, 为此我们定义一个特征折射率 N=(n−1)×106 来凸显折射率的相对变化。由物理实验可知, 对流层特征折射率由干、湿两部分的贡献组成[18]:
(1)
式中, P 为总气压(Pa), T 为温度(K), e 为水蒸气分气压(Pa), K1, K2 和 K3 分别为经验常数, 可由实验测得。这样, 只要通过对雷达波的传播路径进行积分, 即可得到对流层折射率对雷达波相位延迟的贡献 项[19]:
(2)
其中, λ 是雷达波的波长; θ 是雷达波与地面垂向方向之间的夹角, 即入射角; h1 是积分的高度下边界, 与地表高程相关; htop 是积分的高度上边界, 与所在纬度相关。
本研究中, 通过 Matlab 代码包 Toolbox for Re-ducing Atmospheric InSAR Noise (TRAIN)[20], 使用欧洲中期天气预报中心提供的最新一代大气模型ERA5[21]来估计大气噪声。虽然 ERA5 已经是目前时空分辨率最高的大气模型, 但其 31km 的空间分辨率远小于 InSAR 观测结果的几十米空间分辨率。同时, ERA5 长达 1 小时的时间分辨率也难以对快速变化的对流层物理参数进行精确的拟合。为了与InSAR 观测结果的高时空分辨率进行匹配, 计算过程中需要对 ERA5 模型进行时空线性插值处理。这就使得大气模型校正法只能反映大气信号(尤其是降雨较多季节的显著大气信号)的整体特征, 对局部细节特征的拟合效果较差, 降低了大气校正的有效性。与研究区内地面连续 GPS 台站时间序列观测的对比表明, 大气模型校正法将残差的平均值从37.94mm 降至 17.83mm, 去除了约 53.0%的大气噪声项, 故单独使用大气模型校正法对 InSAR 观测结果进行大气校正的有效性不高。
共景叠加法(common scene stacking, CSS)假定地表形变信号为随时间变化的线性信号或弱非线性信号, 大气噪声则是与时间无关的随机信号。通过对具有共同景的前后多个干涉对进行叠加, 就可以抵消形变信号而仅剩大气噪声[22]。叠加的干涉对越多, 对该共景大气噪声项的估计越准确。我们在时间轴上等间隔地采集 2N+1 个 SAR影像, 对于其中的第 i–1、第 i 和第 i+1 景, 它们之间的干涉对形成的相位差如下:
(3)
(4)
式中,
代表第 i 景和第 i–1 景形成干涉对的相位差, αi是第 i 景的大气噪声项,Δτ 是形变项,∈是误差项。在形变信号为随时间变化的线性信号或弱非线性信号假设前提下, 可以认为等时间间隔采样的两个干涉对的形变信号部分是相等的, 即
。忽略误差项, 得到第 i景的大气噪声项:
(5)
在此基础上, 如果继续以第 i 景为共景, 叠加更多的干涉对, 那么可以得到
(6)
在大气噪声是时间无关的随机信号的前提下, 我们可以认为
(7)
由此, 即可得到第 i 景的大气噪声项:
(8)
本研究中, 通过 Matlab 代码包 Matlab Code for Atmospheric Noise Depression through Iterative Sta-cking (MCANDIS)[23], 并在其基础上进行针对广域InSAR 时序分析的并行化改进, 对大气噪声进行估计。共景叠加法中大气噪声的随机信号假设对大气信号波动平缓时期(如降雨较少的季节)的校正效果较好, 对于大气信号剧烈波动时期(如降雨较多的7—9 月), 因季节性趋势相位和非线性相位的存在, 反而会越校正效果越差, 导致出现负校正现象。通过与研究区内地面连续 GPS 台站的时间序列观测值对比, 共景叠加法将残差的平均值从 37.94mm 降至 16.66mm, 去除了约 56.1%的大气噪声项, 故单独使用共景叠加法对 InSAR 观测结果进行大气校正的有效性不高。
单独使用大气模型校正法或共景叠加法对大气噪声进行估计, 各自存在优缺点。大气模型校正法对大气噪声的整体特征(尤其是降雨较多季节的信号)估计效果较好, 但空间细节较差; 共景叠加法对大气噪声的细节特征(尤其是降雨较少的季节)估计效果较好, 但降雨较多的季节会出现负校正。从校正结果来看, 两种校正方法可以相互补充。因此, 本研究使用大气模型校正法与共景叠加法相结合的联合大气校正法对大气噪声进行估计。
联合大气校正方法的步骤如下: 先用大气模型校正法对大气噪声的整体特征进行估计; 在去除大气噪声的整体部分后, 再用共景叠加法对残余大气噪声的细节特征进行估计; 将两个部分的总和作为对大气噪声总的估计值。同样, 通过与研究区内地面连续 GPS 台站的时间序列观测值对比, 联合大气校正法将残差的平均值从 37.94mm 降至 9.77mm, 去除了约 74.3%的大气噪声项, 显著高于大气模型校正法的 53.0%和共景叠加法的 56.1%。
联合大气校正方法可以同时利用两种校正方法的优势, 避免各自的缺陷, 同时较好地实现降雨较多季节和降雨较少季节的大气校正, 是一种较为全面且效果显著的大气校正方法。当然, 联合校正方法也有缺点, 如处理流程繁琐、计算时间较长、继承自共景叠加法的线性形变信号假设和舍弃前后无足够叠加对的影像信息等。我们使用经过并行化改进和算法优化的 MCANDIS 程序来解决联合大气校正方法中对计算资源的需求, 优化处理流程, 降低计算时间。因此, 联合大气校正方法不失为一种有效的大气噪声校正方法。
受长期地下水开采活动影响, 华北地区存在显著的大范围地面沉降现象[24]。为了验证广域InSAR时序分析技术反演地表形变的精度和有效性, 我们选取华北平原为研究区域(图 4), 使用欧洲航天局Sentinel-1A/1B 采集的 2015 年 7 月 30 日到 2019 年 6月 3 日升轨 142 轨道共 92 景 6 帧 SAR 影像数据。使用 GAMMA 软件[25]生成干涉图, 使用并行化 Sta-MPS 方法生成时间序列, 并使用 SRTM[26]提供的 30m 分辨率的 DEM 作为参考高程数据, 估计并移除地形相位效应。
黑色曲线内为华北平原, 青色虚线框内为升轨 142 轨道(A142) SAR 影像数据覆盖区域, 下同; 蓝色线示意主要河流
图4 研究区地理位置
Fig. 4 Geographical location of the study area
图5 干涉对的时空基线
Fig. 5 Spatio-temporal baseline of interference pairs
图6 联合大气校正前的解缠相位时间序列
Fig. 6 Unwrapped phase time series before joint atmospheric correction
图7 联合大气校正后的解缠相位时间序列
Fig. 7 Unwrapped phase time series after joint atmospheric correction
图8 华北平原平均形变速率场
Fig. 8 Mean deformation rate field in the North China Plain
在时空相干性最优的条件下挑选一个共同的参考景, 并以共主景方式形成干涉对, 最终选定的主景是 2017 年 12 月 10 日, 干涉对的时空基线见图 5。通过并行化 StaMPS 方法, 在覆盖区内识别出 3600万个 PS 像元。随后, 使用联合大气校正方法对缠绕结果中的大气噪声项进行压制, 其中前 3 景和后3 景因叠加对数不足而被舍弃。图 6 和 7 分别是剩余 86 景联合大气校正之前和之后的时间序列, 通过对比, 可见联合大气校正方法的校正效果良好, 尤其是 7—9 月降雨季节中强烈的大气噪声项被很好地估计并去除。
经过广域 InSAR 时序分析技术的一系列处理, 我们获得华北平原的 InSAR 时间序列结果。通过简单的最小二乘回归计算, 并以连续 GPS 速度场作为速度参照值进行标定, 获得华北平原 LOS 方向的平均速度场(图 8(a))。
华北平原的构造运动以水平运动为主, 且运动速率较小[27]。考虑到近期华北平原地震活动性较弱, 近 20 年内未发生超过 Mw 5.0 的地震(数据来源于美国地震学研究联合会(Incorporated Research In-stitutions for Seismology, IRIS), 可知华北平原现阶段形变以垂向形变为主导。如图 8(b)所示, 我们使用连续 GPS 台站的水平方向平均速度数据(近东南向, 3～4mm/a), 将其空间二维线性插值视为华北平原的水平方向平均速度场, 结合卫星轨道参数数据, 通过视向线方向的平均速度场, 求解华北平原垂直方向的平均速度场, 结果如图 8(c)所示。
由此, 我们成功地解算出华北平原 2015—2019年的平均垂向速度场。可以看出, 华北平原的地表形变场主要特征为与人类活动直接相关的垂向形变信号, 而非与构造活动相关的水平形变信号。华北平原内部存在大量显著且连续的地面沉降信号, 最大沉降速率高达 100 mm/a, 集中在人口密集的城市周围及大量种植农作物的高用水量区域[28], 其形变成因与地下水开采和南水北调工程回灌密切相 关[29]。同时, 华北平原南部的鲁西平原地区也有部分与城镇群紧密伴随的沉降区域。此外, 华北平原西北方向的太行山山区和东南方向的泰山山区展现广泛分布的 2～3mm/a 的抬升信号。其中, 山区的抬升信号可能由平原内部长期大范围抽取地下水造成的质量卸载和基岩回弹导致[30], 其物理机制值得进一步探究。
本文提出一种并行化 StaMPS 方法和联合大气校正方法组成的广域 InSAR 时序分析技术, 使用华北平原 2015—2019 年的 SAR 卫星影像数据, 成功地获取华北平原内地表大空间尺度、长时间跨度、全空间分辨率和高精度的形变信息, 进而监测到平原内部由长期地下水开采导致的高达 100mm/a 的强烈的大范围沉降信号。
相对于既有算法, 并行化 StaMPS 方法可以大幅度节约计算时间, 实现对广域形变场的有效监测; 联合大气校正方法可以有效地去除大气噪声, 有效性显著高于两种校正方法单独使用的效果, 且通过并行化改进提高了计算效率。广域InSAR时序分析技术可适用于大范围地表形变的高精度连续监测。
致谢 研究工作得到北京大学周仕勇教授的帮助, 中国地震局地震研究所赵斌研究员提供连续GPS观测数据, 对此表示由衷的感谢。
参考文献
[1] 张兆吉, 费宇红, 陈宗宇, 等. 华北平原地下水可持续利用调查评价. 北京: 地质出版社, 2009
[2] 杨会峰, 曹文庚, 支传顺, 等. 近 40 年来华北平原地下水位演变研究及其超采治理建议. 中国地质, 2021, 48: 1142–1155
[3] 雷坤超, 罗勇, 陈蓓蓓, 等. 北京平原区地面沉降分布特征及影响因素. 中国地质, 2016, 43: 2216–2228
[4] Liu P, Li Q, Li Z, et al. Anatomy of subsidence in Tianjin from time series InSAR. Remote Sensing, 2016, 8(3): 266
[5] Jiang L, Bai L, Zhao Y, et al. Combining InSAR and hydraulic head measurements to estimate aquifer parameters and storage variations of confined aquifer system in Cangzhou, North China Plain. Water Resour-ces Research, 2018, 54(10): 8234–8252
[6] Massonnet D, Rossi M, Carmona C, et al. The dis-placement field of the Landers earthquake mapped by radar interferometry. Nature, 1993, 364: 138–142
[7] Zebker H A, Rosen P A, Goldstein R M, et al. On the derivation of coseismic displacement fields using dif-ferential radar interferometry: the Landers earthquake. Journal of Geophysical Research: Solid Earth, 1994, 99(B10): 19617–19634
[8] Li F K, Goldstein R M. Studies of multibaseline spaceborne interferometric synthetic aperture radars. IEEE Transactions on Geoscience and Remote Sen-sing, 1990, 28(1): 88–97
[9] Zebker H A, Villasenor J. Decorrelation in interfero-metric radar echoes. IEEE Transactions on Geoscience and Remote Sensing, 1992, 30(5): 950–959
[10] Ferretti A, Prati C, Rocca F. Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(1): 8–20
[11] Lyons S, Sandwell D. Fault creep along the southern San Andreas from interferometric synthetic aperture radar, permanent scatterers, and stacking. Journal of Geophysical Research: Solid Earth, 2003, 108(B1): 2047
[12] Werner C, Wegmuller U, Strozzi T, et al. Interferome-tric point target analysis for deformation mapping // Proceedings of the IGARSS 2003. Toulouse, 2003: 4362–4364
[13] Hooper A, Zebker H, Segall P, et al. A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers. Geophysical Research Letters, 2004, 31(23): L23611
[14] Hooper A, Segall P, Zebker H. Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to Volcán Al-cedo, Galápagos. Journal of Geophysical Research, 2007, 112(B7): B07407
[15] Hooper A, Bekaert D, Spaans K, et al. Recent advan-ces in SAR interferometry time series analysis for measuring crustal deformation. Tectonophysics, 2012, 514: 1–13
[16] Williams S, Bock Y, Fang P. Integrated satellite in-terferometry: tropospheric noise, GPS estimates and implications for interferometric synthetic aperture ra-dar products. Journal of Geophysical Research: Solid Earth, 1998, 103(B11): 27051–27067
[17] Jolivet R, Agram P S, Lin N Y, et al. Improving InSAR geodesy using global atmospheric models. Journal of Geophysical Research: Solid Earth, 2014, 119(3): 2324–2341
[18] Smith E K, Weintraub S. The constants in the equation for atmospheric refractive index at radio frequencies. Proceedings of the IRE, 1953, 41(8): 1035–1037
[19] Hanssen R F. Radar interferometry: data interpretation and error analysis. Dordrecht: Springer Science & Business Media, 2001
[20] Bekaert D, Walters R, Wright T, et al. Statistical com-parison of InSAR tropospheric correction techniques. Remote Sensing of Environment, 2015, 170: 40–47
[21] Hersbach H, Bell B, Berrisford P, et al. The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society, 2020, 146: 1999–2049
[22] Tymofyeyeva E, Fialko Y. Mitigation of atmospheric phase delays in InSAR data, with application to the eastern California shear zone. Journal of Geophysical Research: Solid Earth, 2015, 120(8): 5952–5963
[23] Wang K, Fialko Y. Observations and modeling of co-seismic and postseismic deformation due to the 2015 Mw 7.8 Gorkha (Nepal) earthquake. Journal of Geo-physical Research: Solid Earth, 2018, 123(1): 761–779
[24] 郭海朋, 白晋斌, 张有全, 等. 华北平原典型地段地面沉降演化特征与机理研究. 中国地质, 2017, 44(6): 1115–1127
[25] Werner C, Wegmüller U, Strozzi T, et al. Gamma SAR and interferometric processing software // Proceedings of the ERS-ENVISAT Symposium. Gothenburg, 2000: 1620
[26] Farr T G, Rosen P A, Caro E, et al. The shuttle radar topography mission. Reviews of Geophysics, 2007, 45 (2): RG2004
[27] 徐杰. 渤海湾盆地构造及其演化. 北京: 地震出版社, 2015
[28] 张雅芳, 郭英, 沈彦俊, 等. 华北平原种植结构变化对农业需水的影响. 中国生态农业学报, 2020, 28(1): 8–16
[29] Li M, Sun J, Xue L, et al. Characterization of aquifer system and groundwater storage change due to south-to-north water diversion project at Huairou ground-water reserve site, Beijing, China, using geodetic and hydrological data. Remote Sensing, 2022, 14: 3549
[30] Amos C B, Audet P, Hammond W C, et al. Uplift and seismicity driven by groundwater depletion in central California. Nature, 2014, 509: 483–486
Wide-area InSAR Time Series Analysis Technique for Monitoring of Surface Deformation in the North China Plain
Abstract Based on the SAR satellite images from 2015 to 2019, this paper proposes a novel wide-area InSAR time series analysis technology for high-precision and continuous monitoring of the surface deformation in the North China Plain. Firstly, standard interferometric processing is carried out on SAR images to obtain the interferogram. Then the parallelized StaMPS method is used to extract all PS pixels in the interferogram and obtain the deformation time series with full resolution. Next the atmospheric noise is estimated and removed using the joint atmospheric correction method, through combination of the atmospheric model correction and common scene stacking. After this processing flow, we successfully obtain a solution of large spatial scale, long-term span, full spatial resolution, and high precision deformation in the North China Plain, and detect the subsidence deformation pattern with a strong signal of up to 100 mm/a in central plain due to groundwater extraction. Compared with existing algorithms, the parallelized StaMPS method saves at least 60% of computing time through real-time distribution among multiple computing nodes. Besides, the joint atmospheric correction method can remove about 74.3% of the atmospheric noise, and its effectiveness is significantly higher than that of using the two correction methods alone. This wide-area InSAR time series analysis technique can effectively realize high-precision continuous monitoring of large-scale surface deformation.
Key words wide-area InSAR time series analysis technique; parallelized StaMPS method; joint atmospheric cor-rection method; surface deformation; North China Plain