面向图像复原的分层贝叶斯局部高斯混合模型*

张墨华，彭建华

1.国家数字交换系统工程技术研究中心，郑州 450000

2.河南财经政法大学 计算机与信息工程学院，郑州 450000

1 引言

信号处理领域的多种方法[1-7]应用到各种图像复原任务中，退化图像Y（以向量化形式）可以建模为：

其中，Y是观测退化图像；X是潜在的真实图像；H是退化算子；V是加性噪声。不同的退化矩阵H对应着不同的逆问题，例如填充、去模糊或是超分辨等。例如，用于图像超分辨[8-9]，H是下采样矩阵；用于图像去模糊[10]，H是模糊矩阵；用于图像去噪时[11]，H是单位矩阵等。噪声项V可能是恒定方差的经典加性高斯噪声或是与信号相关噪声的更复杂的模型[12]。图像复原问题任务就是从观测图像Y来估计原始清晰图像X。

由于这类逆问题固有的不适定性[13]，通常复原问题在变分或是贝叶斯求解时会在图像上施加某种先验，例如总变分（total variation，TV）[14]、小波分解[15]、图像块的稀疏性[16]等各种流行的图像模型。Buades等人提出图像自相似性先验的方法[17]，开启了基于块的图像复原技术的新时代。通过在贝叶斯框架中引入块先验模型，有些方法致力于去噪问题[18-20]，有些提出通用框架来求解各种逆问题[21-22]，包括去噪、填充、去模糊及放大。

在贝叶斯框架下有两个值得注意的方法：一个是Yu等人提出的基于块的贝叶斯方法[22]，即分段线性估计（piecewise linear estimators，PLE）方法；另一个是Lebrun等人提出的非局部贝叶斯（non-local Bayesian，NLB）算法[23]。PLE是在式（1）下解决图像逆问题的一般框架，而NLB是一种去噪方法（H是单位阵）。这两种方法都是利用从图像块学习的高斯块先验通过迭代求解。PLE是通过高斯混合模型（Gaussian mixture model，GMM）来进行建模，它采用相对较少的分量数目（在实验中用的是19），其参数通过图像块进行学习得到。作者在实验中采用图像子区域（大小为128×128）的方法，可以把PLE看成是半局部的方法。而NLB方法中每个块与一个高斯模型相关，其参数是通过其局部近邻（搜索窗口）中选择的相似块计算得到[23-24]。Zoran和Weiss提出的块对数似然期望（expected patch log likelihood，EPLL）算法[21]采用了类似的方法，但其不是从图像块中迭代地更新高斯混合模型，而是利用大量的分量数目，从大型图像数据库中学习得到。Wang和Morel在文献[25]中指出，在去噪问题上，具有较少的分量的方法，性能要高于拥有较多固定分量数量的方法。

上述的这些复原方法都是基于高斯块先验的贝叶斯框架。相比于有限数量的高斯混合模型方法[21-22]，依赖局部先验的方法已经被证明对于图像去噪任务更为准确[23，25]。特别是，由于局部模型估计的原因，非局部贝叶斯方法要比分段线性估计方法在处理图像去噪问题上更为优秀。另一方面，分段线性估计方法在其他应用场景，如图像填充、去模糊及图像放大等应用中，取得了优秀的结果。

因此，本文聚焦于如何将局部块模型应用到更为通用的逆问题求解，而不仅仅是去噪问题。但是这一应用的主要难点在于模型的估计问题，特别是当退化图像中有较高比例缺失像素情况下，估计问题是严重不适定问题。

本文提出根据高斯先验来建模图像块，其参数（均值μ和协方差矩阵Σ）从局部相似块中进行估计。为了处理这一问题，利用分层贝叶斯模型[26]，将模型参数的先验知识（先验参数的概率分布）考虑进来。贝叶斯统计中，μ和Σ称为超参数，超参数的先验称之为超先验。这种模型应用于图像复原领域，最早由Molina[27]提出，近年来应用到图像分割和去卷积[28]及超声图像的复原[29]中。将分层贝叶斯方法应用于局部高斯混合模型的基于块的图像复原建模中，通过超先验来减少由于退化图像中块信息丢失而带来的算法不稳定性。基于L2范数在局部窗口中完成相似块的聚类，利用累加平方图及快速傅里叶变换的数值优化方法，加快相似性度量的计算时间。在从块到图像的聚合过程中，采用基于马式距离的高斯分布相似度的聚合权重，结合图像上的空间域高斯相似度，从而更好地拟合自然图像的统计特性。

