用于取得流动通道的尺度相关信息的方法与系统
技术领域
本发明涉及一种用于根据对流过未知尺度(dimension)的通道的流体进行的测量来确定通道尺度相关信息的方法与系统。本发明允许使用非侵入查询信号(non-invasive interrogating signal)来取得基本上实时的速度测量,并且进而使用这些测量来评估包含该流体的通道的特性,包括瞬时容积流速。
背景技术
在许多情况下,需要取得流动通道的尺度相关信息,而此时却不容易或不希望接近该流动通道或者以其他方式取得直接的尺度测度。当流动通道随时间变化时,取得这些信息尤其富有挑战性。一个例子是取得人体或其他患者的诸如升主动脉之类的血管(即动脉或静脉)的尺度相关信息。人们感兴趣的尺度相关信息可以是流动通道的尺度或者从尺度中推导出的或依赖于尺度的其他信息,诸如定量的流动速度信息(例如容量或质量流速),血管弹性/健康,每次心跳的输送量,射血微量(injection fraction)等等。
人们特别感兴趣的通道尺度相关特性为容积流速(volumetric flow rate,VFR)。流过通道的流体的容积流速(VFR)依赖于流体的速度以及通道的横截面积。当知道在可以被认为是瞬时的短时间段之间的这些值时,可以计算瞬时VFR。确定患者体中流体的VFR对于(例如)评定心脏功能是有用的。在这种应用中,流体为流过诸如主动脉之类的闭通道的液体——例如血液。不幸的是,各个患者之间的变化使人们在不进行测量以确定通道尺度的情况下无法知道患者体内的通道尺度。
当前,用来测量血液流量的主要技术采用时变标记,诸如染剂、冷却生理盐水或者用小型电热器加热的血液,通过导管注入心脏。然后,使用该标记的变化速度来估计总体血液流速。一般地,导管被导入患者的股动脉,并且穿过静脉系统进入并通过心脏进入肺动脉。该方法的一种变体是注入染剂而不用导管,但需要在心房附近注射,并且在诸如耳部之类的地方测量染剂浓度。不幸的是,这些时变标记方法依赖于标记稀释瞬态(marker dilutiontransient),并且提供整体流量数据,而不是通过通道的流速数据。另外,时变标记方法不提供瞬时VFR。
另外一组测量心脏流量的方法称为“菲克(Fick)方法”,其使用动脉与混合静脉血之间氧含量差异以及身体总耗氧量来计算心输出量。典型菲克方法使用动脉与肺动脉导管来测量氧含量。有关方法测量患者气道中氧与二氧化碳数据,以避免使用导管。然而,与时变标记方法类似,菲克方法测量总体心输出量,而不是流速。最后,诸如生物阻抗(bioimpedance)之类的其他方法也可以用来测量心输出量。但是,与菲克方法和时变标记方法类似,这些方法不测量流速,而只测量总体心输出量。
存在某些非侵入确定流体速度的方法,其中多普乐(Doppler)方法是在血流量分析中最常见的。然而,这些方法只测量流体速度,并且如果需要定量流速信息则需要其他测量以提供通道尺度。这些通道尺度测量有时基于与成像或时间延迟相关的方法取得。也可以使用X射线与不透射线的染剂,但是结果很耗时,并且不容易理解。另外,来自超声信号的回声的时间延迟可以用来提供对通道尺度的估计,但是要取得精确的结果需要将超声换能器置于通道的近旁,但除非将这些设备注入患者体内,否则这是不可能的。因此,非侵入确定通道流体VFR存在显著缺陷。
由于测量人为现象以及有关处理,现有的流动通道测度还可能有问题、不确定或不准确。例如,由于从测量人为现象中分离感兴趣信号的困难性,升主动脉的超声测量是有问题的。这些人为现象可以与来自该主动脉之外的生理组织或者其他可能隐藏感兴趣信号的通道的回声、通道内噪声源或者其他干扰成分有关。有些方法已经试图通过将探针置于该主动脉近旁来最小化某些人为现象分量,因此产生了发病率折衷(trade-off)。其他方法试图通过采用非常窄的信号来最小化某些人为现象成分,但由此将对主动脉或者所考虑的其他通道的对准(targeting)复杂化。
发明内容
考虑到以上情况,本发明的一个大的方面是提供一种用于根据定量流速信息(例如速度测量)来取得流动通道的尺度相关信息的方法与系统,所述流动通道包括患者的流动通道,诸如血管。相关的一个方面涉及确定诸如流体流过通道的(多个)容积流速(VFR)之类地定量流速。本发明的另一方面是提供一种用于非侵入确定流体(例如血液)流过活体患者的VFR的方法与系统。本发明的另一方面是提供一种用于确定流体通道尺度并将该流体通道尺度数据与流体速度数据相结合以确定瞬时VFR数据的非侵入方法与系统。具体地,本发明的一个方面是使用同一信号来测量流体速度数据并确定流体通道尺度。本发明的另一方面是非侵入地取得诸如血管等流动通道的基本实时的速度相关以及尺度信息,从而能够得到瞬时定量流速信息,包括对于尺度随时间变化的通道。本发明的另一方面是准确地处理与有关于流动通道的测量相关的人为现象。本发明的另一方面是提供对于根据其他方法在计算上难以处理的公式的快速计算。
在医疗情况下,本发明允许确定流动通道的尺度相关信息,而不用将探头引入患者体内,从而减少了痛苦与不适,并且还降低了诸如可能感染等其他问题的可能性。另外,本发明还能够进行可以用来确定瞬时VRF的测量,从而使门诊医生能够检查时间变化,并且将VFR测量与诸如心电图(EKG)测量之类的其他生命体征测量相结合。
根据本发明的一个方面,患者体内生理组织的流动的流动特性信息用来取得与流动通道的尺度相关的处理后信息。本发明的这一方面可以实现为一过程,其至少部分由处理器进行,并且本发明的这一方面可以实现于(例如)软件产品或者其他逻辑,用于执行此类逻辑的处理单元,或者用于进行相关医疗处理的系统。
在一个实现中,流动特性信息为定量流动特性信息,即,与流动速度或者其导出数有关的信息。例如,可以在相对于通道横截面的一个或更多个点处侵入地或者非侵入地测量流动速度。这些流动速度测量可以用来计算该通道横截面的导出信息,诸如平均流动速度,以及与相对于诸如半径等通道尺度的速度变化有关的信息。这些流动速度测量还可以用来计算表征瞬时速度剖面的其他参数。可以重复这些测量以提供与速度剖面的时间变化有关的导出信息。因此,应该理解流动特性信息可以各种形式提供,或者可以由不同的参数表征。另外,流动特性信息可以基于针对所考虑的流动通道而进行的单一或多个测量。另外,这些流动特性信息可以用于计算的多个步骤,例如,速度剖面信息可以用来导出诸如尺度相关信息等第一处理后信息,并且该第一处理后信息可以与平均速度或者其他流动特性信息相结合,以导出第二处理后信息。
根据所考虑的本发明的特定实现,取得这些信息的方式可能不同。例如,在软件产品或者处理器实现的情况下,这些信息可以模拟或数字信号的形式取得,该信号或者直接从测量设备或者通过中间处理来接收。在医疗系统实现中,通过对患者进行医疗处理,可以取得流动特性信息。在这方面,诸如流动速度测量等信息可以侵入地取得,例如,通过将测量元素引入流动通道中,或者将探头置与患者体内靠近流动通道的地方,或者可以非侵入地取得,例如,通过接收来自通道的信号,诸如在超声模式情况下的反射信号。根据本发明,这些信息的获得可以与感兴趣的生理过程同步。
根据所考虑的应用,使用流动特性信息所获得的处理后信息可能不同。这些信息可能包括:(例如)有关流动通道的尺度信息,诸如半径、主/次轴尺度,横截面积或者表征通道尺度的其他参数;从尺度信息导出的信息,诸如定量流速;或者在其他方面依赖于这些尺度信息的信息(即使尺度信息没有被作为中间步骤确定)。在这方面,可以取得的医疗信息的例子包括:血管的面积、容积流速、压力梯度、感兴趣时间段内的血液容积或者弹性;以及患者的心脏抽运周期、容积输送量或者心室射血微量。
人们特别感兴趣的一种应用涉及确定诸如患者的升动脉等流体通道的容积流速。本发明人认识到可以分析运动流体的速度剖面中的变化来计算流体通道的尺度。进而,可以将流体通道尺度与流体速度剖面数据结合以计算在该通道中流动的流体的VFR。
在这方面,沿通道长度的非稳定层流包含以一定速度运动的流体元素,该速度取决于其与通道壁的距离、通道几何形状、作用在流体上的压力梯度、流体特性以及流体元素的初始速度。直接接触通道壁的流体元素不运动,且具有v=0,其中v为沿通道长度的速度。离开通道壁的流体元素的速度规则地跃迁到取决于与通道壁的距离的速度。当压力梯度不反向时,距离通道壁最远的流体元素的速度最大。流体元素相对于通道几何形状所呈现模式的形状与尺度定义了速度剖面。
本发明人还认识到通道的几何形状可以使用一个或更多个无量纲变量来表征,这些变量将诸如在圆形横截面半径上给定点等尺度值与所考虑尺度的最大范围相关。例如,在通道为圆柱型管的情况下,无量纲半径可以定义为在任意给定点上的半径除以该管的总半径。一个或更多个无量纲变量可以用来表征几何形状,例如,一个无量纲半径表征圆形管,两个无量纲轴表征椭圆形管(一个用于主轴,一个用于次轴)。更复杂的几何形状可以由多变量函数等等表示。无量纲时间可以表征速度剖面从一个形状变为另一形状所需的时间。如下所述,无量纲时间的定义涉及流体的粘度与密度以及时间与总体通道尺度。
在这方面,本发明另一方面针对一种外部测量流过通道的流体的速度的方法。该方法涉及使用诸如查询信号等非侵入手段,测量流过通道的流体的速度。进而,测定速度用来计算通道的面积,例如,横截面积,其然后与测定速度一起用来计算流过通道的流体至少一个VFR。
针对以上方面所述的特征存在各种改进。也可以将其他特征融入以上方面,以形成本发明的多个实施例。从以下描述中显然可以看到这些改进与其他特征,并且这些改进与特征可以独立或者联合地存在。例如,流体速度可以使用超声查询信号来测量。这样的速度测量可以由速度剖面来表征。进而,速度剖面函数可以用来计算第一时间的速度剖面的速度剖面参数。然后,速度剖面参数可以用来计算通过通道的流体的平均速度。可以使用速度剖面参数以及表征速度剖面如何随时间变化的函数关系来计算无量纲时间。进而,无量纲时间被与通道的尺度相关,从而可以计算通道尺度,并且用通道尺度来确定通道的横截面积。通道横截面积可以与平均流体速度一起使用来确定流过通道的流体至少一个VFR。在这方面,可以通过区分从运动流体发出的信号与从周围区域发出的信号,来处理速度测量中的误差,所述周围区域可能产生在其他情况下将扰乱对于速度剖面与流体通道尺度的确定的信号噪声。另外在这方面,可以吸收随机测量误差,从而可以从实际速度剖面中区分出这些误差。
根据本发明的另一方面,提供了一种计算VFR的系统。该系统包括数据处理器,该处理器使用测定流体速度来计算通道面积,并使用该通道面积与测定流体速度来计算流过通道的流体至少一个VFR。该系统可以进一步包括一个或更多个输出模块,以向用户提供指示至少一个容积流速的输出。该系统可以进一步速度测量设备,以测量流体速度并提供测定速度给数据处理器。在这方面,数据处理器可以包括用来使用两个或更多个流动速度剖面或者速度流动分布计算通道尺度的逻辑,以及用于使用通道尺度数据以及所测量的速度数据来计算至少一个VFR的逻辑。
针对以上方面所述的特征存在各种改进。也可以将其他特征融入以上方面,以形成本发明的多个实施例。从以下描述中显然可以看到这些改进与其他特征,并且这些改进与特征可以独立或者联合地存在。例如,这样的系统也可以包括用于从测定数据中区分流动速度分布的逻辑,所述测定数据还包括与有关于在通道内流动的组织的数据不同的数据。这样的系统也可以包括用于计算压力梯度与其他导出参数的逻辑,所述导出参数诸如VFR的变化速度或者将通道尺度与压力梯度相关的通道弹性的测度。数据处理器一般为一个或更多个电子设备,其使用一个或更多个半导体部件,诸如微处理器、微控制器或者存储器装置。数据处理器一般具有数据存储装置,其使用非易失性存储器设备,诸如一个或更多个磁介质(诸如EPROM或EEPROM)。数据存储装置可以用来(例如)存储用来促进计算的数据。可以被存储来促进计算的数据的例子包括但不限于被存储用来避免或减少对计算特殊函数的需求的数据,所述特殊函数诸如贝塞尔函数、Lommel函数、贝塞尔-零函数(Bessel-zero)、改进贝塞尔函数(modified Bessel)、伽马函数、对数-伽马函数(log-Gramma)以及超几何函数函数。
在本系统的一种实施方式中,速度测量设备可以使用查询信号,该查询信号被发送进入包含具有流动流体的通道的区域,所述流动流体的VFR或其他参数待定。这样的查询信号可以使用幅度随时间变化的发射能量,诸如超声波能量或者电磁内量(包括在可见频谱内的光以及频率在可见频谱上或下的电磁辐射)。超声波能量一般使用在50kHz与50MHz范围内的频率,或者更典型地在在500kHz与50MHz范围内。这样的查询信号可以被连续发送或者以脉冲发送。速度测量设备可以使用速度或相位的变化,在查询信号与反向散射能量或者产生反射的运动组织交互作用时发生所述变化。这种实施方式的一个例子包括但不限于多普乐测量。
用于计算VFR的系统可以进一步包括用于根据外部信号对测量进行定时的逻辑。这些外部定时信号包括与流动诱导现象(flow inducing phenomena)有关的信号,诸如与造成压力梯度变化的因素有关的信号。例如,该系统可以使用与心脏的电活动有关的信号,诸如心电图(EKG)信号,从而对何时进行测量进行排列或者控制,或者将测量相关。在这方面,控制何时进行测量可能影响如何解释已经进行的测量。
所述一个或更多个数据输出模块可以是任意适当类型,诸如显示器、可听声音生成装置或者数据输出端口,诸如那些发送电子或电磁信号的端口,以及相关的局域网或广域网端口。该系统也可以包括一个或更多个控制,诸如电源控制、信号敏感度控制或者使信号或所计算的结果根据流动介质的特性而改变的调整。这些调整可以包括那些基于流体特性的调整,并且包括对于标准流体特性的调整。此类标准流体特性的一些例子包括但不限于粘度、密度以及运动粘度。基于流体特性的调整可以包括与标准流体特性相关的特性,并且包括作为相关流体特性的血细胞比容。基于血细胞比容的调整可以包括使用血细胞比容的单一值、血细胞比容的范围或者与血细胞比容有关的患者特有数据,诸如患者种族、年龄或者性别。
根据本发明的另一方面,提供了一种用于计算流动通道的尺度相关信息的软件产品。该软件产品包括在处理器(例如数据处理器)上执行的数据处理器指令,以使用诸如与流体速度相关的信息等测定流动特性来计算尺度相关信息,诸如通道半径或者通道面积。该软件产品可以进一步包括用于以下目的的指令:使用尺度相关信息与测定流动特性一起来计算至少一个其他值,该值可能是进一步的尺度相关值,诸如流过通道的流体的VFR。该软件产品可以进一步包括输出指令,该指令被配置来提供指示流动特性、尺度相关信息以及所述其他值中至少一个的输出给用户。该软件产品可以进一步包括速度测量指令,用来取得流体速度的测量,并且提供测定速度给数据处理器。在这方面,数据处理器指令可以被配置来使用以下来计算通道尺度:两个或更多个流动速度剖面或速度流动分布,以及使用通道尺度数据与测定速度数据来计算至少一个VRF的指令。
针对以上方面所述的特征存在各种改进。也可以将其他特征融入以上方面,以形成本发明的多个实施例。从以下描述中显然可以看到这些改进与其他特征,并且这些改进与特征可以独立或者联合地存在。例如,这样的软件产品也可以包括用于从测定数据中区分流动速度分布的指令,所述测定数据还包括与有关于在通道内流动的组织的数据不同的数据。这样的软件产品也可以包括用于计算压力梯度与其他导出参数的指令,所述导出参数诸如VFR的变化速度或者将通道尺度与压力梯度相关的通道弹性的测度。
用于计算VFR的软件产品可以进一步包括用于根据外部信号对测量进行定时的指令。如上所述,这些外部定时信号包括与流动诱导现象有关的信号,诸如与造成压力梯度变化的因素有关的信号。例如,该软件产品可以使用与心脏的电活动有关的信号,诸如心电图(EKG)信号,从而对何时进行测量进行排列或者控制,或者将测量相关。
根据本发明的另一实施方式,可以基本实时地取得患者流动通道的基本实时的速度测量与基本实时的尺度相关信息。在每种情况下,可以非侵入地取得信息,诸如通过基于同一或多个信号集合的超声处理。例如,如上所述,可以使用超声测量,以取得连续时间上的速度剖面信息,并且与速度剖面时间变化有关的信息可以用来导出尺度相关信息。从速度测量和/或初始尺度相关信息,可以取得其他尺度相关信息,诸如VRF值或其他导出信息。在这方面,可以将连续速度剖面测量在时间上进行得足够紧密,使得诸如通道横截面积与压力梯度等其他流动参数可以假定不变,或者连续测量可以用脉搏周期同步和/或可以考虑到其他流动参数的变化。在任何情况下,可以取得瞬时结果,从而允许向医师提供更及时的展示,以及可能希望的对于生理过程的更好地相关。
根据本发明的另一方面,提供了一种处理信号中人为现象的方法与装置,所述信号被用来确定有关运动中的生理组织的信息,例如,在患者体内流动通道中的生理流体。本发明人已认识到这些人为现象可能与和周围组织相关联的信号部分有关,或者与和从流动通道发出的信号部分相关联的噪声有关。某些试图处理人为现象的现有技术已经认识到(例如)对于与流动通道有关的超声信息的谱分析显示双峰特性,但是这些现有技术假定可以使用频率过滤来孤立感兴趣的谱部分,而不引入不可接收的不准确性。根据本发明,提供了一种表征输入信号为感兴趣信号部分与非希望信号分量的数学模型。然后,使用该数学结构进行分析,以在数学上吸收非希望信号部分,而不用频率过滤以及相关假定。例如,可以使用多参数函数来将输入信号参数化,该输入信号包括非希望信号部分,并且可以使用公知的统计处理来吸收非希望信号分量。使用这种方法,可以吸收通道内与通道外人为现象之一或者两者。这就改进了结果的准确性并且放松了对于将查询信号对准流动通道的要求。在后一方面,可以精确地进行此处所公开的标准超声测量或更先进的测量,而不用将探头置于患者体内通道附近,或者没有与非常窄的信号相关联的对准困难。
附图说明
为了更全面地理解本发明及其进一步的优点,现在将结合附图详细描述本发明,其中:
图1显示已完全展开的层流的无量纲速度剖面(velocity profile);
图2显示正在展开的层流的无量纲速度剖面;
图3显示在从零速度开始后速度剖面的变化;
图4显示在从非零速度开始后的速度剖面的变化,该非零速度大于将在稳态下发生的速度;
图5显示当压力梯度反向时速度剖面的变化;
图6显示查询信号与产生测定信号的区域;
图7显示示例性无量纲速度剖面;
图8描绘圆形区域的累积概率分布;
图9显示无量纲速度(φ)如何随横截面积的微量变化;
图10描绘正在展开的层流的无量纲速度的示例性累积概率函数;
图11描绘已展开的层流的无量纲速度的示例性累积概率函数;
图12描绘无噪声与有噪声的示例性累积概率函数;
图13描绘无噪声与有噪声的示例性概率密度函数;
图14描绘包括缓慢运动区域与噪声的示例性累积概率函数;
图15描绘有与无误差源的示例性累积函数;
图16描绘有与无误差源的示例性密度函数;
图17显示系统的一种实施方式的方框图;
图18显示计算序列的一种实施方式;
图19显示谱累积数据点的例子;以及
图20显示正在进行拟合的谱累积数据点的例子。
具体实施方式
现在将参照附图,这些附图至少有助于描绘本发明各个有关特征。虽然将主要针对确定流过患者的通道(例如动脉)的流体(例如血液)的VFR来描述本发明,但应该理解本发明并不限于这种应用,本发明还可用于多种与在运动生理组织上进行测量有关的应用以及其他流动测量应用,尤其当不容易或者不适合或有困难直接接近流动通道。因此,符合以下描述以及现有技术的修改与变化落入本发明范围之内。此处所述的实施方式还用来解释实施本发明的最佳方式,并且使本领域技术人员能够在这些以及其他实施方式中使用本发明,或者在其他具有本发明的特定应用所需的各种修改的实施方式中使用本发明。
如上所述,沿通道长度的非稳定层流包含以依赖于以下因素的速度移动的流体元素:其距通道壁的距离、通道的几何形状、作用于该流体的压力梯度、流体的特性以及流体元素的初速度。直接接触通道壁的流体元素不移动,并且具有v=0,其中v为沿通道长度的速度。离开通道壁的流体元素的速度规则地过渡到依赖于距通道壁的距离的速度。当压力梯度不反向时,距离壁最远的流体元素的速度移动最快。图1显示在通道为圆柱形管并且施加压力梯度已经长得足以使流体流动已经达到稳定状态的情况下的速度。流体元素速度相对通道几何形状所具有的模式的形状与尺度为速度剖面。
在这方面,通道的几何形状可以使用无量纲变量来表征,该无量纲变量将诸如相对于圆形横截面半径的特定位置的值之类的尺度值与所考虑尺度的最大范围相关。例如,无量纲半径可以定义为在任意点上的半径除以该管的总半径。可以使用一个或更多个无量纲变量来表征几何形状:一个无量纲半径表征圆形管,两个无量纲轴表征椭圆形管(一个用于主轴,一个用于次轴)。在有些情况下,可以使用无量纲变量的替代值,诸如纵横比,其中纵横比为椭圆中主轴对次轴的比率,并且纵横比可以与无量纲主轴结合使用。
速度剖面从一个形状变换为另一形状所需的时间可以用无量纲时间表征。无量纲时间的一种定义使用流体的粘度与密度以及时间与总体通道尺度。在这方面,无量纲时间的一般形式为 τ = μt ρf ( D ) , ]]>其中μ=粘度,ρ=密度,f(D)为表征该通道的尺度的参数的函数。f(D)具有长度平方(L2)的尺度。
无量纲时间τ也称为时间常量,并且可以用来计算使特定事件发生所需多少个时间常量。例如,人体主动脉中血流量遵循由心速确定的周期模式。一般地,心跳周期持续大约0.01与0.1个时间常量(以下将更详细地描述所使用的时间常量),并且在大多数情况下持续在大约0.01到0.06个时间常量之间。本发明的一个方面是其计算对于小的以及大的时间常量的VFR,包括至少与在主动脉心跳周期中所遇到的时间常量一样小的时间常量。
在这方面,本发明使用在一个或更多个时间间隔上发生的、两个或更多个速度剖面之间的变化,以计算对于无量纲时间的值,然后使用无量纲时间的一个或更多个实例的值来计算f(D)中的参数。另外,在这方面,本发明可以使用表征f(D)的参数来计算流体通道的横截面积。另外,在这方面,本发明可以使用一个或更多个所计算的流体通道的横截面积以及一个或更多个速度剖面来计算一个或更多个VFR。
以下的推导与示例使用圆形管作为通道。然而,应该注意本发明不限于该几何形状。圆形管中的非稳定层流是数学推导的出发点。在这方面,沿管长度存在压力梯度,这使该流体加速,直到其达到稳态速度。与以下变量一起,可以使用圆柱坐标:
μ=流体粘度
ρ=流体密度
L=管长度
R=管半径
t=时间
p=压力。注意
为沿该管的压力梯度。
v,vz=沿管长度的速度。
Vmax在时间t=∞处流体的最大速度。在r=0处将产生Vmax,并且 V max = ( p 0 - p L ) R 2 4 μL ]]>为在t=∞处的最大速度。
等式1-3定义分别对于速度、半径以及时间的无量纲变量。在这方面,等式1定义无量纲速度,其为实际速度除以在时间t=∞处产生的最大速度,等式2定义无量纲半径,等式3定义无量纲时间。
φ = 4 v z μL ( p 0 - p L ) R 2 - - - ( 1 ) ]]>
ξ = r R - - - ( 2 ) ]]>
τ = μt ρ R 2 - - - ( 3 ) ]]>
可以对半径R解等式3:
R = μt τρ - - - ( 4 ) ]]>
从作为参考融入此处的参考文献1(R.Byron Bird,Warren E.Steward,andEdwin N.Lightfoot,Transport Phenomena,John Wiley&Sons,Inc.(1960),pages127-130),可以得知与流体通道中力平衡有关的以下关系:
ρ ( ∂ ∂ t v z ) = p 0 - p L L + μ ( ∂ ∂ r r ( ∂ ∂ r v z ) ) r - - - ( 5 ) ]]>
初始与边界条件如下:
对于从0速度开始的情况,初始条件:在t=0处,对于0≤r以及对于r≤R流体具有vz=0,其中R=圆形管的半径。
第一边界条件为:在r=0处,流体将具有vz=有限值。
第二边界条件为:在r=R处,流体具有vz=0。
然后可以将等式(5)乘以
并且在代入无量纲速度、半径和时间之后,变为:
∂ ∂ τ φ = 4 + ∂ ∂ ξ ξ ( ∂ ∂ ξ φ ) ξ - - - ( 6 ) ]]>
可以解等式(6),使得满足以下条件:
对于具有零速度的情况的初始条件为:在τ=0;φ=0
第一边界条件:在ξ=1;φ=0
第二边界条件:在ξ=0;φ=有限值。
当τ=∞时出现稳态。因此,解为稳态项与瞬态项的和:
φ(ξ,τ)=φ∞(ξ)-φt(ξ,τ) (7)
在τ=∞处出现稳态,并且在稳态处 ∂ ∂ τ = 0 . ]]>将此项代入等式(6),并且在φ=0与ξ=1处求解,产生等式8:
φ∞(ξ)=1-ξ2 (8)
参考文献1导出以下与速度剖面等式推导有关的表达式:
φ t ( ξ , τ ) = Σ n = 1 ∞ B n e ( - α ( 0 , n ) 2 τ ) J ( 0 , α ( 0 , n ) ξ ) - - - ( 9 ) ]]>
其中J(n,x)为第一种类型的n次贝塞耳函数,α(0,n)为第一种类型的0次贝塞耳函数的第n个正实根。
初始条件出现在τ=0时,所以
φ ( ξ , 0 ) = 1 - ξ 2 - ( Σ n = 1 ∞ J ( 0 , α ( 0 , n ) ξ ) B n ) - - - ( 10 ) ]]>
目前为止,推导还不依赖于初始速度分布。速度流动剖面等式的标准推导现在将假定初始速度为零。这样的假定大大简化了数学推导,但其没有生成可以用来确定信道尺度的等式。此处所示推导的剩余部分将偏离速度流动剖面等式的标准推导。
在这方面,初始速度分配可以由一般化等式表示。可以使用此类一般化等式的各种形式。当压力梯度不反向时(虽然在单一方向上的大小可以改变),以下形式的等式将正确地表示所有可能的初始速度分布:
φ=a(1-ξk) (11)
其中a=对于所测量的速度剖面的最大φ。
如果初始速度为零,则a=0。当流速从非零初始速度增加时,a<1。如果速度降低,则a>1。如果a>1,则压力梯度下降。同样的函数形式也可用于有量纲速度。在这种情况下,a=对于所测量的速度剖面的最大速度,并且k=速度剖面形状参数。k在闭区间[2,∞]上定义。对于已完全展开的层流,剖面为抛物线形,并且k=2。对于未完全展开的层流,k>2。
将等式(11)代入等式(10)中的φ(ξ,0)得到等式(12):
a ( 1 - ξ k ) = 1 - ξ 2 - ( Σ n = 1 ∞ J ( 0 , α ( 0 , n ) ξ ) B n ) - - - ( 12 ) ]]>
然后,通过将等式(12)乘以J(0,α(0,n)ξ)然后积分,可以对Bn求解等式(12)。结果是等式(13)。
∫ 0 1 J ( 0 , α ( 0 , m ) ξ ) ξa - J ( 0 , α ( 0 , m ) ξ ) ξa ξ k dξ = ∫ 0 1 J ( 0 , α ( 0 , m ) ξ ) ξadξ ]]>
+ ∫ 0 1 - J ( 0 , α ( 0 , m ) ξ ) ξa ξ 3 dξ + ( Σ n = 1 ∞ B n ∫ 0 1 - J ( 0 , α ( 0 , m ) ξ ) ξJ ( 0 , α ( 0 , n ) ξdξ ) - - - ( 13 ) ]]>
只有m=n时贝塞尔函数的正交性才导致贡献。因此,可以消除该和。对于Bm求解导致等式(14)
B m = - 2 aα ( 0 , m ) 2 - aα ( 0 , m ) ( 2 - k ) LommelS 1 ( 1 + k , 0 , α ( 0 , m ) ) - 4 J ( 1 , α ( 0 , m ) ) α ( 0 , m ) 3 - - - ( 14 ) ]]>
其中LommelS1为Lommel函数s。
对于初始条件具有非零速度的情况,对于φ(ξ,τ)的结果等式(包括表征初始速度条件的变量)为
φ ( ξ , τ , a , k ) = 1 - ξ 2 ( Σ n - 1 ∝ ( - 2 ( aα ( 0 , n ) 2 - aα ( 0 , n ) ( 2 - k ) LommelS 1 ( 1 + k , o , α ( 0 , n ) ) - 4 ) c ( - α ( 0 , n ) 2 τ ) J ( 0 , α ( 0 , n ) ξ ) J ( 1 , α ( 0 , n ) ) α ( 0 , n ) 3 ) ) - - - ( 15 ) ]]>
等式(11)描述了速度分布。在测量速度剖面的任意特定时间,其将具有值at,kt,从而导致等式(16)
φ ( t ) = a t ( 1 - ξ k t ) - - - ( 16 ) ]]>
无量纲时间τ为时间的函数。速度分布从初始速度剖面开始展开,该初始速度剖面由等式(11)表征,其中a与k为表征初始速度分布的形状的参数。
结果为等式(17):
a t ( 1 - ξ k t ) = 1 - ξ 2 ( Σ k = 1 ∞ ( - 2 ( aα ( 0 , n ) 2 - aα ( o , n ) ( 2 - k ) LommelS 1 ( 1 + k , o , α ( 0 , n ) ) - 4 ) c ( - α ( 0 , n ) 2 τ ) J ( 0 , α ( 0 , n ) ξ ) J ( 1 , α ( 0 , n ) ) α ( 0 , n ) 3 ) ) - - - ( 17 ) ]]>
从n=1到∞的和可以通过向N求和来近似,如等式18所示:
a t ( 1 - ξ k t ) = 1 - ξ 2 ( Σ k = 1 ∞ ( - 2 ( aα ( 0 , n ) 2 - aα ( o , n ) ( 2 - k ) LommelS 1 ( 1 + k , o , α ( 0 , n ) ) - 4 ) c ( - α ( 0 , n ) 2 τ ) J ( 0 , α ( 0 , n ) ξ ) J ( 1 , α ( 0 , n ) ) α ( 0 , n ) 3 ) ) - - - ( 18 ) ]]>
流动剖面取决于四个变量。其中两个变量a与k定义初始速度剖面。第三个变量τ表征所经过的时间。第四个变量ξ表征管的半径。
在这方面,本发明使用来自速度测量的数据,以计算速度剖面参数,然后根据在已知时间段内速度剖面参数中所发生的变化,确定诸如信道尺度等其他变量的值。以下方法为使用速度剖面变化以确定通道尺度,然后计算VFR以及导出参数的例子。使用以上公式计算通道尺度与VFR所需的一般步骤如下:
(1)在特定时间,使用适当的装置进行速度测量,诸如通过使用诸如超声波束或电磁波等查询信号,并生成速度剖面。众所周知,已知根据信号转移时间以及在该介质中的信号速度,可以确定与所接收的信号相关的范围信息。因此,可以确定对于流体粒子的速度测度以及这些粒子的位置,以取得剖面有关信息。
(2)使用该速度剖面来计算表征速度剖面函数的参数的值。这些参数的例子有等式(11)中的a与k,并且等式(11)为表征速度剖面的函数的例子。如下详述,这些计算可以只使用两个数据点,或者结合诸如最小二乘或最小绝对值方法之类的统计曲线拟合技术,该计算可以使用多于两个数据点。
(3)在经过了已知时间间隔t(例如,以秒测量)后,进行更多的速度测量,并且计算在时间t上表征速度剖面的参数的值。例如,使用等式(11),对于第二个剖面a与k可以被标记为at、kt。具体地,可以分析在短得足够使诸如通道尺度与压力梯度等流动参数的变化可以被安全地假定为可忽略不计的时间段上的速度剖面变化。这些参数随着患者的脉搏周期变化,该脉搏周期一般具有的0.5~4.0Hz的频率。因此,用来分析速度剖面参数变化的时间间隔最好不长于大约0.025秒,并且最好不长于大约0.01秒。然而,应该理解,可替换地,此(多个)时间间隔可以被相对于患者脉搏周期地相位同步,或者可以考虑流动参数的变化。
(4)使用表征速度剖面如何随时间变化的函数关系计算无量纲时间τ。例如,使用等式(17)或(18)以及a、k、at、kt的值来计算τ。
(5)使用将无量纲时间与通道尺度相关的等式来计算通道尺度。例如,使用等式(4)连同t的已知值、粘度以及密度(或者运动粘度,其为粘度与密度的商)来计算管的半径R。
(6)使用现在已知的通道尺度来计算通道的横截面积A。例如,使用半径R来计算横截面积。
(7)使用测定的速度剖面来计算平均速度vm。例如,如下所述,由等式11定义的、通过具有圆形横截面的管的速度剖面的平均值为 v mean = ka 2 + k , ]]>其中a为在测量时间发生的最大有量纲速度。如下所述,在这种情况下速度剖面函数的参数为a与k,可以直接求解,从而在计算平均速度时可以方便地使用它们。
(8)计算VFR=A×vm。
(9)使用定义无尺度速度的关系,计算压力梯度。例如,使用等式(1)重新排列的形式计算压力梯度。
(10)计算在两个不同时间的半径与压力梯度,并且将半径的变化除以压力梯度的变化以计算通道弹性。
如上所述,当初始速度剖面为零时,则在等式(15)中a=0。在稳态下(τ=∞),完全展开的层流剖面为抛物线形,如图1所示。最高速度在管的中心处,并且在管壁处v=0。
随着流动的展开,速度剖面变钝,如图2所示。图2表示当初始速度为零并且已经过足够时间τ=0.1时出现的情况。图3绘出当其从τ=0变到τ=1的速度剖面。当τ=1时速度剖面几乎完全展开。图4绘出当由于压力梯度下降速度减慢时的情况。
图1-4所示情况的所有流动剖面都由等式(11)准确表示。该函数形式为N趋向∞时等式(18)的极限。
图5显示压力梯度已经反向的情况。在图5中,无量纲半径在标记为r/R的轴上。图5还在图的底部显示了非负开始速度剖面,并在图的顶部显示了最终的稳态速度剖面。这两个剖面都符合等式(11)的函数形式。图5中的中间曲线速度剖面显示了当流场在初始速度剖面与最终剖面之间反向时的转换速度剖面。在转换时,管中心处的流体向下流,而边沿处的流体已经反向并向上流。例如,如果主动脉瓣膜具有显著回流,则可能发生这种情况。对于压力梯度反向的情况,需要对于初始速度剖面的更复杂的表示,并且比等式(11)所使用的两个参数包含得更多。
如上所述,本发明的另一方面是对误差来源的调整。误差可能存在于所测定的初始速度剖面中,尤其是当使用查询信号来产生测定数据时。图6示意性地显示了诸如超声波等查询信号被发射到包含流动流体的区域中的情况,该流动流体包含在管中,该管被比该管中流动的大部分组织的运动(与该流动流体同向或者反向)慢得多的组织区域所包围。与该流动流体比较,所述缓慢运动的组织基本静止。为了说明的目的,图6中只显示了包含流体的管的一部分。
图6所显示的情况类似于以下条件下所发生的情况:如果对在升动脉中流动的血液进行超声多普乐测量,并且多普乐换能器置于胸骨上切迹并且向心脏对齐。某些查询信号被反射或者反向散射,这没有在图6中示出,并且这些信号可能被作为返回信号而测量。可以分析返回信号的频谱以生成对于所查询区域的速度剖面的估计。
在这种情况下,存在两个误差来源。在从流动区域所测量的信号中,可能存在第一误差来源,其增加诸如噪声等误差。在从表示缓慢运动的组织区域的信号中,可能存在第二误差来源,其混淆来自流动流体的信号,从而缓慢运动的速度信号被添加到来自的流动流体的信号中。根据本发明,这两个信号源得到处理,并且被防止在所估计的VFR中产生显著误差。一般的解决方案是:通过引入选择性地从测定信号中挑选这些误差的数学结构来吸收误差。以下是对可用来估计测定信号分布或者从测定信号中导出的变量的误差吸收函数的一般与具体例子的推导。测定信号的例子是Doppler测定的频率。导出变量为使用一个或更多个测定信号所计算的变量,诸如从多普乐频率测量计算出来的速度。
可用来估计通道尺度与VFR的一般误差吸收函数的一般形式由以下等式(19)给出:
M(f)=S(f,|xs|)+F(f,|xF|) (19)
其中
M(f)=测定信号
S(f,|xx|)=来自缓慢运动区域的信号,其中|xx|为表征该缓慢运动区域的(待定)参数向量。
F(f,|xF|)=来自运动区域的信号,其中|xF|为表征诸如噪声等误差的(待定)参数向量。
图7显示使用a=1等式(11)而计算的、在通道中运动的流体的速度剖面。无量纲半径的任意特定值ξ*,映射到有且只有一个无量纲速度φ*。对于管来说,每个无量纲半径形成围绕该管中心的圆。所有在特定ξ*处移动通过该管的流体将以同一速度φ*运动。任意随机选择的流体元素具有100%的可能性位于在管中心与ξ=1的最外半径之间的区域内。流体元素位于从最外半径延伸到内侧半径的带内的可能性为管的总横截面积与区间[ξ,1]上ξ的带的半径的比率。因此, F ( ξ ) = π - π ξ 2 π , ]]>其化简为等式(20),
F(ξ)=1-ξ2 (20)
该结果可以被认为是随机流体元素处于在管的半径ξ与外侧半径(ξ=1)之间的带中的概率。例如,随机ξ位于ξ=0.9与ξ=1半径之间的环中的概率为0.19。换而言之,圆圈面积的19%位于半径ξ=0.9之外。
对于在其之外存在所选择的面积部分的ξ值来解决相关问题。在这种情况下,等式(21)可以用来显示在其之外存在19%面积的半径为ξ=0.9。
ξ = - F + 1 - - - ( 21 ) ]]>
计算在其之外存在90%面积的半径得到ξ=0.3162。等式(21)可以被认为是累积概率分布,如图8所示。图8显示了(例如)100%的面积位于0以上的半径,而没有面积位于半径1之外。另外,图8显示面积的19%位于半径ξ=0.9之外,90%的面积位于由ξ=0.3162与ξ=1所限定的区域之内。
无量纲速度与无量纲半径之间的关系,即等式(11),φ=α(1-ξk),显示当ξ变小时φ变大。因此,将等式(21)代入等式(11)得到等式(22)或者当环孔在ξ与1之间延伸时所发生的最大无量纲速度。
φ=-a(-1+(-F+1)(1/2k)) (22)
例如,当a=1,k=6,F=0.5时所发生的最大无量纲速度为φ=0.0875。可以绘制等式22,如图9所示,其显示φ如何随微量(fraction)F的变化而改变。注意图9使用a=0.1与k=6进行计算。
φ的最大范围为区间[0,a]。变量a为在测定速度剖面时所发生的无量纲速度的最大值。
等式(22)中的微量F可以被认为是概率,并且可以针对更常规的等式(23)求解,该等式表示给出小于或等于特定φ的无量纲速度的微量的φ的累积分布。
F ( φ ) = - ( - φ - a a ) ( 2 1 k ) + 1 - - - ( 23 ) ]]>
φ具有闭区间[0,a]上的域,其中a为在特定测量时所发生的最大无量纲速度。因此,该累积函数只是对于[0,a]上的φ进行了定义。例如,如果a=0.1并且k=6,则对φ可能具有的最大范围的绘图导致图10。图10显示当a=0.1并且k=6时,50%的无量纲速度小于或等于φ=0.0875。对于已完全展开的层流a=1且k=2,相应的累积分布具有在0与1之间的定义的φ,如图11所示。因此,相应的等式为F已完全展开的层流(φ)=φ,并且对于已完全展开的层流,具有任意特定无量纲速度的概率是线形的。
通过使用等式(24)可以找到参数化的速度剖面的平均值。
φ meon = ka 2 + k - - - ( 24 ) ]]>
如果无量纲速度被用来对a与k求解,则k的值将与无量纲速度用于计算时相同,并且a将成为在该测定速度剖面中的最大速度的有量纲值。然后,可以使用等式(25)计算平均有量纲速度。
v mean = ka 2 + k - - - ( 25 ) ]]>
等式(25)为累积概率函数。通过对累积概率函数进行微分,找到φ的概率密度函数。结果在等式(26)中给出。与上相同,φ存在于区间[0,a]上。当存在已完全展开的层流(k=2且a=1)时,结果将为f(φ)=1。对于正在展开的层流(包括当压力梯度降低时所发生的)k>2,因此正在展开的层流的概率密度函数总具有小于1的正指数。
f ( φ ) = 2 ( - φ a + 1 ) ( 2 1 k - 1 ) ka - - - ( 26 ) ]]>
测量误差改变所测定的φ为φ的正确值的概率。发生特定φn以及其上被添加了噪声而导致小于或等于φ的总值的概率为具有φn的概率乘以φn+噪声≤φ的概率的乘积,其可以表示为p(φ)p(φn+噪声≤φ)。这两个概率的乘积为总测定值φn≤φ的累积概率。
具有特定φn的概率为在φndφ处计算的密度函数的乘积。设p(φn)表示φ具有特定值φn的概率。一般地,噪声具有与其相关的、由参数向量H定义的分布。对于高斯分布,参数为平均与标准偏差。设对于噪声生成函数的概率密度函数为gnoise(φ,H)。具有特定φn、且围绕其周的噪声小于或等于速度φ的概率为两个概率的乘积:(特定φn的概率)×(x+φn≤φ的概率)。结果是等式(27)给出的联合概率:
p ( φ n , x + φ n ≤ φ ) = p ( φ n ) ∫ - ∞ φ g noise ( x , H ) dx - - - ( 27 ) ]]>
该累积概率为对于所有φn具有(x+φn≤φ)的概率,为对于所有φn的各个概率的和,其与在所有φn上积分相同。将等式(26)作为p(φn)并且设置积分使φn在区间[0,a]上(因为对于φ最大值为a),这将导致等式(28),其形成一般情况:
F noise ( φ ) = ∫ 0 a 2 ( - φ n a + 1 ) ( 2 1 k - 1 ) ∫ - ∞ φ g noise ( x , H ) dx ka d φ n - - - ( 28 ) ]]>
人们特别感兴趣的情况是当噪声为高斯(正态)分布时。当假定许多独立的因素对测量误差有贡献并且这些因素产生相加效应是合理的时,中心极限定理指出误差关于真值具有高斯分布。因此,噪声密度函数以φn为中心,并且φn为平均值,而且噪声函数具有标准偏差σ:
g noise ( φ n ) = 1 2 2 e ( - 1 / 2 ( x - φ n ) 2 σ 2 ) σ π - - - ( 29 ) ]]>
将此等式与等式(27)结合,产生联合概率:
p ( φ n , x + φ n ≤ φ ) = p ( φ n ) ∫ - ∞ φ 1 2 2 e ( - 1 / 2 ( x - φ n ) 2 σ 2 ) σ π dx - - - ( 30 ) ]]>
该累积概率为对于所有φn具有(x+φn≤φ)的概率,为对于所有φn的各个概率的和,其与在所有φn上的积分相同。代入并在区间[0,a](因为对于φ最大值为a)上关于φn积分将导致等式(31):
F noise ( φ ) = ∫ 0 a 2 ( - φ n a + 1 ) ( 2 1 k - 1 ) ∫ - ∞ φ 1 2 2 e ( - 1 / 2 ( x - φ n ) 2 σ 2 ) σ π dx ka d φ n - - - ( 31 ) ]]>
在计算并用x替换φn后,等式(31)成为显式地包含对于实际速度分布的φ的分布以及改变该实际值并产生测定分布的随机噪声的等式,为以下等式:
F noise ( φ ) = ∫ 0 a 2 ( - x a + 1 ) ( 2 1 k - 1 ) ( - 1 2 erf ( 1 2 2 ( - φ + x ) σ ) + 1 2 ) ka dx - - - ( 32 ) ]]>
在这种情况下,无噪声的φ在区间[0,a]上,而测定φ(包含噪声)可以延伸到区间[0,a]之外。图12显示在无量纲速度φ中的高斯分布随机误差如何影响累积分布。这些误差扩展了测定速度的可能值,从而它们既小于也大于在流动流体中所发生的实际值。
等式(28)为具有用来补差测定信号中误差的误差吸收函数的累积谱函数的一般形式。等式(32)为当测定误差函数为高斯型时的特定情况。不用高斯型,也可以使用其他函数形式,诸如指数型、柯西型以及Riceian。
因此相应于等式32的概率密度函数为:
f noise ( φ ) = ∫ 0 a ( - x a + 1 ) ( 2 1 k + 1 ) e ( - 1 / 2 ( - φ + x ) 2 σ 2 ) 2 ka π σ dx - - - ( 33 ) ]]>
图13显示该概率密度函数如何受高斯分布随机误差的存在的影响。应该注意如图13所示的概率密度函数的扩展是测量误差引起的。使用误差吸收函数来降低或去除这些误差,包括那些由随机波动所引起的、一般被称为“噪声”的误差,是被用做本发明的新颖性的一部分的一般技术。然而,应该注意,根据本发明,使用误差吸收函数来降低或去除这些误差特别有利于确定速度剖面(不管是无量纲还是有量纲)以及随后的、导致确定VFR的计算。
概率密度函数为由多普乐测量方法所造成的频移而产生的谱密度函数。当引入测量噪声时,可以求解表征误差吸收函数gnoise(φ,H)的参数,以及表征速度剖面的形状的a与k的值。例如,如果噪声吸收函数为高斯型,则求解标准偏差的值以表征该噪声吸收函数。
如上所述,图9显示查询信号被用作确定流动区域的VFR的过程的一部分的确定速度剖面的系统的一部分。在这方面,以上误差吸收函数可以用来减少或去除在流动区域的信号中的测量误差。然而,图9还显示了另一可能误差来源为来自在该流动区域之外的、缓慢流动或非流动的组织的混淆信号。在这方面,可以如下地降低或去除由缓慢流动或非流动的组织的信号所引起的误差。
累积概率分布的一般形式为
F(φ)=Fslow(φ)+Fflowing(φ) (34)
如下详述,Fslow(φ)表示例如由等式32所表示的、在其已经对于产生具有表示非流动组织(或者不运动或者缓慢运动)的流动Fslow(φ)的信号的部分区域进行了调整之后的流动组织的分布。Fslow(φ)将由非流动速度的分布表示,其一般为高斯型,通常具有非零平均值。然而,应该注意,对于特定应用,最适合的分布可以包括指数型、柯西型、Riceian型或者其他分布。然而,在缺乏与此相反的其他信息时,高斯分布要所使用的适当分布,这是因为假定非流动区域比流动区域的运动慢很多。另外,可以通过多个独立的时变参数来产生非流动组织的速度,而所有这些时变参数都使用中心极限定理。
作为示例,对于以下推导,将使用高斯分布作为误差吸收函数。可以假定缓慢运动的速度具有高斯分布,并且平均值可以具有任意值,包括在进行测量的任意时刻的正值、负值或者零。另外,例如,如果测量的是包括升动脉的区域并且在测量期间所发生的呼吸过程引起离开信号源位置的非流动区域的缓慢运动,则可以发生该非零平均值。在这种情况下,部分速度信号将来自流动区域,其余信号将来自缓慢运动区域。来自这两个区域的信号微量用来对每个区域对总体速度分布的贡献进行加权。例如,在查询信号为多普乐系统的一部分的情况下,缓慢运动的速度产生多普乐反射,其相加于由流动的速度所产生的反射。它们产生了特定速度上的反射总和。在这方面,所有反射被相加,并且总和被用来将每个速度上的计数数目规格化。在这方面,来自特定区域的计数数目取决于其面积,例如,来自缓慢运动区域的计数数目取决于其面积,并且来自流动区域的计数数目取决于其面积。对于每个区域的总面积的微量——其在概率分布中最重要——定义了对于来自流动区域与来自缓慢运动区域的信号给定多少权重。在这方面,如果As=缓慢运动区域的面积,Af=流动区域的面积,缓慢运动区域的微量为 f s = A S A S + Af , ]]>流动区域的微量为 f f = A f A S + A f . ]]>总微量必须等于1,因此fs+ff=1必须为真,这就有fs=1-ff。这种关系导致等式(35),其为对于φ的总累积概率
F noiseslow ( φ ) = ( 1 - f f ) ∫ - ∞ φ 1 2 2 e ( - 1 / 2 ( φ - μ s ) 2 σ s 2 ) σ s π dφ + f f ∫ 0 a 2 ( 1 - x a ) ( 2 1 k - 1 ) ( 1 2 erf ( 1 2 2 ( φ - x ) σ ) + 1 2 ) ka dx - - - ( 35 ) ]]>
具有任意特定φ或更小的总概率为两个概率的和。左边的项为与缓慢运动区域相关的概率,第二项为具有由高斯噪声表示的测量误差的流动速度的概率。
在缓慢运动区域由高斯分布表示的情况下,由缓慢运动区域引入的误差由两个参数表征。一个参数是平均值μs,另一个为标准偏差σs。当进行速度测量时,目标仍为计算a与k。使用误差吸收函数需要计算表征误差的参数,从而去除误差的效果。例如,在等式(35)中,计算以下参数:a、k、ff、σ、μs以及σs。一般地,如此情况所示,可以使用对速度的多次测量,并且使用诸如最小二乘或者最小绝对值之类的统计方法来取得测定数据的最佳拟合。一般地,这些等式是非线性的,并且需要实施非线性拟合技术,诸如非线性最小二乘法。
图14显示对于缓慢运动的平均速度为零的情况基于等式(35)的累积概率函数。
与等式(35)相关的概率密度函数在等式(34)中给出:
F noise , slow ( φ ) = 1 2 ( 1 - f f ) 2 e ( - 1 / 2 ( φ - μ s ) 2 σ s 2 ) σ s π + f f ∫ 0 a ( 1 - x a ) ( 2 1 k - 1 ) e ( - 1 / 2 ( φ - x ) 2 σ 2 ) 2 π σka dx - - - ( 36 ) ]]>
图15与16为对于ff=0.7、μs=0.05、σs=0.04、a=0.4、k=8、σ=0.02的情况的累积函数与密度函数。这些图显示当没有测量噪声或者缓慢运动区域效应以及引入这两个误差源时的无量纲速度分布。这些图还显示了当已知误差吸收函数的参数时,可以恢复成为无误差数据的计算结果。
以上对于密度函数以及累积概率函数的等式也可以用于有量纲速度,而不是无量纲速度。测定的速度数据将具有速度单位,并且在变量φ在等式中出现处使用。在以下的讨论中,当既可以使用有量纲速度也可以使用无量纲速度时,将使用一般名词“速度”。
以上对于密度函数的等式已经被描述为概率密度函数。这种描述便于说明推导。这些密度函数一般被称为谱密度函数,并且定义对于测定流场的速度谱。当进行超声波测量(包括多普乐测量)时,使用标准计算方法,频谱可以转换为速度谱。
在本发明的一种实施方式中,测定数据可能拟合一个或更多个多参数函数,这些函数包含对于速度谱密度以及对于一个或更多个误差吸收函数的参数。在本发明的另一种实施方式中,测定数据可能拟合一个或更多个多参数函数,这些函数包含对于累积速度谱以及对于一个或更多个误差吸收函数的参数。使用累积速度谱而不是密度函数降低了测量误差的效应,这是因为累积速度谱是密度函数的平滑积分。
当对速度谱密度函数或者类似速度谱函数进行拟合时,拟合过程将确定导致在测定数据与函数预测值之间最小误差的参数的值。这些最小误差就是使用任何适当标准所确定的,并且包括最小化最小二乘误差的和或者最小化最小绝对值误差的和。拟合过程可以将所有参数作为一个集合拟合,可以首先估计对于误差吸收函数的值,并且在计算对于定义速度谱的参数的值之前去除误差效应。首先去除误差效应的分步处理导致计算几乎无误差的数据分布。所述无误差数据为对于在通道中流动的组织的速度分布,人们希望得到对于它的速度谱参数。例如这些参数为使用等式(11)计算的a与k,以用于等式(17)或等式(18)来计算τ,从τ计算半径R。
在参数拟合过程中,一般使用迭代处理,其中系统地调整参数值以取得对于速度数据的最佳拟合。在本发明的一种实施方式中,该迭代处理系统地选择误差吸收参数的值,以计算近似无误差的速度分布。然后,使用该近似无误差的速度分布来计算速度分布参数。在函数对于测定速度数据的拟合(Over fit)不够好的情况下,可以使用速度分布参数来计算误差吸收参数的修正值。然后,使用该修正后的误差吸收参数来重新计算速度分布参数的值。继续这种去除误差并计算速度分布参数的迭代处理,直至产生对于数据的适当好的拟合。使用这种在计算速度分布参数之前首先从测定数据中去除误差的实施方式,具有以下优点:在指定数目的迭代中一般会产生对于速度分布参数的更好的估计。应该注意取得对于速度分布参数(诸如等式(18)中的a与k)的准确估计比计算其他参数值更重要,这是因为计算无量纲时间τ所使用的是速度分布参数,无量纲时间τ随后成为计算流动通道尺度的基础。
所描述的用于使用误差吸收函数来降低或者去除误差的方法适用于无量纲速度与有量纲速度两者。当使用有量纲速度v而不是无量纲速度φ时,与速度相关的参数具有量纲。例如,当使用高斯分布时,平均值与标准偏差将具有速度单位。类似地,速度剖面参数可以具有单位,因为(例如)上面所使用的变量a可能具有速度单位。
在已经测量了两个速度分布并且已知其表征参数之后,可知变量a、k、at以及kt的值,并且可以用于等式(17)或等式(18)以求解τ。
然而,应该注意,通过使用不同的等式,求解τ可以在数学上更容易一些。在这方面,测量速度剖面时的峰值速度为at,at为管中心处的速度,此时半径为零并且ξ=0。将此代入等式(18)产生对于at的等式:
a t = 1 - 2 ( Σ n = 1 N ( 4 - aα ( 0 , n ) 2 - aα ( 0 , n ) ( 2 - k ) LommelS 1 ( k + 1 , 0 , α ( 0 , n ) ) ) e ( α ( 0 , n ) 2 τ ) J ( 1 , α ( 0 , n ) ) α ( 0 , n ) 3 ) - - - ( 37 ) ]]>
在这种情况下,等式(38)为从等式(16)开始并求解kt的结果,假定对于可以表示为ξi与φi的、ξ与φ的选定值的已知值φi。
k t = ln ( - φ t + a t a t ) ln ( ξ t ) - - - ( 38 ) ]]>
除了导致不确定操作的极端值ξi=0或者ξi=1,可以为ξt选择任意值,以及对于φi的相应值。因此,ξi与φi可以是速度剖面曲线上的任何有序对。因为我们从等式(16)与(18)知道φ为τ与ξ的函数,所以可以通过代入等式(38)来消除φi。结果为等式(39)
k t = ln ( - ξ t 2 + 2 ( Σ n = 1 N ( - 4 + aα ( 0 , n ) 2 - aα ( 0 , n ) ( 2 - k ) LommelS 1 ( k + 1,0 , α ( 0 , n ) ) ) e ( - α ( 0 , n ) τ ) ( J ( 0 , α ( 0 , n ) ξ t ) - 1 ) J ( 1 , α ( 0 , n ) α ( 0 , n ) 3 ) - 1 + 2 ( Σ n = 1 N ( - ( - 4 + aα ( 0 , n ) 2 - aα ( 0 , n ) ( 2 - k ) LommelS 1 ( k + 1,0 , α ( 0 , n ) ) ) e ( - α ( 0 , n ) 2 τ ) J ( 1 , α ( 0 , n ) ) α ( 0 , n ) 3 ) ) ) ln ( ξ t ) - - - ( 39 ) ]]>
等式39用表征初始速度分布(a,k)、无量纲时间(τ)(其也是时间常量)以及选定ξ(其标记为ξi)的参数给出了kt。
不再使用φ。该等式隐含地用kt给出了τ。因此,当初始速度剖面以及随后的速度剖面都已经被表征后,就可以确定时间常量τ的值。对于这些计算,可以选择开区间(0,1)上的任意ξi。例如,可以做出使用ξi=0.9的判定。在向ξi分配了显式值之后,结果为对于τ的任意指定值给出kt的显式等式。例如,如果ξi=9/10,则结果为等式(40)。在计算之后,等式(40)变为等式(41),当已知a、k以及kt的值时,等式(41)可以对τ迭代求解。如上所述,LommelS1为Lommels函数。
k t = ln ( - 81 100 + 2 ( Σ n = 1 N ( - 4 + aα ( 0 , n ) 2 - aα ( 0 , n ) ( 2 - k ) LommelS 1 ( k + 1,0 , α ( 0 , n ) ) ) e ( - α ( 0 , n ) 2 τ ) ( J ( 0 , 9 10 α ( 0 , n ) ) - 1 ) J ( 1 , α ( 0 , n ) σ ( 0 , n ) 3 ) - 1 + 2 ( Σ n = 1 N ( - ( - 4 + aα ( 0 , n ) 2 - aα ( 0 , n ) ( 2 - k ) LommelS 1 ( k + 1,0 , α ( 0 , n ) ) ) e ( - α ( 0 , n ) 2 τ ) J ( 1 , α ( 0 , n ) ) α ( 0 , n ) 3 ) ) ) ln ( 9 10 ) - - - ( 40 ) ]]>
k t = - 9.491221577 × ]]>
ln ( - 81 + 2 ( Σ n = 1 N ( - 4 + aα ( 0 , n ) 2 - aα ( 0 , n ) ( 2 - k ) LommelS 1 ( k + 1,0 , α ( 0 , n ) ) ) e ( - α ( 0 , n ) 2 τ ) ( J ( 0 , 9 10 α ( 0 , n ) - 1 ) ) J ( 1 , α ( 0 , n ) α ( 0 , n ) 3 ) - 1 + 2 ( Σ n = 1 N ( - ( - 4 + aα ( 0 , n ) 2 - aα ( 0 , n ) ( 2 - k ) LommelS 1 ( k + 1,0 , α ( 0 , n ) ) ) e ( - α ( 0 , n ) 2 τ ) J ( 1 , α ( 0 , n ) ) α ( 0 , n ) 3 ) ) ) ]]>
(41)
计算Lommels函数很耗时,并且不容易准确计算。因此,本发明的一个方面是预先计算用速度剖面表征参数来隐型定义τ的等式中的一个或更多个项的值,然后在用来形成可用于计算时间常量τ的值的内插函数的简化等式中使用这些预先算好的值。以下为如何实现本发明这个方面的例子。在这个例子中,等式(41)可以用来预先计算其求和项的值,并且这些预先算好的值可以用于简化等式,该简化等式形成将ki与τ相关的内插函数。
当τ与k的数值值代入等式(41)时,结果为a表示的kt的等式。例如,如果τ=0.1且k=12,则等式(41)变为等式(42):
k t = - 9.491221577 ln ( . 01000000000 - 27.73984105 - 70.73715877 a - . 3851895036 - . 8293031252 a ) - - - ( 42 ) ]]>
等式(42)具有等式(43)的形式。等式(39)以及诸如等式(41)等基于等式(39)的等式可以转换为等式(43)的形式
k t = A ln ( B + Ca D + Ea ) - - - ( 43 ) ]]>
其中
A = 1 ln ( ξ i ) - - - ( 44 ) ]]>
B = - ξ i 2 - 8 ( Σ n = 1 N e ( - α ( 0 , n ) 2 τ ) ( J ( 0 , α ( 0 , n ) ξ i ) - 1 ) J ( 1 , α ( 0 , n ) ) α ( 0 , n ) 3 ) - - - ( 45 ) ]]>
C = - 2 ( Σ n = 1 N ( - α ( 0 , n ) 2 + α ( 0 , n ) ( 2 - k ) LommelS 1 ( k + 1,0 , α ( 0 , n ) ) ) e ( - α ( 0 , n ) 2 τ ) ( J ( 0 , α ( 0 , n ) ξ i ) - 1 ) J ( 1 , α ( 0 , n ) ) α ( 0 , n ) 3 ) - - - ( 46 ) ]]>
D = - 1 + 8 ( Σ n = 1 N e ( - α ( 0 , n ) 2 τ ) J ( 1 , α ( 0 , n ) ) α ( 0 , n ) 3 ) - - - ( 47 ) ]]>
E = 2 ( Σ n = 1 N ( - α ( 0 , n ) 2 + α ( 0 , n ) ( 2 - k ) LommelS 1 ( k + 1,0 , α ( 0 , n ) ) ) e ( - α ( 0 , n ) 2 τ ) J ( 1 , α ( 0 , n ) ) α ( 0 , n ) 3 ) - - - ( 48 ) ]]>
一旦选定ξi,A总是相同。对于任意选定的τ,B与D为常量。虽然可以在任何时间计算值,但优选方法是在需要测量之前选择τ与k的值,并且在表中预先计算并存储B、C、D以及E的值。例如,这些值可以存储在便于快速访问的查找表中,诸如按τ与k的值索引各个项。使用这种方法,求解τ开始于选择τ的假定值,并且使用包含k的实际值与τ的假定值的查找表项,并且使用a的实际值计算对于kt的计算值。选择τ的一组假定值,并计算kt的其他值。结果是将kt与τ相关的表。在这方面,可以使用适当的内插函数,诸如内插多项式,以将τ与kt相关。然后,可以使用该内插函数来使用kt的实际值求解τ。
上述对于用于有效求解τ的方法的推导与描述使用了ξ与无量纲速度。同样的用于有效求解τ的方法的推导与描述也使用ξ与有量纲速度来应用。使用有量纲速度的结果除了速度剖面参数a包含速度量纲外完全相同。
本发明的另一方面是从一个或更多个VRF以及与用来计算VFR的速度剖面相关的时间范围中发现周期性与非周期性流速行为的存在与性质。用来发现周期性流速行为的存在与性质的手段包括诸如傅立叶分析、人工智能、统计分析以及包含在复杂度理论与分析中所使用的非线性分析之类的数学手段。识别周期性流速行为的存在与性质包括识别揭示流动周期性行为的流速变化,包括来自活体患者的心跳速度的抽运周期(pumping cycle period)。
本发明的另一方面是使用一个或更多个VFR以及在速度测量之间的一个或更多个时间间隔,以计算在一个或更多个时间间隔之间的时间中已经流过的流体容积。当发生周期性流动时,这些计算尤其有用,诸如当血液从心脏流出时。当抽运周期与将流体容积输送量与时间相关的数据结合使用时,这些计算也可以用来确定每个周期所输送的血液容积。用来确定每次心跳容积输送量的心跳速度或者其他周期可以是外生性输入,诸如由用户输入的参数或者来自外部来源的数据,诸如EKG(心电图),或者其可以基于发现周期性与非周期性流速行为的存在与性质,如上所述。容积输送量为导出参数。
本发明的另一方面是使用每次心跳的容积输送量与左心室总容积相结合,以计算射血微量。可以使用诸如查询信号等手段来估计左心室总容积,所述手段计算从查询信号所激励的区域的湍流区段返回的信号的性质。射血微量为导出参数。
为了说明的目的,以下描述针对特定医疗设备应用,然而应该理解本系统对于所述特定医疗设备应用之外的应用也同样有用并可以容易地使用。还应该理解本发明可能采用不同于所描述与所显示的其他配置。例如,在不同于所述的配置中,计算装置可以组合或分离功能单元。作为另一个例子,计算的顺序可以不同于所述地进行。
图17显示使用超声换能器1来将查询信号2引入患者3内的系统。该换能器可以是各种现有超声单元中的任意一种,并且可以以脉冲序列的形式提供信号2。根据本发明,换能器1可以置于患者3之外,例如在胸骨上切迹,并且可以使得信号瞄准升动脉。通过处理在相应于所希望的深度的时间窗期间所接收的返回信号分量可以进行深度瞄准(depthwise targeting)。
查询信号2从患者组织4反向散射,并且返回信号由接收换能器6检测。接收换能器6将返回信号5的机械能(未显示)转换为震荡电信号(未显示),该电信号被引入模拟放大与处理模块7。该模块执行各种功能,包括模拟数字转换、放大、滤波以及其他信号改进。来自模拟放大与处理模块7的输出为输入数字处理模块8的数据的数字化流。数字处理模块8使用诸如快速傅立叶变换等公知的数学处理,将该数字化信号数据转换为频谱,然后,转换为速度数据,该数据成为速度向量,数字处理模块8将该向量在数据处理模块9的存储器位置(未显示)上存储为有量纲数字速度剖面数据。数据处理模块9将该速度剖面数据(未显示)转换为表示逐块连续有量纲速度谱密度函数的数据向量(未显示)。以后将更详细地描述密度函数数据向量的其他处理。
数字处理模块8可以包含一个或更多个微处理器,以及一个或更多个数据存储装置,这使它能够进行确定VRF与导出参数所需的计算。数据存储装置包含用于查找表中的参数值,该查找表有利于避免或减少对于计算特殊函数的需要的计算。另外,诸如通道尺寸、平均流速等特定值可以被预先确定,并且存储在高速缓存或其他存储中,以与以后测量过程所求得的值相结合。图17还图解该系统具有控制模块10,其接收来自用户11的设置参数(未显示)。设置参数可以包括诸如患者参数(例如血细胞比容)以及有关何时应该开始与终止数据读取与处理的导通/断开信号等数据。控制模块10接收来自数据处理模块9的、有关所接收数据(未显示)的状态以及提高还保持高质量测量数据所需的查询信号2的变化的信号。控制模块10可选地接收生物活动12的信号,诸如与有关心脏活动的电信号。生物信号12可以用来将查询信号2同步到生物函数,或者将诸如VFR之类的、所计算的值与这些生物函数相关。当查询信号应该开始与停止时,控制模块10向换能器功率模块13发信号。换能器功率模块13提供驱动信号14给超声换能器1。数据处理模块9生成进入输出模块16的输出数据15。
输出模块16生成可视显示17与可听设备18以通知用户该系统的运行状态以及系统对患者3的查询结果。输出模块16也可以包括用于提供数据给其他设备(诸如EKG)或者处理系统的端口,例如,通过LAN(局域网)或WAN(广域网)。应该理解,以上所述的各种处理模块可以在一个或更多个计算机中实现,这些计算机位于本地或者通过网络对接,用适当的逻辑配置以执行相关算法。另外,可以将特定信号驱动与处理组件融入诸如超声系统之类的查询信号设备。
在本发明的方法的一种实施方式中,VRF与导出参数的计算如图18所示地进行。
步骤1.预先计算减少或者消除对于计算特殊函数的需求的参数。在这方面,步骤1可能包括计算并为以后使用而存储等式(43)所使用的变量A、B、C、D以及E的值。使用等式(44)至(48)计算这些变量的值。对于求和上限N的优选值为大于4的值,并且更优选的是N大于8,再优选的是N大于12。已经发现当使用Waterloo Maple,Inc.,Waterloo,Ontario,Canada所销售的数学软件包Maple预先计算变量的值时,值50能提供适当结果。然后,这些预先计算的值融入计算通道尺度、VFR以及导出参数所需的计算的软件中。
变量的存储值将用来生成将kt与τ相关的函数,而不或者少计算特殊函数。优选的是,使用τ与k的值之间的非等距间隔,τ与k用来计算所存储变量的值,所存储变量用来计算避免或者降低对于计算特殊函数的需求的函数的数据。更优选的是,使用随着τ与k的值增加而增加的τ与k间隔。更优选的是,τ与k间隔从τ与k的最小值开始指数增长。这种指数增长可以来自让每个后续值为前一个值的常数倍。选择所具有的τ的数目以及所具有的k的数目,并且还选择它们的增长因子。为了计算人体升动脉的流动,在计算中所使用的参数值为:
τ的数目:优选的是在2至100的范围内,更优选的是在5至20的范围内,更优选的是在15至25的范围内。
τ的开始值:优选的是大于-1小于1,更优选的值是大于0小于0.1,更优选的开始值在0与0.01的范围内,更优选的开始值具有0.001的量级。
每个后续τ的大小:优选的是前值大小的1.001至10倍的范围内,更优选的是前值大小的1.01至5倍的范围内,更优选的是前值大小的1.1至2倍的范围内,更优选的是前值大小的大约1.5倍。
k的数目:优选的是在2至100的范围内,更优选的是在5至20的范围内,更优选的是在15至25的范围内。
k的开始值:优选的是大于0小于100,更优选的开始值是大于1小于10,更优选的开始值在1.5与2.5的范围内,更优选的值大约为2。
每个后续k的大小:优选的是前值大小的1.001至10倍的范围内,更优选的是前值大小的1.01至5倍的范围内,更优选的是前值大小的1.1至2倍的范围内,更优选的是前值大小的大约1.125倍。
对于τ=0的特殊情况,不需要项,因为在这种情况下,kt=k。用来预选计算变量的τ的最终值取决于流体流动系统的性质。存在需要进行计算的最大的τ。优选的是小于100的最大的τ,这是因为在这个上限之外存在十分小的增加的分辨率。进一步地,限制预先计算的数据的大小允许对于相同数据存储器大小的、在较小τ处的更高分辨率。更优选的是小于10的最大的τ,更优选的是小于2的最大的τ。
K的最小可能值为2,从而k的优选下限为2。优选地,于k的最大值的下限在3与20之间。更优选地,k的最大值的下限在7与20之间。随着k变大,速度剖面变得非常平坦。优选地,k的最大值的上限为1,000,000,更优选地,k的上限为1,000,更优选地,k的上限为100。
步骤2.测量速度数据,并计算速度谱密度函数。
取得测定速度数据,并使用诸如傅立叶分析之类的手段计算速度谱密度函数。在本发明的一个实施例中,可以使用超声查询信号来取得测定速度数据。将包含在谱密度函数中的速度的数目标记为N。将密度函数的N个有序对中每一个标记为(vj,fj)。
步骤3.计算谱累积函数。
通过使用诸如辛普森规则之类的适当数值积分方法,数值地积分在密度函数下的面积,并且通过除以密度函数下的总面积来规格化N个不同的积分结果来使用密度函数数据计算累积函数。结果是N个有序对,其具有定义在区间[0,1]上的范围。将累积函数有序对的每一个标记为(vj,Fj)。图19显示(vj,Fj)数据的图的例子。
步骤4.计算速度剖面表征参数。
本实施方式中的描述假定该流动系统可以使用圆柱型通道来建模,并且压力梯度可以不改变方向地变化。
步骤4.1.计算对于速度剖面表征以及错误吸收函数的开始估计值。
对于本实施方式,需要初始值的参数为ff、μs、σs、σ、a以及k。使用以下步骤生成开始值。
σ:数值地搜索谱密度函数,直至找到相应于右侧众数(right-side mode)的拐点的速度(v或者φ)的两个值,然后采用
通过计算累积概率函数的二阶差分(其为寻找速度谱密度函数的斜率的数值对等过程)来寻找拐点,并且从大φ向小φ反向计算,寻找(1)当二阶差分变得小于0时所发生的最小值,以及(2)当在二阶差分再次取得值<0之前二阶差分大于0时的最大值。φ的这两个值将给出两个拐点,并且峰值将在这两个拐点之间。假定峰值与这些拐点等距离,并且峰的形状近似正态。这些假定意味着这两个拐点分隔开2σ。通过将这两个拐点之间的距离除以2来计算σ。
ff:假定流动速度一般不小于零,并且来自测量流动速度的噪声一般不使测量下降变得小于零,假定所有φ<0都是因为非流动元素。应该注意这种假定只是用来形成对于ff的初始估计,并且在进行最小二乘拟合之后,等式的解不依赖于此假定为真。从φ的最低值开始搜索谱密度函数,并且找到相应于左侧众数峰值的φ值。在这个值上,累积分布的值F(φ)为值(1-ff)的1/2。
μs:使用与左侧众数峰值对应的φ值。
σs:找到左侧众数峰值的拐点,并且计算
这使用了与形成σ初始猜测相同的逻辑。对于被来自流动数据的数据溢出所污染的右拐点可能具有问题,在这种情况下使用
基本的假定是左侧众数的分别可以近似为正态分布。对于正态分布拐点为距离均值的一个标准偏差。
通过进行以下步骤4.2至4.4估计v0与k的初始估计。然后,不改变v0(其将被用做a的值)与k的初始估计,并且从ff、μs、σs以及σ的初始估计开始,改进ff、μs、σs以及σ的估计。这些改进后的值为那些最小化作为实际Fn与由等式(35)所预测的值之间的差异而计算的误差平方和的值。可以使用现有技术中公知的多维数值函数最小化方法。
步骤4.2.调整累积谱函数,以去除来自缓慢运动组织区域的信号的效应。调整Fj的测定值,以去除缓慢运动组织反射的效应。对于估计值(ff,μs,σs),通过将φ=vj以及Fj代入以下等式(49),对缓慢运动组织的效应进行调整,并且计算每个(vj,Fj)的Fno_slowj,
F no _ slow = F j - ( 1 - f f ) ( lerf ( 1 2 ( - φ + μ s ) 2 σ s ) 2 + 1 2 ) f f - - - ( 49 ) ]]>
项 ( 1 - f f ) ( - lerf ( 1 2 ( - φ + μ s ) 2 σ s ) 2 + 1 2 ) ]]>去除了由缓慢运动组织所引起的偏移。除以ff将该数据再次规格化到正确的标度上。结果剩下与由σ所引起的测量误差相结合的、用于流动组织的 F j = F no _ slo w j . ]]>
步骤4.3.调整累积谱函数,以去除测量噪声的效应。
使用 F j = F no _ slo w j , ]]>通过使用等式(50),对于测量噪声的效应进行校正。
F no _ nois e j = 2 F j - 1 2 ( - F j + 1 ) ( ln ( a ) ln ( - φ j = a a ) ) ∫ 0 ( - F j + 1 ) ( ln ( a ) ln ( - φ j + a a ) ) erf ( 1 2 2 ( φ j + i ( ln ( - φ j + a a ) ln ( - F j + 1 ) ) - a ) σ ) + 1 dt ]]>
(50)
将该计算的结果标记为(vn,Fn)。Fn已经被调整去除了缓慢运动组织与测量噪声的效应。
步骤4.4.计算速度剖面参数的值。计算k与v0。在该实施方式中,这些参数为等式(11)的a与k。进行使用(vn,Fn)的最小二乘拟合,以计算k与a的值。
4.4.1.如果最大Fn>1,则再次规格化集合(vn,Fn),并继续将N个有序对的集合标记为(vn,Fn)。
4.4.2.使用该N个数据有序对,进行最小二乘拟合,以计算k与a。在这种情况下,所使用的数据为有量纲速度。因此,所计算的a的值具有单位,并且是在速度测量时所存在的最大速度。将该最大速度标记为v0。以下步骤实现计算k与v0的最小二乘拟合。
4.4.2.1.使用等式(51),以对每个(vn,Fn)计算Gn。
Gn=-Fn+1 (51)
4.4.2.2.通过求出等式(52)的根来求出k。在等式(53)至(56)中定义了变量A、B、C与D。使用数值搜索过程。优选地,该搜索过程将使用k=2的下限,以及最好大于2的上限,以及在大约20至1,000,000范围内的大数来约束该根。
0 = ( D - C ) ( Σ n = 1 N ln ( G n ) ( G n ( 1 2 k ) v n N - 2 G n ( 1 2 k ) v n A + G n ( 1 2 k ) v n B - G n ( 1 2 k ) D + G n ( 1 2 k ) C + G n k D - G n k C ) ) ( N - 2 A + B ) 2 - - - ( 52 ) ]]>
A = Σ n = 1 N G n ( 1 / 2 k ) - - - ( 53 ) ]]>
B = Σ n = 1 N ( G n ( 1 / 2 k ) ) 2 - - - ( 54 ) ]]>
C = Σ n = 1 N G n ( 1 / 2 k ) v n - - - ( 55 ) ]]>
D = Σ n = 1 N v n - - - ( 56 ) ]]>
4.4.2.3.使用等式(57),以使用(vn,Gn)的N个值以及在前面步骤中所计算的k值来计算v0。
v n = - ( Σ n = 1 N v n ) + ( Σ n = 1 N G n ( 1 / 2 k ) v n ) - N + 2 ( Σ n = 1 N G n ( 1 / 2 k ) ) - ( Σ n = 1 N ( G n ( 1 / 2 k ) ) 2 ) - - - ( 57 ) ]]>
4.5.计算误差吸收参数的值。
对于该实施方式,待计算值的参数为ff、μs、σs以及σ。使用在步骤4.4.2.3中计算的v0(其将被用做a的值)与k的值,计算最小化作为实际Fn与等式(35)所预测的值之间的差异而计算的误差的平方和的ff、μs、σs以及σ的值。可以再一次使用现有技术中公知的多维数值函数最小化方法。
4.6.计算误差平方和以及其他退出标准的值。
所使用的多维数值函数最小化方法将要求计算误差的大小以及参数ff、μs、σs、σ、a以及k中的变化量。优选的方式是计算在至少a以及k中的变化。
4.7.验证是否需要重复确定参数值的计算。
如果误差平方和一直持续显著减小,或者如果a以及k的值一直持续变化多于适当的容许值,则应该重复从步骤4.2开始的计算。如果在一次迭代到下一次迭代之间a以及k的值变化大于10%,则最好进行再一次迭代。更优选地是,如果在一次达代到下一次迭代之间a以及k的值变化大于1%,则进行再一次迭代。
图20显示正在向测定数据进行拟合的累积分布可能如何显示。随着使用更多次迭代来改进参数值,该曲线将更紧密地与数据相应。
对于上述实施方式的步骤4的改进是选择N个有序对的子集,并使用这些数据来进行初始计算。与使用数据全集相比,使用数据子集将使计算进行得更快。数据子集将产生正在计算的参数的近似值。通过在随后的计算中使用更多的或者全部数据,可以改进这些近似值。
对步骤4的进一步改进为:
a.对所有参数进行初始估计。
b.使用最小化误差平方和作为标准,通过对每个变量进行一维搜索,确定ff、μs、σs中每一个的更好的值。在这些搜索过程中,不改变a、k或σ的值。
c.使用ff、μs、σs、σ的最佳值,进行上述步骤以调整Fj,并计算k与v0(其与a相同)的新值。
d.使用最小化误差平方和作为标准,确定σ的更好的值。再一次地,不允许改变a或k的值。
e.使用最小化误差平方和作为标准,使用多维搜索来确定对于作为集合的所有参数(k除外)的改进后的值。该值集合不允许k值变化(其使用最近被计算的k值),但是其允许生成所有其他参数(包括a)的外插值。
f.使用最小化误差平方和作为标准,对每个参数使用一维搜索,以确定改进后的值(例如,返回的到该进一步改进方法的步骤b)。
步骤5.使用来自上一次速度测量与来自此次测量的速度剖面参数,计算τ。对于该实施方式,使用来自上一次速度测量的k与v0以及来自当前速度测量的k来计算τ。(如果只进行了一次速度测量,则跳过其余步骤返回步骤2)。
使用等式(42)和在等式(44)至等式(47)中所定义的变量,以及来自上一次速度测量的k与a(对于该实施方式a为v0的值)以及k的当前值(其为分配给kt的值)来计算τ。k与a的值将落入用于步骤1中所进行的计算的范围内。找到在步骤1中所使用的、高于以及低于来自测定速度分布的a的值的a的两个值。还选择在预先计算中所使用的、低于与来自测定速度分布的a的值相邻的值的a的下一个值。因此,将选择步骤1中所使用的三个a的值。选择来自步骤1的、跨越来自测定速度分布的k值的k值。选择步骤1中所使用的、与低于为实际速度分布所计算的a的值相邻的第三个k值。对于步骤1中所使用的a与k的九个组合的每一个,使用A、B、C以及D的值形成从τ值计算kt的九个等式。选择τ的假定值,并计算kt的值。使用内插多项式根据对于测定分布的a与k的实际值来计算kt的估计值。结果为相应于τ的假定值的kt值。假定τ的另一个值,并且使用相同的内插处理来计算kt的第二值。第三次重复此过程,以取得τ与kt的第三组值。这三个数据点将kt与τ相关。使用这些值来计算将kt与τ相关的内插多项式的系数。已知对于最近测定的速度剖面的kt的实际值,所以使用该内插多项式直接求解产生该已知kt值的τ值。
如果假定所使用的τ值没有跨越所计算的值,则重复该过程。如果所假定的τ值与τ的计算值相差大于10%,则重复该过程。最终结果为初始速度分布变为第二速度分布所需的τ值。
步骤6:使用τ(时间常量的值)、ρ(流体密度)、μ(流体粘度)以及t(在流体速度测量之间经过的时间)来计算通道尺度。通过使用这些变量的每一个的已知值,使用等式(4)来计算半径R。在该实施方式中,R为表征通道尺度的唯一参数。使用现在已知的其半径的值来计算通道的横截面积。
可以使用流体粘度或密度的测定值来进行计算。可替换地,这些计算可以基于缺省值,或者基于诸如性别、年龄等患者特性或者与血细胞比容的相关性的值。可以使用运动粘度而不是对于流体粘度与密度的分离的值。运动粘度比粘度变化小,这将导致降低的误差以及与其他诸如血细胞比容等因素的更好的相关性。
步骤7:计算流动流体的平均速度剖面。
对于该实施方式,使用等式(25)以及最近计算的a(其为最近计算的v0)与k的值来计算平均速度。
步骤8:计算容积流速。
将该平均流体速度乘以横截面积,以计算容积流速。
步骤9:计算压力梯度与所导出的参数。
在此时,所有数据都已知,从而可以计算压力梯度,以及其他导出参数,诸如速度变化、作为压力梯度函数的半径变化以及所输送的容积。当可以使用诸如EKG信号等外部数据时,可以计算其他导出参数,诸如将EKG信号与速度与所输送容积相关。
步骤10:显示并存储结果。
步骤11:可以重复取得速度测量并进行计算的整个过程。需要进行搜索以最小化误差平方和或者寻找根的各种计算的初始值可以基于对于参数的最近所计算的值。
上述步骤可以包括存储在存储介质上的指令。这些指令可以由处理系统检索并执行。指令的一些例子有软件、程序代码或者固件。存储介质的一些例子有存储器设备、磁带、盘片、集成电路以及服务器。当被处理系统执行时,这些指令用来指挥该处理系统根据本发明运行。名词“处理系统”指单个的处理设备或者一组可互操作的处理设备。处理系统的一些例子有集成电路与逻辑电路。本领域技术人员应该熟悉指令、处理系统以及存储介质。
本领域技术人员应该理解对于上述实施方式的各种变化落入本发明的范围之内。结果,本发明不限于上述的特定实施例与图例,而只由权利要求限定。