CN112120730A - 一种基于混合子空间投影的广义旁瓣相消超声成像方法 - Google Patents

Abstract

本发明涉及一种基于混合子空间投影的广义旁瓣相消超声成像方法，属于超声成像技术领域。该方法包括：预处理接收的回波信号；划分子阵，对其前后向平滑和对角加载处理，获得样本协方差矩阵；对样本协方差矩阵进行特征分解，并依据样本协方差矩阵的特征值大小以及期望方向向量与样本协方差矩阵的特征矢量相关性，构造混合信号子空间；将广义旁瓣相消器的期望方向向量投影到混合信号子空间的左奇异空间，并基于改进后的期望方向向量构造新的阻塞矩阵；形成广义旁瓣相消波束；将该波束的权值投影到混合信号子空间进行修正；利用修正后的权值对回波信号加权求和，得到合成信号，本发明可以显著提高波束形成算法的稳健性和成像质量。

Description

一种基于混合子空间投影的广义旁瓣相消超声成像方法

技术领域

本发明属于超声成像技术领域，涉及一种基于混合子空间投影的广义旁瓣相消超声成像算法。

背景技术

超声成像具有安全、实时、便捷等优点，因此在医学领域得到广泛应用。数字波束合成是超声成像最为关键的环节。而数字波束合成算法中应用最广泛、最简单的波束形成算法是延时叠加(Delay And Sum,DAS)算法，它是根据阵元通道几何位置关系对所接收的回波信号进行延时量的计算，然后对延时后的数据对齐叠加。传统DAS算法复杂度低，成像速度快，但由于其采用固定窗函数加权导致主瓣宽度增加，分辨率较低。

最小方差无失真波束形成器(Minimum Variance,MV)能够有效提高成像图像分辨率，通过干扰信号和期望信号的方向性差异对干扰及噪声信号进行滤除。然而其稳健性差，不利于实际应用。广义旁瓣相消器(Generalized Sidelobe Canceller,GSC)作为最小方差无失真波束形成器的一种等效结构其对成像的改善效果与其相差不大。近年来，有学者提出将特征空间算法引入广义旁瓣相消器算法(Eigen Space Generalized SidelobeCanceller,ESGSC)中，通过将广义旁瓣相消自适应权矢量投影到基于协方差矩阵特征值大小所构造出的信号子空间中，以提高广义旁瓣消息器的对比度和噪声抑制能力。然而基于特征空间的广义旁瓣相消算法对于复杂环境成像中，往往会造成图像失真。尤其对于强散斑以及强散射点周围的背景进行成像时，该算法会产生较为严重的黑区伪像。这是由于干扰与噪声对应的特征值大于期望信号所对应的特征值，从而导致了传统特征空间广义旁瓣相消器稳健性变差。

综上所述，目前亟需一种既能保持特征空间广义旁瓣相消算法对比度及噪声抑制性能，又能保持强散斑背景不失真、具有强稳健性的波束形成算法。

发明内容

有鉴于此，本发明的目的在于提供一种基于混合子空间投影的广义旁瓣相消超声成像方法，在保持传统特征空间广义旁瓣相消自适应波束形成算法对比度性能基本不变的情况下，有效提高算法稳健性，改善强散斑背景成像效果，从而提高算法综合成像质量。

为达到上述目的，本发明提供如下技术方案：

一种基于混合子空间投影的广义旁瓣相消超声成像算法，具体包括以下步骤：

S1：对超声阵元接收的回波信号进行放大、AD转换和延时聚焦处理，以获得超声回波数据；得到延时聚焦处理之后的信号x_d(k)；

S2：将接收阵列依次划分为具有重叠阵元的子阵，然后对相应接收子阵的回波信号进行前后向平滑和对角加载处理，以获得样本协方差矩阵；

S3：对样本协方差矩阵进行特征分解，并依据样本协方差矩阵的特征值大小以及广义旁瓣相消器的期望方向向量与样本协方差矩阵的特征矢量的相关性，共同构造混合信号子空间；

