工程科学与技术   2021, Vol. 53 Issue (4): 200-208
空域抽样与相干因子融合的超声阵列自适应波束形成算法
宋寿鹏, 李棋     
江苏大学 机械工程学院,江苏 镇江 212013
基金项目: 国家自然科学基金项目(51375217)
摘要: 针对超声成像中自适应波束形成算法计算效率低的问题,提出了一种基于空域抽样与相干因子融合的自适应波束形成算法。该算法通过分析阵列数据波束图,推导出不同阵元数目下的最大抽取因子,将阵列中所有阵元接收的数据按最大抽取因子进行空域等间隔抽样,得到阵元稀疏后的回波数据,减少了采集数据量。将空域抽样数据输入波束形成器,计算其协方差矩阵,再将协方差矩阵构造为Toeplitz矩阵,结合最小方差原理从构造的Toeplitz矩阵计算出抽样数据的自适应加权值,并引入相干因子对自适应加权值进行修正,以突出抽样数据中的有效信息,抑制干扰成份。在使用不对等数据信息和使用空域抽样数据成像的情况下,采用该算法、最小方差算法和融合相干因子的最小方差算法分别对裂纹和横通孔缺陷仿真成像。结果表明:在数据信息不对等时,在成像质量上,该算法介于另外两种算法之间;在成像时间上,与另外两种算法相比,该算法的成像时间平均减少85%以上。在相同的空域抽样数据下,在成像质量上,该算法成像效果优于其他算法;在成像时间上,与另外两种算法相比,该算法的成像时间平均减少65%以上。
关键词: 波束形成算法    空域抽样    相干因子    计算复杂度    
Adaptive Beamforming Algorithm for Ultrasonic Array with the Combination of Spatial Sampling and Coherence Factor
SONG Shoupeng, LI Qi     
School of Mechanical Eng., Jiangsu Univ., Zhenjiang 212013, China
Abstract: In order to solve the problem of low computational efficiency of adaptive beamforming algorithms in ultrasonic imaging, an adaptive beamforming algorithm for ultrasonic array with the combination of spatial sampling and coherence factor was proposed. The maximum decimation factor with different numbers of array elements was deduced according to the beam pattern. The sparse echo data was obtained by spatially sampling the whole array element data using the maximum decimation factor. Therefore, the amount of data used for beamforming was greatly reduced. Taking the spatial sampling data as the input of a beamformer and constructing the covariance matrix as Toeplitz matrix, the adaptive weights of the sampling data were obtained according to the principle of minimum variance. Then, the adaptive weights were modified by introducing the coherence factor to highlight the effective information of the sampling data. Under the case of unequal data and spatial sampling data, the proposed algorithm, minimum variance algorithm and minimum variance algorithm combined with coherence factor were used to simulate the imaging of cracks and cross-drilled holes respectively. The results show that: for unequal data, the imaging quality of the proposed algorithm is between the other two algorithms; in terms of imaging time, compared with the other two algorithms, the average imaging time of the proposed algorithm is reduced by more than 85%. For the same spatial sampling data, the imaging quality of the proposed method is better than the other algorithms; in terms of imaging time, compared with the other two algorithms, the average imaging time of the proposed algorithm is reduced by more than 65%.
Key words: beamforming algorithm    spatial sampling    coherence factor    computational complexity    

自适应波束形成技术是超声成像的关键技术之一,是根据阵列中阵元接收的信号动态地计算权值并对信号进行加权以控制波束进行成像,直接决定着成像的质量[1-3]

自适应算法中典型的最小方差(minimum variance,MV)算法可以在保持期望信号能量的同时,使波束形成器输出的信号功率最小化[4-6]。为了进一步提高MV算法的成像性能,不少学者将波束形成算法与相干因子(coherence factor,CF)相融合,以达到提高成像分辨率的目的[7-9]。Wang等[10]提出将CF与最小方差无畸变(minimum variance distortionless response,MVDR)算法结合应用于高帧率成像的算法,提高成像分辨率和对比度以降低旁瓣等级。在此基础上,Asl等[11]指出该方法可以有效减小聚焦误差,在声速不均匀时表现出良好的鲁棒性。Xu等[12]提出了MV算法与CF和相位相干因子(phase coherence factor,PCF)结合的方法,获得了高对比度和高分辨率的医学超声图像。虽然这些改进和研究可以有效提升MV算法的成像分辨率,但是在计算复杂度上却没有改善。对于一个阵元数为N的线性阵列,自适应MV算法的计算复杂度高达 $ O\left( {{N^3}} \right)$ ,严重影响了超声成像的时效性。

