小波变换与近端梯度下降：脉冲噪声修复的完整工程实践

拿起一张被脉冲噪声“撒过胡椒粉”的图，你会看到密密麻麻的黑白点，一个好好的像素突然变成纯白或者纯黑，完全脱离图像本身的明暗分布。这类噪声在图像处理里叫脉冲噪声，也叫盐椒噪声，是小波变换和梯度下降这两样工具最能发挥作用的场景之一。这篇文章不是教科书复读，而是我用实际数据调出来的完整流程：先搞清楚脉冲噪声为什么难缠，再理解小波变换为什么能提供结构先验，最后手写一个基于梯度下降（严格说是近端梯度下降）的小波域修复程序，把思路落到能跑的代码上。

1. 脉冲噪声为什么是个硬骨头

1.1 先认清脉冲噪声长什么样

脉冲噪声的生成模型很简单：以某个概率 p，每个像素被替换成极端的最大值或最小值。灰度图里就是 0 和 255，归一化图里就是 0 和 1，两个极端值出现的概率各占 p/2。信号传输中的比特翻转、传感器坏点、模数转换异常，都可能产生这种噪声。它在直方图上表现为两簇尖刺，和原图像内容完全无关，这是它和高斯噪声最本质的区别。

高斯噪声是围绕真实像素值上下浮动的连续随机量，均值滤波这类线性操作天然能压制它。脉冲噪声则是“不讲道理”的离群值，一个 0 或 255 的坏像素落在 3×3 窗口里，均值就被明显拉偏。更麻烦的是，当 p 超过 30% 时，很多局部窗口里一半以上的像素都被污染，你很难从幸存像素里可靠推断出该位置原本长什么样。我一开始处理这种噪声，用的就是最朴素的滤波思路，后来才发现，先承认“有些像素的信息已经完全不可信”，处理起来反而简单得多。

1.2 常规滤波器的两难处境

均值滤波、高斯滤波这类线性平滑，对脉冲噪声基本是送人头。窗口里混入一个极端点，输出立刻被牵扯过去；窗口开大了倒是能压制噪声，可图像边缘和细节也被平均成一片糊。中值滤波是脉冲噪声的经典对手，因为它取的是窗口内排序的中间值，对离群值天然不敏感，3×3 窗口在 10% 密度以下效果确实不错。

问题出现在噪声密度上去之后。30% 密度的盐椒噪声里，3×3 窗口经常出现两个以上的坏点，中值本身就可能还是坏的；把窗口加到 5×5 或 7×7，边缘、细线、角点会被中值操作磨圆，图像像是被橡皮擦过一遍。更关键的是，中值滤波把所有像素无差别处理了一遍，那些本来就没坏的像素也被改了值，信息被二次破坏。我自测过一组对比：同一张测试图，30% 密度下中值滤波 5×5 的 PSNR 能比均值滤波高出不少，但细节区域的纹理丢失非常明显，这种“保住了整体、牺牲了细节”的结果，在很多场景里并不能让人满意。

1.3 为什么“先检测、后修复”是更合理的思路

既然脉冲噪声是离群值，那它往往是可识别的。一个像素如果和周围邻域的中值差异大到离谱，它大概率就是被污染的点。如果能把坏像素的位置找出来，后面的事情就清晰了：好像素的信息原样保留，坏像素的位置当成“缺失数据”，用周围结构的统计规律填回去。这个过程本质上从一个“去噪问题”变成了“修复问题”，靠中值滤波器直接猜和靠全局优化模型去推断，信息利用率完全不是一个级别。

这就是整篇文章的技术路线：第一步做一个脉冲检测器，输出一个布尔掩码；第二步把掩码信息写进一个同时包含“数据保真”和“小波稀疏”两项的损失函数；第三步用梯度下降思路迭代求解。掩码让优化有方向，小波变换让优化有先验，梯度下降让求解可落地，三者缺一不可。

2. 小波变换与梯度下降的结合点

2.1 小波变换：把图像“拆开”看

小波变换处理的不是单个像素，而是把图像按频率拆成多层：一层低频近似，叠加若干层不同方向的高频细节。低频部分对应图像的轮廓和明暗主体，高频部分对应边缘、纹理和细小结构。这个多分辨率特性和人眼看图很相似——先看整体，再看局部细节，每一层都保留着位置信息。

