CN106023106A - 一种自适应参数更新的贝叶斯图像反卷积恢复方法 - Google Patents

Abstract

本发明属于图像处理技术领域，公开了一种自适应参数更新的贝叶斯图像反卷积恢复方法，包括：观测模型建模；图像先验建模；求观测模型的下界：利用优化最小化方法获得观测模型的下界函数；求图像先验模型的下界：利用优化最小化方法获得图像先验模型的下界函数；利用最大化后验概率迭代更新观测模型超参数、图像先验超参数和图像；本发明提供的自适应参数更新的贝叶斯图像反卷积恢复方法自适应更新能力好，计算代价低，恢复效果显著，具有去除高斯噪声和椒盐噪声的能力，可应用于在高斯噪声或椒盐噪声污染情形下且受高斯模糊退化图像的恢复。

Description

一种自适应参数更新的贝叶斯图像反卷积恢复方法

技术领域

本发明涉及图像处理技术领域，尤其涉及一种自适应参数更新的贝叶斯图像反卷积恢复方法，用于当图像受到高斯噪声或者椒盐噪声污染的情形下，对图像进行贝叶斯反卷积恢复。

背景技术

图像反卷积是计算机视觉和图像处理中的一个基础和关键问题。在实际成像过程中，由于受湍流、运动、抖动等多种因素的影响，观测图像普遍包含湍流、运动、抖动等模糊退化，如相机在成像过程中由于手的抖动导致的运动模糊。

此外，由于电磁干扰及通信***的故障和缺陷等因素，图像数据常常受到噪声的干扰，例如电视信号会由于大气层的干扰而产生脉冲噪声(椒盐噪声)、电子电路的热效应和光电传感器的光子波动等会产生加性噪声(如高斯噪声)，等等。很多成像应用***的处理结果很大程度上都依赖于反卷积后图像的质量，所以图像反卷积预处理是十分必须的、也是很多图像预处理应用***的一个重要组成部分。

现有变分正则化图像反卷积方法往往能有效去除噪声和模糊，但方法中涉及的关键正则化参数需要手动选取，方法耗时且不实用。近年来，变分贝叶斯方法已经用于图像反卷积，现有贝叶斯图像反卷积方法往往仅仅去除高斯噪声和模糊，性能受限，如D.Babacan等人于2008年提出的变分贝叶斯总变分图像反卷积方法仅仅能处理高斯噪声和模糊退化，对脉冲噪声污染退化无能为力。

发明内容

针对上述问题，本发明的目的在于提供一种自适应参数更新的贝叶斯图像反卷积恢复方法，既能处理高斯噪声又能处理脉冲噪声，解决统一模型下多种噪声和模糊的有效去除问题。

本发明实现的技术思路为：先将观测模型建模为广义高斯分布；再将图像先验模型建模为广义拉普拉斯分布，通过优化最小化方法分别求他们的下界，最后通过最大化后验概率迭代自适应更新超参数和求解图像，具有去除高斯噪声和椒盐噪声和自适应参数更新能力。

为达到上述目的，本发明的实施例采用如下技术方案予以实现。

一种自适应参数更新的贝叶斯图像反卷积恢复方法，所述方法包括如下步骤：

步骤1，获取观测图像y，所述观测图像包含高斯噪声和高斯模糊，或者所述观测图像包含椒盐噪声和高斯模糊；所述观测图像y＝N_im(Hx)+n，其中n表示高斯噪声，N_im表示椒盐噪声，H表示模糊核函数矩阵，x是待恢复的清晰图像，y是观测图像；

步骤2，采用零均值的广义高斯分布对观图像y中的噪声进行建模，从而得到观测图像y的条件概率分布函数p(y|β，x)，p(ο|·)表示已知·时ο的条件概率分布，β为多变量广义高斯观测模型的超参数，β服从均匀分布；

步骤3，对所述待恢复的清晰图像x进行先验建模，得到待恢复的清晰图像x的先验概率分布函数p(x|α)，α为图像先验模型超参数，α服从均匀分布；

步骤4，在分层贝叶斯模型下，构建联合全局概率分布p(α，β，x，y)：

p(α，β，x，y)＝p(α)p(β)p(y|β，x)p(x|α)

步骤5，计算观测图像y的条件概率分布函数p(y|β，x)的下界M₁(β，x，y)：p(y|β，x)≥M₁(β，x，y)；计算待恢复的清晰图像x的先验概率分布函数p(x|α)的下界M₂(α，x)：p(x|α)≥M₂(α，x)；

