CN102499712B - 一种基于特征空间的前后向自适应波束形成方法 - Google Patents

Abstract

本发明公开了一种基于特征空间的前后向自适应波束形成方法，涉及医学超声成像技术领域，首先对接收阵元的多路采样信号进行聚焦延时处理和前后向平滑处理，得到样本协方差矩阵估计，并对其进行对角加载后与方向向量结合，计算得到自适应波束形成权值；然后将对角加载后的前后向协方差矩阵估计进行特征分解，构造信号子空间，将自适应波束形成权值投影到该信号子空间中，得到新的自适应波束形成权值；最后将新的自适应波束形成权值对经过前向平滑处理的多路数据进行加权求和，从而得到一路自适应波束信号；它解决了现有的自适应波束形成算法在提高图像分辨率、对比度以及对方向误差敏感等问题，全面提高了超声成像的整体质量。

Description

一种基于特征空间的前后向自适应波束形成方法

技术领域

本发明属于超声成像技术领域，具体涉及一种基于特征空间的前后向自适应波束形成方法。

背景技术

在医学超声成像中，来自非信源方向上的回波信号引入了广泛的图像杂波，导致图像质量下降。目前，一般采用传统延时叠加方法(DAS)减弱非信源方向上信号的干扰，提高信源方向上信号的强度。该方法对不同通道接收到的数据施以不同的延时，再将其相加，得到期望聚焦点的幅值。DAS算法简单，易实现，但其空间分辨率较低，成像对比度较差，对方向向量误差敏感。为了减小旁瓣信号对成像对比度的影响，人们通常采用幅度变迹技术，对不同的阵元施加不同的加权值(通常为Hanning或Blackman函数)。幅度变迹虽然压制了旁瓣，但代价却是增加了主瓣宽度，降低了成像的空间分辨率。产生问题的根源在于，用于幅度变迹的加权值都是固定的，同空间中散射点的分布以及接收到的数据是无关的。

由此，引出了人们对自适应波束形成的研究，希望根据接收的数据计算出动态的加权值，从而达到提高成像空间分辨率的目的。目前最常用的自适应波束形成算法是最小方差波束形成算法(MV)。该方法是由Capon在1969年提出来的，基本思想是在无失真约束条件下在指定方位或频率上获得最小方差，进而推导出最优权向量。虽然Capon自适应波束形成算法对于提高图像的分辨率有其独特的优势，但是算法的稳健性不如传统的延时叠加算法，如果聚焦方向不够精确，成像效果会严重下降，而且不能显著提高图像的对比度。

综上所述，急需发明一种能够同时提高图像对比度、分辨率，降低方向误差敏感度，并且稳健性好的波束形成算法，以全面提高超声图像的整体质量。

发明内容

有鉴于此，为了解决上述问题，本发明提出一种能够同时提高图像对比度、分辨率以及波束形成算法的稳健性，解决对方向误差比较敏感的问题，有效克服了传统的自适应波束形成算法稳健性低，不能显著提高图像对比度等问题，全面提高了超声图像的整体质量。

本发明的目的是提出一种基于特征空间的前后向自适应波束形成方法；

本发明的目的是这样实现的：

本发明提供的一种基于特征空间的前后向自适应波束形成方法，包括以下步骤：

S1：接收传感器阵列的采样信号；

S2：对采样信号进行前后向平滑处理，得到样本协方差矩阵估计；

S3：对样本协方差矩阵估计进行对角加载，得到空间协方差矩阵估计；

S4：将空间协方差矩阵估计与方向向量结合，计算得到自适应波束形成权值；

S5：对空间协方差矩阵估计进行特征分解，构造信号子空间；

S6：将自适应波束形成权值投影到信号子空间中，得到新的自适应波束形成权值；

S7：将新的自适应波束形成权值对采样信号数据进行加权求和，得到自适应波束信号。

进一步，所述步骤2中进行前后向平滑处理前还包括对采样信号进行聚焦延时处理，得到聚焦延时后的信号x(k)；所述信息x(k)表示为x(k)＝[x₁(k)，x₂(k)，...，x_M(k)]，其中，M表示传感器阵列的阵元个数；