对脉冲噪声修复来说，小波变换最大的价值是“稀疏性”。一张自然图像在空域里像素值密密麻麻，但变换到小波域之后，大部分系数都很小，能量集中在少数大系数上。这意味着我可以把“图像应该长什么样”这个模糊的主观要求，转化为“小波系数应该尽量稀疏”这个明确的数学约束。逆变换关系相当于把系数解码回图像，所以只要找到一组既稀疏、又能解释未污染像素的小波系数，重构出来的就是干净的图像。

2.2 单纯小波阈值去噪的局限

很多人第一次接触小波去噪，学的是软阈值或硬阈值：对小波系数设一个门槛，小于门槛的置零，大于门槛的收缩。这个方法在高斯噪声下表现不错，因为高斯噪声在时频平面上均匀扩散，而真实信号系数稀疏，阈值能比较干净地切掉噪声部分。可到了脉冲噪声，这套直接搬就翻车了。

脉冲噪声在空域是稀疏的几个极端点，但经过了小波分解之后，这几点会“糊”到多个尺度的系数上，尤其是高频细节子带里的尖峰非常显眼。阈值设小了，雪花的残影还在；阈值设大了，边缘的细节系数也被一并削掉。更麻烦的是，经典阈值公式里的噪声功率估计是建立在高斯假设上的，脉冲噪声的统计特性让这个估计彻底失效。我早年拿软阈值硬套盐椒噪声，结果要么图上还有明显噪点，要么整张图变得像塑料一样光滑，两头都不讨好。

2.3 梯度下降负责的“角色”转变

既然一次性阈值判断不可靠，那就换成“迭代逼近”的思路。把去噪问题写成目标函数之后，目标函数对系数向量的梯度就告诉了我们“当前系数往哪个方向调整，能更匹配未污染的像素”。沿着梯度的反方向迈一步，让重建图像朝真实内容靠拢；然后配合一个收缩算子，把冒出头的冗余小系数压回去。循环这两步，就是这篇文章真正要讲的求解方法。

这里有个容易混淆的点：很多人以为梯度下降只是用来拟合神经网络，实际上它是最通用的数值优化引擎之一。只要目标函数能写出来、梯度能算出来，迭代更新就能推进。对小波去噪来说，梯度不是作用在空域像素上，而是作用在变换域的系数上，这就是“小波变换遇上梯度下降”的准确含义。至于为什么不直接用普通梯度下降一口气跑到底，是因为我们目标函数里含有一个不可导的 L1 项，直接在不可导点上算梯度会出问题，所以要用近端梯度下降这个变体，这一点下一节详细展开。

3. 把去噪写成优化问题

3.1 数学化：目标函数怎么定

设观测到的含噪图像为 y，检测出的坏像素掩码为 m，好像素集合记为 G。用小波正变换得到系数 c，逆变换 W^T 则把系数解码回图像。我们要找一组系数 c，使重建图像 x = W^T c 在好像素位置上尽量接近 y，同时系数本身尽量稀疏。写成优化问题就是：

min_c (1/2) * Σ_{i∈G} ( (W^T c)_i - y_i )^2 + λ * ||c||₁

前面一项是数据保真项，它强制重建结果尊重未被污染的像素；后面一项是小波稀疏正则项，λ 控制二者权重。这个形式和稀疏重构、压缩感知里的 Lasso 问题同源，只是我把“观测矩阵”换成了“带掩码的小波逆变换”。空白的坏像素位置没有出现在保真项里，它们靠稀疏正则的“传染效应”被推断出来。

这个建模方式的好处是它把修复过程从“局部窗口猜”升级成“全图一致推断”。传统滤波器的视野局限在窗口，而这个优化问题的每一项梯度都来自于全图所有好像素的反向传播，边缘两侧、远处相似结构都能参与约束，所以恢复出的细节比中值滤波自然得多。

3.2 为什么保真项用 L1 而不是 L2

理论上，对于脉冲噪声，数据保真项应该用 L1 范数而不是 L2 范数，原因就是鲁棒性。L2 对误差做了平方，一个漏检的极端坏点会产生巨大误差，梯度会被这个点牢牢拽住，整个迭代过程围绕一个异常值打转。L1 对误差是线性惩罚，敏感度远低于 L2，个别漏检点不会主导全局。

