CN110208735B - 一种基于稀疏贝叶斯学习的相干信号doa估计方法 - Google Patents

Abstract

本发明提供了一种基于稀疏贝叶斯学习的相干信号DOA估计方法，获取接收阵列的输出信号，网格化观测空间，构造超完备阵列流形，结合稀疏表示的思想，将DOA估计问题转化为稀疏信号重构问题，求解如下稀疏矩阵方程，建立稀疏贝叶斯概率模型，采用变分贝叶斯推断方法计算各隐变量的近似后验分布，计算入射信号的DOA估计值。本发明无需预先估计入射信号个数，且不涉及解相干操作，有效地实现相干信号的DOA估计，针对相干信号，所提DOA估计方法的角度分辨能力好，与现有的其它DOA估计方法相比，该方法具有更好的测向精度。

Description

一种基于稀疏贝叶斯学习的相干信号DOA估计方法

技术领域

本发明涉及信号处理技术领域，尤其是一种相干信号DOA的估计方法。

背景技术

波达方向(DOA，Direction-of-Arrival)估计是阵列信号处理领域中的核心任务之一，该技术对由空间中按照某种特定方式放置的传感器阵列接收到的来波信号进行处理，从而估计出目标信号的个数、来波方向等信号参数，被广泛应用于雷达、声纳、无线通信等领域。DOA估计问题与谐波恢复有着密切联系，而谐波恢复问题可以通过参数化技术来解决，特别是以多重信号分类(MUSIC，Multiple Signal Classification)算法为代表的具有高角度分辨性能的子空间类算法。在利用MUSIC算法进行DOA估计时，对阵列输出数据的协方差矩阵进行特征值分解，得到信号子空间和噪声子空间，利用二者的正交性构造空间方位谱函数，以角度搜索的方式得到DOA估计结果。然而，子空间算法需要预知目标信号个数，在快拍数较少和信噪比较低的情况下，算法的估计性能出现严重下降，并且不能用于处理相干信号(协方差矩阵出现秩亏损，无法正确分离信号子空间和噪声子空间)。若使用子空间类算法处理相干信号，需要进行解相干操作。

不同于子空间类算法，基于稀疏贝叶斯学习的DOA估计算法将DOA估计问题转换为稀疏信号重构问题，在贝叶斯估计框架下对入射信号的统计特性进行估计，进而获得DOA估计结果。基于稀疏贝叶斯学习的DOA估计算法不需要提前预知目标信号个数，且不涉及能够处理相干信号。在实际应用中，由于信号的多径传播特性，阵列接收到的信号不再是单一的不相关信号，因此，研究基于稀疏贝叶斯学习的高分辨、稳健的相干信号DOA估计具有重要的应用价值。

发明内容

为了克服现有技术的不足，本发明提供一种稀疏贝叶斯学习框架下的相干信号DOA估计方法。本发明可以解决现有子空间类DOA估计算法在处理相干信号时需要预知目标个数以及解相干操作这一问题。

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

步骤一：获取接收阵列的输出信号Y；

设定使用M个全指向性传感器组成接收阵列，且假定空间中有N个远场窄带相干信号，分别以角度θ_n入射到接收阵列上，其中，n＝1,2,…,N，利用接收阵列对入射信号进行接收采样，阵列的输出信号为：

Y＝[y(1),y(2),…y(L)] (1)

其中，y(t)(t＝1,…,L)表示阵列在t时刻的输出信号，L表示快拍数；

步骤二：网格化观测空间，构造超完备阵列流形A；

将观测空间角度在范围[-90°,90°]内以1°的角度间隔均匀划分，得到角度网格点集Θ＝{θ₁,…,θ_K}，其中K为网格点总数，且K＞＞N；根据角度网格点集Θ，构造阵列流形：

A＝[a(θ₁)，a(θ₂)，…，a(θ_k)，…，a(θ_K)] (2)

其中，

表示在网格点θ_k上的导向矢量，简记为a_k，d_m为接收阵列中第m个传感器的位置坐标，m＝1,…,M，λ为入射相干信号的波长，j为虚数单位；