本文的主要贡献有以下两点：第一，提出将分层贝叶斯模型引入到图像块的高斯局部先验建模中并用于求解图像复原问题，使得对于缺失像素的逆问题求解更加稳定；第二，利用基于马式距离的高斯分布相似度的聚合权重方法，将局部高斯模型扩展到高斯混合模型框架下。

2 全局高斯混合模型

PLE和EPLL这类方法属于全局高斯混合模型[21-22]，它们只是在初始化和计算聚类权重上有细微差别。通常这类方法假定图像中每个块（n维）独立来自于M个有限多元高斯分布{N(μm,Σm)}1＜m＜M中的一个，该分布的均值为μm，方差为Σm。每个块xi都是从这些有限数量的高斯中的某一个独立地生成，其概率为：

为了最大化所有块的概率分布，需要得到整个图像的有限高斯分布，这可以通过在初始化高斯分布之后，迭代地进行下述步骤得到：

（1）根据估计的(μm,Σm)1＜m＜M，由式（2）度量相似性，将块进行聚类。产生每个块的高斯概率是由{N(μm,Σm)}1＜m＜M确定，这引入了一个基于模型的框架，将恢复的块xi分配给M个估计高斯分布中之一。

（2）更新{N(μm,Σm)}1＜m＜M。每个均值向量和协方差矩阵的估计值，是基于聚类块进行更新的。这些估计可以通过每个聚类中的样本均值和样本协方差得到。

（3）每个块的恢复基于其所属高斯分布，利用维纳滤波而得到。

PLE和EPLL这种全局高斯混合模型方法将多个块从图像的不同部分分配到一个聚类。这种全局聚类限制了利用图像中邻近块的相干性。一些成功的图像去噪算法，例如依赖于平均像素的非局部均值[3，16]，通过设置与像素或块之间的几何距离成反比的平均权重来考虑邻近块的相干性。BM3D（blockmatching and 3D filtering）[2]、NCSR（nonlocally centralized sparse representation）[30]、FNCSR（fast nonlocally centralized sparse representation）[4]等算法也是采用约束在一个有限大小的窗口中分组相似的块，协同进行去噪。文献[31]指出，相比于全局窗口大小，小窗口中的块更可能从高斯分布导出，意味着更易遵循高斯分布。

3 分层贝叶斯的局部高斯混合模型

为了对图像应用空间约束并利用邻域块的相干性，假设邻域中的相似块可以从具有特定均值和协方差的单个多变量高斯概率分布导出，让每个块与一个高斯模型相关，其参数是通过其局部近邻（搜索窗口）中选择的相似块计算得到。非局部贝叶斯算法就是采用这一思路，应用于图像去噪问题[23]。本文关注于将局部块模型如何应用到更为通用的逆问题求解，而不仅仅是去噪问题。这一应用的主要难点在于模型的估计问题，特别是当退化图像中有较高比例缺失像素情况下，估计问题是严重不适定问题[32]。

为此本文提出分层贝叶斯的局部高斯混合模型（hierarchical Bayesian local Gaussian mixture model，HBLGMM），图1给出所提出模型的总体思路。

Fig.1 Hierarchical Bayesian local Gaussian mixture model图1 分层贝叶斯的局部高斯混合模型

