ADMM图像去噪实战：Plug-and-Play框架与MATLAB实现解析

简介：这是一份面向图像处理学习者和科研人员的ADMM图像去噪MATLAB源码包，围绕交替方向乘子方法在图像去噪与去模糊中的应用展开，适合希望掌握优化算法落地实践的读者。资源共19个文件，以15个m脚本为主，涵盖总变分去噪、BM3D、NLM等经典去噪器，以及PSNR计算、图像上下采样、投影函数等配套工具，另含测试图片与效果对比图，便于直接运行验证。压缩包整体仅154KB，结构紧凑，代码模块化程度高，可快速移植到医学影像、遥感图像或视频去噪等场景。已有1056人学习下载，说明其内容对相关方向研究具有参考价值。通过源码和示例，读者能直观理解ADMM的变量更新、拉格朗日乘子迭代与松弛项调整流程，并对照demo_deblur.m等演示脚本，快速搭建自己的去噪实验流程。

1. ADMM 图像去噪：为什么这个优化器比直接滤波更值得调试

拿到一张带噪图像，第一反应通常是中值滤波、高斯滤波或者小波软阈值。这些方法速度快，但遇到纹理密集或边缘锐利的图像时，很容易把细节一并抹掉。ADMM（Alternating Direction Method of Multipliers，交替方向乘子法）走的是另一条路：它把去噪问题拆成一个数据保真项和一个先验正则项，然后交替求解，让两项在迭代中互相妥协。直观地说，它的输出不是“滤”出来的，而是“解”出来的——噪声被建模成优化目标的一部分。这就意味着你可以在同一套框架下换不同的先验，TV、BM3D、NLM，甚至学出来的去噪网络，而迭代逻辑基本不动。适合想深入理解图像反问题、或者在做去模糊和超分辨时需要稳定求解器的人。下面这份 MATLAB 源码包里的PlugPlayADMM_deblur.m和四个wrapper_*文件，正好演示了这套思路。

2. ADMM 迭代框架与 Plug-and-Play 先验：从变量分裂到子问题求解

2.1 变量分裂：为什么要把一个目标函数拆成两个变量

经典的图像去噪或去模糊问题可以写成：

min_x (1/2) || A x - y ||_2^2 + λ R(x)

其中A是观测矩阵，在纯去噪时A是单位阵；在去模糊时A是卷积算子。R(x)是正则项，TV 的话是||D x||_1，BM3D 的话甚至没有一个显式的数学表达式。问题在于，A和R耦合在一起，直接求导或梯度下降都不方便。ADMM 的核心操作是引入辅助变量v，把问题改写成带等式约束的形式：

min_{x, v} (1/2) || A x - y ||_2^2 + λ R(v) s.t. x = v

这样一来，目标函数变成了两个变量块，可以分别求解。这就是变量分裂（variable splitting）。对纯去噪来说，A = I，第一个子问题实际上是一个简单的二次最小化，闭式解就是一个矩阵运算；而第二个子问题只剩正则项R(v)，这一步本质上是在做一个“去噪”操作。Plug-and-Play 的关键洞察就在这里：如果prox_{λR}(z) = argmin_v (1/2)||v - z||_2^2 + λ R(v)这个近端算子可以看作一个去噪器，那么任何现成的去噪算法——BM3D、NLM、甚至深度学习模型——都能直接塞进 ADMM 迭代中，而不用推导它们对应的代价函数。

2.2 三个核心更新步骤与残差收敛判据

ADMM 的每一次迭代按以下三步执行。这里给出一段标准的 MATLAB 骨架，对应源码里PlugPlayADMM_deblur.m的核心循环部分：