步骤三：结合稀疏表示的思想，将DOA估计问题转化为稀疏信号重构问题，求解如下稀疏矩阵方程：

Y＝AX+n (4)

其中，X表示K×L维的信号矩阵，n表示M×L维的加性高斯白噪声矩阵；由于X中仅有N个非零行向量，X为稀疏矩阵；

步骤四：建立稀疏贝叶斯概率模型；

首先，对阵列输出信号Y的每一列向量进行复高斯分布假设，则Y的似然函数表示为：

其中，Y_·i,X_·i分别表示矩阵Y,X的第i列向量，I_M表示单位矩阵，β>0表示噪声精度(precision)，噪声精度是指噪声方差的倒数，对β进行伽马先验分布假设，即：

p(β)＝Gamma(β|c，d) (6)

其中，c,d为伽马分布的参数；

接下来，对信号矩阵X构造分层稀疏先验；

在第一层先验中，引入由隐变量组成的D×L维矩阵Z、K×D维矩阵W及K维向量μ，将X的每列进行线性变换X_·i＝WZ_·i+μ+ε,(i＝1,…,L)，其中Z_·i表示矩阵Z的第i列向量，ε表示0均值的复高斯向量，则X服从如下复高斯分布：

其中，Λ＝diag(γ)，γ＝[γ₁,…,γ_K]^T，diag(·)表示生成对角矩阵运算；超参数γ含有预设网格点方向上入射信号的功率信息；采用公式(7)所示线性高斯先验模型可将信号X的自相关分量包含于矩阵Λ^-1中，互相关分量包含于矩阵W中；

在第二层先验中，分别对超参数W,Z,μ,γ进行先验假设，假设W,Z的每一列以及μ均服从零均值复高斯分布，γ的每个元素独立同分布，均服从伽马分布，即：

其中，W_·i表示W的第i列向量，γ_i表示γ的第i个元素，α＝[α₁,…,α_D]^T，a,b为伽马分布的参数，δ为精度参数，(·)^-1表示矩阵求逆运算；

在第三层先验中，对超参数α的每一个元素进行伽马先验假设：

其中，g,h为伽马分布的参数；

设置参数初值a₀＝b₀＝c₀＝d₀＝g₀＝h₀＝10^-6，δ₀＝10^-3；记ξ＝{X，W，Z，μ，α，β，γ}，称作隐变量集；

步骤五：采用变分贝叶斯推断方法计算各隐变量ξ_i的近似后验分布q(ξ_i)，得到：

①

即X中第i列向量X_·i的后验分布为复高斯分布，其均值

和方差Σ_X分别为：

其中，Ψ＝Λ^-1＝(diag(γ))^-1，Y_i,Z_i分别表示Y,Z的第i列向量，<·>表示求期望运算，(·)^H表示矩阵的共轭转置运算；

②

即W中第j列向量

的后验分布为复高斯分布，其均值

和方差

分别为：

其中，X_ji表示X的第j行第i列元素，γ_j,μ_j分别表示γ,μ的第j个元素；

③

即Z中第i列向量Z_·i的后验分布为复高斯分布，其均值

和方差Σ_Z分别为：

Σ_Z＝<I_D+W^HΛW>^-1 (18)

④

即μ的后验分布为复高斯分布，其均值μ_μ和方差Σ_μ分别为：

Σ_μ＝<LΛ+δI_K>^-1 (20)

⑤

即α中第i个元素α_i的后验分布为伽马分布,该分布的参数g,h为：

g＝g₀+K (21)

其中，||·||₂表示向量的2范数运算；

⑥q(β)＝Gamma(β|c,d)，即β的后验分布为伽马分布，该分布的参数c,d分别为：

c＝c₀+LM (23)

其中，tr(·)表示矩阵的迹；

⑦

即γ中第j个元素γ_j的后验分布为伽马分布，该分布的参数a,b分别为：

a＝a₀+L (25)

其中，W_j·表示W的第j列向量，|·|表示绝对值；

根据{X,W,Z,μ,α,β,γ}的更新公式(13)～(26)，在设置隐变量初值后，进行隐变量的迭代更新直至满足收敛条件，则停止迭代；隐变量初值X₍₀₎＝A^H(AA^H)^-1Y,

