CN111175692B - 基于分层合成Lasso先验模型的离格稀疏贝叶斯DOA估计方法 - Google Patents

Abstract

本发明涉及一种基于分层合成Lasso先验模型的离格稀疏贝叶斯DOA估计方法，首先根据入射信号在空域内稀疏分布的特点建立稀疏信号模型，其次对模型中各变量进行概率假设构建HSL概率模型，继而采用变分推断算法得到概率模型中未知参数的更新公式，在赋予概率模型中未知参数初值后，以迭代更新的方式对阵列的输出数据进行处理，计算得到模型参数的估计值，最后利用得到的参数估计值，采用一维搜索的方法计算DOA估计值。

Description

基于分层合成Lasso先验模型的离格稀疏贝叶斯DOA估计方法

技术领域

本发明属于阵列信号处理领域，特别涉及一种基于分层合成Lasso先验模型的离格稀疏贝叶斯DOA估计方法。

背景技术

在波达方向(DOA,Direction of Arrival)估计领域中，稀疏贝叶斯学习(SBL,Sparse Bayesian Learning)类估计方法属于稀疏重构类估计方法，其具有超分辨、在低信噪比 和小快拍数条件下的估计结果较好、无需预设复杂的模型参数等优点。与经典的子空间类方法MUSIC以及范数优化类稀疏重构方法L1-SVD相比,SBL类估计方法具有更 好的估计性能，在近些年受到了广泛的关注与研究。SBL类估计方法是将DOA估计 问题转换为稀疏信号重构问题来求解，其首先在划分好的空域角度网格点集上建立稀 疏信号模型，再充分挖掘阵列输出数据的统计特性，构建稀疏诱导型概率模型，进而 采用期望-极大(EM,Expectation-Maximization)算法或变分推断算法得到概率模型中 未知参数的估计值。由于概率模型中的某些参数与预设的角度网格点集存在一一对应 的关系，选择一个合适的参数估计值，其极值点所对应的角度网格点即为DOA估计 值。

传统的SBL估计方法只在高斯先验分布假设的基础上(高斯先验、Laplace先验等)来构建稀疏诱导型概率模型，以保证实现稀疏信号的精确逆恢复，例如iRVM(Z. Liu,Z.Huang,and Y.Zhou,An efficient maximum likelihood method for direction-of-arrival estimation via sparse Bayesian learning,in IEEE Trans.WirelessCommun.,2012；11(10):1-11)。虽然基于高斯先验的概率模型能够描述稀疏信号的统计 特性，但其不是最优的。此外，注意到实际的信号入射角度可能并未落在预设的空域 角度网格点上，为了获得精确的DOA估计值，在建立稀疏信号模型时引入离格误差 (即真实入射角度偏离预设角度网格点的角度值)，并在后续参数估计中估计出离格误 差,例如OGSBI方法(Z.Yang,L.Xie,and C.Zhang,Off-grid direction of arrival estimation usingsparse Bayesian inference,in IEEE Trans.Signal Process,2013；61(1):38-43) 和SURE-IR方法(J.Fang,F.Wang,Y.Shen,H.Li,and R.S.Blum,Super-resolutioncompressed sensing for line spectral estimation:an iterative reweightedapproach,in IEEE Trans.Signal Process.,2016；64(18):649-4662)，但这两种方法所采用的信号先验分布在 统计学上不够稀疏，因此在分辨空域临近信源时存在性能局限。因此，为了提升SBL 类DOA估计方法的估计性能，构建怎样的稀疏诱导概率模型，以及当存在离格误差 时如何估计该离格误差是亟待解决的问题。

发明内容

要解决的技术问题

为了提高现有SBL类DOA估计方法的估计精度，本发明提出一种使用分层合成Lasso(HSL,hierarchical synthesis lasso)先验模型的离格稀疏贝叶斯DOA估计方法。

技术方案

1.一种基于分层合成Lasso先验模型的离格稀疏贝叶斯的DOA估计方法，其特征在于步骤如下：

