面向图像先验建模的可扩展高斯混合模型

张墨华,彭建华

(1.国家数字交换系统工程技术研究中心,郑州 450002;2.河南财经政法大学 计算机与信息工程学院,郑州 450002)

0 概述

图像去噪在数学上属于不适定问题,需要使用图像先验进行正则化,从而生成有意义的解[1]。自然图像的一些自身属性可以作为有用的图像先验用于图像去噪,如稀疏性、多尺度自相似与非局部自相似性等。文献[2-4]采用图像的稀疏性作为先验信息,应用小波与曲波等稀疏变换,或是通过学习图像块得到完备字典以及图像块的稀疏表示。文献[5]采用图像像素在空间域和强度域中表现出的相似性作为先验信息,通过相似块进行协同处理。文献[6]利用同一幅图像中相同尺度或不同尺度的相似子块进行去噪。文献[7-8]将图像的非局部自相似性作为先验信息,对相隔较远的相似块进行非局部处理。文献[9-11]则是使用多种图像先验组合的方法进行图像去噪。

除了使用预定义的图像先验,文献[12]还提出了从自然图像中学习先验的方法。生成式图像先验学习方法通常从一组干净的外部图像中学习先验模型,并将其应用于给定的退化图像[13-14],或者从给定的退化图像中学习先验[15]。近年来,判别式图像先验学习方法也开始得到广泛应用。该方法根据清晰-噪声图像对学习去噪模型,可进一步细分为基于深度学习的方法[16]、基于随机场的方法[17]与基于反应扩散的方法[18]。

在上述方法中,文献[12]提出的使用高斯混合模型(Gaussian Mixture Model,GMM)先验进行自然图像块建模的EPLL方法取得了较大成功。这一方法的主要思想是从干净图像块中学习先验,并将其用于求解图像复原问题。但EPLL方法的混合分量数目固定(文献[12]中为200),学习训练完毕后分量不易进行扩展。针对这个问题,文献[19]提出非参贝叶斯方法,为结构化数据(如文本文档、时间序列和图像)的无监督建模提供了灵活框架,其模型的复杂度可以根据数据规模增减进行自适应变化。狄利克雷过程(Dirichlet Process,DP)是最流行的非参贝叶斯方法,其采样结果是一个分布,而不是变量,因此DP混合模型可以看作无穷多个混合分量的混合模型[20]。作为分层贝叶斯框架下的非参先验,DP可以根据具体参数值对该模型生成的数据进行划分,具有自然解释特性,分量的数量是随机的,可以随着观测到的新数据增长。但DP的学习是一个复杂的最优化问题,难以精确求解,只能进行近似求解,常用的方法之一是变分推理。变分推理中基于传统均值场的坐标上升算法按照先后顺序对数据集中所有的局部变分参数与全局变分参数进行更新,因此在大数据集上效率较低,同时由于更新缓慢,容易陷入局部最优解。

针对上述问题，本文构建一种基于狄利克雷过程的可扩展高斯混合模型，从清晰图像数据库中学习外部通用先验，使模型复杂度根据数据规模自适应变化。同时，为改善模型的推理过程，提出一种基于批次更新方式的可扩展变分算法，将数据集分为若干批次，在每个批次访问中对图像的局部变分参数与全局变分参数进行更新，从而求解图像去噪中所有隐变量的变分后验分布，实现通用先验学习。

1 理论定义

定义1(狄利克雷过程) 令G0是集合Θ上的分布,α是正值实数。对于Θ上的有限划分A1,A2,…,Ar,则向量G(A1,A2,…,Ar)是随机的。如果有:

G(A1,A2,…,Ar)～Ddir(αG0(A1),αG0(A2),…,αG0(Ar)

(1)

则称G是基分布G0和集中参数α上的狄利克雷过程(DP)所生成的分布,记作:G～DDP(α,G0)。

假定从G中独立生成N个随机变量ηn,则有:

ηn～G,n∈{1,2,…,N}

(2)

联合分布{η1,η2,…,ηn}服从波利亚坛子模型(Polya’s urn scheme)[21]。

定理1令G～DDP(α,G0),对于任何集合A,有E[G(A)]=G0(A)

Var[G(A)]=Go(A)(1-Go(A))/(α+1)

从定理1可知,随着α值增大,方差逐步减小。这样DP将更多质量集中于均值附近。α值越小,G就更为离散。证明可以参见文献[22]。

定义2(狄利克雷过程混合模型) 如果ηn是第n个观测变量分布F的参数,即有F(x|ηn),则DPMM可以看作为无穷多个混合分量的混合模型,即有:

G|{α,G0}～DDP(α,G0)

ηn|G～G

xn|ηn～F(x|ηn)

(3)

在DP混合模型中,DP作为分层贝叶斯框架下的非参先验,可以根据具体参数值对该模型生成的数据进行划分,因此,DP混合模型拥有自然解释特性,其中分量的数量(划分中原子的数目)是随机的,并且能够随着观测到的新数据进行增长。

文献[23]根据折棍子构造的过程,对DP给出了更为清晰的特征表述。考虑2个独立随机变量无限集合,vi～Bbeta(1,α),ηi～G0,i={1,2,…},它们的折棍子过程表示如下:

(4)

由式(4)可见,G是离散的。G的支持集包括可数无限原子集合,这些原子集合从G0独立生成。混合比例βi(v)通过将单元长度的“棍子”连续折断成无限数目的分段生成,每分段的大小与棍子剩余部分成正比,从Bbeta(1,α)分布独立生成。

在DP混合中,向量βi(v)是一个包含比例的无限向量,{η1,η2,…}是表示混合分量的原子。令zn表示数据点xn所关联的混合分量的分派变量。数据生成过程如下:

1)生成vi|α～Bbeta(1,α),i={1,2,…}。

