CN109298383A - 一种基于变分贝叶斯推断的互质阵波达方向角估计方法 - Google Patents

Abstract

本发明提供了一种基于变分贝叶斯推断的互质阵波达方向角估计方法，利用互质阵对空间信号进行接收采样，根据采样协方差矩阵构造稀疏模型向量，将波达角估计问题转换为线性模型中的稀疏重构问题，绘制幅度谱图，从幅度谱图中找出峰值点所对应的x轴坐标即为所求的波达方向角估计值。本发明避免了传统方法的角度搜索或角度匹配过程，同时采样频率可以低于奈奎斯特采样率，在低快拍数情况下得到较好的估计性能，大大降低了模型维度，在运算量和稀疏重构速度上得到改善，避免了复杂分布形式的直接求解，降低算法复杂度，提高算法收敛速度。

Description

一种基于变分贝叶斯推断的互质阵波达方向角估计方法

技术领域

本发明涉及信号处理领域，特别涉及一种互质阵波达方向角估计方法。

背景技术

信号的波达方向角DOA估计技术是阵列信号处理领域的一个重要内容，它是指利用传感器阵列接收空间信号，再运用现代信号处理方法快速准确地估计出信号源的入射方向，已在包括雷达、声呐、无线通信等众多领域得到广泛应用。

传统波达方向角估计方法主要使用的阵列结构是均匀线列阵，可分辨的目标数受限于阵列的物理阵元数。稀疏阵列能够在阵元数目相同时，获得比传统均匀线列阵更大的孔径，降低了经济成本，因而受到广泛关注。近年来，P.P.Vaidyanathan(P.Pal andP.P.Vaidyanathan.Coprime sampling and the music algorithm[C].IEEE,in DigitalSignal Processing Workshop and IEEE Signal Processing Education Workshop(DSP/SPE),2011. 289-294.)提出的互质阵受到广泛关注。互质阵通常由两个共线均匀线阵组成，阵元数和阵元间距分别为M和Nd，N和Md(其中M，N为互质数，d为半波长)，该阵元数为M+N-1的互质阵最多可分辨O(MN)个信源。

互质阵波达方向角估计方法大致可以划分为：1)基于子空间类方法，2)基于稀疏表示的方法。以基于空间平滑的多重信号分类SS-MUSIC为代表的子空间类方法在低信噪比和小快拍数情形下性能严重下降。近年来出现的基于稀疏表示的波达方向角估计方法利用信号的空域稀疏性进行建模，以贪婪算法、凸优化方法和基于稀疏贝叶斯学习的算法为主要手段而展开目标参数求解。其中贪婪类算法在低信噪比情况下估计性能大幅下降；而凸优化方法运算量大，在低信噪比情况下估计精度不理想。基于稀疏贝叶斯学习的算法对未知参数指定先验概率分布，通过最大后验概率估计获得目标参数的优化解。该算法通常采用期望-最大EM算法进行迭代求解目标参数，但EM 算法不适用于先验分布计算困难的情况，参数初始值影响迭代得到的局部最大值，且具有计算复杂度较高、收敛速度较慢等问题。

发明内容

为了克服现有技术的不足，本发明提供一种基于变分贝叶斯推断的互质阵波达方向角估计方法，以提高在低信噪比、小快拍数情况下波达角度估计的精度以及算法的收敛速度。

本发明解决其技术问题所采用的技术方案包括如下步骤：

步骤1：采用M个传感器按照位置集合{d₁,…,d_M}布放形成互质阵列，并假设有K个独立窄带信号从远场入射到该互质阵列，其中M＝M₁+2M₂-1，且M₁和M₂互为质数， {d₁,…,d_M}＝{M₂m₁d,0≤m₁≤M₁-1}∪{M₁m₂d,1≤m₂≤2M₂-1}，0＜d≤λ/2，λ为入射窄带信号波长，m₁和m₂分别为布放位置间隔为M₂d和M₁d的传感器的编号；

步骤2：利用互质阵对空间信号进行接收采样，得到输出信号y(t)，对该输出信号作时间平均，计算采样协方差矩阵R：

其中，N为采样快拍数，H表示共轭转置运算；

步骤3：根据采样协方差矩阵R构造稀疏模型向量z，构造步骤如下：

步骤3a：构建一个M(M-1)×M²维的选择矩阵J：

J＝[J₁,…,J_m,…,J_M-1]^T

其中，J_m＝[e_(m-1)(M+1)+2,e_(m-1)(M+1)+3,…,e_i,…,e_m(M+1)]是一个M²×M维矩阵，m为一个整数，且1≤m≤M-1，e_i是一个仅第i个元素为1其余元素为0的M²×1维列向量；