S4：将广义旁瓣相消器的期望方向向量投影到混合信号子空间的左奇异空间，并基于改进后的期望方向向量构造新的阻塞矩阵；

S5：利用改进后的期望方向向量和新的阻塞矩阵进行广义旁瓣相消自适应波束形成；

S6：将广义旁瓣相消波束形成后得到的权值投影到混合信号子空间进行修正；

S7：利用修正后的权值对超声阵列回波信号进行加权求和得到最终的自适应波束合成信号。

进一步，在步骤S2中，进行前后向平滑和对角加载处理，得到估计样本协方差矩阵，具体包括以下步骤：

S21：把N个阵元依次划分为阵元数目为L的子阵，并分别计算各个子阵的样本协方差矩阵R_l(k)，然后根据以下公式计算前向估计协方差矩阵

公式中

表示第l个子阵的前向输出向量，

为

的共轭转置。其中，

表示第k个采样时刻下第l个阵元的聚焦后的超声回波数据，[·]^T表示转置运算；

S22：通过下式计算得到后向估计协方差矩阵

公式中

表示第l个子阵的后向输出向量，

表示

的共轭转置；

S23：通过以下计算公式计算前向估计协方差矩阵和后向估计协方差矩阵的求和平均，得到前后向估计协方差矩阵

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

进行对角加载，得到对角加载后的协方差矩阵

其中，

Δ为空间噪声与信号功率之比，

为信号的等效功率，I为单位矩阵。

进一步，在步骤S3中，对样本协方差矩阵进行特征分解，并依据特征值大小以及广义旁瓣相消器的期望方向向量与样本协方差矩阵的特征矢量的相关性，共同构造混合信号子空间，具体包括以下步骤：

S31：通过下式对

进行特征分解：

其中，λ_i为

的特征值，e_i为λ_i对应的特征向量，

为e_i的共轭转置，特征向量矩阵E＝[e₁ … e_N]；E^H为E的共轭转置，特征值矩阵Λ＝diag[λ₁ … λ_N]；

S32：基于样本协方差矩阵特征值大小构造信号子空间的第一部分EA_s：

首先对样本协方差矩阵的特征值大小进行排序：λ_max≥λ₂≥λ_i≥···≥λ_min，其中λ_i(i＝1,2,···,N)，ea_i(i＝1,2,···,N)为特征值λ_i对应的模值归一化后的特征向量。从最大的特征值λ_max开始，依次计算相邻两个特征值的比值λ_i/λ_i+1，当比值大于设定阈值α时，则之前的H个特征值对应的特征向量构成期望阈值信号子空间：EA_s＝[ea₁,ea₂,···,ea_H]；

S33：基于样本协方差矩阵的特征矢量与广义旁瓣相消器的期望方向向量的相关性构造信号子空间的第二部分Ep_s：

首先，分别计算分解后并进行模值归一化后的样本协方差矩阵的特征矢量与广义旁瓣相消器的期望方向向量的投影值pm(i)(i＝1,2,···,N)：

其中，

表示广义旁瓣相消器的期望方向向量，|·|表示求取绝对值，将所有投影模值按照从大到小排列为：

pm_[1] pm_[2] pm_[3] … pm_[N]

此时与pm_[z]对应的模值归一化后的特征矢量依次为：ep_z(z＝1,2,···,N)。

pm(i)≥β·max(pm(i))

然后给定合适的门限值β，利用所有满足上述不等式的Z个特征向量共同构建基于相关性的信号子空间Ep_s＝[ep₁,ep₂,···,ep_Z]。

S34:得到两次分解的信号子空间EA_s和Ep_s后，将二者拼接，得到最终的混合信号子空间E_mix-s：

E_mix-s＝[EA_s Ep_s]＝[ea₁,ea₂,···,ea_H,ep₁,ep₂,···,ep_Z]