for k = 1:maxit % 子问题1: 求解 x，相当于保真项的最小化 % A 是观测矩阵，AT 是 A 的转置，rho 是惩罚参数 x_old = x; x = (AT*y + rho*(v - u)) / (AT*A + rho); % 子问题2: v 更新，等价于用去噪器处理 (x + u) v = denoiser(x + u, sigma); % 对偶变量 u 更新：拉格朗日乘子的梯度上升 u = u + (x - v); % 残差监测 r_prim = norm(x - v, 'fro'); s_dual = rho * norm(v - v_old, 'fro'); if r_prim < tol && s_dual < tol break; end end

逻辑说明：第一行x的更新是在固定辅助变量v和对偶变量u的情况下，求解一个带惩罚项的二次问题。rho越大，x越接近v - u，保真项的影响相对减弱。第二行是 Plug-and-Play 的关键——把x + u输入任意去噪器denoiser，sigma是该去噪器对应的噪声标准差估计。第三行更新对偶变量u，等价于把违反等式约束的程度反馈回下一步迭代。残差r_prim是原始残差，s_dual是对偶残差，两者都收敛才代表真正收敛。

参数说明：maxit一般取 30 到 100 即可，去模糊任务可能需要更多次迭代。rho的典型范围是 0.01 到 10，直接决定收敛速度和解的平滑度。tol设成1e-4到1e-5比较合适，太小会增加不必要计算，太大会提前停止导致伪影残留。上面代码里没有除以矩阵规模，实际判断残差时建议用norm(x - v, 'fro') / norm(y, 'fro') < tol，这样对不同尺寸图像更鲁棒。

2.3 去噪器即近端算子：TV、BM3D、NLM 和 RF 的统一接口

源码包里有四个 wrapper 文件：wrapper_TV.m、wrapper_BM3D.m、wrapper_NLM.m和wrapper_RF.m。它们的作用是把不同的去噪算法包装成统一的denoiser(z, sigma)接口。以wrapper_TV.m为例，它的内部通常调用一个求解全变分去噪的子程序，比如通过 Chambolle 投影算法或 FISTA 求解min_v (1/2)||v - z||_2^2 + lambda * TV(v)。wrapper_BM3D.m则内部调用 BM3D 库，wrapper_NLM.m调用非局部均值方法，wrapper_RF.m可能对应随机场或某种滤波。

这样做的工程价值是：你不需要为每种先验重写整个 ADMM 循环，只需要把denoiser这行换成不同 wrapper。而且denoiser的输入z是x + u，这个量在迭代中会逐渐逼近真实图像，但又始终带着噪声残留，所以 wrapper 内部的sigma参数要跟着迭代动态调整或保持固定。源码包中常见做法是固定sigma，理论分析表明 Plug-and-Play ADMM 在固定近端算子下仍能收敛到不动点，只要该去噪器满足适当的 Lipschitz 条件。

3. MATLAB 源码结构拆解：wrapper_TV.m 到 PlugPlayADMM_deblur.m 的调用链

3.1 压缩包目录与每个脚本的职责

拿到压缩包后，先看目录结构再动手，不然容易在downsample2.m和upsample2.m这两个文件上迷惑。以下是源码包中的关键文件清单及职责：

注意utilities里的constructGGt.m和defGGt.m，这两个函数主要处理卷积算子的 Gram 矩阵。在A是卷积模糊核时，AT*A对应频域里的特征值，constructGGt.m构建循环卷积矩阵的频域表示，defGGt.m可能是定义其快速算子。shepard_initialize.m用于初始化，afun.m则是定义A的线性算子接口，方便传入迭代求解器。proj.m可能用于投影约束，比如非负约束。

3.2 从 demo_deblur.m 到主循环的调用关系

我们先看demo_deblur.m的典型调用方式。这个脚本不会直接调用 ADMM，而是通过一组参数配置交给PlugPlayADMM_deblur.m。常见配置如下：

