CN106842212B - 基于特征空间分解的多重变迹快速自适应波束合成方法 - Google Patents

Abstract

本发明提供一种基于特征空间分解的多重变迹快速自适应波束合成方法：预先在给定的成像环境下采用EIBMV波束合成得到各通道的权值向量，对此权值向量进行主成分分析得到的主成分作为真正的EIBMV权值进行波束合成，波束合成的权值可以采用此预先求知的权值向量，因此简化计算过程。本发明将提高对比度的EIBMV波束合成与降低计算量的PCA、提高图像分辨率的DAX和NSI组合在一起，实现了一维和二维基于特征空间分解的多重变迹快速自适应波束合成方法，克服了MV波束合成得到图像对比度、分辨率不佳等问题。

Description

基于特征空间分解的多重变迹快速自适应波束合成方法

技术领域

本发明属于超声成像领域，具体针对超声成像的波束合成过程，旨在提供一种可保证超声图像质量的同时降低运算过程的复杂度、提高成像效率的基于特征空间分解的多重变迹快速自适应波束合成方法。

背景技术

在现代超声设备中，有两大要素决定了最终成像的质量，其一是设备的硬件水平，其二就是成像方法，而成像过程中最关键的技术就是波束合成。相比成像方法而言，硬件设施的升级更新需要很高的工艺水平，也意味着成本需求的增加，因此新型高效的波束合成方法的研究一直是超声成像领域的热点。

目前应用最为广泛的波束合成方法是延时叠加(Delay and Sum，DAS)波束合成，其基本思路是将各阵元接收到的信号根据其与目标点的距离进行延时，然后将延时后的信号相加得到准确的组织回波信号。但是DAS波束合成方法存在着很多问题：

1)DAS波束合成所形成的波束主瓣宽度较宽，旁瓣水平较高，图像空间分辨率和对比度差；

2)为降低旁瓣水平，DAS波束合成一般采用的是传统的变迹技术，即利用窗函数(例如汉宁窗)给参与合成波束的各接收阵元一个预先给定的权值，但是以增加主瓣宽度，即降低分辨率为代价。

在自适应波束合成中，权值的选择不再与接收信号无关，而是根据各阵元接收信号的特点来调整权值以使得在合成方向上信号能量最大，而尽量压制来自其他方向的信号和干扰，其中最小方差(Minimum Variance，MV)自适应波束合成的应用最为广泛，相比于DAS，MV波束合成在降低旁瓣水平的同时并没有牺牲分辨率，但其中仍然存在一些问题：

1)MV波束合成对于旁瓣和杂波的抑制程度不够，在图像对比度的提高上没有很大作用；

2)MV波束合成得到的图像横向分辨率可以进一步改善；

3)MV波束合成过程中需要很大的计算量，难以用于实时成像。

针对以上不足，BabakMohammadzadehAsl提出了基于特征空间的MV波束合成，即EIBMV(Eigenspace-based Minimum Variance，EIBMV)波束合成，其中各接收阵元权值的获取是通过将MV波束合成中的权值映射到接收信号协方差矩阵的特征空间，使得最终得到的超声图像在保持MV波束合成的高的横向分辨率的同时，有效地压制了旁瓣水平，提高了对比度。零值剪影成像(Null Subtraction Imaging,NSI)是在得到各通道权值之后，通过进一步的变迹方法来降低旁瓣水平同时改善横向分辨率。

而Chi HyungSeo针对旁瓣和杂波抑制提出了新的方法——互相关双重变迹方法(Dual Apodization with cross-correlation，DAX)。对于各通道信号的加权又称为变迹，而在传统的波束合成中权值都是通过一次变迹得到的。在DAX中使用了一对变迹函数。这两个变迹函数应用之后会得到相似的主瓣信号(来自于合成方向上)和存在很大差异的旁瓣信号(来自于其他干扰方向)，因此将两次变迹得到的信号相加并与两者的归一化互相关系数相乘就得到了旁瓣和杂波被有效压制的所需信号。

为降低MV波束合成的计算复杂度，减少计算量，Kyuhong Kim等人提出了快速MV波束合成，通过主成分分析(Principle Component Analysis，PCA)降低信号自相关矩阵的维数来降低波束合成的复杂度。

在目前的波束合成方法中，大都是针对某一方面的不足进行改善。因此，综合解决对比度提高、分辨率改善、计算量减少问题的波束合成方法，将大大增强超声成像的有效性。

发明内容

针对目前的超声图像存在的图像对比度和空间分辨率不佳的问题，本发明的目的在于提供一种基于特征空间分解的多重变迹快速自适应波束合成方法。

为达到上述目的，本发明采用了以下技术方案：