如图1所示，在局部窗口中利用L2范数计算相似块，来充分利用邻域块的相干性，这里假定局部相似图像块满足均值μ和协方差矩阵Σ的高斯先验分布。对均值和协方差参数引入贝叶斯超先验（hyper prior），对图像块均值μ和协方差Σ引入Norm-Wishart先验（参数为μ0和Σ0），利用联合最大后验估计公式用来估计图像块和参数采用交替最小化更新的方式，使高斯统计的局部估计更加稳定；在从块到图像的聚合过程，对于重叠区域利用块对于不同局部高斯分布的相似度，计算基于马式距离（Mahalanobis distance）的高斯权值，从而从局部高斯分布扩展到高斯混合模型。这样可以利用局部模型估计的准确性来处理通用的复原问题，特别是处理丢失像素值（填充问题）。

3.1 块退化模型

观测图像Y分解为I个相互重叠的块{yi}i=1,2,…,I，块大小为。每个块yi∈Rn×1可以看作是随机变量Yi的实现，Yi可以表示为：

其中，Hi∈Rn×n是退化算子；Xi∈Rn×1是需要估计的原始块；Vi∈Rn×1是加性噪声项，通常以高斯分布进行建模。在给定Xi情况下Yi的条件分布为：

假定在每个块的高斯先验模型为p(Xi|μ,Σ)～N(μ,Σ)，其中均值μ和协方差矩阵Σ是未知的。为了简化计算，利用精度矩阵Λ=Σ-1来表示。本文假定(μ,Λ)服从Normal-Wishart先验。Normal-Wishart先验是未知均值和协方差矩阵的多元高斯分布的共轭先验。(μ,Λ)的分布为：

其中，参数包括μ0和Σ0，尺度参数λ＞0和自由度ν＞n-1。

假定已经有观测相似块{Yi}i=1,2,…,M，要恢复的复原块为{Xi}i=1,2,…,M。假定未知{Xi}是独立的且服从相同的高斯模型，计算联合最大后验：

式（6）乘积第一项为式（4）的噪声模型，第二项是图像块集合{Xi}上的高斯先验，第三项为式（5）的超先验。其中第一项在p({Yi}|{Xi},μ,Λ)上忽略了μ和Λ依赖，这是因为这些参数都是由{Xi}确定的。

3.2 优化过程

函数f是分别在变量({Xi},μ)和Λ上的双凸函数。对于给定的超参数μ0和Σ0，可以采用交替凸最小化机制来最小化该函数。即在每轮迭代中，首先固定Λ，相对于({Xi},μ)求最小化，然后进行交替，固定({Xi},μ)，相对于Λ求最小化。

通过f对每个变量求微分，可以得到最优化方程。假定Λ是固定值并且协方差不依赖于{Xi}。函数上是凸的，并且可以得到最小值：

假定变量({Xi},μ)是固定的。函数Λ↦f({Xi},μ,Λ)在上是凸的，并且可以得到最小值：

从式（8）可以看出μ的MAP估计量是两项值的加权平均：先验μ0和利用从相似恢复块得到的均值估计量。参数γ控制着对于先验μ0的置信度。类似地，Λ的MAP估计量也是由三项构成：先验Σ0和μ的协方差，以及来自块Xi所估计的协方差矩阵。

根据上述分析，可以得到f的交替最小化过程，如算法1所示。

算法1f交替最小化

算法1中，迭代次数设置为Iters1。当l→∞时，函数在每个步骤中会逐渐减少，函数f的连续性意味着该序列下端有界，因此会收敛。函数的收敛性意味着序列是有界的。由此可知，它具有至少一个累积点Λ*)，并且存在严格增加的整数序列(lk)，使得收敛到。因此通过交替最小化方案生成的序列至少有一个聚点。实验中也发现算法经过几轮迭代后就会收敛。

3.3 相似块搜索

在当前样例块的搜索窗口中找到与样例块L2范数距离小于指定阈值的相似块。阈值的计算是通过容忍参数ϵ乘以最邻近距离来确定。在所有实验中，搜索窗口设计为32×32（块大小设置为8×8），ϵ设置为1.5。在前一轮迭代的结果图像上进行相似块的比较，块间相似度的计算是基于图像块的L2范数进行的。

受文献[33-34]启发，采用数值优化方法来加速欧式距离的计算。令R是块抽取算子，R(p)∈Rn表示以图像中以像素p=(p1,p2)为中心抽取的块。当且仅当位移t=(i,j)在半径为a（块的半径）的圆范围内（记作），则像素(p1+i,p2+j)在块R(p)中。给定两个块，其L2范数距离为：

