4. 带通滤波参数扫描与频谱验证

为什么先做参数扫描，而不是直接套一组“常用频带”

带通滤波的效果取决于四个拐点频率：低截 f1、低通 f2、高通 f3、高截 f4（单位 Hz）。对 CMP 道集而言，目标不是“滤得干净”，而是保留有效反射的频带、同时压低面波与高频环境噪声。同一工区不同 CMP 的有效频带可能不同，因此本例采用“参数扫描 + 频谱验证”的闭环：先对同一 CMP 道集施加多组频带组合，再逐组比较滤波前后的振幅谱，用可量化的判据选出参数。

下面给出一个可直接运行的 Python 实现（依赖 numpy、scipy、matplotlib），它完成三件事：合成一个带噪声的 CMP 道集（含 3 个有效反射同相轴 + 低频面波 + 高频随机噪声）、对 5 组频带做零相位带通滤波、输出频谱对比图与频带能量占比表。

import numpy as np
from scipy.signal import butter, sosfiltfilt
import matplotlib.pyplot as plt

# ---------- 1. 合成 CMP 道集（示例数据，非实测） ----------
dt, nt, ntr = 0.002, 1000, 24          # 2 ms 采样, 2 s 记录, 24 道
t = np.arange(nt) * dt
offsets = np.linspace(100, 2400, ntr)  # 偏移距 100~2400 m

rng = np.random.default_rng(7)
data = np.zeros((ntr, nt))
for i, x in enumerate(offsets):
    for t0, v, amp in [(0.45, 2200, 1.0), (0.95, 2800, 0.8), (1.45, 3300, 0.6)]:
        tt = t0 + x**2 / (2 * v**2 * t0)          # 双曲线时距（NMO 前）
        data[i] += amp * np.exp(-((t - tt) / 0.018)**2) * np.sin(2*np.pi*35*(t-tt))
    data[i] += 0.9 * np.sin(2*np.pi*6*t) * np.exp(-t/0.6)   # 低频面波
    data[i] += 0.25 * rng.standard_normal(nt)               # 宽带随机噪声

# ---------- 2. 频带组合扫描 ----------
bands = [(5, 10, 60, 90), (8, 12, 55, 80), (10, 15, 50, 70),
         (12, 18, 45, 65), (15, 20, 40, 60)]
fs = 1.0 / dt
results = []
for f1, f2, f3, f4 in bands:
    sos = butter(4, [f1, f2, f3, f4], btype='bandpass', fs=fs, output='sos')
    filt = sosfiltfilt(sos, data, axis=1)      # 零相位，避免同相轴时移
    results.append((f1, f2, f3, f4, filt))

# ---------- 3. 频谱对比 ----------
def amp_spec(seg):
    n = seg.shape[1]
    f = np.fft.rfftfreq(n, dt)
    a = np.abs(np.fft.rfft(seg, axis=1)).mean(axis=0)
    return f, a

f, a_raw = amp_spec(data)
fig, ax = plt.subplots(figsize=(9, 5))
ax.semilogy(f, a_raw, 'k', lw=2, label='原始道集')
colors = ['#1f77b4', '#2ca02c', '#ff7f0e', '#9467bd', '#d62728']
for (f1, f2, f3, f4, filt), c in zip(results, colors):
    _, a_f = amp_spec(filt)
    ax.semilogy(f, a_f, c, lw=1.2, label=f'{f1}-{f2}-{f3}-{f4} Hz')
ax.set_xlim(0, 120); ax.set_xlabel('频率 / Hz'); ax.set_ylabel('平均振幅谱')
ax.legend(fontsize=8); ax.grid(alpha=0.3)
plt.tight_layout(); plt.show()

# ---------- 4. 频带能量占比（判据） ----------
def band_ratio(seg, lo, hi):
    f, a = amp_spec(seg)
    m = (f >= lo) & (f <= hi)
    return a[m].sum() / a.sum()

print('频带组合        有效带(15-50Hz)占比   噪声带(>70Hz)占比')
for f1, f2, f3, f4, filt in results:
    print(f'{f1}-{f2}-{f3}-{f4} Hz   {band_ratio(filt,15,50):.3f}'
          f'                {band_ratio(filt,70,120):.3f}')

运行后应看到：原始谱在 6 Hz 附近有面波能量峰、在 70 Hz 以上有噪声底；随着低截频率从 5 Hz 抬到 15 Hz，低频端被逐步削掉，而 15–50 Hz 的有效反射频带保持平坦。判据是有效带占比尽量高、噪声带占比尽量低，同时目视检查同相轴连续性未被破坏。

参数选择与验证要点

• 零相位滤波：使用 sosfiltfilt 而非单向 sosfilt，避免滤波引入相位延迟导致同相轴时移，否则后续速度谱拾取会出现系统性偏差。
• 拐点不要过陡：4 阶 Butterworth 在 f2 与 f3 之间形成平顶，若阶数过高会在拐点处产生振铃，反而污染弱反射。
• 判据要双约束：只看“噪声带占比低”会诱导把频带收得过窄，导致有效反射高频成分丢失、分辨率下降；因此必须同时看有效带占比。
• 逐 CMP 复核：本例只对一个 CMP 扫描。实际生产中应抽取浅、中、深不同时段的若干 CMP 分别扫描，若最优频带差异明显，则需分时窗或分区域设定参数，而不是全数据统一一组。