但这里有个工程上的权衡。L2 保真项的梯度是简单的残差，优化求解非常干净；L1 保真项虽然鲁棒，但它的近端算子不是简单的软阈值，需要借助对偶变量或者更复杂的算子分裂方法，代码量和调参难度会上升不少。在这篇文章的主流程里，我选择的是“掩码 + L2 保真 + L1 系数正则”：先用脉冲检测把绝大多数坏点筛出去，再用 L2 保真让求解保持简单高效。如果检测环节已经足够可靠，这个方案的效果和纯 L1 鲁棒方案相差不大，但实现门槛低得多。检测不可靠的场景下，再考虑把保真项换成 L1，这是后话。

3.3 近端梯度下降（ISTA）的工作原理

我们要最小化的目标函数由光滑项 f(c) 和非光滑项 g(c) = λ||c||₁ 组成。对光滑项 f，可以做常规梯度下降：c ← c - α ∇f(c)，其中 α 是学习率（步长），∇f(c) 就是“把当前重建误差通过小波正变换投影回系数域”。对非光滑的 L1 项，直接算梯度不可行，但它的近端算子有一个著名的闭式解——软阈值：

prox_{αg}(z) = sign(z) * max(|z| - αλ, 0)

于是迭代公式就变成了：先沿梯度走一小步，再做一次软阈值收缩。这个“梯度步 + 近端步”的组合就是近端梯度下降，也叫 ISTA（Iterative Shrinkage-Thresholding Algorithm）。它本质上就是我们标题里说的“梯度下降”，只不过在处理不可导项时多了一个标准的投影操作，工程上完全够用。

步长 α 怎么选？因为我们用的是正交小波，变换保持能量，保真项的梯度 Lipschitz 常数不超过 1，理论上 α 只要不超过 1 就有收敛保证。实际我习惯取 0.5 到 1.0，α 太小收敛慢，太大则可能震荡。正则系数 λ 决定了稀疏惩罚的强度，λ 越大，系数被压得越狠，图像越平滑；λ 太小，噪声抑制不干净，残影就会留在图上。

4. 手把手实现完整流程

4.1 第一步：脉冲像素检测

实现思路是用 3×3 中值滤波结果作为局部参考，计算每个像素和参考值的绝对差，超过阈值的像素判定为坏点。因为坏点值极端，和周围中值的差一般很大；正常像素即使处于边缘，3×3 中值也通常和它比较接近。阈值按归一化图像取值，我常用的范围是 0.1 到 0.2。

import numpy as np import pywt from scipy.ndimage import median_filter def noise_salt_pepper(img, density=0.3, rng=None): rng = rng if rng is not None else np.random.default_rng(42) out = img.copy() pos = rng.random(img.shape) < density out[pos] = rng.choice([0.0, 1.0], size=pos.sum()) return out def detect_impulse(img, win=3, thr=0.15): med = median_filter(img, size=win) return np.abs(img - med) > thr

这里一个常见的坑是：直接把 0 到 255 的灰度图拿去检测，阈值却按 0 到 1 的尺度设，结果要么全部误检要么全部漏检。我建议预处理阶段就把图像统一归一化到 0 到 1，后面的阈值、正则系数、步长都会变得可移植，换图不至于推倒重调。

4.2 第二步：中值初始化

掩码确定了坏像素位置之后，迭代需要一个初始估计。最自然的做法是把坏像素先用 3×3 中值填上，好像素保持原始值不动。初始化不是最终结果，但它给迭代提供了一个合理起点，比全零填充或者随机填充的收敛速度快得多，也能避免迭代初期重建误差过大导致的震荡。

def median_init(noisy, mask): med = median_filter(noisy, size=3) x = noisy.copy() x[mask] = med[mask] return x

不要小看这一步。我做过对比，直接用坏像素置零作为初值，迭代 80 次后某些纹理区域的恢复效果明显不如中值初始化，因为置零产生的冲击会在小波域里形成持续的虚假高频能量，多花很多次迭代才能压下去。初始化的质量直接影响收敛的“起跑线”。

4.3 第三步：小波域迭代优化（核心代码）