一种基于特征空间分解的快速自适应波束合成方法，该自适应波束合成方法包括以下步骤：在一定的成像条件下，将用于波束合成的对应各通道回波信号的自相关矩阵进行特征空间分解，根据给定阈值挑选出其中较大的特征向量构造信号子空间，将利用MV波束合成得到的权值向量投影到此信号子空间得到新的权值向量；将此权值向量进行主成分分析，将得到的主成分作为用于所述成像条件下进行基于特征空间的最小方差自适应波束合成的权值向量。只要用于波束合成的数据是在完全一致的成像条件下获得，波束合成的所需的权值计算过程可以省略，都可以使用此预先计算得到的权值向量。

进一步的，上述自适应波束合成方法具体包括以下步骤：

1)在给定的成像条件下，进行MV波束合成求得各通道权值向量w_MV，并将其作为样本进行主成分分析；

2)计算权值向量的协方差矩阵R_W：

其中，w_q是第q个MV波束合成权值向量，μ是所有权值向量样本的均值，Q为权值向量样本数，Q大于或等于空间平滑所设定的子阵元数目L；(权值向量是指某通道对应的各采样点的权值组成的向量)

3)对R_W进行特征分解：R_W＝VΛV^H；Λ是对角矩阵，该对角矩阵主对角线上的值由R_W的特征值按照降序组成，V由R_W的特征向量组成，组成V的特征向量彼此正交，且与Λ中的特征值一一对应；

4)根据w_MV＝Vβ，将求最优权值的问题转化为求最优β的问题，即：

minβ^HR₁β，满足β^Hv₁＝1

其中R₁＝V^HRV，R＝E[y·y^H]，R是指用于求解w_MV时的各通道回波信号的协方差矩阵，y为经过延时处理之后用于波束合成的各通道回波信号；v₁＝V^Ha，a为方向向量；

5)根据拉格朗日乘子法求得：

6)将R₁重写为R₁＝E[u·u^H]，其中u＝V^Hy＝[u₁，u₂,…,u_P]^T，然后对R₁进行空间平滑处理和深度方向的样本值平均：

其中P为用于波束合成的阵元总数目，L为空间平滑处理中设定的子阵元数目，L≤P/2，K为深度方向上用于平均处理的点数；然后进行对角加载(Diagonal Loading，DL)，即用代替其中γ₁为加载量，Δ₁＝1/P；选择V能量集中的前I个主成分代表MV波束合成的权值空间 满足I＜L；

7)计算

其中：在中插入一个列，与其他向量正交，且该列的每个元素值均为以避免主成分分析中中心化处理造成进而导致β^Hv₁＝1无法成立的问题；

8)将进行特征分解：Λ_S是对角矩阵，该对角矩阵主对角线上的值由的特征值按照降序组成，V_s由的特征向量组成，组成V_S的特征向量彼此正交，且与Λ_S中的特征值一一对应；

9)构造信号子空间E_S：设定Λ_S中特征值的阈值为最大特征值的δ倍(一般为0.1-0.4)，从而从V_s中选择对应特征值大于此阈值的特征向量来构造信号子空间E_S；

10)计算w_EIBMV：

一种基于特征空间分解的多重变迹快速自适应波束合成方法，包括以下步骤：

1)进行波束合成之前改变换能器的工作模式，即以M个阵元为间隔，将换能器中用于波束合成的全部阵元分为两组，第一组阵元对应通道为第(1…M，2M+1…3M，4M+1…5M，…)通道，第二组阵元对应第(M+1…2M，3M+1…4M，5M+1…6M，…)通道。在用第一组阵元进行波束合成时，将第二组阵元对应的通道信号设为零值，在用第二组阵元进行波束合成时，将第一组阵元对应的通道信号设为零值；然后计算两组阵元波束合成后的数据的互相关系数，通过改变M的值，得到多个互相关系数值(个数小于等于P/2)；

2)在全阵元工作模式下进行波束合成得到输出数据；将输出数据与1)中得到的对应于不同M值的各个互相关系数分别相乘后叠加再按最大值进行归一化处理，得到第一组输出数据集；

3)在全阵元工作模式下进行波束合成，将求得的权值向量在合成波束所用阵元的方向分为左右对称的两部分，左侧阵元的权值保持不变，右侧阵元的权值变为原值的相反数，用形成的新的权值向量对每个通道进行加权，得到第二组输出数据集；

4)将2)和3)中得到的两组输出数据集相减后按最大值进行归一化处理，得到成像射频数据，在经过后续的包络检波、对数压缩等进一步的信号处理步骤后即可用于超声成像；

所述1)—3)中，波束合成均采用上述基于特征空间分解的快速自适应波束合成方法。

进一步的，对于二维面阵超声换能器，将二维面阵分为x轴与y轴(参见中国专利“超声二维面阵的三维宽波束小区域快速空化成像方法”，CN104777485A)，然后分别按照步骤1)至步骤4)在x轴及y轴方向进行波束合成(针对x轴及y轴方向上的每一个线阵)。

