CN102938138A - 一种基于多元统计模型的分形小波自适应图像去噪方法 - Google Patents

Abstract

本发明公开了基于多元统计模型的分形小波自适应图像去噪方法，包括以下步骤：步骤1：对含噪图像进行同态变换；通过同态变换，含乘性噪声的原图像I_B转换为只含加性噪声的图像I_B′；步骤2：先对含噪信号f(k)进行分形小波变换，选择小波基和小波分解层数j，得到相应的小波系数

和

步骤3：选择MGGD多元统计模型自适应求解参数α和β；在对自然图像的小波系数分布情况进行分析后，获得最适合的参数值α和β；步骤4：对分解得到的小波系数

和利用分形小波编码方法对噪声图像进行无噪预测编码；步骤5：利用

和

进行小波重构，得到估计信号和

即为去噪后的图像信号。与其它算法相比较，具有更好的去噪效果和更强的边缘保持能力，特别适用于消除高斯和脉冲混合噪声。

Description

一种基于多元统计模型的分形小波自适应图像去噪方法

技术领域

本发明涉及图像处理技术领域，尤其涉及的是一种基于多元统计模型的分形小波自适应图像去噪方法。

背景技术

由于图像信号在获取、传输和存储过程中，不可避免地受到噪声的干扰，噪声降低了图像质量，淹没了图像的边缘和细节特征，给图像分析和后续处理带来困难。图像噪声的消除是图像处理中的一个重要研究内容，能否有效地滤除噪声直接影响着图像处理后续工作的进行。在进一步进行边缘检测、图像分割、特征提取和模式识别等处理之前，采用适当的方法去除噪声是一个非常重要的预处理步骤^[1]。如何在有效去除噪声的同时保持图像细节的清晰度和图像的对比度成了人们研究的热点。

传统的图像去噪算法是根据图像频谱分布规律，从频率上将图像中的有用信息与噪声分开，例如小波方法去噪^[2-5]。这些方法一般认为噪声的能量集中于图像的高频部分，而图像的有用信息的频谱则分布于图像低频部分的一个有限区域内。然而，在许多情况下，图像的有用信息中也有分布于图像的高频部分的，例如图像边缘，所以采用这些方法在去除噪声的同时也损失了部分有用信息，即缺乏特征保持性。一些改进方法采用了如Contourlet等其它变换代替小波变换^[6]，取得了更好的效果，但其基本假设不变。而噪声在图像的低频部分也有一定的分量，简单地滤除图像的高频成分无法去除这部分噪声分量，即没有有效地将图像有用信息与噪声数据区分开来。如中值滤波去噪是对图像中所有的像素进行滤波，改变了图像中未被脉冲噪声污染的像素点，所以在有效滤除脉冲噪声的同时，会出现对图像的边缘细节和纹理部分过度平滑，造成去噪后的图像清晰度较低。中值滤波算法在去除噪声时仅考虑了邻域内像素的排序信息，忽略了像素的时序信息，因此会在边缘处产生抖动并会删除一些重要的图像细节^[7-9]。空域去噪中的高斯滤波能有效地去除图像平滑区域的噪声，但由于高斯滤波器^[10]是各向同性的，对边缘和细节不加区分，因此该方法容易造成图像边缘和细节模糊。采用软阈值去噪时很难确定一种对所有图像都适用的阈值选取方法，且容易产生图像的伪吉普斯现象^[11-13]。近年来很多研究工作集中于保留图像细节上，但大多方法仍然假设图像分段平滑，这样虽然能保留边缘信息，但对纹理细节的保留贡献不大。如TV及其改进算法通过最小化图像全变分得到分段平滑图像^[14-17]。现有的算法通常将图像各部分分为平滑与边缘两种模式，通过滤除平滑模式图像中的高频分量进行去噪，既将其处理为理想的平滑模式。当图像含有纹理时，纹理也会被当作以上两种模式处理，所以会损失部分细节信息。为了保留图像的纹理信息，一些学者提出了针对纹理自相似性的图像去噪方法。如分形小波去噪算法^[18，19]利用图像分形小波变换中块自相似性来调整尺度因子。通过分形小波预测编码来达到图像去噪的目的。虽然该方法能取得比较良好的去噪效果，但是该方法在图像去噪过程中边缘细节的保持方面还有待于加强。非局部均值方法通过假设相邻区域具有自相似性，搜索领域内的相似图像区域并进行加权平均实现图像去噪^[20]。基于混合线性模型的去噪方法不假设图像分段平滑，仅假设图像具有自相似性，利用图像的相似性区分图像信号与噪声^[21]。通用离散图像去噪算法假设图像具有平稳分布，用非参数估计方法统计图像块的分布，再用最小代价准则来实现图像去噪^[22]。静态小波变换(SWT)利用时间不变性来实现图像去噪^[23]。基于PDE的非线性扩散滤波方法(P-M)^[24]是一种非线性的各向异性去噪方法，各向异性的去噪模型根据图像的梯度值确定扩散的速度，能够兼顾噪声消除和边缘保持两方面的要求。尽管P-M方法在抑制噪声与保留图像重要特征方面取得了一定的效果，但却表现出病态且不稳定。

参考文献：

[1]A.Pizurica，W.Philips.Estimating the probability ofthe presence of a signal of interest in multiresolution single andmultiband image denoising[J].IEEE Trans.Image Proces.，2006，15(3)：654-665.

[2]A.Chambolle，R.A.DeVore，N.Lee，et al.Nonlinear wavelet image processing：variational problems compression andnoise removal through wavelet shrinkage[J].IEEE Trans.Image Process.，1998，7(3)：319-335.

