CN102855623A - 基于经验模态分解的心肌超声造影图像生理参数测量方法 - Google Patents

Abstract

基于经验模态分解的心肌超声造影图像生理参数测量方法，涉及一种心肌超声造影图像的时间-强度曲线分析方法。本发明包括如下步骤：一、对心肌区域按圆周分为六个区域进行灌注信号提取，分别提取每一个心肌分割区域的时间-强度曲线；二、对提取出的时间-平均灰度值强度曲线进行EMD分解，得到扩展的信号本征模态函数，对得到的第一个本征模态函数进行Hilbert变换，得到瞬时谱参数，对每个频率点的幅值进行累加，得到边际谱；三、利用边际谱的能量阈值对边际谱进行阈值判定。本发明为基于心肌超声造影图像的心肌梗塞或心肌微循环诊断提供了一种分析途径，从而提高基于心肌超声造影的生理参数测量的准确性，降低了医生的主观依赖性。

Description

基于经验模态分解的心肌超声造影图像生理参数测量方法

技术领域

本发明涉及一种心肌超声造影图像的时间-强度曲线分析方法，具体涉及一种利用经验模态分解进行时间-强度曲线的特征提取及利用希尔伯特边际谱进行生理参数的测量方法。

背景技术

心脏是人体的重要器官，随着心脏每次收缩，将携带氧气和营养物质的血流经主动脉输送到全身，以供给各组织细胞代谢需要。由于脂质代谢不正常，血液中的脂质沉着在原本光滑的冠动脉内膜上，形成血栓。微血栓一旦形成，会使微血管发生严重狭窄甚至阻塞，产生微循环灌注异常，严重者发生心肌急性缺血、坏死、甚至危及生命。即便采用冠脉血运重建术进行治疗，也大多存在心肌“低灌注”或“无灌注”异常，即心脏在缺血后再灌注时心肌存在微循环障碍。缺血心肌微循环灌注异常将严重影响受累心肌功能的恢复，导致病人心脏功能减退。

心肌超声造影（Myocardial Contrast Echocardiography, MCE）是一种影像新技术，通过将含有微气泡的对比剂注入冠状循环，来观察二维超声心动图上的心肌显影。由于微气泡通过心肌时完全保持在血管内，且微泡的大小及变形性与红细胞相当，因此可将其视作红细胞流动的示踪剂，故MCE可用来在跳动的心脏上估价心肌微循环状况。目前临床上常用定性法分析MCE，该方法虽简便、直观，但很大程度上依赖观察者的经验，缺乏客观定量分析。超声造影图像感兴趣区域的灌注强度曲线可以反映器官、组织的血流灌注情况。该技术可以在不破坏造影剂气泡的前提下，实时动态地显示组织的造影剂灌注，使精确的定量分析成为可能。通过对造影剂灌注的时间-强度曲线进行实时探测，获得反映心肌血流灌注的生理参数，以便为心肌梗塞或心肌微循环诊断提供定量数据。

经验模态分解(Empirical Mode Decomposition，简称EMD)方法是由美国NASA的黄锷博士提出的一种信号分析方法。它不会像傅立叶分解或小波分解那样，预先设定好分解的基函数，而是会根据信号自身的时间尺度特征来进行信号分解。EMD 方法在理论上可以应用于任何类型的信号的分解，尤其在处理非线性、非平稳信号上，具有非常明显的优势。利用心肌超声造影图像提取出的灌注强度曲线的内部时间尺度的变化做能量与频率的解析，将该曲线展开成数个本征模态函数（Intrinsic Mode Function，IMF）的形式，再利用希尔伯特变换（Hilbert Transform，HT）获得IMF的瞬时频率和振幅，从而对提取出的时间-强度信号进行分析，从而获得心肌灌注图像的生理参数，为诊断提供可能。

EMD是HHT（Hilbert-Huang Transform）算法的重要步骤，它认为信号是由许多个不同的模态组合而成的，可以是线性或非线性，分别对应不同过程，而且要求其相对于横坐标轴是对称的。要成为IMF内在模态函数有两个条件：

