CN101530334B - 超声诊断系统 - Google Patents

Abstract

提供一种超声诊断系统，将由超声波探头所接收的对象的压缩前后的反射回波信号作为多个切片面的三维体积数据存储，让置所存储的压缩后的反射回波信号的包络线相对于压缩前的反射回波信号的包络线在三维方向上依次相对移动，求出压缩前后的包络线的相关系数成为最大的三维位置以及在该位置下的所述压缩前后的包络线的相位差；并根据所述三维位置和相位差，求出伴随所述压缩的所述对象的组织的三维位移以及所述对象的组织的三维应变分布。这样，可以在不受到位移量的制约的情况下推定位移及应变分布，缩短计算处理时间，并且对于三维方向上的变位及应变也能应对。

Description

超声诊断系统

本申请是申请号为“03820755.9(PCT/JP03/09731)”，申请日为2003年7月31日，发明名称为“超声诊断系统和应变分布显示方法”之申请的分案申请。 

技术领域

本发明涉及超声诊断系统和应变分布(strain distribution)显示方法，该超声诊断系统和应变分布显示方法允许用户利用超声诊断设备，对组织的硬度进行定量的测量。 

背景技术

超声诊断不仅应用于对组织结构的观察，而且应用于以下领域(超声组织表征)：对诸如声速、阻尼系数等的组织内物理量进行测量，而且根据这样测量的物理量，生成诊断图像。作为该领域的一部分，已知这样的领域，其中对组织的硬度，即弹性特性，进行测量。因为组织的弹性特性与病理状况有着紧密的关系，因此上述领域正在被重点研究。例如，已知受到以下影响的组织展示出比正常组织大的硬度：硬化性肿瘤，诸如乳腺癌、甲状腺癌等；肝硬化；等等。按照惯例，通过触摸来检测组织的硬度。然而，通过触摸的检测具有难以进行客观分析的缺点，需要熟练的外科医生，并具有以下限制：只能检测具有某种尺寸或更大尺寸、且位于体表附近的受影响组织。 