[3]S.G.Chang，B.Yu，M.Vetterli.Adaptive wavelet thresholding for image denoising and compression[J].IEEE Trans.Image Process.，2000，9(9)：1532-1546.

[4]I.Johnstone，B.Silverman.Empirical Bayes selection ofwavelet thresholds，Ann.Stat.，2005，33(4)：1700-1752.

[5]H.Othman，S.E.Qian.Noise reduction of hyperspectral imagery using hyb rid spatial-spectral derivative-domain waveletshrinkage[J].IEEE Trans.Geosci.Remote，2006，44(2)：397-408.

[6]Z.F.Zhou，P.L.Shui.Contourlet based image denoising algorithm using directional windows[J].Electronics Letters，2007，43(2)：92-93.

[7]T.C.Lin.A new adaptive center weighted median filter for suppressing impulsive noise in images[J].InformationSciences，2007，177(4)：1073-1087.

[8]S.Q.Yuan，Y.H.Tan.The solutions of equation based noise detector for an adaptive median filter[J].Pattern Recognition，2006，39(11)：2252-2257.

[9]S.M.Mahbubur Rahman，Md.Kamrul Hasan，Wavelet-domain iterative center weighted median filter for imagedenoising[J].Signal Processing，2003，83(5)：1001-1012.

[10]K.R.Castleman.Digital image processing[M].New Jersey：Prentice-Hall，1979.

[11]K.Q.Huang，Z.Y.Wu，G.S.K.Fung，et al.Color image denoising with wavelet thresholding based on human visualsystem model[J].Signal Processing：Image Communication，2005，20(2)：115-127.

[12]G.Y.Chen，T.D.Bui，A.Krzyak.Image denoising with neighbour dependency and customized wavelet and threshold[J].Pattern Recognition，2005，38(1)：115-124.

[l3]D.L.Donoho.Denoising by softthresholding[J].IEEE Transactions on Imformation Theory，1995，41(3)：613-627.

[14]L.I.Rudin，S.Osher，E.Fatemi.Nonlinear total variation based noise renoval algorithms[J].Physica D：NonlinearPhenomena，1992，60(1-4)：259-268.

[15]Q.Chen，P.Montesinos，Q.S.Sun，et al.Adaptive total variation denoising based on difference curvature[J].lmage andVision Computing，2010，28(3)：298-306.

[16]G.Steidl，J.Weickert.Relationsbetweensoftwayeletshrinkageandtotalvariationdenoising[J].PatternRecognition，Lecture Notes in Computer Science，2002，2449：198-205.

[l7]D.Lazzaro，L.B.Montefusco.Edge-preserving wavelet thresholding for image denoising[J].Journal of Computationaland Applied Mathematics，2007，210：222-231.

[18]Mohsen Ghazel，George H.Freeman，Edward R.Vrscay.Fractal-Wavelet Image Denoising Revisited[J].IEEETRANSACTIONS ON IMAGE PROCESSING，2006，15(9)：2669-2675.

[19]P.Luo，X.P.Gao.Image denoising algorithm based on dual tree com plex wavelettransforn[J].Acta Photonica Sinica，2008，37(3)：604-608.

[20]A.Buades，B.Coll，J.M.Morel.Image denoising by non-local averaging[C].IEEE International Conference onAcoustics，Speech and Signal Processing，Philadelphia，PA，United states，2：1125-1128，2005.

[21]曹扬，罗予频，杨士元.基于混合线性模型的图像去噪[J].计算机学报，2009，32(11)：2260-2263.

[22]T.Weissman，E.Ordentlich，G.Seroussi，et al.Universal discrete denoising：Known channel[J].IEEE Transactions onInformation Theory，2005，51(1)：5-28.

[23]X.H.Wang，S.H.Robert，Y.H.Song.Microarray image enhancement by denoising using stationary wavelet transform[J].IEEE Transactions on Nanobioscience，2003，2(4)：184-189.

[24]B.I.Lee，S.H.Lee，T.S.Kim，et al.Harmonic decomposition in PDE-based denoising techniquefor magnetic resonanceelectrical impedance tomography[J].IEEE Transactions on Biomedical Engineering，2005，52(11)：1912-1920.

[25]G.Davis.Awavelet-based analysisfractal image compression[J].IEEETrans.Image Process，1998，7(2)：141-154.

[26]H.Krupnik，D.Malah，E.Karnin.Fractal representation of images via the discrete wavelettransform[J].IEEE Conv.Elect，1995，15(3)：7-9.

[27]E.R.Vrscay.A generalized class offractal-wavelet transforms for image representation and compression[J].Canad.J.Elect.Comput.Eng.，1998，23(1)：69-84.

[28]A.Benazza Benyahia，J.C.Pesquet.Building robust wavelet estimators for multicomponent images using Steins’principle[J].IEEETrans.Image Proces.，2005，14(11)：1814-1830.

[29]Steve De Backer，Aleksandra

Bruno Huysmans，et al.Denoising of multicomponent images using waveletleast-squares estimators[J]Image andVision Computing，2008，26：1038-1051.

[30]C.S.Tong，M.Wong.Adaptative approximate nearest neighbor search forfractal image compression[J].IEEE Trans.Image Process.，2002，11(6)：605-615.

[31]H.Chipman，E.Kolaczyk，R.McCulloch.Adaptive Bayesian wavelet shrinkage[J].J.Amer.Statist.Assoc.1997，92(440)：1413-1421.

[32]S.J.Roberts，D.Husmeie r，I.Rezek，et al.Bayesian approaches to Gaussian mixture modeling[J].IEEE Trans.Pattern Anal.Machine Intell.，1998，20(11)：1133-1142.

[33]Dongwook Cho，Tien D.Bui.Multivariate statistical modeling for image denoising using wavelet transforms[J].SignalProcessing：Image Communication，2005，20：77-89.

