小波变换与梯度下降结合：脉冲噪声迭代去噪原理与Python实现

脉冲噪声这玩意儿，在真实图像处理里远比教科书上描述的烦人。教科书爱拿高斯噪声举例子，可工程上遇到的往往是传感器坏点、传输信道误码、存储介质损伤带来的椒盐噪点——画面上随机位置突然冒一两个纯白或纯黑的像素，密度高了以后整张图跟筛子似的。

面对这种噪声，很多人第一反应是中值滤波，因为效果直观、实现也快。但如果你把噪声密度拉到0.3甚至0.5，再碰上图像本身有密集的纹理，中值滤波就开始露怯：窗口小了滤不干净，窗口大了细节糊成一片。我后来换了一条路，把去噪明确建模成一个最优化问题：小波变换负责把图像从像素域送进一个更能分离信号与噪声的坐标系，梯度下降负责在这个坐标系里反复迭代，逼近真实的干净图像。

这篇文章不适合只想“调个库就算完事”的人。它更适合你希望搞清楚“脉冲噪声为什么难处理”“正则化为什么有效”以及“小波变换和梯度下降究竟如何协作”的读者。我会从噪声建模、小波原理到完整的Python实现，把整套方案从头到尾走一遍，顺便把我在实际调参过程中踩过的坑都交代清楚。

1. 脉冲噪声的脾气：稀疏的极端值，不是温和的随机扰动

1.1 脉冲噪声的数学模型

脉冲噪声也叫椒盐噪声，数学模型很直白：对于一张图像X，每个像素以概率p被污染，污染后的值要么是极小值（椒，0），要么是极大值（盐，255），原来是什么内容完全不保留。写成公式就是：

Y(i,j) = X(i,j)，概率为 1-p Y(i,j) = 0 或 255，概率为 p，其中0和255各占p/2

这里的p就叫噪声密度。密度0.1意味着每10个像素里平均有1个坏点。可能你觉得10%不算多，但512乘512的图像里就有26000多个坏点，已经足够让整张图像看起来到处是黑洞和白斑。

对比高斯噪声：高斯噪声每个像素都被污染，但每个像素的偏离都不大；脉冲噪声只污染少量像素，但污染值偏离极大。这两种噪声在统计特性上完全不同，直接决定了后面我们为什么要选L1正则而不是L2正则。

1.2 为什么中值滤波在重度噪声下会失效

中值滤波的思路是：图像里某个像素被噪声污染了，那我就用邻域内的中位数替代它。原理简单，对稀疏脉冲也确实有效。不过缺点也是几个层面的。

首先，噪声密度上去以后，一个窗口里可能出现多个坏点，中位数本身就被污染了。比如3x3窗口里如果有3个坏点，中值很可能就是一个坏点；如果窗口开大，虽然抗噪声能力强一点点，但代价是细节丢失更严重。

其次，中值滤波是局部窗口操作，它天然抹平纹理和细线条。图像上的文字边缘、斑马线、织物纹理，只要窗口范围里有结构变化，就会被“团成一团”。我用camera测试图做过实验，3x3中值滤波在10%密度下能压掉大部分噪声，但人物脸部、衣服边缘的细节丢失很明显。

更隐蔽的问题是，中值滤波没有任何全局信息，它不“知道”图像本身有什么规律。我们做去噪时，要么用局部窗口硬扛，要么用全局先验引导恢复。后者就是优化方法的出发点。

1.3 把去噪问题翻译成优化语言

把去噪当成优化问题，就是找一个干净的X，让它同时满足两个要求：一是X要贴近观测到的Y，这是数据保真；二是X要满足我们对“自然图像”的认知，这是正则化。

写成数学形式就是：

min_X 0.5 * ||Y - X||^2 + lambda * R(X)

R(X)是惩罚项，它约束X不能乱来。对高斯噪声，L2保真项是合适的选择；对脉冲噪声，理论上用L1保真项更匹配，但工程上我们仍然可以用L2保真项，只要在正则项上引入足够强的稀疏先验，这套框架照样非常好用。

R(X)怎么选？这就轮到小波变换登场了。小波的职责是给图像一个坐标系，让“自然图像”在这个坐标系里尽量稀疏，让噪声和图像分家。下一节详细说。

2. 小波变换：寻找信号在变换域的稀疏表示

2.1 小波分解到底把图像变成了什么

