CN100515345C - 通过超声造影成像中的补给曲线拟合进行的血液流动估计 - Google Patents

Abstract

用于非侵入性地量化可通过破坏-补给过程获得的组织灌注的系统和方法，通过导出至少一个局部组织灌注值，来提供表示再灌注期间的局部制剂浓度的信号。该系统包括用于实现以下功能的装置：将具有S形特性的时间函数校准或联系到与再灌注期间的局部制剂浓度成比例的信号，并且在具有S形特性的函数的至少一个参数的至少一个值与至少一个局部组织灌注值(如平均速度、平均渡越时间、平均流动、灌注体积)或属性(如血液流动模式、流动分布方差或斜度)之间建立对应关系。

Description

通过超声造影成像中的补给曲线拟合进行的血液流动估计

技术领域

本发明涉及一种包括充气微泡的破坏/补给监控的过程中的血液流动估计技术，其中流动参数是根据对补给动力特性的分析导出的。更具体而言，本发明涉及活体组织中一种非侵入性地量化灌注(perfusion)的方法。本发明还涉及一种用于执行该方法的计算机程序，以及一种包含该程序的产品。此外，本发明还涉及一种用于非侵入性地量化灌注的相应系统，以及用于此系统中的一种装置。

背景技术

利用载液中的气泡的悬浮液作为有效的超声反射体是本领域中公知的。开发这些悬浮液作为增强超声成像的手段遵循了以下早期观测结果：向静脉迅速注射水溶液会导致溶解的气体通过形成气泡而离开溶液。由于血管内的气泡相对于血液的声阻抗差异很大，我们发现这些血管内的气泡是超声的极佳反射体。将载液中的气泡的悬浮液注射到活生物体的血流中大大加强了超声回波描记(echography)成像，从而增强了内部器官的可视化。由于器官和位于深处的组织的成像对于确立医学诊断可能是至关紧要的，因此已投入了大量努力来开发高浓度气泡的稳定悬浮液，这种稳定悬浮液同时又是易于制备和施用的、包含最少量的无活性物种，并且能够长时间存储且易于分发。

但是简单地将自由气泡散布到水介质中实际上引起的兴趣有限，这是因为这些气泡一般没有稳定到足以用作超声造影剂。

因此，人们已对下述方法产生了兴趣：这种方法例如利用乳化剂、油、增稠剂或糖，或者通过在多种系统中混入或封入气体或其前体(precursor)，来稳定气泡，以便用于回波描记和其他超声研究。这种稳定后的气泡在本领域中一般被称为“微泡”，并且可以分成两个主要的类别。第一类稳定后的气泡或微泡在本领域中一般被称为“微气泡”，并且包括这样的水悬浮液：其中气泡被包含表面活性剂(即两性分子材料)的非常薄的包膜限制在气/液分界面处。第二类微泡在本领域中一般被称为“微球”或“微囊”，并且包括这样的悬浮液：其中气泡被天然或合成聚合物所形成的固体材料包膜所包围。另一种超声造影剂包括聚合物或其他固体的多孔微粒的悬浮液，这种悬浮液承载着夹带在微粒的孔中的气泡。

本发明尤其涉及(虽然不限于)将包含微气泡的水悬浮液的超声造影剂(UCA)用于实施一种包含所述UCA的灌注、破坏和补给监视的技术。充气微球也可很方便地用于本技术。

微气泡通常被定义为本质上是由一层两性分子材料所稳定的充气微泡。微气泡的水悬浮液通常是通过如下方法来制备的：将诸如预先形成的冷冻干燥的脂质体或者冻冻干燥或喷雾干燥的磷脂液这样的粉状两性分子材料与空气或其他气体相接触然后与水载液相接触，并且搅动以生成微气泡悬浮液，该微气泡悬浮液必须在制备后短时间内施用。

适当的充气微泡(尤其是微气泡和微球)水悬浮液及其制备的示例例如在以下专利申请中公开：EP 0458745、WO 91/15244、EP0554213、WO 94/09829和WO 95/16467。

1998年，研究者提出了在从成像平面进行破坏之后通过超声成像仪器监视基于微气泡的超声造影剂(UCA)的补给速率(Wei，K.、Jayaweera，A.R.、Firoozan，S.、Linka，A.、Skyba，D.M.和Kaul，S.，“Quantification of Myrocardial Blood Flow With Ultrasound-InducedDestruction of Microbubbles Administered as a Constant VenousInfusion”，Circulation，vol.971998)。这种局部破坏微气泡的可能性本质上具有以下用途：在测量时在器官中向像平面提供所谓的“阴性丸状”制剂，该器官在另外的情况下则基本处于恒定制剂灌注的状态。在连续(即所谓的“实时”)或间歇(即被触发)成像的状态下，观察像平面中UCA的重灌注速率允许了对器官灌注的估计，即对局部流动参数的估计。

此技术已被广泛采用。广泛发表的文献一直在报告利用作为时间函数的补给视频或多普勒信号的最优拟合，该最优拟合利用描述单室体积中的指示剂的稀释动力特性的表达式(增长的一级指数函数)。

例如参见以下出版物：

K.Wei，，Detection and Quantification of Coronary StenosisSeverity With Myocardial Contrast Echocardiograph，Progress inCardiovascular Diseases，44(2)，2001，81-100：图8示出视频强度与脉冲间隔之间的关系，该关系被拟合到一级指数函数。

Kevin Wei、Elizabeth Le、Jian-Ping Bin、Matthew Coggins、Jerrel Thorpe、Sanjiv Kaul，Quantification of Renal Blood Flow WithContrast-Enhanced Ultrasound，J.Am Col1 Cardiol，2001；37：1135-40：图2示出了视频强度与脉冲间隔的一级指数关系。

Kharchakdjian，R.、Burns，P.N.和Henkelman，M.，FractalModeling of Microbubble Destruction-Reperfusion in UnresolvedVessels，IEEE Ultrasonics Symposium，2001：此论文论述了不同生理流动条件下不同类型的再灌注浓度-时间曲线。

Rim，S.-J.、Leong-Poi，H.、Lindner，J.R、Couture，D.、Ellegala，D.、Masson，H.、Durieux，M、Kasse，N.F.和Kaul S.，Quantificationof Cerebral Perfusion with Real-Time Contrast-EnhancedUltrasound，Circulation，vol.104，2001，2582-2587：图2和3示出了通过一级指数函数拟合的声强度-时间图，而所记录的数据被描述为与制剂浓度成比例。

Schlosser等人，Feasibility of the Flash-Replenishment Conceptin Renal Tissue：Which Parameters Affect the Assessment of theContrast Replenishment？，Ultrasound in Med.&Biol.，Vol.27，pp937-944，2001：此文章分析了造影剂补给，并且还用Wei等人所介绍的一级指数函数来应用了非线性曲线拟合。

Murthy TH、Li P、Locvicchio E、Baisch C、Dairywala I、Armstrong WF、Vannan M.，Real-Time Myocardial Blood FlowImaging in Normal Human Beings with the use of MyocardialContrast Echocardiograph，J Am Soc Echocardiogr，2001，14(7)：698-705：图7示出视频强度-时间曲线被用“单相指数关联方程”拟合。

WO 02/102251描述了微气泡破坏/补给，并在其图2b中示出了微脉管视频强度-时间的一级指数函数，根据该函数微脉管流动强度被描述为由一级函数的起始斜率的正切表示。其图2c示出了视频强度-脉冲间隔的一级指数函数(间歇模式)。

本发明的发明人已观察到现有的启发式方法给出了鼓舞人心的结果，这是因为所有回波成像仪器中的回波信号在以视频信号形式提供给观察者之前，都经历严重的非线性压缩(也称为对数压缩)。从而利用一级指数函数拟合视频数据允许了产生与实际局部器官灌注相关的流动估计，这种流动估计到目前为止已被判断为符合要求。

但是，本发明的发明人已观察到已知方法对于用户选择的仪器设置非常敏感，所述仪器设置例如是接收器增益、对数压缩等等。所提取的参数对于每种仪器类型也是特定的，从而无法在使用不同设备或设置的研究者之间进行比较。此外，从现有技术提取的灌注参数只是相对估计，不适用于对流动参数进行绝对量化的估值。

发明内容

本发明是由以下事实所激发的：从指示剂稀释理论(所谓的“一级指数函数”)导出的表达式不能正确描述像平面或断层切片中破坏之后的UCA再灌注动力特性。

本发明通过提供一种新方法来针对解决上述问题，从而允许基本上独立于用户和仪器的流动参数估计，并且就绝对物理而言具有物理意义。

在一个方面中，本发明提供了一种非侵入性地量化活体组织中的灌注的方法，该方法开始于以下步骤：提供指示组织中的成像造影剂的补给的回波信号序列。S形参数时间函数被与回波信号相关联(校准或联系)。然后在S形函数的一个或多个参数的至少一个值与至少一个局部组织灌注值(例如平均速度、平均流量、灌注体积)或属性(例如血液流动模式)之间建立对应关系。

本说明书和权利要求书中定义的S形函数、或具有S形特性的函数是这样的数学函数：其包括具有基本上恒定的初始值的初始部分，具有基本上恒定的最终值的最终部分，以及初始部分和最终部分之间的中央部分，在所述中央部分中所述S形函数从初始值单调地变化到最终值。优选地，所述函数在其中央部分至少具有一个零值二阶导数。此外，所述S形函数在其初始部分和最终部分中优选地具有基本上为零的一阶导数。S形函数的示例是“误差函数”、双曲正切、反曲函数、累积正态分布函数、累积对数正态分布函数或其任何多项式近似。

在本发明的一个实施例中，所述造影剂包括能够反射声能的微泡。所述提供回波信号序列的步骤包括在超声成像装置的成像平面中施加超声脉冲；所述超声脉冲是以高到足以导致该平面内存在的微泡的有效部分被破坏的声压来施加的。然后更多的超声脉冲的序列被施加在所述成像平面中；所述更多的超声脉冲的声压足够低到保持所述微泡的大部分。施加更多的超声脉冲的序列的步骤在预定的后续时刻被重复；然后源自所述平面的由所述更多的超声脉冲产生的回波信号被记录，以便监视所述后续时刻处成像平面内的微泡的补给。

作为进一步的增强，在关联所述S形函数之前所述回波信号被处理。具体而言，使回波信号与微泡的局部浓度成比例；从而产生与成像平面内的任何位置处的造影剂的浓度成比例的处理后的回波信号。基于超声波束几何形状和UCA破坏的程度，在一个实施例中，本发明将上述S形函数的最佳拟合所确定的成像平面上的平均流动速度与达到稳态UCA浓度的一半所需的时间相联系。优选地，S形灌注函数被拟合在UCA所生成的回波信号上，所述回波信号的瞬时幅度与生成这些回波信号的UCA的局部浓度成比例。此比例关系通常是通过对从超声分析获得的最终数据进行适当的线性化来获得的，所述最终数据一般是两种类型的。第一类数据被称为“图像”，并且包括显示为模拟或数字视频信号或者任何其他灰度级或彩色幅度2维图(2D图)的回波信号，这些图像是通过包括对生成的回波信号进行非线性动态范围压缩(例如对数压缩)在内的过程来获得的。图像通常与具有预定幅度的图像元素(象素)相关联。第二类数据被称为“原始”回波信号，并且包括具有与超声回波幅度成比例的幅度的信号，通常是直接从超声装置获得的射频(rf)回波信号。

这里所使用的术语“线性化”或“线性化的信号”适用于被以下述方式处理过的超声回波信号：这种方式使得信号幅度与产生该信号的局部UCA微泡浓度成正比。这反映了随机间隔的散射体群体对声能的散射性质，该散射性质导致回波信号功率与UCA浓度成比例。当处理与声压幅度成比例的rf或解调rf信号时，这种线性化可通过对原始回波信号幅度适当取平方来获得的，而当处理对数压缩后的图像时，这种线性化可通过适当的反对数压缩并对每个象素的幅度值取平方来获得的，从而获得处理后的信号的幅度与UCA浓度之间的比例关系。

或者，S形函数也可被拟合在“未线性化”的回波数据(例如图像)上，即不与局部UCA浓度成比例的数据。在此情况下，所选择的可被拟合在未线性化的数据上的S形函数被与引起UCA浓度和回波数据之间的非线性的过程(例如S形函数的平方根和对数压缩)相同的过程来修改。

