3. 初始模型选取不当引发发散的诊断与替换

问题定位：初始模型如何把迭代“推”向发散

在同一个反演问题上，初始模型 m₀ 决定了第一次线性化的工作点。若 m₀ 远离真解，梯度方向与真实下降方向偏差大，残差可能先增大再振荡；若 m₀ 接近真解，线性化误差小，残差单调下降并快速进入收敛带。诊断的关键不是“初始模型好不好看”，而是看前若干次迭代的残差变化率与模型更新量是否呈现“越走越远”的特征。

下面用一个一维层状介质直流电阻率反演的教学示例说明。设真实模型为三层：ρ = [100, 300, 50] Ω·m，层厚 h = [10, 15] m。反演参数化为每层电阻率的对数 x = ln ρ，正演采用简化的一维解析近似（仅用于演示迭代行为，不作为生产正演）。目标函数为 Φ = ‖d - F(m)‖² + λ‖m - m_ref‖²，用高斯-牛顿法迭代。

诊断判据：从迭代轨迹读出“初始模型问题”

把远、近两个初始模型分别代入同一反演程序，记录每次迭代的残差 Φ_k 与模型更新量 ‖Δm_k‖。判据如下：

• 远初始点特征：Φ₁ > Φ₀（第一步残差反而上升），随后 Φ 在若干值之间来回摆动，‖Δm_k‖ 不随迭代衰减，说明线性化点太远、高斯-牛顿方向不可信。
• 近初始点特征：Φ_k 单调下降，‖Δm_k‖ 逐次减小，通常 5~8 次迭代进入收敛带。
• 排除其他原因：若调整初始模型后仍振荡，才转向步长、线搜索或权重问题；本例只处理“换初始模型即可恢复”的情形。

可运行的教学示例：对比两种初始模型

下面代码用简化的一维正演（三层电阻率的解析近似）构造观测数据，并分别从远、近两个初始模型出发做高斯-牛顿迭代，打印残差轨迹。正演公式为教学简化式，目的是让迭代行为可复现，不用于生产解释。

import numpy as np

# 真实模型：三层电阻率 (Ω·m) 与层厚 (m)
rho_true = np.array([100.0, 300.0, 50.0])
h_true   = np.array([10.0, 15.0])

# 教学用简化正演：视电阻率对层参数的近似响应
def forward(rho, h):
    # 仅用于演示迭代行为，非生产正演
    r1, r2, r3 = rho
    h1, h2 = h
    return np.array([
        r1,
        r1 + (r2 - r1) * h1 / (h1 + h2),
        r2,
        r2 + (r3 - r2) * h2 / (h1 + h2),
        r3
    ])

d_obs = forward(rho_true, h_true)

# 反演参数：三层电阻率的对数
def unpack(x):
    return np.exp(x)

def residual(x):
    return d_obs - forward(unpack(x), h_true)

# 数值雅可比（教学用，步长 1e-4）
def jacobian(x, eps=1e-4):
    J = np.zeros((len(d_obs), len(x)))
    for j in range(len(x)):
        xp = x.copy(); xp[j] += eps
        J[:, j] = (forward(unpack(xp), h_true) - forward(unpack(x), h_true)) / eps
    return J

def gauss_newton(x0, lam=1e-3, max_iter=12):
    x = x0.copy()
    hist = []
    for k in range(max_iter):
        r = residual(x)
        phi = r @ r
        hist.append(phi)
        J = jacobian(x)
        # 带阻尼的正规方程
        A = J.T @ J + lam * np.eye(len(x))
        g = J.T @ r
        dx = np.linalg.solve(A, g)
        x = x + dx
    return x, hist

# 远初始模型：三层都取 20 Ω·m
x_far  = np.log(np.array([20.0, 20.0, 20.0]))
# 近初始模型：接近真解
x_near = np.log(np.array([80.0, 200.0, 70.0]))

x_f, hist_f = gauss_newton(x_far)
x_n, hist_n = gauss_newton(x_near)

print("远初始点残差轨迹:", np.round(hist_f, 2))
print("近初始点残差轨迹:", np.round(hist_n, 2))
print("远初始点反演结果:", np.round(unpack(x_f), 1))
print("近初始点反演结果:", np.round(unpack(x_n), 1))

示例运行后预期看到：远初始点的残差序列在前几步出现上升或平台，随后在相近量级间摆动；近初始点的残差序列单调下降并趋于小值。反演结果上，远初始点可能停在偏离真解的局部极小，近初始点更接近 [100, 300, 50]。这些数值属于教学示例的预期现象，实际运行结果会随正演近似与阻尼设置略有差异。

替换方案设计：用先验信息构造新的 m₀

替换不是“随便换个值再试”，而是把可获得的先验信息转成初始模型与参考模型：

1. 钻孔/测井约束：若有分层信息，直接把层厚与层电阻率赋值给 m₀，并把同一组值作为正则化参考 m_ref，使迭代被“拉”在合理范围内。
2. 视电阻率快速反演：没有钻孔时，先用视电阻率曲线做一次快速一维反演，把其结果作为 m₀；这一步计算量小，但能把初始点从“远离真解”拉到“同一量级”。
3. 多起点试探：构造 3~5 个量级不同的初始模型（如整体偏大、整体偏小、分层递增），各跑少量迭代，取前 5 次残差下降最快者继续完整反演。

替换效果验证清单

完成替换后，把新旧两条残差轨迹画在同一张图上对比，是判断“初始模型是否为发散主因”的最直接证据：若新轨迹单调下降而旧轨迹振荡，即可确认初始模型选取不当是本次不收敛的主因。