CN102708550B - 一种基于自然图像统计特性的盲去模糊方法 - Google Patents

Abstract

本发明公开了一种基于自然图像统计特性的盲去模糊方法，首先，确定模糊图像梯度的先验分布和点扩散函数的先验分布；其次，用近似后验分布来逼近真实后验分布，得到这些近似后验分布的参数更新等式；将图像进行金字塔分解，将低分辨率图层按参数更新等式迭代估计得到当前层的点扩散函数，用双阈值法对点扩散函数进行处理，然后将当前层的点扩散函数进行上采样的结果作为高分辨率图层迭代的初始值，直到最后一层迭代收敛到最佳估计；接着，对模糊图像解卷积得到清晰图像。本发明避免了噪声在更高分辨率图像中的放大，增加了结果的鲁棒性，提出的抑制振铃效应方法简单易行，在保持图像细节信息的同时，一定程度上降低了振铃效应对图像复原的影响。

Description

一种基于自然图像统计特性的盲去模糊方法

技术领域

本发明涉及一种基于自然图像统计特性的盲去模糊方法，可应用于图像采集过程中因相机抖动造成的图像模糊，属于计算机数字图像处理领域。

背景技术

随着现代数字技术的发展以及图像成像设备的普及，数字图像在国家安全，遥感，医学影像，交通监控和人类日常生活等领域中得到了广泛的应用。众所周知，在数字图像的采集过程中，相机等采集设备在快门开启的瞬间会产生不可避免的轻微抖动，这种抖动通常会使我们得到的图像是模糊的。模糊图像会对图像的进一步应用如图像分析、目标提取及识别等带来相当的困难，特别是对于图像获取来说，许多场景如交通监控等只是瞬间发生，无法重现，因此，图像去模糊技术在现代数字图像处理中显得尤为重要。

发明内容

本发明的目的在于提供一种基于自然图像统计特性的盲去模糊方法，

为实现上述目的，本发明所采用的技术方案是：

基于自然图像统计特性的盲去模糊方法包括如下步骤：

1)将相机抖动模糊图像表达为清晰图像与点扩散函数的卷积再加上噪声的形式，

$B=K\⊗L+N---\left(1\right)$

式(1)中，B表示模糊图像，K表示模糊函数，L表示未模糊清晰图像，N表示噪声，表示卷积操作，其中只有B已知，

图像梯度系数的直方图在零点处有大的峰值，而在远离零点处具有长的尾，在梯度域内对模糊图像进行复原，采用零均值高斯混合模型对未模糊清晰图像梯度进行建模，

根据卷积运算的性质，自然图像降质模型在梯度域内表示为

$\▿B=K\⊗\▿L+\▿N---\left(2\right)$

其中和分别表示模糊图像、清晰图像和噪声的梯度；

2)假定点扩散函数具有空间移不变性，即模糊图像的全图受到同一个点扩散函数的影响，选取图像的部分区域来代替全图进行点扩散函数的估计，选取的部分区域记为P，梯度域内的图像降质模型变为：

$\▿P=K\⊗\▿L_p+\▿N---\left(3\right)$

式中L_p表示在模糊图像中选取的区域所对应的清晰图像；

3)根据相机抖动模糊图像点扩展函数统计特性，点扩展函数概率密度分布类似于指数分布，采用混合指数分布对点扩展函数进行建模；

4)根据贝叶斯原理，得到清晰图像梯度点扩散函数K的后验概率：

$p(K,{\▿L}_p|\▿P)=\frac{p(\▿P|K,\▿L_p)p(K,{\▿L}_p)}{p(\▿P)}---\left(4\right)$

采用分布来逼近真实的后验概率分布通过近似分布q和真实后验概率分布之间的Kullback-Leibler散度的最小化来实现近似分布q的优化；

5)根据散度的计算定义贝叶斯变分方法的代价函数，然后根据代价函数和模糊图像梯度的先验分布和点扩散函数的先验分布，用变分贝叶斯期望最大化定理求出各个参数的近似后验概率分布；

6)通过变分期望最大化定理的变分最大化来推导各个参数的近似后验概率函数的分布参数的更新等式；

7)对模糊图像L_p进行金字塔分解，得到S层由分辨率低到高的图像金字塔，令第S层为最高层，即为图像分辨率最高的一层，为第一层选择模型初始化值，通过带入分布参数的更新等式中，反复迭代计算代价函数，直到代价函数收敛于一个设定的阈值，得到各个参数的最优解，由此估计出金字塔第1层的点扩散函数，记为K1，以及当前层清晰图像的梯度值接着采用双阈值方法对估计得到点扩散函数K1进行处理，然后将第一层清晰图像的梯度值利用双线性插值法放大到金字塔第二层大小得到第一层清晰图像梯度放大值同时也将点扩散函数放大到第二层点扩展函数大小得到一个新的点扩展函数K2'，同样的再次使用双阈值法对新的点扩展函数K2'进行处理，将第一层清晰图像梯度放大值和经过双阈值法处理的点扩展函数K2'作为第二层迭代的初始值，再次通过分布参数的更新等式的迭代估计得到金字塔第二层的点扩散函数K2和当前层清晰图像的梯度值以此类推，最终得到金字塔S层点扩展函数K，

双阈值方法特征为：

使用两个阈值t_low和t_high来抑制点扩散函数K中的噪声，其中t_high＞t_low，定义两个掩膜M_low和M_high