步骤1：建立阵列输出数据的稀疏信号模型，具体内容如下：

使用M元均匀线阵对空间中K个分别以角度{θ₁,…,θ_K}入射的远场窄带信号进行采集，阵列输出M×L维的数据矩阵Y，其中L为采样快拍数；选定一个参考坐标，则该 均匀线阵中各阵元与参考坐标之间的距离间隔为{d₁,…,d_M}；接下来均匀划分空域角度 范围[-90°,90°]得到角度网格点集Θ＝[θ₁,…,θ_N]；定义N×L维的稀疏信号矩阵X，M×L维 的噪声信号矩阵N，M×N维的阵列流形A(Θ)，建立阵列输出数据Y的稀疏信号模型：

Y＝A(Θ)X+N

式中，

下面为了方便叙述，将阵列流形A(Θ)简记为A；

步骤2：建立基于HSL先验的概率模型，具体内容如下：

假设噪声服从均值为0方差为σ²的复高斯分布，阵列输出数据Y中第i列向量Y_·i服从均值为AX_·i、协方差矩阵为δ^-1I_M的复高斯分布，则Y的似然函数为多元复高斯分布， 其表达式为：

式中，δ＝σ^-2表示噪声精度，I_M表示M维单位矩阵，||·||₂表示向量的2范数；

假设δ服从参数为c,d的Gamma分布，其表达式为：

式中，

假设稀疏信号矩阵X中第i行向量X_i·服从Laplace先验分布，该分布的参数为η_i，记η＝[η₁,…,η_N]^T，[·]^T表示转置，则X的概率分布表示为：

式中，||·||₂表示向量的2范数；

假设

服从参数为a_i,b_i的Gamma分布，则η服从的概率分布可以表示为：

式中，a＝[a₁,…,a_N]^T，b＝[b₁,…,b_N]^T，[·]^T表示转置；

将X和η的概率分布联合起来并对η求积分，得到参数为a,b的HSL先验分布，其 表达式为：

式中，a＝[a₁,…,a_N]^T，b＝[b₁,…,b_N]^T，||·||₂表示向量的2范数，[·]^T表示转置；

步骤3：使用变分推断算法得到概率模型中各参数的更新公式，具体内容如下：

稀疏信号矩阵X第i列向量X_·i,i＝1,…,L的均值

和协方差矩阵

参数向量η中第i个元素η_i,i＝1,…,N的更新公式分别为：

式中，diag(·)表示对角矩阵，E(·)表示求期望，||·||₂表示向量的2范数，[·]^H表示共轭 转置，[·]^-1表示求逆，υ和ε_i,i＝1,…,N表示变分参数，其计算公式分别为：

式中，{·}_ii表示取括号内矩阵的第i行第i列元素，||·||₂表示向量的2范数，[·]^T表示 转置，[·]^H表示共轭转置，[·]^-1表示求逆；

噪声精度δ的更新公式为：

式中，Tr(·)表示矩阵的迹，||·||₂表示向量的2范数；

步骤4：根据所得到的参数更新公式对数据Y进行处理，得到概率模型参数的估 计值，具体内容如下：

对概率模型中的参数ε,X,η,δ赋初值，同时令c＝d＝10^-6，将阵列输出数据Y代入所得参数ε,X,η,δ的更新公式中，计算其更新值，不断重复这个过程直至满足预设的更新 停止条件，得到最终的参数估计值

步骤5：根据所得到的参数估计值，使用一维搜索方法计算得到DOA的估计值， 所涉及的具体内容如下：

按照如下表达式分别计算一维搜索方法中的三个中间计算变量值：

C＝δ^-1I-Adiag^-1(υ)A^H

式中，C_-i表示矩阵C＝δ-¹I-Adiag^-1(υ)A^H剔除第i列向量后的矩阵，diag(·)表示对 角矩阵，[·]^H表示共轭转置，[·]^-1表示求逆；

