超声剪切波弹性摄影中的自相关引导的互相关的制作方法

本发明涉及跟踪剪切波诱发的运动，并且更具体而言，涉及在运动跟踪中的对相关的使用

背景技术：

超声剪切波弹性摄影是能够测量软组织的诸如剪切模量和剪切黏度的机械性质的新的医学模态。因为组织弹性与病理有关，因此弹性摄影能够提供额外的临床信息以增加诊断置信度。一个范例是通过测量肝硬度的肝纤维化分级。

在基于声辐射力的超声剪切波弹性摄影中，专用脉冲序列包括一个或多个长推动脉冲(通常每个长数百微秒)和和一系列交错跟踪脉冲。由于声辐射力的效应，推动脉冲使得聚焦区中的组织移动远离探头表面，同时建立在垂直于推动波束的方向上传播远离聚焦区域的剪切波。针对聚焦深度处的沿剪切波路径的每个横向位置，由剪切波诱发的组织运动将主要在与推动波束相同的方向上。同时，针对给定位置，跟踪脉冲能够监测这样的动态响应并且将其转化为表示作为时间的函数的组织移动的幅度的位置特异性位移波形。这样的波形能够在沿剪切波传播路径的多个位置处来计算。所述波形能够被输入到用于计算传播的速度的过程中。例如，傅里叶变换能够在这些剪切波形上执行以估计剪切波相位速度。备选地，剪切波幅度峰-峰空间和时间计算也能够产生剪切波传播速度。作为结果，能够确定组织机械性质的绝对值。具体而言并且通过范例，剪切波在组织内部传播的速度由以下项来支配：剪切模量、剪切黏度、组织密度以及通过某些机械模型的剪切波频率。组织越硬，波移动就越快。刚度的量度然后能够被用在例如肝纤维化分级中。

准确、可靠并且有效的运动跟踪是在以任何形式的超声剪切波弹性摄影对组织性质的最终估计中的目标。总体而言，在超声成像中有两种主要的运动跟踪技术：1)通过自相关的相位移位、2)通过互相关的时间移位(或者其他备选，诸如和绝对差(SAD))。

在基于自相关的方法中，感兴趣区结构的位移诱发在由移动的介质背散射的相继的高频超声回波上的相位移位。自相关已经被实施在用于实时彩色流成像的多数商用超声系统上。当超声波传播通过软组织时，谱经历由于频率相关的衰减的下频移。因此，假设在基于多普勒的方法中的恒定中心频率将导致位移估计虚假。已知在快速时间域(即，在超声脉冲的轴向方向上)使用额外的自相关，以估计局部中心频率，并且随后改进位移估计准确性。参见T.Loupas,J.T.Powers和R.W.Gill在IEEE Trans.Ultrason.Ferroelectr.Freq.Control 42.4.，第672–688页，(1995)上的“An Axial Velocity Estimator for Ultrasound Blood Flow Imaging,Based on a Full Evaluation of the Doppler Equation by Means of Two-Dimensional Autocorrelation Approach”一文。

通过互相关的时间移位估计通过使用射频(RF)数据或者从RF数据传达的符合信号从一个脉冲到另一脉冲进行互相关来估计时间延迟。

技术实现要素：

下文在本文中提出的内容解决了上述问题中的一个或多个。

用于位移估计的两种标准方法是自相关和互相关。自相关比互相关更有计算效率。然而，其受混叠影响--最大位移限制是四分之一波长(针对具有中心频率3Hz的标准肝脏图像探头，限制近似为120μm)。