步骤6，根据观测图像y的条件概率分布函数p(y|β，x)的下界、待恢复的清晰图像x的先验概率分布函数p(x|α)的下界，得到所述联合全局概率分布p(α，β，x，y)的下界p(α)p(β)M₁(β，x，y)M₂(α，x)：

p(α，β，x，y)＝p(α)p(β)p(y|β，x)p(x|α)≥p(α)p(β)M₁(β，x，y)M₂(α，x)

步骤7，求解后验分布p(α，β，x|y)＝p(α，β，x，y)/p(y)的下界，得到关于图像先验模型超参数α、多变量广义高斯观测模型的超参数β、待恢复的清晰图像x的泛函；

步骤8，迭代求解关于图像先验模型超参数α、多变量广义高斯观测模型的超参数β、待恢复的清晰图像x的泛函，直到满足迭代停止条件，并将最后一次迭代得到的待恢复的清晰图像作为反卷积恢复的清晰图像。

本发明技术方案的有益效果：

(1)本发明提供的贝叶斯图像反卷积方法，将观测模型和图像先验模型分别建模为广义高斯分布和广义拉普拉斯分布，具有更大兼容性，既能处理高斯噪声又能处理椒盐噪声；(2)本发明提供的贝叶斯图像反卷积方法，利用优化最小化方法分别求观测模型和图像先验模型概率分布的下界，通过最大化后验概率，得到观测模型和图像先验模型超参数的自适应更新公式，相比于现有文献中手动选取参数，本发明提供的方法在迭代过程中自适应更新超参数，具有更强的实用价值；有利于贝叶斯图像反卷积方法走向实用。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的贝叶斯图像反卷积恢复方法的流程示意图；

图2为本发明实施例中的退化观测图像和恢复图像一；

其中，图2(a)和图2(c)是实施例中受高斯噪声污染的退化图像；图2(b)和图2(d)是实施例中恢复的图像；

图3为本发明实施例中的退化观测图像和恢复图像二；图3(a)和图3(c)是实施例中受椒盐噪声污染的退化图像；图3(b)和图像3(d)是实施例中恢复的图像。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

一种自适应参数更新的贝叶斯图像反卷积恢复方法，参考如图1，所述方法包括如下步骤：

步骤1，获取观测图像y，所述观测图像包含高斯噪声和高斯模糊，或者所述观测图像包含椒盐噪声和高斯模糊；所述观测图像y＝N_imp(Hx)+n，其中n表示高斯噪声，N_imp表示椒盐噪声，H表示模糊核函数矩阵，x是待恢复的清晰图像，y是观测图像。

其中n大小：N×1，H大小：N×N，x大小：N×1，y是观测退化图y大小：N×1，N是观测图像的像素个数。

当n≠0，N_imp＝1时，表示观测图像受到高斯噪声污染；

当n＝0，N_imp≠1时，表示观测图像受到椒盐噪声污染。

步骤2，采用零均值的广义高斯分布对观图像y中的噪声进行建模，从而得到观测图像y的条件概率分布函数p(y|β，x)，p(ο|·)表示已知·时ο的条件概率分布，β为多变量广义高斯观测模型的超参数，β服从均匀分布。

采用零均值的广义高斯分布对观测图像y中的噪声进行建模，得到观测图像y的条件概率分布函数p(y|β，x)：

$p\left(y\right|\mathrm{\β},x)=K_1{\mathrm{\β}}^{\frac{{\mathrm{\λ}}_{1}N}{p}}\exp (-\frac{\mathrm{\β}}{p}||Hx-y|{|}_p^p)$

其中，λ₁是正的实数，0＜p≤2，N是观测图像y的像素个数，exp(·)表示以e为底的指数函数，β是多变量广义高斯观测模型的超参数，||·||_p表示p范数处理，表示p范数的p次方处理，p(ο|·)表示已知·时ο的条件概率分布，K₁为正的常数，H表示模糊核函数矩阵，x是待恢复的清晰图像，y是观测图像。

步骤3，对所述待恢复的清晰图像x进行先验建模，得到待恢复的清晰图像x的先验概率分布函数p(x|α)，α为图像先验模型超参数，α服从均匀分布。

对所述待恢复的清晰图像x进行先验建模，得到待恢复的清晰图像x的先验概率分布函数p(x|α)：

