脉冲噪声怎么去？小波稀疏表示+梯度下降实战指南

如果你处理过传感器数据、音频信号或者照片一定见过那种突兀的毛刺原本平滑的曲线中间突然冒出一个尖峰图像上零散分布着白点或黑点。这就是脉冲噪声。它不是每个样本都加随机扰动而是少数样本被严重污染幅度可以很大。传统方法里中值滤波是最常用的招数但窗口尺寸稍不留神就会把有用的细节一起抹掉。前一阵我在折腾小波去噪时发现把脉冲噪声建模成一个优化问题用小波变换做稀疏表示再用梯度下降迭代求解能获得比传统滤波更可控的效果。这篇文章就手把手把这个流程走一遍从噪声模型、目标函数到近端梯度下降的原理再到可以直接跑通的Python代码和调参经验。1. 脉冲噪声为什么值得单独处理1.1 脉冲噪声的数学模型脉冲噪声和我们常说的白噪声、高斯噪声完全不是一回事。高斯噪声是“到处都有、幅度不大”的随机扰动而脉冲噪声恰恰相反它出现的位点很少但一旦出现幅度通常非常大甚至可以远远超过正常信号的动态范围。用数学语言描述观测信号 y 和干净信号 x 之间满足y x n其中 n 是稀疏噪声大部分位置为 0只有少量位置是一个较大的冲击值。以图像里的椒盐噪声为例坏像素直接变成 0 或 255其他像素不受影响以工业传感器为例通信丢包或者电源浪涌会在连续正常采样的曲线上留下几个剧烈的尖峰。这个“稀疏但幅度巨大”的特性决定了我们在设计去噪算法时不能简单地套用最小二乘思想。最小二乘认为每个误差点的成本是误差的平方而一个幅度为 100 的脉冲会给目标函数带来 10000 的惩罚。算法为了降低总惩罚会拼命把重建信号往这个异常值方向拉结果就是脉冲噪声不但没有被清除反而“溶解”进了重建结果里。这是脉冲噪声让人头疼的根源也是必须专门设计算法的原因。1.2 均值滤波和中值滤波的短板面对脉冲噪声很多人第一个想到的是均值滤波。均值滤波的原理是在一个窗口内将所有样本求平均用平均值替代中心点。问题是均值天生对异常值敏感。一个脉冲样本混进窗口后会直接改变平均值的走向窗口越大被拉偏得越明显。所以均值滤波在处理脉冲噪声时效果很差它实际上是把异常值的影响“扩散”到了周围样本。中值滤波要好得多。它把窗口内样本排序取中间值因为脉冲噪声排序后通常落在两端中间值大概率是正常的。这就是它为什么能成为脉冲噪声处理“默认方案”的原因。但中值滤波有两个绕不开的问题。第一个是窗口尺寸不好选窗口太小密集的脉冲去不干净窗口太大原本的阶跃、尖峰和图像边缘这些有用结构会被当成“少数派”杀掉。第二个问题是它没有全局结构意识中值滤波只看到局部窗口完全不知道信号在小波域有哪些稀疏特性。实际使用中中值滤波过后的信号虽然表面“干净”了但在阶跃位置会被修圆在边缘位置会产生锯齿后续做特征提取时会非常难受。1.3 更可控的思路把去噪写成一个优化问题我采用的思路干净信号不再通过滤波器试出来而是把它定义成一个优化问题的解。我们需要同时满足两个目标第一重建信号要尽量贴近未被污染的观测值。因为噪声是稀疏的大多数观测点是可信的所以保真项用鲁棒误差让少数脉冲点可以被“忽略”。第二真实信号在小波域通常是稀疏的也就是说大部分小波系数接近零只有少数系数承载真实结构。因此我们希望在重建结果的小波系数上加稀疏正则让算法自动选出最简洁、最合理的结构来解释信号。于是去噪问题变成了一个目标函数最小化问题min_c 鲁棒保真项(y, W^T c) λ × 稀疏正则项(c)这里的 c 是小波系数W^T c 表示把系数合成回信号。后面的所有讨论都围绕怎么把这个优化问题讲清楚、怎么用梯度下降类算法去求解展开。2. 小波变换给信号一个稀疏的“词典”2.1 小波变换在做什么要理解后面算法为什么有效先得明白小波变换到底在干什么。傅里叶变换把信号拆成正弦波的叠加它能告诉你信号包含哪些频率分量但把“某个频率出现在哪一段时间”这个位置信息丢掉了。小波变换则保留了位置信息它用一组不同尺度的小波基函数去和信号做内积每个尺度相当于一种“放大镜”大尺度看到信号的整体轮廓小尺度看到局部突变和细节。用人话说小波变换像给信号做一个多分辨率拆解。离散小波变换会把一个信号拆成一个近似系数序列和若干个细节系数序列。近似系数是低频骨架细节系数记录不同尺度下的局部变化。如果信号是分段平滑的那么绝大多数细节系数是 0只有少数几个系数承载阶跃、尖峰或边缘。这种“大量接近 0、少数比较大”的性质正好是优化问题里最需要的那一类稀疏结构。2.2 小波分解后的系数长什么样在 Python 里用常见的开源小波工具库做一次分解非常直接import numpy as np import pywt t np.linspace(0, 1, 512) x np.sin(2 * np.pi * 5 * t) 1.5 * (t 0.5) coeffs pywt.wavedec(x, waveletdb4, level3) print(len(coeffs)) # 4 print(coeffs[0].shape) # 近似系数 print(coeffs[1].shape) # 第3层细节系数 print(coeffs[2].shape) # 第2层细节系数 print(coeffs[3].shape) # 第1层细节系数得到的 coeffs[0] 是低频近似系数coeffs[1]、coeffs[2]、coeffs[3] 是从低分辨率到高分辨率的细节系数。阶跃位置会同时出现在多层细节系数上但能量集中在少数几个点随机噪声则会比较均匀地散布到几乎每一个细节系数上。这正是正则项能起作用的原因噪声想要“解释”自己就必须动用大量小系数而稀疏正则恰恰会对这种动作施加惩罚。2.3 为什么选择在系数空间做优化我们不直接优化原始信号而是优化小波系数 c。原因有两个。第一正交小波变换是线性可逆映射系数空间和信号空间是同一个空间的两种坐标表达。在系数空间做迭代更新配合软阈值操作非常自然软阈值就是在小波系数上做收缩。第二真实信号在小波域的“有效维数”比样本数低得多。一个 512 点的信号真正起作用的细节系数可能只有几十个。梯度下降在一个有效维数很低的空间里搜索就不容易被那些少数脉冲点带偏。很多人在去噪时忽略了这一点直接在样本点上一格一格地调整信号很容易陷入局部噪音中出不来。小波变换相当于给了我们一个筛子把有信息的部分和没信息的部分提前分开了。3. 目标函数鲁棒保真项 小波稀疏正则3.1 为什么 L2 保真会栽在脉冲噪声上现在来看目标函数的设计。最直觉的保真项是 L2 范数也就是让重建信号在平方误差意义下贴近观测信号min_x ||y - x||_2^2 正则项这个目标函数对高斯噪声非常合适因为高斯噪声的误差分布本身就是高斯分布平方误差最大似然解就是均值。但脉冲噪声完全不是这样一个幅度为 100 的脉冲平方误差是 10000而一个正常样本的误差 1 平方后只有 1。在总惩罚里异常值占据绝对主导地位优化过程会拼尽全力去追赶这个异常点结果恢复出来的信号仍然带着脉冲的影子。改用 L1 范数之后异常值的惩罚从“平方增长”变成“线性增长”幅度 100 的脉冲贡献是 100不会再压过大量正常点。L1 保真项在统计上的意义是“寻找能让多数样本残差尽可能小的位置估计”它天然抵抗少量极端值。3.2 用 Huber 损失平滑 L1纯 L1 范数在 0 点不可导梯度下降直接求导会卡在奇异点上。我的做法是使用 Huber 损失替代纯 L1。Huber 损失的定义是当 |r| δ 时ρ(r) r^2 / (2δ) 当 |r| δ 时ρ(r) |r| - δ/2在 0 附近是一个窄的二次抛物线过了阈值 δ 之后就变成线性增长。这样既保留了 L1 对异常值的鲁棒性又保证了梯度在整个定义域上都存在。δ 的取值决定了这个近似的“锐利程度”δ 越小越接近纯 L1梯度变化越剧烈δ 越大曲线越光滑但对异常值的鲁棒性会下降。我在实践中的默认值是 δ 0.05 × 信号标准差信号到脉冲的动态范围差异大的时候可以再调大一点。3.3 正则参数 λ 的作用和直观理解目标函数里的 λ 控制稀疏正则的重量λ 越大软阈值会把越多的小波系数压成 0重建信号越平滑但也越容易磨掉真实细节 λ 越小保留的系数越多重建信号越接近含噪观测脉冲毛刺可能残留。实际调参时我会看两个信号特征。如果重建结果里还看得到明显的孤立尖峰说明 λ 偏小应该调大如果阶跃位置被明显抹平高频振动的振幅变小了说明 λ 偏大应该调小。网格搜索把 λ 按对数间隔取值比如 [0.02, 0.05, 0.1, 0.2, 0.5]用小波去噪后的信噪比或者视觉感受选一个最优值。4. 算法实现梯度下降遇上小波系数软阈值4.1 在系数空间计算梯度目标函数由两部分组成保真项和正则项。保真项对小波系数 c 是有梯度的因为 Huber 损失可导。我们把重建信号 x_hat W^T c 代入保真项得到残差 r y - x_hat。Huber 损失对残差的导数是ρ(r) r / δ当 |r| δ ρ(r) sign(r)当 |r| δ根据链式法则保真项对系数 c 的梯度是gradient -A^T ρ(r)其中 A 是小波重构算子A^T 是小波分析算子。因此在负梯度方向上更新系数就是c ← c lr × wavedec(ρ(r), 小波)这里直接用开源库把 ρ(r) 做一次离散小波分解就能得到所有子带的梯度分量。这一步是整个算法里最核心也最容易踩坑的地方不要在样本域更新 x 之后再变回系数而是直接用分析算子把梯度带回系数空间各子带自动对齐代码更简洁计算也更稳。4.2 近端步骤软阈值保真项处理完之后还要处理小波稀疏正则项。λ||c||_1 在 0 点不可导不能直接代入梯度更新。标准做法是“近端梯度下降”把目标函数拆成“可导部分”和“不可导但那部分有简单近端算子”两部分。对于 L1 范数它的近端算子就是软阈值函数soft_threshold(c, t) sign(c) × max(|c| - t, 0)每一次迭代先做保真项的梯度更新再做一次软阈值收缩收缩量为 lr × λ。需要注意近似系数 c[0] 承载的是信号低频轮廓不应该被强制稀疏化所以软阈值只作用于各层细节系数。如果对近似系数也做阈值重建信号的背景水平会被整体压低出现不该有的直流偏置。4.3 一个可以直接运行的 Python 实现把上面的逻辑串起来就是一维信号的完整去噪函数import numpy as np import pywt def huber_grad(r, delta): return np.where(np.abs(r) delta, r / delta, np.sign(r)) def soft_threshold(x, thr): return np.sign(x) * np.maximum(np.abs(x) - thr, 0.0) def wavelet_denoise(y, waveletdb4, level3, lam0.1, deltaNone, maxiter200, lr0.4): if delta is None: delta 0.05 * np.std(y) c list(pywt.wavedec(y, wavelet, levellevel)) n len(y) for _ in range(maxiter): x_hat pywt.waverec(c, wavelet)[:n] r y - x_hat g huber_grad(r, delta) grad_c pywt.wavedec(g, wavelet, levellevel) for i, subband in enumerate(c): subband lr * grad_c[i] if i 0: subband[:] soft_threshold(subband, lr * lam) return pywt.waverec(c, wavelet)[:n] # 构造测试信号 np.random.seed(42) t np.linspace(0, 1, 512) x np.sin(2 * np.pi * 4 * t) 1.2 * (t 0.5) 0.3 * np.sin(2 * np.pi * 30 * t) y x.copy() pulse_positions np.random.choice(512, 40, replaceFalse) y[pulse_positions] np.random.randn(40) * 3 * np.std(x) x_hat wavelet_denoise(y, lam0.08, maxiter200) def snr(clean, est): return 10 * np.log10(np.sum(clean ** 2) / np.sum((clean - est) ** 2)) print(含噪信号 SNR:, snr(x, y)) print(去噪后 SNR:, snr(x, x_hat))我在实际运行中的经验是含噪信号 SNR 大概在 8 dB 左右时这个函数能把 SNR 提升到 15 dB 以上。如果脉冲密度更高适当增大 λ如果细节丢失明显就减小 λ。代码里 lr 取 0.4对正交小波来说是一个比较稳妥的步长不会太震荡收敛速度也还过得去。5. 实验与效果在一维信号上看结果5.1 怎么生成带脉冲的测试信号做实验之前先明确评价标准。我用的是信噪比 SNR定义是信号能量和噪声能量之比的对数形式SNR 10 × log10( sum(x^2) / sum((x - x_hat)^2) )测试信号包含三部分一个 4 Hz 低频正弦、一个 30 Hz 高频正弦还有一个在 0.5 时刻出现的阶跃。这个组合能同时考验算法对平滑结构、高频细节和边缘阶跃的保持能力。给干净信号随机挑选 40 个位置叠加一个大幅值扰动就得到了含脉冲噪声的观测信号。实验里固定随机种子保证结果可复现。这样你就可以自己跑一遍看看每一步输出信号的波形变化比只看最后的数字更有感觉。5.2 去噪效果与量化指标在同样的测试信号上中值滤波窗口取 3 和 5 都能去掉不少脉冲但结果显示30 Hz 的高频分量在窗口内被明显削平。小波加梯度下降的方法在 SNR 指标上的提升通常更明显而且阶跃的上升沿保持得比较陡。我的经验是不要只看总 SNR。脉冲噪声的密度较低时可能一个极端尖峰被压下去之后总能量就上来了但细节可能已经破损。所以要同时观察重建信号和干净信号的误差曲线特别关注阶跃位置附近有没有过冲高频分量有没有被衰减。如果误差曲线在非脉冲位置也有规律摆动说明 λ 太大把真实结构压掉了。5.3 与中值滤波对比时的选择中值滤波的优势是速度快、实现简单对于密度极低且幅度极大的脉冲效果不差。但如果信号本身带高频成分中值滤波就是一个低通性质的局部排序操作它会模糊边缘削平尖峰。小波加梯度下降的优势在于整个优化过程有一个明确的保真项和稀疏先验你可以通过调整 λ 精确控制“干净程度”和“细节保留”之间的平衡。什么时候该用中值滤波当你的任务只要求“把毛刺去掉曲线平滑即可”对细节结构不敏感时中值滤波更省事。什么时候用本文这套方法当信号中有阶跃、尖峰、边缘、高频振动这些需要保留的结构或者你需要一个可以反复调试的全局优化模型时这个方法明显更合适。6. 常见问题排查与调参技巧6.1 恢复结果里还是有脉冲毛刺出现这种情况第一优先级检查 λ。λ 太小软阈值没有把噪声对应的那些小系数压到零残留的幅度就会在重建信号里表现为小毛刺。把 λ 从 0.05 逐步调到 0.2观察毛刺是否消失。第二优先级检查 δ。Huber 损失的 δ 如果特别小梯度在 0 附近接近于不可导迭代过程可能会在几个系数之间来回震荡导致结果不稳定。这时候可以把 δ 提高到 0.1 倍信号标准差。还有个容易被忽略的细节是步长 lr。lr 过大时近端梯度下降可能在最优解附近弹跳看起来像去噪没收敛实际上只是没走稳。把 lr 降到 0.2 到 0.3振荡通常会消失。6.2 恢复结果太光滑细节被磨平这是 λ 偏大的典型表现。细节系数被软阈值压得太狠真实的高频结构和阶跃信息也被当成噪声处理了。先减小 λ同时检查分解层数。层数太少高频细节只有一层系数承载一旦软阈值力度偏大细节就没有退路适当增加到 4 层让高频信息分散到多尺度上每一层只需要做一点轻微收缩整体细节保存得更好。假如换用小波基db4 偏向平滑Haar 或者说 db1 对阶跃更友好。如果你的信号本身是分段常数信号用 Haar 小波做正则往往比 db4 保留更锋利的阶跃。6.3 不收敛或者收敛很慢刚跑完第一次迭代就想看效果肯定不对近端梯度下降通常需要迭代 100 次以上才会比较稳定。如果迭代到两三百次后目标函数还在缓慢下降可以用 FISTA 加速。动量思想的加入很简单在每次迭代时先把当前估计和一个动量项差值叠加再走梯度更新和软阈值步长可以适当加大到 1.0。另一个收敛问题是如果你改变随机种子或真实信号能量很大delta 的计算会变化。默认 delta 设置为 0.05 × std(y)对大多数场景适用。但如果信号本身动态范围极宽std 可能被少数脉冲拉大导致 delta 过大鲁棒性下降。更稳妥的做法是用信号的绝对中位差计算尺度在代码里就是 np.median(np.abs(y - np.median(y))) × 1.4826然后乘以 0.1 到 0.2。6.4 小波基、分解层数和 δ、λ 的选择速查我把调参经验整理成一个速查表参数作用经验值/选择方向调整方向λ稀疏正则力度0.02 到 0.5有毛刺增大细节丢失减小δHuber 平滑参数0.05 到 0.2 倍尺度震荡增大太钝减小lr梯度步长0.2 到 0.5不收敛减小慢可试 FISTAlevel小波分解层数3 到 5信号长/细节多可增加小波基具体小波形态db4、sym6、Haar平滑选 db4阶跃选 Haar7. 扩展图像椒盐噪声也能这样处理7.1 二维小波版本思路一维信号的流程可以直接搬到二维图像。图像里的脉冲噪声就是椒盐噪声像素点变成纯黑或纯白。处理时把一维小波分解换成二维小波分解pywt.wavedec2 和 pywt.waverec2。代码逻辑完全一致先对含噪图像做小波分解得到二维系数列表然后在迭代中反复进行“二维小波重构、计算 Huber 梯度、二维小波分解回系数、对细节子带做软阈值”四步。需要注意二维小波分解的细节子带是三个方向水平、垂直、对角线。软阈值应同时施加在三个方向的细节子带上。图像去噪时 λ 和一维不一样我的经验是取一维的 1/2 到 1/3 左右因为图像小波系数幅值分布更宽过多收缩会让边缘发虚。如果你只有低密度椒盐噪声甚至可以直接先用脉冲检测把坏点标出来然后只对坏点周围的系数做迭代修复能保留更多纹理。7.2 彩色图像和通道处理彩色图像有 RGB 三个通道。如果直接对每个通道分别跑小波去噪通道之间的相关性会被破坏恢复出的图像容易出现彩色噪点或伪影。更推荐的做法是先把图像转到 YCbCr 颜色空间Y 通道承载亮度细节对 Y 通道使用小波加梯度下降Cb 和 Cr 通道因为脉冲噪声通常也发生在色差通道可以只做轻度修复或者不处理。如果必须在 RGB 上处理统一使用同一个 λ并且对三个通道分别做归一化避免某个通道的噪声幅度不同导致恢复色偏。实际项目里还可以先把椒盐噪声点检测出来训练一个简单的掩膜然后用小波优化只填充掩膜位置这对彩色图像尤其有效。8. 写在最后我的实操体会我最初上手时直接用 L2 保真项跑同一个流程结果脉冲毛刺根本没有被抑制恢复信号里留下几个明显的鼓包。换成 Huber 损失之后问题立刻缓解。这件事让我印象很深在脉冲噪声场景里保真项的选择比小波基的选择更关键缺失鲁棒保真项的模型后面正则再精细也会被异常值带偏。如果你手里有真实的传感器信号建议先做一个简单的脉冲检测把明显异常的采样点标记出来再去跑小波优化。这样做有两个好处一是可以让 Huber 损失的 δ 不必调得很小二是可以让算法把精力放在正常点和异常点之间的结构填充上。后面可以把这种思路扩展成自适应 λ 或者加入脉冲密度估计从“手调参数”进一步变成“自动识别”。但基础思路不变小波负责给信号找一个稀疏坐标梯度下降负责在这个坐标里稳定地找最优解二者配合才能在脉冲噪声面前既足够鲁棒又保留住真正有用的细节。