水文动态时移电阻率层析成像反演的自动微分框架

# 面向水文动态时间推移电阻率层析成像反演的自动微分框架
来源：https://arxiv.org/html/2608.14661

Pu Yang¹, Zhengyang Fang¹,², Yuxin Liu¹, Xuan Su³, Deshan Feng¹, Hang Chen²  
¹ 中南大学地球科学与信息物理学院，长沙 410083，中国  
² 爱荷华大学地球、环境与可持续发展学院，爱荷华市，IA 52245，美国  
³ 湖南工商大学管理科学与工程学院，长沙 410205，中国  
[[email protected]](mailto:[email protected])

**摘要。** 时间推移电阻率层析成像（TL-ERT）可提供地下水文变化的空间分布信息。然而，对长时间监测序列进行反演计算要求很高。修改数据失配、正则化、模型参数化或岩石物理转换可能也需要新的梯度推导和单独的实现。本文提出了 AD-TLERT，这是一个基于自动微分的、统一的、经 GPU 加速的时间推移 ERT 反演框架。该框架将模型参数化、可微岩石物理转换、正演模拟、数据失配、正则化和辅助约束整合到一个计算链中。因此，不同的反演公式可以重用相同的偏微分方程导数实现，而无需为每种情况重新推导完整的 ERT 敏感度。与 pyGIMLi 的比较显示，其正演响应、梯度和恢复的电阻率模型高度一致。在所测试的配置下，AD-TLERT 实现了约 51 倍的加速。合成实验表明，反演选择会影响所恢复异常的幅度、形态和时间行为。通过将梯度传播通过嵌入的岩石物理关系，AD-TLERT 实现了直接含水率反演，并且对于所测试的模型，其结果比反演后转换更为准确。一项实际应用进一步展示了如何结合 ERT、温度和土壤水分观测来成像融雪驱动的山坡湿润过程。AD-TLERT 为时间推移 ERT 反演和水文解释提供了一个高效且灵活的框架。

**关键词：** 自动微分；电阻率层析成像；时间推移反演；含水率估计；水文地球物理学

通讯作者：Hang Chen, [email protected]

## 1. 引言
时间推移电阻率层析成像（TL-ERT）广泛应用于监测地下动态过程，包括水分入渗、地下水流动、溶质运移、冻融循环以及热-水力-力学过程（Daily 等，1992；Kemna 等，2002；Cassiani 等，2009；Krautblatter 等，2010；Chambers 等，2014；Uhlemann 等，2017；Chen 等，2024，2025）。重复性测量提供了对孔隙流体含量、盐度、温度和岩性敏感的电性变化的空间分布观测（Singha 等，2015；Binley 等，2015）。这些观测补充了基于点的水文和环境测量（Chen，2022）。然而，在许多应用中，主要关注的并非某一时刻的绝对电阻率，而是与过程相关的变量（如含水率或饱和度）的时空演化（Binley 等，2002；Tso 等，2019）。因此，时间推移 ERT 反演需要恢复电性变化并将其与具有物理意义的地下状态联系起来。在更大的时空尺度上，过程模拟可以进一步指导长期电性监测数据的解释（Robinson 等，2024）。随着监测序列的增长，这项任务变得计算量巨大。ERT 反演是非线性且不适定的，且时间相关信号相对于背景电阻率结构可能很小（LaBrecque 等，1996）。重复性测量还受到测量噪声、与接触条件和电极位置变化相关的系统误差以及环境变异的影响（Tso 等，2017；Oldenborger 等，2005）。

现有的时间推移反演方法大致可分为基于差异的方法和时间耦合方法。基于差异的反演估计相对于基线数据集或模型的变化，并可结合空间协方差先验（LaBrecque 和 Yang，2001；Hermans 等，2016）。时间耦合反演通过连续模型之间的约束联合估计多次测量，包括同时和全序列 4D 公式（Kim 等，2009；Hayley 等，2011；Loke 等，2014）。后来的扩展引入了自适应时间加权、Lp 范数惩罚以及结合空间和时间约束（Karaoulis 等，2011；Kim 等，2013；Karaoulis 等，2014）。对于长时间序列，窗口化反演将这种耦合限制在短的、重叠的子集上，在保持局部时间一致性的同时减小了每个反问题的规模（Wilkinson 等，2022）。然而，全序列和窗口化公式都需要在多次测量、电流源和波数下进行重复的正演和敏感度计算。由此产生的计算成本限制了密集监测数据的反演以及对替代公式的系统评估。所恢复的变化也强烈依赖于反问题的公式化方式。数据失配函数决定了测量误差和大残差如何影响模型更新（Farquharson 和 Oldenburg，1998）。空间和时间正则化影响恢复异常的幅度、形态和平滑度（Constable 等，1987；Portniaguine 和 Zhdanov，1999；Fiandaca 等，2015）。优化方法影响收敛和计算成本，而模型参数化则决定了反演是在电阻率、电导率、相对变化或与过程相关的变量空间中进行。没有一种反演设置是普遍最优的，因为合适的公式取决于数据质量、测量覆盖范围、预期的异常几何形状以及先验信息（Caterina 等，2014；Hermans 等，2016）。

成熟的框架，如 SimPEG，为地球物理模拟和反演提供了模块化组件（Cockett 等，2015）。然而，将时间推移 ERT 实现扩展到新的时空目标、可微岩石物理映射和 GPU 稀疏求解器，仍然可能需要特定问题的导数和实现工作。更改一个组件可能需要新的梯度推导或单独的实现，这使得受控比较变得困难。