其中，λ₂是正的实数，0＜q≤2，α是图像先验模型超参数，exp(·)表示以e为底的指数函数，D₁、D₂、D₃、D₄、D₅、D₆分别表示0°、45°、90°、135°方向一阶差分算子[1-1]、[0 1；-1 0]、[1；-1]、[10；0-1]以及水平和竖直方向二阶差分算子[1 -2 1]、[1；-2；1]对应的差分矩阵，||·||_q表示q范数处理，表示q范数的q次方处理，K₂为正的常数。

需要补充的是，超参数使用扁平非正常超先验建模超参数分布使得其中const.是constant的缩写，表示常数，扁平先验是指该先验在定义域上是均匀分布。

步骤4，在分层贝叶斯模型下，构建联合全局概率分布p(α，β，x，y)：

p(α，β，x，y)＝p(α)p(β)p(y|β，x)p(x|α)

图像反卷积问题的贝叶斯推断是基于未知后验分布p(α，β，x|y)的估计p(α，β，x|y)＝p(α，β，x，y)/p(y)，由于观测模型和图像先验的非凸非平滑特点，直接最大化后验概率p(α，β，x|y)比较困难，因此，本发明实施例使用优化最小化(majorization-minimization，MM)方法来解决上述问题，首先需要求非凸先验的下界。

步骤5，计算观测图像y的条件概率分布函数p(y|β，x)的下界M₁(β，x，y)：p(y|β，x)≥M₁(β，x，y)；计算待恢复的清晰图像x的先验概率分布函数p(x|α)的下界M₂(α，x)：p(x|α)≥M₂(α，x)。

已知不等式其等价表达式为t＞0，z＞0，0＜p≤2；

设因此有

其中，(Hx-y)_i表示向量Hx-y的第i个元素，z_i表示N维变分列向量z的第i个元素，且设有函数：

$M_1(\mathrm{\β},x,y)={\mathrm{\β}}^{\frac{{\mathrm{\λ}}_{1}N}{p}}\exp (-\frac{\mathrm{\β}}{2}\underset{i=1}{\overset{N}{\Σ}}\frac{{(Hx-y)}_i^2+\[(2-p)/p\]z_i}{z_i^{1-p/2}})$

$M_2(\mathrm{\α},x)={\mathrm{\α}}^{\frac{{\mathrm{\λ}}_{2}N}{q}}\exp (-\frac{\mathrm{\α}q}{2}\underset{i=1}{\overset{N}{\Σ}}\frac{\underset{j=1}{\overset{6}{\Σ}}{\left(D_{j}x\right)}_i^2+\[(2-q)/q\]u_i}{u_i^{1-q/2}})$

N是观测图像y或待恢复的清晰图像x的像素个数，u_i表示N维变分列向量u的第i个元素，λ₁是正的实数，0＜p≤2，exp(·)表示以e为底的指数函数，β是多变量广义高斯观测模型的超参数，H表示模糊核函数矩阵，λ₂是正的实数，0＜q≤2，α是图像先验模型超参数，D_j表示差分矩阵；

观测图像y的条件概率分布函数p(y|β，x)的下界为M₁(β，x，y)：

$p\left(y\right|\mathrm{\β},x)=K_1{\mathrm{\β}}^{\frac{{\mathrm{\λ}}_{1}N}{p}}\exp (-\frac{\mathrm{\β}}{p}||Hx-y|{|}_p^p)\≥K_1{M}_1(\mathrm{\β},x,y)$

待恢复的清晰图像x的条件概率分布函数p(x|α)的下界为M₂(α，x)：

$p\left(x\right|\mathrm{\α})=K_2{\mathrm{\α}}^{\frac{{\mathrm{\λ}}_{2}N}{q}}\exp (-\mathrm{\α}|\left|\sqrt{\underset{i=1}{\overset{6}{\Σ}}{\left(D_{i}x\right)}^2}\right|{|}_q^q)\≥K_2{M}_2(\mathrm{\α},x).$

步骤6，根据观测图像y的条件概率分布函数p(y|β，x)的下界、待恢复的清晰图像x的先验概率分布函数p(x|α)的下界，得到所述联合全局概率分布p(α，β，x，y)的下界p(α)p(β)M₁(β，x，y)M₂(α，x)：

p(α，β，x，y)＝p(α)p(β)p(y|β，x)p(x|α)≥p(α)p(β)M₁(β，x，y)M₂(α，x)

步骤7，求解后验分布p(α，β，x|y)＝p(α，β，x，y)/p(y)的下界，得到关于图像先验模型超参数α、多变量广义高斯观测模型的超参数β、待恢复的清晰图像x的泛函。

最大化后验分布p(α，β，x|y)＝p(α，β，x，y)/p(y)的下界等价于：

