CN109407046A - 一种基于变分贝叶斯推断的嵌套阵列波达方向角估计方法 - Google Patents

Abstract

本发明提供了一种基于变分贝叶斯推断的嵌套阵列波达方向角估计方法，采用三层先验分布模型，并利用变分推断寻找后验概率分布的近似分布，求的入射信号波达方向角。本发明采用了嵌套阵阵型结构进行波达方向角估计，克服了现有技术中采用均匀直线阵时可分辨目标数低于阵元数的缺点，提高了在阵元数相同条件下可分辨的目标个数，基于稀疏重构理论，使得待估计稀疏向量的后验概率分布更趋于平缓，更符合信号的稀疏特性，避免了后验概率的直接求解，特别适合后验概率复杂、难以求解的情况，且降低了运算复杂度。

Description

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

技术领域

本发明涉及信号处理技术领域，尤其是一种阵列信号波达方向角估计方法，可用于对信号到达方向角，俯仰角等多种参数的准确估计，从而获得信号源位置。

背景技术

信号的波达方向角DOA估计是阵列信号处理领域的一个重要内容，它利用按一定方式在空间布放的传感器组成的阵列对空间声信号感应接收，再运用现代信号处理方法快速准确地估计信号的波达方向，俯仰角等多种参数。高分辨DOA估计在雷达、声纳、无线通信等领域具有重要应用价值。

目前，超分辨DOA估计方法可以分为子空间类算法和稀疏重构类算法。以多重信号分类MUSIC算法为代表的子空间类DOA估计方法，在低信噪比和小快拍数情况下估计性能严重下降。稀疏重构类算法可以分为凸优化算法和基于稀疏贝叶斯学习SBL的 DOA估计算法。其中，凸优化算法运算复杂度大，在低信噪比下估计精度差；传统的基于SBL准则的DOA估计算法对待估计稀疏向量作高斯-伽马先验两层先验概率分布假设，并通过期望极大EM准则以迭代的方式得到最优最优估计，算法复杂度较高，收敛速度慢。

现有大部分超分辨DOA估计算法聚焦于使用均匀线列阵，一个M元均匀线列阵最多可分辨M-1个目标信号，在目标个数很多的情况下无法满足估计需要。为了突破阵元数给定情况下，均匀线列阵阵型对最大可分辨目标数的限制，一些新的阵型结构相继被提出，其中最具代表性的是嵌套阵列和互质阵列。Piya Pal等人在其发表的论文“Nested arrays：anovel approach to array processing with enhanced degrees of freedom” (《IEEETransactions on Signal Processing》，VOL 58，No.8，August，2010)中提出了嵌套阵列结构及其适用的DOA估计算法——空间平滑多重信号分类SS-MUSIC算法。使用阵元数为M+N嵌套阵可以分辨MN+N-1较好的估计性能，需保证足够的采样快拍数，较高的输入信噪比，入射信号源数目已知。

发明内容

为了克服现有技术的不足，本发明提出一种基于变分贝叶斯推断的嵌套阵波达方向角估计方法。在贝叶斯估计理论框架下，采用三层先验分布模型，并利用变分推断寻找后验概率分布的近似分布，降低算法复杂度，提高估计精度和算法收敛速度。

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

步骤1：使用M个接收传感器形成两级嵌套阵列，详细步骤如下：

步骤1a：用M₁个传感器以间隔d水平布放形成第一均匀直线阵，将每一个传感器称为一个阵元，以第一均匀直线阵的首阵元作为起始阵元，其中，M₁≥1，0＜d≤λ/2，λ为入射窄带信号的波长；

步骤1b：用M-M₁个传感器以间隔(M₁+1)d水平布放形成第二均匀直线阵，其首阵元放置于距参考阵元M₁d位置处，将第一均匀直线阵和第二均匀直线阵共线布放组成两级嵌套阵，其中M≥2；

步骤2：假设有K个窄带信号从远场入射到两级嵌套阵上，使用两级嵌套阵对空间信号接收采样，得到输出信号x(t)，计算采样协方差矩阵R：