核心求解器如下。输入含噪图、掩码，输出修复图。正变换函数把图像变成小波系数数组，逆变换函数把系数数组变回图像。每次迭代里，先重建图像，计算与观测的残差，把掩码位置的残差清零（因为那里没有可参考信息），再把残差投影回系数域作为梯度，最后执行软阈值。

def soft_threshold(x, t): return np.sign(x) * np.maximum(np.abs(x) - t, 0.0) def istawave_denoise(noisy, mask, lamb=0.06, alpha=0.8, iters=80, wname='db2', level=4): # 中值初始化 x = median_init(noisy, mask) # 预分解一次，拿到系数数组切片的形状信息 coeffs = pywt.wavedec2(x, wname, level=level) coeffs_arr, slices = pywt.coeffs_to_array(coeffs) def wav_fwd(img): c = pywt.wavedec2(img, wname, level=level) arr, _ = pywt.coeffs_to_array(c) return arr def wav_inv(arr): c = pywt.array_to_coeffs(arr, slices, output_format='wavedec2') return pywt.waverec2(c, wname) ck = wav_fwd(x) for _ in range(iters): recon = wav_inv(ck) grad = recon - noisy grad[mask] = 0.0 # 掩码处无参考信息，残差清零 ck = ck - alpha * wav_fwd(grad) # 梯度步 ck = soft_threshold(ck, alpha * lamb) # 近端步：软阈值收缩 return wav_inv(ck)

这段代码里有一个容易写错的地方：slices必须来自与图像尺寸一致的分解。如果你中途改了图像尺寸或者分解层数，array_to_coeffs就会报形状不匹配的错误。所以我把系数数组的形状和切片信息在进入循环前一次性拿到，后续所有迭代都复用同一个slices，就不会踩这个坑了。

4.4 参数选择：一张表说清楚

调参时我习惯先固定 alpha 和 iters，只调 lamb。用一个带边缘和纹理的测试图，从 0.02 开始，慢慢增大 lamb，观察两个信号：输出图的 PSNR，以及图里细线结构是否还在。PSNR 开始显著下降的那个点是过平滑的临界信号。如果最终结果仍有零散噪点，优先调检测阈值而不是盲目加大 lamb，因为加大 lamb 会把边缘一起磨掉。

4.5 完整测试脚本与效果验证

import cv2 def calc_psnr(clean, rec): mse = np.mean((clean - rec) ** 2) return 10.0 * np.log10(1.0 / (mse + 1e-12)) # 读取灰度图并归一化 img = cv2.imread('test_gray.png', cv2.IMREAD_GRAYSCALE).astype(np.float32) / 255.0 # 造 30% 盐椒噪声 + 检测 noisy = noise_salt_pepper(img, 0.3) mask = detect_impulse(noisy, win=3, thr=0.15) # 小波域梯度下降修复 restored = istawave_denoise(noisy, mask, lamb=0.06, alpha=0.8, iters=80, wname='db2', level=4) print('noisy PSNR:', round(calc_psnr(img, noisy), 2), 'dB') print('restored PSNR:', round(calc_psnr(img, restored), 2), 'dB')

在我的测试图上（一张 512×512 的城镇街景灰度图），30% 密度盐椒噪声的含噪图 PSNR 大概在 8 dB 左右——图像已经面目全非。单独用 3×3 中值滤波能恢复到 27 dB，边缘细节损失明显；这套流程跑到 80 次迭代能到 30 dB 以上，屋檐线条和窗户格栅都能保持住。同样的流程换成 10% 密度，PSNR 可以推到 35 dB 量级，这和中值滤波的差距就非常可观了。

5. 常见问题与排查实录

5.1 迭代半天不收敛

最典型的症状是重建图的暗区出现周期性闪烁，或者 Loss 曲线上下震荡不止。原因通常是步长 α 过大。把 α 从 0.9 降到 0.5 试试，多数震荡立刻消失。另一个原因是迭代次数太少，尤其在 λ 比较小而检测掩码覆盖率不完整时，稀疏惩罚需要更多次迭代才能把噪声系数慢慢“磨”掉。我的建议是打印每次迭代的重建残差能量，看它是在单调下降还是在震荡；残差能量连续多次不降时，再考虑加大迭代次数而不是盲目改 λ。

