CN114280533B - 一种基于l0范数约束的稀疏贝叶斯DOA估计方法 - Google Patents

Abstract

本发明是一种基于l0范数约束的稀疏贝叶斯DOA估计方法。本发明涉及水下声学探测技术领域，本发明基于水下阵列信号，构建阵列接收的空域稀疏信号模型；构建超参数概率分布模型，通过贝叶斯公式获得待恢复稀疏信号的后验概率分布；构建超参数的目标函数，并在信号功率的目标函数中引入l0范数约束，通过期望最大化算法(EM)对超参数进行迭代更新；通过收敛得到的超参数重构空域稀疏信号，进而获得DOA估计结果迭代终止后，确定空间谱。在稀疏贝叶斯学习结构中，本发明仅加快了超参数γ的收敛速度，对信号模型无特殊需求。

Description

一种基于l0范数约束的稀疏贝叶斯DOA估计方法

技术领域

本发明涉及水下声学探测技术领域，是一种基于l0范数约束的稀疏贝叶斯DOA估计方法。

背景技术

波达方向(DOA)估计技术是通过阵列接收信号方向与阵列轴夹角来估测目标方位信息的技术或方法。作为目标跟踪，定位，识别和检测的前提基础，DOA估计技术在雷达，声纳以及通信领域被广泛应用。常见的DOA估计方法包括波束形成类方法，子空间类方法以及稀疏重构类方法。其中，稀疏重构类方法具有超高分辨率，同时对噪声、有限信号快拍及相关信号有较高鲁棒性，在DOA估计领域应用广泛。在稀疏重构类方法中，常见有贪婪类算法、凸优化类算法、迭代类算法以及稀疏贝叶斯学习(SBL)类算法。

传统的SBL类算法从贝叶斯角度将DOA估计问题转换成了稀疏重构问题，通过对目标信号进行稀疏先验假设来利用其空域稀疏性，具有优秀的DOA估计性能。但是由于SBL计算量庞大，无法保证DOA估计的实时性，因而难以在工程实践中应用。

发明内容

本发明在稀疏贝叶斯学习中引入l0范数约束，能够在保证估计估计精度的前提下实现快速DOA估计，本发明提供了一种基于l0范数约束的稀疏贝叶斯DOA估计方法。本发明提供了以下技术方案：

一种基于l0范数约束的稀疏贝叶斯DOA估计方法，包括以下步骤：

步骤1：基于水下阵列信号，构建阵列接收的空域稀疏信号模型；

步骤2：构建超参数概率分布模型，通过贝叶斯公式获得待恢复稀疏信号的后验概率分布；

步骤3：构建超参数的目标函数，并在信号功率的目标函数中引入l0范数约束，通过期望最大化算法对超参数进行迭代更新；

步骤4：通过收敛得到的超参数重构空域稀疏信号，进而获得DOA估计结果迭代终止后，确定空间谱。

优选地，所述步骤1具体为：

步骤1.1：将K个远场窄带信号s_k(t)，k＝1,…,K，分别从不同的方向θ＝[θ₁,θ₂,…,θ_K]入射到M阵元的均匀线列阵上，确定阵列接收信号模型如下：

y(t)＝A(θ)s(t)+e(t)，

其中，y(t)＝[y₁(t),y₂(t),…,y_M(t)]^T，s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T；

为服从

分布的加性高斯白噪声；A(θ)＝[a(θ₁),a(θ₂),…,a(θ_K)]为阵列流型，a(θ_k)为第k个目标的导向矢量；

步骤1.2：将阵列接收信号模型扩展为空域稀疏信号模型如下：

y(t)＝A(Θ)x(t)+e(t)，

其中，Θ＝[Θ₁,Θ₂,…,Θ_I]为DOA探测范围(0°,180°)离散后的采样网格，I为网格数目；当

为距离目标方位θ_k最近的网格点；

x(t)满足以下条件

步骤1.3：根据多快拍情况，阵列接收的空域稀疏信号多快拍模型通过下式表示：

Y＝A(Θ)X+E，

其中，Y＝[y(t₁),y(t₂),…,y(t_L)]，X＝[x(t₁),x(t₂),…,x(t_L)]，

E＝[e(t₁),e(t₂),…,e(t_L)]。

优选地，所述步骤2具体为：

当x(t)服从

其中，Γ＝diag(γ)为对角矩阵，对角元素来自于信号功率向量γ＝[γ₁,γ₂,…γ_I]；

当X矩阵中的所有列具有相同的稀疏支持，则X的概率密度函数通过下式表示：

p(X；Γ)＝|πΓ|^-Texp[-tr(X^HΓ^-1X)]，

在加性高斯白噪声假设下，似然概率密度函数通过下式表示：

通过贝叶斯公式推导，得到X的后验概率密度函数如下

均值

与方差矩阵Σ分别为：