被执行曲线拟合的回波信号通常是通过在超声装置的成像平面中施加超声脉冲序列来获得的，成像平面内微泡的补给是通过记录源自所述平面中包含的微泡的所述超声回波信号的幅度(作为时间的函数)来监视的。

一个用于提供可指示组织中的成像造影剂补给的回波信号序列的实施例可以包括以下步骤：

-向器官或感兴趣的区域提供恒定UCA供应，或者在破坏-补给方法所要求的几秒钟期间允许充分恒定的UCA灌注的情况下以丸剂方式注射UCA；

-紧挨施加UCA破坏帧之前，记录感兴趣的区域中的回波信号，其中声能级别低于能够破坏所述微泡的预定阈值；

-通过超声脉冲在超声成像装置的成像平面中施加所述UCA破坏帧，所述UCA破坏帧的声能高于所述预定阈值，并且足以引起该平面内存在的微泡的有效部分的破坏；

-紧接施加所述UCA破坏帧之后，记录感兴趣的区域中的回波信号，声能级别低于所述预定的微泡破坏阈值。

以下步骤可在任何可指示组织中的成像造影剂补给的回波信号序列(例如可根据上述方法获得的回波信号序列)上执行，以确定所述组织的灌注参数：

-局部执行所述信号对S形函数(例如“误差函数”家族)的最佳拟合估计，所述最佳拟合估计或者可在用户选择的感兴趣区域(AOI)内在局部逐个象素级别上进行，以便允许生成灌注的参数图像，或者可在象素群组级别上进行，所述象素群组是在通过用相干波束获得超声图像时超声图像的斑点性质所确定的，

-在具有S形特性的函数的至少一个参数的至少一个值与至少一个局部组织灌注值(例如平均渡越时间、平均速度、平均流量、灌注体积)或属性(例如血液流动模式、流动分布方差或斜度)之间建立对应关系。

优选地，在执行最佳拟合估计之前，回波信号被处理(例如线性化)，以获得与局部UCA浓度成比例的处理后的信号，以便对与UCA浓度成比例的数据执行所述最佳拟合。

具有S形特性的函数的参数值与局部组织灌注值之间的对应关系的一个示例是，可通过计算破坏帧所破坏的区域(或切片)的半厚度除以局部参数拟合所确定的平均渡越时间之比率来确定平均流动速度。

另一个示例是作为最佳拟合值找到的幅度一旦被校准，就可被解释为与被分析区域中血液体积成比例的量，并用于根据幅度和流动速度值之积估计流量值。

另一个示例是不同流动分布(例如方差和斜度)的分布范围也可根据最佳拟合值找到，以模拟各个灌注函数之和。

这些流动参数估计可被显示在所使用的AOI内，或者以参数图像的形式显示为二维图。

根据上述方法获得的S形函数也可用于估计再灌注期间微泡的渡越时间或流动速度的概率密度分布(例如不同毛细管中的)。

为此，在本发明的一个实施例中，S形函数被通过小波分解方法分析。

优选地，与局部浓度成比例的回波信号在被通过小波分解方法分析之前被微分两次。

确定用于分解的母小波的一个建议选择包括用于描述单个流动值的S形函数的累积正态分布函数的二阶导数。

在另一个的实施例中，回波信号被通过单步骤或多步骤过程分析，以估计不同流动渡越时间或速度处的贡献的分布。

具体而言，选择第一组流动渡越时间或速度；然后通过多个S形函数的线性组合与回波信号的最佳拟合来进行第一估计。

在优选实现方式中，进行第二估计以定义第二组流动渡越时间或速度；以此情况下，第一估计被用作定义所述第二组的基础。

更优选地，第二估计被用于提供进行第三估计的初始值集合。

第二估计的推荐选择包括使用立方样条外推。

进一步改进方案的方式是使用神经网络分析来进行第三估计。

通常，神经网络由多个权重(用于流动渡越时间或速度)以及多个偏置值(用于加权后的流动渡越时间或速度)所定义。通过反复调整偏置值和权重训练神经网络；优选地偏置值和每个负权重被周期性地重置为零。

在优选实施例中，重置是以50到200次迭代为周期执行的。

作为进一步的增强，第一估计最多是用16个S形函数来进行的，优选地最多是用8个S形函数来进行的。

此外，第二估计至少是用一组8个流动渡越时间或速度来进行的，优选地至少是用一组16个流动渡越时间或速度来进行的.

本发明的另一个方面提供了一种用于执行上述方法的计算机程序。

本发明的另一个方面提供了一种包含此计算机程序的程序产品。

此外，本发明的另一个方面提供了一种用于非侵入性地量化活体组织中的灌注的相应系统。

本发明的另一个方面提供了一种用于该系统中的装置，这种装置具有用于输入所述回波信号的装置。

正如下文中将要证明的：

(1)通过根据本发明的方法从S形函数中提取的参数是独立于所使用的设备和设置的，并且可在使用不同设备或设置的研究者之间进行比较。此外，这些提取的参数适用于绝对量化估值。

(2)本发明提供了一种简单的用于体内血液流动参数的量化方法，尽管实际器官灌注非常复杂，这种方法的效果却惊人的好。即使在存在具有不同速度和来自任意不同方向的叠加的流动贡献的情况下，也可对沿着与扫描平面垂直的方向的平均流动速度进行估计。

(3)根据本发明，曲线拟合可基于这样的参数表达式，这种参数表达式集成了超声探头在垂直方向上的声灵敏度属性的，同时又考虑到了回波信号功率与局部微泡浓度之间的对应关系。于是就足以在对所选择的拟合函数进行任何修改时都能保持这些对应属性。利用此条件，以及对微泡破坏区域在垂直方向上的宽度的了解，可以估计绝对血液流动参数，及其统计分布(例如流动渡越时间或速度)。

附图说明

在以示例方式给出的附图中：

图1a示出从超声换能器发射的声束；

图1b示出超声成像仪器在垂直方向的典型波束灵敏度分布；

图2示出在造影剂再灌注已被破坏的区域时造影剂信号功率的S形函数，这种S形函数被用于根据本发明的方法中；

图3a和3b示出根据本发明的一个实施例的实验记录的破坏之后的回波功率函数，该函数是在体外设置下记录的，并且被用累积对数正态分布函数拟合；

图4a和4b分别示出利用现有技术一级指数拟合进行的流动速度估计以及利用根据本发明的方法进行的相应的估计；

图5a和5b示出流动渡越时间的两个对数概率分布函数；

图6a和6b示出从图5a和5b的流动分布产生的S形回波功率再灌注函数，该函数可用于本发明的另一实施例中；

图7a和7b示出实验记录的破坏之后的回波功率函数，该函数是在体外设置下记录的，并且是根据本发明的另一实施例利用根据渡越时间的S形对数正态分布的再灌注函数之和拟合的；

图8示出被发现是对图7的数据的最优拟合的对数正态渡越时间分布；

图9a-d示出本发明的另一实施例中的再灌注函数的小波分解过程，该分解是就其各个流动贡献的分布而言的；

图10至14示出通过单步骤过程或多步骤过程中的第一步骤进行的回波功率数据分析，该分析用于估计不同流动渡越时间或速度下的贡献分布；

图5示出神经网络；

图16至34示出通过多步骤过程的后续步骤进行的回波功率数据分析，该分析用于估计不同流动渡越时间或速度下的贡献分布；

图35至37示出单步骤或多步骤过程的其他应用；

图38的框图示出典型医学超声成像系统的主要元件；

图39的框图示出根据本发明的参数流动成像的一个特定示例中的流动估计方法/系统的主要功能元件；

图40a和40b示出两个感兴趣区域的补给函数的实验结果，所述两个区域分别处于肾皮层的正常灌注和灌注不足的区域中。

具体实施方式

在以下描述中，所使用的符号的定义如下：

x	像平面内的空间坐标，与波束正交的(横)向
		y	像平面上的空间坐标，与波束正交的(垂直)方向
z	像平面内的空间坐标，沿波束的(深度)方向
		D	垂直(y)方向上微泡破坏的空间范围
F	超声频率
		c	声速
A	超声波长
		m	自然对数的均值
s	自然对数的标准差
		μ	对数正态分布的均值
σ<sup>2</sup>	对数正态分布的方差
		γ	对数正态分布的偏度
τ	从受破坏区域的边缘到中心的流动渡越时间
		τmean	τ的均值
C	微泡的概率密度或相对浓度
		A	参数方程中的幅度因子
O	参数方程中的偏移因子
		a	y方向上的半孔径宽度
v	局部流动速度(v<sub>x</sub>，v<sub>y</sub>，v<sub>z</sub>分量)
		K<sub>Tx</sub>，K<sub>Rx</sub>	发射和接收处的2α/λz参数
K	由K<sup>2</sup>＝K<sup>2</sup><sub>Tx</sub>+K<sup>2</sup><sub>Rx</sub>确定的发射-接收参数
		Y	无单位变量Y＝Ky
Θ	流动方向和像平面的法向之间的角度

β	一级指数函数的“速度”项
		GL	视频信号的灰度级
t	时间
		dt	小波分析中的时间采样
η	累积对数正态概率分布函数中的偏差项
		Γ	任意比例常数

本发明提出了一种器官灌注量化方法，该方法基于对基于微泡的超声造影剂(UCA)的再灌注动力特性的研究。此方法要求施用UCA、等待一定时间以便在器官或感兴趣区域中确立UCA浓度的稳态(通常在5到30秒之间)，在足够高的声压下应用一个或多个制剂破坏帧，以破坏包含像平面的切片中的UCA微泡，并监视再灌注或补给动力特性。达到所需的稳态灌注可按以下方式通过连续滴注UCA或通过丸剂来实现：该方式使得在持续几秒(通常在1到15秒之间)达到相当恒定的滴注速率。

再灌注动力特性

再灌注动力特性是在正比于局部UCA浓度的线性信号的基础上来分析的，所述信号也就是正比于反向散射的声功率的信号。然后可从补给动力特性推出灌注参数。传统的“B-模式”或“2D多普勒”下的超声成像是层析方法，其中组织切片被聚焦声束的迅速扫描所询问。这种成像模式中的空间分辨率主要由三个维度中的每个深度处的发射-接收超声灵敏度分布所决定，所述三个维度是：横向，此方向被定义为像平面内跨过声束宽度的方向；轴向，像平面内沿着波束传播方向的方向；以及垂直方向，即与像平面垂直的方向。

UCA微泡由于施加足够高的声压而被破坏之处的切片厚度由垂直方向上的声束宽度和实际施加的声级别所确定。在补给阶段期间，UCA微泡重新进入破坏后的切片的体积中，并且被回波描记仪器根据其沿波束垂直方向上的空间灵敏度而检测到。申请人发现的是所观察到的作为时间函数的回波功率和实际再灌注动力特性之间的关系由脉冲-回波模式中沿垂直方向的空间发射-接收分布唯一确定。与此相反，对于此关系一般被接受的概念是从指示剂稀释理论借来的，该理论描述了指示剂被随机稀释在均匀介质中时其浓度随时间的演化(例如Wei，K.、Jayaweera，A.R.、Firoozan，S.、Linka，A.、Skyba，D.M.和Kaul，S.，“Quantification of Myrocardial Blood Flow WithUltrasound-Induced Destruction of Microbubbles Administered as aConstant Venous Infusion(利用超声导致以恒定静脉滴注方式施用的微气泡的破坏来量化心肌血液流动)”，Circulation，vol.97 1998)。先前的研究者的方法大多数是基于造影剂补给期间观察到的视频强度级别的，所述视频强度级别是强烈地由回波描记仪器的所谓“对数压缩”所确定的量。这导致了将一级指数函数选择为再灌注动力特性的模型，所述一级指数函数的一般形式如下：

GL(t)＝A(1-e^-βt)，

其中A是稳态幅度，β是一级指数函数的“速度”项，时间原点取在紧随最后的UCA破坏脉冲之后的瞬间。在现有技术中(例如所引证的Wei等人的文章)，值A、β和Aβ一般被解释为正比于被分析的区域内的“血液体积”、“血液速度”和“血液流量”的量。但是，此方法不是基于正比于给定时刻处局部UCA浓度的函数的，并且由于对诸如增益、对数压缩参数等用户设置非常灵敏而受到损害。