W₍₀₎＝1_K×M,Z₍₀₎＝1_M×L,μ₍₀₎＝1_K×1,α₍₀₎＝1_M×1,γ₍₀₎＝1_K×1,其中(·)_(r)表示第r步迭代中的变量，

表示n₁×n₂维的全1矩阵，收敛条件为：||γ_(r)-γ_(r-1)||₂≤10^-4；

步骤六：计算入射信号的DOA估计值；

根据步骤五中得到的各隐变量估计值，以γ为纵坐标Θ为横坐标绘制空间方位谱图，在谱图中峰值点所在的网格点角度区域内，通过一维搜索的方式计算得到入射信号的DOA估计值

计算公式如下：

其中，

表示导向矢量a_j对θ_j的偏导值，Re(·)表示取复数的实部，θ_i表示空间方位谱图中第i个峰值点所在的网格点角度区域。

由于超参数γ含有入射信号的功率信息，空间方位谱图中峰值点所对应的网格点即为入射信号DOA估计值；然而，当入射信号没有落在预设网格点上，即网格失配时，根据谱图峰值点所确定的DOA估计值存在较大的误差，需要进行进一步的精确计算。构造目标函数：

其中，

表示

中与γ_j无关的部分，

表示

中与γ_j有关的部分；将

对γ_j求偏导并令结果为0，得到：

其中，

将

对θ_j求偏导并令结果为0，即：

将式(29)代入式(30)中得到：

根据式(31)即可得到DOA估计值的计算公式(27)。

本发明的有益效果在于由于本发明无需预先估计入射信号个数，且不涉及解相干操作，有效地实现相干信号的DOA估计。本发明采用了多层先验模型，使用非零均值复高斯分布对目标信号的空域稀疏性和信号之间的相干性进行描述，有利于相干信号的稀疏重构。针对相干信号，所提DOA估计方法的角度分辨能力好，与现有的其它DOA估计方法相比，该方法具有更好的测向精度。本发明也能应用于独立信号和混合信号(相干+不相关)的DOA估计问题。

附图说明

图1是本发明DOA估计结果空间方位谱图。

图2是本发明与现有五种DOA估计方法在不同信噪比条件下对两组相干信号DOA估计结果的均方根误差对比图。

图3是本发明与现有五种DOA估计方法在不同快拍数条件下对两组相干信号的DOA估计结果的均方根误差对比图。

具体实施方式

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

本发明的步骤为：

步骤一：获取接收阵列的输出信号Y；

设定使用M个全指向性传感器组成接收阵列，且假定空间中有N个远场窄带相干信号，分别以角度θ_n入射到接收阵列上，其中，n＝1,2,…,N，利用接收阵列对入射信号进行接收采样，阵列的输出信号为：

Y＝[y(1),y(2),…y(L)] (1)

其中，y(t)(t＝1,…,L)表示阵列在t时刻的输出信号，L表示快拍数；

步骤二：网格化观测空间，构造超完备阵列流形A；

将观测空间角度在范围[-90°,90°]内以1°的角度间隔均匀划分，得到角度网格点集Θ＝{θ₁,…,θ_K}，其中K为网格点总数，且K＞＞N；根据角度网格点集Θ，构造阵列流形：

A＝[a(θ₁)，a(θ₂)，…，a(θ_k)，…，a(θ_K)] (2)

其中，

表示在网格点θ_k上的导向矢量，简记为a_k，d_m为接收阵列中第m个传感器的位置坐标，m＝1,…,M，λ为入射相干信号的波长，j为虚数单位；

步骤三：结合稀疏表示的思想，将DOA估计问题转化为稀疏信号重构问题，求解如下稀疏矩阵方程：

Y＝AX+n (4)

其中，X表示K×L维的信号矩阵，n表示M×L维的加性高斯白噪声矩阵；由于X中仅有N个非零行向量，X为稀疏矩阵；

步骤四：建立稀疏贝叶斯概率模型；

首先，对阵列输出信号Y的每一列向量进行复高斯分布假设，则Y的似然函数表示为：