其中，∑(·)表示求和运算，上标H表示共轭转置运算，N为采样快拍数；

步骤3：构造选择矩阵J，根据采样协方差矩阵R，得到无噪协方差向量z，详细步骤如下：

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

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

其中，j＝1，2，…，M-1，J_j＝[e_(m-1)(M+1)+2，…，e_i，…，e_m(M+1)]为一个M²×M维矩阵，e_i是一个仅第i个元素为1其余元素为0的M²×1维向量，上标T表示转置运算；

步骤3b：向量化采样协方差矩阵R，根据选择矩阵J，得到无噪协方差向量z为：

z＝Jvec(R)

其中，vec(·)表示矩阵向量化运算；

步骤4：网格化观测空间，构造超完备基B，详细步骤如下：

步骤4a：将观测空间角度[-90°，90°]均匀划分为G个角度，得到观测空间网格点θ＝[θ₁，θ₂，…，θ_G]，其中，G＞＞M²/2+M-1＞K；

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

A(θ)＝[a(θ₁)，…，a(θ_g)，…，a(θ_G)]

其中，a(θ_g)为θ_g角度对应的导向矢量，g＝1，2，…，G，且：

其中，d_m为嵌套阵第m个阵元布放位置的坐标，m＝1，2，…，M，两级嵌套阵所有阵元位置的坐标集合为{nd，n＝1，…，M₁}∪{nd(M₁+1)，n＝1，…，M-M₁}；

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

B＝J[A*(θ)⊙A(θ)]

其中，上标*表示共轭运算，⊙表示khatri-Rao积运算；

步骤5：根据步骤3和步骤4的结果，将波达方向角估计问题转化为稀疏重构问题，求解关于无噪协方差向量z稀疏矩阵方程：

z＝Bη+ε

其中，η为G×1维未知稀疏向量，ε为M(M-1)×1维估计误差向量；

步骤6：定义超参数γ＝[γ₁，γ₂，…，γ_G]^T、τ＝[τ₁，τ₂，…，τ_G]^T和∑，其中，γ和τ为控制稀疏向量η的三层先验概率分布的未知参数向量，∑为控制无噪采样协方差向量z的条件概率分布的未知参数矩阵；利用贝叶斯变分推断算法得到稀疏向量η以及超参数γ、τ和∑的更新公式，以迭代的方式得到γ的收敛解；

步骤7：以步骤6中得到的收敛解向量γ的幅度值作纵坐标，以观测空间网格点θ＝[θ₁，θ₂，…，θ_G]作横坐标，绘制幅度谱图，在幅度谱图中按幅值从大到小的顺序找到前 K个峰值，前K个峰值所对应的横坐标角度值即为所求的入射信号波达方向角。

本发明的有益效果在于：

1)本发明采用了嵌套阵阵型结构进行波达方向角估计，克服了现有技术中采用均匀直线阵时可分辨目标数低于阵元数的缺点，提高了在阵元数相同条件下可分辨的目标个数。

2)本发明基于稀疏重构理论，在贝叶斯估计框架下，对待估计稀疏向量采用三级分级先验概率分布模型，使其边际概率分布为拉普拉斯分布，与传统单超参数的高斯概率分布相比，拉普拉斯分布曲线更平缓，使得待估计稀疏向量的后验概率分布更趋于平缓，更符合信号的稀疏特性。

3)本发明采用贝叶斯变分推断算法求解波达方向角估计问题中的稀疏矩阵方程，与传统稀疏贝叶斯学习框架下所使用的期望极大EM算法相比，变分贝叶斯推断算法寻求后验概率的近似概率分布，避免了后验概率的直接求解，特别适合后验概率复杂、难以求解的情况，且降低了运算复杂度。

附图说明

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

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

具体实施方式

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

步骤1：使用M个接收传感器形成两级嵌套阵列，详细步骤如下：