$M_{\mathrm{low}}\left(i\right)=\left\{\begin{array}{cc}1& k\left(i\right){\&GreaterEqual;t}_{\mathrm{low}}{k}_{\max }\\ 0& k\left(i\right)<t_{\mathrm{low}}{k}_{\max }\end{array}\right.$

$M_{\mathrm{high}}\left(i\right)=\left\{\begin{array}{cc}1& k\left(i\right){\&GreaterEqual;t}_{\mathrm{high}}{k}_{\max }\\ 0& k\left(i\right)<t_{\mathrm{high}}{k}_{\max }\end{array}\right.$

其中k_max表示PSF中的最大值，在得到两个掩膜后，以掩膜M_high中值为1的元素为中心，观察其8邻域中的点在掩膜M_low中的值，若掩膜M_low中的值为1则相应的使这个点在掩膜M_high中的值也为1，否则这个点在掩膜M_high中的值还是为0，判断完所有掩膜M_high中值为1的点为一次迭代，不断迭代直到掩膜M_high中点的值没有变化；

8)利用点扩散函数K，对模糊图像B用Richardson-Lucy算法进行解卷积得到清晰图像L；

9)利用模糊图像B得到图像的细节区，振铃区和平坦区，其中细节区包含图像细节，振铃区包含图像振铃效应所在区域，平坦区包含基本无细节存在的区域，振铃效应通常出现在细节区域附近的平坦区域，结合形态学图像处理划分得到三个区域掩膜M_r，掩膜M_p，掩膜M_d，分别表示振铃区掩膜，平坦区掩膜，细节区掩膜；

10)利用振铃区掩膜、平坦区掩膜、细节区掩膜，针对不同区域采用不同程度的模糊均值滤波器对图像L进行滤波处理得到最终模糊图像复原结果。

本发明的有益效果如下：

a)针对采集过程中因相机抖动引起的图像模糊，基于自然图像梯度和点扩散函数的先验统计特性，采用了变分贝叶斯方法估计点扩散函数，变分贝叶斯方法是最近几年发展起来的一种贝叶斯近似方法，其原理是用未知变量和参数的近似后验分布来逼近他们的真实分布，使贝叶斯方法能解析实现，能学习模型结构和模型参数，本发明充分利用变分贝叶斯估计方法在学习参数过程中避免过拟合的优点和模型选择的能力，来准确估计图像模糊模型的各参数；

b)为了避免结果收敛到局部最小值，对图像进行金字塔分解，在多尺度空间上对模糊核进行估计，利用低一级的模型估计结果作为高一级的初始值，最终得到最优化解；

c)如果低级模糊核估计不准确会导致低级模糊核中噪声在高级模糊核中放大最终影响图像去模糊的效果，本发明提出了一种双阈值方法对模糊核进行噪声抑制；

d)针对图像去模糊过程中的振铃效应，本发明提出了一种简单易行的减少振铃效应的方法，在保持图像边缘和细节信息的同时，一定程度上降低了振铃效应对图像去模糊的影响。

附图说明

图1是模糊自然图像；

图2是未抑制振铃前去模糊图像；

图3是表示有振铃区的图像；

图4是去模糊结果进行振铃效应抑制后图像。

具体实施方式

基于自然图像统计特性的盲去模糊方法包括如下步骤：

1)将相机抖动模糊图像表达为清晰图像与点扩散函数的卷积再加上噪声的形式，

$B=K\&CircleTimes;L+N---\left(1\right)$

式(1)中，B表示模糊图像，K表示模糊函数，L表示未模糊清晰图像，N表示噪声，表示卷积操作，其中只有B已知，

图像梯度系数的直方图在零点处有大的峰值，而在远离零点处具有长的尾，在梯度域内对模糊图像进行复原，采用零均值高斯混合模型对未模糊清晰图像梯度进行建模，

根据卷积运算的性质，自然图像降质模型在梯度域内表示为

$\&dtri;B=K\&CircleTimes;\&dtri;L+\&dtri;N---\left(2\right)$

其中和分别表示模糊图像、清晰图像和噪声的梯度；

2)假定点扩散函数具有空间移不变性，即模糊图像的全图受到同一个点扩散函数的影响，选取图像的部分区域来代替全图进行点扩散函数的估计，选取的部分区域记为P，梯度域内的图像降质模型变为：

$\&dtri;P=K\&CircleTimes;\&dtri;L_p+\&dtri;N---\left(3\right)$

式中L_p表示在模糊图像中选取的区域所对应的清晰图像；

3)根据相机抖动模糊图像点扩展函数统计特性，点扩展函数概率密度分布类似于指数分布，采用混合指数分布对点扩展函数进行建模；

4)根据贝叶斯原理，得到清晰图像梯度点扩散函数K的后验概率：

$p(K,{\&dtri;L}_p|\&dtri;P)=\frac{p(\&dtri;P|K,\&dtri;L_p)p(K,{\&dtri;L}_p)}{p(\&dtri;P)}---\left(4\right)$

采用分布来逼近真实的后验概率分布通过近似分布q和真实后验概率分布之间的Kullback-Leibler散度的最小化来实现近似分布q的优化；

5)根据散度的计算定义贝叶斯变分方法的代价函数，然后根据代价函数和模糊图像梯度的先验分布和点扩散函数的先验分布，用变分贝叶斯期望最大化定理求出各个参数的近似后验概率分布；

6)通过变分期望最大化定理的变分最大化来推导各个参数的近似后验概率函数的分布参数的更新等式；

7)对模糊图像L_p进行金字塔分解，得到S层由分辨率低到高的图像金字塔，令第S层为最高层，即为图像分辨率最高的一层，为第一层选择模型初始化值，通过带入分布参数的更新等式中，反复迭代计算代价函数，直到代价函数收敛于一个设定的阈值，得到各个参数的最优解，由此估计出金字塔第1层的点扩散函数，记为K1，以及当前层清晰图像的梯度值接着采用双阈值方法对估计得到点扩散函数K1进行处理，然后将第一层清晰图像的梯度值利用双线性插值法放大到金字塔第二层大小得到第一层清晰图像梯度放大值同时也将点扩散函数放大到第二层点扩展函数大小得到一个新的点扩展函数K2'，同样的再次使用双阈值法对新的点扩展函数K2'进行处理，将第一层清晰图像梯度放大值和经过双阈值法处理的点扩展函数K2'作为第二层迭代的初始值，再次通过分布参数的更新等式的迭代估计得到金字塔第二层的点扩散函数K2和当前层清晰图像的梯度值以此类推，最终得到金字塔S层点扩展函数K，