针对MV算法中由于协方差矩阵反演造成的计算复杂度高且计算效率低的问题,不少学者在这一领域开展了研究。Park等[13]提出一种结合FPGA的波束形成方法,通过减少计算量达到降低成像计算复杂度的目的,加快了成像速度。Chen等[14]提出一种基于坐标旋转数字计算机(coordinate rotation digital computer,CORDIC)处理器的低复杂度、高吞吐量的波束形成算法。这些研究虽然有效提高了超声成像效率,但本质上是将软件算法与硬件结合,借助硬件提高成像速度。一些学者进一步改进MV算法中的协方差矩阵,通过降低矩阵反演的计算复杂度等级来提高成像速度。Kim等[15]利用主成分分析(principal component analysis,PCA)实现降维,通过预先计算的常规MV权值离线估计主成分并通过选定的主成分的线性组合逼近MV权值。Park等[16]提出运用QR分解技术将空间协方差矩阵转换为标量矩阵,无需对矩阵反演即可得到自适应加权值和波束形成器的输出。这些方法虽然可以使协方差矩阵反演的计算复杂度降低,但是计算得到的加权值是近似值并非实际值。Nilsen等[17]研究了波束空间(beam space,BS)域中的MV波束形成方法,把数据从阵列空间转换至波束空间,这一操作可将矩阵反演计算复杂度等级由3维降至2维,但存在方向矢量误差时,该方法的性能可能会下降。Asl等[18]通过对采样协方差矩阵作近似,将其为Toeplitz矩阵后再进行矩阵反演。虽然上述方法中矩阵反演的计算复杂度等级降至2维,但由于参与运算的数据量大,因此计算复杂度仍然很高。为此,一些学者开展了减少采样阵元数目的研究以减少处理数据量。

Sakhaei[19]描述了一种抽样MV波束形成器算法,首先结合接收波束的波束模式的分析对全部数据抽样,然后利用全部数据计算出加权系数并对抽样数据进行加权,最后在医学成像上验证了其可行性。在此基础上,Shamsian等[20]提出级联结构的快速波束形成方法,将低通滤波器与最小方差波束形成器级联。其中:第1阶段,通过低通滤波器去除回波信号中的离轴噪声,再对数据进行抽样;第2阶段,将抽样数据作为MV算法的输入以抑制轴上干扰,并提出抽取子波束的MV算法(decimated sub-beam MV,DSMV)。该方法实现了数据量和计算复杂度上的降低,但对阵元数据进行抽样这一过程不可避免会丢失部分信息,进而影响成像质量。

为了在保证成像分辨率的同时,降低算法计算复杂度进而提高计算效率,本文提出了一种空域抽样与相干因子融合的超声阵列自适应波束形成算法(decimated minimum variance combined with coherence factor,DCFMV)。该算法先利用数据空域抽样以减少数据量,然后基于MV原则计算自适应加权值,并融合相干因子对期望信号进行增强。为了进一步简化运算,在对数据处理时,将数据协方差矩阵构造为Toeplitz矩阵的结构。

1 空域抽样与相干因子融合的自适应波束形成算法 1.1 阵列回波信号的空域抽样

由于阵列中阵元接收通道多,数据量大,导致超声波束形成算法计算复杂度高,成像时效性差。为了减少采集数据,将阵列数据进行空域抽样,空域抽样首先需要确定抽取因子(或称抽取间隔) $D$ ,按 $D$ 大小对阵列信号空域等间隔抽样,得到阵列空域抽样数据。此时,阵列数量由 $N$ 变为 $\left\lfloor {N{\rm{/}}D} \right\rfloor $ ,其中, $\left\lfloor \cdot \right\rfloor $ 表示向下取整, $N$ 为阵元数量。如果将空域抽样前的数据称为阵列完备数据,则经空域抽样后的数据即为非完备数据。抽取因子的选取会影响成像效果,若抽样不当,会导致波束图上出现栅瓣,分散主瓣的能量,在成像时会形成伪像。以16阵元线性阵列为例,当抽取因子 $D{\rm{ = }}2$ 时,其完备数据与非完备数据的波束图如图1所示,由图1可以看出,数据抽取会导致波束图中出现额外波瓣,这些波瓣称为栅瓣,能量会从栅瓣泄露。

图1 阵列波束对比 Fig. 1 Contrast of the array beam

对完备阵列数据进行空域抽样,那么其接收波束模式 $B(u)$ 可以表示为:

$B(u) = \frac{{\sin \left( {\dfrac{{{\text{π}}uN}}{{2D}}} \right)}}{{\sin \left( {\dfrac{{{\text{π}}u}}{2}} \right)}}$ (1)

式中, $u$ 为空间归一化频率,取值位于区间 $\left[ {0,1} \right]$ 。当 $u = 2KD{\rm{/}}N$ $K$ 为介于0与 $\left\lfloor {{N/{2D}}} \right\rfloor $ 间的任意正整数,此时 $\sin \left(\dfrac{{{\text{π}}uN}}{{2D}}\right) = 0$ ,第1个零值点 ${u_1} = 2D{\rm{/}}N$