另一方面，已知这样一种方法，其中向组织施加静压力，以便使组织压缩和变形，并利用超声波来测量与施加的压力相对应的组织内应变(strain)，以便估计组织的弹性特性(J.Ophir，I.Cespedes，H.Ponnekanati，Y. Yazdi，and X.Li，“Elastography：A quantitative methodfor 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信号之间的相位差。此外，上述的高速处理方法可以应用于三维配置。此外，可以通过以超声波束的切片间距为单位、在切片方向上改变测量点，来进行上述计算。 

此外，可以利用显示装置来显示这样获得的弹性模量分布。在该情况下，弹性模量计算装置优选地通过，在对象的组织展示出各向同性的弹性和接近不可压缩性的假定下，将对象的组织划分为有限数量的矩形平行六面体单元，来建立三维有限元模型。此外，弹性模量计算装置优选地在对象的组织的每个单元都展示出均匀的弹性模量、均匀应力和均匀应变的假定下，利用弹性方程、根据关于上述应变分布的信息，来计算弹性模量分布。 

如上所述，在组织展示出各向同性的弹性的假定下，进行根据本发明的计算。其原因在于，在施加于组织的外部静压力下，应力与应变之间的关系展示出线性关系。因此，可以在组织用作弹性模型的假定下，进行近似计算。另外，组织展示出各向同性的特性，以及因此，可以在组织用作各向同性的弹性模型的假定下，进行根据本发明的计算。另一方面，关于本发明，在组织展示出接近不可压缩性的假定下，进行计算。其原因在于，如果在组织展示出完全不可压缩性的假定下，进行计算，即利用0.5的组织内恒定Poisson比(泊松比)进行计算，则弹性方程变成特殊情况，造成了不能利用有限元方法进行计算的问题。注意，关于本发明，利用均匀的组织内泊松比进行计算。在该情况下，只应估计Young模量(杨氏模量)，以估计弹性模量分布，由此简化了反演问题。注意，泊松比展示出足够的组织内均匀性，以及因此，关于本发明，优选地利用0.49的泊松比进行计算。根据本发明的弹性模量分布计算能够仅仅根据可以高精度地计算的轴向应变分布，来重构弹性模量分布，由此能够稳定地计算弹性模量分布。 

附图说明

图1描述了超声诊断设备的机理。 

图2(A)、(B)显示了通过施加静压力的组织弹性测量方法的特定例子，以及通过施加静压力的组织弹性测量方法的机理。 

图3显示了空间相关方法的机理。 

图4显示了Doppler(多普勒)方法的机理。 

图5显示了组合自相关方法的机理。 

图6显示了用于执行组合自相关方法的基本算法的电路配置的框 图。 

图7显示了根据本发明实施例的超声诊断系统的示例配置的框图。 

图8显示了三维组合自相关方法的基本算法的流程图。 

图9所示为，根据本发明的超声诊断系统中所采用的三维组合自相关方法的基本算法的流程图，以及图7所示处理的一部分的详细流程图。 

图10显示了，用于详细描述具有提高的计算速度的组合自相关方法的流程图，该组合自相关方法在图9所示的步骤S15中被执行。 

图11显示了，用于执行根据本发明的超声诊断系统中所采用的三维组合自相关方法的基本算法的电路配置的框图。 

图12显示了模拟过程的示意图。 

图13(A)～(C)显示了，在5.0MHz的超声波中心频率的情况下每一点的点扩展函数的例子。 

图14显示了组织模型的示意图。 

图15显示了，利用各估计方法估计的、由于水平方向位移造成的估计应变误差。 

图16(A)～(C)显示了，在0.0mm的水平方向位移的情况下利用各估计方法(组合自相关方法、扩展组合自相关方法和空间相关方法)估计的应变分布。 

图17(A)～(C)显示了，在0.4mm的水平方向位移的情况下利用各估计方法(组合自相关方法、扩展组合自相关方法和空间相关方法)估计的应变分布。 

图18(A)～(C)显示了，用于模拟倾斜方向的压缩的组织模型。 

图19(a)～(d)显示了，通过对轴向的组织模型的简单压缩进行模拟而获得的应变分布的估计结果。 

图20(a)～(d)显示了，通过对倾斜方向的组织模型压缩进行模拟而获得的应变分布的估计结果。 

图21(a)、(b)显示了两种组织模型例子，该两种组织模型例子的每一种都具有三维结构。 

图22(a)～(d)是显示了利用内含物包含模型估计的结果的第一图。 

图23(a)、(b)是显示了利用内含物包含模型估计的结果的第二图。 

图24(a)～(d)是显示了利用层状结构模型估计的结果的第一图。 

图25(a)、(b)是显示了利用层状结构模型估计的结果的第二图。 

图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=\frac{2L}{c}...\left(1\right)$

在此，c表示组织内的声速，在穿过软组织时，c可以被确定为大约1500米/秒的常量。因此，根据从投射超声波直到收到反射回波信号的时间t，来计算探头与反射物之间的距离L。此外，反射回波信号具有关于组织的声学特性的信息，并且因此可以根据反射回波信号，将诸如B模式层析图像的组织信息图像显示在监视器上。 

例如，已知这样一种方法，其中利用超声诊断设备来测量代表组织硬度的弹性特性。上述用于测量弹性特性的方法有一种机理，其中将机械振动作用于组织上，并且基于横波(transverse wave)以高传播速度穿过硬组织、并且以低传播速度穿过软组织这一事实，根据这样产生的横波的传播速度，来估计硬度信息。严格地讲，传播穿过组织的横波的传播速度v取决于组织ρ的密度、剪切模量μ₁、剪切粘度μ₂和振动的角频率ω，如下式(2)所表示。 

$v=\sqrt{\frac{2({{\mathrm{\μ}}_1}^2+{\mathrm{\ω}}^2{{\mathrm{\μ}}_2}^2)}{\mathrm{\ρ}({\mathrm{\μ}}_1+\sqrt{{{\mathrm{\μ}}_1}^2+{\mathrm{\ω}}^2{{\mathrm{\μ}}_2}^2})}}...\left(2\right)$

另一方面，如图2所示，已提出了这样一种方法，其中将静压力施加于组织上，并且根据在所施加的静压力下的组织应变分布，来估 计组织的弹性特性。上述方法是基于以下事实而设计的：这种静压力在硬组织内造成小应变，而在软组织内造成大应变。图2(A)显示了通过施加静压力的组织弹性测量方法的特定例子。图2(B)显示了通过施加静压力的组织弹性测量方法的机理。由图可以清楚地看出，关于上述方法，利用带有超声波探头10的常规超声诊断设备、在不施加压力的情况下，正常地测量组织11的超声回波信号(在不施加压力的情况下的RF信号)。随后，通过超声波探头10轻微地(大约百分之几)按压组织11，并测量在施加压力下的组织11的超声回波信号(在施加压力的情况下的RF信号)。然后，根据在向组织施加压力和没有施加压力的情况下测量的RF信号，来估计代表受压时的组织内的每一点的位移的位移分布。位移分布估计方法的主要例子包括利用空间相关的方法和利用Doppler(多普勒)效应的方法。 

图3显示了空间相关方法的机理。关于该方法，通过利用二维相关函数进行模板匹配，根据在向组织施加压力和没有施加压力的情况下测量的RF信号(或RF信号的包络信号)，来估计在施加压力下的组织内的应变分布。以下将描述该方法的特定处理。首先，分别以i₁(t，x)和i₂(t，x)来表示在向组织施加压力和没有施加压力的情况下测量的RF信号(或RF信号的包络信号)，则上述两个信号之间的交叉相关系数C(t，x；n，m)以下式(3)来表示。 

$C(t,x;n,m)=\frac{\underset{v=-t_0/2}{\overset{t_0/2}{\mathrm{\Σ}}}\underset{w=-x_0/2}{\overset{x_0/2}{\mathrm{\Σ}}}{i}_1(t+v,x+w)i_2(t+v+n{L}_t,x+w+m{L}_x)}{\sqrt{\underset{v=-t_0/2}{\overset{t_0/2}{\mathrm{\Σ}}}\underset{w=-x_0/2}{\overset{x_0/2}{\mathrm{\Σ}}}{i_1}^2(t,x)}\·\sqrt{\underset{v=-t_0/2}{\overset{t_0/2}{\mathrm{\Σ}}}\underset{w=-x_0/2}{\overset{x_0/2}{\mathrm{\Σ}}}{i_2}^2(t+v+n{L}_t,x+w+m{L}_x)}}...\left(3\right)$

在此，t代表超声波束方向(轴向)的坐标点，x代表与超声波束方向正交(水平方向)的坐标点，t₀代表相关窗口的轴向尺寸，x₀代表相关窗口的水平方向尺寸，L_t代表轴向的采样间隔，以及L_x代表水平方向的采样间隔。此外，n和m为整数。在该情况下，如果将显示出最大值的交叉相关系数的组合(n，m)表示为(k，1)，则在测量点(t，x)处的轴向位移u_y和水平方向位移u_x分别以下式来表示。 

u_y＝kL_t

u_x＝1L_x

注意，在水平方向的采样间隔L_x比轴向采样间隔L_t粗糙的情况下进行数据采样，由此导致了水平方向的位移检测精度比轴向位移检测精度差。对每个测量点进行上述处理，由此估计位移分布。空间相关方法的优点在于，使用户能够估计二维位移向量分量。此外，即使组织内出现大的应变(大约5％)，空间相关方法也允许用户估计位移分布。然而，不同于常规超声测量系统，上述方法引起大量的计算，造成实时计算的困难。此外，位移估计精度受采样间隔的限制，由此造成了与Doppler(多普勒)方法相比、精度较差的问题，将在后面描述Doppler(多普勒)方法。 

图4显示了Doppler(多普勒)方法的机理。关于该方法，根据在向组织施加压力和没有施加压力的情况下的RF信号，利用Doppler(多普勒)效应来估计在施加压力下的组织内的位移分布，Doppler(多普勒)效应也用于血流测量中。以下将描述该方法的特定处理。首先，利用如下式(4)所表示的模型，来表示在向组织施加压力和没有施加压力的情况下的RF信号。 

在此，i₁(t)代表在没有施加压力情况下的RF信号，i₂(t)代表在施加压力情况下的RF信号，A(t)代表包络信号，ω₀代表超声波的中心角频率，以及τ代表时移(time shift)。一旦对两个RF信号的每一个执行正交检波，就获得了基带信号，如下式所示。 

s₁(t)＝A(t)e^jθ

$s_2\left(t\right)=A(t-\mathrm{\τ})e^{j({\mathrm{\ω}}_0\mathrm{\τ}+\mathrm{\θ})}...\left(5\right)$

同样，以下式来表示上述两个信号之间的复相关函数R₁₂(t)。 

$R_{12}\left(t\right)={\∫}_{-t_0/2}^{t_0/2}{s}_1(t+v)s_2{(t+v)}^{*}\mathrm{dv}=R_A\left(t\right)e^{-j{\mathrm{\ω}}_0\mathrm{\τ}}...\left(6\right)$

在此，R_A(t)代表包络信号的自相关函数，以及t₀代表超声波束轴向的相关窗口尺寸。此外，星号“*”代表复共轭算子。因此，从相关函数R₁₂(t)的相位φ(t)，获得了由于施加压力而造成的时移τ和轴向位移 u_y，如下式(7)所示。 

$\mathrm{\τ}=-\frac{\mathrm{\φ}\left(t\right)}{{\mathrm{\ω}}_0}$

$u_y=\frac{\mathrm{c\τ}}{2}...\left(7\right)$

注意，c代表组织内的声速，并且假定在组织内c为常量。 

利用Doppler(多普勒)方法，对每一个测量点执行上述处理，并且通过与根据Doppler(多普勒)效应开发的血流测量一样的方式，估计位移分布。因此，Doppler(多普勒)方法具有实时测量的优点。此外，Doppler(多普勒)方法利用相位信息进行计算，由此与空间相关方法相比，可以更高精度地估计位移。然而，Doppler(多普勒)方法的缺点在于：大位移的测量，例如超声波中心频率的波长的四分之一或更大的位移的测量，造成了混叠(aliasing)，导致了不能估计正确位移的问题。此外，Doppler(多普勒)方法的缺点还在于，不能够象可以从上式所推定的那样，来估计除轴向以外的位移。 

为了解决上述问题，本发明人提出了“组合自相关方法(CA法)”，该方法同时具有上述两种方法的优点。图5显示了由本发明人提出的组合自相关方法的机理。关于组合自相关方法，利用RF信号的包络信号的相关进行计算，由此解决了作为Doppler方法的缺点的混叠问题。以下将描述该方法的特定处理。 

首先，如同Doppler方法的情况一样，利用如下式所表示的模型，来表示在向组织施加压力和没有施加压力的情况下的RF信号。 

在此，i₁(t)代表在没有施加压力情况下的RF信号，i₂(t)代表在施加压力情况下的RF信号，A(t)代表包络信号，ω₀代表超声波的中心角频率，以及τ代表时移。一旦对两个RF信号的每一个执行正交检波，就获得了基带信号，如下式所示。 

s₁(t)＝A(t)e^jθ

$s_2\left(t\right)=A(t-\mathrm{\τ})e^{j({\mathrm{\ω}}_0\mathrm{\τ}+\mathrm{\θ})}...\left(9\right)$

然后，利用下式来表示上述两个信号之间的复相关函数R₁₂(t；n)。 

$R_{12}(t;n)={\∫}_{-t_0/2}^{t_0/2}{s}_1(t+v)s_2{(t+n\frac{T}{2}+v)}^{*}\mathrm{dv}=R_A(t;\mathrm{\τ}-n\frac{T}{2})e^{-j{\mathrm{\ω}}_0(\mathrm{\τ}-n\frac{T}{2})}...\left(10\right)$

(n＝...，-2，-1，0，1，2...) 

在此，T代表超声波的周期，R_A(t；τ)代表包络信号的自相关函数，以及t₀代表相关窗口尺寸。此外，星号“*”代表复共轭算子。在此，n代表变量号，并且对不同的n进行每次计算。在由变量号确定的点的周围，计算每个测量点处的位移。注意，在n＝0的情况下，组合自相关函数和如式(6)所示的Doppler(多普勒)方法中的自相关函数相等。也就是，在对超声波波长的四分之一或更大的位移进行测量的情况下，n等于0的组合自相关函数造成了混叠问题。为了解决上述问题，关于组合自相关方法，将包络相关系数C(t；n)定义为如下式(11)所表示。 

$C(t;n)=\frac{\left|R_{12}(t;n)\right|}{\sqrt{\left|R_{11}(t;0)\right|\·\left|R_{22}(t;n)\right|}}...\left(11\right)$

注意，R₁₁(t；0)代表s₁(t)的自相关函数，以及R₂₂(t；n)代表s₂(t+nT/2)的自相关函数。如果将显示出最大值的包络相关系数C(t；n)的n表示为k，则利用n＝k的R₁₂(t；k)的相位φ进行计算。在该情况下，无混叠地计算位移。其原因在于，在二分之一波长的间距下进行包络相关计算。注意，二分之一波长的计算间距是用于计算位移、同时防止混叠的最大间距。从而，利用φ(t；k)来计算时移τ和轴向位移u_y，如下式所示。 

$\mathrm{\τ}=-\frac{\mathrm{\φ}(t;k)}{{\mathrm{\ω}}_0}+k\frac{T}{2}$

$u_y=\frac{\mathrm{c\τ}}{2}...\left(12\right)$

注意，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，并且延迟的回波信号y₁＝y(t-T)被输入到第二受压正交检波(QD)电路1321。同样，第二延迟电路135使已经被第一延迟电路134延迟的回波信号y₁＝y(t-T)延迟超声波的周期T，并且延迟的回波信号y₂＝y(t-2T)被输入到下一个第二受压正交检波(QD)电路1322(未显示)。注意，该电路具有N个延迟电路，回波信号被连续延迟，并且被延迟了整数倍周期T的回波信号以同样的方式被输入到相应的受压正交检波电路。 

第一相关计算电路133根据信号Ix和Qx计算相关值Rxx，并将相关值Rxx输出到第二相关系数计算电路1380至138N。第二相关计算电路1340从受压正交检波电路1320接收正交检波信号Iy(t)和Qy(t)，根据信号I_y和Qy计算相关值Ryy，并将相关值Ryy输出到第二相关计算电路1380。第一相关系数计算电路1350从未受压正交检波电路131接收正交检波信号Ix(t)和Qx(t)，并从第一受压正交检波电路1320接收正交检波信号Iy(t)和Qy(t)，据此计算复基带信号S_R和S_I，并将基带信号S_R和S_I输出到第三相关计算电路1360和相位差计算电路1370。第三相关计算电路1360从第一相关系数计算电路1350接收复基带信号S_R和S_I，据此计算相关值|Rxy|，并将计算的相关值|Rxy|输出到第二相关系数计算电路1380。相位差计算电路1370从第一相关系数计算电路1350接收复基带信号S_R和S_I，并根据复基带信号S_R和S_I计算相位差φ₀(t)。第二相关系数计算电路1380从第一相关 计算电路133接收相关值Rxx，从第三相关计算电路1360接收相关值|Rxy|，以及从第二相关计算电路1340接收相关值Ryy，根据上述相关值计算相关系数C₀(t)，并输出计算的相关系数C₀(t)。 

第二受压正交检波电路(QD)1321接收被第一延迟电路134延迟的回波信号y₁＝y(t-T)，并将正交检波信号Y₁(t)＝Iy₁+jQy₁(Iy₁(t)，Qy₁(t))输出到第一相关系数计算电路1341和第二相关计算电路1351。第二相关计算电路1341从第二受压正交检波电路(QD)1321接收正交检波信号Iy₁(t)和Qy₁(t)，根据信号Iy₁(t)和Qy₁(t)计算相关值Ry₁y₁，并将相关值Ry₁y₁输出到第二相关计算电路1381。第一相关系数计算电路1351从未受压正交检波电路131接收正交检波信号Ix(t)和Qx(t)，以及从第二受压正交检波电路(QD)1321接收正交检波信号Iy₁(t)和Qy₁(t)，据此计算复基带信号S_R1和S_I1，并将复基带信号S_R1和S_I1输出到第三相关计算电路1351和相位差计算电路1371。第三相关计算电路1361从第一相关系数计算电路1351接收复基带信号S_R1和S_I1，根据复基带信号S_R1和S_I1计算相关值|Rxy₁|，并将计算的相关值|Rxy₁|输出到第二相关系数计算电路1381。相位差计算电路1371从第一相关系数计算电路1351接收复基带信号S_R1和S_I1，并根据复基带信号S_R1和S_I1计算相位差φ₁(t)。第二相关系数计算电路1381从第一相关计算电路133接收相关值Rxx，从第三相关计算电路1361接收相关值|Rxy₁|，以及从第二相关计算电路1341接收相关值Ry₁y₁，根据上述相关值计算相关系数C₁(t)，并输出计算的相关系数C₁(t)。 

同样地，从自第一延迟电路135向下的相应延迟电路接收信号的第二受压正交检波电路(QD)1322至132N的每一个，第二相关计算电路1342至134N的每一个，第一相关系数计算电路1352至135N的每一个，第三相关电路1362至136N的每一个，相位差计算电路1372至137N的每一个，以及第二相关系数计算电路1382至138N的每一个，都执行如同如上所述的第一级和第二级电路元件一样的处理，由此输出相关系数C₂(t)至C_N(t)以及相位值φ₂(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收到的水平方向最大相关相对应的位置，计算水平方向位移u_x，并将该位移输出到水平方向应变计算单元96。另一方面，轴向位移计算单元95根据与从二维复相关计算单元93收到的轴向最大相关和相位相对应的位置，计算轴向位移u_y，并将该位移输出到轴向应变计算单元97。水平方向应变计算单元96对从水平方向计算单元94收到的水平方向位移u_x进行空间微分，以便计算水平方向的应变分布ε_x，并将该应变分布输出到量化单元98。另一方面，轴向应变计算单元97对从轴向计算单元95收到的轴向位移u_y进行空间微分，以便计算轴向应变分布ε_y，并将该应变分布输出到量化单元98。量化单元98对水平方向应变分布ε_x和轴向应变分布ε_y进行量化，以便对这些应变分布进行灰阶显示(或彩色显示)，并在显示单元99上显示信息。显示单元99显示这样量化的每一个应变分布。 

接下来，将描述图7所示的超声诊断系统中采用的扩展组合自相关方法的操作。首先，考虑这样一种情况，其中组织被轻微地压缩(即百分之几或更少)。在该情况下，从局部观点看，假定压力造成了组织内的每一点的平行位移。也就是，通过下式所示的模型，来表示向组织施加压力和没有施加压力的情况下的RF信号。 

在此，i₁(t，x)代表没有施加压力情况下的RF信号，i₂(t，x)代表施加压力情况下的RF信号，A(t，x)代表包络信号，ω₀代表超声波的中心角频率，τ代表时移，该时移用作代表轴向位移的时间参数，u_x代表水平方向位移，以及θ代表初始相位。注意，和Doppler(多普勒)方法和组合自相关方法不同，关于该方法，施加压力和没有施加压力情况下的RF信号由考虑到水平方向位移的模型来表示。要在最后阶段获得的参数是轴向位移u_y＝cτ/2(即时移τ)和水平方向位移u_x。注意， c代表组织内的声速，并且假定在组织内c是常量。 

然后，正交检波器92对向组织施加压力和没有施加压力的情况下的RF信号进行正交检波。也就是，将具有与超声波中心频率相同的频率的正弦波和余弦波应用于RF信号，随后进一步对RF信号进行低通滤波，由此获得复基带信号s₁和s₂，如下式(14)所示。 

s₁(t，x)＝A(t，x)e^jθ

$s_2(t,x)=A(t-\mathrm{\τ},x-u_x)e^{j({\mathrm{\ω}}_0\mathrm{\τ}+\mathrm{\θ})}...\left(14\right)$

然后，将s₁(t，x)和s₂(t+nT/2，x+mL)之间的二维复相关函数R₁₂(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\frac{T}{2}+v,x+\mathrm{mL}+w)}^{*}\mathrm{dvdw}$

$=R_A(t,x;t-n\frac{T}{2},u_x-\mathrm{mL})e^{-j{\mathrm{\ω}}_0(\mathrm{\τ}-n\frac{T}{2})}...\left(15\right)$

(n＝-N_min，...，-2，-1，0，1，2，...，N_max) 

(m＝-M_min，...，-2，-1，0，1，2，...，M_max) 

在此，T代表超声波的周期，L代表采样间隔(束线间距)，R_A(t，x；τ，u_x)代表包络信号的自相关函数，t₀代表相关窗口的轴向长度，x₀代表相关窗口的水平方向长度。另一方面，v代表用于积分的时间(τ)轴向的变量值，w代表用于积分的束线方向的变量值，以及星号“*”代表复共轭算子。然后，利用如下式(16)所示的二维复相关函数，来定义二维包络相关系数C(t，x；n，m)。 

$C(t,x;n,m)=\frac{\left|R_{12}(t,x;n,m)\right|}{\sqrt{\left|R_{11}(t,x;\mathrm{0,0})\right|\·\left|R_{22}(t,x;n,m)\right|}}...\left(16\right)$

注意，R₁₁(t，x；0，0)代表S₁(t，x)的自相关函数，以及R₂₂(t，x；n，m)代表S₂(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)，则s₁(t，k)与s₂(t+kT/2，x+lL)之间的时移|τ-kT/2|小于T/2。即，|φ(t，x；k，l)|＝ω₀|t-kT/2|小于π。也就是，关于该方法，利用无混叠的相位φ(t，x；k，l)进行计算，由此在每一个测量点(t，x)处获得正确的时移τ、正确的轴向位移u_y以及正确的水平方向位移u_x，如下式(17)所示。 

$\mathrm{\τ}=-\frac{\mathrm{\φ}(t,x;k,l)}{{\mathrm{\ω}}_0}+k\frac{T}{2}$

$u_y=\frac{\mathrm{c\τ}}{2}$

u_x＝lL                          ...(17) 

注意，c代表组织内的声速(假定为1500m/s的常量，1500m/s是软组织内的正常声速)。关于该方法，对所有测量点进行上述计算，由此获得轴向位移分布u_y(x，y)和水平方向位移分布u_x(x，y)。 

此外，关于该方法，对上述位移分布的每一个进行空间微分，由此获得轴向应变分布ε_y(x，y)和水平方向应变分布ε_x(x，y)，如下式(18)所示。 

${\mathrm{\ϵ}}_y(x,y)=\frac{{\∂u}_y(x,y)}{\∂y}...\left(18\right)$

${\mathrm{\ϵ}}_x(x,y)=\frac{{\∂u}_x(x,y)}{\∂y}$

如上所述，关于该方法，根据在向组织施加压力和没有施加压力的情况下的RF信号，来估计轴向和水平方向的位移(应变)分布。注意，如同可以从上式u_x＝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。否则，流程返回到步骤S 13，并且使扫描方向寄存器中存储的变量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；u₀，v₀，n₀)。 