步骤3b：构建加权矩阵：

其中，T表示转置运算，表示Kronecker积，N为采样快拍数；

步骤3c：向量化采样协方差矩阵R，根据选择矩阵J和加权矩阵W，得到稀疏模型向量z：

z＝W^-1/2Jvec(R)

其中，vec(·)表示向量化运算，W-^1/2为W的逆的厄尔米特平方根；

步骤4：对观测空间进行网格划分，构造超完备基B，将波达角估计问题转换为线性模型中的稀疏重构问题，详细步骤如下：

步骤4a：根据目标信号的空域稀疏特性，将观测空域[-90°,90°]以r度等间隔划分成Q 个角度，得到空域离散角度θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]，θ_q为可能的目标信号来波方向角，q＝1,…,Q，Q＞＞M；

步骤4b：构造一个对应空域角度稀疏化的M×Q维的阵列流形矩阵A(θ)：

A(θ)＝[a(θ,),…a(θ,,),…,a(θ_O)]

其中，a(θ_q)表示对应角度θ_q的导向矢量，且

其中，j为虚数单位；

步骤4c：根据阵列流形矩阵A(θ)和步骤3中得到的选择矩阵J和加权矩阵W，得到超完备基B：

B＝W^-1/2J[A^*(θ)⊙A(θ)]

其中，B＝[b(θ₁),…,b(θ_q),…,b(θ_Q)]，q＝1,2,…,Q，b(θ_q)称为基向量；

步骤5：根据步骤3和步骤4得到的结果，将波达方向角估计问题转化为求解如下稀疏方程：

z＝Bp+ε

其中，p＝[p₁,…,p_q,…,p_Q]^T，是一个Q×1维的未知向量，ε为估计误差；

步骤6：定义两个超参数向量γ＝[γ₁,…,γ_q,…,γ_Q]^T，β＝[β₁,…β_q,…,β_Q]^T，γ_q为控制p分布的未知先验方差，β_q为控制γ_q分布的未知参数，均称为超参数，定义一个超参数δ为控制z条件概率分布的未知方差的逆，δ服从参数为(c,d)的伽马先验分布，采用变分贝叶斯推断求解步骤5中的稀疏方程，得到超参数γ_q、β_q,q＝1,…,Q和δ的估计值的更新公式，以及未知向量p的收敛解，其中，收敛条件为超参数γ_q相邻两次迭代的变化量小于阈值即为收敛，阈值取值为10^-4；

步骤7：以波达方向角范围θ＝[θ₁,…,θ_Q]的值为x轴坐标，以步骤6中得到的向量p的收敛解的幅度值为y轴坐标，绘制幅度谱图，从幅度谱图中按照幅度从高到低的顺序排序，寻找幅度最大的前K个谱峰，前K个谱峰的峰值点所对应的x轴坐标即为所求的波达方向角估计值。

本发明的有益效果在于：

1)充分利用入射信号源的空域稀疏特性结合稀疏表示进行建模，将波达方向角估计问题转化为稀疏重构问题，避免了传统方法的角度搜索或角度匹配过程。同时采样频率可以低于奈奎斯特采样率，在低快拍数情况下得到较好的估计性能。

2)采用矢量化建模，将多测量矢量模型转变为单测量矢量模型，大大降低了模型维度，在运算量和稀疏重构速度上得到改善。

3)采用分层先验模型描述未知参数的统计特性，利用变分贝叶斯学习算法获得参数的后验概率分布的近似分布，避免了复杂分布形式的直接求解，降低算法复杂度，提高算法收敛速度。

附图说明

图1是本发明与基于空间平滑的多重信号分类SS-MUSIC方法的幅度谱图作对比；

图2是本发明与现有五种波达方向角估计方法在不同信噪比条件下的均方根误差对比图；

图3是本发明与现有五种波达方向角估计方法在不同快拍数条件下的均方根误差对比图；

图4是本发明与现有五种波达方向角估计方法在不同波达方向角度间隔条件下的均方根误差对比图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

本发明所用的M元互质阵由两个传感器个数分别为M₁和2M₂-1的均匀线列阵组成，所有传感器共线放置，传感器位置的坐标集合分别为{M₂m₁d,0≤m₁≤M₁-1}和 {M₁m₂d,1≤m₂≤2M₂-1}。假定空间K个远场独立窄带信号入射到该互质阵上，且在信号传播过程中加有均值为零的高斯白噪声。其中，M₁和M₂互为质数，K≥1，0＜d≤λ/2，λ为入射窄带信号波长。