步骤1a：用M₁个传感器以间隔d水平布放形成第一均匀直线阵，将每一个传感器称为一个阵元，以第一均匀直线阵的首阵元作为起始阵元，其中，M₁≥1，0＜d≤λ/2，λ为入射窄带信号的波长；

步骤1b：用M-M₁个传感器以间隔(M₁+1)d水平布放形成第二均匀直线阵，其首阵元放置于距参考阵元M₁d位置处，将第一均匀直线阵和第二均匀直线阵共线布放组成两级嵌套阵，其中M≥2；

步骤2：假设有K个窄带信号从远场入射到两级嵌套阵上，使用两级嵌套阵对空间信号接收采样，得到输出信号x(t)，计算采样协方差矩阵R：

其中，∑(·)表示求和运算，上标H表示共轭转置运算，N为采样快拍数；

步骤3：构造选择矩阵J，根据采样协方差矩阵R，得到无噪协方差向量z，详细步骤如下：

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

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

其中，j＝1，2，…，M-1，J_j＝[e_(m-1)(M+1)+2，…，e_i，…，e_m(M+1)]为一个M²×M维矩阵，e_i是一个仅第i个元素为1其余元素为0的M²×1维向量，上标T表示转置运算；

步骤3b：向量化采样协方差矩阵R，根据选择矩阵J，得到无噪协方差向量z为：

z＝Jvec(R)

其中，vec(·)表示矩阵向量化运算；

步骤4：网格化观测空间，构造超完备基B，详细步骤如下：

步骤4a：将观测空间角度[-90°，90°]均匀划分为G个角度，得到观测空间网格点θ＝[θ₁，θ₂，…，θ_G]，其中，G＞＞M²/2+M-1＞K；

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

A(θ)＝[a(θ₁)，…，a(θ_g)，…，a(θ_G)]

其中，a(θ_g)为θ_g角度对应的导向矢量，g＝1，2，…，G，且：

其中，d_m为嵌套阵第m个阵元布放位置的坐标，m＝1，2，…，M，两级嵌套阵所有阵元位置的坐标集合为{nd，n＝1，…，M₁}∪{nd(M₁+1)，n＝1，…，M-M₁}；

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

B＝J[A*(θ)⊙A(θ)]

其中，上标*表示共轭运算，⊙表示khatri-Rao积运算；

步骤5：根据步骤3和步骤4的结果，将波达方向角估计问题转化为稀疏重构问题，求解关于无噪协方差向量z稀疏矩阵方程：

z＝Bη+ε

其中，η为G×1维未知稀疏向量，ε为M(M-1)×1维估计误差向量；

步骤6：定义超参数γ＝[γ₁，γ₂，…，γ_G]^T、τ＝[τ₁，τ₂，…，τ_G]^T和∑，其中，γ和τ为控制稀疏向量η的三层先验概率分布的未知参数向量，∑为控制无噪采样协方差向量z的条件概率分布的未知参数矩阵；利用贝叶斯变分推断算法得到稀疏向量η以及超参数γ、τ和∑的更新公式，以迭代的方式得到γ的收敛解；

步骤7：以步骤6中得到的收敛解向量γ的幅度值作纵坐标，以观测空间网格点θ＝[θ₁，θ₂，…，θ_G]作横坐标，绘制幅度谱图，在幅度谱图中按幅值从大到小的顺序找到前K个峰值，前K个峰值所对应的横坐标角度值即为所求的入射信号波达方向角。

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

步骤1：使用M个传感器接收机形成一个两级嵌套阵列。

1a)用M₁个传感器以间隔d水平布放形成第一均匀直线阵，每一个传感器称为一个阵元，以第一均匀直线阵的首阵元作为起始阵元，其中，M₁≥1，0＜d≤λ/2，λ为入射窄带信号的波长。

1b)用M-M₁个传感器以间隔(M₁+1)d水平布放形成第二均匀直线阵，其首阵元放置于距参考阵元M₁d位置处，将两个均匀直线阵共线布放组成两级嵌套阵，其中 M≥2，0＜d≤λ/2，λ为入射窄带信号的波长。