Σ＝Γ-ΓA(Θ)^TΩ^-1A(Θ)Γ

且，Ω＝σ²I_M+A(Θ)ΓA(Θ)^H。

优选地，所述步骤3具体为：

步骤3.1：构建超参数的目标函数，并在信号功率的目标函数中引入l0范数约束，通过期望最大化算法对超参数进行迭代更新；

步骤3.2：将γ与σ²合称为超参数Δ＝{γ,σ²}，将X作为隐藏变量，得到超参数的目标函数通过下式表示：

通过将目标函数

关于σ²的梯度归零，确定σ²的更新公式：

l0范数是稀疏性的最佳度量，考虑超参数γ的稀疏性，在目标函数中加入l0范数通过下式表示：

其中，ε为正则化参数，用于平衡新的惩罚项与旧目标函数；将l0范数进行指数近似如下：

其中，δ为近似参数，用于使指数函数逼近l0范数；

将新的目标函数

关于γ的梯度归零，给出γ的更新公式：

步骤3.3：为了减小计算复杂度，γ的更新公式通过下式表示：

其中，

对超参数γ与σ²的更新公式进行迭代，直到迭代结果满足迭代终止条件。

优选地，在每一次更新后，X的后验概率均值u与方差矩阵Σ也相应更新，终止条件包括最大迭代次数N_iter与比例门限ζ，当迭代次数n大于N_iter，或ζ_th＝min(γⁿ)/max(γⁿ)小于ζ时终止迭代，其中γⁿ为第n次迭代后更新的超参数γ。

优选地，所述步骤4具体为：

迭代终止后，将均值

作为X的最优估计，在此基础上，计算空间谱如下：

P_X＝[p₁,p₂,…,p_I]

通过对空间谱搜峰，将峰值所对应网格位置作为目标方位估计结果θ。

本发明具有以下有益效果：

本发明在超参数γ的目标函数中引入了l0范数约束，相当于在γ的迭代更新中加入了零值吸引，加快了γ中的近零参数的收敛速度，进而在不损失DOA估计精度的前提下使得稀疏贝叶斯学习过程计算加速，提高了稀疏贝叶斯学习类算法的实时性。对所引入的l0范数进行了指数函数近似，使其导数具有简单的表达式，不会造成计算负担。在稀疏贝叶斯学习结构中，本发明仅加快了超参数γ的收敛速度，对信号模型无特殊需求。因此，其同时适用于在网模型与离网模型，易于移植到多种稀疏贝叶斯学习类方法中，具有较高的通用性。

附图说明

图1为基于l0范数约束的稀疏贝叶斯DOA估计方法流程图；

图2为收敛性能对比图；

图3为DOA估计精度对比图；

图4为空间谱对比图。

具体实施方式

以下结合具体实施例，对本发明进行了详细说明。

具体实施例一：

根据图1至图4所示，为解决上述技术问题采取的具体优化技术方案是：本发明涉及一种基于l0范数约束的稀疏贝叶斯DOA估计方法：

一种基于l0范数约束的稀疏贝叶斯DOA估计方法，包括以下步骤：

步骤1：基于水声信号中水下阵列信号，构建阵列接收的空域稀疏信号模型。考虑K个远场窄带信号s_k(t)，k＝1,…,K，分别从不同的方向θ＝[θ₁,θ₂,…,θ_K]入射到M阵元的均匀线列阵(ULA)上，给出阵列接收信号模型如下：

y(t)＝A(θ)s(t)+e(t)，

其中，y(t)＝[y₁(t),y₂(t),…,y_M(t)]^T，s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T。

为服从

分布的加性高斯白噪声。A(θ)＝[a(θ₁),a(θ₂),…,a(θ_K)]为阵列流型，a(θ_k)为导向矢量。

接着，将阵列接收信号模型扩展为空域稀疏信号模型：

y(t)＝A(Θ)x(t)+e(t)，

其中，Θ＝[Θ₁,Θ₂,…,Θ_I]为DOA探测范围(0°,180°)离散后的采样网格，I为网格数目。假设

为距离目标方位θ_k最近的网格点，x(t)满足以下条件

考虑多快拍情况，阵列接收的空域稀疏信号多快拍模型如下

Y＝A(Θ)X+E，

其中，Y＝[y(t₁),y(t₂),…,y(t_L)]，X＝[x(t₁),x(t₂),…,x(t_L)]，E＝[e(t₁),e(t₂),…,e(t_L)]。

步骤2：构建超参数概率分布模型，通过贝叶斯公式获得待恢复稀疏信号的后验概率分布。

假设X服从概率密度函数如下

p(X；Γ)＝|πΓ|^-Texp[-tr(X^HΓ^-1X)]，

其中，Γ＝diag(γ)且γ＝[γ₁,γ₂,…γ_I]。进一步地，在加性高斯白噪声假设下，似然概率密度函数为