参照图1，本发明的具体实现步骤如下：

步骤1：计算互质阵接收数据的采样协方差矩阵R。

用互质阵列的M个接收传感器以固定的采样频率对空间信号进行采样，采样点数为N，得到输出信号y(t)，对输出信号y(t)的时间平均，计算互质阵列的采样协方差矩阵R：

其中，N表示采样快拍数，H表示共轭转置运算。

步骤2：根据采样协方差矩阵R构造稀疏模型向量z。

向量化采样协方差矩阵R得到采样协方差向量，对采样协方差向量进行预处理，剔除噪声分量，正则化估计误差使其分布服从渐进标准正态分布，以获得无噪、仅含有限信源分量的正则化表示的稀疏模型向量z，具体预处理步骤如下：

(2a)构造一个M(M-1)×M²维的选择矩阵J，用以剔除采样协方差向量中的噪声分量：

J＝[J₁,…,J_m,…,J_M-1]^T，

其中，J_m＝[e_(m-1)(M+1)+2,…,e_i,…,e_m(M+1)]，e_i是一个仅第i个元素为1其余元素为0的M²×1 维列向量，T表示转置运算。

(2b)构造权矩阵W，用以正则化采样协方差向量的估计误差：

其中，表示Kronecker积。

(2c)向量化采样协方差矩阵R，根据选择矩阵J和权矩阵W，得到稀疏模型向量z：

z＝W^-1/2Jvec(R)，

其中，W-^1/2为W-¹的厄尔米特平方根，vec(·)表示向量化运算。

步骤3：对观测空间离散网格化，构造超完备基B：

基于入射信号空域分布的稀疏性，离散化观测模型的参数空间，将波达角估计问题转换为现行模型中的稀疏重构问题。根据稀疏信号重构理论，任意信号都可以由一个基矩阵线性表示，构造超完备基B矩阵用以线性表示步骤2中得到的稀疏模型向量z，具体构造步骤如下：

(3a)根据入射信号的空域稀疏性，将观测空域[-90°,90°]以r度等间隔划分成Q个角度，得到空域离散角度θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]，θ_q为可能的目标信号来波方向角，q＝1,…,Q， Q＞＞M。

(3b)构造一个对应空域角度稀疏化的M×Q维的阵列流形矩阵A(θ)：

A(θ)＝[a(θ_,),…a(θ_q)，…a(θ_Q)，

其中，a(θ_q)表示对应角度θ_q的导向矢量：

其中，j为虚数单位。

(3c)根据阵列流形矩阵A(θ)和步骤2中得到的选择矩阵J和加权矩阵W，得到超完备基B：

B＝W^-1/2J[A^*(θ)⊙A(θ)]，

其中，B＝[b(θ₁),…,b(θ_q),…,b(θ_Q)]，b(θ_q)称为基向量，q＝1,2,…,Q，⊙⊙表示Khatri-Rao 积。

步骤4：将稀疏模型向量z用超完备基B线性表示，将波达角估计问题转化为稀疏信号重构问题，求解如下稀疏方程：

z＝Bp+ε

其中，p＝[p₁,…,p_q,…,p_Q]^T是一个Q×1维的未知稀疏向量，p中非零元素所对应基向量的角度即为所求的波达方向角，ε为估计误差。

步骤5：采用变分贝叶斯推断算法对上述稀疏方程中的未知稀疏向量p进行参数估计，得到其参数估计的收敛解。

对未知的稀疏向量p指定层先验概率分布，定义一个超参数向量γ＝[γ₁,…,γ_q,…,γ_Q]^T表示空域离散角度θ₁,…,θ_q,…,θ_Q上入射信号功率，p服从零均值协方差矩阵为Γ＝diag(γ)的高斯先验分布。定义一个参数向量β＝[β₁,…β_q,…,β_Q]^T，其第q个元素β_q为控制超参数向量γ第 q个元素γ_q服从的伽马县衙分布的未知参数，且β_q均服从参数为c、d伽马先验分布，其中 q＝1,…,Q。定义一个参数δ，对稀疏模型向量z指定均值为Bp协方差为δ^-1I_M(M-1)条件高斯概率分布。γ，β，δ称为超参数。设置各参数初始值，利用变分贝叶斯学习算法，以迭代更新的方式得到超参数γ_q、β_q和δ的估计值，进而得到稀疏向量p的最大后验概率估计。一旦确定了p的收敛解中非零元素的位置，即可得到入射信号的波达角估计值，其求解步骤如下：

(5a)初始化，p⁽⁰⁾＝(B^HB)^-1B^Hz，参数a、b、c、d设置为很小的值，如a＝b＝c＝d＝10^-6。其中，上角标(i),i＝0,1,…表示第i次迭代，(·)^-1表示矩阵求逆，|·|表示求绝对值，||·||表示求向量2范数。

(5b)根据超完备基B和稀疏模型向量z，计算未知向量p⁽ⁱ⁾的均值向量μ⁽ⁱ⁾和协方差矩阵

μ⁽ⁱ⁾＝Γ^(i-1)B^H((δ^-1)⁽ⁱ⁾I_M(M-1)+BΓ^(i-1)B^H)^-1z，

其中，Γ＝diag(γ)，I_M(M-1)为M(M-1)阶单位矩阵，(·)^-1表示矩阵求逆，diag(·)表示构造对角矩阵。

(5c)根据变分贝叶斯学习准则，更新超参数和δ⁽ⁱ⁾：

其中，μ_q为均值向量μ的第q个元素，(·)^*表示共轭运算，||·||表示向量2范数，tr(·)表示矩阵的迹。

(5d)判断收敛条件||<p⁽ⁱ⁾>-<p^(i-1)>||≤10^-4是否满足，若满足，则停止迭代更新；若不满足，则返回上一步继续迭代更新计算。

步骤6：根据得到的p的收敛解绘制幅度谱图，得到波达方向角的估计值。

以波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]的值为x轴坐标，以收敛解向量p的幅度值为y 轴坐标，绘制幅度谱图，从谱图中按照从高到低的顺序找到前K个幅值较大的谱峰，这些峰值点所对应的x轴坐标点即为所求的入射信号波达方向角。