步骤2：得到嵌套阵输出信号x(t)，并计算采样协方差矩阵R。

假定有K个窄带信号从远场入射到该嵌套阵上，且在信号传播过程中加入高斯白噪声，使用该嵌套阵对空间信号进行接收采样，得到输出信号x(t)，按照下式计算采样协方差矩阵R：

其中，∑(·)表示求和运算，上标H表示共轭转置运算，N为采样快拍数。

步骤3：构造选择矩阵J，根据采样协方差矩阵R，得到无噪采样协方差向量z。

向量化采样协方差矩阵R后的结果，可看作孔径扩展的虚拟差合阵的单快拍输出信号，加性噪声分量仅包含在采样协方差向量的特定位置的元素中，可以通过一个线性变换剔除含噪元素来消除噪声成分，得到无噪采样协方差向量z，具体实现步骤如下：

3a)构造一个M(M-1)×M²维的选择矩阵：

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

其中，j＝1，2，…，M-1，J_j＝[e_(m-1)(M+1)+2，…，e_i，…，e_m(M+1)]为一个M²×M维矩阵，e_i是一个仅第i个元素为1其余元素为0的M²×1维向量，上标T表示转置运算；

3b)向量化采样协方差矩阵，根据选择矩阵J，得到无噪协方差向量z：

z＝Jvec(R)

其中，vec(·)表示矩阵向量化运算。

步骤4：网格化观测空间，构造超完备基B。

根据入射信号在空间稀疏分布的特性，将观测空间角度离散网格化，基于稀疏表示的思想，构造对应空间稀疏化后的超完备基B用以线性表示无噪采样协方差向量，具体实现步骤如下：

4a)将观测空间角度[-90°，90°]以一定的角度间隔均匀划分为G个角度，得到观测空间网格点θ＝[θ₁，θ₂，…，θ_G]，其中，G＞＞M²/2+M-1＞K；

4b)构造一个对应空域稀疏化后的M²×G维的阵列流形A(θ)：

A(θ)＝[a(θ₁)，…，a(θ_g)，…，a(θ_G)]，

其中，a(θ_g)，g＝1，2，…，G为θ_g角度对应的导向矢量：

其中，d_m，m＝1，2，…，M为嵌套阵第m个阵元扽位置坐标，该嵌套阵的阵元位置坐标集合为{nd，n＝1，…，M₁}∪{nd(M₁+1)，n＝1，…，M-M₁}；

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

B＝J[A*(θ)⊙A(θ)]＝[b(θ₁)，b(θ₂)，…，b(θ_G)]，

其中，上标*表示共轭运算，⊙表示khatri-Rao积运算，b(θ_g)，g＝1，…，G称为基向量。

步骤5：根据步骤(3)和步骤(4)的结果，将波达方向角估计问题转化为稀疏重构问题，求解关于无噪协方差向量z稀疏矩阵方程：

z＝Bη+ε；

其中，η为G×1维未知稀疏向量，η中数值大的元素所对应基向量的角度即为所求波达方向角，ε为M(M-1)×1维估计误差向量。

步骤6：利用变分贝叶斯推断算法对稀疏向量η进行参数估计。

首先，对稀疏向量η采用三层分层先验分布模型：定义超参数向量γ＝[γ₁，γ₂，…，γ_G]^T，对η指定均值为0，协方差矩阵为Γ＝diag(γ)的高斯先验分布；定义超参数向量τ＝[τ₁，τ₂，…，τ_G]^T，对γ中元素γ_i，i＝1，2，…，G指定均值为的指数先验分布；对τ中各元素指定参数为a，b的伽马先验分布，其中，diag(·)表示构造对角矩阵运算，上标T表示转置运算，参数a，b的值很小。将η的三层分层先验分布对γ_i，i＝1，2，…，G进行积分后得拉普拉斯概率分布。

