基于多尺度融合和相关性分析的全方向M型心动图优化研究

黄立勤 李 伟 林 强

（福州大学物理与信息工程学院，福州 350108）

引言

LEJ-2型全方向M型心动图系统是福州大学生物医学工程研究所研制开发的产品［1］，该系统是从超声心动图序列图像的研究入手，得到了从一般视频输出序列图像中提取动态信息的一种方法，即全方向灰度（位置）-时间波形图法。该方法应用于心脏超声图像中，称为全方向M型心动图。它能在心脏室壁的任意结构、任意位置的某一条方向线上重建灰度（位置）—时间波形图即组织运动曲线，通过该运动曲线可以体现心脏结构某个位置的运动情况。对组织运动曲线图进行一阶微分，可以得到组织瞬时速度，二阶微分则可以得到瞬时组织加速度。在此基础上作进一步研究，挖掘出速度场、加速度场、功率谱等直接和心脏结构机能和组织机体弹性素质相关的重要信息，并让这些信息可以以心电为时间基准互相同步地与时相相比较显示出来。这些数据对于心血管疾病和血液动力学的研究有着非常重要意义［2］。

由于全方向M型心动图的信息间接来源于二维超声序列图像，图像质量必然受超声图像质量的影响。而超声图像本质上具有模糊性和不均匀性特点，由此产生的图像质量严重阻碍了全方向M型心动图运动曲线的准确提取，如图1所示。如何去除噪声、其他组织干扰等造成的虚假边缘，提取出准确的组织运动曲线信息是全方向M型心动图动态信息分析的基础和前提，它的准确与否将严重影响到动态信息检测的精度。

图1 全方向M型心动图Fig.1 Onmi-directional M-mode Cardiography

几种广泛使用的经典的边缘检测方法有Sobel、Roberts和Laplacian等，这些算法的核心思想是假设边缘点对应于原始图像灰度级梯度的局部极值点。但是，全方向M型心动图含有较大噪声，这些算法对噪声非常敏感，常常会把噪声当作边缘点检测出来，而真正的边缘由于噪声的干扰可能被漏检。目前全方向M型心动图系统所设计使用检测算法是针对全方向M型心动图的构造特点，层进式地在有效搜索范围内运行线状搜索模板［3］，实现目标运动曲线的检测。但该方法对于目标运动曲线的检测中的虚假的边缘点会误检，仍需要较多的人工干预。

上述检测算子对运动曲线检测都在固定的尺度空间进行。固定尺度的边缘检测算子对定位精度和噪声抑制等性能指标难以兼顾。小尺度算子虽然对噪声敏感，但是有利于边缘定位；大尺度算子对噪声的抑制能力强，但对于边缘定位精度较差，有时会丢失某些局部的细节［4］。Marr通过对神经生理学和心理物理学的分析，指出人的视觉前期处理中有多个分辨率的边缘算子在对图像作卷积，如果各边缘检测算子输出进行组合就可以提高定位精度，从而减少噪声干扰［5］。由于小波算子具有良好的时频局域化特性，特别适合于多尺度分析，因此考虑利用小波来构造多尺度边缘检测算子，通过对不同尺度下的多尺度运动曲线进行融合，从而实现运动曲线的检测。

全方向M型心动图是心脏二维超声序列图像中取样线上各像素点因心脏结构在血液动力与心肌弹力的共同作用下引起运动产生的灰度（位置）随时间变化扫描展开而重建出的波型图，由于心脏是人体器官，它的运动前后具有时间、位置、能量相关性的特点。因此我们还可以分析这些相关性信息来结合前面经过多尺度融合得到的数据，进一步限制运动曲线的搜索条件和搜索范围以获得更好的目标运动曲线。

1 算法描述

根据以上思路，处理方法可以分为3个步骤。步骤1：首先通过构建小波函数，在不同的固定尺度空间下进行运动曲线检测。

步骤2：对不同尺度下的运动曲线进行融合。

步骤3：结合心脏运动的相关性信息生成最后需要的组织运动曲线。

