2. 目标函数量级失衡导致的不收敛排查

问题定位：先看两项在目标函数中的实际贡献

目标函数通常写成 Φ(m) = Φ_d(m) + λ·Φ_m(m)，其中 Φ_d 是数据残差项，Φ_m 是模型正则项，λ 是权重。不收敛的常见诱因不是某一项算错，而是两项在数值上相差几个数量级：例如 Φ_d ≈ 1e6、Φ_m ≈ 1e-2，此时梯度几乎完全由数据项主导，正则项形同虚设；反过来若 Φ_m 远大于 Φ_d，迭代会过早被拉向先验模型，残差停滞不降。

排查的第一步不是改代码，而是把每次迭代的 Φ_d、Φ_m、λ·Φ_m 和总目标 Φ 分别打印出来，观察它们的量级比 r = Φ_d / (λ·Φ_m)。经验上 r 落在 1e-2 ~ 1e2 之间时两项能共同起作用；若 r > 1e4 或 r < 1e-4，就应优先怀疑量级失衡，而不是步长或初始模型。

构造一个可复现的失衡算例

下面用一个一维线性反演的最小例子把失衡现象做出来：模型 m 为 40 个网格点的慢度扰动，正演矩阵 G 为 60×40 的平滑核，观测 d = G·m_true + 噪声。目标函数取 Φ = ||Gm − d||² + λ·||Lm||²，L 为一阶差分算子。故意把 λ 设成 1e-8，使正则项贡献远小于数据项，再用固定步长梯度下降迭代 200 次，记录每步的 Φ_d、Φ_m 与 Φ。

import numpy as np

rng = np.random.default_rng(7)
n_data, n_model = 60, 40

# 平滑正演核：每个数据点受邻近若干模型点加权
G = np.zeros((n_data, n_model))
for i in range(n_data):
    c = i * (n_model - 1) / (n_data - 1)
    for j in range(n_model):
        G[i, j] = np.exp(-0.5 * ((j - c) / 3.0) ** 2)
G /= G.sum(axis=1, keepdims=True)

m_true = np.sin(np.linspace(0, 3 * np.pi, n_model)) * 0.5
d = G @ m_true + rng.normal(0, 0.01, n_data)

# 一阶差分正则算子
L = np.zeros((n_model - 1, n_model))
for i in range(n_model - 1):
    L[i, i], L[i, i + 1] = -1.0, 1.0

def objective(m, lam):
    r = G @ m - d
    phi_d = r @ r
    phi_m = (L @ m) @ (L @ m)
    return phi_d, phi_m, phi_d + lam * phi_m

def grad(m, lam):
    return 2 * G.T @ (G @ m - d) + 2 * lam * L.T @ (L @ m)

def run(lam, step, iters=200):
    m = np.zeros(n_model)
    hist = []
    for k in range(iters):
        pd, pm, pt = objective(m, lam)
        hist.append((k, pd, pm, pt))
        m = m - step * grad(m, lam)
    return np.array(hist), m

# 失衡配置：λ 过小
hist_bad, m_bad = run(lam=1e-8, step=2e-3)
for k in (0, 50, 100, 150, 199):
    kk, pd, pm, pt = hist_bad[k]
    print(f"iter {int(kk):3d}  Φd={pd:.3e}  Φm={pm:.3e}  λΦm={1e-8*pm:.3e}  Φ={pt:.3e}")

典型输出（示例数据，实际数值随随机种子略有差异）会显示 Φd 在 1e0 量级、Φm 在 1e-1 量级，而 λ·Φm 只有 1e-9 量级，量级比 r ≈ 1e9。此时 Φ 曲线几乎完全跟随 Φd，正则项对梯度贡献可忽略，模型在数据稀疏区域出现尖峰，残差在 100 次迭代后开始小幅振荡而不继续下降。

归一化与调权：让两项回到同一量级

修正分两步，顺序不能颠倒。第一步做无量纲归一化：把 Φd 除以数据点数 N_d，把 Φm 除以模型参数个数 N_m，消除“项的大小随网格规模变化”的假象。第二步再调 λ，使 λ·Φm/N_m 与 Φd/N_d 处于同一数量级。一个实用的初值估计是 λ₀ = (Φd/N_d) / (Φm/N_m)，在第一次迭代时计算一次即可。

# 归一化目标函数
def objective_norm(m, lam):
    r = G @ m - d
    phi_d = (r @ r) / n_data
    phi_m = ((L @ m) @ (L @ m)) / n_model
    return phi_d, phi_m, phi_d + lam * phi_m

def grad_norm(m, lam):
    return (2 / n_data) * G.T @ (G @ m - d) + (2 * lam / n_model) * L.T @ (L @ m)

# 用初始模型估计 λ0
m0 = np.zeros(n_model)
pd0, pm0, _ = objective_norm(m0, 1.0)
lam0 = pd0 / pm0
print(f"λ0 = {lam0:.3e}")

def run_norm(lam, step, iters=200):
    m = np.zeros(n_model)
    hist = []
    for k in range(iters):
        pd, pm, pt = objective_norm(m, lam)
        hist.append((k, pd, pm, pt))
        m = m - step * grad_norm(m, lam)
    return np.array(hist), m

hist_fix, m_fix = run_norm(lam=lam0, step=2e-3)
for k in (0, 50, 100, 150, 199):
    kk, pd, pm, pt = hist_fix[k]
    print(f"iter {int(kk):3d}  Φd/Nd={pd:.3e}  λΦm/Nm={lam0*pm:.3e}  Φ={pt:.3e}")

调整后 Φd/N_d 与 λ·Φm/N_m 都落在 1e-2 ~ 1e-1 量级，量级比接近 1。此时梯度中两项贡献相当，残差曲线应呈现单调下降并在 150 次迭代附近趋于平台，模型不再出现高频尖峰。

验证收敛路径是否真的改变

只看最终残差不够，要对比调整前后的收敛路径。建议同时记录三条曲线：Φd/N_d、λ·Φm/N_m 和总目标 Φ，并计算相邻两次迭代的相对下降率 |Φ_{k+1} − Φ_k| / Φ_k。失衡时该下降率会在 1e-3 附近反复跳动；归一化调权后应稳定衰减到 1e-6 以下。若调权后残差下降但模型出现新的振荡，说明 λ 调得过大，应回到 λ₀ 附近做小范围扫描（如 0.1λ₀、λ₀、10λ₀ 三档），选残差平台最低且模型最平滑的一档。