本发明公开了一种简单的量化方法，尽管实际器官灌注(微脉管系统结构、流动的随机方向、不同流动值之和等等)非常复杂，这种方法在体内情形下的效果却惊人的好。

通过具有足够高的声强度的脉冲应用UCA破坏帧导致某个体积中的UCA微泡的不足，该体积在垂直方向上由超声波束压力分布所确定，在横向上由超声探头所扫描的区域范围所确定。确定UCA破坏帧的高声强脉冲可以是在超声成像装置的成像平面内的不同方向上应用的单个脉冲或者最好是多个脉冲(脉冲系列)。UCA破坏帧可以是在所述成像平面中顺序应用的单个帧或多个帧。例如，为了实现较深区域处的微泡的有效部分的破坏，可能需要多个破坏帧。如前所述，所应用的声能级别应该高于预定阈值，并且能够破坏UCA微泡。更具体而言，所述级别应该足够高，以便导致成像平面内出现的微泡的有效部分的破坏。“有效部分”的破坏意思是被破坏的微泡的量应该足够高，以便能够检测到从微泡接收到的回波信号在破坏之后紧接着测量的值与到达稳定灌注状态时的值之间有相当大的变化。在实践中，破坏像平面中至少50％的微泡一般就足以获得可接受的实验数据。优选地，所述被破坏的微泡量至少为75％，更优选为至少90％，而在最优选的情形下达到100％。

参见图1a，回波描记图像中的UCA微泡的重新出现(补给其被破坏的区域)速率一方面由图像中的每个位置处的局部血液灌注速率所决定，另一方面由在基本上垂直的方向上超声探头的声灵敏度模式所决定。这种再灌注的速率，以及更一般而言灌注参数的值，是未知变量，对这种未知变量的估计可向临床医生提供有价值的信息，以用于评估局部组织病理。

垂直方向上的声灵敏度模式是以下将论述的参数，这是因为关于它的知识是理解本发明的基础。当在连续波模式中被激励时，在具有矩形几何形状的聚焦孔隙的领域中，垂直方向y上的声压分布大致由以下函数确定：

$p\left(y\right)\≅\mathrm{\Γ}\·\sin c\left(K_{\mathrm{Tx}}y\right),$

其中Γ是任意比例常数，

$K_{\mathrm{Tx}}=\frac{2a}{\mathrm{\λ}\·z},$ $\mathrm{\λ}=\frac{c}{f},$

“sinc”函数表示

$\sin c\left(x\right)\≡\frac{\sin \left(\mathrm{\πx}\right)}{\mathrm{\πx}},$

其中f是超声频率，c是声音在传播介质中的速度，λ是超声波长，a是换能器在垂直方向上的半孔径，z是从换能器探头到感兴趣的深度的距离，y是垂直方向上偏离轴的距离[例如，Kinsler LE、Frey AR等人，Fundamentals of Acoustics(声学基础)，J.Wiley&Sons，1982]。在脉冲激励情形下，正如回波描记成像模式中一般的情形那样，在接近声脉冲波形的中心(或平均)频率的频率处，峰值压强分布的主瓣与连续波情形非常一致。

由于有兴趣使局部检测到的信号幅度与局部UCA浓度之间产生对应关系，因此最好就回波信号强度来表达灵敏度模式。这些模式可由发射和接收分布的组合效果所确定，其中发射和接收分布一般可能是不同的。在发射处，超声换能器域中垂直方向上的声功率分布P_Tx(y)大致上由以上给出的压强分布的平方所确定。从而对于矩形孔隙，它可由以下形式的函数来表达：

$P_{\mathrm{Tx}}\left(y\right)\≅{\sin c}^2\left(K_{\mathrm{Tx}}y\right),$

其中K_Tx由超声发射状况所确定。在实践中，此功率分布可根据以下公式由高斯函数G_Tx(y)来近似：

$P_{\mathrm{Tx}}\left(y\right)\≅G_{\mathrm{Tx}}\left(y\right)=e^{-{(1.94\·K_{\mathrm{Tx}}\·p)}^2}.$

在接收模式中，类似地对P_Rx(y)的近似为：

$P_{\mathrm{Rx}}\left(y\right)\≅G_{\mathrm{Rx}}\left(y\right)=e^{-{(1.94\·K_{\mathrm{Rx}}\·y)}^2}.$

其中K_Rx由超声接收状况所确定。

在脉冲-回波情形中，在一级近似下，超声换能器对于离轴对象的功率灵敏度PE(y)由发射波束分布P_Tx(y)与接收波束分布P_Rx(y)之积确定。从而发射-接收灵敏度模式PE(y)可由高斯型G(y)近似如下：

$\mathrm{PE}\left(y\right)\≅G\left(y\right)\≡G_{\mathrm{Tx}}\left(y\right)\·G_{\mathrm{Rx}}\left(y\right)=e^{-{(1.94\·y)}^2(K_{\mathrm{Tx}}^2+K_{\mathrm{Rx}}^2)},$

从而导致由K²＝K² _Tx+K² _Rx所确定的发射-接收K值的定义

PE(Y)与G(Y)的紧密对应关系在图1b中示出，其中具有任意的矩形几何形状的换能器的主瓣，无单位量Y≡Ky的值的范围为-1至+1。

在任何情况下，为了描述作为本发明的主要目标的流动估计方法，都不需要讨论此分布的确切形状。超声成像和声学领域的任何技术人员都能够通过确定图像的不同深度或部分处的K和D的值，来使本发明适应于任何给定的超声成像系统的实际的波束灵敏度模式和微泡破坏范围。这种确定可通过为每种探头和操作条件进行一次简单的校准过程来执行，如下所述。

实际上，K值可按上述方式理论地确定，或者可通过测量发射-接收波束灵敏度来实验地确定，所述测量是通过以下方式来进行的：在垂直方向上扫描像平面上的小点反射体，并利用垂直方向上空间位移的最佳拟合高斯函数来调整记录的分布图。

在UCA微泡再次进入切片体积这时，首先假设与成像平面垂直的方向上的运动，则与UCA微泡的局部浓度成比例的线性回波信号由被波束灵敏度加权的垂直方向上发射/接收波束模式所截取的微泡的越来越大的比例所确定。对于已重新进入切片直到位置y’≡Y’/K的浓度均一且为单位值的UCA微泡，以波束灵敏度PE(Y)所检测到的回波描记功率信号E(Y’)从数学上而言可表达为积分：

$E\left(Y^{\′}\right)={\∫}_{-\∞}^{Y^{\′}}\mathrm{PE}\left(Y\right)\mathrm{dY}\≅{\∫}_{-\∞}^{Y^{\′}}G\left(Y\right)\mathrm{dY}.$

利用“误差函数”实现本发明

当考虑与高斯函数G(Y)非常接近的PE(Y)函数的实际性质时，申请人发现用“误差函数”erf(q)来表示回波功率信号E(Y’)将会是有利的，该误差函数的定义如下[例如参见：Abramowitz，M.和Stegun，I.A.(eds.)所著Handbook of Mathematical Functions Dover Publications，Inc.，New-York，1972，pp.297-329中Gautschi，W.的“Error Functionand Fresnel Integrals(误差函数和菲涅耳积分)”]：

$\mathrm{erf}\left(q\right)=\frac{2}{\sqrt{\mathrm{\π}}}{\∫}_o^q{e}^{-p^2}\mathrm{dp}.$

此定义符合以下性质：erf(0)＝0，erf(-q)＝-erf(q)且 $\underset{q\→\∞}{\lim }\mathrm{erf}\left(q\right)=1.$ 微泡从一侧(例如Y的负值)对被破坏的切片进行补给的物理情形证明使用所谓的“累积正态分布函数”是正确的，该“累积正态分布函数”在灌注估计的上下文中应当被称为perf(q)，并且定义如下：

$\mathrm{perf}\left(q\right)\≡\frac{1}{\sqrt{\mathrm{\π}}}{\∫}_{-\∞}^q{e}^{-p^2}\mathrm{dp},$

其验证了以下性质： $\underset{q\→-\∞}{\lim }\mathrm{perf}\left(q\right)=0,$ perf(0)＝0.5，并且 $\underset{q\→\∞}{\lim }\mathrm{perf}\left(q\right)=1.$ 此外，perf(q)可以简单地由erf(q)表示为：

perf(q)＝0.5·(1+erf(q)).

利用上述定义，从图形上来说，perf函数是S形函数，其表示从-∞到Y积分时高斯型发射-接收波束G(Y)灵敏度分布图下的能量，如图2的示例所示，其中初始的基本上平坦的部分(平稳段)对应于微泡破坏后再灌注过程的开始，中间的倾斜部分和最后的平稳段对应于完成灌注的稳态。对于根据本发明的任何S形函数，如果万一所记录的数据的持续时间不足以估计最后的稳态值，则应用破坏脉冲之前所测量的稳态值可用于拟合模型中，作为再灌注信号的期望渐近稳态值。

考虑在像平面的两侧对称延伸的厚度为D的受破坏微泡区域，则可理解在负Y和正Y方向上距离探头轴D/2处确定了一个位置。于是新鲜微泡流以流动速度v在Y方向上对受破坏的切片进行补给。从而在非破坏性监视阶段期间，作为时间函数的回波功率恢复函数E(t)也由perf函数所表示，现以一般形式将其表达为：

E(t)＝perf[1.94·Kv(t-D/2v)]＝perf[1.94·Kv(t-τ)]，

其中D/2v表示微泡从受破坏区域的边缘行进到像平面的中央部分所需的时间延迟τ。通过这种方式，可通过测量局部微泡浓度到达其最大(或稳态)值的一半所需的时间延迟，来以实验方式直接估计流动速度v。此延迟等同于“平均渡越时间”，易通过以下方式来估计：利用参数方程对用实验的线性回波功率信号进行最佳拟合分析，其中该参数方程除包括幅度和偏移项外，还包括延迟参数τ。对于垂直于成像平面的恒定均一流动的情形，这种参数拟合方程F(t)可以采取以下形式：

$F\left(t\right)=O+A\·\mathrm{perf}\left[1.94\frac{\mathrm{KD}}{2}(\frac{t}{\mathrm{\τ}}-1)\right],$

其中参数O、A和τ分别代表偏移、幅度和渡越时间延迟最佳似合值，从而允许了利用关系v＝D/2τ来估计所分析的区域中的流动速度。

进行此流动速度分析仅需的先验知识是扫描平面周围受破坏微泡的局部厚度D，所述局部厚度是可以为每种探头类型、回波描记仪器和操作模式以合理近似作为深度的函数而映射的值。不同深度处的D值例如可通过以下方法以实验方式确定：将UCA微泡嵌入在凝胶中，并通过直接光学观察来测量受破坏微泡的范围。或者，该值也可通过以下方法在体内或体外以声学方式确定：在较低声功率下，使用成像平面与第一系统的成像平面垂直的第二超声成像系统，来显现受破坏的UCA微泡的范围。又或者，可以基于传输束分布图以及对UCA微泡破坏的声压阈值的了解来理论地估计D。通过考虑组织衰减可随着深度在D值上应用校正因子。

注意，在实际做法中，UCA微泡从两侧相等地对受破坏的区域进行补给；上述回波功率补给函数E(t)保持有效，因为不论流动方向如何，只有每个时刻波束内微泡的整体浓度才是有关紧要的。

还要注意，在特定情况下，虽然根据本发明的S形函数已被与回波信号相关联以便估计流动参数，但是数据集合或最佳拟合函数可能不会在所分析的数据集合的有限时间间隔中展现用于进行曲线拟合的数学函数的S形特性。这种情况例如可能发生在被超声波束从垂直方向上询问的体积内只有小部分微气泡被破坏时，或者发生在由于存在大流动速度时帧速率有限因而实际微气泡补给被采样不足时。

在微泡速度矢量和扫描平面的法方向之间的角度θ不再是零而是不同于90°的任意值时，所估计的v和最佳拟合参数τ之间的关系就简单地是：

$v=\frac{D}{2\mathrm{\τ}\cos \mathrm{\θ}}.$

在不知道流动方向θ的情况下，假设θ＝0而估计出的流动速度则对应于垂直于像平面的速度分量v_y。

