盲去卷积图像复原：从交替最小化到Python实践

简介：这是一份面向图像处理、光学成像、天文学及医学影像等领域研究者的盲去卷积MATLAB算法实现，主要解决因大气湍流、镜头缺陷、像素响应不均等未知模糊核导致的图像退化问题，可在恢复清晰图像的同时估计点扩散函数（PSF）。压缩包共含1个m文件，包体约2KB，代码精简，聚焦盲去卷积核心流程，并可能结合Tikhonov或全变分正则化策略来提升稳定性。目前已有542人学习下载，适合具备一定MATLAB和图像复原基础的中高级学习者深入研读。文件内通常涵盖图像读取与预处理、迭代去卷积、PSF初值与更新、残差或相似度准则判断等模块，能够帮助读者快速掌握Richardson-Lucy等迭代类算法的工程实现，也可作为科研项目中的基础验证脚本或教学演示代码。

1. 盲去卷积复原：不知道点扩散函数，怎么把模糊图像救回来

手里只有一张严重失焦的老照片，想复原出文字和边缘，却连镜头模糊的半径都不清楚，这是典型的盲去卷积复原任务。图像退化可以看成清晰图与点扩散函数（Point Spread Function，PSF）的卷积再加噪声，普通去卷积假设 PSF 已知，而盲去卷积要把 PSF 和被复原图同时从一张观测图里估计出来。问题从一条方程变成两条未知量，解空间急剧膨胀，局部极小、振铃、核分裂都会冒出来。做图像算法、显微与天文数据处理的工程或科研人员，每天都在和这类问题打交道。后面按“建模 → 交替迭代 → 调参 → 排错”的顺序，把 Blind Deconvolution Algorithm 的常见解法一步步落地。

2. 盲去卷积的数学模型与交替最小化：未知PSF和清晰图一起估计

2.1 退化模型与目标函数：为什么盲去卷积是高度病态的

图像退化可以写成离散卷积形式：y[m,n] = ∑ h[i,j]·x[m-i,n-j] + n[m,n]，简写为 y = h ⊗ x + n。这里 h 是 PSF，x 是潜在清晰图，n 是加性噪声。盲去卷积要求同时估计 h 和 x，但最小二乘目标 min ‖y − h⊗x‖² 在数学上非常别扭。

第一个问题是尺度模糊：若 (x, h) 是一组解，则 (cx, h/c) 也产生同样的 y，相当于图像变亮、核变暗，卷积结果完全不变。第二个问题更麻烦：卷积对两个输入是对称的，h⊗x 和 x⊗h 的结果一样，所以把一个窄核配上一张严重模糊的图，也能解释观测结果。如果我们不做任何约束，退化过程可有无数个解释，普通矩阵求逆直接算只能得到一张纯噪声图。

盲去卷积之所以还“可解”，完全靠先验撑住：清晰自然图像的梯度是重尾分布，而 PSF 非负、有界、支持域有限，两个先验组合起来才能从乱解里挑出一个像样的组合。下表把非盲去卷积和盲去卷积放在一起看差异。

2.2 交替最小化框架：固定PSF求图像，固定图像求PSF

直接对联合目标做梯度下降几乎不收敛，常见做法是交替最小化（Alternating Minimization）。外层迭代每轮做两次子问题求解：

x⁽ᵏ⁺¹⁾ = argminₓ ‖y − h⁽ᵏ⁾⊗x‖² + λ_tv·TV(x)

h⁽ᵏ⁺¹⁾ = argminₕ ‖y − h⊗x⁽ᵏ⁺¹⁾‖² + λ_h‖h‖²，约束 h≥0，∑h=1

固定 h 更新 x，就是一个普通去卷积问题，可以用 Richardson-Lucy 迭代，也可以用梯度下降；固定 x 更新 h，则是一个带约束的核估计问题。两轮交替反复推进，这是 Majority 的 Blind Deconvolution Algorithm 实现的骨架。Richardson-Lucy 的盲版本本质上也落在该框架里，只是 x 更新时塞进了泊松噪声假设下的乘法迭代。

这个框架的实现成本低，因为两块都有成熟算法可复用；缺点是贪心式交替容易掉进局部极小，最终解对初始化相当敏感。工程上为了稳，一般把 h 初始化为一个宽高斯核，x 直接取 y，再配合多尺度策略把局部极小问题减小，后者在第 4 章展开。交替轮数不能太少，通常外层 20~30 轮、内层 3~5 轮起步。

2.3 图像和PSF的正则化与约束项

自然图像的正则项首选全变分（TV），TV(x) = ∑√(∇ₓ² + ∇ᵧ²)，它允许强边缘存在，同时抑制平坦区噪声。核的正则项一般取二范数或平滑约束，因为 PSF 通常不是高频振荡的形状，一个轻微的二范数惩罚能让核估计更干净。