1）在整个函数中，极值点的数目与穿越零点的数目相等或者相差1；

2）在任何时刻，由局部极值包络线所定义的包络线局部均值为零。

其中，第一个条件与传统高斯平稳过程中窄频宽的要求类似。第二个条件将整体性要求改变为局部性要求，使得瞬时频率不会因为不对称波形的存在而导致不必要的晃动。依托这两个条件构建起来的EMD及HHT被认为是强有力地求解非线性、非平稳信号的自适应方法，是近年来对以傅立叶变换为基础的线性及稳态谱分析的重大突破，并得到了广泛的应用。

临床上经常使用拟合曲线对造影剂灌注的心肌参数进行定性测量，但是这样的测量方法非常依赖于医生的经验，而且存在一定的测量误差；同时由于拟合曲线只反映了时间-平均灰度曲线的变化趋势信息，而曲线中蕴含的更多信息却无法通过拟合曲线而表现，所以通常情况下健康的心肌区域和心肌梗塞区域的拟合曲线是基本类似的，测量的生理参数无法表现待诊断区域的性质，因此利用拟合曲线进行生理参数测量的临床实用价值受到一定的限制。

发明内容

为了解决了当前心肌超声造影图像定性理解的主观依赖性，本发明提出了一种基于经验模态分解的心肌超声造影图像生理参数测量方法，利用经验模态分解进行时间-强度曲线的特征提取，通过特征挖掘以及希尔伯特边际谱，来进行心肌灌注情况下的生理参数测量方法，从而为基于心肌超声造影图像的心肌梗塞或心肌微循环障碍诊断提供了另一种分析途径。

本发明是通过以下技术方案实现的：

一、图像信号提取部分：对心肌造影图像的心肌区域按圆周顺时针方向分为六个区域进行灌注信号提取，这六个区域分别代表前壁心肌区域、侧壁心肌区域、后壁心肌区域、下壁心肌区域、后间隔心肌区域和前间隔心肌区域，分别提取每一个心肌分割区域的时间-强度曲线，曲线以横时间为坐标，每个区域心肌的平均灰度值为纵坐标；

二、信号特征提取部分：对提取出的时间-平均灰度值强度曲线进行EMD分解，得到扩展的信号本征模态函数，对得到的第一个本征模态函数进行Hilbert变换，得到瞬时谱参数，包括瞬时幅值、瞬时频率和瞬时相位，对每个频率点的幅值进行累加，得到Hilbert-Huang边际谱；

三、阈值判定部分：利用Hilbert边际谱的能量阈值对边际谱进行阈值判定，如果边际谱能量集中且大于阈值，可以得到造影剂灌注过程中的生理参数，包括呼吸频率和心跳频率参数。

本发明具有如下优点：

1 、对心肌图像进行六区域分割，可以兼顾临床诊断及处理信息量的要求，利用区域内的像素平均值可以反映出该区域内的造影剂灌注情况，从而能够利用时间-平均像素强度曲线来描述心肌超声造影过程所携带的信息。经过该信息提取过程，可以很好地将造影过程中的生理参数信息反映出来。

2、对提取出的时间-强度曲线进行经验模态分解，得到能够反映该曲线特征的本征模态函数，由于EMD会根据信号自身的时间尺度特征来进行信号分解，所以得到的IMF会含有丰富的频率信息。利用HT的瞬时频率概念得到希尔伯特谱，通过对其在时间域积分，利用边际谱来研究在呼吸和心跳不同频率内幅度值的变化情况，从而为心肌灌注的参数提取提供可能。

3、Hilbert边际谱可以反映出信号能量在频率轴上的分布情况，利用Hilbert边际谱的能量阈值对边际谱进行阈值判定，如果边际谱能量集中且大于阈值，可以得到造影剂灌注过程中的呼吸频率和心跳频率参数，为心肌梗塞和心肌微循环提供定量依据。

