CN113673419B

CN113673419B - 适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法

Info

Publication number: CN113673419B
Application number: CN202110953837.8A
Authority: CN
Inventors: 杨益新; 张亚豪
Original assignee: Northwestern Polytechnical University
Current assignee: Northwestern Polytechnical University
Priority date: 2021-08-19
Filing date: 2021-08-19
Publication date: 2024-05-28
Anticipated expiration: 2041-08-19
Also published as: CN113673419A

Abstract

本发明涉及一种适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，采用MVDR‑DL对强干扰进行抑制，构建关于MVDR‑DL波束功率输出线性关系的贝叶斯概率模型，将稀疏贝叶斯算法推广至波束域。在贝叶斯框架下进行DOA估计，避免了超参数选取的问题。本发明在进行DOA估计时，每次迭代仅更新一个网格点对应的信号参数，从而避免了矩阵求逆运算，有效降低了方法的计算量。在强干扰环境下实现了对目标信号方位有效估计的同时，增强了算法的实用性。

Description

适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法

技术领域

本发明属于信号处理等领域，特别涉及一种适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法。

背景技术

基于阵列接收信号的水下目标方位(Direction of arrival，DOA)估计是被动声纳信号处理中一项主要任务。近十几年稀疏重构类DOA估计算法因其对快拍数和信噪比的要求较低受到了广泛的关注。这类算法将空间划分为离散的网格，假设信号在有限个网格点上。通过对每个网格点上的信号参数进行估计，实现DOA估计。根据估计原理不同，算法可分为基于l_p范数的算法和稀疏贝叶斯算法。相比于基于l_p范数的算法来说，稀疏贝叶斯算法无需选取任何超参数，因此在实际信号处理中更易实现。大多数稀疏贝叶斯算法采用期望最大算法或变分贝叶斯算法迭代估计DOA值，需要大量的矩阵求逆运算，算法运算量较大。国外学者M.E.Tipping等人提出一种快速稀疏贝叶斯方法(M.E.Tipping and A.C.Faul,“Fast marginal likelihood maximisation for sparse Bayesian models,”inProc.9th Int.Workshop Artif.Intell.Stat.,2003,vol.1.)，该方法在每次迭代时仅更新一个网格对应的信号参数，从而避免了矩阵求逆，大大提升了计算效率。不同于主动声纳自主发射信号并通过接收反射回波进行目标探测，被动声纳是通过接收舰船辐射噪声来进行目标探测的，因而具有更好的隐蔽性。随着海洋战略和经济的不断发展，水面及水下舰船越来越多，这使得被动声纳的接收信号越来越复杂。当功率较低的目标信号远离接收阵时，距离接收阵较近的大功率水面舰船的辐射噪声等可看作为强干扰信号。强干扰信号的存在将影响对目标信号的DOA估计精度，甚至掩蔽目标信号。由于强干扰的存在，导致稀疏重构类算法对目标信号的方位估计性能存在一定程度的下降。

Yang等人(Y.Yang,Y.Zhang,and L.Yang,“Wideband sparse spatial spectrumestimation using matrix filter with nulling in a strong interferenceenvironment,”J.Acoust.Soc.Am.143(6),3891–3898(2018).)提出一种基于零陷矩阵滤波器的稀疏谱估计方法，使用零陷矩阵滤波器在干扰方向上形成深凹槽抑制强干扰，并利用稀疏谱估计方法进行弱目标信号方位估计。然而，设计零陷矩阵滤波器需求解一个凸优化问题，算法计算量较大。同时，稀疏谱估计算法属于l_p范数类算法，需选取一个合适的超参数保证算法性能，然而该参数选取通常较为困难，造成该算法在实际应用时存在一定局限性。

发明内容

本发明解决的技术问题是：为了高效解决强干扰环境下对目标信号的DOA估计问题，本发明给出一种基于波束功率输出的快速稀疏贝叶斯(Fast sparse Bayesianlearning based on beamformer power outputs，FSBL-BPO)方法。该方法采用基于对角加载的最小方差无失真响应(Minimum variance distortionless response with diagonalloading，MVDR-DL)波束形成器作为预处理器，在强干扰方向形成凹槽充分抑制强干扰信号，并计算波束功率输出。构建适用于MVDR-DL波束功率输出和波束响应线性关系的贝叶斯概率模型。在贝叶斯框架下对模型参数进行迭代更新，每次迭代时仅更新一个网格对应的信号参数，避免大量矩阵求逆计算，实现强干扰环境下对目标信号的快速方位估计。

本发明的技术方案是：适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，包括以下步骤：