双阈值方法特征为：

使用两个阈值t_low和t_high来抑制点扩散函数K中的噪声，其中t_high＞t_low，定义两个掩膜M_low和M_high

$M_{\mathrm{low}}\left(i\right)=\left\{\begin{array}{cc}1& k\left(i\right){\&GreaterEqual;t}_{\mathrm{low}}{k}_{\max }\\ 0& k\left(i\right)<t_{\mathrm{low}}{k}_{\max }\end{array}\right.$

$M_{\mathrm{high}}\left(i\right)=\left\{\begin{array}{cc}1& k\left(i\right){\&GreaterEqual;t}_{\mathrm{high}}{k}_{\max }\\ 0& k\left(i\right)<t_{\mathrm{high}}{k}_{\max }\end{array}\right.$

其中k_max表示PSF中的最大值，在得到两个掩膜后，以掩膜M_high中值为1的元素为中心，观察其8邻域中的点在掩膜M_low中的值，若掩膜M_low中的值为1则相应的使这个点在掩膜M_high中的值也为1，否则这个点在掩膜M_high中的值还是为0，判断完所有掩膜M_high中值为1的点为一次迭代，不断迭代直到掩膜M_high中点的值没有变化；

8)利用点扩散函数K，对模糊图像B用Richardson-Lucy算法进行解卷积得到清晰图像L；

9)利用模糊图像B得到图像的细节区，振铃区和平坦区，其中细节区包含图像细节，振铃区包含图像振铃效应所在区域，平坦区包含基本无细节存在的区域，振铃效应通常出现在细节区域附近的平坦区域，结合形态学图像处理划分得到三个区域掩膜M_r，掩膜M_p，掩膜M_d，分别表示振铃区掩膜，平坦区掩膜，细节区掩膜；

10)利用振铃区掩膜、平坦区掩膜、细节区掩膜，针对不同区域采用不同程度的模糊均值滤波器对图像L进行滤波处理得到最终模糊图像复原结果。

实施例

1)将模糊图像表达为清晰图像与点扩散函数(point spread function PSF)的卷积再加上噪声的形式，

$B=K\&CircleTimes;L+N---\left(1\right)$

由近年关于自然图像统计特性的研究表明，自然图像的梯度基本服从重尾分布，即图像梯度系数的直方图在零点处有较大的峰值，而在远离零点处具有较长的尾，因此在梯度域内对模糊图像进行复原，采用零均值高斯混合模型对未模糊清晰图像梯度进行建模，

根据卷积运算的性质，自然图像降质模型可以在梯度域内表示为

$\&dtri;B=K\&CircleTimes;\&dtri;L+\&dtri;N---\left(2\right)$

假设中各个元素是独立同分布的，则

$p(\&dtri;L(i,j))=\underset{c=1}{\overset{C}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_c^{\left(L\right)}G(\&dtri;L(i,j)|0,v_c)---\left(3\right)$

其中表示原始清晰图像的梯度中元素，G表示零均值高斯分布，即(a是分布的均值，b>0是分布的逆方差或精度)，C表示该混合分布包括零均值高斯成分的总数，和v_c分别表示第c个零均值高斯成分的混合权重和逆方差，混合权重之和为1，即

在图像处理领域，通常假设噪声N为零均值高斯噪声，在噪声N为独立同分布的零均值噪声的前提下，也为独立同分布的零均值高斯噪声，因此，本发明假设噪声为零均值高斯噪声，且中各个元素为独立同分布，则噪声的概率密度分布为：

$p(i,j)=G(\&dtri;N(i,j)|0,\mathrm{\&gamma;})---\left(4\right);$

2)假定PSF具有空间移不变性，即模糊图像的全图受同一个点扩散函数的影响，因此为了提高算法的执行效率，同时保证PSF的估计精度，选取图像的部分区域来代替全图进行PSF的估计，选取的部分区域记为P，则梯度域内的图像降质模型可变为：

$\&dtri;P=K\&CircleTimes;\&dtri;L_p+\&dtri;N---\left(5\right)$

式中L_p表示在模糊图像中选取的区域所对应的清晰图像，大小为I×J，因此，式(3)相应的变为：

$p(\&dtri;L_p(i,j))=\underset{c=1}{\overset{C}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_c^{\left(L\right)}G(\&dtri;L_p(i,j)|0,v_c)---\left(6\right);$

3)根据相机抖动模糊图像点扩展函数统计特性，点扩展函数概率密度分布类似于指数分布，因此采用混合指数分布对点扩展函数进行建模，

假设点扩散函数K中各个元素是独立同分布的，则其概率密度分布为