进一步，所述步骤2中进行前后向平滑处理，具体包括以下步骤：

S21：把M个阵元的数据分为阵元数目为L的子阵，分别计算各个子阵的相关矩阵，然后根据以下公式来计算，得到前向协方差矩阵估计为

$\tilde{R}\left(k\right)=\frac{1}{M-L+1}\underset{l=1}{\overset{M-L+1}{\mathrm{\Σ}}}{x}_d^l\left(k\right)x_d^l{\left(k\right)}^H,$

公式中， $x_d^l\left(k\right)={[x_d^l\left(k\right),x_d^{l+1}\left(k\right),\·\·\·,x_d^{l+L+1}\left(k\right)]}^T$ 为第l个子阵的输出向量，表示的共轭转置；

S22：定义后向重叠向量为 ${\tilde{x}}_d^l\left(k\right)={[x_d^{M-l+1}\left(k\right),x_d^{M-l}\left(k\right),\·\·\·,x_d^{M-l+L+2}\left(k\right)]}^H,$ l＝1，2，...，M；通过以下公式计算，得到后向协方差矩阵估计

${\tilde{R}}^b\left(k\right)=\frac{1}{M-L+1}\underset{l=1}{\overset{M-L+1}{\mathrm{\Σ}}}{\tilde{x}}_d^l\left(k\right){\tilde{x}}_d^l{\left(k\right)}^H,$

其中，表示第l个子阵的后向输出向量，表示的共轭转置；

S23：通过以下公式来计算前后向协方差矩阵估计：

${\tilde{R}}_{\mathrm{FB}}\left(k\right)=\frac{1}{2}(\tilde{R}\left(k\right)+{\tilde{R}}^b\left(k\right)),$

其中，表示前向协方差矩阵估计与后向协方差矩阵估计的求和平均；

进一步，通过下列公式对前后向协方差矩阵进行对角加载得到对角加载后的前后向协方差矩阵估计

${\tilde{R}}_{\mathrm{DL}}={\tilde{R}}_{\mathrm{FB}}+\mathrm{\γI},$

其中，Δ为算法加入的空间噪声与信号功率比，为信号等效功率，I表示单位矩阵；

进一步，所述步骤S4中，通过以下公式来计算自适应波束形成权值：

$w=\frac{{\tilde{R}}_{\mathrm{DL}}^{-1}a}{a^H{\tilde{R}}_{\mathrm{DL}}^{-1}a},$

公式中，a为方向向量，w为自适应波束形成权值，表示的逆矩阵；

进一步，所述步骤S5中，通过以下公式对进行特征分解可得：

${\tilde{R}}_{\mathrm{DL}}=\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{e}_i{e}_i^H=E_s{\mathrm{\Λ}}_s{E}_s^H+E_n{\mathrm{\Λ}}_n{E}_n^H,$

其中，λ₁≥λ₂≥…≥λ_L是相应的L个特征值，其对应的特征向量分别为e_i，i＝1，2，...，L；用大于最大特征值0.5倍的特征值所对应的特征向量组成信号子空间E_s，而其余特征值对应的特征向量组成噪声子空间E_n，表示E_s的共轭转置，表示E_n的共轭转置，num为信号子空间中特征向量的个数，则Λ_s＝diag{λ₁，λ₂，…，λ_num}，Λ_n＝diag{λ_num+1，…，λ_L}，E_s＝[e₁，e₂，…，e_num]，E_n＝[e_num+1，…，e_L]；

进一步，所述步骤S6中新的自适应波束形成权值通过以下公式来计算：

$w_1=E_s{E}_s^{H}w,$

其中，w₁表示将自适应波束形成权值投影到由协方差矩阵特征空间构造的信号子空间E_s后的自适应波束形成权值；

进一步，所述步骤S7中，将投影后的自适应波束形成权值对所述的聚焦延时后的采样信号通过以下公式进行加权求和，计算得到自适应波束信号：

$y\left(k\right)=\frac{1}{M-L+1}\underset{l=1}{\overset{M-L+1}{\mathrm{\Σ}}}{w_1}^H\left(k\right)x_d^l\left(k\right),$