1.1 固定尺度空间的运动曲线检测

对任意的二维图像函数f（x，y）∈L2（R2）进行小波变换，它具有两个方向分量

写成矢量形式为

式中，fk（x，y）是f（x，y）被θk（x，y）平滑后的图像。检测时把连续小波变换改成二进制形式，即取k=2j（j∈Z）。式（4）和式（5）分别代表了图像f（x，y）沿需要x和y方向上的偏导，因此二维小波变换矢量代表的就是梯度［7］。

小波变换在尺度k=2j上梯度矢量的模以及梯度矢量与水平方向的夹角可以分别计算

平滑后的图像f*θ2j（x，y）的突变点对应的是局部极大值的点，所以沿梯度方向检测矢量模的极大值点可以获得全方向M型心动图的运动曲线（边缘）。图像在（x，y）点处的边缘强度可以用式（7）来表示，图像在（x，y）点处的方向矢量可以用式（8）来表示。具体步骤如下：

步骤1：要确定运动曲线（边缘）必须删除在局部范围里面梯度不是最大的点，只保留局部范围梯度最大的点［8］。首先对点（x，y）在梯度图根据式（8）的方向归并到图2所示的4个方向，然后与该方向上的相邻点的梯度比较，如果不是则表明不是边缘点。经过处理得到初步边缘点图fb1（x，y）。

图2 方向图Fig.2 Direction figure

步骤2：对初步边缘点图fb1（x，y）进行进一步分析。首先设定两个上下限阈值Ta和Tb。这两个阈值的设置可以有不同的算法，在本系统实验测试中，首先确定Tb的值，Tb的取值为初步边缘点图fb1（x，y）中梯度直方图中频率最大的梯度，在实验中发现Tb≈2Ta时，实验效果较好。若边缘点梯度大于Tb认为该点一定是边缘点，可以得到确定边缘点图fb2（x，y）；若边缘点梯度小于Tb同时大于Ta，则认为该点可能是边缘点，可以得到可能边缘点图fb3（x，y）。fb2（x，y）的点使用较高的阈值计算出来，有比较多的间断，因此通过在fb2（x，y）曲线边缘处的点在fb3（x，y）的8个相邻像素寻找可以邻接的点，如果存在邻接点，则把该点添加到fb2（x，y），一直不断进行循环计算，最后得到fb2（x，y）是在一个固定尺度下的全方向M型心动图的运动曲线图。

在小尺度下检测获得的运动曲线比较准确，但是有较多的噪声干绕；而在大尺度空间下检测运动曲线的噪声得到抑制，但许多边缘信息被平滑掉了，而且得到的是位移后的结果。因此可以在不同的尺度空间下进行不同的阈值选择。对于较大尺度空间，阈值可以选高一些用来减少噪声产生的影响；在较小尺度空间，阈值选小一些尽可能地得到完整边缘信息。在较大尺度空间下获得的边缘不完整的范围较大。

1.2 不同尺度下的运动曲线融合

在不同尺度空间下检测得到一系列不同尺度下的运动曲线。因为在相邻尺度下的运行曲线最接近，因此运动曲线融合是对相邻的尺度下的运动曲线进行，并且按照尺度大小顺序逐步进行融合。在选择尺度时是先设定最大尺度和最小尺度，再选择尺度的间隔。虽然尺度越小检测的运动曲线越接近真实，但噪声较大，如果噪声淹没了运动曲线，则没有意义。同样尺度越大，越会抑制噪声，但如果边缘失真严重也没有意义［9］。由于不同尺度下的检测算子对全方向M型心动图的运动曲线同一位置的响应并不相同，因此融合不能进行简单的叠加。本方案利用了不同尺度下的运动曲线在位置、强度和方向上的联系分析，依次采用传递、继承、生长来处理运动曲线的融合。

设k1，k2是两个相邻尺度，k1＞k2，定义Fk1，k2（i，j）是尺度k1的3像素×3像素邻域中的像素是尺度k2上局部模极大值点（i，j）的关联域。定义集合Hk1是尺度k1上的局部模极大值点。则相关性可以表示为