$p\left(K(m,n)\right)=\underset{d=1}{\overset{D}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_d^{\left(K\right)}E\left(K(m,n)\right|{\mathrm{\&lambda;}}_d)---\left(7\right)$

其中K(m,n)表示点扩展函数K中的元素，点扩展函数K的大小为(2M_K+1)×(2N_K+1)，E表示指数分布，即 $E\left(x\right|a)=\left\{\begin{array}{cc}{\mathrm{ae}}^{-\mathrm{ax}},& x\&GreaterEqual;0\\ 0,& x<0\end{array}\right.,$ 式(6)中D表示指数混合分布的总数，和λ_d分别表示第d个指数分布的混合权重和尺度因子，混合权重之和为1，即

4)根据贝叶斯原理，得到清晰图像梯度和点扩散函数K的后验概率：

$p(K,{\&dtri;L}_p|\&dtri;P)=\frac{p(\&dtri;P|K,\&dtri;L_p)p(K,{\&dtri;L}_p)}{p(\&dtri;P)}---\left(8\right)$

本发明采用变分贝叶斯方法来解决式(8)的估计问题，用一个易计算的分布来逼近真实的后验概率分布通过近似分布q和真实后验概率分布p之间的Kullback-Leibler(KL)散度的最小化来实现近似分布q的优化，

首先，选取式(4)(6)(7)中各参数的先验分布：

混合高斯模型的权重的先验分布为Dirichlet分布

$p\left({\left\{{\mathrm{\&pi;}}_c^{\left(L\right)}\right\}}_{c=1}^C\right)=D\left({\left\{{\mathrm{\&pi;}}_c^{\left(L\right)}\right\}}_{c=1}^C\right|c^{\left(L\right)}),c=\mathrm{1,2},...,C---\left(9\right)$

混合高斯模型的逆方差的先验分布为Gamma分布

$p\left(v_c\right)=\mathrm{\&Gamma;}\left(v_c\right|a^{\left(L\right)},b^{\left(L\right)}),c=\mathrm{1,2},...,C---\left(10\right)$

混合指数分布权重的先验分布为Dirichlet分布

$p\left({\left\{{\mathrm{\&pi;}}_d^{\left(K\right)}\right\}}_{d=1}^D\right)=D\left({\left\{{\mathrm{\&pi;}}_d^{\left(K\right)}\right\}}_{d=1}^D\right|c^{\left(K\right)}),c=\mathrm{1,2},...D---\left(11\right)$

混合指数分布逆方差的先验分布为Gamma分布

$p\left({\mathrm{\&lambda;}}_d\right)=\mathrm{\&Gamma;}\left({\mathrm{\&lambda;}}_d\right|a_d^{\left(K\right)},b_d^{\left(K\right)}),d=\mathrm{1,2},...D---\left(12\right)$

噪声概率分布中的逆方差的先验分布为Gamma分布

p(γ)＝Γ(γ|a^(γ),b^(γ))   (13)

其中，Γ(·)表示Gamma分布，D(·)表示Dirichlet分布，c^(L)，a^(L)，b^(L)，a^(γ)，b^(γ)是超参数，

其次，将式(6)(7)中的模型参数写成集合的形式，即：

${\mathrm{\&pi;}}^{\left(L\right)}=\{{\mathrm{\&pi;}}_1^{\left(L\right)},{\mathrm{\&pi;}}_2^{\left(L\right)},...,{\mathrm{\&pi;}}_C^{\left(L\right)}\},v=\{v_1,v_2,...,v_C\},$

${\mathrm{\&pi;}}^{\left(K\right)}=\{{\mathrm{\&pi;}}_1^{\left(K\right)},{\mathrm{\&pi;}}_2^{\left(K\right)},...,{\mathrm{\&pi;}}_D^{\left(K\right)}\},\mathrm{\&lambda;}=\{{\mathrm{\&lambda;}}_1,{\mathrm{\&lambda;}}_2,...,{\mathrm{\&lambda;}}_D\},$

变分贝叶斯估计观测数据为模糊图像的梯度未知量用参数集合Θ表示，则参数集的真实后验概率密度用来表示，同理，近似分布q用参数集集合表示为q(Θ)，真实分布和近似后验分布之间的散度为：

$D_{\mathrm{KL}}\left(q\right|\left|p\right)=\underset{\mathrm{\&Theta;}}{\&Integral;}q\left(\mathrm{\&Theta;}\right)\log \left[\frac{q\left(\mathrm{\&Theta;}\right)}{p\left(\mathrm{\&Theta;}\right|\&dtri;p)}\right]\mathrm{d\&Theta;}={\&lang;\log \left[\frac{q\left(\mathrm{\&Theta;}\right)}{p\left(\mathrm{\&Theta;}\right|\&dtri;P)}\right]\&rang;}_q\&GreaterEqual;0---\left(14\right)$

其中〈·〉_q表示概率密度函数q(·)下求期望；

5)根据散度的计算定义贝叶斯变分方法的代价函数，然后根据代价函数和模糊图像梯度的先验分布和PSF的先验分布，用变分贝叶斯期望最大化定理求出各个参数的近似后验概率分布，

由贝叶斯定理

$p\left(\mathrm{\&Theta;}\right|\&dtri;P)=\frac{p(\&dtri;P|\mathrm{\&Theta;})p\left(\mathrm{\&Theta;}\right)}{p(\&dtri;P)}---\left(15\right)$

将(15)式带入(14)式，得：

$D_{\mathrm{KL}}\left(q\right|\left|p\right)=\underset{\mathrm{\&Theta;}}{\&Integral;}q\left(\mathrm{\&Theta;}\right)\log \left[\frac{q\left(\mathrm{\&Theta;}\right)p(\&dtri;P)}{p(\&dtri;p|\mathrm{\&Theta;})p\left(\mathrm{\&Theta;}\right)}\right]\mathrm{d\&Theta;}---\left(16\right)$

p(Θ)的积分与参数集Θ无关，因此式(16)中可将其移到积分号外面，由此得到代价函数：