$\begin{array}{c}\stackrel{\‾}{\mathrm{\Θ}}=\underset{\mathrm{\Θ}}{\arg \min }\{\frac{\mathrm{\β}}{2}\underset{i=1}{\overset{N}{\Σ}}\frac{{(Hx-y)}_i^2+\[(2-p)/p\]z_i}{z_i^{1-p/2}}+\frac{\mathrm{\α}q}{2}\underset{i=1}{\overset{N}{\Σ}}\frac{\underset{j=1}{\overset{6}{\Σ}}{\left(D_{j}x\right)}_i^2+\[(2-q)/q\]u_i}{u_i^{1-q/2}}\\ -\frac{{\mathrm{\λ}}_{1}N}{p}\log \mathrm{\β}-\frac{{\mathrm{\λ}}_{2}N}{q}\log \mathrm{\α}\}\end{array}$

其中Θ＝{α，β，x}，log表示对数运算符，表示求泛函最小时对应的变量Θ的值N是观测图像y或待恢复的清晰图像x的像素个数，u_i表示N维变分列向量u的第i个元素，λ₁是正的实数，0＜p≤2，exp(·)表示以e为底的指数函数，β是多变量广义高斯观测模型的超参数，H表示模糊核函数矩阵，λ₂是正的实数，0＜q≤2，α是图像先验模型超参数，D_j表示差分矩阵。

由上述最小化泛函可看出，利用优化最小化准则将原始观测模型和先验模型中含有的非凸先验用一系列凸泛函代替，并引入变分向量z和u。于是通过交替最小化方法可迭代地求解上述泛函中的未知量α、β、x、z和u。

步骤8，迭代求解关于图像先验模型超参数α、多变量广义高斯观测模型的超参数β、待恢复的清晰图像x的泛函，直到满足迭代停止条件，并将最后一次迭代得到的待恢复的清晰图像作为反卷积恢复的清晰图像。

步骤8具体包括如下子步骤：

(8a)初始化参数p，q，λ₁，λ₂，β⁰，α⁰，i＝1，...，N，迭代次数k初值为0；

(8b)对于第k次迭代的结果z^k，u^k，利用下式求解关于待恢复的清晰图像x的最小化问题：

$\underset{x}{\min }\left\{\frac{{\mathrm{\β}}^k}{2}\underset{i=1}{\overset{N}{\Σ}}\frac{{(Hx-y)}_i^2+\[(2-p)/p\]z_i^k}{{\left(z_i^k\right)}^{1-p/2}}+\frac{{\mathrm{\α}}^{k}q}{2}\underset{i=1}{\overset{N}{\Σ}}\frac{\underset{j=1}{\overset{6}{\Σ}}{\left(D_j{x}^k\right)}_i^2+\[(2-q)/q\]u_i^k}{{\left(u_i^k\right)}^{1-q/2}}\right\}$

可得关于x的更新式为：

$({\mathrm{\βH}}^T{W}_{F}H+\mathrm{\α}\underset{i=1}{\overset{6}{\Σ}}{D}_j^T{W}_R{D}_j)x^{k+1}={\mathrm{\βH}}^T{W}_{F}y$

其中k表示迭代次数，diag(·)表示对角矩阵处理；在迭代的过程中，为了防止z和u取值过小而使得W_F和W_R取到无穷大的不合理值，作如下处理：

其中

这里ε₁＝0.001。

(8c)利用下式求解关于变分向量z的最小化问题：

$\underset{z}{min}\left\{\underset{i=1}{\overset{N}{\Σ}}\frac{{({\mathrm{Hx}}^k-y)}_i^2+\[(2-p)/p\]z_i}{z_i^{1-p/2}}\right\}$

可得关于z的更新式为i＝1，...，N；

(8d)利用下式求解关于变分向量u的最小化问题：

$\underset{u}{\min }\left\{\underset{i=1}{\overset{N}{\Σ}}\frac{\underset{j=1}{\overset{6}{\Σ}}{\left(D_j{x}^k\right)}_i^2+\[(2-q)/q\]u_i}{u_i^{1-q/2}}\right\}$

可得关于u的更新式为i＝1，...，N；

(8e)利用下式更新多变量广义高斯观测模型的超参数β：

$\underset{\mathrm{\β}}{\min }\frac{\mathrm{\β}}{2}\underset{i=1}{\overset{N}{\Σ}}\frac{{({\mathrm{Hx}}^k-y)}_i^2+\[(2-p)/p\]z_i^i}{{\left(z_i^k\right)}^{1-p/2}}-\frac{{\mathrm{\λ}}_{1}N}{p}\log \mathrm{\β}$