若抽取因子 $D$ 取值较小,表明阵列数据信息丢失少,成像效果好,但数据量减少不明显,成像时效性也得不到明显提升;若抽取因子 $D$ 取值较大,则阵列数据丢失的信息较多,虽然成像时效性可以明显提高,但牺牲了成像质量。因此,为解决成像质量与时效性的矛盾,在空域数据抽取中选择合适的抽取因子 $D$ 至关重要。

如果抽取的非完备数据通过栅瓣泄露的功率与完备数据的总功率相比可以忽略不计,那么,最大抽取因子 ${D_{\max }}$ 近似满足[20]

${D_{\max }} = \frac{1}{{{u_1}}}$ (2)

${u_1}$ 代入式(2)可以得到:

$2D \cdot {D_{\max }}{\rm{ = }}N$ (3)

考虑到抽取因子 $D$ 和最大抽取因子 ${D_{\max }}$ 均为正数,则存在:

$N \le {D^2}{\rm{ + }}D_{\max }^2$ (4)

由于 ${D^2}{\rm{ + }}D_{\max }^2 \le 2D_{\max }^2$ ,则:

${D_{\max }} \ge \sqrt {\frac{N}{2}} $ (5)

考虑到阵元抽取间隔只能取整数阵元,因此 $\sqrt {{N / 2}} $ 向下取整,抽取因子 $D$ 可选取区间 $\left[ {1,{\rm{fix}}\left( {\sqrt {{N / 2}} } \right)} \right]$ 上的任一整数,一般情况下选取最大整数,其中, ${\rm{fix}}\left( \cdot \right)$ 表示向零方向取整。

确定抽取因子 $D$ 之后,对接收阵元按抽取因子进行空域等间隔抽样。以 $N{\rm{ = }}16$ $D{\rm{ = }}2$ 为例,空域抽样数据原理示意图如图2所示。由图2可知:发射阵元采用全矩阵发射模式;各阵元发射时,接收阵元按 $D{\rm{ = }}2$ 等间隔抽样,得到空域抽样数据。

图2 空域抽样数据原理示意图 Fig. 2 Schematic diagram of spatial sampling data

1.2 最小方差波束形成算法

得到抽样数据后,采用最小方差波束形成方法使阵列通道所接收到的噪声信号能量和非期望方向上的干扰信号能量达到最小,以此获得最优权值。将空域抽样后的信号作为波束形成器的输入,那么,波束形成器的输出 ${{y}}(k)$ 可以表示为:

${{y}}(k) = {{{w}}^{\rm{H}}}(k){{x}}(k)$ (6)

式中: $k$ 为不同采样点的序列号; ${{w}}$ 为加权列向量; $ {\left[ \cdot \right]^{\rm{H}}} $ 为共轭转置; ${{x}}(k){\rm{ = [}}{x_1}(k{\rm{ - }}{\varDelta _1}),{x_2}(k{\rm{ - }}{\varDelta _2}), \cdots ,{x_i}(k - {\varDelta _i}), \cdots ,$ ${x_M}(k - {\varDelta _M}){]^{\rm{T}}}$ ,为延时后的空域抽样信号,其中, ${x_i}(k)$ 为第 $i$ 个接收阵元在 $k$ 时刻的接收回波, ${\varDelta _i}$ 为第 $i$ 个接收阵元的延时时间, ${\left[ \cdot \right]^{\rm{T}}}$ 为转置,下标 $M$ 为空域抽样后的接收阵元数。

在约束条件为 ${{{w}}^{\rm{H}}}{{a}}{\rm{ = }}1$ 的情况下,通过最小化输出功率即 $\min ({{{w}}^{\rm{H}}}{{Rw}}{\rm{)}}$ 求解最优权向量。利用拉格朗日乘子法求解最优权向量 ${{{w}}_{\rm{opt}}}$ ,得到:

${{{w}}_{{\rm{opt}}}} = \frac{{{{{R}}^{{\rm{ - 1}}}}{{a}}}}{{{{{a}}^{\rm{H}}}{{{R}}^{{\rm{ - 1}}}}{{a}}}}$ (7)

式中: ${{R}}{\rm{ = }}E\left[ {{{x}}(k){{{x}}^{\rm{T}}}(k)} \right]$ $M{\rm{ \times }}M$ 维的空间样本协方差矩阵,其中, $ E[\cdot ]$ 为期望; ${{{R}}^{{\rm{ - 1}}}}$ ${{R}}$ 的逆矩阵; ${{a}}$ 为方向的列向量,一般令 ${{a}}{\rm{ = }}{\left[ {1,1, \cdots ,1} \right]^{\rm{T}}}$