小波变换和傅里叶变换的差别在于：傅里叶变换把信号拆成不同频率的正弦波，每个分量在整条时间或空间轴上展开，它擅长表达平稳信号，却不擅长表达局部的突变；而小波变换同时保留频率和位置信息，就像把图像拆成一系列“带着位置的频率片断”。

对图像做二维离散小波变换（DWT），一次分解得到四个子带：近似分量cA（低频概貌）、水平细节cH、垂直细节cV、对角细节cD。迭代对cA继续分解，就能形成多层金字塔结构。cA是缩略图，细节子带记录各个方向的结构信息。

对一张自然图像，cA的能量最大，细节子带里大多数系数接近0，只有边缘附近才有较大的幅值。这就是“稀疏性”：图像在小波域被压缩了。一张512x512的camera图，三层小波分解后，超过90%的细节系数绝对值小于某个很小阈值，你把这部分系数清零再重构，视觉上几乎看不出区别——图像压缩算法和去噪方法都吃这个红利。

2.2 脉冲噪声在小波域里的表现形式

脉冲噪声跟高斯噪声在小波域里的表现完全不一样。高斯噪声能量均匀铺开，小波系数整体偏大，但不太存在个别极端值；脉冲噪声则是极少数像素上的尖峰，相当于图像里突发的局部突变，做DWT之后，能量会扩散到各层细节子带的对应位置，表现为一批幅度异常大的孤立系数。

有意思的是，脉冲噪声在小波域依然保持着“稀疏性”——只不过它是另一种稀疏：噪声系数数量少但幅度大，真实图像系数数量多但有很多小幅度尾巴。一次软阈值能砍掉幅度小的系数，但要区分“属于边缘的大系数”和“属于噪声的大系数”，单靠一次阈值远远不够。这正好给后续迭代优化留出了空间。

2.3 用Python快速完成正向/逆向小波变换

Python里做DWT最顺手的是PyWavelets库，命令非常简洁。

import pywt import numpy as np from skimage import data, img_as_float from skimage.metrics import peak_signal_noise_ratio as psnr img = data.camera() img = img_as_float(img) # 三层小波分解 coeffs = pywt.wavedec2(img, wavelet='db2', level=3, mode='periodization') # 返回的第一个是低频近似，后面是各层高频细节 cA, (cH3, cV3, cD3), (cH2, cV2, cD2), (cH1, cV1, cD1) = coeffs # 重构 recon = pywt.waverec2(coeffs, wavelet='db2', mode='periodization')

几个要点：wavelet选'db2'、'sym4'这类正交小波基就够用；mode='periodization'让各层尺寸保持对半关系，小波算子近似正交，这对后面梯度下降的步长选择非常重要。如果数据没有归一化到[0,1]，软阈值里的阈值量纲会跟你闹别扭，我习惯先归一化再做所有处理，最后再乘回去。

3. 目标函数：数据保真项与稀疏先验的博弈

3.1 L2保真项：噪声的统计特性决定距离度量

目标函数里最直观的部分是保真项0.5 * ||Y - X||^2。它是Y和X之差的平方和，也叫L2范数。为什么选它？因为高斯噪声的最大似然估计对应L2，最小二乘在这种噪声下是统计最优的。

但脉冲噪声不是高斯噪声，理论上更匹配的保真项是L1或者加权L1，因为脉冲噪声是稀疏的大偏差，L1对离群值更稳健。不过完全抛弃L2会导致问题非线性更强、求解更复杂。这篇走工程路线：在迭代优化框架里，保真项仍然用L2，但在正则项上引入强稀疏先验，已经能在脉冲噪声场景下拿到很不错的结果。如果你想更极致，可以在L2前面加对角权重W，让被污染位置的权重趋近0，这就是鲁棒回归的思路。

记住这个结论即可：L2配对高斯噪声，L1配对脉冲噪声，实际工程往往折中处理。

3.2 L1正则项：为什么稀疏先验比平滑先验更适合脉冲噪声

正则化R(X)有很多选择。最常用的是总变差（TV）正则，它鼓励图像分片光滑，对高斯噪声很有效，但会把纹理也平滑掉。我这里用的小波域L1正则：

R(X) = lambda * ||WX||_1

WX表示图像的小波系数，||·||_1是系数的绝对值之和。这个正则的含义是：希望图像的小波系数尽量稀疏，也就是大部分系数为0或很小。为什么对脉冲噪声有效？因为干净图像的小波系数天然稀疏，而脉冲噪声会把很多噪声系数顶上很大的值；让系数稀疏化，就是同时“抑制噪声系数”和“保留真实边缘系数”。

