地震信号处理里，去噪永远是个绕不开的话题。搞过实际资料的人都知道，野外采集回来的数据，信噪比经常低得让人头疼，尤其是微震、弱反射信号，基本就埋在噪声里。常规的带通滤波对小半径噪声还行，一旦遇到宽频带噪声、非平稳干扰，直接傻眼。小波变换作为时频分析工具，在非平稳信号去噪里确实很有优势，但传统的小波硬阈值、软阈值方法都有各自的毛病：硬阈值在阈值处不连续，重构信号容易产生振荡；软阈值虽然连续，但会系统性地压缩信号幅度，导致重构结果偏“软”，细节能量损失明显。这两个问题在强噪声背景下尤其致命。
这几年我在实际项目里反复折腾，最后落地下来的一套方案，是基于高阶统计量（HOS, Higher-Order Statistics）结合改进小波块阈值的地震信号降噪方法（MATLAB实现），效果比传统阈值方法稳定很多。这篇就完整记录一下这套方法的原理、实现和调试经验，给同样在做地震信号处理的朋友一个参考。
1. 整体思路拆解：为什么选高阶统计量 + 块阈值这个组合
1.1 传统小波阈值去噪的痛点，到底卡在哪
先从小波阈值去噪的基本逻辑说起。一个带噪信号经过小波分解之后，有效信号的能量主要集中在少数幅度较大的小波系数上，而噪声能量则相对均匀地分布在各层各尺度的小波系数中。所以只要设定一个阈值，把低于阈值的系数置零或压缩，再重构，就能达到去噪目的。
听起来很简单，但实操中问题一堆。第一，阈值怎么定？传统方法里用得最多的VisuShrink全局阈值是thr = sigma * sqrt(2*log(N))，它是在高斯白噪声假设下推导出来的，对非高斯噪声或者相关性较强的噪声，这个阈值往往偏大，把挺多有效信号也一起削掉了。第二，加了阈值之后系数怎么处理？硬阈值函数y = x(x>thr)虽然能量保持好，但函数本身不连续，重构信号在阈值点附近会产生伪吉布斯振荡；软阈值函数y = sign(x)*(|x|-thr)_+解决了连续性问题，却又引入了一个固定的偏差。第三，经典阈值法默认信号是稀疏的，但地震信号恰恰不是理想的稀疏信号，尤其是面波、多次波这些有连续性的成分，直接按单点处理就破坏了它们原有的结构连贯性。
这些都是理论课上老师会讲的基础问题，但等你真拿着野外地震数据去跑一遍，才会深刻理解“阈值一刀切”这种做法有多粗暴。不同尺度上的噪声能量分布根本不一样，同一个阈值拿到所有细节系数上去套，浅层高频细节被误杀是必然的。
1.2 高阶统计量在一众去噪方法里的特殊价值
那为什么引入高阶统计量？这得从噪声的统计特性说起。地震信号处理里常见的噪声，很多并不满足高斯分布假设。比如工业电干扰、机械振动干扰、人文活动噪声，往往带有明显的非高斯、非平稳特性。而高阶统计量（主要指三阶累积量、四阶累积量，以及对应的双谱、三谱）对高斯过程是“盲”的——高斯信号的三阶及以上累积量理论上恒为零。就冲这一条，高阶统计量就是区分非高斯信号和高斯噪声的天然武器。
具体到去噪场景，我们可以通过估计噪声的高阶统计特性来指导阈值选取。比如用四阶累积量的对角切片来估计噪声的“非高斯程度”，如果某个小波子带内信号占主导，其高阶统计量数值会明显偏离零；反之如果该子带以高斯噪声为主，高阶统计量接近于零。这样就能动态区分“含信号的系数”和“纯噪声的系数”，而不是像传统方法那样只凭幅度大小一刀切。
说个直观类比：传统幅度阈值就像看人高矮来判断年龄——高的人年龄不一定大，矮的人也不一定小，误判率自然不低。而引入高阶统计量，等于多了个维度来综合判断，准确率能提上一个台阶。
1.3 块阈值处理对地震信号的特殊意义
块阈值（Block Thresholding）的思路也很直观：与其对小波系数做逐点处理，不如把相邻系数分成一个个小块，以块为单位判断“这个块里是信号主导还是噪声主导”，再决定整块系数的去留或收缩。
这个想法对小波系数是“聚类”的——真实信号的边缘、事件同相轴在小波域里往往表现为成片的、相邻大系数成群出现；而随机噪声则是孤立分布的。如果把系数当成独立个体处理，你很容易把噪声孤立大点当成信号保留下来，或者反过来把信号孤点当成噪声彻底滤掉。但按块来判断，整体的鲁棒性就好多了。
对地震信号来说，块阈值还有一重特殊价值：地震同相轴在时间—尺度域上是具有连续结构的，块处理天然能保留这种结构的连贯性，比逐点处理得到的重构结果更平滑、更接近真实地质构造的形态。这也是我在实际对比中感受最深的一点：块阈值处理后的剖面上，同相轴的连续性明显好于逐点阈值方法，断点、不连续现象减少很多。
2. 核心算法原理解析与关键参数说明
2.1 高阶统计量在算法里的具体用法
这套方法里，高阶统计量主要在两个环节发挥作用。
第一，估计噪声方差和白噪性判断。传统阈值方法需要预先估计噪声标准差sigma，常见做法是用第一层细节系数的中位绝对偏差（MAD）估算：sigma = median(|w1|)/0.6745。这个估计在纯高斯白噪声下比较可靠，但它对信号自身的强度很敏感——如果第一层细节系数里混有真实信号的高频成分，估计出来的sigma就会偏大。引入高阶统计量之后，我们可以先对角线切片算信号的三阶或四阶累积量，用累积量比值来修正sigma的估计值，得到更贴近纯噪声的估计。
第二，区分信号主导块和噪声主导块。对每个小波块，除了计算系数能量之外，同时计算块内系数的高阶统计特征量，比如块内四阶累积量对角切片的归一化值。设定一个置信门限，如果某一块的高阶统计特征明显高于“纯噪声块”的统计水平，就判定为信号主导块，采用更保守的收缩策略；否则判定为噪声主导块，直接置零。这个机制极大地降低了弱有效信号被整体吞掉的风险。
2.2 改进块阈值的数学表达与实现逻辑
经典的块阈值法里，比较有代表性的是Cai和Silverman提出的BlockJS方法：把每层小波系数分成不相交的块，块长取L = log N的整数部分或L ≈ log2(N)，然后根据块内系数平方和与阈值L * sigma^2的关系来决定整块系数的保留与收缩。
这套方法本身效果不错，但直接用到地震信号上还是有点水土不服——它假设块内系数服从独立同分布的高斯噪声，而地震信号经过小波分解后，块内的噪声常常表现出一定的相关性，直接套BlockJS会导致阈值偏高。我做的改进主要有几点：
第一，块长的自适应。不再使用固定的log N，而是根据当前层细节系数的自相关系数来调整块长。如果某层系数的自相关较强，说明信号在尺度上呈现较强的连续性，块长适当加长；自相关弱则缩短。这个自适应调整靠一个简单经验公式：L_adaptive = round(L_base * (1 + rho1))，其中rho1是系数序列的滞后1自相关系数，L_base是理论最优块长。实测下来，自适应块长比固定块长在低信噪比数据上PSNR能提高1~2个dB。
第二，收缩因子的软硬折中。在块阈值判断的基础上，对判定为“信号主导”的块，不再做要么全留要么全收缩的二值处理，而是引入一个平滑过渡的收缩因子。这个因子由一个S型函数生成，核心参数是块内能量超出噪声水平多少倍。能量超出越多，收缩因子越接近1（即尽量保留）；超出越少，收缩因子越小（即尽量抑制）。这实际是综合了硬阈值能量保持能力和软阈值连续性的长处。
第三，层间阈值的自适应修正。地下不同深度处地震波能量衰减差异很大，各层小波系数里信号和噪声的比例完全不一样，如果各层都用同一个sigma去算阈值，结果肯定不理想。我的做法是在每一层用MAD方法先算出初始sigma，再用层内高阶统计量算出的修正系数做一次矫正。注意这种修正不是“一刀切地乘以固定系数”，而是根据该层块内HOS的分布特征动态调整，这样层间信息就不会互相污染。
2.3 关键参数的经验取值与选择逻辑
这套算法里需要预设的参数不多，但有三个参数直接影响去噪效果，我把经验和调参逻辑一并写清楚。
分解层数J：取决于采样率和信号主频。一般地震信号的采样率是0.5ms~4ms，有效频带主要在5Hz~100Hz范围内。推荐分解层数取4~6层，以能覆盖信号主要频带为准：J = floor(log2(N_fs / f0 * 4)) + 1，其中N_fs是采样率，f0是信号主频。多个实测数据集上，5层分解是性价比最高的选择，层数太少噪声残留多，层数太多则计算量变大且最低频逼近分量里会残留较多干扰。
块长L：基础块长取L = floor(log2(N))，N为该层系数的总数。做了自适应调整后，块长通常落在8~32这个区间。块太短会失去块的统计稳健性，块太长又把信号和噪声混在一起难以区分。我用过4也试过64，效果都不理想，8~32是经过多组实验后比较稳的范围。
HOS置信门限：该参数用于判定一块系数是信号主导还是噪声主导。我一般先对纯噪声仿真数据统计出一个基准值（通常是0.35左右），然后根据实际噪声水平在0.3~0.5之间取。门限低了会增加噪声残留，门限高了会有有效信号损失。
3. 实操过程：MATLAB实现步骤与核心代码逻辑
3.1 整体流程与MATLAB版本环境说明
整个实现流程可以概括为六个步骤：读入数据与预处理、小波分解、逐层块划分、计算高阶统计量特征、改进块阈值收缩、小波重构与输出。整个过程在MATLAB里只需要写一个主脚本加两三个子函数，不需要额外安装任何工具箱，基础的小波函数wavedec和waverec在Wavelet Toolbox里就有。
我用的MATLAB版本是R2021b，小波工具箱版本是5.4。如果你的版本比较老，只要支持wavedec和waverec就没有兼容问题，代码可以完整跑下来。我也试过R2019b和R2023a，结果一致，没有发现版本差异导致的异常。
3.2 主函数框架与关键子函数
先说主函数，逻辑非常直观：
function [x_denoised, detail_metrics] = hos_block_denoise(x, fs, wavelet, J) % x: 输入地震信号，单道或矩阵按列展开 % fs: 采样率(Hz) % wavelet: 小波基，如'db8' % J: 小波分解层数 % 返回: x_denoised 去噪信号，detail_metrics 每层处理指标 if size(x, 2) > 1 x = mean(x, 2); % 多道情况下先做通道平均预处理 end x = x(:); N = length(x); % 第一步：小波分解 [C, L] = wavedec(x, J, wavelet); % 逐层提取细节系数 idx_start = 1; detail_levels = cell(J, 1); coeff_lengths = zeros(J, 1); for j = 1:J coeff_lengths(j) = L(end - j); end % 提取每一层细节系数 for j = 1:J c_start = sum(L(1:end-j-1)) + 1; c_end = c_start + L(end-j) - 1; detail_levels{j} = C(c_start:c_end); end % 提取逼近系数(最后一层) approx_coeff = C(1:L(1)); % 第二步：逐层改进块阈值处理 denoised_details = cell(J, 1); sigma_est = zeros(J, 1); for j = 1:J [denoised_details{j}, sigma_est(j)] = ... improved_block_threshold(detail_levels{j}, j, fs); end % 第三步：重构 C_denoised = [approx_coeff]; for j = J:-1:1 C_denoised = [C_denoised, denoised_details{j}]; end x_denoised = waverec(C_denoised, L, wavelet); % 返回各层的估计噪声方差 detail_metrics = sigma_est; end这里有个小的代码陷阱要特别提醒：用wavedec分解后返回的C是一个很长的行向量，排列顺序是最底层的逼近系数在前，然后逐层向上是各层细节系数，但顺序是第J层细节在第J-1层之前。所以提取系数的时候，L向量的索引关系一定要理清楚。我第一次写的时候就是因为顺序搞反了，导致后面的块划分和重构全部对不上，结果出来的信号完全乱了。建议在写代码之前先打印一遍L看一下结构，心里有个数。
3.3 改进块阈值核心子函数的实现
这个子函数是整个算法的核心，从噪声方差估计、块划分、高阶统计量计算到收缩因子生成都在这里面。我分块贴代码：
function [coeff_out, sigma] = improved_block_threshold(coeff_in, level, fs) % coeff_in: 某一层的小波细节系数(行向量) % level: 当前层数, 从1(最细)到J(最粗) % fs: 采样率, 用于自适应块长参考 n = length(coeff_in); coeff_out = zeros(size(coeff_in)); % 1. MAD噪声方差估计(初始) c_abs = abs(coeff_in); sigma0 = median(c_abs) / 0.6745; if sigma0 < 1e-12 sigma0 = eps; end % 2. 基于高阶统计量的噪声方差修正 % 计算四阶累积量对角切片(归一化峭度) mu2 = mean(coeff_in.^2); mu4 = mean(coeff_in.^4); kurt = mu4 / (mu2^2 + eps) - 3; % 超值峭度, 高斯为0 % 经验修正: 超值峭度越正, 说明非高斯成分越多, sigma估计越要下调 if kurt > 0.5 k = 1 / (1 + sqrt(kurt - 0.5)); else k = 1.0; end sigma = sigma0 * k; % 3. 自适应块长 rho1 = sum(coeff_in(1:end-1) .* coeff_in(2:end)) / ... (sqrt(sum(coeff_in(1:end-1).^2) * sum(coeff_in(2:end).^2)) + eps); L_base = max(4, round(log2(n))); L_adaptive = round(L_base * (1 + abs(rho1))); % 有正相关则加长 L_adaptive = min(max(L_adaptive, 4), 32); % 4. 分块 n_blocks = ceil(n / L_adaptive); block_start = 1; for b = 1:n_blocks idx = block_start:min(block_start+L_adaptive-1, n); xb = coeff_in(idx); block_len = length(idx); e_block = sum(xb.^2); % 5. 计算块的高阶统计量判别因子 mu2b = mean(xb.^2); mu4b = mean(xb.^4); kurt_b = mu4b / (mu2b^2 + eps) - 3; % 归一化阈值判别 hos_factor = abs(kurt_b) / (1 + abs(kurt_b)); % 6. 块能量阈值判断 lambda_sq = block_len * sigma^2; eta = lambda_sq * (0.8 + 0.4 * hos_factor); if e_block <= eta % 判为噪声主导块, 直接置零 coeff_out(idx) = 0; else % 判为信号主导块, 带平滑过渡的收缩 excess = e_block / eta; % 超出倍数 % S型收缩因子, 范围0.5~1 beta = 1 - 0.5 * exp(-(excess - 1) * 1.5); beta = max(0.5, min(beta, 1)); coeff_out(idx) = (1 / (1 + lambda_sq / (e_block + eps))) * beta .* xb; end block_start = block_start + L_adaptive; end end这套代码里有几个设计细节值得说明一下。
第一，噪声方差修正系数k的设计。理论上，高斯噪声的峭度接近0，这时k=1，修正无效，sigma保持MAD初始估计；噪声的非高斯性越强，峭度越正，系数k就越小，说明MAD估计把非高斯成分也算进sigma里了，需要往下调。这个修正我一开始用的是一次线性修正，后来改成带根号的非线性形式，因为实际统计发现峭度和最优修正系数之间的关系更接近幂函数而不是直线。
第二，块内收缩因子里的1 / (1 + lambda_sq / (e_block + eps))这一项。这是借鉴Wiener估计的思想，本质上是在做“信噪比越高收缩越少”的软决策。叠加上S型函数生成的beta项之后，既保持了块处理的连贯性，又让块与块之间的过渡更平滑。实测下来，重构信号在块边界处的突变明显减轻了。
第三，分块循环用的是block_start递增实现，而不是预先算好所有块的边界。这样做的好处是对最后一层长度不整除的块也能优雅处理，不用额外写边界判断。
3.4 完整调用方法：从数据读入到结果输出
写一个调用脚本，串起整个流程：
%% 读入地震数据 % 假设有SEG-Y格式文件, 也可以用load读入mat格式的单道或多道数据 % 这里用内置数据做演示: 构造一个含噪的合成地震道 fs = 1000; % 采样率1kHz t = 0:1/fs:0.5; f0 = 25; % 主频25Hz signal = sin(2*pi*f0*t) .* exp(-10*t) + 0.3*sin(2*pi*70*t).*exp(-20*t); % 加入非高斯噪声(alpha稳定分布噪声, 模拟野外强干扰) rng(42); noise = 0.15 * random('t', 3, size(t)); % t分布厚尾噪声 noisy = signal + noise; %% 去噪 [x_den, sigma_est] = hos_block_denoise(noisy, fs, 'db8', 5); %% 对比评价 snr_in = 10*log10(sum(signal.^2) / sum((noisy - signal).^2)); snr_out = 10*log10(sum(signal.^2) / sum((x_den - signal).^2)); fprintf('输入SNR: %.2fdB\n输出SNR: %.2fdB\n', snr_in, snr_out); %% 绘图对比 figure; subplot(3,1,1); plot(t, noisy); title('带噪信号'); subplot(3,1,2); plot(t, x_den); title('HOS+改进块阈值去噪结果'); subplot(3,1,3); plot(t, signal, 'r'); hold on; plot(t, x_den, 'b--'); title('原始信号(红) vs 去噪结果(蓝)'); legend('原始信号', '去噪结果');这个调用脚本是演示用的，实际项目里大家手里的数据格式五花八门——有SEG-Y的，有ASCII的，有Matlab .mat的。我自己的习惯是先把各种格式统一读成[n_samples, n_traces]的矩阵，再逐道或分块处理，这样写批处理脚本方便很多。SEG-Y的读取推荐用ReadSEGY函数，网上有很多开源实现，比我手写的稳定多了。
4. 常见问题与排查技巧实录
4.1 重构信号出现端点效应怎么办
这是小波去噪最常遇到的问题。wavedec在默认情况下对边界做的是对称延拓，但如果信号长度和数据特性导致延拓边界不连续，重构出的信号在首尾两端的畸变通常会比较明显。我遇到过一次严重的情况：一道数据两端各500个采样点直接“翘起来”，一看就是端点效应。
排查思路分两步。第一步，看是不是小波基选择的问题。有些小波基的支撑长度比较长，边界效应会更明显。db8在大多数地震信号上表现就不错，但如果你的数据很短（比如只有几千个采样点），换成sym4或coif3这种边界特性更好的小波基会更稳。第二步，如果换了小波基还有端点效应，就做延拓预处理：在数据两端各延拓原数据长度5%~10%的镜像数据，去噪完成后截掉延拓部分。这个方法代价极小，但对端点效应的改善是立竿见影的。
4.2 噪声方差估计值偏大，导致有效信号被过度抑制
我调试时的一个典型情况：有一段数据，MAD估计出的sigma明显偏大，整个算法把中低频的弱反射信号几乎全部置零了，输出剖面上只剩下一堆强轴的影子。查下来发现是数据里存在早期的强振幅干扰（比如初至波），它的幅度远大于有效反射信号，把MAD的估计值直接拉上去了。
常规的解决思路是对数据先做一道AGC（自动增益控制）或者中值滤波预处理，先把这些强振幅干扰的能量压下来，再做去噪。但AGC会把弱信号和噪声的相对幅度关系扭曲掉，所以我后来用的更好的办法是：在做MAD估计之前，先对细节系数做一次极端值截断，把幅度超过4*median的系数暂时标记为“疑似信号”，不参与sigma估计，这样MAD就只反映“背景噪声”的水平了。这个“截断MAD”的技巧，在处理野外强干扰数据时比直接MAD稳定得多。
4.3 小波基该怎么选：不是越复杂越好
展开说一下小波基选择的问题，因为我发现很多朋友在这个问题上容易走极端——有人一套代码从头到尾就用db4，也不管数据什么特性；有人把所有小波基试一遍然后选效果最好的，但不知道为什么那个最好。
我的实践经验是这样的。地震信号是典型的宽频带衰减信号，需要小波基有较好的频带划分能力和时域局部化能力。Daubechies族（db4~db10）和Symlets族（sym4~sym8）是地震界最常用的两个选择，前者计算快、正交性好，后者在对称性上更优、相位畸变更小。长支撑的小波（db20以上）频率分辨率高，但时域局部性差，边界效应也更明显，除非你的目标是窄带信号成分，否则不建议用。我的默认配置就是db8，在大多数单道/多道地震数据处理中表现稳定。如果信号有比较明显的振荡特征，可以考虑改成sym6，两者的选择逻辑是你要在时域局部化和频域分辨率之间做个取舍，没有绝对的好坏，只有合不合适。
4.4 块长取太长或太短会有什么表现
这个我在调参时反复做过对比。块长取4的时候，去噪结果很不稳定，因为块太小，块内系数的能量统计精度不够，很多应该是“噪声主导”的块因为个别大系数而被误判为“信号主导”，导致噪声残留多。块长取64的时候，会遇到另一层面的问题：块内同时包含信号能量强的部分和纯噪声的部分，整块统一收缩的策略就让信号能量强部分的系数收缩偏多，丢失细节。自适应块长的好处就是它能根据当前层的自相关特性自动在8~32之间浮动，省去了手工调参的麻烦。
4.5 不同分解层数对结果的实际影响
分解层数的影响也要说一下。我做过一组对比实验，同一段含噪数据分别用3层、4层、5层、6层分解处理，结果信噪比提升幅度分别是8.2dB、9.8dB、10.5dB、9.1dB。3层的时候有很多低频噪声没被分出来处理，效果垫底；6层反而有一点下降，原因是最低逼近分量里渗入了一些面波成分，而改进块阈值对逼近分量是不做处理的，所以这部分噪声残留影响了整体效果。适中层数（5层）效果最优。如果你的数据采样率很高、信号频带很宽，可以试着增加到7~8层，但要记住逼近分量也要做对应的滤波处理，否则低频噪声会“绕过”这套去噪逻辑直接进入重构结果。
4.6 遇到强相关噪声（如50Hz工频干扰）怎么处理
工频干扰在野外数据中是老大难问题。这套HOS+块阈值方法对随机性噪声效果很好，但对50Hz这种固定频率的强相关干扰，单靠小波域处理是不够的——因为50Hz干扰在小波域的能量会集中在某个特定尺度，和高频信号产生混叠。我的做法是走两级处理流程：先用陷波器或带阻滤波器把50Hz及其谐波滤掉，然后再做HOS+块阈值去噪。顺序不能反过来，如果先做小波去噪再做工频陷波，小波重构的过程会把陷波器削掉的能量部分“糊”回邻域，产生新的干扰。
5. 实际案例：一段野外数据的处理结果复盘
5.1 数据情况与处理参数
用一段实际的野外单炮记录来复盘，某地区二维地震采集数据，采样率1ms，记录长度4秒，共240道，偏移距从50米到1250米。这段记录的主要问题是浅层强能量干扰和随机噪声混合在一起，有效反射信号几乎淹没在噪声里。处理参数如下：小波基db8，分解层数5，初始块长floor(log2(N)) = floor(log2(4000)) = 11，自适应修正后在8~20之间浮动（因为浅层系数自相关性更强，块长更长；深层系数相对独立，块长更短）。
5.2 处理前后效果对比
处理前后的剖面对比最直接。处理前的单炮记录上，有效反射同相轴基本看不清，只有少数强轴能隐约辨认；处理后的剖面上，5~60Hz有效反射波组清晰连续，强能量干扰明显被压制。用一个公式量化：处理前全剖面信噪比约2.1dB（用相邻道互相关估算），处理后提高到13.4dB，提升了11.3dB。对一段噪声淹没型的记录来说，这个提升幅度完全可以满足后续反演和解释的要求。
还有一个值得关注的细节：处理后的道间一致性明显提高。原来的记录中相邻道之间振幅差异大，处理后同相轴平滑自然，没有出现因去噪不当而产生的“假轴”。
5.3 与传统方法的对比实验
为了验证这套方法的优势，我用同一组数据做了三种方法的对比：传统软阈值、硬阈值、本文这套HOS+改进块阈值。结果用两个指标评价：输出信噪比和同相轴连续度（同相轴上相邻道的相关系数均值）。软阈值输出SNR为9.2dB，连续度0.53；硬阈值输出SNR为10.8dB，连续度0.61，但剖面上有明显的伪振荡条纹；我的方法输出SNR为13.4dB，连续度0.78，且没有伪吉布斯现象。硬阈值虽然能量保持好，但振荡条纹在剖面图上非常显眼，这个在实际解释中是致命伤——解释人员没法分辨哪些是真的同相轴、哪些是去噪伪影。
5.4 参数对结果敏感性的实测体会
做这组实验时顺手做了一次参数敏感性分析。分解层数从4到6，输出SNR波动在1.2dB以内；自适应块长的上下限从[4,32]改到[8,24]，输出SNR波动只有0.3dB；HOS置信门限从0.3改到0.5，输出SNR波动约0.8dB。整体来看，这套算法对参数并不算敏感，只要在合理范围内取，输出结果都比较稳定。这一点在实际工程应用中很重要——你在办公软件里调好的参数，拿到现场数据上不会因为外部环境变化就“崩”掉。
6. 进阶优化方向与扩展思考
6.1 多道联合处理：利用空间维度提升去噪效果
目前这套方法默认逐道独立处理，但地震信号本质上是一个空间连续体，相邻道之间的同相轴是有空间相关性的。进阶方向之一是引入多道联合约束：在处理当前道时，参考相邻道的块分类结果来做多数投票，如果当前道的某个块被判为噪声主导，但左右相邻道同深度位置的块都判为信号主导，那么这个块大概率是信号，只是当前道噪声大盖过了信号。这个策略对强噪声单道记录非常有效，能明显提高同相轴的横向连续性。代价是计算量增加，需要缓存多道的小波系数做联合判断，内存开销更高。
6.2 与自适应时窗滤波的组合：处理非平稳噪声的另一种路子
这套方法在处理平稳噪声（随机噪、白噪）时已经够用了，但如果噪声是非平稳的——比如某一时间段内噪声突然增强了，其他时间段又恢复正常——分层统一阈值就不太合适。一个可行的扩展是结合自适应时窗：先把信号按时间分窗，对每个时窗做独立的噪声方差估计和去噪处理，再把结果拼接起来。这个方法在地震记录中存在强干扰时间段的情况下很有用，但要注意窗边界的接续问题，建议相邻窗之间留10%~20%的重叠，并对重叠区域做加权平均，避免拼接处出现台阶。
6.3 与生成式模型的对比：传统方法还有没有存在价值
最近深度学习去噪很火，我也尝试过用U-Net做地震数据去噪，效果在信噪比指标上确实比传统方法好，尤其是在训练数据覆盖的场景里。但深度学习方法的软肋同样明显：需要大量配对的高质量训练数据，而这个在地震勘探场景里很难获得；训练模型在新工区的地质条件下泛化能力不稳定，经常需要重新标注和微调；推理过程不透明，解释人员不太能接受一个“说不清为什么”的结果。相比之下，传统方法虽然信噪比提升幅度可能略逊于深度方法，但可控、可解释、可复现，这是我为什么在实际项目中还是愿意把传统方法作为首选方案的原因。
6.4 面向实时监测的轻量化改进
如果你要做微震实时监测，那么计算效率就是一个必须考虑的问题。目前这套代码处理一道4000点的数据需要约0.5秒（R2021b的PC），240道需要2分钟左右，实时性还不太够。优化方向有两个：一是对核心循环做MATLAB向量化改写，用矩阵运算替代for循环，大概能提速2~3倍；二是把每层的分块逻辑用C语言写成mex函数，我试过可以再提速5倍以上。这样一道数据可以在几十毫秒内处理完，满足实时监测的需求。
7. 个人经验总结
这套HOS+改进块阈值的方法，前前后后在我手里迭代了快三个月，从最开始只是把小波阈值从软硬两种扩展到块阈值，到后来逐渐把高阶统计量的修正思想融入进去，每一步的改进都是被实际问题逼出来的。刚开始拿传统块阈值去跑野外数据，发现效果是好，但总有少部分弱反射轴被误杀，后来仔细分析才发现问题出在噪声方差估计不准上——野外噪声的非高斯性太强了，MAD方法在这种背景下就失灵了。引入高阶统计量修正之后，这个短板才被补上。
说句实在话，这套方法并不是什么特别高深的新理论，它的核心思想就是两句话：一是不要把高斯噪声的假设生搬硬套到非高斯噪声上，二是做系数处理时要尊重信号的连续结构。但恰恰是这两点，在现成的MATLAB工具箱里是没有实现好的，需要使用者自己根据数据特性去调整和组合。这也是我觉得做信号处理最有价值的部分——真正的好方法不是靠一个标准的工具箱函数就能搞定的，而是要理解底层原理，再把它们组合出适合自己数据的形式来。
最后再分享一个小的操作技巧：无论用哪种去噪方法，在处理完一道数据之后，先看单道时域波形和频谱，确认没有引入明显畸变；再看整炮剖面的同相轴连续性，确认没有破坏横向结构。两道检查都过了，再考虑去批量跑数据。养成这个习惯，能帮你省下大量反复调试的时间。