$\begin{array}{c}{C}_{\mathrm{KL}}\left(q\right|\left|p\right)=D_{\mathrm{KL}}\left(q\right|\left|p\right)-\log p(\&dtri;P)=\underset{\mathrm{\&Theta;}}{\&Integral;}q\left(\mathrm{\&Theta;}\right)\log \left[\frac{q\left(\mathrm{\&Theta;}\right)}{p(\&dtri;P|\mathrm{\&Theta;})p\left(\mathrm{\&Theta;}\right)}\right]\mathrm{d\&Theta;}\\ ={\&lang;\log \left[\frac{q\left(\mathrm{\&Theta;}\right)}{p(\&dtri;P|\mathrm{\&Theta;})p\left(\mathrm{\&Theta;}\right)}\right]\&rang;}_q\&GreaterEqual;\log p(\&dtri;p)\end{array}---\left(17\right)$

最小化代价函数C_KL(q||p)等价于最下化KL散度D_KL(q||p)，

变分贝叶斯估计中通常假设各参数集中Θ中的参数是相互独立的，因此q(Θ)可以转变为以下形式：

$\begin{array}{c}q\left(\mathrm{\&Theta;}\right)=q\{{\&dtri;L}_p,K,\mathrm{\&gamma;},{\mathrm{\&pi;}}^{\left(L\right)},v,{\mathrm{\&pi;}}^{\left(K\right)},\mathrm{\&lambda;}\}\\ =q\left({\&dtri;L}_p\right)q\left(K\right)q\left(\mathrm{\&gamma;}\right)q\left({\mathrm{\&pi;}}^{\left(L\right)}\right)q\left(v\right)q\left({\mathrm{\&pi;}}^{\left(K\right)}\right)q\left(\mathrm{\&lambda;}\right)\end{array}---\left(18\right)$

进一步的代价函数式(17)可变为一些简单项之和的形式：

$\begin{array}{c}{C}_{\mathrm{KL}}=\left(q\right|\left|p\right)={\&lang;\log \frac{q(\&dtri;L_p)q\left({\mathrm{\&pi;}}^{\left(L\right)}\right)q\left(v\right)}{p({\&dtri;L}_p,{\mathrm{\&pi;}}^{\left(L\right)},v)}\&rang;}_q+{\&lang;\log \frac{q\left(K\right)q\left({\mathrm{\&pi;}}^{\left(K\right)}\right)q\left(\mathrm{\&lambda;}\right)}{p(K,{\mathrm{\&pi;}}^{\left(K\right)},\mathrm{\&lambda;})}\&rang;}_q\\ +{\&lang;\frac{q\left(\mathrm{\&gamma;}\right)}{p(\&dtri;P,\mathrm{\&gamma;}|\&dtri;L_p,K)}\&rang;}_q\\ =C_{\mathrm{KL}}^{\left(L\right)}\left(q\right|\left|p\right)+C_{\mathrm{KL}}^{\left(K\right)}\left(q\right|\left|p\right)+C_{\mathrm{KL}}^{\left(\mathrm{\&gamma;}\right)}\left(q\right|\left|p\right)\end{array}---\left(19\right)$

利用Jensen不等式进行进一步简化，由于清晰图像梯度的先验概率分布为混合高斯模型，则有：