进一步，在步骤S4中，将广义旁瓣相消器的期望方向向量投影到混合信号子空间的左奇异空间中，并基于改进后的期望方向向量构造新的阻塞矩阵：

首先，对混合子空间E_mix-s进行奇异值分解：E_mix-s＝UDV^H，其中，U是左奇异空间，D为奇异值在对角线上的对角矩阵，V是右奇异矩阵，V^H表示V的共轭转置。

其次，将广义旁瓣相消器的期望方向向量

投影到混合子空间E_mix-s较大特征值所对应的左奇异矢量空间中：

其中U^H表示U的共轭转置。

最后，求取满足方程组

的一组解，以此构造出新的阻塞矩阵B_mix-s，其中0表示“零矩阵”。

进一步，步骤S5中，利用改进后的期望方向向量和新的阻塞矩阵进行广义旁瓣相消自适应波束形成：

w＝w_q-B_mix-sw_a

其中，w表示广义旁瓣相消器的加权矢量，一般情况下w_q＝[1,1,....,1]^T表示长度为L的广义旁瓣相消器主干路的非自适应权矢量，w_a为广义旁瓣相消器辅助支路的自适应权矢量；其中，w_q与w_a的关系为：

其中，B_mix-s ^H为阻塞矩阵B_mix-s的共轭转置，

为样本协方差矩阵，(·)^-1表示矩阵求逆运算。

进一步，在步骤S6中，将广义旁瓣相消波束形成后得到的权值投影到混合信号子空间进行修正：

其中，w_mix-ES表示最终的基于混合子空间投影的广义旁瓣相消算法自适应加权矢量。E_mix-s表示混合信号子空间，

为E_mix-s的共轭转置，B_mix-s表示改进后的阻塞矩阵。

进一步，在步骤S7中，利用修正后的权值对超声阵列回波信号进行加权求和，得到最终的自适应波束合成信号：

其中，y_mix-EGSC(k)表示计算得到的第k个采样时刻下的自适应波束信号，w_mix-ES表示最终的基于混合子空间投影的广义旁瓣相消算法自适应加权矢量，

表示w_mix-ES的共轭转置。

本发明的有益效果在于：本发明相比于现有的特征空间广义旁瓣相消算算法，能够在保证其高对比度性能的同时，有效避免图像失真问题，提高算法稳健性，改善背景成像质量。本发明可以大幅度降低现于特征空间广义旁瓣相消算法的背景区域方差，同时避免黑区伪像的产生，在保证对比度和分辨率质量的同时，大幅度改善强散斑及强散射点周围的背景成像效果。

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

附图说明

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

图1为本发明所述广义旁瓣相消超超声成像方法的流程图；

图2为本发明基于混合信号子空间投影广义旁瓣相消器整体结构框架图；

图3为5种算法点目标成像结果；

图4为点目标成像50mm处5种算法横向分辨率曲线图；

图5为5种算法吸声斑目标成像结果；

图6为5种算法强散斑目标成像结果；

图7为5种算法geabr_0数据成像结果；

图8为geabr_0实验77.5mm处散射点横向截面图。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

请参阅图1～图8，图1为本发明提供的超声成像方法流程图，如图1和图2所示，一种基于混合子空间投影的广义旁瓣相消超声成像算法，包括以下步骤：

步骤S1：对超声阵元接收的回波信号进行放大、AD转换和延时聚焦处理，以获得超声回波数据；得到延时聚焦处理之后的信号x_d(k)，k表示为对应采样深度的采样时刻。

步骤S2：将接收阵列依次划分为具有重叠阵元的子阵，然后对相应接收子阵的回波信号进行前后向平滑和对角加载处理，以获得样本协方差矩阵，具体包括以下步骤：

S21：把N个阵元依次划分为阵元数目为L的子阵，并分别计算各个子阵的样本协方差矩阵R_l(k)，然后根据以下公式计算前向估计协方差矩阵

其中，

表示第l个子阵的前向输出向量，

为

的共轭转置。其中，