根据所得的参数估计值的极大值确定一维搜索的搜索角度范围Ω_i，下标i表示极大 值的序号，则在角度范围Ω_i中搜索得到的DOA估计值的表达式为：

式中，

表示在θ∈Ω_i范围内求极大值，Re(·)表示取复数的实部，[·]^H表示共轭 转置，[·]-¹表示求逆，[·]′表示求导。

有益效果

本发明首先根据入射信号在空域内稀疏分布的特点建立稀疏信号模型，其次对模型中各变量进行概率假设构建HSL概率模型，继而采用变分推断算法得到概率模型中 未知参数的更新公式，在赋予概率模型中未知参数初值后，以迭代更新的方式对阵列 的输出数据进行处理，计算得到模型参数的估计值，最后利用得到的参数估计值，采 用一维搜索的方法计算DOA估计值。当存在离格误差时，通过所设计的一维搜索方 法计算得到的DOA估计值包含有离格误差，而不需要在建立稀疏信号模型和概率模 型时引入离格误差，进而需要在参数估计过程中对离格误差进行估计。本发明所提方 法不仅构建了更优的稀疏诱导型概率模型，而且设计了一个一维搜索过程用于精确计 算DOA的估计结果，获得优于传统SBL类DOA估计方法的估计性能。

本发明的基本原理和实施方案经过了计算机数值仿真的验证，其结果表明：与其他五种DOA估计方法SURE-IR、OGSBI、iRVM,、L1-SVD和MUSIC相比，所提方 法获得更好的DOA估计性能。

附图说明

图1为实施实例中分别利用本发明所提方法与五种经典方法获得的估计结果均方根误差RMSE随信噪比SNR变化的曲线，以及其与克拉美罗下界曲线的对比图。

图2为实施实例中分别利用本发明所提方法与其他五种经典方法获得的估计结果均方根误差RMSE随快拍数变化的曲线，以及其与克拉美罗下界曲线的对比图。

图3为实施实例中分别利用本发明所提方法与其他五种经典方法获得的估计结果均方根误差RMSE随信号入射角度间隔变化的曲线，以及其与克拉美罗下界曲线的对 比图。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

参见图1-图3，本发明的主要内容有：

1.使用M元均匀线阵进行信号采集获得阵列输出数据，建立稀疏信号模型。

2.根据所建立的稀疏信号模型，构造基于HSL先验的概率模型。

3.使用变分推断算法计算得到概率模型中分布参数的估计值，并对所得到的参数估 计值使用一维搜索方法计算得到DOA的估计值。

4.通过计算机数值仿真分别给出了传统方法和本发明所提方法的DOA估计性能对比图，以此证明了本发明所提方法相比传统方法具有更好的DOA估计性能。

结合说明书附图，本发明所采用的技术方案可分为以下5个步骤：

步骤一：使用M元均匀线阵对空间中的K个远场窄带入射信号进行采集，得到阵 列的输出数据Y。其中这K个入射信号的入射角度分别为{θ₁,…,θ_K}，阵列输出数据Y 是一个M×L维矩阵，L是采样快拍数，数据Y中含有入射信号信息与噪声信号信息。

步骤二：建立阵列输出数据的稀疏信号模型，包括以下子步骤：

子步骤一：以角度间隔Δθ均匀划分空域角度范围[-90°,90°]，得到N个离散的角度 网格点组成空域角度网格点集Θ＝[θ₁,…,θ_N]；

子步骤二：定义N×L维的稀疏信号矩阵X，M×L维的噪声矩阵N，M×N维的阵列流形A(Θ)，则稀疏信号模型可表示为：

Y＝A(Θ)X+N

式中，

其中，a(θ_i),i＝1,…,N表示定义在入射角度θ_i方向的导向矢量，{d₁,…,d_M}表示阵列中各 阵元与参考位置之间的距离，λ表示窄带入射信号的中心频率；下面为了方便叙述， 将阵列流形A(Θ)简记为A。

步骤三：根据上一步建立的稀疏信号模型，建立基于HSL先验的概率模型，包括 以下子步骤：

