CN1313054C - 一种多尺度的生物组织位移估计方法 - Google Patents

Abstract

本发明属于超声弹性成像技术领域，涉及一种多尺度的生物组织位移估计方法。该方法包括从组织压缩前的二维射频信号中的第一条扫描线的数据取一小段尺度为T₁的数据d₁，在搜索范围内求该小段数据与组织后的扫描线数据s₂(n)的互相关函数R_1，1(τ)，将数据d₁平均细分成尺度都为T₂的N段数据，分别求各数据的与扫描线数据s₂(n)的互相关函数，再依次对每段数据继续平均细并求各数据的与扫描线数据s₂(n)的互相关函数并进行加权平均，得到其位移t₁，依次从扫描线数据s₁(n)中取一小段尺度为T₁即数据个数为U₁的数据d₂、d₃、…、d_L，即得到第一条扫描线数据s₁(n)对应的组织的位移估计；同样的依次得到各条扫描线数据对应的组织的位移估计。本发明可提高组织位移估计的精度。

Description

一种多尺度的生物组织位移估计方法

技术领域

本发明属于超声弹性成像技术领域，特别涉及生物组织位移估计方法。

背景技术

生物组织弹性模量的变化通常与其病理现象有关。例如，恶性的病理损害，例如乳房硬癌、前列腺癌、甲状腺癌及肝转移等，通常表现为硬的小结。乳房硬癌是乳腺癌的最常见形式，大约占乳腺癌总数的四分之三，由于其基质密度增大而表现为致密的硬块。而其他类型的乳腺癌如导管内癌和乳头状瘤则表现为柔软的组织，良性的乳腺纤维囊性病也很少表现为硬块。

生物组织的弹性模量信息对于疾病的诊断过程具有重要的参考价值。然而，包括X射线成像、超声成像、计算机断层成像(CT)和磁共振成像(MRI)等在内的传统医学成像模态都不能直接提供关于弹性模量这一组织的基本力学属性的信息。1991年，J.Ophir提出超声弹性成像(ultrasound elastography)的方法，对组织的弹性模量分布进行定量估计、成像。目前，超声弹性模量已经成为医学超声成像的一个研究热点，广泛应用于乳房、前列腺、动脉粥样斑块、心肌动力学以及高强度聚焦超声与射频消融引起的损害(lesion)的检测与评估。

超声弹性成像的基本原理为：将超声探头嵌于一块挤压平板中，沿着探头的纵向压缩组织，分别采集组织压缩前、后的射频信号；组织被压缩时，组织内将会产生一个沿压缩方向的应变，如果组织内部弹性模量分布不均匀，组织内的应变分布也会有所差异；弹性模量较大的区域，引起的应变比较小；反之，弹性模量较小的区域，相应的应变比较大。通过一些方法估计出组织内部不同位置的位移，从而计算出组织内部的应变分布情况，用来间接描述组织内部的弹性模量分布，从而描述组织的生理、病理状态。

对于二维超声弹性成像，一般采用线阵的B型超声探头，采集组织压缩前、后的探头每一条扫描线的射频信号，分别进行上面描述的位移估计，从而计算出每一条扫描线对应组织的一维应变分布。最后把所有扫描线对应的一维应变分布按扫描线顺序组成一个二维应变分布，以灰度图或者伪彩图的形式表示，用来间接描述组织内部的弹性模量分布。

一般的超声弹性成像方法包括以下步骤：

1、利用商用B型超声仪器(一般采用线阵探头)得到待测生物组织(一般为人体组织，也可以为动物组织，以下简称组织)压缩前的一幅数字化的二维射频信号(可以采用模拟射频信号输出端接信号放大器，再接高速数据采集卡，得到数字化的二维射频信号；也可以在数字化B型超声仪器上直接得到数字化的二维射频信号)；

2.手持该B型超声仪器的探头或者利用步进电机或者螺旋装置驱动该探头，沿着探头的纵向对该组织施加一个微小的挤压(对组织施加的压缩量一般控制在为1％的数量级)，得到组织压缩后的一幅数字化的二维射频信号；

