CN1586411A - 一种二维综合互相关的生物组织位移估计方法 - Google Patents

Abstract

本发明涉及一种二维综合互相关的生物组织位移估计方法，属于超声弹性成像技术领域，本方法包括从组织压缩前、后的二维射频信号中分别取出第m+1条扫描线的数据，m为二维综合互相关的综合系数，从压缩前的扫描线数据中取多个小段长度为T的数据，求各小段数据与压缩后扫描线数据互相关函数，对各相关函数进行加权平均，得到二维综合互相关函数，该各互相关函数的最大值对应的位置是各段数据对应的位移，同样得到各条扫描线数据对应的组织的位移估计；本发明通过综合考虑相邻的多条扫描线数据的信息，来减小组织横向位移引入的误差，较好地实现抑制由组织横向位移引入的组织纵向位移估计及应变估计误差，从而提高组织纵向位移估计的精度。

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{\τ})}}({\mathrm{\τ}}_1\≤\mathrm{\τ}\≤{\mathrm{\τ}}_2)$

其中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_N，每段数据的序号依次错开V个采样数据(如V＝round(0.4×T×U₀)、V＝round(0.5×T×U₀)等)，直到再错开V个采样数据将超出s₁(n)的范围，按步骤4、5相同的方法依次得到各段数据对应的位移t₂、t₃、…、t_N，其中N为小段数据的总数；则位移序列t₁、t₂、…、t_N为第一条扫描线数据s₁(n)对应的组织的位移估计；

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

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

${\mathrm{\ϵ}}_1=\frac{t_2-t_1}{V},{\mathrm{\ϵ}}_2=\frac{t_3-t_2}{V},\…{,\mathrm{\ϵ}}_{N-1}=\frac{t_N-t_{N-1}}{V}$

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

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

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

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

超声弹性成像中，对组织施加一个小的压缩量，利用互相关分析估计的组织位移是纵向位移，即沿压缩方向的位移。但是，对组织施加一个小的压缩量的时候，组织的运动是很复杂的，受到组织内部弹性模量分布、组织的几何形状、边界条件等因素的影响。组织不仅沿着压缩方向(即线阵探头的纵向)有一个压缩，沿着垂直于压缩方向的方向(包括线阵探头的横向以及垂直于探头扫描平面的方向)也会膨胀。研究表明，沿着垂直于压缩方向的位移会引起压缩前、后信号的互相关函数的幅度降低，也就是说，压缩前、后信号的相似性降低。而超声弹性成像正是利用压缩前、后信号的相似性来跟踪组织位移的，所以，沿着垂直于压缩方向的位移会使得组织位移估计的精度降低。并且，垂直于探头扫描平面方向的位移的影响比横向位移的影响小，所以，减少沿着垂直于压缩方向的位移的影响，主要是减小垂直于探头扫描平面的方向的位移影响。

上述的位移估计方法中是对每一条扫描线的信号进行互相关分析，得到组织的位移估计，组织横向位移引入的影响较大，组织纵向位移估计的精度较差，应变分布估计的信噪比较低。

发明内容

本发明的目的是为克服已有技术的不足之处，提出了一种二维综合互相关的组织位移估计方法，能够充分利用从商用B型超声仪器得到的射频信号的二维信息，通过综合考虑相邻的多条扫描线数据的信息，来减小组织横向位移引入的误差，较好地实现抑制由组织横向位移引入的组织纵向位移估计及应变估计误差，从而提高组织纵向位移估计的精度。

本发明提出的一种二维综合互相关的组织位移估计方法，包括以下步骤：

1.从组织压缩前、后的二维射频信号中分别取出第m+1条扫描线的数据，设为s_1，m+1(n)和s_2，m+1(n)，n表示该两条扫描线上的数据序号，1≤n≤n_max，n的最大值n_max由该B型超声仪器的探查深度、发射的超声波在组织中的传播速度以及射频信号的采样频率决定，m为二维综合互相关的综合系数，m为整数(其取值范围最好为1≤m≤5)；

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

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

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