在实际检测中,由于阵列中各阵元接收到的回波信号来自同一区域内的检测对象的反射或散射,因此,信号具有高相关性使得计算所得协方差矩阵 ${{R}}$ 奇异,从而导致加权值 ${{w}}$ 难以求得。为了解决加权值 ${{w}}$ 求解这一难题,需要对回波信号解相关。解相关就是先将协方差矩阵进行划分,分为 $P{\rm{ = }}M{\rm{ - }}L{\rm{ + }}1$ 个重叠子阵,子阵长度 $L{\rm{ < }}M{\rm{/}}2$ ;然后,计算求得每个子阵的协方差矩阵;最后将每个子阵的协方差矩阵相加取平均值。最终得到平滑处理后的协方差矩阵:

${{R}}{\rm{ = }}\frac{1}{P}\sum\limits_{p = 1}^P {{{{x}}_p}(k){{x}}_p^{\rm{H}}(k)} $ (8)

式中, ${{{x}}_p}(k)$ 为第 $p$ 个子阵的接收数据。

为增强MV算法对建模误差的鲁棒性,对协方差矩阵 ${{R}}$ 进行修正,在协方差矩阵中引入一个加载因子 $\varepsilon $ ,也就是用 ${{R}}{\rm{ + }}\varepsilon {{I}}$ 代替 ${{R}}$ 。其中: ${{I}}$ 为单位阵列,从而起到压缩干扰信号的作用; $\varepsilon $ 的取值遵循[21]

$\varepsilon {\rm{ = }}\gamma \cdot {\rm{tr}}\left\{ {{R}} \right\}$ (9)

式中, ${\rm{tr}}\left\{ \cdot \right\}$ 为矩阵的迹, $\gamma $ 为一恒定常数。

回波信号经解相关和协方差矩阵修正后,波束形成器的输出 ${{y}}(k)$ 可改写为:

${{y}}(k) = \frac{1}{P}\sum\limits_{p = 1}^P {{{w}}_p^{\rm{H}}(k){{{x}}_p}(k)} $ (10)

式中,第 $p$ 个子阵在 $k$ 点的加权值 ${{{w}}_p}(k)$ 为:

${\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{w}}_p}(k) = \frac{{{{({{{R}}_p}{\rm{ + }}\gamma \cdot {\rm tr}\left\{ {{{{R}}_p}} \right\}{{I}})}^{{\rm{ - 1}}}}{{a}}}}{{{{{a}}^{\rm{H}}}{{({{{R}}_p}{\rm{ + }}\gamma \cdot {\rm tr}\left\{ {{{{R}}_p}} \right\}{{I}})}^{{\rm{ - 1}}}}{{a}}}}$ (11)

式中, ${{{R}}_p} = E\left[ {{{{x}}_p}(k){{x}}_p^{\rm{T}}(k)} \right]$ 为第 $p$ 个子阵的协方差矩阵。

1.3 相干因子计算

将空域抽样后的数据作为最小方差波束形成算法的输入,可以有效减少算法中的数据处理量,但是,相比于完备数据,此时的输入信号中不可避免地丢失了部分有用信息。

为增强输入信号中的有效信息,引入相干因子。相干因子是通过计算目标点处相干能量与总能量之比,给相干能量较大的目标点赋予更大的权值 ${C_{\rm{f}}}(k)$ ,从而提高超声回波图像中相干能量,通过提高有效信息的区分度改善超声回波图像质量。相干因子的计算公式为[10]

${\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;C_{\rm{f}}}(k){\rm{ = }}\frac{{{S\!_{{\rm{xg}}}}}}{{{S\!_{{\rm{sum}}}}}}{\rm{ = }}\frac{{{{\left| {\displaystyle\sum\limits_{i = 1}^M {{{{x}}_i}(k{\rm{ - }}{\varDelta _i})} } \right|}^2}}}{{M\displaystyle\sum\limits_{i = 1}^M {{{\left| {{{{x}}_i}(k{\rm{ - }}{\varDelta _i})} \right|}^2}} }}$ (12)

式中, ${S\!_{\rm{xg}}}$ 为空域抽样信号中的相干能量, ${S\!_{{\rm{sum}}}}$ 为空域抽样信号的总能量。 ${C_{\rm{f}}}(k)$ 取值在0~1之间,取值越大,表示回波信号中相干能量占比越高。

对于采样点 $k$ ,空域抽样与相干因子融合的波束形成算法的输出修正为:

${\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{y}}(k){\rm{ = }}\frac{{{C_{\rm{f}}}(k)}}{P}\sum\limits_{p = 1}^P {{{w}}_p^{\rm{H}}(k){{{x}}_p}(k)} $ (13)
2 协方差矩阵反演的计算过程与复杂度分析

理想情况下协方差矩阵 ${{R}}$ 具有Toeplitz矩阵的结构,但实际中通常会出现相干源和低信噪比的情况,因此不满足该结构。在对 ${{R}}$ 反演时,为了降低计算复杂度等级将其构造为Toeplitz矩阵,该方法在有效降低求取加权值的计算复杂度等级的同时,也可以使 ${{R}}$ 接近真实的数据协方差矩阵。