然后，定义一个维的超参数矩阵∑，根据步骤5中的稀疏方程，无噪采样协方差向量z服从高斯条件概率分布，其均值为Bη，协方差矩阵为∑；对∑指定逆威沙特概率分布，其分布的自由度个数为v，其维的标量矩阵为Ψ，其中，

最后，利用贝叶斯变分推断算法，得到η、γ、τ和∑的更新公式，在赋予各超参数初值后，进行迭代计算，当满足收敛条件||γ^(q+1)-γ^(q)||₂/||γ^(q)||≤ξ或达到最大迭代次数，终止迭代过程，此时γ^(new)的值为最优估计值；否则，继续进行迭代过程。其中，||·||₂表示向量2范数运算，上标q表示第q次迭代，q＝1，2，…，ξ表示判决门限，可取10^-4，η、γ、τ和∑的更新公式如下：

6a)稀疏向量η的均值μ和协方差矩阵∑_η的更新公式分别为

μ^(q+1)＝Γ^(q)B^H(∑^(q)+BΓ^(q)B^H)^-1z，

其中，Γ^(q)＝diag(γ^(q))，(·)^-1表示矩阵求逆运算，diag(·)表示对角矩阵；

6b)超参数向量γ中各元素γ_i，i＝1，2，…，G的更新公式为

其中，K(·)表示二阶改进贝塞尔函数，μ_i表示均值μ的第i个元素，τ_i表示超参数向量τ的第i个元素；

6c)超参数向量τ中各元素τ_i，i＝1，2，…，G的更新公式为

其中，γ_i表示超参数向量γ的第i个元素，

6d)超参数矩阵∑更新公式为

其中，上标H表示共轭转置运算。

步骤7：根据超参数向量γ的最优估计值，绘制幅度谱图，得到所求波达方向角。

以观测空间网格点θ＝[θ₁，θ₂，…，θ_G]为横坐标(单位为度)，以γ的最优估计的各元素幅值取10log₁₀(γ_i)，i＝1，2，…，G为纵坐标(单位为分贝dB，)，绘制幅度谱图。从该谱图中按照从大到小的顺序得到前K个峰值，这些峰值所对应的横坐标角度值即为所求的入射信号波达方向角。

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

1.仿真条件：

采用6元两级嵌套阵，组成该嵌套阵的两个均匀直线阵各有3个阵元，其阵元位置坐标为[0，1，2，3，7，11]λ/2，观测空间角度为[-90°，90°]，空间网格化划分间隔为1°。

波达方向角估计的均方根误差RMSE的计算公式如下：

其中，J表示实验次数，J＝500，为第j次实验中第k个入射信号的波达方向角估计值，θ_k为第k个入射信号的真实波达方向角。

2.仿真内容与结果：

仿真1：假设有2个远场窄带独立信号，分别以波达方向角-4°和4°入射到该嵌套阵上，采样快拍数为100，信噪比SNR从-5dB到20dB变化。采用本发明和现有SS-MUSIC、 CSA、SPICE、SBL和L1-SRACV算法分别进行500次独立的波达方向角估计实验DOA 估计，分别计算不同信噪比条件下各算法估计结果的均方根误差，得到均方根误差- 信噪比曲线与克拉美罗下界如图1所示。此外，在信噪比SNR为0dB时，分别计算每种方法的平均运算时间(单位为秒)，结果如表1所示：

表1

算法	CSA	L1-SRACV	SS-MUSIC	SBL	SPICE	本发明
							平均运算时间/s	0.5079	0.3789	0.0026	1.4957	0.0521	0.1140

其中，表1是本发明与现有五种波达方向角算法的平均运算时间对比表，图1中横坐标为信噪比，纵坐标为均方根误差。

从图1可以看出，本发明比其他算法在低信噪比条件下具有更好的估计性能，同时也更接近克拉美罗下界。

从表1可以看出，除SPICE算法外，本发明比其他稀疏重构类算法运算效率更高，尽管SS-MUSIC和SPICE算法运算时间更短，但它们的估计性能从图3中可知差于本发明，因此本算法更具实际意义。