3.从步骤1和2的得到的组织压缩前、后的二维射频信号中分别取出第一条扫描线的数据，设为s₁(n)和s₂(n)，n表示该两条扫描线上的数据序号，1≤n≤n_max，n的最大值n_max由该B型超声仪器的探查深度、发射的超声波在组织中的传播速度以及射频信号的采样频率决定；

4.从该扫描线数据s₁(n)中取一小段长度为T的数据d₁，其数据个数为U，U＝round(T×U₀)，其中，T的单位为mm，U₀代表1mm的组织对应的数据个数，由发射的超声波在组织中的传播速度以及射频信号的采样频率决定，round(.)代表四舍五入的取整操作，该数据d₁的序号从n₁到n₁+U-1，n₁可在1≤n₁≤U的范围内选择；在τ₁到τ₂确定的搜索范围内求该小段数据与扫描线数据s₂(n)的互相关函数R(τ)，计算公式如下

$R\left(\mathrm{\τ}\right)=\frac{\underset{i=n_1}{\overset{n_1+U-1}{\mathrm{\Σ}}}{s}_1\left(i\right)s_2(i-\mathrm{\τ})}{\sqrt{\underset{i=n_1}{\overset{n_1+U-1}{\mathrm{\Σ}}}{s}_1^2\left(i\right)\·\underset{i=n_1}{\overset{n_1+U-1}{\mathrm{\Σ}}}{s}_2^2(i-\mathrm{\τ})}}\left({\mathrm{\τ}}_1{\≤\mathrm{\τ}\≤\mathrm{\τ}}_2\right)$

其中i为计算过程中表示数据序号的循环变量，τ₁为0，τ₂为对组织施加的压缩量，以采样数据的个数表示；为了提高位移估计的精度，一般还需要对计算得到互相关函数进行插值，如抛物线等常规插值方法；

5.确定步骤4得到的互相关函数R(τ)的最大值对应的位置t₁，t₁就是数据d₁在组织压缩后的位移(即s₁(n)中的序号从n₁到n₁+U-1的小段数据d₁的在组织压缩后移动到s₂(n)中的序号从n₁-t₁到n₁+U-1-t₁的位置)；

6.依次从扫描线数据s₁(n)中取一小段长度为T即数据个数为U的数据d₂、d₃、…、d_L，每段数据的序号依次错开V个采样数据，1≤V≤U，直到再错开V个采样数据将超出s₁(n)的范围，按步骤4、5相同的方法依次得到各段数据对应的位移t₂、t₃、…、t_L，其中L为小段数据的总数；则位移序列t₁、t₂、…、t_L为第一条扫描线数据s₁(n)对应的组织的位移估计；

7.利用与步骤3-6相同的方法，依次得到第2、3、…、M条扫描线数据对应的组织的位移估计，其中M为表示探头的扫描线总数，由探头决定；

8.对第一条扫描线数据s₁(n)对应的组织的位移估计序列t₁、t₂、…、t_L求差分，得到组织第一条扫描线s₁(n)对应组织的应变分布，计算公式如下，

${\mathrm{\ϵ}}_1=\frac{t_2-t_1}{V},$ ${\mathrm{\ϵ}}_2=\frac{t_3-t_2}{V},\·\·\·,$ ${\mathrm{\ϵ}}_{L-1}=\frac{t_L-t_{L-1}}{V}$

其中，ε₁、ε₂、…、ε_L-1分别为d₁、d₂、…、d_L-1对应的组织应变；

9.利用与步骤8相同的方法，依次得到组织第2、3、…、M条扫描线数据对应的组织的应变分布；

10.将步骤9得到的M条扫描线数据对应的应变分布，按照扫描线的顺序组合成一个二维数据，并以灰度图或者伪彩图的形式表示出来，就得到组织的二维应变分布图。

在超声弹性成像中，关键的问题在于对组织的位移分布进行估计，也就是上面描述的方法的步骤3-7。互相关函数的值越大，说明压缩前、后的小段数据吻合得越好，互相关函数的最大值位置代表了压缩前的小段数据在压缩后对应的位置，从而可以求出该小段数据的位移，也就是该小段数据对应的组织的位移。