一个相关的限制出现在水文解释中。通常先反演电阻率，然后将其转换为含水率或其他参数（Tso 等，2019）。由于数据拟合和正则化是在电阻率空间中应用的，平滑和幅度偏差可能会非线性地传播到转换后的估计中。互补地球物理数据的联合反演可以改善地下表征（Linde 等，2006），而耦合水文地球物理反演可以通过将过程和岩石物理模型嵌入反演中，更直接地估计水文状态或参数（Hinnell 等，2010；Pleasants 等，2022）。当存在合适的到反演变量的映射时，井中测量、传感器观测和其他辅助数据同样可以作为约束引入，但每个添加的约束都可能需要一个专门的耦合项及其导数。

自动微分（AD）提供了一条通往更灵活实现的途径（Baydin 等，2018；Margossian，2019）。当模型参数化、正演模拟、数据失配、正则化和辅助约束表示为可微操作时，梯度可以通过计算链一致地构建。对于稀疏有限元问题，AD 可以与线性系统解的隐式微分相结合，避免通过内部因式分解步骤进行微分。基于 GPU 的稀疏线性代数可以进一步加速时间推移反演所需的重复正伴随求解。然后，可以在相同的正演和梯度实现中评估替代的失配函数、约束、优化策略和岩石物理参数化（Sambridge 等，2007；Liu 等，2025）。

尽管 AD 最近已被应用于 ERT 反演（Rincón 等，2026），但其与 GPU 加速的 2.5D 时间推移反演的集成、反演选择的系统评估、水文变量的直接反演以及辅助观测的灵活整合仍然有限。在本研究中，我们开发了 AD-TLERT，这是一个用于时间推移 ERT 中水文动态反演的 GPU 加速可微框架。该框架结合了有限元正求解器、隐式微分、GPU 稀疏线性代数、时间耦合、空间和时间正则化、辅助约束以及灵活的模型转换。我们首先将其正演响应、梯度、恢复的电阻率模型和计算性能与 pyGIMLi 进行了基准测试（Rücker 等，2017）。然后，在共同的合成设置下，比较了代表性的数据失配函数、正则化策略和优化算法。进一步将基于已建立的电导率-饱和度公式的可微岩石物理关系（Archie，1942；Waxman 和 Smits，1968）嵌入反演中，以比较直接含水率反演与传统的反演后转换。最后，该框架应用于融雪驱动山坡湿润过程的实际观测，使用温度和土壤水分观测作为辅助约束。它为测试反演假设和解释水文动态提供了一个高效且可扩展的框架。

*参考插图*  
图 1。面向水文动态时间推移 ERT 反演的自动微分框架工作流程。步骤 1 定义地下模型和测量配置。步骤 2 使用 GPU 加速有限元求解器计算合成 ERT 数据。步骤 3 使用所选目标函数评估观测数据与合成数据之间的差异。步骤 4 通过隐式自动微分计算目标函数相对于反演变量的梯度。步骤 5 使用自适应优化方案更新模型。在路线 1 中，直接更新电阻率模型。在可选的路线 2 中，电阻率模型通过岩石物理参数化表示，相应的岩石物理参数在反演过程中更新。

## 2. 方法论
本节介绍 AD-TLERT 底层的数学公式和详细推导。如图 1 所总结，该框架将时间推移反演组织为五个可微阶段：模型定义、GPU 加速有限元正演模拟、失配评估、通过隐式自动微分计算梯度以及自适应模型更新。在这个模块化的计算图中，可以选择和组合替代的目标函数、优化算法以及空间和时间正则化策略，而无需为每种反演公式重新推导完整的 ERT 敏感度。该框架还支持两种更新路线：直接反演电阻率模型，以及可选的岩石物理路线，其中水文或岩石物理参数被映射到电阻率，并通过相同的梯度链进行更新。以下各小节构建了时间推移反问题，推导了 2.5D ERT 正演算子的隐式梯度，并描述了可微岩石物理参数化。

### 2.1. 可微时间推移 ERT 反演
电阻率层析成像（ERT）正演模拟将地下电导率分布映射到相应的测量响应（Dey 和 Morrison，1979）。对于电导率 $\sigma$，电位 $u$ 满足
$$-\nabla\cdot(\sigma\nabla u)=q,$$
其中 $q$ 表示施加的电流源。离散化并投影到接收器位置后，预测数据写为
$$\mathbf{d}_{\mathrm{pred}}=\mathcal{F}(\mathbf{m}),$$
其中 $\mathbf{m}=\{\mathbf{m}_{t}\}_{t=1}^{N_t}$ 包含所有测量时刻的模型参数。时间推移反问题表述为
$$\min_{\mathbf{m}}\ \Phi(\mathbf{m})=\sum_{t=1}^{N_{t}}\Phi_{d}^{(t)}\left(\mathcal{F}(\mathbf{m}_{t}),\mathbf{d}_{\mathrm{obs}}^{(t)}\right)+\alpha_{s}\sum_{t=1}^{N_{t}}\Phi_{s}(\mathbf{m}_{t})+\alpha_{t}\sum_{t=2}^{N_{t}}\Phi_{t}(\mathbf{m}_{t},\mathbf{m}_{t-1})+\alpha_{c}\Phi_{c}(\mathbf{m}),$$
其中 $\Phi_{d}$、$\Phi_{s}$ 和 $\Phi_{t}$ 分别表示数据失配项、空间正则化项和时间正则化项。可选项 $\Phi_{c}$ 纳入了额外约束，如井中或传感器观测、岩石物理先验，或联合反演中使用的耦合项。所有模型转换和 ERT 正