简介：《基于超拉普拉斯先验的快速图像去卷积》论文的Matlab实现包，面向图像处理研究者与工程师，用于复现盲去卷积算法，解决运动模糊、光学模糊等造成的图像退化问题。压缩包共七个文件，包含四个Matlab脚本（主去卷积程序、图像求解函数、测试脚本及信噪比函数）、一张测试图片、一份模糊核数据和一份说明文档，整体大小仅约两兆字节。目前已有1126人浏览学习，可直接运行测试脚本对比恢复前后的图像差异。该实现以超拉普拉斯先验刻画图像梯度的稀疏性，在去噪与边缘保持之间取得平衡；配合预置的模糊核，可验证不同参数下的复原效果，也为后续研究（如深度学习方法）提供了传统基线。通过阅读源码和说明文档，能理解迭代优化与先验建模的关键细节，便于迁移到自己的项目中。

1. Hyper-Laplacian Priors 凭什么成为图像去卷积的默认选项？

深夜处理一段行车记录仪视频，路牌文字糊成一条光带。试过 Wiener 滤波，边缘是绷得住了，但字幕四周一圈一圈的振铃简直像鬼影；换成普通的 L2 正则去卷积，噪声倒是压住了，可文字边缘软得像水彩。最后让我收工的，是标题里这个经典思路——Fast Image Deconvolution using Hyper-Laplacian Priors。它做的事一句话能讲清：自然图像里边缘梯度不是高斯分布，而是重尾的，与其用 α=2 的二次惩罚强行压平滑，不如用一个 α 在 0.5～0.8 的 Hyper-Laplacian 分布去拟合真实梯度，再用查表或加权迭代把求解速度拉回实时区间。这篇笔记想把这条路线拆开：目标函数怎么立、迭代怎么解、参数怎么调、真实落地时在哪几步最容易翻车。

2. 把图像去卷积拆成目标函数：用 Hyper-Laplacian Prior 改写传统 Wiener 滤波

2.1 退化模型与“反卷积不适定”到底在哪

所有去卷积问题都从一个假设出发：观测到的模糊图像 y 是清晰图像 x 与模糊核 k 卷积后再叠加噪声的结果。写成连续形式是 y = k ⊗ x + n。这里 k 是点扩散函数（PSF），可能是失焦产生的圆盘、运动产生的拖尾、或者是光学系统里没法用几何光学描述的那一部分响应。n 是传感器噪声，通常按高斯加性噪声处理，但真实情况下它也包含椒盐噪声和读噪声的混合。

问题在于：给定 y 和 k，直接求 x 是一个不适定问题。卷积在频域里是乘法，理论上只要把 Y(ω)/K(ω) 做个傅里叶逆变换就能得到 X。可 K(ω) 在高频处会趋近于零，于是这一除法把噪声里最微弱的高频分量放大成了满屏的颗粒。这就是为什么教科书里一定会提到“直接逆滤波是不行的”——不是数学上无解，而是解对噪声极其敏感。

正则化就是在目标函数后面加一项惩罚，把解拉向“看起来像自然图像”的集合。传统 Wiener 滤波在频域里把噪声功率和信号功率的比值作为阻尼因子，形式上很优雅，但它假定信号和噪声都是高斯分布。真实照片的高频分量不是高斯噪声，而是由边缘构成的稀疏强响应，用高斯先验去约束它，结果就是边缘被当成需要抹平的东西。这也是 Wiener 去卷积出来的图像总有一种“塑料感”的根源。

2.2 从 MAP 目标函数到“保真项 + 梯度惩罚项”的完整推导

把去卷积看成一个最大后验估计（MAP）问题：需要找的清晰图像 x，使得给定观测 y、模糊核 k 和先验分布 p(x) 条件下，后验概率 p(x | y, k) 最大。取负对数之后，优化目标变成：

J(x) = ||Kx − y||² + λ · Σᵢ ρ(Dᵢx)