$\begin{array}{c}\log p\left({\&dtri;L}_p(i,j)\right)=\log \underset{c=1}{\overset{C}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_c^{\left(L\right)}G(\&dtri;L_p(i,j)|0,v_c)\\ \&GreaterEqual;\underset{c=1}{\overset{C}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{\mathrm{ijc}}^{\left(L\right)}\log \frac{{\mathrm{\&pi;}}_c^{\left(L\right)}G\left({\&dtri;L}_p(i,j)\right|0{,v}_c)}{{\mathrm{\&lambda;}}_{\mathrm{ijc}}^{\left(f\right)}}\end{array}---\left(20\right)$

其中是一组权重参数，且对任意的i和j满足

由于点扩散函数的先验概率分布为混合指数分布，利用Jensen不等式化简之后有：

$\begin{array}{c}\log p\left(K(m,n)\right)=\log \underset{d=1}{\overset{D}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_d^{\left(L\right)}E\left(K(m,n)\right|{\mathrm{\&lambda;}}_d)\\ \&GreaterEqual;\underset{d=1}{\overset{D}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{\mathrm{mnd}}^{\left(K\right)}\log \frac{{\mathrm{\&pi;}}_d^{\left(K\right)}E\left(K(m,n)\right|{\mathrm{\&lambda;}}_d)}{{\mathrm{\&lambda;}}_{\mathrm{mnd}}^{\left(K\right)}}\end{array}---\left(21\right)$

其中是一组权重参数，且对任意的m和n满足

利用变分贝叶斯期望最大化得到Θ中各参数的近似后验概率分布为：

$q(\&dtri;L_p(i,j))=G\left({\&dtri;L}_p(i,j)\right|\&dtri;{\hat{L}}_p(i,j),\&dtri;{\tilde{L}}_p(i,j))---\left(22\right)$

$q\left({\left\{{\mathrm{\&pi;}}_c^{\left(L\right)}\right\}}_{c=1}^C\right)=D\left({\left\{{\mathrm{\&pi;}}_c^{\left(L\right)}\right\}}_{c=1}^C\right|{\left\{{\stackrel{-}{c}}_c^{\left(L\right)}\right\}}_{c=1}^C)---\left(23\right)$

$q\left(v_c\right)=\mathrm{\&Gamma;}\left(v_c\right|{\stackrel{-}{a}}_c^{\left(L\right)},{\stackrel{-}{b}}_c^{\left(L\right)})---\left(24\right)$

$q\left(K(m,n)\right)=G^{*}\left(K(m,n)\right|\hat{K}(m,n),\tilde{K}(m,n))---\left(25\right)$

$q\left({\left\{{\mathrm{\&pi;}}_d^{\left(K\right)}\right\}}_{d=1}^D\right)=D\left({\left\{{\mathrm{\&pi;}}_d^{\left(K\right)}\right\}}_{d=1}^D\right|{\left\{{\stackrel{-}{c}}_d^{\left(K\right)}\right\}}_{d=1}^D)---\left(26\right)$

$q\left({\mathrm{\&lambda;}}_d\right)=\mathrm{\&Gamma;}\left({\mathrm{\&lambda;}}_d\right|{\stackrel{-}{a}}_d^{\left(K\right)},{\stackrel{-}{b}}_d^{\left(K\right)})---\left(27\right)$

$q\left(\mathrm{\&gamma;}\right)=\mathrm{\&Gamma;}\left(\mathrm{\&gamma;}\right|{\stackrel{-}{a}}^{\left(\mathrm{\&gamma;}\right)},{\stackrel{-}{b}}^{\left(\mathrm{\&gamma;}\right)})---\left(28\right)$

其中，(22)式中为均值，为方差的倒数，式(25)中 $G^{*}\left(x\right|a,b)\&Proportional;\left\{\begin{array}{cc}{e}^{-\frac{b}{2}{(x-a)}^2},& x\&GreaterEqual;0\\ 0,& x<0\end{array}\right.;$

6)通过变分期望最大化定理的变分最大化来推导各个参数的近似后验概率函数的分布参数的更新等式，

近似后验概率密度函数的分布参数的更新等式为：

$\&dtri;{\tilde{L}}_p(i,j)=\underset{c=1}{\overset{C}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{\mathrm{ijc}}^{\left(L\right)}{\&lang;v_c\&rang;}_q+\underset{m,n}{\mathrm{\&Sigma;}}{\&lang;\mathrm{\&gamma;}\&rang;}_q{\&lang;K{(m,n)}^2\&rang;}_q---\left(29\right)$

$\begin{array}{c}\&dtri;{\hat{L}}_p(i,j)\&dtri;{\tilde{L}}_p(i,j)=\\ \underset{m,n}{\mathrm{\&Sigma;}}{\&lang;K(m,n)\mathrm{\&gamma;}(\&dtri;P(i+m,j+n)-\underset{m^{\&prime;}\&NotEqual;m,n^{\&prime;}\&NotEqual;n}{\mathrm{\&Sigma;}}K(m^{\&prime;},n^{\&prime;}){\&dtri;L}_p(i+m-m^{\&prime;},j+n-n^{\&prime;}))\&rang;}_q\end{array}---\left(30\right)$

${\stackrel{-}{a}}_c^{\left(L\right)}=a^{\left(L\right)}+\frac{1}{2}\underset{i,j}{\mathrm{\&Sigma;}}{\mathrm{\&lambda;}}_{\mathrm{ijc}}^{\left(L\right)}{\&lang;{\left({\&dtri;L}_p(i,j)\right)}^2\&rang;}_q---\left(31\right)$

${\stackrel{-}{b}}_c^{\left(L\right)}=b^{\left(L\right)}+\frac{1}{2}\underset{i,j}{\mathrm{\&Sigma;}}{\mathrm{\&lambda;}}_{\mathrm{ijc}}^{\left(L\right)}---\left(32\right)$

${\stackrel{-}{c}}_c^{\left(L\right)}=c_c^{\left(L\right)}+\underset{i,j}{\mathrm{\&Sigma;}}{\mathrm{\&lambda;}}_{\mathrm{ijc}}^{\left(L\right)}---\left(33\right)$

${\mathrm{\&lambda;}}_{\mathrm{ijc}}^{\left(L\right)}\&Proportional;\exp \left({\&lang;\log \left({\mathrm{\&pi;}}_c^{\left(L\right)}G(\&dtri;L_p(i,j)|0,v_c)\right)\&rang;}_q\right)---\left(34\right)$

$\tilde{K}(m,n)=\underset{i,j}{\mathrm{\&Sigma;}}{\&lang;\mathrm{\&gamma;}\&rang;}_q{\&lang;\&dtri;P{(i-m,j-n)}^2\&rang;}_q---\left(35\right)$

$\begin{array}{c}\hat{K}(m,n)\tilde{K}(m,n)=\\ \underset{i,j}{\mathrm{\&Sigma;}}{\&lang;\mathrm{\&gamma;}(\&dtri;P(i,j)-\underset{n^{\&prime;}\&NotEqual;j-n}{\underset{m^{\&prime;}\&NotEqual;i-m,}{\mathrm{\&Sigma;}}}K(m^{\&prime;},n^{\&prime;}){\&dtri;L}_p(i-m^{\&prime;},j-n^{\&prime;})){\&dtri;L}_p(i-m,j-n)\&rang;}_q\\ -\underset{d=1}{\overset{D}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{\mathrm{mnd}}^{\left(L\right)}{\&lang;{\mathrm{\&lambda;}}_d\&rang;}_q\end{array}---\left(36\right)$

${\stackrel{-}{a}}_d^{\left(K\right)}=a^{\left(K\right)}+\frac{1}{2}\underset{m,n}{\mathrm{\&Sigma;}}{\mathrm{\&lambda;}}_{\mathrm{mnd}}^{\left(K\right)}{\&lang;K(m,n)\&rang;}_q---\left(37\right)$

${\stackrel{-}{b}}_d^{\left(K\right)}=b^{\left(K\right)}+\frac{1}{2}\underset{m,n}{\mathrm{\&Sigma;}}{\mathrm{\&lambda;}}_{\mathrm{mnd}}^{\left(K\right)}---\left(38\right)$

${\mathrm{\&lambda;}}_{\mathrm{mnd}}^{\left(K\right)}\&Proportional;\exp \left({\&lang;\log \left({\mathrm{\&pi;}}_d^{\left(K\right)}G(K(m,n)|0,v_d)\right)\&rang;}_q\right)---\left(39\right)$

${\stackrel{-}{a}}^{\left(\mathrm{\&gamma;}\right)}=a^{\left(\mathrm{\&gamma;}\right)}+\frac{1}{2}\underset{i,j}{\mathrm{\&Sigma;}}{\&lang;{(P(i,j)-\underset{m,n}{\mathrm{\&Sigma;}}K(m,n)\&dtri;P(i-m,j-n))}^2\&rang;}_q---\left(40\right)$

其中，I×J为大小，-M_K＜i≤I+M_K,-N_K＜j≤J+N_K,-M_K＜m≤M_K,-N_K＜n≤N_K，

根据式(29)到式(41)所述的更新等式通过变分贝叶斯期望最大化不断迭代，直达算法收敛，算法的收敛标准为代价函数从一次迭代到下一次迭代的增量的绝对值小于一个预定的阈值；

7)更新等式收敛到局部最小值，对模糊图像L_p进行金字塔分解，得到S层由分辨率低到高的图像金字塔，

其中φ为点扩散函数的最大尺寸，

令第S层为最高层，即为图像分辨率最高的一层，为第一层选择模型初始化值，通过带入参数更新等式中，反复迭代计算代价函数，直到代价函数收敛于一个设定的阈值，得到各个参数的最优解，估计出金字塔第一层的PSF，记为K1，以及当前层清晰图像的梯度值然后利用双阈值法处理K1，此时阈值t_low和t_high分别为0.01和0.05，进一步的，将利用双线性插值法放大到金字塔第二层大小得到第一层清晰图像梯度放大值同时也将点扩散函数放大到第二层点扩展函数大小得到一个新的点扩展函数K2'，同样的再次使用双阈值法对新的点扩展函数K2'进行处理，此时阈值t_low和t_high分别为0.01和0.04，将和K2'作为第二层迭代的初始值，再次通过参数更新等式的迭代估计得到金字塔第二层的点扩散函数K2和当前层清晰图像的梯度值以此类推，最终得到金字塔S层PSF K，

双阈值方法特征为：

本方法使用两个阈值t_low和t_high来抑制PSF中的噪声，其中t_high＞t_low，定义两个掩膜M_low和M_high

$M_{\mathrm{low}}\left(i\right)=\left\{\begin{array}{cc}1& k\left(i\right){\&GreaterEqual;t}_{\mathrm{low}}{k}_{\max }\\ 0& k\left(i\right)<t_{\mathrm{low}}{k}_{\max }\end{array}\right.$

$M_{\mathrm{high}}\left(i\right)=\left\{\begin{array}{cc}1& k\left(i\right){\&GreaterEqual;t}_{\mathrm{high}}{k}_{\max }\\ 0& k\left(i\right)<t_{\mathrm{high}}{k}_{\max }\end{array}\right.$

其中k_max表示PSF中的最大值，在得到两个掩膜后，以M_high中值为1的元素为中心，观察其8邻域中的点在M_low中的值，若M_low中的值为1则相应的使这个点在M_high中的值也为1，否则这个点在M_high中的值还是为0，本方法认为判断完所有M_high中值为1的点为一次迭代，不断迭代直到M_high中点的值没有变化；

8)利用估计得到的点扩散函数K，对模糊图像B用Richardson-Lucy算法进行解卷积得到清晰图像L(见附图二)；

9)由于得到的清晰图像L有明显的振铃效应，为了减小振铃效应的影响，提出了一种自适应模糊均值滤波器去振铃效应方法，利用模糊图像B得到图像的细节区，振铃区和平坦区，其中细节区包含图像细节，振铃区包含图像振铃效应所在区域，平坦区包含基本无细节存在的区域，具体步骤如下：

①利用canny边缘检测算子检测模糊图像B的边缘得到二值边缘图像I_bw，

②将二值边缘图像I_bw进行膨胀处理，得到二值边缘图像I_bw'，进一步的对二值边缘图像I_bw'进行连通域检测，删除面积小于一定阈值的连通域，使得图像细节丰富的区域变为一个连续的二值块，得到的二值图像M_d'比实际图像的细节区要宽，对图像再次进行形态学腐蚀处理，使得到实际图像的细节区域掩膜M_d，

③振铃效应通常出现在细节区域附近的平坦区域，基于这个特性，利用细节区掩膜M_d进行当前像素是否为振铃区的判定，

从细节区掩膜M_d左上角的像素开始，判断当前像素M_d(i,j)是否属于细节区，若当前像素M_d(i,j)＝0，则当前像素不属于细节区，进一步判断当前像素是否为振铃区，由于振铃效应通常出现在细节区域附近的平坦区域，因此若在离当前像素距离为T(T由预先设置)范围内的像素中有任何一个像素值为1表示此像素在细节区附近，则可判定此像素属于振铃区，否则此像素属于平坦区，判定为振铃区的像素组成振铃区掩膜M_r，判定为平坦区的像素组成平坦区掩膜M_p；

10)利用振铃区、平坦区、细节区掩膜，针对不同区域采用不同程度的模糊均值滤波器清晰图像L进行滤波处理，使得保持图像边缘和细节信息的同时，一定程度上降低振铃效应对图像去模糊的影响，得到最终模糊图像复原结果(见附图4)，

模糊加权均值滤波器的表达式为：

$y_i=\frac{\underset{j=1}{\overset{N_f}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}{x}_j}{\underset{j=1}{\overset{N_f}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}}---\left(42\right)$

w_ij＝exp(-(x_i-x_j)²/k)   (43)

其中N_f为当前滤波窗口内像素的总个数，x_i为当前像素，x_j为其领域像素的灰度值，y_i为当前像素滤波后的灰度值，w_ij为滤波器权重，k为尺度扩散参数，当k较小时，滤波器平滑性能较弱，当k较大时，滤波器平滑性能较强，将细节区尺度参数值选为k_d，振铃区为k_r，平坦区为k_p，k_d≤k_p＜k_r，使得滤波后所得图像在保持图像边缘和细节信息的同时，一定程度上降低振铃效应。

Claims (1)

1.一种基于自然图像统计特性的盲去模糊方法，其特征在于包括如下步骤：

1)将相机抖动模糊图像表达为清晰图像与点扩散函数的卷积再加上噪声的形式，

$B=K\&CircleTimes;L+N---\left(1\right)$

式(1)中，B表示模糊图像，K表示模糊函数，L表示未模糊清晰图像，N表示噪声，表示卷积操作，其中只有B已知，

图像梯度系数的直方图在零点处有大的峰值，而在远离零点处具有长的尾，在梯度域内对模糊图像进行复原，采用零均值高斯混合模型对未模糊清晰图像梯度进行建模，

根据卷积运算的性质，自然图像降质模型在梯度域内表示为

$\&dtri;B=K\&CircleTimes;\&dtri;L+\&dtri;N---\left(2\right)$

其中和分别表示模糊图像、清晰图像和噪声的梯度；

2)假定点扩散函数具有空间移不变性，即模糊图像的全图受到同一个点扩散函数的影响，选取图像的部分区域来代替全图进行点扩散函数的估计，选取的部分区域记为P，梯度域内的图像降质模型变为：

$\&dtri;P=K\&CircleTimes;\&dtri;L_p+\&dtri;N---\left(3\right)$

式中L_p表示在模糊图像中选取的区域所对应的清晰图像；

3)根据相机抖动模糊图像点扩展函数统计特性，点扩展函数概率密度分布类似于指数分布，采用混合指数分布对点扩展函数进行建模；

4)根据贝叶斯原理，得到清晰图像梯度点扩散函数K的后验概率：