本发明的有益效果体现在：

为提高成像速度，预先在给定的相同的成像环境下得到各通道的权值向量。此权值向量的计算采用EIBMV波束合成，即将各通道回波信号的自相关矩阵进行特征空间分解，给定阈值挑选出其中较大的特征向量来构造信号子空间，最终的权值向量是将此信号子空间与MV波束合成得到的权值向量相乘得到。对此权值向量进行主成分分析得到的主成分可以作为真正的EIBMV权值进行波束合成，这在统计学中已得到验证。因此，只要成像环境和相关条件确定不变，波束合成的权值可以采用此预先求知的权值向量，因此简化计算过程。

而为了进一步压制旁瓣和杂波以提高最终得到图像的对比度及分辨率，采用DAX多重变迹方法。在每次进行波束合成之前改变换能器的工作模式，即以M个阵元为间隔，将换能器全部阵元分为两组，第一组阵元对应通道为第(1…M，2M+1…3M，4M+1…5M，…)通道，第二组阵元对应第(M+1…2M，3M+1…4M，5M+1…6M，…)通道。在用第一组阵元进行波束合成时，将第二组阵元对应的通道信号设为零值，在用第二组阵元进行波束合成时，将第一组阵元对应的通道信号设为零值。然后计算两组波束合成后的数据的互相关系数，通过改变M的值，可以得到多个互相关系数值，分别与全阵元工作模式下波束合成的数据相乘后叠加再进行归一化处理，此时的各通道信号的旁瓣和其他干扰已经得到较大抑制。

其次，改善分辨率用到NSI零值剪影成像，通过将每次波束合成中求得的权值分为左右对称的两组，左侧保持不变，右侧变为原来的相反数，用此零均值处理的权值对对应通道信号加权。

最后，将多重变迹DAX与NSI两次处理后的数据相减，归一化之后进行包络检波和对数压缩，用于最终成像。

总之，本发明将提高对比度的EIBMV波束合成与降低计算量的PCA、提高图像分辨率的DAX和NSI组合在一起，实现了一维和二维基于特征空间分解的多重变迹快速自适应波束合成方法，克服了MV波束合成得到图像对比度、分辨率不佳等问题。

附图说明

图1为用于采集数据的超声线阵换能器的结构示意图，其中：(a)超声线阵换能器立体视图；(b)俯视图。

图2为应用DAS波束合成得到的超声图像，其中：(a)点散射子；(b)低回声仿体。

图3为应用本发明设计的基于特征空间分解的多重变迹快速一维自适应波束合成方法得到的超声图像，其中：(a)点散射子；(b)低回声仿体。

图4为DAS和本发明设计的波束合成方法(Proposed)所得图像的横向分辨率的对比。

图5本发明中应用的多重变迹技术中改变线阵换能器工作模式的示意图，其中(a)M＝2；(b)M＝4；(c)M＝8。

图6本发明中基于特征空间分解的多重变迹快速二维自适应波束合成方法流程图，其中：1.时间延迟；2.沿x轴基于特征空间分解的多重变迹快速自适应波束合成；3.沿y轴基于特征空间分解的多重变迹快速自适应波束合成；4.二维波束合成输出。

图7为应用本发明设计的基于特征空间分解的多重变迹快速二维自适应波束合成方法得到的超声仿真图像。

具体实施方式

下面结合附图和实施例对本发明做进一步详细说明。

本发明所述基于特征空间分解的多重变迹快速自适应波束合成方法，仍然基于MV波束合成，实现方案如下：

一、基于特征空间分解的多重变迹快速一维自适应波束合成方法

线阵换能器向成像区域发射平面波，遇到目标后发生散射，产生与发射方向相反的回波；利用数据采集卡对超声回波信号进行采样并转化为射频数据，然后将所有阵元的射频数据存储并传送给计算机；根据目标点位置，计算回波信号的时间延迟，并对接收信号进行时间校正：假设用于波束合成的线阵阵元个数为P，因为成像在二维平面，因此空间目标点坐标设为(x,0,z)，接收阵元的回波信号相对于参考阵元的时间延迟为：

其中，i＝1,...,P，c为声速。

则每个阵元的延迟校正后的回波信号为：

y_i(t)＝x_i(t-Δ_i(t))

其中，x_i(t)表示每个阵元原始的接收信号。

波束合成的目标是对接收信号优化加权以去除来自其他非合成方向的杂波等干扰信号，在MV波束合成中权值的选择是在满足合成波束方向上的增益等于1的条件下使输出信号能量最小，即

min E[|z(t)|²]＝min w^H(t)Rw(t)，满足w^H(t)a＝1

其中E[·]表示求输出信号的期望，w(t)＝[ω₁(t),ω₂(t),…,ω_P(t)]表示各通道的权值，w^H(t)是对w(t)的共轭转置，a是方向向量，因为接收信号已经经过延迟，所以a本质上是单位向量。R是各通道信号的协方差矩阵：