子步骤一：假设噪声为高斯白噪声，且其服从均值为0方差为σ²的复高斯分布； 设定阵列输出数据Y的各列向量Y_·i,i＝1,…,L独立同分布，服从均值为AX_·i、协方差矩阵 为δ^-1I_M的复高斯分布，记为

则阵列输出数据Y的似然函数可以表示 为一个多元复高斯分布：

式中，δ＝σ^-2表示噪声精度，I_M表示M维单位矩阵，||·||₂表示向量的2范数；

进一步，定义δ服从Gamma分布：

式中，c>0为Gamma分布的形状参数，d>0为Gamma分布的尺度参数，

子步骤二：定义稀疏信号矩阵X的每一行向量X_i·,i＝1,…,N独立同分布且均服从Laplace先验分布，则X的概率分布可以表示为：

式中，η＝[η₁,…,η_N]^T表示分布参数向量，其元素η_i,i＝1,…,N为行向量X_i·的分布参数，[·]^T表示转置；

进一步，定义

服从参数为a_i,b_i的Gamma分布，则η服从的概率分布可以表示为：

式中，a＝[a₁,…,a_N]^T，b＝[b₁,…,b_N]^T，[·]^T表示转置；

子步骤三：利用以上定义的各概率分布，通过积分消除对X分布中参数向量η， 得到HSL先验分布：

式中，a＝[a₁,…,a_N]^T，b＝[b₁,…,b_N]^T，||·||₂表示向量的2范数，[·]^T表示转置。

步骤四：利用变分推断算法计算HSL概率模型中的参数，包括以下子步骤：

子步骤一：定义两个变分参数υ和ε_i,i＝1,…,N，其计算公式分别为：

式中，diag(·)表示对角矩阵，{·}_ii表示括号内矩阵的第i行第i列元素，||·||₂表示向量 的2范数，[·]^T表示转置，[·]^H表示共轭转置，[·]^-1表示求逆；

子步骤二：利用变分推断算法，分别得到稀疏信号矩阵X中列向量X_.i,i＝1,…,L的均 值

和协方差矩阵

参数向量η中各元素η_i,i＝1,…,N、噪声精度δ的更新公式：

式中，diag(·)表示对角矩阵，E(·)表示求期望，Tr(·)表示矩阵的迹，||·||₂表示向量的 2范数，[·]^H表示共轭转置，[·]^-1表示求逆；

子步骤三：对ε,X,η,δ赋初值，并选择c,d等于10^-6，将阵列输出数据Y代入上一步得到的公式对ε,X,η,δ进行计算更新，重复这个计算更新过程直至收敛，得到最终的估 计值

步骤五：根据步骤四中得到的参数估计值，通过一维搜索的方式计算得到DOA估计值，第i个DOA估计值的表达式为：

式中，

C_-i表示矩阵C＝δ^-1I-Adiag^-1(υ)A^H剔除第i列向量 后的矩阵，将s_i,q_ij,C称作该一维搜索方法的中间计算变量；

表示在θ∈Ω_i范围内 求极大值的运算，Ω_i表示搜索角度范围，其是根据在步骤四中得到的参数估计值的第i个极大值确定的，Re(·)表示取复数的实部，diag(·)表示对角矩阵，[·]^H表示共轭转置， [·]^-1表示求逆，[·]′表示求导。

利用计算机进行数值仿真，检验本发明所提方法的估计性能。

仿真所使用的均匀线阵的阵元个数为10，阵元之间的间隔等于入射信号的半波长。 以1°为间隔，均匀划分空域角度范围[-90°,90°]得到离散的角度网格点集。为了比较本发 明所提方法和经典方法的估计性能差异，选择SURE-IR、OGSBI、iRVM,、L1-SVD 和MUSIC方法作为比较对象，分别使用本发明所提方法和这五种经典方法在不同的 仿真条件下进行DOA估计，每种方法均进行200次独立实验，比较不同方法在所设 定的条件下获得的DOA估计结果的均方根误差RMSE。分别绘制这六种方法估计结 果的RMSE随信噪比、快拍数和信号入射角度间隔变化的曲线，并与克拉美罗下界 CRLB进行比较。RMSE的计算公式如下：