[34]J.Portilla，V.Strela，M.Wainwright，et al.Image denoising using scale mixtures of Gaussians in the waveletdomain[J].IEEE Trans.Image Process，2003，12(11)：1338-1351.

[35]H.Yan，Y.H.Yu，M.F.Zhao.Mixed statistical model image denoising based on shift-invariant non-aliasing Contourlettransform[J].Optics and Precision Engineering，2010，18(10)：2269-2279.

发明内容

本发明所要解决的技术问题是针对现有技术的不足提供一种基于多元统计模型的分形小波自适应图像去噪方法。

本发明的技术方案如下：

一种基于多元统计模型的分形小波自适应图像去噪方法，包括以下步骤：

步骤1：对含噪图像进行同态变换；通过同态变换，含乘性噪声的原图像I_B转换为只含加性噪声的图像I_B′；

步骤2：先对含噪信号f(k)进行分形小波变换，选择小波基和小波分解层数j，得到相应的小波系数和

步骤3：选择MGGD多元统计模型自适应求解参数α和β；在对自然图像的小波系数分布情况进行分析后，获得最适合的参数值α和β；

步骤4：对分解得到的小波系数

和

利用分形小波编码方法对噪声图像进行无噪预测编码；

步骤5：利用

和

进行小波重构，得到估计信号

和

即为去噪后的图像信号。

所述的方法，所述步骤3的具体方法为：在所述的分析过程中，利用Daubechies20滤波器对图像集进行分形小波分解，寻找最接近每个子带分析的MGGD多元统计模型；确定最适合的参数问题就可以转化为数据拟合问题；如果考虑两个分布函数均方差，残差的L²范数可以用下式来定义：

$R_1={\left|\right|p_2\left(\stackrel{\mathrm{\ρ}}{x}\right|\mathrm{\α},\mathrm{\β})-p_1\left(\stackrel{\mathrm{\ρ}}{x}\right)\left|\right|}_{L_2}^2$

$=\underset{i}{\mathrm{\Σ}}{(p_2\left({\stackrel{\mathrm{\ρ}}{x}}_i\right|\mathrm{\α},\mathrm{\β})-p_1\left({\stackrel{\mathrm{\ρ}}{x}}_i\right))}^2$

为此，利用Matlab的优化工具箱lsqcurvefit()函数来分析，通过最小化R₁得到最接近

的

及其参数α、β；定义对数残差的L²范数为：

$R_2={\left|\right|1n{p}_2\left(\stackrel{\mathrm{\ρ}}{x}\right|\mathrm{\α},\mathrm{\β})-1n{p}_1\left(\stackrel{\mathrm{\ρ}}{x}\right)\left|\right|}_{L_2}^2$

$={\left|\right|1n\frac{p_2\left(\stackrel{\mathrm{\ρ}}{x}\right|\mathrm{\α},\mathrm{\β})}{p_1\left(\stackrel{\mathrm{\ρ}}{x}\right)}\left|\right|}_{L_2}^2$

$=\underset{i}{\mathrm{\Σ}}{(1n{p}_2\left({\stackrel{\mathrm{\ρ}}{x}}_i\right|\mathrm{\α},\mathrm{\β})-{1\mathrm{np}}_1\left({\stackrel{\mathrm{\ρ}}{x}}_i\right))}^2$

由于没有样本时

可能为0和

取值小得不合理时，利用lsqcurvefit()函数可能得到不准确解；风险值R₂和参数α、β之间的关系可以通过观测获得。

本发明提出一种基于多元统计模型的分形小波自适应图像去噪算法。根据分形小波变换具有低熵性、多分辨率、去相关性和选基的灵活性等特点，通过选用扩展的GGD模型MGGD(multivariate generalized Gaussian distribution)模型建立多元统计模型的基础上利用分形小波变换来完成图像去噪。通过最小化残差R₁得到最接近

的

及自适应调整参数α和β；通过使用四叉树分割来实现对噪声图像自适应地预测分形小波无噪图像编码以达到去噪目的。与其它算法相比较，具有更好的去噪效果和更强的边缘保持能力，特别适用于消除高斯和脉冲混合噪声。

附图说明

图1为二维分形小波变换；

图2为Lena图像、噪声图像和几种图像去噪方法去噪后的结果图像；(a)原始图像，(b)添加高斯噪声后的图像，(c)基于小波变换的隐马尔模型方法，(d)边缘保护小波方法，(e)局部双变收缩方法，(f)分形小波方法，(g)多元统计方法，(h)本发明方法；

图3为Lena图像、噪声图像和几种图像去噪方法去噪后的结果图像的灰度直方图；(a)原始图像，(b)添加高斯噪声后的图像，(c)基于小波变换的隐马尔模型方法，(d)边缘保护小波方法，(e)局部双变收缩方法，(f)分形小波方法，(g)多元统计方法，(h)本发明方法；

图4为Lena图像和几种图像去噪方法去噪后的结果图像的边缘检测图；(a)原始图像，(b)基于小波变换的隐马尔模型方法，(c)边缘保护小波方法，(d)局部双变收缩方法，(e)分形小波方法，(f)多元统计方法，(g)本发明方法；

具体实施方式

以下结合具体实施例，对本发明进行详细说明。

1分形小波变换去噪

分形小波变换(Fractal Wavelet Transforms，FWT)是在分形图像压缩过程中，为了减少块效应和计算的复杂性而引入的^[25-29]。分形小波变换操作涉及将小波系数子树缩放和复制到较低的子树，完全类似于分形图像编码器在空间域的操作。分形小波去噪的本质是从噪声图像中预测出无噪声图像的分形编码。