步骤1：假设K_S个目标信号和K_D个干扰信号分别从和方向上入射至M元的均匀线阵列，目标信号和干扰信号之间互不相关；当阵列接收到信号后，将接收信号划分为N段，对每一段进行傅里叶变换后，宽带信号划分为L个子带；第l个子带上第n段阵列接收信号对应的傅里叶变换系数记为采样协方差矩阵计算为

上标“H”为共轭转置运算；

步骤2：建立波束域模型，包括以下步骤：

步骤2.1：通过常规波束形成算法确定目标信号所在方位区域Θ_S＝[Θ_SL,Θ_SR]，其中Θ_SL和Θ_SR分别为区域Θ_S的左右边界；将该区域均匀划分为K_B个网格对于第l个子带，该区域上K_B个MVDR-DL波束形成器组成的波束形成矩阵W_l为

其中为指向φ_k的MVDR-DL波束形成器的加权量，a_l(φ_k)为第l个子带上指向φ_k的阵列流形，为通过求解的(M-K_S-K_D)个小特征值的平均值得到的噪声功率估计值，I_M为M维的单位矩阵，上标“-1”表示矩阵求逆；

通过该矩阵对步骤1得到的采样协方差矩阵进行滤波，得到第l个子带上波束域协方差矩阵为

其中P_l ^S和P_l ^D表示第l个子带上的目标信号和干扰信号的协方差矩阵，σ_l表示第l个子带上的噪声功率，和表示第l个子带上目标信号和干扰信号的阵列流形矩阵，E_l为第l个子带上的扰动误差；

步骤2.2：对步骤2.1得到的协方差矩阵进行按列向量化运算，得到

其中和分别表示第l个子带上目标信号和干扰信号的功率向量，和分别表示和矩阵按列向量化的向量，⊙表示Khatri-Rao乘积，上标“*”为共轭运算；

步骤2.3：定义矩阵对于第m行第n列的元素[J]_mn，若[J]_mn＝1，否则[J]_mn＝0；将步骤2.2得到结果与该矩阵相乘，因MVDR-DL在区域Θ_S外对干扰的波束响应较低，因此第l个子带上波束功率输出的线性关系式表示为：

其中

步骤2.4：将区域Θ_S均匀划分为K_G个网格，网格点组成的向量记为基于该离散网格，步骤2.3中的公式重新表示为

式中为第l个子带上该网格上的阵列流形矩阵。p_l是一个稀疏向量，当p_l第m个元素等于第n个元素，否则为0；

步骤3：在贝叶斯框架下进行迭代计算，最终输出DOA估计值，包括以下子步骤：

步骤3.1：构建贝叶斯概率模型，其中扰动概率模型为：

其中N(·)表示实高斯分布，表示Hadamard积；

信号概率模型表示为：

p(p_l；γ)＝N(0,Γ^-1),l＝1,...,L

其中为信号稀疏参数，上标“T”为转置运算，Γ＝diag(γ)为以γ中元素为对角元素的对角矩阵；

构建好后，给定信号稀疏参数γ，信号功率后验协方差矩阵Σ_l，信号功率后验均值向量μ_l，噪声功率σ_l，以及参数S_l,p，迭代初始值；

步骤3.2：参数迭代求解，分别完成对信号稀疏参数γ，信号功率后验协方差矩阵Σ_l，信号功率后验均值向量μ_l，噪声功率σ_l，以及参数S_l,p，的更新；

在第i次迭代中，首先对γ进行更新：

其中上标(i)表示第i次迭代，表示γ⁽ⁱ⁾第p个元素，

p⁽ⁱ⁾表示使函数

最大时的下标，的值为：当时，否则为∞，∞表示无穷大；|·|表示取绝对值；记存放前(i-1)次迭代中更新下标的活跃集为Φ^(i-1)，对更新分为三种情况：

(1)若时，则Φ⁽ⁱ⁾＝Φ^(i-1)∪{p⁽ⁱ⁾}

(2)若时，则Φ⁽ⁱ⁾＝Φ^(i-1)

(3)若Φ⁽ⁱ⁾为将p⁽ⁱ⁾从Φ^(i-1)中移除的集合；

对应三种不同情况，首先完成对信号功率后验协方差矩阵Σ_l,l＝1,...,L的更新，再进行参数S_l,p，的更新；最终先后完成对信号功率后验均值向量μ_l,l＝1,...,L以及噪声功率σ_l,l＝1,...,L的更新；

更新后的参数值若满足迭代终止条件，则得到所估计的功率谱为

该功率谱中峰值对应的方位即为目标信号的DOA估计值；

若不满足迭代终止条件，则继续进行更新，直至满足迭代终止条件。