本发明的一个令人惊奇地不寻常的特征是：即使存在具有不同速度v_i并且来自不同的任意方向θ_i的叠加的流动贡献的情况下，也可对垂直于扫描平面的平均流动速度v_y进行估计。这一点在下文中说明。当由从N个不同方向和/或以N个不同速度流动的微泡做出贡献时，回波信号功率可表达为：

$E\left(t\right)=\underset{i}{\overset{N}{\mathrm{\Σ}}}{C}_i\·\mathrm{perf}[1.94\·{\mathrm{Kv}}_{i_y}(t-{\mathrm{\τ}}_i)],$ 其中 ${\mathrm{\τ}}_i=\frac{D}{{2v}_{i_y}},$ $v_{i_y}=v_i\cos {\mathrm{\θ}}_i,$

C_i是速度为

的微泡的相对浓度，其中速度

定义为沿y方向的个体相对速度。在这种情况下，对E(t)的形状的影响是它不再是纯perf函数。

利用“logperf”再灌注函数实现本发明

为了对这种实验性的再灌注功率数据(即如上所述由不同流动之和产生)进行成功的曲线拟合，我们发现使用累积对数正态概率分布函数形式的经验性S形参数函数将会是有益的。以下将在目前的上下文中所使用的这一函数称为logperf再灌注函数：

$\log \mathrm{perf}\left(t\right)=O+\frac{A}{2}\·(1+\mathrm{erf}\left(\frac{\ln (t/\mathrm{\τ})}{\mathrm{\η}\sqrt{2}}\right)),$

其中O、A、τ和η是拟合参数。令人惊讶的是，用这种logperf(t)函数来拟合实验性的再灌注回波功率信号允许了对由最佳拟合值τ所给出的平均渡越时间τmean的良好估计。使用logperf参数函数作为拟合函数在图3a和3b中示出。为了获得实验数据，采用了与Veltman等人所述的设置类似的设置(“On the design of a capillary flowphantom for the evaluation of ultrasound contrast agents at very lowflow velocities”，Ultrasound in Med.Biol，28(5)，625-634，2002)。在蠕动泵(Gilson Minipulse 3，Villier等人，法国)的控制下以从1到30mm/s的流动速度值用微气泡悬浮液(磷脂稳定的含全氟丁烷气体的微气泡，平均直径约为2-3μm)灌注一束微纤维(170根内径为240μm的纤维，Hospal AN69，法国)。在蒸馏水中稀释微气泡相当于约每毫升10⁶个气泡的量流经微纤维。

所使用的回波描记仪器是Esatune扫描仪(Esatoe，佛罗伦萨，意大利)，该扫描仪被用于CnTI(S)(实时造影成像)模式中，采用在2.6MHz发射频率和43Hz帧速率下工作的LA532E探头，采用0.07的机械瞄准标线以便进行再灌注监视，并且使用最大发射功率下的一秒钟的破坏帧序列。按以下方式放置探头：产生微纤维束的截面图像，在流动方向和像平面法向之间形成50°的角度θ。对于每个破坏-补给序列，以15秒为周期利用快速rf-抓取器收集实时rf-信号(FEMMINA，Scabia等人，“Hardware and software platform for processing andvisualization of echographic radio-frequency signals(处理和可视化回波描记射频信号的硬件和软件平台)”，IEEE Trans.Ultra.Ferr.Freq.Contr.，49(10)，1444-1452，2002)。然后计算用户定义的AOI内的平均回波功率信号，以获得流动速度值分别为4和8mm/s的情况下破坏之后微纤维中UCA微气泡的补给，如图3a和3b中的虚线数据点所示。

在拟合之前，通过以光学方式测量每毫升约含5·10⁶个微气泡的1.5％琼脂凝胶中的受破坏区域的宽度，确定D值为7.2mm。在Matlab曲线拟合工具箱(MathWorks，Natrick，MA，美国)中实现logperf函数，并且利用可信区域方法计算最佳拟合(Byrd，R.H.、R.B.Schnabel和G.A.Schultz，“Approximate Solution of the Trust RegionProblem by Minimization over Two-Dimensional Subspaces(通过二维子空间上的最小化获得可信区域问题近似解)”，MathematicalProgramming，Vol.40，pp 247-263，1988)。最佳拟合在图3中以实线示出。对于图3a和3b，所找到的τmean值分别为1.421和0.6714，η值分别为0.6197和0.8079，其中A分别为A＝1114和1258。利用这些τmean、θ和D值，对于图3a和3b，所估计的流动速度值则分别是3.9和8.3mm/s，这与实际值4和8mm/s良好吻合。

在图4a和4b中，本发明比起基于一级指数拟合进行的估计来说的主要优点是，利用上述体外设置，越来越大的流动速度值、60°的角度θ以及两个不同的仪器增益值(高增益＝低增益的2倍)。图4a示出利用现有技术一级指数最佳拟合获得的β值。注意到估计的线性良好，但是在所使用的两个增益值之间有显著差异。例如，在20mm/s的流动速度下，对于低增益所估计的β值约为3.3s^-1，而对于高增益约为4.7s^-1。这表明现有技术方法是依赖于系统的。图4b示出利用logperf函数的最佳拟合值的情况下所估计的D/(2τ)值。注意到估计的线性良好，并且对于所测试的两个增益值，估计的绝对值都在真实值的5％之内(对于回归线 $y\≅0.95x$ )。在同样的流动速度为20mm/s的示例中，对于低增益值和高增益值，所估计的D/(2τ)值都约是19mm/s。这表明利用根据本发明的logperf函数所估计的流动速度是独立于系统和绝对物理量的。从而由根据本发明的方法提取的参数独立于所使用的设备或设置，因此可在使用不同设备或设置的研究者之间进行比较。此外，这些提取的参数适用于绝对量化估值。

对于“对数正态”流动分布实现本发明

UCA的补给分析也可用于评估未知组织内的对数正态流动分布的均值、方差和偏度的最佳估计。这可以提供关于微脉管网络的组织的信息。在平均渡越时间为τ的血管的对数正态分布的情况下，组织灌注通常表示为准连续体[Qian H.和Bassingthweighte A.，A Class ofFlow Bifurcation Models with Lognormal Distribution and FractalDispersion(利用对数正态分布和分形散布的一类流动分歧模型)，J.of Theoretical Biology，2000，205，261-268]。当用按与血液相同的速度流动的微泡灌注组织时，已观察到微泡浓度分布C(τ)也可由以下形式的对数正态概率密度分布描述：

$C\left(\mathrm{\τ}\right)=\frac{e^{-\frac{{[\ln \left(\mathrm{\τ}\right)-m]}^2}{{2s}^2}}}{\mathrm{\τs}\sqrt{2\mathrm{\π}}},$

其中m和s分别是τ的自然对数的均值和标准差。对数正态分布的均值η、方差σ²和偏度γ由以下式子给出：

$\mathrm{\μ}=e^{m+\frac{s^2}{2}},$ ${\mathrm{\σ}}^2=e^{s^2+2m}(e^{s^2}-1),$ $\mathrm{\γ}=\sqrt{e^{s^2}-1}(2+e^{s^2}).$

上述渡越时间的概率密度验证了它被规一化为单位1：

$\underset{0}{\overset{\∞}{\∫}}C\left(\mathrm{\τ}\right)\mathrm{d\τ}=1.$

正如申请人所观察到的，对未知组织内的对数正态流动分布的均值、方差和偏度的估计可通过以下方式来实现：将回波功率E(t)表达为由流动渡越时间τ的对数正态分布加权的各个perf函数的组合。此函数被称为lognormperf(t)，由下式给出：

$\log \mathrm{normperf}\left(t\right)=O+A{\∫}_0^{\∞}C\left(\mathrm{\τ}\right)\mathrm{perf}(1.94\·\frac{K\·D}{2\mathrm{\τ}}(t-\mathrm{\τ}))\mathrm{d\τ},$

所有变量的定义如前。

图5a和5b以模拟示例的方式示出了渡越时间的两个对数正态分布，对于其一σ²＝0.5和γ＝1.88，对于第二个σ²＝0.02和γ＝0.34，两者都满足τmean(＝μ)＝1.25。在示例性的参数O＝0，A＝1，D＝7.2mm，θ₁＝0以及K＝1.41mm^-1的情况下，这两个分布给出图6a和6b所示的补给函数lognormperf(t)＝E(t)。注意图6的两个示例之间的差异证实了lognormperf(t)＝E(t)的形状取决于渡越时间的分布的偏度，对于低偏度分布(γ＝0.34)非常接近perf函数，而对于较高的偏度分布(γ＝1.88)是上升较急剧且“肩部”(朝着接近稳态值的弯曲)较缓和的S形。

对于渡越时间的对数正态统计分布的情形，这些后来的考虑在图7a和7b中利用与图3a相同的实验数据以示例方式示出。这里，数据点被用参数表达式lognormperf(t)拟合，其中包括具有m和s参数的流动渡越时间的对数正态分布。通过考虑在发射和接收中相同的测量物理条件，即a＝2.5mm，f＝2.6MHz，z＝45mm，c＝1480mm/s，K值被计算为0.2761mm^-1。对于最佳拟合(D＝7.2mm且偏移O＝0)找到的值如下：对于4和8mm/s的流动速度分别为：μ＝1.53和0.78，σ＝0.47和0.39，γ＝0.96和1.6，A＝1114和1258。流动渡越时间的两个相应的对数正态概率密度分布在图8中示出，在θ＝50°的情况下，相应的估计出的流动速度为3.7和7.2mm/s，又一次分别与实际施加的流动速度4和8mm/s合理地吻合。

图6a和6b的补给函数的显著特征是：尽管两种流动速度条件导致具有很不相同的偏度值的分布，但是相应的估计出的平均渡越时间值μ＝τmean就量化意义而言与实际平均流动速度是吻合的。

从而，根据本发明的在气泡破坏之后的参数流动估计允许了仅基于对每个感兴趣的深度处被破坏的气泡的厚度D的了解，来直接估计垂直于成像平面的方向上的平均流动速度。

在流动速度或渡越时间的对数正态分布的情况下，比如为表示各种器官中的微脉管系统而描述的那些情况，对补给函数E(t)的实际形状的分析允许了估计流动分布的方差和斜度，从而可以评估组织病理。如上所述这，这些估计可按以下方法来进行：利用perf函数或perf函数的线性组合，对从实验线性回波数据导出的功率信号进行曲线拟合。

利用不同S形函数实现本发明

虽然上述perf函数正确地表示了以垂直方向上具有高斯型波束灵敏度的成像平面内的单个流动通道的回波能量再灌注，但是也可使用S形函数的不同的等价参数形式以实际可接受的精度来近似erf函数，所述等价参数形式例如是三角函数或具有多项式项的函数。这些等价的形式可从对erf函数的已知近似导出。可能的示例包括但不限于：

$\mathrm{erf}\left(q\right)\≅\mathrm{sign}\left(q\right)\{1-\frac{1}{1+a_1\left|q\right|+a_2{q}^2+a_3{\left|q\right|}^3+a_4{q}^4}\},$

其中a₁＝0.278393，a₂＝0.230389，a₃＝0.000972，a₄＝0.078108，对于q≥0，sign(q)＝1，对于q＜0，sign(q)＝-1，或者

$\mathrm{erf}\left(q\right)\≅\tanh (1.203\·q).$

注意上述式子还暗示了：

$\mathrm{perf}\left(q\right)\≅0.5[1+\mathrm{erf}\left(q\right)]=\frac{1}{1+e^{-2.406\·q}}\≡\mathrm{sigmoid}(2.406\·q).$

利用这些近似表达式，则可利用通过流动渡越时间或速度的概率分布(例如对数正态分布)加权的perf函数积分来拟合所观察到的再灌注功率函数，以便表示组织灌注。从而可根据对所观察到的数据的最佳拟合估计诸如平均流动速度、流动分布的方差或斜度这样的参数。这些参数可在回波描记图像内用户定界的回波能量区域(“感兴趣的区域”)内计算，该区域或者是由用户手动定界的，或者是通过回波描记图像中的边界描绘的已知方法来自动确定的，或者是根据其他器官解剖特征的。或者，可在各个图像元素内计算参数，以便将这些参数显示为空间参数制图(也称为参数成像)，而不要求用户定义感兴趣的区域，其中所述图像元素例如是数字视频像素(图片元素)、像素群组、局部成像分辨率所确定的个体斑点“颗粒”、像素颗粒群组或其他自动确定的图像区域。