R＝E[y(t)y^H(t)]，其中y(t)＝[y₁(t),y₂(t),…,y_P(t)]

通过拉格朗日乘子法求得w(t)：

实际数据处理过程都是在数字域进行，接收的回波信号都需要经过采样转换为离散样本值，因此协方差矩阵应该表示为：

其中，N_s表示每个通道的总采样点数。

因为各阵元接收的超声回波信号是高度相关的，这与MV波束合成对合成方向与其他方向信号不相关的假设是相悖的，因此为了克服此问题：采用子阵元平均或称为空间平滑的方法，即将用于合成波束的P(与后面的NSI方法选择的阵元集合对应)个阵元划分为P-L+1个互相重叠的子阵，L是子阵长度；同时为了保持散斑统计特性，在合成波束深度方向将相邻的2K+1个样本值进行平均，K一般取为两个波长的采样点数，以对应深度为第k个采样点为例，经过上述两个过程得到的协方差矩阵表示为：

其中y^l(n)＝[y_l(n),y_l+1(n),…,y_l+L-1(n)]^T。在计算过程中，只有上述样本值参与，因此为了增强协方差矩阵的鲁棒性，引入对角加载(Diagonal Loading，DL)技术，即用加载协方差矩阵(I为单位阵)代替其中为加载量，Δ＝1/P。因此：

如果每次波束合成都要重复上述步骤，将带来极大的计算量。因此在成像环境和相关条件不变的情况下，可以利用主成分分析PCA来减少计算时间。具体步骤：

1)在给定的成像条件下，先进行MV波束合成求得各通道权值向量w_MV，并将其作为样本进行主成分分析。

2)计算权值向量的协方差矩阵R_W：

其中，w_q是第q个MV权值向量，μ是所有权值向量样本的均值，Q为权值向量样本数，要求要大于或等于空间平滑所设定的子阵元数目；

3)对R_W进行特征分解：R_W＝VΛV^H；Λ是对角矩阵，主对角线上的值由R_W的特征值按照降序组成，V由R_W的特征向量组成，彼此正交，且与Λ中的特征值一一对应。

4)因为V是满秩矩阵，因此MV波束合成中要计算的w_MV都可以表示为组成V的特征向量的线性组合，即w_MV＝Vβ，进而将求最优权值的问题转化为求最优β的问题，即：

minβ^HR₁β，满足β^Hv₁＝1

其中R₁＝V^HRV，R＝E[y·y^H]是指用于求解w_MV时的各通道信号的协方差矩阵，y为经过延时处理之后用于波束合成的各通道信号；v₁＝V^Ha，a为方向向量；

5)根据拉格朗日乘子法求得：

6)将R₁重写为R₁＝E[u·u^H]，其中u＝V^Hy＝[u₁，u₂,…,u_P]^T。然后对R₁进行空间平滑处理和深度方向的样本值平均：

其中P为用于波束合成的阵元总数目，L为空间平滑处理中设定的子阵元数目，要求L≤P/2，K为深度方向上用于平均处理的点数；然后进行对角加载(Diagonal Loading，DL)，即用代替其中为加载量，Δ₁＝1/P；通过主成分分析发现w_MV的能量主要集中在前几个主成分，因此只需要选择很少数目的主成分即可代表MV波束合成的权值空间，用表示，即 由V的前I个主成分组成，要满足I＜L；

7)计算

其中：但是此时计算得到的是不稳定的，因为的每一列的元素和可能为零值，导致进而无法满足这一问题源于在计算R_W时对样本权值向量进行了中心化处理，这是PCA中很重要的一个步骤，却也因此消除了V中每一列的直流成分，使得约束条件很难满足，而这对于基于MV的自适应波束合成又是基本要求。因此，需要给再加入直流成分，即在中插入一个列，与其他向量正交，且该列的每个元素值均为

在上述步骤，通过对预先在相同成像条件下采集的w_MV权值向量作为样本进行主成分分析来降低进行权值计算时所需的协方差矩阵的维数，可以很大程度减少计算量。而对于MV波束合成存在的超声图像对比度不高的问题，用到基于特征空间的自适应波束合成，EIBMV利用已经计算得到的w_MV，将其投影到协方差矩阵的信号子空间。具体计算：

1)将进行特征分解：Λ_S是对角矩阵，主对角线上的值由的特征值按照降序组成，V_s由的特征向量组成，彼此正交，且与Λ_S中的特征值一一对应；

2)构造信号子空间E_S：设定Λ_S中特征值的阈值为最大特征值的δ倍(一般为0.1-0.4)，从而从V_S中选择对应特征值大于此阈值的特征向量来构造信号子空间E_S，以尽可能保证此特征空间包含了主瓣信号而消除了旁瓣和杂波信号；

