技术领域
本发明涉及超声诊断系统和应变分布(strain distribution)显示方 法，该超声诊断系统和应变分布显示方法允许用户利用超声诊断设备， 对组织的硬度进行定量的测量。
背景技术
超声诊断不仅应用于对组织结构的观察，而且应用于以下领域(超 声组织表征)：对诸如声速、阻尼系数等的组织内物理量进行测量，而 且根据这样测量的物理量，生成诊断图像。作为该领域的一部分，已 知这样的领域，其中对组织的硬度，即弹性特性，进行测量。因为组 织的弹性特性与病理状况有着紧密的关系，因此上述领域正在被重点 研究。例如，已知受到以下影响的组织展示出比正常组织大的硬度： 硬化性肿瘤，诸如乳腺癌、甲状腺癌等；肝硬化；等等。按照惯例， 通过触摸来检测组织的硬度。然而，通过触摸的检测具有难以进行客 观分析的缺点，需要熟练的外科医生，并具有以下限制：只能检测具 有某种尺寸或更大尺寸、且位于体表附近的受影响组织。
另一方面，已知这样一种方法，其中向组织施加静压力，以便使 组织压缩和变形，并利用超声波来测量与施加的压力相对应的组织内 应变(strain)，以便估计组织的弹性特性(J.Ophir，I.Cespedes，H. Ponnekanati，Y.Yazdi，and X.Li，“Elastography：A quantitative method for imaging the elasticity of biological tissue”，Ultrasonic Imaging，Vol. 13，pp.111-134，1991)。常规技术是基于以下事实被开发的：具有大硬 度的组织在压力下展示出小的组织应变，以及另一方面，具有小硬度 的组织在压力下展示出大的组织应变。也就是，关于上述常规方法， 向组织施加静压力，并根据在这样施加的压力下的组织内的应变分布， 来估计组织的弹性特性。
特别是，在不通过超声波探头向组织施加压力的情况下，利用带 有超声波探头的超声诊断设备，来进行超声回波信号(在没有施压的 情况下的射频(RF)信号)的正常测量。随后，外科医生通过超声波 探头、轻微地向组织施加压力(大约百分之几)，此后测量穿过被施压 的组织的超声回波信号(压力下的RF信号)。然后，利用空间相关方 法，根据在向组织施压和没有施压的情况下的RF信号，来估计位移 分布，该位移分布代表由于这样施加的压力而造成的组织的每一点的 位移。
空间相关方法具有以下机理：其中，通过利用二维相关函数的模 板匹配(mask matching)，根据在向组织施加压力和没有施加压力的 情况下的RF信号(或RF信号的包络信号)，来估计在施加的压力下 的组织内的位移分布。即，将具有某种尺寸的二维相关窗口(模板) 应用于和没有向组织施压的情况下的层析成像数据相对应的RF信号 数据，以便通过利用自相关处理、检测在应用了相关窗口的RF信号 数据与向组织施压的情况下的RF信号数据之间的最大相关，来估计 二维表面上的期望测量点的位移。例如，对以网格的形式设置的每个 测量点，都执行上述自相关处理，由此估计应变分布。一般来说，由 于水平方向的采样比轴向的采样粗糙，使得对于利用空间相关方法的 处理，其水平方向(超声波束的扫描方向)的位移检测精度比轴向的 位移检测精度差。如上所述，空间相关方法具有使用户能够估计二维 位移向量的优点。而且，虽然上述空间相关方法具有位移估计精度受 采样间隔的限制的缺点，但是空间相关方法具有以下优点：即使在组 织发生大的变形(例如，大约5％)的情况下，也使用户能够估计位 移分布。然而，和常规超声诊断不同，空间相关方法具有以下缺点： 空间相关处理耗费大量计算，导致了难以实时地处理。
因此，本发明的目的是提供一种用于实时地获得位移分布、应变 分布和弹性模量分布(elastic modulus distribution)的方法。
发明内容
为了解决上述问题，根据本发明的超声诊断系统根据在施加压力 和没有施加压力的情况下、利用超声波探头测量对象(subsject)而获 得的反射回波信号(RF信号)，来获得对象的组织的位移，该超声诊 断系统包括：存储装置，用于存储利用所述超声波探头检测的、诸如 信号的包络信号的信号特性；相关计算装置，用于根据在向对象施加 压力和没有施加压力的情况下所述存储装置中存储的所述特性，计算 在向所述对象施压和没有施压的情况下的所述特性之间的相关系数， 以及在施压和没有施压的情况下的接收信号之间的相位差；以及位移 计算装置，用于根据由相关计算装置这样获得的、在施压与没有施压 情况下的RF信号之间的相关系数和相位差，来计算由于所述施加压 力而造成的每个测量点的位移。此外，根据本发明的超声诊断系统包 括：应变计算装置，用于通过对每个测量点的位移进行空间微分，来 计算对象的组织的应变分布；以及显示装置，用于显示这样获得的应 变分布。
如上所述，关于根据本发明的超声诊断系统，根据在施压与没有 施压情况下的诸如包络信号的特性之间的相关，来计算每个测量点的 位移，由此能够实时地估计位移分布。此外，关于根据本发明的超声 诊断系统，可以通过在所述超声波束方向上、以所述超声信号的二分 之一波长的间距改变所述测量点，来获得在向测量点施压和没有施压 情况下的包络信号之间展示出最大相关系数的每个测量点的位置，由 此解决作为Doppler(多普勒)方法的缺点的混叠(aliasing)问题。
注意，对于这样的相关计算，其中通过沿时间轴改变在施压情况 下的包络信号，以获得在施压的情况下每个测量点的包络信号的自相 关函数，来计算在向测量点施压和没有施压的情况下的包络信号之间 展示出最大相关系数的每个测量点的位置，该相关计算导致了长计算 时间，经常造成不能进行实时计算的问题。
因此，关于根据本发明的超声诊断系统，首先，优选地计算在施 压和没有施压情况下的包络信号的自相关函数。然后，通过改变对应 于预定的测量点变化的、这样获得的自相关函数之间的相位差，例如 通过以超声信号的二分之一波长的间距改变自相关函数之间的相位 差，来获得相关系数。这减少了用于位移计算的计算时间，由此能够 进行高速处理。此外，根据本发明的超声诊断系统可以进一步包括： 存储装置，用于存储正交检波的RF信号的包络信号；相关计算装置， 用于计算在向测量点施加压力和没有施加压力的情况下的包络信号之 间展示出最大相关系数的每个测量点的位置，所述测量点被二维相关 窗口包围；以及位移计算装置，用于根据由相关计算装置这样获得的、 展示相关系数和相位差的每个测量点位置，来计算由于施加压力而造 成的每个测量点的至少二维位移。
也就是，在本说明书中提出了被称为“组合自相关方法(CA法)” 的方法。如上所述，组合自相关方法具有以下优点：在利用相关窗口 的空间相关方法中的二维和三维位移测量；以及在Doppler(多普勒) 方法中的实时和高精度计算。组合自相关方法允许用户估计位移分布， 而在一定程度上与水平方向的位移无关。在该情况下，二维方向可以 包括：超声波束方向，超声波探头在该超声波束方向接收超声信号； 以及超声波束扫描方向。此外，关于根据本发明的超声诊断系统，优 选地，通过在超声波束方向上以超声信号的二分之一波长的间距、以 及在超声波束扫描方向上以超声波束间距，改变测量点，来检测展示 出最大相关系数的每个测量点的位置。注意，虽然根据本发明，在超 声波束方向上的测量点的改变间距不限于超声信号的二分之一波长， 但是优选地采用小于超声信号的二分之一波长的间距。
为了进一步提高位移计算的计算速度，首先，优选地计算在施压 和没有施压的情况下的包络信号的自相关函数。然后，通过在对应于 预定的测量点变化的范围中，改变自相关函数之间的相位差，来获得 展示出最大相关系数的每个测量点的位置。
此外，根据本发明的位移计算不仅可以应用于二维计算，而且可 以应用于三维计算。关于采用具有一维阵列结构的超声波探头的三维 配置，存储装置中存储的帧数据包括由多个帧数据集组成的体积数据， 多个帧数据集的每一个用作切片数据。另一方面，关于采用具有二维 阵列结构的超声波探头的三维配置，数据包含通过在切片方向上扫描 而获得的包络信号。相关计算装置通过在三维方向上、关于体积数据， 改变被三维相关窗口包围的测量点，来检测在向测量点施压和没有施 压情况下的包络信号之间展示出最大相关系数的每个测量点的位置， 所述测量点被三维相关窗口包围。同时，相关计算装置计算在施压和 没有施压的情况下的RF信号之间的相位差。在该情况下，三维方向 可以包括：超声波束方向，超声波探头在该超声波束方向接收超声信 号；超声波束扫描方向；以及与上述两个方向正交的切片方向。此外， 相关计算装置优选地计算，在超声波束方向、超声波束扫描方向、以 及与上述两个方向正交的切片方向上施加压力和没有施加压力的情况 下的RF信号之间的相位差。此外，上述的高速处理方法可以应用于 三维配置。此外，可以通过以超声波束的切片间距为单位、在切片方 向上改变测量点，来进行上述计算。
此外，根据本发明用于获得弹性模量分布的方法可以包括弹性模 量计算装置，该弹性模量计算装置用于：通过将对象划分为有限数量 的单元(element)，来建立至少二维或三维的有限元模型；以及根据 用于建立模型的信息和这样获得的应变分布，来计算弹性模量分布。 此外，可以利用显示装置来显示这样获得的弹性模量分布。在该情况 下，弹性模量计算装置优选地通过，在对象的组织展示出各向同性的 弹性和接近不可压缩性的假定下，将对象的组织划分为有限数量的矩 形平行六面体单元，来建立三维有限元模型。此外，弹性模量计算装 置优选地在对象的组织的每个单元都展示出均匀的弹性模量、均匀应 力和均匀应变的假定下，利用弹性方程、根据关于上述应变分布的信 息，来计算弹性模量分布。
如上所述，在组织展示出各向同性的弹性的假定下，进行根据本 发明的计算。其原因在于，在施加于组织的外部静压力下，应力与应 变之间的关系展示出线性关系。因此，可以在组织用作弹性模型的假 定下，进行近似计算。另外，组织展示出各向同性的特性，以及因此， 可以在组织用作各向同性的弹性模型的假定下，进行根据本发明的计 算。另一方面，关于本发明，在组织展示出接近不可压缩性的假定下， 进行计算。其原因在于，如果在组织展示出完全不可压缩性的假定下， 进行计算，即利用0.5的组织内恒定Poisson比(泊松比)进行计算， 则弹性方程变成特殊情况，造成了不能利用有限元方法进行计算的问 题。注意，关于本发明，利用均匀的组织内泊松比进行计算。在该情 况下，只应估计Young模量(杨氏模量)，以估计弹性模量分布，由 此简化了反演问题。注意，泊松比展示出足够的组织内均匀性，以及 因此，关于本发明，优选地利用0.49的泊松比进行计算。根据本发明 的弹性模量分布计算能够仅仅根据可以高精度地计算的轴向应变分 布，来重构弹性模量分布，由此能够稳定地计算弹性模量分布。
附图说明
图1描述了超声诊断设备的机理。
图2显示了通过施加静压力的组织弹性测量方法的特定例子，以 及通过施加静压力的组织弹性测量方法的机理。
图3显示了空间相关方法的机理。
图4显示了Doppler(多普勒)方法的机理。
图5显示了组合自相关方法的机理。
图6显示了用于执行组合自相关方法的基本算法的电路配置的框 图。
图7显示了根据本发明实施例的超声诊断系统的示例配置的框 图。
图8显示了三维组合自相关方法的基本算法的流程图。
图9所示为，根据本发明的超声诊断系统中所采用的三维组合自 相关方法的基本算法的流程图，以及图7所示处理的一部分的详细流 程图。
图10显示了，用于详细描述具有提高的计算速度的组合自相关方 法的流程图，该组合自相关方法在图9所示的步骤S15中被执行。
图11显示了，用于执行根据本发明的超声诊断系统中所采用的三 维组合自相关方法的基本算法的电路配置的框图。
图12显示了模拟过程的示意图。
图13显示了，在5.0MHz的超声波中心频率的情况下每一点的点 扩展函数的例子。
图14显示了组织模型的示意图。
图15显示了，利用各估计方法估计的、由于水平方向位移造成的 估计应变误差。
图16显示了，在0.0mm的水平方向位移的情况下利用各估计方 法(组合自相关方法、扩展组合自相关方法和空间相关方法)估计的 应变分布。
图17显示了，在0.4mm的水平方向位移的情况下利用各估计方 法(组合自相关方法、扩展组合自相关方法和空间相关方法)估计的 应变分布。
图18显示了，用于模拟倾斜方向的压缩的组织模型。
图19显示了，通过对轴向的组织模型的简单压缩进行模拟而获得 的应变分布的估计结果。
图20显示了，通过对倾斜方向的组织模型压缩进行模拟而获得的 应变分布的估计结果。
图21显示了两种组织模型例子，该两种组织模型例子的每一种都 具有三维结构。
图22是显示了利用内含物包含模型估计的结果的第一图。
图23是显示了利用内含物包含模型估计的结果的第二图。
图24是显示了利用层状结构模型估计的结果的第一图。
图25是显示了利用层状结构模型估计的结果的第二图。
图26显示了超声诊断系统的基本配置。
图27显示了，超声诊断系统中采用的超声扫描仪的特定配置。
具体实施方式
以下将参考附图，来描述根据本发明实施例的超声诊断系统。根 据本发明的超声诊断系统采用一种被称为“扩展组合自相关方法”的 方法，其中利用组合自相关方法、通过一维窗口、根据包络信号进行 相关计算的一维检波处理被扩展为，通过二维相关窗口、利用二维检 波，来处理水平方向的位移。此外，利用根据实施例的扩展组合自相 关方法，只对轴向间距为超声波波长的一半、水平方向间距为束线间 距(beam-line pitch)的网格点，执行包络相关计算，以便减小计算量， 由此使得能够进行高速计算。注意，根据本实施例的扩展组合自相关 方法如同组合自相关方法一样，采用相位信息来提高轴向位移的估计 精度。然而，由于缺少用作载体(carrer)的信号，从而不利用相位信 息来估计水平方向的位移。因此，如同空间相关方法一样，水平方向 位移的估计精度受限于采样间隔(超声波束线间距)。注意，根据本实 施例的扩展组合自相关方法没有利用特殊的机理来提高水平方向位移 的估计精度，因为可以利用后面描述的弹性模量分布重构方法、仅仅 根据轴向应变(位移)分布，来估计弹性模量分布。
在描述根据本实施例的扩展组合自相关方法的特定配置之前，将 参考图1至图6来描述作为扩展组合自相关方法的基础的组合自相关 方法。图1描述了超声诊断设备的原理。由图1可以清楚地看出，用 作超声波探测器的超声波探头具有将电信号转换成超声波以及将超声 波转换成电信号的功能，这允许用户向组织11投射超声波脉冲信号。 穿过组织11的超声波脉冲信号的一部分在具有互不相同的声阻抗的 区域之间的第一界面12上被反射。将被称为“反射回波信号12a”的 反射超声波穿过组织11、投向超声波探头10，其余的超声波穿过第一 界面12。穿过第一界面12的超声波脉冲信号的一部分在具有互不相 同的声阻抗的区域之间第二界面13上被反射。同样地，在第二界面 13上反射的、将被称为“反射回波信号13a”的反射超声波脉冲信号 穿过组织投向超声波探头10，另一方面，其余的超声波脉冲信号穿过 第二界面13。反射超声波回波信号被超声波探头10接收，以便被转 换为反射回波信号，该反射回波信号是电信号。在该情况下，用以下 表达式(1)来表示，从超声波探头10投射超声波脉冲信号、直到收 到从与超声波探头10距离L的反射物14(具有互不相同的声阻抗的 区域之间的界面)反射的回波信号为止的时段t。
t = 2 L c - - - ( 1 ) ]]>
在此，c表示组织内的声速，在穿过软组织时，c可以被确定为大 约1500米/秒的常量。因此，根据从投射超声波直到收到反射回波信 号的时间t，来计算探头与反射物之间的距离L。此外，反射回波信号 具有关于组织的声学特性的信息，并且因此可以根据反射回波信号， 将诸如B模式层析图像的组织信息图像显示在监视器上。
例如，已知这样一种方法，其中利用超声诊断设备来测量代表组 织硬度的弹性特性。上述用于测量弹性特性的方法有一种机理，其中 将机械振动作用于组织上，并且基于横波(transverse wave)以高传播 速度穿过硬组织、并且以低传播速度穿过软组织这一事实，根据这样 产生的横波的传播速度，来估计硬度信息。严格地讲，传播穿过组织 的横波的传播速度v取决于组织ρ的密度、剪切模量μ1、剪切粘度μ2和振动的角频率ω，如下式(2)所表示。
v = 2 ( μ 1 2 + ω 2 μ 2 2 ) ρ ( μ 1 + μ 1 2 + ω 2 + μ 2 2 ) - - - ( 2 ) ]]>
另一方面，如图2所示，已提出了这样一种方法，其中将静压力 施加于组织上，并且根据在所施加的静压力下的组织应变分布，来估 计组织的弹性特性。上述方法是基于以下事实而设计的：这种静压力 在硬组织内造成小应变，而在软组织内造成大应变。图2(A)显示了通 过施加静压力的组织弹性测量方法的特定例子。图2(B)显示了通过施 加静压力的组织弹性测量方法的机理。由图可以清楚地看出，关于上 述方法，利用带有超声波探头10的常规超声诊断设备、在不施加压力 的情况下，正常地测量组织11的超声回波信号(在不施加压力的情况 下的RF信号)。随后，通过超声波探头10轻微地(大约百分之几) 按压组织11，并测量在施加压力下的组织11的超声回波信号(在施 加压力的情况下的RF信号)。然后，根据在向组织施加压力和没有施 加压力的情况下测量的RF信号，来估计代表受压时的组织内的每一 点的位移的位移分布。位移分布估计方法的主要例子包括利用空间相 关的方法和利用Doppler(多普勒)效应的方法。
图3显示了空间相关方法的机理。关于该方法，通过利用二维相 关函数进行模板匹配，根据在向组织施加压力和没有施加压力的情况 下测量的RF信号(或RF信号的包络信号)，来估计在施加压力下的 组织内的应变分布。以下将描述该方法的特定处理。首先，分别以i1(t， x)和i2(t，x)来表示在向组织施加压力和没有施加压力的情况下测量的 RF信号(或RF信号的包络信号)，则上述两个信号之间的交叉相关 系数C(t，x；n，m)以下式(3)来表示。
C ( t , x , n , m ) = Σ v = - t 0 / 2 t 0 / 2 Σ w = - x 0 / 2 x 0 / 2 i 1 ( t + v , x + w ) i 2 ( t + v + n L t , x + w + m L x ) Σ v = - t 0 / 2 t 0 / 2 Σ w = - x 0 / 2 x 0 / 2 i 1 2 ( t , x ) · Σ v = - t 0 / 2 t 0 / 2 Σ w = - x 0 / 2 x 0 / 2 i 2 2 ( t + v + n L t , x + w + m L x ) - - - ( 3 ) ]]>
在此，t代表超声波束方向(轴向)的坐标点，x代表与超声波束 方向正交(水平方向)的坐标点，t0代表相关窗口的轴向尺寸，x0代 表相关窗口的水平方向尺寸，Lt代表轴向的采样间隔，以及Lx代表水 平方向的采样间隔。此外，n和m为整数。在该情况下，如果将显示 出最大值的交叉相关系数的组合(n，m)表示为(k，l)，则在测量点(t，x) 处的轴向位移uy和水平方向位移ux分别以下式来表示。
uy＝kLt
ux＝lLx
注意，在水平方向的采样间隔Lx比轴向采样间隔Lt粗糙的情况 下进行数据采样，由此导致了水平方向的位移检测精度比轴向位移检 测精度差。对每个测量点进行上述处理，由此估计位移分布。空间相 关方法的优点在于，使用户能够估计二维位移向量分量。此外，即使 组织内出现大的应变(大约5％)，空间相关方法也允许用户估计位移 分布。然而，不同于常规超声测量系统，上述方法引起大量的计算， 造成实时计算的困难。此外，位移估计精度受采样间隔的限制，由此 造成了与Doppler(多普勒)方法相比、精度较差的问题，将在后面 描述Doppler(多普勒)方法。
图4显示了Doppler(多普勒)方法的机理。关于该方法，根据 在向组织施加压力和没有施加压力的情况下的RF信号，利用Doppler (多普勒)效应来估计在施加压力下的组织内的位移分布，Doppler (多普勒)效应也用于血流测量中。以下将描述该方法的特定处理。 首先，利用如下式(4)所表示的模型，来表示在向组织施加压力和没 有施加压力的情况下的RF信号。
在此，i1(t)代表在没有施加压力情况下的RF信号，i2(t)代表在施 加压力情况下的RF信号，A(t)代表包络信号，ω0代表超声波的中心 角频率，以及τ代表时移(time shift)。一旦对两个RF信号的每一个执 行正交检波，就获得了基带信号，如下式所示。
s1(t)＝A(t)ejθ
s 2 ( t ) = A ( t - τ ) e j ( ω 0 τ + θ ) - - - ( 5 ) ]]>
同样，以下式来表示上述两个信号之间的复相关函数R12(t)。
R 12 ( t ) = ∫ - t 0 / 2 t 0 / 2 s 1 ( t + v ) s 2 ( t + v ) * dv = R A ( t ) e - j ω 0 τ - - - ( 6 ) ]]>
在此，RA(t)代表包络信号的自相关函数，以及t0代表超声波束轴 向的相关窗口尺寸。此外，星号“*”代表复共轭算子。因此，从相关 函数R12(t)的相位φ(t)，获得了由于施加压力而造成的时移τ和轴向位移 uy，如下式(7)所示。
τ = - φ ( t ) ω 0 ]]>
u y = cτ 2 - - - ( 7 ) ]]>
注意，c代表组织内的声速，并且假定在组织内c为常量。
利用Doppler(多普勒)方法，对每一个测量点执行上述处理， 并且通过与根据Doppler(多普勒)效应开发的血流测量一样的方式， 估计位移分布。因此，Doppler(多普勒)方法具有实时测量的优点。 此外，Doppler(多普勒)方法利用相位信息进行计算，由此与空间相 关方法相比，可以更高精度地估计位移。然而，Doppler(多普勒)方 法的缺点在于：大位移的测量，例如超声波中心频率的波长的四分之 一或更大的位移的测量，造成了混叠(aliasing)，导致了不能估计正 确位移的问题。此外，Doppler(多普勒)方法的缺点还在于，不能够 象可以从上式所推定的那样，来估计除轴向以外的位移。
为了解决上述问题，本发明人提出了“组合自相关方法(CA法)”， 该方法同时具有上述两种方法的优点。图5显示了由本发明人提出的 组合自相关方法的机理。关于组合自相关方法，利用RF信号的包络 信号的相关进行计算，由此解决了作为Doppler方法的缺点的混叠问 题。以下将描述该方法的特定处理。
首先，如同Doppler方法的情况一样，利用如下式所表示的模型， 来表示在向组织施加压力和没有施加压力的情况下的RF信号。
在此，i1(t)代表在没有施加压力情况下的RF信号，i2(t)代表在施 加压力情况下的RF信号，A(t)代表包络信号，ω0代表超声波的中心 角频率，以及τ代表时移。一旦对两个RF信号的每一个执行正交检波， 就获得了基带信号，如下式所示。
s1(t)＝A(t)ejθ
s 2 ( t ) = A ( t - τ ) e j ( ω 0 τ + θ ) - - - ( 9 ) ]]>
然后，利用下式来表示上述两个信号之间的复相关函数R12(t；n)。
R 12 ( t ; n ) = ∫ - t 0 / 2 t 0 / 2 s 1 ( t + v ) s 2 ( t + n T 2 + v ) * dv = R A ( t ; τ - n T 2 ) e - j ω 0 ( τ - n T 2 ) , ( n = . . . , - 2 , - 1,0,1,2 . . . ) - - - ( 10 ) ]]>
在此，T代表超声波的周期，RA(t；τ)代表包络信号的自相关函数， 以及t0代表相关窗口尺寸。此外，星号“*”代表复共轭算子。在此， n代表变量号，并且对不同的n进行每次计算。在由变量号确定的点 的周围，计算每个测量点处的位移。注意，在n＝0的情况下，组合自 相关函数和如式(6)所示的Doppler(多普勒)方法中的自相关函数 相等。也就是，在对超声波波长的四分之一或更大的位移进行测量的 情况下，n等于0的组合自相关函数造成了混叠问题。为了解决上述 问题，关于组合自相关方法，将包络相关系数C(t；n)定义为如下式(11) 所表示。
C ( t ; n ) = | R 12 ( t ; n ) | | R 11 ( t ; 0 ) | · | R 22 ( t ; n ) | - - - ( 11 ) ]]>
注意，R11(t；0)代表s1(t)的自相关函数，以及R22(t；n)代表s2(t+nT/2) 的自相关函数。如果将显示出最大值的包络相关系数C(t；n)的n表示 为k，则利用n＝k的R12(t；k)的相位φ进行计算。在该情况下，无混叠 地计算位移。其原因在于，在二分之一波长的间距下进行包络相关计 算。注意，二分之一波长的计算间距是用于计算位移、同时防止混叠 的最大间距。从而，利用φ(t；k)来计算时移τ和轴向位移uy，如下式所 示。
τ = - φ ( t ; k ) ω 0 + k T 2 ]]>
u y = cτ 2 - - - ( 12 ) ]]>
注意，c代表组织内的声速，并且假定在组织内c为常量。
利用组合自相关方法，对每个测量点都执行上述处理，由此估计 位移分布，它是Doppler(多普勒)方法的扩展方法。因此，组合自 相关方法具有实施测量的优点。此外，和Doppler(多普勒)方法不 同，组合自相关方法具有以下优点：使用户能够利用包络相关来估计 包含大位移(即，超声波波长的四分之一或更大的位移)的位移分布。
图6显示了用于执行组合自相关方法的基本算法的电路配置的框 图。在图6中，在没有施加压力的情况下获得的回波信号x(t)被输入 到未受压正交检波电路(QD)131，用于正交检波，并且这样检测的 正交检波信号Ix(t)和Qx(t)被输入到第一相关计算电路133和第一相关 系数计算电路1350至135N。另一方面，在施加压力的情况下获得的 回波信号y(t)被输入到第一受压正交检波电路(QD)1320，用于正交 检波，并且这样检测的正交检波信号Y(t)＝Iy+jQy(Iy(t)，Qy(t))被输入 到第一相关系数计算电路1340和第二相关计算电路1350。第一延迟 电路134使回波信号y(t)延迟超声波的周期T，并且延迟的回波信号 y1＝y(t-T)被输入到第二受压正交检波(QD)电路1321。同样，第二延 迟电路135使已经被第一延迟电路134延迟的回波信号y1＝y(t-T)延迟 超声波的周期T，并且延迟的回波信号y2＝y(t-2T)被输入到下一个第二 受压正交检波(QD)电路1322(未显示)。注意，该电路具有N个延 迟电路，回波信号被连续延迟，并且被延迟了整数倍周期T的回波信 号以同样的方式被输入到相应的受压正交检波电路。
第一相关计算电路133根据信号Ix和Qx计算相关值Rxx，并将 相关值Rxx输出到第二相关系数计算电路1380至138N。第二相关计 算电路1340从受压正交检波电路1320接收正交检波信号Iy(t)和 Qy(t)，根据信号Iy和Qy计算相关值Ryy，并将相关值Ryy输出到第 二相关计算电路1380。第一相关系数计算电路1350从未受压正交检 波电路131接收正交检波信号Ix(t)和Qx(t)，并从第一受压正交检波电 路1320接收正交检波信号Iy(t)和Qy(t)，据此计算复基带信号SR和 SI，并将基带信号SR和SI输出到第三相关计算电路1360和相位差计 算电路1370。第三相关计算电路1360从第一相关系数计算电路1350 接收复基带信号SR和SI，据此计算相关值|Rxy|，并将计算的相关值 |Rxy|输出到第二相关系数计算电路1380。相位差计算电路1370从第 一相关系数计算电路1350接收复基带信号SR和SI，并根据复基带信 号SR和SI计算相位差φ0(t)。第二相关系数计算电路1380从第一相关 计算电路133接收相关值Rxx，从第三相关计算电路1360接收相关值 |Rxy|，以及从第二相关计算电路1340接收相关值Ryy，根据上述相关 值计算相关系数C0(t)，并输出计算的相关系数C0(t)。
第二受压正交检波电路(QD)1321接收被第一延迟电路134延 迟的回波信号y1＝y(t-T)，并将正交检波信号Y1(t)＝Iy1+jQy1(Iy1(t)， Qy1(t))输出到第一相关系数计算电路1341和第二相关计算电路1351。 第二相关计算电路1341从第二受压正交检波电路(QD)1321接收正 交检波信号Iy1(t)和Qy1(t)，根据信号Iy1(t)和Qy1(t)计算相关值Ry1y1， 并将相关值Ry1y1输出到第二相关计算电路1381。第一相关系数计算 电路1351从未受压正交检波电路131接收正交检波信号Ix(t)和Qx(t)， 以及从第二受压正交检波电路(QD)1321接收正交检波信号Iy1(t)和 Qy1(t)，据此计算复基带信号SR1和SI1，并将复基带信号SR1和SI1输 出到第三相关计算电路1351和相位差计算电路1371。第三相关计算 电路1361从第一相关系数计算电路1351接收复基带信号SR1和SI1， 根据复基带信号SR1和SI1计算相关值|Rxy1|，并将计算的相关值|Rxy1| 输出到第二相关系数计算电路1381。相位差计算电路1371从第一相 关系数计算电路1351接收复基带信号SR1和SI1，并根据复基带信号 SR1和SI1计算相位差φ1(t)。第二相关系数计算电路1381从第一相关计 算电路133接收相关值Rxx，从第三相关计算电路1361接收相关值 |Rxy1|，以及从第二相关计算电路1341接收相关值Ry1y1，根据上述相 关值计算相关系数C1(t)，并输出计算的相关系数C1(t)。
同样地，从自第一延迟电路135向下的相应延迟电路接收信号的 第二受压正交检波电路(QD)1322至132N的每一个，第二相关计算 电路1342至134N的每一个，第一相关系数计算电路1352至135N的 每一个，第三相关电路1362至136N的每一个，相位差计算电路1372 至137N的每一个，以及第二相关系数计算电路1382至138N的每一 个，都执行如同如上所述的第一级和第二级电路元件一样的处理，由 此输出相关系数C2(t)至CN(t)以及相位值φ2(t)至φN(t)。如上所述，图6 所示的用于执行组合自相关方法的基本算法的电路具有以下配置：其 中对于延迟电路134至13N的每一个，都使在施加压力下的回波信号 y(t)延迟了超声波的二分之一周期T/2(二分之一波长)的时段，并且 相应的正交检波电路(QD)1320至132N对这样延迟的每一个回波信 号进行正交检波(quadrature-detection)。
如上所述，通过在向组织施加压力的情况下的估计位移分布的空 间微分，获得了应变分布。应变分布代表组织的相对弹性特性，并且 因此，基于应变分布的诊断展示出与基于弹性模量分布的诊断相似的 效果。然而，在造成整个受影响组织的硬化的肝硬化的情况下，难以 进行如同弹性特性分布一样的组织诊断，弹性特性分布允许外科医生 进行定量估计。因此，近些年来，用于组织弹性模量分布的重构方法 正在被研究。注意，所有这些重构方法都处于研究阶段，并且还没有 建立标准方法。
另一方面，可以基于如上所述的、组织内的应变分布和应力分布， 来获得组织弹性模量分布。然而，难以利用现有技术直接测量应力分 布。因此，尤其是基于应变分布来重放弹性模量分布，以满足在向组 织施加压力情况下的边界条件，即需要解决反演问题。一般来说，难 以解决反演问题，并且只提出了很少的弹性模量重构方法。以下将描 述常规的弹性模量重构方法。
首先，已知这样一种方法，该方法是基于通过一维模型(一维弹 性模型)来表示组织的假定而提出的。也就是，已知这样一种方法， 其中基于通过一维弹性模型来表示组织的假定，来计算弹性模量分布。 基于上述假定，将弹性模量确定为应变的反数(inverse number)。严 格地讲，上述方法不是弹性模量重构方法。关于上述方法，只获得了 应变的倒数，并且因此，如同应变分布的情况一样，只能获得组织的 相对弹性特性。
其次，已知这样一种方法，其中弹性方程被简化为没有应力项(假 定组织显示出各向同性弹性、不可压缩性以及平面应变)。关于上述方 法，所形成的、以代表平面应力状态的弹性方程被简化为没有应力项， 并且利用没有应力项的简化弹性方程、基于应力分布(包括剪切应变 分量的应变张量的所有分量)来重构组织弹性模量分布，以满足边界 条件(体表上的施压分布，或者体表上的位移分布)。注意，该方法需 要这样的区域(参考区)，在该区域已获得了弹性模量。
第三，已知这样一种方法，其中结合了弹性微分方程(假定组织 显示出各向同性弹性、不可压缩性以及平面应变)。关于上述方法，所 形成的、以代表平面应力状态的弹性方程被简化为没有应力项，并且 通过连续地把关于弹性模量的无应力项的简化微分方程与作为参考的 体表附近弹性模量结合、基于应变分布(包括剪切应变分量的应变张 量的所有分量)，来重构组织弹性模量分布。注意，该方法需要这样的 体表附近区域(参考区)，在该体表附近区域已预先获得了弹性模量分 布。此外，上述方法具有以下问题：由于与作为参考的体表附近弹性 模量结合而造成了误差积累，使得离体表越远，计算误差就越大。
第四，已知一种利用扰动法的方法(假定组织显示出各向同性弹 性、不可压缩性以及平面应变)。关于上述方法，利用带有迭代法的扰 动法，基于超声波束方向(轴向)的体表施压分布和体表应变分布， 通过解已被形成以代表平面应力状态的弹性方程，来重构组织弹性模 量分布。
已经描述了基本机理和要解决的特定问题。以下将描述根据本发 明的实施例，以便解决上述问题。图7显示了根据本发明实施例的超 声诊断系统的简要配置框图。超声波探头91包括常规的扇形扫描探头 (扇形调相阵列探头)、线性扫描探头(线性阵列探头)或凸形扫描探 头(凸形阵列探头)，具有向要观察的组织投射超声波、以及接收反射 的超声波的功能。
在向组织施加压力和没有施加压力的情况下获得的RF信号从超 声波探头91输出到正交检波器92。正交检波器92把向组织施加压力 和没有施加压力的情况下的RF信号转换为施加压力和没有施加压力 的复包络信号(IQ信号)，并将IQ信号输出到二维复相关计算单元 93。二维复相关计算单元93计算在向组织施加压力的情况下的RF信 号与没有向组织施加压力的情况下的RF信号之间的二维相关，将显 示出最大相关的位置输出到水平方向位移计算单元94和轴向位移计 算单元95，并将相应的相关函数相位输出到轴向计算单元95。注意， 在轴向、以二分之一超声波中心频率的间距来计算相关，该间距是获 得相位、同时防止混叠的最大间距。以这种间距来计算相关，以便允 许超声诊断系统的实时显示。因此，本发明不限于以二分之一波长的 间距进行计算，而是，可以构造这样的配置，其中高精度地计算相关。
水平方向位移计算单元94根据与从二维复相关计算单元93收到 的水平方向最大相关相对应的位置，计算水平方向位移ux，并将该位 移输出到水平方向应变计算单元96。另一方面，轴向位移计算单元95 根据与从二维复相关计算单元93收到的轴向最大相关和相位相对应 的位置，计算轴向位移uy，并将该位移输出到轴向应变计算单元97。 水平方向应变计算单元96对从水平方向计算单元94收到的水平方向 位移ux进行空间微分，以便计算水平方向的应变分布εx，并将该应变 分布输出到量化单元98。另一方面，轴向应变计算单元97对从轴向 计算单元95收到的轴向位移uy进行空间微分，以便计算轴向应变分 布εy，并将该应变分布输出到量化单元98。量化单元98对水平方向 应变分布εx和轴向应变分布εy进行量化，以便对这些应变分布进行灰 阶显示(或彩色显示)，并在显示单元99上显示信息。显示单元99 显示这样量化的每一个应变分布。
接下来，将描述图7所示的超声诊断系统中采用的扩展组合自相 关方法的操作。首先，考虑这样一种情况，其中组织被轻微地压缩(即 百分之几或更少)。在该情况下，从局部观点看，假定压力造成了组织 内的每一点的平行位移。也就是，通过下式所示的模型，来表示向组 织施加压力和没有施加压力的情况下的RF信号。
在此，i1(t，x)代表没有施加压力情况下的RF信号，i2(t，x)代表施 加压力情况下的RF信号，A(t，x)代表包络信号，ω0代表超声波的中 心角频率，τ代表时移，该时移用作代表轴向位移的时间参数，ux代表 水平方向位移，以及θ代表初始相位。注意，和Doppler(多普勒)方 法和组合自相关方法不同，关于该方法，施加压力和没有施加压力情 况下的RF信号由考虑到水平方向位移的模型来表示。要在最后阶段 获得的参数是轴向位移uy＝cτ/2(即时移τ)和水平方向位移ux。注意， c代表组织内的声速，并且假定在组织内c是常量。
然后，正交检波器92对向组织施加压力和没有施加压力的情况下 的RF信号进行正交检波。也就是，将具有与超声波中心频率相同的 频率的正弦波和余弦波应用于RF信号，随后进一步对RF信号进行低 通滤波，由此获得复基带信号s1和s2，如下式(14)所示。
s1(t，x)＝A(t，x)ejθ
s 2 ( t , x ) = A ( t - τ , x - u x ) e j ( ω 0 τ + θ ) - - - ( 14 ) ]]>
然后，将s1(t，x)和s2(t+nT/2，x+mL)之间的二维复相关函数R12(t，x； n，m)定义为如下式(15)所示。
R 12 ( t , x ; n , m ) = ∫ - x 0 / 2 x 0 / 2 ∫ - t 0 / 2 t 0 / 2 s 1 ( t + v , x + w ) s 2 ( t + n T 2 + v , x + mL + w ) * dvdw ]]>
= R A ( t , x ; t - n T 2 , u x - mL ) e - j ω 0 ( τ - n T 2 ) , ( n = - N min , · · · , - 2 , - 1,0,1,2 , · · · , N max ) , ( m = - M min , · · · , - 2 , - 1,0,1,2 , · · · M max ) - - - ( 15 ) ]]>
在此，T代表超声波的周期，L代表采样间隔(束线间距)，RA(t，x； τ，ux)代表包络信号的自相关函数，t0代表相关窗口的轴向长度，x0代 表相关窗口的水平方向长度。另一方面，v代表用于积分的时间(τ)轴 向的变量值，w代表用于积分的束线方向的变量值，以及星号“*”代 表复共轭算子。然后，利用如下式(16)所示的二维复相关函数，来 定义二维包络相关系数C(t，x；n，m)。
C ( t , x ; n , m ) = | R 12 ( t , x ; n , m ) | | R 11 ( t , x ; 0,0 ) | · | R 22 ( t , x ; n , m ) | - - - ( 16 ) ]]>
注意，R11(t，x；0，0)代表S1(t，x)的自相关函数，以及R22(t，x；n，m) 代表S2(t+nT/2，x+mL)的自相关函数。如同组合自相关方法一样，包络 相关系数用于解决混叠问题。也就是，在每一个测量点(t，x)处，对于 所有变量号n和m，都获得由C(t，x；n，m)和φ(t，x；n，m)构成的组合{C(t，x； n，m)，φ(t，x；n，m)}。关于该方法，假定在足够的范围内，即在足以进行 包络相关的范围内，确定变量号对(n，m)。在该情况下，与显示出最大 包络相关系数的(n，m)＝(k，l)相对应的相位φ(t，x；k，l)匹配无混叠的相 位。其原因在于，如果将显示出最大包络相关系数的变量号对(n，m) 表示为(k，l)，则s1(t，k)与s2(t+kT/2，x+lL)之间的时移|τ-kT/2|小于T/2。 即，|φ(t，x；k，l)|＝ω0|t-kT/2|小于π。也就是，关于该方法，利用无混叠的 相位φ(t，x；k，l)进行计算，由此在每一个测量点(t，x)处获得正确的时移 τ、正确的轴向位移uy以及正确的水平方向位移ux，如下式(17)所 示。
τ = - φ ( t , x ; k , l ) ω 0 + k T 2 ]]>
u y = cτ 2 ]]>
ux＝lL …(17)
注意，c代表组织内的声速(假定为1500m/s的常量，1500m/s是 软组织内的正常声速)。关于该方法，对所有测量点进行上述计算，由 此获得轴向位移分布uy(x，y)和水平方向位移分布ux(x，y)。
此外，关于该方法，对上述位移分布的每一个进行空间微分，由 此获得轴向应变分布εy(x，y)和水平方向应变分布εx(x，y)，如下式(18) 所示。
ϵ y ( x , y ) = ∂ u y ( x , y ) ∂ y - - - ( 18 ) ]]>
ϵ x ( x , y ) = ∂ u x ( x , y ) ∂ y ]]>
如上所述，关于该方法，根据在向组织施加压力和没有施加压力 的情况下的RF信号，来估计轴向和水平方向的位移(应变)分布。 注意，如同可以从上式ux＝lL看出，水平方向的位移检测精度受水平 方向的采样间隔(束线间距)的限制，以及因此，该方法具有水平方 向精度较低的缺点。然而，该方法具有实时观察的优点，由此提高了 实际性能。
上述的扩展组合自相关方法具有以下功能：对于每一次计算，利 用在预定范围内应用的二维相关窗口，来分析在向组织施加压力的情 况下、组织与超声波探头之间的水平方向的相对位移，由此处理组织 的水平方向位移。然而，具有这种功能的该二维扩展组合自相关方法 不能估计，在向组织施加压力的情况下、包含与轴向和水平方向都正 交的方向(与二维超声扫描平面(切片方向)正交的方向)上的位移 的应变分布。为了解决上述问题，很容易地将以上的二维扩展组合自 相关方法扩展为利用应用于三维范围的三维窗口的三维扩展组合自相 关方法，由此使系统能够更稳定。
图9和图10是描述三维组合自相关方法的基本算法的流程图。注 意，图10是对图9所示处理的一部分进行详细描述的流程图。
在步骤S11中，扫描线寄存器1存储“1”，扫描线寄存器1用作 用于执行对于第一扫描线至第L扫描线相同的处理的计数器。扫描线 寄存器1中存储的值在步骤S23的确定处理中被检查。
在步骤S12，对于每一个循环处理，使厚度方向(帧方向)上的 变量在-U至U之间递增。注意，在步骤S18的确定处理中检查计数。
在步骤S13，对于每一个循环处理，使扫描方向(扫描线方向) 上的变量在-V至V之间递增。注意，在步骤S17的确定处理中检查 计数。
在步骤S14，对于每一个循环处理，使距离方向(轴向)上的变 量在0和M之间递增。注意，在步骤S16的确定处理中检查计数。
在步骤S15，利用组合自相关方法，计算距离方向(轴向)的包 络信号的相关系数C(l，t；u，v，n)。注意，常规的组合自相关方法不能展 示出足以进行三维计算的计算速度，以及因此，利用高速组合自相关 方法来计算相关系数C(l，t；u，v，n)。后面将描述高速自相关方法。
在步骤S16，对在上述步骤S14中确定的变量进行处理。即，确 定距离方向寄存器中存储的变量n是否已达到预定最大值M。如果确 定变量n已达到M，则流程前进到步骤S17。否则，流程返回到步骤 S14，并且使距离方向寄存器中存储的变量n递增。
在步骤S17，对在上述步骤S13中确定的变量进行处理。即，确 定扫描方向寄存器中存储的变量v是否已达到预定最大值V。如果确 定变量v已达到V，则流程前进到步骤S18。否则，流程返回到步骤 S13，并且使扫描方向寄存器中存储的变量v递增。
在步骤S18，对在上述步骤S12中确定的变量进行处理。即，确 定厚度方向寄存器中存储的变量u是否已达到预定最大值U。如果确 定变量u已达到U，则流程前进到步骤S19。否则，流程返回到步骤 S12，并且使扫描方向寄存器中存储的变量u递增。
在步骤S19，系统在u＝(-U，…，0…，U)、v＝(-V，…，0…，V)和n＝(0， 1，…，N)的范围内，搜索展示出在步骤S12至S18中计算的最大相关 系数C(l，t；u，v，n)的变量组合(u，v，n)。注意，虽然本实施例假定只有正 压力施加于组织，而采用了这种距离方向范围n＝(0，1，…，N)，但是不 用说，在其中正压力和负压力都被施加于组织上的配置中，应该采用 距离方向范围n＝(-M，…，0，…，N)。
在步骤S20，计算在步骤S19中获得的相关系数C(l，t；u，v，n)的相 位差φ(l，t；u0，v0，n0)。
在步骤S21，相位差n0π+φ(l，t；u0，v0，n0)被计算，作为最终阶段的 相位差。
在步骤S22，利用点(u0，v0，n0)附近的点(u，v，n)处的相关系数C(l，t； u，v，n)，来计算扫描方向的位移v＝v0+Δv和厚度方向的位移u＝u0+Δu。
在步骤S23，检查在上述步骤S11中扫描线寄存器1中存储的变 量。即，确定变量1是否已达到L。如果确定1已达到L，则流程前进 到步骤S24。否则，流程返回到步骤S11，并且使扫描线寄存器1中存 储的变量递增。
在步骤S24，通过对在向组织施加压力的情况下估计的位移分布 进行空间微分，来计算应变分布。
图10是对图9中所示的步骤S15中的上述高速组合自相关方法 进行详细描述的流程图。
在步骤S31，对没有施压和施压情况下的RF信号进行正交检波， 由此获得I信号和Q信号，如下所述。
x(t)→Ix，Qx(其中X(t)＝Ix+jQx)
y(t)→Iy，Qy(其中Y(t)＝Iy+jQy)
在步骤S32，根据下式计算相关Rxy、Rxx和Ryy。
Rxy＝∫X(t+v)·Y*(t+v)dv
Rxx＝∫X(t+v)·X*(t+v)dv
Ryy＝∫Y(t+v)·Y*(t+v)dv
在步骤S33，利用下式，根据这样获得的相关Rxy、Rxx和Ryy 来计算相关系数C0。
C 0 = | Rxy | / Rxx · Ryy ]]>
在步骤S34，将变量号n设置为1。
在步骤S35，计算Yn(t)＝Y(t-nT)ejωnT。
在步骤S36，根据下式计算Rxyn和Rynyn。
Rxyn＝∫X(t+v)·Yn*(t+v)dv
＝∫X(t+v)·Y*(t-nT+v)ejωnTdv
Rynyn＝∫Yn(t+v)·Yn*(t+v)dv
＝∫Y(t-nT+v)·Y*(t-nT+v)dv
在步骤S37，利用下式，根据Rxyn和Rynyn计算相关系数Cn。
C n = | Rxy n | / Rxx · Ry n y n ]]>
在步骤S38，使变量号n递增。
在步骤S39，确定变量号n是否已达到预定最大值M。如果确定 n已达到M，则处理结束。否则，流程返回到步骤S35，并重复相同 处理。
利用图10中流程图所示的方法，在步骤S35中从Y计算Yn，以 便计算Rxyn和Rynyn。这允许简单的电路配置。以下将描述从Y计算 Yn的方法。
首先，用下式来表示没有施压情况下的回波信号x(t)。
x(t)＝u(t)cos(ωt+θ)
另一方面，用下式表示在轴向施加压力的情况下的回波信号y(t)。
y(t)＝x(t+τ)＝u(t+τ)cos(ω(t+τ)+θ)。
则，用下式来表示x、y和yn的正交检波信号。
x(t)→
Ix＝0.5u(t)cosθ
Qx＝-0.5u(t)sinθ
(X(t)＝Ix+jQx＝0.5u(t)e-jθ)
y(t)→
Iy＝0.5u(t+τ)cos(ωτ+θ)
Qy＝-0.5u(t+τ)sin(ωτ+θ)
(Y(t)＝Iy+jQy＝0.5u(t+τ)e-j(ωτ+θ))
yn(t)＝y(t-nT)
＝u(t+τ-nT)cos(ω(t+τ-nT)+θ)
＝u(t+τ-nT)cos(ωt+ω(τ-nT)+θ)
在此，T代表二分之一周期。则，获得如下式所示的Iyn和Qyn。
Iyn＝0.5u(t+τ-nT)cos(ω(τ-nT)+θ)
Qyn＝-0.5u(t+τ-nT)sin(ω(τ-nT)+θ)
(Yn＝Iyn+Qyn＝0.5u(t+τ-nT)e-j(ω(τ-nT)+θ))
因此，根据上式获得了以下关系式。
Yn(t)＝Iyn+Qyn
＝0.5u(t+τ-nT)e-j(ω(τ-nT)+θ)
＝Y(t-nT)ejωnT
如同可以从上式看出，从Y(t)＝Iy+jQy计算Yn(t)。
因此，从X和Y计算Rxyn和Rynyn，如下式所示。
Rxyn＝4∫X(t+v)·Yn*(t+v)dv
＝4∫X(t+v)·Y*(t-nT+v)ejωnTdv
|Rxyn|＝Run
＝∫u(t+v)u(t+τ-nT+v)dv
＝4∫|X(t+v)·Yn*(t+v)|dv
＝4∫|X(t+v)·Y*(t-nT+v)ejωnT|dv
＝4∫|X(t+v)·Y*(t-nT+v)|dv
Rynyn＝∫u(t+τ-nT+v)u(t+τ-nT+v)dv
＝4∫|Yn(t+v)·Yn*(t+v)|dv
＝4∫Y(t-nT+v)·Y*(t-nT+v)dv
在此，星号“*”代表复共轭算子。
图11显示了用于执行三维组合自相关方法的基本算法的电路配 置的框图。关于用于执行如图7所示的组合自相关方法的常规电路配 置，设置了大量的正交检波电路1320至132N，并且这些正交检波电 路1320至132N中的处理需要大量的处理时间，造成了难以进行高速 计算，导致了难以进行实时图像显示。因此，本实施例采用用于执行 如图11所示的上述基本算法的电路配置，由此允许操作用于执行组合 自相关方法的高速电路。
未受压正交检波电路(QD)131接收没有施压情况下的回波信号 x(t)，对收到的回波信号进行正交检波，并将正交检波的信号Ix(t)和 Qx(t)输出到第一相关计算电路133和第一相关系数计算电路1350至 135N。受压正交检波电路(QD)132接收施压情况下的回波信号y(t)， 对收到的回波信号进行正交检波，并将正交检波的信号Y(t)＝Iy+ jQy(Iy(t)，Qy(t))输出到第一相关系数计算电路1350、第二相关计算电 路1340、第一延迟电路134和第二延迟电路135。第一延迟电路134 和第二延迟电路135使正交检波信号Y(t)延迟超声波的周期T，并将 延迟的正交检波信号Y(t-T)输出到第一相关计算电路1351、第三延迟 电路136和第四延迟电路137。第三延迟电路136和第四延迟电路137 的每一个都使正交检波信号Y(t-T)延迟超声波的周期T，并将延迟的 正交检波信号Y(t-2T)输出到随后的第一相关系数计算电路和延迟电 路(未显示)。以同样的方式，利用N个延迟电路的每一个，使正交 检波信号连续地延迟周期T，并将延迟的信号输出到相应的第一相关 系数计算电路。
第一相关计算电路133根据信号Ix和Qx，计算相关值Rxx，并 将相关值Rxx输出到第二相关系数计算电路1380至138N。第二相关 计算电路1340从受压正交检波电路132接收正交检波信号Iy(t)和 Qy(t)，根据信号Iy和Qy计算相关值Ryy，并将相关值Ryy输出到第 二相关系数计算电路1380。第一相关系数计算电路1350从未受压正 交检波电路131接收正交检波信号Ix(t)和Qx(t)，并从受压正交检波电 路132接收正交检波信号Iy(t)和Qy(t)，计算复基带信号SR和SI，并 将复基带信号SR和SI输出到第三相关计算电路1360和相位差计算电 路1370。第三相关计算电路1360从第一相关系数计算电路1350接收 复基带信号SR和SI，根据复基带信号SR和SI计算相关值|Rxy|，并将 相关值|Rxy|输出到第二相关系数计算电路1380。相位差计算电路1370 从第一相关系数计算电路1350接收复基带信号SR和SI，并根据复基 带信号SR和SI计算相位差φ0(t)。第二相关系数计算电路1380从第一 相关计算电路133接收相关值Rxx，从第三相关计算电路1360接收相 关值|Rxy|，以及从第二相关计算电路1340接收相关值Ryy，根据这些 相关值计算相关系数C0(t)，并输出相关系数C0(t)。
第二相关计算电路1341从第一延迟电路134和第二延迟电路135 接收延迟的正交检波信号Iy(t-T)和Qy(t-T)，根据信号Iy(t-T)和Qy(t-T) 计算相关值Ry1y1，并将相关值Ry1y1输出到第二相关系数计算电路 1381。第一相关系数计算电路1351从未受压正交检波电路131接收正 交检波信号Ix(t)和Qx(t)，并从第一延迟电路134和第二延迟电路135 接收延迟的正交检波信号Iy(t-T)和Qy(t-T)，获得复基带信号SR1和SI1， 并将复基带信号SR1和SI1输出到第三相关计算电路1361和相位差计 算电路1371。第三相关计算电路1361从第一相关系数计算电路1351 接收复基带信号SR1和SI1，根据复基带信号SR1和SI1获得相关值|Rxy1|， 并将相关值|Rxy1|输出到第二相关系数计算电路1381。相位差计算电 路1371从第一相关系数计算电路1351接收复基带信号SR1和SI1，并 根据复基带信号SR1和SI1计算相位差φ1(t)。第二相关系数计算电路 1381从第一相关计算电路133接收相关值Rxx，从第三相关计算电路 1361接收相关值|Rxy1|，以及从第二相关计算电路1341接收相关值 Ry1y1，根据这些相关值计算相关系数C1(t)，并输出相关系数C1(t)。
以同样的方式，利用从第三延迟电路135和第四延迟电路136向 下排列的第二相关计算电路1342至134N、第一相关系数计算电路 1352至135N、第三相关计算电路1362至136N、相位差计算电路1372 至137N、以及第二相关系数计算电路1382至138N，对连续延迟的正 交检波信号Iy(t-2T)、…、Iy(t-NT)和Qy(t-2T)、…、Qy(t-NT)的每一 个执行如上所述的相同处理，由此输出相关系数C2(t)至CN(t)和相位 差φ2(t)至φN(t)。
接下来，将描述利用三维有限元方法的弹性模量分布重构方法。 为了简化用于重构弹性模量分布的反演问题，关于本实施例，用模型 来表示组织。这允许用户利用有限元方法来执行本实施例中提出的弹 性模量分布重构方法。特别是，关于本实施例，通过基于下述假定的 模型来表示组织。
首先，假定组织展示出各向同性的弹性。另一方面，对处于外部 静压力下的组织，估计组织应变的分布。注意，对组织施加静压力， 以促使组织的轻微压缩，由此允许用户利用向组织施加压力情况下的 RF信号与没有施加压力情况下的RF信号之间的关系来进行计算。因 此，在该情况下，应力与应变之间的关系显示出线性关系。也就是， 可以假定用弹性模型表示组织，来进行近似的计算。此外，在该情况 下，假定组织是各向同性的，以及因此，假定组织展示出各向同性弹 性。
此外，假定组织展示出接近不可压缩性。一般来说，众所周知， 组织展示出接近不可压缩性(Poisson比(泊松比)v＝0.5)的可压缩 性。关于本实施例，利用0.49的组织内恒定泊松比进行计算。注意， 关于本实施例，不是在组织展示出完全不可压缩性的假定之下进行计 算。其原因在于，如果在组织展示出完全不可压缩性的假定之下进行 计算，即利用0.5的组织内恒定泊松比进行计算，则弹性方程将变成 特殊情况，导致了不能利用本实施例中提出的有限元方法进行计算的 问题。此外，关于本实施例，假定泊松比在组织内是恒定的，以及因 此，只应为弹性模量分布估计Young模量(杨氏模量)，由此简化了 反演问题。注意，一般来说，关于本实施例，利用0.49的组织内恒定 泊松比进行计算。
接下来，用三维有限元模型来表示组织。关于根据本实施例的弹 性模量分布重构方法，利用有限元方法进行计算，以及因此，将组织 划分为有限数量的矩形平行六面体单元。注意，假定每一个单元在其 中展示出恒定的弹性模量、恒定应力以及恒定应变。一般来说，为了 看出解决反演问题的方法，看出与反演问题相对应的正演问题是重要 的。关于本实施例，反演问题是，根据应变分布和边界条件来估计弹 性模量分布。因此，与上述反演问题相对应的正演问题是，根据弹性 模量分布和边界条件来计算应变分布。关于本实施例，利用作为一种 周知的数值分析方法的有限元方法(FEM：有限元方法)，来解决上 述问题。
关于有限元方法，通过由有限数量的单元构成的模型，来表示充 当要估计的连续体的组织，并且利用数值分析来求解代表每个单元的 特性的联立线性方程。注意，后面将描述利用有限元方法的特殊计算。 简而言之，关于有限元方法，根据用作输入值的组织弹性模量分布和 边界条件，来获得用作输出值的应变(位移)分布和应力分布。
关于本实施例，在组织展示出各向同性弹性的假定之下进行近似 计算，以及因此，弹性方程(平衡方程、应变与位移之间的关系、应 力与应变之间的关系)在如下式所示的组织内条件下保持。
用下式(19)来表示平衡方程。
∂ σ ij ∂ x j + f i = 0 , ( i , j = 1,2,3 ) - - - ( 19 ) ]]>
在此，充当i和j的参考数字1、2和3分别代表x、y和z。另一 方面，用下式(20)来表示应变与位移之间的关系。
ϵ ij = 1 2 ( ∂ u i ∂ x j + ∂ u j ∂ x i ) - - - ( 20 ) ]]>
用下式(21)来表示应力与应变之间的关系(通用化的Hooke法 则)
σ ij = E 1 + v ( ϵ ij + v 1 - 2 v δ ij ϵ nn ) - - - ( 21 ) ]]>
利用张量表示上述弹性方程。因此，实际上有三个平衡方程，六 个应变位移关系式，以及六个应力应变关系式。注意，坐标(x1，x2，x3) 代表(x，y，z)。其它参考字符代表如下的特性。
E：杨氏模量(也被称为“弹性模量”)
v：泊松比
εij：应变张量
(εnn＝ε11+ε22+ε33：体积应变)
sij：应力张量
δij：Kronecker delta(克罗内克尔增量)
ui：位移向量
fi：体积力向量
(注意，重力被认为是可以忽略的，以及因此，假定在该情况下 fi为0)。
然后，为εij求解应力与应变之间的关系式。结果，将应变与应力 之间的关系式转变为如下式(22)所示。
ϵ ij = 1 + v E ( σ ij - v 1 + v δ ij σ nn ) - - - ( 22 ) ]]>
在此，snn＝s11+s22+s33。于是，将i＝j＝2代入上式(22)，并 对上式求解杨氏模量E，由此获得下式(23)。
E = σ 22 - v ( σ 11 + σ 33 ) ϵ 22 ]]>
= σ y - v ( σ x + σ z ) ϵ y - - - ( 23 ) ]]>
如同可以从上式(23)看出，可以根据轴向(在本实施例中为超 声波束方向)的应变分量和所有方向的应力分量来计算杨氏模量，即 弹性模量。然而，难以利用当前技术对用于上式的应力分布进行直接 测量。因此，关于本实施例，交替地估计和更新应力分布和弹性模量 分布，以致于所估计的弹性模量分布变得接近于实际分布。如下来执 行用于重构弹性模量分布的特定过程。
首先，假定将均匀分布用作用于估计弹性模量分布的初始分布 {E^0}。第二，利用三维有限元方法，来计算由于初始弹性模量分布{E^0} 所引起的应力分布{S^0}。特别是，将应变位移关系式和应力应变关系 式代入以上的平衡方程中，由此形成新的平衡方程，如下式(24)所 示。将该新的平衡方程应用于组织模型内的每个单元。
∂ ∂ x j [ E 2 ( 1 + v ) ( ∂ u i ∂ x j + ∂ u j ∂ x i ) + Ev ( 1 + v ) ( 1 - 2 v ) δ ij ∂ u n ∂ x n ] + f i = 0 - - - ( 24 ) ]]>
在此，
∂ u n ∂ x n = ∂ u 1 ∂ x 1 + ∂ u 2 ∂ x 2 + ∂ u 3 ∂ x 3 - - - ( 25 ) ]]>
然后，在以下边界条件下、利用Gaussian(高斯)消去法对联立 方程求解位移，由此获得与弹性模量分布{E^0}对应的位移分布{u^0}。
ui|y＝底部＝0
σi|y＝顶部＝pi … (26)
σn|x，z＝侧面＝0
在上式中，pi代表体表上的外部压力向量，以及sn代表与侧面正 交的应力分量。第一式表示底部固定的边界条件，第二式表示体表的 应力分布匹配外部压力分布的边界条件，以及第三式表示每一侧面具 有自由端的边界条件。接下来，将位移分布{u^0}代入应变位移关系式 中，由此获得与弹性模量分布{E^0}对应的应变分布{ε^0}。然后，将应 变分布{ε^0}代入应力应变关系式中，由此获得与弹性模量分布{E^0} 对应的应力分布{S^0}。
第三，利用下式(27)，根据利用三维有限元方法计算的应力分布 以及利用扩展组合自相关方法估计的轴向(y方向)应变分布{εy}，来 更新弹性模量分布{E^k}。
E ^ k + 1 = σ ^ y k - v ( σ ^ x k + σ ^ z k ) ϵ y - - - ( 27 ) ]]>
注意，通过在i＝j＝2的情况下、对上述应力应变关系式求解杨 氏模量E，来获得上式。在上式中，k代表迭代号。
第四，根据这样更新的弹性模量分布和上述边界条件，连续进行 三维有限元分析，由此更新应力分布。
然后，连续执行第三和第四处理，由此弹性模量分布接近实际的 弹性模量分布。注意，如果弹性模量分布满足下式(28)，则确定所估 计的弹性模量分布达到了收敛，并结束估计处理。
1 N Σ l N | E ^ l k + 1 - E ^ l k | < Γ - - - ( 28 ) ]]>
在此，l代表单元号，N代表单元的总数，以及Γ代表阈值。
已经描述了本实施例中提出的利用三维有限元模型的弹性模量分 布重构方法。关于该方法，根据三维平衡方程来估计弹性模量分布。 因此，关于该方法，在比常规方法更接近于实际组织的假定下进行计 算，由此使得能够更精确地估计弹性模量分布。此外，关于该方法， 仅仅根据能够高精度地估计的轴向应变分布，来重构弹性模量分布， 由此使得能够可靠地重构弹性模量分布。注意，关于该方法，为组织 估计三维弹性模量分布，以及因此，需要采用二维阵列超声波探头， 或者需要对切片方向的一维阵列超声波探头进行机械扫描，以便对要 分析的组织进行三维扫描。
以下，将描述根据本实施例的、基于扩展组合自相关方法和三维 有限元模型的弹性模量分布重构方法的优点，利用模拟证实了该方法 的优点。图12显示了模拟过程的示意图。
首先，创建具有用于估计测试的弹性模量分布的组织模型。注意， 该组织模型中包含散射点，这些散射点反射超声波回波信号。第二， 向组织模型施加外部压力，以便在计算机模拟中造成组织模型的压缩。 然后，利用有限元方法等，计算由于压缩而引起的每个散射点的位移。 第三，根据在向组织模型施加压力和没有施加压力的情况下的散射点 位移分布，来模拟在施压和没有施压情况下的RF信号。第四，将扩 展组合自相关方法应用于在施压和没有施压情况下的模拟RF信号， 由此估计组织应变分布。第五，通过利用三维有限元模型的弹性模量 分布重构方法，根据利用扩展组合自相关方法估计的应变分布以及为 模拟组织模型的压缩而确定的边界条件(外部压力分布等等)，来估计 组织弹性模量分布。
虽然在该模拟中，不同的弹性模量分布用于组织模型，但是假定 在此使用的各种弹性模量分布都展示出各向同性的弹性。注意，该模 拟中使用的弹性模量一般匹配与根据本实施例的组织弹性模量分布测 量系统的主要用途相对应的乳房组织的弹性模量。另一方面，每一个 组织模型中都包含散射点，这些散射点用于模拟在向组织施加压力和 没有施加压力的情况下的反射RF信号。特别是，每一个组织模型都 包含平均密度为500点/cm3的散射点。此外，利用均值为1.0、标准差 为0.3的正态伪随机数，来确定在没有施压情况下的散射点位置。然 后，通过根据来自有限元分析的结果计算在没有施压情况下的每个散 射点的位移，来获得在施压情况下的散射点位置。在此，虽然关于实 际组织内的散射点的信息是未知的，但是这样确定散射点的每个参数， 以致于根据模拟的RF信号而生成的B模式图像类似于实际组织的B 模式图像。
关于该方法，通过使施压和没有施压情况下的散射点函数与超声 系统的点扩展函数进行卷积，来计算在向组织施压和没有施压的情况 下的模拟RF信号，如下式(29)所示。
i1(x，y，z)＝∫∫∫h(x-x′，y-y′，z-z′)t1(x′，y′，z′)dx′dy′dz′ … (29)
i2(x，y，z)＝∫∫∫h(x-x′，y-y′，z-z′)t2(x′，y′，z′)dx′dy′dz′
在此，i1(x，y，z)代表没有向组织施压的情况下的RF信号，i2(x，y， z)代表向组织施压的情况下的RF信号，h(x，y，z)代表超声系统的点扩 散函数(脉冲响应)，t1(x，y，z)代表没有向组织施压的情况下的散射点 函数，以及t2(x，y，z)代表向组织施压的情况下的散射点函数。注意， 散射点函数在每一个散射点处代表这样预先确定的散射系数，并且在 除散射点以外的位置代表散射系数0。另一方面，根据由于组织模型 的应变而引起的每个散射点的位移，从没有向组织施压情况下的散射 点函数t1(x，y，z)，来计算向组织施压情况下的散射点函数t2(x，y，z)。 注意，通过对利用有限元分析计算的单元节点处的位移向量进行线性 插值，来计算由于组织的压缩而引起的每个散射点的位移。
关于根据本实施例的模拟，假定采用无超声波阻尼的非聚焦超声 系统。也就是，假定点扩散函数h(x，y，z)对于每一点都恒定。此外， 假定可以将点扩散函数分为如下式(30)所示的、相对于各个方向的 函数。
h(x，y，z)＝hx(x)hy(y)hz(z) … (30)
在此，hy(y)代表超声波束方向的点扩散函数。另一方面，hx(x)和 hz(z)的每一个代表与超声波束方向正交的方向上的点扩散函数。特别 是，假定x方向代表超声层析图像的平面内方向(水平方向)，以及z 方向代表垂直于超声层析图像的方向(切片方向)。注意，根据利用超 声设备、对从线目标(水中延伸的、直径为0.13mm的线(wire))反 射的回波信号进行实际测量而获得的回波信号分布，来建立每个方向 的点扩散函数。图13描述了用于利用5.0MHz的超声波中心频率进行 模拟的每个点扩散函数。图13(A)显示了轴向的点扩散函数hy(y)。通 过使Gaussian(高斯)函数与正弦波相乘，来获得点扩散函数hy(y)， 并且hy(y)用作从线目标反射的实际反射回波分布的近似分布。另一方 面，图13(B)显示了水平方向的点扩散函数hx(x)，图13(C)显示了切片 方向的点扩散函数hz(z)。注意利用Gaussian(高斯)函数来建立上述 点扩散函数的每一个，以便以同样的方式用作从线目标反射的实际反 射回波分布的近似分布。每个函数的参数随中心频率而变，将在后面、 在说明每个模拟时对此进行描述。
接下来，将描述作为根据本实施例的位移(应变)分布估计方法 的扩展组合自相关方法的优点，通过模拟证实了该方法的优点。首先， 将描述估计水平方向的组织位移的优点，这是优于组合自相关方法的 优点。
图14显示了组织模型。组织模型由60mm×60mm(二维尺寸) 的外部尺寸和均匀弹性模量分布构成。然后，模拟组织的压缩，以造 成轴向3％的均匀应变。假定在该模拟中，采用简单的一维弹性模型， 作为用来仅仅评价扩展组合自相关方法的优点的组织模型。此外，利 用水平方向0.0mm至1.4mm的位移模拟组织的轴向压缩，以评价处 理水平方向位移(关于超声波探头的、组织的水平方向相对位移)的 优点。此外，假定模拟水平方向上的组织的简单平行位移，该位移对 应于超声波探头相对于组织滑动的情形。
然后，在具有组织应变和没有组织应变的情况下，为组织模型模 拟RF信号。用于超声系统模拟的参数包括：5.0MHz的中心频率， 0.5mm的脉冲宽度；1.0mm的超声波束宽；0.5mm的扫描线间距；以 及30MHz的采样频率。然后，利用本实施例中提出的扩展组合自相 关方法，根据在施压和没有施压情况下的模拟RF信号，来估计应变 分布。另外，准备利用组合自相关方法估计的应变分布和利用空间相 关方法估计的应变分布，以进行比较。注意，将相同尺寸的相关窗口 应用于具有相同搜索范围的各方法，由此允许用户在它们之间对精度 等进行简单的比较。特别是，扩展组合自相关方法和空间相关方法采 用应用于5.6mm(轴向)×7.5mm(水平方向)的二维搜索范围中的、 且尺寸为1.6mm(轴向)×2.5mm(水平方向)的二维窗口。另一方 面，关于用于分析一维位移的组合自相关方法，将具有相同的1.6mm 轴向长度的一维相关窗口应用于相同的5.6mm轴向范围中。
图15至图17显示了利用各方法的应变分布估计结果。图15显示 了利用各方法估计的水平方向的应变误差。在此，附图标记“◇”代表 利用组合自相关方法的误差，“□”代表利用扩展组合自相关方法的误 差，以及“△”代表利用空间自相关方法的误差。注意，用下式(31) 来表示估计的应变的误差e。
e = Σ i N ( ϵ ^ i - ϵ i ) 2 Σ i N ϵ i 2 - - - ( 31 ) ]]>
在此，ε^i代表所估计的应变，εi代表实际应变(理想值)，i代表 单元号，以及N代表单元的总数。另一方面，图16显示了，利用各 方法(组合自相关方法、扩展组合自相关方法和空间相关方法)估计 的、包含0.0mm水平方向位移的组织的应变分布。图17显示了，利 用各方法(组合自相关方法、扩展组合自相关方法和空间相关方法) 估计的、包含0.4mm水平方向位移的组织的应变分布。注意，图16 和图17显示了相对于沿轴向的每一深度的估计应变的平均值和标准 偏差。
如同可以从模拟结果看出，关于常规的组合自相关方法(CA法)， 水平方向上超过超声波束尺寸(在该情况下为二分之一束宽，即 0.5mm)的组织相对位移导致了应变估计质量的迅速恶化。另一方面， 可以看出，扩展组合自相关方法能够稳定地估计应变，与水平方向位 移的大小无关。此外，可以看出，虽然空间相关方法也能够稳定地估 计应变，而与水平方向位移的大小无关，但是估计结果显示出较差的 精度，即与扩展组合自相关方法相比，显示出两倍或更大的误差。此 外，对上述方法的处理时间进行比较，虽然扩展组合自相关方法需要 的处理时间是组合自相关方法的3.6倍，但是扩展组合自相关方法需 要的处理时间仅是空间相关方法的1/(7.7)倍，如下表所示。如同可以 从上述结果看出，扩展组合自相关方法在某种程度上能够实现实时计 算。
方法 处理时间 规一化的处理时间 组合自相关方法 26秒 1/(3.6) 扩展组合自相关方法 1分34秒 1.0 空间相关方法 12分5秒 7.7
接下来，将描述在倾斜方向上施压的情况下的估计结果。注意， 和利用具有二维均匀结构的简单组织模型的上述模拟不同，关于该估 计，利用具有如同实际组织一样的三维结构的组织模型进行估计。注 意，在理想的测量中，应该利用超声波探头在超声波束方向(轴向) 上向组织施加压力。现在，将描述受倾斜方向上的组织压缩的影响的 估计结果。
图18显示了用于估计倾斜方向上的组织压缩的影响的组织模型。 如图18(A)所示，该组织模型具有外部尺寸为60mm×60mm×60mm 的三维结构，并且包含直径为15mm、长度为60mm、高硬度的圆柱 内含物。假定包围内含物的材料具有10Kpa的弹性模量(杨氏模量)， 并且内含物的弹性模量是上述材料的3倍，即具有30Kpa的弹性模量。 注意，根据乳房组织的弹性模量和乳腺癌的弹性模量，来确定上述弹 性模量，乳腺癌的诊断是本发明的主要目的。然后，在以下两种情形 下，模拟组织模型的压缩。一种情形对应于，由于沿轴向、从顶面向 组织模型施加200Pa的均匀外部压力而造成的2％的轴向组织模型压 缩，如图18(B)所示。另一种情形对应于，由于沿倾斜方向、从顶面 向组织模型施加均匀外部压力(轴向200Pa和水平方向30Pa)而造成 的倾斜方向的组织模型压缩，如图18(C)所示。
然后，对于上述两种情形，模拟在施压和没有施压情况下的RF 信号，接着利用扩展组合自相关方法估计应变分布。注意，也利用组 合自相关方法和空间相关方法估计应变分布，以进行比较。在此，对 于各方法，将具有相同尺寸的相关窗口应用于相同搜索范围内的计算， 由此允许简单的比较。注意，如同上述的模拟一样，相关窗口具有相 同尺寸。此外，将用于模拟RF信号的其它参数确定为如同上述模拟 的情况一样的值。也就是，其它参数包括：5.0MHz的中心频率，0.5mm 的脉冲宽度；1.0mm的水平方向超声波束宽；2.0mm的切片方向超声 波束宽；0.5mm的扫描线间距；以及30MHz的采样频率。
图19和图20显示了上述模拟的模拟结果。图19显示了，其中组 织模型在轴向被压缩的简单情形下的应变分布的估计结果。图20显示 了，其中组织模型在倾斜方向被压缩的情形下的应变分布的估计结果。 注意，对于各情形，把通过三维有限元分析获得的轴向应变分布用作 利用理想方法估计的应变分布，该轴向应变分布充当实际的应变分布。 注意，图19和图20是沿组织模型的中线取的组织模型横截面视图， 图19和图20显示了估计结果。在图20中，利用理想方法估计的应变 分布展示出左右不对称性。其原因在于，压力是在倾斜方向施加的。 特别是，在该情况下，压力是在附图中的右下方向上施加的，导致了 附图中右上部分的水平方向上的大位移。
首先，已证实了，如同组合自相关方法一样，对于包含轴向压缩 的应变分布，扩展组合自相关方法显示出一般相同的估计结果。另一 方面，在倾斜方向上向组织施加压力的情形下，也已证实了，虽然由 于水平方向的大位移，而造成不能利用组合自相关方法估计某些区域， 但是扩展组合自相关方法能够稳定地估计应变分布，而与水平方向的 位移大小无关，如同以上的上述模拟所描述的那样。另一方面，还已 证实了，虽然空间相关方法能够稳定地进行估计，而与水平方向的位 移大小无关，但是空间相关方法的估计精度比扩展组合自相关方法差 得多。除上述模拟结果以外，已证实了，扩展组合自相关方法具有应 变分布估计的优点。
如上所述，已证实了，可以利用上述扩展组合自相关方法、高速、 高精度地估计组织应变分布。现在将描述，通过对本实施例中提出的、 利用三维有限元模型的弹性模量分布重构方法进行模拟而获得的估计 结果。注意，弹性模量分布重构方法是根据应变分布来估计弹性模量 分布的方法，该方法用作在组织弹性成像系统的第二级中执行的方法。
本实施例中提出的弹性模量分布重构方法的主要功能是，根据三 维动力学平衡方程来估计弹性模量分布。现在，通过对根据本实施例 的弹性模量分布重构方法与二维重构方法进行比较，来证实本实施例 中提出的弹性模量分布重构方法的优点，该二维重构方法除了利用二 维动力学平衡方程进行计算以外，具有同样的处理。注意，假定在组 织中出现平面应变，利用二维重构方法进行计算。在该模拟中，采用 两种模型作为组织模型，该两种模型的每一种具有如同实际组织一样 的三维结构，如图21所示。也就是，图21显示了具有三维结构的两 种组织模型例子。图21(a)显示了包含充当乳房肿瘤模型的内含物的内 含物包含模型。特别是，该内含物包含模型具有100mm×100mm× 100mm的外部尺寸，并包含直径为20mm的坚硬内含物。假定该内含 物具有30kPa的弹性模量，并且包围该内含物的材料具有10kPa的弹 性模量。以如同上述模拟一样的方式，根据实际乳房组织的弹性模量 来确定上述弹性模量。另一方面，内含物和包围内含物的材料都展示 出接近不可压缩性，以及因此，都被假定为具有0.49的相同泊松比。 图21(b)显示了用于模拟诸如肌肉的层状结构组织的层状结构模型。层 状结构模型具有100mm×100mm×100mm的外部尺寸，并包含厚度 为20mm的硬层。假定硬层具有30kPa的弹性模量，并且包围该硬层 的材料具有10kPa的弹性模量。注意，层状结构模型也具有0.49的均 匀泊松比。
然后，在图21(a)所示的内含物包含模型的情况下，在计算机上模 拟，在沿轴向、从模型的顶面施加的100Pa均匀外部压力下的压缩。 另一方面，在图21(b)所示的层状结构模型的情况下，在计算机上模拟， 在沿轴向、从模型的顶面施加的150Pa均匀外部压力下的压缩。注意， 对上述两种模型模拟不同外部压力下的压缩，以致于对每种模型模拟 大约1％的相同应变。然后，对于每种组织模型，模拟在施压和没有 施压的情况下的RF信号，并利用扩展组合自相关方法估计轴向的应 变分布。随后，利用三维弹性模量分布重构方法，根据估计的轴向应 变分布以及为压缩模拟而确定的边界条件，来估计弹性模量分布。同 样，利用二维重构方法，根据相同的轴向应变分布和相同边界条件， 来估计弹性模量分布，以进行比较。在此，用于模拟RF信号的参数 包括：3.75MHz的中心频率，0.75mm的脉冲宽度；2.0mm的水平方 向超声波束宽；2.0mm的切片方向超声波束宽；以及2.0mm的扫描线 间距。另一方面，用于利用扩展组合自相关方法进行计算的参数包括： 3.2mm(轴向)×4.0mm(水平方向)的相关窗口尺寸，以及11.2mm (轴向)×14.0mm(水平方向)的搜索范围。另一方面，关于利用三 维有限元模型的弹性模量分布重构方法，各组织模型由50000矩形平 行六面体单元构成，每一个矩形平行六面体单元都具有2mm(轴向) ×2mm(水平方向)×5mm(切片方向)的尺寸。
图22至图25显示了模拟结果。图22和图23显示了内含物包含 模型的估计结果。另一方面，图24和图25显示了层状结构模型的估 计结果。注意，虽然可以利用三维重构方法来估计三维弹性模量分布， 但是图24和图25的每一幅图都显示了沿模型的中线取得的模型横截 面视图。其原因在于，利用二维重构方法只能估计二维弹性模量分布。 因此，图24和图25显示了沿模型中线取得的估计结果的横截面视图， 由此允许用户对它们进行比较。另一方面，在下表中，列出了利用各 组织模型获得的三维重构方法和二维重构方法的估计值。
包围模型核 心的区域中 的弹性模量 误差[％] 包围模型核 心的区域中 的标准偏差 [％] 模型核心中 的对比度误 差[％] 内含物包含 模型 三维重构方法 3.5 15.5 11.0 二维重构方法 30.9 17.9 35.9 层状结构模 型 三维重构方法 8.5 26.8 3.1 二维重构方法 24.9 22.1 43.5
在此使用的估计参数包括：包围模型核心的区域中的弹性模量误 差eS；包围模型核心的区域中的标准偏差SDS；以及模型核心区域中 的对比度误差eC，它们由下式来定义。
e S = 1 N S Σ i N S | E ^ i - E i | E ‾ S ]]>
SD S = 1 N S Σ i N S ( E ^ i - E ‾ S ) 2 E ‾ S ]]>
e C = | ( E ^ C - E ‾ S ) - ( E C - E S ) | E C - E S - - - ( 32 ) ]]>
注意，在上式中，S代表包围核心的区域中的和，E^代表估计的 弹性模量，E代表实际弹性模量，NS代表包围核心的区域中的单元总 数，E-S代表包围模型核心的区域中的弹性模量的平均值，E^C代表所 估计的模型核心区域中的弹性模量，EC代表模型核心区域中的实际弹 性模量，以及ES代表包围模型核心的区域中的实际弹性模量。
如同可以从上述模拟结果看出，已证实了，与基于组织内出现平 面应力的假定而形成的二维重构方法相比，本实施例中提出的、利用 三维有限元模型的弹性模量分布重构方法具有更高精度地估计弹性模 量分布的优点。虽然三维重构方法能够对弹性模量分布进行精确估计， 但是利用二维重构方法估计的弹性模量分布显示出比实际弹性模量分 布小的值。不用说，其原因在于，关于二维重构方法，在与二维计算 平面正交的方向上的计算受到限制。在该评价中，已经清楚地证实了， 利用与实际组织对应的三维计算的弹性模量分布重构方法适用于对实 际组织进行分析。
接下来，将描述采用上述扩展组合自相关方法和利用三维有限元 模型的上述弹性模量分布重构方法的超声诊断系统的特定配置。图26 显示了超声诊断系统的基本配置。超声诊断系统包括：三维超声扫描 仪281，超声诊断设备282，个人计算机283，脉冲电动机控制器284， 脉冲电动机285，以及压力计286等。三维超声扫描仪281具有向组 织投射超声脉冲信号、以及接收从组织反射的超声回波信号的功能。 注意，该系统采用利用三维有限元模型的弹性模量重构方法，以及因 此，该系统需要组织内的三维数据。因此，该超声诊断系统具有利用 三维超声扫描仪281进行三维扫描的配置。超声诊断设备282具有控 制三维超声扫描仪281、以及显示实时超声B模式图像的功能，由此 允许用户确定要测量的区域的位置。注意，可以在不进行更改的情况 下，将常规的超声诊断设备用作超声诊断设备282。该超声诊断设备 包括全数字装置，该全数字装置具有用于临时存储测量的RF信号的 帧存储器。个人计算机283接收由超声诊断设备282测量的RF信号， 执行处理(上述本实施例中提出的处理)以估计组织弹性特性，以及 显示估计结果。脉冲电动机285具有控制施加于组织上的压力的功能。 脉冲电动机285被固定在固定臂的顶端，并且三维超声扫描仪281被 安装在脉冲电动机285的活动部分上。该系统具有以下配置，其中脉 冲电动机控制器284控制脉冲电动机285，以调节超声扫描仪281的 垂直方向位置，由此允许用户调节施加于组织上的压力，该压力高精 度地促使组织发生百分之几的轻微压缩。压力计286具有测量体表上 的压力的功能；压力用作为重构弹性模量分布所需的边界条件。注意， 压力计286位于超声扫描仪281与体表之间。在此，假定通过超声扫 描仪281的组织压缩造成了体表上的均匀压力分布，则将利用压力计 286测量的压力用作体表上的压力。
图27显示了超声诊断系统中包括的超声扫描仪281的特定配置。 三维超声扫描仪281不具有其中二维地部署了超声换能器的二维阵列 配置，而是具有三维扫描配置，在该三维扫描配置中，在水中、在切 片方向上进行二维凸形扫描探头的机械扫描。
图26所示的超声诊断系统主要被设计用于乳腺癌的诊断，以及因 此，针对乳腺癌的诊断来确定超声扫描仪的特性参数。特别是，在此 使用的超声扫描仪的主要特性参数包括：7.5MHz的超声波中心频率， 30MHz的采样频率，71的扫描线数，44的帧数，30°的换能器扫描角， 以及0.5秒的三维扫描周期的持续时间。注意，通过上述换能器扫描 角，在切片方向上进行凸形探头的扫描。此外，上述的帧数表示在凸 形探头的一个扫描周期内获取的扫描图像(帧)的数量。另一方面， 通过利用水中的线(wire)进行实际测量获得的超声脉冲的特性包括： 大约0.5mm的脉冲宽度；大约1.5mm的水平方向束宽；以及大约 2.6mm的切片方向束宽。
接下来，将描述图26所示的超声诊断系统进行弹性特性测量的操 作。首先，用户移动被安装在臂上的三维扫描仪281，以便将超声扫 描仪281定位在要测量的期望部分，同时监测通过超声诊断设备282 获得的实时B模式图像。注意，在对超声扫描仪281进行定位的时候， 不进行超声扫描仪281的三维扫描(即，不进行凸形探头的机械扫描)， 而只有与扫描仪的中心面相对应的B模式图像被显示在超声诊断设备 282上。其原因在于，关于在此使用的超声诊断设备282，不能在进行 三维扫描的同时、显示实时的B模式图像。不用说，可以采样这样一 种超声诊断设备，该设备具有能够在进行三维扫描的同时、显示实时 的B模式图像的功能，由此允许用户在进行三维扫描的同时、对三维 超声扫描仪281进行定位。在把超声扫描仪281定位在要测量的期望 位置之后，用户固定臂的活动部分，以便固定超声扫描仪281。然后， 系统对没有向组织施压的情况下的三维RF信号进行测量。注意，一 旦用户按压用于三维扫描的按钮，三维扫描就自动地进行。在此，一 个三维扫描周期的持续时间仅为0.5秒。没有施压情况下的测量RF 数据被存储在超声设备内的帧存储器中。随后，一旦用户按压用于操 作脉冲电动机控制器284的按钮一次，以控制组织的压缩，被固定在 臂上的脉冲电动机285就使超声扫描仪281移动预定距离，由此通过 超声扫描仪281压缩组织。随后，脉冲电动机285停止，同时维持组 织的压缩。在该状态下，用户再次按压用于三维扫描的按钮，由此测 量在向组织施压的情况下的RF数据。注意，以如同没有施压情况下 的RF数据一样的方式，把向组织施压情况下的RF数据存储在超声设 备282中包括的帧存储器中。此外，利用被安装在超声扫描仪281的 末端的压力计，来测量施加于组织上的压力。然后，测量结束，并且 释放施加于组织上的压力，此后释放对象。
在释放对象之后，系统通过个人计算机283访问超声诊断设备282 中包括的帧存储器，并且在向组织施压和没有施压的情况下的RF数 据被存储在个人计算机283中包括的硬盘上。执行这种处理是因为， 超声设备中包括的帧存储器具有临时存储数据的功能，即该帧存储器 只有用于存储一个测量周期的数据的容量。个人计算机283执行这样 的程序，该程序利用扩展组合自相关方法和采用三维有限元模型的弹 性模量分布重构方法进行计算，以便根据向组织施压和没有施压的情 况下的RF数据，来估计应变分布和弹性模量分布。然后，个人计算 机283通过执行显示程序，在监视器上显示B模式图像、应变图像和 弹性模量图像。
根据本发明的、利用应变分布显示方法和弹性模量分布显示方法 的超声诊断系统具有以下优点：能够仅仅根据超声波束方向(轴向) 的应变分布，来重构弹性模量分布；以及能够与水平方向位移无关地， 估计位移分布。
注意，虽然已经描述了采用包络信号的配置，但是本发明不限于 上述配置，而是代表包括振幅、波高和波数的波动性之间的关系的任 何参数都可以被采用。