相关申请的交叉引用
本申请要求于2011年8月19日提交的序列号为6,1525,378的美国临 时申请的优先权，其全部内容通过引用合并到本申请中。
技术领域
本发明总体上涉及医学成像，具体地涉及使用例如组织的对象的薄体 积(thinvolume)的超声成像来确定对象的粘弹性属性的弹性成像。
背景技术
弹性成像通常涉及对组织的机械属性如弹性、粘度、松弛时间、剪切 模量、孔隙率等进行成像的方法。通常用以下4个步骤来执行弹性成像：
1.通过在组织中引起一些形式的变形或运动来对组织进行激励；
2.使用医学成像设备来观察并且记录描绘组织中不同的位置在一段 时间内的运动的一系列图像；
3.根据这一系列图像来估计在不同位置和不同时刻组织的位移；以及
4.根据所估计的位移来估计组织的机械属性如弹性和粘度。
在本领域中通过组合不同类型的激励已经提出了具有不同成像形式 的许多弹性成像系统。已知的成像形式包括超声和磁共振成像(MRI， magneticresonanceimaging)、以及光学相干断层成像(OCT，optical coherencetomography)和X射线计算机断层成像(CT，x-raycomputed tomography)。也已经提出了用于估计位移、以及根据所估计的位移来估 计弹性和粘度的不同方法。
尽管也已经研究了在磁共振弹性成像(MRE，magneticresonance elastography)中使用的瞬态激励，然而本领域中的大多数MRE方法还 是使用稳态激励。尽管也已经研究了稳态激励，然而本领域中的大多数超 声弹性成像方法还是使用瞬态激励。
与MRI、OCT或CT相比，超声成像具有某些优势，如更低的成本、 更轻的重量和更容易的操作。然而，提供组织属性的成像的现有实时超声 弹性成像系统使用计算密集型技术，该技术需要复杂且昂贵的计算硬件或 者要求仅在2D成像平面内采集图像。这样的2D测量引入了误差或者降 低了测量弹性的绝对值的能力，而是仅测量整个图像的相对变化。因此， 现有的使用仅在2D成像平面内采集数据的探针的超声弹性成像系统易于 产生不准确的测量结果。
发明内容
根据本发明的一个方面，提供了一种用于对在组织的期望横截面中的 组织的至少一个机械属性进行成像的弹性成像方法。该方法包括：向组织 施加稳态振动以在组织中产生组织位移；通过在一段时间内采集表示组织 的薄体积的一组图像数据来在上述一段时间内对薄体积进行超声成像，薄 体积包括所述期望横截面和与所述期望横截面相邻的至少两个平面；对于 薄体积中的多个空间点计算地估计在上述一段时间内的多个相位同步的 位移，其中，在期望横截面和与所述期望横截面相邻的至少两个平面中的 每个平面上包括至少一个空间点；通过使用多个相位同步的位移来计算地 确定期望横截面上的组织的至少一个机械属性。与所述期望横截面相邻的 至少两个平面的数量选择成使得能够在实时刷新率范围内计算地确定所 述组织的至少一个机械属性，实时刷新率范围可以定义为每5秒至少1 个新帧。至少一个机械属性可以包括从以下属性中选择的任意一个或更多 个属性：绝对弹性、绝对剪切模量、绝对横波速度和绝对粘度。另外，通 过使用有限元方法、本振频率估计、行波扩展和直接转换(direct inversion)中的任意一个或组合根据多个相位同步的位移来计算至少一个 机械属性。
根据本发明的另一个方面，存在一种用于对在组织的期望横截面中的 所述组织的至少一个机械属性进行成像的弹性成像系统。该系统包括：配 置成生成稳态振动的至少一个振动源；配置成在一段时间内采集表示组织 的薄体积的一组图像数据的超声探针，该薄体积包括期望横截面和与所述 期望横截面相邻的至少两个相邻平面；与超声探针通信以接收来自超声探 针的图像数据并且包括处理器的电路，该处理器执行上述弹性成像方法； 以及与电路通信以接收并且显示期望横截面的一个或更多个图像以及在 实时刷新率范围内确定的组织的机械属性的显示设备。
超声探针可以是3D超声探针。更具体地，3D超声探针可以是包括 具有外部运动台的2D探针的机械3D探针。可选地，3D超声探针可以是 包括换能器的二维矩阵的多维探针。
振动源可以是电磁音圈并且可以安装在可调节臂上。臂可以安装在超 声机、病床或便携杆上。可选地，振动源可以配置为手持设备。同样，振 动源和超声探针可以集成在一起成为单个设备。
附图说明
图1是根据本发明的实施例的振动源的透视图。
图2是描绘了根据本发明的实施例的用于使用外部探针对对象进行 成像的配置的示意图，其中，振动源安装在连接至超声机的可调节臂上。
图3是描绘了根据本发明的实施例的用于使用内腔探针对对象进行 成像的配置的示意图，其中，振动源安装在连接至超声机的可调节臂上。
图4是描绘了根据本发明的实施例的用于使用外部探针对对象进行 成像的配置的示意图，其中，振动源安装在连接至病床的可调节臂上。
图5是描绘了根据本发明的实施例的用于使用外部探针对对象进行 成像的配置的示意图，其中，振动源安装在连接至便携杆的可调节臂上。
图6是描绘了根据本发明的实施例的用于使用外部探针对对象进行 成像的配置的示意图，其中，振动源直接使病床振动。
图7是描绘了根据本发明的实施例的用于使用外部探针对对象进行 成像的配置的示意图，其中，使用了手持式振动源。
图8是根据本发明的实施例的集成的振动源和超声探针的示意性视 图，其中，振动源位于探针的侧面。
图9是根据本发明的实施例的集成的振动源和超声探针的示意性视 图，其中，振动源配置成使探针的超声换能器振动。
图10是根据本发明的实施例的集成的振动源和超声探针的示意性视 图，其中，振动源配置成使惯性块振动。
图11是描绘了根据本发明的实施例的用于估计两个超声RF线之间 的位移的互相关方法的示图。
图12是描绘了根据本发明的实施例的在不同时间间隔对正弦信号的 采样的示图。
图13是描绘了根据本发明的实施例的使用超声探针从组织的一部分 依次采集RF线的示图。
图14是描绘了根据本发明的实施例的以两个不同的采样频率采样的 正弦信号的示图。
图15是描绘了根据本发明的实施例的当通过从组织的一部分依次采 集RF线来对组织进行成像时从所述一部分中的不同的点采集数据的时间 线的示图。
图16是描绘了根据本发明的实施例的机动的外部3D探针的示图。
图17是描绘了根据本发明的实施例的机动的内腔3D探针的示图。
图18是描绘了根据本发明的实施例的集成有振动源的矩阵阵列3D 探针的示图，其中，探针一次采集一个平面的3D数据。
图19是描绘了根据本发明的实施例的集成有振动源的机动的3D探 针的示图。
图20是描绘了激励与同步信号、以及图像采集与平面切换信号之间 的相关性的时序图。
图21是根据本发明的实施例的与外部2D探针一起使用以生成3D数 据的外部运动台的透视图。
图22是根据本发明的实施例的与内腔2D探针一起使用以生成3D数 据的外部运动台的透视图。
图23是描绘了根据本发明的实施例的用于执行对象的实时弹性成像 的超声弹性成像系统的框图。
图24是描绘了由弹性成像系统中的GPU实现的基于相关性的运动 估计算法的框图。
图25描绘了由弹性成像系统获取的示例性绝对弹性图像。
具体实施方式
方向术语如“顶部”、“底部”、“向上”、“向下”、“竖直地”和“横向 地”在下面的描述中仅用于提供相对参考的目的，而并不意在暗示对任意 元素在使用过程中或相对于环境如何显示进行任何限制。
本文中所描述的实施例总体上涉及一种弹性成像方法和系统，用于在 某一段时间内获取受激励的组织的超声图像，然后在实时刷新率范围内计 算地确定组织的一个或更多个机械属性。因为仅对受激励的组织的薄体积 进行成像并且处理，所以该方法可以实时地执行弹性成像。该薄体积包括 组织的期望横截面和与期望横截面相邻的至少两个平面。相邻平面的最大 数量选择成使得计算机系统能够在实时刷新率范围内计算地确定机械属 性。在此背景中，“相邻”表示与期望横截面分隔开并且在期望的横截面 旁边，具体地，合适的相邻平面可以紧邻着期望横截面或者在其与期望横 截面之间具有一个或更多个其他平面。同样在此背景中，实时刷新率定义 为每5秒至少一个新帧。
该方法总体上涉及：使用能够提供稳态振动的振动源来激励待成像的 组织；在所选择的一段时间内在至少两个不同的时刻采集受激励的组织的 以射频(RF)数据形式的超声图像；根据RF数据来估计受激励的组织 的位移，其包括针对超声图像中的多个空间点估计在一段时间内的相位同 步位移；以及根据相位同步位移来确定组织的机械属性如绝对弹性和粘 度。下面更详细地描述这些步骤中的每个步骤。
对体积使用位移测量是有利的，这是因为其允许在全部三个空间方向 上执行空间波长的测量。成像处理的速度与体积的大小成比例，所以较小 的体积由于其实现了更快的成像速度而是有利的。计算空间波长所需的在 任意一个空间方向上的数据点的最小数量是3个——2个数据点仅定义直 线并且因此不能用于估计由正弦波构成的波形的空间波长。3个数据点允 许进行曲率测量。换句话说，3个数据点允许计算可以用于测量空间波长 的二阶空间导数。
激励生成
现在参照图1至图10，提供了可以生成稳态激励以激励对象组织的 振动源和用于向对象施加振动的装置的不同实施例，它们用于与弹性成像 方法一起使用并且构成了弹性成像系统的一部分。
在图1示出的实施例中，通过能够生成在2Hz至1000Hz范围内的 稳态激励的机械振动源101来生成对象组织的激励。振动源101可以是电 磁音圈，或者是气动或液压驱动的。杆102和垫103用于将振动传递给组 织。
图2至图6示出了用于把振动源101施加于对象并且使用超声探针对 对象进行成像的不同的配置。在这些实施例中，振动源101安装在可调节 臂201上。可调节臂201上设置有手柄202。临床医生可以通过使用手柄 202并且将振动垫103布置成在期望的位置处以期望的角度与患者205接 触来调节臂201。在图2和图3示出的实施例中，臂201安装在超声机203 上。在图4示出的实施例中，臂201安装在病床204上。在图5示出的实 施例中，臂201安装在便携杆206上。图6示出了另一实施例，其中，振 动源101直接通过在其上有患者的床204来将振动传递给患者。
图7示出了又一实施例，其中，振动源101是手持式设备。操作者使 用一只手401握住振动源使得振动垫103与患者接触，并且使用另一只手 402握住超声探针403。可选地，一个操作者可握住振动源101，而另一 操作者使用超声探针来执行扫描。
图8是根据本发明的另一实施例的集成的振动源和超声探针的示意 图，其中，振动源101连接至超声探针403并且定位在探针403的侧面或 顶部。可选地，两个或更多个振动源可以连接至探针的两个或更多个侧面。
图9示出了根据本发明的另一实施例的集成的振动源和超声探针。在 这个实施例中，至少一个振动源101直接连接至换能器元件阵列502。所 生成的振动通过用作图1的垫103的在阵列502中的一个或更多个换能器 元件传递至组织。
图10示出了根据本发明的又一实施例的集成的振动源和超声探针。 在这个实施例中，振动源101连接至探针403，并且块503附接至杆102。 当通过振动源101在方向504上移动块503时，反作用力使超声探针403 在相反的方向505上移动。当探针403开始与组织接触时，探针403在方 向505上的运动传递至组织。块503可以通过弹簧或弹簧-阻尼系统连接 至探针403，弹簧或弹簧-阻尼系统能够优化，以在具体的频率范围内实现 更大的运动，如使用具体的中心频率和品质因数来调谐弹簧-块-阻尼系统 以共振。可选地，振动源101的杆102可以连接至探针403而不是其本体 101，实际上使用振动源本体作为反应块并且潜在地降低了总的探针质量。
在另一实施例中，振动源101的可调节臂201连接至探针403而不是 超声机203。
在本发明的又一实施例中，通过使用来自超声机的声辐射力来在组织 中内部生成激励。
为了增加组织中的运动的幅值，可以同时使用本文中描述的激励源中 的两个或更多个任意组合。
稳态激励
弹性成像方法的本实施例利用稳态激励。在这种类型的激励中，激励 信号的每个频率分量的幅值和相位随时间保持恒定，并且在瞬态效应已经 消失之后研究组织。
在弹性成像方法的一些实施例中，施加给振动源101的激励信号v(t) 可以是频率为f的纯正弦波的形式：
v(t)＝asin(2πft)(1)
当激励器垫103或探针403以这个频率振动，开始与组织接触并且在空间中保持固定时，在短时间段之后组织中的振动图案会达到稳态。在这个稳态下，组织中的每个点可以按正弦模式以同样的频率f振动，但是在每个位置处具有不同的相位和幅值。如果对象不移动，并且探针和激励器的位置在空间中固定，则相位和幅值将不随着时间改变，因而称为稳态。在这个设置中，每个组织点的位移可以用其幅值和其相位来表示为复数。用表示的复数表示定义为在点x处频率为f的位移的相位复矢量(phasor)。
在弹性成像方法的一些其他实施例中，激励信号v(t)可以是不同频率 的至少两个正弦波的总和：
v ( t ) = Σ i b i s i n ( 2 πf i t + θ i ) - - - ( 3 ) ]]>
在稳态下，组织中的每个点的位移u(t)也是与包括在激励中的频率相 同的频率的正弦波的总和：
并且每个具体的频率fi的位移的相位复矢量定义为：
位移估计
参照图11，弹性成像方法的所述实施例采用根据通过超声机在两个 不同的时刻收集的射频(RF)数据602和603来估计组织位移的处理。 这样的组织位移估计可以基于本领域中已知的处理，如由AManduca、R Muthupillai、PJRossman、JFGreenleaf和RLEhman在“Local wavelengthestimationformagneticresonanceelastography”(Proceedings InternationalConferenceonImageProcessing,1996,Vo1.3,pp.527- 530.)中描述的处理。
在不同的时刻采集两组RF数据602和603，要估计这两组RF数据之间的位移。RF数据602和603被分成块604，并且根据相似度测量执行对匹配的块的搜索。每个块604表示位于组织中典型的点607处的小体积的组织。
在弹性成像方法的一种实施例中，使用归一化相关性608作为相似度测量。相关性605的峰值示出了块之间的匹配(高相似度)，并且假设匹配的块之间的错位606为典型的点607的位移。在图12中示出了这个位移701，作为对于正弦激励的时间的函数。结果是在离散时间tk与tk+1之间对在处的小体积组织的位移的估计。
在位移跟踪的一种实现中，相对于给定的时刻t1来计算点的所有位移。这种方法被称为绝对位移跟踪，并且对于不同的时刻t1,t2,t3,...,tK生成位移
u ( x ‾ , ( t 1 , t 1 ) ) , u ( x ‾ , ( t 1 , t 2 ) ) , u ( x ‾ , ( t 1 , t 3 ) ) , ... , u ( x ‾ , ( t 1 , t K ) ) - - - ( 6 ) ]]>
在另一实现中，计算连续采集时间的数据之间的位移。这种方法被称 为相对位移跟踪，并且对于不同的时刻t1,t2,t3,...,tK生成位移
u ( x ‾ , ( t 1 , t 2 ) ) , u ( x ‾ , ( t 2 , t 3 ) ) , u ( x ‾ , ( t 3 , t 4 ) ) , ... , u ( x ‾ , ( t K - 1 , t K ) ) - - - ( 7 ) ]]>
可以通过取相对位移的累积总和来得到绝对位移，并且可以通过取绝 对位移的差来得到相对位移。
相位复矢量位移估计
为了使弹性成像方法的本实施例使用上述位移数据来确定组织的机械属性，位移或相位复矢量必须是相位同步的。可以执行相位复矢量位移估计步骤，以根据所计算的位移在激励频率fi求得组织中的每个点607的相位同步的相位复矢量(“同步相位复矢量”)这会是具有挑战性的计算，这是因为其中通过超声机来采集组织的超声图像的处理不是瞬时处理。由于声音在组织中有限的传播速度(大约1540m/s)，根据成像的深度、视场、分辨率、脉冲序列方案等，采集单一图像的处理可能需要数百微秒到数百毫秒。因此超声图像不是组织运动在单个时间瞬间的真正快照。这表示不同点的运动的采样时间不是相同的，而是分布在整个时间轴。不考虑采样时间的相互关系和激励频率的对相位复矢量的盲估计导致不同步的相位复矢量。在使用不同步的相位复矢量来估计组织弹性之前，必须对不同步的相位复矢量进行同步。
如下面将要讨论的，使用对位移进行采样的时刻tk与激励频率fi之间 的相互关系来计算组织中的不同的空间点的相位同步相位复矢量。
计算相位复矢量
如图12所示，激励具有单个频率f。在这个实施例中，点607的测量位移701接近频率为f的正弦波。患者和持有探针的医生的手的运动、以及测量噪声也可能促成点607的测量位移701。通过估计激励频率为f的相位复矢量，过滤了这些外来分量的影响。假设通过超声采集处理，在时刻t1,t2,...,tK对点607进行了成像，如图12所示。
在弹性成像方法的一种实施例中，超声机记录其采集的每个RF线的采集的时间戳。通过添加由x的深度确定的用以说明由超声发出的超声脉冲到达点对该点进行采样、并且然后返回超声探针所需要的时间的固定延迟，根据包含点的RF线的时间戳来求得对点进行扫描的时刻，即t1至tk。
可以对所估计的位移数据通过模型的最小方差拟合(leastsquares errorfitting)来求得位移的幅值和相位。例如，在绝对位移的情况下，模 型是：
其中，需要常量c来说明位移全部是相对于点在时刻t1的位移而言的情形：
c = - u ( x ‾ , t 1 ) - - - ( 9 ) ]]>
通过使用本领域中已知的最小二乘法来对下面的方程组进行求解，求 得幅值和相位：
1 s i n ( 2 π f ( t 1 - t 1 ) ) c o s ( 2 π f ( t 1 - t 1 ) ) 1 s i n ( 2 π f ( t 2 - t 1 ) ) c o s ( 2 π f ( t 2 - t 1 ) ) . . . . . . . . . 1 s i n ( 2 π f ( t K - t 1 ) ) c o s ( 2 π f ( t K - t 1 ) ) c α β = u ( x ‾ , ( t 1 , t 1 ) ) u ( x ‾ , ( t 1 , t 2 ) ) . . . u ( x ‾ , ( t 1 , t K ) ) - - - ( 10 ) ]]>
注意，相对于t1来估计相位，这表示假设在t1处相位为零。
在弹性成像方法的另一实施例中，激励包含多个频率，并且可以通过 下式求得不同频率的相位复矢量中的每个相位复矢量的幅值和相位：
对于不同的时间写出这个公式，得到：
a i = α i 2 + β i 2 , ]]>
同样相位都相对于t1而言的。
图13示出了另一实施例，其中，对点x进行扫描的时间点t1,...,tK在时 间上可以等间隔：
t2-t1＝t3-t2＝...＝tK-tK-1＝T(13)
许多超声机以按照901、902、903、904的顺序逐线扫描组织，直到 采集到待成像的整个扇区905，并且从第一条线重复同样的序列906。在 这种情况下，超声机不需要记录所有采集的RF线的时间戳。知道时间间 隔T就足够了。由N和τ来确定时间间隔T，其中，N是扇区中线的数量， τ是扫描每条线的时间：
T＝Nτ(14)
扫描每条线的时间τ必须足以让超声脉冲到达在深度d处的待成像的 最深的点908并且返回探针：
τ > 2 d 1540 m / s - - - ( 15 ) ]]>
总之，可以通过知道T来根据矩阵方程(10)或(12)确定相位复 矢量。
在弹性成像方法的另一实施例中，使用基于傅里叶的方法来计算相位 复矢量，连续扫描之间的时间间隔T或激励的频率fi调整成使得：
1)激励的频率共享公共周期f-1(注意，仅在这种情况下稳态激励变 为周期性的，周期为f-1)：
fi＝kfk∈{1,2,3,...}(16)
2)时间间隔T是激励的公共周期f-1的有理分式：
T = m n f - 1 , m , n ∈ { 1 , 2 , 3 , 4 , ... } - - - ( 17 ) ]]>
换句话说，从m个激励的周期得到n个样本。
为了这个目的，应当使用以下技术中的一个或组合来调整激励频率 fi或时间间隔：
-改变该扇区中线的数量N；
-改变成像的深度d、或点908；
-在采集序列中增加等待时间；
-使对线的采集与具有可调周期的外部脉冲源同步。这个脉冲源 可以是生成激励的同一个源。
如上所述，T与fi之间的这样的关系使得能够使用基于傅里叶的方法 来计算相位复矢量。图14示出了示例702，其中，T是激励周期的1/5。 在这种情况下，还可以采用位移的离散傅里叶变换来计算相位复矢量。在 这种情况下取得相位复矢量的另一种数学上等价的方法是乘以exp(j2πft) 并且关于时间求和：
U f = 1 5 Σ i = 1 5 exp ( j 2 π f · i T ) u ( x , ( t 1 , t i ) ) - - - ( 18 ) ]]>
图14还示出了示例703，其中，T是激励周期的五分之六。在这种 情况下，以与在702的情况下完全相同的相对相位但是在运动的不同周期 来对运动进行采样。然而，在这种情况下也可以使用同样的技术(离散傅 里叶变换、以及乘以exp(j2πft))来计算相位复矢量。即使在采样率低于激 励频率703时也能计算相位复矢量的情形是频率域中的激励的窄带宽的 结果。在本领域中这项技术被称为“带通采样”，并且RSinkus、J Lorenzen、DSchrader、MLorenzen、MDargatz和DHolz在 “High-resolutiontensorMRelastographyforbreasttumourDetection” (Phys.Med.Biol.45,2000,pp.1649-1664)中解释了这项技术。
对相位复矢量进行相位同步
因为相对于点在时刻t1901的位移取得每个点607的位移(如在图13中示例性地示出)，所以相位复矢量的所计算的相位是相对于t1而言的。对于两个不同的点和参考时刻通常是不同的：
因此相位差不表示所述两个点的运动的实际相位差。其还包括由参考时刻的差异t1-t′1引起的相位差。可以通过对这个外加的相位差进行补偿来使相位复矢量同步(相位同步)：
在记录由超声机采集的每个RF线的时间戳的弹性成像方法的实施例 中，使用时间戳来计算参考时间的差异并且使用上述公式使相位复矢量同 步。
在逐线901、902、903采集图像并且一条线和下一条线的采集之间的 延迟为常量τ的弹性成像方法的实施例中，如图13和图15所示，执行相 位补偿不需要知道RF线的绝对时间戳。
如图15所示，可以绘制采集的时间线图，根据采集的时间线图可以计算任意两个点的采集之间的延迟。为了说明的目的，在扫描线序列和时间线图二者上标记了9个点1005至1013。例如对于在两个连续的线上的相同深度处的两个点1005和1008，时间延迟等于τ，并且这两个所述点的对应的相位复矢量可以通过下式来同步：
在一些特殊的情况下，两个点之间的相位差可以变为零，换句话说， 点已经同步并且不需要在他们之间进行相位补偿。为了做到这一点，使得：
exp(j2πf(t1-t′1))＝1
f(t1-t′1)＝mm∈{0,1,2,3,...}(22)
t1-t′1＝mf-1m∈{0,1,2,3,...}
为了利用这个属性，调整连续的线的扫描之间的时间间隔τ和激励的 频率fi，使得：
1)激励的频率共享公共周期f-1；
2)时间间隔τ是激励的公共周期f-1的整数倍。
估计位移的两个分量
在本领域中已知使用常规的声束偏转来在不同角度采集RF线以计算 位移的两个分量(轴向和横向)的方法，例如由RZahiri-Azar、ABaghani、 SESalcudean和RRohling在“2DHighFrameRateDynamic ElastographyUsingDelayCompensatedandAngularlyCompounded MotionVectors:PreliminaryResults”(IEEEtrans.Ultrason.,Ferroelect., Freq.contr,Vol57,No11,Nov.202010,pp.2421-2436)中讨论的。
在弹性成像方法的一些实施例中，至少在两个不同的角度采集RF线， 并且沿两个或更多个方向计算的位移的相位复矢量相结合以计算轴向和 横向位移。
根据相位复矢量位移计算绝对弹性和粘度
通过激励源生成并且传递至组织的在2Hz至1000Hz范围内的机械 振动作为横波传播通过组织。假设为纯粹的弹性模型的这些波的传播的控 制方程是：
ρ ∂ 2 ∂ t 2 u ( x , y , z , t ) = μ ( ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 + ∂ 2 ∂ z 2 ) u ( x , y , z , t ) - - - ( 23 ) ]]>
其中，ρ是密度，μ是剪切模量，点x的坐标用(x,y,z)来表示。因为生物 组织是几乎不可压缩的，所以其弹性模量E等于3μ，并且因此弹性和剪 切模量的测量是等同的。对于频率为f的稳态激励，控制频率为f的位移 的相位复矢量的方程变成：
- ρ ( 2 π f ) 2 U f ( x , y , z ) = E 3 ( ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 + ∂ 2 ∂ z 2 ) U f ( x , y , z ) - - - ( 24 ) ]]>
实际上生物组织不完全是弹性的。对组织的粘弹性性质进行建模的一 种方法是考虑μ为复数，其是频率的函数，
μ＝μ(f)＝μRE(f)+jμIM(f)(25)
其中，μRE与剪切弹性相关联，并且μIM与剪切粘度有关。
用于根据相位复矢量来计算弹性和粘度的绝对值的方法在本领域中 是已知的，并且可以在本实施例的弹性成像方法中实现以计算、但不限于 计算弹性、剪切模量、剪切波速和粘度中的任意一个。
例如，波动方程的直接反演的方法，如由RSinkus、JLorenzen、D Schrader、MLorenzen、MDargatz和DHolz在“High-resolutiontensor MRelastographyfor20breasttumourdetection”(Phys.Med.Biol.45, 2000,pp.1649-1664)中所描述的，可以用于通过转换波动方程(24)来 求得弹性，
E = - 3 ρ ( 2 π f ) 2 U f ( ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 + ∂ 2 ∂ z 2 ) U f - - - ( 26 ) ]]>
可以对Uf进行滤波，例如，可以在被用在公式(26)中之前进行低通滤 波以去除噪声，或进行高通滤波以去除零频率分量，或者进行上述两种处 理。
在本振频率估计(LFE，localfrequencyestimation)中，以在现有 技术中已知的方式如由AManduca、RMuthupillai、PJRossman、JF Greenleaf和RLEhman在“Localwavelengthestimationformagnetic resonanceelastography”(ProceedingsInternationalConferenceonImage Processing,1996,Vo1.3,pp.527-530)中所描述的，对相位复矢量进行定 向滤波并且使相位复矢量通过一组滤波器组。假设本振频率在所使用的滤 波器对的带宽范围内，每对滤波器的输出的比率给出对信号的本振频率的 估计。使用该估计的加权总和作为对本振频率的测量。本振(空间角)频 率k通过下式与弹性相关：
E = - 3 ρ ( 2 π f ) 2 1 - k 2 ( x , y , z ) - - - ( 27 ) ]]>
此外，作为算法的一部分，Uf可以在用在滤波器组中之前进行适当 的滤波。
在行波扩展(TWE，travellingwaveexpansion)中，使用具有不同 的幅值和相位的不同方向上的行波的总和作为用本振(空间角)频率作为 模型参数的模型。然后通过求得最佳本振频率来使模型本振地拟合到相位 复矢量。可以如在(27)中以在本领域中已知的方式根据所估计的本振(空 间角)频率来计算弹性。
在有限元方法(FEM，finiteelementmethod)中，使正演模型(forward model)迭代地拟合到相位复矢量，每次调整机械属性的假设分布，直到 得到机械属性的实际绝对值；由HEskandari、SESalcudean、RRohling 和JOhayon在“ViscoelasticCharacterizationofSoftTissuefrom DynamicFiniteElementModels”(PhysicsinMedicineandBiology,Vol. 53,No.22,pp.6569-6590,Nov.2008)中描述了示例方法。
绝对弹性成像的体积成像
如公式(24)所描述的，弹性E的准确估计需要在3D体积中测量位 移相位复矢量。换句话说，振动实际上在生物组织的3D体积中传播。为 了测量这样的振动的空间频率，需要在3D体积中对其进行测量。因此， 在弹性成像方法的本实施例中，使用单个3D超声探针从组织的3D体积 收集超声数据(体积成像)。
一般来说，合适的3D超声探针将声波发送至限定待成像的组织的3D 体积中。来自反射的声波的接收数据创建了解剖结构(anatomy)的体积 数据集(通常缩写为“体积”)，这与创建截面的图像的2D超声探针不同。 至少可以通过以下两种已知方法来实现实时3D超声成像：
1)机械扫描：通过将2D探针与用于快速移动2D探针的机动机构组合使 得2D图像重复地扫描通过感兴趣的体积来构建专用3D探针。通常以摆 动方式来实现重复扫描，其中每次摆动产生3D体积。因为探针运动是受 控的，所以由于每次摆动的2D图像集之间的空间关系是已知的，并且图 像重建成3D笛卡尔坐标体积。在下文中该设备称为机械3D探针；
2)多维阵列：在不使用机动机构而是使用换能器的二维阵列来扫描感兴 趣的3D体积的情况下创建专用探针。体积采集的速度通常高于机械探针， 但是探针的复杂度有所增加，并且图像质量可能较差。该探针被称为多维 探针。
在弹性成像方法的本实施例中使用的3D超声探针可以是本领域中已 知的或在下面描述的机械3D探针或多维3D探针。
合适的已知机械3D探针的示例是用于通用电气公司(GE医疗集团， 查尔方特-圣贾尔斯，英国(GEHealthcare,ChalfontSt.Giles,United Kingdom))的Voluson730超声机的RAB2-5H46701M。合适的已知多 维探针的示例是用于飞利浦iU22超声机(飞利浦医疗保健，安多佛，马 萨诸塞州，美国(PhilipsHealthcare,Andover,Massachusetts,USA))的 X7-2。
图16至图22中示出了根据不同实施例的机械3D探针和多维3D探 针。图16示出了外部3D机动线性阵列超声探针。内腔3D机动超声探针 在图17中示意性地示出并且能够根据对多个平面的成像来采集3D体积 图像。图18示出了在探针内集成有振动源的矩阵阵列3D超声探针。图 19示出了根据本发明的实施例的也集成有振动源的机动的3D超声探针。 图21和图22示出了使用具有外部运动台的2D探针的一对机械3D超声 探针。无论在弹性成像方法的本实施例中使用哪种3D超声探针，都通过 以下步骤来采集体积数据：
(a)将组织的体积分为多个平面。
(b)在一段时间内对每个平面进行成像使得多次扫描平面中的每个 点；以及
(c)对所有平面重复该处理。
如上所述，3D超声探针采集包括期望横截面和至少两个相邻平面的 组织的薄体积的图像。例如，在图16、图17和图18中，正在成像的薄 体积是所有平面的凸形外壳，期望横截面示出为1204，相邻平面示出为 1203、1202和1205、1206，1203、1202和1205、1206都与横截面1204 相邻，平面1203和1205紧邻横截面1204。在一些情况下，需要调整相 邻平面之间的间隔，以更准确地捕捉在立面方向上的相位复矢量的变化。
如图18所示的3D超声探针包括用于操纵超声成像平面并且采集多个平面1202、1203、1204、1205、1206以形成体积的换能器的矩阵阵列1201。在这个探针中，电子地切换成像平面，因此与其他类型的探针相比，能够更快地执行采集，并且能够达到更高的整体帧速率。图18的实施例中使用的振动机构是与图10中示出的振动机构基本相同的块503。具有换能器的矩阵阵列的该超声探针可以使用合成孔径波束形成，以重建覆盖3D体积的RF线。由ArendtJensen、SvetoslavIvanovNikolov、KimGammelmark和MortenPedersen在“Syntheticapertureultrasoundimaging”(Ultrasonics44,2006,pp.e5-eI5)中描述了这样的探针的示例。可选地，具有换能器的矩阵阵列的该超声探针可以使用微波束形成，以重建覆盖3D体积的RF线。由ChristopherHall在“4-DimensionalUltrasonicImaging”(AdvancesinHealthcareTechnologyCareShapingtheFutureofMedical,PhilipsResearchBookSeries,2006,Volume6,No.2,pp.99-116)中描述了这样的探针的示例。
图19中示出的超声探针包括线性换能器阵列1401和用于使阵列扫过 一扇形角的电机1402。在这个探针中，探针从一个面到下一个面的运动 在超声RF线的采集处理中产生了新的时延。这是探针内的机械运动停止 所需要的时间。如果不同平面上所采集的RF线的时间戳是可用的，则可 以如公式(20)中一样使用该信息来使相位复矢量同步。振动机构是与图 10中示出的振动机构基本相同的块503。图20示出了根据本发明的实施 例的激励和同步信号、图像采集和平面切换信号之间的相关性。如图20 所示，激励(V)的频率共享公共周期f-1，并且以同样的周期f-1生成同 步信号(S)。同步信号S用于使每个成像面上的采集的开始1504与振动 源同步。图20示出了用于这个目的的一个示例时序图。在这个图中，同 步信号S与振动信号V是同步的。通过采集信号A来控制RF线的采集， P是例如前往电机1402的控制器的平面切换信号。在已经针对一个平面 采集了数据1501之后，命令电机移动至新的平面1502。在电机已经稳定 在新的位置1503之后，系统等待同步信号S上的下一个上升沿，以启动 下一个平面的采集1504。
可选地，包括具有外部运动台的2D探针的机械3D探针可以用于对 3D体积进行成像。例如，图21中示出了根据一种可选实施例的与外部 2D探针一起使用以生成3D数据的外部运动台1601，而图22示出了根据 另一实施例的包括与内腔2D探针一起使用以生成3D数据的外部运动台 1603的机械3D探针。外部运动台1601、1603移动固定至所述运动台的 2D探针403，以采集在图21中标记为1602以及在图22中标记为1604 的多个平面。在这些实施例中，不同平面的采集时间可以使用时间戳来记 录或者与激励同步，以便使不同平面的位移相位复矢量同步。
不像位于平行的平面内等间距的笛卡尔位置的MRE中的位移测量的 规则间距，使用超声的位移测量通常不在等间距的笛卡尔位置处。从图 16和图17中的平面的不平行间隔可以看到这一点。在这样的超声数据的 情况下，因此需要适应测量的具体间隔，使得在笛卡尔坐标系中执行绝对 弹性和粘度的计算，如在公式(23)、(24)、(26)和(27)中所示。假设 相位复矢量的相位同步包括测量的位置，则相位同步步骤应该包含超声波 束和平面的实际空间位置。
实时地执行计算
现在参照图23，提供了一种执行上述弹性成像方法的弹性成像系统。 该系统包括耦接至超声探针403的超声机1702；耦接至振动源101的激 励信号发生器1704(或者称为波形发生器)；与超声机通信并且包括具有 存储器的处理单元1708(GPU)的电路，存储器上编程有由GUP1708 可执行的实施弹性成像方法的步骤和指令；以及用于显示成像的组织和该 组织的由弹性成像方法确定的机械属性的显示器1713。
GPU1708可以具有并行处理能力，并行处理能力可以用于相对于RF 数据的采集时间降低用于估计机械属性的计算时间。超声机1702可以从 探针403采集RF数据1707，并且将RF数据1707发送至GPU1708。 GPU1708具有以在本领域中已知的方式生成B模式图像1715并且把B 模式图像1715发送至显示器1713的B模式生成软件模块1709。GPU1708 还具有计算至少一个轴向上的组织位移的位移估计模块1710、相位复矢 量生成软件模块1711和计算组织的机械属性的软件模块1712。在一种实 现中，GPU1708将机械属性的图像发送至显示器1713，机械属性的图像 叠加在B模式图像上并且向用户显示。波形发生器1704生成用于振动源 101的信号1705，并且还生成用于超声机1702的触发信号1706。
并行处理GPU可以包括数百个处理单元。这样的GPU能够同时执行数千个线程(顺序程序)。弹性成像方法中定义的计算是高度并行的，因此适于由并行处理GPU来执行。按照NVIDIA公司的“NVIDIACUDAC编程指南”和“用于CUDA架构的OpenCL编程”中指定的方法，由提供的CUDATM库可以用于编译GPU的多线程程序。由VictorPodlozhnyuk在“Imageconvolutionwithcuda”(Technicalreport,NVIDIACorporation,2007)中提供了如何使用CUDATM库来进行图像滤波的详细的示例。
图24是示出了执行弹性成像方法的弹性成像系统的另一实施例的框 图，该实施例执行基于相关性的运动估计算法以执行弹性成像方法。在该 系统中，在单独的线程中执行计算每个对应的两个RF线1803与1804之 间的位移估计的互相关的每个步骤。两个内核1805和1806在主机处理器 (CPU)1700上运行。互相关内核1805将RF线1803和1804从CPU 存储器1802复制1807到GPU1708的全局存储器1809。数据被分成窗口。 然后，根据两个RF线1801和1808求得多个相邻的窗口(例如4个窗口) 的位移所需的数据被复制1810到GPU1708中的GPU块1811，并且对 其他组的四个窗口重复同样的到其他块的过程，使得在GPU1708的块 1811至块1820之间分配求得互相关的任务。在每个块中，128个线程1813、 1814、1815计算每个窗口的不同变化量的互相关。所有的线程并行运行 并且在计算所有块的互相关的时刻被同步。互相关被写入全局存储器 1809。
第二内核1806——或称为峰值查找内核——在块之间分配互相关并 且求得每个窗口的互相关的峰值，在结束时同步以确保已经求得所有峰 值，并且写回GPU1708上的全局存储器1809。然后，最后的位移被传送 回CPU存储器1802。
在另一实施例中，使用CULATM库来执行快速多线程矩阵代数。使 用CULATM库通过对公式(12)进行求解来根据时间位移计算相位复矢 量。在一个矩阵中把不同点的时间位移分为一组(公式12的右手边)，使 得对CUDATM求解程序的单一调用可以同时提供所有点的相位复矢量。 优选地，由内核来使不同点的相位复矢量同步。内核在单独的线程中执行 每个点的相位补偿。
另外，可以使用由提供的CUFFTTM来执行LFE滤波所需的傅里叶变换。可选地，可以使用CULATM库来对FEM的矩阵方程进行求解。
也可以以有效的方式在现场可编程门阵列(FPGA， field-programmablegatearray)上或者在数字信号处理器(DSP，digital signalprocessor)上执行上述计算中的一些。
薄体积的实时采集
有意义的是，实时地或者以每五秒至少一帧以及优选地接近或者在每 1/30秒一帧的高帧速率沿着组织的横截面来采集描绘组织的机械属性(例 如弹性)的2D图像。放射科医师、超声技术人员、医生、护士和超声机 的其他用户通常使用实时2D成像来定位恶性肿瘤和进行诊断。本实施例 实现了弹性成像方法的实时计算处理并且从而使用户可以在缓慢地移动 探针时观察所显示的图像的变化。
在弹性成像方法的所描述的实施例中，从组织的期望横截面和与期望 横截面相邻的至少两个其他平面采集数据。这样的数据的体积称为薄体 积，这是因为体积在期望平面上的长度是体积在第三(平面外)维度上的 长度的至少3倍。薄体积可比更大的体积更容易进行实时采集。然后在所 述期望横截面上使用数据来估计例如弹性的机械属性，并且向用户示出作 为2D图像的结果。
在体积内使用位移相位复矢量的测量是有益的，这是因为其允许在全 部三个空间方向上执行空间波长的测量。成像处理的速度与体积的大小成 比例。较小的体积是有益的，这是因为其能够实现更快的成像率。计算空 间波长所需的在空间方向中任意一个空间方向上的数据点的最小数量是 三个，这是因为三个数据点允许计算可以用于测量空间波长的二阶空间导 数。因此，使用薄体积和同步采集以及高效的相位矫正，本发明的实施例 能够实现实时的或高帧速率的超声弹性成像。
图25示出了通过弹性成像系统的实施例收集的示例绝对弹性图像。 该系统用于使用以下成像参数对组织模型进行成像：扫描每条线的时间τ 等于125微秒，125微秒使得能够进行深度为80mm的成像。扇区中的 线的数量N等于32，线间距为1.2mm。所得到的扇区宽度等于38.4mm。 每个平面扫描30次。平面的采集时间是30×32×125μs＝120ms。采 集相邻平面之间间距为1.2mm的七个平面，其需要7×120ms＝0.84秒。 使用将晶体从一个平面移动至下一个平面需要20ms的机动3D探针(图 19)。因此每个帧的总采集时间是：0.84s+7×20ms＝0.98秒。数据的 处理没有对采集时间(小于100ms)造成严重的系统开销。在该示例中 薄体积的尺寸在感兴趣的平面中是38.4mm乘以80mm，对于第三(平 面外)维度是8.4mm。当数据变得可用时对位移和相位复矢量进行估计。 感兴趣的平面是中央平面，并且所采集的七个平面的中央平面中的弹性使 用LFE方法来估计并且实时地示出。所得到的帧速率是每秒一帧。图25 并排示出了叠加在B模式图像上的所得到的弹性图像的一帧和B模式图 像。因而可以在不到一秒钟内完成对在通常使用的大约80mm乘以40 mm的合理的感兴趣区域中的组织弹性的绝对值的成像。这意味着弹性成 像系统可以以大于每秒1帧的刷新率显示绝对弹性图像。
标准超声检查中所使用的其他图像尺寸涉及最多200mm的深度以 及每扇区64线的线密度。如以前一样假设7个平面的采集，弹性图像刷 新率将近似是每5秒一个新帧。
薄体积中所使用的平面的数量和所得到的成像帧速率之间总会存折 衷。确定立面方向上相位复矢量的变化的精确度取决于所使用的点的数 量。三个平面是最小值并且将产生最高弹性成像帧速率。更多的平面会提 供另外的所测量的相位复矢量位置并且因此会提高弹性图像精确度。已经 发现，使用标准的弹性模型，七个平面提供了图像精确度与成像帧速率之 间的良好折衷。
尽管已经描述并且示出了具体的实施例，然而这样的实施例应当认为 仅仅是说明性的而不作为对本发明的限制。