本发明进一步的技术方案是；所述步骤3.1中的各参数迭代初始值为：

信号稀疏参数的初始值γ⁽⁰⁾：

其中表示γ⁽⁰⁾第p个元素；

表示的第p列，p⁽⁰⁾表示使函数

最大时对应的下标，的值为：当Δ＞0时，否则

活跃集的初始值Φ⁽⁰⁾：

Φ⁽⁰⁾＝{p⁽⁰⁾}

第l个子带上信号功率后验协方差矩阵的初始值

第l个子带上信号功率后验均值向量的初始值

第l个子带上噪声功率的初始值

参数S_l,p，的初始值：

本发明进一步的技术方案是：所述步骤3.2中对应第一种情况时，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为

其中仅包含活跃集Φ^(i-1)中记录的“激活”网格对应的阵列流形，

和更新为

其中

本发明进一步的技术方案是：所述步骤3.2中对应第二种情况时，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为

其中为的第p⁽ⁱ⁾列，表示在第p⁽ⁱ⁾行第p⁽ⁱ⁾列的元素，为除了第p⁽ⁱ⁾个元素为1其余元素为0的向量；

和更新为

本发明进一步的技术方案是：所述步骤3.2中对应第三种情况时，将该公式代入第二种情况中进行计算，即可完成更新。

本发明进一步的技术方案是：所述步骤3.2中完成信号功率后验协方差矩阵Σ_l以及参数S_l,p，的更新后，对第l个子带上信号功率后验均值向量μ_l的更新为

最后，第l个子带上噪声功率σ_l的更新为

本发明进一步的技术方案是：所述步骤3.2中的迭代终止条件为：当迭代满足||γ⁽ⁱ⁾-γ^(i-1)||₂/||γ^(i-1)||₂≤10^-3，或者迭代次数大于Iter_max＝1000时，其中||·||₂表示l₂范数，迭代终止。

发明效果

本发明的技术效果在于：采用MVDR-DL波束形成器代替零陷矩阵滤波器来充分抑制强干扰信号，避免了强干扰信号对后续目标信号方位估计影响的同时，降低了计算量；计算MVDR-DL波束功率输出，构建了关于波束功率输出与波束响应的线性关系的贝叶斯概率模型，将稀疏贝叶斯类算法推广至波束域中，提升了稀疏贝叶斯算法在强干扰环境下对目标信号方位估计的稳健性；在贝叶斯框架下进行参数迭代更新，避免了超参数选取的问题；每次迭代时仅更新一个网格点对应的信号参数，避免大量矩阵求逆计算，从而在有效解决强干扰环境下对目标信号方位估计问题的同时，实现了快速DOA估计。

附图说明

图1基于波束功率输出的快速稀疏贝叶斯方法的总体流程框图

图2 FSBL-BPO迭代流程

图3采用常规波束形成算法估计目标所在区域结果

图4快速贝叶斯方法的方位估计结果

图5基于零陷矩阵滤波器的稀疏谱估计方法方位估计结果

图6 FSBL-BPO方法的方位估计结果

具体实施方式

理解的是，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”、“顺时针”、“逆时针”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

参见图1-图6，本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤1：假设K_S个目标信号和K_D个干扰信号分别从和方向上入射至M元的均匀线阵列，目标信号和干扰信号之间互不相关。当阵列接收到信号后，将接收信号划分为N段，对每一段进行傅里叶变换后，宽带信号划分为L个子带。第l个子带上第n段阵列接收信号对应的傅里叶变换系数记为采样协方差矩阵计算为

上标“H”为共轭转置运算。

步骤2：建立波束域模型，包括以下步骤：

步骤2.1：通过常规波束形成算法确定目标信号所在方位区域Θ_S＝[Θ_SL,Θ_SR]，其中Θ_SL和Θ_SR分别为区域Θ_S的左右边界。将该区域均匀划分为K_B个网格对于第l个子带，该区域上K_B个MVDR-DL波束形成器组成的波束形成矩阵W_l为

其中为指向φ_k的MVDR-DL波束形成器的加权量，a_l(φ_k)为第l个子带上指向φ_k的阵列流形，为通过求解的(M-K_S-K_D)个小特征值的平均值得到的噪声功率估计值，I_M为M维的单位矩阵，上标“-1”表示矩阵求逆。

通过该矩阵对步骤1得到的采样协方差矩阵进行滤波，得到第l个子带上的波束域协方差矩阵为

其中P_l ^S和P_l ^D表示第l个子带上的目标信号和干扰信号的协方差矩阵，σ_l表示第l个子带上的噪声功率，和表示第l个子带上目标信号和干扰信号的阵列流形矩阵，E_l为第l个子带上的扰动误差。