在超声弹性成像中，选择的压缩前信号中的小段数据用来跟踪对应的小段组织的位移，称为跟踪波段。其长度称为跟踪波段长度，或者尺度。选取合适的尺度在超声弹性成像的位移估计中非常重要。当对组织施加的压缩量比较小(如小于1％)的时候，尺度越大，包含的信息越多，位移估计的精度越高，对随机噪声的干扰也越不敏感。当对组织施加的压缩量比较大(如1-5％)的时候，尺度越大，由于组织被压缩，组织压缩前、后信号重合的部分越小，利用互相关函数最大值位置进行位移估计的精度越低，同时，尺度越大，又越不容易受到随机噪声的干扰；而当尺度较小的时候，组织压缩前、后信号重合的部分较大，位移估计的精度较高，然而此时的位移估计又对随机噪声的干扰比较敏感。

为了增加组织压缩前、后信号重合的程度，一般采用线性插值的方法，把被“压缩”的信号“拉伸”(stretching)成与原信号同长，从而增加压缩前后波形的重合度。在波形拉伸中，拉伸系数与对组织施加的压缩量或应变(以百分比表示)有关。当拉伸系数接近组织的应变时，波形的重合度最大，从而互相关函数最大值最大。但是组织的应变在时延估计前是未知的，拉伸系数的选择只能根据施加的位移量和组织的深度作大致判断。并且由于组织的不均匀性，组织内部应变的分布也不均匀，因此特定的拉伸系数在组织的某些部分可能拉伸不足，在某些部分又拉伸过量，这反过来又引入误差。因此，有人提出一种对不同位置的波段进行自适应拉伸系数的方法，使得拉伸后的波段与压缩前波段的互相关函数最大值最大，然后利用拉伸系数直接估计出该位置组织的应变。但是这种自适应拉伸的方法计算时间较长，在实际应用中可能受到限制。

从实际应用的角度，希望对组织施加的压缩量大一些好，这样可以更能“看清”组织内部弹性模量差异引起的应变分布差异的细节。一般采用多次压缩(multi-compression)的方法，将大位移分解为小位移之和，逐次对组织施压一个小的位移量，计算出相应的应变分布。最后将多幅应变分布叠加，得到大位移时的应变分布。然而，该方法增加了数据采集时间。

发明内容

本发明为了解决对组织施加的压缩量比较大的时候大尺度和小尺度的矛盾，提出一种多尺度的组织位移估计方法；在组织压缩量比较大时，将两种或两种以上不同尺度对应的互相关函数进行加权平均，从而综合了多种不同尺度的信息进行组织位移估计，可对只选择单一的大尺度或小尺度进行组织位移估计的优缺点实现扬长避短，提高组织位移估计的精度。

本发明提出的一种多尺度的组织位移估计方法，包括以下步骤：

1.从组织压缩前、后的二维射频信号中分别取出第一条扫描线的数据，设为s₁(n)和s₂(n)，n表示该两条扫描线上的数据序号，1≤n≤n_max，n的最大值n_max由该B型超声仪器的探查深度、发射的超声波在组织中的传播速度以及射频信号的采样频率决定；

2.从该扫描线数据s₁(n)中取一小段尺度为T₁的数据d₁，其数据个数为U₁，U₁＝round(T₁×U₀)，其中，T₁的单位为mm，U₀代表1mm的组织对应的数据个数，由发射的超声波在组织中的传播速度以及射频信号的采样频率决定，round(.)代表四舍五入的取整操作，该数据d₁的序号从n₁到n₁+U₁-1，n₁可在1≤n₁≤U₁的范围内选择；在τ₁到τ₂确定的搜索范围内求该小段数据与扫描线数据s₂(n)的互相关函数R_1，1(τ)，计算公式如下：

$R_{\mathrm{1,1}}\left(\mathrm{\τ}\right)=\frac{\underset{i=n_1}{\overset{n_1+U_1-1}{\mathrm{\Σ}}}{s}_1\left(i\right)s_2(i-\mathrm{\τ})}{\sqrt{\underset{i=n_1}{\overset{n_1+U_1-1}{\mathrm{\Σ}}}{s}_1^2\left(i\right)\·\underset{i=n_1}{\overset{n_1+U_1-1}{\mathrm{\Σ}}}{s}_2^2(i-\mathrm{\τ})}}\left({\mathrm{\τ}}_1{\≤\mathrm{\τ}\≤\mathrm{\τ}}_2\right)$

