基于二维集合经验模式分解的医学MR图像特征提取方法.pdf

一种基于二维集合经验模式分解的医学MR图像特征提取方法：步骤一，进行初始化；步骤二，向原始信号中加入白噪声序列；步骤三，通过筛分提取第i个BIMF分量；步骤四，加入不同的白噪声序列重复步骤三；步骤五，更新余量；步骤六，若极值点个数小于给定阈值或分解达到指定层数时算法停止，否则转到步骤二进行下层分解；步骤七，将BIMF的集成平均值作为最终的分解结果；步骤八，对得到的BIMF分解结果进行Hilbert变换，构建解析信号，得到所有BIMF图像的局部特征信息，同时对第一个BIMF图像进行相位一致性的边缘特征提取；最后综合分析图像的所有特征信息。本发明通过提高MR图像处理的精确度，从而提高诊断正确性。

权利要求书
1.  一种基于二维集合经验模式分解的医学MR图像特征提取方法，其特征在于，对于二维图像信号I(x,y)进行如下处理：
步骤一，初始化r0(x,y)＝I(x,y)，定义求解BIMF分量次数的控制变量为i，并且首次求解BIMF分量的i＝1；
步骤二，向原始信号中加入白噪声序列；
步骤三，通过筛分提取第i个BIMF分量，具体操作如下：
(1)初始化：h0(x,y)＝ri-1(x,y)，定义对应求解第i个BIMF分量时的筛选控制变量为j，并且首次求解第i个BIMF分量的j＝1；
(2)利用极值点选择算法找出hj-1(x,y)中所有的极大值点与极小值点；
(3)使用曲面插值算法分别对(2)中的极大值点和极小值点进行插值，得到上包络面emax(x,y)和下包络面emin(x,y)；
(4)计算上下包络面的均值mj-1(x,y)＝[emax(x,y)]+emin(x,y)]/2；
(5)得到筛选函数hj(x,y)＝hj-1(x,y)-mj-1(x,y)；
(6)判断：当hj(x,y)满足筛分停止条件，则BIMFi(x,y)＝hj(x,y)；否则j＝j+1，转到(2)；
步骤四，加入不同的白噪声序列重复步骤三；
步骤五，更新余量ri(x,y)＝ri-1(x,y)-BIMFi(x,y)；
步骤六，如果ri(x,y)的极值点个数小于给定的阈值或分解达到指定层数时则算法停止；否则，转到步骤二继续进行下一层的分解；
步骤七，将每次得到的对应BIMF分量的集成平均值作为最终的分解结果；
步骤八，对得到的BIMF分解结果进行Hilbert变换，构建解析信号，得到所有BIMF图像的局部特征信息，同时对第一个BIMF图像进行相位一致性的边缘特征提取得到图像的边缘特征信息；
步骤九，综合分析图像的所有特征信息，完成图像特征提取。

2.  根据权利要求1所述的基于二维集合经验模式分解的医学MR图像特征提取方法，其特征在于，所述的步骤三的第(6)步中筛分停止条件为满足BIMF分解结果的基本条件，即通过相邻两次筛选结果的标准差进行判断，标准偏差SD的计算公式为：
SD=Σm=1MΣn=1N[|hi,j-1(m,n)-hi,j(m,n)|2hi,j-12(m,n)];]]>
其中，hi,j-1(m,n)和hi,j(m,n)分别是第i个BIMF分解过程中两个连续处理结束的时间序列，二维EMD的标准偏差SD的阈值η设定在0.1～0.3之间，M*N等于二维图像信号I(x,y)的像素大小。

3.  根据权利要求1所述的基于二维集合经验模式分解的医学MR图像特征提取方法，其特征在于：所述的步骤六中ri(x,y)的极值点个数小于给定的阈值或分解达到指定层数时，其中阈值设为2个，得到3～6个BIMF分量以及1个余量，则算法停止。

4.  根据权利要求1所述的基于二维集合经验模式分解的医学MR图像特征提取方法，其特征在于，所述的步骤八中构建解析信号的具体操作方法为：
a.将一维信号x(t)通过Hilbert变换公式得到
b.构造复信号z(t)=x(t)+jx^(t)=a(t)ejθ(t);]]>
则z(t)称为x(t)的解析信号，其中：
a(t)=x2(t)+x^2(t),θ(t)=arctan(x^(t)x(t));]]>
分别表示信号x(t)的局部振幅和局部相位。

5.  根据权利要求1所述的基于二维集合经验模式分解的医学MR图像特征提取方法，其特征在于：所述的Hilbert变换采用二维欧式空间向量值扩展后的Riesz变换。

6.  根据权利要求5所述的基于二维集合经验模式分解的医学MR图像特征提取方法，其特征在于，所述的Riesz变换求解局部特征的具体过程为：
a.根据Riesz核的空域表达式：
(Rx(x),Ry(x))=(x2π|x|3,y2π|x|3),x=(x,y)∈R2;]]>
上式中x为信号分析时域，即图像域，进行Fourier域的变换后，得到表达式：
u＝(u,v)∈R2，u为信号分析频域；
b.根据上述公式，对于二维输入信号，其单演信号为：
fM(x,y)＝(f(x,y),Rx{f}(x,y),Ry{f}(x,y))＝(f,Rx*f,Ry*f)；
其中，*表示卷积运算；
c.分别通过以下公式得到局部振幅lA，局部相位lp，局部方向lo：
lA=|fM(x,y)|=f2+Rx2{f}+Ry2{f};]]>
lp=atan2(Rx2{f}+Ry2{f},f)p∈[0,π];]]>
lo=arctan(Ry{f}Rx{f})o∈[0,π];]]>
通过局部相位，进一步得到局部频率lf为：
lf=(∂lp∂x)2+(∂lp∂y)2.]]>

