一种应用于超声成像的稳健宽带自适应波束形成方法

摘要

本发明涉及一种应用于超声成像的稳健宽带自适应波束形成方法，包括：将所得到的超声回波数据由时域转换至频域，然后进行频带的划分；计算所述频带划分后的各子频带的协方差矩阵；根据得到的所述协方差矩阵，利用特征空间分解方法将假定的导向矢量投影在信号子空间中；通过凸优化方法估计真实的导向矢量；利用最小方差方法的基本原理，得到各子频带空间的波束形成处理结果；对所述各子频带空间的波束形成处理结果进行叠加，并转换到时域，得到用于最终成像的输出结果。本发明充分利用超声回波信号的宽带性特点，在频域内对真实环境下的导向矢量进行重构，使其能在相位畸变的条件下具备较为稳健的性能，并同时提高非均匀介质下的超声成像质量。

说明书

技术领域

本发明涉及医学超声诊断及超声无损检测等领域，尤其涉及到一种应用于超声成像的稳健宽带自适应波束形成方法。

背景技术

传统的医学超声诊断及超声无损检测系统中普遍使用了延时叠加的波束形成方法，这种方法简单可靠，易于实现，但性能却比较一般。为提高超声成像的质量，最小方差自适应波束形成方法受到了广泛的关注，然而这种方法面临如下的两个主要问题：(1)经典的最小方差波束形成方法是根据窄带信号的模型所提出的，而超声成像领域的回波信号一般属于宽带信号的范畴，所以直接将最小方差波束形成方法应用于超声成像中并未能充分展现出原方法的优势；(2)对于非均匀介质的成像过程中往往存在着相位畸变，导致最小方差波束形成方法中的协方差矩阵和导向矢量的估计存在一定误差，继而影响了该方法的稳健性，难以在较为复杂的应用环境中实现理想的性能。

目前存在的稳健自适应波束形成方法，大都围绕在协方差矩阵的重构及导向矢量的估计上，如在最小方差波束形成方法中，常使用对角加载技术来提高协方差矩阵估计的稳健性，这种方法简单易行，而对角加载因子这一参数的选取则对性能有着重要的影响。采用前后向相结合的空间平滑方式对超声回波信号去相干，无需使用对角加载即可在一定程度上保证协方差矩阵估计的稳健性，同时也能够提高最终成像的对比度。但上述方法都集中在时域内处理，并未充分考虑到超声回波信号的宽带性特点，也未能对相位畸变条件下失配的导向矢量进行有效的估计，导致最终的性能提升较为有限。

发明内容

本发明针对现有技术的不足，对最小方差自适应波束形成方法进行改进，充分利用超声回波信号的宽带性特点，在频域内对真实环境下的导向矢量进行重构，使其能在相位畸变的条件下具备较为稳健的性能，并同时提高非均匀介质下的超声成像质量。

为实现上述目的，本发明提供了一种应用于超声成像的稳健宽带自适应波束形成方法，应用于等间距的M个阵元的超声换能器阵列中，该方法主要使用在对接收的回波数据进行波束形成处理的阶段，包括以下的步骤：

将所得到的超声回波数据由时域转换至频域，然后进行频带的划分；

计算所述频带划分后的各子频带的协方差矩阵；

根据得到的所述协方差矩阵，利用特征空间分解方法将假定的导向矢量投影在信号子空间中；

通过凸优化方法估计真实的导向矢量；

利用最小方差方法的基本原理，得到各子频带空间的波束形成处理结果；

对所述各子频带空间的波束形成处理结果进行叠加，并转换到时域，得到用于最终成像的输出结果。

进一步地，所述将所得到的超声回波数据由时域转换至频域，然后进行频带的划分步骤，包括：

对B模式超声成像下每一根扫描线相应的各个阵元所接收到回波数据进行短时傅立叶变换，将变换后所得到的频域数据划分为K个子频带。

更进一步地，所述将所得到的超声回波数据由时域转换至频域，然后进行频带的划分步骤包括：

根据发射的脉冲宽度预先设计出所述短时傅立叶变换的窗函数的宽度，然后进行回波数据空间的划分；

对所述回波数据空间划分后的各段数据进行“补零”处理，然后由时域转换至频域；

对所述频域数据进行频带的划分。

进一步地，所述计算所述频带划分后的各子频带的协方差矩阵步骤，包括：

对各子频带的数据进行空间平滑处理，其中以L个阵元为一组来计算协方差矩阵，共分为M-L+1组，计算的公式如下：

