基于小波变换的图像去噪方法研究报告附MATLAB程序

基于小波变换的图像去噪方法的研究

（附送程序，见上传者“我的文档”）

摘要

图像在生成和传输的过程中会受到各种噪声的干扰,对信息的处理、传输和存储造成极大的影响。寻求一种既能有效地减小噪声,又能很好地保留图像边缘信息的方法,是人们一直追求的目标。小波分析是局部化时频分析,它用时域和频域联合表示信号的特征,是分析非平稳信号的有力工具。它通过伸缩、平移等运算功能对信号进行多尺度细化分析,能有效地从信号中提取信息。

本文对基于小波变换的图像去噪方法进行了研究分析, 详细介绍了其去噪原理和算法,分析了去噪过程中参数的选取问题,并给出了一些选取依据;详细介绍了小波系数相关性去噪方法的原理和算法。最后对均值滤波、中值滤波和维纳滤波方法在高斯噪声下进行了分析比较,并给出了仿真实验结果，结果证明小波去噪十分有效，其结果好于其它3种滤波。

关键词：小波变换，图像去噪，阈值，阈值函数

1．引言

数字图像在我们日常生活中起着非常重要的作用，它与我们的日常生活息息相关，例如在卫星、电视、核磁共振、计算机视觉、地球信息系统以及天文学中应用非常广泛。但是一般情况下采集到的数字图像是含有噪声的。噪声[1]可以理解为“妨碍人们感觉器官对所接收的信源信息理解的因素”。图像在生成和传输的过程中灰受到各种噪声的干扰，对信号的处理、传输和存储造成极大的影响。数字图像之所以含有噪声这是因为在图像的采集、获取、编码和传输的过程中，所有的图像均不同程度地被可见或不可见的噪声“污染”。对于这种“污染”，如果信噪比(SNR)低于一定水平，就会影响图像场景内容的表示，直接导致图像质

量的下降。除了视觉质量上下降外，噪声还可能掩盖一些重要的图像细节，使图像的熵增大，从而对于图像数据的有效压缩起到了一定的妨碍作用。对于图像在采集、获取过程造成的“污染”，我们虽然尽量提高硬件设备以获取质量更高的图像，但图像传感器的截止频率总是有一定的，受硬件水平和价格的限制，且图像在编码和传输过程中造成的“污染”，必需采取有效的降噪技术才能提高图像的质量。

对图像进行去噪最初主要是在空域内进行的，图像空域去噪方法很多，主要是通过各种滤波器对图像进行去噪。例如均值滤波器、顺序统计滤波器、维纳滤波器等。为了进一步提高去噪的效果，在变换域中进行降噪处理成为有效的方法，图像变换域去噪就是对图像进行某一种变换，然后将图像从时域变换到变换域中，再对变换域中的图像变换系数按照某种方法进行处理，最后再对处理后的系数按照某种方法进行反变换，这样就实现了将图像去除图像噪声的目的。将图像从时域转换到变换域的变换方法很多，例如傅立叶变换、小波变换等等。

小波变换是在短时傅立叶变换的基础上发展起来的一种新型的变换方法。小波变换具有多分辨率分析的特点，在时域、频域都具有较强的表征信号局部特征的能力，因此基于小波分析的图像去噪技术已成为图像去噪的一个重要方法。本文采用小波阈值去噪的方法，从去噪的效果上比较了多种去噪方法的优劣，实验证明小波去噪在图像噪声处理中起到很好的效果。

2．小波变换概述

2.1 小波变化去噪技术研究现状

上个世纪八十年代Mallet 提出了MRA(Multi_Resolution Analysis)，并首先把小波理论运用于信号和图像的分解与重构，利用小波变换模极大值原理进行信号的奇异性检测，提出了交替投影算法用于信号重构，为小波变换用于图像处理奠定了基础[1]。后来，人们根据信号与噪声在小波变换下模极大值在各尺度上的不同传播特性，提出了基于模极大值去噪的基本思想。1992年，Donoho和Johnstone 提出了“小波收缩”，它较传统的去噪方法效率更高。“小波收缩”被Donoho 和Johnstone证明是在极小化极大风险中最优的去噪方法，但在这种方法中最重要的就是确定阈值。1995年，Stanford大学的学者D.L.Donoho和I.M.Johnstone 提出了通过对小波系数进行非线性阈值处理来降低信号中的噪声[2]。从这之后的小波去噪方法也就转移到从阈值函数的选择或最优小波基的选择出发来提高去噪的效果。影响比较大的方法有以下这么几种：Eero P.Semoncelli和Edward H.Adelson提出的基于最大后验概率的贝叶斯估计准则确定小波阈值的方法[3]；Elwood T.Olsen等在处理断层图像时提出了三种基于小波相位的去噪方法：边缘跟踪法、局部相位方差阈值法以及尺度相位变动阈值法；学者Kozaitis结合小波变换和高阶统计量的特点提出了基于高阶统计量的小波阈值去噪方法[4]；