$p(K,{\&dtri;L}_p|\&dtri;P)=\frac{p(\&dtri;P|K,\&dtri;L_p)p(K,{\&dtri;L}_p)}{p(\&dtri;P)}---\left(4\right)$

采用分布来逼近真实的后验概率分布通过近似分布q和真实后验概率分布之间的Kullback-Leibler散度的最小化来实现近似分布q的优化；

5)根据散度的计算定义贝叶斯变分方法的代价函数，然后根据代价函数和模糊图像梯度的先验分布和点扩散函数的先验分布，用变分贝叶斯期望最大化定理求出各个参数的近似后验概率分布；

6)通过变分期望最大化定理的变分最大化来推导各个参数的近似后验概率函数的分布参数的更新等式；

7)对模糊图像B进行金字塔分解，得到S层由分辨率低到高的图像金字塔，令第S层为最高层，即为图像分辨率最高的一层，为第一层选择模型初始化值，通过带入分布参数的更新等式中，反复迭代计算代价函数，直到代价函数收敛于一个设定的阈值，得到各个参数的最优解，由此估计出金字塔第1层的点扩散函数，记为K1，以及当前层清晰图像的梯度值接着采用双阈值方法对估计得到点扩散函数K1进行处理，然后将第一层清晰图像的梯度值利用双线性插值法放大到金字塔第二层大小得到第一层清晰图像梯度放大值同时也将点扩散函数放大到第二层点扩展函数大小得到一个新的点扩展函数K2'，同样的再次使用双阈值法对新的点扩展函数K2'进行处理，将第一层清晰图像梯度放大值和经过双阈值法处理的点扩展函数K2'作为第二层迭代的初始值，再次通过分布参数的更新等式的迭代估计得到金字塔第二层的点扩散函数K2和当前层清晰图像的梯度值以此类推，最终得到金字塔S层点扩展函数K，