3)计算w_EIBMV：

以上完成波束合成的权值确定过程。而为了进一步压制旁瓣和杂波同时减小主瓣宽度，提高分辨率，采用NSI零值减影成像技术和DAX多重变迹技术，具体步骤：

1)进行波束合成之前改变换能器的工作模式，即以M个阵元为间隔阵元，将换能器中用于波束合成的全部阵元分为两组，第一组阵元对应通道为第(1…M，2M+1…3M，4M+1…5M，…)通道，第二组阵元对应第(M+1…2M，3M+1…4M，5M+1…6M，…)通道。在用第一组阵元进行波束合成时，将第二组阵元对应的通道信号设为零值，在用第二组阵元进行波束合成时，将第一组阵元对应的通道信号设为零值；然后计算两组阵元波束合成后的数据的互相关系数，通过改变M的值，得到多个互相关系数值；

2)在全阵元工作模式下进行波束合成得到输出数据；与1)中得到的对应于不同M值的各个互相关系数相乘后叠加再进行归一化处理，得到第一组输出数据集；

3)在全阵元工作模式下进行波束合成，将求得的权值向量在合成波束所用阵元的方向分为左右对称的两部分，左侧的权值保持不变，右侧的权值变为原值的相反数，用形成的新的权值向量对每个通道进行加权，得到第二组输出数据集；

4)将2)和3)中得到的两组数据集相减后进行归一化处理，在经过后续的包络检波、对数压缩等进一步的信号处理步骤后即可用于超声成像。

二、基于特征空间分解的多重变迹快速二维自适应波束合成方法

将二维面阵分为x轴与y轴；首先按照x轴对时延后的信号进行基于特征空间分解的多重变迹快速一维自适应波束合成方法中的快速超声自适应波束合成的计算，随后进行基于特征空间分解的多重变迹快速一维自适应波束合成方法中的多重变迹叠加；在此基础上，再按照y轴重复上述步骤，得到最终的二维基于特征空间分解的多重变迹快速自适应波束合成的结果，具体的：

2.1沿x轴波束合成

1)在给定的成像条件下，先按照x轴进行MV波束合成求得权值向量w_xMV，并将其作为样本进行主成分分析。要保证预先求得的权值向量的个数不小于空间平滑中设定的沿x轴子阵长度L；

2)计算权值向量的协方差矩阵R_xW：

其中，w_xq是沿x轴第q个MV权值向量，μ_x是所有权值向量样本的均值，Q为权值向量样本数，要求要大于或等于空间平滑所设定的子阵元数目；

3)对R_xW进行特征分解：R_xW＝V_xΛ_xV_x ^H；V_x是对角矩阵，主对角线上的值由R_xW的特征值按照降序组成，V_x由R_xW的特征向量组成，彼此正交，且与Λ_x中的特征值一一对应；

4)将MV波束合成中要计算的w_xMV可以表示为组成V_x的特征向量的线性组合，即w_xMV＝V_xβ_x，进而将求最优权值的问题转化为求最优β_x的问题，即：

minβ_x ^HR_x1β_x，满足β_xβ_x ^Hv_x1＝1

其中R_x1＝V_x ^HR_xV_x，R_x＝E[y_x·y_x ^H]是指用于求解w_xMV时的各通道信号的协方差矩阵，y_x为经过延时处理之后用于x方向波束合成的各通道信号；v_x1＝V_x ^Ha，a为方向向量；

5)根据拉格朗日乘子法求得：

6)将R_x1重写为其中同样对R_x1进行空间平滑处理和深度方向的样本值平均：

其中然后进行对角加载DL，即代替其中为加载量，Δ₁＝1/P。通过主成分分析发现w_xMV的能量主要集中在前几个主成分，因此只需要选择很少数目的主成分即可代表MV波束合成的权值空间，用表示，由V_x的前I个主成分组成，要满足I＜L；

7)计算其中：并给加入直流成分；

8)按照特征分解的方式将进行特征分解，并按照前述构造信号子空间E_S的方式构造信号子空间E_xS；

9)计算w_xEIBMV：

以上完成沿x方向波束合成的权值确定过程。而为了进一步压制旁瓣和杂波同时减小主瓣宽度，提高分辨率，采用NSI零值减影成像技术和DAX多重变迹技术，具体步骤：

1)进行波束合成之前改变换能器的工作模式，即沿x轴以M个阵元为间隔阵元，将换能器中用于波束合成的全部阵元分为两组，第一组阵元对应通道为沿x方向第(1…M，2M+1…3M，4M+1…5M，…)通道，第二组阵元对应第(M+1…2M，3M+1…4M，5M+1…6M，…)通道。在用第一组阵元进行波束合成时，将第二组阵元对应的通道信号设为零值，在用第二组阵元进行波束合成时，将第一组阵元对应的通道信号设为零值；然后计算两组阵元波束合成后的数据的互相关系数，通过改变M的值，可以得到多个互相关系数值；