G.P.Nason等利用原图像和小波变换域中图像的相关性用GCV(general cross- validation)法对图像进行去噪；Hang.X和Woolsey等人提出结合维纳滤波器和小波阈值的方法对信号进行去噪处理[5]，Vasily Strela等人将一类新的特性良好的小波(约束对)应用于图像去噪的方法[6]；同时，在19世纪60年代发展的隐马尔科夫模型(Hidden Markov Model)，是通过对小波系数建立模型以得到不同的系数处理方法；后又有人提出了双变量模型方法[7]，它是利用观察相邻尺度间父系数与子系数的统计联合分布来选择一种与之匹配的二维概率密度函数。这些方法均取得了较好的效果，对小波去噪的理论和应用奠定了一定的基础。

小波具有低墒性、多分辨率、去相关性、选基灵活性等特点。

2.2 连续小波变换

1. 小波基函数

所谓小波(Wavelet)，即存在于一个较小区域的波。小波函数的数学定义是：

设()t ψ为一平方可积函数，即()()R L t 2∈ψ，若其傅立叶变换()w ψ

?满足： ()

∞=?

dw C R

w

w 2

ψψ（2.1）

时，则称()t ψ为一个基本小波或小波母函数，并称上式是小波函数的可容许条件。

根据小波函数的定义，小波函数一般在时域具有紧支集或近似紧支集，即函数的非零值定义域具有有限的范围，这即所谓“小”的特点;另一方面，根据可容许性条件可知()00==w w ψ，即直流分量为零，因此小波又具有正负交替的波动性。

将小波母函数()t ψ进行伸缩和平移，设其伸缩因子(亦称尺度因子)为a ，平移因子为b ，并记平移伸缩后的函数为()t b a ,ψ，则:

()()0;,,1

,≠∈=--a R b a a t t b a τ

ψψ（2.2）

并称()t b a ,ψ为参数a 和b 小波基函数。由于a 和b 均取连续变换的值，因此又称为连续小波基函数，它们是由同一母函数()t ψ经伸缩和平移后得到的一组函数系列。

2. 连续小波变换

将()R L 2空间的任意函数()t f 在小波基下进行展开，称其为函数()t f 的连续小波变换CWT ，变换式为：

()()()dt t f f b a WT a

b t R

a

b a f -?>=

=

,,,（2.3）

当小波的容许性条件成立时，其逆变换为：

()()()?

?

+∞∞

--+∞

∞-?=

db b a WT t f a b t f a da C ψψ

,2

1（2.4）

其中()

∞=?

dw C R

w

w 2

ψψ为()t ψ的容许性条件。

另外，在小波变换过程中必须保持能量成比例，即:

()()???

=R

R

f R

a da

dx x f C db b a WT 2

2

,2

ψ（2.5）

由CWT 的定义可知，小波变换和傅立叶变换一样，也是一种积分变换，其中

()b a W T f ,为小波变换系数。

可见小波变换对函数()t f 在小波基上的展开具有多分辨率的特性，这种特性正是通过缩放因子a 和平移因子b 来得到的。

一个一维函数()t f 的连续小波变换是一双变量的函数，变量比()t f 多一个，因此称连续小波变换是超完备的，因为它要求的存储量和它代表的信息量都显著增加了。对于变量超过一个的函数来说，这个变换的维数也将增加。

若()t f 是一个二维函数，则它的连续小波变换是：

()()(

)d xdy t t f b b a WT b y b x a

y x f y x

--+∞∞-+∞

∞

-??

=

,

,,,211

ψ

（2.6）

其中，x b ，y b 表示在两个维度上的平移，二维连续小波逆变换为：

()()(

)y

x

a

b y a

b x y x f a da

C db db b b a WT y x f y x

??

?

+∞∞---+∞∞

-+∞

=

,

,,,0

13ψ

ψ

（2.7）

同样的方法可以推广到两个或两个以上的变量函数上。

2.3 离散小波变换

计算机中的图像信息是以离散信号形式存放的，所以需要将连续小波变换离散化。而最基本的离散化方法就是二进制离散，一般将这种经过离散化的小波及其变换叫做二进小波和二进变换。需要注意的是这里的离散化都是针对连续的尺度因子a 和连续平移因子b 的，而不是针对时间t 的。这儿限制尺度因子a 总是正数。