其中i为计算过程中表示数据序号的循环变量，τ₁为0，τ₂为对组织施加的压缩量，以采样数据的个数表示(为了提高位移估计的精度，最好对计算得到互相关函数进行插值处理，如抛物线等常规插值方法)；最后得到一个T₁尺度下的互相关函数R_1，1(τ)；

3.将数据d₁平均细分成N段数据d_2，1，d_2，2，…，d_2，N，其尺度都为T₂，T₂＝T₁/N，其中N是尺度T₁细分的段数，即大尺度T₁与小尺度T₂的比值，N≥2，一般可取2≤N≤8；每段小尺度T₂数据的个数为U₂，U₂＝round(T₂×U₀)，第j(1≤j≤N)小段数据d_2，J的序号从n_2，j到n_2，j+U₂-1，n_2，j＝n₁+round((j-1)×T₂×U₀)，在τ₁到τ₂的搜索范围内分别求数据d_2，1，d_2，2，…，d_2，N与扫描线数据s₂(n)的互相关函数R_2，j(τ)(1≤j≤N)，计算公式如下

$R_{2,j}\left(\mathrm{\τ}\right)=\frac{\underset{i=n_{2,j}}{\overset{n_{2,j}+U_2-1}{\mathrm{\Σ}}}{s}_1\left(i\right)s_2(i-\mathrm{\τ})}{\sqrt{\underset{i=n_{2,j}}{\overset{n_{2,j}+U_2-1}{\mathrm{\Σ}}}{s}_1^2\left(i\right)\·\underset{i=n_{2,j}}{\overset{n_{2,j}+U_2-1}{\mathrm{\Σ}}}{s}_2^2(i-\mathrm{\τ})}}\left({\mathrm{\τ}}_1{\≤\mathrm{\τ}\≤\mathrm{\τ}}_2\right)$

最后得到N个T₂尺度下的互相关函数R_2，j(τ)(1≤j≤N)；

4.利用与步骤3相同的方法，分别将步骤3得到的N段数据d_2，1，d_2，2，…，d_2，N中的每段数据继续平均细分成N段数据，这样总共有N²段数据，其尺度为T₃，T₃＝T₂/N；在τ₁到τ₂的搜索范围内分别求这N²段数据与扫描线数据s₂(n)的互相关函数R_3，j(τ)(1≤j≤N²)，得到N²个T₃尺度下的互相关函数R_3，j(τ)(1≤j≤N²)；利用同样的方法，将该N²段小段数据中的每段数据平均细分成N段数据，继续依次再对细分的每段数据平均细分成N段，其尺度依次减小到前一尺度的1/N，即T₄＝T₃/N，…，T_m＝T_m-1/N；m为采用的不同尺度总的个数，2≤m≤6(一般取m＝2或者m＝3即可)，依次计算得到N³，…，N^m-1个在尺度T₄，…T_m的互相关函数；

5.对步骤2-4得到的1，N，…，N^m-1个T₁，T₂，…，T_m尺度下的互相关函数进行加权平均，得到综合了不同尺度信息的复合的互相关函数，公式如下：

$R\left(\mathrm{\τ}\right)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\underset{l=1}{\overset{N^{k-1}}{\mathrm{\Σ}}}{\mathrm{\α}}_{k,l}{R}_{k,l}\left(\mathrm{\τ}\right)$

其中α_k，l为加权系数，α_k，l≥0，并且 $\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\underset{l=1}{\overset{N^{k-1}}{\mathrm{\Σ}}}{\mathrm{\α}}_{k,l}=1,$ 一般可取 ${\mathrm{\α}}_{k,l}=\frac{1}{1+N+\·\·\·+N^{m-1}},$ 即采用算术平均，即各个互相关函数的加权系数都相等。