通过贝叶斯公式，给出X的后验概率密度函数：

其中，均值

与方差矩阵Σ分别为

Σ＝Γ-ΓA(Θ)^TΩ^-1A(Θ)Γ

且，Ω＝σ²I_M+A(Θ)ΓA(Θ)^H。

步骤3：构建超参数的目标函数，并在信号功率的目标函数中引入l0范数约束，通过期望最大化算法(EM)对超参数进行迭代更新

将γ与σ²合称为超参数Δ＝{γ,σ²}，同时把X看作隐藏变量，可以得到超参数的目标函数如下

通过将目标函数

关于σ²的梯度归零，给出σ²的更新公式：

考虑超参数γ的稀疏性，在其目标函数中加入l0范数近似函数如下

式中ε为正则化参数，用于平衡新的惩罚项与旧目标函数。δ为近似参数，用于使指数函数逼近l0范数。将新的目标函数

关于γ的梯度归零，给出γ的更新公式：

其中，

对超参数γ与σ²的更新公式进行迭代，直到迭代结果满足迭代终止条件。在每一次超参数更新后，X的后验概率均值

与方差矩阵Σ也相应更新。

步骤4：通过收敛得到的超参数重构空域稀疏信号，进而获得DOA估计结果

迭代终止后，将

作为X的最优估计，在此基础上计算空间谱：

P_X＝[p₁,p₂,…,p_I]

通过对空间谱P_X进行搜峰，将峰值所对应网格位置作为目标方位估计结果θ。

采用仿真对本发明所设计的基于l0范数约束的稀疏贝叶斯DOA估计方法进行验证，并对结果进行说明。

在仿真结果中同时给出了本发明在在网模型下的应用(后文简称l0-SBL)与在离网模型下的应用(后文简称l0-OGSBI)，且分别与传统的SBL方法，离网稀疏贝叶斯学习(OGSBI)方法进行了性能对比。

仿真1用于验证本发明对超参数γ收敛速度的提升，实验结果由200次蒙特卡洛实验下统计得到的γ的均方根误差(RMSE-γ)给出。仿真各参数如下：水下布放8阵元等阵元间距均匀线列阵(ULA)，DOA探测范围被划分为网格集Θ＝[0°,1°,…,180°].考虑单个窄带远场目标信号入射，入射角度在DOA探测范围内随机生成，信号快拍数为200，信噪比为5dB.相关参数设置为δ＝8，ε＝1×10^-4，N_iter＝500和ζ＝1×10^-8。

各对比方法的RMSE-γ随迭代次数的变化如图2所示。分析图2可知，各对比方法中的超参数γ均可随迭代次数增加而逐渐收敛。本发明达到稳态所需的迭代次数明显少于传统SBL与OGSBI方法，验证了本发明对超参数γ收敛速度的提升。

仿真2用于验证本发明对DOA估计精度的影响，实验结果由200次蒙特卡洛实验下统计的，估计方位θ的均方根误差(RMSE-θ)给出。仿真各参数如下：水下布放8阵元等阵元间距均匀线列阵(ULA)，DOA探测范围被划分为网格集Θ＝[0°,1°,…,180°].考虑单个窄带远场目标信号入射，入射角度在DOA探测范围内随机生成，信号快拍数为200。信噪比由-5dB变化到10dB.相关参数设置为δ＝8，ε＝1×10^-4，N_iter＝500和ζ＝1×10^-8。

各对比方法的RMSE-θ随信噪比的变化如图3所示。分析图3可知，各对比方法的RMSE-θ均随信噪比增加而逐渐降低。本发明保持了与传统SBL，OGSBI相似的RMSE-θ，验证了本发明不会损失DOA估计精度。

仿真3各参数如下：水下布放8阵元等阵元间距均匀线列阵(ULA)，DOA探测范围被划分为网格集Θ＝[0°,1°,…,180°].考虑独立窄带远场双目标信号分别从θ₁＝61.25°和θ₂＝118.34°入射，信号快拍数为200，信噪比为5dB.相关参数设置为δ＝8，ε＝1×10^-4，N_iter＝500和ζ＝1×10^-4。

表1算法计算时间

各对比方法的空间谱结果如图4所示，表1给出了各算法的计算时间(单位：秒)。分析图4可知，本发明可得到与传统方法相似的空间谱，且本发明得到的空间谱拥有更尖锐的谱峰与更平坦的谱背景。分析表1可知，本发明的计算时间明显小于传统计算方法。由上可知，本发明在保证在优秀DOA估计性能的前提下具有更快的计算速度，提高了稀疏贝叶斯学习类算法的实时性。