本发明的效果可以通过以下仿真说明：

1.仿真条件与方法：

采用6个传感器形成互质阵，传感器的位置集合为{0,2,3,4,6,9}倍的入射信号半波长，观测空域角度范围为[-90°,90°]，空间网格离散化间隔为1°。

2.仿真内容与结果：

仿真1：假设有7个独立窄带信号分别以角度-50°、-35°、-10°、5°、20°、30°和-45°入射到互质阵，信噪比SNR为0dB，快拍数为100，利用本发明和现有的SS-MUSIC算法分别进行波达方向角估计试验，结果如图1所示，图中横坐标表示观测角度值，纵坐标表示归一化的幅度谱值。

从图1中可以看出，在相同条件下，本发明能够正确估计入射信号的波达方向，且本发明具有比SS-MUSIC算法更尖锐的谱峰。

仿真2：假设有2个独立窄带信号分别以角度-3°和3°入射到互质阵，信噪比由-10dB 增加到20dB，快拍数为100，利用本发明和现有的L1-SRACV、SS-MUSIC、SBL和SPICE 算法分别进行500次独立的波达方向角估计试验，分别计算不同信噪比条件下六种方法的均方根误差RMSE，均方误差计算公式如下：

其中，L表示试验次数，L＝200，表示第l次试验的DOA估计值，θ_k表示对应的真实入射信号波达方向角。试验结果及克拉美罗下界曲线如图2所示，图中横坐标表示信噪比值，单位dB，纵坐标表示均方根误差，单位为度。

从图2中可以看出，本发明在不同信噪比下的均方根误差小于其他算法，在低信噪比情况下更加接近克拉美罗下界。

仿真3：假设有2个独立窄带信号分别以角度-3°和3°入射到互质阵，信噪比为0dB，快拍数从50到500变化，利用本发明和现有的L1-SRACV、SS-MUSIC、SBL和SPICE算法分别进行500次独立的波达方向角估计试验，分别计算不同信噪比条件下六种方法的均方根误差RMSE，试验结果及克拉美罗下界曲线如图3所示，图中横坐标表示采样快拍数，纵坐标表示均方根误差。

从图3中可以看出，本发明在不同快拍数条件下的均方根误差均小于其他算法，且在小快拍数下也有良好的估计性能。

仿真4：假设有2个独立窄带信号的入射到互质阵，两个信号的入射角度间隔从3°到13°变化，信噪比为0dB，快拍数为100，利用本发明和现有的L1-SRACV、SS-MUSIC、SBL 和SPICE算法分别进行500次独立的波达方向角估计试验，分别计算不同信噪比条件下六种方法的均方根误差RMSE，试验结果及克拉美罗下界曲线如图4所示，图中横坐标表示两个入射信号方向角的角度间隔，纵坐标表示均方根误差。

从图4中可以看出，本发明在不同入射角度间隔下的均方根误差小于其他算法，能分辨出以小角度间隔入射的信号。