3.依次从组织压缩前、后的二维射频信号中取出第m+1条扫描线的前m条和后m条扫描线的数据，即第1到m条扫描线和第m+2到2m+1条扫描线的数据，压缩前、后的扫描线数据分别设为s_1，1(n)、s_1，2(n)、…s_1，m(n)、s_1，m+2(n)、s_1，m+3(n)、…、s_1，2m+1(n)和s_2，1(n)、s_2，2(n)、…、s_2，m(n)、s_2，m+2(n)、s_2，m+3(n)、…、s_2，2m+1(n)，依次取与d_1，m+1同样长度(T)和同样序号(序号从n₁到n₁+U-1)的一小段数据d_1，1(n)、d_1，2(n)、…、d_1，m(n)、d_1，m+2(n)、d_1，m+3(n)、…、d_1，2m+1(n)，利用与步骤2相同的计算方法，在τ₁到τ₂确定的搜索范围内分别求小段数据d_1，1(n)、d_1，2(n)、…、d_1，m(n)、d_1，m+2(n)、d_1，m+3(n)、…d_1，2m+1(n)与对应的压缩后的扫描线数据的互相关函数，即求d_1，1(n)与s_2，1(n)的互相关函数R₁(τ)，求d_1，2(n)与s_2，2(n)的互相关函数R₂(τ)，…，求d_1，m(n)与s_2，m(n)的互相关函数R_m(τ)，求d_1，m+2(n)与s_2，m+2(n)的互相关函数R_m+2(τ)，求d_1，m+3(n)与s_2，m+3(n)的互相关函数R_m+3(τ)，…，求d_1，2m+1(n)与s_2，2m+1(n)的互相关函数R_2m+1(τ)；

4.对步骤2-3得到的互相关函数R₁(τ)、R₂(τ)、…、R_m(τ)、R_m+1(τ)、R_m+2(τ)、…、R_2m+1(τ)进行加权平均，得到包含了射频信号二维信息的复合的互相关函数，即二维综合互相关函数R_m+1′(τ)，计算公式如下

${R_{m+1}}^{\′}\left(\mathrm{\τ}\right)=\underset{k=-m}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_k{R}_{m+1+k}\left(\mathrm{\τ}\right)$

其中，α_k代表各扫描线数据对应互相关函数的权重，α_k＞0(-m≤k≤m)，并且所有扫描线数据对应互相关函数的权重之和为1，即 $\underset{k=-m}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_k=1;$ (为了提高位移估计的精度，一般还需要对该二维综合互相关函数R_m+1′(τ)进行插值，如抛物线插值)；

5.确定该互相关函数R_m+1′(τ)的最大值对应的位置t₁，t₁就是数据d_1，m+1在组织压缩后的位移，即s_1，m+1(n)中序号从n₁到n₁+U-1的小段数据d_1，m+1在组织压缩后移动到s_2，m+1(n)中的序号从n₁-t₁到n₁+U-1-t₁的位置；

6.依次从该扫描线数据s_1，m+1(n)中取一小段长度为T即数据个数为U的数据d_2，m+1、d_3，m+1、…、d_N，m+1，每段数据的序号依次错开V个采样数据(如V＝round(0.4×T(1)×U₀)、V＝round(0.5×T(1)×U₀)等)，直到再错开V个采样数据将超出s_1，m+1(n)的范围，按步骤2-5相同的方法依次得到各段数据对应的位移t₂、t₃、…、t_N，其中N为小段数据的总数；则位移序列t₁、t₂、…、t_N为第m+1条扫描线数据s_1，m+1(n)对应的组织的位移估计；

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

本发明的原理：

本发明方法是在计算第j(m+1≤j≤M-m)条扫描线的压缩前数据中的某个跟踪波段的位移时，不仅计算该跟踪波段与第j扫描线的压缩后数据的互相关函数，还计算和第j条扫描线相邻的2m条扫描线的压缩前数据中的相同长度、相同位置的跟踪波段和对应的压缩后数据的互相关函数，即计算了从第j-m到第j+m的2m+1条扫描线的压缩前数据中的相同长度、相同位置的跟踪波段和对应的压缩后数据的互相关函数，将计算得到的2m+1个互相关函数进行加权平均，得到包含了射频信号二维信息的复合的互相关函数，即二维综合互相关函数。

如图1所示，在计算第j(m+1≤j≤M-m)条扫描线的压缩前数据中的某个跟踪波段的位移时，采用从第j-m到第j+m的2m+1条扫描线11，计算每条扫描线的压缩前数据中的相同长度、相同位置的跟踪波段和对应的压缩后数据12之间的互相关函数13，并将计算得到的2m+1个互相关函数13进行加权平均14，得到二维综合互相关函数15，利用二维综合互相关函数15进行组织纵向位移的估计。