2)生成ηi|G0～G0,i={1,2,…}。

3)对于第n个数据点:

(1)生成zn|{v1,v2,…}～Mmult(β(v))。

(2)生成xn|zn,ηn～F(xn|ηzn)。

在图像去噪任务中,每个图像具有其特定的聚类频度wi={wi1,wi2,…,wiK,…},通过施加分层DP先验,聚类频度可以看作一个有限的狄利克雷分布:

[wi1,wi2,…,wiK,wi>K]～

Ddir(αβ1,αβ2,…,αβK,αβ>K)

(5)

从式(5)可知,wi的均值为β,方差由集中参数α所决定,下标>K表示除了前K个分量外其余分量累加之和。本文约定观测数据从指数家族分布生成,DP的基分布是对应的共轭先验。

2 图像生成模型

使用基于DP的高斯混合模型对图像的生成过程进行建模。自然图像可以看作图像块组合,每个块包含D个像素,并且有其对齐基点b,b的取值在D中,即:b～Ccat(1/D,…,1/D)。一幅图像由所有对齐基点下各个块重叠生成,每个对齐基点下,可能有部分块处于部分观测的状态,即部分参与最终图像的生成,部分不参与。

参照文献[12]的方法,将每个图像块去平均化,对齐基点下每个图像块gibn看成由高斯混合分布生成,即每个块表示为均值为0、精度矩阵为S的高斯混合,如下式所示:

(6)

其中,i表示第i幅图像(共N幅),b表示对齐基点(共D个),n表示图像i在对齐基点b中的第n个块(共Nib个)。精度矩阵S的先验为Wishart分布,满足:

S～Wwis(μ,V)

(7)

图像块均值mibn满足高斯分布:

mibn～Nnorm(τ,σ)

(8)

每个图像块可以指派为高斯混合模型中某个分量(聚类)中的一个,用zibn满足分类分布,表示为:

zibn～Ccat(w1,w2,…,wk)

(9)

给定通过对齐基点b生成的、块均值为mibn的块gibn,图像xi的采样可以通过式(10)生成:

(10)

在图像去噪任务中,观测图像yi是干净图像xi的退化版本,是由向xi中加入高斯白噪声后所得。其生成定义如下:

yi～Nnorm(yi|xi,σ2I)

(11)

其中,σ2为噪声方差。

整个图像生成过程如图1(a)所示,其变分参数的推理过程如图1(b)所示,第3节将给出详细过程。

图1 图像生成有向图模型及推理过程Fig.1 Directed graph model of image generation andits inference process

3 图像先验学习

3.1 变分分布

本文通过变分分布的参数推理实现图像先验模型的学习。给定干净观测图像xi,需要学习模型中的各种隐变量,如全局图像变量与局部图像变量。全局图像变量包括折棍子过程的折断比例vk,与每个聚类的精度矩阵Sk;局部图像变量包括聚类频度wi、图像块gibn、块均值mibn与块聚类指派zibn。直接计算后验分布难度较大,因此常用变分推理进行近似求解。变分推理是一种近似似然和后验的确定性方法[24],采用优化方法尝试寻找一个合适变分分布q,它与真实的后验在KL散度上无限接近。

定义待学习的隐变量集合W={g,S,wi,mi,zi},观测图像变量x={x1,x2,…,xN},超参数θ={α,μ,V}。隐变量的后验分布与其近似分布q的KL散度为:

D(q(W)‖p(W|x,θ))=Eq[lbq(W)]-

Eq[lbp(W,x|θ)]+lbp(x|θ)

(12)

式(12)中KL散度的最小化求解可以转换为对数边缘似然下界(ELBO)的最大化求解[25]。如下所示:

lbp(x|θ)≥Eq[lbp(W,x|θ)]-

Eq[lbq(W)]≜L