其中，Y_·i,X_·i分别表示矩阵Y,X的第i列向量，I_M表示单位矩阵，β>0表示噪声精度(precision)，噪声精度是指噪声方差的倒数，对β进行伽马先验分布假设，即：

p(β)＝Gamma(β|c，d) (6)

其中，c,d为伽马分布的参数；

接下来，对信号矩阵X构造分层稀疏先验；

在第一层先验中，引入由隐变量组成的D×L维矩阵Z、K×D维矩阵W及K维向量μ，将X的每列进行线性变换X_·i＝WZ_·i+μ+ε,(i＝1,…,L)，其中Z_·i表示矩阵Z的第i列向量，ε表示0均值的复高斯向量，则X服从如下复高斯分布：

其中，Λ＝diag(γ)，γ＝[γ₁,…,γ_K]^T，diag(·)表示生成对角矩阵运算；超参数γ含有预设网格点方向上入射信号的功率信息；采用公式(7)所示线性高斯先验模型可将信号X的自相关分量包含于矩阵Λ^-1中，互相关分量包含于矩阵W中；

在第二层先验中，分别对超参数W,Z,μ,γ进行先验假设，假设W,Z的每一列以及μ均服从零均值复高斯分布，γ的每个元素独立同分布，均服从伽马分布，即：

其中，W_·i表示W的第i列向量，γ_i表示γ的第i个元素，α＝[α₁,…,α_D]^T，a,b为伽马分布的参数，δ为精度参数，(·)^-1表示矩阵求逆运算；

在第三层先验中，对超参数α的每一个元素进行伽马先验假设：

其中，g,h为伽马分布的参数；

设置参数初值a₀＝b₀＝c₀＝d₀＝g₀＝h₀＝10^-6，δ₀＝10^-3；记ξ＝{X，W，Z，μ，α，β，γ}，称作隐变量集；

步骤五：采用变分贝叶斯推断方法计算各隐变量ξ_i的近似后验分布q(ξ_i)，得到：

①

即X中第i列向量X_·i的后验分布为复高斯分布，其均值

和方差Σ_X分别为：

其中，Ψ＝Λ^-1＝(diag(γ))^-1，Y_i,Z_i分别表示Y,Z的第i列向量，<·>表示求期望运算，(·)^H表示矩阵的共轭转置运算；

②

即W中第j列向量

的后验分布为复高斯分布，其均值

和方差

分别为：

其中，X_ji表示X的第j行第i列元素，γ_j,μ_j分别表示γ,μ的第j个元素；

③

即Z中第i列向量Z_·i的后验分布为复高斯分布，其均值

和方差Σ_Z分别为：

Σ_Z＝<I_D+W^HΛW>^-1 (18)

④

即μ的后验分布为复高斯分布，其均值μ_μ和方差Σ_μ分别为：

Σ_μ＝<LΛ+δI_K>^-1 (20)

⑤

即α中第i个元素α_i的后验分布为伽马分布,该分布的参数g,h为：

g＝g₀+K (21)

其中，||·||₂表示向量的2范数运算；

⑥

即β的后验分布为伽马分布，该分布的参数c,d分别为：

c＝c₀+LM (23)

其中，tr(·)表示矩阵的迹；

⑦

即γ中第j个元素γ_j的后验分布为伽马分布，该分布的参数a,b分别为：

a＝a₀+L (25)

其中，W_j·表示W的第j列向量，|·|表示绝对值；

根据{X,W,Z,μ,α,β,γ}的更新公式(13)～(26)，在设置隐变量初值后，进行隐变量的迭代更新直至满足收敛条件，则停止迭代；隐变量初值X₍₀₎＝A^H(AA^H)^-1Y,

W₍₀₎＝1_K×M,Z₍₀₎＝1_M×L,μ₍₀₎＝1_K×1,α₍₀₎＝1_M×1,γ₍₀₎＝1_K×1,其中(·)_(r)表示第r步迭代中的变量，

表示n₁×n₂维的全1矩阵，收敛条件为：||γ_(r)-γ_(r-1)||₂≤10^-4；