式（12）的最后一项可以看作是两个块的卷积。因此可以借助快速傅里叶变换在频域进行乘法操作来计算。利用上述算法，对一幅图像每个像素只需要处理一次。通过SSI，可以以常量时间得到图像任意矩阵中每个像素的平方累加和。

3.4 块聚合过程

把估计块放回其原始位置可以重建整个图像。为了改善恢复性能，如BM3D这类算法对重叠的评估块利用加权平均的方法来重构恢复图像[2]。本文采用高斯核的方法，这种方法已经在双边滤波和非局部均值[4]等方法中应用来获取平均权重，这些权重的非归一化形式用来表示块p与q的相似度，可以表示为：

其中，ξ用于设定可能权值合适的尺度，d是度量块p和块q的距离。在非局部均值方法中，权值基于中心像素为p和q的块距离来计算。归一化通过除以像素的权值之和实现。受高斯核用于平均权值的启发，本文提出基于式（14）中的高斯核的块的聚合权重，其中d采用马式距离度量，该距离已被用于基于GMM的聚类方法，作为测量聚类的高斯分布相似性的标准。因此第r个高斯分布导出的块的权重由式（14）计算得到。

这些权重分配给整个块而不是中心像素。与其他基于高斯核的方法（如非局部均值法）相比，本文方法权重是基于块与模型的相似度来测量的，而不是块与块相似度。通过使用马氏距离，与估计高斯分布越相似的块就有更高的权值。

虽然使用L2范数距离的相似性度量，聚类的方法似乎没有使得聚类块与高斯分布较高概率拟合，但是通过使用基于与高斯分布相似度的聚合权重，可以得到图像上的空间域高斯相似度。因此，随着迭代的进行，具有更高相似性的高斯分布的块通过L2范数聚类分组在一起，这表明使用这种权重进行平均的重要性。

3.5 初始化

由于通过迭代过程解决非凸问题，因此良好的初始化是至关重要的。本文利用Yu等人提出的方法[22]，该方法是通过从不同方向边缘的合成图像以及DCT基来表示各向同性的模式，学习K个分量的GMM协方差矩阵（PLE中设置为19），从而完成PLE算法的初始化过程。文献[22]中指出在字典学习中，最突出的原子代表局部边缘，这些边缘可用于表示和恢复轮廓，因此这种初始化有助于恢复受损的块。PLE算法第一轮的输出即是初始的理想图像Co。

3.6 图像复原算法

分层贝叶斯的局部高斯混合模型完整图像复原算法见算法2，迭代次数设置为Iters2。算法需要两个阶段：通过算法1最小化f；估计超参数μ0和Σ0。为了估计这些参数，需要依赖于第一阶段估计所有块的聚合计算得到的图像。

算法2HBLGMM图像复原算法

3.7 计算复杂度

算法的计算复杂度主要集中于算法1中计算(Xl,μl)的步骤，矩阵的逆运算利用了Cholesky分解，这样总的时间复杂度为。其中NSP是相似块的数目，n是块中像素数目，NP是要恢复的块的总数目。算法可以对图像进行子区域划分，利用多核架构采用并行方式进行并行计算来加速。

4 实验及分析

本章给出图像去噪和图像填充两个图像复原问题求解结果。通过与一些优秀的方法进行对比进而验证本文提出的模型。在填充问题上考虑两种情况的插值：一种是随机观测的像素；一种是放大问题（可以看成是均匀观测像素的插值问题）。对比方法的代码来自作者公开代码。本文实验用图一部分来自经典的合成图像，如Lena、Barbara等，另一部分选取BSDS数据集中的测试图像[35]。部分实验用图如图2所示，第一行为经典合成图像，第二行为BSDS图像。

采用峰值信噪比（peak signal to noise ratio，PSNR）作为算法的客观度量，PSNR定义为：