% demo_deblur.m img = im2double(imread('data/test_256.png')); h = fspecial('gaussian', [9 9], 3); % 模糊核 y = imfilter(img, h, 'circular'); y = y + 0.01 * randn(size(img)); % 加噪 opt.method = 'TV'; % 可选 'BM3D' / 'NLM' / 'RF' opt.rho = 0.15; opt.sigma = 1/255; % 去噪器内部噪声参数 opt.maxit = 60; opt.tol = 1e-4; [x_out, psnr_curve] = PlugPlayADMM_deblur(y, h, opt);

PlugPlayADMM_deblur.m内部首先根据opt.method选择对应的 wrapper，把函数句柄赋给denoiser变量。然后在迭代循环中，x子问题的求解涉及afun、constructGGt和defGGt。如果A是归一化的卷积算子，AT*A可以用defGGt快速计算，避免显式存储矩阵。这就是源码中downsample2.m和upsample2.m出现的原因——某些多尺度实现里，先下采样再卷积，再上采样恢复原始尺寸，shepard_initialize.m则是用来做初始化插值。

3.3 为什么需要 Utilities 里的频域算子加速

直接构造A矩阵在 256x256 图像上有 65536 维，显式矩阵是 65536×65536，根本存不下。所以afun.m采用“算子即函数”的思路，afun(x, mode)根据mode返回A*x或A'*y。在循环卷积假设下，A可以在频域对角化，constructGGt.m构造AT*A的频域对角元素，然后x子问题的解变为：

% 用频域方法快速求解 x 子问题 num = AT*y + rho * (v - u); den = GGt + rho; % GGt 是 AT*A 的频域特征 x = real(ifft2(fft2(reshape(num, sz)) ./ den));

这里GGt来自constructGGt.m的输出，den在每个频点上是标量，因此逐点相除即完成求解。downsample2.m和upsample2.m在这里的作用是处理多尺度模糊核或边界条件，如果模糊核不是循环卷积，而是空间变化或降采样后的算子，需要这两步把图像和核对齐到同一尺度。实际运行时会发现，这套频域加速方案比直接解线性方程快一到两个数量级，这是值得保留的工程细节。

4. 关键参数与去噪效果验证：psnr.m、test_256.png 与真实噪声场景

4.1 PSNR 计算与迭代曲线绘制

源码包里的psnr.m是一个实用函数，用来评估恢复图像与原始图像的峰值信噪比。它的公式是：

PSNR = 10 * log10(MAX^2 / MSE)

其中MAX是像素最大值，对灰度图通常是 1（双精度）或 255（uint8）。为了验证算法是否在迭代中确实改善图像，我通常会在PlugPlayADMM_deblur.m的循环里记录每一步的 PSNR，然后画出收敛曲线：

% 在 ADMM 循环中记录每步的 PSNR true_img = img; % 仅在测试时使用 for k = 1:maxit % ... 迭代更新 ... psnr_curve(k) = psnr(x, true_img); end plot(1:maxit, psnr_curve); xlabel('Iteration'); ylabel('PSNR (dB)'); grid on;

逻辑说明：这个曲线能直观反映两个问题——算法是否收敛，以及是否过拟合。如果 PSNR 先上升后缓慢下降，说明rho或sigma设置得过大，导致迭代后期往噪声上拟合。如果 PSNR 单调上升但迟迟不平稳，说明maxit不够或者tol太紧。psnr函数在 MATLAB 的 Image Processing Toolbox 或自带的psnr.m中都可用，但这个包里的版本是自己实现的，输入输出可能与官方略有差异，建议查看源码确认输入顺序是psnr(original, noisy)还是反过来。

4.2 参数rho、sigma和maxit的交互影响

这三个参数不是独立的。rho是惩罚参数，控制变量分裂的刚性；sigma是辅助变量更新时去噪器的强度参数。两者有直接的耦合关系。下面是一组在 256x256 灰度图上、模糊核为 9x9 高斯核时常用的参数组合：