式中，

和θ⁽ⁱ⁾分别表示第i次实验中真实的DOA值和估计得到的DOA值，I表 示独立实验的次数，K表示入射信号个数，||·||₂表示向量的2范数，[·]^1/2表示开平方。

1)比较本发明所提方法和其他五种方法在不同信噪比条件下的估计性能

假设存在3个非相干窄带入射信号，其分别以角度-10°+ζ、3°+ζ和11°+ζ入射到均匀线阵上，其中ζ为[-0.5°,0.5°]范围内的随机角度。设置快拍数L＝50，3个入射信号 的信噪比SNR相等，均从-5dB增大到20dB，分别计算本发明所提方法和其他五种方 法在不同SNR条件下的RMSE以及CRLB。图1为在不同SNR的条件下，不同方法 的RMSE曲线和CRLB曲线的对比图。从图1中可知，所有方法的RMSE随着SNR 的增加而减小，即所有方法的DOA估计精度均随着SNR的增加而改善。在低SNR 条件下，除开SURE-IR和MUSIC方法，其余四种稀疏重构类方法均可以成功分辨3 个入射信号的DOA。本发明所提方法在不同SNR条件下的RMSE均是最小的，并且 与CRLB曲线的值非常接近，这表明所提方法在不同SNR条件的DOA估计精度是最好的。

2)比较本发明所提方法和其他五种方法在不同快拍数条件下的估计性能

在上一实验的基础上，固定SNR＝0dB，设置快拍数L由20增大至100，分别计算 本发明所提方法和其他五种方法在不同快拍数条件下的RMSE以及CRLB。图2为在 不同快拍数的条件下，不同方法的RMSE曲线和CRLB曲线的对比图。从图2中可知， 所有方法的RMSE基本上均随着快拍数的增大而减小，但是MUSIC方法在快拍数为 50左右时才能够将这三个入射信号分辨开来，并且其DOA估计精度要低于本发明所 提方法。对于SURE-IR方法来说，即便快拍数达到100，其也无法成功分辨这三个信 号。从图2中明显可以看出，本发明所提方法在不同快拍数条件下的DOA估计精度 是所有方法中最高的。

3)比较本发明所提方法和其他五种方法在不同信号入射角度间隔下的估计性能

考虑两个非相干窄带信号入射到均匀线阵上，信号的SNR均为0dB，快拍数L＝50。设置信号入射角度间隔的入射角度间隔的绝对值|θ₂-θ₁|从5°增大到15°，这两个信号的 入射角度分别为-|θ₂-θ₁|/2+ζ和|θ₂-θ₁|/2+ζ，其中ζ的设置与之前的实验相同。分别计 算本发明所提方法和其他五种方法在不同信号入射角度间隔条件下的RMSE以及 CRLB。图3为在不同信号入射角度间隔的条件下，不同方法的RMSE曲线和CRLB 曲线的对比图。从图3中可知，

当两个信号的入射角度间隔小于10°时，SURE-IR方法已不能分辨入射信号，其余四种 经典方法可分辨的信号入射角度间隔阈值分别为6°(L1-SVD)、7°(OGSBI)、8°(iRVM)和9°(MUSIC)，而本发明所提方法的可分辨的信号入射角度间隔阈值分别为6°，此时 所提方法的DOA估计精度接近CRLB，并且本发明所提方法的DOA估计精度是所有 方法中最高的。

Claims (1)

1.一种基于分层合成Lasso先验模型的离格稀疏贝叶斯的DOA估计方法，其特征在于步骤如下：

步骤1：建立阵列输出数据的稀疏信号模型，具体内容如下：

