简介针对地震资料处理中随机噪声干扰严重、弱信号易被淹没等问题这份约8KB的MATLAB程序包以小波分析与D-S证据理论融合为技术主线面向地震勘探、信号处理与地球物理反演方向的科研人员和工程师提供了一套可用于实际资料处理的去噪参考实现。压缩包内仅含1个m脚本但功能覆盖了数据读取、小波基选择与多尺度分解、软硬阈值去噪、D-S证据融合、井曲线辅助校正、小波重构以及去噪前后剖面对比等完整流程可帮助使用者系统理解从噪声分离到多源信息整合的各个环节并掌握MATLAB中小波工具箱与证据理论的结合方式。资源体积精简仅8KB已吸引215人学习下载适合有一定MATLAB基础的学习者直接运行验证、替换自己的地震数据测试或在此基础上按需改造阈值策略与融合参数针对复杂地质条件下的噪声压制问题进行针对性优化是开展地震数据去噪研究的一手脚本参考。1. 地震数据去噪为什么绕不开小波拿到一批含噪地震记录时最先看到的往往是有效反射被面波、工业电干扰和随机噪声糊成一片。做过几年处理的人的第一反应是带通滤波但滤波后剖面变“干净”了同相轴也变“肉”了高频细节和弱反射一起被削掉。这正是地震数据去噪和小波去噪这类工具在工区处理流程里始终占有一席之地的原因地震信号本身是非平稳的有效信号与噪声在频带上大面积重叠单纯做频率切除解决不了问题。小波变换的时频局部化特性能把信号和噪声按不同尺度拆开再针对每个尺度单独压制噪声从而在保幅和去噪之间多出一个可调的维度。pengbeng_v82.zip 这一类在圈内流传的处理脚本包核心就是把“小波分解 → 阈值处理 → 重构”这套流程封装成可直接运行的模块省去自己从头写边界处理和系数选择的功夫。这篇文章就顺着这条技术线把原理、参数、可复现代码和实际踩过的坑一次讲完。2. 小波去噪用于地震数据的基本原理与选型2.1 为什么带通滤波压不掉与有效波同频的噪声地震记录中的随机噪声并不总是白噪很多时候它的频谱和有效反射波的频谱高度重合。比如面波的主频集中在 5–15 Hz而深层反射信号同样落在这一频段50 Hz 工业干扰也会和部分有效波重叠。传统带通滤波只能按频率一刀切切掉噪声的同时必然损失同频段的有效信号。小波变换的不同之处在于它同时保留了频率和位置信息一次分解把原始道集按频带拆成近似分量和若干细节分量每个细节分量还带有时间位置信息。噪声如果只在某个短时间内占据某个尺度就可以在那个尺度的时间段内单独压制而不影响其他时刻的有效波。对地震数据做小波去噪本质上是在每个分解尺度上估计噪声强度然后对该尺度的细节系数做非线性收缩最后用处理过的系数重构信号。这个过程对随机噪声、尖峰脉冲和部分相干噪声都有较好的压制效果。实际处理中工区不同、激发接收条件不同噪声类型也不同小波基和分解层数的选择会显著影响最终剖面质量。2.2 小波基怎么选db4、db8 与 sym8 的取舍小波基的选择直接决定分解系数的能量集中程度。地震数据通常更看重线性相位和光滑性sym8 和 db8 是处理地震数据时最常见的两个选择。下表列出几组常用小波基的对比方便在不同任务里快速做决定。小波基消失矩对称性紧支撑适用场景db44近似对称是常规随机噪声压制计算量小db88近似对称是同相轴较平滑的地震数据保幅较好sym88近似对称是兼顾对称性与平滑性工区默认首选coif510近似对称是对振幅保真要求较高的 AVO 道集bior4.44严格对称是可逆性好但正交性差不常用于去噪选择小波基时主要看消失矩和对称性。消失矩越高越能压制信号中的多项式趋势但支撑长度也随之变长边界效应更明显。对地震数据来说层面反射在时间域上是较光滑的脉冲状信号db8 或 sym8 在多数情况下能在光滑度和定位精度之间取得较好的平衡。如果数据采样率较低或深层信号能量弱可以适当降到 db4避免小波变换本身引入过多振荡。2.2.1 分解层数与采样率的关系分解层数的确定和采样率直接相关。每分解一层频带对半划分第 j 层细节分量对应的频带约为 [fs/2^(j1), fs/2^j]近似分量对应 [0, fs/2^(j1)]。以 1 ms 采样fs1000 Hz为例第 1 层细节对应 250–500 Hz第 2 层对应 125–250 Hz第 3 层对应 62.5–125 Hz。地震有效波的频率范围通常在 5–80 Hz所以分解到第 4 层31.25–62.5 Hz或第 5 层15.625–31.25 Hz就已经覆盖主要有效频段。一般经验是让最低的细节分量下边界接近有效波最低频率的一半这样既不会把有效波的低频成分拆掉又能把高于有效频带的噪声单独分出去。2.3 阈值函数中的软阈值、硬阈值与折中方案拿到细节系数后需要先做阈值收缩。硬阈值把绝对值小于阈值的系数置零大于阈值的保持不变。它的优点是保幅能力好缺点是系数在阈值附近不连续重构信号容易在局部出现振荡也就是常说的伪吉布斯现象。软阈值把系数往零方向收缩一个阈值连续性更好压制噪声更彻底但会系统性低估大系数导致反射振幅偏弱。地震数据去噪最怕振幅失真所以很多处理流程选择折中方案。常见做法是使用 garrote 阈值函数系数绝对值大于阈值时按 x - threshold²/x 收缩既保持连续性又比软阈值更接近原始系数。对后续要做 AVO 分析的数据这个差异值得关注逐道输出振幅对比是必要的。3. 用 Python 在本地跑通地震数据小波去噪的最小流程3.1 准备数据从 SEGY 读出单道地震记录地震数据最常用的存储格式是 SEG-Y处理前先把道数据读进内存转成 numpy 数组。用 ObsPy 的函数读取最省事以下代码读取一个 SEG-Y 文件并取出前 10 道验证波形形状。import numpy as np import matplotlib.pyplot as plt from obspy import read st read(raw_data.sgy, formatSEGY) trace st[0].data.astype(np.float32) print(每道采样点数:, len(trace)) print(采样率(Hz):, st[0].stats.sampling_rate) plt.plot(trace[:500]) plt.title(First 500 samples of raw trace) plt.show()这里st[0]取出第一道地震记录sampling_rate通常由 SEG-Y 二进制头中的采样间隔字段解析而来。绘图时只取前 500 个采样点方便先看一眼噪声水平和有效波的相对位置。若数据是道集形式后续直接沿 axis1 逐道处理即可。3.2 用 PyWavelets 对单道数据做小波分解与阈值去噪去噪核心代码用 PyWavelets 完成。先pywt.wavedec做 5 层分解再用pywt.threshold对各层细节系数做阈值收缩最后pywt.waverec重构。完整代码如下。import pywt def denoise_trace(trace, waveletsym8, level5, modesoft): 对单道地震数据做小波去噪 coeffs pywt.wavedec(trace, wavelet, levellevel, modeperiodization) sigma np.median(np.abs(coeffs[-1])) / 0.6745 threshold sigma * np.sqrt(2 * np.log(len(trace))) coeffs_th [coeffs[0]] for c in coeffs[1:]: c_th pywt.threshold(c, threshold, modemode) coeffs_th.append(c_th) return pywt.waverec(coeffs_th, wavelet, modeperiodization) denoised denoise_trace(trace, waveletsym8, level5, modegarrote)参数modeperiodization表示信号边界按周期延拓能避免标准延拓带来的边界振荡。sigma用最高阶细节系数绝对值的 MEDIAN 估计噪声标准差除以 0.6745 是经验常数对应正态分布的中位数与标准差关系。threshold采用 Donoho 的通用阈值公式噪声越强、采样点越多阈值越高。pywt.threshold的mode支持soft、hard、garrote地震数据建议直接用garrote。3.3 批量处理整个道集并输出去噪结果实际工区里一个炮集通常有数百道逐道调用上面的函数即可。读入全部道数据后沿道方向循环处理。data np.stack([tr.data for tr in st]).astype(np.float32) data_den np.zeros_like(data) for i in range(data.shape[0]): data_den[i] denoise_trace(data[i], waveletsym8, level5, modegarrote) from obspy import Stream, Trace new_traces [] for i in range(data_den.shape[0]): tr st[i].copy() tr.data data_den[i] new_traces.append(tr) Stream(new_traces).write(denoised_data.sgy, formatSEGY)批量处理时要注意道头信息的保留特别是炮点坐标、道号等关键信息。上面代码用st[i].copy()复制原始 Trace 对象只替换数据部分可以最大程度保留 SEG-Y 头字段。写出的新文件可以直接进后续速度分析或叠加流程。4. 阈值策略与地震道集批量去噪的工程化4.1 全局阈值与逐尺度阈值的差异全局阈值方法计算简单但实际地震记录的噪声水平在不同尺度上并不一致。浅层高频噪声多深层以低频噪声为主单一阈值会顾此失彼阈值定高了深层有效细节系数被压制定低了浅层高频随机噪声留存过多。更稳健的做法是逐尺度估计噪声方差对各层细节系数单独设阈值。最高层细节分量几乎全是噪声可以直接用它的中位绝对偏差来估计全局噪声水平但对较低层要用该层自身的中位绝对偏差估计局部噪声。def denoise_trace_per_level(trace, waveletsym8, level5): coeffs pywt.wavedec(trace, wavelet, levellevel, modeperiodization) coeffs_th [coeffs[0]] for c in coeffs[1:]: sigma_j np.median(np.abs(c)) / 0.6745 th_j sigma_j * np.sqrt(2 * np.log(len(trace))) coeffs_th.append(pywt.threshold(c, th_j, modegarrote)) return pywt.waverec(coeffs_th, wavelet, modeperiodization)逐尺度阈值处理后低频段噪声压制相对温和有效波低频成分保留更完整。处理深层弱信号时要格外关注这一点因为深层信号本身能量弱全局阈值容易把这些微弱振幅连同事域上的高频抖动一起抹掉。4.2 大规模批量处理时的文件组织与参数管理把去噪流程放进正式处理流程时靠命令行参数管理比改脚本要高效得多。一份用于批量处理多个 SEG-Y 文件的脚本目录组织建议按以下结构拆分raw/放原始数据denoised/放输出params.json集中管理小波基、分解层数和阈值模式。处理脚本读 JSON 参数不支持逐文件改代码。python denoise_batch.py \ --input raw/line01.sgy \ --output denoised/line01_den.sgy \ --wavelet sym8 \ --level 5 \ --threshold-mode garrote \ --segypath data/*.sgy参数的意义分别对应小波基类型、分解层数、阈值函数模式以及通配符匹配一批文件。把level从 4 调到 6效果差异会很明显建议先在一个炮集上测试再批量跑工区。工区覆盖次数高、叠加道数多的数据去噪参数可以保守一些避免过度处理损伤有效信号覆盖次数低的数据则需要更强的去噪参数但阈值上限以不出现明显振幅凹槽为准。4.3 噪声水平估计的稳健性地震记录里的野值尖峰噪声会严重干扰噪声方差估计。中位绝对偏差比标准差稳健符合地震道存在零星尖峰脉冲的实际分布因此优先使用中位绝对偏差而不是np.std。若数据中强振幅异常道较多建议先用中值滤波剔除单道尖峰再做小波分解否则最高层细节系数会被少数大值主导导致阈值虚高。5. 去噪效果评估与常见参数误区5.1 两组指标SNR 提升与波形相关系数判断去噪好坏不能只看剖面“干净不干净”。量化指标要落到信噪比提升和有效波保真度两个维度。对合成记录或已知子波的模型数据可以直接计算去噪前后的 SNR。def compute_snr(clean, noisy): noise clean - noisy snr 10 * np.log10(np.sum(clean**2) / np.sum(noise**2) 1e-12) return snr def compute_cc(clean, denoised): a clean - np.mean(clean) b denoised - np.mean(denoised) return np.sum(a*b) / (np.sqrt(np.sum(a*a) * np.sum(b*b)) 1e-12)compute_snr以未加噪的干净数据为参考衡量去噪后残差的大小。compute_cc计算干净信号与去噪信号的相关系数越接近 1 说明波形保持得越好。实际工区没有干净参照时可以取未受噪声干扰的某一时窗计算去噪前后波形的相关系数作为替代指标。每日处理速报里建议同时列出这两个指标单看 SNR 会掩盖振幅被削平的问题。5.2 边界效应与 Gibbs 振荡的规避小波变换的边界处理是生产中最常被忽略的环节。使用默认的对称延拓时信号两端会被重复镜像重构后头部几十个毫秒常出现异常振幅压制了浅层初至的真实形态。改用modeperiodization可减轻边界问题但要求信号长度为 2 的整数次幂内部会自动截断到合适长度。更稳妥的办法是对每道信号两端各延拓一段零值或线性趋势处理完后切除延拓段。注意切除长度不能小于小波基支撑长度的一半否则延拓效果会打折扣。5.3 “细节系数全置零”是最常见的错误用法很多刚接触地震数据小波处理的人以为去噪就是把所有细节系数清零只保留近似分量。这个做法相当于做了一次极端低通滤波有效波中含高频成分的波形细节会全部消失同相轴变圆滑但分辨率明显下降。小波去噪的精髓是“保留大于阈值的系数”而不是丢弃所有细节系数。处理时用直方图查看各层系数分布凡是统计分布明显偏离高斯、出现长尾巴的细节层都保有大量有效信号成分不能一刀切置零。6. 用子带能量占比自适应确定分解层数的一个小技巧分解层数选多少最合适在实际生产中经常靠试。这里给出一个替代做法对单道数据逐层计算细节分量的能量占比找到能量占比较高的频带边界依此确定最大分解层数。具体实现是对数据做 8 层分解计算每层细节系数的均方根能量找出能量占比低于 1% 的最低层把该层序号减 1 作为实际分解层数。def find_adaptive_level(trace, waveletsym8, max_level8): coeffs pywt.wavedec(trace, wavelet, levelmax_level, modeperiodization) energies [] for c in coeffs[1:]: energies.append(np.sqrt(np.mean(c**2))) energies np.array(energies) ratio energies / (np.sum(energies) 1e-12) for i, r in enumerate(ratio): if r 0.01: return max(1, i) return max_level这个函数从最高频的细节层往低频方向扫描找到第一个能量占比不足 1% 的层。该层以上的高频分量对有效信号贡献极小主要成分是噪声。选择该层序号作为分解层数既保证有效频带内的细节被保留又避免把纯噪声层纳入阈值处理引入不必要的计算。用这个方法处理实际工区 1 ms 采样的数据时结果通常在 4 到 5 层之间跳变。对比同一炮集在固定 5 层与自适应层数下的去噪输出多数情况下两者差异很小但自适应方法在遇到深部弱反射数据时能自动减少分解层数避免把弱有效信号的能量在过深的分解中逐层摊薄。可以把这段逻辑封装成一个小工具在批量处理前先对每个炮集抽几道做层数探测输出推荐的层数范围再人工确认一次即可。这样做既保留了人工质控的环节又不用在每批数据上都重复试参。本文还有配套的精品资源点击获取 
长河君 · 2026/9/15 3:16:45 · 阅读 1.2k 地震数据小波去噪实战：从原理到Python实现
简介针对地震资料处理中随机噪声干扰严重、弱信号易被淹没等问题这份约8KB的MATLAB程序包以小波分析与D-S证据理论融合为技术主线面向地震勘探、信号处理与地球物理反演方向的科研人员和工程师提供了一套可用于实际资料处理的去噪参考实现。压缩包内仅含1个m脚本但功能覆盖了数据读取、小波基选择与多尺度分解、软硬阈值去噪、D-S证据融合、井曲线辅助校正、小波重构以及去噪前后剖面对比等完整流程可帮助使用者系统理解从噪声分离到多源信息整合的各个环节并掌握MATLAB中小波工具箱与证据理论的结合方式。资源体积精简仅8KB已吸引215人学习下载适合有一定MATLAB基础的学习者直接运行验证、替换自己的地震数据测试或在此基础上按需改造阈值策略与融合参数针对复杂地质条件下的噪声压制问题进行针对性优化是开展地震数据去噪研究的一手脚本参考。1. 地震数据去噪为什么绕不开小波拿到一批含噪地震记录时最先看到的往往是有效反射被面波、工业电干扰和随机噪声糊成一片。做过几年处理的人的第一反应是带通滤波但滤波后剖面变“干净”了同相轴也变“肉”了高频细节和弱反射一起被削掉。这正是地震数据去噪和小波去噪这类工具在工区处理流程里始终占有一席之地的原因地震信号本身是非平稳的有效信号与噪声在频带上大面积重叠单纯做频率切除解决不了问题。小波变换的时频局部化特性能把信号和噪声按不同尺度拆开再针对每个尺度单独压制噪声从而在保幅和去噪之间多出一个可调的维度。pengbeng_v82.zip 这一类在圈内流传的处理脚本包核心就是把“小波分解 → 阈值处理 → 重构”这套流程封装成可直接运行的模块省去自己从头写边界处理和系数选择的功夫。这篇文章就顺着这条技术线把原理、参数、可复现代码和实际踩过的坑一次讲完。2. 小波去噪用于地震数据的基本原理与选型2.1 为什么带通滤波压不掉与有效波同频的噪声地震记录中的随机噪声并不总是白噪很多时候它的频谱和有效反射波的频谱高度重合。比如面波的主频集中在 5–15 Hz而深层反射信号同样落在这一频段50 Hz 工业干扰也会和部分有效波重叠。传统带通滤波只能按频率一刀切切掉噪声的同时必然损失同频段的有效信号。小波变换的不同之处在于它同时保留了频率和位置信息一次分解把原始道集按频带拆成近似分量和若干细节分量每个细节分量还带有时间位置信息。噪声如果只在某个短时间内占据某个尺度就可以在那个尺度的时间段内单独压制而不影响其他时刻的有效波。对地震数据做小波去噪本质上是在每个分解尺度上估计噪声强度然后对该尺度的细节系数做非线性收缩最后用处理过的系数重构信号。这个过程对随机噪声、尖峰脉冲和部分相干噪声都有较好的压制效果。实际处理中工区不同、激发接收条件不同噪声类型也不同小波基和分解层数的选择会显著影响最终剖面质量。2.2 小波基怎么选db4、db8 与 sym8 的取舍小波基的选择直接决定分解系数的能量集中程度。地震数据通常更看重线性相位和光滑性sym8 和 db8 是处理地震数据时最常见的两个选择。下表列出几组常用小波基的对比方便在不同任务里快速做决定。小波基消失矩对称性紧支撑适用场景db44近似对称是常规随机噪声压制计算量小db88近似对称是同相轴较平滑的地震数据保幅较好sym88近似对称是兼顾对称性与平滑性工区默认首选coif510近似对称是对振幅保真要求较高的 AVO 道集bior4.44严格对称是可逆性好但正交性差不常用于去噪选择小波基时主要看消失矩和对称性。消失矩越高越能压制信号中的多项式趋势但支撑长度也随之变长边界效应更明显。对地震数据来说层面反射在时间域上是较光滑的脉冲状信号db8 或 sym8 在多数情况下能在光滑度和定位精度之间取得较好的平衡。如果数据采样率较低或深层信号能量弱可以适当降到 db4避免小波变换本身引入过多振荡。2.2.1 分解层数与采样率的关系分解层数的确定和采样率直接相关。每分解一层频带对半划分第 j 层细节分量对应的频带约为 [fs/2^(j1), fs/2^j]近似分量对应 [0, fs/2^(j1)]。以 1 ms 采样fs1000 Hz为例第 1 层细节对应 250–500 Hz第 2 层对应 125–250 Hz第 3 层对应 62.5–125 Hz。地震有效波的频率范围通常在 5–80 Hz所以分解到第 4 层31.25–62.5 Hz或第 5 层15.625–31.25 Hz就已经覆盖主要有效频段。一般经验是让最低的细节分量下边界接近有效波最低频率的一半这样既不会把有效波的低频成分拆掉又能把高于有效频带的噪声单独分出去。2.3 阈值函数中的软阈值、硬阈值与折中方案拿到细节系数后需要先做阈值收缩。硬阈值把绝对值小于阈值的系数置零大于阈值的保持不变。它的优点是保幅能力好缺点是系数在阈值附近不连续重构信号容易在局部出现振荡也就是常说的伪吉布斯现象。软阈值把系数往零方向收缩一个阈值连续性更好压制噪声更彻底但会系统性低估大系数导致反射振幅偏弱。地震数据去噪最怕振幅失真所以很多处理流程选择折中方案。常见做法是使用 garrote 阈值函数系数绝对值大于阈值时按 x - threshold²/x 收缩既保持连续性又比软阈值更接近原始系数。对后续要做 AVO 分析的数据这个差异值得关注逐道输出振幅对比是必要的。3. 用 Python 在本地跑通地震数据小波去噪的最小流程3.1 准备数据从 SEGY 读出单道地震记录地震数据最常用的存储格式是 SEG-Y处理前先把道数据读进内存转成 numpy 数组。用 ObsPy 的函数读取最省事以下代码读取一个 SEG-Y 文件并取出前 10 道验证波形形状。import numpy as np import matplotlib.pyplot as plt from obspy import read st read(raw_data.sgy, formatSEGY) trace st[0].data.astype(np.float32) print(每道采样点数:, len(trace)) print(采样率(Hz):, st[0].stats.sampling_rate) plt.plot(trace[:500]) plt.title(First 500 samples of raw trace) plt.show()这里st[0]取出第一道地震记录sampling_rate通常由 SEG-Y 二进制头中的采样间隔字段解析而来。绘图时只取前 500 个采样点方便先看一眼噪声水平和有效波的相对位置。若数据是道集形式后续直接沿 axis1 逐道处理即可。3.2 用 PyWavelets 对单道数据做小波分解与阈值去噪去噪核心代码用 PyWavelets 完成。先pywt.wavedec做 5 层分解再用pywt.threshold对各层细节系数做阈值收缩最后pywt.waverec重构。完整代码如下。import pywt def denoise_trace(trace, waveletsym8, level5, modesoft): 对单道地震数据做小波去噪 coeffs pywt.wavedec(trace, wavelet, levellevel, modeperiodization) sigma np.median(np.abs(coeffs[-1])) / 0.6745 threshold sigma * np.sqrt(2 * np.log(len(trace))) coeffs_th [coeffs[0]] for c in coeffs[1:]: c_th pywt.threshold(c, threshold, modemode) coeffs_th.append(c_th) return pywt.waverec(coeffs_th, wavelet, modeperiodization) denoised denoise_trace(trace, waveletsym8, level5, modegarrote)参数modeperiodization表示信号边界按周期延拓能避免标准延拓带来的边界振荡。sigma用最高阶细节系数绝对值的 MEDIAN 估计噪声标准差除以 0.6745 是经验常数对应正态分布的中位数与标准差关系。threshold采用 Donoho 的通用阈值公式噪声越强、采样点越多阈值越高。pywt.threshold的mode支持soft、hard、garrote地震数据建议直接用garrote。3.3 批量处理整个道集并输出去噪结果实际工区里一个炮集通常有数百道逐道调用上面的函数即可。读入全部道数据后沿道方向循环处理。data np.stack([tr.data for tr in st]).astype(np.float32) data_den np.zeros_like(data) for i in range(data.shape[0]): data_den[i] denoise_trace(data[i], waveletsym8, level5, modegarrote) from obspy import Stream, Trace new_traces [] for i in range(data_den.shape[0]): tr st[i].copy() tr.data data_den[i] new_traces.append(tr) Stream(new_traces).write(denoised_data.sgy, formatSEGY)批量处理时要注意道头信息的保留特别是炮点坐标、道号等关键信息。上面代码用st[i].copy()复制原始 Trace 对象只替换数据部分可以最大程度保留 SEG-Y 头字段。写出的新文件可以直接进后续速度分析或叠加流程。4. 阈值策略与地震道集批量去噪的工程化4.1 全局阈值与逐尺度阈值的差异全局阈值方法计算简单但实际地震记录的噪声水平在不同尺度上并不一致。浅层高频噪声多深层以低频噪声为主单一阈值会顾此失彼阈值定高了深层有效细节系数被压制定低了浅层高频随机噪声留存过多。更稳健的做法是逐尺度估计噪声方差对各层细节系数单独设阈值。最高层细节分量几乎全是噪声可以直接用它的中位绝对偏差来估计全局噪声水平但对较低层要用该层自身的中位绝对偏差估计局部噪声。def denoise_trace_per_level(trace, waveletsym8, level5): coeffs pywt.wavedec(trace, wavelet, levellevel, modeperiodization) coeffs_th [coeffs[0]] for c in coeffs[1:]: sigma_j np.median(np.abs(c)) / 0.6745 th_j sigma_j * np.sqrt(2 * np.log(len(trace))) coeffs_th.append(pywt.threshold(c, th_j, modegarrote)) return pywt.waverec(coeffs_th, wavelet, modeperiodization)逐尺度阈值处理后低频段噪声压制相对温和有效波低频成分保留更完整。处理深层弱信号时要格外关注这一点因为深层信号本身能量弱全局阈值容易把这些微弱振幅连同事域上的高频抖动一起抹掉。4.2 大规模批量处理时的文件组织与参数管理把去噪流程放进正式处理流程时靠命令行参数管理比改脚本要高效得多。一份用于批量处理多个 SEG-Y 文件的脚本目录组织建议按以下结构拆分raw/放原始数据denoised/放输出params.json集中管理小波基、分解层数和阈值模式。处理脚本读 JSON 参数不支持逐文件改代码。python denoise_batch.py \ --input raw/line01.sgy \ --output denoised/line01_den.sgy \ --wavelet sym8 \ --level 5 \ --threshold-mode garrote \ --segypath data/*.sgy参数的意义分别对应小波基类型、分解层数、阈值函数模式以及通配符匹配一批文件。把level从 4 调到 6效果差异会很明显建议先在一个炮集上测试再批量跑工区。工区覆盖次数高、叠加道数多的数据去噪参数可以保守一些避免过度处理损伤有效信号覆盖次数低的数据则需要更强的去噪参数但阈值上限以不出现明显振幅凹槽为准。4.3 噪声水平估计的稳健性地震记录里的野值尖峰噪声会严重干扰噪声方差估计。中位绝对偏差比标准差稳健符合地震道存在零星尖峰脉冲的实际分布因此优先使用中位绝对偏差而不是np.std。若数据中强振幅异常道较多建议先用中值滤波剔除单道尖峰再做小波分解否则最高层细节系数会被少数大值主导导致阈值虚高。5. 去噪效果评估与常见参数误区5.1 两组指标SNR 提升与波形相关系数判断去噪好坏不能只看剖面“干净不干净”。量化指标要落到信噪比提升和有效波保真度两个维度。对合成记录或已知子波的模型数据可以直接计算去噪前后的 SNR。def compute_snr(clean, noisy): noise clean - noisy snr 10 * np.log10(np.sum(clean**2) / np.sum(noise**2) 1e-12) return snr def compute_cc(clean, denoised): a clean - np.mean(clean) b denoised - np.mean(denoised) return np.sum(a*b) / (np.sqrt(np.sum(a*a) * np.sum(b*b)) 1e-12)compute_snr以未加噪的干净数据为参考衡量去噪后残差的大小。compute_cc计算干净信号与去噪信号的相关系数越接近 1 说明波形保持得越好。实际工区没有干净参照时可以取未受噪声干扰的某一时窗计算去噪前后波形的相关系数作为替代指标。每日处理速报里建议同时列出这两个指标单看 SNR 会掩盖振幅被削平的问题。5.2 边界效应与 Gibbs 振荡的规避小波变换的边界处理是生产中最常被忽略的环节。使用默认的对称延拓时信号两端会被重复镜像重构后头部几十个毫秒常出现异常振幅压制了浅层初至的真实形态。改用modeperiodization可减轻边界问题但要求信号长度为 2 的整数次幂内部会自动截断到合适长度。更稳妥的办法是对每道信号两端各延拓一段零值或线性趋势处理完后切除延拓段。注意切除长度不能小于小波基支撑长度的一半否则延拓效果会打折扣。5.3 “细节系数全置零”是最常见的错误用法很多刚接触地震数据小波处理的人以为去噪就是把所有细节系数清零只保留近似分量。这个做法相当于做了一次极端低通滤波有效波中含高频成分的波形细节会全部消失同相轴变圆滑但分辨率明显下降。小波去噪的精髓是“保留大于阈值的系数”而不是丢弃所有细节系数。处理时用直方图查看各层系数分布凡是统计分布明显偏离高斯、出现长尾巴的细节层都保有大量有效信号成分不能一刀切置零。6. 用子带能量占比自适应确定分解层数的一个小技巧分解层数选多少最合适在实际生产中经常靠试。这里给出一个替代做法对单道数据逐层计算细节分量的能量占比找到能量占比较高的频带边界依此确定最大分解层数。具体实现是对数据做 8 层分解计算每层细节系数的均方根能量找出能量占比低于 1% 的最低层把该层序号减 1 作为实际分解层数。def find_adaptive_level(trace, waveletsym8, max_level8): coeffs pywt.wavedec(trace, wavelet, levelmax_level, modeperiodization) energies [] for c in coeffs[1:]: energies.append(np.sqrt(np.mean(c**2))) energies np.array(energies) ratio energies / (np.sum(energies) 1e-12) for i, r in enumerate(ratio): if r 0.01: return max(1, i) return max_level这个函数从最高频的细节层往低频方向扫描找到第一个能量占比不足 1% 的层。该层以上的高频分量对有效信号贡献极小主要成分是噪声。选择该层序号作为分解层数既保证有效频带内的细节被保留又避免把纯噪声层纳入阈值处理引入不必要的计算。用这个方法处理实际工区 1 ms 采样的数据时结果通常在 4 到 5 层之间跳变。对比同一炮集在固定 5 层与自适应层数下的去噪输出多数情况下两者差异很小但自适应方法在遇到深部弱反射数据时能自动减少分解层数避免把弱有效信号的能量在过深的分解中逐层摊薄。可以把这段逻辑封装成一个小工具在批量处理前先对每个炮集抽几道做层数探测输出推荐的层数范围再人工确认一次即可。这样做既保留了人工质控的环节又不用在每批数据上都重复试参。本文还有配套的精品资源点击获取