式（9）可以用来表示在k2尺度上的点（i，j）与k1尺度的相关性。其中φk1（i，j）和φk2（i，j）代表尺度k1，k2上极大值点（i，j）的梯度方向，α表示它们方向差所设定的阈值。

1）运动曲线的传递

由式（9）我们可以进一步得到由k1尺度传递到k2尺度运动曲线。

尺度k2上的极大值点（i，j）如果与k1相关，则认为点（i，j）由k1传承下来。Tk1，k2代表的是在k1和k2不同尺度空间反映全方向M性心动图中同一的运动曲线。

2）运动曲线的继承

由式（9）还可以得到由k1尺度继承到k2尺度运动曲线。

Pk1，k2保留了k1尺度空间上的不与k2尺度存在相关性的点，代表的是k1到k2的继承关系。

定义Dk1，k2和φk1，k2为通过上述两个步骤对k1和k2尺度进行融合后的边缘增强图和梯度方向图。

3）运动曲线的生长

前面介绍的运动曲线的传递仅在3像素×3像素邻域中小窗口内进行，这使得运动曲线信息无法传递到较远的地方。如果增大了窗口尺寸，有可能有相关性的局部模极大值并不对应于同一运动曲线［10］。另外，噪声在大尺度空间存在较少，即使大尺度邻域内存在噪声，根据算法顺序，优先采用运动曲线的传递，其次才是运动曲线的继承。因此算法是优先采用传递的方法选择小尺度空间上正确曲线上的点，不会采用继承的方法得到噪声点。得到运动曲线的传递和继承后，还需要进行运动曲线的生长。本方案采用小窗口迭代的方法实现运动曲线的生长，具体方法如下：

不同尺度下的运动曲线融合的整个过程是在从大尺度下运动曲线附近寻找小尺度下运动曲线的相关性，若存在相关性，则用小尺度下运动曲线代替大尺度下运动曲线即完成传递；否则大尺度下运动曲线保留下来即继承；然后在小尺度下完成运动曲线的迭代扩展。因为扩展是在3像素×3像素的邻域内逐步进行，而且对方向有限制，所以不容易扩展到噪声里面。整个流程如图3所示，图3中尺度k0＞k1＞k2＞…＞kn-1＞kn。

图3 运动曲线融合流程Fig.3 Flow chart of motion curves fusion

为去除各种干扰对全方向M型运动曲线融合结果的影响，本方案按照8邻域细化算法对融合得到的运动曲线进行细化处理［11］，对边缘轮廓进行8邻域跟踪搜索，判断边缘轮廓是否长度大于阈值。如果边缘轮廓长度小于阈值，则认为是孤立点噪声造成的假边缘，在图像中将其剔除掉。由此可以看出，该方法可以很好地消除超声图像的斑点噪声影响，且预处理算法容易丢失精度信息，因此本实验对要分析的全方向M型心动图不进行降噪预处理。

1.3 全方向M型心动图基于相关性分析的运动曲线检测

全方向M型心动图系统经过了前面所阐述的多个尺度下进行融合检测后，可能包含需要的和不需要的运动曲线，并且仍然可能出现断点情况。根据全方向M型心动图的具体特点，对其进行进一步的相关性分析，才能最终获得符合要求的运动曲线图。

1）全方向M型心动图信息相关性

