2026/9/13 15:23:27
FEATURED · 精选文章 
简介面向地震全波形反演与优化算法研究者资源给出了基于SEISCOPE优化工具箱的截断牛顿算法示例程序对应Metivier等人2013年发表于SIAM Journal on Scientific Computing的经典方法。资源共22个文件主体为15个Fortran 90源文件实现核心优化迭代逻辑另含2个数据文件用于输入输出以及Visual Studio工程与解决方案文件便于在Windows环境下编译调试整体仅31KB轻量易读。已有239人学习适合正在学习全波形反演、拟牛顿或截断牛顿法的高年级研究生与科研人员。通过这份代码读者可对照论文理清截断牛顿算法在海量地震数据反演中的具体实现步骤包括迭代求解、梯度与Hessian向量积计算等关键环节作为自己工具箱开发或实验对比的参考起点。1. TRN.ZIP 里的截断牛顿优化解决的是波形反演哪一环地震波全波形反演FWI本质上是一个带多频率正演模拟的大规模非线性最小二乘问题。模型参数动辄几十万到上千万目标函数在速度—密度空间里既不凸又存在周期性相位差带来的局部极小值。工程上大多数实现停留在梯度类方法或者 Gauss-Newton 变体上但一旦介质出现强散射、速度对比超过 20%梯度类方法的收敛半径就明显不够用。截断牛顿优化Truncated Newton简称 TRN通过在海森矩阵方向上进行截断共轭梯度求解在每一步外迭代里近似得到一个牛顿方向既能保留牛顿法对二阶曲率的利用又不必显式组装海森矩阵这让它成为 seiscope 这类教学与科研用反演框架里“往深处走一步”的常见选择。适合已经跑通过一次最简 FWI、被局部极值和慢收敛卡住想弄明白优化器内部到底在干什么的那类使用者。2. 从梯度到牛顿截断牛顿优化在波形反演里的原理2.1 梯度方向为什么不够用Gauss-Newton 又差在哪儿FWI 目标函数写成常见形式J(m) 1/2 · Σ_f || P·u_f(m) - d_f ||²其中m是模型向量u_f是频率f时的正演波场P是检波器采样算子d_f是观测数据。普通梯度下降只使用一阶信息更新公式是m_{k1} m_k - α_k·∇J(m_k)。梯度方向在远离极小值时下降较快但进入波场干涉严重的区域后不同参数维度之间的耦合让梯度更新产生大量之字形震荡。Gauss-Newton 方法把牛顿方向中的海森矩阵近似为Re(J^H J)其中J是 Jacobian 矩阵。这个近似忽略了数据残差对模型参数的二阶导数项优点是正定且计算代价可控缺点是在强散射介质中残差项本身携带重要的二阶信息忽略之后容易高估更新步长导致迭代发散。完整牛顿法使用真实海森H_k · p_k -∇J(m_k)但海森矩阵的规模是模型参数量的平方对三维 FWI 完全不现实。把这两个事实放在一起结论就很清楚需要一个不用显式构造海森、又比 Gauss-Newton 更接近完整牛顿方向的方法。这就是 TRN 出现的位置。2.2 “截断”指的是内层共轭梯度迭代被提前终止TRN 的核心思路是把海森系统H·p -g的求解交给共轭梯度CG法而 CG 只迭代有限步不等到完全收敛。这个提前终止动作就是“截断”。外层是通常的非线性迭代框架for k 0, 1, 2, ...: 计算梯度 g_k 用 CG 迭代求解 H_k·p_k -g_kCG 次数 N_cg(k) 由收敛判据决定 沿 p_k 做线搜索得到步长 α_k m_{k1} m_k α_k·p_kCG 迭代终止条件一般写成||r_j|| min(η_k, ||g_k||^θ) · ||g_k||r_j是当前 CG 残差η_k是用户设定的容差θ通常取 0.9 左右。η_k设在 0.1 到 0.5 之间既保留牛顿方向质量又避免在数据噪声很大时去精确求解一个本身就不准的线性系统。截断带来的副产物是每次外迭代的计算量不再是固定值而是随梯度范数和模型复杂度动态变化。这个动态性在实测数据上看起来非常直观前几次外迭代 CG 次数多越往后越少。2.3 一次 Hessian-向量积需要四次波场模拟CG 迭代内部只需要海森矩阵与向量的乘积H·v这是 TRN 能用于大规模问题的关键。在频域波形反演中H·v可以通过波场积分算子展开而不用显式生成海森矩阵。假设正演算子为L(m)u s伴随波场为λ那么H·v的常见计算路径是u solve(L, src) # 正向场 λ solve(L^*, P^H (P u - d)) # 伴随场 δu -solve(L, (∂L/∂m·v)·u) # 扰动正向场 δλ solve(L^*, P^H P δu - (∂L/∂m·v)^H λ) # 扰动伴随场 H·v Re[ ((∂L/∂m)^H δλ (∂L/∂m·δu)^H λ) ]每频率需要两次正向求解和两次伴随求解通常笼统叫“四次正演”。这是 TRN 比 L-BFGS 贵的主要原因但相比显式海森已经是数量级的节省。实现时注意∂L/∂m·v这一项表示模型扰动对波场方程的影响在有限差分离散中等于对差分系数的线性扰动不需要额外存储 Jacobian。seiscope 框架里正演和伴随算子已经有了成熟接口TRN 部分只需要把上述流程接进目标函数和梯度计算之后。3. 把截断牛顿优化接进 seiscope 波形反演主循环3.1 从 TRN.ZIP 到可运行的工程目录拿到TRN.ZIP后第一件事不是读文档而是先确认包的结构。常见做法是解压后目录里包含一个优化器核心、一个线搜索模块、以及若干示例模型脚本。可以先用下面命令确认文件布局unzip TRN.ZIP -d trn_work cd trn_work tree -L 2输出里如果看到类似trn_solver、line_search、example/这样的目录说明结构是统一的。我一般会把源码放到src/把模型声明和观测数据放到单独目录避免运行示例时把生成文件写进源码树。无论包怎么组织最终反演代码至少要暴露三个接口forward_solver正演、gradient梯度、hessian_vector_productHv 积。如果包没直接提供 Hv需要自己补后面接 seiscope 的建模器时也用得上。运行环境方面seiscope 系代码通常依赖 NumPy 和 SciPy某些示例还需要matplotlib来做结果可视化。不需要 GPU在二维小模型上 CPU 已经能在一小时内完成一个几十次外迭代的 TRN 实验。3.2 反演主循环的最小 Python 骨架下面是一个不依赖具体 seiscope API 的 TRN 外层循环骨架重点在于表达控制流把优化器与正演解耦。实际使用时应把forward_op替换成 seiscope 的正演封装。import numpy as np class TruncatedNewtonFWI: def __init__(self, forward_op, adjoint_op, model0, params): self.fwd forward_op # 正演 self.adj adjoint_op # 伴随 self.m model0.copy() self.eta params.get(eta, 0.2) # CG 截断容差 self.cg_max params.get(cg_max, 30) # 内层最大迭代次数 def compute_gradient(self, freq): u self.fwd(self.m, freq) r u - self.data[freq] lam self.adj(self.m, freq, r) g self._apply_dLdm_adj(self.m, u, lam, freq) return g.real def compute_Hv(self, v, freq): u self.fwd(self.m, freq) du self._solve_perturbation(self.m, u, v, freq) lam self.adj(self.m, freq, self.data[freq] - u) dlam self._solve_adj_perturbation(self.m, u, du, lam, v, freq) Hv ( self._apply_dLdm_adj(self.m, u, dlam, freq) self._apply_dLdm_adj(self.m, du, lam, freq) ) return Hv.real def outer_step(self, freq): g self.compute_gradient(freq) p self._truncated_cg(g, freq) # 内层 CG 求解 H p -g alpha self._line_search(g, p, freq) self.m alpha * p return np.linalg.norm(g)逻辑说明compute_gradient先正演得到波场残差再走一次伴随得到梯度compute_Hv实现的是上一章的四波场流程。_truncated_cg内部运行共轭梯度每轮只需要调用compute_Hv不需要知道海森矩阵的显式形式。_line_search采用回溯线搜索目标函数下降就接受步长下降过少就减半。参数说明里最重要的三个量eta越小内层 CG 解越精确但单次外迭代成本越高cg_max是内层 CG 次数上限防止在大模型上无限制迭代freq是当前参与反演的频率频率延续策略下每次外迭代可以换一个频率。3.3 和 seiscope 模型对接先做梯度校验再谈收敛接 seiscope 时最常见的错误是梯度方向错了但是反演仍然在前几次迭代里“看着在下降”后边突然发散。为此正式跑 TRN 前一定要做一次有限差分梯度校验。方法是对模型m加入小扰动δ·v比较解析梯度与差分梯度的点积def check_gradient(fwi, m, v, eps1e-6, freq5.0): g_analytic fwi.compute_gradient(freq) f_plus fwi.objective(m eps * v, freq) f_minus fwi.objective(m - eps * v, freq) g_fd (f_plus - f_minus) / (2 * eps) ratio np.dot(g_analytic, v) / g_fd print(ratio:, ratio) # 接近 1.0 说明梯度实现正确eps取1e-6到1e-7太大则差分误差主导太小则浮点噪声主导。ratio在0.9 ~ 1.1之间视为通过。seiscope 的建模器可能需要额外的吸收边界项参与梯度这会导致解析梯度和差分梯度系统性偏离所以要先把边界条件固定住再测。这一步通过之后TRN 的外层方向计算才有意义。4. 影响截断牛顿优化的参数与收敛控制4.1 内层迭代容差、线搜索与正则化的三件套TRN 实际调参时先看一眼收敛曲线是哪种形态。下表是参数影响速查表按重要性排列参数常见范围调小的影响调大的影响内层 CG 容差eta0.1 ~ 0.5方向更接近牛顿方向单步成本高方向质量下降近似退化为梯度法内层 CG 最大迭代数10 ~ 50每个频率计算快收敛慢计算量线性增加可能过度拟合海森近似误差回溯线搜索衰减因子0.5 ~ 0.8步长收缩快迭代次数变多步长收缩慢容易越过下降区间模型正则化权重1e-4 ~ 1e-2反演结果更容易出现高波数噪声模型过于平滑低波数信息被压制eta的设定有一个局部最优区间。把它设成 0.05 时内外层迭代路径看起来非常干净但每次外迭代的波场模拟次数可能翻倍设成 0.6 以上时截断后的方向里还残留明显的负曲率分量线搜索被迫频繁缩短步长。实际数据上更推荐在反演前期用较大的eta0.3 左右快速进入目标函数主谷后期把eta调小到 0.1 收敛局部细节。这个策略和频率延续的节奏配合起来非常顺手。回溯线搜索里还有一层容易被忽略的参数就是初始步长。TRN 的搜索方向尺度受海森近似影响并不是单位化的所以初始步长建议用α0 -g·p / (p·H·p)来估计这个值可以直接从 CG 最后一步的结果中取不需要额外正演。4.2 频率延续策略对 TRN 收敛的帮助由于 FWI 目标函数非凸TRN 和梯度法一样需要低频先行。常见做法是从 3~5 Hz 开始每个频率做若干次外迭代后向高频推进。seiscope 的示例里通常把频率组分拆成数组TRN 外循环在相邻频率间采用上一频率的解作为新频率的初始模型。这个 Warm Start 对 TRN 尤其重要因为牛顿方向容易把模型推入一个对高频数据拟合更好、但低频信息已经丢失的局部极小值里。低频阶段建议低频迭代每频率只做 5~10 次外迭代。过多反而会让模型被低频数据里的振幅误差带偏。高频阶段再把这组频率的eta调小让模型精细尺度被逐步锁住。如果内存允许可以把多个频率同时放进目标函数里求和梯度这会让 Hv 积的波场模拟次数按频率数倍增所以更常见的做法是每次只激活一个频率。4.3 三个常见坑内层跑满、阻尼不当、边界污染第一个坑是把cg_max设得很大比如 100却不提高eta的精度要求。结果内层 CG 在后半段其实是在追海森近似本身的噪声尤其在使用有限差分正演且数据含噪时CG 残差会停滞不降白白烧算力。解决方法是让cg_max与eta联动eta0.2时cg_max设 20~30 就够。第二个坑是正则化权重设上头。TRN 对海森矩阵的正定化依赖比梯度法更强因为牛顿方向要求 H 正定才能得到下降方向当正则化过小时 CG 可能遇到不收敛或负曲率方向过大时反演分辨率急剧下降。建议用自适应正则化每次外迭代检查 CG 是否在少量迭代内产生非下降搜索方向若是则把正则化权重乘上 1.5否则逐步恢复。第三个坑是吸收边界条件在 Hv 积里被忽略。全波形反演的梯度公式假定正演算子完美满足边界条件但实际有限差分实现里边界系数同样依赖于模型参数如果不把边界对模型参数的导数包含进∂L/∂m·v最终的近似海森方向会在近边界区域产生五到十道快速变化的噪声条纹。排查方法是把目标梯度集中到模型边界附近的三个网格点做有限差分校验如果比值在边界处明显偏离 1.0就说明边界导数项被遗漏了。5. 实战验证用合成数据判断 TRN 是否真正收敛5.1 设计一个能区分 TRN 与梯度法的小模型用一个 64×64 网格的二维声波走时模型做验证。背景速度 2500 m/s在中心位置放置一个半径 6 个网格、速度 2800 m/s 的圆形异常体观测来自 8 个炮点、每炮 32 个检波器。低频从 2 Hz 开始这个设置在几百次迭代内能产生肉眼可辨认的收敛差异。运行 TRN 时建议每隔几次外迭代把模型和收敛曲线同时画出来。正常收敛的曲线形态是前三次外迭代目标函数快速下降随后进入一个平缓台阶。对比使用 L-BFGS 的结果TRN 在台阶处能继续向下走出一个小幅但明确的下降这正是海森方向里的二阶信息在低速异常体边界处起作用的表现。5.2 从三张图判断问题出在优化器还是正演器第一张图是目标函数曲线。如果曲线出现“锯齿形”的反复震荡优先怀疑线搜索的初始步长估算不对回到上一章α0 -g·p / (p·H·p)的做法重新计算。第二张图是梯度范数曲线。梯度范数平稳下降是最好的信号如果梯度范数下降但目标函数不降说明搜索方向里有显著的非线性残差项需要把eta调小如果梯度范数先降后飙升这通常是 CG 迭代超限内层求解出来的p根本不是可信下降方向。第三张图是模型剖面沿垂直方向的切片。TRN 在正确收敛时会把速度异常体的边缘锐化切片的过渡带只有两到三个网格而 L-BFGS 得到的过渡带往往更宽且含有旁瓣振荡。看到这种差异才说明你真正把截断牛顿优化的能力用上了。最后运行一次check_gradient并在脚本开头打印ratio如果该值不是接近 1.0优化器本身调得再好也无济于事问题一定出在建模器到优化器的接口这一层。本文还有配套的精品资源点击获取
截断牛顿优化：从原理到工程实践的FWI应用指南
简介面向地震全波形反演与优化算法研究者资源给出了基于SEISCOPE优化工具箱的截断牛顿算法示例程序对应Metivier等人2013年发表于SIAM Journal on Scientific Computing的经典方法。资源共22个文件主体为15个Fortran 90源文件实现核心优化迭代逻辑另含2个数据文件用于输入输出以及Visual Studio工程与解决方案文件便于在Windows环境下编译调试整体仅31KB轻量易读。已有239人学习适合正在学习全波形反演、拟牛顿或截断牛顿法的高年级研究生与科研人员。通过这份代码读者可对照论文理清截断牛顿算法在海量地震数据反演中的具体实现步骤包括迭代求解、梯度与Hessian向量积计算等关键环节作为自己工具箱开发或实验对比的参考起点。1. TRN.ZIP 里的截断牛顿优化解决的是波形反演哪一环地震波全波形反演FWI本质上是一个带多频率正演模拟的大规模非线性最小二乘问题。模型参数动辄几十万到上千万目标函数在速度—密度空间里既不凸又存在周期性相位差带来的局部极小值。工程上大多数实现停留在梯度类方法或者 Gauss-Newton 变体上但一旦介质出现强散射、速度对比超过 20%梯度类方法的收敛半径就明显不够用。截断牛顿优化Truncated Newton简称 TRN通过在海森矩阵方向上进行截断共轭梯度求解在每一步外迭代里近似得到一个牛顿方向既能保留牛顿法对二阶曲率的利用又不必显式组装海森矩阵这让它成为 seiscope 这类教学与科研用反演框架里“往深处走一步”的常见选择。适合已经跑通过一次最简 FWI、被局部极值和慢收敛卡住想弄明白优化器内部到底在干什么的那类使用者。2. 从梯度到牛顿截断牛顿优化在波形反演里的原理2.1 梯度方向为什么不够用Gauss-Newton 又差在哪儿FWI 目标函数写成常见形式J(m) 1/2 · Σ_f || P·u_f(m) - d_f ||²其中m是模型向量u_f是频率f时的正演波场P是检波器采样算子d_f是观测数据。普通梯度下降只使用一阶信息更新公式是m_{k1} m_k - α_k·∇J(m_k)。梯度方向在远离极小值时下降较快但进入波场干涉严重的区域后不同参数维度之间的耦合让梯度更新产生大量之字形震荡。Gauss-Newton 方法把牛顿方向中的海森矩阵近似为Re(J^H J)其中J是 Jacobian 矩阵。这个近似忽略了数据残差对模型参数的二阶导数项优点是正定且计算代价可控缺点是在强散射介质中残差项本身携带重要的二阶信息忽略之后容易高估更新步长导致迭代发散。完整牛顿法使用真实海森H_k · p_k -∇J(m_k)但海森矩阵的规模是模型参数量的平方对三维 FWI 完全不现实。把这两个事实放在一起结论就很清楚需要一个不用显式构造海森、又比 Gauss-Newton 更接近完整牛顿方向的方法。这就是 TRN 出现的位置。2.2 “截断”指的是内层共轭梯度迭代被提前终止TRN 的核心思路是把海森系统H·p -g的求解交给共轭梯度CG法而 CG 只迭代有限步不等到完全收敛。这个提前终止动作就是“截断”。外层是通常的非线性迭代框架for k 0, 1, 2, ...: 计算梯度 g_k 用 CG 迭代求解 H_k·p_k -g_kCG 次数 N_cg(k) 由收敛判据决定 沿 p_k 做线搜索得到步长 α_k m_{k1} m_k α_k·p_kCG 迭代终止条件一般写成||r_j|| min(η_k, ||g_k||^θ) · ||g_k||r_j是当前 CG 残差η_k是用户设定的容差θ通常取 0.9 左右。η_k设在 0.1 到 0.5 之间既保留牛顿方向质量又避免在数据噪声很大时去精确求解一个本身就不准的线性系统。截断带来的副产物是每次外迭代的计算量不再是固定值而是随梯度范数和模型复杂度动态变化。这个动态性在实测数据上看起来非常直观前几次外迭代 CG 次数多越往后越少。2.3 一次 Hessian-向量积需要四次波场模拟CG 迭代内部只需要海森矩阵与向量的乘积H·v这是 TRN 能用于大规模问题的关键。在频域波形反演中H·v可以通过波场积分算子展开而不用显式生成海森矩阵。假设正演算子为L(m)u s伴随波场为λ那么H·v的常见计算路径是u solve(L, src) # 正向场 λ solve(L^*, P^H (P u - d)) # 伴随场 δu -solve(L, (∂L/∂m·v)·u) # 扰动正向场 δλ solve(L^*, P^H P δu - (∂L/∂m·v)^H λ) # 扰动伴随场 H·v Re[ ((∂L/∂m)^H δλ (∂L/∂m·δu)^H λ) ]每频率需要两次正向求解和两次伴随求解通常笼统叫“四次正演”。这是 TRN 比 L-BFGS 贵的主要原因但相比显式海森已经是数量级的节省。实现时注意∂L/∂m·v这一项表示模型扰动对波场方程的影响在有限差分离散中等于对差分系数的线性扰动不需要额外存储 Jacobian。seiscope 框架里正演和伴随算子已经有了成熟接口TRN 部分只需要把上述流程接进目标函数和梯度计算之后。3. 把截断牛顿优化接进 seiscope 波形反演主循环3.1 从 TRN.ZIP 到可运行的工程目录拿到TRN.ZIP后第一件事不是读文档而是先确认包的结构。常见做法是解压后目录里包含一个优化器核心、一个线搜索模块、以及若干示例模型脚本。可以先用下面命令确认文件布局unzip TRN.ZIP -d trn_work cd trn_work tree -L 2输出里如果看到类似trn_solver、line_search、example/这样的目录说明结构是统一的。我一般会把源码放到src/把模型声明和观测数据放到单独目录避免运行示例时把生成文件写进源码树。无论包怎么组织最终反演代码至少要暴露三个接口forward_solver正演、gradient梯度、hessian_vector_productHv 积。如果包没直接提供 Hv需要自己补后面接 seiscope 的建模器时也用得上。运行环境方面seiscope 系代码通常依赖 NumPy 和 SciPy某些示例还需要matplotlib来做结果可视化。不需要 GPU在二维小模型上 CPU 已经能在一小时内完成一个几十次外迭代的 TRN 实验。3.2 反演主循环的最小 Python 骨架下面是一个不依赖具体 seiscope API 的 TRN 外层循环骨架重点在于表达控制流把优化器与正演解耦。实际使用时应把forward_op替换成 seiscope 的正演封装。import numpy as np class TruncatedNewtonFWI: def __init__(self, forward_op, adjoint_op, model0, params): self.fwd forward_op # 正演 self.adj adjoint_op # 伴随 self.m model0.copy() self.eta params.get(eta, 0.2) # CG 截断容差 self.cg_max params.get(cg_max, 30) # 内层最大迭代次数 def compute_gradient(self, freq): u self.fwd(self.m, freq) r u - self.data[freq] lam self.adj(self.m, freq, r) g self._apply_dLdm_adj(self.m, u, lam, freq) return g.real def compute_Hv(self, v, freq): u self.fwd(self.m, freq) du self._solve_perturbation(self.m, u, v, freq) lam self.adj(self.m, freq, self.data[freq] - u) dlam self._solve_adj_perturbation(self.m, u, du, lam, v, freq) Hv ( self._apply_dLdm_adj(self.m, u, dlam, freq) self._apply_dLdm_adj(self.m, du, lam, freq) ) return Hv.real def outer_step(self, freq): g self.compute_gradient(freq) p self._truncated_cg(g, freq) # 内层 CG 求解 H p -g alpha self._line_search(g, p, freq) self.m alpha * p return np.linalg.norm(g)逻辑说明compute_gradient先正演得到波场残差再走一次伴随得到梯度compute_Hv实现的是上一章的四波场流程。_truncated_cg内部运行共轭梯度每轮只需要调用compute_Hv不需要知道海森矩阵的显式形式。_line_search采用回溯线搜索目标函数下降就接受步长下降过少就减半。参数说明里最重要的三个量eta越小内层 CG 解越精确但单次外迭代成本越高cg_max是内层 CG 次数上限防止在大模型上无限制迭代freq是当前参与反演的频率频率延续策略下每次外迭代可以换一个频率。3.3 和 seiscope 模型对接先做梯度校验再谈收敛接 seiscope 时最常见的错误是梯度方向错了但是反演仍然在前几次迭代里“看着在下降”后边突然发散。为此正式跑 TRN 前一定要做一次有限差分梯度校验。方法是对模型m加入小扰动δ·v比较解析梯度与差分梯度的点积def check_gradient(fwi, m, v, eps1e-6, freq5.0): g_analytic fwi.compute_gradient(freq) f_plus fwi.objective(m eps * v, freq) f_minus fwi.objective(m - eps * v, freq) g_fd (f_plus - f_minus) / (2 * eps) ratio np.dot(g_analytic, v) / g_fd print(ratio:, ratio) # 接近 1.0 说明梯度实现正确eps取1e-6到1e-7太大则差分误差主导太小则浮点噪声主导。ratio在0.9 ~ 1.1之间视为通过。seiscope 的建模器可能需要额外的吸收边界项参与梯度这会导致解析梯度和差分梯度系统性偏离所以要先把边界条件固定住再测。这一步通过之后TRN 的外层方向计算才有意义。4. 影响截断牛顿优化的参数与收敛控制4.1 内层迭代容差、线搜索与正则化的三件套TRN 实际调参时先看一眼收敛曲线是哪种形态。下表是参数影响速查表按重要性排列参数常见范围调小的影响调大的影响内层 CG 容差eta0.1 ~ 0.5方向更接近牛顿方向单步成本高方向质量下降近似退化为梯度法内层 CG 最大迭代数10 ~ 50每个频率计算快收敛慢计算量线性增加可能过度拟合海森近似误差回溯线搜索衰减因子0.5 ~ 0.8步长收缩快迭代次数变多步长收缩慢容易越过下降区间模型正则化权重1e-4 ~ 1e-2反演结果更容易出现高波数噪声模型过于平滑低波数信息被压制eta的设定有一个局部最优区间。把它设成 0.05 时内外层迭代路径看起来非常干净但每次外迭代的波场模拟次数可能翻倍设成 0.6 以上时截断后的方向里还残留明显的负曲率分量线搜索被迫频繁缩短步长。实际数据上更推荐在反演前期用较大的eta0.3 左右快速进入目标函数主谷后期把eta调小到 0.1 收敛局部细节。这个策略和频率延续的节奏配合起来非常顺手。回溯线搜索里还有一层容易被忽略的参数就是初始步长。TRN 的搜索方向尺度受海森近似影响并不是单位化的所以初始步长建议用α0 -g·p / (p·H·p)来估计这个值可以直接从 CG 最后一步的结果中取不需要额外正演。4.2 频率延续策略对 TRN 收敛的帮助由于 FWI 目标函数非凸TRN 和梯度法一样需要低频先行。常见做法是从 3~5 Hz 开始每个频率做若干次外迭代后向高频推进。seiscope 的示例里通常把频率组分拆成数组TRN 外循环在相邻频率间采用上一频率的解作为新频率的初始模型。这个 Warm Start 对 TRN 尤其重要因为牛顿方向容易把模型推入一个对高频数据拟合更好、但低频信息已经丢失的局部极小值里。低频阶段建议低频迭代每频率只做 5~10 次外迭代。过多反而会让模型被低频数据里的振幅误差带偏。高频阶段再把这组频率的eta调小让模型精细尺度被逐步锁住。如果内存允许可以把多个频率同时放进目标函数里求和梯度这会让 Hv 积的波场模拟次数按频率数倍增所以更常见的做法是每次只激活一个频率。4.3 三个常见坑内层跑满、阻尼不当、边界污染第一个坑是把cg_max设得很大比如 100却不提高eta的精度要求。结果内层 CG 在后半段其实是在追海森近似本身的噪声尤其在使用有限差分正演且数据含噪时CG 残差会停滞不降白白烧算力。解决方法是让cg_max与eta联动eta0.2时cg_max设 20~30 就够。第二个坑是正则化权重设上头。TRN 对海森矩阵的正定化依赖比梯度法更强因为牛顿方向要求 H 正定才能得到下降方向当正则化过小时 CG 可能遇到不收敛或负曲率方向过大时反演分辨率急剧下降。建议用自适应正则化每次外迭代检查 CG 是否在少量迭代内产生非下降搜索方向若是则把正则化权重乘上 1.5否则逐步恢复。第三个坑是吸收边界条件在 Hv 积里被忽略。全波形反演的梯度公式假定正演算子完美满足边界条件但实际有限差分实现里边界系数同样依赖于模型参数如果不把边界对模型参数的导数包含进∂L/∂m·v最终的近似海森方向会在近边界区域产生五到十道快速变化的噪声条纹。排查方法是把目标梯度集中到模型边界附近的三个网格点做有限差分校验如果比值在边界处明显偏离 1.0就说明边界导数项被遗漏了。5. 实战验证用合成数据判断 TRN 是否真正收敛5.1 设计一个能区分 TRN 与梯度法的小模型用一个 64×64 网格的二维声波走时模型做验证。背景速度 2500 m/s在中心位置放置一个半径 6 个网格、速度 2800 m/s 的圆形异常体观测来自 8 个炮点、每炮 32 个检波器。低频从 2 Hz 开始这个设置在几百次迭代内能产生肉眼可辨认的收敛差异。运行 TRN 时建议每隔几次外迭代把模型和收敛曲线同时画出来。正常收敛的曲线形态是前三次外迭代目标函数快速下降随后进入一个平缓台阶。对比使用 L-BFGS 的结果TRN 在台阶处能继续向下走出一个小幅但明确的下降这正是海森方向里的二阶信息在低速异常体边界处起作用的表现。5.2 从三张图判断问题出在优化器还是正演器第一张图是目标函数曲线。如果曲线出现“锯齿形”的反复震荡优先怀疑线搜索的初始步长估算不对回到上一章α0 -g·p / (p·H·p)的做法重新计算。第二张图是梯度范数曲线。梯度范数平稳下降是最好的信号如果梯度范数下降但目标函数不降说明搜索方向里有显著的非线性残差项需要把eta调小如果梯度范数先降后飙升这通常是 CG 迭代超限内层求解出来的p根本不是可信下降方向。第三张图是模型剖面沿垂直方向的切片。TRN 在正确收敛时会把速度异常体的边缘锐化切片的过渡带只有两到三个网格而 L-BFGS 得到的过渡带往往更宽且含有旁瓣振荡。看到这种差异才说明你真正把截断牛顿优化的能力用上了。最后运行一次check_gradient并在脚本开头打印ratio如果该值不是接近 1.0优化器本身调得再好也无济于事问题一定出在建模器到优化器的接口这一层。本文还有配套的精品资源点击获取 