2)在全阵元工作模式下进行波束合成得到输出数据集；与1)中得到的对应于不同M值的各个互相关系数相乘后叠加再进行归一化处理，得到第一组输出数据集。此时的各通道信号的旁瓣和其他干扰已经得到较大抑制；

3)在全阵元工作模式下进行波束合成，将每次波束合成求得的w_xEIBMV在合成波束所用阵元的方向分为左右对称的两部分，左侧的权值保持不变，右侧的权值变为原值的相反数；用形成的新的权值向量对每个通道进行加权，得到第二组输出数据集；

4)将2)和3)中得到的两组数据集相减后得到沿x轴方向波束合成结果。

2.2沿y轴波束合成

1)在给定的成像条件下，先按照y轴进行MV波束合成求得权值向量w_yMV，并将其作为样本进行主成分分析。要保证预先求得的权值向量的个数不小于空间平滑中设定的沿y轴子阵长度L；

2)计算权值向量的协方差矩阵R_yW：

其中，w_yq是沿y轴第q个MV权值向量，μ_y是所有权值向量样本的均值，Q为权值向量样本数，要求要大于或等于空间平滑所设定的子阵元数目；

3)对R_yW进行特征分解：R_yy＝V_yΛ_yV_y ^H；V_y是对角矩阵，主对角线上的值由R_yw的特征值按照降序组成，V_y由R_yw的特征向量组成，彼此正交，且与Λ_y中的特征值一一对应；

4)将MV波束合成中要计算的w_yMV都可以表示为组成V_y的特征向量的线性组合，即w_yMV＝V_yβ_y，进而将求最优权值的问题转化为求最优β_y的问题，即：

minβ_y ^HR_y1β_y，满足β_yβ_y ^Hv_y1＝1

其中R_y1＝V_y ^HR_yV_y，R_y＝E[y_y·y_y ^H]是指用于求解w_yMV时的各通道信号的协方差矩阵，y_y为经过延时处理之后用于xy方向波束合成的各通道信号；v_y1＝V_y ^Ha，a为方向向量。

5)根据拉格朗日乘子法求得：

6)将R_y1重写为其中同样对R_y1进行空间平滑处理和深度方向的样本值平均：

其中然后进行对角加载DL，即代替其中为加载量，Δ₁＝1/P。通过主成分分析发现w_yMV的能量主要集中在前几个主成分，因此只需要选择很少数目的主成分即可代表MV波束合成的权值空间；

7)计算其中：并给加入直流成分；

8)按照特征分解的方式将进行特征分解，并按照前述构造信号子空间E_S的方式构造信号子空间E_yS；

9)计算w_yEIBMV：

以上完成波束合成的权值确定过程。而为了进一步压制旁瓣和杂波同时减小主瓣宽度，提高分辨率，采用NSI零值减影成像技术和DAX多重变迹技术，具体步骤：

1)进行波束合成之前改变换能器的工作模式，即沿y轴以M个阵元为间隔阵元，将换能器中用于波束合成的全部阵元分为两组，第一组阵元对应通道为沿y方向第(1…M，2M+1…3M，4M+1…5M，…)通道，第二组阵元对应第(M+1…2M，3M+1…4M，5M+1…6M，…)通道。在用第一组阵元进行波束合成时，将第二组阵元对应的通道信号设为零值，在用第二组阵元进行波束合成时，将第一组阵元对应的通道信号设为零值；然后计算两组阵元波束合成后的数据的互相关系数，通过改变M的值，可以得到多个互相关系数值；

2)在全阵元工作模式下进行波束合成得到输出数据集；与1)中得到的对应于不同M值的各个互相关系数相乘后叠加再进行归一化处理，得到第一组输出数据集。此时的各通道信号的旁瓣和其他干扰已经得到较大抑制；

3)在全阵元工作模式下进行波束合成，将每次波束合成求得的w_yEIBMV在合成波束所用阵元的方向分为左右对称的两部分，左侧的权值保持不变，右侧的权值变为原值的相反数；用形成的新的权值向量对每个通道进行加权，得到第二组输出数据集；

4)将2)和3)中得到的两组数据集相减后得到沿y轴波束合成结果。

2.3在两个方向波束合成后最终得到的数据分别进行归一化处理，再经过后续的包络检波、对数压缩等进一步的信号处理步骤后即可用于超声三维成像(三维成像参见中国专利“超声二维面阵的三维宽波束小区域快速空化成像方法”，CN104777485A)。

三、发明效果验证

1、采用基于特征空间分解的多重变迹快速一维自适应波束合成方法