这里有一个经常踩的坑：L2正则（岭回归）也惩罚大系数，但它是平方惩罚，边缘处的真实大系数很容易被一道压得糊掉。L1是线性惩罚，具备“对小系数推得更狠、对大系数惩罚相对温和”的特性，所以L1在保留边缘上天然优于L2。这也是稀疏表示领域的标准正则通常选L1而不是L2的根本原因。

3.3 正则系数lambda的选取策略

lambda是保真项和正则项之间的天平。lambda太小，噪声残留，图像还是花的；lambda太大，图像变得过平滑，细节全没了。没有一劳永逸的公式，但有几个实用策略。

第一，数据归一化后，lambda在0.05到0.5之间做一次网格搜索，用PSNR或SSIM挑最优。第二，观察迭代收敛后的残差分布：如果残差的方差接近你估计的噪声方差，这个lambda基本合理。第三，经验上高噪声密度用更大的lambda，低噪声密度用小lambda。

提醒一句：lambda的具体数值依赖小波基的范数和分解层数，换小波基、换图像尺度，效果都会漂移。所以最好是代码里写成可配置参数，跑完看指标再调，不要指望一次性调到最优。

4. 梯度下降进场：从一次软阈值到迭代软阈值

4.1 对不光滑目标函数，普通梯度下降为什么棘手

把目标函数写完整，在小波系数域上更简洁：

min_w 0.5 * ||Y - W^T w||^2 + lambda * ||w||_1

这里w是小波系数向量，W^T是逆小波变换。第一个部分是光滑的，第二个部分含绝对值，在w=0处不可导。普通梯度下降要求目标函数可导，遇到不可导点就没法算梯度，直接强上会振荡、发散。

解决办法是使用近端梯度下降，也叫ISTA（Iterative Shrinkage-Thresholding Algorithm）。它将目标拆成光滑部分+非光滑部分：光滑部分正常走梯度，非光滑部分用“近端算子”处理。对L1范数来说，近端算子恰好就是软阈值函数。直觉上就是把系数朝0方向压缩，压缩量固定为lambda，压缩到0就停手。这个软阈值操作本身，就等价于对小波系数做了一次最优的稀疏化。

4.2 近端梯度下降（ISTA）的推导与更新公式

ISTA的迭代式特别简单：

1. 用当前系数w重构出图像x = W^T w；
2. 计算残差r = x - Y；
3. 对残差做小波变换，得到梯度grad = W r；
4. 系数沿梯度负方向走一步：z = w - alpha * grad；
5. 对细节系数做软阈值稀疏化：w_new = soft_threshold(z, alpha * lambda)。

这样迭代100到300次，系数w会逐步收敛到一个平衡点：既让重构图像贴近观测Y，又让小波系数保持稀疏。

还有一个很容易被忽视的点：ISTA的第一次迭代，因为初始w取的是噪声图像的小波系数，残差为0，梯度为0，所以第一步做的就是经典的一次性小波软阈值。之后的每一次迭代，都在对它做“重构—算残差—修正—再收缩”。所以把ISTA看作不同层次的一次软阈值的重复精修，是一点毛病没有的。

4.3 步长选择：Lipschitz常数与正交小波的特殊性

步长alpha是ISTA能否收敛的关键。条件是不能超过光滑部分梯度的Lipschitz常数L，否则迭代发散。本文目标函数光滑部分的L等于小波变换算子W的谱范数平方，对正交小波来说L=1，所以alpha=1理论上没问题。

但实际用pywt时，非正交小波或边界延拓会导致算子不严格正交，alpha=1有时会在前几步出现轻微振荡。我个人的工程习惯是：先用alpha=1跑，如果PSNR曲线振荡，就把alpha减到0.8或0.5；再不行就检查是否用了非正交小波，或者mode没设成periodization。对90%的场景，alpha=1配合db2/sym4与periodization模式，收敛都比较稳。

至于FISTA加速，就是在系数更新时加上前一步的外推动量，收敛速度通常能快一个数量级，工程落地时再考虑，初学者先跑通普通ISTA最重要。

5. 完整实操：基于Python的脉冲噪声去除全流程

5.1 环境准备与数据获取

依赖库是numpy、scipy、scikit-image、PyWavelets。安装命令：

pip install numpy scipy scikit-image PyWavelets

测试图用skimage自带的camera图，方便复现。如果手头有其他灰度图也可以，建议先用灰度图试通；彩色图可以逐通道处理，但要注意通道间相关性可能影响指标，最好转成YCbCr后只在Y通道上做去噪，Cb和Cr通道可以用轻量处理。

