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

1、(10)申请公布号 (43)申请公布日 (21)申请号 201410660753.5 (22)申请日 2014.11.18 G06T 7/00(2006.01) G06F 5/00(2006.01) G06K 9/46(2006.01) A61B 5/055(2006.01) (71)申请人 陕西师范大学 地址 710062 陕西省西安市长安区陕西师范 大学 (72)发明人 范虹 刘晓杰 毛玉龙 (74)专利代理机构 西安通大专利代理有限责任 公司 61200 代理人 陆万寿 (54) 发明名称 基于二维集合经验模式分解的医学 MR 图像 特征提取方法 (57) 摘要 一种基于二维集合经验模式

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

3、图像的所有特 征信息。本发明通过提高 MR 图像处理的精确度， 从而提高诊断正确性。 (51)Int.Cl. (19)中华人民共和国国家知识产权局 (12)发明专利申请 权利要求书3页 说明书14页 附图9页 (10)申请公布号 CN 104392444 A (43)申请公布日 2015.03.04 CN 104392444 A 1/3 页 2 1. 一种基于二维集合经验模式分解的医学 MR 图像特征提取方法， 其特征在于， 对于二 维图像信号 I(x,y) 进行如下处理 ： 步骤一， 初始化 r0(x,y) I(x,y)， 定义求解 BIMF 分量次数的控制变量为 i， 并且首次 求解 BI

4、MF 分量的 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

