二维EMD的纹理分析及图像瞬时频率估计

1、第17卷第10期2005年10月计算机辅助设计与图形学学报JOURNALOFCOMPUTER2AIDEDDESIGN&COMPUTERGRAPHICSVol117,No110Oct1,2005二维EMD的纹理分析及图像瞬时频率估计沈滨崔峰彭思龙(中国科学院自动化研究所国家专用集成电路设计工程研究中心北京100080)(bin1shenia1ac1cn)摘要将一维经验模式分解(EMD)方法和Hilbert2Huang变换(),的二维EMD算法,AM2FM形式的瞬时频率1关键词;变换;自适应;AM2FM表示;直接法BidimensionalEMDforTextureAnalysisandE

2、stimationoftheInstantaneousFrequenciesofanImageShenBinCuiFengPengSilong(NationalASICDesignEngineeringCenter,InstituteofAutomation,ChineseAcademyofSciences,Beijing100080)AbstractAnewalgorithmbasedonbi2dimensionalempiricalmodedecomposition(BEMD)isdevelopedtoextractfeaturesoftextureatmultiple2scalespat

3、ialfrequencies1Hilbert2Huangtransform(HHT)anddi2rectmethodareappliedtoestimatetheinstantaneousfrequenciesofthemulti2componentAM2FMimagere2spectively1Allalgorithmsaredemonstratedintheexperimentwitheithersyntheticornaturalimages1KeywordsBEMD;textureanalysis;Hilbert2Huangtransform(HHT);adaptive;AM2FMre

4、presenta2tion;directmethod1引言纹理分析主要有4个研究方向1:纹理分割,纹理分类,纹理合成和根据纹理的三维重建1纹理建模主要有以下方法:统计模型,几何模型,随机场模型,分形模型和信号处理模型1纹理图像的信号处理模型包括许多方法:Gabor滤波器2,小波变换,AM2FM表示3等1图像的AM2FM模型将非平稳纹理分析方法已经应用于许多领域,如遥感图像分析,文档图像分析,医学图像分析,自动检测等1针对非线性、非平稳性信号的处理,Huang等4提出了一种自适应的、非监督的、数据驱动的多分辨分解方法经验模式分解(EmpiricalModeDe2composition,EMD),

5、该方法将信号分解为局部窄带的各分量内蕴模函数(IntrinsicModeFunction,IMF)之和,并对分解后得到的各分量IMF进行Hilbert变换,获得分量的瞬时频率和振幅,即Hilbert谱1EMD方法是一种自适应的信号处理方法,已用于一维(地震和生物医学)信号处理1Nunes的、多分量的图像在空域2频域上表示为局部窄带的各分量之和,从各分量中提取出图像的局部信息,瞬时频率是其中重要的信息1收稿日期:2004-07-05;修回日期:2005-01-13基金项目:国家自然科学基金(60272042,10171007)等5将一维EMD思想应用于二维图像处理中12346计算机辅助设计与图形

6、学学报2005年美国NASA宇航中心将这种形式的Hilbert变换称为Hilbert2Huang变换(Hilbert2HuangTransform,HHT)1由于EMD方法具有很好的局部性和自适应性,我们将其应用于纹理的多尺度分解1本文将一维EMD方法和HHT直接应用于对纹理图像的处理,提出了一种新的分析方法BEMD(BidimensionalEMD),并提出了两种估计图像的瞬时频率的方法1抽取,提取出图像的二维IMF1我们针对二维图像提出了一种新的二维EMD方法BEMD,将图像的每一行和每一列看成是一维信号,用一维EMD方法分别对它们进行分解,从而提取出图像的水平和垂直方向的纹理结构1因此本

7、文的BEMD水平和垂直方向的纹理图像11=,I;2基于EMD的纹理分析211一维EMDEMD,值,从而得到residue的行均值矩阵IStep31确定Irowmean;rowmean的每一列的极大值点和极小值点,自适应地分解为频率由高到低的、局部窄带的各分量,即IMF;该分解算法称为筛分过程1再对分解后得到的各个IMF进行Hilbert变换,得到各分量的Hilbert谱1IMF具有两个特点:(1)函数的过零点数目与用三次样条插值得到每一列的上下包络及两包络的均值,从而得到residue的均值矩阵Imean;mean;meanStep41residue=residue-IStep51重复Step