表示第k个采样时刻下第l个阵元的聚焦后的超声回波数据，[·]^T表示转置运算；

S22：通过下式计算得到后向估计协方差矩阵

其中，

表示第l个子阵的后向输出向量，

表示

的共轭转置；

S23：通过以下计算公式计算前向估计协方差矩阵和后向估计协方差矩阵的求和平均，得到前后向估计协方差矩阵

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

进行对角加载，得到对角加载后的协方差矩阵

其中，

Δ为空间噪声与信号功率之比，

为信号的等效功率，I为单位矩阵。

步骤S3：对样本协方差矩阵进行特征分解，并依据样本协方差矩阵的特征值大小以及广义旁瓣相消器的期望方向向量与样本协方差矩阵的特征矢量的相关性，共同构造混合信号子空间，具体包括以下步骤：

S31：通过下式对

进行特征分解：

其中，λ_i为

的特征值，e_i为λ_i对应的特征向量，

为e_i的共轭转置，特征向量矩阵E＝[e₁ … e_N]；E^H为E的共轭转置，特征值矩阵Λ＝diag[λ₁ … λ_N]；

S32：基于样本协方差矩阵特征值大小构造信号子空间的第一部分EA_s，具体包括：

首先对样本协方差矩阵的特征值大小进行排序：λ_max≥λ₂≥λ_i≥···≥λ_min，其中λ_i(i＝1,2,···,N)，ea_i(i＝1,2,···,N)为特征值λ_i对应的模值归一化后的特征向量。从最大的特征值λ_max开始，依次计算相邻两个特征值的比值λ_i/λ_i+1，当比值大于设定阈值α时，则之前的H个特征值对应的特征向量构成期望阈值信号子空间：EA_s＝[ea₁,ea₂,···,ea_H]；

S33：基于样本协方差矩阵的特征矢量与广义旁瓣相消器的期望方向向量的相关性构造信号子空间的第二部分Ep_s，具体包括：

首先，分别计算分解后并进行模值归一化后的样本协方差矩阵的特征矢量与广义旁瓣相消器的期望方向向量的投影值pm(i)(i＝1,2,···,N)：

其中，

表示广义旁瓣相消器的期望方向向量，|·|表示求取绝对值，将所有投影模值按照从大到小排列为：

pm_[1] pm_[2] pm_[3] … pm_[N]

此时与pm_[z]对应的模值归一化后的特征矢量依次为：ep_z(z＝1,2,···,N)。

pm(i)≥β·max(pm(i))

然后给定合适的门限值β，利用所有满足上述不等式的Z个特征向量共同构建基于相关性的信号子空间Ep_s＝[ep₁,ep₂,···,ep_Z]。

S34:得到两次分解的信号子空间EA_s和Ep_s后，将二者拼接，得到最终的混合信号子空间E_mix-s：

E_mix-s＝[EA_s Ep_s]＝[ea₁,ea₂,···,ea_H,ep₁,ep₂,···,ep_Z]

步骤S4：将广义旁瓣相消器的期望方向向量投影到混合信号子空间的左奇异空间，并基于改进后的期望方向向量构造新的阻塞矩阵：

首先，对混合子空间E_mix-s进行奇异值分解：E_mix-s＝UDV^H，其中，U是左奇异空间，D为奇异值在对角线上的对角矩阵，V是右奇异矩阵，V^H表示V的共轭转置。

其次，将广义旁瓣相消器的期望方向向量

投影到混合子空间E_mix-s较大特征值所对应的左奇异矢量空间中：

其中U^H表示U的共轭转置。

最后，求取满足方程组

的一组解，以此构造出新的阻塞矩阵B_mix-s，其中0表示“零矩阵”。

步骤S5：利用改进后的期望方向向量和新的阻塞矩阵进行广义旁瓣相消自适应波束形成：

w＝w_q-B_mix-sw_a

