高阶统计量+小波块阈值：低信噪比地震信号去噪实战

前几个月处理一批微震监测数据时，遇到了让我头疼的地震信号降噪问题。剖面信噪比低到只有0 dB左右，有效波在随机噪声里若隐若现；常规小波阈值去噪跑完，噪声是淡了，但波峰被磨圆、同相轴断断续续，压根没法用于后续拾取和反演。试了好几版方案后，真正把问题解决的是一条组合思路：高阶统计量（HOS）加上改进的小波块阈值。这篇文章就把这套方法从头讲透，包括为什么二阶统计量不够用、怎么用HOS在小波域里区分有效信号和噪声，以及MATLAB里的完整实现和调参陷阱。适合做信号处理、地球物理资料处理，或者刚接触小波去噪但对效果不满意的同学参考。

1. 地震信号降噪的痛点：低信噪比下常规小波阈值为什么失效

1.1 地震记录里的噪声不是简单的高斯白噪声

很多教材里讲去噪，总假设地震噪声是平稳高斯白噪声，这在实际资料中站不住脚。地震记录中的噪声成分非常杂：环境微震、工业电干扰、风噪、仪器自振、多次波残留等。整体看，随机背景噪声可以近似为平稳高斯过程，但局部会出现强振幅脉冲、相关干扰等非高斯成分。有效的地震反射波、折射波或微震P波/S波，则是典型的非平稳、非高斯信号。

这个“非高斯”特性非常关键。高斯过程只需要均值、方差（二阶统计量）就完全描述；而地震有效波含有相位耦合、突变和子波形态信息，仅看幅度或功率谱根本区分不开。尤其当地震信号强度接近噪声水平时，两者在幅值上可能完全重叠。

1.2 逐点阈值去噪的三个问题

传统小波阈值去噪，是先把信号做小波分解，然后把低于阈值的系数置零或压缩，再重构。实现很简单，但实际用在低信噪比地震资料上，有三个绕不开的坑。

第一个问题是假设独立。逐点阈值把每个小波系数当成独立样本处理，但有效地震事件在小波域里是有结构的，相邻尺度和相邻时刻的系数会一起出现强振幅。点处理会把这种“成团”的结构打散，结果就是同相轴连续性变差。

第二个问题是阈值选择的两难。固定阈值按噪声方差σ²和信号长度N估计，常见的是VisuShrink阈值σ√(2lnN)。在噪声很强时，部分噪声系数幅度会超过阈值、留在结果里；而部分弱信号系数反而低于阈值被削掉。没有任何判别信息去区分“该不该保留”。

第三个问题是只用二阶统计量。阈值去噪本质上是在振幅维度上做分类，可信号和噪声在小波域里的区别不只在振幅，还在分布的“形状”。逐点阈值看不到这一点。

1.3 块阈值：把小波系数从“点”看成“块”

后来研究者提出块阈值，比如NeighBlock、BlockJS，思路是把相邻小波系数分成块，用块内总能量做判决。如果某一块的能量显著超过噪声本底，就保留这一块，否则整块置零。这样利用了系数的邻域相关性，对瞬态信号和地震事件更友好。

但经典块阈值仍然有一个隐含假设：噪声是高斯白噪声。它用块能量与噪声方差比较，对强脉冲干扰和有色噪声很容易误判。于是就有了本文的出发点——在块判决中引入高阶统计量，用“块内非高斯程度”来做二次分类，把小波块阈值从“能量门限”升级为“结构感知门限”。

2. 高阶统计量：给信号噪声分类找一把更准的尺子

2.1 为什么二阶统计量分不开信号和噪声

在地震去噪里，信噪比低的时候，信号和噪声在振幅分布上大量重叠。你用方差做标准，看到的只是“谁的能量大”。但有效地震波在统计上不是高斯的，它含有子波、反射系数序列带来的非对称或重尾特征；而高斯随机噪声的高阶累积量理论上为零。

高阶累积量正是捕捉这种差异的工具。三阶累积量反映分布不对称性，四阶累积量反映分布“尖峰重尾”程度。由于地震信号和很多非高斯干扰都有明显的重尾特征，所以四阶累积量在实际中表现最稳定。三阶累积量对某些对称分布信号不敏感，而地震子波经反射系数调制后，在小波域里的系数分布往往是对称重尾，用三阶反而容易漏检。

2.2 峰度：工程上最顺手的HOS指标

四阶累积量单独用不好选尺度，工程上一般用标准化峰度（Kurtosis）：

K = E[(X - μ)^4] / (E[(X - μ)^2])^2 - 3

对于高斯分布，K = 0；对于拉普拉斯分布或地震子波叠加信号，K明显大于0。随机脉冲噪声的峰度可能更高。