7.  根据权利要求6所述的基于二维集合经验模式分解的医学MR图像特征提取方法，其特征在于，所述的Riesz变换获取局部特征信息的具体步骤如下：
步骤1：分解原始图像得到K个BIMF；
步骤2：计算第i个BIMF的Riesz变换得到相应的单演信号mi，其中，1≤i≤K；
步骤3：分别计算出mi对应的局部振幅，局部相位，局部方向和局部频率；
步骤4：重复步骤2及步骤3，直到所有BIMF的局部特征都得到求解。

8.  根据权利要求1所述的基于二维集合经验模式分解的医学MR图像特征提取方法，其特征在于，所述的对第一个BIMF进行相位一致性的边缘特征提取的方法为：
a.对输入图像进行BEEMD分解；
b.取第一个BIMF并计算其相位一致性度量PC；
c.对PC进行非极大值抑制得到边缘图像；
d.对边缘图像进行滞后阈值处理得到最终的二至边缘图。

9.  根据权利要求8所述的基于二维集合经验模式分解的医学MR 图像特征提取方法，其特征在于：所述相位一致性度量PC的计算公式为：PC(x,y)=F(x,y)2+H(x,y)2Σm*nAm*n(x,y)+ϵ;]]>
其中，从余弦和正弦小波的卷积导出小波的响应和分别计算出尺度m*n的偶分量em*n(x)和奇分量om*n(x)，变换[em*n(x,y),om*n(x,y)]=[f(x,y)*Mm*ne,f(x,y)*Mm*no]]]>的结果振幅得到局部能量：Am*n(x,y)=em*n(x,y)2+om*n(x,y)2;]]>对偶分量和奇分量求和得到：F(x,y)=Σm*nem*n(x,y)2,]]>H(x,y)=Σm*nem*n(x,y)2,]]>ε为机器能够显示的最小值。