估计渡越时间或速度的概率密度分布

根据上述方法获得的S形函数也可用于估计再灌注期间微泡的渡越时间或速度的概率密度分布(例如在不同毛细管中)。正如已经知道的，在健康组织(即无异常的组织)中，微泡的流动速度与血液相同；在此情况下，渡越时间或速度的分布构成对数正态分布。因此，只要通过将所估计的分布与对数正态分布相比较，就可将其用于检测生理异常。

已开发出了计算概率密度分布函数的新方法，在以下情况下该方法用在所提出的方法中尤其有利：在所述具有S形特性的函数的参数的至少一个值与局部组织灌注值或属性之间建立了对应关系。但是下述新方法对于在其他的灌注分析上下文(特征例如在于回波功率作为时间的函数而线性增大或指数衰减)中提供概率密度分布函数计算也是有价值的，所述其他的灌注分析上下文例如是UCA的丸剂施用的情形。

通过小波分析实现本发明

可以在考虑到再灌注回波数据是各个perf函数的未知之和的结果、而不是产生于渡越时间的对数正态分布的情况下来分析该数据。在这种情况下，发现应用在灌注函数的二阶导数上或灌注函数的平滑或低通滤波后的版本上的连续小波分解提供了对流动渡越时间的概率分布的有用估计。这里，用于小波比例和延迟系数估计的“母小波(mother wavelet)”(或模型小波)无疑可选择为速度为v_mw的perf函数的二阶导数。这种函数满足小波定义的所有要求。在这里的气泡破坏区域的UCA补给的上下文中，在数字形式上，连续小波变换的比例sc和延迟del系数取为：

$\mathrm{sc}=\frac{2\mathrm{\τ}\·v_{\mathrm{mw}}}{D\·\mathrm{dt}}$ $\mathrm{del}=\frac{\mathrm{\τ}}{\mathrm{dt}},$

其中dt表示采样时间间隔。

图9示出通过对再灌注曲线的二阶导数使用连续小波变换来估计渡越时间的任意分布的过程。要分析的是图9a的理论再灌注示例，它实际上是由均值分别以0.56和1.35秒为中心的渡越时间的双峰型分布给出的。此函数的二阶导数按相对幅度比例在图9b中示出。选中来执行小波分解的母小波在图9c中示出，它实际上被取为初等perf函数的二阶导数。在图9d中，将小波分解的结果与用于生成图9a的再灌注函数的原始的渡越时间双峰型分布相比较，并发现二者具有可接受的良好吻合度，这证实了小波分析在本发明的流动分布估计的上下文中的实用性。利用付立叶、Radon、希尔伯特、Z-或任何其他积分变换可执行类似的再灌注函数分析以估计渡越时间分布，在每种情况下，都基于以下认识：灌注渡越时间的个体值生成个体的功率贡献，该贡献是时间的函数，具有此文献中所公开的perf函数的类型。

对于利用上述方法来拟合观察到的再灌注数据，有以下替换方案：在某些情况下，在执行曲线拟合之前向实验数据应用任何已知形式的平滑或时间平均，例如移动平均、低通滤波、中值滤波等，或者这些平滑或时间平均的组合，将会是有利的。另一种替换方案可以是利用本发明中公开的S形拟合函数对“非线性”回波数据(例如对数压缩的数据)执行曲线拟合，其中所述S形拟合函数被与导致回波信号的非线性响应的过程相同的过程(例如S形函数的对数压缩)而修改。

利用神经网络分析重建流动速度的概率密度分布

在本发明的不同的实施例中，再灌注回波数据改为由单步骤或多步骤过程所分析。

分析的起始点是上述回波功率补给信号E(t)。

为清晰起见重复有关公式，垂直方向上的脉冲-回波声灵敏度模式PE(y)可由高斯型G(y)近似：

$\mathrm{PE}\left(y\right)\≅G\left(y\right):=e^{-{(1.94\·y)}^2\·K^2}---\left(1.1\right),$

有利地，补给函数E(t)被表示为：

E(t)＝perf(t)，(1.2)

其中perf(t)：＝O+A·φ(1.94·Kv(t-τ))，(1.3)

φ(q)＝0.5·(1+erf(q))

并且 $\mathrm{erf}\left(q\right)=\frac{1}{\sqrt{\mathrm{\π}}}\underset{0}{\overset{q}{\∫}}{e}^{-p^2}\mathrm{dp}$

O是偏移因子，A是幅度因子， $\mathrm{\τ}=\frac{D}{2v}$ (1.4)是微气泡从受破坏区域的边缘行进到中央所需的渡越时间。

在微气泡速度矢量v与y方向之间的角度θ_v不再为零而是任意的不等于90°的值的情况下，perf(t)的定义(1.3)变为

perf(t)：＝O+A·φ(1.94·K′v(t-τ))，(1.5)

其中K′：＝K·cosθ_v并且

$\mathrm{\τ}=\frac{D^{\′}}{2v},$

其中 $D^{\′}:=\frac{D}{\cos {\mathrm{\θ}}_v}.$

令v_y为v沿着y方向的分量，即

v_y＝v·cosθ_v；

则

perf(t)＝O+A·φ(1.94·K·v_y(t-τ_y))，(1.6)

其中 ${\mathrm{\τ}}_y=\frac{D}{{2v}_y}.$

注意(1.6)等价于定义(1.3)，其中v_y取代了v，τ_y取代了τ。Perf(t)(1.6)可以单独就v_y而言或单独就τ_y而言更恰当地写成：

${\mathrm{perf}}_{v_y}\left(t\right):=O+A\·\mathrm{\φ}(1.94\·K\·v_y(t-\frac{D}{{2v}_y})),---\left(1.7\right)$

${\mathrm{perf}}_{{\mathrm{\τ}}_y}\left(t\right):=O+A\·\mathrm{\φ}(1.94\·K\frac{D}{2{\mathrm{\τ}}_y}(t-{\mathrm{\τ}}_y)).---\left(1.8\right)$

在感兴趣的区域中有许多毛细管，每个毛细管的特征在于有一个且仅有一个再灌注速度。从而每个毛细管的补给函数由(1.7)或(1.8)给出(个体perf函数)：

$E\left(t\right)={\mathrm{perf}}_{v_y}\left(t\right)$

或者等价地

$E\left(t\right)={\mathrm{perf}}_{{\mathrm{\τ}}_y}\left(t\right).$

如果P(V_y)是感兴趣的区域中速度沿y方向的分量的概率密度分布，则在这种区域中的回波功率补给函数E(t)可表达为由每个毛细管的影响的概率所加权的个体perf函数的组合：

$E\left(t\right)=\∫P\left(v_y\right)\·{\mathrm{perf}}_{v_y}\left(t\right){\mathrm{dv}}_y.---\left(1.9\right)$

(1.9)的离散版本为：

$E\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}P\left(v_{i_y}\right)\·(v_{{i+1}_y}-v_{i_y})\·{\mathrm{perf}}_{v_{i_y}}\left(t\right).---\left(1.10\right)$

概率密度分布P(v)满足归一化属性：

∫P(v_y)dv_y＝1，(1.11)

以使得

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}P\left(v_{i_y}\right)\·(v_{{i+1}_y}-v_{i_y})\≅1.---\left(1.12\right)$

就τ_y而言：

$E\left(t\right)=\∫P\left({\mathrm{\τ}}_y\right)\·{\mathrm{perf}}_{{\mathrm{\τ}}_y}\left(t\right)d{\mathrm{\τ}}_y,$

其离散版本为：

$E\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}P\left({\mathrm{\τ}}_{i_y}\right)\·({\mathrm{\τ}}_{{i+1}_y}-{\mathrm{\τ}}_{i_y})\·{\mathrm{perf}}_{{\mathrm{\τ}}_{i_y}}\left(t\right),---\left(1.13\right)$

其中是P(τ_y)是感兴趣的区域中τ_y的概率密度分布；从而，类似于(1.11)和(1.12)，

∫P(τ_y)dτ_y＝1

并且

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}P\left({\mathrm{\τ}}_{i_y}\right)\·({\mathrm{\τ}}_{{i+1}_y}-{\mathrm{\τ}}_{i_y})\≅1.---\left(1.14\right)$

为了简化符号，以下将省略下标Y；从而将用符号v和τ取代v_y和τ_y，它们将分别被称为速度和渡越时间。

当微气泡以与血液流动相同的速度流动时(表现健康组织(无异常)的特征)，P(v)是lognormal概率密度函数：

$P\left(v\right)=\frac{e^{\frac{-{(\ln \left(v\right)-m)}^2}{2s^2}}}{\mathrm{vs}\sqrt{2\mathrm{\π}}},$

其中m和s分别是v的自然对数的均值和标准差。

类似地，就τ而言，健康组织的特征在于

$P\left(\mathrm{\τ}\right)=\frac{e^{\frac{-{(\ln \left(\mathrm{\τ}\right)-m)}^2}{{2s}^2}}}{\mathrm{\τs}\sqrt{2\mathrm{\π}}},$

其中m和s分别是τ的自然对数的均值和标准差。

图10示出对数正态情形的补给函数的示例。为构造此曲线使用了公式(1.10)，其中：

-对数正态概率密度分布参数：m＝1.5，s＝0.45；

-速度矢量：

v₁＝v_min·b⁰，...，v_n＝v_min·b^n-1，(1.15)

其中 $b=e^{\frac{\ln (v_{\min }/v_{\max })}{n-1}},$ v_min＝1，v_max＝25，n＝32。

方程(1.10)和(1.13)将补给函数表达为(个体)perf函数的线性组合。实际回波功率补给信号可表达为

E′(t)＝E(t)+N(t)，

其中E(t)是(理论)补给函数，N(t)表示噪声。

从而在补给阶段的时刻t₁，...，t_j所记录的信号向量可写作

E′：＝(E′(t₁)，...，y(t_j))＝(E(t₁)，...，E(t_j))+(N(t₁)，...，N(t_j)).(1.16)

目标是在不假设关于P(v)或P(τ)的形式的任何信息的情况下，从E’中任取i，找到P(v_i)的良好近似，或者等价地说是找到P(τ_i)的良好近似。

在下文中，问题陈述同样会按照v和τ来处理；尤其是会同样地描述P(v)或P(τ)的重建。已开发了一种新的单步骤或多步骤方法，有利地，可通过将

优化工作箱用于第一步骤，将

神经网络工具箱用于第二和更多的步骤，来实现此方法。

第一步骤的实现

令E’由(1.16)定义，并令sum_perf(p，t)为由以下线性组合定义的函数：

$\mathrm{sum}\_\mathrm{perf}(p,t)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i\·{\mathrm{perf}}_{v_i}\left(t\right),$

其中p＝(p₁，...，p_n)是加权因子向量，n是为进行分析所考虑的流动速度数目。令sum_perf(p)为由以下式子定义的向量

sum_perf(p)：＝(sum_perf(p，t₁)，...，sum_perf(p，t_j))..(1.17)

如果sum_perf被视为p的函数，则p的函数f(误差函数)可定义为

$f\left(p\right)=\underset{k=1}{\overset{j}{\mathrm{\Σ}}}|\mathrm{sum}\_\mathrm{perf}(p,t_k)-E^{\′}\left(t_k\right)|.$

从而目标是找到受约束非线性多变量函数的最小值

在p_i≥0 $\∀i=1...n$ 的条件下的

(1.18)

实际上，如果可以找到受约束函数(1.18)的解p_min，则可以通过以下简单规则找到P(v_i)的近似：

$P\left(v_i\right)\≅\frac{{\left(p_{\min }\right)}_i}{v_{i+1}-v_i},\∀i=1...n.---\left(1.19\right)$

此外补给函数的近似可为：

$E\left(t\right)\≅\mathrm{sum}\_\mathrm{perf}(p_{\min },t).---\left(1.20\right)$

编程语言包含称为优化工具箱的工具箱，这种工具箱对于解决最小化问题尤其有用。尤其是