以图1所示超声线阵换能器的阵列为例，其单阵元尺寸为0.27mm×0.27mm，换能器阵元间距为0.3mm，阵元与通道数为128。发射中心频率为7MHz，焦点距离线阵几何中心20mm。

多重变迹技术中超声换能器的工作模式如下：

进行波束合成之前，将线阵换能器的阵元以M个阵元为间隔分为两组，再各自进行波束合成。具体在得到图3的超声图像时，设定了3个M的值，分别为2、4、8(如图5所示，其中ON表示的为两组中的一组)。第一次将阵元以2为间隔分组，将得到的合成后的输出数据进行互相关计算去掉二者共有的杂波信号而保留各自独有的主瓣信号特征，给得到的系数矩阵设定一个阈值ρ，即大于ρ的系数保持不变，小于ρ的系数设定其值等于ρ，然后与全阵元工作模式下得到的波束合成数据相乘；改变M的值，重复以上操作，将得到的三组数据相加后进行归一化处理，这就是多重变迹技术在本发明中的应用。

上述基于特征空间分解的多重变迹快速一维自适应波束合成方法与DAS波束合成的对比如图2、图3所示，用图1中的线阵换能器发射平面波，将采集到的两组数据分别用DAS波束合成和基于特征空间分解的多重变迹快速自适应波束合成方法进行处理。得到的超声图像表明，本发明实现的自适应波束合成方法在图像对比度、分辨率方面要明显优于DAS波束合成方法。

如图4所示，对以上两种波束合成方法得到的超声图像的横向分辨率进行量化对比，选取深度为35mm处的图像信息。结果表明用本发明的自适应波束合成方法处理得到的图像横向分辨率相比DAS波束合成方法提高了大约2倍。

2、基于特征空间分解的多重变迹快速二维自适应波束合成方法

所采用二维面阵超声换能器(举例说明)，其二维面阵规格32×32，单阵元尺寸为0.3mm×0.3mm，换能器阵元间距为0.1mm。发射中心频率为7MHz。

方法流程如图6所示，首先按照x轴对时延后的信号进行一维基于特征空间分解的多重变迹快速自适应波束合成方法中的快速超声自适应波束合成的计算，随后进行一维基于特征空间分解的多重变迹快速自适应波束合成方法中的多重变迹叠加；在此基础上，再按照y轴重复上述步骤，得到最终的二维基于特征空间分解的多重变迹快速自适应波束合成的结果。

仿真结果如图7所示，其结果能够高分辨显示空间中各个目标的分布。

Claims (6)

1.一种基于特征空间分解的快速自适应波束合成方法，其特征在于：该自适应波束合成方法包括以下步骤：

在一定的成像条件下，将用于波束合成的对应各通道回波信号的自相关矩阵进行特征空间分解，根据给定阈值挑选出其中较大的特征向量构造信号子空间，将MV波束合成得到的权值向量投影到此信号子空间得到新的权值向量；将此权值向量进行主成分分析，将得到的主成分作为用于所述成像条件下进行基于特征空间的最小方差自适应波束合成的权值向量；

所述自适应波束合成方法具体包括以下步骤：

1)在给定的成像条件下，进行MV波束合成求得各通道权值向量w_MV；

2)计算权值向量的协方差矩阵R_W：

其中，w_q是第q个MV波束合成权值向量，μ是所有权值向量样本的均值，Q为权值向量样本数，Q大于或等于空间平滑所设定的子阵元数目L；

3)对R_W进行特征分解：R_W＝VΛV^H；Λ是对角矩阵，该对角矩阵主对角线上的值由R_W的特征值按照降序组成，V由R_W的特征向量组成，组成V的特征向量彼此正交，且与Λ中的特征值一一对应；

4)根据w_MV＝Vβ，将求最优权值的问题转化为求最优β的问题，即：

minβ^HR₁β，满足β^Hν₁＝1

其中R₁＝V^HRV，R＝E[y·y^H]是指用于求解w_MV时的各通道信号的协方差矩阵，y为经过延时处理之后用于波束合成的各通道信号；v₁＝V^Ha，a为方向向量；

5)根据拉格朗日乘子法求得：

6)将R₁重写为R₁＝E[u·u^H]，其中u＝V^Hy＝[u₁,u₂,…,u_P]^T，然后对R₁进行空间平滑处理和深度方向的样本值平均：

其中P为用于波束合成的阵元总数目，L为空间平滑处理中设定的子阵元数目，L≤P/2，K为深度方向上用于平均处理的点数；然后用代替其中为加载量，Δ₁＝1/P；选择V能量集中的前I个主成分代表MV波束合成的权值空间 满足I＜L；

7)计算

其中：在中插入一个列，与其他向量正交，且该列的每个元素值均为

8)将进行特征分解：Λ_S是对角矩阵，该对角矩阵主对角线上的值由的特征值按照降序组成，V_S由的特征向量组成，组成V_S的特征向量彼此正交，且与Λ_S中的特征值一一对应；