5.2 脉冲噪声注入与评价指标

先写注入函数。skimage的util.random_noise也能做，但自己实现更清楚：

def salt_pepper(img, prob, salt=1.0, pepper=0.0): out = img.copy() mask = np.random.rand(*img.shape) < prob val = np.random.rand(*img.shape) out[mask] = np.where(val[mask] < 0.5, pepper, salt) return out np.random.seed(2024) noisy = salt_pepper(img, 0.2)

评价指标用PSNR和SSIM。PSNR高不代表视觉一定好，最好把降噪前后的局部放大图贴出来比对。为了公平，指标都在同样坐标范围[0,1]下计算。

5.3 对照方案：中值滤波与小波软阈值

中值滤波直接用scipy.ndimage的median_filter。注意scipy.signal里的medfilt2d在新版本Scipy中已经越来越边缘化，用ndimage版本更稳。

from scipy.ndimage import median_filter res_med = np.clip(median_filter(noisy, size=3), 0, 1)

一次小波软阈值去噪的实现，先定义软阈值函数：

def soft_threshold(x, thr): return np.sign(x) * (np.abs(x) - thr).clip(min=0) def wavelet_soft_threshold(noisy, thr=0.1, level=3, wavelet='db2'): coeffs = pywt.wavedec2(noisy, wavelet, level=level, mode='periodization') new_coeffs = [coeffs[0]] # 低频保留，不收缩 for i in range(1, len(coeffs)): cH, cV, cD = coeffs[i] new_coeffs.append((soft_threshold(cH, thr), soft_threshold(cV, thr), soft_threshold(cD, thr))) return pywt.waverec2(new_coeffs, wavelet, mode='periodization')

这个对比的意义在于：把小波软阈值和ISTA放在一起看，才能看出“迭代优化”带来的增益。一次性软阈值就像只看一眼图纸就下锤子，ISTA则是不断测量、不断修正。

5.4 核心代码：ISTA迭代去噪

下面是ISTA的核心实现。我在系数域做迭代，低频近似分量只做梯度更新、不做软阈值收缩，这样能保住图像的整体亮度，避免细节被过度压掉。

def ista_denoise(noisy, lambd=0.1, level=3, wavelet='db2', iterations=200, alpha=1.0): coeffs = pywt.wavedec2(noisy, wavelet, level=level, mode='periodization') w = [c.copy() for c in coeffs] for _ in range(iterations): # 1. 由当前小波系数重构图像 x = pywt.waverec2(w, wavelet, mode='periodization') # 2. 计算残差 residual = x - noisy # 3. 残差的小波变换就是目标函数光滑部分的梯度 grad = pywt.wavedec2(residual, wavelet, level=level, mode='periodization') # 4. 梯度步 + 软阈值收缩 w_new = [w[0] - alpha * grad[0]] # 低频只走梯度步，不收缩 thr = alpha * lambd for i in range(1, len(w)): if isinstance(w[i], tuple): cH, cV, cD = w[i] gH, gV, gD = grad[i] w_new.append((soft_threshold(cH - alpha * gH, thr), soft_threshold(cV - alpha * gV, thr), soft_threshold(cD - alpha * gD, thr))) else: w_new.append(soft_threshold(w[i] - alpha * grad[i], thr)) w = w_new return pywt.waverec2(w, wavelet, mode='periodization')

调用方式：

result_ista = ista_denoise(noisy, lambd=0.1, level=3, wavelet='db2', iterations=200)

整个流程就四步：读图归一化、加噪声、三套方法各自去噪、算PSNR和SSIM对比。如果想观察收敛过程，可以在循环里每隔一定次数存一次重构结果，画成一行小图。

5.5 实验结果对比与效果分析

我实测的camera 512x512灰度图，噪声密度0.2，db2三层小波，lambda=0.1，ISTA迭代200次，结果大致如下：

从这张表能看出三个关键结论：第一，同样是软阈值思路，迭代的ISTA比一次性软阈值能高出2.4dB左右，视觉上细节轮廓更干净；第二，中值滤波在中低噪声密度下并不差，但纹理区域糊掉的问题在SSIM上体现得很清楚；第三，高密度脉冲噪声下，比如密度拉到0.5以上，ISTA的优势更明显，中值滤波基本报废。

操作时可以把残差图也画出来：残留的亮点就是算法没处理干净的地方。对ISTA来说，残差里剩下的大多是随机分布的孤立点，细节结构处的分布已经很少，这是收敛良好的信号。

