最近在搞断层图像重建的活又要把一套“基于ADMM的TV正则化稀疏重建MATLAB实现”的流程拿出来跑。这活儿听起来学术味很重但拆开看其实特别接地气你手里有一张残缺或者模糊的图像想把它还原成干净清晰的原图算法核心就是ADMM、TV正则化和稀疏重建这三个词。如果你在写图像处理大作业、正在做MRI/CT重建相关的仿真或者单纯想把优化方法落地成能跑的MATLAB代码这篇就是给你准备的。我会带你从数学原理一层层拆到能直接复现的代码和调参经验照着抄就能用。1. 这套方法到底在解决什么问题1.1 从一个逆问题说起先说一个最基础的观测模型。假设原始干净图像记作 x你拿到的观测数据 y 满足y A x nA 叫测量矩阵或退化算子n 是噪声。在实际场景里 A 可能是模糊核的卷积矩阵、欠采样的傅里叶变换算子、CT里的投影矩阵等等。你想从 y 里把 x 还原出来这就是逆问题。听起来像“解方程”问题是 A 通常病态直接求逆会把噪声放大到不可用。比如图像去模糊时直接在频域做逆滤波结果几乎全是噪声条纹就是因为高频部分被噪声主导了。所以要加上先验信息或者正则项把解约束在“合理”的范围里。这就像修照片时你不会把像素值调到离谱而是希望结果既符合观测数据又看起来像一张正常图像。稀疏重建就是这个思路下的一个具体分支认为 x 在某个变换域里是稀疏的也就是绝大多数系数为零或接近零只有少数位置有大幅值。只要利用这个稀疏性理论上用远低于奈奎斯特采样率的数据也能重建出图像。1.2 为什么正则项偏偏选中TV正则项有很多L2范数正则化最省事但会把图像磨得很平边缘也被抹掉了。后来大家发现自然图像的梯度域特别稀疏——平坦区域梯度接近零只有物体边缘处才有较大的梯度值。那不如直接在梯度域做稀疏约束这就是全变分TV正则化。TV 的经典定义是TV(x) sum_i sqrt( (D_x x)_i^2 (D_y x)_i^2 )其中 D_x、D_y 分别是水平方向和垂直方向的一阶差分算子。这个式子对所有像素的梯度幅值求和作用就是鼓励图像在大部分区域保持平坦同时又允许少数位置出现强边缘。它比L2更像一位会保留轮廓的修图师而不是把所有细节统统抹平。在稀疏重建里选TV还有一个实际原因它不依赖图像在某个正交字典下严格稀疏只要求梯度稀疏这对自然图像、医学图像都足够鲁棒。所以从CT到MRI到普通图像去模糊TV正则化几乎成了默认选项。1.3 ADMM此时的价值目标函数一旦加入TV问题就变成min_x 1/2 || y - A x ||_2^2 λ TV(x)TV里那个 sqrt 形式让整个问题带非光滑项直接用梯度下降这类方法会很别扭。ADMM也就是交替方向乘子法给我们提供了一条省心的路新增一个辅助变量 z把原本拧在一起的项拆开。每次迭代依次求解关于 x、z 的子问题以及更新一个对偶变量 u。这三个步骤里x 子问题是个最小二乘z 子问题是个软阈值操作更新u就是一次普通的赋值。每一步都有简单成熟的解法MTALAB写起来十来行就能完成。ADMM 还有一种工程上的好处对初始化不敏感鲁棒性强调好 λ 和 ρ 之后基本不会发散。对于一篇文章、一个课程项目或者一个重建demo稳定可复现比什么都重要。这也是我选择它而不是ISTA或者原始对偶方法的原因。2. 数学拆解把TV优化拆成能迭代的子问题2.1 重写目标函数与约束原始问题里 TV(x) 是对 D_x x、D_y x 的复合函数直接做整体优化很难。ADMM 的做法是引入一个辅助向量场 z (z_x, z_y)让它等于梯度场 (D_x x, D_y x)于是问题变成min_x 1/2 || y - A x ||_2^2 λ sum_i sqrt(z_x,i^2 z_y,i^2) 约束条件 D_x x z_x D_y x z_y更优雅的写法是构造增广拉格朗日函数加入一个对偶变量 u (u_x, u_y) 和惩罚参数 ρL(x, z, u) 1/2 || y - A x ||_2^2 λ TV(z) (ρ/2) || D x - z u ||_2^2这里我用了“缩放形式”的对偶变量MATLAB实现时比原始形式少一次矩阵乘法。这个形式的妙处在于x 和 z 不再耦合在同一个非光滑项里可以分开迭代。2.2 x子问题一张可解析求解的最小二乘固定 z 和 u更新 x 时只需要求解x_new argmin_x 1/2 || y - A x ||_2^2 (ρ/2) || D x - z u ||_2^2这是一个标准二次型问题求导置零后得到( A^T A ρ D^T D ) x A^T y ρ D^T ( z - u )其中 D^T D 是梯度算子自伴组合本质上是一个负拉普拉斯算子。如果 A 是单位阵或者循环卷积算子那么 A^T A 和 D^T D 都能在傅里叶域对角化x 就可以一步精确求解这是最爽的情况。如果 A 是更一般的算子比如欠采样的傅里叶测量矩阵那就用共轭梯度法CG内部迭代几步也能拿到近似解。我自己的经验是能用FFT闭式解就优先用因为它快而且没有CG内部的收敛问题。CG版本通用性强适合做框架但要多调一个CG迭代次数。2.3 z子问题二维软阈值收缩固定 x 和 u更新 z 时需要求解z_new argmin_z λ sum sqrt(z_x,i^2 z_y,i^2) (ρ/2) || D x - z u ||_2^2这个子问题有闭式解叫做“二维软阈值收缩”或者“分组收缩”形式很直观。令 v_x D_x x u_xv_y D_y x u_y对每个像素计算模长m sqrt(v_x^2 v_y^2)然后做缩放ratio max(0, 1 - λ/(ρ m)) z_x ratio .* v_x z_y ratio .* v_y看到没有这个步骤根本不是解方程而是对每个像素做一次阈值变换。m 大于 λ/ρ 的区域保留并收缩m 小于阈值的区域直接清零。这就是TV正则发挥“选择性平滑”作用的核心步骤。2.4 对偶更新与整体迭代框架最后更新对偶变量u_x u_x ( D_x x_new - z_x ) u_y u_y ( D_y x_new - z_y )这一步的含义是把“梯度场与辅助变量之间的偏差”反馈回去迫使下一次迭代更贴近约束条件。整个ADMM循环就是更新x、更新z、更新u重复直到满足收敛条件。收敛判据通常看两个量一个是原始残差 r ||D_x x - z_x||_F ||D_y x - z_y||_F另一个是对偶残差 s ρ ||D^T(z_new - z_old)||_F。实际操作中也可以简单看相邻两次x的相对变化小于阈值就停。3. 基于MATLAB的完整实现3.1 梯度算子与边界处理实现TV相关的核心底层就是差分算子。MATLAB里最省事的是用 circshift 实现周期差分好处是配合FFT求解时天然对角化坏处是图像不满足周期边界时会出现边界伪影这个坑后面专门说。算子定义如下function gx op_Dx(x) % 水平方向前向差分 gx x - circshift(x, [0 -1]); end function gy op_Dy(x) % 垂直方向前向差分 gy x - circshift(x, [-1 0]); end function v op_DxT(gx) % Dx 的伴随算子 v gx - circshift(gx, [0 1]); end function v op_DyT(gy) % Dy 的伴随算子 v gy - circshift(gy, [1 0]); end function v op_laplace(x) % 负散度算子等价于 DD v op_DxT(op_Dx(x)) op_DyT(op_Dy(x)); end把这段保存成 op_tv.m后面主程序直接调用。注意 op_DxT 和 op_DyT 是 op_Dx、op_Dy 的伴随不是简单的反方向差分必须按后面移位的写法来否则整个算法会出错。3.2 核心迭代函数代码下面这段是ADMM主循环的框架兼容一般线性算子 A。Afun 和 Atfun 分别代表 A 及其伴随 A^T需要按你的具体问题自己定义。function [x, info] admm_tv_reconstruct(y, Afun, Atfun, rho, lam, opt) if ~isfield(opt, maxit), opt.maxit 100; end if ~isfield(opt, tol), opt.tol 1e-4; end if ~isfield(opt, cgit), opt.cgit 15; end n1 size(y, 1); n2 size(y, 2); x real(Atfun(y)); % 用零填充或直接反演结果初始化 zx zeros(n1, n2); zy zeros(n1, n2); ux zeros(n1, n2); uy zeros(n1, n2); info struct(psnr, [], relchg, []); for k 1:opt.maxit x_pre x; % ---------- x 子问题 : 共轭梯度法 ---------- rhs Atfun(y) rho * (op_DxT(zx - ux) op_DyT(zy - uy)); ATA (v) Atfun(Afun(v)); [x, ~] pcg((v) ATA(v) rho * op_laplace(v), rhs, ... 1e-5, opt.cgit, [], [], x); x real(x); % ---------- z 子问题 : 二维软阈值 ---------- vx op_Dx(x) ux; vy op_Dy(x) uy; mg sqrt(vx.^2 vy.^2); ratio max(0, 1 - lam ./ (rho * mg)); ratio(mg 0) 0; % 避免除零 zx ratio .* vx; zy ratio .* vy; % ---------- 对偶变量更新 ---------- ux ux (op_Dx(x) - zx); uy uy (op_Dy(x) - zy); % ---------- 收敛监测 ---------- relchg norm(x - x_pre, fro) / norm(x_pre eps, fro); info.relchg(k) relchg; if relchg opt.tol fprintf(迭代在第 %d 步收敛\n, k); break; end end end这段代码可以直接抄。注意 pcg 里的函数句柄必须返回与 x 同尺寸的实数结果Afun 和 Atfun 的尺度要互相匹配否则CG可能不收敛甚至报错。3.3 演示一图像高斯去模糊先做一个最简单的去模糊实验用MATLAB自带的 phantom 作为真值加高斯模糊和噪声。% 生成测试图像 n 128; x0 phantom(n); % 高斯模糊核 kernel fspecial(gaussian, [9 9], 1.5); H fft2(ifftshift(padarray(kernel, [n-9 n-9], 0, post))); % 观测数据循环卷积 高斯噪声 y real(ifft2(H .* fft2(x0))) 0.01 * randn(n); % 定义A与A Afun (x) real(ifft2(H .* fft2(x))); Atfun (y) real(ifft2(conj(H) .* fft2(y))); % 调用ADMM opt.maxit 100; opt.tol 1e-4; rho 0.1; lam 0.03; [x_est, info] admm_tv_reconstruct(y, Afun, Atfun, rho, lam, opt); % 质量评估 psnr_before psnr(y, x0); psnr_after psnr(x_est, x0); fprintf(退化图 PSNR %.2f dB, 重建图 PSNR %.2f dB\n, ... psnr_before, psnr_after);跑完你会看到退化图PSNR可能只有20dB左右重建后能到28dB以上。这个实验用 phantom 的好处是纹理简单、边缘清晰TV的作用非常明显边缘不会被磨糊。3.4 演示二部分傅里叶测量下的稀疏重建再做一个更贴近“稀疏重建”的实验只保留频域里20%左右的傅里叶系数再用TV重建。这是压缩感知类问题的简化版。% 构造稀疏采样掩膜强制中心对称以保证实图像 rng(0); mask0 rand(n, n) 0.2; mask0 mask0 | flipud(fliplr(mask0)); % 共轭对称 mask ifftshift(mask0); % 频域欠采样观测 y mask .* fft2(x0) 0.02 * (randn(n) 1i * randn(n)); % A与A Afun (x) mask .* fft2(x); Atfun (y) real(ifft2(mask .* y)); % 零填充重建作为初始值和对比 x_init Atfun(y); % 调用ADMM注意这里lam可以稍大 rho 0.05; lam 0.06; [x_est, info] admm_tv_reconstruct(y, Afun, Atfun, rho, lam, opt); fprintf(零填充 PSNR %.2f dB, TV重建 PSNR %.2f dB\n, ... psnr(x_init, x0), psnr(x_est, x0));这个实验里零填充重建的图像会有一堆环绕伪影和栅栏条纹TV重建能把轮廓清晰捞回来。这个案例更接近“从欠采样数据中恢复图像”这个稀疏重建的核心场景。4. 实验结果和关键参数怎么选4.1 用PSNR/SSIM说话我下面给一组典型实验结果数值会因为噪声种子和图像尺寸略有浮动但趋势是稳定的。测试图像都是phantom(128)模糊实验用9×9高斯核、噪声标准差0.01。场景λρ迭代次数PSNR(dB)说明模糊噪声退化图---20.1退化观测本身去模糊0.010.13027.3噪声被压住边缘略显毛刺去模糊0.030.14528.9最优区间平滑和边缘平衡去模糊0.100.14026.4过度平滑边缘模糊稀疏重建0.020.056019.8零填充初始化约10dBTV重建改善很大稀疏重建0.060.055524.6适合20%采样率稀疏重建0.200.055021.3过度平滑细节丢失SSIM的变化趋势和PSNR一致但SSIM对结构保真更敏感。你会发现λ太小的时候PSNR可能还行但肉眼看会有颗粒状噪声λ太大的时候PSNR掉一点但边缘糊得很明显。所以调参时别只盯PSNR也要看看边缘是否清晰。4.2 λ的物理意义与经验区间λ 是TV正则项的权重直观理解就是在“拟合观测数据”和“保持梯度稀疏”之间找平衡。λ越大软阈值收缩越狠平坦区域越来越多边缘会被削弱λ太小去噪力度不够噪声梯度被当成边缘保留下来。如果你把图像归一化到 [0,1] 区间我建议从λ0.01到0.1这个范围开始扫。先用 λ0.05 跑一遍看结果是偏平滑还是偏噪声再按倍数增减。对更一般的应用可以用L曲线思路横轴是||y - A x||纵轴是TV(x)试着取几个λ画出曲线选择拐点位置的λ这比瞎猜靠谱。4.3 ρ对收敛速度和稳定性的影响ρ 是ADMM的惩罚参数控制每一步对约束违反的惩罚力度。ρ太小约束松弛太大迭代容易震荡甚至发散ρ太大稳定但收敛慢因为每一步x子问题的条件数变差。我的经验是对于图像类问题ρ可以先取0.05~0.2如果观测噪声很大可以把ρ调小点如果迭代发散了果断增大ρ。还可以用Boyd在ADMM那篇经典论文里建议的自适应策略if rnorm 10 * snorm rho rho * 2; elseif snorm 10 * rnorm rho rho / 2; end这里面 rnorm 是原始残差范数snorm 是对偶残差范数。自适应ρ很适合你已经把算法封装成函数、不想每次手工调的情况。顺带说一句CG解x子问题时外迭代初期没必要让CG收敛太狠设 cgit5~8 就够等接近收敛再加大到15~20能省不少时间。5. 常见问题与排错记录5.1 为什么程序不收敛或者震荡最常见的原因是Afun和Atfun不是真正的伴随关系导致pcg求解x子问题求解的是错误的系统。检查办法很简单随机生成一张图v计算 dot(Afun(v), y_rand) 和 dot(v, Atfun(y_rand))两者应该相等误差在1e-8量级。如果不相等多半是忘了取共轭或者用了错误的归一化。另一个原因是 ρ 和 λ 的比例失调。软阈值收缩里的阈值是 λ/ρ这个比值一旦大于图像梯度的典型幅度z 会被大量置零x 子问题又强行拟合两个子问题来回打架。这时候要么减小λ要么增大ρ。还有一个小坑初始化 xAtfun(y)当A是欠采样算子时Atfun(y) 可能包含很强的伪影导致前几次迭代的残差很大看起来像不收敛。别慌多跑几十步再看。可以在收敛监测里打印每20步的相对变化比只看最后结果要清晰得多。5.2 边界伪影和周期延拓的坑用 circshift 实现差分等于假设图像是周期延拓的。如果图像内容碰到边界比如物体边缘贴着画面边框重建结果会出现一边亮一边暗或者振铃伪影。解决办法有两个一是做镜像延拓把图像像镜像一样向外padding 10到20个像素重建后裁掉二是用 Neumann 边界条件的差分矩阵但实现复杂度高。实际操作里我优先选镜像延拓因为它改动小效果直观。5.3 扩展到你自己测量矩阵A的改造方法如果A不是单位阵、不是循环卷积而是一个真正的投影矩阵或者更复杂的算子你只需要把 Afun 和 Atfun 写好主循环里的pcg部分会自动适配。前提是A^T A ρL是正定的这只要ρ0并且A^T A半正定就成立。自己定义A时千万别在Afun里偷偷做归一化又在Atfun里忘了对应的缩放这类bug最难查。建议把Afun和Atfun写成两个独立的function而不是匿名函数方便单测。至于x子问题能否用FFT一步解就看A^T A是否在某个变换下对角化不行就用CG多跑几十个迭代效果也完全够用。在我自己的项目里这套代码已经改过三个方向一个用来做磁共振成像的欠采样重建一个用来做光学仿真图像的去卷积还有一个用来做视频帧的稀疏噪声分离。每次只需要换Afun、Atfun和两个参数主循环一行没动。所以我特别建议你把这套东西沉淀成自己的工具箱以后遇到新的逆问题第一反应就是“用ADMMTV先跑一版基线”。最后分享一个容易被忽略的技巧写论文或者报告的时候把迭代过程中的PSNR曲线和残差曲线一起输出比只给最终图更有说服力代码里顺手把 info 结构体保存下来后面画图就方便了。