由于受到各种因素的影响，压缩前、后的信号必然包含有信号成分和噪声成分。这样，压缩前、后射频信号之间的互相关函数也包含了信号成分和噪声成分。

本发明提出的二维综合互相关的方法的基本思路就是对多条扫描线对应的压缩前、后射频信号之间的互相关函数进行加权平均，得到包含了射频信号二维信息的复合的互相关函数，从而突出互相关函数中的信号成分，而互相关函数中由组织横向位移影响引入的噪声成分抵消了一部分。这样，组织横向位移引入的噪声减少，加权平均后利用互相关函数来估计组织的位移精度更高。同时，其他噪声也得到一定程度的抵消，从而进一步提高了组织纵向位移估计的精度。

如果相邻的2m+1条扫描线的压缩前数据中的相同长度、相同位置的跟踪波段的位移相同，在没有噪声的理想情况下，他们与对应的压缩后数据的互相关函数的最大值位置相同，而在含有噪声例如组织横向位移引入的噪声的情况下，这些互相关函数的最大值位置存在一定的偏差。本发明提出的二维综合互相关的组织位移估计方法，就是将这些互相关函数进行加权平均，从而减小由于噪声例如组织横向位移引入的噪声干扰导致的偏差，从而提高组织纵向位移估计的精度。

如果相邻的2m+1条扫描线的压缩前数据中的相同长度、相同位置的跟踪波段的位移存在一些差异，则本发明提出的二维综合互相关的组织位移估计方法的得到的综合互相关函数的为相邻的2m+1条扫描线的压缩前数据中的相同长度、相同位置的跟踪波段与对应的压缩后数据之间的互相关函数的平滑化结果，对组织纵向位移估计也取到一定的平滑作用，从而消除噪声的干扰，提高组织纵向位移估计的精度。同时，该方法也导致超声弹性成像的横向分辨率的减小。一般情况下，采用商用B型超声仪器的探头(一般为线阵探头)的扫描线之间间距较小(如0.1mm的数量级)，当采用不太多的若干条扫描线(m≤5表示扫描线数目不超过11条)的压缩前数据中的相同长度、相同位置的跟踪波段与对应的压缩后数据之间的互相关函数进行加权平均后，对成像的横向分辨率的影响并不明显，而又能提高组织纵向位移估计的精度。

影响二维综合互相关方法的参数包括综合系数m与各扫描线信号互相关的权重α_k(-m≤k≤m)。在实际的应用中，如果m的取值过大，将会减小超声弹性成像的横向分辨率；取值过小，则二维综合互相关方法的效果不明显。一般m的取值为1，或2。权重α_k(-m≤k≤m)的取值有多种方法，简单的，可以令各扫描线数据对应的互相关函数的权重相等，均为1/(2m+1)，也可以为其他的满足α_k＞0(-m≤k≤m)和 $\underset{k=-m}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_k=1$ 的加权方式。

如果m＝0，或者α_-m＝α_-m+1＝…＝α_-1＝α₁＝α₂＝…＝α_m＝0，α₀＝1，则

                              R_j+1′(τ)＝R_j+1(τ)

表示这时候的二维综合互相关函数退化成为一般方法采用的互相关函数，由于该方法没有充分利用射频信号的二维信息，容易受到组织横向位移的影响，从而影响组织纵向位移估计的精度。

因此，本发明提出的二维综合互相关的组织位移估计方法，要求满足m≥1和α_k＞0(-m≤k≤m)。

本发明提出的二维综合互相关的组织位移估计方法，提高了第m+1条到第M-m条扫描线对应的组织位移估计的精度。对于第1到第m条扫描线以及第M-m+1到第M条扫描线，可以不进行组织位移估计，这将减小超声弹性成像的横向范围(即扫描线数目)，但是由于数目较小(m≤5表示扫描线数目不超过11条)，在最后的二维应变分布图上并不明显；也可以采用一般的方法进行组织位移估计，从而不减小超声弹性成像的横向范围；虽然这些扫描线位于探头的边缘，组织横向位移的影响较大，但是由于数目较小，在最后的二维应变分布图上也不明显。

本发明的特点：

1)综合相邻的多条扫描线的压缩前数据中的相同长度、相同位置的跟踪波段与对应的压缩后数据之间的多个互相关函数进行组织纵向位移估计；

2)将上述的多个互相关函数进行加权平均，得到包含了射频信号二维信息的复合的互相关函数，即二维综合互相关函数；