函数fmincon从初始估计开始找到几个变量的标量函数的受约束最小值。fmincon函数使用一种子空间可信区域方法，该方法基于以下文献中描述的内部反射牛顿方法：Coleman，T.F.和Y.Li，“An Interior，Trust Region Approachfor Nonlinear Minimization Subject to Bounds(受限非线性最小化的内部可信区域方法)”，SIAM Journal on Optimization，Vol.6，pp.418-445，1996以及Coleman，T.F.和Y.Li，“On the Convergence ofReflective Newton Methods for Large-Scale Nonlinear MinimizationSubject to Bounds(大型受限非线性最小化的反射式牛顿方法的收敛)”，Mathematical Programming，Vol.67，Number 2，pp.189-224，1994。每次迭代涉及利用预处理共轭梯度(PCG)方法求大型线性系统的近似解。

f的梯度的解析公式被提供给fmincon，以便梯度不是数值地计算的，算法变得更精确。选择赋予fmincon的初始估计是相当重要的。实际上(受约束)函数f具有许多局部极小值，这些局部极小值允许了找到补给函数E(t)的良好近似，但不是P(v_i)的良好近似。

在N(t_k)＝0 $\∀k=1,...,j$ 的情形下，通过使用以下初始估计获得非常好的结果：

· ${\left(p_{\mathrm{initial}}\right)}_i=\frac{1}{n},\∀i=1...n---\left(1.21\right)$

·(p_initial)_i＝(v_i+1-v_t)·(max(v₁，...，v_n)-min(v₁，...，v_n))，i＝1...n(1.22)

给定的感兴趣的间距中的速度值可根据其算术或几何级数来选择。

图11示出根据图10所示的补给函数E(t)的采样重建P(v)。应用了具有速度向量(1.15)(与用于构造图10的采样相同，即n＝32)和初始化(1.22)的算法。

图12示出双峰型概率密度分布情形下的重建

$P\left(v\right)=\frac{2}{3}\log \mathrm{normal}(v,m_1,s_1)+\frac{1}{3}\log \mathrm{normal}(v,m_2,s_2),$

其中m₁＝1.5，s₁＝0.45，m₂＝1.7·m₁， $s_2=\frac{s_1}{2}.$

为了理解在实际实验信号的情形下该算法是否有效，它被应用到以下述方式构造的采样：

E(t)＝E(t)+N(t)(1.23)

其中N(t)是白噪声时间信号，与瞬时信号相比其具有10％的能量幅度。

图13示出在对数正态P(v)情形下利用(1.23)构造的“有噪”采样，图14是从这个有噪采样中找到的对数正态P(v)的重建。

图14清楚示出在噪声相对较高的情形下可能损失精确度。因此这种“单步骤”方法对于具有低噪声的采样尤其有用。对于有噪采样，按下述方式以第二和第三步骤完成该方法并且减少第一步骤中考虑的速度数目n将会是有利的。

实现第二和第三步骤：(神经网络)

神经网络的基本对象是神经元。简单神经元的模型可描述如下：标量输入u被传输经过一个将其强度乘以标量权重p的连接，以形成积p·u，此积也是标量；然后p·u和偏置b的和是转移函数g的自变量，该函数产生标量输出a。通常转移函数是阶跃函数或反曲函数或恒定函数(id)。

在具有向量输入的神经元中，输入u是向量u＝(u_i，...，u_n)，p是权重向量；内积p·u被与偏置b相加，以形成转移函数g的自变量，该函数产生标量输出a，如图15所示。

两个或更多个如前所述的神经元可被组织在一层中，一个神经网络可包含一个或多个这种层。但是不必详细考虑这些更复杂的网络；实际上，使用具有向量输入的仅由一个神经元形成的网络就足够了。令net为此类型的网络，并令g：＝id为其转移函数。则net可用

函数来创建。

通过调整权重和偏置值可训练神经网络，以使得特定输入导致特定输出。从而，对于网络net的训练过程，需要输入向量u₁，...，u_n和相应的目标输出向量a＝(a₁，...a_n)。在训练期间，反复调整网络的权重和偏置，以最小化网络性能函数(net.performFcn)。用newlin创建的网络的默认性能函数是均方误差(mse)，此误差是网络输出w＝g(p·u+b)＝id(p·u+b)＝p·u+b和目标输出a之间的平均平方误差：

$\mathrm{mse}=\frac{1}{n}{\left|\right|w-a\left|\right|}^2.$

令u_i为由下式定义的列向量

$u_i\left(k\right):={\mathrm{perf}}_{v_i}\left(t_k\right),\∀i=1,...,n,\∀k=1,...,j;---\left(1.24\right)$

令a为由下式定义的向量

a：＝E′..(1.25)

由于

$p\·u=\mathrm{sum}\_\mathrm{perf}(p,t_k),\∀k=1,...,j,$

于是

$\mathrm{mse}=\frac{1}{j}\underset{k=1}{\overset{j}{\mathrm{\Σ}}}{(\mathrm{sum}\_\mathrm{perf}(p,t_k)+b\left(t_k\right)-E^{\′}\left(t_k\right))}^2,$

其中sum_perf(p)由(1.17)定义。

如上所述，在训练期间，网络的权重和偏置被反复调整以最小化mse：如果偏置值被初始化为零并且在所有训练过程期间它都被保持为非常接近零，则此最小化问题与上述第一步骤中的问题(1.18)非常类似。

除了偏置值问题外，还必须解决其他问题：如何限制权重为正值或零，以及如何选择权重的初始估计。采用了以下方法：

1.偏置值：偏置值被初始化为0，并且每50或100次迭代它就被重置回0；

2.初始权重估计：利用均匀分布的只有8个速度的向量 $v_{\mathrm{eight}}=(v_{{\mathrm{eight}}_1},...{,v}_{{\mathrm{eight}}_8})$ 来应用前一节中描述的优化算法；通过这种方式，通过(1.19)，为分析的第一步骤确定了

的8个近似。通过利用立方平滑样条(函数csaps)在速度域拟合这8个值，通过在同样是均匀分布的v₁，...，v_n处对此样条估值，最后通过用v_i+1-v_i乘以每个值，则为分析的第二步骤获得了n个初始权重估计，其中n可以是比为第一步骤所选择的值大的任何值(例如在图16-18的情形中n＝32)；一般而言，已观察到在分析第一步骤中要考虑的速度的初始数目是从4到16，优选为从6到10，而对于第二步骤，要考虑的速度的数目为从8到64，优选为从16到48；

3.权重的非负约束：每50或100次迭代所有负权重由零替换。

通过使用这些规则和训练

函数traingdx，获得了非常好的结果。traingdx根据批梯度下降动量和自适应学习速率更新权重和偏置值。

·梯度下降：梯度下降算法在性能函数mse减小最快的方向上更新网络权重和偏置，该方向也就是与梯度相反的方向。此算法的一次迭代可写为

W_i+1＝W_I-α_Ig_I”

其中是W_I当前权重或当前偏置的向量，g_I是当前梯度，a_I是学习速率。

·批模式：在批模式中，只有当整个训练集合已被应用到网络之后，网络的权重和偏置才被更新；

·随动量下降的梯度：动量是可添加到梯度下降算法的技术，它允许网络只响应于局部梯度，还响应于误差平面中最近的趋势。动量起的作用与低通滤波器类似，它允许网络忽略误差平面中的小特征；

·自适应学习速率：如果在训练过程期间允许改变学习速率，则可改进梯度下降算法的性能。

此外，为了降低噪声的负面影响，在应用新算法之前，向有噪向量E’应用了中值滤波器(

函数medfilt1)。

结果在下文中给出。

在描述结果之前，有必要描述仅每50或100次迭代就用零替换所有的负权重和当前偏置的优点以及描述停止标准的基本原理。

在训练迭代期间，性能函数mse的值减小；但是当算法被停止以便用零替换所有的负权重和当前偏置时，训练过程所找到的值被修改，因此mse增大。从而迭代从性能函数的较高值再次开始。

如果停止之前的迭代次数J太小，则mse在增大之前减小得不够，从而算法可能不那么有效。

根据经验已观察到在图16-18所示的特定示例中，为获得最佳结果，J优选地应为50或更大，一般而言，J的值最少为10较好，至少为25更佳，至少为50则更好，一直到例如约200，通常为100。

对于所应用的停止标准，所述停止在mse的相对变化下降到指定容限tol之下时介入。注意mse的相对变化不是在每次迭代之后测量的，而是在负权重和当前偏置被零替换时测量的。于是迭代总次数始终是J的倍数。

虽然已就v而言描述了第一和第二步骤，但是也可就τ而言来描述它们。唯一的差别时使用了τ的向量而不是v的向量来构造采样E’，并利用这些算法来分析它。

图16、17和18分别示出在n＝32和J＝25，50，100情况下，从相同的有噪E’重建对数正态P(τ)(m＝-0.2191，s＝0.45)。在任何情况下，算法都由于mse的相对变化下降到容限tol＝0.01之下而停止，但是在J＝25的情况下mse的最终绝对值和在J＝50，100的情况下mse的最终绝对值之间有可感测的差异。表1示出mse的初始值和最终值。很明显在特定示例中，选择J＝25不允许mse减小到足以良好地重建P(τ)。

表1mse的初始值和最终值

	初始mse	最终mse
			J＝25	111.434	12.5768
J＝50	111.434	0.5925
			J＝100	111.434	0.4729

当通过方程(1.23)构造合成数据E’时，使用了偏移因子O的给定值和幅度因子A的给定值。在此情况下，在应用算法时可使用这些已知值。否则，在处理实际数据时，可假设O＝0，但幅度值未知。通过假设E’的渐近值可给出A的一级近似，并将此值用于分析。但是，这种方式可能不是处理该问题的完全符合要求的方式；从而，为了提供在实际情况下该方法的效率，该方法可按以下方式来修改：

令

为perf_τ(t)的修改，其由下式定义

以使得补给函数可表达为

从而可用新的

取代perf_τ(t)来应用该算法。从现在开始算法的新版本将被称为“最终算法”。与旧版本唯一的差别是找到了A·P(τ_i)而不是P(τ_i)的近似[A·P(τ_i)]_est。

在实际幅度已知的合成情况下，可将[A·P(τ_i)]_est与实际A·P(τ_i)相比较。此外，由于概率密度函数满足方程(1.14)，因此它遵循

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}A\·P\left({\mathrm{\τ}}_i\right)\·({\mathrm{\τ}}_{i+1}-{\mathrm{\τ}}_i)\≅A.$

从而可找到A的近似如下：

$A_{\mathrm{esi}}:=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[A\·P\left({\mathrm{\τ}}_i\right)]}_{\mathrm{est}}\·({\mathrm{\τ}}_{i+1}-{\mathrm{\τ}}_i).$

图19、20分别示出在J＝50，100的情况下通过最终算法从图16、17、18的同样的E’重建P(v)；tol为0.01。幅度的实际值为A＝100。

如上所述，目标是在不假设关于P(v)的形式的任何信息的情况下找到P(v)的良好近似，或者等价地说找到的P(τ)良好近似。利用上述最终算法找到的几个概率密度分布的重建在下文中描述，其中始终有10％的白噪声。

图21示出在均值μ＝12.5、方差σ²＝11，即均值A·P(v)；P(v)为高斯型；J＝50；tol＝0.01的情况下高斯型P(v)的重建

图22表示在P(v)为高斯型；J＝100；tol＝0.01的情况下A·P(v)的重建。

图23表示在P(v)为梯形；J＝50；tol＝0.01的情况下梯形P(v)的重建，即A·P(v)的重建。

图24表示在P(v)为梯形；J＝100；tol＝0.01的情况下A·P(v)的重建。

图25表示由下式定义的瑞利型P(v)的重建：

$P\left(v\right)=\frac{2v}{b}{e}^{-{\left(\frac{v}{b}\right)}^2},$

其中b＝6。即这是在P(v)为瑞利型；J＝50；tol＝0.01的情况下A·P(v)的重建。

图26是在P(v)为瑞利型；J＝100；tol＝0.01的情况下A·P(v)的重建。

图27是由下式定义的双峰型P(v)的重建：

$P\left(v\right)=\frac{1}{2}\log \mathrm{normal}(v,m_1,s_1)+\frac{1}{2}\log \mathrm{normal}(v,m_2,s_2),$