其中，w表示广义旁瓣相消器的加权矢量，一般情况下w_q＝[1,1,....,1]^T表示长度为L的广义旁瓣相消器主干路的非自适应权矢量，w_a为广义旁瓣相消器辅助支路的自适应权矢量；其中，w_q与w_a的关系为：

其中，B_mix-s ^H为阻塞矩阵B_mix-s的共轭转置，

为样本协方差矩阵，(·)^-1表示矩阵求逆运算。

步骤S6：将广义旁瓣相消波束形成后得到的权值投影到混合信号子空间进行修正：

其中，w_mix-ES表示最终的基于混合子空间投影的广义旁瓣相消算法自适应加权矢量。E_mix-s表示混合信号子空间，

为E_mix-s的共轭转置，B_mix-s表示改进后的阻塞矩阵。

步骤S7：利用修正后的权值对超声阵列回波信号进行加权求和，得到最终的自适应波束合成信号：

其中，y_mix-EGSC(k)表示计算得到的第k个采样时刻下的自适应波束信号，w_mix-ES表示最终的基于混合子空间投影的广义旁瓣相消算法自适应加权矢量，

表示w_mix-ES的共轭转置。

验证实验：

Field II是丹麦理工大学基于声学原理开发的一款超声实验仿真平台，其在理论研究上获得了广泛的认可和使用。为验证本发明算法的有效性，利用Field II对超声成像中常用的点散射目标和吸声斑目标以及强散斑目标进行成像，并利用实际实验数据geabr_0进行成像对比实验。在点目标仿真实验中，在纵向深度为40mm、45mm、65mm、70mm处分别置一散射点，横向距离位于中心0mm处的。并且在纵向55mm处，设置横向相邻间隔为±2mm的三个散射点。在50mm、和60mm处分别设置2个间隔为3mm的散射点，以此观察各个算法的横向分辨率强弱。采用发射定点聚焦和接收动态聚焦方式，发射聚焦为50mm处，并设置图像的成像动态范围为60dB。同时，设一中心在40mm，半径为3mm的圆形区域吸声斑，外部随机分布着100000个散射点，采用5种算法分别进行成像，设定成像动态范围为60dB。另外，设置强散斑成像实验：首先在30mm-50mm深度随机分布100000个幅值为区间[-1,1]中随机数的散射点。其次，将中心点位于直角系坐标(0mm,40mm),半径为3mm的斑状囊肿内的幅值提高至原先的40倍，设定成像动态范围为80dB。实验所采用的阵元中心频率为3.33MHz，阵元数目为64个，间距为0.2413mm，采样频率为17.76MHz，声速为1500m/s，设成像动态范围为60dB。

对上述四个实验目标采用延时叠加算法(DAS)、传统广义旁瓣相消算法(CF)、基于特征空间分解的广义旁瓣相消算法(ESGSC)、基于相关性特征空间分解的广义旁瓣相消算法(CEDGSC)以及本发明的基于混合信号子空间的广义旁瓣相消算法(MixEGSC)进行对比成像实验。

图3给出了上述5种算法点目标成像结果，从图3可以看出DAS算法其横向伪影最多、分辨率最低。GSC算法通过阻塞矩阵对消干扰噪声，有效抑制了旁瓣伪像的形成，提高对比度。ESGSC在GSC的基础上进一步提高分辨率和对比度，对于50mm处散射点成像几乎没有伪影黏连，这是由于ESGSC算法通过特征空间投影改进阻塞矩阵的阻塞能力以及期望信号的准确性，使得广义旁瓣相消器更充分的发挥出噪声对消的作用。CEDGSC算法基于期望方向向量与特征矢量相关性分解构造信号子空间，能提高ESGSC的稳健性，但是在成像分辨率上介于GSC与ESGSC之间。MixEGSC算法将ESGSC与CEDGSC算法的空间分解方法相结合。能在提高算法稳健性情况下，不削弱ESGSC算法原先的高对比度优势。

