在石油勘探和地质研究中，测井曲线就像是给地层做的CT扫描，能反映出不同深度地层的物理性质，比如放射性（伽马射线）、电阻率、密度、孔隙度等。但这些曲线往往夹杂着噪音、异常值，且不同物理曲线反映的地层特征尺度不同，有的反映厚层趋势（如泥岩层），有的反映薄层或突变（如砂岩与泥岩的交界面）。传统的分析方法很难把这些不同尺度的特征同时提取清楚。
小波变换就像一把数学显微镜，它可以把原始信号分解成不同频率的成分：低频部分对应地层的大尺度变化（比如从泥岩到砂岩的缓慢过渡），高频部分对应小尺度细节（如薄互层、冲刷面）。特别是小波包分解，它比普通小波分解更精细，能把高频部分也继续分解，从而更准确地定位那些薄而关键的地质体，比如泥岩冲刷面、砂岩体。
提出方法对伽马射线（GR）曲线做小波包分解，然后看各个频带能量分布，找到对地质解释最有用的那几个频带（比如低频A5反映大套泥岩，高频D5反映薄层突变），再结合电阻率、密度等曲线，就能把地层中的关键界面和岩性变化自动识别出来，为地质建模和储量评价提供定量依据。

这个算法其实就是用数学上的小波包分解，把原本混在一起的测井信号拆成不同频率的成分，然后挑出其中最能代表大套岩性变化和薄层突变的那几个成分，再结合其他测井曲线（电阻率、密度等）相互印证，从而自动找到地层中的关键界面（比如泥岩冲刷面）和有利储层位置（砂岩体）。
好处在于，它不用人工去一条一条曲线地找，而是让计算机根据信号本身的频率特征去判断，既快又客观。而且通过把多种曲线（伽马、电阻率、密度）叠在一起看，能避免单条曲线误判，最后得到的地质解释更可靠。

如果你对信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测有疑问，或者需要论文思路上的建议，欢迎咨询

担任《MSSP》《中国电机工程学报》《宇航学报》《控制与决策》等期刊审稿专家，擅长领域：信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测

算法步骤

数据准备与清洗读取原始测井数据（LAS文件），把深度从3100米到4075米这一段取出来，剔除掉那些关键曲线（伽马、密度、孔隙度、电阻率等）有缺失的深度点。

异常值去除用两种方法去噪音：一种是Hampel滤波器（窗口40个点，阈值1.8倍中位数绝对偏差），主要针对突变尖峰；另一种是滚动IQR（四分位距）法，用于平稳去除那些明显偏离周围趋势的孤立点。 对伽马、密度、孔隙度、电阻率这些关键曲线分别处理，得到干净的曲线。

泥质含量计算用伽马曲线算出泥质含量Vsh（Larionov公式），然后用这个Vsh去校正孔隙度曲线（APLC），去掉泥质对孔隙度的影响，得到校正后的中子孔隙度NPHI_shcorr。

小波包分解把清洗后的伽马曲线做小波包分解（选db7小波，分解3层），得到8个频带（从很低频到很高频）。计算每个频带的能量占比，看看能量主要集中在哪里——能量高的频带说明该频带对地层变化贡献大。

关键频带重建与解释用小波包分解得到的系数，只保留A5（最低频）和D5（较高频）进行信号重建。

A5反映大尺度趋势，用来识别“砂体”（SB），也就是大套砂岩发育的地方。

D5反映小尺度突变，用来识别“泥岩冲刷面”（MFS），也就是沉积环境突然变化的位置。

联合多曲线验证把A5、D5的识别结果与电阻率、密度、孔隙度曲线放在一起看，确定最终的MFS和SB位置，形成最终的地质解释结果。

def wavelet_packet_decompose(signal, wavelet='db7', level=3):
    """
    对信号进行小波包分解，返回各频带信号及能量占比
    """
    # 小波包分解
    wp = pywt.WaveletPacket(data=signal, wavelet=wavelet, mode='symmetric', maxlevel=level)
    nodes = wp.get_level(level, order='freq')  # 获取第level层的所有节点，按频率排序
    signals = [n.data for n in nodes]

    # 统一长度（取最短的）
    min_len = min(len(s) for s in signals)
    signals = [s[:min_len] for s in signals]

    # 计算各频带能量占比
    energies = np.array([np.sum(s**2) for s in signals])
    energies_norm = energies / np.sum(energies)

    return signals, energies_norm

def reconstruct_band(coeffs, wavelet, band_name):
    """
    根据小波包系数重建指定频带
    band_name: 例如 'a5', 'd5' 等
    """
    coeffs_zero = [np.zeros_like(c) for c in coeffs]
    # 假设 coeffs 的顺序为 [a5, d5, d4, d3, d2, d1]
    index_map = {'a5': 0, 'd5': 1, 'd4': 2, 'd3': 3, 'd2': 4, 'd1': 5}
    coeffs_zero[index_map[band_name]] = coeffs[index_map[band_name]]
    return pywt.waverec(coeffs_zero, wavelet)

# ---------- 主流程 ----------
# 1. 准备干净的伽马曲线（已归一化）
gr_clean = log_clean2['HSGR_clean'].values  # 已清洗、去噪
gr_norm = (gr_clean - np.mean(gr_clean)) / np.std(gr_clean)

# 2. 小波包分解（三层）
signals, energies = wavelet_packet_decompose(gr_norm, wavelet='db7', level=3)

# 3. 用普通小波分解得到系数（用于重建）
coeffs = pywt.wavedec(gr_norm, 'db7', level=5)

# 4. 重建关键频带
rec_a5 = reconstruct_band(coeffs, 'db7', 'a5')   # 低频趋势 → 砂体趋势
rec_d5 = reconstruct_band(coeffs, 'db7', 'd5')   # 高频细节 → 冲刷面信号

# 5. 裁剪至统一深度
min_len = min(len(rec_a5), len(rec_d5))
depth_aligned = depth[:min_len]
rec_a5 = rec_a5[:min_len]
rec_d5 = rec_d5[:min_len]

# 6. 识别砂体（SB）和泥岩冲刷面（MFS）
#   砂体：低频信号（A5）的低谷区
sb_signal = -rec_a5   # 负向为砂体响应
sb_idx, _ = find_peaks(sb_signal, distance=80)   # 距离阈值80个采样点

#   冲刷面：高频信号（D5）与低频信号（A5）组合的强振幅区
mfs_signal = np.abs(rec_d5) + np.abs(rec_a5)
mfs_signal = (mfs_signal - mfs_signal.min()) / (mfs_signal.max() - mfs_signal.min())
mfs_idx, _ = find_peaks(mfs_signal, height=0.6, distance=80)

# 7. 输出结果
print("砂体（SB）位置（采样点）:", sb_idx)
print("泥岩冲刷面（MFS）位置（采样点）:", mfs_idx)

如果你对信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测有疑问，或者需要论文思路上的建议，欢迎咨询

担任《MSSP》《中国电机工程学报》《宇航学报》《控制与决策》等期刊审稿专家，擅长领域：信号滤波/降噪，机器学习/深度学习，时间序列预分析/预测，设备故障诊断/缺陷检测/异常检测