第一项叫保真项，衡量重构图像经过同一个模糊核后再投影回来与观测的差异；第二项是正则项，其中 Dᵢx 是图像在某个方向的差分（水平梯度、垂直梯度，有时也加上二阶差分），ρ 是惩罚函数。λ 是两者的平衡权重。K 是模糊核构成的卷积矩阵，实际写成矩阵形式只是为了推导方便，实现时都是直接做卷积或 FFT。

传统选择 ρ(t) = t²，也就是 L2 范数，对应高斯先验。它的优点是目标函数变成二次型，求导后是线性方程组，一次迭代就能收敛；缺点前面说了，对边缘过度平滑。后来有人换 ρ(t) = |t|，即 L1 范数，对应拉普拉斯先验。L1 在零点不可导，迭代要用次梯度或近端算子，但它对大幅度的梯度惩罚相对温和，能更大概率保留边缘。我对这两种选择都跑过实验，L1 比 L2 好一截，但中等强度的纹理区域会出现一种“磨砂”效果，像锐化滤镜拉过头的打印品。问题出在自然图像梯度的经验分布既不是高斯也不是标准拉普拉斯，而是介于两者之间、尾部更重的形态。

2.3 Hyper-Laplacian 的 α 次范数：为什么 0.5~0.8 才是常见工作的区间

Hyper-Laplacian Prior 的基本思想是：先从大量自然图像里统计梯度幅值的概率分布，得到一个重尾的经验曲线；这个曲线介于高斯和拉普拉斯之间，用一个更灵活的形式 ρ(t) = |t|^α 来拟合，α 可以从 0.3 取到 1.2。实验观察是：自然图像的梯度幅值分布在对数坐标下近似一条直线，斜率对应的 α 一般在 0.5 到 0.8 之间。α 取 2 时退化为 L2，α 取 1 时退化为 L1，而真正的拟合点落在二者中间。

为什么这个看似微小的改动影响这么大？当 α < 1 时，惩罚函数在零点附近的导数是无穷大，这会让优化过程强烈抑制微小的梯度波动——噪声就是一堆小幅度的随机梯度，正好被压住；而在梯度幅度较大的边缘区域，惩罚的增长是次线性的，也就是说边缘不会被“砍掉”。这正好同时处理了振铃和过度平滑两个矛盾。另一个数学上的麻烦是：α < 1 时整个目标函数是非凸的，严格求解会落入局部极小。论文里给出的应对是查表法加变量分裂，实际工程里也更常配合多尺度金字塔去逼近全局最优。

这也是我特别想提醒的一点：Hyper-Laplacian 不是“把 α 设小就万事大吉”，它把去卷积从凸优化问题换成了一个非凸问题。α=0.8 时凸性相对来说还能接受，α=0.4 时你得额外接受结果对初值敏感这个代价。我一般在小纹理区域比较多的图像上用 0.7～0.8，纯文字或图纸渲染图可以用 0.5，后者边缘极稀疏，更容易收敛。

3. 用 IRLS 做快速求解：最小可复现的 Hyper-Laplacian 图像去卷积流程

3.1 IRLS 的基本套路：把非二次的 α 次范数二次化

直接优化 J(x) 时，α<1 使梯度惩罚项不可导。一个稳定的求解方案是 Iteratively Reweighted Least Squares（IRLS），核心思想是：把 |t|^α 写成 |t|^(α−2) · t²，然后每次迭代用当前的 t 值算出权重 w = α · max(|t|, ε)^(α−2)，把非二次项“冻结”成一个已知的二次项。这里 ε 是一个极小的正数，防止梯度为零时除零。这样每一轮迭代都只需要解一个线性系统——问题从“非线性优化”降格成“反复解线性方程”。