5、(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 分量的集成平均值作为最终的分

6、解结果 ； 步骤八， 对得到的 BIMF 分解结果进行 Hilbert 变换， 构建解析信号， 得到所有 BIMF 图 像的局部特征信息， 同时对第一个 BIMF 图像进行相位一致性的边缘特征提取得到图像的 边缘特征信息 ； 步骤九， 综合分析图像的所有特征信息， 完成图像特征提取。 2.根据权利要求1所述的基于二维集合经验模式分解的医学MR图像特征提取方法， 其 特征在于， 所述的步骤三的第 (6) 步中筛分停止条件为满足 BIMF 分解结果的基本条件， 即 通过相邻两次筛选结果的标准差进行判断， 标准偏差 SD 的计算公式为 ： 其中， hi,j-1(m,n) 和 hi,j(m,n) 分别

7、是第 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图像特征提取方法， 其 特征在于， 所述的步骤八中构建解析信号的具体操作方法为 ： 权

8、 利 要 求 书 CN 104392444 A 2 2/3 页 3 a. 将一维信号 x(t) 通过 Hilbert 变换公式得到 b. 构造复信号 则 z(t) 称为 x(t) 的解析信号， 其中 ： 分别表示信号 x(t) 的局部振幅和局部相位。 5.根据权利要求1所述的基于二维集合经验模式分解的医学MR图像特征提取方法， 其 特征在于 ： 所述的 Hilbert 变换采用二维欧式空间向量值扩展后的 Riesz 变换。 6.根据权利要求5所述的基于二维集合经验模式分解的医学MR图像特征提取方法， 其 特征在于， 所述的 Riesz 变换求解局部特征的具体过程为 ： a. 根据 Riesz

9、核的空域表达式 ： 上式中 x 为信号分析时域， 即图像域， 进行 Fourier 域的变换后， 得到表达式 ： u (u,v) R2， u 为信号分析频域 ； b. 根据上述公式， 对于二维输入信号， 其单演信号为 ： fM(x,y) (f(x,y),Rxf(x,y),Ryf(x,y) (f,Rx*f,Ry*f) ； 其中， * 表示卷积运算 ； c. 分别通过以下公式得到局部振幅 lA， 局部相位 lp， 局部方向 lo： 通过局部相位， 进一步得到局部频率 lf为 ： 7.根据权利要求6所述的基于二维集合经验模式分解的医学MR图像特征提取方法， 其 特征在于， 所述的 Riesz 变换获

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

11、 a. 对输入图像进行 BEEMD 分解 ； b. 取第一个 BIMF 并计算其相位一致性度量 PC ； c. 对 PC 进行非极大值抑制得到边缘图像 ； d. 对边缘图像进行滞后阈值处理得到最终的二至边缘图。 9.根据权利要求8所述的基于二维集合经验模式分解的医学MR图像特征提取方法， 其 特征在于 ： 所述相位一致性度量 PC 的计算公式为 ： 其 中，从 余 弦 和 正 弦 小 波 的 卷 积 导 出 小 波 的 响 应和 分 别 计 算 出 尺 度 m*n 的 偶 分 量 em*n(x) 和 奇 分 量 om*n(x)，变 换 的 结 果 振 幅 得 到 局 部 能 量 ：对 偶 分

12、量 和 奇 分 量 求 和 得 到 ： 为机器能够显示的最小值。 权 利 要 求 书 CN 104392444 A 4 1/14 页 5 基于二维集合经验模式分解的医学 MR 图像特征提取方法 技术领域 0001 本发明属于医学图像处理技术领域， 具体涉及一种基于二维集合经验模式分解的 医学 MR 图像特征提取方法。 背景技术 0002 随着 MRI 技术的发展， 产生了大量的医学图像， 这些医学图像是医务工作者进行 疾病的诊断、 疾病的跟踪、 手术的流程、 术后康复的重要依据材料， 因此对这些图像数据的 处理成为临床应用的最大难题。目前， 这些图像的临床分析主要通过医生对图像的定性评 价来完

13、成。 由于缺乏图像特征的定量度量， 人们视觉感知的差异， 不同特征和诊断标准的使 用， 以及人工观察的主观因素和医生的经验等， 导致不同医生诊断结论不同， 多位医生之间 难以形成统一的、 正确的诊断结果。这样， 良性或恶性的诊断结果都需要依靠活检进行确 认， 大量的阴性活检造成了病人不必要的痛苦和费用， 使病人的身心受到伤害。 最为关键的 是如此大量的图像信息若依常规方式逐层解读， 诊断工作量大， 容易造成医生的疲劳， 导致 误诊或漏诊率的上升。随着计算机技术和图像处理技术突飞猛进的发展， 为这些诊断数据 的处理提供了新的方法， 计算机辅助诊断技术应运而生。 0003 计算机辅助诊断技术以医学

14、影像技术为基础， 结合计算机的处理分析功能， 帮助 影像医师发现病灶区域以提高临床诊断的效率。其中医学图像的特征提取能够获取图像 中独特的有别于其他图像的特征信息， 进而实现医学图像目标的自动分类、 匹配、 识别等功 能， 是实现计算机辅助诊断的一个关键步骤。 0004 医学图像数据的多样性和重要性要求运用高效准确的图像处理算法提取图像中 的特征信息。 有效的特征提取算法， 可以准确的分离出人体组织和病灶区域， 这样就可以减 轻医生的负担， 提高医生的诊断效率， 给医生的临床诊断提供更多的参考信息。目前， 常用 的图像特征提取方法主要有梯度算子、 拉普拉斯算子等， 这些算法通用性较差， 对噪声

15、也非 常敏感 ； 小波分析在图像处理中应用广泛， 该方法具有多尺度、 多分辨率的特性， 但小波基 的选择是关键， 采用不同的小波基， 分解效果会有所不同 ； 分形理论作为一种非线性的处理 方法， 也是图像特征提取中比较热门的理论。 这些算法大多都是直接在原图像上进行处理， 效果受噪声影响较大。而 Snake 算法、 人工神经网络方法、 统计方法等， 大多针对辅助诊断 或模式识别应用， 对图像特征区域划分较细， 因此算法一般过于复杂、 需要人工干预且时间 复杂度很大， 难以满足临床诊断的需求。 0005 经验模式分解 (Empirical Mode Decomposition， 简称 EMD)

16、方法由美国华裔科学 家 HUANG 于 1998 年提出， 方法扩展了 Hilbert 变换的应用， 突破了传统数据分析方法只能 分析线性或平稳数据的局限， 开创了一种处理非线性非平稳数据的有效方法， 能自适应处 理非线性非平稳数据的特性使其迅速在诸多领域得到了广泛的应用。近年来， 二维经验模 式分解 (Bi-Dimensional Empirical Mode Decomposition， BEMD) 在图像去噪、 图像增强、 特征提取图像分割、 图像融合和图像压缩等图像处理各个方面的应用研究也正处于不断深 入的过程之中， EMD 在图像处理领域的应用已成为学者们研究的热点。 说 明 书 C

17、N 104392444 A 5 2/14 页 6 0006 虽然， EMD 方法在各领域的应用已经取得了很大的成功， 但由于该方法提出不久， 其理论研究远未成熟， 在应用该方法时存在着端点效应、 模态混叠等各类问题。针对这些 问题不少学者提出了相应的解决方法。其中， 值得一提的是， 2009 年， HUANG 等人在分析 了白噪声统计特性的基础上提出了一种新的噪声辅助数据分析方法， 集合经验模式分解 (Ensemble Empirical Mode Decomposition， 简称 EEMD)。EEMD 方法是对原始 EMD 方法的 巨大改进， 这种方法通过给信号加入极小幅度白噪声， 利用白

18、噪声频谱均衡分布的特点， 用 白噪声来均衡噪声的特性， 较为理想的解决了模态混叠问题， 因而非常适合用于医学 MR 图 像的特征提取。 0007 发明目的 0008 本发明的目的在于针对上述现有技术中的缺陷， 提供一种能够解决大多数特征提 取算法对 MR 图像噪声较为敏感， 需要人工干预， 时间复杂度大， 在医学 MR 图像的特征提取 上效果不理想的问题， 进而提高 MR 图像特征提取效率及精度的基于二维集合经验模式分 解的医学 MR 图像特征提取方法。 0009 为了实现上述目的， 本发明采用的技术方案为， 对于二维图像信号 I(x,y) 进行如 下处理 ： 0010 步骤一， 初始化 r0

19、(x,y) I(x,y)， 定义求解 BIMF 分量次数的控制变量为 i， 并且 首次求解 BIMF 分量的 i 1 ； 0011 步骤二， 向原始信号中加入白噪声序列 ； 0012 步骤三， 通过筛分提取第 i 个 BIMF 分量， 具体操作如下 ： 0013 (1) 初始化 ： h0(x,y) ri-1(x,y)， 定义对应求解第 i 个 BIMF 分量时的筛选控制变 量为 j， 并且首次求解第 i 个 BIMF 分量的 j 1 ； 0014 (2) 利用极值点选择算法找出 hj-1(x,y) 中所有的极大值点与极小值点 ； 0015 (3) 使用曲面插值算法分别对 (2) 中的极大值点和

20、极小值点进行插值， 得到上包 络面 emax(x,y) 和下包络面 emin(x,y) ； 0016 (4) 计算上下包络面的均值 mj-1(x,y) emax(x,y)+emin(x,y)/2 ； 0017 (5) 得到筛选函数 hj(x,y) hj-1(x,y)-mj-1(x,y) ； 0018 (6)判断 ： 当hj(x,y)满足筛分停止条件， 则BIMFi(x,y)hj(x,y) ； 否则jj+1， 转到 (2) ； 0019 步骤四， 加入不同的白噪声序列重复步骤三 ； 0020 步骤五， 更新余量 ri(x,y) ri-1(x,y)-BIMFi(x,y) ； 0021 步骤六， 如

21、果 ri(x,y) 的极值点个数小于给定的阈值或分解达到指定层数时则算 法停止 ； 否则， 转到步骤二继续进行下一层的分解 ； 0022 步骤七， 将每次得到的对应 BIMF 分量的集成平均值作为最终的分解结果 ； 0023 步骤八， 对得到的BIMF分解结果进行Hilbert变换， 构建解析信号， 得到所有BIMF 图像的局部特征信息， 同时对第一个 BIMF 图像进行相位一致性的边缘特征提取得到图像 的边缘特征信息 ； 0024 步骤九， 综合分析图像的所有特征信息， 完成图像特征提取。 0025 所述的步骤三的第(6)步中筛分停止条件为满足BIMF分解结果的基本条件， 即通 过相邻两次筛

22、选结果的标准差进行判断， 标准偏差 SD 的计算公式为 ： 说 明 书 CN 104392444 A 6 3/14 页 7 0026 0027 其中， hi,j-1(m,n) 和 hi,j(m,n) 分别是第 i 个 BIMF 分解过程中两个连续处理结束的 时间序列， 二维 EMD 的标准偏差 SD 的阈值 设定在 0.1 0.3 之间， M*N 等于二维图像信 号 I(x,y) 的像素大小。 0028 所述的步骤六中 ri(x,y) 的极值点个数小于给定的阈值或分解达到指定层数时， 其中阈值设为 2 个， 得到 3 6 个 BIMF 分量以及 1 个余量， 则算法停止。 0029 所述的步骤

23、八中构建解析信号的具体操作方法为 ： 0030 a. 将一维信号 x(t) 通过 Hilbert 变换公式得到 0031 b. 构造复信号 0032 则 z(t) 称为 x(t) 的解析信号， 其中 ： 0033 0034 分别表示信号 x(t) 的局部振幅和局部相位。 0035 所述的 Hilbert 变换采用二维欧式空间向量值扩展后的 Riesz 变换。 0036 所述的 Riesz 变换求解局部特征的具体过程为 ： 0037 a. 根据 Riesz 核的空域表达式 ： 0038 0039 上式中 x 为信号分析时域， 即图像域， 进行 Fourier 域的变换后， 得到表达式 ： 004

24、0 u 为信号分析频域 ； 0041 b. 根据上述公式， 对于二维输入信号， 其单演信号为 ： 0042 fM(x,y) (f(x,y),Rxf(x,y),Ryf(x,y) (f,Rx*f,Ry*f) ； 0043 其中， * 表示卷积运算 ； 0044 c. 分别通过以下公式得到局部振幅 lA， 局部相位 lp， 局部方向 lo ： 0045 0046 0047 0048 通过局部相位， 进一步得到局部频率 lf为 ： 说 明 书 CN 104392444 A 7 4/14 页 8 0049 0050 所述的 Riesz 变换获取局部特征信息的具体步骤如下 ： 0051 步骤 1 ： 分解

25、原始图像得到 K 个 BIMF ； 0052 步骤 2 ： 计算第 i 个 BIMF 的 Riesz 变换得到相应的单演信号 mi， 其中， 1 i K ； 0053 步骤 3 ： 分别计算出 mi对应的局部振幅， 局部相位， 局部方向和局部频率 ； 0054 步骤 4 ： 重复步骤 2 及步骤 3， 直到所有 BIMF 的局部特征都得到求解。 0055 所述的对第一个 BIMF 进行相位一致性的边缘特征提取的方法为 ： 0056 a. 对输入图像进行 BEEMD 分解 ； 0057 b. 取第一个 BIMF 并计算其相位一致性度量 PC ； 0058 c. 对 PC 进行非极大值抑制得到边缘

26、图像 ； 0059 d. 对边缘图像进行滞后阈值处理得到最终的二至边缘图。 0060 所述相位一致性度量 PC 的计算公式为 ： 0061 其 中，从 余 弦 和 正 弦 小 波 的 卷 积 导 出 小 波 的 响 应和 分 别 计 算 出 尺 度 m*n 的 偶 分 量 em*n(x) 和 奇 分 量 om*n(x)，变 换 的 结 果 振 幅 得 到 局 部 能 量 ：对 偶 分 量 和 奇 分 量 求 和 得 到 ： 为机器能够显示的最小值。 0062 与现有技术相比， 本发明具有如下有益效果 ： 本发明针对目前临床上人工逐层解 读大量 MR 图像信息不但诊断工作量大， 而且容易因为医生

27、的疲劳而导致误诊或漏诊率上 升的状况， 研究基于二维集合经验模式分解 BEEMD 的医学 MR 图像特征提取问题， 首先通过 EMD 方法提取复杂信号在每个时间中的局部震荡模式， 将复杂信号分解为一系列有限个从 高频到低频排列的固有模态函数之和， 然后对每个 IMF 做 Hilbert 变换， 计算每个 IMF 的能 量， 即瞬时频率和振幅， 从而形成时间 - 频率 - 振幅谱， 即 Hilbert 谱， 提取信号的瞬时频率 特征。 本发明将集合经验模式分解算法由一维推广到更加复杂的二维空间得到二维集合经 验模式分解算法， 以此来抑制 BEMD 算法中出现的模态混叠现象， 从而获取到图像的局部

28、幅 值、 局部相位、 局部方向和局部频率 ； 在 BEEMD 分解的基础上再结合相位一致性算法得到 MR 图像的边缘特征， 为下一步特征分类和识别提供依据。 结合临床应用背景， 本发明重点克服 了现有的大多数特征提取算法对 MR 图像噪声较为敏感， 需要人工干预且时间复杂度大， 在 医学MR图像的特征提取上效果不理想的问题， 同时提高了MR图像特征提取的效率及精度。 通过对 lena 图像以及临床 MR 乳腺和大脑图像进行分解和特征提取， 提取结果和 BEMD 的提 取结果进行对比， 对比结果说明本发明在性能能上明显优于 BEMD， 能客观地提取和分析 MR 图像， 为 MR 图像的处理提供新

29、的理论与方法， 使医疗专家摆脱繁重的人工观察和诊断 ； 同 说 明 书 CN 104392444 A 8 5/14 页 9 时能提供更精确的辅助诊断数据， 降低医生诊断的主观性， 从而提高诊断正确性和客观性。 0063 进一步的， 本发明在筛分停止判断中二维标准偏差SD的阈值设定在0.10.3 之间， 由于二维标准偏差SD目前还没有明确标准， 很大程度上要靠经验值， 若取值过小， 停止条件则过于严格， 会导致分解的迭代次数过多、 计算量大， 图像过度分解结果不够理 想 ； 若 取值偏大， 停止条件宽松， 造成分解结果尺度特征分离不明显。实验中 我们取 0.2， 一方面它的大小合适， 可以作为一

30、个较准确的筛选值 ； 另一方面这个值不至于过小， 可 以降低迭代次数和计算量。 附图说明 0064 图 1 本发明图像特征提取方法的整体流程图 ； 0065 图 2 合成信号的 EMD 分解及其 Hilbert 谱图 ： 其中， 图 2(a) 为由频率为 20Hz 和 100Hz 的余弦信号叠加而成的信号的 EMD 分解图 ； 图 2(b) 为频率为 20Hz 和 100Hz 的余弦 信号叠加而成的信号的 EMD 分解图对应的 Hilbert 谱图 ； 0066 图 3(a) 为待合成的两个模拟信号 x1 与 x2 的波形图 ； 0067 图 3(b) 为 x1 与 x2 的 EMD 算法分解

31、图 ； 0068 图 3(c) 为 x1 与 x2 的 EEMD 算法分解图 ； 0069 图 4 本发明 BEEMD 的算法流程图 ； 0070 图 5 本发明结合 BEEMD 和相位一致性算法的图像边缘特征提取流程图 ； 0071 图 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 图 ；

32、 图 6(g) 为余量图 ； 图 6(h) 为 BIMF1 5 的重建 图 ； 0072 图 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 的重建图 ； 0073 图 8 为乳腺 MR 图像的原图 ； 0074 图 9 为乳腺 MR

33、 图像的 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 的重建图 ； 0075 图 10 为乳腺 MR 图像的 BEMD 分解图 ： 其中， 图 10(a) 为 MR 图像原图 ； 图 10(b) 为 BEMD 的 BIMF1 图 ； 图 10

34、(c) 为 BEMD 的 BIMF2 图 ； 图 10(d) 为 BEMD 的 BIMF3 图 ； 图 10(e) 为 BEMD 的 BIMF4 图 ； 图 10(f) 为 BEMD 的 BIMF5 图 ； 图 10(g) 为余量图 ； 图 10(h) 为 BIMF1 5 的重建图 ； 0076 图 11 为大脑 MR 图像的原图 ； 0077 图 12 为大脑 MR 图像的局部特征提取 BEMD 与 BEEMD 对比图， 其中， 第一行为基于 BEMD 的提取结果， 第二行为基于 BEEMD 的提取结果， 图 12(a) 为局部幅值对比图 ； 图 12(b) 为局部相位对比图 ； 图 12(

35、c) 局部方向对比图 ； 图 12(d) 为局部频率对比图 ； 0078 图 13 为乳腺 MR 图像的局部特征提取 BEMD 与 BEEMD 对比图， 其中， 第一行为基于 说 明 书 CN 104392444 A 9 6/14 页 10 BEMD 的提取结果， 第二行为基于 BEEMD 的提取结果， 图 13(a) 为局部幅值对比图 ； 图 13(b) 为局部相位对比图 ； 图 13(c) 局部方向对比图 ； 图 13(d) 为局部频率对比图 ； 0079 图 14 为乳腺 MR 图像的边缘特征提取结果图 ： 图 14(a) 为 Prewitt 算子检测结果 图 ； 图 14(b) 为 R

36、oberts 算子检测结果图 ； 图 14(c) 为 Sobel 算子检测结果图 ； 图 14(d) 高 斯拉普拉斯算子检测结果图 ； 图14(e)为Canny算子检测结果图 ； 图14(f)本发明检测结果 图 ； 0080 图 15 为大脑 MR 图像的边缘特征提取结果图 ： 图 15(a) 为 Prewitt 算子检测结果 图 ； 图 15(b) 为 Roberts 算子检测结果图 ； 图 15(c) 为 Sobel 算子检测结果图 ； 图 15(d) 高 斯拉普拉斯算子检测结果图 ； 图15(e)为Canny算子检测结果图 ； 图15(f)本发明检测结果 图。 具体实施方式 0081 为

37、了使本领域技术人员更好地理解本申请中的技术方案， 下面结合附图对本申请 技术方案做进一步的详细说明， 所描述的实施例仅仅是本申请的部分实施例， 而不是全部 的实施例。基于本申请中的实施例， 本领域普通技术人员在没有做出创造性劳动的前提下 所获得的所有其它实施例， 都应当属于本发明保护的范围。 0082 参见图 1， 本发明基于二维集合经验模式分解的医学 MR 图像特征提取方法， 对于 二维图像信号 I(x,y) 进行如下处理 ： 0083 步骤一， 初始化 r0(x,y) I(x,y)， 定义求解 BIMF 分量次数的控制变量为 i， 并且 首次求解 BIMF 分量的 i 1， 对 BIMF

38、分量进行 N 求解 ； 0084 步骤二， 向原始信号中加入白噪声序列 ； 0085 步骤三， 通过筛分提取第 i 个 BIMF 分量， 具体操作如下 ： 0086 (1) 初始化 ： h0(x,y) ri-1(x,y)， 定义对应求解第 i 个 BIMF 分量时的筛选控制变 量为 j， 并且首次求解第 i 个 BIMF 分量的 j 1 ； 0087 (2) 利用极值点选择算法找出 hj-1(x,y) 中所有的极大值点与极小值点 ； 0088 (3) 使用曲面插值算法分别对 (2) 中的极大值点和极小值点进行插值， 得到上包 络面 emax(x,y) 和下包络面 emin(x,y) ； 008

39、9 (4) 计算上下包络面的均值 mj-1(x,y) emax(x,y)+emin(x,y)/2 ； 0090 (5) 得到筛选函数 hj(x,y) hj-1(x,y)-mj-1(x,y) ； 0091 (6)判断 ： 当hj(x,y)满足筛分停止条件， 则BIMFi(x,y)hj(x,y) ； 否则jj+1， 转到 (2) ； 筛分停止条件为满足 BIMF 分解结果的基本条件， 即通过相邻两次筛选结果的标 准差进行判断， 标准偏差 SD 的计算公式为 ： 0092 0093 其中， hi,j-1(m,n) 和 hi,j(m,n) 分别是第 i 个 BIMF 分解过程中两个连续处理结束 的时间

40、序列， 二维 EMD 的标准偏差 SD 的阈值 设定在 0.1 0.3 之间， 实施例 取 0.2， M*N 等于二维图像信号的像素大小 ； 0094 步骤四， 加入不同的白噪声序列重复步骤三 ； 说 明 书 CN 104392444 A 10 7/14 页 11 0095 步骤五， 更新余量 ri(x,y) ri-1(x,y)-BIMFi(x,y) ； 0096 步骤六， 如果 ri(x,y) 的极值点个数小于给定的阈值或分解达到指定层数时， 得到 3 6 个 BIMF 分量以及 1 个余量， 则算法停止 ； 极值点个数给定的阈值为 2 个 ； 0097 步骤七， 将每次得到的对应 BIMF

41、 的集成平均值作为最终的分解结果 ； 0098 步骤八， 对得到的 BIMF 分解结果进行 Riesz 变换， 构建解析信号， 得到局部振幅、 相位、 频率等所有BIMF图像的局部特征信息， 同时对第一个BIMF图像进行相位一致性的边 缘特征提取得到图像的边缘特征信息 ； 0099 步骤九， 综合分析图像的所有特征信息， 完成图像特征提取。 0100 1.EMD 算法基本原理为 ： 0101 经验模式分解 (EMD) 是一种直观的， 直接的、 先验的、 自适应的分解方法， 方法根据信号的特征时间尺度将信号分解为一组固有模态函数 (Intrinsic Mode Function,IMF)， 或称

42、基本模式分量的和。每个 IMF 分量必须满足两个条件 ： 0102 1) 在整个信号序列中， 极值点 ( 即大值点与极小值点 ) 的个数 Ne和过零点数目 Nz 必须相等或最多相差不超过一个， 即 0103 (Nz-1) Ne (Nz+1) (1) 0104 2) 任意时间点 ti上， 信号序列局部最大值所确定的上包络线 Smax(t) 与局部最 0105 小值所确定的下包络线 Smin(t) 关于时间轴局部对称， 即均值为零 ： 0106 (Smax(t)+Smin(t)/2 0 (2) 0107 IMF表征了信号内在的波动模式。 由上述定义可知， 在IMF每一个波动周期(两零 点之间 )

43、只有一个单纯的波动模式， 没有其它波动叠加上去， 因此它可以作为分解信号的 基本量， 也就是说任意信号可以被分解为若干个 IMF。EMD 算法就是这样一个 “筛选” ( 的过 程， 具体由以下几步组成 ： 0108 步骤 1 ： 找出信号 x(t) 的所有局部极大值和局部极小值 ； 0109 步骤 2 ： 利用三次样条插值方法分别计算极小值插值和极大值插值并得到对应信 号包络 emin(t) 和 emax(t) ； 0110 步骤 3 ： 计算局部均值 m(t) (emin(t)+emax(t)/2。 0111 步骤 4 ： 将原始输入信号减去局部均值得到振荡信号 ： 0112 h(t) x(

44、t)-m(t)。 0113 步骤 5 ： 当 h(t) 满足 IMF 的条件时， h(t) 就成为一个 IMF ； 否则， 用 h(t) 替换步骤 1 中的 x(t) 并从步骤 1 开始重复上述过程。 0114 步骤 6 ： 令 c1 h(t)， 则 c1为第一个 IMF， 对应的余量 r1 x(t)-c1。 0115 当迭代终止条件不满足时， 用 r1替换原信号并重复上述所有步骤并得到第二个 IMF 分量 c2及余量 r2， 如此反复， 最后得到第 n 个 IMF、 cn和 rn。当满足如下预先约定的条 件时， 停止筛分过程 ： cn或 rn小于预先设定的值， 或者 rn变成一个单调函数。

45、0116 经过上面的分解过程， x(t) 最终被分解为 n 个 IMF 分量和一个余量 rn。原始信号 x(t) 可以表示为对分解后得到的各 IMF 分量分别求 Hilbert 变换， 得到 对应的瞬时频率和瞬时振幅， 进而得到 Hilbert 谱和 Hilbert 边缘谱。图 2(a) 为由频率为 20Hz 和 100Hz 的余弦信号叠加而成的信号 x1(t) cos(220t)+cos(2100t) 说 明 书 CN 104392444 A 11 8/14 页 12 的 EMD 分解结果， 图 2(b) 为其对应的 Hilbert 谱。 0117 2.EEMD 算法实现 0118 原始 E

46、MD 算法的主要缺点之一就是模态混叠的频繁出现。模态混叠有两种表现形 式 ： (1)单个IMF中包含不同的时间尺度成份 ； (2)同一尺度出现在不同的IMF中。 模态混叠 的出现会导致产生严重错假的时频分布也使IMF失去物理意义。 集合经验模式分解(EEMD) 主要解决模态混叠的问题， 其分解步骤如下 ： 0119 步骤 1 ： 向原始信号中加入白噪声序列 ； 0120 步骤 2 ： 将添加了白噪声的信号通过 EMD 算法分解为一系列的 IMF ； 0121 步骤 3 ： 重复步骤 1、 步骤 2， 但每次加入不同的白噪声序列 ； 0122 步骤 4 ： 将每次得到的对应 IMF 的集成平均值

47、作为最终的分解结果。 0123 EEMD 方法需要两个重要的参数 ： 所添加白噪声的幅值和集成平均次数。WU 和 HUANG 通过实验总结得出只要添加的噪声经过频率调制而且集成平均次数足够大， 那么噪 声幅值和集成平均次数的增加几乎就不会改变分解结果， 并建议大多数情况下添加噪声的 幅值应该满足其大约为 0.2 个数据信号的标准差。EEMD 方法是对原始 EMD 方法的巨大改 进， 这种方法通过给信号加入极小幅度白噪声， 利用白噪声频谱均衡分布的特点， 用白噪声 来均衡噪声的特性， 较为理想的去除了模态混叠。 0124 原始经验模式分解与 EEMD 分解的结果对比如图 3(a) (c) 所示。

48、其中， 图 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在 抗模态混叠方面的有效性。 0125 3.BEEMD 算法实现 0126 1).BEMD 算法实现 ： 0127 将一维 EMD 算法推广到更加复杂的二维空间便得到了二维经验模式分解算法 BEMD。与一维经验模式分解算法相似， 其分解步骤也是通过筛分过程来完成， 每次筛分 完成后得到一个对应的模态函数。分解完成后得到数量有限的几个二维固有模态函数 (Bi-Dimensional Intrinsic Mode Function， 简称 BIMF) 及一个余量。对二维图像信号 I(x,y) 的