9)构造信号子空间E_S：设定Λ_S中特征值的阈值为最大特征值的δ倍，从而从V_s中选择对应特征值大于此阈值的特征向量来构造信号子空间E_S；

10)计算权值向量w_EIBMV：

2.根据权利要求1所述一种基于特征空间分解的快速自适应波束合成方法，其特征在于：所述δ为0.1-0.4。

3.一种基于特征空间分解的多重变迹快速自适应波束合成方法，其特征在于：包括以下步骤：

1)进行波束合成之前改变换能器的工作模式，即以M个阵元为间隔，将换能器中用于波束合成的全部阵元分为两组，在用第一组阵元进行波束合成时，将第二组阵元对应的通道信号设为零值，在用第二组阵元进行波束合成时，将第一组阵元对应的通道信号设为零值；然后计算两组阵元波束合成后的数据的互相关系数，通过改变M的值，得到多个互相关系数值；

2)在全阵元工作模式下进行波束合成得到输出数据；将输出数据与1)中得到的对应于不同M值的各个互相关系数分别相乘后叠加再按最大值进行归一化处理，得到第一组输出数据集；

3)在全阵元工作模式下进行波束合成，将求得的权值向量在合成波束所用阵元的方向分为左右对称的两部分，左侧阵元的权值保持不变，右侧阵元的权值变为原值的相反数，用形成的新的权值向量对每个通道进行加权，得到第二组输出数据集；

4)将2)和3)中得到的两组输出数据集相减后按最大值进行归一化处理，得到成像射频数据；

所述1)、2)以及3)中，波束合成均采用基于特征空间分解的快速自适应波束合成方法，该自适应波束合成方法在一定的成像条件下，将用于波束合成的对应各通道回波信号的自相关矩阵进行特征空间分解，根据给定阈值挑选出其中较大的特征向量构造信号子空间，将MV波束合成得到的权值向量投影到此信号子空间得到新的权值向量；将此权值向量进行主成分分析，将得到的主成分作为用于所述成像条件下进行基于特征空间的最小方差自适应波束合成的权值向量。

4.根据权利要求3所述一种基于特征空间分解的多重变迹快速自适应波束合成方法，其特征在于：所述自适应波束合成方法具体包括以下步骤：

1)在给定的成像条件下，进行MV波束合成求得各通道权值向量w_MV；

2)计算权值向量的协方差矩阵R_W：

其中，w_q是第q个MV波束合成权值向量，μ是所有权值向量样本的均值，Q为权值向量样本数，Q大于或等于空间平滑所设定的子阵元数目L；

3)对R_W进行特征分解：R_W＝VΛV^H；Λ是对角矩阵，该对角矩阵主对角线上的值由R_W的特征值按照降序组成，V由R_W的特征向量组成，组成V的特征向量彼此正交，且与Λ中的特征值一一对应；

4)根据w_MV＝Vβ，将求最优权值的问题转化为求最优β的问题，即：

minβ^HR₁β，满足β^Hν₁＝1

其中R₁＝V^HRV，R＝E[y·y^H]是指用于求解w_MV时的各通道信号的协方差矩阵，y为经过延时处理之后用于波束合成的各通道信号；v₁＝V^Ha，a为方向向量；

5)根据拉格朗日乘子法求得：

6)将R₁重写为R₁＝E[u·u^H]，其中u＝V^Hy＝[u₁,u₂,…,u_P]^T，然后对R₁进行空间平滑处理和深度方向的样本值平均：

其中P为用于波束合成的阵元总数目，L为空间平滑处理中设定的子阵元数目，L≤P/2，K为深度方向上用于平均处理的点数；然后用代替其中为加载量，Δ₁＝1/P；选择V能量集中的前I个主成分代表MV波束合成的权值空间 满足I＜L；

7)计算

其中：在中插入一个列，与其他向量正交，且该列的每个元素值均为

8)将进行特征分解：Λ_s是对角矩阵，该对角矩阵主对角线上的值由的特征值按照降序组成，V_s由的特征向量组成，组成V_s的特征向量彼此正交，且与Λ_S中的特征值一一对应；

9)构造信号子空间E_S：设定Λ_s中特征值的阈值为最大特征值的δ倍，从而从V_s中选择对应特征值大于此阈值的特征向量来构造信号子空间E_S；

10)计算权值向量w_EIBMV：

5.根据权利要求4所述一种基于特征空间分解的多重变迹快速自适应波束合成方法，其特征在于：所述互相关系数值的个数小于等于P/2。

6.根据权利要求3所述一种基于特征空间分解的多重变迹快速自适应波束合成方法，其特征在于：对于二维面阵超声换能器，将二维面阵分为x轴与y轴，然后分别按照步骤1)至步骤4)在x轴及y轴方向进行波束合成。