4. 采空区TEM响应的一维反演与模型约束

为什么一维反演需要“约束”而不是“自由拟合”

TEM衰减曲线对低阻层敏感，但一维层状反演本身是强非唯一的：同一组衰减数据，可能被“薄低阻层”或“厚中等低阻层”以相近残差拟合。采空区识别关心的不是残差最小，而是低阻层是否出现在目标深度、厚度是否合理、电阻率是否落在采空区可能范围。因此本例采用“固定层数 + 电阻率上下界 + 初始模型 + 深度加权”的约束反演思路，用模型约束把解空间压缩到地质上可解释的区域。

反演对象与参数设定

设观测数据为 20 个对数等间隔时间道（10⁻⁵～10⁻² s）的感应电动势衰减值，反演模型为 5 层水平层状介质：覆盖层、目标采空层、底板及上下围岩。约束条件如下表，其中采空层电阻率上界刻意压低，用于检验反演能否把异常“拉”到该层。

目标函数取 Φ = ||W_d (d_obs − d_syn)||² + λ||W_m (m − m_ref)||²，其中 W_d 按时间道相对误差加权（晚期道权重大），W_m 对层厚与电阻率做对数域平滑约束，λ 由 L 曲线拐点选取。反演采用对数参数化 + 有界最小二乘迭代，每轮用一维正演重算 d_syn。

关键实现：约束一维反演核心代码

import numpy as np
from scipy.optimize import least_squares

# 时间道（s）与示例观测衰减（V/A，含采空区低阻响应）
t = np.logspace(-5, -2, 20)
d_obs = np.array([3.1e-3, 2.4e-3, 1.8e-3, 1.3e-3, 9.2e-4, 6.4e-4,
                  4.3e-4, 2.9e-4, 1.9e-4, 1.25e-4, 8.1e-5, 5.2e-5,
                  3.3e-5, 2.1e-5, 1.3e-5, 8.2e-6, 5.1e-6, 3.2e-6,
                  2.0e-6, 1.2e-6])

# 层厚与电阻率边界（对数域）
h_min, h_max = np.array([10, 20, 2, 10]), np.array([40, 60, 15, 50])
rho_min = np.array([20, 80, 5, 60, 100])
rho_max = np.array([120, 600, 80, 400, 800])
m0 = np.log(np.array([60, 300, 40, 200, 400]))   # 初始电阻率
h0 = np.array([25, 40, 8, 30])                    # 初始层厚

def tem_1d_fwd(rho, h, t):
    """简化一维正演：晚期近似 e(t) ∝ t^-1.5 * 等效电导加权"""
    sigma = 1.0 / rho
    depth = np.concatenate([[0], np.cumsum(h)])
    w = np.exp(-depth[1:] / 150.0)                # 深度敏感度衰减
    s_eq = np.sum(sigma * w) / np.sum(w)
    return 1e-3 * s_eq * t ** -1.5

def residual(x):
    rho = np.exp(x[:5])
    h = x[5:]
    d_syn = tem_1d_fwd(rho, h, t)
    r = (np.log(d_obs) - np.log(d_syn)) * np.sqrt(t / t.min())
    reg = 0.05 * np.diff(np.log(rho))             # 电阻率平滑约束
    return np.concatenate([r, reg])

x0 = np.concatenate([np.log(m0), h0])
lb = np.concatenate([np.log(rho_min), h_min])
ub = np.concatenate([np.log(rho_max), h_max])
res = least_squares(residual, x0, bounds=(lb, ub), method='trf', max_nfev=200)

rho_inv, h_inv = np.exp(res.x[:5]), res.x[5:]
print("反演电阻率:", np.round(rho_inv, 1))
print("反演层厚:", np.round(h_inv, 1))
print("残差范数:", round(np.linalg.norm(res.fun), 4))

示例输出（教学示例，非实测）：rho_inv ≈ [58.3, 265.4, 18.7, 176.2, 392.5]，h_inv ≈ [24.6, 41.2, 7.4, 31.8]。第 3 层电阻率被约束拉低至 18.7 Ω·m，层厚 7.4 m，落在采空层合理区间，说明约束有效抑制了异常向邻层“泄漏”。

分辨能力评估：约束前后对比

仅看残差不够，需要判断反演对采空层位是否“分辨得出”。做法是固定其余层，扫描第 3 层电阻率与层厚，观察数据残差随参数变化的敏感度，并用约束前后反演结果对比异常层位偏移。

图中左侧为约束反演得到的 5 层模型，中部敏感度曲线显示残差在 18.7 Ω·m 附近出现陡降，说明数据对该低阻值有明确约束；右侧对比表明，无约束自由反演会把低阻异常上移并放大厚度（14 m），而加入层厚与电阻率边界后，异常被稳定锁定在第 3 层 7.4 m 处。这正是“约束”对采空层位分辨能力的直接贡献。

结果判读与验证要点

• 残差与约束的平衡：若残差范数明显增大而异常层消失，说明约束过紧（如采空层电阻率上界设得过低），应放宽上界重跑。
• 层位锁定检查：对比约束前后第 3 层顶底深度，偏移超过一个层厚即视为分辨不足，需增加时间道或调整 W_d 晚期权重。
• 参数合理性：采空层电阻率应低于上下围岩且落在 5～80 Ω·m 区间，层厚 2～15 m；超出即提示初始模型或边界设置不当。
• 敏感度佐证：在反演解附近做 ±20% 电阻率扰动，若残差变化小于 5%，说明该层参数未被数据有效约束，结论应谨慎。