三个必须保留的硬约束：h≥0，∑h=1，x≥0。h 非负保证 PSF 物理可解释；和归一化消掉尺度模糊的残余自由度；x 非负对多数光学成像成立，防止反卷积把图推出负亮度。约束的加入方式很简单，每次更新后做 clip 和归一化投影即可。参数方面，如果图像先归一化到 0~1，λ_tv 经验范围在 1e-3 到 1e-2，λ_h 在 1e-4 到 1e-3。λ_tv 太大会把图磨成卡通画，太小则振铃和噪声一起放大，第 4 章会给出调参顺序。

3. 用Python复现Blind Deconvolution Algorithm：交替迭代代码与参数整定

3.1 构造退化样本与PSF初始化

先用合成数据跑通流程，才能判断算法行为。下面的代码生成高斯核和模糊观测图，全程用 NumPy 与 SciPy，没有额外依赖。

import numpy as np from scipy.signal import fftconvolve def make_psf(size, sigma): ax = np.linspace(-(size // 2), size // 2, size) xx, yy = np.meshgrid(ax, ax) psf = np.exp(-(xx**2 + yy**2) / (2 * sigma**2)) return psf / psf.sum() def degrade(img, psf, noise=0.005): y = fftconvolve(img, psf, mode='same') y += np.random.normal(0, noise, y.shape) return np.clip(y, 0, 1)

make_psf 生成主瓣宽度由 sigma 控制的高斯核，并做归一化；degrade 用 fftconvolve 完成快速卷积，再叠加高斯白噪。输入图像请预先归一化到 0~1 的灰度图，否则后续正则参数要跟着缩放。初始 PSF 我一般用 sigma=2、尺寸 11 的高斯核，x 直接用观测图 y，不用做任何预处理就能跑交替迭代。

3.2 交替迭代核心函数 blind_deconvolution

def tv_grad(u): gx = np.gradient(u, axis=0) gy = np.gradient(u, axis=1) gnorm = np.sqrt(gx**2 + gy**2 + 1e-8) return -np.gradient(gx / gnorm, axis=0) - np.gradient(gy / gnorm, axis=1) def blind_deconvolution(y, h0, outer=30, inner=5, lr_x=0.25, lr_h=0.05, lam_tv=0.008, lam_h=0.0005): x = y.copy() h = h0.copy() h /= h.sum() for _ in range(outer): # 固定 h，更新 x：数据保真项梯度 + TV 正则梯度 for _ in range(inner): res = fftconvolve(x, h, mode='same') - y grad_x = fftconvolve(res, h[::-1, ::-1], mode='same') x = x - lr_x * (grad_x + lam_tv * tv_grad(x)) x = np.clip(x, 0, None) # 固定 x，更新 h：对 h 求梯度，再做非负与归一化投影 for _ in range(inner): res = fftconvolve(x, h, mode='same') - y grad_h = fftconvolve(res, x[::-1, ::-1], mode='same') h = h - lr_h * (grad_h + lam_h * h) h = np.clip(h, 0, None) h /= h.sum() return x, h

grad_x 的计算用的是翻转卷积核再做卷积，等效于卷积算子的转置，可以把它理解为把残差“反传”回图像空间；tv_grad 是 TV 项的负梯度近似，用散度算子实现。h 更新时 grad_h 同样用 x 的翻转卷积完成，加 lam_h*h 相当于给核加了一个轻微的二范数惩罚，让核更容易收敛成平滑山包而不是噪声碎片。内层迭代设 5 次，是为了每次外层交替前让变量充分松弛，内层太小会抖，太大等于反复解子问题，耗时成倍增加却不一定更好。x 和 h 每次更新后都做非负投影，h 再归一化到和为 1，这对应 2.3 节的两个硬约束。

3.3 关键参数与收敛诊断

跑完一轮可以在循环里记录当前数据残差和图像 PSNR。如果观测图 y 有真实清晰图 x_true，可以打印 psnr(x, x_true)：正常情况会先上升，之后可能缓慢下降，这往往是振铃开始占上风，选 PSNR 峰值对应的轮次即可，不必跑满迭代。如果 h 出现了明显的多峰分裂，说明 lam_h 或 lr_h 不合适，优先降低 lr_h。模型完全失衡时，x 上会出现高频棋盘纹理，此时先检查约束投影是否生效，再查 h 是否被 clip 成全部为 0。

4. 盲去卷积复原实战：多尺度策略、评价指标与调参顺序

4.1 coarse-to-fine多尺度：绕过局部极小的常用策略

单尺度交替最小化对核初始尺寸非常敏感，核尺寸给大几个像素，h 就容易把能量摊成一个平盘。实际工程中很少只跑一层，而是用金字塔从粗到细推进：先把观测图降采样两层，在最粗尺度上估计一个粗核和粗图像，然后插值放大作为下一层初始值。粗尺度上核相对图更小，优化曲面更平滑，局部极小更少。

from scipy.ndimage import zoom def multiscale_blind(y, psf_size=15, levels=3, **opts): pyr = [y] for _ in range(levels - 1): pyr.append(zoom(pyr[-1], 0.5)) h = make_psf(max(3, psf_size // (2 ** (levels - 1)) + 1), 1.0) x = pyr[-1] for lev in range(levels - 1, -1, -1): if lev < levels - 1: x = zoom(x, 2.0) h = zoom(h, 2.0, order=1) h /= h.sum() x, h = blind_deconvolution(pyr[lev], h, **opts) return x, h

核心逻辑是先把观测图缩成金字塔，在最底层用一个较小的高斯核启动，每向上一层就把上一轮的 x 和 h 双线性放大。zoom 默认阶数为 3，核上采样时装成 order=1 更稳，避免核插值产生负值。粗层迭代不必像细层那样严格，内层轮数可以少一半，因为它的任务只是提供一个靠谱的初始支撑域。多尺度跑完后 h 的尺度是相对于最终原图的，无需再缩放。

4.2 复原质量怎么评：PSNR、SSIM与核误差

合成实验里 ground truth 齐全，评价要同时看图像质量和核精度两个维度。

核误差老被低估：FFT 卷积的周期平移会让估计核整体位移几个像素，数值上看误差很大，但对复原图影响很小。所以比较前用 ndimage.shift 把 h_est 的质心挪到 h_gt 质心位置再算。现实任务没有真实图时，可以看数据残差 ‖y − h⊗x‖/‖y‖，残差平稳下降说明代数上收敛了，再配合可视化看边缘是否振铃。

4.3 调参顺序与三个误用

我一般按“先核后图”的顺序调参：第一步固定 psf_size，用测试结果判断估计核是否光滑、支持域是否合理；第二步调 lam_tv，看边缘过冲被压住且纹理不过度磨皮；最后才动 lr_h 和 inner。把 lam_tv 从 1e-3 逐步加到 5e-2，每次只翻一倍，观察 PSNR 和振铃强度。噪声大的观测图，先做一次轻微的高斯预滤波，否则算法会把噪点当作精细结构交给图像通道。

三个常见误用要避开。第一个是用真实核做初始化去“测算法”，得到的收敛轨迹完全没有参考意义，盲的起点和真值不能太接近。第二个是把彩色图像三个通道分开做盲去卷积，三个核必然漂移出不同结果，正确做法是转 YUV 只处理 Y 亮度通道，或者假设三个通道共享同一核，联合估计一个核之后再通道独立做非盲去卷积。第三个是只用平滑纹理图测试，自然图梯度稀疏性没有被激活，算法很容易逃回模糊解，测试图里必须混入文字、边缘或随机纹理。

5. 盲去卷积进阶：去振铃、边界处理与自相关估计PSF尺寸

5.1 振铃的来源：边界截断与噪声放大

振铃是盲去卷积最常见的失败形态，原因有两个。一是频域卷积默认循环边界，图像外侧像素与对侧发生假相关，边界突变被当成高频振荡放大；二是反卷积本质是高频放大器，噪声里那点能量会被成倍拉高。边界问题最直接的解是边缘锥化（edge taper）：把图像四周渐变为均值，或者用反射边界扩展后裁剪掉外圈，让卷积窗口内没有硬跳变。正规则项方面，振铃越严重，优先增大 lam_tv 而不是减少 lr_x。限制 outer 次数同样有效，PSNR 曲线见顶后继续迭代只会往振铃方向漂移，用早停选中间结果比盲目多跑几轮更实际。

5.2 用自相关粗估PSF支持域：先定尺寸再进循环

psf_size 选错比正则权重选错后果更严重，所以动手前需要用一种廉价方式估一下。对高细节图像，清晰图自相关近似为尖峰，观测图的自相关就近似等于 PSF 的自相关；自相关主瓣的负波瓣半径能给出核支持的粗略上限。

def autocorr2d(a): a = a - a.mean() fa = np.fft.fft2(a) c = np.real(np.fft.ifft2(fa * np.conj(fa))) return np.fft.fftshift(c) c = autocorr2d(y) idx = np.unravel_index(np.argmin(c), c.shape) radius = int(np.hypot(idx[0] - c.shape[0] // 2, idx[1] - c.shape[1] // 2)) print('suggested psf radius:', max(3, radius))

原理是自相关峰旁边那一圈负旁瓣位置与核尺度正相关，argmin 找到最负的点，它到中心的距离就是核半径的粗略估计。这个方法在轻微失焦图上能给出 3~8 像素的合理范围，在强运动模糊图上会偏大，但至少可以排除把 psf_size=20 直接扔进去的盲目做法。得到 radius 后，把 psf_size 设为 2*radius+1，正好覆盖核支持域。

5.3 验证流程：先跑仿真，再上真实图像

盲去卷积里最容易翻车的环节是“自我感觉良好”，因为人眼对清晰化图像有天然偏好。标准验证流程是：用已知核生成合成退化图，跑同一条完整链路，分别计算 PSNR 和核相对误差；核误差压到 0.1 以下，再换用多张不同纹理的测试图，确认不是过拟合到单张图。真实图像上没有 ground truth，只能核对残差平稳性和核形态，所以仿真验证这一步不能省。把模糊图自相关跑一遍，先用负瓣半径把 psf_size 定下来，再进交替循环，是现在最省事的盲去卷积起步做法。

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