在步骤S21，相位差n₀π+φ(l，t；u₀，v₀，n₀)被计算，作为最终阶段的相位差。 

在步骤S22，利用点(u₀，v₀，n₀)附近的点(u，v，n)处的相关系数C(l，t；u，v，n)，来计算扫描方向的位移v＝v₀+Δv和厚度方向的位移u＝u₀+Δu。 

在步骤S23，检查在上述步骤S11中扫描线寄存器1中存储的变量。即，确定变量1是否已达到L。如果确定1已达到L，则流程前进到步骤S24。否则，流程返回到步骤S11，并且使扫描线寄存器1中存储的变量递增。 

在步骤S24，通过对在向组织施加压力的情况下估计的位移分布进行空间微分，来计算应变分布。 

图10是对图9中所示的步骤S 15中的上述高速组合自相关方法进行详细描述的流程图。 

在步骤S31，对没有施压和施压情况下的RF信号进行正交检波，由此获得I信号和Q信号，如下所述。 

x(t)→Ix，Qx(其中X(t)＝Ix+jQx) 

y(t)→I y，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来计算相关系数C₀。 

$C_0=\left|\mathrm{Rxy}\right|/\sqrt{\mathrm{Rxx}}\·\sqrt{\mathrm{Ryy}}$

在步骤S34，将变量号n设置为1。 