3)利用上述的二维综合互相关函数进行组织纵向位移估计，从而减小组织横向位移的影响。

附图说明

图1为本发明提出的二维综合互相关的组织位移估计的的示意图。

图2为本实施例的计算机仿真的组织模型；

图3为利用有限元分析计算得到的组织应变分布的理想结果；

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

图5为综合系数m＝1的二维综合互相关的组织位移估计方法得到的组织应变分布的计算机仿真结果；

图6为综合系数m＝2的二维综合互相关的组织位移估计方法得到的组织应变分布的计算机仿真结果。

具体实施方式

本发明提出的二维综合互相关的组织位移估计方法结合具体实施例及附图详细

说明如下：

实施例1利用计算机程序和一般的超声散射模型仿真得到一块模拟的组织在压缩前和压缩后的二维射频信号。模拟的组织结构如图2所示，组织21大小60×60mm²，组织内分布有3个弹性模量较大的圆形异物22、23、24，它们的弹性模量是组织1的2倍，它们的直径均为5mm；组织压缩比为1％，即压缩量为0.6mm；探头中心频率为3.5MHz，-3dB带宽为2.0MHz，探头扫描线宽度和间隔分别为2mm和0.4mm，探头宽度和组织宽度一致，也是60mm，因此总共有151条扫描线，即M＝151，并且探头中央部位和探头边缘部位分别对应组织中央部位和组织边缘部位；射频信号的采样频率为20MHz，假设超声波在组织内的传播速度为1540m/s，因此1mm的组织长度对应1mm/(1540×10³mm/s×1/20×10⁶/2)Hz≈26个数据，因为组织深度为60mm，所以，每一条扫描线的数据为60×26＝1560个，对组织施加的压缩量以采样数据的个数表示为60×1％×26≈16个采样数据。

图3表示利用美国MSC公司的MARC软件进行有限元分析，计算得到的该实施例采用的组织模型的理想的应变分布。横向和纵向分别表示组织的横向位置和纵向位置(即组织深度)，灰度表示计算出来的理想应变的大小，灰度值越大(即颜色越亮或越白)，表示应变越大，灰度值越小(即颜色越暗或越黑)，表示应变越小，31为灰度值与应变大小的对照关系。图3中，较暗的区域(即32-34)与弹性模量较大的组织层(即图2中的22-24)对应，说明弹性模量较大的区域应变较小。

本实施例中，m的取值为1，令各扫描线数据对应的互相关函数的权重α_k(-m≤k≤m)相等，均为1/(2m+1)，即为1/3；跟踪波段的长度T取3mm，即跟踪波段的数据个数为26×3＝78个，相邻的跟踪波段错开26个采样数据，即V＝26。

利用本发明提出的二维综合互相关的组织位移估计方法，估计出第m+1条到第M-m条扫描线对应的组织位移，从而得到对应的应变分布；而对于第1到第m条扫描线以及第M-m+1到第M条扫描线，不进行组织位移估计，减小超声弹性成像的横向范围(即扫描线数目)，但是由于数目较小(m≤5表示扫描线数目不超过11条)，影响不大。

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

1.从组织压缩前、后的二维射频信号(计算机仿真得到)中分别取出第一条扫描线的数据，设为s_1，m+1(n)和s_2，m+1(n)，1≤n≤1560，m为二维综合互相关方法的综合系数，m＝1；

2.从该扫描线数据s_1，m+1(n)中取一小段长度为T的数据d_1，m+1，T＝3mm，其数据个数为U，U＝78，该数据的序号从13到90；在0到16的搜索范围内求该小段数据与扫描线数据s_2，m+1(n)的互相关函数R_m+1(τ)，计算公式如下

$R_{m+1}\left(\mathrm{\τ}\right)=\frac{\underset{i=n_1}{\overset{n_1+U-1}{\mathrm{\Σ}}}{s}_{1,m+1}\left(i\right)s_{2,m+1}(i-\mathrm{\τ})}{\sqrt{\underset{i=n_1}{\overset{n_1+U-1}{\mathrm{\Σ}}}{s}_{1,m+1}^2\left(i\right)\·\underset{i=n_1}{\overset{n_1+U-1}{\mathrm{\Σ}}}{s}_{2,m+1}^2(i-\mathrm{\τ})}}(0\≤\mathrm{\τ}\≤16)$

其中i为计算过程中表示数据序号的循环变量；