协方差矩阵 ${{R}}$ $M \times M$ 维矩阵,将 ${{R}}$ 下三角和上三角部分各条对角线上元素取平均得到Toeplitz矩阵 ${{{R}}_{\rm{T}}}$ ${{{R}}_{\rm{T}}}$ 每一条对角线上的任意两个值都是相等的,任意元素的计算公式如下:

${\;\;\;\;\;\;\;\;\;\;\;\left\{\!\!\!\! {\begin{array}{*{20}{c}} {{r_{{\rm{ - }}q}} = \dfrac{1}{{M - q}}\displaystyle\sum\limits_{i = 1}^{M - q} {{r_{i,(i + q)}}} ,} \\ {{r_q} = \dfrac{1}{{M - q}}\displaystyle\sum\limits_{i = 1}^{M - q} {{r_{(i + q),i}}} ,} \end{array}} \right.0 \le q \le M - 1}$ (14)

式中, ${r_{i,(i{\rm{ + }}q)}}$ 为矩阵 ${{R}}$ 中第 $i$ $(i + q)$ 列的元素, ${r_{(i{\rm{ + }}q),i}}$ 为矩阵 ${{R}}$ 中第 $(i + q)$ $i$ 列的元素, ${r_{{\rm{ - }}q}}$ 为矩阵 ${{{R}}_{\rm{T}}}$ 上三角部分第 $q$ 条次对角线上的任一元素的值, ${r_q}$ 为矩阵 ${{{R}}_{\rm{T}}}$ 下三角部分第 $q$ 条次对角线上的任一元素的值。

递推求出 ${{{R}}_{\rm{T}}}$ 的逆矩阵 ${{R}}_{\rm{T}}^{{\rm{ - 1}}}$ 的各列,具体流程如下:

1)根据Zohar算法[22-23],求解线性方程组:

$\left\{\!\!\!\! {\begin{array}{*{20}{c}} {{{{R}}_{\rm{T}}}{{x}} = {{f}}{\rm{,}}} \\ {{{{R}}_{\rm{T}}}{{y}} = {{{e}}_M}} \end{array}} \right.$ (15)

式中, ${{{e}}_M}{\rm{ = }}{\left( {0,0, \cdots ,1} \right)^{\rm{T}}}$ ${{f}}{\rm{ = }}({r_1}{\rm{ - }}{r_{1{\rm{ - }}M}},$ ${r_2}{\rm{ - }}{r_{2{\rm{ - }}M}}, \cdots ,{r_{M{\rm{ - }}1}}{\rm{ - }}{r_{{\rm{ - }}1}}, $ $ 0{)^{\rm{T}}}$ 。求得方程组(15)的唯一解为 ${{x}}$ ${{y}}$ ,这一步骤共需 $4{M^2}{\rm{ + }}4$ 次乘除运算、 $4{M^2}{\rm{ - }}5M{\rm{ + }}8$ 次加减运算。

2)根据Toeplitz矩阵定理, ${{R}}_{\rm{T}}^{ - 1}$ 的列向量 ${{{z}}_j}(j{\rm{ = }}1, $ $ 2, \cdots ,M)$ 满足如下关系[24]

$\left\{\!\!\!\! {\begin{array}{*{20}{l}} {{{{z}}_M} = {{y}}{\rm{,}}} \\ {{{{z}}_j} = {{Q}}{{{z}}_{j + 1}} + {{{y}}_{M - j}}{{x}} - {{{x}}_{M - j}}{{y}}{\rm{,}}} \end{array}} \right.j = M - 1, \cdots ,2,1$ (16)

式中, $ {{Q}} = \left[ {\begin{array}{*{20}{c}} 0&1& \cdots &0 \\[-2pt] \vdots & \vdots &{}& \vdots \\ 0&0& \cdots &1 \\ 1&0& \cdots &0 \end{array}} \right]$

递推求出矩阵的逆的各列元素,需要 ${M^2}{\rm{ - }}M$ 次乘除运算、 ${M^2}{\rm{ - }}M$ 次加减运算。

综合步骤1)和2),一共需要 $5{M^2}{\rm{ - }}M{\rm{ + 4}}$ 次乘除运算、 $5{M^2}{\rm{ - }}M{\rm{ + }}4$ 次加减运算。因此,对于 $M{\rm{ \times }}M$ 维矩阵 ${{R}}$ ,经过Toeplitz结构化后反演的计算复杂度为 $O({M^2})$ ,即 $O({N^2}{\rm{/}}{D^2})$

将本文提出的空域抽样与相干因子融合的最小方差波束形成方法(decimated minimum variance combined with coherence factor,DCFMV)与传统MV算法及融合相干因子的最小方差波束形成算法[11](minimum variance combined with coherence factor,CFMV)计算协方差矩阵上的计算复杂度进行对比,结果如表1所示。

表1 协方差矩阵反演的计算复杂度对比 Tab. 1 Comparison of computational complexity of covariance matrix inversion

