论文信息：Liu, L., Yang, B., Zhang, Y., Xu, Y., Peng, Z., & Yang, D. (2024). Calculating sensitivity or gradient for geophysical inverse problems using automatic and implicit differentiation. Computers & Geosciences, 193, 105736. https://doi.org/10.1016/j.cageo.2024.105736

在地球物理反演里，灵敏度矩阵与目标函数梯度是几乎所有确定性优化算法的"燃料"，它们的准确与否直接决定反演能否收敛、收敛得多快。然而，对每个算子手工推导、编码并调试其导数公式，往往是耗时且极易出错的一环。本文（第一作者刘连，来自南方科技大学；浙江大学团队共同完成，张壹参与其中）提出并系统检验了一套自动微分 + 隐式微分（Automatic Differentiation with Implicit Differentiation，简称 ADID）的方案，能够自动而精确地算出灵敏度与梯度，并将其接入 3D 电磁法与直流电阻率反演。下面围绕这篇 Computers & Geosciences 论文，梳理其动机、方法与数值结果。

灵敏度/梯度为什么是反演绕不开的"硬骨头"

地球物理反演是非线性的、且因观测覆盖不足与模型空间存在零空间而严重欠定。求解这类问题要么走确定性路线（如 Occam、非线性共轭梯度、L-BFGS），要么走随机路线（如贝叶斯 / 转维抽样）。无论哪种，确定性框架下都要反复计算两类量：正演响应与灵敏度/梯度——后者正是目标函数对模型参数的偏导数。

对灵敏度/梯度的计算，传统上有三条路径（McGillivray & Oldenburg, 1990）：

  • 扰动法（perturbation）：对每个模型参数额外求解一个全新的正演问题，代价随参数个数线性暴涨。
  • 灵敏度方程法（sensitivity equation）：计算量与扰动法相当。
  • 伴随方程法（adjoint equation，AE）：只需解少量额外的子正演来拼出灵敏度矩阵，或多次解一个伴随正演来得到梯度向量，是过去数十年里电磁反演（MT、DC、CSEM）的首选。

伴随法与灵敏度方程法都利用了线性方程组系数矩阵可复用的特性，但它们的痛点在于：需要手工推导、编码并调试相关的导数公式。尤其是当代码"看起来能跑、结果不离谱"时，隐藏在导数实现里的错误格外难揪出来。

方法：让自动微分"接管"难啃的线性求解部分

自动微分（AD）本质上是链式法则 + 编程技术：复杂的函数被拆成一系列初等代数运算（指数、对数、三角等），它们的导数都是良定义的，再用链式法则拼起来即可得到任意复合函数的导数，而且保证精确、无需手工推导。自 Sambridge 等（2007）把 AD 引入地球物理以来，它已广泛用于地震全波形反演等优化问题。

但 AD 在地球物理中一直"卡壳"在一点：地球物理正演（如 MT/DC 的电磁场）通常要求解线性方程组 Ku = s，其中 K 是系数矩阵、u 是物性场。线性求解器本身是迭代过程，不属于初等代数函数，AD 无法显式地对它的解求导。常见的两个补丁都不理想：

  • 摊开求解器（unroll，如把 BiCGstab、QMR 的迭代全程放进 AD 图里）：要求用 AD 库重写求解器，开销巨大，且 BiCGstab 这类线性求解器自身的数值误差会在 AD 过程中被放大，导致导数偏差；
  • 自定义 AD 算子：即本文采用的思路——用隐式微分（ID）把线性方程组的解视为由 F(K,u,s)=Ku−s=0 定义的隐函数。由全微分与隐函数定理可严格导出解 u 对 s 和 K 的导数，从而只需多解一次与伴随方程结构相同的线性方程组（利用 K 通常对称的性质），就能把 ∂y/∂s 与 ∂y/∂K 全部拿到手，再交给 AD 的链式法则完成剩余步骤。

这两者的数学计算量与经典 AE 完全一致，唯一差别是执行顺序。因此在理论上，ADID 不会有任何额外的计算误差，也无需重写求解器——这正是它的核心价值：精度有保证，人力开销几乎为零。

实现与算例验证

