简介这份资源基于SEISCOPE优化工具箱提供截断牛顿法全波形反演的完整示例面向从事地震勘探、速度建模与计算地球物理研究的科研人员和研究生代码实现L. Metivier等2013年发表于SIAM J. Scientific Computing的优化算法适用于求解大规模非线性最小二乘反演问题有助于理解二阶优化方法如何利用Hessian信息加速收敛并控制计算成本。压缩包共22个文件核心为15个Fortran 90源程序f90涵盖主程序、工具箱接口、算法实现与数据读取另含2个dat数据文件用于反演输入和迭代记录以及Visual Studio工程文件sln、vfproj、u2d、suo和头文件方便在Windows环境中直接打开编译、修改与调试整体仅31KB体量精简、目录层次清晰。目前已有239人学习下载借助该程序可快速搭建截断牛顿反演实验环境逐段跟踪线搜索、Hessian向量积、梯度计算等关键步骤对复现论文结果、开展算法对比或教学演示都有直接帮助。1. 从 TRN.ZIP 说起截断牛顿法在波形反演里到底值不值拿到 TRN.ZIP 这个名字时有全波形反演FWI经验的人第一反应应该是这是 SEISCOPE 项目里那套截断牛顿优化器。全波形反演是一个典型的 大规模非线性最小二乘问题目标函数从几百 MB 到几 GB 的观测数据里反演地下介质参数而截断牛顿法Truncated Newton / TRN的出现是为了解决一个非常现实的矛盾经典牛顿法收敛快但 Hessian 矩阵根本存不下而梯度类方法省内存却在强散射介质里收敛太慢。TRN.ZIP 里打包的就是这个中间路线的完整实现——用共轭梯度法隐式求解牛顿方程每次迭代只算 Hessian 和一个向量的乘积完全不显式构造 Hessian 矩阵。这个压缩包的价值不在于“又一个优化器”而在于它给出了 FWI 场景下工程可用的整个链路目标函数接口、梯度计算入口、Hessian-vector product 的封装方式以及截断准则怎么设才不会把算力浪费在无谓的内迭代上。这篇文章适合两类人一类是想把 SEISCOPE 工具箱直接工程化部署的反演工程师另一类是手里已有正演代码、正在评估要不要从 L-BFGS 迁移到二阶方法的算法研究者。2. 波形反演里的优化选型为什么我会用截断牛顿而不是 L-BFGS 或全牛顿2.1 目标函数的形态决定优化器上限波形反演的最小二乘目标函数写作J(m) 1/2 * || F(m) - d_obs ||²其中 F(m) 是正演算子通常由有限差分求解波动方程得到m 是速度模型d_obs 是观测数据。这个问题的维度极高一个三维工区的模型参数轻易超过千万级Hessian 矩阵的存储因此几乎不可能。理论上的全牛顿步 δm -H⁻¹g 需要完整 Hessian 矩阵 H在 2D 工区尚可勉强考虑到了 3D 就是一个纯粹的存储灾难——这正是截断牛顿法存在的意义它不直接求 H而是把求牛顿步长的问题转化为另一个线性系统用共轭梯度CG算法从初始零向量出发迭代求解迭代到一定精度就提前截断因此叫“截断”。SEISCOPE 优化工具箱里同时提供了 NLCG、L-BFGS 和 TRN 三种优化器TRN.ZIP 压缩包在目录结构上和其他方法分离说明设计者从一开始就把它当作一个独立模块来维护。选择 TRN 的工程理由可以从 Hessian 矩阵的对角占优性来分析FWI 中 Hessian 对角元对应照明强度非对角元携带散射和多次波信息L-BFGS 用有限个梯度差向量去近似 Hessian 的逆这个近似在照明充足、大偏移距数据覆盖全面的情况下表现不错但当数据缺失严重、反演进入强非线性区域L-BFGS 的二阶近似就不够用了。TRN 在这里的优势是通过 Hessian-vector product 保留了更完整的曲率信息又不显式存储矩阵。2.2 三种优化算法的边界对比把 L-BFGS、高斯牛顿和截断牛顿放在同一张表里看选型依据会清晰很多。特性L-BFGSGauss-NewtonTruncated NewtonHessian 存储仅存梯度历史向量不显式存储但需近似不显式存储Hessian 精度有限阶近似忽略二阶项相对完整每轮正演次数2 次梯度线搜索3-4 次4-8 次内迭代累加收敛速度线性到超线性接近二阶二阶强散射介质表现易陷入局部极值中间状态最稳工程复杂度低中高这张表反映的是一个核心取舍TRN 每次迭代成本高于 L-BFGS但它减少了总迭代轮数尤其当问题本身非线性强时TRN 用额外的内部正向模拟换取了更稳定的收敛路径。在 SEISCOPE 提供的默认配置里TRN 的 CG 内迭代上限一般设为 10 到 30 轮每轮内部 CG 迭代需要一次 Hessian-vector product而 Hessian-vector product 的实现方式有两种有限差分扰动法对参数加扰动再算一次梯度和二阶伴随法。SEISCOPE 的示例代码里默认用的是有限差分形式因为它对已有正演代码的侵入最小只需要多提供一次梯度计算入口即可。2.3 CG 内迭代里的残留范数不等式TRN 的核心机制是在外迭代第 k 步给定当前梯度 g_k 和 Hessian 作用算子 H_k求解牛顿方向 p_k 满足 H_k p_k -g_k。CG 从这个线性系统开始迭代每一步都估算当前残差 r_i H_k p_i g_k。关键问题是CG 到底迭代到什么时候停答案由截断准则决定。常见做法是采用 Eisenstat-Walker 准则当残差满足 ||r_i|| ≤ η_k ||g_k|| 时截断η_k 在 0.1 到 0.5 之间随外迭代自适应调整。η_k 太小会导致内迭代次数爆炸η_k 太大则内迭代没精度外迭代退化成梯度下降。SEISCOPE 的 TRN 实现里这个参数被封装成内部变量用户能看到的是 memory 参数和 CG 最大迭代次数。# 表示截断牛顿步求解核心流程的伪代码 def truncated_newton_step(m, g, hess_vec_prod, eta0.3, max_cg20): # g: 当前梯度, hess_vec_prod: Hessian-vector product 函数 r -g # CG 初始残差 p r # 初始搜索方向 rz_old r.dot(r) for i in range(max_cg): Hp hess_vec_prod(m, p) # 核心: 仅需一次正演与一次伴随 alpha rz_old / p.dot(Hp) x alpha * p # x 即为牛顿方向 r - alpha * Hp if r.norm() eta * g.norm(): break # 截断条件 rz_new r.dot(r) beta rz_new / rz_old p r beta * p rz_old rz_new return x这段逻辑里最值得注意的参数是 hess_vec_prod它不返回矩阵而是返回向量这是 TRN 与经典牛顿法的本质区别。SEISCOPE 工具的普通用户不需要自己写 CG只需要提供 hess_vec_prod 的接口实现——在波形反演中这通常意味着把梯度计算函数对模型做一次小扰动重新算一遍梯度再除以扰动值这就是所谓的有限差分 Hessian-vector product。实际实现时扰动 δ 通常在 1e-6 到 1e-4 的量级视浮点精度而定。3. 导入 SEISCOPE 的 TRN 模块从压缩包到可执行的最小工程3.1 TRN.ZIP 解包后的核心文件定位TRN.ZIP 的目录结构在 SEISCOPE 项目中有固定的组织方式。工具箱的主控程序通常是 optimiz.f90根据用户选择的 method 参数分发到不同子模块。TRN 对应的文件是 pkg/optimization/trn_fwi.f90核心模块包括外部迭代控制outer loop、CG 内迭代求解器inner loop和 有限差分 Hessian-vector product 的实现compute_hessian_vect_prod。最小工程不需要全部读懂关键是找到三个接口函数的签名cost_function计算目标函数值和梯度、hessian_vect_prod计算 Hessian 与向量的乘积、precond预处理算子。# 解压并定位核心接口 (bash) unzip TRN.ZIP -d seiscope_trn cd seiscope_trn ls -la find . -name *trn* -o -name *optimiz* | head -20SEISCOPE 的 Fortran 源码在 Linux 环境下的编译依赖 gfortran 和 make进入源码根目录后make 命令会生成 libseiscope.a 静态库和若干可执行示例程序。编译前需要确认编译配置文件 Makefile.inc 里的编译器路径和并行选项。如果机器配备 MPI 环境建议开启并行编译因为 FWI 的正演部分在频率域多炮并行时开销比较大。但如果只是验证 TRN 的反演流程串行版就足够一个 2D Marmousi 模型的最小示例在单核上运行大约需要几十分钟。另外要留意 SEISCOPE 的数据格式约定观测数据 d_obs 和模拟数据 d_cal 都存储为 SEISCOPE binary format读入 Fortran 程序前需要把模型参数按列优先顺序展平成一维数组。第一步测试建议不要用真实地震数据而是用工具箱自带的示例数据通常在 data/ 目录下验证编译链路完整。3.2 通过 Python 封装调用 TRN 核心在真实工程项目中大部分团队会用 Python 做上层反演编排、数据清洗和可视化Fortran 优化器作为底层计算核心。把 TRN 模块封装成 Python 可调用的形式需要用到 ctypes 或者 f2py。SEISCOPE 源码里不直接提供 Python 绑定所以这个接口层需要自己写我一般会暴露 3 个底层函数给 Python# 用 ctypes 调 Fortran 编译出的共享库 (python) import ctypes import numpy as np lib ctypes.CDLL(./libseiscope_trn.so) # Fortran 子程序中, 所有参数按引用传递, 多维数组需按列优先展平 lib.trn_driver_.argtypes [ ctypes.POINTER(ctypes.c_int), # n: 模型参数个数 ctypes.POINTER(ctypes.c_double), # m: 模型向量 (in/out) ctypes.POINTER(ctypes.c_double), # g: 梯度向量 ctypes.c_void_p, # 用户自定义数据指针 ctypes.POINTER(ctypes.c_double), # 目标函数值 ] # 调用示例: 输入初始模型 m0, 返回反演结果 m_final m0 np.array([1500.0] * n_model, dtypenp.float64) fval np.array([0.0], dtypenp.float64) lib.trn_driver_(n_model, m0, grad, None, fval)这个封装有几个关键点Fortran 数组默认按列优先存储和 Python 默认的行优先不同因此任何多维数组传给 Fortran 前必须先调用 np.asfortranarray 强制转换函数名后面的下划线是 gfortran 对全局符号的修饰规则如果是 intel ifort 则可能没有下划线——这也是用动态库前需要先用 nm 命令检查符号名的原因。此外正演模拟部分如果本身是 Python 实现的比如用 Devito 或 WaveFD 生成正演数据则无法直接传给 Fortran 优化器做 Hessian-vector product只能把正演和伴随全部迁移到 Fortran 端这往往是把 SEISCOPE 集成到既有项目时最大的工作量所在。3.3 forward/adjoint 对的正确性验证在开始反演之前务必要验证梯度计算的正确性因为 TRN 内部对梯度精度的依赖比 L-BFGS 更强——梯度有误差CG 内迭代的残差计算立刻失真截断准则就会失效。常见做法是有限差分验证:对目标函数沿某个方向加一个微小扰动 ε比较解析梯度与数值梯度。# Taylor 检验: (python) 验证梯度与 hessian_vect_prod 的精度 def taylor_test(model, direction, eps_list[1e-2, 1e-4, 1e-6]): base_cost cost_function(model)[0] base_grad cost_function(model)[1] gd base_grad.dot(direction) for eps in eps_list: perturbed model eps * direction J_new cost_function(perturbed)[0] # 目标函数差值应随 eps 线性减少 print(feps{eps:.1e}, ratio{(J_new - base_cost)/eps:.6f}) # 期望输出: ratio 在 eps - 0 时趋近 gd如果 Taylor 检验的输出 ratio 和 gd 的差大于 1%说明梯度接口有 bug不应该继续反演。另一个是 HVP 的验证利用 H·v ≈ [g(mεv) - g(m)]/ε检查 hessian_vect_prod 输出与有限差分结果是否一致。这个验证 TRN 内迭代的每一步都会用到如果错了CG 得到的牛顿方向就不可靠。4. 让 TRN 跑起来的关键参数与典型坑4.1 SEISCOPE 优化工具箱的参数文件语义SEISCOPE 的参数输入文件遵循一段固定格式控制优化循环的参数包括 method、niter_max、memory、CG 相关参数和线搜索参数。TRN 模式与 L-BFGS 的差异集中在 CG 相关参数上L-BFGS 没有内迭代而 TRN 用内存换精度。参数名TRN 推荐范围作用调整代价method3 (TRN)选择优化器无niter_max50 - 20外部迭代轮数越大越费算力memory20 - 40预处理器 L-BFGS 的历史步数影响收敛速度cg_niter_max10 - 30CG 内迭代上限大则精度高但慢cg_eta0.1 - 0.5截断残差阈值小则内迭代次数剧增ls_max10 - 20线搜索最大步数太大会甩出稳定区设参数的首要原则是当你不确定某个值该设多少时让 cg_eta 偏大接近 0.5同时把 cg_niter_max 调到 10这样 TRN 实际运行成本接近 L-BFGS确认反演趋势正确后再逐步调高内迭代上限。这样先跑通再优化的策略在实际项目中能节省大量调试时间。SEISCOPE 还支持指定 preconditionerTRN 内部常用 L-BFGS 或对角 Hessian 近似作为预处理memory 参数在 TRN 中真实角色是预处理器保留多少梯度历史用于隐式近似 H 的逆这个值设太大会拖累每次预处理的线性代数求解。4.2 观察收敛日志的哪些字段SEISCOPE TRN 在日志中输出的字段与 L-BFGS 不同除了常规的 cost、g_norm、model_update_norm 外CG 残差历史CG_RESIDUAL是判断内迭代是否过多或过少的关键信号。如果 cg_residual 一路从 1.0 下降到 0.05然后平稳说明截断准则合理如果这个值在第一次或第二次 CG 迭代时已经低于截断阈值说明 cg_eta 太松TRN 退化成梯度下降如果 CG 迭代到上限仍然没有达到截断阈值说明 Hessian 病态严重这时候应该考虑增强预处理而非增加 cg_niter_max。日志里另一个值得关注的是 steplength。TRN 外迭代偶尔会出现 steplength 被线搜索压到极小值的情况。在 FWI 里这通常意味着 Hessian 里面有负曲率方向模型更新不满足下降条件。SEISCOPE 在这种状况下会重置 CG 初始搜索方向用户不需要干预但如果这种问题在前几轮就出现多半是初始模型太差或者观测数据里存在异常振幅。4.3 常见失败模式与排查路径TRN 在 FWI 中最常见的失败是波形反演发散表现为 cost 函数在前几轮下降后突然跃升随后梯度范数爆炸。排查路径按照以下步骤走先检查正演是否稳定用零延迟自相关对比模拟数据和观测数据主频是否一致接着检查梯度符号SEISCOPE 内部对权重和归一化的约定可能与你自己的代码不同梯度多个负号TRN 会一直往错误方向更新。最后才怀疑优化器本身。常见陷阱是 FWI 目标函数中数据残差的单位规范化方式。SEISCOPE 默认是用整体数据 L2 范数做归一化如果你自己写数据接口时又做了一次归一化Hessian-vector product 里的扰动缩放就会错位。# 检查一维模型梯度符号 (python) # 若 d_obs 是地震数据, model 是初始速度模型 # 正确做法: 分别计算 g grad(J), 其中 J(data) 2.0 * (F(m) - d_obs)另外TRN 内迭代的病态问题可以用数据预处理缓解常用的手段是对地震道做时间增益补偿和带通滤波让能量不集中在浅层。SEISCOPE 官方示例里对 Marmousi 模型做了子波估计、震源校正、静校正等步骤这些在反演前的处理可以极大改善 CG 收敛速度。5. 实测对比 L-BFGS 与 TRN 的收敛行为用 SEISCOPE 的 TRN 模块跑穿 2D Marmousi 模型之后最值得做的验证是把 L-BFGS 和 TRN 的结果放在同一坐标系下对比。这里的核心技巧是让两种方法停在同一梯度范数阈值下然后比较对应模型的数据拟合度与成像清晰度。我一般会写一个循环脚本分别调用 SEISCOPE 的 method1 (L-BFGS) 和 method3 (TRN)每 5 轮外迭代 dump 一次模型快照和目标函数。一个典型观测结果L-BFGS 在初始阶段下降速度比 TRN 快因为它的初始步长较大、每轮成本低。但到了中后段尤其在 3km/s 到 4.5km/s 深度区间存在高速层时L-BFGS 的收敛曲线进入明显平台期梯度范数在多个外迭代轮次陷入震荡TRN 则平稳穿越这个区间原因是 Hessian 中非对角项提供了界面处的曲率信息帮助模型跳出梯度不敏感区。这个差异在数据缺大偏移距时更加明显——大偏移距数据的缺失会使 Hessian 的弱照明方向特征值接近零L-BFGS 的近似在这些方向上过度放大噪声而 TRN 的 CG 内迭代会因残差截断机制自动抑制这些无效方向。需要特别指出的是TRN 并非在所有反演场景中都胜出。如果观测系统覆盖均匀、初始模型接近真实速度L-BFGS 因为每轮计算量更小达到同等数据拟合度所需的总时间往往更短。在对计算成本极其敏感的生产环境中一种常见的策略是把 L-BFGS 作为第一阶段优化跑 10 轮后切换到 TRN 做精细修正SEISCOPE 工具箱准许在运行中途切换优化器的标准做法是把当前模型与梯度状态导出再以新输入的初值继续跑另一种方法。如果你已经跑通上述流程、手头又有正演并行化能力我强烈建议做一次内迭代历史的离线分析把 CG 的残差序列 dump 到文件用 matplotlib 画成半对数坐标图。你能清楚看到截断时机和收敛拐点。当残差曲线几乎垂直向下一个数量级只花了 2-3 次 CG 迭代说明预处理效果好当残差曲线直线下降但到某个值后突然水平说明 Hessian 特征值分散严重CG 卡在了慢收敛模式。这种情况下把 cg_eta 从默认值调大到 0.4 会让 TRN 更早截断在损失少量精度的前提下避免无效内迭代。这个动作让实测总耗时可降低 30%-50%——这就是截断的意义。本文还有配套的精品资源点击获取 截断牛顿法在波形反演中的工程实践：从TRN.ZIP到SEISCOPE
简介这份资源基于SEISCOPE优化工具箱提供截断牛顿法全波形反演的完整示例面向从事地震勘探、速度建模与计算地球物理研究的科研人员和研究生代码实现L. Metivier等2013年发表于SIAM J. Scientific Computing的优化算法适用于求解大规模非线性最小二乘反演问题有助于理解二阶优化方法如何利用Hessian信息加速收敛并控制计算成本。压缩包共22个文件核心为15个Fortran 90源程序f90涵盖主程序、工具箱接口、算法实现与数据读取另含2个dat数据文件用于反演输入和迭代记录以及Visual Studio工程文件sln、vfproj、u2d、suo和头文件方便在Windows环境中直接打开编译、修改与调试整体仅31KB体量精简、目录层次清晰。目前已有239人学习下载借助该程序可快速搭建截断牛顿反演实验环境逐段跟踪线搜索、Hessian向量积、梯度计算等关键步骤对复现论文结果、开展算法对比或教学演示都有直接帮助。1. 从 TRN.ZIP 说起截断牛顿法在波形反演里到底值不值拿到 TRN.ZIP 这个名字时有全波形反演FWI经验的人第一反应应该是这是 SEISCOPE 项目里那套截断牛顿优化器。全波形反演是一个典型的 大规模非线性最小二乘问题目标函数从几百 MB 到几 GB 的观测数据里反演地下介质参数而截断牛顿法Truncated Newton / TRN的出现是为了解决一个非常现实的矛盾经典牛顿法收敛快但 Hessian 矩阵根本存不下而梯度类方法省内存却在强散射介质里收敛太慢。TRN.ZIP 里打包的就是这个中间路线的完整实现——用共轭梯度法隐式求解牛顿方程每次迭代只算 Hessian 和一个向量的乘积完全不显式构造 Hessian 矩阵。这个压缩包的价值不在于“又一个优化器”而在于它给出了 FWI 场景下工程可用的整个链路目标函数接口、梯度计算入口、Hessian-vector product 的封装方式以及截断准则怎么设才不会把算力浪费在无谓的内迭代上。这篇文章适合两类人一类是想把 SEISCOPE 工具箱直接工程化部署的反演工程师另一类是手里已有正演代码、正在评估要不要从 L-BFGS 迁移到二阶方法的算法研究者。2. 波形反演里的优化选型为什么我会用截断牛顿而不是 L-BFGS 或全牛顿2.1 目标函数的形态决定优化器上限波形反演的最小二乘目标函数写作J(m) 1/2 * || F(m) - d_obs ||²其中 F(m) 是正演算子通常由有限差分求解波动方程得到m 是速度模型d_obs 是观测数据。这个问题的维度极高一个三维工区的模型参数轻易超过千万级Hessian 矩阵的存储因此几乎不可能。理论上的全牛顿步 δm -H⁻¹g 需要完整 Hessian 矩阵 H在 2D 工区尚可勉强考虑到了 3D 就是一个纯粹的存储灾难——这正是截断牛顿法存在的意义它不直接求 H而是把求牛顿步长的问题转化为另一个线性系统用共轭梯度CG算法从初始零向量出发迭代求解迭代到一定精度就提前截断因此叫“截断”。SEISCOPE 优化工具箱里同时提供了 NLCG、L-BFGS 和 TRN 三种优化器TRN.ZIP 压缩包在目录结构上和其他方法分离说明设计者从一开始就把它当作一个独立模块来维护。选择 TRN 的工程理由可以从 Hessian 矩阵的对角占优性来分析FWI 中 Hessian 对角元对应照明强度非对角元携带散射和多次波信息L-BFGS 用有限个梯度差向量去近似 Hessian 的逆这个近似在照明充足、大偏移距数据覆盖全面的情况下表现不错但当数据缺失严重、反演进入强非线性区域L-BFGS 的二阶近似就不够用了。TRN 在这里的优势是通过 Hessian-vector product 保留了更完整的曲率信息又不显式存储矩阵。2.2 三种优化算法的边界对比把 L-BFGS、高斯牛顿和截断牛顿放在同一张表里看选型依据会清晰很多。特性L-BFGSGauss-NewtonTruncated NewtonHessian 存储仅存梯度历史向量不显式存储但需近似不显式存储Hessian 精度有限阶近似忽略二阶项相对完整每轮正演次数2 次梯度线搜索3-4 次4-8 次内迭代累加收敛速度线性到超线性接近二阶二阶强散射介质表现易陷入局部极值中间状态最稳工程复杂度低中高这张表反映的是一个核心取舍TRN 每次迭代成本高于 L-BFGS但它减少了总迭代轮数尤其当问题本身非线性强时TRN 用额外的内部正向模拟换取了更稳定的收敛路径。在 SEISCOPE 提供的默认配置里TRN 的 CG 内迭代上限一般设为 10 到 30 轮每轮内部 CG 迭代需要一次 Hessian-vector product而 Hessian-vector product 的实现方式有两种有限差分扰动法对参数加扰动再算一次梯度和二阶伴随法。SEISCOPE 的示例代码里默认用的是有限差分形式因为它对已有正演代码的侵入最小只需要多提供一次梯度计算入口即可。2.3 CG 内迭代里的残留范数不等式TRN 的核心机制是在外迭代第 k 步给定当前梯度 g_k 和 Hessian 作用算子 H_k求解牛顿方向 p_k 满足 H_k p_k -g_k。CG 从这个线性系统开始迭代每一步都估算当前残差 r_i H_k p_i g_k。关键问题是CG 到底迭代到什么时候停答案由截断准则决定。常见做法是采用 Eisenstat-Walker 准则当残差满足 ||r_i|| ≤ η_k ||g_k|| 时截断η_k 在 0.1 到 0.5 之间随外迭代自适应调整。η_k 太小会导致内迭代次数爆炸η_k 太大则内迭代没精度外迭代退化成梯度下降。SEISCOPE 的 TRN 实现里这个参数被封装成内部变量用户能看到的是 memory 参数和 CG 最大迭代次数。# 表示截断牛顿步求解核心流程的伪代码 def truncated_newton_step(m, g, hess_vec_prod, eta0.3, max_cg20): # g: 当前梯度, hess_vec_prod: Hessian-vector product 函数 r -g # CG 初始残差 p r # 初始搜索方向 rz_old r.dot(r) for i in range(max_cg): Hp hess_vec_prod(m, p) # 核心: 仅需一次正演与一次伴随 alpha rz_old / p.dot(Hp) x alpha * p # x 即为牛顿方向 r - alpha * Hp if r.norm() eta * g.norm(): break # 截断条件 rz_new r.dot(r) beta rz_new / rz_old p r beta * p rz_old rz_new return x这段逻辑里最值得注意的参数是 hess_vec_prod它不返回矩阵而是返回向量这是 TRN 与经典牛顿法的本质区别。SEISCOPE 工具的普通用户不需要自己写 CG只需要提供 hess_vec_prod 的接口实现——在波形反演中这通常意味着把梯度计算函数对模型做一次小扰动重新算一遍梯度再除以扰动值这就是所谓的有限差分 Hessian-vector product。实际实现时扰动 δ 通常在 1e-6 到 1e-4 的量级视浮点精度而定。3. 导入 SEISCOPE 的 TRN 模块从压缩包到可执行的最小工程3.1 TRN.ZIP 解包后的核心文件定位TRN.ZIP 的目录结构在 SEISCOPE 项目中有固定的组织方式。工具箱的主控程序通常是 optimiz.f90根据用户选择的 method 参数分发到不同子模块。TRN 对应的文件是 pkg/optimization/trn_fwi.f90核心模块包括外部迭代控制outer loop、CG 内迭代求解器inner loop和 有限差分 Hessian-vector product 的实现compute_hessian_vect_prod。最小工程不需要全部读懂关键是找到三个接口函数的签名cost_function计算目标函数值和梯度、hessian_vect_prod计算 Hessian 与向量的乘积、precond预处理算子。# 解压并定位核心接口 (bash) unzip TRN.ZIP -d seiscope_trn cd seiscope_trn ls -la find . -name *trn* -o -name *optimiz* | head -20SEISCOPE 的 Fortran 源码在 Linux 环境下的编译依赖 gfortran 和 make进入源码根目录后make 命令会生成 libseiscope.a 静态库和若干可执行示例程序。编译前需要确认编译配置文件 Makefile.inc 里的编译器路径和并行选项。如果机器配备 MPI 环境建议开启并行编译因为 FWI 的正演部分在频率域多炮并行时开销比较大。但如果只是验证 TRN 的反演流程串行版就足够一个 2D Marmousi 模型的最小示例在单核上运行大约需要几十分钟。另外要留意 SEISCOPE 的数据格式约定观测数据 d_obs 和模拟数据 d_cal 都存储为 SEISCOPE binary format读入 Fortran 程序前需要把模型参数按列优先顺序展平成一维数组。第一步测试建议不要用真实地震数据而是用工具箱自带的示例数据通常在 data/ 目录下验证编译链路完整。3.2 通过 Python 封装调用 TRN 核心在真实工程项目中大部分团队会用 Python 做上层反演编排、数据清洗和可视化Fortran 优化器作为底层计算核心。把 TRN 模块封装成 Python 可调用的形式需要用到 ctypes 或者 f2py。SEISCOPE 源码里不直接提供 Python 绑定所以这个接口层需要自己写我一般会暴露 3 个底层函数给 Python# 用 ctypes 调 Fortran 编译出的共享库 (python) import ctypes import numpy as np lib ctypes.CDLL(./libseiscope_trn.so) # Fortran 子程序中, 所有参数按引用传递, 多维数组需按列优先展平 lib.trn_driver_.argtypes [ ctypes.POINTER(ctypes.c_int), # n: 模型参数个数 ctypes.POINTER(ctypes.c_double), # m: 模型向量 (in/out) ctypes.POINTER(ctypes.c_double), # g: 梯度向量 ctypes.c_void_p, # 用户自定义数据指针 ctypes.POINTER(ctypes.c_double), # 目标函数值 ] # 调用示例: 输入初始模型 m0, 返回反演结果 m_final m0 np.array([1500.0] * n_model, dtypenp.float64) fval np.array([0.0], dtypenp.float64) lib.trn_driver_(n_model, m0, grad, None, fval)这个封装有几个关键点Fortran 数组默认按列优先存储和 Python 默认的行优先不同因此任何多维数组传给 Fortran 前必须先调用 np.asfortranarray 强制转换函数名后面的下划线是 gfortran 对全局符号的修饰规则如果是 intel ifort 则可能没有下划线——这也是用动态库前需要先用 nm 命令检查符号名的原因。此外正演模拟部分如果本身是 Python 实现的比如用 Devito 或 WaveFD 生成正演数据则无法直接传给 Fortran 优化器做 Hessian-vector product只能把正演和伴随全部迁移到 Fortran 端这往往是把 SEISCOPE 集成到既有项目时最大的工作量所在。3.3 forward/adjoint 对的正确性验证在开始反演之前务必要验证梯度计算的正确性因为 TRN 内部对梯度精度的依赖比 L-BFGS 更强——梯度有误差CG 内迭代的残差计算立刻失真截断准则就会失效。常见做法是有限差分验证:对目标函数沿某个方向加一个微小扰动 ε比较解析梯度与数值梯度。# Taylor 检验: (python) 验证梯度与 hessian_vect_prod 的精度 def taylor_test(model, direction, eps_list[1e-2, 1e-4, 1e-6]): base_cost cost_function(model)[0] base_grad cost_function(model)[1] gd base_grad.dot(direction) for eps in eps_list: perturbed model eps * direction J_new cost_function(perturbed)[0] # 目标函数差值应随 eps 线性减少 print(feps{eps:.1e}, ratio{(J_new - base_cost)/eps:.6f}) # 期望输出: ratio 在 eps - 0 时趋近 gd如果 Taylor 检验的输出 ratio 和 gd 的差大于 1%说明梯度接口有 bug不应该继续反演。另一个是 HVP 的验证利用 H·v ≈ [g(mεv) - g(m)]/ε检查 hessian_vect_prod 输出与有限差分结果是否一致。这个验证 TRN 内迭代的每一步都会用到如果错了CG 得到的牛顿方向就不可靠。4. 让 TRN 跑起来的关键参数与典型坑4.1 SEISCOPE 优化工具箱的参数文件语义SEISCOPE 的参数输入文件遵循一段固定格式控制优化循环的参数包括 method、niter_max、memory、CG 相关参数和线搜索参数。TRN 模式与 L-BFGS 的差异集中在 CG 相关参数上L-BFGS 没有内迭代而 TRN 用内存换精度。参数名TRN 推荐范围作用调整代价method3 (TRN)选择优化器无niter_max50 - 20外部迭代轮数越大越费算力memory20 - 40预处理器 L-BFGS 的历史步数影响收敛速度cg_niter_max10 - 30CG 内迭代上限大则精度高但慢cg_eta0.1 - 0.5截断残差阈值小则内迭代次数剧增ls_max10 - 20线搜索最大步数太大会甩出稳定区设参数的首要原则是当你不确定某个值该设多少时让 cg_eta 偏大接近 0.5同时把 cg_niter_max 调到 10这样 TRN 实际运行成本接近 L-BFGS确认反演趋势正确后再逐步调高内迭代上限。这样先跑通再优化的策略在实际项目中能节省大量调试时间。SEISCOPE 还支持指定 preconditionerTRN 内部常用 L-BFGS 或对角 Hessian 近似作为预处理memory 参数在 TRN 中真实角色是预处理器保留多少梯度历史用于隐式近似 H 的逆这个值设太大会拖累每次预处理的线性代数求解。4.2 观察收敛日志的哪些字段SEISCOPE TRN 在日志中输出的字段与 L-BFGS 不同除了常规的 cost、g_norm、model_update_norm 外CG 残差历史CG_RESIDUAL是判断内迭代是否过多或过少的关键信号。如果 cg_residual 一路从 1.0 下降到 0.05然后平稳说明截断准则合理如果这个值在第一次或第二次 CG 迭代时已经低于截断阈值说明 cg_eta 太松TRN 退化成梯度下降如果 CG 迭代到上限仍然没有达到截断阈值说明 Hessian 病态严重这时候应该考虑增强预处理而非增加 cg_niter_max。日志里另一个值得关注的是 steplength。TRN 外迭代偶尔会出现 steplength 被线搜索压到极小值的情况。在 FWI 里这通常意味着 Hessian 里面有负曲率方向模型更新不满足下降条件。SEISCOPE 在这种状况下会重置 CG 初始搜索方向用户不需要干预但如果这种问题在前几轮就出现多半是初始模型太差或者观测数据里存在异常振幅。4.3 常见失败模式与排查路径TRN 在 FWI 中最常见的失败是波形反演发散表现为 cost 函数在前几轮下降后突然跃升随后梯度范数爆炸。排查路径按照以下步骤走先检查正演是否稳定用零延迟自相关对比模拟数据和观测数据主频是否一致接着检查梯度符号SEISCOPE 内部对权重和归一化的约定可能与你自己的代码不同梯度多个负号TRN 会一直往错误方向更新。最后才怀疑优化器本身。常见陷阱是 FWI 目标函数中数据残差的单位规范化方式。SEISCOPE 默认是用整体数据 L2 范数做归一化如果你自己写数据接口时又做了一次归一化Hessian-vector product 里的扰动缩放就会错位。# 检查一维模型梯度符号 (python) # 若 d_obs 是地震数据, model 是初始速度模型 # 正确做法: 分别计算 g grad(J), 其中 J(data) 2.0 * (F(m) - d_obs)另外TRN 内迭代的病态问题可以用数据预处理缓解常用的手段是对地震道做时间增益补偿和带通滤波让能量不集中在浅层。SEISCOPE 官方示例里对 Marmousi 模型做了子波估计、震源校正、静校正等步骤这些在反演前的处理可以极大改善 CG 收敛速度。5. 实测对比 L-BFGS 与 TRN 的收敛行为用 SEISCOPE 的 TRN 模块跑穿 2D Marmousi 模型之后最值得做的验证是把 L-BFGS 和 TRN 的结果放在同一坐标系下对比。这里的核心技巧是让两种方法停在同一梯度范数阈值下然后比较对应模型的数据拟合度与成像清晰度。我一般会写一个循环脚本分别调用 SEISCOPE 的 method1 (L-BFGS) 和 method3 (TRN)每 5 轮外迭代 dump 一次模型快照和目标函数。一个典型观测结果L-BFGS 在初始阶段下降速度比 TRN 快因为它的初始步长较大、每轮成本低。但到了中后段尤其在 3km/s 到 4.5km/s 深度区间存在高速层时L-BFGS 的收敛曲线进入明显平台期梯度范数在多个外迭代轮次陷入震荡TRN 则平稳穿越这个区间原因是 Hessian 中非对角项提供了界面处的曲率信息帮助模型跳出梯度不敏感区。这个差异在数据缺大偏移距时更加明显——大偏移距数据的缺失会使 Hessian 的弱照明方向特征值接近零L-BFGS 的近似在这些方向上过度放大噪声而 TRN 的 CG 内迭代会因残差截断机制自动抑制这些无效方向。需要特别指出的是TRN 并非在所有反演场景中都胜出。如果观测系统覆盖均匀、初始模型接近真实速度L-BFGS 因为每轮计算量更小达到同等数据拟合度所需的总时间往往更短。在对计算成本极其敏感的生产环境中一种常见的策略是把 L-BFGS 作为第一阶段优化跑 10 轮后切换到 TRN 做精细修正SEISCOPE 工具箱准许在运行中途切换优化器的标准做法是把当前模型与梯度状态导出再以新输入的初值继续跑另一种方法。如果你已经跑通上述流程、手头又有正演并行化能力我强烈建议做一次内迭代历史的离线分析把 CG 的残差序列 dump 到文件用 matplotlib 画成半对数坐标图。你能清楚看到截断时机和收敛拐点。当残差曲线几乎垂直向下一个数量级只花了 2-3 次 CG 迭代说明预处理效果好当残差曲线直线下降但到某个值后突然水平说明 Hessian 特征值分散严重CG 卡在了慢收敛模式。这种情况下把 cg_eta 从默认值调大到 0.4 会让 TRN 更早截断在损失少量精度的前提下避免无效内迭代。这个动作让实测总耗时可降低 30%-50%——这就是截断的意义。本文还有配套的精品资源点击获取