使用M元均匀线阵对空间中K个分别以角度{θ₁,…,θ_K}入射的远场窄带信号进行采集，阵列输出M×L维的数据矩阵Y，其中L为采样快拍数；选定一个参考坐标，则该均匀线阵中各阵元与参考坐标之间的距离间隔为{d₁,…,d_M}；接下来均匀划分空域角度范围[-90°,90°]得到角度网格点集Θ＝[θ₁,…,θ_N]；定义N×L维的稀疏信号矩阵X，M×L维的噪声信号矩阵N，M×N维的阵列流形A(Θ)，建立阵列输出数据Y的稀疏信号模型：

Y＝A(Θ)X+N

式中，

下面为了方便叙述，将阵列流形A(Θ)简记为A；

步骤2：建立基于HSL先验的概率模型，具体内容如下：

假设噪声服从均值为0方差为σ²的复高斯分布，阵列输出数据Y中第i列向量Y_.i服从均值为AX_.i、协方差矩阵为δ^-1I_M的复高斯分布，则Y的似然函数为多元复高斯分布，其表达式为：

式中，δ＝σ^-2表示噪声精度，I_M表示M维单位矩阵，||·||₂表示向量的2范数；

假设δ服从参数为c,d的Gamma分布，其表达式为：

式中，

假设稀疏信号矩阵X中第i行向量X_i·服从Laplace先验分布，该分布的参数为η_i，记η＝[η₁,…,η_N]^T，[·]^T表示转置，则X的概率分布表示为：

式中，||·||₂表示向量的2范数；

假设

服从参数为a_i,b_i的Gamma分布，则η服从的概率分布可以表示为：

式中，a＝[a₁,…,a_N]^T，b＝[b₁,…,b_N]^T，[·]^T表示转置；

将X和η的概率分布联合起来并对η求积分，得到参数为a,b的HSL先验分布，其表达式为：

式中，a＝[a₁,…,a_N]^T，b＝[b₁,…,b_N]^T，||·||₂表示向量的2范数，[·]^T表示转置；

步骤3：使用变分推断算法得到概率模型中各参数的更新公式，具体内容如下：

稀疏信号矩阵X第i列向量X_.i,i＝1,…,L的均值

和协方差矩阵

参数向量η中第i个元素η_i,i＝1,…,N的更新公式分别为：

式中，diag(·)表示对角矩阵，E(·)表示求期望，||·||₂表示向量的2范数，[·]^H表示共轭转置，[·]^-1表示求逆，υ和ε_i,i＝1,…,N表示变分参数，其计算公式分别为：

式中，{·}_ii表示取括号内矩阵的第i行第i列元素，||·||₂表示向量的2范数，[·]^T表示转置，[·]^H表示共轭转置，[·]^-1表示求逆；

噪声精度δ的更新公式为：

式中，Tr(·)表示矩阵的迹，||·||₂表示向量的2范数；

步骤4：根据所得到的参数更新公式对数据Y进行处理，得到概率模型参数的估计值，具体内容如下：

对概率模型中的参数ε,X,η,δ赋初值，同时令c＝d＝10^-6，将阵列输出数据Y代入所得参数ε,X,η,δ的更新公式中，计算其更新值，不断重复这个过程直至满足预设的更新停止条件，得到最终的参数估计值

步骤5：根据所得到的参数估计值，使用一维搜索方法计算得到DOA的估计值，所涉及的具体内容如下：

按照如下表达式分别计算一维搜索方法中的三个中间计算变量值：

C＝δ^-1I-Adiag^-1(υ)A^H

式中，C_-i表示矩阵C＝δ^-1I-Adiag^-1(υ)A^H剔除第i列向量后的矩阵，diag(·)表示对角矩阵，[·]^H表示共轭转置，[·]^-1表示求逆；

根据所得的参数估计值的极大值确定一维搜索的搜索角度范围Ω_i，下标i表示极大值的序号，则在角度范围Ω_i中搜索得到的DOA估计值的表达式为：

式中，

表示在θ∈Ω_i范围内求极大值，Re(·)表示取复数的实部，[·]^H表示共轭转置，[·]^-1表示求逆，[·]′表示求导。

Images