IRLS 的收敛性和 α 大小有关。α 越接近 1，收敛越稳；α 越小，权重变化越剧烈，需要更多外迭代。实际跑的时候，我看的不只是目标函数下降曲线，还会监控梯度惩罚项的变化量，如果连续两轮的变化小于 1%，就停止。额外加了一条保险：如果迭代过程中 PSNR 连续两轮下降超过 0.2dB，就回退到上一轮的 x，再调小 α 或者减小步长重启。这个回退机制在α<0.5时几乎每次都会触发，后来变成一个固定习惯。

3.2 组装保真项与正则项的线性系统

把 IRLS 用到具体目标函数上，每一轮需要解：

(KᵀK + λ · Dᵀ W D) x = Kᵀ y

其中 D 是由水平差分 D_x 和垂直差分 D_y 堆叠成的梯度算子矩阵，W 是每轮更新的对角权值矩阵，每个对角线元素对应一个像素在某个梯度方向上的权重。这个系统有两个实现难点：一是 KᵀK 和 DᵀWD 的维度是图像像素数的平方，显式构造矩阵不现实；二是正规方程有条件数问题，需要做共轭梯度或预处理。

工程上我的做法是：显式构造稀疏梯度算子（换行列向量很简单），共轭梯度迭代不显式拿矩阵，而是在每一步计算矩阵向量乘积。保真项部分用卷积算子。这样做的好处是避免使用大矩阵求逆，共轭梯度只跑 20～50 步即可，因为每一轮 IRLS 都只是对上个解的微调。边界处理与循环托普利兹近似在第五章专门讲，这里先按内部图像区域来。

3.3 落地到灰度图像：可直接改的最小 Python 实现

下面这段代码是我通常用来验证先验效果的骨架，它不追求极致性能，但结构清楚，能直接跑通一张 64×64 的灰度图。清晰图像 x0 先用一个已知模糊核 k 卷积加高斯噪声得到观测 y，再用 IRLS 还原。代码里的核心逻辑是外层 IRLS 更新权重、内层共轭梯度解方程。

import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla from scipy.signal import convolve2d from scipy.sparse import spdiags def gradient_matrices(h, w): """构造水平与垂直一阶差分算子，行数 = 像素数""" n = h * w idx_r = np.arange(n) idx_cr = np.arange(n) val_r = np.ones(n) rows, cols, vals = [], [], [] # 水平差分：x(r,c) - x(r,c+1)，列方向右边界不约束 for r in range(h): for c in range(w - 1): i = r * w + c rows += [i, i] cols += [i, i + 1] vals += [1.0, -1.0] Dx = sp.coo_matrix((vals, (rows, cols)), shape=(h*w, h*w)).tocsr() # 垂直差分：x(r,c) - x(r+1,c) rows, cols, vals = [], [], [] for r in range(h - 1): for c in range(w): i = r * w + c rows += [i, i] cols += [i, i + w] vals += [1.0, -1.0] Dy = sp.coo_matrix((vals, (rows, cols)), shape=(h*w, h*w)).tocsr() return Dx, Dy def solve_one_step(Kx, Dx, Dy, y_flat, lam, W_flat): """给定当前梯度计算权重，解一次加权最小二乘""" n = y_flat.shape[0] # 共轭梯度法求解 (KtK + λ Dt W D) x = Kt y # 这里 Kx 表示模糊核卷积，卷积转置做到算子层面 def matvec(x): x_mat = x.reshape(y_shape) kx = convolve2d(x_mat, kernel, mode="same", boundary="fill") ktx = convolve2d(kx, kernel_rev, mode="same", boundary="fill") dx = Dx @ x dy = Dy @ x reg = lam * (Dx.T @ (W_flat[:, 0] * dx) + Dy.T @ (W_flat[:, 1] * dy)) return ktx.ravel() + reg A = spla.LinearOperator((n, n), matvec=matvec) y_mat = y_flat.reshape(y_shape) b = convolve2d(y_mat, kernel_rev, mode="same", boundary="fill").ravel() x_est, info = spla.cg(A, b, x0=y_flat, rtol=1e-4, maxiter=40) return x_est.reshape(y_shape)