图4为上述5种算法点目标50mm处横向分辨率对比图，可以直观的看到ESGSC对对比度的提高最高，而MixEGSC的图像对比度仅次于ESGSC。为了更加直观的对比5种算法成像分辨率，表2给出了不同深度下5种算法-10dB的半峰值宽度(FWHM)数据对比，从表2可以看出，DAS分辨率最低，而GSC算法分辨率最高。三种基于不同特征空间分解的广义旁瓣相消算法在分辨率上较为接近。MixEGSC构造出的混合子空间投影法仍然可以保持ESGSC高分辨率、高对比度的优势。

表2.不同深度下5种算法-10dB的FWHM对比

图5为上述5种算法吸声斑目标成像结果。从图5可以看出，5种算法对于暗斑轮廓都能较为准确清晰地成像，然而DAS的斑内部存在大量旁瓣伪像，这是由于其分辨率与对比度不足导致的。相对于DAS算法有了很大的改善，GSC算法也存在较明显伪影，这是因为GSC算法虽然能够有效提高分辨率，但是在提升对比度方面仍然不足。ESGSC算法能够有效提高算法的对比度，因此他对于黑斑内的伪像抑制能力最佳。然而ESGSC存在明显缺点在于算法稳健性差，易引起图像失真，因此提出CEDGSC算法对其稳健性进行修正。但是CEDGSC算法会降低ESGSC算法的对比度，因此斑内伪影抑制效果变差。MixEGSC的暗斑成像效果和伪影抑制能力能够保持几乎与ESGSC相同，这是由于它结合了CEDGSC与ESGSC的优点，通过构造混合信号子空间的方式优化GSC阻塞矩阵与权值，因此能在提高ESGSC算法稳健性的同时，保持其高对比度、强抑噪能力的优势不受影响。

表3.不同算法暗斑背景质量指标比较

表3提供了不同算法暗斑背景成像质量指标，可以看出DAS虽然具有优秀的背景成像质量，然而其内部平均功率过高，即斑内伪像抑制能力弱，因此成像效果不佳。GSC算法的优势在于内部平均功率得到了改善，但是外部平均功率不佳，导致CR降低。ESGSC虽然具有最高的CR，但是其背景区域方差数值最大，表示其背景平滑程度差，背景成像质量低。因此提出CEDGSC算法对其稳健性进行修正。但是CEDGSC算法会降低ESGSC算法的CR值。而MixEGSC算法能够在对ESGSC背景均匀性进行修正的同时，最大程度的保持ESGSC在CR方面的优势。

为了验证CEDGSC算法与MixEGSC算法对ESGSC算法稳健性及背景成像质量的改进效果。图6给出了5种算法的强散斑成像效果对比，从图6中可以看出，ESGSC算法背景质量最差，因为存在明显的图像失真，DAS斑内部亮度效果最好，背景均匀性好。而GSC和CEDGSC算法虽然不存在图像失真，但是其背景均匀性仍有待提高。MixEGSC算法图像没有失真，并且具有良好的背景质量。

表4.不同算法强散斑背景质量指标比较

表4给出了不同算法强散斑背景质量指标的对比，以更直观的反映MixEGSC及CEDGSC对ESGSC算法的背景成像质量改善效果。从表中可以看出CEDGSC降低了ESGSC的背景区域方差，但是其CR值也遭到了大幅度的削弱，而MixEGSC算法大幅度改善ESGSC背景区域方差的同时保持优秀的CR值，可以看到其CR数值在5种算法当中仅次于ESGSC，并且在CNR指标上相对于ESGSC也得到了明显的改善。