在步骤S35，计算Y_n(t)＝Y(t-nT)e^jωnT。 

在步骤S36，根据下式计算Rxy_n和Ry_ny_n。 

Rxy_n＝∫X(t+v)·Y_n*(t+v)dv 

＝∫X(t+v)·Y*(t-nT+v)e^jωnTdv 

Ry_ny_n＝＝∫Y_n(t+v)·Y_n*(t+v)dv 

＝∫Y(t-nT+v)·Y*(t-nT+v)dv 

在步骤S37，利用下式，根据Rxy_n和Ry_ny_n计算相关系数C_n。 

$C_n=\left|{\mathrm{Rxy}}_n\right|/\sqrt{\mathrm{Rxx}}\·\sqrt{{\mathrm{Ry}}_n{y}_n}$

在步骤S38，使变量号n递增。 

在步骤S39，确定变量号n是否已达到预定最大值M。如果确定n已达到M，则处理结束。否则，流程返回到步骤S35，并重复相同处理。 

利用图10中流程图所示的方法，在步骤S35中从Y计算Y_n，以便计算Rxy_n和Ry_ny_n。这允许简单的电路配置。以下将描述从Y计算Y_n的方法。 

首先，用下式来表示没有施压情况下的回波信号x(t)。 

x(t)＝u(t)cos(ωt+θ) 