逻辑说明：外层 IRLS 在每次迭代时计算每个像素两个方向的梯度幅值 gx = |Dx·x|、gy = |Dy·x|，然后按 Hyper-Laplacian 的导数公式更新权重 w = α * max(g, ε)^(α−2)。权重矩阵 W 是对角的，只存两个长度等于像素数的向量即可。代码里 matvec 函数返回 Kᵀ(Kx) + λ Dᵀ(WDx)，把卷积算子当成线性系统的主对角块，用共轭梯度解方程。

参数说明：kernel 是归一化后的 PSF，kernel_rev 是其 180° 旋转后的核，用于计算 Kᵀ 的卷积。lam 是正则项系数，数值上对应噪声逆方差；α 典型取 0.6~0.8；ε 作为防除零保护，常见做法是取梯度中位数的 1%～5%。共轭梯度迭代次数不要给太大，我一般设 30~50 次；外迭代 10 轮内基本稳定，再往下跑收益很小。这段代码的性能瓶颈在 Python 层的循环构图，换到生产环境时我会把差分算子做成双线性滤波形式的矩阵乘法，或者直接用傅里叶域循环矩阵，速度能提升两个数量级。

4. 参数怎么调才不翻车：α、λ、迭代数与金字塔策略

4.1 α 与 λ 的常用区间和调节方向

Hyper-Laplacian 里有三个参数直接影响输出：α、λ、以及 IRLS 的外迭代次数。很多人第一次上手会把 α 当重点，其实 λ 才是决定成败的那一个。我用表格列出常用区间和调节方向，这些数值来自我在一批公开测试集上的回归经验，不是拍脑袋。

λ 的调法我向来有固定的顺序：先用 1e-2 起步，若输出出现噪声颗粒就加大，若振铃明显就减小。如果这两个症状同时出现，说明不是 λ 不对，而是模糊核估计不准，或者噪声模型和图像内容不匹配。λ 跟图像灰度范围强相关，如果你把输入从 0～255 缩放成 0～1，λ 数值要除以 255²，很多新人在这里翻车。

α 的调节逻辑相对独立：固定 λ，α 越接近 0.5，边缘越“脆”，但同时中低纹理区域会逐渐出现不自然的阶梯状色块；α 越接近 1，图像整体越平滑，接近 L1 的效果。一个保险做法是对同一张图跑 α=0.6 和 0.8 两个版本，对比边缘锐度，而不是盲信单一数值。

4.2 多尺度金字塔：模糊核尺寸超过 15 像素时的必然选择

当模糊核尺寸大于图像短边的 3% 以上，直接对全分辨率图做 IRLS，结果大概率会掉进局部极小。原因是目标函数非凸，误差函数带有大量局部极小盆，直接优化容易收敛到一个振铃严重的局部解。解决方法是构造图像金字塔：从低分辨率开始估计大尺度结构，然后依次上采样作为下一层的初值。

具体尺寸参数我一般这样定：金字塔层数等于 log2(核大小/5) 取整数，最高层尺寸保证核至少占 20～30 个像素。每层做固定轮数的 IRLS，得到结果后上采样到下一层。核本身在每一层按比例缩放，同一金字塔里核长度保持一致的比例。这个小改动让盲去卷积的成功率提升非常明显，我险些以为是自己代码写错了——实际上只是初始解离全局最优太远。

4.3 收敛判定：看相对残差和正则项变化量，不要只看 PSNR

我见过不少工程实现以“外迭代 20 次”作为固定终止条件，这其实浪费计算时间。更好的做法是同时监控三个量：观测误差 ||Kx−y||²、正则项 Σρ(Dx)、以及相邻两轮 x 的相对变化。前两者是目标函数的两个部分，后者反映解是否还在移动。

实现时常见的做法是为了避免计算完整目标函数。正则项可以在 IRLS 更新权重时顺手算出来，代码几乎免费。收敛标准我写这颗：相对变化小于 1% 或正则项变化小于 0.5% 即停止。同时检查目标函数值是否持续下降，若连续两轮上升，说明权重更新过猛，应退回上一轮解并减小 update step（可以把权重更新公式改成加权平均：w_new = 0.7·w_new + 0.3·w_old，能消掉大部分振荡）。