注意sigma的取值范围取决于图像像素值是否归一化到[0,1]。若图像是 uint8，sigma需要除以 255，否则去噪强度会被低估。图中result.jpg是一个参考输出，我建议你跑完自己的配置后对比一下，差异如果超过 1 dB，优先检查rho和sigma的量级是否匹配。

4.3 真实噪声场景下的自适应策略

原始图像data/test_256.png是无噪声的地面真值，通常用于模拟实验。真实场景中噪声水平未知，可以用一个简单的估计器——先对图像做小波分解，取最细对角线子带的绝对中位差作为sigma估计。常见做法是：

% 估计噪声标准差 function sigma = estimate_noise(img) if size(img,3) == 3 img = rgb2gray(img); end y = im2double(img); [cA, cH, cV, cD] = dwt2(y, 'db2'); sigma = median(abs(cD(:))) / 0.6745; % MAD 估计 end

这个估计值可以直接传给wrapper_TV或wrapper_BM3D。但在 Plug-and-Play ADMM 中，sigma同时承担了正则强度的角色，所以如果把估计的噪声标准差直接当sigma用，rho需要相应调大，否则恢复结果过于平滑。一个实用技巧是先用sigma_est跑 20 次迭代，看x的 PSNR 变化，如果迭代后期 PSNR 下降，就适当缩小sigma或增大rho。对于纹理丰富的图像，将sigma乘以 0.8 通常能保留更多细节。

5. 把 TV 换成 BM3D 或 NLM 之前要懂的收敛性检查技巧

5.1 不同去噪器的行为差异与停机准则调整

当opt.method从'TV'切换到'BM3D'或'NLM'时，ADMM 迭代的行为会发生显著变化。TV 去噪器是严格的近端算子，它对应的代价函数是凸的，所以 ADMM 有理论收敛保证。BM3D 和 NLM 不是简单近端算子，它们内部包含分组、变换、阈值和聚合等非线性操作，不满足标准近端映射的定义。实际表现是：PSNR 曲线可能在前几步快速上升，然后在小范围内波动，甚至出现周期性振荡。

我一般会在切换去噪器时进行两个检查。第一个是画原始残差和对偶残差的双纵轴图。如果对偶残差震荡而原始残差单调下降，通常说明rho取小了，尝试把rho放大 2 到 5 倍。第二个检查是固定denoiser的输入输出差异，即norm(v - z)的轨迹。对于稳定的近端算子，这个差异应该随迭代单调递减；对于非近端去噪器，它可能保持在一个平台。如果平台太高，说明sigma太大，去噪器过度修改辅助变量，导致x和v始终不能靠拢。

5.2 一个实用的自动调参技巧：残差比值监控

这里分享一个针对 Plug-and-Play ADMM 的调参技巧，适用所有 wrapper。在每次迭代时计算一个比值ratio(k) = r_prim / (rho * norm(v - x_old))，这个比值反映原始残差与对偶残差之间的相对大小。理想的收敛轨迹是ratio在 1 附近小幅摆动并最终趋于 0。如果ratio始终大于 5，说明rho太大，需要减小；如果始终小于 0.2，说明rho太小，需要增大。在 MATLAB 中实现：

% 在循环中输出残差比值 ratio = r_prim / (rho * norm(v - x_old, 'fro') + eps); fprintf('Iter %03d r_prim=%.2e ratio=%.2f\n', k, r_prim, ratio);

这个方法只用两次范数计算，代价极小，但能提前暴露参数失配。当我从wrapper_TV切换到wrapper_BM3D时，通常会发现同样rho下ratio翻了几倍，这正是因为 BM3D 的平滑强度与 TV 不同，所以需要重新调rho才能让两个残差平衡。最后提一个小细节：源码包里的wrapper_NLM.m可能会一次性计算整个图像块的权重矩阵，内存占用是随图像尺寸平方增长的，测大图之前先算一遍whos看变量大小，别等内存爆炸再回头。

本文还有配套的精品资源，点击获取