可得关于β的迭代更新式：

${\mathrm{\β}}^{k+1}=\frac{{\mathrm{\λ}}_{1}N}{\underset{i=1}{\overset{N}{\Σ}}|{\mathrm{Hx}}^k-y{|}_i^p}$

(8f)利用下式更新图像先验模型超参数α：

$\underset{\mathrm{\α}}{\min }\frac{\mathrm{\α}q}{2}\underset{i=1}{\overset{N}{\Σ}}\frac{\underset{j=1}{\overset{6}{\Σ}}{\left(D_j{x}^k\right)}_i^2+\[(2-q)/q\]u_i^k}{{\left(u_i^k\right)}^{1-q/2}}-\frac{{\mathrm{\λ}}_{2}N}{q}log\mathrm{\α}$

可得关于α的迭代更新式：

${\mathrm{\α}}^{k+1}=\frac{{\mathrm{\λ}}_{2}N}{q\underset{i=1}{\overset{N}{\Σ}}{\[\underset{j=1}{\overset{6}{\Σ}}{\left(D_{j}x\right)}_i^2\]}^{\frac{q}{2}}}$

(8g)令迭代次数k加1，并依次重复子步骤(8b)至(8f)，直到||x^k-x^k+1||_F/||x^k||_F≤ε；或迭代次数k达到预设的最大迭代次数；

示例性的，迭代终止条件ε＝10^-6。

(8h)获取最后一次迭代得到的待恢复的清晰图像，并将其作为反卷积恢复的清晰图像。

在本实施例中，参数λ₁和λ₂需要去逼近与先验分布p(y|β，x)和p(x|α)相应的划分函数。但是，与先验分布p(y|β，x)和p(x|α)相应的划分函数在解析上难以处理，因此，有必要近似处理先验分布p(y|β，x)和p(x|α)的划分函数。近似地确定先验分布p(y|β，x)和p(x|α)的划分函数后，就可以近似地确定λ₁和λ₂的值。且λ₁和λ₂的值一旦选定，就保持不变。

示例性的，

当n≠0，N_imp＝1时，表示图像受高斯噪声污染。此时，可以设置p＝2，q＝1，λ₁＝1.2，λ₂＝0.4。

当n＝0，N_imp≠1时，表示图像受椒盐脉冲噪声污染。此时，可以设置p＝1，q＝1，λ₁＝1.2，λ₂＝0.4。

在本实施例中，部分参数的初始化值如下：迭代开始时，首先迭代初始化图像x⁰为输入的退化图像y。初始化i＝1，...，N；β⁰＝10^-3，α⁰＝10^-3；最大迭代次数maxIter＝500；上述参数初始化值对所有测试图像均适用。

为了验证本发明提供的贝叶斯图像反卷积方法的有效性，图2中的图2(a)和图2(c)表示受高斯模糊和高斯噪声污染的退化图像，提出方法应用到图2(a)和图2(c)获得相应的恢复图像为图2(b)和图2(d)，其中dB是decibel的缩写，表示分贝，BSNR是blursignal-to-noise ratio的缩写，表示模糊信噪比，PSNR是peak signal-to-noise ratio的缩写，表示峰值信噪比，BSNR和PSNR的单位均是dB。图3中的图3(a)和图3(c)表示受高斯模糊和椒盐噪声污染的退化图像，提出方法应用到图3(a)和图3(c)获得相应的恢复图像为图3(b)和图3(d)。

从图2和图3的这2组退化图像可看出，图像被模糊和噪声退化后，图像本身细节基本被淹没，从视觉上难以分辨，这为图像后续处理带来很大的困难。从图2(b)和图2(d)以及图3(b)和图3(d)可看出，本发明提供的贝叶斯反卷积方法有效地去除图像模糊，较好地保留了图像本身的细节结构，也最大程度地抑制图像强噪声。这从主观视觉验证了本发明提供的方法的有效性。

表1列出了本发明给出的方法应用到四种不同噪声级别污染的退化图像恢复的PSNR结果对比。从这个表可看出，提出的方法恢复结果的PSNR相比退化图像的PSNR获得了明显的提高，此外，本发明提出的方法的PSNR也高于2008年Babacan等人提出的方法，且Babacan等人提出的方法很难去除椒盐噪声。这再一次从客观评价指标上验证了本发明方法的有效性。