$>\hat{R}\left({\omega }_i\right)=\frac{1}{M-L+1}\underset{l=1}{\overset{M-L+1}{\Sigma }}{x}_d^l\left({\omega }_i\right)x_d^l{\left({\omega }_i\right)}^H$>

其中，是协方差矩阵，是频域上的回波数据，是的共轭转置，i＝1,2，…，K，K为子频带的个数。

更进一步地，所述利用特征空间分解方法将假定的导向矢量投影在信号子空间中的步骤，包括：

按下式对得到的协方差矩阵在特征空间进行分解：

$>\hat{R}={\mathrm{U\Lambda U}}^H=U_s{\Lambda }_s{U_s}^H+U_p{\Lambda }_p{U_p}^H$>

其中Λ＝diag[λ₁,λ₂,...,λ_num,...,λ_L]，且特征值(λ₁≥λ₂≥...λ_num≥...≥λ_L)，U＝[u₁,u₂,...,u_num,...u_L]是与之相对应的特征矢量，选取参数β介于0.1～0.5，满足λ_num≥βλ₁；

将整个特征空间划分为信号子空间Us及噪声子空间Up两部分。

利用得到的信号子空间Us，计算假定的导向矢量在信号子空间内的投影：

$>a=U_s{U_s}^H\bar{a}$>

其中，Us为信号子空间，U_s^H为Us的共轭转置，假定的导向矢量中的元素均取为“1”。

更进一步地，所述通过凸优化方法估计真实的导向矢量，包括：

利用凸优化方法估计真实情况下的导向矢量a，计算公式如下：：

$>\underset{a}{\min }{a}^H{\hat{R}}^{-1}a$>

$>\left(\begin{array}{cc}s.t.& \left|\right|a-\tilde{a}\left|\right|\le ϵ\end{array}\right)$>

其中，是协方差矩阵的逆，导向矢量是假定导向矢量在信号子空间内的投影，参数ε是所设定的误差范围，一般介于之间。

更进一步地，所述利用最小方差方法的基本原理，得到各子频带空间的波束形成处理结果，包括：

根据得到的协方差矩阵及真实情况下的导向矢量，利用最小方差波束形成的基本原理计算出一组最优权向量；

具体地，按如下公式进行计算：

$>w_{MV}\left({\omega }_i\right)=\frac{{\hat{R}}^{-1}\left({\omega }_i\right)a}{a^H{\hat{R}}^{-1}\left({\omega }_i\right)a}$>

其中w_MV(ω_i)是最优向量，是得到的协方差矩阵的逆，导向矢量a是估计的真实情况下的结果；

根据得到的最优权向量以及如下公式计算出各子频带数据的波束形成的输出结果：

$>y\left({\omega }_i\right)=w_{MV}^H\left({\omega }_i\right)·x_d\left({\omega }_i\right)$>

其中，是w_MV(ω_i)的共轭转置，y(ω_i)是第i个子频带的波束形成的输出结果，x_d(ω_i)是频域上的回波数据。

进一步地，所述对所述各子频带数据的波束形成处理结果进行叠加，得到频域空间的波束形成输出结果步骤，包括：

对每一段数据中的各子频带空间的波束形成的处理结果进行累加，得到频域空间的波束形成输出结果，计算的公式如下：

$>Y\left(\omega \right)=\underset{i=1}{\overset{K}{\Sigma }}y\left({\omega }_i\right)$>

其中，Y(ω)是频域空间的波束形成输出结果，y(ω_i)是第i个子频带的波束形成的输出结果。

本发明的优点在于：

本发明通过短时傅立叶变换将时域的回波数据转换到各个子频带中处理，充分利用了超声回波信号的宽带特点，能够取得更为理想的性能；同时在各个频带内将假定的导向矢量投影在信号子空间内，并采用凸优化方法进行重构，提高了算法在相位畸变条件下的稳健性。

附图说明

图1为本发明实施例提供的流程示意图；

图2为本发明实施例提供的对原始时域上的回波数据进行时延预处理的示意图；

图3为本发明实施例提供的得到最终输出结果的原理框图；

图4为本发明实施例提供的各子频带空间进行波束形成的处理过程。

具体实施方式

下面结合附图和实施例，对本发明的技术方案做进一步的详细描述。

