多尺度变换（通用7篇）
多尺度变换 篇1
现实中的图像,无论是在在采集、转换和传输中,都难以避免收到成像设备或者外部环境噪声干扰等影响,降低了图像的细腻行,因此很多现实图像都是被噪声污染的图像,这些噪声对后续的图像分析、图像理解造成很大的影响,因此对图像进行降噪是整个图像处理过程中首要任务。图像处理主要分空域处理和频域处理,而频域处理具有空域处理无法比拟的优越性,因而在实际应用得更广泛。早期的频域处理主要傅里叶变换、沃什变换、小波变换等,由于小波变换具有多分辨率的特性,可以对信号进行局部分析,因此,当前主流的频域降噪算法都是基于小波变换。但是小波变换只能最优的逼近一维信号,二维小波只是对一维简单的张成,方向性选择比较差,不能有效的表示二维信号的轮廓信息,而图像信号边缘以及纹理特征往往比较复杂,各个方向都可能表现异相,二维小波在图像降噪的过程中,会消除图像在纹理以及边缘细节信息,造成模糊。2002年M.N.Do和Vetterli M.提出了一种“真正”的二维图像稀疏表达方法—Contourlet变换,这种变换能够很好地表征二维信号的各向异性特征,具有多分辨特性、局域性、方向性的优点。由于Contourlet变换比小波变换更稀疏的表示方式,更好地表示图像的细节,在图像降噪的过程中能更好地保护图像的纹理以及边缘,使得降噪后的图像在细节方面更突出。
1 Contourlet 变换
Contourlet变换先进行多分辨率分解然后在进行多方向分解,首先用拉普拉斯金字塔(LP,Laplacian Pyramid)变换进行多尺度分解,并捕捉图像信号的奇异点,接着使用方向性滤波器组(DFB,Directional Filter Bank)进行多方向分解,将同个方向上的奇异点合成,捕捉二维信号的高频成分。图1给出了Contourlet变换的流程图。LP变换能够有效的避免由方向性滤波器所造成的低频分量的“泄漏”。
Contourlet变换的核心LP和DFB,因此Contourlet变换也称“塔式方向滤波器组”。而Contourlet变换的最终是用类似于轮廓段 (contour segment)的基结构来逼近表示二维图像信息。图2表示一个Contourlet变换可能分解的频率。由于拉普拉斯金字塔和方向性滤波器组具备的完全重构特性,因此其组合PDFB也具备能完全重构。但是由于拉普拉斯金字塔本身具有冗余性,因此Con-tourlet变换具有冗余性。
图像经Contourlet分解之后绝大多数系数幅值接近零,系数表现相当稀疏,而在原图像的边缘细节以及纹理特征部分稀疏幅值比较大,且在他们在各尺度上具有有一定的相关性。图3给出了lena图的Contourlet分解的例子,图像经3级LP分解,以及8个方向的分解,从图我们可以看到最大多系数都接近零,而在图像的轮廓和边缘部分系数的幅值比较大。
2 基于层结构的多阈值Contourlet降噪算法
阈值降噪是图像降噪处理中最常用的方法,其关键是选取一个适当的阈值把图像噪声系数和信号系数区分开来,然后对噪声系数进行置零或者修改其值,随后通过反变换重构图像,达到图像降噪目的。
文献[8]中提出了基于层结构的Contourlet多阈值图像降噪算法。文献中提出“算法使用Contourlet变换来代替小波变换,Con-tourlet的基函数能有效地对分段光滑的线段进行表示,而对奇异点则影响不大”,在降噪处理的过程中,即使将噪声系数误认为图像细节系数给予保留,或者相反,将图像细节系数误认噪声给予删除,也不会使降噪重构的图像出现显著数值较大的孤立点,从而造成图像边缘细节或纹理特征稀疏的丢失。其方法如下:
1)首先,对带噪声的图像进行CT变换,确定对图像进行分解的级数J和各个级数J上所要分解的方向数s。
2)根据下面式子计算出各个尺度和各个方向上噪声方差估计值,
wgi,j(s) 为带噪声图像高频系数。
3)对图像CT变换后的系数方差进行估计
4) 阈值调整:对低分解层的高频系数进行调整。
5) Contourlet反变换:将处理后的Contourlet变换系数进行反变换,可以重构经过降噪之后的图像。
3 基于Contourlet的相关尺度降噪算法
文献[10]提出了一种基于Contourlet变换的相关降噪新算法,其思想主要是将图像的Contourlet系数跟含噪声图像的Contourlet变换系数进行比较,随着分辨率越高,图像的细节部分对应的Contourlet变换系数越稳定,而噪声则衰减越大,即代表图像纹理特征或边缘细节的Contourlet系数相关比较强,而代表噪声的Contourlet变换系数相关比较弱或不相关。采用了利用小波系数相关量来确定图像边缘信息,并且与阈值函数相关结合。其过程如下:
4 基于Contourlet相关尺度多阈值降噪算法
采用硬阈值函数进行降噪处理后的信号虽然可以比较好地保留图像边缘纹理等细节信息,都由于这种理想滤波器会造成图像出现振铃或者伪吉布斯效应等视觉失真,而基于相关尺度的降噪方法,有时候易将过多噪声系数误判为有用信号,为了取得比较好的降噪效果,该文对以上两种方法进行综合,采用了相关尺度多阈值降噪方法,将阈值函数与相关尺度结合起来。在基于层机构的阈值基础上,对于判断为噪声的Contourlet系数进行相关量的计算,具体算法如下:
1)首先对带有噪声的图像进行CT变换,一般情况下确定对图像进行分解的级数J为3级和各个级数J上所要分解的方向数s为8。
2)运用文献[8]的方法,计算3个尺度以及8个方向上噪声方差估计值,并对对图像的Contourlet系数的方差进行估计最后得到阈值函数:
5 结果及分析
本文采用标准的大小为512*512像素测试图象Lena, Peppers为实验图象,进行仿真实验,比较了多阈值降噪、相关降噪和本文提出的多阈值相关降噪效果的比较,分解级数为4,4,8;表1为各个降噪方法后的信噪比(SNR),图4是lena图各种方法降噪后的结果。
从表1可以看出虽然尺度相关多阈值降噪的SNR不是最高,但在局部图上可以看出此方法在边缘细节的保持上效果更好。
6 结论
在Contourlet变换多阈值和相关性降噪的基础上,给出了一种结合尺度相关多阈值降噪算法,利用了不同尺度间Contourlet系数的相关性以及噪声信号的不相关性,保留了更完整的图像边缘细节信息,与层机构多阈值方法相结合,在图像降噪的过程中能更好地保护图像的纹理以及边缘,使得降噪后的图像在细节方面更突出。对各种图像的实验结果也表明了这种方法在降噪上有效可取。
多尺度变换 篇2
近年来,计算机视觉技术不断发展,临床医学对疾病诊断的准确性要求越来越高,三维医学图像可以为医生提供准确的特征参数及病理数据,是计算机辅助诊断的基础。医学图像边缘检测是三维重建的第一步,能够准确快速提取图像边缘是众多学者的研究热点。对图像区域构造边缘检测算子[1]是经典的边缘检测算法原理,常见边缘检测算子有Roberts算子、Prewitt算子、Sobel算子、Canny算子等。这些都是基于导数的边缘算子,采用2×2、3×3的模板为核与原始图像进行卷积和运算,然后选取合适的阈值作为边缘提取的依据。其中,Roberts算子采用2×2模板计算灰度函数对x和y的偏导数,检测速度快,算法简单,但易丢失细节,且对噪声敏感。Prewitt算子采用3×3模板计算灰度函数对x和y方向的偏导数,算法简单,有方向性,但是检测不精确,间断点较多。Marr[2]等提出的LOG算子,用高斯滤波器对图像进行平滑,再用Laplace算子计算二阶导数[3],对图像进行平滑处理,有效减少了噪声对Laplace求导的影响,但LOG算子检测方向性效果差。Canny算子同样采用高斯滤波器平滑图像,与LOG算子不同的是Canny算子在对图像进行高斯平滑、梯度计算后求出梯度幅值的非极大值抑制,双阈值检测并确定图像边缘。Canny算子是一个具有滤波、增强和检测功能的优化算子[3],边缘检测精度高,误判率低,但计算复杂,对噪声敏感。
小波变换[4]在图像边缘检测中由于其定位准确和噪声抑制能力强受到国内外学者的广泛重视。1992年Mallat利用二阶B样条小波变换实现了多尺度边缘提取,为小波边缘提取奠定了基础[5]。陆艳飞等[6]提出基于B样条二进小波的多尺度边缘检测算法,能有效地做平滑处理和抑制伪边缘;王森等[7]根据小波多分辨率特性进行多尺度相似矩阵重构,运算速度快,边缘提取精度高。基于小波多尺度分析的医学图像边缘检测,选择高斯函数[8]一阶导数为小波函数,结合数学形态学[9]知识,对医学图像进行边缘提取,并同经典边缘提取算法进行比较,实验结果表明,该算法的性能要优于经典边缘检测算法的性能,可准确、有效地检测图像边缘。
1 小波多尺度变换边缘检测原理
医学图像边缘检测有如下标准[10]:(1)尽可能减少边缘检测漏检率;(2)检测算法抗噪声干扰能力强,受边缘方向影响较小;(3)检测算法简单快速,精确度高。
实际上,图像处理和分析都是面向某种具体应用,上述评价标准间互相矛盾,难以兼顾。医学图像受噪声干扰比较严重,如果边缘检测不准确,会对后期三维重建等图像处理和分析工作造成很大影响。本文以脑部、腹部CT图像为例,采用多尺度小波变换法进行边缘检测,并同经典边缘检测算法结果进行比较。实验结果表明:该算法检测边缘连续性好,抗噪性好,边缘定位准确。本文主要从多尺度小波算子、小波函数的选取以及尺度的选择3个方面提高边缘提取的精度并减少噪声干扰。
1.1 小波函数的选取
在图像的边缘检测中,小波函数种类很多,也可以自己定义,小波函数直接关系到检测结果的好坏,小波变换相当于对尺度确定图像进行带通滤波。高斯函数的一阶导数和中心B样条函数导数是2种常见、高效的小波函数,中心B样条函数导数在计算的快速性和有效性方面远远优于高斯函数的一阶导数截尾造成的影响,但是边缘提取效果较差,故本文选择高斯函数一阶导数作为小波函数。一个Gaussian正态分布表达式如公式(1)所示:
可以证明,该函数满足下列条件:
1.2 多尺度小波算子
如1.1所述,θ(x,y)满足,是二维平滑函数。本文用到多尺度小波算子数学模型如下:
其中,公式(1)为高斯平滑函数θ(x,y)表达式;f(x,y)为待处理图像;公式(2)为高斯函数满足的积分条件;公式(3)、(4)为2个基本小波,分别为高斯平滑函数θ(x,y)在x和y方向的偏导数;公式(5)、(6)、(7)为引入尺度参数后公式(1)、(3)、(4)的变形;公式(8)、(9)为原始图像经高斯变换后沿x和y方向的小波变换表达式;公式(10)为原始图像沿x和y方向的小波变换矢量表达式;公式(11)、(12)分别为图像在不同尺度小波变换下的幅值和相角表达式;公式(13)为多尺度表达式,s为尺度参数。由以上推导公式可知,图像边缘为沿Asf(x,y)方向的Msf(x,y)极大值点的集合。
1.3 尺度参数的选择
在小波变换分析中,小尺度算子对噪声敏感但定位准确;大尺度算子会丢失部分细节,定位不准但抗噪声干扰能力强,有助于对大结构弱边界的检测。因此,滤波尺度参数s的选择取决于当前像素所处的区域是平滑区还是边缘区。若是边缘区,则选取最小尺度smin以精确定位图像边缘;若是平滑区,则选取最大尺度smax以减少边缘检测噪声。若无法确定,选取尺度smin+i,i的值根据实际图像噪声大小而定,本文i可选范围为i=1,2,3。
2 算例仿真
本文以2幅CT图像为例在MATLAB平台上仿真,一幅是脑部CT(960像素×720像素),第二幅是腹部CT(980像素×827像素)。算例中对这2幅CT图像分别采用多尺度小波变换和经典边缘检测算法进行边缘检测。
2.1 多尺度小波变换程序流程
本文算法的程序流程图如图1所示。
2.2 实验结果及分析
采用本文的检测算法和传统的边缘检测算子分别对2幅CT图像边缘进行检测,用MATLAB进行实验仿真,实验效果如图2~5所示。图2、3为脑部CT图像边缘检测结果,图4、5为腹部CT图像边缘检测结果。
从上述实验结果中不难发现:(1)如图2、4所示,Roberts算子边缘定位不准,丢失大量细节,检测效果较差;Soble算子、Prewitt算子检测质量有所改善,但边缘点不连续。此3种算子对噪声比较敏感,因此Roberts算子、Soble算子、Prewitt算子仅适合含噪声较少且简单的图像。(2)如图2、4所示,Canny算子在抗噪声干扰和精确定位之间提供较佳的折中方法,Canny算子对急剧变化方向上的边缘特别敏感,但计算复杂,存在着大量的噪声,而且边缘连续性及光滑性也很不理想;LOG算子较Canny算子受噪声干扰较小,但检测边缘点不连续,检测效果不佳。(3)图3、5采用多尺度小波变换边缘检测,分别显示了s=2、4、8的边缘检测结果。基于小波变化的边缘检测算法边缘定位准确、细腻,没有产生伪边缘,且抗噪性能优势明显。但随着尺度的增大,图像的主要轮廓渐显现出来,但边缘模糊,边缘细节丢失。
3 结语
医学图像边缘检测在三维重构中至关重要,本文针对医学图像三维重构中的边缘检测问题,利用小波变换的特点,采用基于多尺度的医学图像边缘检测算法,建立数学模型,选择合适的小波函数和尺度参数后以脑部、腹部CT图像为例,进行边缘检测。实验结果表明,该方法简单实用,检测效果明显优于已有的传统边缘检测算法,得到的边缘图像保留丰富的细节,抗噪性能好,边缘连续的同时也保证了边缘定位的准确性。边缘检测效果明显优于经典的边缘检测算子。
摘要：目的 :提出基于多尺度小波变换的边缘检测算法,能够准确有效地解决医学图像三维重构中的边缘提取问题。方法:结合小波变换理论,在经典边缘检测算法基础上采用一种多尺度小波变换法对脑部、腹部CT图像进行边缘提取,建立数学模型,并同经典边缘检测算法结果进行比较。结果:该算法检测边缘连续性好,抗噪性好,边缘定位准确。结论:该算法性能优于经典边缘检测算法的性能,可实现快速、准确的医学图像边缘检测。
关键词：小波变换,多尺度变换,医学图像,边缘检测
参考文献
[1]王钧铭,赵力.一种基于数学形态学的车牌图像分割方法[J].电视技术,2007,31(10):84-86.
[2]MARR D,HIDRETH E.Theory of edge detection[J].Processings of the Royal Society,1980:187-217.
[3]蔡胜仓.基于数学形态学的边缘检测及其在医学图像处理中的应用[D].青岛:青岛大学,2007:35-47.
[4]FENG X C,GAN X B,SONG G X.Wavelet theory and functional analysis[M].Xi′an,China:Electron Technological University Press,2003.
[5]MALLAT S,HWANG W L.Singularity detection and processing with wavelets[J].IEEE Trans on Information Theory,1992,38(2):617-643.
[6]陆艳飞,汪艳丽.基于B-样条二进小波的多尺度边缘检测[J].新疆师范大学学报(自然科学版),2015,34(2):58-61.
[7]王森,伍星,刘韬,等.基于二进小波变换的多尺度图切割方法[J].计算机工程与应用,2015,51(13):9-14.
[8]刘飞飞,董绍华.基于小波多尺度分析的图像边缘检测[J].图像处理,2009,25(3):271-274.
[9]林雯.基于数学形态学的线粒体电镜图像边缘检测算法[J].九江学院学报(自然科学版),2012,4(4):30-34.
多尺度变换 篇3
多尺度小波变换图像边缘检测算法存在的问题是大尺度下图像边缘细节信息会有丢失,导致检测到边缘点定位不准确而偏离实际边缘;小尺度下图像边缘点又易受噪声影响导致检测到伪边缘[2]。选取合适的尺度并采用合理的边缘细化算法是平衡多尺度下边缘检测优/缺点的主要思路。笔者在多尺度小波变换得到的边缘灰度图像中选择最优尺度灰度图像,再选择合适的灰度门限,将图像转换为二值图像后,用边缘细化算法将边缘准确定位到单像素宽度,以期较传统边缘检测算子在噪声容限和单像素边缘定位方面有所提高。
1 多尺度边缘检测①
对于一幅普通图像而言,边缘就是图像颜色或灰度发生突变之处,边缘检测的目的就是检测出图像中突变发生的位置。传统的边缘检测算子对图像进行边缘检测的基本原理就是求微分运算,找到图像中的突变点作为图像边缘,原理简单、操作容易,但是存在的最大问题就是对噪声敏感。
较好地解决噪声敏感问题的方法之一是对图像进行多尺度边缘检测,这需要首先对图像进行多尺度表达。对图像进行小波变换是获得图像多尺度表达的一种方法,小波变换能够把图像分解成多种尺度成分,并对大小不同的尺度成分采用相应的时域或空域取样步长,从而不断地聚焦到图像的任意微小细节、间断点、奇异点和边缘[3]。小波变换所独有的这种多尺度特性恰好可以用于多尺度下图像的边缘检测[4]。人在观察目标时,随着分辨率的增加,越走近目标观察获得的信息就越丰富;反之,获得的信息越少[5]。在噪声抑制与图像细节表现之间寻求平衡是多尺度分析的主要任务。由于小波变换在各尺度上都提供了图像的边缘信息,所以称为多尺度边缘,小波变换是提取多尺度边缘的基础。
2 图像小波变换
小波变换是将一个函数表示为一个尺度成分和这个尺度下的一个小的波动成分的一种运算。
2.1 连续小波变换
θ(t)是一个平滑函数,,小波ψ(t)是θ(t)的一阶导数,即,s(s>0)为尺度参数。待分析信号f(t)∈L2(R)在尺度s下的连续小波变换为:
在连续小波变换中,ψ为实函数,小波变换的模极大值|Wf(s,u)|就是f(t)经磨光后的函数的一阶导数的极大值,对应信号f(t)的突变点。
在实际问题中,经常处理的是离散数字图像。计算离散二维小波变换的常用方法是二进小波变换。适当选择小波及其滤波器,对数字图像信号直接利用离散快速二进小波变换算法求出各尺度下的小波系数。这种方法对所用小波有较大限制,但是计算速度快,是边缘检测的常用方法[6]。二维离散小波变换可以通过滤波器的卷积运算实现,设θ(n,m)是一个二维平滑实函数,该函数满足
其中W12jf(n,m)和W22jf(n,m)分别是尺度2j下图像f(n,m)的水平方向与垂直方向的小波变化系数(高频信息)。以f(n,m)二维离散小波的小波变换系数计算每一点(n,m)的模值:
对每个像素点(n,m)计算相角Af(2j,n,m)的正切值:
2.2 求边界点
确定阈值T>0,对于Mf(2j,n,m),如果Mf(2j,n,m)≥T,Mf(2j,n,m)取得局部最大值,即此时的(n,m)为模极大值点。因为梯度的局部极大值对应着图像f(n,m)的锐变处,而且图像的边缘也处于图像的锐变处,所以图像f(n,m)的梯度局部极大值点就对应着图像的边缘点[7]。
局部梯度幅值最大就是图像的局部高频信息,图像的噪声也是局部高频信息,经过小波变换后也可能产生边缘,这种边缘为伪边缘。图像进行小波变换后无论在哪个尺度上,尖锐边缘都有很大的信号值;相反,噪声点的信号值会随着尺度的增大而衰减[8]。
从图1所示的实验结果可以看出:随着尺度的增大,图像的细节逐渐减少消失,边缘逐渐平滑。这正是多尺度小波变换对高频信息进行平滑滤波的体现。但同时也存在图像细节丢失的缺陷,直观感觉就是图像变模糊。
3 边缘细化
经小波边缘提取得到的图像边缘比较粗,难以达到单像素的精度,需要对图像的边缘进行细化。图像边缘细化是图像处理中的基本技术,它要求完整地保存图像的拓扑结构,以便于代替原始图像进行识别和处理。将一个图像的主要边缘清晰、完整地提取出来,得到细化的、完整的边缘,将为图像检索、目标分割及识别等后续处理带来极大的便利[9]。
此处的边缘细化算法采用改进的Zhang-suen算法,其基本思想是在细化过程中不断移动3×3的模板,使它与图像中的各点重合。在扫描图像的过程中,定义待删除的像素为P1,P1周围的近邻像素位置关系如图2所示。
定义n(P1)和s(P1)两个参数,其中n(P1)是与P1相邻的非零像素个数,s(P1)是沿着P2→P3→P4→P5→P6→P7→P8→P9→P2的顺序由0过渡到1的总次数。对于考察点P1=1,如果同时满足2≤n(P1)≤6、s(P1)=1、P2P4P8=0且P2P6P8=0这4个条件,则可将P1删除。细化过程是将满足条件的像素点予以删除。反复迭代直到再也没有像素点满足上述条件为止,此时完成检验,得到了细化后的图像边缘。这样,所有检测到的边缘可以定位到单像素宽度,完整地保留了图像的边缘信息[10]。
4 实验结果比较分析
图3给出了不存在噪声时,传统边缘检测算子检测的结果。
由图3可以看出,传统的边缘检测算子在没有噪声的情况下可以得到较好的边缘检测效果。实际的数字图像总是存在一定功率的噪声,如椒盐噪声及高斯噪声等,其中椒盐噪声可以通过中值滤波较好地去除,高斯噪声是图像噪声容限方面要考虑的主要噪声。图4给出了存在高斯噪声(μ=0,σ=0.01)时传统的边缘检测算子的检测的结果。
可以看出,在有高斯噪声存在的情况下,传统的边缘检测方法检测的图像边缘出现缺失,伪边缘被检测出大范围出现,对边缘检测结果造成了极大破坏。采用基于二进小波变换的多尺度边缘细化检测方法检测得到的边缘如图5所示。
从视觉直观感受来看,加入高斯噪声的传统边缘检测方法的检测结果受噪声干扰很大;而采用笔者提出的方法后,图像边缘的连续性和准确性都比传统方法的检测结果有较大改善。
在此,定量地来探讨图像边缘检测结果的比较。边缘检测过程会产生3个主要类型的误差,即边缘点丢失、边缘点定位偏离以及将噪声波动误认为是边缘点。以Canny最优边缘检测准则为参考依据,Canny最佳边缘检测三准则具体如下:
a.最优检测。对于真实存在的边缘不漏检,当然也不会把非边缘点检出,使得输出信噪比最大。
b.最优检测精度。所得边缘点的位置与实际边缘点的位置最近。
c.检测点与边缘点一一对应。每个实际存在的边缘点与检测到的边缘点一一对应。
设f(x)为用于边缘检测的滤波器,边缘发生处的x=0,信号中的噪声是加性高斯白噪声n(x),其方差为n02,则边缘检测函数的3个性能指标———信噪比SNR、检测精度L和伪边界平均距离M的计算式分别为[11]:
信噪比SNR越大,能够准确检测到的真实边角和较少含有伪边界的可能性也就越大。检测精度L为检测到的边界与真实边界之间倒数的数学期望,检测精度越高,测量误差越小。伪边界平均距离M为随机噪声与检测函数卷积之后伪边界出现的平均距离,伪边界平均距离越长测量结果中出现伪边界的个数就越少[12]。
以主要性能指标信噪比SNR作为边缘检测评价标准,将没有噪声存在时的Canny边缘检测算子检测的边缘作为实际边缘,与存在(μ=0,σ=0.01)高斯噪声时传统的最优边缘检测算法Canny边缘检测和笔者所提方法进行比较,Canny边缘检测的SNR=1.005 0,小波多尺度边缘检测的SNR=1.534 8。
5 结束语
基于尺度不变特征变换的图像检索 篇4
基于内容的图像检索(CBIR)是计算机视觉在图像检索中的一个重要应用。常见的基于内容的检索利用全局内容特征,例如颜色、纹理以及形状中的一种或多种特征来进行图像的检索。这种检索往往忽略了图像内容的局部细节部分,而不能达到很好的检索性能。基于局部特征的检索是通过提取图片中的一些关键点,并通过这些点的相似度来决定两幅图片的相似度。局部特征已经被普遍应用于图像检索中[1,2,3]。因为它们不仅容易计算,具有一定的图像变换不变性,而且对于部分被遮挡的物体也能较好地识别。
1 SIFT描述子构造
SIFT是一种对尺度、旋转、仿射以及亮度变化都有很好的不变性的特征[4,5,6]。SIFT的提出是为了解决Harris的角点检测不具有尺度不变性的问题。 SIFT特征描述子的生成主要包括极值点检测、关键点定位、关键点梯度方向计算以及SIFT描述子4个步骤。
1.1 极值点检测
第一步是在尺度空间上检测极值点。SIFT中使用的是高斯核差分的尺度空间,它近似于归一化的拉普拉斯尺度空间。对于要处理的图片,SIFT先用不同的核对其进行高斯平滑。I(x,y)表示原图,G(x,y,σ)为高斯变换,L(x,y,σ)则为原图经过高斯变换后的图片。如下式:
L(x,y,σ)=G(x,y,σ)*I(x,y)。 (1)
相邻尺度的高斯平面相减后,得到不同尺度的高斯核差分(DoG)的尺度空间。在DoG尺度空间上进行的极值点检测能使SIFT特征具有较好的尺度不变性。DoG可直接通过相邻尺度的高斯面相减而得到。如图1中,左边的是经过不同高斯平滑后的高斯平面,右边是相邻尺度的高斯面相减后得到的高斯差分面(DoG)。计算出高斯差分平面之后,对于其上的某一点,将它与周围的8近邻点以及上下相邻尺度的对应位置的9近邻点(共26个点)进行比较,若该点都大于或都小于这26个点,那么这个点将被选取为极值点。
1.2 关键点准确定位
DoG空间中选取出来的极值点是灰度变化的极值点,包含着图片的结构信息。但这些点还需要进一步的处理,去掉低对比度的点,因为低对比度的点对噪音比较敏感。首先,对DoG尺度空间函数D(x,y,σ)进行Taylor展开;其次,将检测到的极值点
由于DoG对边缘有很强的响应,因此检测到的边缘点也要去掉。DoG在沿着边缘方向上的主曲率值会很大,而在与边缘方向垂直方向的主曲率值会较小。令α为沿着边缘方向的主曲率,β为与边缘垂直方向的主曲率。α=γβ,则当γ大于某个阈值(10)时,则作为边缘点舍弃。
过滤掉对比度低的点及边缘的点后,剩下的点即为关键点。由于关键点分布在不同尺度空间的不同位置上,这样可保证SIFT的尺度缩放不变性。
1.3 关键点梯度方向计算
在确定关键点后,为了保证特征的旋转不变性,SIFT通过关键点周围区域的梯度方向分布来表示关键点的方向。在尺度空间的高斯图片L( x,y,σ)中,关键点 ( x,y ) 周围区域的梯度模值m( x,y )和方向θ( x, y )如下:
然后对关键点周围区域的梯度方向进行直方图统计,以每10°为一个桶,平面上360°角可分为36个桶。根据梯度方向值将其关联到对应的桶中,权值由该梯度的模值大小来确定。SIFT用以关键点为圆心的高斯加权窗口来对梯度模值进行加权。
对梯度方向的直方图统计后,每个关键点的方向被定义为它的梯度方向的直方图统计中最大值的方向以及大小在最大值的80%以内的方向。
1.4 SIFT特征描述子构造
前面的步骤确定了关键点的尺度、位置以及方向,保证了其尺度缩放不变性及旋转不变性。关键点的检测已完成,现在将构造SIFT的特征描述子。SIFT特征是以它周围区域的像素点的梯度来表示。
梯度模值和方向的计算与第3步中计算关键点方向时的梯度计算的方式类似,只是这里选取的桶个数是8。对于以关键点为中心的周围区域16*16的像素点,计算其梯度方向和模值大小,然后以4*4的区域为单位进行方向的直方图统计。于是,原来的16*16区域的像素点被划分成4*4个4*4的小区域,每个小区域由对应的8个方向的梯度模值表示。因此,最终特征向量为:4*4*8 = 128维的。这便是SIFT的特征描述子。
为了保证SIFT特征对亮度变化的不变性,还需要对提取出来的特征向量进行归一化。这样,便得到了一个对于尺度缩放、旋转和亮度变换等具有不变性的特征——SIFT特征。
2 基于SIFT相关特征的图像检索
SIFT特征计算的简单性,对于尺度,旋转等变换具有很好的不变性,使它被广泛应用到图像检索中[7,8,9]。基于SIFT特征的图像检索的QBE过程,实际上是一种特征提取以及特征匹配的过程。它包含2个步骤:特征的提取跟特征相似度的计算。
2.1 特征提取
SIFT特征的提取包括2个阶段:离线的操作和在线的操作。对于庞大的图片集,先离线提取它们的SIFT特征存到数据库中,作为图像检索时的检索数据使用。在线图像检索时,对用户给出的检索图片提取SIFT特征后,与数据集中的特征进行匹配。然后,根据图片匹配的点数进行降序排序,匹配数目越多说明图片越相似。于是,图片的匹配转换成了图片中关键点特征的匹配。
将SIFT特征的提取分成离线和在线2个阶段是很有意义的,因为耗时的SIFT对于在线检索的实时性是很不利的。而如果先离线提取大量图像的SIFT特征,在线检索时只需要处理一幅图片的特征提取以及相似特征的检索及匹配问题。这样能很大地提高检索速度。
2.2 特征相似度计算
对于提取后的SIFT特征,检索过程要计算示例图片的特征向量与数据集中的特征向量的相似性。对于特征空间中的点,2点间相似度度量有2种常用的方法:一种是直接计算点间的距离,距离越小,相似度越大,常用的方法有曼哈顿距离和欧氏距离;另一种是计算给定点的最近距离和次近距离,当它们的比值越小,相似度越大[10]。第2种方法相对而言更准确,因为若2个特征向量匹配时,则它们间的相似度应该远大于它们和其他点间的相似度。因此,文中对于示例图片中特征向量的相似特征的检索采用的是第2种方法。
3 实验结果分析
实验中将基于SIFT的检索应用于3种类型的图片集,主要是通过对不同类型的图片比较SIFT特征进行匹配的效果。图片集包含600幅图片,其中建筑物图、室内图以及商标的图片各200幅。3种图片集的选取出于这样的考虑:室内的图片是为了观察对于室内的物体多且杂乱时,SIFT特征的性能如何。对同一建筑物的不同视角,以此对比同一物体的SIFT特征匹配的效果。对于企业徽标进行的实验主要是讨论SIFT对于相似的抽象图案的匹配结果。表1给出了利用SIFT特征对3类数据集的检索结果。
实验结果总结SIFT特征的特点如下:首先,SIFT特征是一种具有尺度、旋转不变性的局部特征,它对于视角变化(±35°)的图片也能较好的匹配。其次,SIFT特征提取的特征数目与图片的具体内容以及图片的分辨率成正比。最后,SIFT特征很适合用于相似图案的检索,它对同一场景或同一物体不同大小或小视角变化的检索结果也较好。从实验结果中可以看到,SIFT特征仍有它的局限性:一是SIFT的特征向量维数太高,对于较大的图片SIFT特征提取的特征数目很大,使得匹配速度慢;二是SIFT特征仅依靠点来描述物体,对位置等全局信息包含得还不够。
4 结束语
SIFT特征是一种对尺度,旋转变换等具有较好的不变性的特征。它很好地描述了局部信息,因此它在计算机视觉相关的邻域有很好的应用。同时,它也有特征向量维度高、匹配时间开销大以及缺少对全局信息的描述,只适用于对同一物体的匹配的局限性。高维特征向量的问题,可以通过对SIFT特征进行降维减少时间开销,或者通过采用分布式计算,从而提高检索速度。缺少对全局信息的描述,可以把SIFT与其他不同的特征(如颜色、形状和位置信息等)结合起来,提高全局信息的检索精度。
参考文献
[1]SCHMID C,MOHR R.Local Grayvalue Invariants forImage Retrieval[J].IEEE Transactions on PatternAnalysis and Machine Intelligence,1997,19(5):530-534.
[2]MIKOLAJCZYK K,SCHMID C.Indexing Based on ScaleInvariant Interest Points[C]//In Proc.of InternationalConference on Computer Vision,2001:525–531.
[3]KE Y,SUKTHANKAR R.Efficient Near duplicateDetection and Sub-image Retrieval[C]//InternationalMultimedia Conference Proc.of ACM,2004,60(2):869-876.
[4]MORTENSEN E N,DENG Hong-li,SHAPIRO L.A SIFTDescriptor with Global Context[C]//Computer Visionand Pattern Recognition,IEEE Computer SocietyConference,2005:184-190.
[5]EGENHOFER M.IEEE Symposium on Visual Languages[C]//IEEE Computer Society,1996:60-67.
[6]BRUNELLI R,MICH O.Image retrieval by examples[J].IEEE Trans.Multimedia,2000,2(3):164-171.
[7]EAKINS J P M,GRAHAM M E.Content-based ImageRetrieval:A Report to the JISC Technology ApplicationsProgramme[R].In:Institute for Image Data Research.University of Northumbria at Newcastle,Newcastle uponTyne,1999.
[8]汪道寅,胡访宇.基于改进SIFT算法的视频序列图像配准[J].无线电工程,2011,41(2):16-18.
[9]李向阳,庄越挺,潘云鹤.基于内容图像检索技术与系统[J].计算机研究与发展,2001,38(3):344-354.
多尺度变换 篇5
滚珠丝杠副是数控机床的关键功能部件,其质量的优劣直接影响数控机床整体性能[1],而振动是影响滚珠丝杠加工表面粗糙度的重要因素。在滚珠丝杠精加工中,磨削是一种常用的加工方法。磨削中的振动对丝杠加工表面粗糙度有着重要的影响,为了实现在线检测及故障诊断,就必须对振动信号进行实时分析。然而在磨削过程中,往往由于振动信号测试设备或是人为因素的干扰导致振动信号呈现低信噪比及野值的特征,从而使得振动信号的可观测性大大地降低。长期以来,国内外学者在振动问题的研究上,大多数集中在机床振动和切削振动方面,滚珠丝杠磨削振动研究相对较少[2]。特别是在多传感器数据进行融合测量与处理[3,4,5]时,各传感器报送测量数据中不可避免地会出现数据突变点,亦即野值。融合处理对采样数据中包含野值点反应较为敏感,如果不将这些野值预先剔除,就将给数据处理带来很大的误差。
小波多尺度分析是时间频率分析的一种新技术,在信号处理和图像分析领域得到了广泛的应用,多尺度变换是小波变换中的重要内容。小波分析或多分辨分析是继Fourier分析之后纯粹数学和应用数学完美结合的典范。小波分析是一种信号的时间尺度(时间-频率)的分析方法,具有多分辨率分析的特点,在时频两域都具有表征信号局部特征的能力,是一种窗口大小固定但形状可变,时间窗和频率窗都可以改变的时频局部变化分析方法。即在信号低频部分具有较高的频率分辨率和较低的时间分辨率,在高频部分具有较高的时间分辨率和较低的频率分辨率。利用小波的这些优点,可以对获取的微弱的信号进行预处理,在强背景噪声下显现并增强信号的目标特性,提高信噪比,同时小波变换可以用于各种奇异信号的检测[6,7,8]。
本文建立在仿真试验及实验数据分析的基础之上,提出了基于q带正交小波多尺度变换的振动信号处理方法,研究了滚珠丝杠磨削加工中低信噪比、野值的振动信号特征。通过小波多尺度变换技术,提高了滚珠丝杠振动信号的信噪比,并且得到了去野值的振动信号的本质特征,为磨削滚珠丝杠实时监测和故障诊断提供了有效的理论基础。
1 q带正交小波多尺度分析
正交小波变换能够将信号投影到一组正交的基函数而并不改变信号本身。如果信号的关键特征能够用少量的基函数表达,那么,通过将信号投影到这些正交的基函数,利用小波较好的空间域和频率域定位特性,信号的主要特征将很容易检测出来。q带正交小波变换对于信号的分解是将同一频率带等分为M段,相对于二进制小波对高频部分的划分就更加细,从而对高频处理上优于二进制小波,而且实验证明它对普通信号处理的效果也比较理想[9,10]。
定义1:若满足
则称φ(t)是q尺度分析的尺度函数,并对应于(q-1)个小波φ(i)(t),满足
令
并记尺度空间子Vj和小波空间子wj分别为:
则对于观测信号Z(k)在一定的尺度/分辨率下,其分解和重构如式(4)。
式(5)中:AjZ(k)为信号Z(k)的离散逼近,)为信号的细节信息。因此用多尺度分析工具可以将信号Z(k)分解为近似信息与细节信息之和,再通过重构可以对原始信号进行有效的滤波处理,提高信号的信噪比。
2 q带正交小波多尺度的处理性能分析
为了研究q带正交小波多尺度的处理性能,本文通过设计一个模拟信号来验证q带正交小波多尺度的处理性能。设模拟振动信号为:
式(6)中r是均值为零、正态分布的噪声,单位:加速度(g),Q为模拟信号中野值,模拟振动信号的采样频率为1 000 Hz。如图1所示,为噪声与模拟信号的时域波形,其中在0.4 s及0.8 s附近模拟信号中出现了野值。
根据上述q带正交小波多尺度理论,令q=2,采用’db2’小波,对模拟信号式(8)进行了尺度因子为5的分解,分解结果如图2所示,图中‘cd‘表示为模拟信号被分解出来的细节信号。图3为模拟信号的时频分布展开图。
由图2及图3可知,根据模拟信号的时域特征可知,模拟信号包含三个正弦信号及两处的野值,而从模拟信号多尺度分解后的细节信号图中可知,细节信号cd3,cd4,cd5为周期正弦信号,而cd1,cd2含有野值信号,因此可以在重构模拟信号时,有选择地选择细节信号进行原始模拟信号的近似。如图4所示,为重构前后的模拟信号对比图,由图可知,重构后的信号不仅提高了信噪比,而且还有效的抑制了野值现象。
3 滚珠丝杠磨削振动实验
3.1 滚珠丝杠磨削实验设计
根据文献[2]可知,滚珠丝杆在磨削加工时,头架处的振动频率最大,本文为了处理在实际磨削环境下测试所得的最大振动信号,因此对滚珠丝杠在加工时头架处的振动信号进行了实验测试。实验中,安装于头架上的传感器如图5所示。所用加速度传感器型号为TG—320。采用VXI数据采集仪来采集振动信号,实验中振动信号采样频率为1 000 Hz。分别在机床x,y,z三个方向进行了测试,其中x为垂直于滚珠丝杠轴向,并指向砂轮面为正向,y为滚珠丝杠轴向方向,并指向头架向为正向,z向与x,y组成右手法则。
3.2 振动测试信号分析
通过上述实验设备,分别在x,y,z三个方向上对滚珠丝杠在磨削加工时的头架处振动信号进行了测试,测试结果如图6—图8所示。为了更好地研究振动信号内部的特征,根上述正交小波多尺度分解理论,对三个方向的振动信号进行了多尺度分解,并且分别在时域和频率上进行了对比研究。
由图6—图8可知,头架在x方向上的振动能量最强,在z方向上的振动最小,x,y,z三个方向的振动信号经过多尺度分解后呈现出多频率的特征,说明滚珠丝杠在磨削加工中的振动信号为非平稳信号,且在高频阶段从信号的时域可知,呈现出高斯噪声的特性,说明振动信号中存在高频噪声。由幅值谱能够看到多尺度分解后细节信号的频率成分按照由大到小的顺序被逐渐分解,同时也较好地保持了原始振动信号的特征。在振动信号分析中,由于滚珠丝杠磨削加工过程中产生的不同信号成分,可以根据频率的成分来选择有效的振动信号,也可直接根据分解的细节信号时域特性来判断振动的严重程度,进行故障诊断及决策。
4 结论
本文提出了基于小波多尺度分解的振动信号滤波方法,通过仿真模拟信号及实验测试滚珠丝杠在磨削加工过程中的振动信号,验证了正交小波多尺度分析对具有低信噪比、高野值及非平稳特性的振动信号分解的有效性,实现了对滚珠丝杠振动信号平稳性分析,揭示了滚珠丝杠磨削加工振动信号的本质特征,为滚珠丝杠在线检测及故障诊断提供了理论基础。
参考文献
[1]王兆旦,朱继生,张瑞,等.滚珠丝杠副性能指标的测试与研究.制造技术与机床,2009;6:128—134
[2]于艳美,李春梅,冯虎田.滚珠丝杠磨削加工中的振动测试与分析.机床与液压,2011;3(39):17—21
[3] Waltz E,Lilnas J.Multisensor data fusion.Bos-ton:ArtechHouse,1990
[4]何友,王国宏,彭应宁.多传感器信息融合及应用.北京:电子工业出版社,2000:38—45
[5]林云.多尺度小波变换在野值剔除中的应用.航天电子对抗,2009;6(29):54—58
[6]许郴,郑链,王克勇,等.基于信号奇异性分析的小目标检测方法.红外技术,2005;27(3):245—249
[7]李利荣,张桂林.基于小波变换的小目标检测.测控技术,2004;23(6):67—69
[8]程正兴.小波分析算法与应用.西安:西安交通大学出版社,1998
[9] Daubechies I.Ten lectures on wavelet.philadephia,Pennsylvania:So-ciety for Industrial and Appl math,1992
分层多尺度建模-计算方法 篇6
1 基本假设
我们所考虑的材料结构都被假定为宏观上足够均匀, 但在微观上是不均匀 (异质) 的 (以可区分的组成物来说, 例如夹杂、晶粒、界面及空洞) , 如图1所示。
由于此细观尺寸比分子尺寸要大得多, 所以我们可以将每种组成相当做连续介质来看待。与此同时, 以尺度区分的准则来说, 细观尺寸又要比宏观试样的特征长度或宏观下载荷的波长要小得多。
目前大部分的多尺度方法都假定材料具有周期性的细观结构, 即认为整个宏观结构件都由空间单胞重叠而成。目前, 有些学者采用了更加合理的假设:局部周期假设。即在宏观上不同的点相应的细观结构可以不同, 而在宏观点周围一小块部分细观形态是重复的。局部和整体周期性的概念示意图如图2所示。笔者采用的是第一种假设。即整体周期性假设。
在此方法中, 每一个宏观材料点的变形 (梯度) 张量FM需要先计算出来。这个材料点的变形张量FM参与构建RVE的边界条件。然后求解RVE的边界值问题, 宏观的应力张量PM就可以通过相应RVE的应力场的体积平均来得到。因此, 在宏观材料点上的变形 (应变) -应力数值关系就很容易得到了。方法流程如图3所示。
2 细观尺度上的问题说明
材料细观结构的物理和几何特征由代表性体积单元RVE来确定。一个典型的二维RVE如图4所示。实际上对RVE的选取是一项十分复杂的工作:
RVE尺寸必须足够大来代表材料的细观结构, 同时RVE尺寸又要足够小来进行更加效率的建模和计算。在文献[2,3,4,5]中详细说明了代表性体积单元RVE的概念和构建。这里假定一个合适的RVE已经建好。
给宏观材料点指定的RVE, 已知RVE的初始状态向量为X (在参考体积V0范围内) , 当前位置向量为x (在当前体积V) , 细观结构变形梯度张量表示为:Fm= (∇0mx) c, 其中∇0m是关于所参考微观结构的梯度算子;c表示共轭。
如图4所示, 此RVE处于一个状态, 在数学上反映为关于柯西应力张量σm或者关于第一Piola-Kirchhoff应力张量Pm=det (Fm) σm (Fmc) -1的平衡方程式, 其表达式为 (不考虑体力) :
∇m·σm=0 V中 或
∇0m·Pmc=0 在V0中 (1)
其中∇m是关于当前单胞细观结构形状的梯度算子。
各细观组成物的力学特征由各自的本构关系来描述, 现对各细观组成物指定时间相关的应力-应变关系:
Pm (α) (t) =RP (α) {Fm (α) (τ) , τ∈[0, t]} (2)
其中t表示当前时间, undefined为可区分的细观组成物的数量 (例如基体、夹杂等) 。
然后指定位移边界条件, 变形状态下RVE上一个点的位置向量可表示为:
x=FM·X (X在Γ0上) (3)
其中Γ0为RVE上未产生变形的边界。
3 宏-细观尺度耦合
使用的是位移边界条件, 假定宏观变形张量FM为细观变形张量Fm的体积平均:
undefined
现在验证位移边界条件 (3) 是否满足式 (4) , 把式 (3) 代入式 (4) , 然后使用散度定理∇0mX=I。
undefined
=FM (5)
4 应用实例
镍基超合金广泛应用于航空、发电站等领域, 尤其是发动机涡轮叶片等热端部件。
4.1 缺口试样拉伸实验的有限元建模
缺口试样拉伸实验主要用于对金属材料塑性损伤和断裂的研究。圆形缺口试样示意图如图5所示, 单位为mm。
使用分层多尺度方法进行此材料拉伸实验的有限元模拟, 考虑到试样是绕中心线轴对称的, 所以使用轴对称单元进行二维分析。试样的二维几何模型及网格划分如图6所示。
由于试样采用镍基超合金, 此材料为多晶体材料, 其细观结构的几何模型和网格划分如图7所示。二维多晶体的几何模型构建参见文献[5,6]。
单元类型上, 选择一阶常规单元, 不使用减缩积分。最后的模型中, 二维网格单元总数为343。关于单元的详细信息见表1。
4.2 材料定义和边界条件设置
缺口试样采用某镍基合金的材料数据, 宏观材料参数见表2。
微观结构采用晶体粘塑性本构模型。下面简要列出其弹性本构方程,
undefined
其中左端是以中间构形为基准状态的Kirchhoff应力张量τ的Jaumann导数。L为刚度张量。
镍基合金各组成相得流动法则和硬化规律方程参见文献[7]的晶体塑性理论部分。
4.3 计算结果
宏观试样分析和细观结构分析的有限元分析结果如图8所示, 图8 (a) 为宏观试样在拉伸载荷下的最大主应变分布, 图8 (b) 是在细观尺度下RVE的最大主应变分布。由图可以看出, 试样在缺口底部的应变量最大, 其值为0.3789。细观结构进行位移继承后得到的最大主应变结果为0.3871, 两者差异2.16%。这在计算缺口部位应变能和疲劳强度的时候相差是巨大的。
引起细观结构应变结果比宏观分析要大的原因可能是由于细观尺度下晶粒大小、晶体学取向等的影响[8]。因此, 为了要获得这种影响的机理进行详细的多尺度分析是必要的。
5 结论
首先提出了多尺度建模方法, 然后讨论了此方法的可行性, 最后将此方法应用到圆形缺口试样拉伸实验的有限元模拟。完成了缺口试样从宏观到细观的跨尺度计算。得到了如下启示:
(a) 宏观 (b) 细观
1) 从宏观试样的分析结果中我们知道在缺口底部会出现应力应变集中, 一般把此区域称为危险区域, 选用此区域的单元作为宏-细观跨尺度计算的连接点;
2) 传统的宏观计算由于没有考虑细观因素的影响, 可能会导致错误的结果, 所以进行考虑细观结构影响的宏-细观的计算是必要的。此方法已经成功应用于某镍基合金的跨尺度分析。
参考文献
[1]V.Kouznetsova, Computational homogenization for themulti-scale analysis of multi-phase materials[D], Eindhoven University of Technology, 2002.
[2]J.Zemanand, M.Sejnoha, Numerical evaluation of effec-tive elastic properties of graphite fiber towimpregnated bypolymer matrix[J], J.Mech.Phys.Solids, 2001, 49:69-90.
[3]Z.ShanandA.M.Gokhale, Representative volume elementfor non-uniform micro-structure[J], Comput.Mater.Sci, 2002, 24:361-379.
[4]C.Huet, Application of variational concepts to size effectsin elastic heterogeneous bodies[J], J.Mech.Phys.Sol-ids, 1990, 38 (6) :813-841.
[5]黎倪.材料微观组织结构的有限元模拟[D].兰州理工大学材料科学与工程学院, 2005.
[6]刘玉振.基于Abaqus的材料微结构有限元计算的后处理研究[D].兰州理工大学材料科学与工程学院, 2008.
[7]王自强, 段祝平.塑性细观力学[M].北京:科学出版社, 1995.
多尺度变换 篇7
标准Kalman滤波是用于线性系统的最小均方意义上的最优状态估计。Kalman滤波对非平稳信号具有较强的估计能力,但标准最优Kalman滤波的缺点和局限性是要求知道系统的数学模型和噪声统计特性。所以,在单一尺度上,如果由于多种不确定因素干扰,导致使用不精确的模型和噪声统计,设计出的Kalman滤波器会导致滤波器性能变坏[1]。
多分辨率的数据进行多级别、多层次的处理能够获得更有价值的信息,而这种信息是单一尺度上所无法获得的。通过多尺度分解,在不同的尺度空间上分别刻画目标的特征属性再加以融合,可以有效地降低不确定因素的干扰,从而提高滤波器的滤波性能[2,3,4]。小波变换的多尺度特点非常适合多尺度信号的处理,特别是依靠小波变换的重构算法对低分辨率补充一定信息可以获得高分辨率近似补充。文献[5]中提出基于小波变换,寻找基于全局信息的一步预测值和估计误差协方差矩阵,通过综合不同尺度的信息来降低Kalman滤波对数学模型和噪声统计特性的依赖,这种方法比起单一尺度上的Kalman滤波效果是明显的,但计算量较大,由于算法侧重于横向数据更新,全局最优预测值和误差矩阵对滤波效果的影响变得较大,容易受到干扰。本文就此提出一种基于多尺度Kalman的融合滤波方法。充分利用初始估计序列来有效地降低向量和矩阵维数,减少运算量,把标准Kalman滤波只在单一尺度和时间轴上,对状态估计值和误差协方差进行数据更新,改进为基于小波变换在尺度轴和时间轴上,通过不同尺度的观测数据,纵向和横向双向数据更新,算法侧重于纵向数据更新,充分利用多尺度数据进行多层次融合处理,从而获得更好的滤波效果。
本文主要思想可以概括为:首先,在最细尺度N上通过标准Kalman滤波得到初始估计序列,然后通过小波变换将状态估计值和状态估计误差协方差的分块序列分解到最粗尺度,并根据相应尺度上的观测数据进行数据更新,然后通过小波重构,实现不同尺度上的数据更新,当纵向回归到原始尺度即最细尺度时,再跳至下一个数据块进行,直至滤波结束。在实际中,也可以利用此方法解决多传感器多分辨率数据的融合,更好地利用不同分辨率数据的互补信息,达到更佳的融合效果。
2 多尺度动态系统的描述及Kalman最优滤波
通常情况,动态系统状态方程在尺度N上的描述为
其中:k为离散时间变量,x(N,k)是状态向量,Φ(N,k)是系统状态转移矩阵,w(N,k)是系统噪声,且服从均值为零,方差为Q(N,k)的正态分布,初始值x(N0,)为一随机变量,其均值和方差是E[x(N,0)]=x0和E{[x(N,O)-x0][x(N,O)-x0]T}=P0。
动态系统量测方程在尺度N上的描述为
其中:C(N)系统观测阵,z(N,k)是对x(N,k)的观测值,观测噪声v(N,k)服从均值为零,方差为R(N,k)的正态分布,初始值x(N0,)与观测噪声v(N,k)以及系统噪声w(N,k)之间是互不相关的。
在尺度N上标准Kalman最优滤波的基本方程,如式(3)~式(7)所示:
3 基于小波变换的多尺度动态系统描述
我们由小波理论可知[6],通过一个脉冲响应为h(l)的低通滤波器可以从尺度i上获得粗尺度i-1上的平滑信号xL(i-,1k),脉冲响应为g(l)的高通滤波器可获得细节信号xH(i-,1k):
通过小波变换将状态方程从尺度i分解到粗尺度i-1,我们可以发现不同尺度之间的参数矩阵之间的内在联系,为了计算方便,一般我们取状态方程中的G(N,k)为单位阵,则i-1尺度上系统状态方程为
其中上标i表示结果由尺度i上小波分解得到的,比较式(10)和(11)可得式(12)~式(13):
由式(13)可推导出:
同样,用小波变换将量测方程从尺度i分解到尺度i-1,则i-1尺度上系统量测方程为
比较式(15)和式(16)可得式(17)~式(18):
从式(12)∼式(13)和式(17)∼式(18)可以看出,通过尺度i上的系统状态和量测方程的参数矩阵,可以很容易地获得尺度i-1上的相应的参数矩阵,研究标准Kalman最优滤波的基本方程可知,在尺度i-1上进行标准Kalman滤波估计还需要和Pk+1/k(i-1),我们将在多尺度Kalman滤波算法中给出计算方法。
4 基于小波变换的多尺度Kalman滤波
首先,假设一长度为T的信号,我们在单一尺度N上用标准Kalman滤波根据观测值进行估计。由于多种不确定因素干扰,滤波效果并不会很理想,在尺度N上,经Kalman滤波后,我们将得到估计序列和相应的估计误差协方差阵序列Pk/k(N),序列长度均为T,我们把得到的这一系列估计序列定义为初始估计序列,将这些序列分割成长度为M的数据块,其中T和M均设定为2的整数次幂,方便后续的小波变换。这样分割后的序列可表示为Xm(N)和Pm(N),则:
对于数据块Xm(N)中的每一个数据元素均来自于初始估计序列的估计值,由于我们要在尺度轴上对此估计值进一步修正,所以我们可以将此值看作一步预测状态估计值,为了保持与标准Kalman最优滤波方程中的符号的一致性,我们用符号来代替数据块Xm(N)的表示,此时是由初始状态估计值组成的数据块。
由式(8)~式(9)可知:经过从尺度N到尺度i的小波塔式分解,信号可分解为尺度i上的平滑信号和相应的各尺度l(i≤l≤N-)1的细节信号,这样我们把初始状态估计序列通过小波变换向尺度i分解,可得到尺度i上的平滑信号和相应各尺度l上的细节信号,如式(21)~式(22)所示:
对于初始状态估计序列,分解后不同尺度上的序列集合为
接下来,我们要获得相应尺度上的预测误差协方差阵。由估计误差协方差阵的定义可知:
其中:为尺度N上的一步预测误差,由于我们是对长度为M的数据块操作,所以此时估计误差协方差为一矩阵块。根据初始估计误差协方差阵序列,我们可以获得数据块Xm(N)的估计误差协方差阵为,如果用表示第i行第j列的元素,则:
同样,为了保持与标准Kalman最优滤波方程中的符号的一致性,我们用符号Pk/k-1(N)来代替数据块表示。这样对Pk/k-1(N)做小波变换实际上是一个二维离散小波变换[6],通过小波变换我们可以得到尺度i上的Pk/k-1(i)。
若在单一尺度上测量,我们可以通过下采样Zk(N)获得尺度i上的测量序列Zk(i),若是多传感器多分辨率系统,则尺度i上的测量值Zk(i)可以通过不同传感器的观测获得。在尺度i上,我们利用Kalman滤波方程对数据块进行数据更新。由式(17)、(18)、(23)和(26),我们可以得到尺度i上的Ck(i),Rk(i),和Pk/k-1(i),并通过式(27)~式(31)在尺度i上Kalman滤波:
数据更新结束后,通过小波重构我们可以得到尺度i+1上的以及Pk/k-1(i+1),在通过尺度i+1上观测值Zk(i+1)用同样的方法对序列进行数据更新,以此类推,直至尺度N,我们便得到了尺度N上的。同时在小波重构时,还可根据实际的研究情况,结合一定的小波域值去噪方法来去除部分噪声影响[7,8]。
5 仿真算例
实验一
系统在尺度4上的状态模型和尺度j(j=4,3,2,1)上的量测模型,如式(32)~式(33)所示:
其中:、w(,4k)~N[,0Q(,4k)]、v(j,k)~N[,0R(j,k)];式中Q(,4k)=diag1(1,),测量矩阵C(j,k)是单位矩阵,测量误差协方差矩阵R(j,k)为
尺度4为最细尺度,初始值x0=[10,]0T和P0=diag(4,4),选用d B4小波分解和重构。图1(a)信号在尺度4上的量测值(T=512,含两个状态分量)。图1(b)为最细尺度(i=4)上Kalman滤波结果。图1(c)为融合4层信息后的多尺度Kalman滤波结果(M=32),可以很直观地看出,多尺度Kalman滤波的有效性。
从表1中可以看出,融合的尺度越多,测量值的改进越明显,但当融合尺度数目过多时,测量值的改进程度明显减少且运算量随之增大,这是由于测量值提供的有用信息变少了,在粗尺度上的数据更新也就相应变少。从表2中可以看出,不同数据块大小对滤波效果也有明显影响,当数据块选取过大时状态估计误差协方差矩阵较大,数据之间的影响变大,导致滤波性能降低。所以,我们可以确定实际需求,确定适当融合尺度数目和数据块大小有效地提高滤波性能和效率。
实验二
继续采用式(32)∼式(33)的系统状态模型和量测模型,尺度4为最细尺度。假定此系统为多传感器多分辨率系统,我们用三个分辨率不同的传感器来观测系统信号,并进行数值仿真,三个传感器对信号观测的结果如图2(a)∼图2(c)所示,其信噪比(SNR)分别为13.4d B、9.2d B、6.7d B,图中为信号的第一个状态分量,初始值x0=[10,0]T和P0=diag(4,4)。这样式(27)的尺度i上的测量值Zk(i)可以通过不同传感器的观测获得,我们选用d B4小波分解和重构,传感器1、2、3的量测分别对应尺度2、3、4上的观测,将三个尺度上的数据融合后的结果图如图2(d)所示,其信噪比为19.3d B,可以看出基于多尺度Kalman的多传感器数据融合滤波算法具有良好的滤波效果,可用于多分辨率多传感器数据融合。
6 结论
本文通过深入分析基于小波变换的动态系统模型,提出一种基于多尺度Kalman的数据融合滤波的方法,利用小波的多尺度特点,把初始估计序列多尺度分解,并进行分层Kalman滤波估计,通过小波重构进行估计融合。该算法将小波多尺度分解和Kalman滤波结合起来,同时还有效地利用了多分辨率的数据信息。通过实验可以看出,该算法对实际中含较强噪声的动态系统状态估计效果较好,同时也能用于多分辨率多传感器数据融合。
摘要：本文通过分析基于小波变换的动态系统模型,提出一种基于小波多尺度的Kalman数据滤波方法,本文利用小波的多尺度特点,把初始估计序列多尺度分解,并在不同尺度层上进行Kalman滤波估计,再利用小波重构来融合各层的估计信息,把标准Kalman滤波只在单一尺度和时间轴上对状态估计值和误差协方差进行数据更新,改进为基于小波变换的尺度轴和时间轴上的双向数据更新,该算法将小波多尺度分解去噪和Kalman滤波相结合,对实际中含较强噪声的动态系统的状态估计效果较好。算法也可用于多分辨率多传感器数据融合。
关键词：多尺度,Kalman滤波,小波变换,数据融合
参考文献
[1]邓自立,许燕.基于Kalman滤波的白噪声估值理论[J].自动化学报,2003,29(1):21-23.DENG Zi-li,XU Yan.White Noise Estimation Theory Based on Kalman Filtering[J].Acta Automatica Sinica,2003,29(1):21-23.
[2]HONG L,CHEN G R,CHUI C K.A Filter-bank-based Kalman Filtering Technique for Wavelet Estimation and Decomposition of Random Signal[J].IEEE Trans.On Circuits and systems II:Analog and Digital Signal Processing,1998,45(2):237-241.
[3]WEN X B,TIAN Z,LIN W.Robust Kalman Filter and Smoothing Recursive Estimator for Multiscale Autoregressive Process[C]//ICSP’04.7th International Conf.Beijing:IEEE,2004,1(1):364-367.
[4]雷明,韩崇昭,元向辉,等.非线性系统的小波分频的扩展Kalman滤波[J].光电工程,2006,33(6):68-72.LEI Ming,HAN Chong-zhao,YUAN Xiang-hui,et al.Wavelet Decomposing based Subbandwise Extended Kalman Filtering for the Nonlinear System[J].Opto-Electronic Engineering,2006,33(6):68-72.
[5]文成林.多传感器单模型动态系统多尺度数据融合[J].电子学报,2001,29(3):341-345.WEN Cheng-lin.Multiscale Data Fusion for Multi-Sensor Single Model Dynamic Systems[J].Acta Electronica Sinica,2001,29(3):341-345.
[6]MALLAT S.A Wavelet Tour of Signal Processing(Second Edition)[M].San Diego:Academic Press,1999.
[7]彭玉华.一种改进的小波变换阈值去噪方法[J].通信学报,2004,25(8):119-123.PENG Yu-hua.An Improved Thresholding Method in Wavelet Transform Domain for Denosing[J].Journal on Communications,2004,25(8):119-123.