一种血管内超声弹性成像二维多级混合位移估计方法

摘要

本发明公开了一种血管内超声弹性成像二维混合位移估计方法，本发明提出的二维混合位移估计方法，结合了数据点的多级搜索和单级跟踪的方法，体现了合理平衡计算量和计算精度的优点。在算法初级阶段，利用二维多级搜索的方法来获取轴向和横向偏移量，而避免在高级阶段对横向方向进行搜索。在高级阶段，为最大程度减小计算量，则利用质量指导和零点相位差跟踪方法来估计轴向位移。这种算法可以很好地适用于血管内超声成像的实时成像要求，以及成像精度。

说明书

技术领域

本发明涉及超声成像技术领域，具体涉及的是一种血管内超声弹性成像二维多级混合位移估计方法。

背景技术

血管内超声成像（intravascular ultrasound，IVUS）是20世纪80年代末迅速发展起来的将无创性的超声技术和有创性的心导管技术相结合的一种新的介入式超声成像技术，是目前唯一商业化用于临床检测的可以实时提供患者冠状动脉血管横截面图像的检查手段。具体方法是在心导管检查过程中，经导丝将超声探头送至心血管腔内进行探测，继而回撤超声导管，可观测到管腔大小、管壁结构的厚度和形态等，观测其收缩和舒张的功能变化，还能对病变进行定性分析甚至是精确测量分析。不仅能反映血管内腔的变化，同时也能反映含有斑块的血管横断面结构以及斑块的性状等等，提供了冠状动脉造影技术所不能提供的重要信息。目前IVUS成像技术已可对冠状动脉甚至更细小血管进行血管内成像，可定性、定量地提供动脉壁微结构灰度图像，对于冠状动脉粥样硬化与狭窄等心血管疾病的诊断与治疗具有重要意义。临床应用经验表明该技术具有直观、准确等优点，被认为是诊断冠心病新的“金标准”。

目前应用于IVUS系统的主要技术有虚拟组织学成像（Virtual Histology IVUS,VH-IVUS）和复合散斑成像(Integrated Backscatter IVUS,IB-IVUS)等，均是基于超声回波强度差别来区分富含脂质的低回声斑块即软斑块和富含纤维成分的高回声斑块即硬斑块两种，但对于脂类或者半脂类的斑块不能够很好的进行定性，而且存在导管位置偏离轴心和运动不均匀等伪影影响。

1991年，Ophir和其研究小组初次提出了超声成像技术——弹性成像（Elastography）以来，逐步成为当前国际医学超声领域最具价值的热点课题之一。弹性成像方法不仅可以提供组织硬度的图像，也就是关于病变的组织特征的信息，以此来定位病变，而且可以根据不同组织间弹性系数不同，在受到外力压迫后组织发生变形的程度不同，将受压前后回声信号移动幅度的变化转化为实时彩色图像。将超声弹性成像方法运用到血管内超声成像技术上，可以分辨高低应变区域，从而识别冠状动脉血管内不同斑块的相对生物力学特征，辅助判断斑块破裂难易程度，进一步的提高心血管疾病的检测手段。