其中，y(k)表示计算得到的自适应波束信号，w₁ ^H表示w₁的共轭转置，表示第l个子阵的输出向量；

进一步，所述子阵阵元数目L的取值为M/2；

进一步，所述对角加载因子γ中Δ的值取为1/100L，所述方向向量a为全1的矢量。

本发明的优点在于：本发明采用一种基于特征空间的前后向自适应波束形成方法，利用前后向平滑技术来计算样本协方差矩阵，提高其精度，利用对角加载提高算法的稳健性；然后利用特征空间法将自适应波束形成权值投影到由样本协方差矩阵特征空间构造的信号子空间中，提高图像的对比度，降低旁瓣等级；解决了现有的自适应波束形成算法在提高图像分辨率、对比度以及对方向误差比较敏感等问题，全面提高了超声成像的整体质量，克服了传统的自适应波束形成算法稳健性低，不能显著提高图像对比度等问题。

本发明的其它优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导与启发。本发明的目标和其它优点可以通过下面的说明书，权利要求书，以及附图中所特别指出的结构来实现和获得。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步的详细描述，其中：

图1是本发明算法的流程图；

图2是前后向平滑算法示意图；

图3是传统的延时叠加(DAS)；

图4是最小方差波束形成(MV)；

图5是本发明所提的方法(EIB+FBMV)；

图6是本发明算法生成的点目标成像图中在深度40mm横向截面图与现有算法相比；

图7是本发明算法生成的点目标成像图中在深度50mm横向截面图与现有算法相比；

图8是传统的延时叠加(DAS)；

图9是最小方差波束形成(MV)；

图10是本发明所提的方法。

具体实施方式

以下将结合附图，对本发明的优选实施例进行详细的描述；应当理解，优选实施例仅为了说明本发明，而不是为了限制本发明的保护范围。

图1是本发明算法的流程图，如图所示：本发明提供的一种基于特征空间的前后向自适应波束形成方法，包括以下步骤：

S1：假设有M个阵元的传感器阵列，接收传感器阵列的采样信号；对采样信号进行聚焦延时处理，得到聚焦延时后的信号x(k)；所述信息x(k)表示为x(k)＝[x₁(k)，x₂(k)，...，x_M(k)]，其中，M表示传感器阵列的阵元个数。

S2：对聚焦后的采样信号进行前后向平滑处理，得到样本协方差矩阵估计；图2是前后向平滑算法示意图；给出了前后向空间平滑算法示意图。

具体包括以下步骤：

S21：把M个阵元的数据分为阵元数目为L的子阵，本发明的实施例中子阵阵元数目L的取值为M/2，L取值的上限为M/2，当L＝M/2时，图像的分辨率最高，稳健性较差。考虑到本发明已采用前后向平滑滤波方法提高算法的稳健性，因此取L＝M/2。分别计算各个子阵的相关矩阵，然后根据以下公式来计算，得到前向协方差矩阵估计为

$\tilde{R}\left(k\right)=\frac{1}{M-L+1}\underset{l=1}{\overset{M-L+1}{\mathrm{\Σ}}}{x}_d^l\left(k\right)x_d^l{\left(k\right)}^H,$

公式中， $x_d^l\left(k\right)={[x_d^l\left(k\right),x_d^{l+1}\left(k\right),\·\·\·,x_d^{l+L+1}\left(k\right)]}^T$ 为第l个子阵的输出向量，表示的共轭转置；

S22：定义后向重叠向量为 ${\tilde{x}}_d^l\left(k\right)={[x_d^{M-l+1}\left(k\right),x_d^{M-l}\left(k\right),\·\·\·,x_d^{M-l+L+2}\left(k\right)]}^H,$ l＝1，2，...，M。通过以下公式计算，得到后向协方差矩阵估计

${\tilde{R}}^b\left(k\right)=\frac{1}{M-L+1}\underset{l=1}{\overset{M-L+1}{\mathrm{\Σ}}}{\tilde{x}}_d^l\left(k\right){\tilde{x}}_d^l{\left(k\right)}^H,$

其中，表示第l个子阵的后向输出向量，表示的共轭转置；

S23：通过以下公式来计算前后向协方差矩阵估计：