在验证阶段用的指标是 PSNR，但 PSNR 在某些情况下不准确——比如振铃集中在边缘附近时，像素误差并不大，视觉观感却很差。所以在调参时我从来都输出一张对比图：原图、观测图、去卷积结果并排放在一张画布里，用眼睛看比用 PSNR 快得多。后来加上结构相似性（SSIM）作为辅助指标，才把客观评价从黑匣子状态里拽出来。

5. 图像去卷积高频踩坑：背光亮条、颗粒感、彩边与迭代不收敛

5.1 现象：输出出现规则的“洗衣板”振铃纹

这是去卷积里最常见也最容易让人怀疑算法本身的失败模式，现象是锐利边缘两侧出现等间距的一圈一圈高亮纹，间距正好等于模糊核尺寸的一半左右。原因是边界处理。多数实现默认用零填充做卷积，高频分量的能量在图像边界处被突变截断，产生虚拟的高频响应；反卷积时这个响应被当成真实信号放大，形成沿边界传播的振铃。

解决方法是换边界模型。我一般先试镜像填充（reflect padding），再配合边界权重：给靠近图像边缘 1 倍核半径的区域分配更低的保真项权重，因为那里本来就没有可靠信息。另一个治标的方法是先把图像四周裁剪掉核半径大小的区域，去卷积完成后再拼接回来，但这样会损失画幅。

5.2 现象：噪声被放大成白色颗粒

现象是去卷积后原本看起来干净的背景区域出现密集的白色小颗粒，像老照片加了一层噪点纹理。原因是保真项权重太大（λ 太小），或者观测噪声方差估计低估了。IRLS 的权重更新天然会放大高频，如果保真项在频域没有足够阻尼，高频噪声就会被成倍放大。

解决的通常做法有两种。其一是前置去噪：在进去卷积之前用双边滤波或 BM3D 先压一层噪声，但注意不要过度去噪导致细节丢失。其二是后置处理：把 λ 调大，然后接受边缘锐度略降；或者把 IRLS 的权重下界设置高一些，防止某个梯度极小的平坦区域被赋予极大的权值。注意平坦区域像素梯度的绝对值很小，max(g, ε)^(α−2) 会很大，相当于让正则项在整个平坦区域内施加过强约束，这个效应会把噪声压成一种“瓷砖状”的马赛克。

5.3 现象：彩色图像边缘出现青/品红晕

当直接对 RGB 三通道独立做去卷积，输出里会出现彩色边缘——物体的亮侧偏青，暗侧偏红。原因是自然图像在 RGB 通道之间的高频分量具有强相关性，独立惩罚梯度会把这些相关性破坏掉，等效于在通道间引入无中生有的色差。

解决方法是把图像转到 YCbCr 或 Lab 空间，对亮度通道做完整去卷积，对色度通道只做轻量平滑或者保真项直接不做梯度惩罚。这样既保留彩色信息，又避开了跨通道的伪彩色。要注意 YCbCr 转换时如果输入是 sRGB，需要先做反伽马变换再转，否则亮度通道的梯度统计和自然图像先验不匹配。

5.4 现象：loss 曲线一直抖动，迭代不收敛

现象是目标函数在前几轮下降后开始来回振荡，或者干脆上一轮 PSNR 高、下一轮低，像在解空间里来回弹跳。原因是 IRLS 权重更新时，某些像素梯度为 0，导致 max(g, ε)^(α−2) 直接冲到无穷大。另一个常见原因是外迭代权重更新太激进，权重值变化率超过了解自身的收敛速度。