在小波域里，把某一块内的小波系数看成一组样本，估计它的峰度。纯噪声块由于高斯性，峰度在0附近波动；包含有效波的块，因为有相关结构，系数分布会出现异常大的尾部，峰度显著偏高。这样我们就把“块里有没有信号”转化成一个假设检验问题：块峰度超过多少，才能判定不是纯噪声。

2.3 小块峰度估计的方差问题

这里有个工程上必须注意的点：峰度估计需要足够样本才稳定。一个块如果只有4、5个系数，K的方差极大，就算纯噪声也可能算出很高的峰度。根据统计理论，纯高斯样本的峰度估计方差近似为24/L，L是样本数。所以块长不能太小，否则HOS判决本身就是噪声。

我在实际中一般取块长L=16~32。块长16时，24/16=1.5，标准差大约1.22，判为信号块的峰度门限至少要放到2.5以上；块长32时，标准差约0.87，门限可以降到2.0以下。这一条后面参数讨论还会细说。

3. 改进块阈值的设计：块内HOS感知的自适应收缩策略

3.1 经典块阈值方法的短板

经典块阈值的判决函数基本是：
if 块能量 > λ * 块长 * σ²，保留块；else 置零。
λ的选取需要知道噪声方差，而且对所有尺度统一。对于地震信号这种窄带、瞬态信号，有些有效弱信号块的能量可能并没有远大于噪声本底，经典方法容易漏掉；有些强脉冲噪声块的能量足够大，又会被误判为信号块。

还有一个容易被忽略的问题：块能量判决在数学上等价于假设“信号出现会让局部方差变大”，但这只在信号能量比噪声大时才成立。低信噪比下，弱信号叠加在噪声上，能量变化很小，能量判决的灵敏度会急剧下降。

3.2 用HOS给块分类，再分别处理

我的改进思路很简单：分两步。
第一步，还是按块计算小波系数峰度K_b。设定一个峰度门限T_k，这是核心。

• K_b > T_k：判为信号块，块内很有可能含有有效事件或强干扰。
• K_b ≤ T_k：判为噪声块，大概率只有随机噪声。

第二步，按块类别使用不同的阈值策略：

• 对噪声块，直接置零或使用较激进的收缩阈值，把它处理干净。
• 对信号块，用较小的软阈值或保持原系数，避免把有效波细节磨掉。
  在信号块内，也存在残余噪声，所以完全不做处理也不行。我一般对信号块用0.5~0.8倍σ作为软阈值，只把明显的噪声尾巴削掉。

这比单一能量门限聪明的地方在于：即使某个块能量不是特别高，但块内系数分布明显非高斯，也值得保留；反过来，一个能量很高的纯高斯噪声块（比如强白噪声），峰度却接近0，可以放心置零。

3.3 阈值判定与收缩函数怎么配合

具体实现时，我用一个小块滑动窗口，步长可以取块长的一半，这样处理完系数后重叠部分取平均，能减少分块边界处的突变。对每个窗口块：

1. 提取小波系数向量d；
2. 计算块峰度k_b；
3. 计算块内方差并估计噪声水平；
4. 根据k_b是否超过T_k决定用lambda_noise还是lambda_signal；
5. 软阈值处理，结果累加到输出数组。

公式上，lambda_signal = beta_signal * sigma（beta_signal = 0.5~0.8），lambda_noise = beta_noise * sigma（beta_noise = 2~3）。对于判断为噪声块的，也可以直接用零。不过完全置零会把弱边缘搞得更碎，所以推荐用大阈值软阈值，而不是一刀切。

这样一来，我们既保留了经典块阈值的邻域相关性优势，又加入了HOS对信号结构的识别能力，在低信噪比下更稳。

4. MATLAB实现：完整流程和核心代码

4.1 数据准备：合成地震记录和加噪

为了验证和复现，先用一个简单的合成地震道做测试。用Ricker子波和随机反射系数褶积，然后加入高斯白噪声，让SNR=0dB左右。

matlab复制rng(1);
fs = 1000;      % 采样率 Hz
t = 0:1/fs:1;   % 1秒记录
N = length(t);

% Ricker子波
f0 = 30;                       % 主频
ricker = (1 - 2*(pi*f0*t).^2) .* exp(-(pi*f0*t).^2);
% 反射系数（稀疏，模拟地震事件）
refl = zeros(size(t));
refl(150) = 1.0;
refl(400) = -0.6;
refl(750) = 0.8;

% 合成干净地震道
clean = conv(ricker, refl, 'same');
clean = clean / max(abs(clean)) * 0.5;

% 加高斯白噪声，SNR=0dB
sigma_noise = std(clean) / (10^(0/20));
noisy = clean + sigma_noise * randn(size(clean));