8、1Step4,直至I=01上述筛分算法可以提取出图像I的一个二维IMF,继续筛分可以得到更多的IMF1极值点数目相等或至多相差1;(2)函数关于局部零均值对称1212HHT由于二维图像每一行(列)的数据与一维信号相比要短很多,更容易受到BEMD方法边界效应的影响,因此需要进一步讨论BEMD方法边界效应的抑制1214BEMD方法边界效应的抑制用EMD方法将数据分解成IMF和趋势项后,对各个IMF进行Hilbert变换,得到各自的瞬时频率和瞬时振幅,瞬时振幅在时间-频率平面上的分布就是Hilbert谱1与Fourier谱相比,Hilbert谱具有更高的分辨率,避免了Fourier谱中许多多余的频率

9、成分1将一信号通过EMD方法和Hilbert变换后,可以表示成AM2FM形式,信号的每一个IMF就是信号的一个AM2FM分量;而且,对信号每进行一次筛分过程,提取出的IMF都是当前信号中局部最高频率的振荡1213二维EMD方法在一维EMD方法的筛分过程中,构成上下包络的三次样条函数在数据序列的两端会出现发散现象,使边界产生较大误差,而且,这种误差随着筛分过程的不断进行而向内传播,从而污染整个数据序列1另外,在进行Hilbert变换时,信号的两端也会出现边界效应1对于较长的一维数据序列,可以根据极值点的情况不断舍弃两端的数据,以确保分解的精度,但对于二维图像来说,这种方法完全不可行1Huang等

10、为了抑制EMD方法的边界效应,提本文的主要目的就是将对一维信号进行处理的EMD方法和Hilbert变换推广到对二维图像的分析出用特征波(CharacteristicWave,CW)对原始信号进行延拓的方法,而且将该方法在美国申请了专利,但没有给出确定合适的CW的具体方法1针对EMD方法边界效应的抑制问题,邓拥军等6提出了基于神经网络技术对数据序列进行延拓的方法,张郁山等7提出了基于AR模型的Levinson2Durbin算法对数据序列进行线性预测的中去,将图像表示成AM2FM形式1作为一维IMF的推广,图像的每一个二维IMF都是图像的一个零均值的二维AM2FM分量1而且,对图像每进行一次筛分过

11、程,提取出的二维IMF都具有当前图像中局部最高的空间振荡频率,也就是当前图像的纹理特征1Nunes等5利用数学形态学中的形态重建算子方法1本文采用最大熵谱估计法,即Burg方法来进行边界延拓,这是由于Burg方法尤其适用于对短数和径向基函数提出了一个二维筛分过程以进行纹理10期沈滨等:二维EMD的纹理分析及图像瞬时频率估计2347据情况的预测1Burg方法采用最小化正向和反向两个预测误差的和,利用Levinson递推关系来决定预测滤波器的参数1本文方法如下:在用EMD方法对行(列)数据进行筛分前,首先用Burg方法对数据两端分别进行延拓,直至数据两端各获得一个极大值点和一个极小值点,对延拓数据

12、得到的两个极大(小)值点和数据自身的极大(小)值点进行三次样条插值,得到上(下)包络1由此来抑制EMD方法的边界效应,尤其i:R2R1称ai(x)是第i个分量的幅度调制函数,i(x)是该分量的频率调制函数1t(1jti(x)对于具有K个分量的图像,AM2FM表示为KRet(x)=i=1(x)(),给定一个初ii,:幅度xai),i(x)=是抑制EMD方法的边界效应对低频响,这是由于EMD响较小,1()yx1i(x)=yi式(1)通过对一个单分量图像si(x)进行二维Hilbert变换得到,ti(x)=si(x)+jHsi(x),()x和垂直频率5xi3多分量2图像的瞬时频率估计信号的AM2FM

13、表示可以有效地刻画非平稳信号,估计其瞬时频率1由于多分量信号在某些时刻其瞬时频率不惟一,需要将其分解为瞬时频率惟一的单分量信号1纹理图像一般都是非平稳的多分量信号1为了能将纹理图像分解为单分量信号之和,Havlicek等3其中,H 表示si(x)在ex=1,0T方向的二维Hilbert变换Hsi(x)=Rsi(x-yex)y(xT)si(x)3,xTexey=0,11T提出用复小波滤波器组对纹理图像进行分解,滤波器组的设计原则是对于图像的每一像素点,滤波器组的每一通道至多只对一个分量响应1显然这种方法需要纹理图像的先验知识,不能满足非平稳信号处理的自适应性,而且滤波器组结构复杂1由于每一个IM