仿真2：在仿真1的基础上，固定信噪比SNR为0dB，将采样快拍数从20到400变化，采用本发明和现有SS-MUSIC、CSA、SPICE、SBL和L1-SRACV算法分别进行500次独立的波达方向角估计实验DOA估计，分别计算不同快拍数条件下各算法的均方根误差，得到均方根误差-快拍数曲线如图2所示，图中横坐标为快拍数，纵坐标为均方根误差。

从图2中可以看出，本发明在小快拍数条件下与其他算法相比具有更好的估计性能。

Claims (1)

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

步骤1：使用M个接收传感器形成两级嵌套阵列，详细步骤如下：

步骤1a：用M₁个传感器以间隔d水平布放形成第一均匀直线阵，将每一个传感器称为一个阵元，以第一均匀直线阵的首阵元作为起始阵元，其中，M₁≥1，0＜d≤λ/2，λ为入射窄带信号的波长；

步骤1b：用M-M₁个传感器以间隔(M₁+1)d水平布放形成第二均匀直线阵，其首阵元放置于距参考阵元M₁d位置处，将第一均匀直线阵和第二均匀直线阵共线布放组成两级嵌套阵，其中M≥2；

步骤2：假设有K个窄带信号从远场入射到两级嵌套阵上，使用两级嵌套阵对空间信号接收采样，得到输出信号x(t)，计算采样协方差矩阵R：

其中，∑(·)表示求和运算，上标H表示共轭转置运算，N为采样快拍数；

步骤3：构造选择矩阵J，根据采样协方差矩阵R，得到无噪协方差向量z，详细步骤如下：

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

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

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

步骤3b：向量化采样协方差矩阵R，根据选择矩阵J，得到无噪协方差向量z为：

z＝Jvec(R)

其中，vec(·)表示矩阵向量化运算；

步骤4：网格化观测空间，构造超完备基B，详细步骤如下：

步骤4a：将观测空间角度[-90°,90°]均匀划分为G个角度，得到观测空间网格点θ＝[θ₁,θ₂,…,θ_G]，其中，G＞＞M²/2+M-1＞K；

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

A(θ)＝[a(θ₁)，…，a(θ_g)，…，a(θ_G)]

其中，a(θ_g)为θ_g角度对应的导向矢量，g＝1,2,…,G，且：

其中，d_m为嵌套阵第m个阵元布放位置的坐标，m＝1,2,…,M，两级嵌套阵所有阵元位置的坐标集合为{nd,n＝1,…,M₁}∪{nd(M₁+1),n＝1,…,M-M₁}；

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

B＝J[A^*(θ)⊙A(θ)]

其中，上标*表示共轭运算，⊙表示khatri-Rao积运算；

步骤5：根据步骤3和步骤4的结果，将波达方向角估计问题转化为稀疏重构问题，求解关于无噪协方差向量z稀疏矩阵方程：

z＝Bη+ε

其中，η为G×1维未知稀疏向量，ε为M(M-1)×1维估计误差向量；

步骤6：定义超参数γ＝[γ₁,γ₂,…,γ_G]^T、τ＝[τ₁,τ₂,…,τ_G]^T和Σ，其中，γ和τ为控制稀疏向量η的三层先验概率分布的未知参数向量，Σ为控制无噪采样协方差向量z的条件概率分布的未知参数矩阵；利用贝叶斯变分推断算法得到稀疏向量η以及超参数γ、τ和Σ的更新公式，以迭代的方式得到γ的收敛解；

步骤7：以步骤6中得到的收敛解向量γ的幅度值作纵坐标，以观测空间网格点θ＝[θ₁,θ₂,…,θ_G]作横坐标，绘制幅度谱图，在幅度谱图中按幅值从大到小的顺序找到前K个峰值，前K个峰值所对应的横坐标角度值即为所求的入射信号波达方向角。