在常规的、即无引导的互相关方法中，混叠限制被消除。然而，其耗费更多计算时间，所述计算时间与预定义搜索区成比例地相关。在大的区上的搜索使得标准互相关在计算上是昂贵的。然而，如果搜索区被设置太小，则互相关函数的真实峰将被遗漏；如果其太大，则由于在信号中包含的噪声而将出现跳峰。在这样的情况下，互相关能够产生超过真实峰的虚假相关峰。虚假峰将在位移估计中引入不连续和误差。遗憾的是，在实践中，所需要的搜索区的尺寸将是变化的。在临床应用中，超声剪切波弹性摄影诱发的软组织运动大约为十微米(10-μm)。上限主要由美国食品药品监督管理局(FDA)声学输出安全限制所约束。针对一个剪切波推动脉冲，通常在10-20毫秒(ms)的时段内由超声跟踪脉冲来观察所诱发的剪切波运动。如果使用多个推动脉冲，则观察能够增加到50-100ms。尽管通常建议患者在剪切波弹性摄影期间屏息几秒钟，但是在50-100ms的时间段期间，患者运动(心脏运动和呼吸)和背景运动有时能超过四分之一波长并且达到比真实剪切波运动高一个量级以遮蔽真实剪切波信号。在这样的情况下，如果单单利用常规互相关，则必须使用非常大的搜索尺寸以便定位最大互相关函数。这导致更长的计算时间并且潜在地导致虚假跳峰。如果单单使用自相关，则将可能发生相位混叠。

在本文中提出的两步运动跟踪方法对背景运动更为鲁棒。

为了对软组织中的剪切波传播的更准确和有效的跟踪，该新的方法使用自相关和互相关两者。这例如对于超声剪切波弹性摄影的体内应用而言是有用的。两步方法可表征为自相关引导的互相关。在第一步骤中，相位根据复信号的自相关来导出并且被适当地被展开。然后，使用恒定中心频率根据相位来计算初始位移。在第二步骤中，所述初始位移被量化并且被用作开始点以向互相关馈送减小的搜索窗，从而获得具有高精度的最终位移。通过将自相关与互相关整合，所提出的方法非常适于实时应用。在即时建议中的通过互相关的计算比标准互相关快得多。其还减小了位移估计和偏置的变化，而不需要RF信号的平均频率的明确估计。

该新的方法相对于常规自相关的性能改进包括：1)混叠被移除；并且2)不需要考虑中心频率深度衰减和频率随机波动的效应。存在后者的益处是因为使用基于时间移位的互相关的位移估计不受频率改变的影响。此外，该新的方法相对于常规互相关的性能改进包括：1)显著降低的搜索尺寸，以克服跳峰(虚假最大值或最小值)，从而得到改进的估计准确性；以及2)更快速的计算。

在本文中所提出的内容以这样的方式组合自相关和互相关估计器：利用每个的优点，以改进总体运动跟踪性能。例如，由于在步骤一中通过相位展开未移除的任何混叠造成的误差通过接下来的互相关来校正或减轻。作为另一范例，互相关的相对高的计算负担通过为互相关提供由自相关导出的开始点显著缩短的搜索来缓解。相比于任一个体协议，该混合协议能够获得针对剪切波弹性摄影的更准确并且鲁棒的位移估计。

根据所提出的技术的一方面，一种超声运动估计算法，包括：发出在剪切波的传播方向上彼此间隔开的多个超声脉冲，以跟踪由剪切波引起的轴向运动。所述剪切波已经由轴向方向的超声或者机械诱发的推动脉冲诱发。基于跟踪脉冲回波，所述方法使用自相关来估计轴向位移。所述方法然后使用所述估计作为在基于时间域的运动跟踪算法中的开始点，以用于修改所述估计，从而产生经修改的位移。所述修改在运动跟踪准确性方面能够构成在估计上的改进。所述发出可以对应地从多个声学窗(多个超声成像探头分别经由所述窗进行成像)发生。所述自相关和算法能够具体地操作在经由用于跟踪由推动(所述推动是单个推动)诱发的波所引起的运动中所使用的脉冲采集的成像上。所述算法可以涉及经受图像匹配准则影响的递增地增加的在搜索区上的互相关。针对流程的采集能够扩展在多个推动上。

在剪切波弹性摄影中的新颖的、互相关引导的基于时间域的运动跟踪的细节在下文借助如下附图进一步阐述，所述附图未按比例绘制。

附图说明

图1是根据本发明的示范性超声运动估计设备和无噪声剪切波的概念性描绘的示意图；

图2是提供与图1的设备的操作有关的概念的范例的概念图；并且

图3是说明根据本发明的用于在剪切波弹性摄影中的互相关引导的基于时间域的运动跟踪的可能操作的一组流程图。

具体实施方式