1. 尺度与位移的离散化

对连续小波基函数()t b a ,ψ尺度因子a 和平移因子b 进行离散化可以得到离散小波变换()b a W T f ,，从而减少小波变换系数的冗余度。在离散化时通常对尺度因

子a 和平移因子b 按幂级数进行离散化，即取m

m b b a a 0

0,==（m 为整数，,10≠a 但一般都假定10 a ），得到离散小波函数为：

()(

)()0011

,0

00

nb t a t m a a b na t a n m m

m -=

=

--ψψ

ψ（2.8）

其对应系数为：

()()()dt t t f t f C n m n m n m ,,,,ψψ?+∞

∞

->==<（2.9）

2. 二进制小波变换

二进小波变换是一种特殊的离散小波变换，特别地令参数20=a ，10=b ，则有()n t m n m m -=--22,ψψ。该二进尺度分解的原理在二十世纪三十年代由

Littlewood 和 Paley 在数学上进行了研究证明。

离散小波变换为：

()()()dt t t f n m n m WT n m f ?+∞

∞

->==<,,,ψ（2.10）

离散二进小波变换为：

()()()dt t t f n m n m WT n m f ?+∞

∞

->==<,,,ψ（2.11）

二维离散小波变换：

我们考虑二维尺度函数是可分离的情况，也就是：

()()()2121,x x x x φφφ?=（2.12）

设()i x ψ是与()i x φ对应的一维小波函数，则有：

()()()21211,x x x x ψφψ=（2.13）

()()()21212,x x x x φψψ=（2.14） ()()()21213,x x x x ψψψ=（2.15）

以上三式就建立了二维小波变换的基础。

3. 图像的小波变换及Mallat 算法