6.确定步骤5得到的复合的互相关函数R(τ)的最大值对应的位置t₁，则t₁就是数据d₁在组织压缩后的位移(即s₁(n)中的序号从n₁到n₁+U₁-1的小段数据d₁的在组织压缩后移动到s₂(n)中的序号从n₁-t₁到n₁+U₁-1-t₁的位置)：

7.依次从扫描线数据s₁(n)中取一小段尺度为T₁即数据个数为U₁的数据d₂、d₃、…、d_L，每段数据的序号依次错开V个采样数据，1≤V≤U₁，直到再错开V个采样数据将超出s₁(n)的范围，按步骤2-6相同的方法依次得到各段数据对应的位移t₂、t₃、…、t_L，其中L为小段数据的总数；则位移序列t₁、t₂、…、t_L为第一条扫描线数据s₁(n)对应的组织的位移估计；

8.利用与步骤1-7相同的方法，依次得到第2、3、…、M条扫描线数据对应的组织的位移估计，其中M为表示探头的扫描线总数，由探头决定。

本发明的原理：

从误差分析的角度来看，位移估计的误差可以分成相位误差和伪峰误差。相位误差表现为压缩前、后射频扫描线数据之间的互相关函数的最大值位置在对应组织的实际位移值的附近扰动，误差较小。而伪峰误差表现为该互相关函数的最大值位置已经严重偏离对应组织的实际位移值，此时的组织位移估计已经不包含足够的信息量。当对组织施加的压缩量比较大的时候，尺度越大，越不容易受随机噪声的干扰，越不容易出现伪峰误差，但是由于组织压缩前、后信号重合的部分越小，此时的相位误差较大；而尺度较小时，相位误差得到减小，但是此时的伪峰误差较大。

本发明提出的多尺度的方法，综合了多种不同尺度的信息，利用大尺度来减小伪峰误差，利用小尺度来减小相位误差，并且采用加权平均的方法，即使某些跟踪波段计算得到的互相关函数的最大值位置误差较大，对最后的复合的互相关函数的影响也不大，从而减小了误差的干扰，提高了组织位移估计的精度。

本发明的特点：

1)采用两种或两种以上不同尺度的跟踪波段，分别计算它们与压缩后的扫描线数据的互相关函数；

2)首先采用大尺度的跟踪波段，计算它与压缩后的扫描线数据的互相关函数；

3)接着将大尺度的跟踪波段平均细分成分成若干段小尺度的跟踪波段，分别计算它们与压缩后的扫描线数据的互相关函数；

4)如果需要，继续对小尺度的跟踪波段平均细分成若干段更小尺度的跟踪波段，分别计算它们与压缩后的扫描线数据的互相关函数；

5)然后将用不同尺度计算得到的互相关函数进行加权平均，得到综合了多种不同尺度信息的复合的互相关函数；

6)最后复合的互相关函数的最大值位置进行组织位移估计。

附图说明

图1为一般的组织位移估计方法得到的组织应变分布的计算机仿真结果；

图2为本发明提出的多尺度的组织位移估计方法的实施例得到的组织应变分布的计算机仿真结果。

具体实施方式

本发明提出的多尺度的组织位移估计方法结合具体实施例及附图详细说明如下：