以上所述仅是一种基于l0范数约束的稀疏贝叶斯DOA估计方法的优选实施方式，一种基于l0范数约束的稀疏贝叶斯DOA估计方法的保护范围并不仅局限于上述实施例，凡属于该思路下的技术方案均属于本发明的保护范围。应当指出，对于本领域的技术人员来说，在不脱离本发明原理前提下的若干改进和变化，这些改进和变化也应视为本发明的保护范围。

Claims (3)

1.一种基于l0范数约束的稀疏贝叶斯DOA估计方法，其特征是：包括以下步骤：

步骤1：基于水下阵列信号，构建阵列接收的空域稀疏信号模型；

所述步骤1具体为：

步骤1.1：将K个远场窄带信号s_k(t)，k＝1,…,K，分别从不同的方向θ＝[θ₁,θ₂,…,θ_K]入射到M阵元的均匀线列阵上，确定阵列接收信号模型如下：

y(t)＝A(θ)s(t)+e(t)，

其中，y(t)＝[y₁(t),y₂(t),…,y_M(t)]^T，s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T；

为服从

分布的加性高斯白噪声；A(θ)＝[a(θ₁),a(θ₂),…,a(θ_K)]为阵列流型，a(θ_k)为第k个目标的导向矢量；

步骤1.2：将阵列接收信号模型扩展为空域稀疏信号模型如下：

y(t)＝A(Θ)x(t)+e(t)，

其中，Θ＝[Θ₁,Θ₂,…,Θ_I]为DOA探测范围(0°,180°)离散后的采样网格，I为网格数目；当

为距离目标方位θ_k最近的网格点；

x(t)满足以下条件

步骤1.3：根据多快拍情况，阵列接收的空域稀疏信号多快拍模型通过下式表示：

Y＝A(Θ)X+E，

其中，Y＝[y(t₁),y(t₂),…,y(t_L)]，X＝[x(t₁),x(t₂),…,x(t_L)]，

E＝[e(t₁),e(t₂),…,e(t_L)]；

步骤2：构建超参数概率分布模型，通过贝叶斯公式获得待恢复稀疏信号的后验概率分布；

所述步骤2具体为：

令x(t)服从

其中，Γ＝diag(γ)为对角矩阵，对角元素来自于信号功率向量γ＝[γ₁,γ₂,…γ_I]；

当X矩阵中的所有列具有相同的稀疏支持，则X的概率密度函数通过下式表示：

p(X；Γ)＝|πΓ|^-Texp[-tr(X^HΓ^-1X)]，

在加性高斯白噪声假设下，似然概率密度函数通过下式表示：

通过贝叶斯公式推导，得到X的后验概率密度函数如下

均值

与方差矩阵Σ分别为：

Σ＝Γ-ΓA(Θ)^TΩ^-1A(Θ)Γ

且，Ω＝σ²I_M+A(Θ)ΓA(Θ)^H；

步骤3：构建超参数的目标函数，并在信号功率的目标函数中引入l0范数约束，通过期望最大化算法对超参数进行迭代更新；

所述步骤3具体为：

步骤3.1：构建超参数的目标函数，并在信号功率的目标函数中引入l0范数约束，通过期望最大化算法对超参数进行迭代更新；

步骤3.2：将γ与σ²合称为超参数Δ＝{γ,σ²}，将X作为隐藏变量，得到超参数的目标函数通过下式表示：

通过将目标函数

关于σ²的梯度归零，确定σ²的更新公式：

l0范数是稀疏性的最佳度量，考虑超参数γ的稀疏性，在目标函数中加入l0范数通过下式表示：

其中，ε为正则化参数，用于平衡新的惩罚项与旧目标函数；将l0范数进行指数近似如下：

其中，δ为近似参数，用于使指数函数逼近l0范数；

将新的目标函数

关于γ的梯度归零，给出γ的更新公式：

步骤3.3：为了减小计算复杂度，γ的更新公式通过下式表示：

其中，

对超参数γ与σ²的更新公式进行迭代，直到迭代结果满足迭代终止条件；

步骤4：通过收敛得到的超参数重构空域稀疏信号，进而获得DOA估计结果迭代终止后，确定空间谱。

2.根据权利要求1所述的一种基于l0范数约束的稀疏贝叶斯DOA估计方法，其特征是：在每一次更新后，X的后验概率均值

与方差矩阵Σ也相应更新，终止条件包括最大迭代次数N_iter与比例门限ζ，当迭代次数n大于N_iter，或ζ_th＝min(γⁿ)/max(γⁿ)小于ζ时终止迭代，其中γⁿ为第n次迭代后更新的超参数γ。

3.根据权利要求2所述的一种基于l0范数约束的稀疏贝叶斯DOA估计方法，其特征是：所述步骤4具体为：

迭代终止后，将均值

作为X的最优估计，在此基础上，计算空间谱如下：

P_X＝[p₁,p₂,…,p_I]

通过对空间谱搜峰，将峰值所对应网格位置作为目标方位估计结果θ。

Images