综上所述，本发明公开了一种分层贝叶斯图像反卷积恢复方法，首先对观测模型和图像进行先验建模，然后再对超参数进行建模；其次为了克服观测模型和图像先验模型中非凸泛函数值上难以求解的难题，本发明利用优化最小化方法将上述问题转化为易于求解的非凸先验的下界，然后再最大化后验分布结合交替最小化的思想以迭代的方式得到想要的解。上述实验结果表明本发明提出的自适应参数更新的方法不需要任何人工手动调整参数，而且可以有效地抑制高斯噪声和椒盐噪声同时去除模糊。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (7)

1.一种自适应参数更新的贝叶斯图像反卷积恢复方法，其特征在于，所述方法包括如下步骤：

步骤1，获取观测图像y，所述观测图像包含高斯噪声和高斯模糊，或者所述观测图像包含椒盐噪声和高斯模糊；所述观测图像y＝N_imp(Hx)+n，其中n表示高斯噪声，N_imp表示椒盐噪声，H表示模糊核函数矩阵，x是待恢复的清晰图像，y是观测图像；

步骤2，采用零均值的广义高斯分布对观测图像y中的噪声进行建模，从而得到观测图像y的条件概率分布函数p(y|β,x)，β为多变量广义高斯观测模型的超参数，β服从均匀分布；

步骤3，对所述待恢复的清晰图像x进行先验建模，得到待恢复的清晰图像x的先验概率分布函数p(x|α)，α为图像先验模型超参数，α服从均匀分布；

步骤4，在分层贝叶斯模型下，构建联合全局概率分布p(α,β,x,y)：

p(α,β,x,y)＝p(α)p(β)p(y|β,x)p(x|α)

步骤5，计算观测图像y的条件概率分布函数p(y|β,x)的下界M₁(β,x,y)：p(y|β,x)≥M₁(β,x,y)；计算待恢复的清晰图像x的先验概率分布函数p(x|α)的下界M₂(α,x)：p(x|α)≥M₂(α,x)；

步骤6，根据观测图像y的条件概率分布函数p(y|β,x)的下界、待恢复的清晰图像x的先验概率分布函数p(x|α)的下界，得到所述联合全局概率分布p(α,β,x,y)的下界p(α)p(β)M₁(β,x,y)M₂(α,x)：

p(α,β,x,y)＝p(α)p(β)p(y|β,x)p(x|α)≥p(α)p(β)M₁(β,x,y)M₂(α,x)

步骤7，求解后验分布p(α,β,x|y)＝p(α,β,x,y)/p(y)的下界，得到关于图像先验模型超参数α、多变量广义高斯观测模型的超参数β、待恢复的清晰图像x的泛函；

步骤8，迭代求解关于图像先验模型超参数α、多变量广义高斯观测模型的超参数β、待恢复的清晰图像x的泛函，直到满足迭代停止条件，并将最后一次迭代得到的待恢复的清晰图像作为反卷积恢复的清晰图像。

2.根据权利要求1所述的一种自适应参数更新的贝叶斯图像反卷积恢复方法，其特征在于，步骤1中所述观测图像y＝N_imp(Hx)+n：

当n≠0，N_imp＝1时，表示观测图像包含高斯噪声；

当n＝0，N_imp≠1时，表示观测图像包含椒盐噪声。

3.根据权利要求1所述的一种自适应参数更新的贝叶斯图像反卷积恢复方法，其特征在于，步骤2具体为：

采用零均值的广义高斯分布对观测图像y中的噪声进行建模，得到观测图像y的条件概率分布函数p(y|β,x)：

$p\left(y\right|\mathrm{\β},x)=K_1{\mathrm{\β}}^{\frac{{\mathrm{\λ}}_{1}N}{p}}\exp (-\frac{\mathrm{\β}}{p}||Hx-y|{|}_p^p)$

其中，λ₁是正的实数，0<p≤2，N是观测图像y的像素个数，exp(·)表示以e为底的指数函数，β是多变量广义高斯观测模型的超参数，||·||_p表示p范数处理，表示p范数的p次方处理，K₁为正的常数，H表示模糊核函数矩阵，x是待恢复的清晰图像。

4.根据权利要求1所述的一种自适应参数更新的贝叶斯图像反卷积恢复方法，其特征在于，步骤3具体为：

对所述待恢复的清晰图像x进行先验建模，得到待恢复的清晰图像x的先验概率分布函数p(x|α)：