图1通过说明性和非限制性范例描绘了超声运动估计设备100，超声运动估计设备100能用于在自相关引导的基于时间域的运动跟踪中的剪切波弹性摄影。设备100包括图像采集电路104、运动跟踪电路108、控制器112以及数据存储设备116。图像采集电路被配置用于从各自横向位置124的间隔开120的图像采集。设备100还包括多个超声成像探头128。从各自声学窗132执行图像采集。声学窗是探头128的成像表面上与被检查的对象或患者(人或动物)接触、交换用在成像中的超声的区域。在轴向方向136上，发射超声脉冲140，并且从患者回波返回射频(RF)数据144。从被采样的多个成像深度148被动态地回波返回脉冲140，以询问患者内的介质。

尽管图1中的不同标记标注不同的探头128，但是这是出于解释目的。探头128通常可以在对介质进行成像的过程中被相同地操作。它们可以同时并行地被操作，尽管在本文中所提出的内容并不限于这样的操作。

横向方向152是由轴向方向的推动生成的剪切波的传播的方向。所述推动由被聚焦到具体成像深度的超声来生成以递送足够的声辐射力，或者可以是机械诱发的推动。

传播的剪切波主要是在轴向方向136上在该成像深度处将介质变型和位移。图1中的曲线图示出了表示由剪切波对身体组织的位移的三个波形。所述三个波形对应于远离推动脉冲的焦点0mm、5mm和10mm的横向距离。从对应于曲线图的原点的时间处开始，所述推动脉冲持续数百微秒。组织位移被跟踪35ms。出于说明目的，它们被归一化到在0mm处的峰位移。不同横向距离(即，0mm、5mm、10mm)处的剪切波诱发的位移或者变形158表明传播剪切波的位移效应。它们的波形根据组织机械性质和距剪切波原点(即，推动脉冲焦点位置)的横向传播距离而变化。

如本文上面所提到的，自相关是在用于跟踪轴向位移的两步方法中的第一步骤。在移动介质中的感兴趣结构的位移诱发在由移动介质背散射的相继的高频超声回波上的相位移位。

在慢时间中由复信号(复解析或者基带信号)的1-D自相关来估计相位移位。基带数据包括同相(I)和正交(Q)分量，其能在接收-波束形成后通过解调导出，以移除载波频率。从RF数据对复解析信号的导出也是公知的，并且在Elfataoui等人的美国专利No.7873686中有所描述。

在剪切波跟踪中，G₀(x,y)是在推动脉冲的激励之前被采集为参考信号的复合图像。通过其采集参考信号G₀(x,y)的参考脉冲是运动跟踪的初始部分。轴向方向136沿着x轴，并且横向方向152沿着y轴。每个点或者“像素”(x_p,y_q)表示在对应于x_p的成像深度148处并且在对应于y_q的横向位置124处的图像样本。

G_n(x,y)是在剪切波跟踪期间的推动后的第n图像采集。所述图像采集能够经由探头128例如并行地重复，其中，n的范围是在单个推动后的达到例如N＝30或者更多(N由跟踪脉冲重复频率和跟踪持续时间来确定)。但是，跟踪可以延伸超过第二或者甚至随后的推动，其在本文中所提出的内容的期望的范围内。

针对复基带信号或复解析信号，用于计算在像素(x,y)处的1-D相关性R_k,n的一般化的公式被表达如下：

其中，(2M+1)是如图2中看到的相关核203的在轴向方向136上的尺寸；并且R_k,n是快时间(轴向方向)中的k滞后和慢时间(跟踪时间)中的n滞后处的相关系数。G_n^* 是G_n的复共轭。

针对第n剪切波跟踪采集的在(x,y)处的位移D被初始地计算为：

其中，c是声音的速度，∠R_0,n(x,y)是快时间中的零滞后和慢时间中的n滞后处的1-D自相关系数的角，并且f_c是中心频率。因此，在下文被讨论的量化之后，D_n(x_p,y_q)充当介质的轴向位移的估计，所述估计特异于当前像素(x_p,y_q)。角或者“相位”可从如下公式导出：