4、通过对SD大鼠的前壁心肌区域进行结扎后，通过心肌超声造影进行观察，并利用上述所给分析方法进行信号提取及特征提取。分析所得边际谱可以看出，对应正常心肌区域的灌注生理参数，即呼吸频率范围为1~2Hz、心跳频率为7~8Hz，与正常心肌情况下大鼠的基本生理参数一致。对于心肌梗塞区域，从灌注过程中得不到明确的心肌活动生理参数，与心肌坏死的前提一致，从而证明了本发明方法进行心肌灌注的生理参数测量的正确性。

附图说明

图1为基于经验模态分解进行大鼠心肌超声造影图像时间-强度曲线的分析方法流程图；

图2为经验模态分解流程图；

图3为心肌的六区域分割图；

图4为大鼠六个心肌区域的时间-强度曲线；

图5为心肌区域1的时间-强度曲线的EMD分解结果；

图6为心肌区域3的时间-强度曲线的EMD分解结果；

图7为心肌区域1的IMF1中对应的各个瞬时参数；

图8为心肌区域3的IMF1中对应的各个瞬时参数；

图9为区域1的边际谱；

图10为区域3的边际谱；

图11为区域1的边际谱阈值判断图；

图12为区域3的边际谱阈值判断图。

具体实施方式

具体实施方式一：如图1所示，本实施方式的基于经验模态分解的心肌超声造影图像生理参数测量方法包括如下步骤：首先对心肌造影图像的心肌区域按顺时针分为6个区域，分别提取每一个区域在心肌灌注过程中的时间-强度曲线；利用EMD方法对时间-强度曲线进行分解，得到对应不同频率的本征模态函数；分别对本征模态函数做Hilbert变换，得到各种瞬时特征参数，并最终得到Hilbert边际谱，根据能量阈值判断，可以得到心肌造影图像所对应的心肌灌注情况下的生理参数，包括呼吸频率和心跳频率范围，从而为基于心肌超声造影的心肌梗塞或心肌微循环诊断提供定量依据。具体步骤如下：

步骤一：图像信号提取部分，为了满足临床上对于心肌梗塞诊断的要求，将心肌按圆周分为六个区域进行灌注分析，这六个区域分别代表心肌的不同部分。六个心肌区域按顺时针分别标号为Region-1至Region-6，如图3所示。Region-1，即顺时针1点到3点的部分（1点为时钟在1点时的位置，依次类推）对应实际心肌的前壁心肌区域；Region-2，即顺时针3点到5点的部分对应心肌的侧壁心肌区域；Region-3，即顺时针5点到7点的部分对应心肌的后壁心肌区域；Region-4，即7点到9点的部分对应心肌的下壁心肌区域； Region-5，即9点到11点的部分对应心肌的后间隔心肌区域； Region-6，即11点到1点的部分对应心肌的前间隔心肌区域。

心肌区域分割的目的是为了对每个区域进行分别处理，每个区域分别代表了心肌的不同部分，在临床中通常仅对不同区域的心肌进行微循环灌注分析，通过计算每一个区域的平均灌注强度信息，来尽量反映该区域内的心肌灌注情况。实践证明，将心肌进行分区处理能够明显提高灌注生理参数测量的准确性。

时间-强度曲线提取是将分割出来的心肌区域进行时间-强度曲线绘制，由于心肌超声造影图像的帧顺序反映了时间信息，而每一帧中的每个心肌区域中的平均像素强度反映了在该时刻上造影剂的灌注情况，所以以时间为横坐标，每个区域心肌的平均灰度值为纵坐标可以做出时间-强度曲线，该曲线可以反映心肌内造影剂的灌注情况。直接使用平均灰度值绘制曲线的方法，能够完全保留心腔区域灰度变化的信息，而且没有引入心腔区域灰度变化的干扰、简单直观，后面的研究都是基于这种平均灰度值曲线的绘制方法进行的。