作者用 TensorFlow（内置 AD 模块 + 自定义 ID 模块）构建了全套正演与求导框架，源码开源在 GitHub（Geo-LianLiu/ADID）。论文依次用四个算例检验 ADID：

  1. 透明简单样例：构造一个含三角函数、含线性方程组的玩具正演，一步步展示 ADID 如何通过链式法则把 ∂y/∂s、∂y/∂K 算出来。
  2. 合成 2D MT：一个 100 Ω·m 半空间中含 10 Ω·m 低阻块与 1000 Ω·m 高阻块，21 个测点、21 个频点。ADID 计算的梯度与灵敏度同经典 AE 完全吻合，相对误差均小于 10⁻⁴（图1）。
  3. 3D MT：三层模型，网格达 114×114×114。ADID 的灵敏度与解析灵敏度高度一致，正演阻抗与解析解相对误差小于 0.5%，灵敏度相对误差小于 1.0%。
  4. 3D DC 电阻率：两层模型、pole-pole 观测，网格 121×121×110，用开源 SimPEG 提供 AE 参照。ADID 求得的灵敏度与解析解吻合，相对误差同样小于 1.0%（图2）。
ADID与AE计算的梯度对比
图1　合成 2D MT 算例中，ADID（a）与经典伴随方程法 AE（b）计算的梯度对比（引自论文 Figure 2）

两者没有可见差异；文中报告其相对误差均小于 10⁻⁴。这直观说明了 ADID 与 AE 在数学上等价、仅执行顺序不同。

3D DC灵敏度与解析解对比
图2　3D DC 电阻率算例中，ADID 求出的视电阻率灵敏度（实线）与解析解（叉号）对比（引自论文 Figure 9）

(a) 对第一层电阻率、(b) 对第二层电阻率的灵敏度。虚线为相对误差，均小于 1.0%，说明 ADID 与解析解高度一致。

讨论：ADID 的优势与适用边界

作者指出，AE 几十年来一直是灵敏度/梯度计算的默认首选，但 ADID 在特定场景下明显更有优势：电磁与直流问题里，∂K/∂m、∂s/∂m 这类偏导数很容易写错、且实现费时费力（比如常见省掉 ∂s/∂m 的"取巧"会引入误差），而 AD 能自动、精确地处理它们。若用到更复杂、非可微易错的吸收边界（如 GPR 模拟中的 PML），ADID 的优势会更加突出。

另一方面，ADID 的代价是要求用户用 AD 库（如 TensorFlow）来写正演代码，不过这些库的语法与 Numpy、Eigen 等常用库相似，学习成本不高；对已经用 AE 写好导数代码的研究者，ADID 也是一个现成的校验与增强工具。作者还提醒，AD 只擅长初等算子，对 abs(x) 在不可导点（x=0）等情形需要自定义导数；TensorFlow 对 loop、if-else 等控制流也做了良好处理。理论上 ADID 与 AE 计算速度与内存开销相当，大尺度 3D 问题尤为接近；只是小尺度问题中两者的开销分配略有不同。

结论与展望

这项工作把 ADID 成功用于地电与地电磁两大典型反演问题（MT 与 DC），示范了如何自动、精确地计算灵敏度和梯度：ADID 与 AE 执行相同的数学运算、只是顺序不同，因此精度有保障、人力成本大幅下降；其隐式微分模块纯粹基于数学原理、不依赖具体地球物理应用，因而对 Helmholtz 方程（MT）、Poisson 方程（DC）之外的大量相似地球物理问题同样适用。作者认为，ADID 不仅能显著解放研发人员在求导、编码与排错上的精力，更在打造工业级地球物理反演软件中扮演关键角色——让导数这一基础环节变得可靠、可复用、可验证。

延伸阅读与引用

  • Sambridge, M., Rickwood, P., Rawlinson, N., & Sommacal, S. (2007). Automatic differentiation in geophysical inverse problems. Geophysical Journal International, 170, 1–8.
  • McGillivray, P. R., & Oldenburg, D. W. (1990). Methods for calculating Fréchet derivatives and sensitivities for the non-linear inverse problem: a comparative study. Geophysical Prospecting, 38, 499–524.
  • Baydin, A. G., Pearlmutter, B. A., Radul, A. A., & Siskind, J. M. (2018). Automatic differentiation in machine learning: a survey. arXiv:1502.05767.
  • Blondel, M., et al. (2021). Efficient and modular implicit differentiation. arXiv:2105.15183.
  • 本文代码（ADID，Python/TensorFlow）：https://github.com/Geo-LianLiu/ADID

本文为对该论文的中文解读，图片均引自论文原文。