这里SNR的公式按均方根幅度定义。合成信号有效波只有三处，其余都是噪声，比较接近微震数据场景。

4.2 小波分解和分块

用db4小波，分解层数为4。小波分解的近似系数不去动，只处理细节系数。

matlab复制wname = 'db4';
level = 4;
[C, S] = wavedec(noisy, level, wname);

% 提取各层细节系数
det = cell(1, level);
for k = 1:level
    det{k} = detcoef(C, S, k);   % 从第一层（高频）到第level层
end

噪声水平估计用第一层细节的MAD：

matlab复制sigma_est = median(abs(det{1})) / 0.6745;

分块函数可以单独写成m函数。以每层细节系数为输入，块长为blockLen，步长为step（取blockLen/2）。为了抑制边界效应，使用重叠平均。

matlab复制function d_clean = hos_block_threshold(d, blockLen, step, sigma, T_k, beta_signal, beta_noise)
    d = d(:)';
    n = length(d);
    % 结果累加器和权重
    acc = zeros(1, n);
    wgt = zeros(1, n);

    lambda_signal = beta_signal * sigma;
    lambda_noise  = beta_noise * sigma;

    startIdx = 1;
    while startIdx <= n
        seg = startIdx:min(startIdx+blockLen-1, n);
        segLen = length(seg);
        if segLen < 4
            break;
        end
        block = d(seg);
        % 块峰度
        mb = mean(block);
        vb = mean((block - mb).^2);
        if vb == 0
            k_b = 0;
        else
            k_b = mean((block - mb).^4) / (vb^2) - 3;
        end

        if k_b > T_k
            out = wthresh(block, 's', lambda_signal);
        else
            out = wthresh(block, 's', lambda_noise);
        end
        acc(seg) = acc(seg) + out;
        wgt(seg) = wgt(seg) + 1;

        startIdx = startIdx + step;
    end
    d_clean = acc ./ max(wgt, eps);
end

末尾不足一个块时我直接break，如果信号尾部有有效事件会漏掉。更稳妥的做法是补零到blockLen，处理完之后再把结果截断回原始长度，代码上多几行，但值得。

4.3 计算块HOS并设定峰度门限

T_k需要事先确定。可以用一个简单的经验公式：T_k = 2 * sqrt(24 / blockLen)。块长16时，T_k≈2.45；块长32时，T_k≈1.73。但真实资料最好用噪声模拟标定，后文会说。

matlab复制blockLen = 16;
step = 8;
T_k = 2 * sqrt(24 / blockLen);   % 约2.45
beta_signal = 0.6;
beta_noise = 2.5;

det_clean = cell(1, level);
for k = 1:level
    det_clean{k} = hos_block_threshold(det{k}, blockLen, step, sigma_est, T_k, beta_signal, beta_noise);
end

4.4 重构与评价指标

处理完细节系数后，把干净的细节系数放回原来的C结构，再用waverec重构。关键在于C的排列顺序：前面是近似系数，然后依次是第level层到第1层的细节系数。

matlab复制C_clean = C;
idxC = S(1) + 1;   % 跳过近似系数部分
for k = level:-1:1
    d = det_clean{k};
    dlen = length(d);
    C_clean(idxC:idxC+dlen-1) = d;
    idxC = idxC + dlen;
end
clean_recon = waverec(C_clean, S, wname);

有了合成干净数据，可以算SNR提升和RMSE：

matlab复制denoised = clean_recon;
snr_before = 20*log10(norm(clean)/norm(noisy-clean));
snr_after = 20*log10(norm(clean)/norm(denoised-clean));
rmse_after = sqrt(mean((clean - denoised).^2));

如果使用真实数据，没有clean参考，可以对比去噪前后的频谱、同相轴连续性以及处理后的残差剖面。

5. 参数讨论：块长度、分解层数、判定阈值怎么调

5.1 块长度：检测灵敏度和分辨率之间的博弈

块长度直接决定峰度估计的稳定性和检测分辨率。块太小，峰度门限只能取得非常大，否则大量噪声块会误判成信号块；块太大，一个块内可能既包含有效事件又包含一堆噪声，峰度会被平均掉，弱信号就漏检了。

我试过的经验：

• 块长8：分辨率最好，但峰度门限要到3以上，效果很差，误检率很高。
• 块长16：折中，微震和反射波效果都不错。
• 块长32：峰度估计更稳定，但两三个采样点的短促事件会被磨掉。

5.2 小波基和分解层数

小波基我推荐db4或sym8。它们时域紧支撑，和地震子波形态接近，能得到更稀疏的细节系数。Haar太块状，不适合光滑子波；biorthogonal在相位上可能引入偏移，不太适合保幅处理。

