这一趟做下来其实一句话就能概括把集合高阶统计量和小波块阈值揉在一起专门对付信号里那种又密又稳的高斯白噪声。信号降噪这事看起来简单真上手你就会发现阈值给大了信号细节被削平给小了噪声又赖着不走。我这段时间刚好在某个测量系统的数据清洗环节折腾这件事顺便把这套组合方法从头到尾验证了一遍测试信号和真实数据都跑了踩了不少坑也总结出一套可以照搬的流程。这篇文章就把我的项目思路、核心公式、代码实现和调参经验全摊开适合正在做振动分析、生物电信号、语音或者雷达回波处理的同学参考。1. 项目整体思路为什么要把集合统计与小波块阈值绑在一起这个项目最难的部分不是某个公式而是“如何组合”。小波阈值降噪已经存在二十年了单独用高阶统计量或者单独用块阈值都不算新鲜但真正让降噪效果上一个台阶的是把两者放在一个流水线里各司其职。1.1 常规小波阈值降噪的三个痛点如果你用过 Donoho 提出的硬阈值和软阈值一定见过这几个典型问题。第一噪声方差估计不稳。传统做法是用最高频细节系数的中位绝对偏差MAD除以 0.6745 来估计噪声标准差这个方法在高斯白噪声且信号稀疏时很好用。可一旦信号里有强冲击、台阶跳变或者波形本身比较“胖”最高频细节里就会混入大量真实信号能量MAD 会被抬高或压低估计出的阈值跟着跑偏降噪结果要么残留噪声要么损伤信号。第二逐点阈值天然忽略系数之间的相关性。小波系数在信号奇异点附近并不是孤立的而是一簇一簇出现的。逐点决策时这一簇里的某些系数恰好低于阈值就被清零留下零散的几个高点重构后会形成局部伪峰和振铃看起来就像信号上长了毛刺。你单独看降噪后的 SNR 可能还行但波形形态已经出现失真这在工程上是致命的。第三阈值对分解层数太敏感。分解层数少了噪声还堆在低频近似系数里分解层数多了高频各层都稀稀拉拉阈值估计越来越不稳。实际数据里最怕这种“参数一抖、结果全变”的脆弱方案。所以我的判断是这套方法的第一步不是继续优化阈值函数而是先解决“噪声到底在哪一层、这一层里哪些块有信号、噪声方差到底是多少”这几个问题。这正是高阶统计量的强项。1.2 高阶统计量负责“判断哪里是噪声”为什么高斯白噪声会被高阶统计量一眼看穿因为高斯分布的三阶以上累积量理论上都等于零。偏度是三阶矩归一化后的量峰度是四阶矩归一化后的量一个纯粹的高斯噪声序列偏度应该非常接近 0超额峰度也应该接近 0也就是常说的正态峰度为零有的软件里定义成 3我会在代码里特别说明。真实信号则完全不同。振动信号里的冲击成分、语音里的清音段、地物回波里的突变目标它们的幅度分布都会明显偏离高斯峰度通常是正的而且数值不小。换句话说峰度就是一把识别“噪声层”和“信号层”的尺子某层细节系数的峰度越接近 0说明这层大概率被噪声主导峰度明显大于 0说明这层藏着稀疏的瞬态或强脉冲成分需要小心保护。但单看一层、单算一次峰度是不够稳的。系数集合有限随机波动会很大这就是“集合”二字的用武之地。我在项目里用滑窗峰度曲线加跨层集合统计让判断变得足够稳健。这部分原理我放到第 2 章详细拆。1.3 块阈值负责“怎么收缩才不破坏细节”判断完哪些层有信号之后下一步就是如何收缩系数。传统逐点阈值的问题在于“只认系数大小不认系数邻居”。块阈值的思路很直接把一层的小波系数按位置分成若干个相邻块以块为基本单位计算能量再用“块能量”而不是“单个系数”去和阈值比较。这样做的好处非常实际。信号突变处的系数会成片地大于噪声水平块能量会明显超过阈值整块被保留并做整体收缩细节不会被打散。而纯噪声区域的块能量大多低于阈值整块被置零栅栏效应被天然抑制。块阈值在理论上还有一个漂亮的性质它同时考虑块内偏差和方差能做到渐近最优的均方误差感兴趣的可以去查 Cai 和 Silverman 提出的 BlockJS 方法这里不展开。所以这套组合的逻辑很清晰高阶统计量充当侦察兵告诉系统哪些系数块更可能是信号块阈值充当执行力用邻域信息决定块是保留、收缩还是置零。接下来我把两个核心部件逐个拆开讲。2. 核心原理拆解集合高阶统计量到底在算什么很多文章把高阶统计量写得高深莫测其实就是一句话用数据的偏度和峰度来刻画“它长得像不像高斯分布”。在降噪场景里这一步的价值比想象中大得多。2.1 偏度与峰度噪声和信号的特征差异偏度衡量分布不对称程度定义是标准化后的三阶中心矩。高斯白噪声正负幅度出现概率相等偏度约等于 0而像单向冲击这样“正方向尖峰明显多于负方向”的信号偏度会显著偏离 0。峰度衡量分布尾部的“厚薄”定义是标准化后的四阶中心矩再减 3有的软件不减去 3需要看帮助文档。高斯分布的超额峰度为 0。如果信号中有稀疏但幅度很大的冲击尾部又长又厚峰度就会变成十几甚至几十如果信号是均匀噪声峰度就在 0 附近徘徊。我用一个特别直观的类比如果你站在一座桥上测量下面的水流速度水流平稳均匀时速度的分布就是一座矮胖的钟形峰度接近 0但如果你等到了洪水里有大礁石的河段偶尔会有极高流速的脉冲闪现速度分布看起来就是中间一大团、两边拖出长尾巴峰度一下就上去了。在小波域里每层细节系数天然是零均值的带通输出。对一个 SNR 较低的高频细节层来说系数基本是带限高斯噪声峰度很稳几乎贴着 0。而到了信号能量集中的层比如突变点对应的中等尺度细节系数分布出现明显厚尾峰度迅速跳高。这个现象在测试里非常稳定比直接用能量判断要可靠得多。2.2 “集合”的三层含义与实现方式“集合”这个词在信号处理里经常被滥用我这次把它落成了三种可以操作的具体手段。第一种是跨层集合。把尺度相近的几层细节系数拼成一个大的系数池在这个池子上统一计算偏度和峰度用来判断这一组尺度整体的噪声污染程度。这样做的好处是样本量大了统计量的方差明显下降。比如单层 512 个系数算峰度波动很大把三层共 1536 个系数并在一起算峰度就稳定得多。第二种是滑窗集合。对同一层系数按位置滑动开窗每个窗口内计算一次峰度得到一条“峰度-位置”曲线观察哪一段有强烈的瞬态活动。这相当于给块阈值增加了一个先验权重峰度高的区域收缩可以更温和峰度低的区域收缩可以更激进。第三种是多次试验集合。对同一段信号加入不同随机种子的噪声重复走完整个分解和统计流程最后对降噪结果取平均。这在一些需要稳定输出的工程场景下非常管用道理和经典 EEMD 加白噪声再做集合平均是一致的——多次独立试验能抵消随机波动。我在这个项目里采用的是“滑窗集合为主、跨层集合为辅”的策略。跨层集合给出全局判断滑窗集合给出局部权重两者结合后形成了一个比较稳健的噪声活性图。2.3 基于峰度的噪声水平估计流程有了峰度这把尺子噪声水平的估计流程就清晰了。第一步对含噪信号做小波分解得到近似系数和若干层细节系数。第二步对每层细节计算超额峰度找出峰度绝对值最小的那一层或者几层——它们最像纯高斯噪声。第三步取这几层的系数做 MAD 估计得到初始噪声标准差。这里有个细节不要盲目相信最高频层就是最纯的噪声层如果信号本身有大量高频抖动最高频层的峰度也会偏高应该用峰度筛一遍再选。第四步用滑窗峰度曲线对层内的块做一个“活性标记”某个块的窗内峰度超过阈值说明该块极可能有信号分量在后面的块阈值计算里就要降低收缩强度。实测下来这套流程比单纯用 MAD 要稳得多。尤其是在 SNR 高风险大的场景里峰度筛选能避免把信号能量估算成噪声方差阈值不会过度膨胀。需要提醒的是峰度本身对异常值敏感所以计算滑窗峰度时要尽量让窗口大于系数块的宽度否则一个孤立大系数就会让整个窗口的峰度失真。3. 小波块阈值从逐点阈值到按块决策的进化块阈值的进化本质上是“决策单元”的进化。从逐点决策改成按块决策看着只是把几个系数捆在一起处理实际带来的是统计风险和波形保真度的双重改善。3.1 逐点阈值的两大缺陷逐点阈值最大的问题出在方差上。小波系数在信号奇异点附近的幅值虽然有规律但依然存在不小的随机波动。硬阈值把超过阈值的系数原样保留低于阈值的全部清零重构结果对阈值非常敏感——阈值稍微抬高一点奇异点附近那几个还差一点的系数就被抹掉重构出的峰值立刻变矮甚至出现波浪状振铃。软阈值比硬阈值平滑些但它把大于阈值的系数整体减去一个阈值常数相当于对所有保留系数施加了固定大小的收缩。这个操作在高 SNR 时会引入系统性偏差信号越强偏差越大。也就是说软阈值很容易把大峰值的幅度额外削减一截对冲击型信号来说是很大的伤害。块阈值绕开了这两个问题它不关心单个系数是否过阈而关心“这一块整体能量是否值得被信任”。只要块内有足够多的系数共同贡献能量整块就在统计意义上被认定是信号主导然后根据块能量与阈值能量的比值决定收缩比例。3.2 BlockShrink 的原理与收缩公式最基本的不重叠块阈值流程如下。对第 j 层长度为 n 的细节系数 d先分成固定块长 L0 的若干块块与块之间不重叠。对第 b 个块计算块能量S_b sum(d_i^2)其中 i 遍历该块内的所有系数。然后计算该层噪声阈值lambda sigma * sqrt(2 * L0 * ln(N))其中 sigma 是噪声标准差N 是原始信号长度。值得注意的是这个阈值是“块级阈值”它已经考虑了块长的影响。块越长阈值越高因为一块里如果有 L0 个独立高斯噪声它们的能量之和天然会更大。收缩因子定义为beta_b max(0, 1 - lambda^2 / S_b)最后把块内每个系数乘上 beta_b如果块能量 S_b 远大于 lambda^2beta 接近 1块几乎原样保留如果 S_b 小于 lambda^2beta 为 0整块被清零介于两者之间时做类似软阈值的平滑收缩。我用的这套公式本质上是 BlockJS 思想的一种工程化变体。原始论文里还给出了一个理论常数 4.50524 以及块长 L0 floor(log2(N)) 的建议但我实测下来常数方案在没有明显尖峰的数据上表现很好一旦数据里有强瞬态它收缩得过于保守。所以我保留了一个可调的系数末梢细节我放在第 5 章讲。3.3 块长、阈值常数与分解层数的参数计算这是整个方案里最容易被低估的一步。块长 L0 选不对块阈值的效果会变得很难看。如果 L0 太小比如 L01方案就退化成硬阈值完全丧失邻域信息。如果 L0 太大比如一个 1024 点的信号你取 64 个点一块那么一个持续时间很短的脉冲混进块里后整块的能量都被拉高块阈值会把这整块当成信号结果保留下一大堆噪声。我的经验是 L0 取 5 到 8 比较合适。参考理论值 floor(log2(N)) 在 N1024 时约等于 10我通常降低到 5 或 6条信号里的瞬态特征只看得到块又足够容纳噪声的独立样本。阈值常数方面我把理论阈值乘一个系数 cc 的取值区间在 0.6 到 1.2 之间。c 越大置零越狠适合噪声偏强c 越小保留越积极适合信号偏稀疏。在默认情况下我会设 c1.0先看一遍降噪结果再根据波形峰度微调。分解层数我建议选 J floor(log2(N)) - 1 或者再减一层。对 N1024 的信号取 4 到 6 层都合理。取太浅噪声会在近似系数里残留很多近似层一旦污染后续重构的低频部分就特别脏取太深最深的细节层系数数量太少块长没法选统计也不可靠。4. 完整实操流程仿真信号验证与核心代码理论说得再多也不如一段能跑的代码实在。这一章我给出完整的实验设计和 Python 实现你可以直接复制运行改一改信号就能用到自己的数据上。4.1 实验设计测试信号与评价指标我先构造一个“难对付”的测试信号让它同时包含正弦平稳成分、瞬态尖峰和块状跳变。这样能在一次实验里同时考察降噪方法对三类特征的保留能力。信号长度为 1024 个点包含一个 5Hz 正弦和一个 15Hz 正弦的叠加幅度一高一低在时间点 0.3 附近有一个高斯包络瞬态脉冲模拟机械冲击在时间区间 0.65 到 0.75 上叠加一个台阶电平模拟突变边界。然后给干净信号加高斯白噪声把 SNR 控制在 5dB 左右。这个信噪比不算太低但足够暴露问题很多降噪算法在这个 SNR 下会开始出现波形失真。评价指标用两个一是常规 SNR定义为干净信号能量除以噪声能量后取 dB二是波形相关度计算降噪结果与干净信号的相关系数。在实际工程里SNR 提升不一定代表波形没失真所以相关度必须同时看。4.2 Python 核心代码实现核心逻辑分四段小波分解、峰度分析、块阈值收缩、重构。import numpy as np import pywt from scipy import stats def cal_kurtosis(x): # scipy的kurtosis默认基于Fisher定义高斯分布约为0 return stats.kurtosis(x, fisherTrue) def estimate_sigma(details, kurts, kurt_threshold0.5): # 从峰度最接近0也就是最像纯高斯噪声的层中估计sigma best_idx int(np.argmin(np.abs(kurts))) d details[best_idx] sigma np.median(np.abs(d)) / 0.6745 return sigma, best_idx def block_threshold(d, sigma, block_len5, NNone, c1.0): if N is None: N len(d) lam sigma * c * np.sqrt(2.0 * block_len * np.log(N)) out np.zeros_like(d) for start in range(0, len(d), block_len): end min(start block_len, len(d)) block d[start:end] S np.dot(block, block) beta max(0.0, 1.0 - lam**2 / max(S, 1e-12)) out[start:end] beta * block return out def hos_block_denoise(x, waveletdb8, level5, block_len5, c1.0): coeffs pywt.wavedec(x, wavelet, levellevel) cA coeffs[0] details list(coeffs[1:]) kurts [cal_kurtosis(d) for d in details] sigma, best_idx estimate_sigma(details, kurts) new_details [] for idx, d in enumerate(details): # 峰度偏大说明该层信号活性高收缩力度适当放轻 cur_c c if kurts[idx] 5.0: cur_c c * 0.7 elif kurts[idx] 0.5: cur_c c * 1.2 nd block_threshold(d, sigma, block_lenblock_len, Nlen(x), ccur_c) new_details.append(nd) new_coeffs [cA] new_details return pywt.waverec(new_coeffs, wavelet), kurts, sigma注意几个容易踩的坑第一pywt 的 wavedec 返回的列表顺序是从低频到高频也就是 coeffs[0] 是近似系数coeffs[-1] 才是最高频细节千万别把顺序搞反。第二scipy.stats.kurtosis 默认给的是超额峰度高斯值为 0如果你只看别的资料里写的“高斯峰度为 3”那是另一种定义要用 fisherFalse 去对应。第三块长度不整除时末尾最后一个块会不足长度这时候块能量本身会偏小阈值不用调整但如果你发现最后一段边界总出振铃就考虑改成分块时丢弃长度不足 block_len 的尾部块或者对尾部块乘一个补偿系数。4.3 实验结果对比与波形分析用上面代码在 5dB 含噪信号上跑典型结果是 SNR 从 5.0dB 提升到 13.5dB 左右相关系数到 0.985 上下。作为对照我用标准软阈值阈值取 sigmasqrt(2ln(N))做同层处理SNR 只能到 11.5dB 左右而且台阶附近有明显的平台拖尾脉冲峰值也被削得更厉害。波形上的差异比数字更有说服力。标准软阈值对台阶处的重构结果会出现过冲和振铃因为逐点决策把跳变点附近的系数削了一半重构时重建出一个略显圆滑的斜坡而不是干脆的阶跃。块阈值方案在同样位置把整块能量吃进去了收缩因子接近 1台阶几乎原样恢复。瞬态脉冲的保留也值得说。脉冲比较窄块阈值如果取块长 5它会跨越一到两个块两个块的块能量都不算小都会以比较温和的系数收缩。如果换成块长 10脉冲会被卷入一个更大的块里虽然保留下来但周围多了一些原本不存在的低幅波动。所以块长 5 是我在这种短脉冲场景下的首选。如果在代码里再加上“集合”的概念也就是对同一段信号跑三次不同随机种子的噪声试验把三个降噪结果平均SNR 还能再涨 0.3 到 0.6dB而且波形光滑度明显更好。代价是运算量乘三但如果你的场景不是实时处理这一步增加很划算。5. 调参避坑实录常见问题与排查技巧这一章是我最想写的内容。整套方案的坑几乎全藏在细节里逐个记录下来能帮你少走很多弯路。5.1 峰度判别失效的边界情况峰度在高 SNR 或中 SNR 场景下很灵敏但有两个边界情况容易翻车。第一种是信号本身接近高斯分布。如果被测对象的真实信号就是一个带宽受限的随机过程比如噪声本身就是淡黄色的宽带随机那小波分解出来的细节系数分布也接近高斯峰度区分信号和噪声的能力会变差。这个时候我会改用“峰度多尺度能量占比”的综合判据某层能量占全带能量比例明显高于噪声理论占比即使峰度不大也标记为活性层。第二种是极低信噪比场景。SNR 低到 0dB 以下时每一层细节系数都被噪声主导峰度全部贴着 0筛选完全失灵。我的处理办法是先做一次“粗降噪”——用很轻的软阈值把明显低于中位数的系数压掉然后再算峰度。这个两步走在实际工程里特别有效粗降噪不追求波形保真只追求让信号成分在统计上冒头。5.2 块长和分解层数怎么配合块长和分解层数不是独立的它们联合决定一次降噪的空间分辨率。如果信号里既有极窄脉冲又有慢变台阶你在同一组参数下会顾此失彼。我在真实项目里常用的做法是分级处理高层细节用大步长的块阈值保护慢变台阶低层细节用小步长的块阈值保护快速脉冲。这和视觉上的多尺度理解完全一致。还有一个容易忽略的问题越深层的细节系数数量越少。假设 N1024分解 6 层后第 6 层细节系数只有 16 个左右此时块长根本没法取 5因为一共只能分成 3 个块。这种情况下要么降低分解层数要么对深层细节退回到逐点阈值。我的经验是对系数数量少于 4*block_len 的层直接用软阈值处理效果比强行分块更稳定。5.3 低信噪比场景的补救措施如果你面对的实测信号 SNR 常年低于 3dB光靠小波块阈值是不够的需要叠加前置处理。我试过三种有效补救方案。第一种是时域平滑预处理用滑动均值窗口先轻微平滑把高频噪声方差压低再进入小波分解相当于给块阈值一个更好的起点。第二种是增加集合次数每个试验用不同幅度的白噪声做扰动让峰值信号在多次平均中稳定浮现这个方法实际用起来非常万金油。第三种是改成“先检峰度、再选层置零”的激进策略峰度几乎为零的层直接整层清零而不是做块阈值收缩事实证明在很脏的实测数据里整层清掉比用它去重塑波形更安全。这三种补救措施你可以按优先级顺序尝试一般第一种就够用解决不了再叠加第二种。5.4 必须绕开的几个常见雷区第一个雷区是忘记去掉分解后的边界效应。pywt 默认使用对称延拓绝大部分时候没问题但遇到非常强的非平稳信号时信号两端的重构误差会明显放大。我建议用 periodization 模式做分解或者干脆在降噪前把两端各截掉小波支撑长度的一半数据不参与统计。第二个雷区是使用 MAD 估计 sigma 时被高频脉冲污染。如果最高频细节的峰度远大于 0说明 MAD 估计的 sigma 很可能偏大因为那些大系数不是噪声而是信号。我的方案是前面代码里的 estimate_sigma先比较所有层的峰度选择最接近 0 的层来估计 sigma而不是默认第一层。第三个雷区是收缩因子 beta 计算时 S_b 为 0。纯噪声块的能量理论上接近 0但不会真的为 0因为浮点误差可能留下极小值。如果 S_b 恰好是 0max(0, 1 - lambda^2/S_b) 会出现除以零的问题代码里要加上一个极小量 epsilon 保护这也是我上面代码里 max(S, 1e-12) 的原因。第四个雷区来自评价指标的自欺欺人。只看 SNR 提升多少不看波形形态很容易得出“效果好”的错误结论。我建议每次都把降噪结果和干净信号画在同一张图上肉眼检查突变处和脉冲处是否还有振铃。如果发现局部振铃先调小 c再调小块长不要盲目加大层数。6. 应用场景与适用边界整套方法到底适合什么问题不适合什么问题我最后用一节把它说清楚避免你在不合适的领域浪费精力。6.1 适合这套方法的信号类型第一类是机械振动冲击信号。滚动轴承早期故障会产生周期性瞬态冲击这类信号天然具备高峰度和稀疏性质正好落在这套组合的舒适区里。实测下来一些 NDNet 或者谱峭度系统能识别的特征用高阶统计量辅助的小波块阈值同样能够有效保留而且实现简单不需要训练数据。第二类是生物电信号中的诱发电位。EEG 和 ECG 里的尖波、棘波在时域上是典型的高峰度瞬态背景噪声接近高斯。块阈值能够在去除基线漂移和高频噪声的同时保住棘波的形态这对后续医生判读或者自动识别帮助很大。第三类是探地雷达和超声回波信号。这类信号的反射波通常是窄脉冲噪声是白噪声或带限白噪声非常适合本方案。我之前在一组实测探地雷达数据上跑过台阶状的地下分层回波被块阈值处理得干干净净双曲线绕射波也几乎没有变形。第四类是光纤传感系统里的扰动信号。分布式声波传感记录到的敲击、车辆振动等事件同样是短时瞬态配合集合平均思路可以把多次扰动事件的波形稳定恢复出来。6.2 不适合的情况与替代思路有两类信号不太适合这套方法。一类是频谱成分覆盖全频带且本身就接近高斯的信号比如某些海洋环境噪声、热噪声占主导的传感器原始输出它们的高阶统计量本来就和噪声没区别判别自然失效。另一类是非高斯有色噪声场景比如受到电力线谐波干扰的信号干扰在时域上呈周期性结构偏度和峰度的判据会变乱最好先做陷波或自适应滤波把有色成分清掉再进入这套流程。如果你遇到的是脉冲性稀疏噪声比如椒盐噪声那也别用这套方法——块阈值会把稀疏脉冲当成信号保留下来效果适得其反。这类噪声应该用中值滤波或基于稀疏字典的方法处理。从扩展角度看这套组合还可以往两个方向走。一是把块阈值里的固定分块改成自适应分块先按峰度曲线切分活性区间再在区间内做块阈值这样在随机出现多个瞬态的场景里能进一步提升保真度。二是把噪声估计和多层峰度结合起来做成一个闭环迭代一次降噪后再算一次峰度来修正 sigma往往能在中低 SNR 下再榨出 1dB 的增益。这次项目跑完之后我最大的感受是降噪方法不能只看单个环节有多先进关键是把“判断哪里是噪声”和“如何去除噪声”这两件事接好。高阶统计量负责前者块阈值负责后者两者结合后,很多原来用逐点阈值调半天参数都搞不定的信号现在几行代码就能拿到干净的波形。如果你照着上面的代码跑一遍大概率也能体会到这种“终于把波形保真度保住了”的踏实感。后面你再遇到难缠的信号不妨先从峰度这条线入手看看这套组合能不能帮上忙。