并且，在代入到等式(2)中之前，经受相位展开以缓解或避免混叠。函数Re和Im分别提取复值自相关系数R_n的实部和虚部。关于相位展开，在慢时间中的展开的相位移位(即，在n与n+1之间，其中，n＝1，…，N-1)被希望为平滑而没有中断。参考图2，针对采集n的卷绕相位206跟随有采集n+1的卷绕相位208。针对采集n+2的相位210卷绕回去，以2π为模，在约束区间[π，-π]内。这导致卷绕相位212，以及在幅度上接近2π并且不表示介质中的实际轴向运动的相位不连续214。相位展开解决了这一问题，产生展开的相位210。不管相位展开，如果发生中断，则其通常指示显著外部引起的运动，即，起源于剪切波的平面。

等式(1)具有被用在针对“x”的特定值的自相关的核中的2M+1采样深度。针对任意给定采集n，从2M+1采样深度中的每个采样深度和额外的采样深度来采集样本。所述额外的采样深度允许等式(1)被运行多次，每次针对“x”的不同的值。例如，在一次运行中，要针对像素(x_p,y_q)确定D_n(x_p,y_q)；在下一运行中，要针对像素(x_p+1,y_q)确定D_n(x_p+1,y_q)，其中x_p和x_p+1在对应不同的采样深度处。

等式(1)不仅重复地被运行，每次针对其位移估计正被计算的给定像素；该组运行被重复，每次针对不同的采集n(1≤n≤N)。

此外，等式(1)的上述运行中的全部被重复，每次针对不同的横向位置。每个横向位置124对应于“y”的值，因此，y_q的横向位置124不同于y_q+1的横向位置。

在图2中，在给定横向位置124处的采样深度中的两个连续采样深度被指代为218、220。

在给定横向位置124处的给定采集n中，通过在接收窗期间针对所有采样深度，从采样深度218到采样深度220地采集数据。这是所有重做采集到采集。继而，重复的采集(即，每个是采样深度218到采样深度220)是所有重做横向位置到横向位置。

整组采集被完成一次，并且被记录在数据存储设备116中，该数据然后在步骤一和二中重复地被调出。

等式(2)中的轴向位移D(或者“D_AC”)指示针对两步方法中的第二步骤的开始点。

在等式(2)中的D_AC，如从当前像素(x_p,y_q)轴向偏离，通常将具有被设置在采样位置之间的结束点。

替代使用D_AC自身作为开始点，基于由通过箭头224在图2中所示的轴向方向中的像素间隔单元222分离的连续采样深度218、220来量化D_AC。相同的像素间隔或者另外细化的间隔(通过对RF或复信号上采样)可以被用在第二步骤中。在后者情况下，应当在上采样之前和之后使用像素间隔的比率来调节D_AC。

量化，其可以是上或者下，并且通过基于接近度通过舍入或者通过缺省，导致量化的位移[D_AC]226。

这充当身体组织的轴向位移的像素特异性(即，特异于(x_p,y_q))估计。

通过[D_AC]偏移的像素(x_p,y_q)是针对步骤二的开始点。其是精细调谐来自第一步骤的估计226的第二步骤--这是经由逐图像的匹配来完成的，跟随有峰搜索以及多项式拟合或者相位过零检测，以用于精细调谐。

由于步骤二中的逐图像的匹配是基于时间域的，因此不需要在步骤一中考虑针对采集发出的超声的中心频率衰减，并且不需要在步骤一中校正针对中心频率衰减的估计226。有利地，中心频率衰减校正开销被避免。衰减从采样深度228到采样深度230积累231。图2中的“x’s”交叉箭头表示不需要考虑在跟踪脉冲140的传播期间发生的中心频率衰减。

针对作为基于时间域的运动跟踪算法的部分的互相关，开始点234被提供为当前考虑的像素(x_p,y_q)的从在图2中被指代为“x_p”轴向位置的偏移226。偏移或者“定位”226对应于相对于被用在自相关中的相关性核203的参考图像227的移位。被移位的核236被用在互相关中。

基于时间域的互相关通过将来自一个脉冲的接收的RF回波(或者复解析或基带信号)互相关到另一脉冲以估计时间延迟。

下文在本文中立即描述二维互相关，因为步骤一和步骤二可以被用于确定二维偏移。这将为了更大的准确性，以额外的计算为代价。