步骤六：计算入射信号的DOA估计值；

根据步骤五中得到的各隐变量估计值，以γ为纵坐标Θ为横坐标绘制空间方位谱图，在谱图中峰值点所在的网格点角度区域内，通过一维搜索的方式计算得到入射信号的DOA估计值

计算公式如下：

其中，

表示导向矢量a_j对θ_j的偏导值，Re(·)表示取复数的实部，θ_i表示空间方位谱图中第i个峰值点所在的网格点角度区域。

由于超参数γ含有入射信号的功率信息，空间方位谱图中峰值点所对应的网格点即为入射信号DOA估计值；然而，当入射信号没有落在预设网格点上，即网格失配时，根据谱图峰值点所确定的DOA估计值存在较大的误差，需要进行进一步的精确计算。构造目标函数：

其中，

表示

中与γ_j无关的部分，

表示

中与γ_j有关的部分；将

对γ_j求偏导并令结果为0，得到：

其中，

将

对θ_j求偏导并令结果为0，即：

将式(29)代入式(30)中得到：

根据式(31)即可得到DOA估计值的计算公式(27)。

本发明的实施例的步骤如下：

步骤一：得到接收阵列的输出信号Y。

假定有N个远场相干窄带信号以角度

入射到阵元数为M的接收阵列上，设阵列在t(t＝1,2,…,L)时刻的输出信号为y(t)，则阵列的输出信号Y＝[y(1),y(2),…y(L)]，其中，L为快拍数。

步骤二：网格化观测空间，构造超完备阵列流形A。

基于入射信号在空域角度方位内是有限的、稀疏分布的这一特性，将观测空间角度范围[-90°,90°]以某一角度间隔均匀划分成K(K＞＞N)个网格点，得到网格点集Θ＝{θ₁,…,θ_K}，构造超完备阵列流行A＝[a(θ₁),a(θ₂)，…,a(θ_K)]，其中，

表示在网格点θ_k(k＝1,…,K)上的导向矢量，简记为a_k，d_m(m＝1,…,M)为接收阵列中第m个传感器的位置坐标，λ为入射相干信号的波长，j为虚数单位，(·)^T表示矩阵转置运算。

步骤三：结合稀疏表示的思想，将DOA估计问题转化为稀疏信号重构问题，即求解稀疏矩阵方程Y＝AX+n，其中，X表示K×L维的信号矩阵，其为仅有N个非零行向量的稀疏矩阵，n表示M×L维的加性高斯白噪声矩阵。

步骤四：根据稀疏贝叶斯学习理论，建立稀疏贝叶斯概率模型。

指定阵列输出信号Y的每一列向量服从复高斯分布，引入超参数β>0表示高斯白噪声的精度，则Y的似然函数表示为

其中，Y_·i,X_·i分别表示矩阵Y,X的第i列向量，I表示单位矩阵；对β进行伽马先验分布假设，即p(β)＝Gamma(β|c,d)。

构建分层稀疏先验分布，在第一层先验中，引入超参数W,Z,μ,γ，W,Z分别为K×D维和D×L维的矩阵，μ＝[μ₁,…,μ_K]^T，γ＝[γ₁,…,γ_K]^T，指定信号矩阵X的每一列向量服从复高斯分布，

其中，Z_·i表示Z的第i列向量，Λ＝diag(γ)，diag(·)表示生成对角矩阵运算；在第二层先验中，对超参数W,Z,μ,γ进行先验假设，再引入超参数α＝[α₁,…,α_D]^T，指定W的每一行向量均服从零均值、方差为(diag(α))^-1的复高斯分布，指定Z的每一列向量均服从零均值、方差为I_D的复高斯分布，指定μ服从零均值、方差为δ^-1I_K的复高斯分布，指定γ的每一个元素均服从参数为a,b的伽马分布；在第三层先验中，对超参数α的每一个元素指定参数为g,h的伽马分布。

步骤五：根据稀疏贝叶斯学习理论，采用变分贝叶斯推断方法求解隐变量X,W,Z,μ,α,β,γ的后验分布，得到各自统计特性(均值、方差等)的更新公式；设置合适的隐变量初始值，迭代更新隐变量直至收敛，得到隐变量的最优估计值，各隐变量的更新公式如下：