图1为本发明实施例提供的一种自适应波束形成方法流程示意图，本发明实施例应用于B模式下的超声成像。等间距的M个阵元的超声换能器阵列工作在B模式成像线性扫描方式下，对所接收到的回波数据进行波束形成处理，每一根扫描线对应着M个阵元的接收信号(即M为阵元数目)，并且每个阵元的采样数目为N，在B模式超声成像下每一根扫描线对应的所有阵元接收的回波数据可以用M*N的矩阵来描述。需要说明的是，这些时域上的回波数据已经经过时延的预处理，图2给出了时延预处理的示意图，以下步骤中即以预处理得到的x_c作为默认的处理对象，同时需要指出的是由于非均匀介质的影响，回波数据并不能保证得到完美的“对齐”，而本发明也无需对回波数据的时延预处理结果做出额外的假设。图3是本发明实施例提供的得到最终输出结果的原理框图，由图1和图3所示，可知本发明实施例的的波束形成处理方法包括步骤101-103：

步骤101，将所得到的超声回波数据由时域转换至频域，然后进行频带的划分；

具体的，在B模式超声成像下每一根扫描线中的所有阵元接收的回波数据可以用M*N的矩阵来描述，其中M为阵元数目，N为采样数目，需要说明的是，该回波数据已经完成了时延预处理，对该回波数据进行短时傅立叶变换，并将变换后所得到的频域数据划分为K个子频带，为在各个子空间内的波束形成完成必要的准备。

下面进一步详细介绍所述步骤101的处理方法：

步骤1011，根据发射的脉冲宽度预先设计出所述短时傅立叶变换的窗函数的宽度，然后进行回波数据空间的划分；

具体地，根据发射的脉冲宽度设计出合适的窗函数的宽度Nl，将整个数据空间划分为Ng段；其中，N为采样数目，Ng＝N/Nl并取其整数部分，即每一段对应的数据空间大小为M×Nl。需要进一步说明的是，为保证最终成像的纵向分辨率，N1一般不会超过发射脉冲的宽度。

步骤1012，对所述回波数据空间划分后的各段数据进行“补零”处理，然后由时域转换至频域；

具体地，在每个阵元所对应的时域上的回波数据后添加(Nf-Nl)个“0”，然后对每一段数据进行傅里叶变换处理，即将时域上的回波数据转换到频域，经过短时傅里叶变换后的回波数据将有Ng段；

在处理过程中，为了短时傅立叶变换的便捷计算，Nf一般取2的指数幂，并满足Nf>2Nl，从而每一段数据空间的大小变为M*Nf，即变换后得到的每一段频域数据都是M*Nf的矩阵；

步骤1013，对所述频域数据进行频带的划分；

具体地，上述步骤1012共得到Ng段M*Nf的频域数据，对得到的频域数据，进行频带划分，将每一段频带划分为K个子频带，其中，Nl<K<Nf/2。

步骤102，对所述频带划分后的各子频带空间进行波束形成处理，得到所述各子频带空间的波束形成处理结果；

具体地，如图4所示，在所划分的各个子频带中，在获得各子频带的协方差矩阵后，先通过特征空间分解的方法将假定的导向矢量投影在信号子空间内，然后通过凸优化方法估计实际情形下的导向矢量，继而利用最小方差方法的基本原理，得到对各个子频带空间波束形成输出结果。

下面以第i(i＝1，2，…，K)个子频带为例，详细介绍步骤102的处理过程：

步骤1021，计算各子频带的协方差矩阵；

具体地，第i个子频带所对应的数据大小为M×1，首先对维度是M×1的数据进行空间平滑处理，其中以L个阵元为一组来计算协方差矩阵，共分为M-L+1组，计算的公式如下：

$>\hat{R}\left({\omega }_i\right)=\frac{1}{M-L+1}\underset{l=1}{\overset{M-L+1}{\Sigma }}{x}_d^l\left({\omega }_i\right)x_d^l{\left({\omega }_i\right)}^H$>

其中，是协方差矩阵，是频域上的回波数据，是的共轭转置，i＝1,2，…，K，K为子频带的个数。

步骤1022，对步骤1021中得到的协方差矩阵在特征空间进行分解；

具体地，按照如下的公式进行分解：

$>\hat{R}={\mathrm{U\Lambda U}}^H=U_s{\Lambda }_s{U_s}^H+U_p{\Lambda }_p{U_p}^H$>

其中Λ＝diag[λ₁,λ₂,...,λ_num,...,λ_L]，且特征值按照由大到小的顺序排列，即(λ₁≥λ₂≥...λ_num≥...≥λ_L)，而U＝[u₁,u₂,...,u_num,...u_L]则是与之相对应的特征矢量，选取合适的参数β，满足λ_num≥βλ₁，继而将整个特征空间划分为信号子空间Us及噪声子空间Up两部分。为保证算法的稳健性能，选择的参数β介于0.1～0.5；

