简介：一份面向地球物理全波形反演与优化算法研究者的示例工程，演示如何在SEISCOPE优化工具箱中调用截断牛顿（Truncated Newton）算法进行波形反演计算。程序实现对应Metivier等人2013年发表在SIAM Journal on Scientific Computing上的经典方法，适合需要复现该算法或开展反演实验的科研与工程人员参考。包内共22个文件，核心为15个Fortran 90源文件，涵盖算法主程序与子模块；另含两个数据文件、工程配置文件及Visual Studio相关文件，整体压缩包仅31KB，结构精简便于快速阅读与调试。该资源已有239人学习浏览，虽是小型示例，但完整呈现了截断牛顿法在全波形反演中的落地流程，有助于理解工具箱调用方式、数据组织与迭代求解细节，可作为二次开发或算法对比的起点。
1. 为什么全波形反演需要截断牛顿法，而不是多跑几轮 L-BFGS
FWI 的迭代优化器选择，往往比正演精度更容易被低估。很多人在拿到一个波动方程正演代码后，第一反应是用 L-BFGS 往上堆迭代次数，结果在强散射或深部目标上，梯度方向被几何扩散和多重散射拖慢，几十次迭代后成本函数仍在平台期。截断牛顿法（Truncated Newton）把一个完整的牛顿步拆成外迭代和内迭代：外层做线搜索、更新模型，内层用共轭梯度近似求解牛顿方程，全程不需要显式组装海森矩阵，只额外要求你提供一次海森向量积。这正是 SEISCOPE OPTIMIZATION TOOLBOX 里 TRN 示例的价值。如果你想在自己的 2D/3D 全波形反演流程里换一个收敛更扎实的优化器，这个压缩包里的代码和运行日志值得从头到尾拆一遍。适合的读者是已经有正演和伴随代码、但对优化器原理还停留在“梯度下降加个拟牛顿”的人。
2. TRN.ZIP 项目解剖：SEISCOPE 的算法骨架与两个迭代记录文件
2.1 压缩包里的工程骨架，先区分可编译与可忽略
TRN.ZIP不是传统意义上的规范源码包，更像开发者直接打包的工作目录。TRN.sln和TRN.vfproj是 Windows 上 Visual Studio + Intel Visual Fortran 的工程入口；src与trn_src保存核心源码；common放公共数据结构和常量；根目录的iterate_TRN.dat与iterate_TRN_CG.dat是上一次运行留下的外层与内层迭代日志。TRN.u2d、TRN.suo都是 IDE 生成的会话文件，不参与编译，可以忽略；tets从命名上看是 tests 的笔误，多半是作者做回归验证的脚本目录。这种布局的好处是你直接能看到程序的运行产物，坏处是构建顺序和依赖项必须自己判断。我拿到这类项目后第一件事是按“可编译、可运行、可忽略”三类重新归档，否则很容易在无关的调试文件上浪费时间。
2.2 optim_type.h：优化算法和用户代码之间的协议
optim_type.h是理解 SEISCOPE 优化框架的钥匙，它把优化算法和具体反演问题解耦：算法层只认结构体，不关心你的正演是用时域有限差分还是频率域 Born 近似。简化开来看，结构体长这样：
/* 从 SEISCOPE 工具箱接口风格简化的定义 */ typedef struct optim_type { int verbose; /* 是否打印每轮日志 */ int nouter_max; /* 外层截断牛顿最大迭代次数 */ int ncg_max; /* 内层 CG 最大迭代次数 */ double cg_tol; /* 内层 CG 相对残差容差 */ double c1, c2; /* Wolfe 线搜索参数 */ void (*cost)(const double *m, double *fcost, void *ctx); void (*grad)(const double *m, double *g, void *ctx); void (*hess_vec)(const double *m, const double *v, double *Hv, void *ctx); } optim_type;在 SEISCOPE 实际版本里，字段名和组织方式不一定完全一样，但你要准备的三样东西不会变：计算目标函数、计算梯度、给一个向量 v 返回 Hv。这个设计是截断牛顿法能落地的关键。传统牛顿法需要显式组装 N×N 的海森矩阵，FWI 里 N 是模型参数数量，2D 网格到 10^5 已经令密矩阵存储压力很大，3D 更不可能；而用二阶伴随或散射正演，只需要一到两次等效正演就能得到 Hv。TRN 外层通过线搜索保证下降，内层通过 CG 控制牛顿方向精度，两者都只依赖 Hv，所以这套接口能撑起大参数规模的反演。
2.3 外层 iterate_TRN.dat 与内层 iterate_TRN_CG.dat 分别记录什么
这两个文件不是下载占位符，而是算法主动写出的诊断数据，比终端输出可靠。iterate_TRN.dat记录外层截断牛顿迭代，典型列包括：迭代号、目标函数值、梯度 L2 范数、线搜索步长、本轮 CG 次数。iterate_TRN_CG.dat记录某轮外迭代中内层 CG 的每次迭代信息，常见列是 CG 迭代号、残差范数或方向修正量。我一般先看这两个文件，不看终端，因为终端会滚动且通常只打印摘要。快速验证是否正常可以用：
# 查看外层目标函数与梯度范数的变化趋势 head -5 iterate_TRN.dat tail -5 iterate_TRN.dat # 查看内层 CG 记录有多少条，是否有异常长迭代 wc -l iterate_TRN_CG.dat这段命令里head/tail是 Linux/macOS 命令，Windows 上对应 PowerShell 的Get-Content -Head。目标函数应大致单调下降，梯度范数的数量级应从初始值下降几个量级；wc -l统计行数用于判断内层 CG 是否频繁突破预期迭代上限。
2.4 从 Debug 目录和 tets 目录反推运行习惯
Debug目录说明作者使用了 Visual Studio 默认的 Debug 配置构建，根目录出现iterate_TRN_CG.dat则说明程序的工作目录就是工程根目录，输入输出都用相对路径。复现时尽量不要修改这个约束，否则程序很容易找不到模型文件。tets目录提醒我们：TRN 这类与 CG 残差、Wolfe 条件纠缠的算法，必须有几个可重复的小模型做回归用例。你改一行容差，可能让一个在 A 模型上收敛良好的测试在 B 模型上悄悄退化，没有回归用例很难意识到。
2.5 TRN 与 L-BFGS 的本质差异：从近似逆到近似解
L-BFGS 的本质是用过去的梯度差和模型差去逼近海森矩阵的逆，好处是每次迭代只需要额外的向量运算，坏处是逼近质量依赖历史迭代的多样性和问题非线性程度。TRN 则不同，它在每个外层步用 CG 去直接求解牛顿方程，CG 的每次迭代都调用精确的 Hv 算子，因此只要问题本身海森是正定的，它获得的搜索方向比 L-BFGS 的“替换品”更接近真实牛顿方向。这也是为什么在强散射介质、长偏移距数据下，TRN 能把 FWI 从几十次迭代压到十几次，而 L-BFGS 常常卡在梯度平台区。代价是每次外层迭代的正演次数显著增加，所以它不是替代 L-BFGS，而是在预算允许时更高阶的选择。
3. 从 TRN.sln 到反演结果：Windows 构建、跨平台编译与迭代日志判读
3.1 在 Windows 上把 .sln 跑通的关键不是 VS，而是 Intel Fortran
TRN.vfproj是 Intel Visual Fortran 工程文件，打开TRN.sln后如果项目显示为空，通常是因为本机没有安装 Intel Visual Fortran 的 VS 集成。安装 Intel oneAPI 时选择 Fortran 编译器与 Visual Studio 集成组件即可。构建时 Debug/Release 的选择也影响性能：Release 下编译优化不会被调试符号拖慢，FWI 正演热区在 Debug 下可能慢一个量级。若出现无法解析的外部符号，先检查项目属性的库目录是否指向 MKL 的对应架构路径。命令行构建可以这样：
call "%ONEAPI_ROOT%\compiler\latest\env\vars.bat" intel64 msbuild TRN.sln /p:Configuration=Release /p:Platform=x64%ONEAPI_ROOT%是你的 oneAPI 安装根目录；/p:Configuration选择构建配置；/p:Platform要与本机架构及 MKL 库一致。注意vars.bat后需要加intel64或ia32，这取决于终端架构，如果混用，编译链接阶段往往会出现环境变量冲突。
3.2 跨平台手动编译兜底方案
如果只有 Linux 或者不想依赖 Visual Studio，Fortran 源码本身可以跨平台，只差一份构建规则。常见做法是把源码列表拉出来，然后直接编译：
find common src trn_src \( -name "*.f90" -o -name "*.F90" -o -name "*.f" -o -name "*.F" \) > sources.txt ifort -O3 -I common -I src -mkl -o trn_run $(cat sources.txt)find命令收集常见 Fortran 扩展名文件到sources.txt，括号需要转义防止被 shell 解释；-I common -I src把模块目录和公共头文件目录交给编译器；-mkl是 Intel 编译器连接 MKL 的快捷开关，gfortran 下不存在，需要改成-lblas -llapack；$(cat sources.txt)把文件列表展开成多个源文件参数。如果编译报错说找不到某个模块，多半是-I路径顺序不对，Fortran 模块依赖要求编译顺序按照 use 关系从底层往上排。最稳妥的做法是让编译器自动加载.mod文件，而不是手动调源文件顺序。
3.3 主程序流程：TRN 掉进优化器之前要准备好什么
从src/trn_src的代码组织看，主程序并不是把优化器函数一股脑塞进去，而是先准备好四个输入：初始模型、观测数据、目标函数接口、梯度接口。简化流程如下：
! 主程序流程示意 call read_model(m0) ! 初始速度/参数模型 call load_observed_data(dobs) ! 观测波形数据 call misfit_initial(m0, fcost, grad) ! 初始误差和伴随梯度 optim_type%cost => cost_wrapper optim_type%grad => grad_wrapper optim_type%hess_vec => hess_vec_wrapper call TRN_run(m0, fcost, grad, optim_type)这里cost_wrapper、grad_wrapper和hess_vec_wrapper是你自己的 Fortran 函数，TRN 求解器通过函数指针回调。注意hess_vec_wrapper接收的扰动向量 v 来自 CG 迭代，内容没有任何物理含义，可能包含棋盘格噪声；如果你的正演代码对 v 做了平滑或截断，就会返回一个不一致的 Hv 给 CG。这是很多 TRN 实现跑起来比理论慢的原因——Hv 没按“未处理的原向量”计算。
3.4 跑完以后先看数据，而不是看终端
程序退出码正常不代表反演正确。我通常立刻对比两个迭代文件：
# 外层是否单调下降 awk '{print $1, $2, $3}' iterate_TRN.dat | head -30 # 内层 CG 残差是否降到设定容差以下 awk '$NF < 1e-6 {print NR, $NF}' iterate_TRN_CG.dat | head第一条awk命令取前三列，分别对应迭代号、目标函数值、梯度范数；如果第三列的数量级在下降，说明外层搜索方向是有效的。第二条命令把内层 CG 记录中最后一列小于 1e-6 的行打印出来，用来确认内层迭代确实在收敛，而不是被ncg_max硬截断。如果打印为空，说明 CG 基本没达到默认容差，需要去检查预条件和 Hessian 向量积。
注意：不同版本的程序写出文件的列数和排列可能不同，先用
head -1看表头，再决定 awk 的列号，不要照搬这里的列索引。
4. 调参截断牛顿：线搜索、CG 容差与 Hv 的一致性
4.1 外层线搜索：Wolfe 条件怎么给，才算不过分保守
TRN 的线搜索目标不是精确找到极小点，而是找到一个满足充分下降和曲率条件的步长，即 Wolfe 条件。SEISCOPE 示例里的默认参数通常取 c1=1e-4，c2=0.7 到 0.9。c2 太小会让线搜索把每一步试满，成本飙升；c2 太大则容易接受过大的步长，导致下一轮 CG 起点变差。判断线搜索状态的指标是目标函数曲线中是否出现“连续小步长密集下降”。如果步长序列大量落在 0.01 以下，但目标函数还在缓慢下降，这通常不是 c1、c2 需要调，而是梯度方向包含严重的数值噪声；这时要检查正演是否用了足够的模板精度或吸收边界，而不是继续压线搜索参数。
4.2 内层 CG 容差：跟着外层梯度走，不要设死值
截断牛顿的“截断”就是指内层 CG 没有真正解到机器精度。在 L. Metivier、R. Brossier、J. Virieux、S. Operto 于 2013 年发表在 SIAM Journal on Scientific Computing 上的论文中，核心讨论之一就是如何控制内层 CG 的残差。常见做法是采用 Eisenstat-Walker 策略：初始给一个较大容差，随着外层梯度范数减小逐步收紧。用伪代码表达就是：
double eta = 0.5; // 初始容差 if (outer_iter > 0) { eta = fmin(0.5, sqrt(grad_norm / grad_norm_initial)); eta = fmax(eta, 0.1); // 下限保护 } cg_solve(..., eta);这里的eta是内层 CG 的相对残差目标。开始迭代时，模型距离最优解很远，方向稍微粗糙点也能接受，所以容差可以放到 0.1~0.3；等到外层梯度缩小到初始的百分之一，如果还保持 0.3 的容差，CG 给出的方向可能离牛顿方向偏移很大，所以要让容差按梯度下降比例收缩。fmax下限 0.1 是为了避免在数值噪声主导时把 CG 逼到无意义的精细求解。我一般在这个框架基础上保留ncg_max作为硬上限，因为即使容差很小，正演算子的舍入误差也会限制 CG 的实际可达精度。
4.3 预条件与模型维度：TRN 不是免预条件的
很多用户把 TRN 直接套用到时间域 FWI，完全不加预条件，然后发现内层 CG 永远在极限迭代，算法表现比 L-BFGS 还差。问题不在截断牛顿，而在于海森条件数太差。常见做法是给内层 CG 加一个对角预条件子，用模型空间振幅归一化或平滑算子近似海森对角项。要注意，模型网格超过 10^7 量级时，每次 Hv 至少需要一次正演一次伴随，3D 时间域代价仍然很高。因此 TRN 最适合的区间是“一次正演成本适中、但梯度方向需要更准确曲率”的中间规模问题。在超大网格上，我宁可把 CG 内迭代次数压到 5 次以下，让 TRN 退化成一种含自适应正则的梯度法，也比硬上完整牛顿方向划算。
4.4 常见失败模式与排查顺序
| 现象 | 可能原因 | 优先排查 |
|---|---|---|
| 外层 cost 上升且步长趋近 0 | 梯度或 Hv 不一致 | 梯度有限差分校验 |
| 内层 CG 残差长期不变 | Hv 实现与梯度不同源 | 单独测试 Hv |
| cost 下降但梯度范数不降 | 正则项压掉了梯度信息 | 检查目标函数归一化 |
| 单轮时间几乎全在 CG | cg_tol 过严或没有预条件 | 放宽容差、加对角预条件 |
| 换网格尺寸后收敛性突变 | 数据残差未按网格体积缩放 | 检查观测误差归一化 |
这张表是我在实际 FWI 中总结出来的排查顺序。第二行要特别强调：Hv 必须和梯度使用同一套正演与伴随实现，任何网格差异、边界条件差异或近似差分阶数不一致，都会让 CG 残差像噪声一样停在某个水平线上。这类问题不会在梯度检验中发现，需要手动构造一个小扰动方向来自测。
注意：调参顺序应当是先验证梯度与 Hv、再调线搜索参数、最后动 CG 容差；反过来会浪费大量计算。
5. 在自定义 FWI 代码里复用 TRN 思想：梯度与 Hv 的一致性命门
5.1 梯度有限差分模板
在把 TRN 接入自己的 FWI 流程前，先做一次梯度验证。中心差分是常规做法：
def check_gradient(m, misfit, grad, eps=1e-6): mplus = m.copy(); mplus[0] += eps mminus = m.copy(); mminus[0] -= eps num_g = (misfit(mplus) - misfit(mminus)) / (2 * eps) rel_err = abs(num_g - grad[0]) / (abs(num_g) + 1e-30) print(f"rel_err={rel_err:.3e}")这个函数只检验第一个参数，实际应当随机抽一批分量，并保证扰动不超过正演离散误差允许的范围。eps选太大会引入截断误差，选太小则浮点对消开始支配结果。我一般先在 1e-7、1e-6、1e-5 三档测试，选取相对误差稳定且不随eps剧烈变化的一档作为正式测试配置。
5.2 Hv 的 Taylor 展开检验
梯度验证通过后，继续验证 Hessian 向量积。思路是对扰动向量 v，定义 φ(t)=g(m+t·v)，那么 φ'(0)=Hv，可以用差分近似侧的数值导数来对比你的hess_vec。示例代码如下：
import numpy as np t = 1e-6 g_plus = compute_gradient(m + t * v) g_minus = compute_gradient(m - t * v) num_Hv = (g_plus - g_minus) / (2 * t) ana_Hv = hess_vec(m, v) rel = np.linalg.norm(num_Hv - ana_Hv) / np.linalg.norm(ana_Hv) print(rel)这里的v必须取随机方向，不能取网格基向量；随机方向的覆盖更广泛，能抓出类似“只有某些分量有符号错误”的 bug。另一个细节是，如果正演代码中有边界吸收条件，扰动向量如果碰到边界区域，数值导数会出现来源不明的异常，那不代表 Hv 错，而是边界类实现没有对扰动向量做相同处理。测试时把扰动限制在计算域内部，缩小到远处。
5.3 量纲与归一化：正式接入时最值得注意的一件事
当梯度与 Hv 检验都通过后，TRN 才能真正被信任。实际接入时最容易引发问题的不是优化器内部，而是模型物理量纲。以速度模型为例，单位取 km/s 时，速度值在 1.5 到 6.0 之间，数据残差如果是地震振幅，数量级可能在 1e-3 到 1e2 之间。两者相乘产生的梯度和 Hessian 条件数，会让内层 CG 的残差在真正收敛之前就触碰浮点极限。建议在主循环开始前把模型向量的最大值归一化到 1，反演迭代完成后把增量乘回原始量纲，这一步往往比任何 cg_tol 和 c2 都更能影响收敛速度。
本文还有配套的精品资源，点击获取