解决要分两步：首先要给权重加显式截断，我一般把 w 限制在 [1e-4, 1e3] 之间，防止极端值主导方程；其次是把权重更新改成带惯性的形式，w_new = 0.5·w_new + 0.5·w_old，这个系数就是阻尼因子，能有效抑制振荡。另外内层共轭梯度的收敛阈值不能放松太多，否则每一轮的外迭代都从有误差的解起步，误差不断累积。把 rtol 固定在 1e-4 以下，并且把上一轮的共轭梯度解作为下一轮的初值，这点能省一半迭代。

5.5 现象：同一组参数在另一张图上完全失效

现象是同一份代码、同一套参数，换一张图之后输出质量骤然下降，要么振铃，要么噪声爆炸。原因往往不是参数问题，而是模糊核和图像内容的比例关系变了。模糊核是固定大小，但画面的聚焦主体、纹理密度、噪声水平在不同照片之间差异很大，正则强度的“最优值”自然也变。

解决办法是做归一化预处理：把图像和核都归一化到相同的能量水平，噪声方差先估计出来，再按噪声方差去缩放 λ。这样能把参数调整范围压缩到与图像内容无关的相对尺度上。我自己有个固定习惯：处理任何图之前先打印三行信息：图像尺寸、模糊核尺寸、估计噪声方差。这三行信息决定了我把 λ 的起点设在哪。跑多了你会发现大部分“参数玄学”其实只是没做这一步跨图归一化。

6. 再走一步：把 Hyper-Laplacian 先验接进盲去卷积与真实噪声场景

6.1 从非盲到盲：交替估计模糊核与潜像

非盲去卷积假设模糊核已知，但真实照片里我们通常是不知道模糊核的。常见做法是把 Hyper-Laplacian 去卷积作为盲去卷积的内层模块：外层用迭代法交替估计模糊核 k 和潜像 x。初始化一个简单核（比如高斯盘或直线），然后跑几轮 Hyper-Laplacian 去卷积得到 x，再用这个 x 去估计新的 k，再回来跑一次，交替到模糊核变化不大。注意盲去卷积里初始核的形状会影响结果收敛方向，我一般用对称的圆形核起步，比矩形核更容易稳定。

6.2 真实拍摄照片：先估计噪声电平，再决定 λ 数量级

真实抓拍和合成图像的去卷积差别比想象中大。合成模糊是在干净图像上加已知卷积和噪声，而真实照片里的模糊核可能包含光学系统的畸变、传感器 CFA 插值、运动模糊的非线性。拿到真实照片时，我第一步会在图像的纯色区域（比如墙面、天空）采样一些像素块，用它们的方差近似估计噪声电平，再按 λ 与噪声方差成正比这一关系定初始值。另一个常用检查项是把模糊核画出来，如果它带有明显的不对称或旋转特征，说明可能是运动模糊而非失焦，那种情况下 Hyper-Laplacian 先验一样适用，但多尺度金字塔尤其重要。

6.3 验证方法：合成模糊 + PSNR/SSIM 回归，别用一张图下结论

每次调完参数，我的验证流程是固定的：拿 8~12 张不同内容的清晰图，每张做 3 种不同模糊核（高斯、直线运动、失焦圆盘）和 3 种噪声等级，合成出模糊观测，再去卷积，最后统计 PSNR 和 SSIM 的平均值而不是单张最佳结果。这组回归不仅筛掉了盲目调参，也能暴露某些参数组合“只对特定图有效”的隐蔽问题。把测试图固定在一张漂亮的渲染图上做调参，几乎一定会得到一张好看但对其他图完全失效的参数组合——这个坑我踩过不止一次。

现在接手新的去卷积任务时，我已经不再问“用什么先验”，而是先问“噪声多大、核多大、图像里有没有大面积平坦区域”。这三个答案定下来，Hyper-Laplacian 的 α、λ、迭代策略基本就有一个合适的出手点了。如果你也刚开始跑这条路线，我建议先停掉对参数的精调，花半天时间把退化模型和边界处理做成日志打印出来，每一次失败都记录了当时的三项输入，你会比我更早找到稳定复现的手感。希望帮到你。

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