简介：面向图像去模糊与盲去卷积研究者的经典论文配套Matlab实现，突出超拉普拉斯先验在清晰图像恢复中的作用，适用于相机抖动、光学散焦等模糊场景的算法验证与二次开发。资源包共七个文件，包括四个m格式源码文件（覆盖主去卷积流程、图像求解、信噪比评估与测试入口）、一张真实场景示例照片、一份预置模糊核的mat数据文件，外加一个readme说明文档，整体仅约2.05MB，下载与部署成本很低。算法通过最小化带有超拉普拉斯先验的损失函数，交替估计模糊核与清晰图像，在恢复边缘细节和抑制噪声方面表现突出；源码结构清晰，可直接运行复现论文中的实验结果，也可作为改进实验的起始基线。对于需要构建非深度学习去模糊方法、或深入理解经典优化策略的研究者与工程师，这份代码同时提供了一个完整的示例流程，便于对比不同先验或核估计方案的性能。目前该资源已有1126人学习使用。
1. 一次运动模糊，让我把 Hyper-Laplacian Priors 从论文拉回工程
做图像复原时，我最早愿意把数学真正搬进工程的一次，就是 Fast Image Deconvolution using Hyper-Laplacian Priors 这个方案。它解决的问题很具体：一张因为手抖或对焦不准变糊的照片，在已知或近似已知模糊核时，怎么用尽量少的迭代次数把锋利边缘还原出来。和动不动要训练几小时、推断时还要考虑显存和 batch 的深度方法比，这类传统优化方法在几十次迭代内就能出可用结果，而且不依赖成对的清晰/模糊训练数据。
这套方法适合两类人：一类是在做摄影后期、监控图像增强、显微图像复原的工程师，手里经常有单张模糊图和一块事先估计好的核；另一类是刚接触去卷积的同学，想理解“先验”到底如何改变优化问题的形状。它不需要 GPU，CPU 上跑几百乘几百的图像也只需要几秒到十几秒。
下面我会按“退化模型为什么病态 → 超拉普拉斯先验为什么有效 → IRLS 怎么落地 → 最小可复现代码 → 踩坑记录 → 扩展思路”的顺序走一遍。目标是让你看完后能自己写一个可用的去卷积脚本，并且知道参数调到什么程度该停。
2. 为什么普通去卷积会振铃：退化模型与先验选择的由来
2.1 退化方程：卷积、加性噪声和病态的逆算子
图像去卷积的起点是一行很简单的式子：g = h ⊗ x + n。x是清晰图，h是点扩散函数（PSF），也就是模糊核，n是传感器或量化带来的噪声，g是观察到的模糊图。在频域里，卷积变成乘法，于是最直觉的做法是直接算X = G / H，把除法做完再反变换回来。
但这里有个致命的工程问题：运动模糊核的频谱H会在某些频率接近零，尤其是高频方向。相机抖动、失焦、大气扰动都会让高频信息在成像时被衰减到接近零。直接做除法，等于把那些频率上的微小噪声无限放大，结果就是整张图上出现一条条规则振铃。你甚至不需要加多少噪声，只要图像是 8bit 量化，逆滤波的结果就没法看。
所以去卷积本质上不是一个“反卷积”问题，而是一个“在噪声与模糊之间做权衡”的估计问题。我们要找的x不是让h⊗x无限接近g，而是在逼近g的同时，符合自然图像的基本规律。这就是先验进入的地方。
2.2 高斯先验救不了长尾梯度
最常见也最容易想到的正则项是 L2，也就是让λ||∇x||²尽量小。最小二乘加 L2 正则，在概率上等价于假设图像的梯度服从高斯分布。高斯分布的特点是中间很高、尾部衰减极快，这意味着它认为“大的梯度值”几乎不可能出现。
但你去统计任何一张真实照片，会发现自然图像的梯度直方图是一条尖峰加长尾的曲线：大部分像素确实接近零梯度，但边缘处依然存在相当大的梯度值，而且出现频率远高于高斯分布预测。L2 惩罚对大的梯度处罚太重，优化时会拼命把边缘压平来降低代价，最后得到一张“边缘被磨掉”的过平滑图。这也是很多基础去模糊教程里，维纳滤波做完总觉得图像蒙了一层雾的原因。
换句话说，高斯先验把“噪声”和“真实边缘”混为一谈。它看到一个大梯度，第一反应不是“这是细节”，而是“这是异常值”。于是细节连同噪声一起被抹掉。
2.3 超拉普拉斯先验：一个 α 参数同时控制稀疏度和形状
把惩罚从||∇x||²改成||∇x||^α，其中α小于 2，就得到 Hyper-Laplacian 先验。α=2时是高斯，α=1是全变分，α在 0.5 到 0.8 之间时，概率密度更接近自然图像梯度的长尾形状。
具体到目标函数，可以写成：
min_x ||h⊗x - g||² + λ (||gx||^α + ||gy||^α)
其中gx、gy分别是水平、垂直方向的一阶差分。为什么把水平和垂直分开而不是用sqrt(gx² + gy²)？因为真实图像里水平边缘和垂直边缘的统计特性并不完全一样，分开惩罚在恢复斜边和纹理时更稳。这一点在后面实现里会体现为两个独立权重数组。
需要说明的是，α<2会让目标函数变成非凸，全局最优不一定找得到。但在图像去卷积这个任务里，非凸带来的“边缘保持”收益远大于局部最优风险。真正决定速度的不是α的值，而是求解策略。
3. 让 Fast 名副其实：IRLS 和变量分裂怎么把非凸优化压成线性求解
3.1 用迭代重加权把非凸项变成局部二次型
|s|^α在α<2时没有全局统一的二阶近似，直接上牛顿法很容易在s=0附近出现除零和发散。IRLS（迭代重加权最小二乘）的思路是：每次迭代先把非凸惩罚固定成一组权重，再解一个加权最小二乘子问题。
具体来说，在当前位置x^k，用w_i = |s_i|^(α-2)作为权重，把|s_i|^α替换成w_i s_i²。由于α-2是负数，所以梯度接近零的像素权重大，梯度大的像素权重小。这个行为正好符合稀疏先验的直觉：我允许少数大梯度存在，但我希望绝大多数位置梯度尽量小。
每轮外循环更新权重，内层解一次线性系统：
(H^T H + λ Σ D_i^T W_i D_i) x = H^T g
其中D_i是梯度算子，W_i是由当前x计算出的对角权重。外循环 10 到 30 次，内层用共轭梯度法迭代 5 到 20 步，通常就能收敛到不错的结果。
3.2 权重更新的数值细节：eps、分方向和归一化
权重公式|s|^(α-2)在s=0时没有定义，实现时必须加一个小的eps做保护。常见做法是：
wx = max(|gx|, eps)^(α-2)
eps不能给太大，一般取图像动态范围的千分之一到万分之一。如果图像像素值在 0 到 1 之间，eps=1e-3比较稳；如果图像是 0 到 255，建议先除以 255，否则权重量级会差好几个数量级。
另一个工程细节是把水平和垂直方向的权重分开算。很多人图省事，算一个sqrt(gx²+gy²)然后两个方向共用同一个权重。这种做法会让对角方向的梯度被重复惩罚，恢复出来的斜边容易变成阶梯状。
最后，我会对权重做归一化：wx = wx * wx.size / (wx.sum() + 1e-12)。这样每个迭代的权重均值大致为 1，λ的数值含义在不同图像、不同模糊核之间基本可迁移。不归一化的话，同样一个λ=0.02，在一张暗图上可能完全不起作用，换到亮图上又可能过强。
3.3 内层求解：为什么不能全程 FFT 一把梭
如果权重是常数，正则项在频域是对角的，整个线性系统可以用两次 FFT 直接解。但 IRLS 的权重随像素位置变化，它在频域不是对角矩阵，没法直接做除法。这时候有两个选择：一个是上共轭梯度，每次迭代只做一次 FFT 和一次逆 FFT；另一个是论文里常用的变量分裂加半二次优化，把梯度项先用辅助变量替换，得到一个可以逐像素闭式求解的子问题和一个可以 FFT 闭式求解的子问题。
我在工程里更习惯先用 IRLS 加共轭梯度跑通逻辑，因为它的代码最短，也最容易检查每一项对不对。等确认目标函数没问题、参数语义清楚了，再替换成变量分裂方案加速。下面的最小实现就用了这个更直观的写法。
4. 用 Python 复现一个最小实现：从运动模糊到迭代去卷积
4.1 准备测试数据：合成模糊核与噪声
没有干净的数据集时，先用合成模糊自检。我习惯生成一张带矩形和孤立亮点的测试图，用水平运动模糊核做 FFT 卷积，再加一点高斯噪声。这样每一步都有标准答案，能直接算 PSNR。
4.2 主循环：IRLS 权重更新与共轭梯度求解
下面是一个可运行的最小实现。它依赖 NumPy 和 SciPy，输入灰度图范围建议在 0 到 1：
import numpy as np from numpy.fft import fft2, ifft2 from scipy.sparse.linalg import LinearOperator, cg def psf2otf(psf, shape): """把 PSF 转到 OTF，并保证中心在零频位置。""" h, w = psf.shape if h % 2 == 0 or w % 2 == 0: raise ValueError("PSF 尺寸需要是奇数") pad = np.zeros(shape) pad[:h, :w] = psf # 将 PSF 中心移到 (0,0)，负方向的部分自然绕到右侧 return fft2(np.roll(pad, (-(h // 2), -(w // 2)), axis=(0, 1))) def deconv_hyperlap(input_img, psf, lam=0.02, alpha=0.8, outer_iter=25, inner_iter=20): H = psf2otf(psf, input_img.shape) HtH = np.abs(H) ** 2 rhs = ifft2(np.conj(H) * fft2(input_img)).real # 维纳初始化：用一个小常数压住 H 接近零的频率 x = ifft2(np.conj(H) * fft2(input_img) / (HtH + 1e-6)).real for _ in range(outer_iter): # 当前图像的水平、垂直一阶差分（循环边界） gx = np.roll(x, -1, axis=1) - x gy = np.roll(x, -1, axis=0) - x # IRLS 权重，alpha-2 为负数，小梯度给大权重 eps_w = 1e-3 wx = np.maximum(np.abs(gx), eps_w) ** (alpha - 2) wy = np.maximum(np.abs(gy), eps_w) ** (alpha - 2) # 把权重均值归一到 1，让 lam 的含义跨图像稳定 wx *= wx.size / (wx.sum() + 1e-12) wy *= wy.size / (wy.sum() + 1e-12) def matvec(z_flat): z = z_flat.reshape(input_img.shape) # 保真项部分：H^T H z az = ifft2(HtH * fft2(z)).real # 正则项部分：D^T (w * D z) zx = np.roll(z, -1, axis=1) - z zy = np.roll(z, -1, axis=0) - z az += lam * (np.roll(wx * zx, 1, axis=1) - wx * zx) az += lam * (np.roll(wy * zy, 1, axis=0) - wy * zy) return az.ravel() n = input_img.size A_lin = LinearOperator((n, n), matvec=matvec, dtype=np.float64) x_hat, _ = cg(A_lin, rhs.ravel(), x0=x.ravel(), maxiter=inner_iter, tol=1e-3, atol=1e-6) x = x_hat.reshape(input_img.shape) return x下面做一个自检：
rng = np.random.default_rng(3) # 生成一块带边缘的测试图 x_true = np.zeros((96, 96)) x_true[24:72, 24:72] = 0.8 x_true[44:52, 44:52] = 0.1 x_true[10, 10] = 1.0 # 水平运动模糊核 psf = np.zeros((9, 9)) psf[4, 2:7] = 1.0 / 5 H = psf2otf(psf, x_true.shape) g = ifft2(H * fft2(x_true)).real g += rng.normal(0, 0.005, size=g.shape) rec = deconv_hyperlap(g, psf, lam=0.02, alpha=0.8) def psnr(a, b): mse = ((a - b) ** 2).mean() return 10 * np.log10(1.0 / mse) print("blurred PSNR:", psnr(g, x_true)) print("recovered PSNR:", psnr(rec, x_true))4.3 逻辑说明与参数含义
psf2otf是这套代码的地基。PSF 先填充到和原图一样大，再把中心搬到左上角原点。这么做的原因是 FFT 的卷积约定里，原点在图像的左上角，如果中心位置不对，复原结果会整体平移，看起来像“对不齐”。
IRLS 部分的关键在wx、wy两个权重。alpha=0.8时指数是-1.2，所以梯度接近零的像素会拿到很大的权重，梯度大的边缘却只受轻微惩罚。这就实现了“鼓励平坦、容忍边缘”的稀疏先验。lam=0.02是正则化强度，它控制你愿意为清晰边缘牺牲多少保真度。lam太大，结果会变得过于平滑；lam太小，权重保护不住噪声，结果会出现颗粒和振铃。
内层cg的maxiter=20并不追求完全收敛。IRLS 外层的权重本来就在变化，内层解太精确没有意义，反而浪费时间。比较合理的组合是外层 25 次，内层 10 到 20 次。如果发现结果在最后几次迭代还在明显变化，可以加外层次数；如果结果已经稳定但速度不够快，优先降内层次数。
4.4 代码边界说明
这段代码刻意用了循环差分和 FFT 卷积，边界条件隐藏在全图周期性里。它对 PSF 是全尺寸且居中核有效；如果你的 PSF 是裁剪过的局部核，比如只有左上角一块，请先自动补零并做循环移位。另外，代码没有做任何下采样加速，图像超过 1024×1024 后内层 FFT 次数会明显增加。这时建议先把图像缩小到 1/4 尺寸调好参数，再回全分辨率跑少量迭代。
5. Hyper-Laplacian 去卷积常见问题排查：五个让你返工的坑
5.1 λ 调了十倍，结果却纹丝不动
现象是无论把lam从 0.01 改成 0.1，输出图像几乎没有区别，或者某一刻突然从振铃变成一团糊。
原因是权重没有归一化，λ的真实作用量取决于梯度量级。图像像素范围是 0 到 255 还是 0 到 1，wx的绝对值会差十的三次方量级。你也可能在 0.02 时权重过小，0.2 时又瞬间过强，中间没有平滑变化。
解决方式是把输入图像归一化到 0 到 1，并像上面代码那样把wx、wy的均值压到 1 附近。之后lam就可以当作一个 0.01 到 0.1 之间线性可调的旋钮，而不是玄学参数。
5.2 图片四周出现一圈亮边或暗边
现象是中心区域恢复得不错，但四条边附近明显比中间亮，像加上了一个亮边框。
原因是梯度算子和卷积算子的边界条件不一致。FFT 卷积假设图像是周期延拓的，所以差分算子应该用循环差分；如果用了convolve2d配合mode='same'的零填充，边界处会出现一个人为的大梯度，正则项为了保护这个大梯度，会把边界单独“抬起来”或“压下去”。
解决方式有两种：一是全流程都用 FFT 卷积和np.roll差分，保持周期边界；二是在进入算法前，先把图像做镜像扩展 8 到 16 个像素，处理完再裁掉。第二种更适合真实照片，因为真实照片并不满足周期边界。
5.3 横平竖直的边缘不错，斜边却变成阶梯
现象是 45 度斜线恢复成一节一节的小方块，看起来像像素化加重的结果。
原因是把两个方向的梯度合并成了一个共同权重，相当于用sqrt(gx²+gy²)替代了|gx|^α + |gy|^α。对斜边来说，水平和垂直分量同时被惩罚，优化器宁可把能量分摊到两步阶梯，也不愿意留下一个大斜梯度。
解决方式是始终维护wx和wy两个独立权重数组。代码里已经这么做了，如果你参考的是早期版本实现，建议重点检查这里。
5.4 迭代到一半 PSNR 不升反降，外循环在振荡
现象是前三轮恢复效果快速变好，第五轮之后开始出现越来越多细碎噪点，PSNR 掉头向下。
原因是α太小，IRLS 权重更新过于激进。α越接近 0，正则项的非凸性越强，权重在“极小平坦区”和“极强边缘区”之间跳变，整个子问题容易在两个局部最优之间反复横跳。
解决方式是把α控制在 0.5 到 0.8 之间，并且不要在一开始就用最终参数。常见做法是先α=1跑十轮，再用α=0.8接着跑，相当于用全变分结果做热启动。这个技巧在论文和工程实现里都很常见。
5.5 已知核稍微不准，结果就出现梳子状条纹
现象是复原图里出现一排排平行细纹，方向通常和运动方向垂直，看起来像印刷网点。
原因是图像估计和核估计没有交替迭代。Hyper-Laplacian 只负责在给定核的情况下复原图像，但盲去卷积里核往往是从模糊图上估出来的，初始核误差会在高频区域被放大成周期条纹。
解决方式是不要一次把核用到底。每轮先固定核去卷积几张图，再用复原结果反过来更新核，并把核的支撑域约束在较小范围，最后再做一次核精修。这个多尺度交替方案我会在下一章展开。
6. 从已知核走向盲去卷积：多尺度验证与参数迁移
6.1 用合成模糊做最小验证：PSNR 不是唯一标准
我每调一次算法，第一步都是用已知核和合成模糊做回归测试，计算 PSNR。但 PSNR 只能告诉你整体误差，不能告诉你振铃在哪。实际看结果时，我会刻意找三个区域：高光边缘、平坦墙面、细纹理。高光边缘看有没有 overshoot 亮的白边，墙面看有没有把人眼看不到的噪声放出来，细纹理看是不是被磨平。
如果合成测试里这三个区域都合格，再拿真实模糊图去试。真实模糊没有 GT，判断标准变成了“文字边缘锐利但不带白边”和“肤色区域干净但没有塑料感”。这一步没有捷径，只能靠肉眼对比。
6.2 把 Hyper-Laplacian 接进盲去卷积的多尺度框架
盲去卷积的核心是交替更新核和图像。常见做法是先构建高斯金字塔，在最粗尺度估计一个模糊核，然后逐层上采样并精修。每一层的图像复原都用 Hyper-Laplacian 去卷积，但正则化参数λ要随着尺度缩放：下采样到 1/4 时，λ通常要乘以 1.5 到 2 倍，因为下采样本身已经平滑掉了部分噪声，保真项可以放得更松。
在核估计阶段，我会对核加一个||k||²或||k||¹的小正则，并每轮把核的非负约束和能量归一化重新加上。忽略这两个约束是梳状条纹的常见来源。
最后留一个我自己的习惯：把所有可调参数暴露成配置文件或命令行参数，不要写在函数里硬编码。因为α、λ、outer_iter在每一轮尺度上的最优值并不相同，交互式调参比每次改源码快得多。这个经验也是我在远程处理一张高噪声模糊图时得来的教训。希望帮到你。
本文还有配套的精品资源，点击获取