$p\left(x\right|\mathrm{\&alpha;})=K_2{\mathrm{\&alpha;}}^{\frac{{\mathrm{\&lambda;}}_{2}N}{q}}\exp (-\mathrm{\&alpha;}|\left|\sqrt{\underset{j=1}{\overset{6}{\&Sigma;}}{\left(D_{j}x\right)}^2}\right|{|}_q^q)$

其中，λ₂是正的实数，0<q≤2，α是图像先验模型超参数，N是观测图像y的像素个数，D_j表示差分矩阵，||·||_q表示q范数处理，表示q范数的q次方处理，K₂为正的常数。

5.根据权利要求1所述的一种自适应参数更新的贝叶斯图像反卷积恢复方法，其特征在于，步骤5具体为：

已知不等式其等价表达式为t>0，z>0，0<p≤2；

设因此有

其中，(Hx-y)_i表示向量Hx-y的第i个元素，z_i表示N维变分列向量z的第i个元素，且设有函数：

$M_1(\mathrm{\&beta;},x,y)={\mathrm{\&beta;}}^{\frac{{\mathrm{\&lambda;}}_{1}N}{p}}\exp (-\frac{\mathrm{\&beta;}}{2}\underset{i=1}{\overset{N}{\&Sigma;}}\frac{{(Hx-y)}_i^2+\&lsqb;(2-p)/p\&rsqb;z_i}{z_i^{1-p/2}})$

$M_2(\mathrm{\&alpha;},x)={\mathrm{\&alpha;}}^{\frac{{\mathrm{\&lambda;}}_{2}N}{q}}\exp (-\frac{\mathrm{\&alpha;}q}{2}\underset{i=1}{\overset{N}{\&Sigma;}}\frac{\underset{j=1}{\overset{6}{\&Sigma;}}{\left(D_{j}x\right)}_i^2+\&lsqb;(2-q)/q\&rsqb;u_i}{u_i^{1-q/2}})$

N是观测图像y或待恢复的清晰图像x的像素个数，u_i表示N维变分列向量u的第i个元素，λ₁是正的实数，0<p≤2，exp(·)表示以e为底的指数函数，β是多变量广义高斯观测模型的超参数，H表示模糊核函数矩阵，λ₂是正的实数，0<q≤2，α是图像先验模型超参数，D_j表示差分矩阵；

观测图像y的条件概率分布函数p(y|β,x)的下界为M₁(β,x,y)：

||·||_p表示p范数处理，表示p范数的p次方处理；

待恢复的清晰图像x的条件概率分布函数p(x|α)的下界为M₂(α,x)：

$p\left(x\right|\mathrm{\&alpha;})=K_2{\mathrm{\&alpha;}}^{\frac{{\mathrm{\&lambda;}}_{2}N}{q}}\exp (-\mathrm{\&alpha;}|\left|\sqrt{\underset{j=1}{\overset{6}{\&Sigma;}}{\left(D_{j}x\right)}^2}\right|{|}_q^q)\&GreaterEqual;K_2{M}_2(\mathrm{\&alpha;},x)$

||·||_q表示q范数处理，表示q范数的q次方处理。

6.根据权利要求1所述的一种自适应参数更新的贝叶斯图像反卷积恢复方法，其特征在于，步骤7具体为：

最大化后验分布p(α,β,x|y)＝p(α,β,x,y)/p(y)的下界等价于：