在二维超声成像中，位移向量(u,v)的两个分量能够使用2-D散斑跟踪流程来估计。因此，“u”在轴向方向136上，并且“v”在横向方向152上。具体而言，在组成参考图像G₀(x,y)的初始相位敏感信号或图像(即，RF信号、复基带或者复解析信号)中的每个像素(x,y)处，近似等于一个散斑的空间范围的2D相关性核被定义在所述像素周围。在尺寸上对应于下文公式(4)中的M和J的散斑被定义为在初始复图像G₀(x,y)的二维自相关函数的两个维度中的半高宽。该核然后与跟随变形、即推动后的复图像G_n(x,y)互相关。作为2D滞后(k,l)的函数，在第n采集处的像素(x,y)处的得到的2D加权的互相关系数ρ’_k,l被计算为：

在该表达式中，W_ij是在(2M+1)x(2J+1)点相关性核上的样本二维加权函数。所述加权，通过加权平均的权重，降低了在搜索最优(或“最大”)滞后中的跳峰和总体误差的可能性。能够通过使用在尾部处平滑减小到零的加权函数(诸如汉宁窗)来降低高频噪声。

相关系数是具有±K和±L的范围的滞后(k，l)的单位归一化复函数。

如在等式(1)至(3)中，当仅1D相关性和1D搜索被用于跟踪轴向运动时，以上公式被简化到L＝0和J＝0，从而降低到：

对应的相关性核236的尺寸为2M+1。针对G_n(x,y)(0＜n≤N)，等式(5)的右侧的值“x”由开始点234、即(x_p+[D_AC],y_q)替换。滞后k，并且因此范围[-K,K]，在等式(5)的迭代运行上增大。因此，搜素区的范围是相关性核236的中心的任一轴向侧上的M+k像素间隔单元222。与结果为相关系数ρ’_k,n(x,y)的最大值的内容相关联的滞后被用在内插中，以找到相关系数曲线的峰和对应或者“最大”的滞后k_maxI。因为k_maxI在轴向方向136上，因此其为向量。来自步骤一的估计[D_AC]226类似为在相同方向136上的向量。两个向量(或者，在1D情况下等价地，标量)被相加以导出针对像素(x_p,y_q)的经修改的位移，即，位移的经改进的测量结果。在2D情况下，相加的两个向量都是二维的。作为使用开始点234作为针对G_n(x,y)(0＜n≤N)的“x”的值的替换的备选，(x_p-[D_AC],y_q)可以被用作针对G₀(x,y)中的“x”的替换。

操作性地并且参考图3中的示范性流程300，采集计数器被初始化到零(步骤S302)。

探头128被定位为间隔开120并且与患者或对象接触(步骤S304)。

采集参考图像227，即，进行成像采集(步骤S306)。

具体而言，在运行步骤S306过程中，图像采集子流程310被调用，并且被如下地运行。横向位置计数器被初始化到零(步骤S308)。成像深度计数器类似地被初始化到零(步骤S310)。超声脉冲140由在当前横向位置124中的探头128发射(步骤S312)。从源自于当前成像深度148的超声回波采集RF信号144(步骤S314)。如果其不是最后深度(步骤S316)，则深度计数器增大(步骤S318)并且返回到样本接收步骤S314。相反，如果其是最后深度(步骤S316)，则深度计数器被清理(步骤S319)。如果当前位置124不是最后位置(步骤S320)，则位置计数器增大(步骤S322)并且处理分支回到脉冲发射步骤S312。另一方面，如果当前位置是最后位置(步骤S320)，则复解析信号从采集的RF数据来导出(步骤S324)。所导出的解析信号与存储的RF信号对应地被存储(步骤S326)。备选地，如果例如经由多项式拟合技术完成相关性滞后内插，则可以不生成或者不需要复解析数据。类似地，复基带数据可以充当步骤一和步骤二中的RF数据的角色。