另一方面，用下式表示在轴向施加压力的情况下的回波信号y(t)。 

y(t)＝x(t+τ)＝u(t+τ)cos(ω(t+τ)+θ)。 

则，用下式来表示x、y和y_n的正交检波信号。 

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(ωτ+θ)) 

y_n(t)＝y(t-nT) 

＝u(t+τ-nT)cos(ω(t+τ-nT)+θ) 

＝u(t+τ-nT)cos(ωt+ω(τ-nT)+θ) 

在此，T代表二分之一周期。则，获得如下式所示的Iy_n和Qy_n。 

Iy_n＝0.5u(t+τ-nT)cos(ω(τ-nT)+θ) 

Qy_n＝-0.5u(t+τ-nT)sin(ω(τ-nT)+θ) 

(Y_n＝Iy_n+Qy_n＝0.5u(t+τ-nT)e^{-j(ω(τ-nT)+θ)}) 

因此，根据上式获得了以下关系式。 

Y_n(t)＝Iy_n+Qy_n

＝0.5u(t+τ-nT)e^{-j(ω(τ-nT)+θ)}

＝Y(t-nT)e^jωnT

如同可以从上式看出，从Y(t)＝Iy+jQy计算Y_n(t)。 

因此，从X和Y计算Rxy_n和Ry_ny_n，如下式所示。 

Rxy_n＝4∫X(t+v)·Y_n*(t+v)dv 

＝4∫X(t+v)·Y*(t-nT+v)e^jωnTdv 

|Rxy_n|＝Ru_n

＝∫u(t+v)u(t+τ-nT+v)dv 

＝4∫|X(t+v)·Y_n*(t+v)|dv 

＝4∫|X(t+v)·Y*(t-nT+v)e^jωnT|dv 

＝4∫|X(t+v)·Y*(t-nT+v)|dv 

Ry_ny_n＝∫u(t+τ-nT+v)u(t+τ-nT+v)dv 

＝4∫|Y_n(t+v)·Y_n*(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)，计算复基带信号S_R和S_I，并将复基带信号S_R和S_I输出到第三相关计算电路1360和相位差计算电路1370。第三相关计算电路1360从第一相关系数计算电路1350接收复基带信号S_R和S_I，根据复基带信号S_R和S_I计算相关值|Rxy|，并将相关值|Rxy|输出到第二相关系数计算电路1380。相位差计算电路1370 从第一相关系数计算电路1350接收复基带信号S_R和S_I，并根据复基带信号S_R和S_I计算相位差φ₀(t)。第二相关系数计算电路1380从第一相关计算电路133接收相关值Rxx，从第三相关计算电路1360接收相关值|Rxy|，以及从第二相关计算电路1340接收相关值Ryy，根据这些相关值计算相关系数C₀(t)，并输出相关系数C₀(t)。 

第二相关计算电路1341从第一延迟电路134和第二延迟电路135接收延迟的正交检波信号Iy(t-T)和Qy(t-T)，根据信号Iy(t-T)和Qy(t-T)计算相关值Ry₁y₁，并将相关值Ry₁y₁输出到第二相关系数计算电路1381。第一相关系数计算电路1351从未受压正交检波电路131接收正交检波信号Ix(t)和Qx(t)，并从第一延迟电路134和第二延迟电路135接收延迟的正交检波信号Iy(t-T)和Qy(t-T)，获得复基带信号S_R1和S_I1，并将复基带信号S_R1和S_I1输出到第三相关计算电路1361和相位差计算电路1371。第三相关计算电路1361从第一相关系数计算电路1351接收复基带信号S_R1和S_I1，根据复基带信号S_R1和S_I1获得相关值|Rxy₁|，并将相关值|Rxy₁|输出到第二相关系数计算电路1381。相位差计算电路1371从第一相关系数计算电路1351接收复基带信号S_R1和S_I1，并根据复基带信号S_R1和S_I1计算相位差φ₁(t)。第二相关系数计算电路1381从第一相关计算电路133接收相关值Rxx，从第三相关计算电路1361接收相关值|Rxy₁|，以及从第二相关计算电路1341接收相关值Ry₁y₁，根据这些相关值计算相关系数C₁(t)，并输出相关系数C₁(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)的每一个执行如上所述的相同处理，由此输出相关系数C₂(t)至C_N(t)和相位差φ₂(t)至φ_N(t)。 