其中m₁＝1.5，s₁＝0.45，m₂＝1.7·m₁， $s_2=\frac{s_1}{2}.$

图27是在P(v)为双峰型；J＝50；tol＝0.01的情况下A·P(v)的重建。

图28是在P(v)为双峰型；J＝100；tol＝0.01的情况下A·P(v)的重建。

图27所示的重建不如图28的好；实际上已观察到，在比对数正态或高斯型分布更复杂的概率密度分布的情况下，通过使用较小的停止标准容限，例如tol＝0.001，可改进J＝50的算法的性能：

图29是在P(v)为双峰型；N＝50；tol＝0.001的情况下A·P(v)的重建。

图30是由下式定义的双峰型P(τ)的重建：

$P\left(\mathrm{\τ}\right)=\frac{3}{5}\log \mathrm{normal}(\mathrm{\τ},m_1,s_1)+\frac{2}{5}\log \mathrm{normal}(\mathrm{\τ},m_2,s_2),$

其中m₁＝-0.2191，s₁＝0.45，m₂＝-3·m₁， $s_2=s_1^2.$

从而图30是在P(τ)为双峰型；J＝50；tol＝0.01的情况下A·P(τ)的重建。

图31是在P(τ)为双峰型；J＝100；tol＝0.01的情况下A·P(τ)的重建。

可观察到在图30的特定示例中，J＝50的算法不允许精确重建此双峰型分布；如同先前的情况那样，如图32所示，通过使用较小的停止标准容限，可改进算法的性能，该图是在P(τ)为双峰型；J＝50；tol＝0.001的情况下A·P(τ)的重建。

在此情况下，小容限对于改进J＝50的算法非常有用，但是一般而言该问题是非常微妙的：在简单分布情况下使用tol＝0.001可能导致精确度损失，所述简单分布优选以tol＝0.01来重建。例如，图25示出通过用J＝50和tol＝0.01获得的瑞利型分布的非常好的重建。图33示出通过用J＝50和tol＝0.001从相同的E’获得的相同瑞利型分布的重建。

与J＝50的算法的情况不同，如果使用J＝100，则对于任何种类的分布来说容限tol＝0.001都太小了。实际上也观察到了双峰型分布的重建精确度的恶化。例如，对于图30-32的双峰型，通过使用J＝100和tol＝0.001，获得了图34中重建，它是P(τ)为双峰型；J＝100；tol＝0.001的情况下A·P(τ)的重建。

从而实现了以下目的：在不假设关于P(v)的形式的任何信息的情况下，对于任何种类的分布，找到P(v)的良好近似，或者等价地说，找到P(τ)的良好近似。从而可以通过将所找到的重建与作为健康组织的特征的对数正态分布相比较，来了解感兴趣的区域中是否有异常。

由于上述最终算法，在具有10％的白噪声的合成数据的情况下，获得了良好结果。

特别地，利用这些新方法，只要通过将所找到的重建与作为健康组织特征的速度分布相比较，就可以检测生理异常。利用这些新方法，不必区分补给信号。特别地，这些新方法即使对于具有相当大的噪声的数据也可提供非常好的结果。

修改

可以对所描述的方法和系统实施例进行许多修改，并且所述方法和系统除用于确定所述的灌注值或属性外，也可用于确定许多其他灌注值或属性。

此外，以上所述的采用曲线拟合然后进行神经网络训练的多步骤分析可应用到除S形函数外的其他灌注函数，这种灌注函数是在UCA滴注下在破坏帧序列之后获得的。一个示例是到达被研究区域的UCA的丸剂注射。在存在多个流动渡越时间的情况下，期望的回波功率丸剂函数E(t)可被视为回波功率函数E_i(t)之和，每个E_i(t)由其流动渡越时间τ_i所确定。这种基本丸剂函数E_i(t)可以是以下类型的：

$E_i\left(t\right)=e\frac{t}{{\mathrm{\τ}}_i}{e}^{\frac{-t}{{\mathrm{\τ}}_i}}.$

在图35所示的示例中，考虑了流动渡越时间C(τ_i)的对数正态分布，其中m和s分别等于0.313和0.4。以下描述渡越时间的四个部分子分布，分别称为C1_part(τ_i)至C4_part(τ_i)。考虑由部分子分布所确定的四个回波功率函数E1_part(t)至E4_part(t)之和所产生的回波功率函数E(t)，可应用分析方法以提供对对数正态渡越时间分布C(τ)的估计，如图36所示。

因此，根据本发明，分析回波功率函数允许了对原始对数正态分布的估计，如图37所示。

超声成像系统

图38的框图示出了可用于根据本发明的方法中的典型医学超声成像系统的主要元件。在中央处理单元1和中央计时电路2的的控制下，发射波束形成器和脉冲发生器3被用于通过适当的Tx/Rx(发射/接收)复用器4向多元素超声探头5施加激励信号。当在传播介质内发生超声反射时，回波信号通过前置放大器和前置时间增益补偿(TGC)6而被处理。通常，这些接收到的信号接下来被A/D转换器从模拟电压转换为数字值，并且通过接收波束形成器7被组合成聚焦的接收信号。通过经由数字滤波器8(例如带通滤波器)和其他信号调节器9(例如波束形成后TGC)进行的处理，数字射频(rf)信号10被解调器11所处理以生成与回波包络幅度成比例的信号，并且在被写入扫描转换器13之前被进一步非线性地处理(例如通过对数压缩预处理器12)，以便考虑到探头几何形状和波束扫描序列的因素。扫描转换后的信号被任选地再次压缩(通过后处理器14)，然后被转换成视频标准信号。然后视频信号可被显示在视频显示单元15上和/或存储在存储单元16上，或者采取模拟形式(在此情况下存储单元16将会是视频记录设备)，或者采取数字形式，例如计算机文件(在此情况下存储单元16可以是硬盘)。然后这样获得的图像被根据本发明处理，如下所述。

当使用这种超声成像系统来实现根据本发明的方法时，用于根据本发明的再灌注动力特性分析的信号优选地是在幅度敏感模式(例如谐波B-模式、脉冲反转或功率调制多脉冲模式或其他编码激励或解调单脉冲或多脉冲模式)下根据对比度特定检测方案获得的，所述检测方案例如是利用微泡的非线性声响应的方案。或者，用于再灌注动力特性分析的信号可以在多普勒模式(例如谐波功率多普勒、去相关成像，或其他利用响应于连续激励脉冲而发生的响应的变化的模式，所述变化是由于微泡的位置或声响应的变化而引起的)下根据对比度特定检测方案来获得。

用于流动速度估计的系统/方法

图39的框图示出在参数流动成像的一个特定示例中，根据本发明的用于实现流动估计方法的系统的主要功能元件。从包含B-模式图像(包括破坏之后UCA的补给)序列的计算机文件21(例如来自图38的数字文件存储装置16)开始，基于所使用的超声设备特定的对数压缩或其他非线性增益数据23，图像值被适当的线性化22处理，以产生与局部UCA浓度成比例的数据。然后，根据本发明，时间序列数据被最佳拟合方法24调整为描述期望的灌注定律的参数方程，如上所述。这些参数的物理解释25要求来自超声设备的输入数据26，例如Tx/Rx波束灵敏度数据以及对气泡破坏区域的厚度D的估计。最佳拟合参数允许了计算平均流动速度数据、平均渡越时间分布估计或其他流动特性，然后这些数据、估计或特性可通过已知的灰度级或伪彩色编码而被映射到扫描转换存储器27中，并且可选地被与B-模式或2维多普勒数据重迭，并且被显示在视频显示单元28上或以数字或模拟形式存储(未示出)。

根据本发明的系统的功能元件可被存储在一个程序上，该程序在被加载到计算机(CPU 1或与单元11-15相关联的计算)中时，可操作以用于在自动或在手动控制下将和再灌注期间的局部制剂浓度成比例的信号与相应的具有S形特性的时间函数相关。所使用的数学函数例如可用Matlat曲线拟合工具箱(Mathworks，Natrick，MA，美国)来实现。此程序优选地还包含用于处理所记录的回波信号的幅度以使所述信号与局部微泡浓度成比例的装置，并且在被加载到系统的计算机中时可操作以自动地或在手动控制下实现所述处理，然后处理后的信号被校准到或联系到S形参数函数。该程序优选地还包括用于实现以下功能的装置：自动地或在手动控制下，在具有S形特性的函数的至少一个参数的至少一个值与至少一个局部组织灌注值或属性之间建立对应关系。此拟合例如可利用上述可信区域方法来计算。这两个主函数都可以自动或半自动地实现，并且具有允许用户控制的适当用户接口。

对选中的感兴趣区域内的组织灌注的估计可通过用户执行研究来确定，或通过结合在系统程序中的自动装置基于自动边界限定或其他解剖学器官特征来确定。当在个体二维图像元素或象素内估计灌注值，这些灌注值可以按参数图像的形式被显示在查看监视器(15)上，其中象素灰度级强度或色彩渲染是根据局部估计的个体参数的值或其组合来编码的。系统程序也可估计二维图像元素或象素群组内的灌注值，这些群组是通过局部成像分辨率以以下方式来确定的：对于每个象素群组基本上只获得每个象素的一个值。通常，当在象素群组内估计灌注值时，这些群组基本上可以由回波描记仪器的局部斑点模式(自动地)确定。象素群组的大小可根据局部回波描记图像的二维空间付立叶分析确定，其范围与局部最大有效空间频率成反比，这些都是由系统程序自动计算的，或在用户控制下计算的。该程序还可根据最后的破坏超声脉冲与S形函数达到以下幅度的时刻之间的时延来估计流动参数，该幅度等价于紧挨施加破坏超声脉冲之前的浓度值的一半相等的制剂浓度。

系统程序还可计算平均局部血液速度，该速度被表达为以下比率：微泡被破坏处的切片的厚度的一半除以到达1/2最大浓度的时间。为此目的，还可在与成像平面垂直的方向上，从发射波束垂直宽度，提供对微泡被破坏处的切片的厚度值的近似。

对于破坏脉冲之后的大时间值，系统程序可将稳态制剂浓度级别计算为S形函数拟合的渐近值，或者此浓度级别可取为紧挨破坏脉冲之前的值。

系统程序可将S形函数计算为：累积正态概率分布函数；或反曲函数；或累积正态概率分布函数的多项式近似；或包含多项式形式的有限展开的参数表达式；或者累积对数正态概率分布函数的参数表达式，这些都在上文中更详细描述。在一个实施例中，它将S形函数计算为：由流动速度或渡越时间的对数正态概率分布加权的累积正态概率分布函数的和或积分的参数表达式，其最佳拟合参数值代表器官灌注物理量(例如流动速度或渡越时间)及其均值、方差和斜度。

流动速度值可通过用再灌注函数的平滑后或滤波后版本估计到达半稳态回波功率幅度的时间来估计。

与局部浓度成比例的信号在被参数拟合方程调整之前，有利地通过低通滤波器或其他平滑函数被处理。与局部浓度成比例的信号不是被线性化，而是可以被用参数拟合函数来调整，所述参数拟合函数是通过与回波描记仪器中应用的那些具有相同性质的非线性函数修改过的。

正如以上参考图9a-d中所描述的，回波功率数据可通过小波分解方法来分析，以估计不同流动渡越时间或速度处的贡献的分布。为此，与局部浓度成比例的信号在被通过小波分解方法分析之前，优选地首先被低通滤波或平滑；它们也可在被小波分解方法分析之前被微分两次。在小波方法中，用于分解的母小波例如是用于描述单个流动值的再灌注函数的累积正态分布函数的二阶时间导数。

额外，或作为替换地，回波功率数据通过单步骤或多步骤过程而被分析，以估计不同流动渡越时间或速度下的贡献的分布。在第一步骤中，选择第一组流动渡越时间或速度，并且通过多个S形函数的线性组合与回波功率数据的最佳拟合来进行第一估计。然后对于第二组流动渡越时间或速度进行第二估计，其中利用第一估计作为定义第二组的基础。然后第二估计可用于提供进行第三估计的初始值集合。