其中，x是原始的退化图像；是估计图像；m是图像的像素数目。估计图像的质量越高，会有更高的PSNR值。通过每个实验的10次实验平均值得到最终PSNR值。

4.1 参数选择

算法涉及的主要参数包括计算μ和Σ需要设置Norm-Wishart分布的4个参数γ、ν、先验均值μ0和先验协方差Σ0以及计算高斯权重的参数ξ。

4.1.1 γ 和ν 的确定

Fig.2 Part of evaluation image图2 部分实验用图

从式（7）可见，均值μ的计算涉及到相似块的均值估计和先验均值μ0。参数γ关联于先验μ0的置信度，因此它的值可以在先验准确性的置信和相似块提供的信息间进行折衷。随着相似块数目NSP的增多和块中已知像素数NKP的增加，相似块提供的信息会得到改善。因此，γ值的计算如下：

其中，TH是阈值。对于ν可以采用类似的做法：

式（16）加上n是因为Normal-Wishart先验所要求的ν＞n-1。

通过实验证明上述简单的规则是一种有效的方法。上述方法主要针对填充问题，但是在推广到更普通的其他应用中，需要进一步研究。

4.1.2 μ0 和Σ0 的确定

μ0和Σ0可以通过经典的最大似然估计的方法从Co的相似块集合中进行计算：

4.1.3 ξ 的设置

为了获得良好的复原性能，对式（14）中的ξ进行小幅的调整是必要的。例如在去噪实验中，针对不同的噪声级别设置不同的ξ值：在σ≤40时，ξ=0.015；在σ＞40时，ξ=0.01在填充实验中，设置ξ=0.01。通过ξ的小幅调整，PSNR值取得一定的提升。

4.2 图像去噪

在实验中，原始图像利用标准差分别为30、40、50、60、70的加性高斯噪声进行污染，在去噪实验中，由于没有未知像素要插值，因此H为单位阵。在表1中，本文方法（标记为HBLGMM）与BM3D[2]、LSSC（learned simultaneous sparse coding）[16]、NLB[23]、EPLL[21]等优秀的去噪算法进行了对比。表中数值为每个噪声标准差下所有测试图像的平均PSNR值（单位：dB），最好的结果用加粗字体表示。

在所有噪声级别下，本文方法比没有应用局部约束的全局GMM的方法（EPLL），以及基于GMM的方法（NLB），从平均来看都要优秀。与其他全局GMM方法类似，组中所有观测块减掉就是每个聚类的均值向量，完成MAP的估计之后再加上，这样做的主要目的是加速算法的执行。之前基于GMM的方法如PLE、EPLL用于图像去噪的性能总体低于基于稀疏的方法。本文方法的去噪结果从平均结果来看，在噪声标准差于30～60范围比基于稀疏的方法效果更好。LSSC方法在高噪声下PSNR占优，与其内部字典训练的优势有关，但值得一提的是LSSC方法的在线字典训练时间过长，实际可用性较差。对比的示例图像如图3所示，对Barbara图加入标准差为30的噪声，图中第一行为本文方法与其他优秀方法去噪结果对比，第二行为还原图像与原始图像的差值图，第三行为局部区域的放大图对比。从放大区域清楚可见，本文方法在裤子纹理的还原上更加自然。

Table 1 Comparison of denoising results表1 去噪结果对比

图4为高噪声下BSDS集合中测试图像去噪结果对比（噪声标准差为70）。从对比图中可见，BM3D和NLB方法在这一数据集中图像观感相对较差，人脸还原不够自然；BM3D方法在帽子区域有较多波纹，EPLL方法则也在帽子区域不够平滑，LSSC尽管PSNR值更优，但人物脸部、背景及帽子区域有大量的伪像，整体观感较低，而本文方法无论在PSNR还是视觉观感上都较为理想。

4.3 图像填充

图像填充实验中设计两种应用场景：一种是随机遮蔽算子（用于随机填充）；一种是均匀网格下的下采样遮蔽算子（用于图像放大）。