接下来，将描述利用三维有限元方法的弹性模量分布重构方法。为了简化用于重构弹性模量分布的反演问题，关于本实施例，用模型来表示组织。这允许用户利用有限元方法来执行本实施例中提出的弹 性模量分布重构方法。特别是，关于本实施例，通过基于下述假定的模型来表示组织。 

首先，假定组织展示出各向同性的弹性。另一方面，对处于外部静压力下的组织，估计组织应变的分布。注意，对组织施加静压力，以促使组织的轻微压缩，由此允许用户利用向组织施加压力情况下的RF信号与没有施加压力情况下的RF信号之间的关系来进行计算。因此，在该情况下，应力与应变之间的关系显示出线性关系。也就是，可以假定用弹性模型表示组织，来进行近似的计算。此外，在该情况下，假定组织是各向同性的，以及因此，假定组织展示出各向同性弹性。 

此外，假定组织展示出接近不可压缩性。一般来说，众所周知，组织展示出接近不可压缩性(Poisson比(泊松比)v＝0.5)的可压缩性。关于本实施例，利用0.49的组织内恒定泊松比进行计算。注意，关于本实施例，不是在组织展示出完全不可压缩性的假定之下进行计算。其原因在于，如果在组织展示出完全不可压缩性的假定之下进行计算，即利用0.5的组织内恒定泊松比进行计算，则弹性方程将变成特殊情况，导致了不能利用本实施例中提出的有限元方法进行计算的问题。此外，关于本实施例，假定泊松比在组织内是恒定的，以及因此，只应为弹性模量分布估计Young模量(杨氏模量)，由此简化了反演问题。注意，一般来说，关于本实施例，利用0.49的组织内恒定泊松比进行计算。 

接下来，用三维有限元模型来表示组织。关于根据本实施例的弹性模量分布重构方法，利用有限元方法进行计算，以及因此，将组织划分为有限数量的矩形平行六面体单元。注意，假定每一个单元在其中展示出恒定的弹性模量、恒定应力以及恒定应变。一般来说，为了看出解决反演问题的方法，看出与反演问题相对应的正演问题是重要的。关于本实施例，反演问题是，根据应变分布和边界条件来估计弹性模量分布。因此，与上述反演问题相对应的正演问题是，根据弹性模量分布和边界条件来计算应变分布。关于本实施例，利用作为一种周知的数值分析方法的有限元方法(FEM：有限元方法)，来解决上 述问题。 

关于有限元方法，通过由有限数量的单元构成的模型，来表示充当要估计的连续体的组织，并且利用数值分析来求解代表每个单元的特性的联立线性方程。注意，后面将描述利用有限元方法的特殊计算。简而言之，关于有限元方法，根据用作输入值的组织弹性模量分布和边界条件，来获得用作输出值的应变(位移)分布和应力分布。 

关于本实施例，在组织展示出各向同性弹性的假定之下进行近似计算，以及因此，弹性方程(平衡方程、应变与位移之间的关系、应力与应变之间的关系)在如下式所示的组织内条件下保持。 

用下式(19)来表示平衡方程。 

$\frac{{\∂\mathrm{\σ}}_{\mathrm{ij}}}{\∂x_j}+f_i=0$ (i，j＝1，2，3)...(19) 

在此，充当i和j的参考数字1、2和3分别代表x、y和z。另一方面，用下式(20)来表示应变与位移之间的关系。 

${\mathrm{\ϵ}}_{\mathrm{ij}}=\frac{1}{2}(\frac{{\∂u}_i}{{\∂x}_j}+\frac{{\∂u}_j}{{\∂x}_i})...\left(20\right)$

用下式(21)来表示应力与应变之间的关系(通用化的Hooke法则) 

${\mathrm{\σ}}_{\mathrm{ij}}=\frac{E}{1+v}({\mathrm{\ϵ}}_{\mathrm{ij}}+\frac{v}{1-2v}{\mathrm{\δ}}_{\mathrm{ij}}{\mathrm{\ϵ}}_{\mathrm{nn}})...\left(21\right)$

利用张量表示上述弹性方程。因此，实际上有三个平衡方程，六个应变位移关系式，以及六个应力应变关系式。注意，坐标(x₁，x₂，x₃)代表(x，y，z)。其它参考字符代表如下的特性。 

E：杨氏模量(也被称为“弹性模量”) 

v：泊松比 

ε_ij：应变张量 

(ε_nn＝ε₁₁+ε₂₂+ε₃₃：体积应变) 

s_ij：应力张量 

δ_ij：Kronecker delta(克罗内克尔增量) 

u_i：位移向量 

f_i：体积力向量 

(注意，重力被认为是可以忽略的，以及因此，假定在该情况下f_i为0)。 

然后，为ε_ij求解应力与应变之间的关系式。结果，将应变与应力之间的关系式转变为如下式(22)所示。 

${\mathrm{\ϵ}}_{\mathrm{ij}}=\frac{1+v}{E}({\mathrm{\σ}}_{\mathrm{ij}}-\frac{v}{1+v}{\mathrm{\δ}}_{\mathrm{ij}}{\mathrm{\σ}}_{\mathrm{nn}})...\left(22\right)$

在此，s_nn＝s₁₁+s₂₂+s₃₃。于是，将i＝j＝2代入上式(22)，并对上式求解杨氏模量E，由此获得下式(23)。 

$E=\frac{{\mathrm{\σ}}_{22}-v({\mathrm{\σ}}_{11}+{\mathrm{\σ}}_{33})}{{\mathrm{\ϵ}}_{22}}...\left(23\right)$

$=\frac{{\mathrm{\σ}}_y-v({\mathrm{\σ}}_x+{\mathrm{\σ}}_z)}{{\mathrm{\ϵ}}_y}$

如同可以从上式(23)看出，可以根据轴向(在本实施例中为超声波束方向)的应变分量和所有方向的应力分量来计算杨氏模量，即弹性模量。然而，难以利用当前技术对用于上式的应力分布进行直接测量。因此，关于本实施例，交替地估计和更新应力分布和弹性模量分布，以致于所估计的弹性模量分布变得接近于实际分布。如下来执行用于重构弹性模量分布的特定过程。 

首先，假定将均匀分布用作用于估计弹性模量分布的初始分布{E^⁰}。第二，利用三维有限元方法，来计算由于初始弹性模量分布{E^⁰}所引起的应力分布{S^⁰}。特别是，将应变位移关系式和应力应变关系式代入以上的平衡方程中，由此形成新的平衡方程，如下式(24)所示。将该新的平衡方程应用于组织模型内的每个单元。 

$\frac{\∂}{{\∂x}_j}[\frac{E}{2(1+v)}(\frac{\∂u_i}{\∂x_j}+\frac{{\∂u}_j}{\∂x_i})+\frac{\mathrm{Ev}}{(1+v)(1-2v)}{\mathrm{\δ}}_{\mathrm{ij}}\frac{{\∂u}_n}{\∂x_n}]+f_i=0...\left(24\right)$

在此， 