图7给出了5种算法geabr_0数据成像结果。从图7中可以看出，DAS存在明显的暗斑旁瓣伪像，图像受到噪声干扰严重，GSC虽然提高了图像分辨率，但其对比度不佳，ESGSC算法在GSC的基础上大幅度提高了算法的对比度，然而其造成了验证的图像失真问题，主要反映在强散射点以及强散斑周围的图像背景严重抑制，形成黑斑。CEDGSC图像虽然解决了背景失真问题，但是仍存在明显旁瓣伪像，这表明它虽然ESGSC的基础上提高了算法稳健性，但是也极大的损失了对比度性能。最后，MixEGSC算法在图像背景上得到了有效的提高，它的亮斑成像相较于其他算法更加清晰，并且对斑内的旁瓣伪像进行有效的抑制，图像干净清晰，成像质量有显著的提高。它兼顾了CEDGSC对ESGSC在成像背景质量稳健性的改善作用，同时也很好的保持了ESGSC算法的性能。

图8给出了5种算法在77.5mm处散射点的横向分辨率对比曲线。从图中可以看出DAS的分辨率及对比度最低。GSC在分辨率上相对于DAS得到了有效的改善。而ESGSC算法在GSC的基础上进一步地提高了算法对比度。CEDGSC是为了ESGSC算法背景成像稳健性而设计的，然而它会造成ESGSC在对比度上的优势被削弱，而MixEGSC很好地兼顾了两者的优化，在提高ESGSC算法稳健性的同时，能够很大程度的保持该算法在对比度上的优势，实现良好的综合成像能力。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (9)

1.一种基于混合子空间投影的广义旁瓣相消超声成像方法，其特征在于，该方法具体包括以下步骤：

S1：对超声阵元接收的回波信号进行放大、AD转换和延时聚焦处理，得处理后的超声回波信号；

S2：将接收阵列依次划分为具有重叠阵元的子阵，然后对相应子阵的回波信号进行前后向平滑和对角加载处理，获得样本协方差矩阵；

S3：对样本协方差矩阵进行特征分解，并依据样本协方差矩阵的特征值大小以及广义旁瓣相消器的期望方向向量与样本协方差矩阵的特征矢量相关性，共同构造混合信号子空间；

S4：将广义旁瓣相消器的期望方向向量投影到混合信号子空间的左奇异空间，并基于改进后的期望方向向量构造新的阻塞矩阵；

S5：利用改进后的期望方向向量和新的阻塞矩阵进行广义旁瓣相消自适应波束形成；

S6：将广义旁瓣相消波束形成后得到的权值投影到混合信号子空间进行修正；

S7：利用修正后的权值对超声阵列回波信号进行加权求和，得到最终的自适应波束合成信号。

2.根据权利要求1所述的广义旁瓣相消超声成像方法，其特征在于，步骤S2中，得到样本协方差矩阵的具体步骤包括：

S21：把N个阵元依次划分为阵元数目为L的子阵，并分别计算各个子阵的样本协方差矩阵R_l(k)，然后根据以下公式计算前向估计协方差矩阵

其中，

表示第k个采样时刻下第l个子阵的聚焦后的超声回波数据，即为第l个子阵的前向输出向量；

的共轭转置，[·]^T表示转置运算；

S22：通过下式计算得到后向估计协方差矩阵

其中，

表示第l个子阵的后向输出向量，

表示

的共轭转置；

S23：通过以下计算公式计算前向估计协方差矩阵和后向估计协方差矩阵的求和平均，得到前后向估计协方差矩阵

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

进行对角加载，得到对角加载后的样本协方差矩阵

其中，

Δ为空间噪声与信号功率之比，

为信号的等效功率，I为单位矩阵。

3.根据权利要求1所述的广义旁瓣相消超声成像方法，其特征在于，步骤S3中，构造混合信号子空间的具体步骤包括：

S31：通过下式对

进行特征分解：

其中，λ_i为

的特征值，e_i为λ_i对应的特征向量，

为e_i的共轭转置，特征向量矩阵E＝[e₁…e_N]；E^H为E的共轭转置，特征值矩阵Λ＝diag[λ₁…λ_N]；

S32：基于样本协方差矩阵特征值大小构造期望阈值信号子空间EA_s；

S33：基于样本协方差矩阵的特征矢量与广义旁瓣相消器的期望方向向量的相关性构造基于相关性的信号子空间Ep_s；