(13)

将L定义为待优化求解的目标函数,通过参数的学习过程,尽量使L值最大化,从而尽可能减小变分近似后验q与真实后验的差值。采用均值场方法[25],假定后验完全独立,q可以分解为指数家族密度函数的乘积。

(14)

其中,随机度量的随机性通过无限集合{v1,v2,…}和{S1,S2,…}体现,而随机性的生成过程则通过经典的截断折棍子过程表示[21]。假设固定截断位置为T,则混合比例wt(v)在t>T时等于0。T是变分参数,可以自由设置,并非先验模型的一部分。

下面具体给出所有变量的变分后验分布,这些分布参数称为变分参数,是要进行优化求解的对象。

在局部变量因子中,图像块g、图像块均值m、块分量指派z、图像聚类频度w、图像x的变分后验分布分别为:

(15)

在全局变量因子中,折断比例变量v的变分分布为:

(16)

另一个全局变量因子精度矩阵S的变分分布为:

(17)

对于自然图像,所有的对齐基点具有相似的值,设计对齐基点B的变分后验为

q(B)=Ccat(1/D,…,1/D)

(18)

该设计简化了更新操作,同时仍然可以避免不重叠块所造成的伪影情况。

3.2 变分参数的更新

图像通用先验的学习是基于外部干净图像数据库进行的,此时图像x是可观测变量。根据式(12),结合3.1节所述各变分分布,目标函数L可以分解为3个部分:

L≜Limage+Lentropy+Lweight

其中:

Limage≜Eq[lbp(x|g,m,b)+

(19)

(20)

(21)

(22)

(23)

其中,期望Eq[β]通过式(4)进行计算,Cik是对于图像i中所有图像块使用聚类k的情况的统计量,定义为:

(24)

(25)

(26)

从上式中抽取出两个统计量Ck和Fk,定义为:

(27)

(28)

这两个统计量是训练图像集中所有分派给聚类k的数据统计量,将在3.3节中详述。

(29)

(30)

(31)

(32)

其中,Oibn表示在对齐基点b下第n个块中可观测的像素的数目。

(33)

αEq[βk]Tk

(34)

其中,Tk定义为:

(35)

(36)

3.3 通用先验学习算法

传统的坐标上升算法先对数据集中所有局部变分参数进行更新,再对全局变分参数进行更新,因此在大数据集上效率较低。同时,由于全局参数在第一次参与局部参数更新时是初始化值,计算意义不大。此外,传统算法更新较为缓慢,容易陷入局部最优解。针对这些问题,本文利用式(27)、式(28)、式(35)中所得3个统计量的可累加性,采用批次更新方式完成坐标上升算法。将数据集分成若干批次,对所有批次完成访问称为一轮。在每轮批次访问中,首先对批次中每个图像的局部变分参数进行更新,然后进行全局变分参数的更新。

在每个批次中都对式(27)、式(28)、式(35)所得的3个统计量的值进行记录(表示为Ubt),然后累积到全局统计量Uo中,其定义如下:

其中,B表示总批次数。在访问完每个批次bt之后,执行增量更新来累积统计量,从而反映新批次的汇总情况,并删除该批次bt之前的值。

(37)

其中,ISbt是图像数据集IS的批次数据集。综上,通用先验训练过程如算法1所示。

算法1通用先验训练的变分推理算法

输入图像数据集IS

1.Repeat

2.For each batch ISbtin IS do

12.End for

13.Until ELBO收敛

3.4 聚类分量的新增与归并机制

3.4.1 新增机制

在通用先验训练过程中,通过分量新增机制可以增加新的有用分量,帮助规避局部最优解。但即使这些分量得到整个图像数据库的支持,每个批次也无法提供缺失分量所需要的足够样本,因此,整个新增过程需要两轮数据访问。第一轮采集目标数据样本,第二轮生成新的分量,并且使用扩展的模型对每个批次进行更新。具体步骤如下:

2)在第2轮数据访问之前生成新的分量。通过运行有限次数的变分推理,将数据集IS′拟合为包括K′(设置为10)个混合分量的模型。从而将整个模型扩展为K+K′个分量。期间并没有对这些新分量所产生的ELBO的变化进行评估,依赖随后的归并操作来移除不需要的分量。

3)在第2轮中访问每个批次,并且对扩展后的K+K′个混合模型执行局部和全局参数更新。通过对目标数据集X′分析得到统计量U′,并入全局统计量Uo中。数据集X′有两种分派结果:分派给最初的分量(最有可能是k′),或者分派给全新的分量。第2轮结束后,从Uo中去掉U′,使得Uo和全局参数与数据集IS一致。

一次新增多个分量有助于规避局部最优解。尽管新增分量可能会因为增加不必要的分量而导致ELBO小幅下降,但随后归并操作会拒绝较差的新增分量。