自弹性成像提出以来，就有许多方法不断发展起来，总体来说可以分为两大类型：一是一维方法，主要包括两种：1.基于梯度运算的应变估计算法，即先对组织的位移进行估计，然后对位移估计做差分处理，得到组织的应变分布，如时域互相关算法、相位检测法、过零点跟踪法等；2.直接的应变估计算法，不以组织位移为中间估计值，直接得到组织的应变分布，如自适应伸展算法、谱相关算法等。具体方法是利用超声探头，采集组织压缩前和压缩后的探头每一条扫描线的射频信号，再分别利用上述提到的算法对每一条扫描线进行对应组织的应变分布。除分别具有自身的优缺点外，一维方法仅考虑换能器每一条扫描线接收到的组织轴向方向运动的数据，而没有考虑组织在横向的运动，因此并没有包括组织运动的全部信息。而现有的二维方法，优点是可以同时测量运动物体的轴向和横向运动分量，不像一维运动估计算法那样受限于测量角度等问题。但是现有的二维算法中，如图像相减算法（Hein et al., 1993），就是用一帧图像与另一帧图像作相减运算，如果两幅图像中相同位置处的像素的亮度没有变化，则该像素在相减运动结果图像中亮度为零，表示组织没有发生运动；而亮度有变化的区域，在相减运算结果图像中具有不同亮度，表示组织在这个区域有运动。这种方法对存在较大差别的两幅图像的检测效果还可以，但亮度级数太少，对比度不高。优点是具有实时性，缺点是不能得到运动的定量信息，检测不出运动方向。或者二维区域匹配算法，以及基本相关算法，本质上是利用图像序列之间的一种对应关系，找到目标区域的变化，但这种算法计算量非常大，十分耗时，不能很好的应用于实时的血管内超声成像系统，不能满足其实时性要求。

发明内容

本发明的目的在于克服现有技术存在的以上问题，提供一种血管内超声弹性成像二维多级混合位移估计方法，很好地解决平衡计算时效和精度的问题。

为实现上述技术目的，达到上述技术效果，本发明通过以下技术方案实现：

    一种血管内超声弹性成像二维混合位移估计方法，包括以下步骤：

    步骤1）利用商用超声换能器得到待测生物组织的压缩前的一帧二维射频信号；

    步骤2）手持所述超声换能器沿着探头纵向对所述生物组织施加一个微小的挤压，得到压缩后的一帧二维射频信号；

    步骤3）分别对从步骤1和步骤2中得到的组织压缩前、后的射频信号进行计算，得到基带解析数据，具体步骤为先对射频数据进行下采样减少数据点，为避免下采样时组织位移信息丢失，压缩前、后的信号采样频率应尽可能高，而下采样时，在不违背麦奎斯特采样定律的前提下，下采样频率也应该尽可能高，再进行希尔伯特变换，乘以exp^(-jw₀^t),其中 w₀为超声换能器的中心频率，则得到基带解析信号pre(x,y)和post(x,y), 其中x表示沿探头轴向方向采样点，y表示沿探头横向方向采样点；

    步骤4）对步骤3中得到的压缩前、后基带解析信号进行信号处理，分别得到压缩前、后信号的幅值a₁(x,y)和a₂(x,y)；

步骤5）确定第一级窗口的数量及大小，可取9个，每行每列各3个；

    （1）压缩前的目标区域第一级窗口大小为轴向或横向采样点的五分之一；

    （2）压缩后的搜索区域第一级窗口大小略大于轴向或横向采样点的五分之一，可取四分之一；

步骤6）利用步骤5中确定的参数，对步骤4中得到的压缩前、后信号的幅值a₁(x,y)和a₂(x,y)进行第一级窗口搜索，搜索时对压缩前、后帧进行相关系数最大值匹配，相关系数公式如式（1）所示；

（1）在最低分辨率下，压缩前和压缩后帧中的每一个窗口中，每四个采样点作为一个对象点，图2中的上一图中的方形点为压缩前窗口的中心，对齐到压缩后的搜索区域进行相关系数计算，找到最佳匹配的矩阵形区域并做记录；

（2）将分辨率提高一级，每两个采样点作为处理对象点，上一步中找到的外圆方形点（图2中的下一图中的点）周围的x方向三个点作为这一步的匹配候选点，计算相关系数，找到最佳匹配点，并做记录；

（3）分辨率最高的一级，不间隔点进行搜索，将上一步方形点周围的三个点作为候选点，且计算相关系数时不再间隔任何采样点，找到最佳匹配点，完成第一级搜索，只有在最低分辨率时，候选匹配位置尺寸与搜索区域的大小相关；在下一步分辨率提高的情况下，对相关系数的要求高于计算量，候选匹配区域定为三个，得到9个独立的轴向-横向位移估计值；

                                                        (1)