在随机填充实验中，设置20%、30%、40%、50%、60%、70%、80%等像素缺失率的随机遮蔽算子应用于原始图像中。所提出方法的填充性能与MCA（morphological component analysis）[36]、FOE（field of exports）[37]、KR（kernel regression）[38]、BP（beta process）[39]、PLE[22]等方法进行了对比，所有方法图像块的大小均使用的是8×8的大小。用于计算累积权重的式（14）中参数ξ设置为0.01。

Fig.3 Comparison of visual quality of denoising results(σ=30)图3 去噪视觉质量对比(σ=30)

Fig.4 Comparison of visual quality of denoising results(σ=70)图4 去噪视觉质量对比(σ=70)

表2给出本文方法与其他的图像填充方法对比的结果，表中数值为不同数据丢失率下所有测试图像的平均PSNR值。从表2可见，本文方法在所有的丢失率下比对比方法都要优秀。

图5给出局部Barbara图像在随机丢失80%像素情况下还原对比示例，重点是比较图像的平滑和纹理区域的填充结果。从图可见，本文方法比其他方法无论在纹理还是平滑区域，效果更为理想。

图像填充的一种特殊情况是图像分辨率放大问题，可以看成是均匀采样图像的插值问题。然而，图像放大问题要比随机观测像素的填充更具有挑战性。因为这种规则的采样，目前提出的很多算法未能很好地复原图像中潜在纹理。如文献[40]中所讨论的，随机采样要比均匀采样可以得到更好的结果。本文算法对于由均匀采样带来的复原错误更为鲁棒。

Table 2 Comparison of random inpainting results表2 随机填充结果对比

Fig.5 Comparison of visual quality of random inpainting(data lost ratio=80%)图5 80%丢失率下随机填充视觉质量对比

实验中对真实的高分辨图像（high resolution，HR）采用2倍的下采样得到低分辨图像（low resolution，LR），对LR上采样后，再施加均匀采样遮罩（uniform sample masking），得到退化后的图像用于图像填充复原，本文方法与EPLL[21]、Patch Nonlocal[40]、NEDI（new edge-directed interpolation）[41]、NARM（nonlocal autoregressive modeling）[42]、Bicubic等方法进行了对比。

图6给出对比结果，从图可见本文方法比其他方法具有更为清晰的重建结果，能够较好地找到真实的纹理。

4.4 分析

在利用高斯混合模型先验建模下，无法使用一组固定的聚类来准确表示图像中很少出现的块，例如边缘或纹理。像PLE这样的方法实际上是半局部的（128×128区域），并不能解决这个问题。EPLL则是通过从200万个自然图像块中学习得到200个分量的高斯混合模型[21]，这要比PLE方法有更多的混合分量。如文献[25]所指出的，EPLL的策略在去噪任务中要比PLE方法更加低效。在所有上面的实验中，无论在去噪上，还是在填充和放大问题上HBLGMM都要比EPLL和PLE方法更为优秀。从视觉质量来看，图像的细节得到了较好的重建，图像更加清晰。本文方法从局部窗口中的相似块执行估计，通过考虑局部邻域来强化自相似性的假设，这种策略对于模型进行邻域限制，因此使得估计更加稳健。块的相似性是通过高斯分布的相似性进行度量的，聚合权重考虑了图像块与所估计高斯分布的相似性，结合图像空间域的相似性，从而会将更为相似的高斯块通过L2范数度量聚类汇集到一起。受益于聚类中对于每个块使用不同的权重，这在改善图像复原质量上扮演重要的角色。

Fig.6 Comparison of image zooming result图6 图像放大结果对比图

5 结束语

本文提出将分层贝叶斯模型引入到图像块的高斯局部先验建模中并用于求解图像复原问题，使得对于缺失像素的逆问题求解更加稳定，并从局部窗口中进行相似块估计，利用基于马式距离的高斯分布相似度的聚合权重方法，将局部高斯模型扩展到高斯混合模型框架。利用上述方法，图像复原问题可以在统一的框架中进行，将提出的模型应用于图像去噪、图像填充等问题中，取得了良好的效果。对于图像块尝试使用其他的多元分布，比如多元拉普拉斯，以及采用不同的几何距离度量方法将是未来研究的课题。