S34：得到两次分解的信号子空间EA_s和Ep_s后，将二者拼接，得到最终的混合信号子空间：E_mix-s＝[EA_s Ep_s]＝[ea₁,ea₂,…,ea_H,ep₁,ep₂,…,ep_Z]。

4.根据权利要求3所述的广义旁瓣相消超声成像方法，其特征在于，步骤S32中，构造期望阈值信号子空间EA_s，具体包括：首先对样本协方差矩阵的特征值大小进行排序：λ_max≥λ₂≥λ_i≥…≥λ_min，其中ea_i为特征值λ_i对应的模值归一化后的特征向量；从最大的特征值λ_max开始，依次计算相邻两个特征值的比值λ_i/λ_i+1，当比值大于设定阈值α时，则之前的H个特征值对应的特征向量构成期望阈值信号子空间：EA_s＝[ea₁,ea₂,...,ea_H]。

5.根据权利要求3所述的广义旁瓣相消超声成像方法，其特征在于，步骤S33中，构造基于相关性的信号子空间Ep_s，具体包括：首先计算模值归一化后的样本协方差矩阵的特征矢量与广义旁瓣相消器的期望方向向量的投影值pm(i)：

其中，

表示广义旁瓣相消器的期望方向向量，|·|表示求取绝对值，将所有投影模值按照从大到小排列为：

pm_[1] pm_[2] pm_[3]…pm_[N]

此时与pm_[z]对应的模值归一化后的特征矢量依次为：ep_Z，Z＝1,2,…,N；

pm(i)≥β·max(pm(i))

然后选择门限值β，利用所有满足上述不等式的Z个特征向量共同构建基于相关性的信号子空间：Ep_s＝[ep₁,ep₂,…,ep_Z]。

6.根据权利要求1所述的广义旁瓣相消超声成像方法，其特征在于，步骤S4中，构造新的阻塞矩阵，具体包括：

首先，对混合子空间E_mix-s进行奇异值分解：E_mix-s＝UDV^H，其中U是左奇异空间，D为奇异值在对角线上的对角矩阵，V是右奇异矩阵，V^H表示V的共轭转置；

其次，将广义旁瓣相消器的期望方向向量

投影到混合子空间E_mix-s大特征值所对应的左奇异矢量空间中：

其中U^H表示U的共轭转置；

最后，求取满足方程组

的一组解，以此构造出新的阻塞矩阵B_mix-s，其中0表示“零矩阵”。

7.根据权利要求1所述的广义旁瓣相消超声成像方法，其特征在于，步骤S5中，形成广义旁瓣相消自适应波束，具体包括：

w＝w_q-B_mix-sw_a

其中，w表示广义旁瓣相消器的加权矢量，一般情况下w_q＝[1,1,....,1]^T表示长度为L的广义旁瓣相消器主干路的非自适应权矢量，w_a为广义旁瓣相消器辅助支路的自适应权矢量；其中，w_q与w_a的关系为：

其中，B_mix-s ^H为阻塞矩阵B_mix-s的共轭转置，

为样本协方差矩阵，(·)^-1表示矩阵求逆运算。

8.根据权利要求1所述的广义旁瓣相消超声成像方法，其特征在于，步骤S6中，将广义旁瓣相消波束形成后得到的权值投影到混合信号子空间进行修正：

其中，w_mix-ES表示最终的基于混合子空间投影的广义旁瓣相消算法自适应加权矢量，E_mix-s表示混合信号子空间，

为E_mix-s的共轭转置，B_mix-s表示改进后的阻塞矩阵。

9.根据权利要求1所述的广义旁瓣相消超声成像方法，其特征在于，步骤S7中，利用修正后的权值对超声阵列回波信号进行加权求和，得到最终的自适应波束合成信号：

其中，y_mix-EGSC(k)表示计算得到的第k个采样时刻下的自适应波束信号，

表示w_mix-ES的共轭转置。

Images