专利名称：采用两种尺度的生物组织位移估计方法
技术领域：
本发明属于超声弹性成像技术领域，特别涉及生物组织位移估计方法。
背景技术：
生物组织弹性模量的变化通常与其病理现象有关。例如，恶性的病理损害，例如乳房硬癌、前列腺癌、甲状腺癌及肝转移等，通常表现为硬的小结。乳房硬癌是乳腺癌的最常见形式，大约占乳腺癌总数的四分之三，由于其基质密度增大而表现为致密的硬块。而其他类型的乳腺癌如导管内癌和乳头状瘤则表现为柔软的组织，良性的乳腺纤维囊性病也很少表现为硬块。
生物组织的弹性模量信息对于疾病的诊断过程具有重要的参考价值。然而，包括X射线成像、超声成像、计算机断层成像(CT)和磁共振成像(MRI)等在内的传统医学成像模态都不能直接提供关于弹性模量这一组织的基本力学属性的信息。1991年，J.Ophir提出超声弹性成像(ultrasound elastography)的方法，对组织的弹性模量分布进行定量估计、成像。目前，超声弹性模量已经成为医学超声成像的一个研究热点，广泛应用于乳房、前列腺、动脉粥样斑块、心肌动力学以及高强度聚焦超声与射频消融引起的损害(lesion)的检测与评估。
超声弹性成像的基本原理为将超声探头嵌于一块挤压平板中，沿着探头的纵向压缩组织，分别采集组织压缩前、后的射频信号；组织被压缩时，组织内将会产生一个沿压缩方向的应变，如果组织内部弹性模量分布不均匀，组织内的应变分布也会有所差异；弹性模量较大的区域，引起的应变比较小；反之，弹性模量较小的区域，相应的应变比较大。通过一些方法估计出组织内部不同位置的位移，从而计算出组织内部的应变分布情况，用来间接描述组织内部的弹性模量分布，从而描述组织的生理、病理状态。
对于二维超声弹性成像，一般采用线阵的B型超声探头，采集组织压缩前、后的探头每一条扫描线的射频信号，分别进行上面描述的位移估计，从而计算出每一条扫描线对应组织的一维应变分布。最后把所有扫描线对应的一维应变分布按扫描线顺序组成一个二维应变分布，以灰度图或者伪彩图的形式表示，用来间接描述组织内部的弹性模量分布。
一般的超声弹性成像方法包括以下步骤1.利用商用B型超声仪器(一般采用线阵探头)得到待测生物组织(一般为人体组织，也可以为动物组织，以下简称组织)压缩前的一幅数字化的二维射频信号(可以采用模拟射频信号输出端接信号放大器，再接高速数据采集卡，得到数字化的二维射频信号；也可以在数字化B型超声仪器上直接得到数字化的二维射频信号)；2.手持该B型超声仪器的探头或者利用步进电机或者螺旋装置驱动该探头，沿着探头的纵向对该组织施加一个微小的挤压(对组织施加的压缩量一般控制在为1％的数量级)，得到组织压缩后的一幅数字化的二维射频信号；3.从步骤1和2的得到的组织压缩前、后的二维射频信号中分别取出第一条扫描线的数据，设为s1(n)和s2(n)，n表示该两条扫描线上的数据序号，1≤n≤nmax，n的最大值nmax由该B型超声仪器的探查深度、发射的超声波在组织中的传播速度以及射频信号的采样频率决定；4.从该扫描线数据s1(n)中取一小段长度为T的数据d1，其数据个数为U，U＝round(T×U1)，其中，T的单位为mm，U1代表1mm的组织对应的数据个数，由发射的超声波在组织中的传播速度以及射频信号的采样频率决定，round(·)代表四舍五入的取整操作，该数据d1的序号从n1到n1+U-1，n1可在1≤n1≤U的范围内选择；在τ1到τ2确定的搜索范围内求该小段数据与扫描线数据s2(n)的互相关函数R(τ)，计算公式如下R(τ)=Σi=n1n1+U-1s1(i)s2(i-τ)Σi=n1n1+U-1s12(i)·Σi=n1n1+U-1s22(i-τ)(τ1≤τ≤τ2)]]>其中i为计算过程中表示数据序号的循环变量，τ1为0，τ2为对组织施加的压缩量，以采样数据的个数表示；为了提高位移估计的精度，一般还需要对计算得到互相关函数进行插值，如抛物线等常规插值方法；5.确定步骤4得到的互相关函数R(τ)的最大值对应的位置t1t1就是数据d1在组织压缩后的位移(即s1(n)中的序号从n1到n1+U-1的小段数据d1的在组织压缩后移动到s2(n)中的序号从n1-t1到n1+U-1-t1的位置)；6.依次从扫描线数据s1(n)中取一小段长度为T即数据个数为U的数据d2、d3、…、dN，每段数据的序号依次错开V个采样数据，1≤V≤U，直到再错开V个采样数据将超出s1(n)的范围，按步骤4、5相同的方法依次得到各段数据对应的位移t2、t3、…、tN，其中N为小段数据的总数；则位移序列t1、t2、…、tN为第一条扫描线数据s1(n)对应的组织的位移估计；7.利用与步骤3-6相同的方法，依次得到第2、3、…、M条扫描线数据对应的组织的位移估计，其中M为表示探头的扫描线总数，由探头决定；8.对第一条扫描线数据s1(n)对应的组织的位移估计序列t1、t2、…、tN求差分，得到组织第一条扫描线s1(n)对应组织的应变分布，计算公式如下，ϵ1=t2-t1V,ϵ2=t3-t2V,···,ϵN-1=tN-tN-1V]]>其中，ε1、ε2、…、εN-1分别为d1、d2、…、dN-1对应的组织应变；9.利用与步骤8相同的方法，依次得到组织第2、3、…、M条扫描线数据对应的组织的应变分布；10.将步骤9得到的M条扫描线数据对应的应变分布，按照扫描线的顺序组合成一个二维数据，并以灰度图或者伪彩图的形式表示出来，就得到组织的二维应变分布图。
在超声弹性成像中，关键的问题在于对组织的位移分布进行估计，也就是上面描述的方法的步骤3-7。互相关函数的值越大，说明压缩前、后的小段数据吻合得越好，互相关函数的最大值位置代表了压缩前的小段数据在压缩后对应的位置，从而可以求出该小段数据的位移，也就是该小段数据对应的组织的位移。
在超声弹性成像中，选择的压缩前信号中的小段数据用来跟踪对应的小段组织的位移，称为跟踪波段。其长度称为跟踪波段长度，或者尺度。选取合适的尺度在超声弹性成像的位移估计中非常重要。当对组织施加的压缩量比较小(如小于1％)的时候，尺度越大，包含的信息越多，位移估计的精度越高，对随机噪声的干扰也越不敏感。当对组织施加的压缩量比较大(如1-5％)的时候，尺度越大，由于组织被压缩，组织压缩前、后信号重合的部分越小，利用互相关函数最大值位置进行位移估计的精度越低，同时，尺度越大，又越不容易受到随机噪声的干扰；而当尺度较小的时候，组织压缩前、后信号重合的部分较大，位移估计的精度较高，然而此时的位移估计又对随机噪声的干扰比较敏感。
为了增加组织压缩前、后信号重合的程度，一般采用线性插值的方法，把被“压缩”的信号“拉伸”(stretching)成与原信号同长，从而增加压缩前后波形的重合度。在波形拉伸中，拉伸系数与对组织施加的压缩量或应变(以百分比表示)有关。当拉伸系数接近组织的应变时，波形的重合度最大，从而互相关函数最大值最大。但是组织的应变在时延估计前是未知的，拉伸系数的选择只能根据施加的位移量和组织的深度作大致判断。并且由于组织的不均匀性，组织内部应变的分布也不均匀，因此特定的拉伸系数在组织的某些部分可能拉伸不足，在某些部分又拉伸过量，这反过来又引入误差。因此，有人提出一种对不同位置的波段进行自适应拉伸系数的方法，使得拉伸后的波段与压缩前波段的互相关函数最大值最大，然后利用拉伸系数直接估计出该位置组织的应变。但是这种自适应拉伸的方法计算时间较长，在实际应用中可能受到限制。
从实际应用的角度，希望对组织施加的压缩量大一些好，这样可以更能“看清”组织内部弹性模量差异引起的应变分布差异的细节。一般采用多次压缩(multi-compression)的方法，将大位移分解为小位移之和，逐次对组织施压一个小的位移量，计算出相应的应变分布。最后将多幅应变分布叠加，得到大位移时的应变分布。然而，该方法增加了数据采集时间。
发明内容
本发明为了解决对组织施加的压缩量比较大的时候大尺度和小尺度的矛盾，提出一种采用两种尺度的组织位移估计方法。在组织压缩量比较大时，将选择大尺度或小尺度进行组织位移估计的两种方法结合起来，可对这两种组织位移估计的优缺点实现扬长避短，提高组织位移估计的精度。
本发明提出的一种采用两种尺度的组织位移估计方法，包括以下步骤1.从组织压缩前、后的二维射频信号中分别取出第一条扫描线的数据，设为s1(n)和s2(n)，n表示该两条扫描线上的数据序号，1≤n≤nmax，n的最大值nmax由该B型超声仪器的探查深度、发射的超声波在组织中的传播速度以及射频信号的采样频率决定；2.从该扫描线数据s1(n)中取一小段尺度为Ta的数据d1，其数据个数为Ua，Ua＝round(Ta×U1)，其中，Ta的单位为mm，U1代表1mm的组织对应的数据个数，由发射的超声波在组织中的传播速度以及射频信号的采样频率决定，round(·)代表四舍五入的取整操作，该数据d1的序号从n1到n1+Ua-1，n1可在1≤n1≤Ua的范围内选择；在τ1到τ2的搜索范围内求该小段数据与扫描线数据s2(n)的互相关函数R(τ)，计算公式如下R′(τ)=Σi=n1n1+U-1s1(i)s2(i-τ)Σi=n1n1+U-1s12(i)·Σi=n1n1+U-1s22(i-τ)(τ1≤τ≤τ2)]]>其中i为计算过程中表示数据序号的循环变量，τ1为0，τ2为对组织施加的压缩量，以采样数据的个数表示；为了提高位移估计的精度，最好还需要对计算得到互相关函数进行插值处理，如抛物线插值等常规插值方法；3.确定步骤2得到的互相关函数R′(τ)的最大值对应的位置t′1，t′1为利用尺度Ta进行位移估计时数据d1在组织压缩后的位移(即s1(n)中的序号从n1到n1+Ua-1的小段数据d1的在组织压缩后移动到s2(n)中的序号从n1-t′1到n1+Ua-1-t′1的位置)；4.从该扫描线数据s1(n)中的数据段d1的中央位置再取尺度为Tb、序号从n1+round(Ua-Ub2)]]>到n1+round(Ua-Ub2)+Ub-1]]>的一小段数据e1，该数据个数为Ub，Ub＝round(Tb×U1)，其中，Tb的单位为mm，并且满足Tb＜Ta，最好取116Ta≤Tb≤12Ta;]]>在τ3到τ4确定的搜索范围内求该小段数据与扫描线数据s2(n)的互相关函数R″(τ)，计算公式如下R′′(τ)=Σi=m1m1+Ub-1s1(i)s2(i-t′1-τ)Σi=m1m1+Ub-1s12(i)·Σi=m1m1+Ub-1s22(i-t′1-τ)(τ3≤τ≤τ4)]]>其中τ3和τ4为组织位移精细调整的范围，满足如下关系τ4＜τ2，τ3＝-τ4，最好取116τ2≤τ4≤12τ2;]]>为了提高位移估计的精度，最好还对计算得到互相关函数进行插值，如抛物线插值等常规插值方法；
5.确定步骤4得到的互相关函数R″(τ)的最大值对应的位置t″1，t″1为利用尺度Tb进行位移估计时数据e1的位移相对于t′1的偏移；6.利用尺度Tb估计得到的位移偏移t″1对尺度Ta估计得到的位移t′1进行调整，得到的小段数据d1的精细调整后的位移t1，即t1＝t′1+t″17.依次从扫描线数据s1(n)中取一小段尺度为Ta即数据个数为Ua的数据d2、d3、…、dN，每段数据的序号依次错开V个采样数据，1≤V≤Ua，直到再错开V个采样数据将超出s1(n)的范围，按步骤2-6相同的方法依次得到各段数据对应的位移t2、t3、…、tN，其中N为小段数据的总数；则位移序列t1、t2、…、tN为第一条扫描线数据s1(n)对应的组织的位移估计；8.利用与步骤1-7相同的方法，依次得到第2、3、…、M条扫描线数据对应的组织的位移估计，其中M为表示探头的扫描线总数，由探头决定。
本发明的原理本发明采用两种尺度的跟踪波段进行位移估计，先用大尺度的跟踪波段Ta进行位移的初估计，然后用小尺度的跟踪波段Tb进行细调。这种采用两种尺度的方法的思想模仿了人眼观察事物的方式先用大尺度粗看，再用小尺度细看，也就是先整体后局部的思想。
本发明利用大的尺度来消除随机噪声的影响，从而确定了组织位移的大致范围，然后利用小的尺度来对组织位移进行精确估计。因为利用小尺度的跟踪波段进行位移估计的时候，是在大尺度跟踪波段的估计结果的基础上进行，可以在一个小的范围内计算互相关函数，所以随机噪声的影响比较小。这种方法打破了尺度固定不变的思路，很好地解决了对组织施加的压缩量较大的时候大尺度和小尺度的矛盾，提高了组织位移估计的精度，又不容易受随机噪声的影响。
从误差分析的角度来看，位移估计的误差可以分成相位误差和伪峰误差。相位误差表现为压缩前、后射频扫描线数据之间的互相关函数的最大值位置在对应组织的实际位移值的附近扰动，误差较小。而伪峰误差表现为该互相关函数的最大值位置已经严重偏离对应组织的实际位移值，此时的组织位移估计已经不包含足够的信息量。当对组织施加的压缩量比较大的时候，尺度越大，越不容易受随机噪声的干扰，越不容易出现伪峰误差，但是由于组织压缩前、后信号重合的部分越小，此时的相位误差较大；而尺度较小时，相位误差得到减小，但是此时的伪峰误差较大。
本发明提出的采用两种尺度的方法，就是先利用大尺度来减小伪峰误差，然后利用小尺度来减小相位误差，从而综合了两种尺度的优点，减小了误差的干扰，提高了组织位移估计的精度。
本发明的特点1采用两种不同尺度的跟踪波段，利用先粗略估计再精细调整、先整理后局部的思路进行位移估计；2先用大尺度的跟踪波段进行位移的粗略估计，然后用小尺度的跟踪波段进行精细调整；3利用精细调整的结果对粗略估计的结果进行修正，从而得到更为精确的组织位移估计。
图1为一般的组织位移估计方法得到的组织应变分布的计算机仿真结果；图2为本发明提出的采用两种尺度的组织位移估计方法的实施例得到的组织应变分布的计算机仿真结果。
具体实施例方式
本发明提出的采用两种尺度的组织位移估计方法结合具体实施例及附图详细说明如下本发明的实施例利用计算机程序和一般的超声散射模型仿真得到一块模拟的组织在压缩前和压缩后的一维射频信号。组织厚度为60mm，其中0-30mm处的弹性模量是30-60mm处的弹性模量的3倍，采用简单的一维模型，因此0-30mm处的应变是30-60mm处的应变的1/3。对组织施加的压缩量为3.33％，即1.98mm；探头中心频率为3.5MHz，-3dB带宽为2.0MHz，因为采用简单的一维模型，所以只有一条扫描线；射频信号的采样频率为20MHz，假设超声波在组织内的传播速度为1540m/s，因此1mm的组织对应1mm/(1540×103mm/s/(20×106Hz)/2≈26个数据，因为组织深度为60mm，所以，每一条扫描线的数据为60×26＝1560个，对组织施加的压缩量以采样数据的个数表示为60×3.33％×26≈52个采样数据。
本实施例采用的大尺度Ta为4mm，搜索范围为0到52；采用的小尺度Tb为1mm(也可采用2mm或在1mm至2mm范围内的任何尺度)，搜索范围为-10到10(也可为-20到-20或满足10≤τ4≤26和τ3＝-τ4的从τ3到τ4的任何范围)。
本实施例的具体步骤如下1.设组织组织压缩前、后的一维射频信号(即扫描线数据)分别为s1(n)和s2(n)；n表示该两条扫描线上的数据序号，1≤n≤1560；2.从该扫描线数据s1(n)中取一小段尺度为Ta的数据d1，Ta＝4mm，其数据个数为U，U＝104，该数据的序号从32到135；在0到52的搜索范围内求该小段数据与扫描线数据s2(n)的互相关函数R(τ)，计算公式如下R′(τ)=Σi=n1n1+U-1s1(i)s2(i-τ)Σi=n1n1+U-1s12(i)·Σi=n1n1+U-1s22(i-τ)(0≤τ≤52)]]>3.确定该互相关函数R′(τ)的最大值对应的位置t′1，t′1就是利用尺度Ta进行位移估计时数据d1在组织压缩后的位移(即s1(n)中的序号从32到135的小段数据d1的在组织压缩后移动到s2(n)中的序号从32-t′1到135-t′1的位置)；
4.取s1(n)中的序号从71到96的一小段数据e1，该数据e1位于数据d1的中央位置；该数据e1的尺度为Tb，Tb＝1mm，数据个数为Ub，Ub＝26；在-10到10的搜索范围内求该小段数据与扫描线数据s2(n)的互相关函数R″(τ)，计算公式如下R′′(τ)=Σi=m1m1+Ub-1s1(i)s2(i-t′1-τ)Σi=m1m1+Ub-1s12(i)·Σi=m1m1+Ub-1s22(i-t′1-τ)(-10≤τ≤10)]]>5.确定步骤4得到的互相关函数R″(τ)的最大值对应的位置t″1，t″1就是利用尺度Tb进行位移估计时数据e1的位移相对于t′1的偏移；6.利用尺度Tb估计得到的位移偏移t″1对尺度Ta估计得到的位移t′1进行调整，得到的小段数据d1的精细调整后的位移t1，即t1＝t′1+t″17.依次从扫描线数据s1(n)中取一小段尺度为Ta即数据个数为Ua的数据d2、d3、…、dN，每段数据的序号依次错开32个采样数据，直到再错开32个采样数据将超出s1(n)的范围，按步骤2-6相同的方法依次得到各段数据对应的位移t2、t3、…、tN，其中N为小段数据的总数，N＝45；则位移序列t1、t2…、tN为扫描线数据s1(n)对应的组织的位移估计。
本实施例与一般方法的位移估计效果比较如下图1和图2为本实施例的结果与一般方法的结果的比较。图1和图2中，横坐标表示组织的不同位置，纵坐标表示应变的大小；图1中的虚线11和图2中的虚线21表示组织的实际应变分布；图1中的12为采用一般方法估计得到的组织应变分布，图2中的22为采用本发明提出的采用两种尺度的方法估计得到的组织应变分布。
由图1和图2可见，本发明提出的采用两种尺度的方法采用两种不同尺度的跟踪波段进行位移估计，结果比采用尺度固定不变的跟踪段进行位移估计的一般方法的结果要好，位移估计的精度更高，特别是在组织应变比较大的区域(30-60mm位置)更为明显。
权利要求
1.一种采用两种尺度的组织位移估计方法，包括以下步骤1)从组织压缩前、后的二维射频信号中分别取出第一条扫描线的数据，设为s1(n)和s2(n)，n表示该两条扫描线上的数据序号，1≤n≤nmax，nmax为n的最大值；2)从该扫描线数据s1(n)中取一小段尺度为Ta的数据d1，其数据个数为Ua，Ua＝round(Ta×U1)，其中，Ta的单位为mm，U1代表1mm的组织对应的数据个数，round(·)代表四舍五入的取整操作，该数据d1的序号从n1到n1+Ua-1，1≤n1≤Ua；在τ1到τ2的搜索范围内求该小段数据与扫描线数据s2(n)的互相关函数R(τ)，计算公式如下R′(τ)=Σi=n1n1+U-1s1(i)s2(i-τ)Σi=n1n1+U-1s12(i)·Σi=n1n1+U-1s22(i-τ)(τ1≤τ≤τ2)]]>其中i为计算过程中表示数据序号的循环变量，τ1为0，τ2为对组织施加的压缩量，以采样数据的个数表示；3)确定步骤2)得到的互相关函数R′(τ)的最大值对应的位置t′1，则t′1为利用尺度Ta进行位移估计时数据d1在组织压缩后的位移(即s1(n)中的序号从n1到n1+Ua-1的小段数据d1的在组织压缩后移动到s2(n)中的序号从n1-t′1到n1+Ua-1-t′1的位置)；4)再从该扫描线数据s1(n)中的数据段d1的中央位置取尺度为Tb、序号从n1+round(Ua-Ub2)]]>到n1+round(Ua-Ub2)+Ub-1]]>的一小段数据e1，该数据个数为Ub，Ub＝round(Tb×U1)，其中，Tb的单位为mm，并且满足Tb＜Ta；在τ3到τ4的搜索范围内求该小段数据与扫描线数据s2(n)的互相关函数R″(τ)，计算公式如下R′′(τ)=Σi=m1m1+Ub-1s1(i)s2(i-t′1-τ)Σi=m1m1+Ub-1s12(i)·Σi=m1m1+Ub-1s22(i-t′1-τ)(τ3≤τ≤τ4)]]>其中τ3和τ4为组织位移精细调整的范围，满足如下关系τ4＜τ2；5)确定步骤4得到的互相关函数R″(τ)的最大值对应的位置t″1，则t″1为利用尺度Tb进行位移估计时数据e1的位移相对于t′1的偏移；6)利用尺度Tb估计得到的位移偏移t″1对尺度Ta估计得到的位移t′1进行调整，得到的小段数据d1的精细调整后的位移t1，即t1＝t′1+t″17)依次从扫描线数据s1(n)中取一小段尺度为Ta即数据个数为Ua的数据d2、d3、…、dN，每段数据的序号依次错开V个采样数据，1≤V≤Ua，直到再错开V个采样数据将超出s1(n)的范围，按步骤2-6相同的方法依次得到各段数据对应的位移t2、t3、…、tN，其中N为小段数据的总数；则位移序列t1、t2、…、tN为第一条扫描线数据s1(n)对应的组织的位移估计；8)利用与步骤1-7相同的方法，依次得到第2、3、…、M条扫描线数据对应的组织的位移估计，其中M为表示探头的扫描线总数。
2.如权利要求1所述的，其特征在于，还进一步包括在所述第2)、4)步骤中对所述计算得到的互相关函数进行插值处理。
3.如权利要求1所述的，其特征在于，所述第4)步骤中的尺度Tb的取值范围为116Ta≤Tb≤12Ta;]]>所述搜索范围的τ4取值为116τ2≤τ4≤12τ2.]]>
全文摘要
本发明属于超声弹性成像技术领域，涉及采用两种尺度的生物组织位移估计方法，该方法包括从组织压缩前二维射频信号中第一条扫描线的数据取出尺度为T
文档编号A61B8/08GK1586409SQ20041005690
公开日2005年3月2日 申请日期2004年8月20日 优先权日2004年8月20日
发明者白净, 罗建文 申请人:清华大学