3.依次从组织压缩前、后的二维射频信号中取出第m+1条扫描线的前m条和后m条扫描线的数据，即第1到m条扫描线和第m+2到2m+1条扫描线的数据，压缩前、后的扫描线数据分别设为s_1，1(n)、s_1，2(n)、…s_1，m(n)、s_1，m+2(n)、s_1，m+3(n)、…、s_1，2m+1(n)和s_2，1(n)、s_2，2(n)、…、s_2，m(n)、s_2，m+2(n)、s_2，m+3(n)、…、s_2，2m+1(n)，依次取与d_1，m+1同样长度(3mm)和同样序号(序号从13到90)的一小段数据d_1，1(n)、d_1，2(n)、…、d_1，m(n)、d_1，m+2(n)、d_1，m+3(n)、…、d_1，m+1(n)，利用与步骤2相同的计算方法，在0到16的的搜索范围内分别求小段数据d_1，1(n)、d_1，2(n)、…、d_1，m(n)、d_1，m+2(n)、d_1，m+3(n)、…d_1，2m+1(n)与对应的压缩后的扫描线数据的互相关函数，即求d_1，1(n)与s_2，1(n)的互相关函数R₁(τ)，求d_1，2(n)与s_2，2(n)的互相关函数R₂(τ)，…，求d_1，m(n)与s_2，m(n)的互相关函数R_m(τ)，求d_1，m+2(n)与s_2，m+2(n)的互相关函数R_m+2(τ)，求d_1，m+3(n)与s_2，m+3(n)的互相关函数R_m+3(τ)，…，求d_1，2m+1(n)与s_2，2m+1(n)的互相关函数R_2m+1(τ)；

4.对步骤2-3得到的互相关函数R₁(τ)、R₂(τ)、…、R_m(τ)、R_m+1(τ)、R_m+2(τ)、…、R_2m+1(τ)进行加权平均，得到包含了射频信号二维信息的复合的互相关函数，即二维综合互相关函数R_m+1′(τ)，计算公式如下

${R_{m+1}}^{\′}\left(\mathrm{\τ}\right)=\underset{k=-m}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_k{R}_{m+1+k}\left(\mathrm{\τ}\right)$

其中，α_i代表各扫描线数据对应互相关函数的权重，α_k＝1/(2m+1)(-m≤k≤m)；(为了提高位移估计的精度，一般还需要对该二维综合互相关函数R_m+1′(τ)进行插值，如抛物线插值)；

5.确定该互相关函数R_m+1′(τ)的最大值对应的位置t₁，t₁就是数据d_1，m+1在组织压缩后的位移(即s_1，m+1(n)中序号从13到90的小段数据d_1，m+1在组织压缩后移动到s_2，m+1(n)中的序号从13-t₁到90-t₁的位置)；

6.依次从该扫描线数据s_1，m+1(n)中取一小段长度为3mm即数据个数为78的数据d_2，m+1、d_3，m+1、…、d_N，m+1，每段数据的序号依次错开26个采样数据(即V＝26)，直到再错开26个采样数据将超出s_1，m+1(n)的范围，按步骤2-5相同的方法依次得到各段数据对应的位移t₂、t₃、…、t_N，其中N为小段数据的总数，N＝60；则位移序列t₁、t₂、…、t_N为第m+1条扫描线数据s_1，m+1(n)对应的组织的位移估计；

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

在计算第j(m+1≤j≤M-m)条扫描线的压缩前数据中的某个跟踪波段的位移时，计算了从第j-m到第j+m的2m+1条扫描线的压缩前数据中的相同长度、相同位置的跟踪波段和对应的压缩后数据的互相关函数，将计算得到的2m+1个互相关函数进行加权平均，得到二维综合互相关函数，利用二维综合互相关函数进行组织纵向位移的估计。

实施例2的组织模型和参数设计和实施例1相同，只是m的取值为2，各扫描线数据对应的互相关函数的权重α_k(-m≤k≤m)均为1/5，具体步骤也同实施例1相同，只是在m和α_k的取值不同。

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

图4为一般的组织位移估计方法得到的组织应变分布的计算机仿真结果，即综合系数m为0时的结果；图5为实施例1的得到的组织应变分布的计算机仿真结果，即采用综合系数m为1的二维综合互相关的组织位移估计方法得到的结果；图6为实施例2的得到的组织应变分布的计算机仿真结果，即采用综合系数m为2的二维综合互相关的组织位移估计方法得到的结果；图4-6中，横向和纵向分别表示组织的横向位置和纵向位置(即组织深度)，灰度表示估计出来的应变大小，灰度值越大(即颜色越亮或越白)，表示应变越大，灰度值越小(即颜色越暗或越黑)，表示应变越小，41、51和61分别为图4、图5和图6的灰度值与应变大小的对照关系。