14、F在任意时刻都只有一个瞬时频率4,即IMF是单分量信号,因此用EMD方法对非平稳信号进行分解的过程就是将一个多分量信号分解成单分量信号之和的过程;又由于EMD方法具有良好的局部性和自适应性,因此用它处理非平稳信号具有很高的效率1在对多分量的纹理图像进行AM2FM建模过程中,我们用BEMD方法首先将多分量的纹理图像分解成单分量的二维IMF之和,再对二维IMF进行二维Hilbert变换,得到多分量纹理图像的AM2FM表示1最后,我们提出了两种方法来估计单分量的二维IMF的水平瞬时频率和垂直瞬时频率1多分量AM2FM建模技术将图像表示为局部一致的复值单分量之和,第i个单分量具有如下形式(1)ti(x

15、)=ai(x)expji(x),其中,(x)=x1,x2T,ti:R2C,ai:R20,显然,也可以在ey方向上定义不同的二维Hilbert变换,用不同方向的二维Hilbert变换对si(x)进行处理,一般情况下会得到不同的水平瞬时频率和垂直瞬时频率1为了得到惟一的水平瞬时频率和垂直瞬时频率,我们提出了两个估计非平稳的多分量图像的水平频率和垂直频率的方法1311用一维Hilbert变换估计图像瞬时频率用一维Hilbert变换估计瞬时频率的具体步骤如下:Step11用BEMD方法将非平稳的多分量图像分解成数量很少的平稳的单分量二维IMF1Step21估计每一个二维IMF的水平频率:a1初始化1r

16、ow=1(IMF的行索引);b1data=二维IMF第row行数据;c1求出data上下包络和均值means,newdata=data-means;d1对newdata进行Hilbert变换,求得第row行的水平频率;e1row=row+1;f1重复be,直到二维IMF的最后一行1Step31估计每一个二维IMF的垂直频率:a1初始化1column=1(IMF的列索引);由式(1)可得ai(x)=|ti(x)|,i(x)=2348计算机辅助设计与图形学学报2005年b1data=二维IMF第column列数据;c1求出data上下包络和均值means,newdata=data-means;导,

17、有(t)=-(t)sin(t)-(t)2cos(t)=S-2S(t),=(t)S(t)|=S(t)|1Sd1对newdata进行Hilbert变换,求得第column列的垂直频率;e1column=column+1;f1重复be,直到二维IMF的最后一列1通过以上分析,我们提出了估计信号s(t)=a(t)cos(t)的瞬时频率,步骤如下:Step11s(t),用三次在实际应用中,由于BEMD方法对图像的筛分存在误差,因此对于一个二维IMF,其均值矩阵只能是近似于零矩阵,也就是说一个二维IMF行和每一列不一定是一维14,IMFs(t)的上(t);(t)幅值a(t)=max(|upper(t)|,

18、|(t)|);Step31令S(t)=s(t)a(t);Step41确定信号s(t)瞬时频率,IMF,然1312用直接法估计瞬时频率=(t)|S1-S2(t),当|S(t)S(t)|1对信号进行Hilbert变换的目的是要获得信号的幅度和相位,从而得到信号的瞬时振幅和瞬时频率1但是Hilbert变换对有限长的数字信号存在边界效应,并且在信号的频率发生变化时也会产生波动1为了克服Hilbert变换的上述缺点,我们提出了一个新方法直接法,来获得信号的瞬时振幅和瞬时频率1设s(t)是一个IMF,获得s(t)的幅度a(t)和相位(t)的传统方法是对s(t)进行Hilbert变换,使得s(t)=a(t)

19、cos(t),并将s(t)的瞬时频率定义为其相位关于时间的导数1我们将a(t)定义为s(t)的上包络,这个上包络就是EMD方法中所用到的上包络1在实际应用中,由于EMD方法将信号分解成IMF时会存在误差,使得分解结果不完全关于局部零均值对称,为了避免出现|cos(t)|>1的情况,将a(t)的定义修改如下:记s(t)的上下包络分别为upper(t)和lower(t),规定s(t)的瞬时振幅为a(t)=max(|upper(t)|(t)S|,当|=11从直接法可以看出,利用信号本身可以直接估计出它的瞬时频率,而无需再对信号进行Hilbert变换,它完全是一种数据驱动的方法,这就是直接法这个

20、名称的来历14结果与讨论我们将BEMD方法应用于合成纹理图像和自然纹理图像的纹理抽取,以验证BEMD方法的有效性和实用性1411纹理抽取从BEMD方法的分解过程可知:纹理图像分解后得到的每一个IMF都包含着不同的空间尺度信息1本文方法可以从图像中提取出不同尺度的纹理,即使不同纹理之间的尺度相差较大,也可以方便地提取出这些纹理1本文所选用的实验图像取自Brodatz图像库8中部分图像(200×200灰度图像),合成纹理图像来自文献5(400×400灰度图像)1首先,我们将BEMD方法应用于合成纹理图像的分解(如图1所示)1该图像由三个不同幅度和频率的正弦分量之和构成,其水平频

21、率分别为h1=40,h2=14和h3=012;垂直频率分别为v1=60,v2=17和v3=013;相应的振幅分别为ah1=40,ah2=30,ah3=190,av1=40,av2=50和av3=1901|,|lower(t)|)1当a(t)0时,令S(t)=s(t)a(t),于是得到S(t)=cos(t),对其两边求导,有(t)=-(t)sin(t)=-sin(t)S(2)当sin(t)0,即S(t)1时,由的非负性,可由式(2)得出=-S(t)(t)|-S2(t)sin(t)=|S;当sin(t)=0,即|S(t)|=1时,对式(2)两边再求我们用BEMD方法将上述纹理进行分解后得到的结果如

