宇航学报 
JOURNAL OF ASTRONAUTICS
1998年 第19卷 第1期 No1 Vol.9 1998科技期刊

具有理论阈值的子波域滤波方法
潘　泉　戴冠中　张洪才
(西北工业大学自动控制系*西安*710072)
　　摘　要　本文在对多尺度滤波和子波域滤波综述的基础上,用奇异性和自相似性理论对信号与噪声的不同性质进行了分析,揭示了信号与噪声的多尺度数学机理,并对基于此机理的Mallat去噪声算法和子波域波波(WTDF)算法进行了研究,首次给出了WTDF的阈值选取方法和理论证明,完善和推广了WTDF,通过仿真计算表明,具有理论阈值的WTDF较之传统的试凑阈值WTDF性能有明显的改善。
　　主题词　子波域滤波　阈值　子波分解与重构
A THRESHOLD SELECTION METHOD 
FOR WAVELET TRANSFORM DOMAIN FILTER
Pan Quan　Dai Guanzhong　Zhang Hongcai
(Dept.of Automatic Control Northwestern Polytechnical University*Xi'an*710072)
　　Abstract　In this paper we describe the signal processing algorithms which use wavelet transform theory.Then a method to select threshold for wavelet transform domain filter(WTDF) is developed.The threshold selection method for WTDF brings the algorithm to completion.The simulation shows that the new algorithm is better than WTDF without the new threshold.
　　Key words　Wavelet transform　Multiscale analysis　Threshold
1　引　言
　　信号处理是子波变换最重要的工程应用领域,是子波理论发展的最重要工程动力之一。在早期的多尺度信号处理工作中,人们就已注意到信号和噪声在不同尺度上有不同的特征表现,并试图有效地利用这些特征,子波变换的出现为这一思想提供了一个自然而完善的工具［1,5,6］,使信号的多尺度处理技术向前推进了一大步。
　　自八十年代Morlet［1］将子波变换用于信号处理以来,这一领域的研究空前活跃。1992年Mallat［2］等人提出了基于信号奇异性(Singularity)的信号和图像多尺度边缘表示法及其基于此法的图像信号滤波方法。Mallat的方法是目前对信号和图像及噪声的多尺度数学特性刻画得最精细,数学基础最完善的理论之一。Mallat方法存在原理复杂,算法不易操作,编程较困难等不足,但其发展和应用潜力很大。Lu［3］等人直接将子波变换理论与传统的多尺度信号处理方法相结合,给出了一种性能优良的子波滤波器,并在医学图像处理领域得到应用。Lu的方法(子波域滤波器WTDF)对信号和图像去噪的同时,有相当好的保留信号图像边缘特征的能力,但WTDF未能给出多尺度滤波过程叠代求解时的阈值选取方法或选取原则,对滤波性能有所影响,算法在理论上尚欠完善。耶鲁大学的Coifman和斯坦福大学的Donoho另辟蹊径,提出了直接对信号和图像取阈值处理［4］。C-D算法同样没有给出阈值的选取方法和原则,关于C-D方法的研究我们另有专文报告。
　　本文将首先研究信号和噪声的多尺度数学机理(第2节),然后分析研究子波域滤波算法(第3节),最后给出WTDF阈值选取的证明和方法,并通过仿真计算比较了Lu与新算法的滤波性能。
2　信号和噪声的多尺度数学机理
　　在有背景干扰环境下检测信号,不仅与信号的形式和干扰性质有关,也与信号处理的方式本身密切相关。对不同类型的信号寻找最佳的处理技术一直是信号处理及检测的主要问题。
　　传统的谱分析技术,当信号是时间平稳,且有明显区别于噪声的谱特性时是十分有效的。然而实际中经常遇到非平稳信号的情况,由于这类信号的谱特性扩展得很宽,利用傅氏变换的基函数很难与其匹配。以子波变换为基础的处理方法在获得高的信噪比增益的同时,能够保持对信号突变信息的良好分辨,因此对非平稳信号处理有其自身的优越性。
　　首先我们通过奇异性和自相似性理论来分析信号与噪声的不同性质。
　　我们称无限次可导的函数是光滑的或没有奇异性的,若函数在某处有间断或某阶导数不连续,则称该函数在此处有奇异性。傅氏变换一直是研究函数奇异性的基本工具,由一函数的傅氏变换趋于零的快慢可以推断该函数是否有奇异性及奇异性的大小。但傅氏变换缺乏空间局部性,因而一函数的傅氏变换只能确定一函数奇异的整体性质,而难以确定其奇异点在空间的具体分布。
　　设实子波ψ(x)满足连续可微条件,并具有n阶消失矩(n∈N),f(x)∈L2(R),δx0为x0的任一开邻域,则x∈δx0,f(x)在x0处的奇异性为α等价为下述条件成立:
|Wf(s,x)|≤ksa
(1)
其中Wf(s,x)是函数f(*)在尺度s上的子波变换。(1)式表明,
|Wf(s,x)|=O(sa)
(2)
对于奇异点x0,由于其正则度小于邻域δx0内其余的Lipschitz指数,所以当s充分接近于零时,x0处的子波变换模值衰减最慢,从而δx0内的点x收敛到x0且成为模极大值。
　　上述特性可以用来区分噪声和信号。
　　设有突变信号f(x),x0为信号f(x)的局部突变点(奇异点),则在该点处f(x)的子波变换取得模极大值,(1)式变为
log2|W2jf(x)|≤log2k+α*j
(3)
其中W2jf(x)为f(x)的二进子波变换。若信号在x处奇异性大于0,由(3)式知,随尺度j的增大,子波变换模极大值的对数也变大。若信号有负的奇异性,则情况相反。负的奇异性意味着信号具有比不连续(有界,α=0)更差的奇异性,如白噪声的情况。对于白噪声可以证明它是一个几乎处处奇异的随机分布,且具有负的奇异指数α=--ε,ε＞0。
　　除上述的奇异特性以外,噪声的自相似性特性也是重要的特征。如对云图一类的自然信号,可用Browian运动描述的自相似分形信号描述。布朗运动F(x)的局部微发|F(x)-F(x+Δx)|/‖Δx‖H是一个高斯概率分布。这些随机过程是自相似的。即r＞0,F(x)与rHF(rH)是统计一致的,F(x)的一次实现在所有尺度和分辨率上看起来是相似的。
　　奇异性和自相似性是近年来极受瞩目的信号和噪声特性,这两种方法成为信号处理研究走向近代的代表性思路之一。
3　子波域滤波算法
　　Witkin首先将尺度空间方法和信号的子带分解方法结合,设计了一套复杂的由粗尺度到精尺度搜寻信号边缘的方法［2］。1992年Mallat使用子波变换而非传统的子带分解手段,利用信号和图像的边缘(即局部极大值点)来重构信号和图像,随后Mallat又将这种算法用于信号和图像的去噪处理［2］。Mallat算法的基本思路是:①对信号(图像)求其子波变换;②除了最粗尺度外,对其余的细节尺度信号求极大值;③利用信号和噪声奇异性之差异,去除噪声信号所对应的模值极大点;④利用Mallat重构算法——投影迭代法重构信号,得到滤波后的信号或图像。上述方法的最大缺点是算法复杂,且搜索和迭代投影将耗费太多的计算资源。
　　Rosenfield曾指出,在进行数字图像处理时,直接将相邻频带上数据相乘,可以有效且准确地确定信号边缘［3］。基于上述思想,将信号在不同尺度上的子波变换值按相邻相乘处理,利用阈值检验,从噪声中检出重要的信号边缘,并移去噪声,再重构信号。上述子波域滤波方法(Wavelet Transform Domain Filter——WTDF)［3］与Mallat算法相比稍嫌粗糙,但它具有直观性、易实现等特点。
　　WTDF的一个基本出发点是:尖锐信号边缘的子波变换值在子波变换的各尺度上都有较大幅值且幅值变化不大(Mallat也指出这一性质),但噪声的子波变换幅值随着尺度的增加而快速减小。直接将信号子波变换值在几个相邻尺度上进行相关运算(记为Corr1(m,n)),可以相当准确地确定信号边缘或其他重要特征的位置。

(4)
其中N表示离散信号的点数,n表示时间,m表示尺度,W表示信号的子波变换。若M是总的尺度数,应该有1＜M-m+1。直接将信号数据进行尺度空间相关,将锐化和增强主要边缘,同时噪声和小的边缘将受到压缩。当某个局部存在信号边缘或其他重要特征时,此法易将背景噪声除去。
　　WTDF的计算如下:取1=2,令Wf为滤波后子波值,初始化时全部元素置0。
　　①求取各尺度与相邻尺度的Corr2(m,n),得到增强的信号和变弱的噪声;
　　②对Corr2(m,n)重新归一划,使其归一化到W(m,n)的能量上去:

(5)

(6)
归一化后的相关值Corr2为:

(7)
　　③若|Corr2(m,n)|≥|W(m,n)|,则认为n点处子波变换是由信号引发,将W(m,n)赋给Wf的相应位置。并将W(m,n)置0;否则,认为是由噪声引发,W(m,n)保留;
　　④返回步骤①,重复②～③,直到PW(m)满足一个与噪声能量水平有关的阈值比th(m)为止。
　　在给出WTDF算法时,未能给出th(m)阈值的选择原则或方法,而我们的研究表明,阈值的选取对算法的理论完整性和性能有重要的影响。
4　WTDF的阈值选取方法
　　注意到WTDF已检出全部信号,W(m,n)只剩下噪声时,是噪声子波变换值方差的一个无偏估计,并且是最优的。因而可设在尺度m上噪声子波变换值的方差为阈值th(m)，并将阈值条件改为。下面推导噪声序列V(n)～N(0,δ2)在m尺度上的方差。
　　在实际子波计算中,一般将原始的离散数据序列V视为一连续函数在尺度2°上的离散子波值,即S(0,n)=V(n)。其中S(0,n)表示第0级尺度信号。
　　记W(m,n)的傅氏变换为W^(M,ω)=W^(2mω)。由离散子波分解和重构定理我们知:

(8)

(9)
其中H^是尺度函数的共轭滤波器(低通滤波器)的傅氏变换,为高通滤波器的傅氏变换。
　　由(8)、(9)式可得:

(10)
求上式功率谱的数学期望为:

(11)
可求出th(m)为

(12)
其中*表示卷积,Hi表示H的2i尺度膨胀。Gm-1表示G的2m-1尺度扩张。‖*‖F为范数。
5　仿真计算
　　图1～2给出了一个尖形信号(二维图像边缘的一个剖面)的WTDF滤波结果,图1为原始信号和加有高斯白噪声(水平为0.3)的信号,图2为滤波后的重构结果。滤波前后信噪比为:滤波前信噪比为11.97dB;用无理论阈值WTDF滤波,经过十几次的试凑,信噪比增益为6.54dB,用本文的理论阈值WTDF滤波,信噪比增益值达9.74dB。具有理论阈值的WTDF的性能改善是显而易见的。


图1　原始信号和受噪声污染的尖形信号


图2　WTDF滤波重构后的信号
6　结束语
　　在子波域滤波方法的研究中,我们完善和推广了WTDF算法,并用奇异性和自相似性理论对信号与噪声的不同性质进行了分析。在WTDF中,滤波的效果与阈值的选取密切相关,当选取的阈值与实际系统噪声方差偏离时,应当考虑对子波域滤波和阈值决策滤波的自适应能力,这对WTDF在处理多批次数据时能自动调整阈值大小也是十分必要的。
参考文献
1　Yver Meyer.Wavelets:algorithm and application.SIAM Philadelphia,1993
2　Mallat S.Singularity detection and processing with wavelets.IEEE IT 1992,(38)
3　J.Lu.Contrast enhancement of medical images using multiscale edge representation.Optical Engineering,1994,(33)7
4　D.Donoho.Denoising by soft-thresholding.IEEE Trans.on IT,1995,(41)3
5　张贤达.现代信号处理.清华大学出版社,1994
6　刘贵忠.小波分析及应用.西安电子科技大学出版社,1993
本文于1996年4月18日收到
　　*国防科学基金N0.94J 1A.4.100资助项目