全方向M型心动图是心脏二维超声序列图像中取样线上各像素点由于心脏结构在血液动力与心肌弹力的共同作用下，引起运动产生的灰度（位置）随时间变化扫描展开而重建出的波型图如图4所示。全方向M型心动图相邻两列像素代表二维超声序列图像相邻两帧取样线（投射）部位心脏结构的灰度像素点集合。由于心脏结构运动具有速度有限性，一般认为最大峰值速度不超过15cm/s，而相邻两帧二维超声图像的时间间隔一般在20～40ms，因此，反映在全方向M型心动图的列像素集合之间的灰度信息就具有非常丰富的相关性。具体体现为：（1）时间（水平坐标）相关性 各列像素集合代表不同时刻取样线上的灰度信息，且从左至右对应时间上的连续性；（2）位置（垂直坐标）相关性 心脏结构运动速度的有限性与帧时间间隔决定了取样线上所代表的心脏结构的位置是密切相关的，反映在运动曲线的检测上即相邻列运动曲线点不可能出现幅度跳变；（3）灰度相关性又称为能量相关性 待检测的运动曲线在一定程度上近似代表心脏结构内同一质点的运动状态，在未发生回波失落的情况下，反映在运动曲线的检测上即相邻列运动曲线点不可能出现灰度跳变［12］。

图4 二维超声序列图像与全方向M型心动图的联系Fig.4 The relationship between two-dimensional ultrasound image sequences and onmi-directional M-mode cardiography

全方向M型心动图的构造特点及其信息相关性的存在，为检测运动曲线的搜索条件及其搜索范围的设计提供了依据。

2）基于相关性的运动曲线检测算法

心脏在幅度上的运动具有有限性，即在单位时间内（按抽样间隔算）超声图像在相邻帧之间心脏各结构的运动幅度是有限的，而且是和前一帧该结构的运动位置密切相关的，这一重要信息反映在边缘的提取上即为相邻时间点的边界不可能有大幅度跳变，表示如下：

（1）相邻时间点的边界位置是相关的。

（2）边界搜索的区域是有限的。

因此，在边缘检测的算法中也应用了这一特点：

（1）找到当前时间点的边界后，下一时间点对应的边界的搜索范围在该点的附近（向上和向下若干个点）进行。

（2）下一时间点边界搜索的范围是动态改变的，即首先在最小的区域（例如向上3个像素点和向下3个像素点）进行搜索，如果搜索不到满足条件的边界点，再逐级扩大搜索区域，进行下一级搜索。每一级搜索区域垂直方向上的步进值为3。

（3）最大搜索范围是有限的，这个限度的设置也很重要。如果设置太大可能会找到一些干扰点，而且也是系统资源的浪费；设置太小，可能找不到满足条件的边界点。全方向M型心动图是125帧/5s的序列图像进行分析得到的，帧间隔为40ms，估计心脏结构的最大速度为30cm/s，一般心脏结构的速度都小于这个值。根据公式：最大运动范围=最大速度×帧间隔时间，可得心脏各结构在这个时间间隔内的最大运动范围为1.2cm左右，一般心脏结构在这个时间间隔（40ms）内的运动范围均小于该距离。因此，设定最大搜索范围为上下30个像素点，这个距离相当于1～1.5cm。实践验证，这个最大搜索边界设置得比较合理。

（4）如果在这个区域内仍然找不到满足条件的点，即认为该处边界是断点或缺口（连续多个断点）。断点，是指一两个没有找到边界点的时间点，直接用其前后两个时间点的边界位置的平均值替代，即采用直线插值法。缺口处理，可在缺口位置附近依据位置相关性，运用前、后时间段的边界位置信息判断缺口边界所在的大致位置。

（5）根据全方向M型心动图系统检测上边缘或下边缘的需求，从上到下搜索或从下到上搜索获得上边缘或下边缘。

（把此后的内容按给的标题进行“填空”）

2 系统实验

本实验是在福州大学生物医学工程研究所开发的LEJ-2型全方向M型心动图系统上进行改进实验。LEJ-2型全方向M型心动图系统是用Visual C++6.0编写，运行环境是Windows操作系统，CPU Intel Pentium Dual E2180 2.0G，RAM 2G。实验分两个部分。

2.1 多尺度融合的全方向M型运动曲线检测

对一幅全方向M型心动图图像，经过尺度分别为k=1、2、3、4的4级小波变换后，运用本研究提出的多尺度融合算法得到边缘与传统算子边缘检测的结果进行比较。

2.2 结合多尺度融合和相关性分析的全方向M型心动图检测