$\frac{{\∂u}_n}{\∂x_n}=\frac{{\∂u}_1}{{\∂x}_1}+\frac{{\∂u}_2}{{\∂x}_2}+\frac{{\∂u}_3}{{\∂x}_3}...\left(25\right)$

然后，在以下边界条件下、利用Gaussian(高斯)消去法对联立 方程求解位移，由此获得与弹性模量分布{E^⁰}对应的位移分布{u^⁰}。 

u_i|_y＝底部＝0 

σ_i|_y＝顶部＝p_i                ...(26) 

σ_n|_{x，z＝侧面}＝0 

在上式中，p_i代表体表上的外部压力向量，以及s_n代表与侧面正交的应力分量。第一式表示底部固定的边界条件，第二式表示体表的应力分布匹配外部压力分布的边界条件，以及第三式表示每一侧面具有自由端的边界条件。接下来，将位移分布{u^⁰}代入应变位移关系式中，由此获得与弹性模量分布{E^⁰}对应的应变分布{ε^⁰}。然后，将应变分布{ε^⁰}代入应力应变关系式中，由此获得与弹性模量分布{E^⁰}对应的应力分布{S^⁰}。 

第三，利用下式(27)，根据利用三维有限元方法计算的应力分布以及利用扩展组合自相关方法估计的轴向(y方向)应变分布{ε_y}，来更新弹性模量分布{E^^k}。 

${\hat{E}}^{k+1}=\frac{{{\widehat{\mathrm{\σ}}}_y}^k-v({{\widehat{\mathrm{\σ}}}_x}^k+{{\widehat{\mathrm{\σ}}}_z}^k)}{{\mathrm{\ϵ}}_y}...\left(27\right)$

注意，通过在i＝j＝2的情况下、对上述应力应变关系式求解杨氏模量E，来获得上式。在上式中，k代表迭代号。 

第四，根据这样更新的弹性模量分布和上述边界条件，连续进行三维有限元分析，由此更新应力分布。 

然后，连续执行第三和第四处理，由此弹性模量分布接近实际的弹性模量分布。注意，如果弹性模量分布满足下式(28)，则确定所估计的弹性模量分布达到了收敛，并结束估计处理。 

$\frac{1}{N}\underset{l}{\overset{N}{\mathrm{\&Sigma;}}}|{{\hat{E}}_l}^{k+1}-{{\hat{E}}_l}^k|<\mathrm{\&Gamma;}...\left(28\right)$

在此，1代表单元号，N代表单元的总数，以及Γ代表阈值。 

已经描述了本实施例中提出的利用三维有限元模型的弹性模量分布重构方法。关于该方法，根据三维平衡方程来估计弹性模量分布。因此，关于该方法，在比常规方法更接近于实际组织的假定下进行计算，由此使得能够更精确地估计弹性模量分布。此外，关于该方法，仅仅根据能够高精度地估计的轴向应变分布，来重构弹性模量分布，由此使得能够可靠地重构弹性模量分布。注意，关于该方法，为组织估计三维弹性模量分布，以及因此，需要采用二维阵列超声波探头，或者需要对切片方向的一维阵列超声波探头进行机械扫描，以便对要分析的组织进行三维扫描。

以下，将描述根据本实施例的、基于扩展组合自相关方法和三维有限元模型的弹性模量分布重构方法的优点，利用模拟证实了该方法的优点。图12显示了模拟过程的示意图。 

首先，创建具有用于估计测试的弹性模量分布的组织模型。注意，该组织模型中包含散射点，这些散射点反射超声波回波信号。第二，向组织模型施加外部压力，以便在计算机模拟中造成组织模型的压缩。然后，利用有限元方法等，计算由于压缩而引起的每个散射点的位移。第三，根据在向组织模型施加压力和没有施加压力的情况下的散射点位移分布，来模拟在施压和没有施压情况下的RF信号。第四，将扩展组合自相关方法应用于在施压和没有施压情况下的模拟RF信号，由此估计组织应变分布。第五，通过利用三维有限元模型的弹性模量分布重构方法，根据利用扩展组合自相关方法估计的应变分布以及为模拟组织模型的压缩而确定的边界条件(外部压力分布等等)，来估计组织弹性模量分布。 

虽然在该模拟中，不同的弹性模量分布用于组织模型，但是假定在此使用的各种弹性模量分布都展示出各向同性的弹性。注意，该模拟中使用的弹性模量一般匹配与根据本实施例的组织弹性模量分布测量系统的主要用途相对应的乳房组织的弹性模量。另一方面，每一个组织模型中都包含散射点，这些散射点用于模拟在向组织施加压力和没有施加压力的情况下的反射RF信号。特别是，每一个组织模型都包含平均密度为500点/cm³的散射点。此外，利用均值为1.0、标准差为0.3的正态伪随机数，来确定在没有施压情况下的散射点位置。然后，通过根据来自有限元分析的结果计算在没有施压情况下的每个散射点的位移，来获得在施压情况下的散射点位置。在此，虽然关于实际组织内的散射点的信息是未知的，但是这样确定散射点的每个参数， 以致于根据模拟的RF信号而生成的B模式图像类似于实际组织的B模式图像。 

关于该方法，通过使施压和没有施压情况下的散射点函数与超声系统的点扩展函数进行卷积，来计算在向组织施压和没有施压的情况下的模拟RF信号，如下式(29)所示。 

i₁(x，y，z)＝∫∫∫h(x-x′，y-y′，z-z′)t₁(x′，y′，z′)dx′dy′dx′ 

                                                                       ...(29) 

i₂(x，y，z)＝∫∫∫h(x-x′，y-y′，z-z′)t₂(x′，y′，z′)dx′dy′dz′ 

在此，i₁(x，y，z)代表没有向组织施压的情况下的RF信号，i₂(x，y，z)代表向组织施压的情况下的RF信号，h(x，y，z)代表超声系统的点扩散函数(脉冲响应)，t₁(x，y，z)代表没有向组织施压的情况下的散射点函数，以及t₂(x，y，z)代表向组织施压的情况下的散射点函数。注意，散射点函数在每一个散射点处代表这样预先确定的散射系数，并且在除散射点以外的位置代表散射系数0。另一方面，根据由于组织模型的应变而引起的每个散射点的位移，从没有向组织施压情况下的散射点函数t₁(x，y，z)，来计算向组织施压情况下的散射点函数t₂(x，y，z)。注意，通过对利用有限元分析计算的单元节点处的位移向量进行线性插值，来计算由于组织的压缩而引起的每个散射点的位移。 

关于根据本实施例的模拟，假定采用无超声波阻尼的非聚焦超声系统。也就是，假定点扩散函数h(x，y，z)对于每一点都恒定。此外，假定可以将点扩散函数分为如下式(30)所示的、相对于各个方向的函数。 

h(x，y，z)＝h_x(x)h_y(y)h_z(z)               ...(30) 