①X中第i(i＝1,…,L)列向量X_·i的均值

和方差Σ_X：

其中，Ψ＝Λ^-1＝(diag(γ))^-1，Y_i,Z_i分别表示Y,Z的第i列向量，<·>表示求期望运算，(·)^H表示矩阵的共轭转置运算，(·)^-1表示矩阵求逆运算；

②W中第j(j＝1,…,D)列向量

的均值

和方差

为：

其中，X_ji表示X的第j行第i列元素，γ_j,μ_j分别表示γ,μ的第j个元素；

③Z中第i(i＝1,…,L)列向量Z_·i的均值

和方差Σ_Z：

Σ_Z＝<I_D+W^HΛW>^-1

④μ的均值μ_μ和方差Σ_μ为：

Σ_μ＝<LΛ+δI_K>^-1

⑤α中第i(i＝1,…,D)个元素α_i的后验分布为伽马分布,分布的参数g,h_i：

g＝g₀+K

其中，||·||₂表示向量的2范数运算；

⑥β的后验分布为伽马分布，分布的参数c,d：

c＝c₀+LM

其中，tr(·)表示矩阵的迹；

⑦γ中第j(j＝1,…,K)个元素γ_j的后验分布为伽马分布，分布的参数a,b：

a＝a₀+L

其中，W_j·表示W的第j行向量，|·|表示绝对值。

步骤6：计算目标信号的DOA估计值。

以γ最优估计值取10倍的以10为底的对数值为纵坐标(单位为分贝dB)，以网格点集Θ为横坐标，绘制空间方位谱图；由于γ是一个稀疏向量，其大多数元素趋于0，仅包含N个非零值，即谱图中的峰值，取第i(i＝1,…,N)个峰值对应网格点所在的角度区域θ_i为搜索区间，计算该搜索区间内存在的目标信号DOA估计值：

其中，

表示导向矢量a_j对θ_j的偏导值，Re(·)表示取复数的实部。

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

1.仿真条件：

假设有2组共4个远场相干窄带信号(每组由2个信号组成，同一组信号相干，不同组信号不相关)入射到接收阵列上。观测空间角度为[-90°,90°]，空间网格化划分间隔为1°。将所提方法记为CASBL，与MFOCUSS、SS-MUSIC、OGSBI、cRVM、L1-SVD和SS-MUSIC这五种方法进行性能比较，以蒙特卡罗曲线为参考，比较各方法在不同条件下的均方根误差(RMSE,average root mean-square error)曲线。

RMSE的计算公式表示为：

其中，Q表示实验次数，Q＝200，

为第q次实验中第n个入射信号的DOA估计值，

为第n个入射信号的真实DOA值。

2.仿真内容与结果：

仿真1：假设这两组相干信号的入射角度分别为[-20.8°,-12.6°]和[12.3°,19.5°]，相关系数分别为[-0.0349+j0.9994,-0.6490+j0.2622]^T和[0.7092+j0.5541,0.7999+j0.0140]^T。使用12元非均匀线阵进行DOA估计，阵列的阵元位置坐标为[-9,-6,-5,-3,-2,-1,1,2,3,5,6,9]λ/2，其中λ为入射窄带信号的波长。快拍数为100，信噪比(SNR,signal-to-noise ratio)为0dB。采用本发明进行DOA估计，得到的空间方位谱图如图1所示，图中横坐标为角度值，纵坐标为归一化幅度谱值(各元素值除以最大元素值)。

从图1可以看出，本发明能够正确分辨这两组相干信号，且谱图中的峰值较为尖锐，表明所提方法的角度分辨能力好；由于这两组信号之间是不相关的，仿真结果说明该方法不仅能够处理相干信号，也能处理独立信号和混合信号(相干+不相关)的DOA估计问题。

仿真2：假设两组相干信号的入射方向分别为

和

随机变量