通过两个实验例子将结合多尺度融合和相关性分析的全方向M型心动图检测实验结果，与目前已在使用的产品LEJ-2型全方向M型心动图系统中使用的基于线状模板搜索算法得到的检测结果进行比较。

3 结果

图5中（a）是一幅全方向M型心动图图像，对图像进行小波变换，（c）～（f）是经过4级小波变换后获得的高于给定门限局部模极大值点位置图，即各尺度下的边缘图，尺度分别为k=1、2、3、4。从边缘图可以看出，随着尺度的增大，噪声逐渐减少，边缘逐渐平滑。在尺度1时，噪声的影响非常大，边界比较破碎；尺度k=4为最大尺度，提取的边缘体现了原图中的主要边缘，基本不受噪声的影响，但是边缘失真比较严重，且提取的边缘不完整，对于某些连续的边缘，只检测出其中的一段，但是在小尺度空间可以较完整地检测出来。因此需要利用最大尺度空间提供的位置信息，融合各尺度的信息，合成精确的边缘.运用本研究提出的多尺度融合算法，结果如图5（b）所示。通过逐层融合。将小轮廓长度剔除原来断裂的边缘连接起来，而且边缘位置越来越贴近实际边缘位置，但是仍然可能有断点存在。

图6中（a）为Roberts算子（b）为Sobel算子（c）为Prewitt算子（d）为Kirsch算子检测出的边界。通过与图5（b）比较可以看出：对于噪声多、边缘模糊的全方向M型心动图图像，用传统算子检测出的边缘较模糊，去噪效果差，在定位精度、精确检测等方面都不如本研究中所采用的方法。与其他边缘检测算子的比较可以获得类似结果。

图5 各尺度下边缘检测结果。（a）原图；（b）融合之后图；（c）k=1；（d）k=2；（e）k=3；（f）k=4Fig.5 The results of edge detection on different scale.（a）Original image；（b）Image after fusion；（c）k=1；（d）k=2；（e）k=3；（f）k=4

图6 传统算子边缘检测结果。（a）Roberts；（b）Sobel；（c）Prewitt；（d）KirschFig.6 The results of the traditional edge detection operator.（a）Roberts；（b）Sobel；（c）Prewitt；（d）Kirsch

图7显示了通过实验1采用所提出算法与LEJ-2型全方向M型心动图系统中使用的基于线状模板搜索算法的比较结果。对比图7中（b）和（c）可以看出，结合多尺度融合和相关性分析的全方向M型心动图检测在质量较差的情况下，可以检测到的全方向M型心动图清晰的运动边界，同时滤掉了第1、2波峰周围虚假的边缘，而传统的LEJ-2全方向M型自动检测则需要人工修正。

图7 实验1检测结果与原系统检测结果比较。（a）全方向M型心动图1；（b）基于线状模板搜索LEJ-2全方向M型系统自动分析结果；（c）结合多尺度融合和相关性分析的全方向M型心动图检测自动分析结果Fig.7 Experimental results 1 compared with the original system.（a）Onmi-directional M-mode Cardiography 1；（b）LEJ-2 Omnidirectional M-mode Echocardiography System analysis result using linear template；（c）Result using multi-scale integration and correlation analysis

图8显示了通过实验2采用所提出算法与LEJ-2型全方向M型心动图系统中使用的基于线状模板搜索算法的比较结果。对比图8中（b）和（c）可以看出，图8（b）波峰上升的边缘检测到的是错误的边缘，而图8（c）检测到正确的边缘。

4 讨论与结论

多尺度融合的全方向M型运动曲线检测的实验表明由于小波变换有多尺度的特点，可以利用多尺度特性，通过细节和粗节进行逼近，强于其他经典算法；在边缘和噪声的取舍中，由于二者均为高频信号，很难用频带划分。使用小波变换的方法，使得可在大尺度下抑制噪声，小尺度下得到边缘的真实位置；而传统的和经典的边缘检测算法则在此问题上不能提供有效的解决办法。不论选用怎样的小波函数，都可以利用上述算法进行多尺度边缘融合。因此该方法可以有效抑制噪声的干扰，同时保证融合边界定位的准确性。