双阈值方法特征为：

使用两个阈值t_low和t_high来抑制点扩散函数K中的噪声，其中t_high＞t_low，定义两个掩膜M_low和M_high

$M_{\mathrm{low}}\left(i\right)=\left\{\begin{array}{cc}1& k\left(i\right){\&GreaterEqual;t}_{\mathrm{low}}{k}_{\max }\\ 0& k\left(i\right)<t_{\mathrm{low}}{k}_{\max }\end{array}\right.$

$M_{\mathrm{high}}\left(i\right)=\left\{\begin{array}{cc}1& k\left(i\right){\&GreaterEqual;t}_{\mathrm{high}}{k}_{\max }\\ 0& k\left(i\right)<t_{\mathrm{high}}{k}_{\max }\end{array}\right.$

其中k_max表示PSF中的最大值，在得到两个掩膜后，以掩膜M_high中值为1的元素为中心，观察其8邻域中的点在掩膜M_low中的值，若掩膜M_low中的值为1则相应的使这个点在掩膜M_high中的值也为1，否则这个点在掩膜M_high中的值还是为0，判断完所有掩膜M_high中值为1的点为一次迭代，不断迭代直到掩膜M_high中点的值没有变化；

8)利用点扩散函数K，对模糊图像B用Richardson-Lucy算法进行解卷积得到清晰图像L；

9)利用模糊图像B得到图像的细节区，振铃区和平坦区，其中细节区包含图像细节，振铃区包含图像振铃效应所在区域，平坦区包含基本无细节存在的区域，振铃效应通常出现在细节区域附近的平坦区域，结合形态学图像处理划分得到三个区域掩膜M_r，掩膜M_p，掩膜M_d，分别表示振铃区掩膜，平坦区掩膜，细节区掩膜；

10)利用振铃区掩膜、平坦区掩膜、细节区掩膜，针对不同区域采用不同程度的模糊均值滤波器对图像L进行滤波处理得到最终模糊图像复原结果。