其中分别是横向和轴向位移，分别是压缩前和压缩后的图像幅值，分别是的内插窗口平均值，T是窗口大小；

    步骤7）确定第二级的窗口数量及大小，其数量与算法参数无关，其中有九个为第一级9个窗口的子集，因此第一级搜索得到的轴向和横向位移估计就可以传递到第二级，第二级其余的窗口则均匀分布整个帧，第二级运动的方法和第一级相同，但间隔因子初始即取每两个采样点作为一个对象点，压缩前的目标区域窗口大小为轴向或横向采样点的十五分之一，而压缩后的搜索区域窗口大小则略大于目标区域，可取轴向或横向采样点的十分之一；

    步骤8）重复步骤6中的（1）、（2）、（3）步，但第二级的匹配过程由动态跟随来完成，首先第二级的9个窗口由第一级的输出确定一个初始点p1(x0,y0)，其余的窗口初始值则由第一级的结果大致定向，定位于每一个窗口中心，将p1存储于点集S中，以及一个队列L集中，第二步，计算p1的四个邻点，并将其加入到S集中， S2 = {p1,p2, p3, p4, p5}， p2 = (x0 +1, y0)， p3 = (x0 ?1, y0)，p4 = (x0, y0 +1)，p5 = (x0, y0 ?1)；此时S集中有四个点待处理，根据式（1），选出下一个处理点P2：P2 = max (C(p1),C(p2), C(p3), C(p4),C(p5))，其中C即为互相关系数最大，此时用P2将队列L集更新，紧接着搜索当前点P2的邻点，如此经过连续不断的递进搜索，L集队列中永远是相关系数最大值，直到判断条件S集不再有更新，则该区域搜索完毕，得到更普遍分布的第二级轴向-横向位移估计值；

    步骤9）步骤8完成后，利用Savitzky-Golay滤波方法，对得到的横向估计进行平滑；

    步骤10）第三级搜索之前，对第二级得到的横向位移估计值利用双线性插值法，得到第三级的横向位移估计，因为横向位移对噪声敏感，窗口越小，这种影响越大，因此，第三级不再对横向位移估计做更进一步的精确；

    步骤11）确定第三级窗口数量及大小，第三级需要完成整幅帧的搜索，因此数量由窗口大小确定，窗口大小在轴向方向可取四十分之一，第二级获得的轴向位移估计作为第三级窗口的子集，进行结果传递；

    步骤12）第三级因为窗口较小，此时默认每个窗口的弹性模量为常数，位移可以近似看做深度的线性函数，窗口的位移可以假定参考窗口中心点的位移，但由于散斑和背向散射可能会引起小窗口中心位移估计局部产生巨大变差，因此利用对数压缩方法，对基带信号进行处理，以此来减小幅值的变化，如式2和式3所示；

       (2)

             (3)

其中,c表示压缩因子。    步骤13）因为第三级窗口的数量较多，不再利用互相关系数作为匹配准则，而是改用压缩前、后信号相关性最大时，其相位差趋近于零点的特点，进行牛顿迭代快速收敛的方法，如式4所示，其中arg(pre_log(x,y))和arg(post_log(x,y))为步骤12中得到的压缩前、后对数压缩后的相值，W(x,y)为采样点对位移估计的贡献值，即权函数，表示前一步得到的位移值，为当前所求的位移值，w₀是换能器的中心频率，而且因为信号的采样点为离散点，且经过下采样处理，因此实际中很难完全趋近于零点，可以通过限制迭代次数来达到精确位移估计的目的，同时不会进入死循环，迭代次数可以由具体的信号采样点数量及所要求的精度最终确定；

 

 

      (4)

    步骤14）步骤13中的W(x,y)，可取1，也可取步骤4得到的压缩前、后信号的幅值a₁(x,y)和a₂(x,y)对应窗口的和，如式5所示；

                (5)

    步骤15）利用步骤8中的跟踪方法结合步骤13中提到的相位差方法，将轴向位移估计传递到剩下的每一个第三级窗口，完成算法的最后一步搜索，即可得到整幅帧图的位移估计值。

    进一步的，所述步骤2中的微小的挤压控制在1%的数量级。

