2. 地震资料品质分析与频谱特征提取

为什么先做品质分析，而不是直接上属性

地震相识别与河道砂体预测对资料品质高度敏感：主频决定可分辨的砂体厚度下限，频带宽度决定薄层调谐与波形分类的稳定性，信噪比决定相干、曲率、波形聚类等属性是否可用。若在低信噪比或窄频带数据上直接做地震相划分，得到的“相带”往往只是噪声团块或调谐假象。因此本例的任务是：在例1建立的工区中，锁定目的层时窗，完成频谱分析与信噪比评估，输出可量化、可复核的品质参数表，作为后续方法选择的输入条件。

分析对象与关键约定

• 分析时窗：以目的层顶界 T0 为基准，向上取 20 ms、向下取 40 ms，形成 60 ms 的滑动分析时窗；若目的层厚度变化大，则改为“顶界—底界”层间时窗，并在成果表中注明时窗类型。
• 空间抽样：全区按 500 m × 500 m 网格抽取分析道集，另在河道可能发育区加密至 250 m × 250 m，避免只用单一井旁道代表全区。
• 振幅约定：统一使用相对振幅保真处理后的数据体，分析前不做 AGC；若必须做增益，需在成果表中记录增益类型，因为增益会显著改变频谱形态。
• 频率单位：统一用 Hz；采样率 dt 以 ms 表示，Nyquist 频率 = 1000/(2·dt)。例如 dt = 2 ms 时 Nyquist 频率为 250 Hz，频谱分析上限不应超过该值。

频谱特征提取的实现路径

核心做法是对时窗内地震道做加窗傅里叶变换，得到振幅谱后统计主频、频带宽度与低频/高频能量占比。下面给出可直接运行的分析脚本（Python + NumPy/SciPy），输入为按道存储的时窗数据矩阵，输出为逐道频谱参数与全区统计量。

import numpy as np
from scipy.signal import windows

def spectrum_params(traces, dt_ms, f_low_ratio=0.5, f_high_ratio=0.5):
    """
    traces: shape (n_traces, n_samples)，时窗内地震道，相对振幅
    dt_ms : 采样间隔(ms)
    f_low_ratio / f_high_ratio: 以峰值振幅的该比例作为频带边界
    返回: 每道的主频、低频端、高频端、带宽、峰值振幅
    """
    n_tr, n_s = traces.shape
    fs = 1000.0 / dt_ms                 # 采样频率 Hz
    freqs = np.fft.rfftfreq(n_s, d=1.0/fs)
    win = windows.hann(n_s)             # 汉宁窗抑制截断泄漏
    out = np.zeros((n_tr, 5))
    for i in range(n_tr):
        x = traces[i] - np.mean(traces[i])   # 去直流
        amp = np.abs(np.fft.rfft(x * win))
        peak_idx = np.argmax(amp)
        peak_amp = amp[peak_idx]
        f_peak = freqs[peak_idx]
        # 从峰值向两侧搜索，首次低于阈值处作为频带边界
        lo = peak_idx
        while lo > 0 and amp[lo] >= peak_amp * f_low_ratio:
            lo -= 1
        hi = peak_idx
        while hi < len(amp)-1 and amp[hi] >= peak_amp * f_high_ratio:
            hi += 1
        f_lo, f_hi = freqs[lo], freqs[hi]
        out[i] = [f_peak, f_lo, f_hi, f_hi - f_lo, peak_amp]
    return out

# 示例：3 道、时窗 60 ms、dt = 2 ms
dt_ms = 2.0
n_s = int(60 / dt_ms) + 1
rng = np.random.default_rng(7)
t = np.arange(n_s) * dt_ms
# 构造示例道：主频 28 Hz 的 Ricker 子波 + 不同强度噪声
def ricker(f, t_ms):
    tt = (t_ms - 30.0) / 1000.0
    return (1 - 2*(np.pi*f*tt)**2) * np.exp(-(np.pi*f*tt)**2)

clean = ricker(28, t)
traces = np.vstack([
    clean + 0.02*rng.standard_normal(n_s),   # 高信噪比
    clean + 0.15*rng.standard_normal(n_s),   # 中等信噪比
    clean + 0.45*rng.standard_normal(n_s),   # 低信噪比
])
res = spectrum_params(traces, dt_ms)
for i, r in enumerate(res):
    print(f"道{i+1}: 主频={r[0]:.1f} Hz, 频带={r[1]:.1f}-{r[2]:.1f} Hz, "
          f"带宽={r[3]:.1f} Hz, 峰值振幅={r[4]:.3f}")

该脚本的关键点在于：加汉宁窗减少时窗截断造成的频谱泄漏；去直流避免 0 Hz 能量抬高低频端；频带边界用“峰值振幅的 50%”定义（即 -6 dB 点），这一约定必须在成果表中写明，否则不同人用 -3 dB 或 -12 dB 会得到不同的带宽数值。示例数据为合成道，仅用于说明算法行为，不代表实际工区实测结果。

信噪比评估的两种可操作方法

信噪比没有唯一算法，工程上常用以下两种，建议同时计算并交叉验证：

• 频谱能量比法：在目的层时窗内，取信号主频频带（如 f_lo—f_hi）内能量与高频端噪声频带（如 0.7—0.9 倍 Nyquist 频率）能量的比值，适合快速批量评估。
• 相邻道相似法：计算分析道与周围 4 道的归一化互相关最大值，取平均作为信噪比的相对度量。河道边界处横向变化快，该值会自然降低，需与频谱法结合判读，避免把地质边界误判为低信噪比。

品质参数与后续方法选择的对应关系

成果表与判读要点

建议按“分区 + 参数”的形式输出成果表，而不是只给全区一个平均值。下面是一份可直接填写的示例结构（数值为示例，需用实际工区数据替换）：

完成本例后，应得到三样可直接使用的东西：目的层时窗内的主频与带宽数值、分区信噪比等级、以及一份注明约定（时窗、增益、频带边界比例）的品质参数表。这三项将直接决定例4中属性提取的时窗设置与属性类型取舍。