1.1分形小波变换图像编码简介

二维信号(图像)的离散小波变换(DWT)系数是一个标准的矩阵式排列图。假设二维小波基函数是通过一维尺度函数和小波函数的合适的张量积来构造的，则可以得到图1所示的二维分形小波变换图^[18]。每个块

和

0≤k≤K，对应的小波系数分别为

和

是由水平、垂直和对角三个块系数构成的四叉树。在这个矩阵和唯一的子树中，任何小波系数

γ∈{h，v，d}，都是以

为根元素。分形小波变换图像编码可以用“拼贴编码”来完成。

“拼贴编码”程序对图像产生分形小波编码的过程如下：首先考虑一套固定的父子子树的等级值且

对每个未编码子子树

找到其对应的父子树

及相应的尺度系数s_{i，j，i′，j′}，使用(1)计算出来的“拼贴距离”最小。

${\mathrm{\Δ}}_{i,j,i^{\′},j^{\′}}^{\mathrm{\γ}}=\left|\right|X_{k_2^{*}}^{\mathrm{\γ}}-s_{i,j,i^{\′},j^{\′}}{X}_{k_1^{*},i^{\′},j^{\′}}^{\mathrm{\γ}}\left|\right|---\left(1\right)$

分形小波编码的结果由三部分组成：1)父子子树的等级值

2)总数为

个小波系数

这些小波系数在分形小波变换过程中保持不变；3)总数为

个尺度系数和

个父子子树索引(i^γ，j^γ)。

在严格执行图像保真度的限制条件下，标准分形小波编码方案中的三个基本子带可以使用共同的父带和尺度系数。分形小波编码通常是先从已存储小波系数的小波系数树开始。然后通过分形小波缩放和复制的迭代方式生成一个逼近原始图像的“固定点”的小波系数矩阵。“拼贴距离”越小，逼近原始图像的效果越好。

给定一个小波系数树R，假设子子树为

用向量

来表示；父子树为

用向量

来表示。在实践中，父子子树的拼贴编码应使L²范数的误差

最大限度地减少。然后，把图像的小波变换作为一个随机信号。这样，小波系数

和

就可以被看作是从随机变量中抽取出来的随机样本，并分别代表父子树及其对应的子子树小波系数分布。最佳最小二乘尺度系数可以写成^[18]：

$s^{*}=\frac{\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\stackrel{\mathrm{\ρ}}{x}}_k{\stackrel{\mathrm{\ρ}}{x}}_k}{\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left({{\stackrel{\mathrm{\ρ}}{x}}_k}^{\′}\right)}^2}=\frac{\underset{k=1}{\overset{k}{\mathrm{\Σ}}}{\stackrel{\mathrm{\ρ}}{x}}_k{\stackrel{\mathrm{\ρ}}{x}}_k}{\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left({{\stackrel{\mathrm{\ρ}}{x}}_k}^{\′}\right)}^2}---\left(2\right)$

由于统计的样本量有限，严格地说，上面的表达式只是对随机变量的近似统计量。

1.2无噪图像的分形小波预测编码

如果把图像的分形小波变换看作为一个随机信号，那么分形小波编码过程可以归纳为均方差(MSE)的估计问题：对每个未编码的子子树

来说，找到最佳的父子树

及其对应的尺度系数时的均方差可以用(3)来计算：

$\mathrm{MES}=E\left[{({\stackrel{\mathrm{\ρ}}{x}}_k-s_{i,k}^{*}{{\stackrel{\mathrm{\ρ}}{x}}_i}^{\′})}^2\right]$

$=E\left[{\stackrel{\mathrm{\ρ}}{x}}_k^2\right]+s_{i,k}^{*2}E\left[{\left({{\stackrel{\mathrm{\ρ}}{x}}_i}^{\′}\right)}^2\right]-2s_{i,k}^{*}E\left[{\stackrel{\mathrm{\ρ}}{x}}_k{{\stackrel{\mathrm{\ρ}}{x}}_i}^{\′}\right]---\left(3\right)$

其中，尺度系数的估计是通过式(2)完成的。

然而，在实际中必须对噪声图像进行无噪预测编码。用

和分别表示噪声图像的子子树和父子树，则采用正交小波基时，噪声图像和无噪图像对应的小波系数之间满足如下关系：