22、图2所示110期沈滨等:二维EMD的纹理分析及图像瞬时频率估计2349将图1与图2进行对比可以看出,BEMD的分解结果与理想结果非常相似1模式1对应于最高的空间频率,模式2对应于中间频率,残差对应于低频1残差与理想的低频成分之间由于BEMD的边界效应,仍然存在误差,但是应用Burg方法已经显著地降低了误差1为了进一步说明BEMD方法的有效性,我们比较了原始的合成纹理图像及其各频率分量的密度分布图1图3所示为合成纹理图像及其各频率分量在主对角线上的密度分布图1图4所示为合成纹理图像及其用BEMD方法分解后得到的各分量在主对角线上的密度分布图1从中可以看出,对BEMD方法的边界效应的抑制显著地降低

23、了图像分解的误差,尤其是降低了BEMD方法的边界效应对低频分量的影响1图3合成纹理(频率成分)图4合成纹理(BEMD)2350计算机辅助设计与图形学学报2005年图5所示为用BEMD方法分解Brodatz图像库中的纹理图像后得到的结果,从中可以看出纹理结构从细到粗的变化1412EMD方法的边界效应对信号的低频分量影响最大,1a进,分别用Hilbert变换和直接法估计了合成纹理图像各分量的水平频率和垂直频率,结果如图6,7所示1为了更清楚地比较两种方法所取得的结果,我们从两种方法所估计的残差频率分量中各取出相对应的一行和一列数据,对比显示如图8所示1通过构造合成纹理图像,其理想的频率值应图8b,

24、8d表明,在两种方法的边界效应的影响下,对信号低频分量频率的估计更是得出了负频率的结果,完全无法得到对理想值的近似估计1用直接法估计瞬时频率时,首先需要拟合信号的上下包络,而确定信号极值点的准确与否将直接影响信号的上下包络的准确性1信号极值点的确定与信号采样率有直接关系,采样率越大,离散信号的极值点就越接近理想值,所得到的频率估计值也就越近似于理想值1该是:水平频率分别是h1=40,h2=14,h3=012;垂直频率分别是v1=60,v2=17,v3=0131从图6,8b,8d可以看出,Hilbert变换的边界效应十分明显1图6对合成图像用Hilbert变换所得到的瞬时频率估计10期沈滨等:二

25、维EMD的纹理分析及图像瞬时频率估计2351图7对合成图像用直接法所得到的瞬时频率估计图8用直接法和Hilbert变换估计残差瞬时频率结果比较一步的研究15结论我们将一维EMD方法扩展到二维,即BEMD方法,并将其应用于纹理图像的分析中,提取出图像中不同的空间尺度和空间频率信息,并基于BEMD提出了两种估计空间频率的方法1BEMD具有优越的局部分析能力,是一种自适应的、完全非监督的方法1与Nunes等提出的BEMD方法相比,本文提出的抑制EMD边界效应的方法得到的分解结果更为准确,将EMD方法进一步应用于瞬时频率的估计,获得了良好的结果1在实际计算中,BEMD方法运算速度更快,收敛迅速1本文方

26、法也可以应用于纹理压缩等领域1本文BEMD方法对于水平和垂直结构的纹理比较有效,稍加修改,也可以将该方法推广至具有显著方向性的纹理,对于更复杂的纹理情况还需要进21由于EMD还是一种比较新的信号分析方法,有些方面还需要更深的研究和探讨,如对于EMD方法的收敛性还需要进一步研究,标准偏差的定义与选择对分解结果会产生重要影响,因此也需要更深入的工作,但是其对非线性系统和非平稳数据优越的处理能力,一定会在更多领域得到应用1参考文献TuceryanM,JainAK1TextureAnalysisM1In:HandbookofPatternRecognitionandComputerVision12nded1Singa2pore:WorldScientificPublishingCo1,19981207248TurnerMR1TexturediscriminationbyGaborfunctionsJ1BiologicalCybernetics,1986,55(2):71823Hav