${\tilde{R}}_{\mathrm{FB}}\left(k\right)=\frac{1}{2}(\tilde{R}\left(k\right)+{\tilde{R}}^b\left(k\right)),$

其中，表示前向协方差矩阵估计与后向协方差矩阵估计的求和平均。

S3：考虑到当协方差矩阵估计是奇异矩阵时，计算自适应波束形成权值可能会造成误差。因此，采用对角加载对协方差矩阵进行处理，用代替 对样本协方差矩阵估计进行对角加载，得到空间协方差矩阵估计，通过下列公式对前后向协方差矩阵进行对角加载得到对角加载后的前后向协方差矩阵估计

${\tilde{R}}_{\mathrm{DL}}={\tilde{R}}_{\mathrm{FB}}+\mathrm{\γI},$

其中，Δ为算法加入的空间噪声与信号功率比，Δ＝1/100L为一恒定常数，为信号等效功率，I表示单位矩阵。

本发明实施例中的对角加载因子γ中Δ的值取为1/100L。

S4：将空间协方差矩阵估计与方向向量结合，计算得到自适应波束形成权值；通过以下公式来计算自适应波束形成权值：

$w=\frac{{\tilde{R}}_{\mathrm{DL}}^{-1}a}{a^H{\tilde{R}}_{\mathrm{DL}}^{-1}a},$

公式中，a为方向向量，w为自适应波束形成权值，表示的逆矩阵。本发明实施例中的方向向量a为全1的矢量。对于超声成像中的宽带信号，方向向量a不能再简单的表示为一组与信号频率有关的复数。为了解决这个问题，一般设置在超声发射阶段采用定点聚焦，接收阶段采用动态聚焦，这样来自焦点的响应可以近似成一个平面波入射阵元，从而得到一个与频率无关的方向向量，因此，本专利中方向向量a表示为全1的矢量。

S5：对空间协方差矩阵估计进行特征分解，构造信号子空间，通过以下公式对进行特征分解可得：

${\tilde{R}}_{\mathrm{DL}}=\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{e}_i{e}_i^H=E_s{\mathrm{\Λ}}_s{E}_s^H+E_n{\mathrm{\Λ}}_n{E}_n^H,$

其中，λ₁≥λ₂≥…≥λ_L是相应的L个特征值，其对应的特征向量分别为e_i，i＝1，2，...，L；用大于最大特征值0.5倍的特征值所对应的特征向量组成信号子空间E_s，而其余特征值对应的特征向量组成噪声子空间E_n，表示E_s的共轭转置，表示E_n的共轭转置，num为信号子空间中特征向量的个数，则Λ_s＝diag{λ₁，λ₂，…，λ_num}，Λ_n＝diag{λ_num+1，…，λ_L}，E_s＝[e₁，e₂，…，e_num]，E_n＝[e_num+1，…，e_L]。

信号子空间中特征向量个数num的选取直接决定了保持主瓣信号和降低旁瓣等级的能力。反之，特征向量个数的选取也与主瓣信号和旁瓣信号有关，会随着环境的变化而变化。主瓣信号的能量主要集中在较大特征值所对应的特征向量，旁瓣信号集中在小特征值所对应的特征向量。因此，一般用大于最大特征值δ倍的特征值所对应的特征向量组成信号子空间E_s。δ在0到1之间取值，本发明专利中δ取0.5。

S6：将自适应波束形成权值投影到信号子空间中，得到新的自适应波束形成权值；因为考虑到在理想情况下，信号的自适应波束形成权值应该位于信号子空间中，存在于噪声子空间中的分量应该为零；所述新的自适应波束形成权值通过以下公式来计算：

$w_1=E_s{E}_s^{H}w,$

其中，w₁表示将自适应波束形成权值投影到由协方差矩阵特征空间构造的信号子空间中后的自适应波束形成权值。

S7：将投影后的自适应波束形成权值对所述的聚焦延时后的采样信号通过以下公式进行加权求和，计算得到自适应波束信号：

$y\left(k\right)=\frac{1}{M-L+1}\underset{l=1}{\overset{M-L+1}{\mathrm{\Σ}}}{w_1}^H\left(k\right)x_d^l\left(k\right),$