综上所述，本发明降低了DOA估计的运算量，降低了在低信噪比、小快拍情况下的角度估计的估计误差，具有更好的小角度目标分辨性能。

Claims (1)

1.一种基于变分贝叶斯推断的互质阵波达方向角估计方法，其特征在于包括下述步

骤：

步骤1：采用M个传感器按照位置集合{d₁,…,d_M}布放形成互质阵列，并假设有K个独立窄带信号从远场入射到该互质阵列，其中M＝M₁+2M₂-1，且M₁和M₂互为质数，{d₁,…,d_M}＝{M₂m₁d,0≤m₁≤M₁-1}∪{M₁m₂d,1≤m₂≤2M₂-1}，0＜d≤λ/2，λ为入射窄带信号波长，m₁和m₂分别为布放位置间隔为M₂d和M₁d的传感器的编号；

步骤2：利用互质阵对空间信号进行接收采样，得到输出信号y(t)，对该输出信号作时间平均，计算采样协方差矩阵R：

其中，N为采样快拍数，H表示共轭转置运算；

步骤3：根据采样协方差矩阵R构造稀疏模型向量z，构造步骤如下：

步骤3a：构建一个M(M-1)×M²维的选择矩阵J：

J＝[J₁,…,J_m,…，Jm，…，J,J_M-1]^T

其中，J_m＝[e_(m-1)(M+1)+2,e_(m-1)(M+1)+3,…,e_i,…,e_m(M+1)]是一个M²×M维矩阵，m为一个整数，且1≤m≤M-1，e_i是一个仅第i个元素为1其余元素为0的M²×1维列向量；

步骤3b：构建加权矩阵：

其中，T表示转置运算，表示Kronecker积，N为采样快拍数；

步骤3c：向量化采样协方差矩阵R，根据选择矩阵J和加权矩阵W，得到稀疏模型向量z：

z＝W-^1/2Jvee(R)

其中，vec(·)表示向量化运算，W^-1/2为W的逆的厄尔米特平方根；

步骤4：对观测空间进行网格划分，构造超完备基B，将波达角估计问题转换为线性模型中的稀疏重构问题，详细步骤如下：

步骤4a：根据目标信号的空域稀疏特性，将观测空域[-90°,90°]以r度等间隔划分成Q个角度，得到空域离散角度θ＝[θ₁,θ₂,…,θ_q,…，θ_q，…，θ_q,θ_Q]，θ_q为可能的目标信号来波方向角，q＝1,…,Q，Q＞＞M；

步骤4b：构造一个对应空域角度稀疏化的M×Q维的阵列流形矩阵A(θ)：

A(θ)＝[a(θ₁)，…a(θ_q)，…，a(θ_Q)]

其中，a(θ_q)表示对应角度θ_q的导向矢量，且

其中，j为虚数单位；

步骤4c：根据阵列流形矩阵A(θ)和步骤3中得到的选择矩阵J和加权矩阵W，得到超完备基B：

B＝W^-1/2J[A^*(θ)⊙A(θ)]

其中，B＝[b(θ₁),…,b(θ_q),…,b(θ_Q)]，q＝1,2,…,Q，b(θ_q)称为基向量；

步骤5：根据步骤3和步骤4得到的结果，将波达方向角估计问题转化为求解如下稀疏方程：

z＝Bp+ε

其中，p＝[p₁,…,p_q,…,p_Q]^T，是一个Q×1维的未知向量，ε为估计误差；

步骤6：定义两个超参数向量γ＝[γ₁,…,γ_q,…,γ_Q]^T，β＝[β₁,…β_q,…,β_Q]^T，γ_q为控制p分布的未知先验方差，β_q为控制γ_q分布的未知参数，均称为超参数，定义一个超参数δ为控制z条件概率分布的未知方差的逆，δ服从参数为(c,d)的伽马先验分布，采用变分贝叶斯推断求解步骤5中的稀疏方程，得到超参数γ_q、β_q,q＝1,…,Q和δ的估计值的更新公式，以及未知向量p的收敛解，其中，收敛条件为超参数γ_q相邻两次迭代的变化量小于阈值即为收敛，阈值取值为10^-4；

步骤7：以波达方向角范围θ＝[θ₁,…,θ_Q]的值为x轴坐标，以步骤6中得到的向量p的收敛解的幅度值为y轴坐标，绘制幅度谱图，从幅度谱图中按照幅度从高到低的顺序排序，寻找幅度最大的前K个谱峰，前K个谱峰的峰值点所对应的x轴坐标即为所求的波达方向角估计值。