本发明的有益效果是:

本发明提出的二维混合位移估计方法，结合了数据点的多级搜索和单级跟踪的方法，体现了合理平衡计算量和计算精度的优点。在算法初级阶段，利用二维多级搜索的方法来获取轴向和横向偏移量，而避免在高级阶段对横向方向进行搜索。在高级阶段，为最大程度减小计算量，则利用质量指导和零点相位差跟踪方法来估计轴向位移。这种算法可以很好地适用于血管内超声成像的实时成像要求，以及成像精度。

附图说明

图1是算法总体流程示意图;

图2为第一级的搜索策略示意图图；

图3为第二级、第三级的搜索策略流程图。

具体实施方式

下面将参考附图并结合实施例，来详细说明本发明。

    一种血管内超声弹性成像二维混合位移估计方法，包括以下步骤：

    步骤1）利用商用超声换能器得到待测生物组织的压缩前的一帧二维射频信号；

    步骤2）手持该超声换能器沿着探头纵向对该组织施加一个微小的挤压（一般控制在1%的数量级），得到压缩后的一帧二维射频信号；

步骤3）分别对从步骤1和步骤2中得到的组织压缩前、后的射频信号进行计算，得到基带解析数据，具体步骤为先对射频数据进行下采样减少数据点，为避免下采样时组织位移信息丢失，压缩前、后的信号采样频率应尽可能高，而下采样时，在不违背麦奎斯特采样定律的前提下，下采样频率也应该尽可能高，再进行希尔伯特变换，再乘以exp^(-jw₀^t),其中 w₀为超声换能器的中心频率，则得到基带解析信号pre(x,y)和post(x,y)，其中x表示沿探头轴向方向采样点，y表示沿探头横向方向采样点；

    步骤4）对步骤3中得到的压缩前、后基带解析信号进行信号处理，分别得到压缩前、后信号的幅值a₁(x,y)和a₂(x,y)；

    步骤5）确定第一级窗口的数量及大小，可取9个，每行每列各3个，（1）压缩前的目标区域第一级窗口大小为轴向或横向采样点的五分之一；（2）压缩后的搜索区域第一级窗口大小略大于轴向或横向采样点的五分之一，可取四分之一；

    步骤6）利用步骤5中确定的参数，对步骤4中得到的压缩前、后信号的幅值a₁(x,y)和a₂(x,y)进行第一级窗口搜索，搜索时对压缩前、后帧进行相关系数最大值匹配，相关系数公式如式（1）所示，搜索策略可参照图1，（1）在最低分辨率下，压缩前和压缩后帧中的每一个窗口中，每四个采样点（间隔3个点）作为一个对象点，图2的上一图中的方形点为压缩前窗口的中心，对齐到压缩后的搜索区域（图2的下一图）进行相关系数计算，找到最佳匹配的矩阵形区域（图2的下一图中外圆方形点为其中心点）并做记录；（2）将分辨率提高一级，每两个采样点作为处理对象点，上一步中找到的外圆方形点周围的x方向（轴向）三个点作为这一步的匹配候选点，计算相关系数，找到最佳匹配点（图2中的下二图中外圆方形点），并做记录；（3）分辨率最高的一级，不间隔点进行搜索，将上一步方形点周围的三个点作为候选点，且计算相关系数时不再间隔任何采样点，找到最佳匹配点（图2中的下三图中外圆方形点），完成第一级搜索。只有在最低分辨率时，候选匹配位置尺寸与搜索区域的大小相关；在下一步分辨率提高的情况下，对相关系数的要求高于计算量，候选匹配区域定为三个，得到9个独立的轴向-横向位移估计值；

        (1)