其中，y(k)表示计算得到的自适应波束信号，w₁ ^H表示w₁的共轭转置，表示第l个子阵的输出向量。

为了验证该算法的有效性，利用Field II对医学成像中常用的点散射目标和斑散射目标进行成像。Field II是基于线性系统空间响应的原理，它的仿真效果与实际成像很接近，已被国际上广泛认同为仿真超声系统的标准。点散射目标的目标散射点共11个，分布在深度为30～60mm，宽度为10mm的区域内；斑散射目标的成像深度为32～45mm，散射目标中心是位于36mm深度，半径为3mm的圆形吸声体，散射体呈高斯分布。两种仿真模型的仿真参数设置：有效阵元数目M为64，发射信号频率f₀为3.5MHz，系统采样频率f_s为50MHz，阵元中心间隔为半个波长，声速c为1540m/s，子阵长度L为M/2，对角加载量Δ＝1/100L。采用定点发射和动态接收聚焦的模式，信号仿真时加入了60dB的高斯白噪声，成像的动态范围为60dB。

图3至图5给出了本发明算法生成的点目标成像图与现有算法生成的点目标成像图的对比。图3传统的延时叠加(DAS)；图4最小方差波束形成(MV)；图5本发明所提的方法(EIB+FBMV)。图6和图7给出了本发明算法生成的点目标成像图中在深度45mm和50mm处的横向截面图与现有的算法相比。从图3至图7可知，DAS算法分辨率和对比度均为最差；MV算法提高了图像的分辨率，但图像的对比度没有得到很大的改善；而本发明的算法不仅提高了图像的对比度，而且消除了MV算法形成图像中点目标拖着的亮线，算法稳健性得到提高，但是分辨率略低于MV算法。因此，综合考虑图像分辨率、对比度以及稳健性，本发明的方法效果最好。

图8至图10给出了本发明算法生成的斑散射目标成像图与现有的算法相比。图8传统的延时叠加(DAS)；图9最小方差波束形成(MV)；图10发明所提的方法。为了直观的评估波束形成的对比分辨率，引入对比度(CR)，其定义为中心圆的平均功率与外部背景区域的平均功率之差。表1列出各种成像方法的对比度。从表1可知，本发明的方法所得图像的对比度比MV和DAS分别高13.68dB和11.58dB。说明本发明的方法大幅度提高斑散射目标的对比度。

表1吸声斑散射目标体对比度

以上所述仅为本发明的优选实施例，并不用于限制本发明，显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (5)

1.一种用于超声成像的基于特征空间的前后向自适应波束形成方法，其特征在于：包括以下步骤：

S1：接收传感器阵列的采样信号；

S2：对采样信号进行前后向平滑处理，得到样本协方差矩阵估计；

S3：对样本协方差矩阵估计进行对角加载，得到空间协方差矩阵估计；

S4：将空间协方差矩阵估计与方向向量结合，计算得到自适应波束形成权值；

S5：对空间协方差矩阵估计进行特征分解，构造信号子空间；

S6：将自适应波束形成权值投影到信号子空间中，得到新的自适应波束形成权值；

S7：将新的自适应波束形成权值对采样信号数据进行加权求和，得到自适应波束信号；

所述步骤S2中进行前后向平滑处理前还包括对采样信号进行聚焦延时处理，得到聚焦延时后的信号x(k)；所述信号x(k)表示为x(k)＝[x₁(k),x₂(k),…,x_M(k)]，其中，M表示传感器阵列的阵元个数；

所述步骤S2中进行前后向平滑处理，得到样本协方差矩阵估计，具体包括以下步骤：

S21：把M个阵元的数据分为阵元数目为L的子阵，分别计算各个子阵的相关矩阵，然后根据以下公式来计算，得到前向协方差矩阵估计为

$\tilde{R}\left(k\right)=\frac{1}{M-L+1}\underset{l=1}{\overset{M-L+1}{\mathrm{\Σ}}}{x}_d^l\left(k\right)x_d^l{\left(k\right)}^H,$