$\left\{\begin{array}{c}{\stackrel{\mathrm{\ρ}}{x}}_i={\stackrel{\mathrm{\ρ}}{x}}_i+W_{{\stackrel{\mathrm{\ρ}}{x}}_i}\\ {{\stackrel{\mathrm{\ρ}}{x}}^{\′}}_k={{\stackrel{\mathrm{\ρ}}{x}}^{\′}}_k+W_{{{\stackrel{\mathrm{\ρ}}{x}}^{\′}}_k}\end{array}\right.---\left(4\right)$

其中，和

可认为是统计独立的相同分布的AWGN处理^[11]。

也可被认为是统计独立的。可以通过确保子和父子树不重叠来达到独立。并假设图像和噪声信号都是独立的，这种假设是合理可行的。由上述假设可以很容易地得到：

$\left\{\begin{array}{c}E\left[{\stackrel{\mathrm{\ρ}}{x}}_i^2\right]=E\left[{\stackrel{\mathrm{\ρ}}{x}}_i^2\right]-{\mathrm{\σ}}_W^2\\ E\left[{{\stackrel{\mathrm{\ρ}}{x}}^{\′}}_i^2\right]=E\left[{{\stackrel{\mathrm{\ρ}}{x}}^{\′}}_i^2\right]-{\mathrm{\σ}}_W^2\\ E\left[{\stackrel{\mathrm{\ρ}}{x}}_i{{\stackrel{\mathrm{\ρ}}{x}}^{\′}}_k\right]=E\left[{\stackrel{\mathrm{\ρ}}{x}}_i{{\stackrel{\mathrm{\ρ}}{x}}^{\′}}_k\right]\end{array}\right.---\left(5\right)$

在文献[13]中，噪声方差MES被错误地加入，而不是从噪声二阶矩估计中减去。这可以说明Donoho等人在图像去噪过程中遇到的困难。文献[18]在分形小波编码去噪方面取得了一些进步。由上述推导，提出以下重要观测：

从对噪声的观测过程中，可以对原始无噪声图像的统计进行估计。当这些估计具有鲁棒性时，它们可以用来估计无噪声图像的分形小波编码。根据对噪声观测的统计来估计无噪声图像拼贴方差MES。对于给定的子树，选择最小的父子子树方差MES估计。当含噪父子树和子子树的能量远远大于噪声方差时，统计估计量的鲁棒性得以实现，即：

$\left\{\begin{array}{c}E\left[{\stackrel{\mathrm{\ρ}}{x}}_i^2\right]>>{\mathrm{v\σ}}_W^2\\ E\left[{{\stackrel{\mathrm{\ρ}}{x}}^{\′}}_i^2\right]>>{\mathrm{v\σ}}_W^2\end{array}\right.---\left(6\right)$

参数v＞＞1。通过对各种测试图像的观测，v在1.5至2.5之间能够获得最佳分形小波编码去噪效果，在本发明实验中选取v＝2。当(6)的鲁棒性标准不满足时，(5)可能会产生负的二阶矩估计，(2)的预测尺度系数可能是极大的。因此，该预测方法不能适用于稀疏信号信息。对于这种特殊的噪声子树，减小噪声尺度系数的幅度有助于抑制一些噪声。采用(7)来减小噪声尺度系数的幅度：

$s_{i,k}^{*}\≈\min (\frac{E\left[{{\stackrel{\mathrm{\ρ}}{x}}^{\′}}_k^2\right]}{{\mathrm{v\σ}}_W^2},\frac{E\left[{\stackrel{\mathrm{\ρ}}{x}}_k^2\right]}{{\mathrm{v\σ}}_W^2})\×{\hat{s}}_{i,k}^{*}---\left(7\right)$

虽然上述算法是对标准分形小波编码方法的概述，但可以将它推广到其它分形小波编码方案，如使用拼贴误差分解标准的四叉树分形小波分解编码方法。使用子和父三个子带树的分形小波去噪的目的是为使编码有足够大的规模。要不然，如果对局部统计的估计不好可能导致去噪效果不佳。

2.基于多元统计模型的分形小波自适应图像去噪

虽然分形小波预测编码能够在图像去噪过程中得到比较清晰的去噪图像，但是去噪同时丢失了部分边缘信息。由于图像的小波系数的分布近似于高斯分布，为了得到保留良好边缘和纹理细节的去噪图像，本发明利用多元统计模型来优化分形小波编码。

2.1多元统计模型的贝叶斯估计

设I_A为不带噪声自然图像，I_B为带噪声图像，它们之间的关系可以用式(8)来表示：

I_B＝I_A+σC                            (8)

其中，C表示零均值高斯白噪声，C～N(0，1)；σ表示噪声方差。

对噪声图像I_B进行多分辨率分形小波分解后得到第j层第i个水平小波系数

垂直小波系数

和对角小波系数由小波变换的线性关系，可以得出：

$\left\{\begin{array}{c}{y}_{i,j}^h=x_{i,j}^h+{\mathrm{\σz}}_{i,j}^h\\ y_{i,j}^v=x_{i,j}^v+{\mathrm{\σz}}_{i,j}^v\\ y_{i,j}^d=x_{i,j}^d+{\mathrm{\σz}}_{i,j}^d\end{array}\right.---\left(9\right)$

其中，

和

分别表示图像I_A的水平、垂直和对角小波系数；

和

分别表示噪声C的水平、垂直和对角小波系数。

设是一个d维小波系数向量，

其中x₁是在去噪过程中必须考虑的小波系数，(x₂，Λ，x_d)是在去噪过程中可以考虑的相关小波系数(如邻域、父子小波系数)。为了简化计算式，用单下标小波系数x_k，y_k，z_k分别取代

和

噪声图像和噪声对应的小波系数向量分别为

和

则有：

$\stackrel{\mathrm{\ρ}}{y}=\stackrel{\mathrm{\ρ}}{x}+\mathrm{\σ}\stackrel{\mathrm{\ρ}}{z}---\left(10\right)$

在计算过程中，关注的焦点是未知小波系数向量

的估计值。而

的估计值的计算又依赖于噪声图像I_B所对应的小波系数向量

本发明利用最大后验概率(maximum aposteriori MAP)算子最大化概率来估计

可以通过(11)来计算。

由于式(11)中

只是一个已知常量，不影响计算结果。则在最小概率误差下，

最佳值可以通过和

来估计。

首先，由于高斯噪声的每个向量是独立且相等分布的，

满足多元高斯分布

因此，

可以通过(12)来计算：

$1\mathrm{np}\left(\stackrel{\mathrm{\ρ}}{y}\right|\stackrel{\mathrm{\ρ}}{x})=1\mathrm{np}\left(\stackrel{\mathrm{\ρ}}{z}\right)$

$=1n\frac{1}{{\left(2\mathrm{\π}\right)}^{d/2}{\left|{\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{z}}\right|}^{1/2}}\exp \{-\frac{{(\stackrel{\mathrm{\ρ}}{y}-\stackrel{\mathrm{\ρ}}{x})}^T{\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{z}}^{-1}(\stackrel{\mathrm{\ρ}}{y}-\stackrel{\mathrm{\ρ}}{x})}{2}\}$

$=1n\frac{1}{{\left(2\mathrm{\π}\right)}^{d/2}{\left|{\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{z}}\right|}^{1/2}}\exp \{-\frac{{(\stackrel{\mathrm{\ρ}}{y}-\stackrel{\mathrm{\ρ}}{x})}^T(\stackrel{\mathrm{\ρ}}{y}-\stackrel{\mathrm{\ρ}}{x})}{2{\mathrm{\σ}}^2}\}---\left(12\right)$

$=-\frac{d}{2}1n\left(2\mathrm{\π}{\mathrm{\σ}}^2\right)-\frac{{(\stackrel{\mathrm{\ρ}}{y}-\stackrel{\mathrm{\ρ}}{x})}^T(\stackrel{\mathrm{\ρ}}{y}-\stackrel{\mathrm{\ρ}}{x})}{{2\mathrm{\σ}}^2}$

其次，必须为

建立合适的统计模型。为此，对样本图像的小波系数进行检测，发现它们的分布近似于高斯分布。并结合文献[31-35]得出高斯混合模型最适合组建最佳模型。因此，选用扩展的GGD模型的MGGD模型可表示为：

$p\left(\stackrel{\mathrm{\ρ}}{x}\right)=\mathrm{vexp}\{-{\left(\frac{{(\stackrel{\mathrm{\ρ}}{x}-\mathrm{\μ})}^T{\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{x}}^{-1}(\stackrel{\mathrm{\ρ}}{x}-\mathrm{\μ})}{\mathrm{\α}}\right)}^{\mathrm{\β}}\}---\left(13\right)$

其中，α和β为模型的球形参数；v是α、β和协方差矩阵

的归一化常数。如果把式(11)中的

定义为未知函数

则由(11)和(12)可得：

其中，

为方括号内的部分。假设连续、可微，如果存在

满足

则最大化

可改由(15)来计算。

$\▿F\left(\stackrel{\mathrm{\ρ}}{x}\right)=\frac{\∂F\left(\stackrel{\mathrm{\ρ}}{x}\right)}{\∂\left(\stackrel{\mathrm{\ρ}}{x}\right)}=0---\left(15\right)$

最后，利用(15)对(14)进行化简可得：

$\▿F\left(\stackrel{\mathrm{\ρ}}{x}\right)=-\frac{\stackrel{\mathrm{\ρ}}{x}-\stackrel{\mathrm{\ρ}}{y}}{{\mathrm{\σ}}^2}+\▿f\left(\stackrel{\mathrm{\ρ}}{x}\right)=0\⇔\stackrel{\mathrm{\ρ}}{x}=\stackrel{\mathrm{\ρ}}{y}+{\mathrm{\σ}}^2\▿f\left(\stackrel{\mathrm{\ρ}}{x}\right)---\left(16\right)$

这样，假设μ＝0，利用MGGD模型可以得出(16)的更明确的计算为：

$\▿f\left(\stackrel{\mathrm{\ρ}}{x}\right)=-\frac{2\mathrm{\β}}{{\mathrm{\α}}^{\mathrm{\β}}}{\left({\stackrel{\mathrm{\ρ}}{x}}^T{\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{x}}^{-1}\stackrel{\mathrm{\ρ}}{x}\right)}^{\mathrm{\β}-1}{\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{x}}^{-1}\stackrel{\mathrm{\ρ}}{x}---\left(17\right)$

由(16)和(17)可得：

$\stackrel{\mathrm{\ρ}}{x}=\stackrel{\mathrm{\ρ}}{y}-\frac{2{\mathrm{\σ}}^2\mathrm{\β}}{{\mathrm{\α}}^{\mathrm{\β}}}{\left({\stackrel{\mathrm{\ρ}}{x}}^T{\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{x}}^{-1}\stackrel{\mathrm{\ρ}}{x}\right)}^{\mathrm{\β}-1}{\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{x}}^{-1}\stackrel{\mathrm{\ρ}}{x}$

$={(I+\frac{2{\mathrm{\σ}}^2\mathrm{\β}}{{\mathrm{\α}}^{\mathrm{\β}}}{\left({\stackrel{\mathrm{\ρ}}{x}}^T{\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{x}}^{-1}\stackrel{\mathrm{\ρ}}{x}\right)}^{\mathrm{\β}-1}{\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{x}}^{-1}\stackrel{\mathrm{\ρ}}{x})}^{-1}\stackrel{\mathrm{\ρ}}{y}---\left(18\right)$

$=({\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{x}}+\frac{2{\mathrm{\σ}}^2\mathrm{\β}}{{\mathrm{\α}}^{\mathrm{\β}}}{\left({\stackrel{\mathrm{\ρ}}{x}}^T{\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{x}}^{-1}\stackrel{\mathrm{\ρ}}{x}\right)}^{\mathrm{\β}-1}I){\mathrm{\Σ}}_{\stackrel{\mathrm{\ρ}}{x}}\stackrel{\mathrm{\ρ}}{y}$

为了解决(18)中没有通解的问题，可以通过定义α、β及协方差矩阵

为特殊值或数值来达到求解的目的。本发明在实验中采用最小二乘法来自适应求解参数α和β。

2.2算法描述

利用基于多元统计模型的分形小波自适应图像去噪算法进行去噪的具体步骤描述如下：

步骤1：对含噪图像进行同态变换。通过同态变换，含乘性噪声的原图像I_B转换为只含加性噪声的图像I_B′。

步骤2：先对含噪信号f(k)进行分形小波变换，选择合适的小波基和小波分解层数j，得到相应的小波系数

和

例如，当噪声功率较低时，选择db1和harr小波去噪 效果较好，若选用sym4和bior4等小波基做分解，则会使信号的奇异点变得圆滑，去噪 信号有较大的畸变，信噪比也较低；当噪声功率较大时，小波基的洗择影响不大。而对 于类似正弦的比较平滑的信号，不管噪声功率大小，选择sym8，demy，bior4和coif3 等小波基去噪效果较好，而用db1和harr几乎不能恢复信号原来的形状。分解层数一 般为3-5层，一般5层和4层分解信噪比比3层分解大，平滑段的视觉效果较好，但跳 跃处有较大的毛刺。本发明分解层数选4层。

步骤3：选择MGGD多元统计模型自适应求解参数α和β。在对自然图像(即不 带噪声的图像)的小波系数分布情况进行分析后，通过下面的方法来获得最适合的参数值α和β：

众所周知，自然图像的小波系数的详细分布看起来像零均值高斯分布，如GGD^[32]。选用20幅大小为512×512的测试图像进行样本系数提取。在分析过程中，利用Daubechies20滤波器对图像集进行分形小波分解，寻找最接近每个子带分析的MGGD多元统计模型。确定最适合的参数问题就可以转化为数据拟合问题。如果考虑两个分布函数均方差，残差的L²范数可以用(19)式来定义：

$R_1={\left|\right|p_2\left(\stackrel{\mathrm{\ρ}}{x}\right|\mathrm{\α},\mathrm{\β})-p_1\left(\stackrel{\mathrm{\ρ}}{x}\right)\left|\right|}_{L_2}^2$

$=\underset{i}{\mathrm{\Σ}}{(p_2\left({\stackrel{\mathrm{\ρ}}{x}}_i\right|\mathrm{\α},\mathrm{\β})-p_1\left({\stackrel{\mathrm{\ρ}}{x}}_i\right))}^2---\left(19\right)$

为此，利用Matlab的优化工具箱lsqcurvefit()函数来分析，通过最小化R₁得到最接近的

及其参数α、β。由于数目较少的大小波系数比数目较多的小小波系数在计算中更重要，本发明定义对数残差的L²范数为：

$R_2={\left|\right|1n{p}_2\left(\stackrel{\mathrm{\ρ}}{x}\right|\mathrm{\α},\mathrm{\β})-1n{p}_1\left(\stackrel{\mathrm{\ρ}}{x}\right)\left|\right|}_{L_2}^2$

$={\left|\right|1n\frac{p_2\left(\stackrel{\mathrm{\ρ}}{x}\right|\mathrm{\α},\mathrm{\β})}{p_1\left(\stackrel{\mathrm{\ρ}}{x}\right)}\left|\right|}_{L_2}^2---\left(20\right)$

$=\underset{i}{\mathrm{\Σ}}{(1n{p}_2\left({\stackrel{\mathrm{\ρ}}{x}}_i\right|\mathrm{\α},\mathrm{\β})-{1\mathrm{np}}_1\left({\stackrel{\mathrm{\ρ}}{x}}_i\right))}^2$

由于没有样本时

可能为0和

取值小得不合理时，利用lsqcurvefit()函数可能得到不准确解。这时，风险值R₂和参数α、β之间的关系可以通过观测获得。

步骤4：对分解得到的小波系数

和

利用1.2部分介绍的分形小波编码方法对噪声图像进行无噪预测编码。

步骤5：利用

和

进行小波重构，得到估计信号

和

即为去噪后的图像信号。

3实验结果分析

为验证本发明提出算法的去噪效果并和最近相关文献进行对比研究，本发明选用灰度级为255，像素为(512×512)Lena(图2(a)所示)、Boat和Peppers256图像进行去噪实验。由于许多实际噪声可以近似为高斯分布的白噪声，本发明通过在图像中叠加高斯白噪声来研究和比较图像去噪效果。加入高斯噪声后的图像如图2(b)所示。采用本发明提出的去噪方法和文献中提出的去噪方法得到的结果图像分别如图2(c)-2(h)所示。从图2可以看出，本发明提出的去噪算法得到的图像，不但能够有效地去除噪声，而且原图像中的许多边缘和纹理细节特征也得以很好地保留(如帽羽、瞳孔等)，具有较高的视觉质量。

本发明采用定性和定量的比较方法对本发明提出的去噪方法和相关文献中提出的去噪方法的去噪效果进行比较。由于在去噪过程中不可避免地丢失一些图像信息和残留部分噪声信息，这样图像中一些像素的灰度值及图像的边缘特征会发生相应的变化。因此，可以根据图像的灰度直方图和图像的边缘检测结果图来比较各种去噪方法的去噪效果。图3(a)为原始Lena图像的灰度直方图；图3(b)为加入高斯噪声后图像的灰度直方图；采用本发明提出的去噪方法和文献中提出的去噪方法得到的结果图像的灰度直方图分别如图3(c)-3(h)所示。和原始图像比较可以看出，利用文献中提出的去噪算法得到的图像的像素点的灰度值变化比较大，而利用本发明提出的去噪算法得到的图像的像素点的灰度值变化不明显。说明本发明方法的去噪效果更好。图4(a)为原始Lena图像的边缘检测结果图；采用本发明提出的去噪方法和文献中提出的去噪方法得到的结果图像的边缘检测结果图分别如图4(b)-4(g)所示。可以看出，图中分形小波去噪图像轮廓细节的清晰度较好；多元统计模型去噪图像的边缘特征保持得比较好；综合分形小波和多元统计模型的优点的本发明去噪方法既能够保持图像原有的清晰边缘又可以提高去噪图像的清晰度。

为了定量地衡量去噪效果，通过定义均方误差百分比(RMSE)和峰值信噪比(PSNR)两项指标对去噪效果进行统计分析。对于一幅像素为N×M的图像，定义均方误差百分比为：

$\mathrm{RMES}=\sqrt{\frac{1}{N\×M}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{m=0}{\overset{M-1}{\mathrm{\Σ}}}{[\hat{f}(n,m)-f(n,m)]}^2\×100\%}---\left(21\right)$

定义去噪后图像的峰值信噪比为：

$\mathrm{PSNR}=101g{\left[\frac{f_{\max }}{\mathrm{RMES}}\right]}^2---\left(22\right)$

其中，f_max＝maxf(n，m)，n∈[0，1，Λ，N-1]，m∈[0，1，Λ，M-1]，f(n，m)和分别是原始图像和去噪后的图像在(n，m)位置上的像素点的灰度值。表1、表2和表3分别列出了Lena、Boat和Peppers256图像在不同噪声率下使用不同去噪算法的峰值信噪比和均方误差百分比。由结果可见，本发明所提出的算法的去噪效果明显优于其它算法的去噪效果，尤其在平滑高斯噪声和有脉冲噪声在内的混合噪声效果更显著。总的来说，无论从主观视觉效果还是PSNR和RMSE的客观评价方面来看，本发明算法相对于其它几种算法去噪效果都有明显的改进，既能够很好地消除噪声，又能够较好地保持图像边缘和纹理细节。

需要说明的是，不同的图像及不同的含噪度，其最终采用的参数也会不一样。因此，在使用本发明提出的算法去噪时，应不断地自适应调整参数α和β的值，以便达到使图像最终处理的结果最优。

表1不同等级高斯噪声下去噪性能比较(Lena)

表2不同等级高斯噪声下去噪性能比较(Boat)

表3不同等级高斯噪声下去噪性能比较(Peppers256)

本发明提出的基于多元统计模型的分形小波自适应图像去噪算法，在去噪过程中，首先建立了一个多元统计模型，该模型能够更准确地估计各种相关信息，且模型参数改善比较灵活。然后通过与分形小波去噪方法结合可以选择高品质的图像空间。在适度的噪声方差下可以根据拼贴距离在最好的子树域中找到近优父子树。最后通过从噪声图像中预测出无噪声的图像分形小波编码，从而达到优化去噪的目的。

实验结果表明，本发明提出的去噪方法明显优于作为实验结果中列出的现有算法。该方法在去除噪声的同时，能有效地保持图像的边缘及纹理特征，很好地保留图像的精细结构，取得了良好的去噪效果。在图像复原、增强等诸多领域内都具有重要的应用价值。由于采用了预测小波分形编码，优化了算法结构，算法的处理速度比较快。因此，完全可以达到实时图像处理过程中的去噪预处理对处理速度的要求。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (2)

1.一种基于多元统计模型的分形小波自适应图像去噪方法，其特征在于，包括以下步骤：

步骤1：对含噪图像进行同态变换；通过同态变换，含乘性噪声的原图像I_B转换为只含加性噪声的图像I_B′；

步骤2：先对含噪信号f(k)进行分形小波变换，选择小波基和小波分解层数j，得到相应的小波系数

和

步骤3：选择MGGD多元统计模型自适应求解参数α和β；在对自然图像的小波系数分布情况进行分析后，获得最适合的参数值α和β；

步骤4：对分解得到的小波系数

和

利用分形小波编码方法对噪声图像进行无噪预测编码；

步骤5：利用

和

进行小波重构，得到估计信号

和

即为去噪后的图像信号。

2.根据权利要求1所述的方法，其特征在于，所述步骤3的具体方法为：在所述的分析过程中，利用Daubechies20滤波器对图像集进行分形小波分解，寻找最接近每个子带分析的MGGD多元统计模型；确定最适合的参数问题就可以转化为数据拟合问题；如果考虑两个分布函数均方差，残差的L2范数可以用下式来定义：

$R_1={\left|\right|p_2\left(\stackrel{\mathrm{\ρ}}{x}\right|\mathrm{\α},\mathrm{\β})-p_1\left(\stackrel{\mathrm{\ρ}}{x}\right)\left|\right|}_{L_2}^2$

$=\underset{i}{\mathrm{\Σ}}{(p_2\left({\stackrel{\mathrm{\ρ}}{x}}_i\right|\mathrm{\α},\mathrm{\β})-p_1\left({\stackrel{\mathrm{\ρ}}{x}}_i\right))}^2$

为此，利用Matlab的优化工具箱lsqcurvefit()函数来分析，通过最小化R₁得到最接近的

及其参数α、β；定义对数残差的L²范数为：

$R_2={\left|\right|1n{p}_2\left(\stackrel{\mathrm{\ρ}}{x}\right|\mathrm{\α},\mathrm{\β})-1n{p}_1\left(\stackrel{\mathrm{\ρ}}{x}\right)\left|\right|}_{L_2}^2$

$={\left|\right|1n\frac{p_2\left(\stackrel{\mathrm{\ρ}}{x}\right|\mathrm{\α},\mathrm{\β})}{p_1\left(\stackrel{\mathrm{\ρ}}{x}\right)}\left|\right|}_{L_2}^2$

$=\underset{i}{\mathrm{\Σ}}{(1n{p}_2\left({\stackrel{\mathrm{\ρ}}{x}}_i\right|\mathrm{\α},\mathrm{\β})-{1\mathrm{np}}_1\left({\stackrel{\mathrm{\ρ}}{x}}_i\right))}^2$

由于没有样本时

可能为0和

取值小得不合理时，利用lsqcurvefit()函数可能得到不准确解；风险值R₂和参数α、β之间的关系可以通过观测获得。

Images