图像是二维信号，二维多分辨率分析与一维类似，但空间)(2R L 变成)*(2R R L ，一维中引入的尺度函数)(x φ变为),(y x φ。

设{

}Z j j V ∈是)(2R L 的一个多分辨率分析，则可以证明，张量空间{}

Z

j j V ∈2：

j j j V V V ⊕=2构成)*(2R R L 的一个多分辨率分析，并且二维多分辨率分析{}

Z

j j V ∈2

的二维尺度函数),(y x φ为

)()(),(y x y x φφφ= (2.16)

式中：)(x φ是{

}Z j j V ∈尺度函数(一维)。 式(2.28)说明了二维尺度函数的可分离性。对于每一个Z j ∈,函数系

}),()()(),({2,,,,Z m n y x y x m j n j m n j ∈=φφφ构成{}

Z

j j V ∈2

的规范正交基，这里下标j ，n ，

m 的含义是：

)0(),2()2(2),(,,≥--=j m y n x y x j j j m n j φφφ（2.17）

我们将{}

Z

j j V ∈2

称为)*(2

R R L 的可分离多分辨率分析。因)(x φ、)(y φ都是低通的尺度函数，所以{}

Z

j j V ∈2

平滑的低通空间。

如果)(x ψ是一维多分辨率分析{

}Z j j V ∈的正交小波基，则二维多分辨率分析的三个小波函数为：

)

()(),()()(),()

()(),(321y x y x y x y x y x y x ψψψφψψψφψ===（2.18） 对于每一个Z j ∈，它们的整数平移系为：

)

()(),()()(),()

()(),(,,,,3,,,,2,,,,1y x y x y x y x y x y x m j n j m n j m j n j m n j m j n j m n j ψψψφψψψφψ===（2.19） 注意这里的上标只是索引而不是指数。它们构成了Z j j W ∈}{2

的规范正交基。

因此以上的三个正交基中都至少包含一个带通的)(x ψ或)(y ψ，所以它们都是带通的。也就是说这三部分反应的都是细节信息。

具体来说，函数系}2,2(2{)},({,,m y n x y x j j j m n j --=εεψψ，3,2,1,0=≥εj 是

)*(2R R L 的正交归一基，其中均为整数，ε=1,2,3分别对应于水平、垂直和对角三个方向。对于任一二维图像信号)*(),(2R R L y x f ∈，在分辨率j -2下有：

2

,,332,,222

,,112,,),(),(),,(),(,),(),,(),(),(),,(),(,),(),,(Z m n y x y x f f D Z m n y x y x f f D Z m n y x y x f f D Z m n y x y x f f A m n j j m n j j m n j j m n j j ∈=∈=∈=∈=ψψψφ （2.20）

上式表明，在j -2分辨率上将图像分解成f A j 、f D j 1、f D j 2和f D j 3四个子图，其中f A j 代表原图像在j -2分辨率上的近似（即图像的低频部分，不妨用LL 来表示），f D j ε则代表这种近似的误差（即图像的高频部分或“细节”部分）；f

D j 1

对应于垂直方向的高频成分，即水平的边缘（细节）信息（不妨用LH 表示）；f D j 2对应于水平方向的高频成分，即垂直的边缘（细节）信息（不妨用HL 表示）；f D j 3则对应于对角方向的高频成分（不妨用HH 表示）。图2.1形象地表示了二维图像的多分辨率小波分解。图中符号的上标表示图像的小波分解层数，图中示意了图像的2级小波分解。可以看到，在每一分解层上，图像均被分解为LL, LH, HL 和HH 的四个频带；下一层的分解仅对低频分量LL 进行分解。

图2.1二维图像的小波分解

按照Mallat 的快速算法，图像的小波分解算法如图2.2所示：

图2.2图像的小波分解算法

图2.2示意了图像的一步小波分解过程，可以看到：二维图像的小波分解可以对图像依次按行、按列与一维的低通(H)和高通(G)滤波器作卷积来实现，在卷积之后进行相应的降2采样。图像小波分解的重构算法如图2.3所示。

图2.3图像的小波重构算法

图2.3示意了小波分解图像的一步重构过程，小波分解图像的重构是先对列或行进行升2采样(在相邻列或行间插入一零列或零行)，然后再按行、按列与一维的低通或高通滤波器进行卷积，这样递推下去便可重构原图像。二维图像的这种行、列可分离性简化了图像的小波变换。

3．小波阈值去噪

3.1阈值函数的选取

在阈值去噪中，阈值函数体现了对超过和低于阈值的小波系数模的不同处理策略以及不同估计方法。常用的阈值函数有硬阈值函数和软阈值函数两种[8-9]。硬阈值策略保留大于闽值的小波系数，而把小于阈值的小波系数都设定为零。软阈值策略把小于闽值的小波系数置零，把大于阈值的小波系数的绝对值减去阈值以除去噪声的影响。设A 为元素小波系数，t 为阈值，两种阈值函数的表达式分别如下：

（1） 硬阈值函数：

()????

?≥==t

A t

A A t A T x h 0,（3.1）

（2） 软阈值函数：

()????

?≥-==t

A t A t

A t A T x s 0

,（3.2）

硬阈值法得到的小波系数的连续性较差，重构信号可能出现突变或振荡现象；如图3.2（a ）所示。软阈值法的到的小波系数的连续性好，但当小波系数较大时，得到的处理后的小波系数和实际的小波系数有一定的偏差，会导致重构结果的误差。如图3.2（b ）所示。

图3.2 两种阈值方法

3.2阈值的选取

小波阈值去噪方法除了阈值函数的选取，另一个关键因素是对阈值的具体估计。如果阈值较小，去噪后的图像信号与输入比较接近，但是残留了较多噪声。若阈值较大，则得到较多为零的小波系数，对于软阈值策略重建图像变得模糊。在小波域阈值去噪中，阈值的选取直接影响滤波效果。本文用了2种的阈值估计方法。

（1） Visushrink 阈值

也称通用阈值，是Donoho 和Johnstone 提出的，针对多维独立正态变量联合分布，在维数趋于无穷是得出的结论，是基于最小最大估计得出的最优阈值[10]。阈值t 的选择满足

A

A

（a ）硬阈值

（b ）软阈值

()N l t n n 2σ=（3.3）

其中，N 为信号长度，n σ为噪声标准差。n σ的估计公式如下：

6745

.0)(ij n Y Median =

σ（3.4）

（2） Surehrink 阈值

也称Stein 无偏风险阈值，是一种基于Stein 的无偏似然估计的自适应阈值选择[s53]，它是针对软阈值函数得出的结论。SureShrink 阈值的具体计算过程如下：

①求取信号的长度N ；

②将某一层的小波系数的绝对值由小到大排列，得打一个新的向量，X=[x 1,x 2,x 3,…,x N ]，其中x 1

③计算风险向量，R=[r 1,r 2,r 3,…,r N ]，其中

N

x x N i N r i

k i

i i ∑=+-+-=

1)1(2（3.5）

④以R 中最小元素r B 作为风险值，由r B 的相应位置B 求出对应的x B ，则SureShrink 阈值为B n x t σ=。

3.3 小波阈值去噪过程

阈值去噪法就是通过对图像进行小波变换，得到小波变换系数。因为信号对应的小波系数包含有重要的信息，其数据较少，幅值变化较大，而噪声对应的小波系数的分布则恰好相反，通过设定特定的阈值对小波系数进行取舍，就可以得到小波系数估计值，最后通过估计小波系数进行小波重构，就得到去噪后的图像[11]。阈值去噪法实现简单，计算量小，在实际中有着广泛的应用。经过阈值处理后，得到的处理后的小波系数多，因此可以直接对其进行小波重构，如图3.3所示。

图3.3小波阈值去噪过程

在实际中，小波分解层数值J一般取3-5。实际中，J取越大，则噪声和信号表现不同特征越明显，越有利于信噪分离；另一方面，对重构来讲，分解层越多，则失真越大，即重构误差越大。最大分解层数应该与信号的信噪比有关，根据实验表明，对一般的信号而言，若信噪比大于20，则取J=3，否则，取J=4。

本文采用了小波基为coif2的小波分解法，以及以上2中阈值去噪法，过程如下：Array图3.4本文小波去噪流程图

4．实验仿真及结论

本文为说明小波去噪的有效性和优越性，对含有高斯白噪声的Lenna图像进

行消噪处理，其中噪声平均为0，方差2σ为0.01，0.02，0.03，0.04，0.05，0.06。

为了客观评价噪声图像的噪声污染度，以及各种去噪方法的效果，本文引用了峰值信噪比PSNR ：

∑∑-=--=101

22

10

)),('),((1255log 10m i n j

j i I j i I mn PSNR （4.1）

式中I ’(i,j)和I(i,j)分别是有噪图像和原始理想图像在点（i,j ）处的灰度值，m 、n 分别表示图像的行数和列数。

由于在各方差噪声下，噪声图像的峰值信噪比均小于20，所以小波尺度取J =4。本文对方差为0.01噪声图像进行小波系数分解，得到下图：

图4.1小波系数分解图

本文分别用小波去噪、均值去噪、中值去噪和维纳去噪对8位的255灰度值的Lena 图像进行去噪处理，结果图如下所示：

图4.2去噪效果图

此外，本文也对多个方差的噪声进行去噪处理，并得出PSNR 值进行对比，如下表所示。

表4.1 各种去噪方法的PSNR 值比较

通过以上结果分析可知，小波阈值去噪的效果较均值滤波、中值滤波和维纳滤波效果都好。

参考文献

[1] 杨福生.小波变换的工程分析与应用.北京:科学出版社,1999：150-176.

原图像噪声图小波去噪图

中值滤波均值滤波维纳滤波

[2] 阮秋琦.数字图像处理学[M].第二版.北京.电子工业出版社.2007.1-2

[3] D.L.Donoho, I.M.Johnstone. Wavelet Shrinkage: Asymptopia J. R. Stat. Soc. B.1995,

57：301-369.

[4] D.L.Donoho, I.M.Johnstone. Ideal Time-Frequency Denoising. Standford

University,Technical Report, Dept. of Statistics,1994：3397-3415.

[5] 张晔，黄秀明.小波变换及在图像处理中的小波特性分析[J].中国图像图形学报，

1997,2(7):480-484.

[6] Elwood T.Olsen, Biquan Lin. A Wavelet Phase Filter for Emission Tomography.

SPIE,1996，2491：829-839.

[7] Kozaitis, Basuhail. Adaptive Wavelet Threshold Selection Using Higher Order

Statistics for Signal Denoising. SPIE,1998，3391：68-74.

[8] Tu Dan and Shen Jianjun.The design of wavelet domain wiener filter and its

application in inage denoising.Systems Engineering and Electronics.2011.23.4-7 [9] Shubhankar R,Bani K.A Bayesian transformation model for wavelet shrinkage. IEEE

Trans. on IP.2003.12(12).1512-1520

[10] Donoho D L,Johnstone I.Adapting to unknown smoothness via wavelet shinkage.

Journal of Amercan Stat.Assoc.1995.90.1201-1225

[11] 张磊，潘泉，张洪才，戴冠中，小波域滤波阈值参数c的选取，电子学报，

2001，3，29(3).

[12] Xu Chen, Zhao Ruizhen, Gan Xiaobing. ApplicationAlgorithm of Wavelet

Analysis[M]. Beijing：Science Press，2004.

[13] Chao Rui, Zhang Ke, Li Yanjun. A Wavelet TransformBased Image Fusion

Method[J].Acta Electronica Sinica，2004，32(5)：750-753.

[14] Gai Li-ping,Wang Gui-lian. Technology of filtering in medical image

processing[J].ChineseMedical EquipmentJournal,2007,(06):50-51.

[15] Gao Qing-wei,Li Bin. An image de-noising method based on stationary wavelet

transform[J].Journal of Computer Research and Development,2002,(12):1689-1693.