公式中 $x_d^l\left(k\right)={[x_d^l\left(k\right),x_d^{l+1}\left(k\right),...,x_d^{l+L-1}\left(k\right)]}^T$ 为第l个子阵的输出向量，表示的共轭转置；

S22：定义后向重叠向量为 ${\tilde{x}}_d^l\left(k\right)={[x_d^{M-l+1}\left(k\right),x_d^{M-1}\left(k\right),...,x_d^{M-l-L+2}\left(k\right)]}^H,l=\mathrm{1,2},...,M;$ 通过以下公式计算，得到后向协方差矩阵估计

${\tilde{R}}^b\left(k\right)=\frac{1}{M-L+1}\underset{l=1}{\overset{M-L+1}{\mathrm{\Σ}}}{x}_d^l\left(k\right)x_d^l{\left(k\right)}^H,$

其中，表示第l个子阵的后向输出向量，表示的共轭转置；

S23：通过以下公式来计算前后向协方差矩阵估计：

${\tilde{R}}_{\mathrm{FB}}\left(k\right)=\frac{1}{2}(\tilde{R}\left(k\right)+{\tilde{R}}^b\left(k\right)),$

其中，表示前向协方差矩阵估计与后向协方差矩阵估计的求和平均；

所述步骤S3中，通过下列公式对前后向协方差矩阵进行对角加载得到对角加载后的前后向协方差矩阵估计

${\tilde{R}}_{\mathrm{DL}}={\tilde{R}}_{\mathrm{FB}}+\mathrm{\γI},$

其中，Δ为算法加入的空间噪声与信号功率比，为信号等效功率，I表示单位矩阵；

所述子阵阵元数目L的取值为M/2；

所述对角加载因子γ中Δ的值取为1/100L。

2.根据权利要求1所述的用于超声成像的基于特征空间的前后向自适应波束形成方法，其特征在于：所述步骤S4中，通过以下公式来计算自适应波束形成权值：

$w=\frac{{\tilde{R}}_{\mathrm{DL}}^{-1}a}{a^H{\tilde{R}}_{\mathrm{DL}}^{-1}a},$

公式中，a为方向向量，w为自适应波束形成权值，为的逆矩阵。

3.根据权利要求2所述的用于超声成像的基于特征空间的前后向自适应波束形成方法，其特征在于：所述步骤S5中，通过以下公式对进行特征分解可得：

${\tilde{R}}_{\mathrm{DL}}=\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{e}_i{e}_i^H=E_s{\mathrm{\Λ}}_s{E}_s^H+E_n{\mathrm{\Λ}}_n{E}_n^H,$

其中，λ₁≥λ₂≥…≥λ_L是相应的L个特征值，其对应的特征向量分别为e_i,i＝1,2,…,L；用大于最大特征值0.5倍的特征值所对应的特征向量组成信号子空间E_s，而其余特征值对应的特征向量组成噪声子空间E_n，表示E_s的共轭转置，表示E_n的共轭转置，num为信号子空间中特征向量的个数，则Λ_s＝diag{λ₁，λ₂，…，λ_num}，Λ_n＝diag{λ_num+1,…，λ_L}，E_s＝[e₁,e₂,…，e_num]，E_n＝[e_num+1,…，e_L]。

4.根据权利要求3所述的用于超声成像的基于特征空间的前后向自适应波束形成方法，其特征在于：所述步骤S6中新的自适应波束形成权值通过以下公式来计算：

$w_1=E_s{E}_s^{H}w,$

其中，w₁表示将自适应波束形成权值投影到由协方差矩阵特征空间构造的信号子空间E_s后的自适应波束形成权值。

5.根据权利要求4所述的用于超声成像的基于特征空间的前后向自适应波束形成方法，其特征在于：所述步骤S7中，将投影后的自适应波束形成权值对所述的聚焦延时后的采样信号通过以下公式进行加权求和，计算得到自适应波束信号：

$y\left(k\right)=\frac{1}{M-L+1}\underset{l=1}{\overset{M-L+1}{\mathrm{\Σ}}}{w}_1^H\left(k\right)x_d^l\left(k\right),$

其中，y(k)表示计算得到的自适应波束信号，w₁ ^H表示w₁的共轭转置，表示第l个子阵的输出向量。