步骤二：信号特征提取部分，对提取出的时间-平均灰度值强度曲线进行EMD分解，得到扩展的信号本征模态函数。对经过提取的信号进行经验模态分解(EMD)，得到扩展的特征曲线，用本征模态函数(IMF)表示。要成为IMF有两个条件： (1)极值及过零点相差不可以超过1个；(2)上下包络线和为0。要完全严格的满足这两条中的第(2)条是不可能的，所以只能要求其尽量小。可以通过如下的分解方法得到IMF。经验模态分解的流程图如图2所示，具体描述如下：

为了得到分解的本征模态函数IMF，定义                                               

，其中l=1,…,L.代表了IMF的序号，k=1,…K.代表了筛选次数。筛选过程从输入信号开始，定义

，为了得到下一个IMF，我们将前一个IMF函数的余项作为输入信号，

。对于给定信号

的IMF筛选过程，包括以下步骤：

(1)程序初始化，并搜索信号的振幅极大值及振幅极小值。

(2)使用三次样条插值得到上下包络线，分别记为

和

。并计算上下包络的均值：

；

(3)通过前一级信号与包络均值计算出后一级信号：

；

检查步骤(3)得到的信号

是否满足IMF的条件，即当包络均值足够接近零时迭代终止：

；

(4)通过判定准则来判断信号均值是否足够趋近于零，如果不符合则重复步骤1至3。当满足判定准则时，IMF函数

定义为步骤3最后一次的结果，即

；

定义余项：

                                 ；

(5)计算下一个IMF需要重新初始化，将余项作为输入信号，通过步骤1进行重复筛选，

；

对于所有的重复步骤1到5，当余项不包含极值点是EMD分解结束。输入信号可以表示成IMF函数和余项的和：

                              

  。

IMF1中含有时间-强度曲线的丰富频率信息，因此利用Hilbert变换对得到的IMF1进行处理，得到IMF1对应的瞬时谱参数，包括瞬时幅值、瞬时频率和瞬时相位，通过对这些瞬时参数的分析可以获得信号的瞬时特性。频域一定程度上比时域更加能体现出信号的特征，尤其利用Hilbert变换，首先构造原信号的解析信号，在对其进行变换，可以得到瞬时频谱信息。

经过EMD分解得到的IMF 为

，该IMF仍然为时域信号，记为

，其Hilbert变换为

：

；

组合可以得到一个复数

：

；

其中

和

分别表示信号的瞬时幅值和瞬时相位。

；

瞬时频率

便是对相位求导，即：

。

对于一个含有复杂频率的信号，无法使用这种方法进行有效的分析，所以引入经验模态分解方法将信号分解为不同的IMF，然后分别对不同的IMF进行Hilbert变换，就可以得到希尔伯特谱和边际谱。

Hilbert-Huang边际谱表达了每个频率在全局上的幅度(或能量)贡献，它代表了在统计意义上的整组数据在每个频率点的累积幅值分布，也就是说边际谱可以代表信号在每一个频率处的能量累积的情况，这样就可以看出信号能量在频率轴上的分布情况，并可以据此观察是否有能量的分布规律。对EMD分解后的IMF做Hilbert变换得到：

其中

就是希尔伯特谱，其中

和

分别为前面计算出的信号瞬时幅值和瞬时相位，Re表示求取该复数的实部。得到了希尔伯特谱之后，对其在时间域积分，便得到边际谱：

。

步骤三：利用阈值进行判定，Hilbert边际谱可以反映出信号能量在频率轴上的分布情况，并可以据此观察是否有能量的分布规律。对每个心肌区域的IMF分量的瞬时幅值在频率轴上进行累加，得到各个区域的边际谱。心率和呼吸频率是对心肌造影图像影响最重要的生理参数，且这两个参数的频率范围受受试者心肌超声造影灌注过程的影响，考虑到呼吸和心跳对整个灌注过程的作用，所以在边际谱中会反映出该生理参数信息。

结合边际谱的能量阈值，可以得到心肌灌注过程的呼吸和心跳参数，总结如下：

1）边际谱有两个大的能量幅值，且超过给定阈值之上，则该两个能量阈值范围对应的频率点分别对应呼吸频率点和心跳频率点范围。