由于抽取因子 $D \ge 1$ ,因此从表1中可直观看出DCFMV方法比MV和CFMV算法计算复杂度明显降低。

3 仿真与讨论

为了研究DCFMV算法的性能,利用Field Ⅱ仿真软件对本文提出的空域抽样与相干因子融合的最小方差波束形成方法(DCFMV)、传统MV算法、融合相干因子的最小方差波束形成算法(CFMV)[11]这3种方法进行成像对比。仿真试验中选用裂纹和横通孔两种缺陷,结构和分布示意图分别如图3(a)(b)所示。试块上设置裂纹尺寸为 $0.2 \times 10 \times 20\;{\rm{ m}}{{\rm{m}}^3}$ ,横通孔为4个直径为2 mm的通孔缺陷,缺陷之间间隔为20 mm,后续分析中从上至下依次命名为孔1、孔2、孔3和孔4。

图3 缺陷分布示意图 Fig. 3 Schematic diagram of the defect distribution

试验中,超声换能器为线性阵列,阵元数 $N = 64$ ,各阵元中心距为0.6 mm,阵元宽度为0.5 mm,中心频率为5 MHz,采样频率为100 MHz,换能器采用垂直入射方式,超声波在被测试件中的声速为5 900 m/s。子阵列长度 $L = N{\rm{/}}2$ ,对角加载系数 $\gamma = 1{\rm{/}}10L$ ,抽取因子 $D = 4$ ,采用动态孔径的成像方法,成像深度与孔径的比值为定值2。设置超声图像动态显示范围为60 dB。使用的计算机配置为Intel(R)Core(TM)i5–8250U,运行内存8 GB。

3.1 数据信息不对等时成像性能分析

在DCFMV、MV、CFMV3种算法下,裂纹和横通孔缺陷成像结果分别如图45所示,其中,MV和CFMV这两种算法采用完备数据,DCFMV算法采用空域抽样后的非完备数据。

图4 数据不对等时裂纹缺陷成像效果对比 Fig. 4 Comparison of imaging effects of crack defects under unequal data

图5 数据不对等时横通孔缺陷成像效果对比 Fig. 5 Comparison of imaging effects of cross-drilled hole defects under unequal data

图4可知,CFMV算法成像时裂纹与实际缺陷尺寸误差较大,成像结果的量化精度没有MV和DCFMV算法准确。

用阵列成像质量性能指标(array performance indicator,API)表征成像分辨率。API是点扩散函数空间大小的无量纲度量,定义为点扩散函数从其最大值向下大于–6 dB的区域 ${A_{{\rm{ - 6\;dB}}}}$ 与波长的平方之比,计算式[25]如下:

${\rm API} = \frac{{{A_{{\rm{ - 6\;dB}}}}}}{{{\lambda ^2}}}$ (17)

图5可知,3种方法中,CFMV的成像效果最好,MV的成像效果最差。通过式(17)计算图5中横通孔缺陷图像的API值,发现CFMV算法API值最小,DCFMV算法次之,MV算法效果最差。对比图5中3种成像算法的主瓣宽度,发现CFMV算法主瓣最窄,DCFMV算法次之,MV算法主瓣最宽,但在数值上相差不明显。

为了对比DCFMV、MV、CFMV这3种算法的计算复杂度,以成像耗时作为衡量指标,分别对3种成像算法的成像时间进行对比,15个样本取平均,结果如表2所示。

表2 数据不对等时算法成像时间对比 Tab. 2 Comparison of imaging time of algorithms under unequal data

表2可知:通过参考15个样本取平均,DCFMV算法成像时间大幅降低;对于裂纹缺陷,DCFMV算法的成像时间相较于MV算法和CFMV算法分别减小了85.12%和85.24%;对于横通孔缺陷,DCFMV算法的成像时间相比于MV算法和CFMV算法分别缩短了86.11%和86.47%。

3.2 使用空域抽样数据成像性能分析

利用空域抽样数据分别对DCFMV、MV、CFMV这3种成像算法进行了对比,其中抽取因子 $D = 4$ ,成像结果分别如图67所示。

图6 空域抽样数据下裂纹缺陷成像效果对比 Fig. 6 Comparison of imaging effects of crack defects under spatial sampling data

图7 空域抽样数据下横通孔缺陷成像效果对比 Fig. 7 Comparison of imaging effects of cross-drilled hole defects under spatial sampling data

图67可知,在3种算法下,MV算法成像效果最差。为了更好地比较3种成像算法的成像质量,以API作为成像质量衡量指标,在对裂纹缺陷成像时,MV算法、CFMV算法及DCFMV算法计算得到的API值分别为3.53、1.20、1.11,DCFMV算法的API值分别比MV算法和CFMV算法降低了68.77%和8.34%,表明在空域抽样数据下DCFMV算法的成像质量好于MV算法和CFMV算法;在对通孔缺陷成像时,DCFMV算法的API值都大幅度降低,4个孔类缺陷从上至下的API值与MV算法相比分别降低了70.91%、61.86%、66.62%、66.88%,与CFMV算法相比分别降低了17.77%、7.74%、30.07%、11.99%。综上所述,DCFMV算法下的图像质量优于其他两种算法。

