简介这份MATLAB代码包围绕基于正则化的图像超分辨重建展开面向具备一定图像处理基础、希望理解并复现总变分正则化去卷积算法的高校学生与研究人员。它针对噪声与模糊造成的图像退化问题提供可运行的算法框架、测试图像与结果输出便于开展课程实验和小型研究验证。压缩包共22个文件约882KB以6个m脚本承载核心算法与示例流程另有14张jpg、2张png测试图像及输出结果方便直接运行、观察和对比。代码包含经典总变分正则化及L1、L2范数变体的去卷积实现分别侧重边缘保持、稀疏解与平滑效果同时提供去噪、去模糊示例和PSNR质量评价函数数据与输出目录便于复现实验、比较不同正则化策略的差异。已有1845人学习适合作为超分辨重建与图像恢复方向的入门实践和课程实验参考。1. 一张 200×150 的截图要放大四倍正则化超分到底在补什么手上只有一张 200×150 的监控截图要放大四倍放进报告里imresize(...,bicubic)出来的结果边缘发虚、车牌数字糊成一团——插值只是在已有像素之间猜中间值它不会凭空造出高频。正则化图像超分辨重建走的是另一条路先把「高分辨率图像怎么变成这张低分辨率图」写成一个退化方程再把这个方程加上一个先验约束变成一个可以迭代求解的最优化问题用先验把丢掉的高频“猜”回来。这套基于正则化的图像超分辨重建 matlab 代码把退化算子、Tikhonov 与 TV 两种先验、共轭梯度与 ADMM 两套求解器都拆成了可直接调用的函数。它适合已经会用 MATLAB 做图像处理、但被“插值放大不能用”卡住的人也适合想把课本上那几个公式真正跑出图、并且知道参数往哪调的人。2. 退化模型拆解为什么你的超分结果总是比双三次还软绝大多数人第一次写超分代码都是直接对低分辨率图做一次反卷积结果要么糊得更厉害要么满屏振铃。原因不在求解器在于模型少了一项。真实成像链路上高分辨率场景先被光学系统模糊再被传感器采样成低分辨率网格中间还叠了噪声——这三步少写任何一步解出来的东西都会偏向某个错误的方向。2.1 单帧超分是病态问题y DHx n 里四个符号的分工把观测写成y D H x n逐项拆开x是待求的高分辨率图像H是模糊算子由点扩散函数 PSF 决定的卷积D是下采样算子每隔 s 个像素取一个n是噪声。以 256×256 放大 4 倍为例y只有 64×64方程个数 4096未知数 65536欠定十几倍。更麻烦的是D H的零空间非空存在大量不同的x经过同样的模糊和下采样后得到完全一样的y。这就是病态ill-posed的来源。插值法本质上是人为规定「零空间里的分量全取零」所以高频一律丢失正则化则是换一条规则——在所有能解释y的x里挑一个让先验项最小的比如「梯度尽量稀疏」或「二阶导数尽量小」。选哪条规则决定了结果偏软还是偏锐。所以判断一份超分代码靠不靠谱第一个要看的就是它有没有显式建模H和D。只写y Hx而不管下采样那叫去模糊不叫超分把D写成简单裁剪而不做抗混叠放大后会出现明显的锯齿和摩尔纹。2.2 PSF 怎么定三种常见核与一个反直觉结论H的核心是 PSF。MATLAB 里最省事的构造方式是fspecial但选哪种、参数给多少直接决定结果锐度。核类型生成方式典型参数用错时的表现高斯核fspecial(gaussian, k, sigma)k9~13sigma1.2~2.0sigma 偏小 → 振铃偏大 → 过平滑运动模糊fspecial(motion, len, ang)len7~15ang 按实际方向方向估反 → 边缘出现斜向拖影圆盘/散焦fspecial(disk, r)r2~4半径偏小 → 高光处出现光晕环反直觉的地方在于PSF 估计偏小比偏大更容易翻车。核给小了模型认为观测到的模糊很轻求解器就会用很强的锐化去补偿高频被过度放大噪声和振铃同时爆出来视觉上比插值还难看。核给大一点结果只是软至少不发散。我一般的做法是先用一条竖直边缘估个初值——取边缘剖面拟合它的 10%→90% 上升宽度除以 2.56 就近似是高斯 sigma——然后在初值附近做一轮 0.8x 到 1.5x 的扫描看 PSNR 曲线。2.3 边界与伴随写不对就是发散没有第二种可能迭代求解器尤其是共轭梯度要求你提供的伴随算子Aᵀ严格满足内积恒等式⟨Ax, y⟩ ⟨x, Aᵀy⟩。很多人A用imfilter(..., symmetric)Aᵀ却随手写成imfilter(..., replicate)边界上内积对不上迭代会缓慢发散或者卡在一个明显错误的解上而残差曲线看起来还挺正常——这是最难查的一类坑。写完之后先跑一遍自检别急着看图% 伴随算子内积自检A 与 At 必须是严格的一对 x randn(N); y randn(size(A(x))); lhs sum(A(x).*y, all); rhs sum(x.*At(y), all); fprintf(内积相对误差 %.3e\n, abs(lhs - rhs) / abs(lhs));正常应该在1e-12量级。如果跑出来是1e-2别往下调参了先把伴随算子改对。对于对称核加对称边界填充imfilter得到的矩阵本身就是对称的Aᵀ可以直接复用同一段代码一旦换成circular或非对称核就必须显式做核翻转加裁剪这个细节后面第 5 章还会再提一次。2.4 正则项选型Tikhonov、TV、稀疏先验各自的天花板正则项不是越复杂越好它要和你的图像内容、噪声水平匹配。三个主流选项的取舍先验数学形式保边能力求解方式适用场景TikhonovL2λ‖Lx‖₂²弱边缘会被一起压平正规方程 共轭梯度闭式可解噪声大、只求稳定不求锐利全变分 TVL1λ‖∇x‖₁强能保住阶跃边缘ADMM / 原始-对偶需迭代纹理少、边缘硬的目标文字、车牌、建筑稀疏表示λ‖α‖₁, x Φα取决于字典字典学习 OMP/LARS有同域训练样本、纹理重复一个实用结论纹理丰富的自然图像TV 会把细节“油画化”草地、毛发这种地方会变成一块块色斑。这时候要么降 λ要么换非局部相似或稀疏先验。而文字、条形码、工业检测里的规则边缘TV 的表现可以明显压过 Tikhonov。选型不要只看论文里的 PSNR 表要拿你自己那类图跑一遍。3. Tikhonov 版本跑通从算子构造到共轭梯度求解的六十行先把最简单的版本跑通是整个流程里性价比最高的一步。Tikhonov 版本有闭式的正规方程能验证你的退化算子、伴随算子、数据对齐是不是都对等这一版输出的 PSNR 明显高于双三次再上 TV 才有意义。直接上 ADMM 而模型又是错的你会在几十个参数里迷路。3.1 目标函数与正规方程写成最小二乘加二次惩罚min_x 0.5‖DHx - y‖₂² (λ/2)‖Lx‖₂²对x求梯度令其为零得到正规方程(HᵀDᵀDH λLᵀL) x HᵀDᵀ y左端是一个巨大的对称正定矩阵但它是稀疏带状卷积 下采样 局部差分所以永远不要把矩阵显式建出来——256×256 的图矩阵就是 65536×65536双精度要 32 GB。正确做法是把它写成一个函数句柄让共轭梯度每次只做「乘一次」内存占用就是几张图的大小。这是 MATLAB 做这类问题最关键的一条工程习惯。L取离散拉普拉斯五点差分它的作用是把解往「二阶导数为零」的方向拉也就是偏向平滑。这也解释了为什么 Tikhonov 结果看起来偏软——先验本身就是奖励平滑的。3.2 造自测数据imresize造 LR 有个必须知道的偏差没有真实 HR 参照就没法算 PSNR所以自测阶段要用 HR 人工造 LR。这里有个坑imresize默认带抗混叠它的降采样滤波器和你代码里的D不是同一个东西。如果你的A里用的是纯隔点采样而 LR 是imresize造的模型和数据的退化不匹配PSNR 会莫名其妙卡在一个低值上怎么调 λ 都没用。两种处理方式任选一种但必须一致一是自测时把 LR 也用你自己的degrade()生成形成闭环二是正式跑真实图时把D换成带抗混叠的版本。前者适合验证代码正确性后者适合真实场景。3.3 完整脚本主流程加四个辅助函数%% demo_tikhonov.m —— 高斯模糊 4 倍下采样的 Tikhonov 超分 clear; clc; scale 4; psf fspecial(gaussian, 9, 1.6); % 9x9 高斯核sigma1.6 N [256 256]; % 高分辨率网格尺寸 lambda 0.02; % 正则化强度先给个小的 maxit 120; tol 1e-5; % 1) 造一对 HR/LR仅用于自测真实场景跳过此步 Xtrue double(imread(cameraman.tif)); Xtrue Xtrue(1:N(1), 1:N(2)); Y degrade(Xtrue, psf, scale); Y Y 1e-3 * randn(size(Y)); % 加一点点噪声别加多 % 2) 前向算子 A 与伴随算子 At都用函数句柄不建矩阵 A (x) degrade(x, psf, scale); At (y) adegrade(y, psf, scale, N); % 3) 正规方程左端 M AA lambda * LLL 为离散拉普拉斯 b At(Y); M (v) At(A(v)) lambda * lap(lap(v)); % 4) 共轭梯度初值取双三次上采样比零初值快一倍以上 x0 imresize(Y, scale, bicubic); [x, relres] cg_solve(M, b, x0, maxit, tol); fprintf(双三次 PSNR %.2f dB\n, psnr8(Xtrue, x0)); fprintf(Tikhonov PSNR %.2f dB\n, psnr8(Xtrue, x)); figure; imshow([x0, x], []); title(左双三次 右正则化重建);四个辅助函数全部放在同目录下function Y degrade(X, psf, s) % 前向退化对称边界卷积 - 隔点下采样 B imfilter(X, psf, symmetric, same, conv); Y B(1:s:end, 1:s:end); end function X adegrade(Y, psf, s, N) % 伴随退化零填充上采样 - 核翻转后卷积 - 裁剪回 HR 尺寸 Z zeros(N); Z(1:s:end, 1:s:end) Y; X imfilter(Z, rot90(psf, 2), symmetric, same, conv); X X(1:N(1), 1:N(2)); end function y lap(x) % 五点离散拉普拉斯对称边界核自身对称故可复用为 L K [0 1 0; 1 -4 1; 0 1 0]; y imfilter(x, K, symmetric, same, conv); end function [x, relres] cg_solve(Afun, b, x0, maxit, tol) % 标准共轭梯度Afun 必须对称正定 x x0; r b - Afun(x); p r; rs r(:)*r(:); rs0 rs; relres zeros(maxit,1); for k 1:maxit Ap Afun(p); alpha rs / (p(:)*Ap(:)); x x alpha*p; r r - alpha*Ap; rs_new r(:)*r(:); relres(k) sqrt(rs_new/rs0); if relres(k) tol, break; end p r (rs_new/rs)*p; rs rs_new; end relres relres(1:k); end function v psnr8(ref, est) % 统一按 8bit 动态范围算 PSNR mse mean((ref(:) - est(:)).^2); v 10*log10(255^2 / mse); endlambda 0.02不是普适值它和‖A‖²的量级强相关。粗调经验是先跑lambda 0此时退化为最小二乘解噪声被放大得很厉害然后按 10 倍步长往上试1e-3 / 1e-2 / 1e-1 / 1取 PSNR 最高那档再在该档附近做 2 倍细扫。maxit给 120 通常够relres曲线如果在 20 步内就降到1e-5说明lambda偏大、问题被正则项主导了。3.4 彩色图别三个通道各跑一遍走 YCbCr把 RGB 三个通道分别超分再拼回去是最常见的错误之一。三个通道各自的最优lambda不一样而且人眼对色度通道的分辨率远低于亮度通道分开跑不但慢三倍还会在边缘产生彩色镶边。正确做法是转 YCbCr只对 Y 通道做正则化重建Cb、Cr 直接双三次放大img double(imread(peppers.png)); ycc rgb2ycbcr(img ./ 255); % 归一化到 [0,1] 再转色空间 Ych ycc(:,:,1) * 255; % 亮度通道单独拿出来 % ... 对 Ych 做上面的重建得到 Yhr ... ycc_hr imresize(ycc, scale, bicubic); % 色度走插值 ycc_hr(:,:,1) min(max(Yhr,0),255) / 255; % 亮度用重建结果 out ycbcr2rgb(ycc_hr) * 255;注意rgb2ycbcr期望输入是[0,1]还是[0,255]取决于 MATLAB 版本和输入类型换算错了颜色会整体偏绿。跑完先imshow看一眼肤色是否正常比查文档快。4. TV 正则加 ADMM把边缘从“糊”拉回“锐”的三步迭代Tikhonov 跑通之后你会发现一个稳定的现象PSNR 涨了但视觉上还是有点闷尤其是物体轮廓。因为 L2 惩罚对大梯度罚得太狠先验在「压平边缘」。换成 L1 就完全不同——L1 对大梯度的惩罚是线性的只要边缘足够强它宁愿保留边缘、把其他地方的梯度压到零。4.1 从 L2 到 L1TV 保边的直觉解释TV 先验写成λ‖∇x‖₁其中∇x是水平和竖直两个方向的差分。它的效果可以这么理解把全图的梯度总量当成成本优化器会优先把「小幅度的、噪声引起的」梯度砍掉因为砍它们省成本而保留「大幅度的、真实边缘的」梯度砍它们收益低。这就是 TV 能去噪同时保边的原因。代价是阶梯效应平滑渐变区域会被压成一块块常数平台也就是前面提到的“油画感”。另外 TV 的‖·‖₁不可微正规方程那套闭式解不能用了必须换成能处理非光滑项的算法——ADMM 是最工程化的选择它的每一步都有闭式解。TV 还分各向异性|∇ₕx| |∇ᵥx|和各向同性√(∇ₕx² ∇ᵥx²)。前者对轴向边缘更锐后者对斜向边缘不产生锯齿。图像超分一般用各向同性下面的代码也按这个来。4.2 差分算子与它的转置一步写错全盘皆输∇用前向差分边界置零。它的转置是后向差分符号相反。这两个必须严格配对否则 ADMM 一定发散——而且发散的形式是缓慢漂移不是一下子 NaN很容易误判成「参数还没调好」。function d dh(x) % 水平前向差分最后一列置零 d [diff(x, 1, 2), zeros(size(x,1), 1)]; end function x dht(d) % 上面算子的转置 x zeros(size(d)); x(:,1) -d(:,1); x(:,2:end-1) d(:,1:end-2) - d(:,2:end-1); x(:,end) d(:,end-1); end function d dv(x) % 竖直前向差分最后一行置零 d [diff(x, 1, 1); zeros(1, size(x,2))]; end function x dvt(d) % 竖直差分的转置 x zeros(size(d)); x(1,:) -d(1,:); x(2:end-1,:) d(1:end-2,:) - d(2:end-1,:); x(end,:) d(end-1,:); end同样可以用内积自检sum(dh(x).*z,all)应该等于sum(x.*dht(z),all)。这一步花三十秒能省掉后面两小时的困惑。4.3 ADMM 三步更新x 子问题、z 子问题、对偶更新引入分裂变量z ∇x把约束优化转成增广拉格朗日x 步解 (AᵀA ρ∇ᵀ∇) x Aᵀy ρ∇ᵀ(z - u) z 步z soft(∇x u, λ/ρ) —— 逐像素软阈值 u 步u u ∇x - zx步又是一个对称正定系统直接复用第 3 章的cg_solve而且因为有了ρ∇ᵀ∇项条件数比原来更好内迭代 20~30 步就够。z步是逐点闭式解各向同性版本是“组软阈值”把(∇ₕ, ∇ᵥ)当成二维向量整体按模长收缩。%% tv_admm.m —— 接在第 3 章的数据与算子上 lambda 0.03; rho 0.8; % TV 权重 / 惩罚参数 outer 60; inner 25; % 外迭代 / 内 CG 步数 x imresize(Y, scale, bicubic); zh zeros(size(x)); zv zh; % 分裂变量 uh zh; uv zv; % 对偶变量缩放形式 AtY At(Y); for t 1:outer % --- x 子问题 --- rhs AtY rho * (dht(zh - uh) dvt(zv - uv)); Mx (v) At(A(v)) rho * (dht(dh(v)) dvt(dv(v))); x cg_solve(Mx, rhs, x, inner, 1e-4); % --- z 子问题各向同性软阈值 --- gh dh(x) uh; gv dv(x) uv; mag sqrt(gh.^2 gv.^2) eps; shrink max(0, 1 - (lambda/rho) ./ mag); zh gh .* shrink; zv gv .* shrink; % --- 对偶更新 --- uh uh dh(x) - zh; uv uv dv(x) - zv; if mod(t, 10) 0 fprintf(iter %2d PSNR %.2f dB\n, t, psnr8(Xtrue, x)); end end x min(max(x, 0), 255); % 出图前一定要裁剪到合法范围rho的物理含义是「约束z ∇x有多硬」它同时也是x子问题的步长尺度。经验值 0.5~2 之间太小则收敛慢50 次外迭代都未必稳太大则x子问题条件数变差、CG 内迭代不够。lambda/rho才是真正控制稀疏强度的量所以调参时先固定 rho只扫 lambda扫出量级之后再微调 rho。4.4 ρ 与 λ 的联合扫描别用单一 PSNR 做决策两个参数扫起来是二维的全网格扫太慢。实用流程是固定rho 0.8lambda取[0.005 0.01 0.02 0.05 0.1 0.2]六档各跑 40 次外迭代记录 PSNR 和耗时选出 PSNR 最高的一档再在它左右各取一档配合rho取[0.4 0.8 1.5]做个 3×3 的小网格。判断标准不要只看 PSNR。TV 的 PSNR 峰值往往出现在lambda稍大的位置那时候图已经很“平”了视觉上反而不好看。我的习惯是同时看三样PSNR、收敛曲线是否在 40 次内走平、以及放大 2 倍后边缘有没有出现阶梯块。三者都过关才定参数。5. 排查清单结果发灰、振铃、发散、比插值还差的五种现场这一节是从项目实践中总结出的五种常见问题。共同点是现象看起来像“参数没调好”实际根因往往在模型或实现层。按顺序排查不要一上来就狂调 λ。5.1 现象重建结果的 PSNR 反而低于双三次原因三种可能按概率排序——第一lambda太小甚至为 0最小二乘解把噪声放大;第二A与Aᵀ不满足内积恒等式CG 收敛到了一个错误的解;第三用于自测的 LR 是imresize生成的而A里用的是纯隔点采样退化模型与数据不匹配。解决先跑 2.3 节的内积自检误差不在1e-10量级就直接修算子;通过之后把lambda从1e-3按 10 倍往上试看 PSNR 曲线是否出现单峰;如果曲线一直单调上升或单调下降基本可以确定是模型不匹配把 LR 的生成方式换成自己的degrade()再测一次。5.2 现象图像四周出现一圈亮边或暗边越靠边越明显原因边界填充方式不一致。A用symmetricAᵀ用了默认的零填充卷积核在边界处的“质量”对不上于是AᵀA在边界附近不等于AAᵀ迭代会把误差往边界推。另一种可能是adegrade里裁剪的索引写错采到了错位一个像素的区域。解决把A和Aᵀ的imfilter边界选项写成同一个用rot90(psf,2)显式翻转核而不是依赖conv的自动行为;再做一次单点测试——把x设成中心一个冲激看Aᵀ(A(x))的峰值位置是否回到冲激原位偏移就说明索引错位。5.3 现象迭代几十次后数值爆炸图变成全白或全黑原因ρ给得过大导致x子问题条件数极差而inner内迭代次数不够CG 没解准就把误差带进了下一步;或者z子问题里的软阈值公式写成了1 - t/mag但mag没有加eps遇到零梯度时除零产生Inf一路传播。解决先加eps这是最廉价的一步;然后把inner提到 40、rho降到 0.5 重跑。如果还在涨把每次外迭代的max(abs(x(:)))打印出来看它是从第几步开始脱离[0,255]区间——突变点通常就对应某个子问题的公式错误。出图前务必min(max(x,0),255)这个“后悔药”能救回大部分可视化。5.4 现象PSNR 很高但放大看边缘还是糊的原因PSNR 对模糊非常宽容。一张轻微模糊的图和一张锐利的图只要整体均方误差接近PSNR 可能只差 0.3 dB。如果你的目标是视觉锐度而不是数值指标PSNR 会误导你。解决引入梯度域指标或者直接做边缘剖面检查——在图上取一条穿过强边缘的水平线把 HR、双三次、重建三条剖面画在一起看 10%→90% 上升宽度。TV 结果通常比 Tikhonov 窄 30% 左右如果两者剖面几乎重合说明 TV 的lambda给大了先验没起作用实际退化成了 Tikhonov。5.5 现象换一张图效果就崩同一组参数完全不通用原因参数是在单张图上扫出来的过拟合了那张图的纹理统计。纹理密集的图需要小lambda保细节平坦的图需要大lambda压噪声。另外噪声水平变了最优lambda也会跟着变二者大致成反比。解决建立一个小型验证集5~10 张覆盖不同内容类型的图参数在这上面取平均最优而不是单图最优。噪声水平未知时用平坦区域的标准差估一下 sigma然后按lambda ∝ 1/sigma做归一化。这套做法看着笨但比每次换图重新扫参快得多。6. 验收与推进用 PSNR/SSIM 卡线再用学习型先验抬上限最后说验收。我一般分两轨有真值的自测轨和无真值的真实轨。自测轨用第 3 章的方式造 HR/LR 对跑 PSNR 和 SSIM同时把双三次作为基线钉在表里——重建结果比基线高不到 0.8 dB就先别急着上更复杂的先验回头查模型。真实轨没有真值就靠边缘剖面、噪声平坦区标准差、以及“缩小回原尺寸后与原始 LR 的差异”这三个可量化的替代指标。参数扫描的记录表大概长这样可以直接套用lambdarho外迭代PSNR(dB)SSIM耗时(s)边缘宽度(px)0.010.84028.60.81212.42.10.030.84029.40.84613.11.60.100.84028.90.82113.52.80.031.54029.10.83818.71.7一个具体的进阶技巧用 TV 的解作为学习型方法的初始化或后处理。哪怕是简单的去噪自编码器或即插即用Plug-and-Play先验把 TV 结果当x₀喂进去收敛所需迭代数通常能砍掉一半以上而且不容易掉进局部极小。ADMM 框架天然支持这种替换——把z子问题的软阈值换成一个预训练的去噪网络前向推理就行x子问题一行都不用改。还有一个我改不掉的习惯每次跑长迭代之前先把x0、Y、psf、lambda、rho一起save(ckpt.mat)存一份外迭代每 10 次再覆盖存一次。TV-ADMM 调到第 50 次才发现参数理解错了、要从头再来那种感觉很不值。有了断点文件改完参数从第 40 次接着跑五分钟就能验证想法。希望帮到你。本文还有配套的精品资源点击获取