$\begin{array}{c}\stackrel{\&OverBar;}{\mathrm{\&Theta;}}=\underset{\mathrm{\&Theta;}}{\arg \min }\{\frac{\mathrm{\&beta;}}{2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{{\left(Hx-y\right)}_i^2+\&lsqb;\left(2-p\right)/p\&rsqb;z_i}{z_i^{1-p/2}}+\frac{\mathrm{\&alpha;}q}{2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{\underset{j=1}{\overset{6}{\mathrm{\&Sigma;}}}{\left(D_{j}x\right)}_i^2+\&lsqb;\left(2-q\right)/q\&rsqb;u_i}{u_i^{1-q/2}}\\ -\frac{{\mathrm{\&lambda;}}_{1}N}{p}\log \mathrm{\&beta;}-\frac{{\mathrm{\&lambda;}}_{2}N}{q}\log \mathrm{\&alpha;}\}\end{array}$

其中Θ＝{α,β,x}，log表示对数运算符，表示求泛函最小时对应的变量Θ的值N是观测图像y或待恢复的清晰图像x的像素个数，u_i表示N维变分列向量u的第i个元素，λ₁是正的实数，0<p≤2，exp(·)表示以e为底的指数函数，β是多变量广义高斯观测模型的超参数，H表示模糊核函数矩阵，λ₂是正的实数，0<q≤2，α是图像先验模型超参数，D_j表示差分矩阵。

7.根据权利要求6所述的一种自适应参数更新的贝叶斯图像反卷积恢复方法，其特征在于，步骤8具体包括如下子步骤：

(8a)初始化参数p，q，λ₁，λ₂，β⁰，α⁰，i＝1，…，N，迭代次数k初值为0；

(8b)对于第k次迭代的结果z^k，u^k，利用下式求解关于待恢复的清晰图像x的最小化问题:

$\underset{x}{\min }\left\{\frac{{\mathrm{\&beta;}}^k}{2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{{\left(Hx-y\right)}_i^2+\&lsqb;\left(2-p\right)/p\&rsqb;z_i^k}{{\left(z_i^k\right)}^{1-p/2}}+\frac{{\mathrm{\&alpha;}}^{k}q}{2}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{\underset{j=1}{\overset{6}{\mathrm{\&Sigma;}}}{\left(D_j{x}^k\right)}_i^2+\&lsqb;\left(2-q\right)/q\&rsqb;u_i^k}{{\left(u_i^k\right)}^{1-q/2}}\right\}$

可得关于x的更新式为:

$({\mathrm{\&beta;H}}^T{W}_{F}H+\mathrm{\&alpha;}\underset{i=1}{\overset{6}{\mathrm{\&Sigma;}}}{D}_j^T{W}_R{D}_j)x^{k+1}={\mathrm{\&beta;H}}^T{W}_{F}y$

其中k表示迭代次数，diag(·)表示对角矩阵处理；

(8c)利用下式求解关于变分向量z的最小化问题:

$\underset{z}{min}\left\{\underset{i=1}{\overset{N}{\&Sigma;}}\frac{{({\mathrm{Hx}}^k-y)}_i^2+\&lsqb;(2-p)/p\&rsqb;z_i}{z_i^{1-p/2}}\right\}$

可得关于z的更新式为i＝1，…，N；

(8d)利用下式求解关于变分向量u的最小化问题:

$\underset{u}{min}\left\{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{\underset{j=1}{\overset{6}{\mathrm{\&Sigma;}}}{\left(D_j{x}^k\right)}_i^2+\&lsqb;\left(2-q\right)/q\&rsqb;u_i}{u_i^{1-q/2}}\right\}$

可得关于u的更新式为i＝1，…，N；

(8e)利用下式更新多变量广义高斯观测模型的超参数β：

$\underset{\mathrm{\&beta;}}{\min }\frac{\mathrm{\&beta;}}{2}\underset{i=1}{\overset{N}{\&Sigma;}}\frac{{({\mathrm{Hx}}^k-y)}_i^2+\&lsqb;(2-p)/p\&rsqb;z_i^i}{{\left(z_i^k\right)}^{1-p/2}}-\frac{{\mathrm{\&lambda;}}_{1}N}{p}log\mathrm{\&beta;}$

可得关于β的迭代更新式:

${\mathrm{\&beta;}}^{k+1}=\frac{{\mathrm{\&lambda;}}_{1}N}{\underset{i=1}{\overset{N}{\&Sigma;}}|{\mathrm{Hx}}^k-y{|}_i^p}$

(8f)利用下式更新图像先验模型超参数α：

$\underset{\mathrm{\&alpha;}}{\min }\frac{\mathrm{\&alpha;}q}{2}\underset{i=1}{\overset{N}{\&Sigma;}}\frac{\underset{j=1}{\overset{6}{\&Sigma;}}{\left(D_j{x}^k\right)}_i^2+\&lsqb;(2-q)/q\&rsqb;u_i^k}{{\left(u_i^k\right)}^{1-q/2}}-\frac{{\mathrm{\&lambda;}}_{2}N}{q}log\mathrm{\&alpha;}$

可得关于α的迭代更新式：

${\mathrm{\&alpha;}}^{k+1}=\frac{{\mathrm{\&lambda;}}_{2}N}{q\underset{i=1}{\overset{N}{\&Sigma;}}{\&lsqb;\underset{j=1}{\overset{6}{\&Sigma;}}{\left(D_{j}x\right)}_i^2\&rsqb;}^{\frac{q}{2}}}$

(8g)令迭代次数k加1，并依次重复子步骤(8b)至(8f)，直到||x^k-x^k+1||_F/||x^k||_F≤ε；或迭代次数k达到预设的最大迭代次数；

(8h)获取最后一次迭代得到的待恢复的清晰图像，并将其作为反卷积恢复的清晰图像。