在[-0.5°,-0.5°]上均匀分布。使用15元均匀线阵进行DOA估计，阵元间隔为入射窄带信号的半波长。采样快拍数为50，信噪比从-5dB增大到20dB，分别采用本发明和其它五种方法进行200次独立的DOA估计实验，计算不同信噪比条件下各方法估计结果的RMSE，得到均方根误差-信噪比曲线如图2所示。图2中横坐标为信噪比，纵坐标为DOA估计结果的RMSE。

从图2可以看出，CASBL在不同的信噪比条件下均具有最小的RMSE，且最贴近蒙特卡罗曲线，即最接近理想估计结果；特别是在低信噪比的情况下，本发明具有最佳的DOA估计性能。

仿真3：在仿真2的基础上，固定信噪比为0dB，将快拍数从30增大到120，分别采用本发明和其它五种方法进行200次独立的DOA估计实验，计算各方法在不同快拍数条件下RMSE，得到均方根误差-快拍数曲线如图3所示，图中横坐标为快拍数，纵坐标为DOA估计结果的RMSE。

从图3中可以看出，本发明在不同的快拍数条件下均具有小的RMSE，即具有最优的估计性能。

Claims (2)

1.一种基于稀疏贝叶斯学习的相干信号DOA估计方法，其特征在于包括下述步骤：

步骤一：获取接收阵列的输出信号Y；

设定使用M个全指向性传感器组成接收阵列，且假定空间中有N个远场窄带相干信号，分别以角度θ_n入射到接收阵列上，其中，n＝1,2,…,N，利用接收阵列对入射信号进行接收采样，阵列的输出信号为：

Y＝[y(1),y(2),…y(L)] (1)

其中，y(t)(t＝1,…,L)表示阵列在t时刻的输出信号，L表示快拍数；

步骤二：网格化观测空间，构造超完备阵列流形A；

将观测空间角度在范围[-90°,90°]内以1°的角度间隔均匀划分，得到角度网格点集Θ＝{θ₁,…,θ_K}，其中K为网格点总数，且K＞＞N；根据角度网格点集Θ，构造阵列流形：

A＝[a(θ₁)，a(θ₂)，…，a(θ_k)，…，a(θ_K)] (2)

其中，

表示在网格点θ_k上的导向矢量，简记为a_k，d_m为接收阵列中第m个传感器的位置坐标，m＝1,…,M，λ为入射相干信号的波长，j为虚数单位；

步骤三：结合稀疏表示的思想，将DOA估计问题转化为稀疏信号重构问题，求解如下稀疏矩阵方程：

Y＝AX+n (4)

其中，X表示K×L维的信号矩阵，n表示M×L维的加性高斯白噪声矩阵；由于X中仅有N个非零行向量，X为稀疏矩阵；

步骤四：建立稀疏贝叶斯概率模型；

首先，对阵列输出信号Y的每一列向量进行复高斯分布假设，则Y的似然函数表示为：

其中，Y_·i,X_·i分别表示矩阵Y,X的第i列向量，I_M表示M×M维的单位矩阵，β>0表示噪声精度(precision)，噪声精度是指噪声方差的倒数，对β进行伽马先验分布假设，即：

p(β)＝Gamma(β|c，d) (6)

其中，c,d为伽马分布的参数；

接下来，对信号矩阵X构造分层稀疏先验；

在第一层先验中，引入由隐变量组成的D×L维矩阵Z、K×D维矩阵W及K维向量μ，将X的每列进行线性变换X_·i＝WZ_·i+μ+ε,(i＝1,…,L)，其中Z_·i表示矩阵Z的第i列向量，ε表示0均值的复高斯向量，则X服从如下复高斯分布：

其中，Λ＝diag(γ)，γ＝[γ₁,…,γ_K]^T，diag(·)表示生成对角矩阵运算；超参数γ含有预设网格点方向上入射信号的功率信息；采用公式(7)所示线性高斯先验模型可将信号X的自相关分量包含于矩阵Λ^-1中，互相关分量包含于矩阵W中；

在第二层先验中，分别对超参数W,Z,μ,γ进行先验假设，假设W,Z的每一列以及μ均服从零均值复高斯分布，γ的每个元素独立同分布，均服从伽马分布，即：