比较图4-6中的组织边缘部位42、43、52、53、62、63，可见，采用一般的组织位移估计方法时，组织横向位移引入的影响较大，组织纵向位移估计的精度较差，应变分布估计的信噪比较低。而采用m＝1和m＝2的二维综合互相关的组织位移估计方法，可以减小由于组织横向位移引入的影响，提高纵向位移估计的精度，提高应变分布估计的信噪比。并且，m越大，二维综合互相关的方法的效果越明显。

Claims (2)

1、一种二维综合互相关的组织位移估计方法，包括以下步骤：

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

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

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

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

3)依次从组织压缩前、后的二维射频信号中取出第m+1条扫描线的前m条和后m条扫描线的数据，即第1到m条扫描线和第m+2到2m+1条扫描线的数据，压缩前、后的扫描线数据分别设为s_1，1(n)、s_1，2(n)、…s_1，m(n)、s_1，m+2(n)、s_1，m+3(n)、…、s_1，2m+1(n)和s_2，1(n)、s_2，2(n)、…、s_2，m(n)、s_2，m+2(n)、s_2，m+3(n)、…、s_2，2m+1(n)，依次取与d_1，m+1同样长度和同样序号的一小段数据d_1，1(n)、d_1，2(n)、…、d_1，m(n)、d_1，m+2(n)、d_1，m+3(n)、…、d_1，2m+1(b)，利用与步骤2相同的计算方法，在τ₁到τ₂确定的搜索范围内分别求小段数据d_1，1(n)、d_1，2(n)、…、d_1，m(n)、d_1，m+2(n)、d_1，m+3(n)、…d_1，2m+1(n)与对应的压缩后的扫描线数据的互相关函数，即求d_1，1(n)与s_2，1(n)的互相关函数R₁(τ)，求d_1，2(n)与s_2，2(n)的互相关函数R₂(τ)，…，求d_1，m(n)与s_2，m(n)的互相关函数R_m(τ)，求d_1，m+2(n)与s_2，m+2(n)的互相关函数R_m+2(τ)，求d_1，m+3(n)与s_2，m+3(n)的互相关函数R_m+3(τ)，…，求d_1，2m+1(n)与s_2，2m+1(n)的互相关函数R_2m+1(τ)；

4)对步骤2-3得到的互相关函数R₁(τ)、R₂(τ)、…、R_m(τ)、R_m+1(τ)、R_m+2(τ)、…、R_2m+1(τ)进行加权平均，得到包含了射频信号二维信息的复合的互相关函数，即二维综合互相关函数R_m+1′(τ)，计算公式如下

${R_{m+1}}^{\′}\left(\mathrm{\τ}\right)=\underset{k=-m}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_k{R}_{m+1+k}\left(\mathrm{\τ}\right)$

其中，α_k代表各扫描线数据对应互相关函数的权重，α_k＞0(-m≤k≤m)，并且所有扫描线数据对应互相关函数的权重之和为1，即 $\underset{k=-m}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\α}}_k=1;$

5)确定该互相关函数R_m+1′(τ)的最大值对应的位置t₁，t₁就是数据d_1，m+1在组织压缩后的位移，即s_1，m+1(n)中序号从n₁到n₁+U-1的小段数据d_1，m+1在组织压缩后移动到s_2，m+1(n)中的序号从n₁-t₁到n₁+U-1-t₁的位置；

6)依次从该扫描线数据s_1，m+1(n)中取一小段长度为T即数据个数为U的数据d_2，m+1、d_3，m+1、…、d_N，m+1，每段数据的序号依次错开V个采样数据，直到再错开V个采样数据将超出s_1，m+1(n)的范围，按步骤2-5相同的方法依次得到各段数据对应的位移t₂、t₃、…、t_N，其中N为小段数据的总数；则位移序列t₁、t₂、…、t_N为第m+1条扫描线数据s_1，m+1(n)对应的组织的位移估计；

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

2、如权利要求1所述的二维综合互相关的组织位移估计方法，其特征在于，所述的二维综合互相关的综合系数1≤m≤5。

Images