其中和分别表示第l个子带上目标信号和干扰信号的功率向量，和分别表示W_l ^HW和W_l ^HE_lW矩阵按列向量化的向量，⊙表示Khatri-Rao乘积，上标“*”为共轭运算。

其中

式中为第l个子带上该网格上的阵列流形矩阵。p_l是一个稀疏向量，当p_l第m个元素等于第n个元素，否则为0。

步骤3：在贝叶斯框架下进行迭代计算，最终输出DOA估计值，包括以下子步骤：步骤3.1：构建贝叶斯概率模型，其中扰动概率模型为：

其中N(·)表示实高斯分布，表示Hadamard积。信号概率模型表示为：

p(p_l；γ)＝N(0,Γ^-1),l＝1,...,L

其中为信号稀疏参数，上标“T”为转置运算，Γ＝diag(γ)为以γ中元素为对角元素的对角矩阵。

构建好后，给定信号稀疏参数γ，信号功率后验协方差矩阵Σ_l，信号功率后验均值向量μ_l，噪声功率σ_l，以及参数S_l,p，迭代初始值，如下信号稀疏参数的初始值γ⁽⁰⁾：

其中表示γ⁽⁰⁾第p个元素，

表示的第p列，p⁽⁰⁾表示使函数

最大时对应的下标，的值为：当Δ＞0时，否则

活跃集的初始值Φ⁽⁰⁾：

Φ⁽⁰⁾＝{p⁽⁰⁾}

第l个子带上信号功率后验协方差矩阵的初始值

第l个子带上信号功率后验均值向量的初始值

第l个子带上噪声功率的初始值

参数S_l,p，的初始值：

步骤3.2：参数迭代求解，分别完成对信号稀疏参数γ，信号功率后验协方差矩阵Σ_l，信号功率后验均值向量μ_l，噪声功率σ_l，以及参数S_l,p，的更新。

在第i次迭代中，首先对γ进行更新：

其中上标(i)表示第i次迭代，表示γ⁽ⁱ⁾第p个元素，

p⁽ⁱ⁾表示使函数

最大时的下标，的值为：当时，否则为∞，∞表示无穷大。|·|表示取绝对值。记存放前(i-1)次迭代中更新下标的活跃集为Φ^(i-1)，对更新分为三种情况：

(1)若时，则Φ⁽ⁱ⁾＝Φ^(i-1)∪{p⁽ⁱ⁾}

(2)若时，则Φ⁽ⁱ⁾＝Φ^(i-1)

(3)若Φ⁽ⁱ⁾为将p⁽ⁱ⁾从Φ^(i-1)中移除的集合。

针对第一种情况，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为

和更新为

其中

针对第二种情况，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为

其中为的第p⁽ⁱ⁾列，表示在第p⁽ⁱ⁾行第p⁽ⁱ⁾列的元素，为除了第p⁽ⁱ⁾个元素为1其余元素为0的向量。

和更新为

针对第三种情况，将该公式代入第二种情况中进行计算，即可完成更新。

完成信号功率后验协方差矩阵Σ_l以及参数S_l,p，的更新后，对第l个子带上信号功率后验均值向量μ_l的更新为

最后，第l个子带上噪声功率σ_l的更新为

更新后的参数值若满足||γ⁽ⁱ⁾-γ^(i-1)||₂/||γ^(i-1)||₂≤10^-3，其中||·||₂表示l₂范数，或者迭代次数大于Iter_max＝1000时，迭代终止，所估计的功率谱为

该功率谱中峰值对应的方位即为目标信号的DOA估计值；若不满足迭代终止条件，则继续进行更新，直至满足迭代终止条件。

为了验证利用本发明给出的方法在强干扰环境下对目标信号方位估计的有效性，设计仿真实验如下：假设接收阵列为M＝20元均匀线列阵，阵元间距为4m。两个远处航船从-10°和-7°方向上入射至阵列，功率为0dB，视为目标信号，即K_S＝2；一个近处水面舰从10°方向上入射至阵列，功率为20dB，大于远处航船，视为干扰信号，即K_D＝1。所考虑的信号频段为[90，180]Hz。噪声为高斯白噪声，其功率在所考虑的信号频段上为0dB。

将接收信号均匀划分为N＝50段，在每一段上进行傅里叶变换，将信号划分至L＝46个子带。图3为采用常规波束形成算法对目标信号所在区域估计的结果，虚线为目标信号所在区域范围Θ_S的边界Θ_SL＝-14°和Θ_SR＝-2°。故Θ_S为[-14,-2]°。