步骤1023，利用步骤1022得到的信号分量特征空间，计算假定的导向矢量在信号子空间内的投影；

具体地，按如下公式进行计算：

$>a=U_s{U_s}^H\bar{a}$>

其中，Us为信号子空间，U_s^H为Us的共轭转置，假定的导向矢量中的元素均取为“1”；

步骤1024，利用凸优化方法估计真实情况下的导向矢量a；

具体地，按如下公式进行计算：

$>\underset{a}{\min }{a}^H{\hat{R}}^{-1}a$>

$>\left(\begin{array}{cc}s.t.& \left|\right|a-\tilde{a}\left|\right|\le ϵ\end{array}\right)$>

其中是由步骤1021中所得到的协方差矩阵的逆，导向矢量是步骤1023中得到的假定导向矢量在信号子空间内的投影，参数ε是所设定的误差范围，其取值一般介于之间。

步骤1025，根据得到的协方差矩阵及重构的导向矢量，利用最小方差波束形成的基本原理计算出一组最优权向量；

具体地，按如下公式进行计算：

$>w_{MV}\left({\omega }_i\right)=\frac{{\hat{R}}^{-1}\left({\omega }_i\right)a}{a^H{\hat{R}}^{-1}\left({\omega }_i\right)a}$>

其中w_MV(ω_i)是最优权向量，是步骤1021中得到的协方差矩阵的逆，导向矢量a是步骤1024中得到的估计结果。

步骤1026，根据得到的最优权向量以及如下公式计算出第i个子频带空间内的波束形成的输出结果：

$>y\left({\omega }_i\right)=w_{MV}^H\left({\omega }_i\right)·x_d\left({\omega }_i\right)$>

其中，是w_MV(ω_i)的共轭转置，y(ω_i)是第i个子频带的波束形成的输出结果，x_d(ω_i)是频域上的回波数据。

步骤1027，按照上述的方法依次计算出B模式成像下每一根扫描线中各子频带(i＝1，2，3，…，K)的波束形成输出结果；

步骤103，将各子频带空间得到的波束形成输出结果进行累加，并转换到时域，得到用于最终成像的输出结果。

具体地，所述步骤103的处理方法包括：

步骤1031，将每一段数据中的K个子频带的波束形成的输出结果进行累加，得到频域的波束形成输出结果，计算公式如下：

$>Y\left(\omega \right)=\underset{i=1}{\overset{K}{\Sigma }}y\left({\omega }_i\right)$>

其中，Y(ω)是频域空间的波束形成输出结果，y(ω_i)是第i个子频带的波束形成的输出结果。

得到频域空间的波束形成输出结果，还需要进行后续的处理，用于最终的显示成像：

步骤1032，将频域空间的波束形成输出结果转换到时域；

具体地，将每一段数据中的频域空间的波束形成输出结果进行傅立叶逆变换，即将频域空间的波束形成输出结果转换至时域，得到最终的输出结果，在处理过程中，将每段数据的大小取为M×Nl。

步骤1033，重复上述处理，计算得到所有Ng段的输出结果，使得每一根扫描线对应的波束形成输出的数据大小为1×N，用于最终的显示成像。

本发明涉及一种应用于超声成像的稳健宽带自适应波束形成方法，该方法的要点包括：利用最小方差波束形成方法的基本思想，结合短时傅立叶变换、特征空间分解以及导向矢量估计方法来提高波束形成器的性能。首先将所得到的回波信号由时域转换至频域，并根据超声信号的频带特点将其划分为若干子带，接着在各个子带上分别做波束形成处理，即通过特征空间分解的方法将假定的导向矢量投影在信号子空间内，然后通过凸优化方法估计实际情形下的导向矢量，继而利用最小方差方法的基本原理，得到对各个子频带空间的处理结果，并进行叠加得到整个频域的波束形成输出结果，最终利用傅立叶逆变换转为时域上的数据用于最后的显示成像。本发明能够提高在非均匀介质环境下波束形成方法的稳健性，并同时显著改善了超声成像的质量，可以应用于医学超声诊断及超声无损检测等实际需求中。

专业人员应该还可以进一步意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

结合本文中所公开的实施例描述的方法或算法的步骤可以用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器(RAM)、内存、只读存储器(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。