第二估计例如可以利用立方样条外推来进行，而第三估计有利地是用神经网络分析来进行的；在此情况下，神经网络的偏置值和每个负权重优选地被周期性地重置为零(例如在对其进行调整的每50-100次迭代中)。第一估计可用相对小数量的S形函数来进行，通常最多用16个S形函数，优选为最多用8个S形函数。然后第二估计可用较大数量的流动渡越时间或速度来进行，例如用至少8个一组，优选为用至少16个流动渡越时间或速度。

与局部浓度成比例的信号也可从射频回波信号(即图38中的10处)导出，导出方式例如是通过取解调回波信号的平方。或者，与局部浓度成比例的信号是通过以下方式根据视频回波信号导出的：在对数压缩(图38中14处)之后，应用与对数压缩的效果相逆的反向对数压缩，然后取这些信号的平方。在一个实施例中，与局部浓度成比例的信号是通过对比度特定检测方案来获得的，该方案在幅度敏感模式(例如谐波B-模式，脉冲反相或功率调制多脉冲模式，或其他编码激励或解调单脉冲或多脉冲模式)中利用微泡的非线性声响应。

系统程序可在一个或多个模块中提供，这些模块被加载到图像系统的一个或多个计算机中，从而允许在用户的控制下自动或半自动地实现本发明的方法。所有用于实现所述过程步骤的数学函数都可以自动实现。

该程序通常在CD-ROM、DVD或任何其他计算机可读介质上提供；或者，该程序被预先加载到硬盘上，通过网络被发射到计算机，被广播，或者更一般地说被以任何其他可直接加载到计算机的工作存储器中的形式提供。但是，根据本发明的方法导致其本身能被用硬件结构(例如集成在半导体材料的芯片中)实现，或者用软件和硬件的组合来实现。

此外，应当注意上述程序适合于甚至以独立产品的形式被投放到市场上，以便用在现有的超声成像系统中；例如，该程序可用在系统中集成的计算机上，或者用在接收由单独的扫描器提供的回波信号序列(例如通过可移动盘、通信线路或无线连接)的外部计算机上。

示例

为了描述所公开的发明对于体内情形的实用性，利用充气型动脉内气球在小种猪的肾中产生人工闭合(狭窄)。所使用的造影剂是

(Bracco Imaging S.p.A.)，该造影剂被以1mL/分钟的恒定滴注速率静脉内施用。所使用的超声成像设备是来自ATL/Philips(Bothell，WA)的HDI-5000，其中在脉冲反相对比度特定B模式成像中使用了C5-2凸型探头。UCA破坏/补给图像序列被记录和传送到计算机，以供再灌注分析。利用允许灰度级线性化的视频密度计量程序在肾皮层中绘出两个感兴趣的区域(AOI)。第一个AOI-1是在正常灌注的肾皮层中绘出的，而另一个AOI-2是在皮层的灌注不足区域内绘出的。图40a和图40b中分别以时间函数形式示出应用气泡破坏帧之后最初10秒中AOI-1和AOI-2的线性化视频密度计量数据点。最佳拟合logperf能量函数E(t)被示为相应图上的实线。对于这两个AOI找到的最适当的最佳拟合参数为：对于AOI-1和AOI-2分别有A＝59.5和A＝0.62，τ＝1.25s和τ＝3.1s。对于两个区域的D值6和8mm，所确定平均横向流动速度v_y分别为2.4mm/s和0.97mm/s。相对体积流动为A和v_y之积，即对于AOI-1和AOI-2值分别为142.8和0.6，从而表明了两个皮层区域之间非常大的流动差别。

无疑，为了满足本地和特定要求，本领域的普通技术人员可向上述解决方案应用许多修改和变化，但是这些修改和变化都被包括在以下权利要求书所限定的本发明的保护范围内。

Claims (38)

1.一种非侵入性地量化活体组织中的灌注的方法，该方法包括以下步骤：

提供指示所述组织中的成像造影剂的补给的回波信号序列，

将一个S形参数时间函数与所述回波信号相关联，所述S形函数包括具有基本上恒定的初始值的初始部分，具有基本上恒定的最终值的最终部分，以及所述初始部分和所述最终部分之间的中央部分，在所述中央部分中所述S形函数从所述初始值单调地变化到所述最终值，其中所述S形函数在初始部分和最终部分中具有基本上为零的一阶导数；并且

在所述S形函数的至少一个参数的至少一个值与至少一个局部组织灌注值或属性之间建立对应关系，所述属性表示不同流动贡献的分布范围。

2.如权利要求1所述的方法，其中所述造影剂包括能够反射声能的微泡，通过以下步骤获得所述回波信号序列：

在超声成像装置的成像平面中施加第一超声脉冲，所述超声脉冲的声压高到足以导致该成像平面内存在的微泡的有效部分被破坏；

在所述成像平面中施加第二超声脉冲序列，所述第二超声脉冲的声压足够低以保持所述微泡的大部分；并且

在预定的后续时刻重复所述施加第二超声脉冲序列的步骤，并且记录源自所述成像平面的由所述第二超声脉冲产生的回波信号，以监视在所述后续时刻所述成像平面内的微泡的补给。

3.如权利要求2所述的方法，还包括在关联所述S形函数之前执行以下步骤：

处理所述回波信号，以使所述回波信号与所述微泡的局部浓度成比例，从而产生与所述成像平面内的任何位置处的造影剂浓度成比例的处理后的回波信号。

4.如权利要求3所述的方法，其中所述回波信号是以二维图像的形式获得的，所述二维图像与和所述回波信号的幅度成比例的图像元素相关联，并且所述处理包括所述回波信号的线性化。

5.如权利要求3所述的方法，其中所述回波信号的幅度正比于超声回波幅度，并且所述处理包括所述回波信号的线性化。

6.如权利要求5所述的方法，其中所述回波信号是射频信号或解调后的射频信号。

7.如权利要求1所述的方法，其中所述组织灌注值或属性是在选中的感兴趣的区域内被估计。

8.如权利要求1所述的方法，其中所述灌注值或属性是在各个二维图像元素内被估计的，然后被以参数图像的形式在观看监视器上显示为象素灰度级强度或色彩渲染，其中所述象素灰度级强度或色彩渲染是根据局部估计的各个参数值或其组合来编码的。

9.如权利要求1所述的方法，其中所述灌注值或属性是在二维图像元素的群组内被估计的，所述群组是通过局部成像分辨率以以下方式来确定的：对于每个图像元素群组基本上只获得每个象素的一个值。

10.如权利要求2所述的方法，其中所述灌注值或属性是在图像元素的群组内被估计的，所述图像元素群组的大小基本上是所述超声装置的回波描记仪器的局部斑点模式确定的。

11.如权利要求10所述的方法，其中所述图像元素群组的大小是根据局部回波描记图像的二维空间付立叶分析确定的，其范围与局部存在的最大有效空间频率成反比。

12.如权利要求1所述的方法，其中所述对应关系是与从下述值选出的至少一个灌注值建立的：平均渡越时间、平均速度、平均流量和灌注体积。

13.如权利要求1所述的方法，其中所述对应关系是与从下述属性选出的至少一个灌注属性建立的：血液流动模式、流动分布方差和流动斜度。

14.如权利要求2所述的方法，其中所述灌注值或属性是根据所述第一超声脉冲与所述S形函数达到下述值的时刻之间的时延来估计的：该值等价于一个与紧挨着施加所述第一超声脉冲之前存在的浓度值的一半相等的造影剂浓度。

15.如权利要求14所述的方法，其中所述灌注值是平均局部血液速度，该速度被表达为以下比率：所述微泡被破坏的切片的厚度的一半与到达1/2最大浓度的时间之比。

16.如权利要求15所述的方法，其中所述微泡被破坏的切片的厚度是在垂直于所述成像平面的方向上、由所施加的第一超声脉冲的宽度值来近似的。

17.如权利要求1所述的方法，其中所述S形函数是：累积正态概率分布函数，反曲函数，或累积正态概率分布函数的多项式近似。

18.如前述权利要求17所述的方法，其中所述S形函数是包含多项式形式的有限展开的参数表达式，或者累积对数正态概率分布函数的参数表达式。

19.如权利要求17所述的方法，其中所述S形函数是：由流动速度或渡越时间的对数正态概率分布加权的累积正态概率分布函数之和或积分的参数表达式，所述S形函数具有的最佳拟合参数值代表器官灌注的物理量。

20.如权利要求1所述的方法，其中在关联所述S形函数之前，所述回波信号被通过平滑函数处理。

21.如权利要求1所述的方法，其中所述回波信号被通过小波分解方法分析，以估计不同流动渡越时间或速度处的贡献分布。

22.如权利要求21所述的方法，其中与局部浓度成比例的回波信号在被通过所述小波分解方法分析之前被微分两次。

23.如权利要求22所述的方法，其中用于所述分解的母小波是用于描述单个流动值的S形函数的累积正态分布函数的二阶时间导数。

24.如权利要求1所述的方法，其中所述回波信号被分析以估计不同流动渡越时间或速度处的贡献分布，所述分析包括以下步骤：

选择第一组流动渡越时间或速度；并且

通过多个S形函数的线性组合与所述回波信号的最佳拟合来进行第一估计。

25.如权利要求24所述的方法，其中所述第一组流动渡越时间或速度是从4至16个S形函数的线性组合获得的。

26.如权利要求24所述的方法，其中所述分析还包括以下步骤：

进行第二估计以定义第二组流动渡越时间或速度，其中利用所述第一估计作为定义所述第二组的基础。

27.如权利要求26所述的方法，其中所述第二组流动渡越时间或速度是从8至64个S形函数的线性组合获得的。

28.如权利要求26所述的方法，其中所述第二估计被用于提供对所述第二组流动渡越时间或速度进行第三估计的初始值集合。

29.如权利要求26所述的方法，其中所述第二估计是利用立方样条外推来进行的。

30.如权利要求28所述的方法，其中所述第三估计是用神经网络分析来进行的。

31.如权利要求30所述的方法，其中所述神经网络由用于所述第二组流动渡越时间或速度的多个权重以及用于加权后的所述第二组流动渡越时间或速度的多个偏置值所定义，所述第三估计包括：

反复调整所述偏置值和所述权重；并且

周期性地将所述偏置值和每个负权重重置为零。

32.如权利要求31所述的方法，其中所述重置步骤是以等于10至200次迭代之间的周期执行的。

33.如权利要求24所述的方法，其中所述第一估计最多是用8个S形函数来进行的。

34.如权利要求26所述的方法，其中所述第二估计是用一组至少16个流动渡越时间或速度来进行的。

35.一种非侵入性地量化活体组织中的灌注的系统，该系统包括：用于提供指示所述组织中的成像造影剂的补给的回波信号序列的装置；用于将一个S形参数时间函数与所述回波信号相关联的装置，所述S形函数包括具有基本上恒定的初始值的初始部分，具有基本上恒定的最终值的最终部分，以及所述初始部分和所述最终部分之间的中央部分，在所述中央部分中所述S形函数从所述初始值单调地变化到所述最终值，其中所述S形函数在初始部分和最终部分中具有基本上为零的一阶导数；以及用于在所述S形函数的至少一个参数的至少一个值与至少一个局部组织灌注值或属性之间建立对应关系的装置，所述属性表示不同流动贡献的分布范围。

36.如权利要求35所述的系统，其中所述造影剂包括能够反射声能的微泡，所述用于提供回波信号序列的装置包括：用于在超声成像装置的成像平面中施加第一超声脉冲的装置，所述超声脉冲的声压高到足以导致该成像平面内存在的微泡的有效部分被破坏；用于在所述成像平面中施加第二超声脉冲的序列，所述第二超声脉冲的声压足够低到保持所述微泡的大部分；以及用于在预定的后续时刻记录源自所述成像平面的、由所述第二超声脉冲产生的回波信号以监视在所述后续时刻所述成像平面内的微泡的补给的装置。

37.如权利要求36所述的系统，还包括用于在关联所述S形函数之前处理所述回波信号以使所述回波信号与所述微泡的局部浓度成比例，从而产生与所述成像平面内的任何位置处的造影剂的浓度成比例的处理后的回波信号的装置。

38.一种用于如权利要求35至37中任何一项所述的系统中的设备，该设备包括用于输入所述回波信号的装置；所述用于将所述S形参数时间函数与所述回波信号相关联的装置；以及所述用于建立对应关系的装置。

Images