以通孔类缺陷为例,分别计算采用DCFMV、MV、CFMV这3种算法的4个不同深度孔类缺陷主、旁瓣的分布,结果如图8所示。

图8 4个通孔缺陷主、旁瓣分布 Fig. 8 Main and side lobes of four cross-drilled hole defects

图8中可以看出,DCFMV算法的主瓣宽度最小。

为了比较DCFMV、MV、CFMV这3种算法在空域抽样数据下的成像时间,分别计算了3种成像算法的运行时间,15个样本取平均,结果如表3所示。

表3 空域抽样数据下算法成像时间对比 Tab. 3 Comparison of imaging time of algorithms under spatial sampling data

表3可以看出:CFMV算法耗时最长,MV算法次之,DCFMV算法耗时最少;对于裂纹缺陷,DCFMV算法相较于MV算法和CFMV算法成像时间分别减小了65.90%和66.10%;对于横通孔缺陷,DCFMV算法相较于MV和CFMV算法成像时间分别缩短了67.49%和68.61%。

4 结 论

为了降低波束形成算法缺陷成像时的计算复杂度,提高成像运行时效性,本文提出了一种超声阵列波束形成方法。该方法根据阵列阵元数确定最大抽取因子,对阵元数据空域抽样后,数据量得到大幅减少。将协方差矩阵 ${{R}}$ 构造为Toeplitz矩阵对其进行修正,算法计算复杂度降低;融合相干因子提高抽样数据中有效信息占比,进而提高成像质量。DCFMV算法的计算复杂度由 $O({N^3})$ 降为 $O({N^2}{\rm{/}}{D^2})$ 。通过在试块中裂纹和通孔缺陷的仿真成像,DCFMV算法的成像质量较好,与MV算法和CFMV算法相比,DCFMV算法时间效率提高分别提高85%和86%以上。数据空域抽样后,DCFMV算法下的成像API最小,成像质量最好。该算法减小了成像算法对硬件性能的过度依赖,可为超声波束形成及缺陷成像提供理论参考。下一步将结合等间隔抽样的信息损失评估,相干因子融合的效果评估,在矛盾制约条件下挖掘出最佳抽样相干融合策略。