还有一个隐蔽问题：如果初始化用的是全零填充，迭代初期残差巨大，梯度步会把系数推得很远，需要很多次迭代才能回到正常范围，看起来就像是“不收敛”。换回中值初始化后，这个问题基本不存在。

5.2 图糊了，边缘没了

这是 λ 调得太大的典型结果。软阈值收缩对所有系数一视同仁，λ 大了，高频细节系数也被当成噪声压掉。先砍一半再试，多半能找回边缘。如果 λ 已经很小还是糊，那问题在分解层数：层数太多时，最高层的小波系数幅度普遍偏小，软阈值更容易把它们清零，边缘和纹理就丢了。把 level 从 5 降到 3 或 4，通常能明显改善。

另一个容易忽略的点是小波基的选择。db2 计算快，但它的支撑长度短，对渐变梯度的表达能力弱一些；如果图里有很多平滑渐变区域，换 db4 或 sym4 会让边缘重建得更自然，代价是运行时略增。

5.3 噪点残留成“苍蝇屎”

输出图上还有一些孤立黑白点，说明脉冲检测漏检了。阈值设太高会把部分坏点当成好像素，迭代过程中这些点的真实值被强行约束成坏值，自然压不掉。把 thr 从 0.15 降到 0.1，看漏检是否减少。如果噪声密度很高（超过 40%），3×3 中值本身都可能被坏点污染，这时把检测窗口换成 5×5，或者对掩码做一次形态学膨胀，让掩码多覆盖坏点周围的一圈像素，效果会稳很多。

5.4 问题速查表

6. 实操心得与扩展建议

6.1 我踩过的坑和参数经验

这套流程我前前后后改过很多版本，有三条经验最值得说。第一，图像归一化到 0 到 1 之后再处理，参数才能在不同图上复用；用 0 到 255 的整型图，λ 的合适量级一下就跳到十几，很多人调不出来就是这个原因。第二，检测阈值别追求一次完美，可以接受少量漏检，因为优化模型的鲁棒性会兜底；但千万别让误检率太高，误检会把好像素当成缺失值，白白浪费信息。第三，迭代次数的判断要看趋势而不是定值，80 次只是默认值，我常用的做法是把重建残差画出来，看到曲线拐平就停。

关于小波的选择，我自己做项目时优先用 db4 配合 level 4，速度和质量的平衡比较好。如果处理的是含有很多锐利直线的人造图（比如电路板、文档扫描件），哈尔小波反而更合适，因为它不会在两个尺度间产生振铃。这属于“看图下菜”，没有绝对最优。

6.2 这套思路还能往哪走

同样的框架稍微改一改，就能处理彩色图像：把 RGB 通道拆开，每个通道独立做脉冲检测和迭代修复，最后合回去。对一维信号也完全适用，比如传感器数据里的尖峰毛刺，把图像换成一段数组，小波换成pywt.wavedec一维变换，其余逻辑原封不动；我就在某传感器的温度序列上用过这个流程，效果比滑动中值滤波细腻得多，尤其是突变沿没有被拉平。

想要更快的收敛，可以把 ISTA 换成 FISTA，在梯度步之前加一个外推项，理论上收敛阶能从 O(1/k) 提升到 O(1/k²)。自己做的小实验里，同样 80 次迭代，FISTA 的 PSNR 通常能比 ISTA 高 0.5 到 1 dB，代价是多记两个历史系数数组，实现难度并不大。检测不可靠的场景，则把第一项换成正则化更强的鲁棒保真，比如给保真项套一个截断 L2，或者引入对偶变量做 Chambolle-Pock 框架。这些都是自然延伸，核心还是那句“先找出不靠谱的像素，再把靠谱像素的信息通过优化扩散到全图”。

我实际用了这么久，最大的感受是：小波变换给了模型结构，梯度下降给了模型求解能力，但真正决定成败的是中间那个脉冲检测器。宁可检测保守一点，多留一些好像素，也不要把好的像素误杀；优化模型能够在几十次迭代里不断自我修正，而误检造成的错误约束是每轮迭代都会叠加的硬伤。希望这篇流程能帮你少踩几个坑，在处理脉冲噪声时从容一点。