再次参考主流程300，发射推动脉冲(步骤S328)。此处，假设针对整个流程300发射单个推动脉冲，但是备选地之后能够存在额外的一个或多个推动脉冲，以允许在需要额外的跟踪的时间段期间以来自患者运动的更多全局背景噪声为可能代价的额外跟踪。当前，如果图像采集仍然继续(步骤S330)，则采集计数器增大(步骤S332)，图像采集子流程310被调用以用于跟踪(步骤S334)，并且返回步骤S330。然而，如果图像采集不再继续(步骤S330)，则处理指向针对其位移要被测量的采集的数据中的第一横向位置124(步骤S335)。处理指向当前横向位置处的第一像素(x_p,y_q)。根据等式(1)的自相关在零滞后处执行(步骤S338)。等式(3)被用于提取复值自相关系数R_n的相位(步骤S340)。相位基于先前采集的相位确定而被展开(步骤S342)。等式(2)被用于估计针对当前像素的轴向位移D_AC(步骤S344)。D_AC被量化以产生量化的位移[D_AC]266(步骤S346)。被用在自相关中的相关性核203现在被移位，从而得到针对随后的互相关的移位的核236(步骤S348)。相关性滞后k的绝对值基于例如被预期为校正来自步骤一的估计226所需的最小搜索的内容而被初始化(步骤S350)。例如，典型100μm预定义互相关搜索区能够被降低到30μm，这是所提出的两步方法的益处。当滞后k在-K与K之间变化时，根据等式(5)的互相关迭代到迭代地利用复解析信号(步骤S352)。这不是对步骤一的限制—RF、复解析或者复基带信号的任意组合可以已经被使用。从在迭代中生成的一组系数ρ’_k,n(x,y)，选择最大系数ρ_最大，并且其相位240被确定，如针对自相关系数所做的(步骤S354)。选择具有相反极性的相位242的迭代间最邻近系数，并且内插被用于根据零交叉244确定内插的相关性滞后_maxI以及对应的内插的最大互相关系数ρ_maxI(步骤S356)。如果ρ_maxI并不满足图像匹配阈值或“准则”T_IM(步骤S358)，则滞后k增加(步骤S360)。进行关于滞后k现在是否大于最大滞后阈值T_ML的查询(步骤S362)。如果滞后k不大于T_ML(步骤S362)，则针对由步骤S360的仅先前实例引入的k的中间值中的每个重复等式(5)的互相关。然后，返回到步骤S354。另一方面，如果ρ_maxI满足T_IM(步骤S358)或者滞后k超过T_ML(步骤S362)，则执行向量相加，将量化的位移[D_AC]226与内插的相关性滞后_maxI加和(步骤S366)。如果更多像素(x_p,y_q)仍然要被处理(步骤S368)，则像素指针增大(步骤S370)并且返回到自相关步骤S338。否则，如果没有更多像素仍然要被处理(步骤S368)，但是更多空间位置要被处理(步骤S372)，则位置指针增大(步骤S374)并且处理返回到步骤S336。

已经使用体内临床数据验证了以上在本文中所提出的方法。相比于其他算法，在后处理中，针对两步方法找到了经改进的运动跟踪成功率。更高的运动跟踪成功率意味着在相同检查时间下针对更好的诊断性能的更有效的临床数据。

尽管本发明的方法能够有利地被应用于提供针对人或动物对象的医学诊断，但是本发明的范围不这样被限制。更宽泛地，本文公开的技术针对体内或体外的经改进的基于剪切波的运动跟踪。

尽管已经在附图和前面的描述中详细图示和描述了本发明，但是这样的图示和描述应当被认为是说明性或示范性的，而非限制性的；本发明不限于所公开的实施例。

例如，备选的基于时间域的运动跟踪算法的范例是块匹配最小差异方法。

本领域技术人员通过研究附图、公开内容以及权利要求书，在实践请求保护的本发明时能够理解并且实现对所公开的实施例的其他变型。在权利要求书中，词语“包括”不排除其他元件或步骤，并且词语“一”或“一个”不排除多个。权利要求中的任何附图标记不应被解释为对范围的限制。

计算机程序能够临时、暂时或者针对更长时间段存储在诸如光学存储介质或固态介质的适合的计算机可读介质上。这样的介质仅在不是暂态传播信号的意义上是非暂态的，但是包括诸如寄存器存储器、处理器高速缓存、RAM的其他形式的计算机可读介质。

单个处理器或其他单元可以履行权利要求中记载的若干项目的功能。尽管某些措施被记载在互不相同的从属权利要求中，但是这并不指示不能有利地使用这些措施的组合。