图8 实验2检测结果与原系统检测结果比较。（a）全方向M型心动图2；（b）基于线状模板搜索LEJ-2全方向M型系统自动分析结果；（c）结合多尺度融合和相关性分析的全方向M型心动图检测自动分析结果Fig.8 Experimental results 2 compared with the original system.（a）Onmi-directional M-mode Cardiography 2；（b）LEJ-2 Omnidirectional M-mode Echocardiography System analysis result using linear template；（c）Result using multi-scale integration and correlation analysis

结合多尺度融合和相关性分析的全方向M型心动图检测的实验表明：本算法对于全方向M型心动图的虚假边缘和错误边缘的检测效果较好。另外，检测的运动曲线比较平滑，能较好地滤除噪声干扰，而传统的LEJ-2全方向M型自动检测对于噪声比较敏感；经实验测试，计算出运动曲线的时间小于1s，可以符合系统应用实际要求。

由于全方向M型心动图的信息间接来源于二维超声序列图像，全方向M型图像质量必然受超声图像质量的影响。而超声图像本质上具有模糊性和不均匀性特点严重阻碍了边缘的准确提取。本研究设计了一种基于多尺度融合和相关性分析的全方向M型心动图检测算法。通过该算法能自动去除全方向M型心动图的虚假边缘、无关噪声等干扰，大大减轻了LEJ-2全方向M型系统的人工干预程度。本方案试验分析是在福州大学生物医学工程研究所开发的LEJ-2型全方向M型心动图系统基础上进行改进实验，并提供详细实验对比数据和结果分析。该设计实现了对国家发明专利“全方向M型心动图方法及其系统98 125713.5”系统的优化设计。全方向M型心动图系统具有完全的自主知识产权，改进后较大程度地减少了误差，提高了精确度，这是加快我国医疗卫生行业信息化程的非常价值的尝试。

［1］林强，石江宏.超声心动图的一种动态信息—全方向M型心动图［J］.仪器仪表学报，2005，26（4）：437-440.

［2］郭薇，陈斌，叶振盛.全方向M型超声心动图左室径向心肌运动速度梯度评价局部心肌功能的应用价值［J］.中华超声影像学杂志，2010，19（7）：565-568.

［3］Lin Qiang，Wu Wenji，Huang Liqin，et al.An omnidirectional M-mode echocardiography system and its clinical application［J］.Computerized Medical Imaging and Graphics.2006，30：333-338.

［4］阮松，陈松灿，王敏.采用多尺度多级组合分类器快速定位乳腺X片中的感兴趣区域［J］.中国生物医学工程学报，2008，28（5）：675-685.

［5］Marr.视觉计算理论［M］.姚国正，刘磊，汪云九，译.北京：科学出版社，1988：256-260.

［6］张玉叶，周晓东，王春歆.小波估计图像棱边分布的边缘保持规整化复原［J］.计算机辅助设计与图形学学报，2009，21（1）：125-135.

［7］Yitzhaky Y，Peli E.A method for objective edge detection evaluation and detector parameter selection［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence.2003，25（8）：1027-1033.

［8］Canny J.A computational approach to edge detection［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1986，8（6）：679-698.

［9］Konishi S，Yuile AL，Coughlan JM，et al.Statistical edge detection：Learning and evaluating edge cues［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2003，25（1）：54-57.

［10］汪荣贵，朱静，杨万挺，等.基于照度分割的局部多尺度Retinex算法［J］.电子学报，2010，38（5）：1181-1186.

［11］Bao Paul，Zhang Lei.Noise reduction for magnetic resonance image via adaptive multiscale products thresholding［J］.IEEE Transactions on Medical Imaging，2003，22（9）：1089-1099.

［12］Chalana V，Linker DT，Haynor DR，et al.A multiple active contour model for cardiac boundary detection on echocardiographic sequences［J］.IEEE Transactions on Medical Imaging，1996，15（3）：290-298.