本发明的实施例利用计算机程序和一般的超声散射模型仿真得到一块模拟的组织在压缩前和压缩后的一维射频信号。组织厚度为60mm，其中0-30mm处的弹性模量是30-60mm处的弹性模量的3倍，采用简单的一维模型，因此0-30mm处的应变是30-60mm处的应变的1/3。对组织施加的压缩量为4.8％，即2.88mm；探头中心频率为3.5MHz，-3dB带宽为2.0MHz，因为采用简单的-维模型，所以只有-条扫描线；射频信号的采样频率为20MHz，假设超声波在组织内的传播速度为1540m/s，因此1mm的组织对应1mm/(1540×10³mm/s/(20×10⁶Hz)/2≈26个数据，因为组织深度为60mm，所以，每一条扫描线的数据为60×26＝1560个，对组织施加的压缩量以采样数据的个数表示为60×4.8％×26≈75个采样数据。

本实施例采用三种不同的尺度，分别为T₁＝4mm，T₂＝2mm，T₃＝1mm，跟踪波段细分的段数N＝2；搜索范围为0到75。

本实施例的具体步骤如下：

1.设组织组织压缩前、后的一维射频信号(即扫描线数据)分别为s₁(n)和s₂(n)；n表示该两条扫描线上的数据序号，1≤n≤1560；

2.从该扫描线数据s₁(n)中取一小段尺度为T₁的数据d₁，T₁＝4mm，其数据个数为U₁，U₁＝104，该数据的序号从32到135；在0到75的搜索范围内求该小段数据与扫描线数据s₂(n)的互相关函数R_1，1(τ)，计算公式如下

$R_{\mathrm{1,1}}\left(\mathrm{\τ}\right)=\frac{\underset{i=n_1}{\overset{n_1+U_1-1}{\mathrm{\Σ}}}{s}_1\left(i\right)s_2(i-\mathrm{\τ})}{\sqrt{\underset{i=n_1}{\overset{n_1+U_1-1}{\mathrm{\Σ}}}{s}_1^2\left(i\right)\·\underset{i=n_1}{\overset{n_1+U_1-1}{\mathrm{\Σ}}}{s}_2^2(i-\mathrm{\τ})}}(0\≤\mathrm{\τ}\≤75)$

3.将数据d₁平均细分成2段数据d_2，1，d_2，2，其尺度都为T₂，T₂＝2mm，数据个数为U₂，U₂＝52，数据d_2，1，d_2，2的序号分别从32到83和从84到135，即n_2，1＝32，n_2，2＝84；在0到75的搜索范围内分别求数据d_2，1，d_2，2与扫描线数据s₂(n)的互相关函数R_2，1(τ)，R_2，2(τ)，计算公式分别如下

$R_{2,1}\left(\mathrm{\τ}\right)=\frac{\underset{i=n_{\mathrm{2,1}}}{\overset{n_{\mathrm{2,1}}+U_2-1}{\mathrm{\Σ}}}{s}_1\left(i\right)s_2(i-\mathrm{\τ})}{\sqrt{\underset{i=n_{\mathrm{2,1}}}{\overset{n_{\mathrm{2,1}}+U_2-1}{\mathrm{\Σ}}}{s}_1^2\left(i\right)\·\underset{i=n_{\mathrm{2,1}}}{\overset{n_{\mathrm{2,1}}+U_2-1}{\mathrm{\Σ}}}{s}_2^2(i-\mathrm{\τ})}}(0\≤\mathrm{\τ}\≤75)$

$R_{\mathrm{2,2}}\left(\mathrm{\τ}\right)=\frac{\underset{i=n_{\mathrm{2,2}}}{\overset{n_{\mathrm{2,2}}+U_2-1}{\mathrm{\Σ}}}{s}_1\left(i\right)s_2(i-\mathrm{\τ})}{\sqrt{\underset{i=n_{\mathrm{2,2}}}{\overset{n_{\mathrm{2,2}}+U_2-1}{\mathrm{\Σ}}}{s}_1^2\left(i\right)\·\underset{i=n_{\mathrm{2,2}}}{\overset{n_{\mathrm{2,2}}+U_2-1}{\mathrm{\Σ}}}{s}_2^2(i-\mathrm{\τ})}}(0\≤\mathrm{\τ}\≤75)$

4.利用与步骤3相同的方法，分别将步骤3得到的2段数据d_2，1，d_2，2中的每段数据继续平均细分成2段数据，这样总共有4段数据，其尺度为T₃，T₃＝1mm；利用与步骤3相同的方法，在0到75的搜索范围内分别求这4段数据与扫描线数据s₂(n)的互相关函数R_3，1(τ)，R_3，2(τ)，R_3，3(τ)，R_3，4(τ)；

5.对步骤2-4得到的1，2，4个T₁，T₂，T₃尺度下的互相关函数进行加权平均，得到综合了三种不同尺度信息的复合的互相关函数，公式如下

$R\left(\mathrm{\τ}\right)=\frac{1}{7}(R_{1,1}\left(\mathrm{\τ}\right)+R_{\mathrm{2,1}}\left(\mathrm{\τ}\right)+R_{\mathrm{2,2}}\left(\mathrm{\τ}\right)+R_{\mathrm{3,1}}\left(\mathrm{\τ}\right)+R_{\mathrm{3,2}}\left(\mathrm{\τ}\right)+R_{\mathrm{3,3}}\left(\mathrm{\τ}\right)+R_{\mathrm{3,4}}\left(\mathrm{\τ}\right))$

即采用算术平均，各个互相关函数的加权系数都相等，等于1/7；

6.确定步骤5得到的复合的互相关函数R(τ)的最大值对应的位置t₁，t₁就是数据d₁在组织压缩后的位移(即s₁(n)中的序号从n₁到n₁+U₁-1的小段数据d₁的在组织压缩后移动到s₂(n)中的序号从n₁-t₁到n₁+U₁-1-t₁的位置)；

7.依次从扫描线数据s₁(n)中取一小段尺度为T₁即数据个数为U₁的数据d₂、d₃、…、d_L，每段数据的序号依次错开32个采样数据，直到再错开32个采样数据将超出s₁(n)的范围，按步骤2-6相同的方法依次得到各段数据对应的位移t₂、t₃、…、t_L，其中L为小段数据的总数，L＝45；则位移序列t₁、t₂、…、t_L为扫描线数据s₁(n)对应的组织的位移估计。

本实施例与一般方法的位移估计效果比较如下：

图1和图2为本实施例的结果与一般方法的结果的比较。图1和图2中，横坐标表示组织的不同位置，纵坐标表示应变的大小；图1中的虚线11和图2中的虚线21表示组织的实际应变分布；图1中的12为采用一般方法估计得到的组织应变分布，图2中的22为采用本发明提出的多尺度的方法估计得到的组织应变分布。

由图1和图2可见，对组织施加的压缩量为4.8％的时候，一般方法已不能准确地估计出组织的位移分布，计算得到的应变分布偏离实际值太远；而本发明提出的多尺度的方法仍然能够较为准确地估计出组织的位移分布，计算得到的应变分布接近实际值。从而说明，本发明对于较大的组织压缩量，能够提高组织位移估计的精度。

Claims (4)

1、一种多尺度的组织位移估计方法，包括以下步骤：

1)从组织压缩前、后的二维射频信号中分别取出第一条扫描线的数据，设为s₁(n)和s₂(n)，n表示该两条扫描线上的数据序号，1≤n≤n_max，n_max为n的最大值；