将该区域以2°为间隔均匀划分为7个网格点，得到MVDR-DL波束指向角，即K_B＝7。为了进行DOA估计，将该区域以1°为间隔均匀划分为13个网格点，即K_G＝13，在该网格上进行DOA估计。图4、图5和图6分别为快速贝叶斯方法、基于零陷矩阵滤波器稀疏谱估计方法和FSBL-BPO方法的方位估计结果，图中的虚线为目标信号的真实方位。可以看出，快速贝叶斯方法(图4)由于受到了强干扰信号的影响，没有办法估计出两个信号。本发明给出的方法(图6)和基于零陷矩阵滤波器的稀疏谱估计方法(图5)避免了强干扰信号的影响，可以很好地估计出两个目标信号，说明本方法成功将稀疏贝叶斯算法推广至波束域中，增强了稀疏贝叶斯算法在强干扰环境下对目标信号方位估计的稳健性。除此之外，本发明给出的方法进行一次DOA估计的运行时间为0.2s，而基于零陷矩阵滤波器的稀疏谱估计方法的运行时间为700s。可以看出，本发明给出的方法大大提升了计算效率，实现了快速DOA估计。

Claims

1.适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，包括以下步骤：

上标“H”为共轭转置运算；

步骤2：建立波束域模型，包括以下步骤：

其中和分别表示第l个子带上目标信号和干扰信号的功率向量，和分别表示W_l ^HW和W_l ^HE_lW矩阵按列向量化的向量，⊙表示Khatri-Rao乘积，上标“*”为共轭运算；

步骤2.3：定义矩阵对于第m行第n列的元素[J]_mn，若n＝K_B(m-1)+m，[J]_mn＝1，否则[J]_mn＝0；将步骤2.2得到结果与该矩阵相乘，因MVDR-DL在区域Θ_S外对干扰的波束响应较低，因此第l个子带上波束功率输出的线性关系式表示为：

其中

式中为第l个子带上该网格上的阵列流形矩阵，p_l是一个稀疏向量，当p_l第m个元素等于第n个元素，否则为0；

步骤3.1：构建贝叶斯概率模型，其中扰动概率模型为：

其中N(·)表示实高斯分布，表示Hadamard积；

信号概率模型表示为：

p(p_l；γ)＝N(0,Γ^-1),l＝1,...,L

在第i次迭代中，首先对γ进行更新：

其中上标(i)表示第i次迭代，表示γ⁽ⁱ⁾第p个元素，

p⁽ⁱ⁾表示使函数

最大时的下标，的值为：当时，否则为∞，∞表示无穷大；|·|表示取绝对值；

记存放前(i-1)次迭代中更新下标的活跃集为Φ^(i-1)，对更新分为三种情况：

(1)若时，则Φ⁽ⁱ⁾＝Φ^(i-1)∪{p⁽ⁱ⁾}

(2)若时，则Φ⁽ⁱ⁾＝Φ^(i-1)

(3)若Φ⁽ⁱ⁾为将p⁽ⁱ⁾从Φ^(i-1)中移除的集合；

该功率谱中峰值对应的方位即为目标信号的DOA估计值；

2.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.1中的各参数迭代初始值为：

信号稀疏参数的初始值γ⁽⁰⁾：

其中表示γ⁽⁰⁾第p个元素；

表示的第p列，p⁽⁰⁾表示使函数

最大时对应的下标，的值为：当Δ＞0时，否则活跃集的初始值Φ⁽⁰⁾：

Φ⁽⁰⁾＝{p⁽⁰⁾}

第l个子带上信号功率后验协方差矩阵的初始值

第l个子带上信号功率后验均值向量的初始值

第l个子带上噪声功率的初始值

参数S_l,p，的初始值：

3.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.2中对应第一种情况时，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为

和更新为

其中

4.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.2中对应第二种情况时，第i次迭代中第l个子带上信号功率后验协方差矩阵更新为

和更新为

5.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.2中对应第三种情况时，将该公式代入第二种情况中进行计算，即可完成更新。

6.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.2中完成信号功率后验协方差矩阵Σ_l以及参数S_l,p，的更新后，对第l个子带上信号功率后验均值向量μ_l的更新为

最后，第l个子带上噪声功率σ_l的更新为

7.如权利要求1所述的适用于强干扰环境的波束域快速稀疏贝叶斯方位估计方法，其特征在于，所述步骤3.2中的迭代终止条件为：当迭代满足||γ⁽ⁱ⁾-γ^(i-1)||₂/||γ^(i-1)||₂≤10^-3，或者迭代次数大于Iter_max＝1000时，其中||·||₂表示l₂范数，迭代终止。