3.4.2 归并机制

归并操作有助于优化全局数据目标。由于推理的代价随分量数目K的增加而增加,因此为了保持较小的K,设计归并操作时可以将两个分量归并为一个分量。

具体来说,通过随机选择分量ka和kb,比较现有模型q和候选模型q′的ELBO值。由于L(q′)的计算中除了Lentropy项,其他都是U统计量的线性函数,因此可以提前对所有可能的归并对计算Lentropya,b项,每个批次最多需要保存K(K-1)/2个标量。这种预计算会加快后续的归并操作,有助于改善模型质量。通过在每轮数据访问中执行一次新增操作和若干次归并操作,几轮访问后可以提升先验模型的质量,使其结构更加紧凑。

3.5 图像去噪算法

(38)

(39)

在完成{ν,S}的更新后,固定全局变分因子,对所有局部变分变量参数进行多轮更新,直至最大化变分目标L′:

L′≜Eq[lb(y,x,IS,.)-lbq(x,.)]

(40)

综上,算法2给出图像去噪算法。

算法2图像去噪算法

输入噪声图像y,通用先验模型

1.Repeat

7.Until ELBO收敛

8.Return

4 实验结果与分析

本文从BSDS数据集[28]中挑选的200张图像中均匀采样出200万个图像块作为训练集,用来学习通用外部先验知识。实验测试用图选取部分经典合成图像及BSDS图像集中部分图像,如图2所示。

图2 部分测试用图Fig.2 Part of images for testing

利用最大似然估计来计算图像块均值m的参数τ和σ,类似地,使用计算块的经验协方差矩阵来估计精度矩阵S中的超参数n和V,并利用变分学习算法使模型能根据观测数据进行自适应调整。

采用峰值信噪比(PSNR)作为算法的客观度量,PSNR计算过程如下:

(41)

图像块大小在模型中扮演重要的角色。过大会导致边缘模糊,丢失纹理细节,过小将产生锯齿效果。本文通过实验对比,选择8×8的块大小,能较好地适配不同的噪声水平。

使用标准差分别为30、40、50、75的加性高斯噪声对原始图像进行污染。将本文方法的去噪结果与BM3D[4]、LSSC[15]、EPLL[12]、PGPD[10]、NL-Bayes[8]、EPPGIC[11]等优秀去噪算法进行对比,结果如表1所示。表中数值为每个噪声标准差下所有测试图像的平均PSNR值,加粗数据为最高值。

表1 去噪性能对比Table 1 Comparison of denoising performance dB

对于所选的经典合成测试图片,许多优秀算法经过调优取得了出色结果。而本文方法在效果上优于大部分算法,并且在更大的BSDS测试集合中取得了更好的结果,这显示了本文模型对于大图片集上非参学习的应用价值。

EPLL可以看作是本文模型的简化,其模型的分量数目是固定的。对于所有的噪声级别和数据集,本文模型内外部先验结合的方法均优于EPLL,在性能上有所提升,表明贝叶斯非参学习方法能够得到合适的模型聚类分量数目。

图3为各方法对噪声标准差为30时的Barbara图像的去噪结果对比。本文方法在PSNR值上优于LSSC、EPLL 、NL-Bayes、PGPD方法,相比BM3D和EPPGIC方法在PSNR值上略弱,但是在局部恢复细节方面有着突出的表现,如图3(e)所示,原始Barbara图中腿部的两处黑点区域,在上述两种方法的复原图中被过渡平滑了,而本文方法在恢复纹理的同时,上述两处黑点仍保持清晰可见。EPLL方法尽管也复原了黑点,但是腿部条纹恢复不太理想。

图3 去噪视觉效果对比1Fig.3 Comparison 1 of visual effects after denoising

在更高噪声级别下,各种方法的复原都出现了一些伪像。如图4所示,BM3D、PDGD背景有大量的波纹,EPPGIC的背景有较多混杂的伪像,NL-Bayes方法在蛇身周边有较多伪像,EPLL方法背景平滑度较弱,本文方法在背景有少量伪影,但无论在PSNR还是视觉观感上都具有一定优势。图4噪声标准差为75。

图4 去噪视觉效果对比2Fig.4 Comparison 2 of visual effects after denoising

5 结束语

本文构建一种基于狄利克雷过程的可扩展高斯混合模型用于图像去噪。该模型从干净图像数据库中学习通用先验,借助聚类分量新增及归并机制和模型中统计量的可累加性,使模型复杂性可以自适应于训练观测图像的变化。实验结果表明,相比传统去噪模型,该模型能取得更高的峰值信噪比,复原效果更佳。后续将从退化图像中学习内部先验,以此捕获图像内部自相似性,实现内外部先验结合,使模型可用于图像填充、去模糊与超分辨等任务。