2)从该扫描线数据s₁(n)中取一小段尺度为T₁的数据d₁，其数据个数为U₁，U₁＝round(T₁×U₀)，其中，T₁的单位为mm，U₀代表1mm的组织对应的数据个数，round(·)代表四舍五入的取整操作，该数据d₁的序号从n₁到n₁+U₁-1，1≤n₁≤U₁；在τ₁到τ₂的搜索范围内求该小段数据与扫描线数据s₂(n)的互相关函数R_1，1(τ)，计算公式如下：

$R_{\mathrm{1,1}}\left(\mathrm{\τ}\right)=\frac{\underset{i=n_1}{\overset{n_1+U_1-1}{\mathrm{\Σ}}}{s}_1\left(i\right)s_2(i-\mathrm{\τ})}{\sqrt{\underset{i=n_1}{\overset{n_1+U_1-1}{\mathrm{\Σ}}}{s}_1^2\left(i\right)\·\underset{i=n_1}{\overset{n_1+U_1-1}{\mathrm{\Σ}}}{s}_2^2(i-\mathrm{\τ})}}$ (τ₁≤τ≤τ₂)

其中i为计算过程中表示数据序号的循环变量，τ₁为0，τ₂为对组织施加的压缩量，以采样数据的个数表示；

3)将数据d₁平均细分成尺度都为T₂的N段数据d_2，1，d_2，2，…，d_2，N，T₂＝T₁/N，N≥2；每段数据个数为U₂，U₂＝round(T₂×U₀)，第j(1≤j ≤N)小段数据d_2，j的序号从n_2，j到n_2，j+U₂-1，n_2，j＝n₁+round((j-1)×T₂×U₀)，在τ₁到τ₂的搜索范围内分别求数据d_2，1，d_2，2，…，d_2，N与扫描线数据s₂(n)的互相关函数R_2，j(τ)(1≤j≤N)，计算公式如