分解层数一般取4~6。如果采样率1000Hz、主频30Hz，4层就能把主要有效频带和低频背景分开。层数太少，噪声和信号混在一起；层数太多，近似系数会包含一些有意义的低频信号，而你只处理细节系数，会导致低频部分没有去噪。

5.3 峰度判定阈值的蒙特卡洛标定

理论公式T_k=2√(24/L)只是把噪声块峰度近似看作N(0,24/L)下的2σ界，实际上块长度小时分布不一定正态。更可靠的办法是拿纯噪声跑一遍：生成足够多的高斯白噪声序列，做同样的小波分解和分块，统计每个块峰度的95%分位数，用它作为T_k。

matlab复制blockLen = 16;
% 模拟500道纯噪声
K_all = [];
for iter = 1:500
    noise = randn(1, 1024);
    [C_tmp, S_tmp] = wavedec(noise, level, wname);
    for k = 1:level
        d = detcoef(C_tmp, S_tmp, k);
        nseg = floor(length(d)/blockLen);
        for b = 1:nseg
            blk = d((b-1)*blockLen+1 : b*blockLen);
            mb = mean(blk); vb = mean((blk-mb).^2);
            K_all(end+1) = mean((blk-mb).^4)/(vb^2)-3;
        end
    end
end
T_k = quantile(K_all, 0.95);

这样标定的T_k比经验公式更贴你的小波基和分解方式。实测下来，db4、level=4、blockLen=16时，T_k大约在1.8~2.2之间，比理论2.45略低一点。原因是各层噪声小波系数并不是完全独立高斯，存在层间相关性。

5.4 阈值系数的调整

beta_signal和beta_noise也需要调。beta_noise取2~3很稳；beta_signal如果取0.5，信号保真度好，但残留噪声稍多；取0.8，波形更干净，但弱振幅会被压低。如果后续要做振幅反演，beta_signal最好不超过0.6。

6. 实测对比和避坑记录

6.1 合成测试：三种方法效果对比

我做了个对比实验：同一段含噪合成记录，分别用逐点软阈值（全局阈值）、经典NeighBlock能量块阈值、本文HOS块阈值处理。参数如下：

• 逐点软阈值：阈值=σ√(2lnN)
• NeighBlock：块长16，步长8，基于块能量
• 本文：块长16，步长8，T_k=2.2，beta_signal=0.6，beta_noise=2.5

结果用SNR和RMSE比较：

SNR提升只是一个方面。看波形，逐点软阈值在事件起跳处有压缩，NeighBlock在强脉冲噪声附近会误留，HOS块阈值在这几处明显更干净。

6.2 真实资料里的三个坑

第一个坑：强脉冲干扰的峰度比有效波还高。有些环境干扰，比如敲击、风吹缆，在小波域表现为少数极大系数，块峰度非常高，会被HOS判为“信号块”而保留。解决的办法是加一道“能量上限”检查：如果一个块的峰度极高但同时能量也异常高，它更可能是异常噪声而不是地震反射，可以用3~5倍全局平均能量把它单独剔除。

第二个坑：不同层级的细节系数峰度分布不一样。高频第一层受随机噪声影响大，纯噪声块的峰度方差也大；低频细节层噪声相关性更强，纯噪声块的峰度可能略偏移。所以T_k不要所有层统一，最好按层分别标定，或者每层做自己的纯噪声统计。

第三个坑：重叠分块的步长不能设得太小。步长太小会让相邻块高度相关，峰度结果平滑，但计算量变大；步长太大会在边界留下接缝。我用blockLen/2作为步长，图个简单省事，效果也够好。如果追求极致，可以用Hann窗给每个块加权，重叠相加，能进一步抑制边界效应。

6.3 实用技巧：从“看到波形”到“敢拿去用”

最后说几个很多人容易忽略的点。第一，评价去噪效果不要只看SNR，对真实地震资料尤其如此。要拉出单道和剖面对比，看同相轴是否连续、振幅是否保真、高频弱信号有没有被抹掉。第二，处理真实资料前，先拿一段纯噪声（比如初至前的记录）标定T_k，比拍脑袋调参数靠谱得多。第三，如果后续要做走时拾取，可以适当增大beta_signal保留弱信号；如果做AVO或反演，则尽量用较小beta_signal，避免引入人为波形畸变。

这个方法看起来不算复杂，核心就是在块阈值里加了一个“高阶统计量分类器”。但正是这一步，让我那几个晚上反复调阈值的经历变成了历史。如果你也卡在低信噪比地震信号去噪上，不妨照着这个流程试一下：先跑通合成记录，再标定T_k，最后再上真实数据。参数一次调不理想太正常了，但把块长、层数、峰度门限这三个旋钮分开去试，会比从前盲调好得多。