其中分别是横向和轴向位移，分别是压缩前和压缩后的图像幅值，分别是的内插窗口平均值，T是窗口大小；

    步骤7）确定第二级的窗口数量及大小，其数量与算法参数无关，其中有九个为第一级9个窗口的子集，因此第一级搜索得到的轴向和横向位移估计就可以传递到第二级，第二级其余的窗口则均匀分布整个帧，第二级运动的方法和第一级相同，但间隔因子初始即取每两个采样点（间隔一个采样点）作为一个对象点，压缩前的目标区域窗口大小为轴向或横向采样点的十五分之一，而压缩后的搜索区域窗口大小则略大于目标区域，可取轴向或横向采样点的十分之一；

    步骤8）重复步骤6中的（1）、（2）、（3）步，但第二级的匹配过程由动态跟随来完成，首先第二级的9个窗口由第一级的输出确定一个初始点p1(x0,y0)，其余的窗口初始值则由第一级的结果大致定向，定位于每一个窗口中心，将p1存储于点集S中，以及一个队列L集中，第二步，计算p1的四个邻点，并将其加入到S集中， S2 = {p1,p2, p3, p4, p5}， p2 = (x0 +1, y0)， p3 = (x0 ?1, y0)，p4 = (x0, y0 +1)，p5 = (x0, y0 ?1)；此时S集中有四个点待处理，根据式（1），选出下一个处理点P2：P2 = max (C(p1),C(p2), C(p3), C(p4),C(p5))，其中C即为互相关系数最大，此时用P2将队列L集更新，紧接着搜索当前点P2的邻点，如此经过连续不断的递进搜索，L集队列中永远是相关系数最大值，直到判断条件S集不再有更新，则该区域搜索完毕，得到更普遍分布的第二级轴向-横向位移估计值，具体流程可参考图3；

    步骤9）步骤8完成后，利用Savitzky-Golay滤波方法，对得到的横向估计进行平滑；

    步骤10）第三级搜索之前，对第二级得到的横向位移估计值利用双线性插值法，得到第三级的横向位移估计，因为横向位移对噪声敏感，窗口越小，这种影响越大，因此，第三级不再对横向位移估计做更进一步的精确；

    步骤11）确定第三级窗口数量及大小，第三级需要完成整幅帧的搜索，因此数量由窗口大小确定，窗口大小在轴向方向可取四十分之一，第二级获得的轴向位移估计作为第三级窗口的子集，进行结果传递；

    步骤12）第三级因为窗口较小，此时默认每个窗口的弹性模量为常数，位移可以近似看做深度的线性函数，窗口的位移可以假定参考窗口中心点的位移，但由于散斑和背向散射可能会引起小窗口中心位移估计局部产生巨大变差，因此利用对数压缩方法，对基带信号进行处理，以此来减小幅值的变化，如式2和式3所示；

       (2)

             (3)

其中,c表示压缩因子。

     步骤13）因为第三级窗口的数量较多，不再利用互相关系数作为匹配准则，而是改用压缩前、后信号相关性最大时，其相位差趋近于零点的特点，进行牛顿迭代快速收敛的方法，如式4所示，其中arg(pre_log(x,y))和arg(post_log(x,y))为步骤12中得到的压缩前、后对数压缩后的相值，W(x,y)为采样点对位移估计的贡献值，即权函数，表示前一步得到的位移值，为当前所求的位移值，w₀是换能器的中心频率，而且因为信号的采样点为离散点，且经过下采样处理，因此实际中很难完全趋近于零点，可以通过限制迭代次数来达到精确位移估计的目的，同时不会进入死循环，迭代次数可以由具体的信号采样点数量及所要求的精度最终确定；

 

 

      (4)

    步骤14）步骤13中的W(x,y)，可取1，也可取步骤4得到的压缩前、后信号的幅值a₁(x,y)和a₂(x,y)对应窗口的和，如式5所示；

                (5)

    步骤15）利用步骤8中的跟踪方法结合步骤13中提到的相位差方法，将轴向位移估计传递到剩下的每一个第三级窗口，完成算法的最后一步搜索，即可得到整幅帧图的位移估计值。