2）边际谱只含有一个超过能量阈值的频率段，则该能量阈值范围对应的频率点为呼吸频率范围。

3）边际谱中不含有超过能量阈值的频率段，由于心跳频率与呼吸频率都是通过血流的变化在边际谱中反映出来的。心跳频率的体现是通过心脏向心肌中输送血液产生的，呼吸频率的体现是通过胸腔内压力随着呼吸变化而对心肌中的血管挤压而产生的。如果边际谱中不含有超过能量阈值的频率段，则表明待测的生理参数在灌注过程中未被体现，也由此可以得出心肌梗塞严重的结论。

通过对边际谱的分析可以对心肌灌注过程中的生理参数进行半定量测量：当心肌灌注过程的Hilbert边际谱有两处超过能量阈值时，可以得到对应1）的两个生理参数信息，即呼吸频率和心跳频率范围。当Hilbert边际谱有一处超过能量阈值时，可以得到对应2）的一个生理参数信息，即呼吸频率范围。当Hilbert边际谱未有能够超过能量阈值时，即符合3）所描述的情况，表明灌注过程未能体现任何生理参数信息。该生理参数的测量结果可以为判定正常心肌或心肌梗塞区域提供诊断依据。

具体实施方式二：为了评估利用经验模态分解进行心肌超声图像时间-强度曲线分析，及利用边际谱进行心肌灌注过程生理参数的测量研究，本实施方式以SD大鼠的心肌超声造影图像为例，对上述过程进行具体实施。

执行步骤一：对SD大鼠的心肌超声造影图像进行六区域分割，如图3所示，分别表示大鼠心肌的前壁心肌区域、侧壁心肌区域、后壁心肌区域、下壁心肌区域、后间隔心肌区域和前间隔心肌区域。

提取大鼠心肌超声造影图像的时间-强度曲线。以时间为横坐标（每一个超声图像帧），六个分割区域的平均灰度值为纵坐标，绘制时间-强度曲线，如图4所示。图中绘制出的六个区域的时间-强度曲线，从曲线上来看数据杂乱，并不能分析出该区域的心肌灌注状况，也无法直接从时间-强度曲线中得到心肌灌注过程的生理参数。

执行步骤二：对SD大鼠的心肌分割区域1和3（分别代表大鼠的前壁心肌区域和后壁心肌区域）进行EMD分解，得到每个区域对应的六个IMFs及残差，如图5和6所示。EMD分解得到的IMFs分别代表信号的不同的频率特征，相当于经过EMD分解将信号中的不同频率成分分解出来，而最后的残差即余项代表的是信号的趋势。虽然每一个心肌区域的曲线分解得到的是不同的IMF分量，但是直接从IMF分量中无法直接体现心肌灌注过程的生理参数特征和规律，因此需要借助Hilbert变换来进一步得到瞬时特征参数。

图7和8分别给出了区域1和区域3中对应的IMF1经过Hilbert变换后得到的各种瞬时参数，其中A表示瞬时幅值，IP表示瞬时相位，IF表示瞬时频率。

利用Hilbert-Huang边际谱理论对区域1和区域3分别进行能量累加，得到边际谱，如图7所示。

执行步骤三：设定能量阈值为0.2，将图9和10所示的SD大鼠的前壁心肌区域和后壁心肌区域的边际谱进行阈值判定，如图11和12所示，其中水平的红线为阈值范围，纵向的红线即表明该边际谱是否包含超过阈值的频率段。图12中的后壁心肌区域有四条红色的纵向线，表明后壁心肌区域被能量阈值限定后，边际谱的能量主要集中在1~2Hz和7~8Hz两个频率范围内。因为正常情况下呼吸频率会低于心跳频率，因此可以得到心肌灌注过程中的两个生理参数范围，即1~2Hz的频率范围为呼吸频率范围，7~8Hz为心跳频率范围。心跳和呼吸分别以心跳频率和呼吸频率对心肌超声造影过程产生影响，这种影响蕴含在超声造影图像上的心肌区域的灰度变化中，通过心肌区域灰度变化曲线经过EMD后得到的边际谱经过阈值限定后，可以确定其灌注过程的呼吸和心跳频率参数。