$R_{2,j}\left(\mathrm{\τ}\right)=\frac{\underset{i=n_{2,j}}{\overset{n_{2,j}+U_2-1}{\mathrm{\Σ}}}{s}_1\left(i\right)s_2(i-\mathrm{\τ})}{\sqrt{\underset{i=n_{2,j}}{\overset{n_{2,j}+U_2-1}{\mathrm{\Σ}}}{s}_1^2\left(i\right)\·\underset{i=n_{2,j}}{\overset{n_{2,j}+U_2-1}{\mathrm{\Σ}}}{s}_2^2(i-\mathrm{\τ})}}$ (τ₁≤τ≤τ₂)得到N个T₂尺度下的互相关函数R_2，j(τ)(1≤j≤N)；

4.)利用与步骤3)相同的方法，分别将步骤3)得到的N段数据d_2，1，d_2，2，…，d_2，N中的每段数据继续平均细分成N段数据，共有N²段数据，其尺度为T₃，T₃＝T₂/N；在τ₁到τ₂的搜索范围内分别求这N²段数据与扫描线数据s₂(n)的互相关函数R_3，j(τ)(1≤j≤N²)，得到N²个T₃尺度下的互相关函数R_3，j(τ)(1≤j≤N²)；再将该N²段小段数据中的每段数据平均细分成尺度为T₄的N段数据，继续依次再对细分的每段数据平均细分成N段，其尺度依次减小到前一尺度的1/N，即T₄＝T₃/N，…，T_m＝T_m-1/N；m为采用的不同尺度的总个数，2≤m≤6，依次计算得到N³，…，N^m-1个在尺度T₄，…T_m的互相关函数；

5)对步骤2)-4)得到的1，N，…，N^m-1个T₁，T₂，…，T_m尺度下的互相关函数进行加权平均，得到综合了不同尺度信息的复合的互相关函数，公式如下：

$R\left(\mathrm{\τ}\right)=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\underset{l=1}{\overset{N^{k-1}}{\mathrm{\Σ}}}{\mathrm{\α}}_{k,l}{R}_{k,l}\left(\mathrm{\τ}\right)$

其中R_k，l(τ)为尺度T_k下的第l个窗口的互相关函数，α_k，l为加权系数，并且

${\mathrm{\α}}_{k,l}\≥0,\underset{k=1}{\overset{m}{\mathrm{\Σ}}}\underset{l=1}{\overset{N^{k-1}}{\mathrm{\Σ}}}{\mathrm{\α}}_{k,l}=1;$

6)确定步骤5)得到的复合的互相关函数R(τ)的最大值对应的位移t₁；

7)依次从扫描线数据s₁(n)中取一小段尺度为T₁即数据个数为U₁的数据段d₂、d₃、…、d_L，每段数据的序号依次错开V个采样数据，1≤V≤U₁，直到再错开V个采样数据将超出s₁(n)的范围，按步骤2)-6)相同的方法依次得到各段数据对应的位移t₂、t₃、…、t_L，其中L为小段数据的总数；则位移序列t₁、t₂、…、t_L为第一条扫描线数据s₁(n)对应的组织的位移估计；

8)采用与步骤1)-7)相同的方法，依次得到第2、3、…、M条扫描线数据对应的组织的位移估计，其中M为表示探头的扫描线总数。

2、如权利要求1所述的多尺度的组织位移估计方法，其特征在于，所述步骤2)中还包括对得到的互相关函数R_1，1(τ)再进行插值处理。

3、如权利要求1所述的多尺度的组织位移估计方法，其特征在于，所述步骤3)中的大尺度T₁与小尺度T₂的比值N取值为2≤N≤8。

4、如权利要求1所述的多尺度的组织位移估计方法，其特征在于，所述步骤5)中复合的互相关函数的加权系数 ${\mathrm{\α}}_{k,l}=\frac{1}{1+N+\·\·\·+N^{m-1}}.$