在此，h_y(y)代表超声波束方向的点扩散函数。另一方面，h_x(x)和h_z(z)的每一个代表与超声波束方向正交的方向上的点扩散函数。特别是，假定x方向代表超声层析图像的平面内方向(水平方向)，以及z方向代表垂直于超声层析图像的方向(切片方向)。注意，根据利用超声设备、对从线目标(水中延伸的、直径为0.13mm的线(wire))反射的回波信号进行实际测量而获得的回波信号分布，来建立每个方向的点扩散函数。图13描述了用于利用5.0MHz的超声波中心频率进行模拟的每个点扩散函数。图13(A)显示了轴向的点扩散函数h_y(y)。通 过使Gaussian(高斯)函数与正弦波相乘，来获得点扩散函数h_y(y)，并且h_y(y)用作从线目标反射的实际反射回波分布的近似分布。另一方面，图13(B)显示了水平方向的点扩散函数h_x(x)，图13(C)显示了切片方向的点扩散函数h_z(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=\sqrt{\frac{\underset{i}{\overset{N}{\mathrm{\&Sigma;}}}{({\widehat{\mathrm{\&epsiv;}}}_i-{\mathrm{\&epsiv;}}_i)}^2}{\underset{i}{\overset{N}{\mathrm{\&Sigma;}}}{{\mathrm{\&epsiv;}}_i}^2}}...\left(31\right)$

在此，ε^_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显示了沿模型中线取得的估计结果的横截面视图，由此允许用户对它们进行比较。另一方面，在下表中，列出了利用各组织模型获得的三维重构方法和二维重构方法的估计值。 

在此使用的估计参数包括：包围模型核心的区域中的弹性模量误差e_S；包围模型核心的区域中的标准偏差SD_S；以及模型核心区域中的对比度误差e_C，它们由下式来定义。 

$e_S=\frac{\frac{1}{N_S}\underset{i}{\overset{N_S}{\mathrm{\&Sigma;}}}|{\hat{E}}_i-E_i|}{{\stackrel{\&OverBar;}{E}}_S}$

${\mathrm{SD}}_S=\frac{\sqrt{\frac{1}{N_S}\underset{i}{\overset{N_S}{\mathrm{\&Sigma;}}}{({\hat{E}}_i-{\stackrel{\&OverBar;}{E}}_S)}^2}}{{\stackrel{\&OverBar;}{E}}_S}$

$e_C=\frac{|({\hat{E}}_C-{\stackrel{\&OverBar;}{E}}_S)-(E_C-E_S)|}{E_C-E_S}...\left(32\right)$

注意，在上式中，S代表包围核心的区域中的和，E^代表估计的弹性模量，E代表实际弹性模量，N_S代表包围核心的区域中的单元总数，E^- _S代表包围模型核心的区域中的弹性模量的平均值，E^_C代表所估计的模型核心区域中的弹性模量，E_C代表模型核心区域中的实际弹性模量，以及E_S代表包围模型核心的区域中的实际弹性模量。 

如同可以从上述模拟结果看出，已证实了，与基于组织内出现平面应力的假定而形成的二维重构方法相比，本实施例中提出的、利用三维有限元模型的弹性模量分布重构方法具有更高精度地估计弹性模量分布的优点。虽然三维重构方法能够对弹性模量分布进行精确估计，但是利用二维重构方法估计的弹性模量分布显示出比实际弹性模量分布小的值。不用说，其原因在于，关于二维重构方法，在与二维计算平面正交的方向上的计算受到限制。在该评价中，已经清楚地证实了，利用与实际组织对应的三维计算的弹性模量分布重构方法适用于对实际组织进行分析。 

接下来，将描述采用上述扩展组合自相关方法和利用三维有限元模型的上述弹性模量分布重构方法的超声诊断系统的特定配置。图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模式图像、应变图像和弹性模量图像。 

根据本发明的、利用应变分布显示方法和弹性模量分布显示方法的超声诊断系统具有以下优点：能够仅仅根据超声波束方向(轴向)的应变分布，来重构弹性模量分布；以及能够与水平方向位移无关地，估计位移分布。 

注意，虽然已经描述了采用包络信号的配置，但是本发明不限于上述配置，而是代表包括振幅、波高和波数的波动性之间的关系的任何参数都可以被采用。 

Claims (8)

1.一种超声诊断系统，包括：

超声波探头，与对象之间进行超声信号的收发；

存储装置，存储由所述超声波探头所接收的所述对象的压缩前后的反射回波信号作为多个切片面的三维体积数据；

相关计算装置，让由所述存储装置所存储的压缩后的所述反射回波信号的包络线相对于压缩前的所述反射回波信号的包络线在三维方向上依次相对移动，求出压缩前后的包络线的相关系数成为最大的三维位置以及在该位置下的所述压缩前后的包络线的相位差；

位移应变计算装置，根据由所述相关计算装置所求出的所述三维位置和所述相位差，求出伴随所述压缩的所述对象的组织的三维位移以及所述对象的组织的三维应变分布；

显示装置，显示所述三维应变分布。

2.根据权利要求1所述的超声诊断系统，其特征在于，

所述三维方向是由所述超声波探头所接收的超声波束的轴方向、该超声波束的扫描方向和与这些方向正交的切片方向。

3.根据权利要求2所述的超声诊断系统，其特征在于，

在所述超声波束的轴方向上，以所述超声信号的2分之一波长间隔进行移动。

4.根据权利要求2所述的超声诊断系统，其特征在于，

在所述超声波束的轴方向上，以所述超声信号的2分之一波长间隔进行移动；在所述超声波束的扫描方向以所述超声波束间距进行移动；在所述切片方向上以所述超声波束的切片间距进行移动。

5.根据权利要求1所述的超声诊断系统，其特征在于，

所述相关计算装置，在让所述压缩后的所述反射回波信号的包络线相对于所述压缩前的所述反射回波信号的包络线在三维方向上相互依次移动时，在预先设定的三维相关窗口的范围内在三维方向上依次移动。

6.根据权利要求5所述的超声诊断系统，其特征在于，

所述相关计算装置，在求取属于所述三维相关窗口内的所述压缩前后的所述包络线的相关系数成为最大的三维位置时，求出压缩后的所述包络线的自相关函数，让该自相关函数相对于所述压缩前的所述包络线在三维方向上依次移动，求出所述相关系数。

7.根据权利要求1所述的超声诊断系统，其特征在于，

还包括弹性模量计算装置，通过将所述对象划分为有限数量的单元，来建立至少三维有限元模型；根据用于建立所述模型的信息和所述应变分布，来计算弹性模量分布；

在所述显示装置上显示所述弹性模量分布。

8.根据权利要求7所述的超声诊断系统，其特征在于，

所述弹性模量计算装置，假定所述对象的组织为各向同性的弹性体和接近不可压缩性，将所述对象的组织划分为有限数量的矩形平行六面体单元，来建立三维有限元模型，并且假定在各所述单元内弹性模量、应力和应变都均匀，将所述应变分布用于弹性方程式，来计算弹性模量分布。

Images