表1给出了SD大鼠的基本生理参数，即在未做心肌超声造影，也未对大鼠做任何心脏结扎手术时的正常状态下的生理参数。其中心跳每分钟测量结果在475次/分，转换为秒是7.9次/秒，即正常情况下SD大鼠的心跳频率为7.9Hz左右。同样，受测SD大鼠的呼吸频率为85.5次/分，折合成秒为1.5次/秒，即SD大鼠的呼吸频率为1.5Hz。而且该SD大鼠在做心肌超声造影时，对它的前壁心肌区域进行了结扎手术，也就是该前壁心肌区域存在严重的心肌梗塞，但是其他区域心肌活动正常。

表1 SD大鼠的基本生理参数

通过与大鼠的基本生理参数相比较，在心肌超声造影过程中，通过对造影图像进行信号提取并进行特征挖掘，得到的灌注过程的心跳频率和呼吸频率范围（1~2Hz和7~8Hz）与SD大鼠的基本生理参数基本一致，这个结果也与未对大鼠的后壁心肌区域的做任何结扎，其灌注活动正常的条件一致，说明经过分析得到的参数是正确的。

同样的过程，对大鼠前壁心肌区域的边际谱进行阈值限定，如图11所示，由于该心肌区域边际谱比较平坦，且整体的幅值很小，加上阈值之后未见有与阈值相交的频率范围，因此说明在心肌造影灌注过程中，该心肌未有正常的生理活动，不能得到灌注过程中的生理参数信息。与造影前对大鼠的结扎手术前提比较，由于对大鼠前壁心肌区域进行了结扎，前壁心肌区域有严重的心肌梗塞，导致前壁心肌坏死，心肌灌注的像素灰度值几乎没有太大波动，边际谱中无法体现灌注过程中的大鼠心率与呼吸频率的信息，与所设条件一致。

Claims (5)

1.基于经验模态分解的心肌超声造影图像生理参数测量方法，其特征在于所述方法包括如下步骤：

一、图像信号提取部分：对心肌造影图像的心肌区域按圆周顺时针方向分为六个区域进行灌注信号提取，这六个区域分别代表前壁心肌区域、侧壁心肌区域、后壁心肌区域、下壁心肌区域、后间隔心肌区域和前间隔心肌区域，分别提取每一个心肌分割区域的时间-强度曲线，曲线以横时间为坐标，每个区域心肌的平均灰度值为纵坐标；

二、信号特征提取部分：对提取出的时间-平均灰度值强度曲线进行EMD分解，得到扩展的信号本征模态函数，对得到的第一个本征模态函数进行Hilbert变换，得到瞬时谱参数，包括瞬时幅值、瞬时频率和瞬时相位，对每个频率点的幅值进行累加，得到Hilbert-Huang边际谱；

三、阈值判定部分：利用Hilbert边际谱的能量阈值对边际谱进行阈值判定，如果边际谱能量集中且大于阈值，可以得到造影剂灌注过程中的生理参数，包括呼吸频率和心跳频率参数。

2.    根据权利1所述的基于经验模态分解的心肌超声造影图像生理参数测量方法，其特征在于所述瞬时幅值                                               

、瞬时相位

和瞬时频率

的表达式为：

；

。

3.    根据权利1所述的基于经验模态分解的心肌超声造影图像生理参数测量方法，其特征在于所述边际谱表达式为：

。

4.根据权利1所述的基于经验模态分解的心肌超声造影图像生理参数测量方法，其特征在于所述边际谱有两个大的能量幅值，且超过给定阈值之上，则该两个能量阈值范围对应的频率点分别对应呼吸频率点和心跳频率点范围。

5.根据权利1所述的基于经验模态分解的心肌超声造影图像生理参数测量方法，其特征在于所述边际谱只含有一个超过能量阈值的频率段，则该能量阈值范围对应的频率点为呼吸频率范围。

Images