说明书基于二维集合经验模式分解的医学MR图像特征提取方法
技术领域
本发明属于医学图像处理技术领域，具体涉及一种基于二维集合经验模式分解的医学MR图像特征提取方法。
背景技术
随着MRI技术的发展，产生了大量的医学图像，这些医学图像是医务工作者进行疾病的诊断、疾病的跟踪、手术的流程、术后康复的重要依据材料，因此对这些图像数据的处理成为临床应用的最大难题。目前，这些图像的临床分析主要通过医生对图像的定性评价来完成。由于缺乏图像特征的定量度量，人们视觉感知的差异，不同特征和诊断标准的使用，以及人工观察的主观因素和医生的经验等，导致不同医生诊断结论不同，多位医生之间难以形成统一的、正确的诊断结果。这样，良性或恶性的诊断结果都需要依靠活检进行确认，大量的阴性活检造成了病人不必要的痛苦和费用，使病人的身心受到伤害。最为关键的是如此大量的图像信息若依常规方式逐层解读，诊断工作量大，容易造成医生的疲劳，导致误诊或漏诊率的上升。随着计算机技术和图像处理技术突飞猛进的发展，为这些诊断数据的处理提供了新的方法，计算机辅助诊断技术应运而生。
计算机辅助诊断技术以医学影像技术为基础，结合计算机的处理分析功能，帮助影像医师发现病灶区域以提高临床诊断的效率。其中医学图像的特征提取能够获取图像中独特的有别于其他图像的特征 信息，进而实现医学图像目标的自动分类、匹配、识别等功能，是实现计算机辅助诊断的一个关键步骤。
医学图像数据的多样性和重要性要求运用高效准确的图像处理算法提取图像中的特征信息。有效的特征提取算法，可以准确的分离出人体组织和病灶区域，这样就可以减轻医生的负担，提高医生的诊断效率，给医生的临床诊断提供更多的参考信息。目前，常用的图像特征提取方法主要有梯度算子、拉普拉斯算子等，这些算法通用性较差，对噪声也非常敏感；小波分析在图像处理中应用广泛，该方法具有多尺度、多分辨率的特性，但小波基的选择是关键，采用不同的小波基，分解效果会有所不同；分形理论作为一种非线性的处理方法，也是图像特征提取中比较热门的理论。这些算法大多都是直接在原图像上进行处理，效果受噪声影响较大。而Snake算法、人工神经网络方法、统计方法等，大多针对辅助诊断或模式识别应用，对图像特征区域划分较细，因此算法一般过于复杂、需要人工干预且时间复杂度很大，难以满足临床诊断的需求。
经验模式分解(Empirical Mode Decomposition，简称EMD)方法由美国华裔科学家HUANG于1998年提出，方法扩展了Hilbert变换的应用，突破了传统数据分析方法只能分析线性或平稳数据的局限，开创了一种处理非线性非平稳数据的有效方法，能自适应处理非线性非平稳数据的特性使其迅速在诸多领域得到了广泛的应用。近年来，二维经验模式分解(Bi-Dimensional Empirical Mode Decomposition，BEMD)在图像去噪、图像增强、特征提取图像分割、图像融合和图 像压缩等图像处理各个方面的应用研究也正处于不断深入的过程之中，EMD在图像处理领域的应用已成为学者们研究的热点。
虽然，EMD方法在各领域的应用已经取得了很大的成功，但由于该方法提出不久，其理论研究远未成熟，在应用该方法时存在着端点效应、模态混叠等各类问题。针对这些问题不少学者提出了相应的解决方法。其中，值得一提的是，2009年，HUANG等人在分析了白噪声统计特性的基础上提出了一种新的噪声辅助数据分析方法，集合经验模式分解(Ensemble Empirical Mode Decomposition，简称EEMD)。EEMD方法是对原始EMD方法的巨大改进，这种方法通过给信号加入极小幅度白噪声，利用白噪声频谱均衡分布的特点，用白噪声来均衡噪声的特性，较为理想的解决了模态混叠问题，因而非常适合用于医学MR图像的特征提取。
发明目的
本发明的目的在于针对上述现有技术中的缺陷，提供一种能够解决大多数特征提取算法对MR图像噪声较为敏感，需要人工干预，时间复杂度大，在医学MR图像的特征提取上效果不理想的问题，进而提高MR图像特征提取效率及精度的基于二维集合经验模式分解的医学MR图像特征提取方法。
为了实现上述目的，本发明采用的技术方案为，对于二维图像信号I(x,y)进行如下处理：
步骤一，初始化r0(x,y)＝I(x,y)，定义求解BIMF分量次数的控制变量为i，并且首次求解BIMF分量的i＝1；
步骤二，向原始信号中加入白噪声序列；
步骤三，通过筛分提取第i个BIMF分量，具体操作如下：
(1)初始化：h0(x,y)＝ri-1(x,y)，定义对应求解第i个BIMF分量时的筛选控制变量为j，并且首次求解第i个BIMF分量的j＝1；
(2)利用极值点选择算法找出hj-1(x,y)中所有的极大值点与极小值点；
(3)使用曲面插值算法分别对(2)中的极大值点和极小值点进行插值，得到上包络面emax(x,y)和下包络面emin(x,y)；
(4)计算上下包络面的均值mj-1(x,y)＝[emax(x,y)]+emin(x,y)]/2；
(5)得到筛选函数hj(x,y)＝hj-1(x,y)-mj-1(x,y)；
(6)判断：当hj(x,y)满足筛分停止条件，则BIMFi(x,y)＝hj(x,y)；否则j＝j+1，转到(2)；
步骤四，加入不同的白噪声序列重复步骤三；
步骤五，更新余量ri(x,y)＝ri-1(x,y)-BIMFi(x,y)；
步骤六，如果ri(x,y)的极值点个数小于给定的阈值或分解达到指定层数时则算法停止；否则，转到步骤二继续进行下一层的分解；
步骤七，将每次得到的对应BIMF分量的集成平均值作为最终的分解结果；
步骤八，对得到的BIMF分解结果进行Hilbert变换，构建解析信号，得到所有BIMF图像的局部特征信息，同时对第一个BIMF图像进行相位一致性的边缘特征提取得到图像的边缘特征信息；
步骤九，综合分析图像的所有特征信息，完成图像特征提取。
所述的步骤三的第(6)步中筛分停止条件为满足BIMF分解结果的基本条件，即通过相邻两次筛选结果的标准差进行判断，标准偏差SD的计算公式为：
SD=Σm=1MΣn=1N[|hi,j-1(m,n)-hi,j(m,n)|2hi,j-12(m,n)];]]>
其中，hi,j-1(m,n)和hi,j(m,n)分别是第i个BIMF分解过程中两个连续处理结束的时间序列，二维EMD的标准偏差SD的阈值η设定在0.1～0.3之间，M*N等于二维图像信号I(x,y)的像素大小。
所述的步骤六中ri(x,y)的极值点个数小于给定的阈值或分解达到指定层数时，其中阈值设为2个，得到3～6个BIMF分量以及1个余量，则算法停止。
所述的步骤八中构建解析信号的具体操作方法为：
a.将一维信号x(t)通过Hilbert变换公式得到
b.构造复信号z(t)=x(t)+jx^(t)=a(t)ejθ(t);]]>
则z(t)称为x(t)的解析信号，其中：
a(t)=x2(t)+x^2(t),θ(t)=arctan(x^(t)x(t));]]>
分别表示信号x(t)的局部振幅和局部相位。
所述的Hilbert变换采用二维欧式空间向量值扩展后的Riesz变换。
所述的Riesz变换求解局部特征的具体过程为：
a.根据Riesz核的空域表达式：
(Rx(x),Ry(x))=(x2π|x|3,y2π|x|3),x=(x,y)∈R2;]]>
上式中x为信号分析时域，即图像域，进行Fourier域的变换后，得到表达式：
(Hu(u),Rv(u))=(-iu|u|,-iv|u|),u=(u,v)∈R2,]]>u为信号分析频域；
b.根据上述公式，对于二维输入信号，其单演信号为：
fM(x,y)＝(f(x,y),Rx{f}(x,y),Ry{f}(x,y))＝(f,Rx*f,Ry*f)；
其中，*表示卷积运算；
c.分别通过以下公式得到局部振幅lA，局部相位lp，局部方向lo：
lA=|fM(x,y)|=f2+Rx2{f}+Ry2{f};]]>
lp=atan2(Rx2{f}+Ry2{f},f),p∈[0,π);]]>
lo=arctan(Ry{f}Rx{f}),o∈[0,π);]]>
通过局部相位，进一步得到局部频率lf为：
lf=(∂lp∂z)2+(∂lp∂y)2.]]>
所述的Riesz变换获取局部特征信息的具体步骤如下：
步骤1：分解原始图像得到K个BIMF；
步骤2：计算第i个BIMF的Riesz变换得到相应的单演信号mi，其中，1≤i≤K；
步骤3：分别计算出mi对应的局部振幅，局部相位，局部方向和局部频率；
步骤4：重复步骤2及步骤3，直到所有BIMF的局部特征都得到求解。
所述的对第一个BIMF进行相位一致性的边缘特征提取的方法为：
a.对输入图像进行BEEMD分解；
b.取第一个BIMF并计算其相位一致性度量PC；
c.对PC进行非极大值抑制得到边缘图像；
d.对边缘图像进行滞后阈值处理得到最终的二至边缘图。
所述相位一致性度量PC的计算公式为：PC(c,y)=F(x,y)2+H(x,y)2Σm*nAm*n(x,y)+ϵ;]]>
其中，从余弦和正弦小波的卷积导出小波的响应和分别计算出尺度m*n的偶分量em*n(x)和奇分量om*n(x)，变换[em*n(x,y),om*n(x,y)]=[f(x,y)*Mm*ne,f(x,y)*Mm*no]]]>的结果振幅得到局部能量：Am*n(x,y)=em*n(x,y)2+om*n(x,y)2;]]>对偶分量和奇分量求和得到：F(x,y)=Σm*nem*n(x,y)2,H(x,y)=Σm*nom*n(x,y)2;]]>ε为机器能够显示的最小值。
与现有技术相比，本发明具有如下有益效果：本发明针对目前临床上人工逐层解读大量MR图像信息不但诊断工作量大，而且容易因为医生的疲劳而导致误诊或漏诊率上升的状况，研究基于二维集合经验模式分解BEEMD的医学MR图像特征提取问题，首先通过EMD方法提取复杂信号在每个时间中的局部震荡模式，将复杂信号分解为 一系列有限个从高频到低频排列的固有模态函数之和，然后对每个IMF做Hilbert变换，计算每个IMF的能量，即瞬时频率和振幅，从而形成时间-频率-振幅谱，即Hilbert谱，提取信号的瞬时频率特征。本发明将集合经验模式分解算法由一维推广到更加复杂的二维空间得到二维集合经验模式分解算法，以此来抑制BEMD算法中出现的模态混叠现象，从而获取到图像的局部幅值、局部相位、局部方向和局部频率；在BEEMD分解的基础上再结合相位一致性算法得到MR图像的边缘特征，为下一步特征分类和识别提供依据。结合临床应用背景，本发明重点克服了现有的大多数特征提取算法对MR图像噪声较为敏感，需要人工干预且时间复杂度大，在医学MR图像的特征提取上效果不理想的问题，同时提高了MR图像特征提取的效率及精度。通过对lena图像以及临床MR乳腺和大脑图像进行分解和特征提取，提取结果和BEMD的提取结果进行对比，对比结果说明本发明在性能能上明显优于BEMD，能客观地提取和分析MR图像，为MR图像的处理提供新的理论与方法，使医疗专家摆脱繁重的人工观察和诊断；同时能提供更精确的辅助诊断数据，降低医生诊断的主观性，从而提高诊断正确性和客观性。
进一步的，本发明在筛分停止判断中二维标准偏差SD的阈值η设定在0.1～0.3之间，由于二维标准偏差SD目前还没有明确标准，很大程度上要靠经验值，若η取值过小，停止条件则过于严格，会导致分解的迭代次数过多、计算量大，图像过度分解结果不够理想；若η取值偏大，停止条件宽松，造成分解结果尺度特征分离不明显。 实验中η我们取0.2，一方面它的大小合适，可以作为一个较准确的筛选值；另一方面这个值不至于过小，可以降低迭代次数和计算量。
附图说明
图1本发明图像特征提取方法的整体流程图；
图2合成信号的EMD分解及其Hilbert谱图：其中，图2(a)为由频率为20Hz和100Hz的余弦信号叠加而成的信号的EMD分解图；图2(b)为频率为20Hz和100Hz的余弦信号叠加而成的信号的EMD分解图对应的Hilbert谱图；
图3(a)为待合成的两个模拟信号x1与x2的波形图；
图3(b)为x1与x2的EMD算法分解图；
图3(c)为x1与x2的EEMD算法分解图；
图4本发明BEEMD的算法流程图；
图5本发明结合BEEMD和相位一致性算法的图像边缘特征提取流程图；
图6为lena图像BEEMD的分解图：图6(a)为lena图像原图；图6(b)为BEEMD的BIMF1图；图6(c)为BEEMD的BIMF2图；图6(d)为BEEMD的BIMF3图；图6(e)为BEEMD的BIMF4图；图6(f)为BEEMD的BIMF5图；图6(g)为余量图；图6(h)为BIMF1～5的重建图；
图7为lena图像的BEMD的分解图：图7(a)为lena图像原图；图7(b)为BEMD的BIMF1图；图7(c)为BEMD的BIMF2图；图7(d)为BEMD的BIMF3图；图7(e)为BEMD的BIMF4图； 图7(f)为BEMD的BIMF5图；图7(g)为余量图；图7(h)为BIMF1～5的重建图；
图8为乳腺MR图像的原图；
图9为乳腺MR图像的BEEMD分解图：其中，图9(a)为MR图像原图；图9(b)为BEEMD的BIMF1图；图9(c)为BEEMD的BIMF2图；图9(d)为BEEMD的BIMF3图；图9(e)为BEEMD的BIMF4图；图9(f)为BEEMD的BIMF5图；图9(g)为余量图；图9(h)为BIMF1～5的重建图；
图10为乳腺MR图像的BEMD分解图：其中，图10(a)为MR图像原图；图10(b)为BEMD的BIMF1图；图10(c)为BEMD的BIMF2图；图10(d)为BEMD的BIMF3图；图10(e)为BEMD的BIMF4图；图10(f)为BEMD的BIMF5图；图10(g)为余量图；图10(h)为BIMF1～5的重建图；
图11为大脑MR图像的原图；
图12为大脑MR图像的局部特征提取BEMD与BEEMD对比图，其中，第一行为基于BEMD的提取结果，第二行为基于BEEMD的提取结果，图12(a)为局部幅值对比图；图12(b)为局部相位对比图；图12(c)局部方向对比图；图12(d)为局部频率对比图；
图13为乳腺MR图像的局部特征提取BEMD与BEEMD对比图，其中，第一行为基于BEMD的提取结果，第二行为基于BEEMD的提取结果，图13(a)为局部幅值对比图；图13(b)为局部相位对比图；图13(c)局部方向对比图；图13(d)为局部频率对比图；
图14为乳腺MR图像的边缘特征提取结果图：图14(a)为Prewitt算子检测结果图；图14(b)为Roberts算子检测结果图；图14(c)为Sobel算子检测结果图；图14(d)高斯拉普拉斯算子检测结果图；图14(e)为Canny算子检测结果图；图14(f)本发明检测结果图；
图15为大脑MR图像的边缘特征提取结果图：图15(a)为Prewitt算子检测结果图；图15(b)为Roberts算子检测结果图；图15(c)为Sobel算子检测结果图；图15(d)高斯拉普拉斯算子检测结果图；图15(e)为Canny算子检测结果图；图15(f)本发明检测结果图。
具体实施方式
为了使本领域技术人员更好地理解本申请中的技术方案，下面结合附图对本申请技术方案做进一步的详细说明，所描述的实施例仅仅是本申请的部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其它实施例，都应当属于本发明保护的范围。
参见图1，本发明基于二维集合经验模式分解的医学MR图像特征提取方法，对于二维图像信号I(x,y)进行如下处理：
步骤一，初始化r0(x,y)＝I(x,y)，定义求解BIMF分量次数的控制变量为i，并且首次求解BIMF分量的i＝1，对BIMF分量进行N求解；
步骤二，向原始信号中加入白噪声序列；
步骤三，通过筛分提取第i个BIMF分量，具体操作如下：
(1)初始化：h0(x,y)＝ri-1(x,y)，定义对应求解第i个BIMF分量 时的筛选控制变量为j，并且首次求解第i个BIMF分量的j＝1；
(2)利用极值点选择算法找出hj-1(x,y)中所有的极大值点与极小值点；
(3)使用曲面插值算法分别对(2)中的极大值点和极小值点进行插值，得到上包络面emax(x,y)和下包络面emin(x,y)；
(4)计算上下包络面的均值mj-1(x,y)＝[emax(x,y)]+emin(x,y)]/2；
(5)得到筛选函数hj(x,y)＝hj-1(x,y)-mj-1(x,y)；
(6)判断：当hj(x,y)满足筛分停止条件，则BIMFi(x,y)＝hj(x,y)；否则j＝j+1，转到(2)；筛分停止条件为满足BIMF分解结果的基本条件，即通过相邻两次筛选结果的标准差进行判断，标准偏差SD的计算公式为：
SD=Σm=1MΣn=1N[|hi,j-1(m,n)-hi,j(m,n)|2hi,j-12(m,n)];]]>
其中，hi,j-1(m,n)和hi,j(m,n)分别是第i个BIMF分解过程中两个连续处理结束的时间序列，二维EMD的标准偏差SD的阈值η设定在0.1～0.3之间，实施例η取0.2，M*N等于二维图像信号的像素大小；
步骤四，加入不同的白噪声序列重复步骤三；
步骤五，更新余量ri(x,y)＝ri-1(x,y)-BIMFi(x,y)；
步骤六，如果ri(x,y)的极值点个数小于给定的阈值或分解达到指定层数时，得到3～6个BIMF分量以及1个余量，则算法停止；极值点个数给定的阈值为2个；
步骤七，将每次得到的对应BIMF的集成平均值作为最终的分解 结果；
步骤八，对得到的BIMF分解结果进行Riesz变换，构建解析信号，得到局部振幅、相位、频率等所有BIMF图像的局部特征信息，同时对第一个BIMF图像进行相位一致性的边缘特征提取得到图像的边缘特征信息；
步骤九，综合分析图像的所有特征信息，完成图像特征提取。
1.EMD算法基本原理为：
经验模式分解(EMD)是一种直观的，直接的、先验的、自适应的分解方法，方法根据信号的特征时间尺度将信号分解为一组固有模态函数(Intrinsic Mode Function,IMF)，或称基本模式分量的和。每个IMF分量必须满足两个条件：
1)在整个信号序列中，极值点(即大值点与极小值点)的个数Ne和过零点数目Nz必须相等或最多相差不超过一个，即
(Nz-1)≤Ne≤(Nz+1)                  (1)
2)任意时间点ti上，信号序列局部最大值所确定的上包络线Smax(t)与局部最
小值所确定的下包络线Smin(t)关于时间轴局部对称，即均值为零：
(Smax(t)+Smin(t))/2＝0                (2)
IMF表征了信号内在的波动模式。由上述定义可知，在IMF每一个波动周期(两零点之间)只有一个单纯的波动模式，没有其它波动叠加上去，因此它可以作为分解信号的基本量，也就是说任意信号可以被分解为若干个IMF。EMD算法就是这样一个“筛选”(的过程， 具体由以下几步组成：
步骤1：找出信号x(t)的所有局部极大值和局部极小值；
步骤2：利用三次样条插值方法分别计算极小值插值和极大值插值并得到对应信号包络emin(t)和emax(t)；
步骤3：计算局部均值m(t)＝(emin(t)+emax(t))/2。
步骤4：将原始输入信号减去局部均值得到振荡信号：
h(t)＝x(t)-m(t)。
步骤5：当h(t)满足IMF的条件时，h(t)就成为一个IMF；否则，用h(t)替换步骤1中的x(t)并从步骤1开始重复上述过程。
步骤6：令c1＝h(t)，则c1为第一个IMF，对应的余量r1＝x(t)-c1。
当迭代终止条件不满足时，用r1替换原信号并重复上述所有步骤并得到第二个IMF分量c2及余量r2，如此反复，最后得到第n个IMF、cn和rn。当满足如下预先约定的条件时，停止筛分过程：cn或rn小于预先设定的值，或者rn变成一个单调函数。
经过上面的分解过程，x(t)最终被分解为n个IMF分量和一个余量rn。原始信号x(t)可以表示为对分解后得到的各IMF分量分别求Hilbert变换，得到对应的瞬时频率和瞬时振幅，进而得到Hilbert谱和Hilbert边缘谱。图2(a)为由频率为20Hz和100Hz的余弦信号叠加而成的信号x1(t)＝cos(2π·20·t)+cos(2π·100·t)的EMD分解结果，图2(b)为其对应的Hilbert谱。
2.EEMD算法实现
原始EMD算法的主要缺点之一就是模态混叠的频繁出现。模态 混叠有两种表现形式：(1)单个IMF中包含不同的时间尺度成份；(2)同一尺度出现在不同的IMF中。模态混叠的出现会导致产生严重错假的时频分布也使IMF失去物理意义。集合经验模式分解(EEMD)主要解决模态混叠的问题，其分解步骤如下：
步骤1：向原始信号中加入白噪声序列；
步骤2：将添加了白噪声的信号通过EMD算法分解为一系列的IMF；
步骤3：重复步骤1、步骤2，但每次加入不同的白噪声序列；
步骤4：将每次得到的对应IMF的集成平均值作为最终的分解结果。
EEMD方法需要两个重要的参数：所添加白噪声的幅值和集成平均次数。WU和HUANG通过实验总结得出只要添加的噪声经过频率调制而且集成平均次数足够大，那么噪声幅值和集成平均次数的增加几乎就不会改变分解结果，并建议大多数情况下添加噪声的幅值应该满足其大约为0.2个数据信号的标准差。EEMD方法是对原始EMD方法的巨大改进，这种方法通过给信号加入极小幅度白噪声，利用白噪声频谱均衡分布的特点，用白噪声来均衡噪声的特性，较为理想的去除了模态混叠。
原始经验模式分解与EEMD分解的结果对比如图3(a)～(c)所示。其中，图3(a)为待合成的两个模拟信号x1与x2，图3(b)和图3(c)中的输入信号为x1与x2的累加合成信号。图3(b)为EMD算法分解的结果，图3(c)为EEMD算法分解的结果图。从图 3(b)中可以看出，EMD算法分解的IMF1和IMF2中都出现了两个模态混叠的现象，基本上没有区分出信号x1和x2，而图3(c)中EEMD算法基本上很清晰地将x1与x2分解出来，证明了EEMD在抗模态混叠方面的有效性。
3.BEEMD算法实现
1).BEMD算法实现：
将一维EMD算法推广到更加复杂的二维空间便得到了二维经验模式分解算法BEMD。与一维经验模式分解算法相似，其分解步骤也是通过筛分过程来完成，每次筛分完成后得到一个对应的模态函数。分解完成后得到数量有限的几个二维固有模态函数(Bi-Dimensional Intrinsic Mode Function，简称BIMF)及一个余量。对二维图像信号I(x,y)的分解算法如下：
步骤1：初始化r0(x,y)＝I(x,y),i＝1。
步骤2：筛分，提取第i个BIMF分量：
(1)初始化：h0(x,y)＝ri-1(x,y),j＝1。
(2)利用极值点选择算法找出hj-1(x,y)中所有的极大值点与极小值点。
(3)使用曲面插值算法分别对(2)中的极大值点和极小值点进行插值，得到上包络面emax(x,y)和下包络面emin(x,y)。
(4)计算上下包络面的均值mj-1(x,y)＝[emax(x,y)]+emin(x,y)]/2。
(5)得到筛选函数hj(x,y)＝hj-1(x,y)-mj-1(x,y)。
(6)判断，当hj(x,y)满足筛分停止条件，则BIMFi(x,y)＝hj(x,y)； 否则j＝j+1，转到(2)。
步骤3：更新余量ri(x,y)＝ri-1(x,y)-BIMFi(x,y)。
步骤4：如果ri(x,y)的极值点个数小于给定的阈值或分解达到指定层数，则停止算法；否则，转到步骤2继续进行下一层的分解。
原始图像I(x,y)在经过BEMD算法处理后，可以得到一系列BIMF分量和一个余量图像而且它们之间满足：
I(x,y)=Σi=1nBIMFi(x,y)+rn(x,y)---(3)]]>
2)BEEMD实现过程:
BEMD算法也存在着一维EMD中存在的模态混叠问题，将EEMD算法中的噪声辅助方法应用于BEMD中，便得到二维集合经验模式分解算法BEEMD。与EEMD类似，BEEMD通过N次迭代实现，每次迭代的输入信号中加入白噪声，对加入噪声的信号进行BEMD算法分解，N次迭代分解得到的各BIMF的平均值为最后的分解结果，BEEMD的算法流程图如图4所示。
4.局部特征计算：
在一维信号处理中，解析信号(Analytic Signal)是一种重要的复数表示方法，并广泛地应用于信息编码(振幅调制、相位调制)、语音识别、基于雷达的目标检测、地震数据的处理等领域。一维信号的解析函数通过Hilbert变换公式得到，如(4)式所示
x^(t)=1π∫-∞∞x(t)t-t′dt′---(4)]]>
其中为x(t)的Hilbert变换。
构造复信号z(t)：
z(t)=x(t)+jx^(t)=a(t)ejθ(t)---(5)]]>
则z(t)称为x(t)的解析信号，其中
a(t)=x2(t)+x^2(t),θ(t)=arctan(x^(t)x(t))---(6)]]>
分别表示信号x(t)的局部振幅(Local Amplitude)和局部相位(Local Phase)，局部振幅和局部相位满足不变性(Invariance)和同变性(Equivariance)。信号局部能量的变化不会引起局部相位的变化，但信号局部结构的变化将引起局部相位的变化；局部结构的变化不会引起局部振幅的变化，局部振幅代表着局部能量。在未混合不同尺度、不同局部相位信号的信号中，信号能量和结构是独立的信息。为了保持不变性、同变性，信号必须通过带通滤波以去除部分信号。
在一维经验模式分解中，用一维的Hilbert变换来构造输入信号对应的解析信号进而可以得到输入信号的频谱信息。为了进一步分析二维固有模态函数的局部特征必须先得到对应的二维解析信号。要求解析信号必须满足以下特征：
1)变换后的信号与输入信号必须正交。
2)忽略直流分量，解析信号的能量为输入信号能量的两倍。
为了满足上述条件，Michael Felsberg等人基于Riesz变换和向量场理论提出了一种新的二维解析信号并将其命名为单演信号(Monogenic Signal)。单演信号基于一阶Riesz变换，是Hilbert变换在二维欧式空间的向量值扩展，Riesz核的空域表达式如下：
(Rx(x),Ry(x))=(x2π|x|3,y2π|x|3),x=(x,y)∈R2---(7)]]>
其在Fourier域的变换函数为：
(Hu(u),Rv(u))=(-iu|u|,-iv|u|),u=(u,v)∈R2---(8)]]>
对于二维输入信号f(x,y)，其单演信号定义为：
fM(x,y)＝(f(x,y),Rx{f}(x,y),Ry{f}(x,y))＝(f,Rx*f,Ry*f)       (9)
其中，*表示卷积运算。于是，局部振幅lA，局部相位lp，局部方向lo分别通过以下公式得到：
lA=|fM(x,y)|=f2+Rx2{f}+Ry2{f}---(10)]]>
lp=atan2(Rx2{f}+Ry2{f},f),p∈[0,π)---(11)]]>
lo=arctan(Ry{f}Rx{f}),o∈[0,π)---(12)]]>
通过局部相位，可进一步得到局部频率lf为
lf=(∂lp∂z)2+(∂lp∂y)2---(13)]]>
5.基于BEEMD的局部特征提取：
从经验模式分解的过程可以看出，BEEMD算法将图像分解为数个BIMF分量，第一个分量BIMF1包含几乎所有的图像细节信息，第二个分量BIMF2只含有少量的细节信息，以此类推，余量几乎不含任何细节信息。
经过BEEMD分解后得到的BIMF由于其只与输入图像有关本身即可作为图像的一种特征，结合Riesz变换，可以进一步得到局部特征信息，具体步骤如下：
步骤1：使用BEEMD算法分解原始图像image得到K个BIMF。
步骤2：计算第i(1≤i≤K)个BIMF的Riesz变换得到相应的单演信号mi。
步骤3：根据公式(10)，(11)，(12)，(13)分别计算出mi对应的局部振幅，局部相位，局部方向和局部频率。
步骤4：重复步骤2、步骤3直到所有BIMF的局部特征都得到求解。
通过对局部振幅，局部相位等低层局部特征的分析和使用，可以得到更高层次的特征。
6.边缘特征提取：
1)相位一致性
基于梯度的边缘检测方法(Sobel算子，Canny算子等)对图像亮度的变化都很敏感，亮度梯度的影响因素包括光照强度、模糊度、放大倍数等，当将图像放大到两倍时，任何基于梯度的边缘检测算法都必须适当地修改阈值。但是，通常图像的亮度和放大倍数事先并不知晓，因此，有效边缘对应的梯度值通常根据经验来设定。1986年Morrone等人提出了一种新的特征感知模型—局部能量模型，该模型假设图像的特征位于最大相位一致性的像素点，与人类视觉系统关系密切。考虑一维信号f(x)，其Fourier级数展开表达式为：
f(x)=Σnancos(nωx+φn),n≥0,an>0---(14)]]>
其中，an，φn分别为Fourier分量的振幅和初始相位，ω为常数，一般为2π，则相位一致性函数定义为：
PC(x)=maxθ∈[0,22π]Σnancos(nωx+φn-θ)Σnan---(15)]]>
θ为某点处所有Fourier分量的加权平均相位角。显然，PC的取值范围在0与1之间而且是一个无量纲的归一化量。该式表示的是由不同变量的权值与所有分量总振幅的比值，除了难以实现外，其还对噪声敏感，不易于定位。Venkatesh和Owens证明了相位一致性与局部能量成正比，因此，可以用寻找局部能量最大值的方法来替代相位一致性函数的直接计算，而且局部能量可以弥补在有噪声时相位检测的灵敏度。因此，Kovesi提出了一个基于小波的度量以提高噪声下的性能，一维信号f(x)在尺度对于一组小波的响应可以从余弦和正弦小波的卷积导出，分别记为和分别计算出尺度的偶分量en(x)和奇分量on(x)：
[en(x),on(x)]=[f(x)*Mne,f(x)*Mno]---(16)]]>
该变换结果振幅即为局部能量：
An(x)=en(x)2+on(x)2---(17)]]>
在每个点上都有一组向量对应各个尺度的滤波器。由于只对大范围频率产生的相位一致性感兴趣，只需把这组小波滤波器设计成使相邻分量重叠。对偶分量和奇分量求和，得到：
F(x)=Σnen(x)2,H(x)=Σnon(x)2---(18)]]>
因而总能量为：
ΣnAn(x)≈Σnen(x)2+on(x)2---(19)]]>
最后，得到的相位一致性度量为：
PC(x)=F(x)2+H(x)2ΣnAn(x)+ϵ---(20)]]>
Kovesi在二维情况下对相位一致性的扩展方法为先对小波滤波器和图像进行卷积计算，接着计算出平均滤波响应和单个滤波响应之间的差值，从而确定相位一致性。本发明使用Kovesi的相位一致性度量计算方法。相位一致性模型是一个特征检测算子，其度量是一个归一化的无量纲量，具有可以检测大范围的特征，对光照变化、图像的对比度具有不变性。
本发明由此延伸至二维相位一致性度量PC的计算公式为：PC(c,y)=F(x,y)2+H(x,y)2Σm*nAm*n(x,y)+ϵ;]]>
其中，从余弦和正弦小波的卷积导出小波的响应和分别计算出尺度m*n的偶分量em*n(x)和奇分量om*n(x)，变换[em*n(x,y),om*n(x,y)]=[f(x,y)*Mm*ne,f(x,y)*Mm*no]]]>的结果振幅得到局部能量：Am*n(x,y)=em*n(x,y)2+om*n(x,y)2;]]>对偶分量和奇分量求和得到：F(x,y)=Σm*nem*n(x,y)2,H(x,y)=Σm*nom*n(x,y)2;]]>ε为机器能够显示的最小值。
2)结合BEEMD和相位一致性的边缘特征提取
由于经验模式分解算法得到的第一个BIMF几乎包含了原图像中的所有细节信息，因此，本发明仅对第一个BIMF执行相位一致性算法可得到图像的边缘信息，其算法流程图如图5所示。
7.算法性能验证：
1)BEEMD算法验证
经验模式分解根据输入信号的不同而自动地选择相应的基函数，能自适应地对数据进行处理，分解后得到从高频到低频的数个固有模态函数，高频部分有效地提取出原信号的细节信息，低频部分反映了原图像的总体能量趋势的分布。
本发明以MATLAB 7.11.0为实验平台，操作系统为Windows 8，在此基础上实现了BEEMD算法，选取筛分终止条件SD＝0.2，分解5层得到5个BIMF和一个余量。对Lena图像和如图8所示临床MR图像，用BEEMD算法的分解结果如图6(a)～(h)、图9(a)～(h)所示，图7(a)～(h)和图10(a)～(h)分别给出了由BEMD算法分解的结果。实验中所采用的MR图像为乳腺癌诊断核磁共振扫描图像，图片由西门子公司的全体核磁共振成像系统生成，图片大小为512×512(像素)。
图6，7，9，10中，(a)为原始图像，(b)～(f)依次为分解所得的BIMF1～BIMF5，图(g)为余量。图(h)为按公式(21)重建后的图像。
Ir=Σi=1nBIMFi---(21)]]>
其中n为所选取的BIMF的个数。
对比图6(a)～(h)和7(a)～(h)、图9(a)～(h)和10(a)～(h)可以看出：不论是BEEMD算法还是BEMD算法，分解结果中从第一个BIMF到第五个BIMF再到余量图像，BIMF包含的原图像的细节信息(即高频部分)越来越少，且原图的细节信息大部分都集中在BIMF1中；由BIMF1～5重构后的图像几乎与原图一样，进一步说明分解所得的BIMF几乎包含了原图中的所有细节信息，而余量图像中只有原图像 的整体趋势信息；同时体现出了算法将输入信号随频率由高到低分解的滤波特性。但是由BEEMD算法分解得到的BIMF比BEMD算法的结果包含了更加丰富、更加完整的细节信息，说明了EEMD原理扩展到二维之后的有效性。
2)图像特征提取算法验证
为了验证本发明提取图像特征的性能，在Matlab平台上实现上述算法，并使用临床大脑MR图像和乳腺MR图像做测试数据，如图11，8所示，运行环境配置如表1所示：
表1运行环境配置
配置项配置值电脑型号Lenovo G470操作系统Windows 8企业版64位处理器Intel Core i5-2460M CPU@2.4GHz内存容量4GBMatlab版本Version 7.11.0.58432-bit
图12(a)～(d)和图13(a)～(d)分别给出了大脑MR图像和乳腺MR图像的局部特征提取结果，图(a)～(d)中各列依次为所提取局部幅值、相位、方向和频率；第一行为用BEMD算法分解提取的特征；第二行为用BEEMD算法分解后提取的结果。从图中可以看出，这些局部特征清晰地表达了原始图像的细节信息。其中，局部振幅体现了图像的强度信息，展示了整体能量的分布；局部相位信息体现出原图像的结构信息，目标的边界信息都能在局部相位图像中看到；局部方向体现了图像在水平垂直两个方向上的变化信息；局部频率体现了图像的区域边界信息。这些特征都从不同角度体现了图像的特征信息，是原 始图像的最低层特征信息，通过这些低层特征进而提取更高层特征是特征提取中常用的方法之一。
对比两幅图上下两行，可以看出用BEEMD算法分解后获取的特征信息更加细腻、清晰，尤其相位特征更能体现出两种算法的差异，用BEEMD算法所提取的相位特征更明显地体现了图像的结构信息，为后续图像特征的分类奠定了一定的基础，表明BEEMD算法在图像特征提取方面具有一定的优势。
另外，图8和图11边缘检测结果分别如图14(a)～(f)、图15(a)～(f)所示。
由于经验模式分解算法在理论方面的不完善，目前，其分解结果的好坏还没有客观的评价方法，只能通过主观来判断。从图14(a)～(f)及图15(a)～(f)中可以看出，在对乳腺MR图像的边缘提取中，视觉上Prewitt，Roberts和Sobel算子得到的边缘图似乎要更好，它们将心脏部分与乳房部分区了开，但原图中存在的很多细节性边缘都没有检测出来；而包括本发明在内的其他三种算法(拉普拉斯算子、Canny算子以及本发明)得到了更加细致的边缘线。在对大脑MR图像的边缘提取中，Prewitt，Roberts和Sobel算子提取到的边缘信息较其他三种算法明显要差。在直观视觉上，本发明提取到的边缘信息不如Canny算子与高斯拉普拉斯算子检测到边缘那么连续、直观，但仔细观察可以发现本发明提取到的边缘包含了更多的边缘细节，证明了本发明的有效性、可行性。
事实上本发明对某些图像边缘检测效果不如Canny算子(如上大 脑图像)应有几方面的影响。原因一，可能是由于经验模式分解算法的不成熟而导致，因为经验模式分解得到的第一个BIMF提取的细节信息的多少决定了本发明边缘检测的结果；原因二，相位一致性模型的建立过程所造成的误差，都将直接影响到边缘的检测。因此，尽可能地完善经验模式分解算法和相位一致性模型是申请人下一步要研究的内容。