3. TEM数据预处理：坏点剔除与衰减曲线平滑

为什么坏点剔除必须先于平滑

TEM衰减曲线在早期道（约10⁻⁵～10⁻⁴ s）信号强、信噪比高，晚期道信号衰减到nV级，容易受工频、天电、电极接触不良等干扰出现孤立尖跳或整段漂移。若直接做平滑，尖跳会被"抹"进邻道，形成虚假的缓变异常，反而掩盖采空区低阻体引起的真实衰减斜率变化。因此处理顺序应为：先按规则识别并剔除坏点，再对保留点做保形平滑。

坏点识别规则设计

针对单测点衰减曲线，采用三条可量化规则联合判定，任一命中即标记为坏点：

• 规则A（局部跳变）：对第 i 道，计算二阶差分 |V[i-1] − 2V[i] + V[i+1]| / V[i]，超过阈值 τ_jump（示例取 0.35）判为尖跳。
• 规则B（单调性破坏）：TEM感应电动势在无极性反转的正常衰减段应单调递减，若 V[i] > V[i-1]·(1+τ_mono)（τ_mono 示例取 0.05）判为异常回升。
• 规则C（邻点偏离）：用左右各2点的中值作参考 ref[i]，若 |V[i] − ref[i]| / ref[i] > τ_dev（示例取 0.5）判为离群。

剔除后对空缺道用相邻有效道在对数域线性插值补齐，避免平滑时出现断点。

保形平滑：对数域滑动中值 + 三点加权

衰减曲线跨数个数量级，直接在线性域平滑会压制晚期弱信号。做法是：先取 y = ln(V)，在 y 上做窗口宽 3 的滑动中值抑制残余尖点，再做 [0.25, 0.5, 0.25] 三点加权平均降低随机噪声，最后 exp 还原。中值+加权的组合对单调衰减段近似保形，不会像大窗口均值那样把采空区引起的"衰减变缓"抹平。

可运行的Python实现

下面代码对一条含噪衰减曲线完成坏点标记、插值补齐与保形平滑，并输出平滑前后晚期道斜率对比，用于验证异常特征是否保留。

import numpy as np

# 示例数据：20道衰减曲线（单位 nV），含2个坏点
t = np.logspace(-5, -2, 20)          # 时间道 s
V = 1200 * np.exp(-t / 3e-3) + 8     # 理论衰减 + 偏置
rng = np.random.default_rng(7)
V_noisy = V * (1 + 0.03 * rng.standard_normal(20))
V_noisy[5] *= 1.9                    # 注入尖跳坏点
V_noisy[12] *= 1.6                   # 注入回升坏点

def detect_bad(V, tau_jump=0.35, tau_mono=0.05, tau_dev=0.5):
    n = len(V)
    bad = np.zeros(n, dtype=bool)
    for i in range(1, n - 1):
        # 规则A 二阶差分
        d2 = abs(V[i-1] - 2*V[i] + V[i+1]) / max(V[i], 1e-12)
        if d2 > tau_jump:
            bad[i] = True
        # 规则B 单调性破坏
        if V[i] > V[i-1] * (1 + tau_mono):
            bad[i] = True
        # 规则C 邻点中值偏离
        lo, hi = max(0, i-2), min(n, i+3)
        ref = np.median(np.delete(V[lo:hi], i-lo))
        if abs(V[i] - ref) / max(ref, 1e-12) > tau_dev:
            bad[i] = True
    return bad

def fill_and_smooth(V, bad):
    Vf = V.copy()
    idx = np.arange(len(V))
    good = ~bad
    # 对数域线性插值补齐坏点
    Vf[bad] = np.exp(np.interp(idx[bad], idx[good], np.log(V[good])))
    # 对数域滑动中值（窗口3）
    y = np.log(Vf)
    ym = y.copy()
    for i in range(1, len(y) - 1):
        ym[i] = np.median(y[i-1:i+2])
    # 三点加权 [0.25, 0.5, 0.25]
    ys = ym.copy()
    for i in range(1, len(ym) - 1):
        ys[i] = 0.25*ym[i-1] + 0.5*ym[i] + 0.25*ym[i+1]
    return np.exp(ys), bad

bad = detect_bad(V_noisy)
V_smooth, _ = fill_and_smooth(V_noisy, bad)

# 晚期道斜率对比（最后5道，对数域线性拟合）
late = slice(-5, None)
k_raw = np.polyfit(t[late], np.log(V_noisy[late]), 1)[0]
k_sm  = np.polyfit(t[late], np.log(V_smooth[late]), 1)[0]
print("坏点道号:", np.where(bad)[0])
print(f"晚期道斜率 原始={k_raw:.1f}  平滑={k_sm:.1f}  相对变化={abs(k_sm-k_raw)/abs(k_raw)*100:.1f}%")

示例运行预期输出：坏点道号为 [5, 12]；晚期道斜率相对变化通常在数个百分点以内。若相对变化超过约5%，说明平滑窗口过宽或阈值过松，已开始削弱真实衰减特征，应缩小窗口或收紧 τ 值。