参考文献
[1]
Sakhaei S M. Optimum beamforming for sidelobe reduction in ultrasound imaging[J]. IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control, 2012, 59(4): 799-805. DOI:10.1109/TUFFC.2012.2257
[2]
Shang Qiufeng,Zheng Guoqiang,Li Yan. NUFFT two-dimensional imaging algorithm combined with adaptive beamforming[J]. Optoelectronic Technology, 2020, 40(2): 94-99. [尚秋峰,郑国强,李炎. 结合自适应波束形成的NUFFT二维成像算法[J]. 光电子技术, 2020, 40(2): 94-99. DOI:10.19453/j.cnki.1005-488x.2020.02.004]
[3]
Sadeghi M,Mahloojifar A. A novel adaptive apodization to improve the resolution of phased subarray imaging in medical ultrasound[J]. Journal of Medical Ultrasonics, 2020, 47(1): 13-24. DOI:10.1007/s10396-019-00970-2
[4]
Synnevag J F,Austeng A,Holm S. Benefits of minimum-variance beamforming in medical ultrasound imaging[J]. IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control, 2009, 56(9): 1868-1879. DOI:10.1109/TUFFC.2009.1263
[5]
Wang Ping,Jiang Jinyang,Li Fang,et al. Signal to noise ratio dependent postfilter combined with eigenspace-based minimum variance algorithm for ultrasound imaging[J]. Acta Acustica, 2019, 44(1): 136-144. [王平,江金洋,李昉,等. 信噪比后滤波与特征空间融合的最小方差超声成像算法[J]. 声学学报, 2019, 44(1): 136-144. DOI:10.15949/j.cnki.0371-0025.2019.01.015]
[6]
Wang Ping,Xu Qin,Fan Wenzheng,et al. Eigenspace-based forward-backward minimum variance beamforming applied to ultrasound imaging[J]. Acta Acustica, 2013, 38(1): 65-70. [王平,许琴,范文政,等. 超声成像中基于特征空间的前后向最小方差波束形成[J]. 声学学报, 2013, 38(1): 65-70. DOI:10.15949/j.cnki.0371-0025.2013.01.016]
[7]
Li Paichi,Li Menglin. Adaptive imaging using the generalized coherence factor[J]. IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control, 2003, 50(2): 128-141. DOI:10.1109/TUFFC.2003.1182117
[8]
Liu Haolin,Yi Zongrui,Liu Dongquan. Amplitude and phase estimation combined coherence factor applied to medical ultrasound imaging[J]. Journal of Sichuan University(Engineering Science Edition), 2015, 47(Supp2): 125-129. [刘昊霖,易宗锐,刘东权. 相干系数与幅度相位估计融合的医学超声成像算法[J]. 四川大学学报(工程科学版), 2015, 47(增刊2): 125-129. DOI:10.15961/j.jsuese.2015.s2.019]
[9]
Lan Zhengfeng,Jin Liu,Feng Shuai,et al. Joint generalized coherence factor and minimum variance beamformer for synthetic aperture ultrasound imaging[J]. IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control, 2021, 68(4): 1167-1183. DOI:10.1109/TUFFC.2020.3035412
[10]
Wang Shunli,Li Paichi.High frame rate adaptive imaging using coherence factor weighting and the MVDR method[C]//Proceedings of the 2008 IEEE Ultrasonics Symposium.Beijing:IEEE,2008:1175–1178.
[11]
Asl B M,Mahloojifar A. Minimum variance beamforming combined with adaptive coherence weighting applied to medical ultrasound imaging[J]. IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control, 2009, 56(9): 1923-1931. DOI:10.1109/TUFFC.2009.1268
[12]
Xu Mengling,Chen Yimin,Ding Mingyue,et al.Adaptive minimum variance beamforming combined with phase coherence imaging for ultrasound imaging[C]//Proceedings of Medical Imaging 2012:Ultrasonic Imaging,Tomography,and Therapy.San Diego:SPIE,2012,8320:83200E.
[13]
Park J H,Yoon C,Chang J H,et al.A real-time synthetic aperture beamformer for medical ultrasound imaging[C]//Proceedings of the 2010 IEEE International Ultrasonics Symposium.San Diego:IEEE,2010:1992–1995.
[14]
Chen Kuanting,Ma W H,Hwang Y T,et al. A low complexity,high throughput DoA estimation chip design for adaptive beamforming[J]. Electronics, 2020, 9(4): 641. DOI:10.3390/electronics9040641
[15]
Kim K,Park S,Kim J,et al. A fast minimum variance beamforming method using principal component analysis[J]. IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control, 2014, 61(6): 930-945. DOI:10.1109/TUFFC.2014.2989
[16]
Park J,Wi S M,Lee J S. Computationally efficient adaptive beamformer for ultrasound imaging based on QR decomposition[J]. IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control, 2016, 63(2): 256-265. DOI:10.1109/TUFFC.2016.2515260
[17]
Nilsen C I C,Hafizovic I. Beamspace adaptive beamforming for ultrasound imaging[J]. IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control, 2009, 56(10): 2187-2197. DOI:10.1109/TUFFC.2009.1301
[18]
Asl B M,Mahloojifar A. A low-complexity adaptive beamformer for ultrasound imaging using structured covariance matrix[J]. IEEE Transactions on Ultrasonics,Ferroelectrics,and Frequency Control, 2012, 59(4): 660-667. DOI:10.1109/TUFFC.2012.2244
[19]
Sakhaei S M. A decimated minimum variance beamformer applied to ultrasound imaging[J]. Ultrasonics, 2015, 59: 119-127. DOI:10.1016/j.ultras.2015.02.005
[20]
Shamsian S E,Sakhaei S M. Fast adaptive beamforming through a cascade structure for ultrasound imaging[J]. Journal of Medical Ultrasonics, 2019, 46(3): 287-296. DOI:10.1007/s10396-019-00930-w
[21]
Li Jian,Stoica P,Wang Zhisong. On robust Capon beamforming and diagonal loading[J]. IEEE Transactions on Signal Processing, 2003, 51(7): 1702-1715. DOI:10.1109/TSP.2003.812831
[22]
Zohar S. The solution of a toeplitz set of linear equations[J]. Journal of the ACM, 1974, 21(2): 272-276. DOI:10.1145/321812.321822
[23]
Min Rui.Research on airborne SAR three-dimensional imaging theory and key technology[D].Chengdu:School of Electronic Engineering,2012.
闵锐.机载SAR三维成像理论及关键技术研究[D].成都:电子科技大学,2012.
[24]
Lu Quan,Xu Zhong,Ye Zhenglin. A new expression and a fast algorithm for the inversion of toeplitz matrix[J]. Journal of Numerical Methods and Computer Applications, 2005, 26(3): 191-197. [陆全,徐仲,叶正麟. Toeplitz矩阵之逆矩阵的新分解式及快速算法[J]. 数值计算与计算机应用, 2005, 26(3): 191-197. DOI:10.3969/j.issn.1000-3266.2005.03.004]
[25]
Holmes C,Drinkwater B W,Wilcox P D. Post-processing of the full matrix of ultrasonic transmit-receive array data for non-destructive evaluation[J]. NDT & E International, 2005, 38(8): 701-711. DOI:10.1016/j.ndteint.2005.04.002