多分辨分析方法（精选7篇）
多分辨分析方法 篇1
图像融合是指将多源信道所采集到的关于同一目标的图像数据经过图像处理和计算机技术等,最大限度地提取各自信道中的有利信息,最后综合成高质量的图像,以提高图像信息的利用率、改善计算机解译精度和可靠性、提升原始图像的空间分辨率和光谱分辨率,利于监测[1]。图像融合技术可以分为3个层次:像素级融合、特征级融合和决策级融合。像素级融合是在信息的最底层进行的,但是由于它保留了最原始的图像信息,所以信息量最大,进度最高[2]。
目前像素级图像融合领域的研究方法常用的有线性加权平均法、小波变换法、IHS变换融合法、拉普拉斯金字塔变换法、PCA变换融合法等[3]。笔者分析了两种图像融合方法,基于塔形分解的多分辨率分析,对多聚焦图像进行融合,并比较两种图像的融合效果。
1 拉普拉斯金字塔变换
1.1 拉普拉斯金字塔的分解与重构
拉普拉斯金字塔的分解及重构包括3个步骤。一是图像的高斯金字塔分解。图像的拉普拉斯金字塔变换是建立在高斯金字塔变换基础上的。高斯金字塔是一个在尺寸上逐层减半的一组图像序列。序列中的每一级图像均为其前一级图像低通滤波后作隔行隔列降采样,即
式中,G0为源图像作为高斯金字塔的底层,Gl为第l层高斯金字塔图像;Cl和R1分别为第l层子图像的列数和行数;金字塔的总层数为N+1;ω(m,n)是一个二维分离的5×5窗口函数。
二是图像的拉普拉斯金字塔分解。将Cl内插放大,得到放大的图像Cl*,使Cl*的尺寸与Cl-1的尺寸相同,内插图像窗口表示为
式中
为整数时Cl*的尺寸与Cl-1的尺寸是相同的,但实际上二者并不相同。从上式可以看出,通过对原有像素灰度值的加权平均从而确定在原有像素间内插的新像素的灰度值。由于Cl是对Cl-1进行低通滤波得到的,即Cl是模糊化,降采样的Cl-1,所以Cl*所包含的细节信息少于Cl-1。
令
由LP0,LP1,LPl,LPN构成的金字塔即为拉普拉斯金字塔,这个计算过程相当于带通滤波,它的每一层图像是高斯金字塔本层图像与其高一层图像经放大算子放大后图像的差。因此,拉普拉斯金字塔亦可称为带通金字塔[4]。
三是由拉普拉斯金字塔重建原图像。
由式(3)可得
式(4)说明,拉普拉斯金字塔的各层图像经逐步内插放大到与下一层图像一样大,然后再相加就可精确重建原图像,这说明图像的拉普拉斯金字塔分解是原图像的完整表示[5]。
1.2 图像融合实验
图1选用两幅多聚焦图像进行融合,两幅图像为同一场景,大小为640*480,窗口大小都选为3×3,窗口选得太大会偏离中心像素太远,从而使融合效果不敏感,且计算量也会增加。实验平台MATLAB 7.0,结果见图2。
图2是先对两幅图像做三层的拉普拉斯金字塔变换,然后对应层对应点灰度取平均值,再做金字塔反变换所得融合图像。从融合图像可以看出,方法很好地保留了边缘和细节信息,几乎没有产生振铃效应,融合效果良好。
2 小波变换
2.1 小波变换的分解与重构
对于经过配准的两幅图像利用小波变换的基本过程见图3。
分解时使用的小波为bion4.2.小波分解的层数越多,融合时选择的频率范围就越多,但分解层数越多,得到的局部空间越小,系数数量按指数增加,处理的计算量越大,分解的层数不宜过多,一般取2～5层即可。利用先验知识确定所需处理的信息不在某一频段,相应的频段可以不分解[6]。
2.2 图像融合步骤
基于小波变换的图像融合具体步骤如下。一是对源图像进行预处理。图像滤波:对失真的图像直接进行融合,然后便会将图像中的噪声融入到融合图像中,这样就会影响融合图像的效果,因此需要对源图像进行滤波以消除噪声。图像配准:为了综合使用多种成像模式和多聚焦以提供全面的信息,所以笔者需要对有效的信息进行整合,使多幅源图像在空间域中的几何位置完全对应。二是对已配准的源图像进行小波分解,相当于使用一组高低通滤波器进行滤波,分别得到图像的水平高频分量、低频分量、对角分量和垂直高频分量。三是对各分解层分别进行融合处理,采取不同的融合策略,根据低频和高频分量的特点,在各自的变换域进行融合。四是对融合的小波金字塔进行二维离散小波逆变换,重建图像即可得到融合图像[7]。
2.3 图像融合实验
图4选用两幅多聚焦图像进行融合,两幅图像为同一场景,大小为640*480,窗口大小都选为3×3,窗口选得太大会偏离中心像素太远,从而使融合效果不敏感,且计算量也会增加。实验平台MATLAB 7.0,结果见图5。
在基于小波变换的图像融合中存在一个最优分解层数,且最优分解层数与融合策略是有关的。对融合图像质量的评价应从定量和定性两个方面来考虑。笔者得到的评价指标值大多数都随着分解层数的增加而增加,但从视觉上看融合图像会在某一分解层时达到一个最佳结果。
3 图像融合的评价
3.1 主观对比
对比图2和图4可以直接看出基于小波变换的融合结果要比基于金字塔拉普拉斯融合结果的效果好很多,相比之下基于小波变换的融合结果图像更清晰,目标更明确,更容易识别。
3.2 客观对比
3.2.1 多向粗糙度
1)计算图像中大小为2k×2k的活动窗口中像素的平均值,k=0,1,2,…,5。
2)对于每个像素分别计算它在垂直、水平、对角线方向上互不重叠的窗口之间的平均值差。在三个平均值差中选择能使其值最大的k值,相应的粗糙尺寸为Sc(x,y)=2k
3)整幅图像的多向粗糙度为
多向粗糙度的值越大,说明图像的细节信息越多。
3.2.2 标准差
标准差反映了图像灰度相对于灰度平均值地离散情况。
1)用G(x,y)=g表示图像中第(x,y)个像素的灰度为g,L为图像的灰度级数,图像的尺寸表示为M×N则灰度图像的均值表示为
2)设一幅图像的灰度分布为p={p(0),p(1),p(L-1)等于g的像素数与图像总的像素数的比值,且则标准偏差的定义为:
从表1和表2指标结果可以看出,基于小波变换的图像融合方法的指标最大,这也和主观评价的结果是一致的。
4 结论
小波变换融合法与金字塔图像融合法相比,具有很多优点。首先,小波变换数据的大小与图像相同,而金字塔的数据大小是图像的4/3,小波变换更为紧凑。其次,由于可以选择正交小波核,因此不同分辨率包含的信息是唯一的,而金字塔分解在两个不同的尺度之间含有冗余。再次,小波表达式提供了方向信息,而金字塔表达式未将空间方向选择性引入分解过程。最后,金字塔的重构过程可能具有不稳定性,特别是两幅图像存在明显差异区域时,而小波变换图像融合没有类似的问题。小波变换具有不同分辨率包含的信息唯一、重构过程稳定等在信号处理领域的优良特性,因此基于小波分析的融合方法正日益受到人们的重视。
参考文献
[1]郭雷,李晖晖,鲍永生.图像融合[M].北京:电子工业出版社,2008.
[2]闫敬文,屈小波.图像融合研究最新进展[J].厦门理工学院学报,2007,4(15):44-49.
[3]玉振明,高飞.基于金字塔变换的图像融合方法[J].计算机应用研究,2004,21(10):128-130.
[4]李建林,俞建成,孙胜利.基于梯度金字塔图像融合的研究[J].科学技术与工程,2007,7(22):5818-5822.
[5]黄文涛,毕笃彦,毛柏鑫,等.基于中值变换和金字塔分解的图像去噪方法[J].电子与信息学报,2004,26(11):1686-1692.
[6]刘坤,郭雷,李晖晖.一种基于静态小波变换的图像融合算法[J].计算机工程与应用,2007,43(12):59-61.
[7]玉振明,毛士艺,袁运能,等.基于边缘检测的小波变换的图像融合研究[J].电子学报,2005,13(8):35-39.
多分辨分析方法 篇2
表面形貌是指零件在加工过程中残留下来的微观几何形态,可用粗糙度、波纹度、形状误差等来表征。表面形貌客观地反映了表面生成机制,影响工程表面的功能特性[1]。为了有效地分离粗糙度、波纹度、形状误差等表面形貌特征,国内外许多学者对表面滤波技术进行了深入研究[2,3]。目前,表面滤波的方法主要有多项式拟合法、2RC滤波法、高斯滤波法、稳健高斯回归滤波法、样条滤波法和小波滤波法等。多项式拟合法能保证数学处理速度,但受函数形式和多项式次数的制约,其拟合精度有限;2RC滤波器能有效分离出高频和低频部分,但存在非线性相移的缺陷;高斯滤波法能有效分离表面成分且无相移,但存在边界效应和稳健性差的问题;稳健高斯回归滤波法解决了高斯滤波法存在的边界效应和稳健性差的问题[4];样条滤波法能有效提取大曲率曲线且没有边界效应,但其频率传输特性较差[5,6,7]。上述滤波法都缺乏多尺度性能,小波滤波法则具有多尺度性能但存在最优小波函数选择的问题[8]。针对目前表面滤波法所存在的问题,本文提出了一种基于多分辨率奇异值分解(MSVD)的表面滤波方法。该方法不存在边界效应和最优小波函数选择的问题,并且具有多尺度性。
1 一维多分辨率奇异值分解原理
对于任何一维信号X=(x1,x2,…,xN),构造矩阵
对矩阵H进行奇异值分解[9](singular value decomposition,SVD)得到:
其中,δ1≥δ2,u1、u2、v1、v2分别为SVD分解后得到的列向量,H=δ1u1v1T+δ2u2v2T。记H1=δ1u1v1T和H2=δ2u2v2T,H1对应的是大奇异值,反映的是H的主体成分;H2对应的是小奇异值,反映的是H的细节成分。令Hj为j尺度下的构造矩阵,Hj,1为Hj的主体成分,Hj,2为Hj的细节成分;Aj为j尺度下信号的主体成分,Dj为j尺度下信号的细节成分。由式(2)可知:
一维MSVD递推分解方式如图1所示[10]。
其中,Aj→Aj+1、Dj+1算法如图2所示。
由式(2)~式(4)可知:Aj=Aj+1+Dj+1,Aj和Aj+1都为N维向量,所以一维MSVD不会存在边界效应。对一维离散信号进行MSVD可知,MSVD与小波分析具有一定的差别,小波分析的分解层数是有限的,而MSVD不受分解层数的限制。文献[11]证明了{Dj}构成近似公比为0.5的等比数列,所以{Aj}是收敛的,则信号通过SVD迭代分解最终得到的结果是稳定的。
2 二维多分辨率奇异值分解的构造
根据一维MSVD分析理论,本文结合二维小波分析中Mallat算法来构造二维MSVD,其算法如图3所示。其中,Ajx为对矩阵Aj进行行SVD分解得到的主体成分,Djx为对矩阵Aj进行行SVD分解得到的细节成分,AjxAjy为对矩阵Ajx进行列SVD分解得到的主体成分,AjxDjy为对矩阵Ajx进行列SVD分解得到的细节成分,DjxAjy为对矩阵Djx进行列SVD分解得到的主体成分,DjxDjy为对矩阵Djx进行列SVD分解得到的细节成分。
Aj→Aj+1、D1j+1、D2j+1、D3j+1的算法如图4所示。其中,Aj为j尺度下信号的主体成分,通过一维H矩阵的构造方法来构造Aj进行行递推分解的构造矩阵Hjx和列递推分解的构造矩阵Hjy;Hjx,1、Hjx,2分别为Hjx的主体成分和细节成分。Hjy1、Hjy2分别为通过行递推主体部分和细节部分的列构造矩阵。
令
式中,M为二维离散信号的行数;N为二维离散信号的列数。
构造Aj进行行递推分解的构造矩阵Hjx,并进行SVD分解如下:
其中,int()为向下取整函数。且i=1,2,…,2 M;k=1,2,…,N-1。
对Hjx进行SVD分解如下:
令
令
其中,i=1,2,…,2 M;k=1,2,…,N-1。
通过式(3)和式(4)得出Ajx和Djx:
根据Axj和Dxj构造进行列递推分解的,并分别进行SVD分解如下:
分别令
其中,i=1,2,…,M-1;k=1,2,…,2 N。
分别令
其中,i=1,2,…,M-1;k=1,2,…,2 N。
通过式(3)和式(4)得出Aj+1、D1j+1、D2j+1和D3j+1如下:
由式(2)~式(4)和图3、图4可知:
由对一维离散信号的MSVD分解结果可知,一维离散信号通过SVD迭代分解最终得到的结果是稳定的,且不存在边界效应。二维离散信号的SVD迭代分解实质上是通过行分解和列分解迭代的形式进行的,因此二维离散信号MSVD最终得到的结果是稳定的,且不存在边界效应。
3 基于MSVD的表面滤波法
3.1 表面形貌滤波的MSVD模型
令z(x)为二维表面轮廓信号,z1(x)为二维表面形状误差,z2(x)为二维表面波纹度,z3(x)为二维表面粗糙度;z(x,y)为三维表面轮廓信号,z1(x,y)为三维表面形状误差,z2(x,y)为三维表面波纹度,z3(x,y)为三维表面粗糙度。二维表面形貌滤波的MSVD模型可表示为
其中,j为分解层数,z1(x)+z2(x)为表面评定基准线,Aj、Dk的求解见图2。三维表面形貌滤波的MSVD模型可表示为
其中,z1(x,y)+z2(x,y)为表面评定基准面,Aj、Dk1、Dk2、Dk3的求解见图4。
3.2 MSVD模型分解层数的确定
MSVD具有多尺度性,因此可在不同尺度下对表面特征进行分离和提取。对于具有多尺度性的小波分析应用于提取表面形貌评定基准时,存在小波分解次数的确定问题,同理,MSVD用于提取表面形貌评定基准也存在分解层数确定的问题。为此,本文提出了一种基于最大临近层差值条件的方法来确定分解层数。定义二维表面和三维表面最大一维相邻层差值:
其中,Δ(j)为二分递推的第j层Aj与第j+1层Aj+1的最大差值。设定最大邻层差值为Δ,当j≥k时,满足Δ(j)≤Δ,则确定递推分解终止层数为k。因为{Aj}是收敛的,所以必然存在k,使得j≥k时,Δ(j)≤Δ。
3.3 基于MSVD滤波法的计算复杂度
设一维测量信号为X1×N,T11、T12、T13分别为Aj→Hj、Hj→Hj,1和Hj,1→Aj+1的计算复杂度,k为分解层数,则一维MSVD滤波法的计算复杂度T1为
设二维测量信号为XM×N,T21、T22、T23、T24、T25、T26分别为Aj→Hjx、Hjx→Hjx,1、Hjx,1→Ajx、Ajx→Hjy1、Hjy1→Hjy1,1、Hjy1,1→Aj+1的计算复杂度,k为分解层数,则二维MSVD滤波法的计算复杂度T2为
4 仿真计算和实验验证
4.1 一维MSVD的实测实验
通过光切显微镜测量长方形铜块表面轮廓,该实验的数据评定长度为ln=7.5mm,采样间隔为Δx=0.05mm,采样点数为150,如图5所示。对该实验数据进行一维MSVD,设定迭代终止条件Δ=0.05,由图6可知分解层数可选17。
为了验证该方法的可行性与正确性,本文分别采用了高斯滤波法和高斯回归滤波法来提取实验数据的表面基线(图5、图7)。图8为基于一维MSVD滤波法得到的表面基线。从图5、图7、图8可知三种滤波方式得到的表面基线有较好的一致性,并且本文算法不会存在高斯滤波法所存在的边界数据丢失的问题。
4.2 二维MSVD仿真实验
为了说明二维MSVD对三维表面滤波法的有效性,本文以多峰函数peaks(200)作为一个采样点数为200×200的三维基准面,同时在该基准面上加白噪声作为模拟表面,如图9所示。运用基于二维MSVD滤波法对模拟表面进行滤波处理,设置迭代终止条件Δ=0.05,得到的基准面如图10所示。同时本文还分别采取了高斯回归滤波法和双树复小波滤波法来提取基准面,如图11、图12所示。
通过对仿真数据进行对比发现:从基准面的形貌来说,三种滤波结果具有较好的一致性。不过滤波性能还是有一定的差异:采用高斯回归滤波时,得到的基准面不够自然平滑,与理想基准面有一定的差异。采用双树复小波滤波法和基于二维MSVD滤波法得到的基准面自然光滑且不存在边界效应。
为了进一步评价不同滤波法对仿真信号的滤波效果,本文采用信噪比和最大误差两个评价指标来评定。其中二维信噪比RSNR为理想基准面功率和滤波后基准面的噪声功率之比:
最大误差Δmax为理想基准面与滤波后基准面的最大误差:
其中,z为标准基准面的高度矩阵,z1为滤波后基准面的高度矩阵。通过不同滤波法提取基准面的评价指标如表1所示。
从表1可以看出,高斯回归滤波法信噪比最小和最大误差最大,从而评价结果最差;基于二维MSVD滤波模型得到基准面的信噪比最大并且最大误差最小,从而评价结果最优。
4.3 二维MSVD实测数据分析
采用英国Taylor-Hobson公司的CCI型白光干涉式表面测量仪进行数据分析,其最小采样间距为0.078μm,垂直分辨率为0.01nm,实验样品为某磨削工件,采样点数为256×256,采样间距Δx=1μm,截止波长选择为λc=80μm,实测表面如图13所示。为了验证方法的可行性与正确性,本文分别采取高斯回归滤波法、双树复小波滤波法和二维MSVD滤波法对磨削样品进行滤波处理,结果如图14~图16所示。
由图14~图16可知,三种滤波法能有效地提取表面基准面,并且基准面的形貌具有较好的一致性。但高斯回归滤波结果受异常信号的影响有明显的凸峰和凹谷,而双树复小波滤波法和基于二维MSVD滤波法具有多尺度性,可通过选取较大尺度来抑制凸峰和凹谷,但双树复小波滤波法相邻尺度滤波结果差异性较大,并且分解尺度有限。为保证滤波结果的准确性,双树复小波分解尺度可根据截止波长或能量守恒法来确定[11,12,13],由采样间距和截止波长计算双树复小波分解层数为5,其滤波结果见图15。基于二维MSVD滤波法因细节成分是公比小于1的等比数列,所以相邻尺度滤波结果差异随分解尺度的增大而减小。可通过改变最大差值终止条件来微调表面基准面,这是双树复小波滤波法所不具有的。本文设置最大差值终止条件为0.05,其滤波结果如图16所示,从图16可知基于二维MSVD滤波得到基准面存在不明显凸峰和凹谷,在一定程度上能抑制异常信号的影响。
5 结束语
多分辨分析方法 篇3
股票市场是涉及金融、经济、政治、社会以及股民心理等诸多影响因素的复杂的动力学系统, 其变化过程具有非线性、混沌性、长期记忆性等特点[1~2]。Peters等[3]指出金融市场, 包括股票市场, 是由不同投资时间水平的交易者组成, 如短期、中期和长期交易者等。
一、EMD方法的引入
近年来, 小波变换 (Wavelet Transformation, WT) 理论在股票市场系统变量的多时间尺度分析与建模中取得了丰富的成果[4,5]。小波变换在时域和频域都具有良好的多分辨率分析能力, 被誉为数学显微镜。但小波变换实质上是一种窗口可调的傅立叶变换, 其小波窗内的信号必须是平稳的, 因而没有从根本上摆脱傅立叶分析的局限, 小波变换虽然能够在频域和时域内同时得到较高的分辨率, 但仍然存在一定的限制, 这种限制通常会造成很多虚假的谐波, 且小波基函数的选择对小波分解结果有显著的影响[6]。针对小波变换的不足, 1998年, Huang等人提出了一种全新的多分辨率信号分析方法—经验模态分解[7] (Empirical Mode Decomposition, EMD) 。EMD是基于信号局部特征时间尺度, 从原信号中提取本征模态函数 (Intrinsic Mode Function, IMF) 。在线性框架下基于EMD得到的Hilbert谱与小波谱具有相同的表现特性, 而Hilbert谱在频域和时域内的分辨率都远高于小波谱, 依此得到的分析结果可以更准确地反映系统原有的物理特性。由于EMD方法比小波变换有更强的局部表现能力, 所以在处理非线性、非平稳信号时, EMD方法是一种更有效的方法[8], 而金融时间序列 (如股价、股价指数、收益率等) 就是一类典型的非线性、非平稳时间序列。
二、EMD的基本理论和方法
EMD方法中的本征模态函数 (IMF) 要满足两个条件: (1) 在整个数据范围内, 极值点和过零点的数量必须相等或至多相差1; (2) 在任何点处, 所有极大值点形成的上包络线和所有极值点形成的下包络线的平均值始终为零。EMD分解基于如下前提: (1) 被分解的信号至少有两个极值点:一个极大值点和一个极小值点; (2) 局部特征时间尺度定义为信号中两临近极大值点或极小值点的时间间隔; (3) 若信号中没有极值点, 但包含一些拐点, 可以先对信号进行若干次微分, 使极值点显露出来, 最后对分解得到的分量进行积分得到最后结果。
EMD分解方法的基本思想是:加入待分解数据序列的极大值或极小值数目比上跨零点 (或下跨零点) 的数目多2个或2个以上, 则该数据序列就需要进行平稳化处理。首先, 利用三次样条函数把序列x (t) 的局部极大值和局部极小值点分别拟合成x (t) 的上包络线和下包络线, 然后计算两包络线的均值m1。再从原数据序列x (t) 减去m1, 即可得到一个移除低频的新数据序列:h1=x (t) -m1 (1)
通常, h1并不是IMF分量, 为此需要对h1重复以上处理过程以进行k次筛选, 直到所得到的平均包络趋于零为止, 此时所得数据为:h1k-h1 (k-1) -m1k (2)
式中, h1k为第k次筛选所得的数据, h1 (k-1) 为第k-1次筛选所得的数据。利用限制标准差SD的值来判断每次筛选结果是否为IMF分量, SD定义为:
其中T为数据序列的长度。
EMD分解时的限制标准差SD的值一般取在0.2~0.3之间, 即满足0.2
当h1k满足SD的要求时, 令c1=h1k即可得到信号x (t) 的第一个IMF分量, 它代表了原始信号序列中最高频率的组成成分。从原始数据序列x (t) 中减去第一个IMF分量c1, 就得到一个移除高频组分的差值数据序列:r1-x (t) -c1。若r1中仍包含x (t) 的较长的局部特征时间尺度信息, 可将r1作为待分解信号, 再重复式 (1) ~式 (3) 的过程, 直至所剩信号r1中的信息对所研究的内容意义很小或已是单调函数时停止分解运算, 此时的rn代表原始数据序列的趋势或均值。至此, 便得到了信号x (t) 的一系列IMF分量:c1, c2, …, cn且r1-c2=r2, r2-c3=r3, …, rn-1-cn=rn。原始数据序列即可由这些IMF分量及一个均值或趋势项表示：
由于每一个IMF分量是代表一组特征尺度 (频率) 的数据序列, 因此EMD分解实际上就把原始数据序列分解为各种不同特征波动的叠加, 每一个IMF分量既可以是线性的也可以是非线性的, 且每个IMF分量都有实际的物理背景相对应。
三、日收益率时序的EMD分解
(一) 基本资料
沪深300指数于2007年4月8日正式发布, 是由上海和深圳证券市场中选取300只A股作为样本, 其中沪市有179只, 深市121只。样本选择标准为规模大、流动性好的股票。沪深300指数样本覆盖了沪深市场六成左右的市值, 具有良好的市场代表性, 能够反映沪深市场股价变动情况和主流投资动向。本文分析选用的基础数据为沪深300指数的每日收盘指数时序, 样本区间为2005年4月8日至2007年9月6日日, 共计589个交易日数据, 如图1所示。
指数日收益率Rt定义为:,
其中Pt+1为第t+1个交易日的收盘指数, Pt为第t个交易日的收盘指数, 上述收盘指数时序对应的日收益率时序如图2所示。
(二) 日收益率时序的EMD分解
对图2所示的日收益率时序进行EMD分解, 取SD的值为0.25, 并采用边界延拓法来处理分解时的边界问题[9], 结果见图3至图10, 其中分别包含有7个IMF分量 (图3~图9) 和一个趋势项Res分量 (图10) 。从中, 可以得出以下重要结论: (1) 日收益率时序可以分解为7个具有不同波动周期的分量和1个趋势分量, 反映了沪深300指数日收益率变化的复杂的多时间尺度性、多层次性; (2) 对日收益率时序来说, 第一个本征模态函数IMF1是振幅最大, 频率最高, 波长最短的一个波动, 依次下去的其他本征模态函数振幅逐渐减小, 频率逐渐降低, 波长逐渐变大; (3) IMF1分量以准2个交易日 (以下简称日) 为波动的主要周期, 次要波动周期有准3日、准4日和准5日等, 且越长的周期, 其出现次数越少;IMF2分量的波动周期为准4~5日;IMF3分量的波动周期主要以准15日为主;IMF4分量的波动周期以准28日波动周期为主, 以准40日波动周期为辅;IMF5分量以准70日波动周期为主, 嵌套着准80日波动周期;IMF6的波动周期为准140~190日;IMF7分量波动周期为准240日;Res分量反映的是日收益率时序的整体变化趋势, 自2005年4月初至2007年9月初, 沪深300指数所反映的中国股市的整体收益率呈上升趋势, 上升幅度达19.19%, 这与同期中国股市的牛市状态是相符合的。 (4) 从图3~图6可以看出, IMF1~IMF4分量的波动幅度具有较明显的“正的持续性”, 即具有长期记忆性, 反映在日收益率序列波动幅度的变化上, 用平均的观点来看, 表明IMF1分量过去的一个较大幅度的波动趋势意味着将来的一个较大幅度的波动趋势, 过去的一个较小幅度的波动趋势意味着将来的一个较小幅度的波动趋势, 除非发生改变这种趋势的关键事件, 同时也暗示IMF1波动分量序列表现出一定的非高斯性 (即非随机性) 。已有相关研究证实了我国股市收益率存在长期记忆性[10]。
结语
经验模态分解 (Empirical Mode Decomposition, EMD) 作为一种全新的信号分析理论, 是近年来对以傅立叶变换为基础的线性和稳态谱分析的一个重大突破, 比小波变换具有更强的局部特性和适应性。本文运用EMD方法对根据沪深300指数的日收盘指数时序计算得到的其日收益率时序进行了多分辨率 (多时间尺度) 分解, 结果表明其波动分量分别具有准2日、准4~5日、准15日、准28日、准70日、准140~190日、准240日等波动周期, 日收益率的长期记忆性主要表现在其短期波动分量 (准2日、准4~5日、准15日和准28日等) 中, 这揭示了中国股市系统变量运动所具有的复杂的多时间尺度性。必须指出的是, 股市系统是一个具有很强时变性的动力学系统, 本文的结论只是就观测时段内数据计算分析所得。
基于精确的多分辨率分析能力, EMD理论作为一种崭新的非线性、非平稳信号分析理论在金融系统变量时间序列的分析、预测和建模中具有广阔的应用前景。
参考文献
[1]李亚静, 何跃, 朱宏泉.中国股市收益率与波动性长记忆性的实证研究[J].系统工程理论与实践, 2003, 23, (1) :9-14.
[2]徐梅, 张世英.基于小波变换的长记忆随机波动模型估计研究方法[J].中国管理科学, 2006, 14, (1) :7-14.
[3]Peters E.E.Fractal market analysis:applying chaos theory to investment and economics[M].New York:John Wiley Sons, 1994.
[4]宿成建, 刘星, 刘礼培, 等.应用小波分析方法研究沪深股市的溢出效应[J].系统工程学报, 2004, 19, (1) :99-103.
[5]侯守国, 张世英.基于小波分析的股市高频互相关研究[J].中国管理科学, 2006, 14, (3) :1-6.
[6]Tewfiki.A.H.On the optimal choice of a wavelet for signal representation[J].IEEE Trans Information Theory, 1992, (2) :747-765.
[7]Norden, E.H., Shen, Z., Long S.R., et a1.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-station ary time series analysis[J].Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, 1998:454-955.
[8]Huang N.E.A new view of nonlinear water waves-the Hilbert spectrum[J].Ann.Rev.Fluid.Mech., 1999, 31 (1) :417-457.
[9]瞿伟廉, 程磊.应用径向基函数神经网络处理EMD方法中的边界问题[J].华中科技大学学报:城市科学版, 2006, 23, (4) :1-4.
多分辨分析方法 篇4
关键词：多小波变换,图像融合,超声成像,图像空间分辨率,插值
无损检测成像技术在现代工业中有很重要的用途,超声成像技术是其中一种令人瞩目的新技术[1]。目前,在数据量巨大,时效性的要求下,超声成像技术应用存在的主要问题是成像图像分辨率较低,而图像的空间分辨率是评价图像质量的一项关键性指标[2],直接决定其应用价值。
超声成像无损检测目的是确定对象内部缺陷。改变超声图像不清晰,分辨率不高的缺点,主要技术实现采用两种途径:一是重新提升改进硬件(如采用电子聚焦、超声聚焦、连续焦点图像技术、变孔径技术等方法来提高技术指标),但受限于技术发展,同时开发周期长,成本高;二是利用数字信号处理方法,提高算法精度。文献[3][4]利用分形理论双毯分析法和正则化变分法及偏微分方程滤波技术提高超声图像的分辨率。采用图像融合提高图像分辨率常见于遥感图像处理[5],对于材料内部缺陷及焊缝质量检测的超声图像研究目前还很少。文献[6]采用基于局部协方差特征的边缘自适应图像插值方法来提高C扫描图像的分辨率。文献[7]采用基于插值和小波图像融合相结合区域能量最大法提高超声图像的分辨率。
针对超声图像的应用特点,利用多小波理论在图像处理的优势,采用多小波变换图像融合来增强超声图像分辨率,分解后低频高亮各子子带采用加权平均法,高频各子子带采用区域能量加权融合方式。通过与小波变换区域能量最大法[7]比较及多项指标结果进行评价,结果表明,该方法对改善超声图像分辨率有一定增强效果。
1 正交离散多小波变换[8]
在数字图像处理中,正交性能保持信号能量,对称性(线性相位)使信号的边界易于处理,但在实数域中同时满足正交性、对称性、紧支撑的单小波基是不存在的。1994年T.N.T.Goodman,S.L.lee等在小波理论基础上,突破性提出了多小波概念。此后,J.S.Geronimo,D.P.Hardin,P.R.Massopust构造了著名的GHM多小波基。多小波保持了单小波良好的时域和频域局部特征,同时克服了单小波的缺点,将正交性、对称性、光滑性、紧支撑性、高阶消失矩结合在一起。
正交离散多小波变换中,
离散多小波变换对应一个矢值多滤波器组,有r×r矩阵“子带”,它需要r路数据流,所以变换前需要对输入数据进行预处理[9]。该文中采用4系数对称多滤波器组,“低通”滤波器为4个2×2的矩阵G(k)组成,选取r=2,预处理采用基于逼近的“appro”方法[8],边界不处理。
Jiang Q.T于1998年构造出SA4多小波基,它是一组对称—反对称多小波组[10],它的逼近阶数2、支撑长度为[0、3],同时保持正交性和对称性,光滑性好于GHM多小波基和CL多小波基。
其中四个系数矩阵:G0;G1;G2;G3值见文献[12]
其中四个系数矩阵:H0;H1;H2;H3值见文献[12]。
图像多小波变换在第一级分解前首先完成一个预处理步骤。不同于单小波变换,多小波由于存在多个尺度函数,单小波变换生成的一个子带,在多小波变换后且是r×r个子子带。一般而言,L级变换可得到r2(3L+1)个子子带。图1是二维离散多小波分解示意图。
2 超声成像图像融合技术
2.1 提高分辨率的图像融合原理
图像融合是对多源信道采集到的同一目标的图像数据进行计算机图像处理。首先将目标两个或多个图像进行空间定位配准,然后应用某种融合算法将各图有机结合起来,最大限度的提取各个图的有利信息,以提高图像信息的利用率、提升原始图像的空间分辨率和光谱分辨率为目的。
从理论上讲,如果图像采集各图只含一幅低分辨率的图像的信息,要获得更高分辨是不可能的,因为信源无变化就无法获得更多信息。如果已经得到了同一场景下的变化的一组观测图像,并已知成像参数。如果这一组观测图像之间存在较小的相对位移,而且图像足够丰富,就是信源有变化,可以认为它从整体上包含更高分辨率图像的信息。这一组图像之间所含信息是互补的,又是冗余的。因此,基于序列图像来提高图像分辨率本质是一种图像数据的融合问题,属于数据级融合范畴[11],也称像素级融合。
2.2 基于多小波变换的图像融合[12]
采用像素级的融合方法可以保留更多的信息。基于多小波变换的图像融合的方法首先是对已配准空间的图像分别进行预滤波,然后作多小波分解,将图像分解成频域内不同频率段的分解系数,构成一系列子图像(矩阵),代表原图像的各个特征分量。融合处理可以根据不同的特征与细节采用不同的融合方法以达到最佳效果。图2给出了此图像融合算法的示意图(L=1)。
2.3 超声图像融合前插值
在一定采样率的条件下,为获得更多信息,超声成像技术应用微位移采样技术,使图像之间存在较小的相对位移,如C超扫描图像,微位移半个像素距离内采样,移动精度到达1/10单个像素距离,这样图像足够丰富,从整体上记录了更高分辨率图像的信息。
微位移半个像素距离采样实际使采样率提高一倍,对应图像采样也是放大一倍。所以在微位移采样获得的图像进行图像融合前,需要对各单幅图像进行放大,即横向纵向放大一倍。在保证信息熵的要求下,图像放大算法采用数学插值。常用方法有帧内插值、三次样条插值、双线性插值等。该文采用双立方插值法[13]。双立方插值法被插值点的值是由4x4个邻近的点求的,比双线性插值更平滑。
3 区域能量加权超声图像融合方法
3.1 融合流程
根据区域能量融合的思想,多小波变换后,低频高亮4子子带部分采用加权平均法,高频各子子带采用区域能量匹配加权融合方法。融合流程图见图2所示。融合步骤如下:
1)应用双立方插值法,对图像A’,B’分别进行图像插值,得到插值图像A,B;
2)将图像A,B作“appro”预滤波和SA4基多小波一级分解,得到不同尺度的低频分量和高频分量;
3)对获取的低频高亮4个子子带系数采用平均算子,融合值通过对应位置的值加权平均获得;
上式:Ak(2j,x,y),AAk(2j,x,y),ABk(2j,x,y)分别代表融合结果图像、源图像A和源图像B在基带点(x,y)位置上的值,k1,k2是两个权值,不同的数据源采用不同的权值。在本方法中权值采取k1=k2=0.5。
4)获取的高频各子子带系数采用区域能量加权融合方法。
5)对融合后的各高频子带和各低频子带系数进行SA4基多小波逆变换,得到融合结果图像。
3.2 区域能量加权融合规则
用L×K矩阵窗口对图像求能量运算,计算结果作为窗口中心像素选择的方法:判断区域能量匹配度,加权计算像素的值作为融合图像的对应点的像素值。
高频各子子带计算在中心像素近邻区域L×K(该文选3×3)的能量E,见式(7),用M来表示两图像在该区域的匹配度,其式(8);当M大于给定阈值a时,则进行加权运算融合各原图,权值W由式(9)确定,否则,直接选择区域能量大的区域的中心像素,最后的融合值按式(10)计算。
其中Ej,A,Ej,B按式(7)计算。
定义一个匹配度阈值a(a通常取0.5~1,该文取0.75)
式中j表示2j尺度,D(n,m)表示对应子子带(n,m)点的多小波系数,DA和DB表示待融合的源图像,DC表示融合后的图像。由于源图像取自超声成像,可以保证观察缺陷对象集中在中央,多小波重构时所有子子带不再作边沿处理。
4 试验结果及分析
试验选择在Matlab平台上完成,试验图像是实际的缺陷超声图像,由加拿大RD/TECH公司Omniscan超声相控阵探伤仪扫描获得,如下图5所示,图5a中红色等深色、类似椭圆形区域是C-扫描显示的空洞缺陷,其降倍采样图像如图5b和5c所示,其中图5c与图5b在x,y方向上的相对位移都为+0.5像素。
试验1:对图5b和5c采用双立方插值放大2倍,采用SA4多小波一级分解,预处理采用“appro”方法,低频高亮4个子子带系数采用平均算子,高频子子带计算能量E选用近邻区域的窗口大小为3×3,处理后的融合结果分别如图4所示。
试验2:对图5b和5c采用双立方插值放大2倍,采用9/7单小波一级分解,方法见文献[7]。处理后的融合结果分别如图5所示。
从主观评价(目视效果)看,图4好于图5,接近原图,但无法准确判断。客观评价指标采用信息熵(H)、均方误差(MSE)、光谱扭曲度(D)和峰值信噪比(PSNR)度量,计算结果如表1所示。
从客观评价指标上分析上,应用本方法的均方误差和光谱扭曲度减小,同时峰值信噪比增大,说明更多地保留了原图像的有用信息,优于单小波方法。信息熵基本无变化可通过多级多小波分解进一步改善。
5 结论
1)基于插值和多小波变换图像融合的方法对于改善超声图像的分辨率具有一定增强效果。有利于进一步精确定位缺陷边沿。
2)利用多小波变换理论的优势,更良好的时频局部化特性,同时正交性、对称性、光滑性、紧支撑等集一起,应用到超声图像处理效果优于单小波。
多分辨分析方法 篇5
多分辨分析的思想就是先在能量有限函数空间的某个子空间中建立基底, 然后利用简单的伸缩与平移变换, 把子空间的基底扩充到中。
L2R中的一列闭子空间{Vj}j∈Z称为多分辨分析, 如果{Vj}j∈Z满足下列条件:
(5) 存在g (x) ∈V0使得{g (x-n) M∈Z}构成V0的Riesz基, 即
b) 存在0πA≤Bπ∞, 使得对任意{cn}n∈Z∈l2 (Z) 有:
其中, g (x) 称为尺度函数。
条件 (l) 称为一致单调性, 表明子空间列{Vj}j∈Z是嵌套的。条件 (2) 称为渐近完全性, 表明通过增大j, L2 (R) 中的每一个函数能够用它在Vj中的投影Pjf非常接近所希望的逼近;反之, 通过减小j, 投影Pjf能够具有任意小的能量。条件 (3) 称为伸缩规则性, 表明空间列{Vj}j∈Z中任一空间Vj的基可由其中任意空间Vj的基经过伸缩与平移变换得到。条件 (4) 称为平移不变性, 表明函数f (x) 在同一子空间Vj中波形在平移后不变化, 即f (x) ∈Vj圯f (x-2jk) ∈Vj坌j, k∈Z。条件 (5) 称为形Riesz基的存在性, 结合条件 (3) 、 (4) 可知, {g (2-jx-k) |k∈Z}构成Vj的Riesz基。因此又称g (x) 为多分辨分析的生成元。可以证明, 存在函数谆 (x) ∈V0使得{谆 (x-k) |k∈Z}构成V0的标准正交基。数列{hk}或与之等价的函数H (覣) 可以完全决定一个多尺度分析。
2 Mallat分解与重构算法
1988年Mallat受塔式算法的启发, 在多分辨分析的指导下建立了Mallat算法, 它的作用可与FFT在Fourier变换中的作用相提并论, 具体算法可描述如下:
设V0是给定的多分辨分析尺度空间, 尺度因子m=0, 相应的尺度函数为:
小波函数为
由V0=V1⊕W1, 若f (t) ∈V0, 则正交小波分解为
当m=1时, 则:
当尺度因子m=1转到m=2时, 相应于f0→f1, 需补充细节信息d1, 则f0=f1+d1, 类似地则有fm-1=fm+dm。
由此f0=f1+d1可得
为建立φ (t-n) 与φ (t-j) 的系数之间的关系, 代换角标, n用k代之, j用n代之, 则有
由此可得
由此推广至系数表达的一般式
该式就是由小波系数anm和bnm重构函数f (t) 的正交展开式系数anm-1的公式。
函数f (t) 的小波展开系数〈f, φm, n〉与〈f, ψm, n〉用双尺度差分公式表示可按下式计算
同理:
由上述计算可知, Mallat算法本质上不需要知道尺度函数φ (t) 和小波函数ψ (t) 的具体结构, 只由系数就可以实现f (t) 的分解与重构, 因此, 称为快速小波变换。
令:
则可以把滤波处理表示成简洁的形式, 即
由此, 完全重构a0的条件是HH*+GG*=1, H和G是一对共轭正交滤波器组 (Conjugate Quadrature Filters) 。
Mallat算法描述如图1所示。利用小波进行图像压缩的例子如图2所示。
参考文献
[1]Mallat S.A theory for multiresolution signal decomposition:the wavelet representation[J].IEEE Trans on PAMI, 1989 (7) .
[2]Mallat S.Multifrequency channel decomposition of images and wavelets[J].IEEE Trans.ASSIP, 1989 (2) .
[3]Mallat S.A Wavelet Tour of Signal Processing[M].London:Academic press, 1998.
[4]刘贵忠, 邸双亮.小波分析及其应用[M].西安:西安电子科技大学出版社, 2000.
多分辨分析方法 篇6
2014 年世界卫生组织( WHO) 发布的《世界癌症报告》显示,肺癌的发病率居首位,是最普遍和最致命的癌症[1]。对于肺癌来说,PET图像对肺癌的早期诊断有较高的灵敏性,但由于图像空间分辨率较低不能精确定位病灶; CT图像的高解析度可以为肺癌患者的临床诊断和分期问题提供病灶的精确解剖结构,但因软组织分辨率较低在定性诊断上有很大限制。PET/CT融合图像可同时呈现不同模态图像的信息,提高病变部位的可辨识度[2],同时进行病灶的精确定位以及放疗靶区的勾画[3],对肺癌的早期诊断、临床分期、病灶定位、制定治疗方案以及疗效评估等方面都有重要的临床应用价值。目前,PET/CT图像融合的研究主要集中于临床诊断方面,Hari Mukundan等[4]比较了MRI和PET / CT对患者头颈部鳞状细胞癌治疗后的诊断性能,结果表明PET/CT对术后评估的效果更好。
针对PET/CT图像融合算法的讨论相对较少,且以二维小波变换为主。在小波变换的基础上Do等[5]提出了轮廓波变换( Contourlet transform) 理论,Contourlet变换是一种多分辨率、局域性和多方向的稀疏表示方法,文献[6]已经证明相比于小波变换以及改进的分层离散余弦变换,Contourlet变换能更稀疏地表示图像。但由于Contourlet变换不具有平移不变性,会产生频谱混叠现象,而基于非下采样Contourlet变换( NSCT) 理论的提出解决了这一问题[7]。NSCT去除了下采样的环节,具备良好的方向性及平移不变性,可有效表达出图像的纹理特征,避免了Gibbs现象的发生[8,9]。但是NSCT在图像融合中的计算复杂度高,高频图像在多方向上的分解会产生大量的冗余数据,本文将多分辨率变换NSCT与压缩感知相结合应用与PET / CT图像融合,由于多分辨率NSCT分解后的高频子带满足压缩感知理论中稀疏先验条件,可以通过少量的观测数据恢复出原始信号[10],因此本文将医疗成像设备获取到的源图像经过多分辨率NSCT变换后进行压缩测量再做融合处理,提高了PET/CT解剖型和功能型图像融合的质量,降低了医学图像融合在网络传输过程中的数据量,为移动医疗、智慧医疗提供技术支持。
本文提出一种基于多分辨率变换NSCT和压缩感知的肺癌PET/CT图像融合算法,针对移动医疗背景下医学图像融合信息交互的局限性问题,为克服异地传输对带宽的影响,考虑到PET与CT在成像机理上的不同,对低高频子带应结合各自特点采用相应的规则进行融合。本文在NSCT变换域对多方向的高频系数采用基于压缩感知的图像融合方法,低频系数的融合规则采用自适应局部区域能量加权的方法。实验结果表明,将多尺度NSCT变换和压缩感知的理论运用到PET / CT图像融合的工作中,降低了高频信号采样工作量,融合后的图像有较好的视觉效果,并且在客观指标评价上较其他算法有一定的提高。
2 非下采样Contourlet变换
Cunha等[7]在2006 年提出了非下采样Contourlet变换理论( Non - subsampled Contourlet Transform,NSCT) 。在多分辨率、多方向、各向异性的Contourlet变换基础上具有了平移不变性,并且消除了Gibbs现象,是一种冗余的超完备图像稀疏表示方法[11]。非下采样Contourlet变换是由非下采样金字塔滤波器组( NSP) 和非下采样方向滤波器组( NSDFB) 共同组成[12],源图像首先经过塔形滤波器进行多尺度,多分辨率的变换,得到与源图像大小一致的1 个低频子带和1个高频子带,高频子带再经过方向滤波器进行i级方向变换( 方向级数可自行设置) ,得到2i个高频图像,这是NSCT变换层级k = 1 的分解过程。随着层级的增加低频子带会在NSP上进行迭代操作,若图像经过层级为k的NSCT变换后将得到1 个低频图像和个多方向的高频图像,其中k为变换层级,2ik为对应层级的方向分解数。
3 压缩感知理论
图像的可稀疏表示是压缩感知的基础条件[13],从数学角度而言,假设一大小为N × N的图像f ( f ∈RN × N) ,变换基为 Ψ,对f进行变换,则f可表示为
式中: f是图像在时域的表示; α 是图像在 Ψ 域的表示。若 α 的非零值K << N,或者 α 经排序后呈指数级衰减并趋近于零,则认为图像是稀疏的。
在稀疏先验条件下,构造一个与变换基 Ψ 不相关的测量矩阵 Φ 对信号进行线性投影得到感知测量值y,实现信号从N维降为M维,用以精确的重构信号或者图像[14],采样过程可表示为
式中: Φ 为M × N大小的测量矩阵; y是图像f在测量矩阵 Φ 下的线性投影值,且M < N。测量矩阵的构造要使得测量值的数量尽可能少而且要与稀疏变换矩阵不相关,图像经过稀疏变换和测量矩阵之后得到测量值实现降维,还要利用重构算法恢复到高维图像,其数学表达式为
重构算法是图像是否能精确重构的关键,并且该算法要求在可以精确恢复信号的同时,处理过程稳定、计算复杂度低、观测数量少[15]。
4 基于多分辨率变换和压缩感知的肺癌PET/CT图像融合算法
4. 1 算法思想
1) 首先对肺癌患者的PET及CT图像分别进行单层多分辨率NSCT变换,分解后分别得到1 个低频图像和8 个方向的高频图像。
2) 变换后的低频图像大多是逼近信号且稀疏性较差,可直接进行融合以保持图像大量背景信息。
3) 变换后的高频图像主要呈现原始图像的细节信息,相比于低频图像的稀疏度,高频图像的稀疏性更好,对高频系数进行压缩采样能得到更好的重构精度。因此,采用高斯随机矩阵对高频图像进行压缩测量,根据符合测量值物理特性的融合规则进行融合,然后利用正交匹配追踪算法重构融合后的高频图像。
4) 最后,重构后的高频融合图像与融合后的低频图像经过NSCT逆变换得到最终的融合图像。融合框架如图1 所示。
4. 2 关键技术
4. 2. 1 低频融合规则
源图像经过NSCT变换后的低频图像主要包含源图像的近似分量,稀疏性很小,如果使用随机测量矩阵对低频系数进行观测,会破坏低频系数之间的相关性,影响重构精度[16],本文选择直接进行融合。由于PET与CT在成像机理上的不同,PET与CT图像的灰度差较大,通常出现互斥特性,代谢旺盛的恶性病变组织在PET图像中表现为较暗区域,而骨骼及脏器的分布在CT图像中能得到清晰的呈现,如果采用基于像素的简单平均融合算法会降低融合图像的对比度,使重要目标信息淡化[17]。本文中,首先对两幅源图像PET( A)和CT( B) 经过单层NSCT变换的低频图像进行取极小值融合,融合后的图像再与PET的低频图像进行融合
式中: LA和LB分别为PET和CT的低频图像分解系数。在此基础上选择模糊集理论中的高斯型隶属度函数计算融合算子,通过自适应加权的方法对低频子带进行融合。高斯型函数的表现形式为
式中:c和σ2分别代表低频分解系数的均值和方差。根据式(5)计算出PET的低频图像LA与取极小值后低频图像LB的隶属值GAi,j和GBi,j(i∈m,j∈n),以此求出自适应融合算子如下
式中: ωA和 ωB分别为低频图像LA与LB的融合算子,即权重值; LF为最终融合后的低频子带系数。
4. 2. 2 高频融合规则
PET与CT图像经过NSCT变换后得到8 个方向的高频子带,包含了源图像在不同方向上的细节信息且稀疏性较好。由于压缩感知对图像信息的压缩采样是线性变换过程,高频系数转换为观测值以后仍然具有观测值越大包含图像信息量越大的特点,只是随机测量矩阵和NSCT变换矩阵之间存在不相关性,因此,高频系数的融合规则不能根据像素间的关系进行设置[18]。本文中,对高频子带分别进行压缩测量,尽可能保留源图像的边缘特征,采用线性加权的方法基于平均梯度和区域能量计算权值: 平均梯度包含了图像细节信息和纹理特征的变化特征,更适合用于高频的细节分量; 计算区域能量的方法能体现出高频测量系数的局部特征,因此,高频融合规则将测量值的共同信息和互补信息进行综合,使得重构恢复后的高频融合图像信息更加丰富。
首先构造高斯随机测量矩阵对两幅源图像的高频子带系数HA和HB进行线性测量,得到相应子带系数的测量值YA和YB。
其次,根据平均梯度与区域能量计算公式
得到高频子带的平均梯度为GradI和能量EI( I =A,B) ,根据高频子带的平均梯度计算高频子带测量值YA、YB的权重 ωYA和 ωYB分别为
最后,根据高频子带测量值的融合规则表达式
对PET和CT图像分解后的高频子带测量值进行融合得到融合后的高频子带系数YF,再采用正交匹配追踪算法对YF进行重构得到高频融合图像HF。
5 仿真实验及分析
5. 1 实验环境
硬件环境: 仿真硬件平台为Pentium( R) Dual -Core CPU E6700,3. 2 GHz,2. 0 Gbyte内存,操作系统为Windows 7。
软件环境: 软件Matlab R2012b。
实验数据: 采用一组肺癌经配准后的PET及CT图像,如图2 所示,均为灰度图像,大小均为256 × 256。
NSCT变换参数设置: 滤波层级为1,方向级数为3,其中NSP结构采用双正交小波分解,NSDFB采用梯形滤波器。
5. 2 实验结果及分析
1) 实验一: 与其他压缩感知图像融合方法的对比
为了验证非下采样Contourlet变换的优越性,将本文算法与其他压缩感知图像融合方法进行对比,对比方法分别为基于小波变换的压缩感知图像融合( Wavelet-CS) 和基于Contourlet变换的压缩感知图像融合( CT-CS) 。Wavelet-CS试验中稀疏变换矩阵选择小波正交基,融合规则是加权平均的方法; CT-CS实验中高低频融合规则与本文算法相同,3 种方法的实验中测量矩阵均为高斯随机矩阵,重构算法均为正交匹配追踪算法,采样率分别设置为30% ,50% ,70% 。表1 为本文算法( NSCT + CS) 和Wavelet-CS、CT-CS算法下得到的融合图像在客观评价指标下的数据。
融合图像的评价方法一般包括主观评价和客观评价。主观评价对于图像质量检验来说最为可靠,尤其是在医学图像融合中,对临床医生的辅助作用很关键,但是主观评价需要设备、人员组织和严格的环境条件协同作用,在实际应用中较难实现。因此,本文选择客观评价指标来度量融合图像的质量和性能,包括标准差、平均梯度、空间频率、峰值信噪比、图像清晰度、互信息、信息熵。
从表1 中的数据可以看出,在不同的采样率下,本文算法融合图像的标准差、平均梯度、空间频率、图像清晰度、互信息均高于其他两种算法,CT + CS算法下的融合图像的峰值信噪比和信息熵优于本算法与Wavelet-CS算法,可能是由于Contourlet变换会产生频谱叠加现象的原因,图像的灰度变化和边缘波动程度较大,使得图像信息量变大; 而本文算法和CT-CS方法的指标数据均优于Wavelet-CS算法,可知多尺度、多分辨率的稀疏变换能使图像较小波变换更加稀疏,且保留多方向的细节信息。
2)实验二:与其他多分辨率图像融合方法的对比
为检验压缩感知理论在图像融合中的优势,本文算法将分别与NSCT变换图像融合和Contourlet变换图像融合进行比较,本文算法的采样率定为50%,NSCT变换与Contourlet变换的低频融合规则均为加权平均法,高频融合规则均为局部区域能量最大法。如图3为融合结果图像,表2为3种融合方法在客观指标上的比较,图4为客观指标测试数据柱状图。
从融合图像的视觉效果可以看出本文融合算法将PET与CT图像各自的特征综合呈现在一幅融合图像中,患者肺部横断面的骨骼和脏器组织较为清晰,尤其是癌细胞代谢旺盛的病灶部位得到了更完整的保留,提高了病灶部位与毗邻组织的对比度,为病灶的精确定位提供高质量的影像依据。但是,相比于NSCT图像融合方法和Contourlet图像融合方法所得到的融合图像,本文算法的融合图像中骨骼亮度并不高,而且图像左上角信息标注的地方数字有些模糊,融合算法仍可做进一步改进以得到更好的视觉效果。
从表2、图3 以及图4 中可以看出,本文算法得到的融合图像的标准差为23. 279 7、平均梯度为4. 205 1、空间频率为11. 721 3、图像清晰度为4. 995 4,均远大于另两种多尺度的图像融合方法得到的结果; 互信息为1. 994 7,相对于另外两种方法提高的并不是十分明显;而峰值信噪比和信息熵的值均小于NSCT图像融合方法和Contourlet图像融合方法的结果,与实验二同理,Contourlet产生频谱混叠现象,所以峰值信噪比和信息熵为三种方法中最大; 本文算法中结合了多尺度变换和压缩感知理论,多尺度变换使图像表达更加稀疏,压缩感知对高频图像的处理更加具体细致,在综合原始图像信息能力上有较好的优势,仅用高频图像50% 的数据就能融合出高质量的PET/CT图像。
6 结论
多分辨分析方法 篇7
随着电力系统的发展,特别是非线性负载的广泛应用,电力系统的谐波情况越来越复杂,既存在频率是工频整数倍的谐波,还存在大量的非整数次谐波[1],即间谐波,对电能造成严重的污染,降低了电力系统的可靠性。因而,间谐波的分析与检测对于电力系统的监控与保护具有十分重要的意义[2]。
当前主要的谐波检测方法都是以离散傅立叶变换(DFT)为基础的,但是用DFT方法只能精确地检测出整数倍的基波频率,对于非整数次的间谐波却无能为力[3]。此外,DFT方法还存在频谱泄漏问题[4]。
绝大部分非线性负载产生的间谐波,其幅值和频率随时间变化,即具有时变特性[5]。小波是一个时间函数,它的傅立叶变换呈现为带通滤波器的频率特性,即它在时域和频域内都是局部化的[6],克服了傅立叶变换仅有频域局部化而无时域局部化特性的缺点,适用于间谐波的分析和时变谐波的检测,在电力系统谐波检测中逐步得到广泛应用。本文采用基于小波多分辨率分析方法来实现对间谐波的检测,并通过仿真来证明本文所采用方法是可行的。
2 小波变换的定义
小波变换是将信号与一个时域和频域均具有局部化性质的平移伸缩小波基函数进行卷积,将信号分解成位于不同频带—时段上的各个成分。
设ψ(t)是一个小波母函数,将ψ(t)进行伸缩和平移,记
式中:a为尺度因子;b为平移因子;ψa,b(t)为由小波母函数ψ(t)生成的依赖参数a和b的连续小波基函数。
对任意f(t)∈L2(R),f(t)的连续小波变换定义为:
在实际应用中,为了方便用计算机进行分析、处理,信号f(t)都要离散化为离散序列,a和b也必须离散化,成为离散小波变换(DWT)。通常将a和b的离散化公式分别取作a=a0j,b=ka0jb0,j∈Z。则f(t)的离散小波变换定义为:
当a0=2,1时,有b0=1时，有、称为二进小波变换,此时,对应的为:WTf
3 多分辨率分析
多分辨率分析理论是建立在函数空间概念上的,它为正交小波变换提供了数学上的理论基础。
如果对于满足多分辨率分析的一系列闭子空间{V j}j∈Z,Vj⊂Vj-1,φj,k(t)是由{V j}j∈Z产生的尺度函数,它的补空间为{W j}j∈Z,Vj=Vj+1⊕Wj+1,则{W j}j∈Z构成小波空间,相邻空间的剩余系数cj,k和小波系数dj,k满足Mallat算法:
式(4)和式(5)可看作j-1尺度空间的剩余系数cj-,1k经滤波器系数h0(n)、h1(n)进行加权求和得到j尺度空间的剩余系数cj,k和小波系数dj,k。
对于任意函数f(t)∈Vj-1,按照Mallat算法分解一次投影到Vj,Wj空间,可得到剩余系数cj,k和小波系数dj,k;将Vj空间的剩余系数cj,k进一步分解下去,可分别得到Vj+1,Wj+1空间的剩余系数cj+,1k和小波系数dj+,1k。按照这样的方法一直分解下去,可以得到任意尺度空间的分量。类似于信号的分解过程,可以从剩余系数和小波系数重建原信号,式(6)是小波变换系数的重建公式:
4 间谐波检测原理
多分辨率分析通过不断地滤除频率相对较高频带上的细节信号分量,同时保存这些分量进行信号重构。这就是多分辨率思想用于谐波检测的原理[7,8]。
由于小波变换具有良好的时频局部化特征,故可以根据不同尺度的小波变换系数的幅值来测量谐波的频率。如果能够使不同频率的谐波位于不同的频带中,就能够把包括整数次/非整数次的不同频率的谐波分离出来。因此,利用小波变换可以实现整数次和非整数次的谐波含量的测量。
与标准傅立叶变换相比,小波函数ψa,b(t)具有不唯一性。为避免小波变换时基波和谐波发生混叠,在工程应用中,必须选择一个最优小波基。结合电网中谐波的特点和谐波检测的目的,最优小波基应具有以下特征[9]:
(1)正交性。保证基波和谐波的完全分离;
(2)紧支性。保证检测的精度;
(3)对称性。避免信号分解与重构时的失真;
(4)高阶消失矩。保证有限的计算长度和计算量。
综合以上的特征,本文采用Daubechies(dbn)小波系的函数作为仿真用的小波函数。
5 间谐波检测的实现
本文将小波变换应用于电力系统的间谐波测量,对含有间谐波的信号进行正交小波分解。利用多分辨率的概念,将高尺度上的结果看作不含谐波的基波分量。为分析简单,只考虑信号中含有一个间谐波的情况。
设原始电流信号为:
其中,基频ω=50Hz。采样频率f s=4.8k Hz。很明显,原始信号中含有一个4.7次的间谐波。选取数据采样点为前1024点,仿真结果如图1所示。
图1所示波形为借助MATLAB小波分析工具箱的图形用户接口(GUI)对原始信号s(t)进行仿真所得,其中选取的小波函数是具有紧支集的正交小波Danbechies10(db10),对信号进行了5层分解。由图1可以看出,选取该小波函数便可取得比较好的结果。
在图1中,s为含有间谐波的原始信号波形图,5a为小波分解的第五层的近似信号,d1~d5为各层的细节信号。其中,4.7次间谐波体现在第五层的细节信号d5中,随着尺度的增加,时间分辨率降低,间谐波信号在低时间分辨率上被很好地重构。
小波分解后的基波分量包含在近似信号a5中,将原始信号中的基波分量s1与近似信号a5进行比较,如图2所示。可以看出,重构基波分量和原始信号中所含的基波分量的误差比较小,而此微小误差是由于db10正交镜像滤波器组幅频特性曲线尾部小部分重叠产生的[10,11]。
从仿真结果可以看出,原始信号通过小波变换可以分解为基波分量和各次谐波分量。在上述仿真中,不同的尺度具有不同的时间和频率分辨率,进而小波分解能将原始信号的不同频率成分分开,所以谐波/间谐波可以被较好地分解出来。
6 结语
本文中提出了基于小波多分辨率分析的间谐波检测方法,仿真结果表明了该方法是可行的。小波函数虽然具有良好的时频局部化性质,但在具体应用时仍存在着固有的缺陷[12],主要体现在小波窗口能量不集中或分频不到位,从而造成频谱混叠现象的发生,影响信号分析的效果。总体来说,小波变换对信号的分析灵敏度高,并且比较精确。随着小波理论的不断发展和完善,小波变换必将在电力系统间谐波信号处理中发挥更大的作用。
摘要：间谐波是电力系统中一种特殊的谐波,大量存在于电网中,对电网的危害很大,对它进行检测和分析具有重要的意义。针对傅立叶(DFT)方法只能检测出整数倍的基波频率,而对间谐波无能为力的缺点,文章提出了一种基于小波多分辨率分析的间谐波检测方法。仿真结果表明该方法具有较好的分辨率,可以很好地对间谐波进行分析,在电力系统间谐波分析方面具有经典DFT方法无可比拟的优点。
关键词：间谐波,多分辨率分析,间谐波检测
参考文献
[1]郝江涛,刘念,幸晋渝.电力系统间谐波分析[J].电力自动化设备,2004,24(12):23～25
[2]Yacamimi R.Power system harmonics.IV.Interharmonics[J].Power Engineering Journal,1996(4):185～193
[3]崔锦泰,程正兴.小波分析导论[M].西安:西安交通大学出版社,1995:23～181
[4]李圣清,彭玉楼,周有庆.一种改进型自适应谐波电流检测方法的研究[J].高电压技术,2002,28(12):325
[5]向东阳,王公宝,马伟明,等.基于FFT和神经网络的非整数次谐波检测方法[J].中国电机工程学报,2005,23(9):35～39
[6]廖瑞金,宋杨,杜林.两种数字化介损监测方法的软件仿真分析[J].高电压技术,2001,27(2):11～13
[7]高成.Matlab小波分析与应用[M](第二版).北京:国防工业出版社,2007:29～195
[8]柴旭峥,习文山,关根志,等.一种高精度的电力系统谐波分析算法[J].中国电机工程报,2003,23(9):67～70
[9]任震,等.小波分析及其在电力系统中的应用[M].北京:中国电力出版社,2003:202～231
[10]任震,黄群古,等.基于多频带小波变换的电力系统谐波分析新方法[J].中国电机工程学报,2000,20(12):38～46
[11]武小红.基于小波变换的农村电网谐波检测[J].中国农村水利水电,2007,2(2):99～101