其中，W_·i表示W的第i列向量，γ_i表示γ的第i个元素，α＝[α₁,…,α_D]^T，a,b为伽马分布的参数，δ为精度参数，(·)^-1表示矩阵求逆运算；

在第三层先验中，对超参数α的每一个元素进行伽马先验假设：

其中，g，h为伽马分布的参数；

设置参数初值a₀＝b₀＝c₀＝d₀＝g₀＝h₀＝10^-6，δ₀＝10^-3；记ξ＝{X，W，Z，μ，α，β，γ}，称作隐变量集；

步骤五：采用变分贝叶斯推断方法计算各隐变量ξ_i的近似后验分布q(ξ_i)，得到：

①

即X中第i列向量X_i的后验分布为复高斯分布，其均值

和方差∑_x分别为：

∑_x＝[<β>A^HA+<Ψ>^-1]^-1 (14)

其中，Ψ＝Λ^-1＝(diag(γ))^-1，Y_i，Z_i分别表示Y，Z的第i列向量，<·>表示求期望运算，(·)^H表示矩阵的共轭转置运算；

②

即W中第j列向量

的后验分布为复高斯分布，其均值

和方差

分别为：

其中，X_ji表示X的第j行第i列元素，γ_j，μ_j分别表示γ，μ的第j个元素；

③

即Z中第i列向量Z_i的后验分布为复高斯分布，其均值

和方差∑_z分别为：

∑_z＝<I_D+W^HΛw>^-1 (18)

④

即μ的后验分布为复高斯分布，其均值μ_μ和方差∑_μ分别为：

∑_μ＝<LΛ+δI_K>^-1 (20)

⑤

即α中第i个元素α_i的后验分布为伽马分布,该分布的参数g,h为：

g＝g₀+K (21)

其中，||·||₂表示向量的2范数运算；

⑥q(β)＝Gamma(β|c,d)，即β的后验分布为伽马分布，该分布的参数c,d分别为：

c＝c₀+LM (23)

其中，tr(·)表示矩阵的迹；

⑦

即γ中第j个元素γ_j的后验分布为伽马分布，该分布的参数a,b分别为：

a＝a₀+L (25)

其中，W_j·表示W的第j列向量，|·|表示绝对值；

根据{X,W,Z,μ,α,β,γ}的更新公式(13)～(26)，在设置隐变量初值后，进行隐变量的迭代更新直至满足收敛条件，则停止迭代；隐变量初值X₍₀₎＝A^H(AA^H)^-1Y,

W₍₀₎＝1_K×M,Z₍₀₎＝1_M×L,μ₍₀₎＝1_K×1,α₍₀₎＝1_M×1,γ₍₀₎＝1_K×1,其中(·)_(r)表示第r步迭代中的变量，

表示n₁×n₂维的全1矩阵，收敛条件为：||γ_(r)-γ_(r-1)||₂≤10^-4；

步骤六：计算入射信号的DOA估计值；

根据步骤五中得到的各隐变量估计值，以γ为纵坐标Θ为横坐标绘制空间方位谱图，在谱图中峰值点所在的网格点角度区域内，通过一维搜索的方式计算得到入射信号的DOA估计值

计算公式如下：

其中，

表示导向矢量a_j对θ_j的偏导值，Re(·)表示取复数的实部，θ_i表示空间方位谱图中第i个峰值点所在的网格点角度区域。

2.根据权利要求1所述的一种基于稀疏贝叶斯学习的相干信号DOA估计方法，其特征在于：

由于超参数γ含有入射信号的功率信息，空间方位谱图中峰值点所对应的网格点即为入射信号DOA估计值；然而，当入射信号没有落在预设网格点上，即网格失配时，根据谱图峰值点所确定的DOA估计值存在较大的误差，需要进行进一步的精确计算；构造目标函数：

其中，

表示

中与γ_j无关的部分，

表示

中与γ_j有关的部分；将

对γ_j求偏导并令结果为0，得到：

其中，

将

对θ_j求偏导并令结果为0，即：

将式(29)代入式(30)中得到：

根据式(31)即可得到DOA估计值的计算公式(27)。

Images