6. 参数调优与踩坑记录

6.1 lambda、分解层数与迭代次数的影响

先说lambda。它是最影响效果的超参数。我对同一个含噪图分别取lambda=0.03、0.1、0.3做了200轮ISTA，PSNR结果如下：

lambda偏小时，软阈值压不住脉冲造成的大系数；lambda偏大时，真实边缘系数也被压掉，图像看起来有“塑料感”。我调参的顺序是：先固定level=3、迭代200次，lambda从0.02到0.5按对数网格试，选出PSNR最优的，再微调level和迭代次数。不要一上来就同时调三个参数，容易陷入一团乱麻。

分解层数level建议取3到4。层数太少，低频部分太粗糙，收缩不到位；层数太多，最高频子带分辨率低，噪声和边缘混在一起难以分离。对512x512的图，3层基本够用，4层也行但计算量不划算。选小波基时，db2、db4、sym4在图像去噪里的效果差异不大，但db1也就是Haar小波会带来明显的方块伪影，建议避开。

迭代次数也有讲究。我的经验是50轮已经能出肉眼可用的结果，200轮基本收敛，500轮以上边际收益很小。可以每隔20轮打印一次PSNR，观察曲线是否还在涨。如果涨得极慢，说明alpha偏小；如果曲线先涨后跌，说明alpha太大，已经被噪声主导了。

6.2 边界效应与大梯度伪影

用小波处理图像时，边界是最容易出问题的地方。图像左右边缘在DWT里会跟对侧边界发生卷积混合，如果延拓模式不合理，边缘处会产生一圈亮边或暗带。我一般用mode='periodization'，它对周期信号是完美的，对自然图像也足够稳。如果发现边缘伪影明显，可以切到mode='reflect'或'symmetric'对比一下。

另一个问题是迭代过程中高频子带的“振铃”。软阈值收缩系数后，重构图像在强边缘处可能出现类似吉布斯现象的抖动。这个问题的本质是：L1稀疏正则倾向把小的边缘系数误杀，边缘处能量不连续，重构时产生振荡。缓解方法有两个：一是选择更平滑的小波基，比如sym4或coif2；二是给最低频的近似分量也设一个很小的阈值，避免低频里残留噪声被放大。当然，基础版本保留低频不动，在大多数情况下已经够用。

顺带提一个调试技巧：把ISTA每次迭代产生的图像按帧存下来，合成一张横向拼接的进度图，能直观看到“先整体逼近、后细节收缩”的过程。这不是花架子，它帮你判断算法到底在哪一步开始过拟合噪声，对调参帮助非常大。

6.3 扩展方向：联合TV正则、ADMM与工程落地

ISTA这套“小波变换+梯度下降”框架是稀疏去噪的基地型结构，真正落地时通常会再叠几层。

一是联合TV正则。小波L1正则擅长稀疏化系数，但不直接约束图像梯度；加一个TV项可以进一步压制方块伪影。代价是目标函数多了不可导项，软阈值不动，TV子问题要额外用投影算子处理，复杂度直接上一个台阶。

二是把近端梯度换成ADMM。ISTA处理“小波域稀疏+像素域保真”时，如果小波基不是正交的，迭代到后期偶尔会出现收敛抖动；ADMM通过引入分裂变量，把原来耦合的问题变成交替求解子问题，稳健性更好。实现难度更高，但工程上经常用来替换ISTA，尤其适合小波基不完全正交、或者需要叠加多个正则项的场合。

三是应用场景延伸。我做图像增强的预清洗时，会先用这套去噪管线处理低照度图上的脉冲点，再去提亮增强，效果比直接增强再修噪点干净得多。做水印鲁棒性测试时，在嵌入前对载体做一遍脉冲噪声去噪，也可以模拟真实传输链路里的信道误码情况。如果做视频去噪，把单帧ISTA改成“时间维滑窗+空间维小波稀疏”也能扩展。

我个人的体会是：小波变换和梯度下降的组合，核心价值不在某个具体函数，而在于它教会你用“先验+优化”的眼光去拆问题。脉冲噪声只是这个问题的一个载体，你掌握了建模、正则化、近端算子和调参这套动作，换成斑点噪声、条纹噪声、甚至缺失像素修复，路数都是一样的。最后再多说一句：把这套方法用到生产环境时，宁可把lambda设得稍小一点，多跑几十轮迭代，也不要追求一步到位把参数拧得冒烟；稳定、可解释、能随时翻回去调参，才是工程里真正重要的东西。