CN107703477A

CN107703477A - 基于块稀疏贝叶斯学习的准平稳宽带阵列信号波达方向估计方法

Info

Publication number: CN107703477A
Application number: CN201710810635.1A
Authority: CN
Inventors: 段惠萍; 张新月; 梁瀚明; 马姗姗; 方俊
Original assignee: University of Electronic Science and Technology of China
Current assignee: University of Electronic Science and Technology of China
Priority date: 2017-09-11
Filing date: 2017-09-11
Publication date: 2018-02-16
Anticipated expiration: 2037-09-11
Also published as: CN107703477B

Abstract

本发明公开了一种基于块稀疏贝叶斯学习的准平稳宽带阵列信号波达方向估计方法。本发明利用准平稳宽带信号频谱的帧内相关性和帧间独立性为信号设置相应的块稀疏先验分布模型，通过块稀疏贝叶斯模型对稀疏信号进行估计，从而获得精度更高的估计结果。即首先对阵列接收信号进行适当的分帧处理，对每帧信号分别做傅里叶变换并在频域建立其各自的块稀疏贝叶斯模型，然后在各帧信号独立的假设下，联合所有帧的信息建立一个总的贝叶斯模型，并通过超参数向量来控制所有帧待重构信号的稀疏性，最后利用期望最大化准则获得该超参数向量的迭代更新式。本发明充分利用准平稳宽带阵列信号的短时平稳特性来建立块稀疏模型，能够获得更高的波达方向估计精度。

Description

基于块稀疏贝叶斯学习的准平稳宽带阵列信号波达方向估计方法

技术领域

本发明属于阵列信号处理领域，具体涉及一种基于块稀疏贝叶斯学习的准平稳宽带阵列信号波达方向估计方法。

背景技术

波达方向估计是阵列信号处理的一个重要研究方向，随着宽带信号的广泛应用，如何有效地实现对宽带信号波达方向的高精度、高分辨率估计成为当前的一个研究热点。现实生活中准平稳宽带阵列信号的DOA估计具有广泛的应用，如麦克风阵列对声源的定位、航空站系统对视频信号的定位等。在已有的宽带阵列信号DOA估计方法中，子空间类方法如双边相关变换(Two-side Correlation Transformation，TCT)算法，需要精确地预估信源方位，在预估不准确的情况下算法性能会急剧下降；稀疏重构类方法如分布式压缩感知同步正交匹配追踪(Distributed Compressive Sensing-SimultaneousOrthogonalMatching Pursuit，DCS-SOMP)算法虽然有较低的运算复杂度，但其重构性能较差，另外，将L1-SVD(Singular Value Decomposition)方法扩展到宽带信号处理的系列算法因为采用了奇异值分解，需要已知信源个数，而且在低信噪比的情况下性能不理想。

发明内容

本发明的发明目的在于：针对上述存在的问题，提供一种高估计精度的基于块稀疏贝叶斯学习的准平稳宽带阵列信号波达方向估计方法。

本发明的基于块稀疏贝叶斯学习的准平稳宽带阵列信号波达方向估计方法包括：(1)对阵列接收信号进行分帧处理，对每帧信号分别做离散时间傅里叶变换，利用空域稀疏性构建出稀疏模型；(2)分帧构建块稀疏贝叶斯学习模型；(3)在各帧信号独立的假设下，联合所有帧的信息建立一个总的贝叶斯模型，并将控制每帧待重构信号稀疏性的超参数向量设成一致，利用期望最大化准则得到超参数集合的迭代更新公式；(4)通过迭代求解超参数向量从而获得DOA估计结果。

当利用N个阵元组成的均匀线性阵列对K个远场宽带信号进行DOA估计时，本发明的具体步骤为：

步骤1：将阵列接收信号均匀划分成L帧，并对每帧进行F点离散时间傅里叶变换，得到每帧的阵列接收信号向量y_l＝[y_1,l,y_2,l,...,y_F,l]，其中l＝1,2,...,L，f＝1,2,...,F；

将空间角度平均划分成M个网格，每个网格的角度表示为θ_m，其中m＝1,2,...,M。

基于空域稀疏性构建关于y_f,l的稀疏模型：

y_f,l＝A_f(θ)x_f,l+n_f,l f＝1,2,...,F,l＝1,2,...,L

其中A_f(θ)是方向矩阵，x_f,l是待重构信号源向量，n_f,l是噪声(高斯白噪声)向量的傅里叶变换。

步骤2：分帧构建块稀疏贝叶斯学习模型：

y_l＝Φx_l+n_l

其中x_l表示每帧待重构信号源向量，n_l表示每帧的噪声向量，Φ表示N×F的稀疏贝叶斯学习模型的感知矩阵，每帧的阵列接收信号向量y_l的条件概率密度函数为复高斯分布p(y_l|x_l；λ)～CN(Φx_l,λI_NF)，λ是噪声协方差，每帧待重构信号源向量x_l的先验分布假设为复高斯分布p(x_l；Γ_l,B_l)～CN(0,Σ_l)，其中Γ_l＝diag(γ_1,l,...,γ_M,l)控制着x_l的稀疏性，B_l揭示第l帧信号频谱的帧内相关性，符号表示Kronecker积运算。因为每帧的x_l具有相同的稀疏性，即Γ_l非零元素的位置是一样的，所以将每帧的Γ_l约束成同一矩阵Γ来控制这个相同的稀疏性。即Γ＝diag(γ₁,γ₂,...,γ_M)，对应的有

步骤3：联合每帧的阵列接收信号向量y_l，构建总稀疏贝叶斯学习模型：y＝Ψx+n，其中其中I_L表示L×L的单位矩阵，阵列接收信号向量y＝vec([y₁,y₂,...,y_L])，信号源向量x＝vec([x₁,x₂,...,x_L])，噪声向量n＝vec([n₁,n₂,...,n_L])，其中vec(·)为向量化算符。

因信号源向量x的先验分布为复高斯分布Σ₀＝diag(Σ₁,...,Σ_L)，阵列接收信号向量y的条件概率密度函数为复高斯分布I_NFL是维数为NFL×NFL的单位矩阵。

根据贝叶斯公式，利用条件概率密度函数和先验分布函数可以得到信号源向量的后验分布p(x|y；Θ)～CN(μ_x,Σ_x)，则其后验均值μ_x＝Σ₀Ψ^H(λI_NFL+ΨΣ₀Ψ^H)^-1y，后验协方差矩阵(·)^H表示矩阵的共轭转置操作。

因此，可得超参数集合Θ＝{Γ,B₁,...,B_L,λ}，根据期望最大化准则得到各超参数的迭代更新公式：

符号Tr表示求矩阵的迹运算，J＝M(l-1)+m，采用MATLAB的记法定义：

Φ^m＝Φ(:,(m-1)F+1:mF)；从矩阵Φ的第(m-1)F+1列到第mF列的这个子矩阵

从向量μ_x的第(J-1)F+1个元素到第JF个元素的这个子向量

从矩阵Σ_x的第(J-1)F+1行到第JF行，第(J-1)F+1列到第JF列的这个子矩阵。

基于γ_m、B_l和λ预设的初始值进行迭代更新，当满足迭代停止条件时，停止对γ_m、B_l和λ的迭代更新，并执行步骤4。

其中，优选的初始化方式为：γ_m初始化为1，B_l初始化为单位矩阵。而λ则取决于应用的场景。

步骤4：由M个γ_m中大于0的索引下标m所对应的角度θ_m获得信源波达方向估计结果。

其中迭代停止条件可以是：迭代次数达到指定次数；或则连续两次迭代估计信号的误差小于指定门限，其中表示基于当前迭代更新后的参数γ_m、B_l和λ计算得到的信号源向量的后验均值，表示基于当前迭代的上一次迭代后的参数γ_m、B_l和λ计算得到的信号源向量的后验均值。

综上所述，由于采用了上述技术方案，本发明的有益效果是：利用准平稳宽带信号功率谱的短时平稳特性，建立一个具有分块对角协方差矩阵的块稀疏贝叶斯学习模型，探索信号短时谱的帧内相关性，并联合利用各帧信息来估计揭示信源DOA的超参数向量，从而在不需要知道信源个数的前提下，获得信源波达方向估计结果，且在低信噪比、低快拍数、来波方向较为密集的情况下也有较好的估计精度。

附图说明

图1是本发明的处理流程图；

图2是本发明(提出算法)与其他两种算法(SBL-SML(stochastic maximumlikelihood，SML)算法、忽略相关法)的均方根误差随信噪比变化的比较图；

图3是本发明(提出算法)与忽略信号频谱相关性的算法、SBL-SML算法的均方根误差随帧数变化的比较图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图，对本发明作进一步地详细描述。

参见图1，在将本发明的估计方法用在由N个阵元组成的均匀线性阵列对K个方向为的远场宽带信号进行DOA估计时，具体实现步骤如下：

步骤1、将阵列接收信号划分成L帧，并对每帧进行F点离散时间傅里叶变换，则阵列接收信号可以表示成

其中方向矩阵s_f,l是信号源向量的傅里叶变换，n_f,l是高斯白噪声向量的傅里叶变换，是阵列对方向上的信号在频率点f处的导向矢量，v表示信号传播速率。

为构建稀疏模型，将空间角度平均划分成M个网格，构造新的方向矩阵A_f(θ)＝[a_f(θ₁),...,a_f(θ_M)]，则阵列接收信号可以重新表示成：

y_f,l＝A_f(θ)x_f,l+n_f,l f＝1,2,...,F,l＝1,2,...,L

因为M＞＞K，所以每个待重构向量x_f,l均是稀疏的，且具有相同的支撑集，支撑集代表着信号源方向。

步骤2、分帧构建块稀疏贝叶斯学习模型。

为综合利用每个频谱的信息，探索其相关性，分别给每帧信号构造联合矩阵Y_l＝[y_1,l,y_2,l,…,y_F,l]，X_l＝[x_1,l,x_2,l,...,x_F,l]，N_l＝[n_1,l,n_2,l,...,n_F,l]，则有

vec(Y_l)＝diag(A₁,A₂,...,A_F)vec(X_l),l＝1,...,L

为构造块稀疏结构向量，令vec(X_l)和之间实际上是行变换的关系，所以令D∈R^MF×MF是初等行变换矩阵，且满足：

令y_l＝vec(Y_l)∈C^NF×1，n_l＝vec(N_l)∈C^NF×1，G＝diag(A₁,A₂,...,A_F)，则有

y_l＝Gvec(X_l)+n_l＝GDx_l+n_l＝Φx_l+n_l l＝1,...,L

假设所有噪声向量相互独立，且服从复高斯分布，方差均为λ，则阵列接收信号的条件概率密度函数也服从复高斯分布：

p(y_l|x_l；λ)～CN(Φx_l,λI_NF),l＝1,...,L

I_NF是维数为NF×NF的单位矩阵。

假设不同方向上的信号源是相互独立的，即X_l的行向量是相互独立的，也服从复高斯分布：

超参数γ_m,l≥0控制着X_l的行稀疏性，当γ_m,l＝0时，为零，当γ_m,l＞0时，意味着对应方向θ_m有来波信号。B_m,l是正定矩阵，用来表示信号频谱之间的相关性，为避免过拟合，采取T-SBL(temporal extension of SBL)算法中的措施，将不同行的约束成相同矩阵B_l。

根据的分布以及高斯分布的特性，可以得到向量x_l的先验分布也服从复高斯分布：

其中Γ_l＝diag(γ_1,l,...,γ_M,l)，是协方差矩阵，符号表示Kronecker积运算。

步骤3、联合利用每帧信号，构建总的稀疏贝叶斯模型。

因为每帧待重构向量x_l具有相同的稀疏结构，利用这一先验知识，将所有约束成同一个Γ来控制所有x_l的块稀疏性，在此，假设是相互独立的，构造矩阵X＝[x₁,x₂,...,x_L]，则X的概率分布函数为：

其中Σ₀＝diag(Σ₁,...,Σ_L)是协方差矩阵，估计矩阵X等价于估计其向量化算符x＝vec(X)，

p(x；Γ,B₁,...,B_l)＝p(X；Γ,B₁,...,B_l)～CN(0,Σ₀)

定义Y＝[y₁,y₂,...,y_L]，令y＝vec(Y)，则y的条件概率密度函数是：

其中I_L是维数为L×L的单位矩阵。

根据贝叶斯公式，利用条件概率密度函数和先验分布函数可以得到信号源向量的后验分布：

p(x|y；Θ)～CN(μ_x,Σ_x)

其中，超参数集合Θ＝{Γ,B₁,...,B_L,λ}，后验均值μ_x＝Σ₀Ψ^H(λI_NFL+ΨΣ₀Ψ^H)^- ¹y，后验协方差矩阵

要得到超参数集合Θ的估计值，利用EM(期望最大化准则)方法，将x当作隐藏变量，最大化代价函数：

Θ^(old)表示前一次迭代中超参数的估计值，E表示对x的后验分布logp(x|y；Θ^old)求期望。

在估计Γ和B_l时，只需要最大化代价函数中含有Γ和B_l的项，因此Q(Θ)可以简化为

然后对Q(Γ,B₁,...,B_L)分别求γ_m(m＝1,...,M)的偏导数，得

令上式为0，即可获得γ_m的迭代更新公式为：

符号Tr表示求矩阵的迹运算，其中J＝M(l-1)+m，采用MATLAB的记法定义：

对Q(Γ,B₁,...,B_L)分别求B_l(l＝1,...,L)的偏导数

同理可得到B_l的迭代更新公式：

估计噪声方差λ时，只需要最大化代价函数中含有λ的项，因此Q(Θ)可以简化为

利用T-SBL算法的简化步骤，可以得到噪声方差λ的估计式：

采用MATLAB的记法定义：Φ^m＝Φ(:,(m-1)F+1:mF)。

步骤4、利用上述各超参数的迭代更新公式进行迭代，直至连续两次迭代估计信号的误差小于指定门限或者迭代次数达到指定次数，根据迭代停止时的输出结果Γ即可获得信源波达方向估计结果。

本发明的效果通过仿真图2和图3进行说明：

仿真实验中，信号源为语音信号，取其700Hz-1200Hz的频谱，四个信号源的真实方向是[-10°,-5°,3°,10°]，均匀线性阵列的阵元个数是28，阵元间距是半波长，分帧处理之后，每帧信号取了9个频点。蒙特卡洛实验次数为W＝200，用均方根误差(RMSE)来比较每个算法的性能，其定义为是第w次实验估计的第k个信源的方向。

图2是帧数取5时，本发明与忽略信号频谱相关性的算法、SBL-SML算法的均方根误差随信噪比SNR变化的比较图，SNR取值范围是[4:12]dB，间隔为2dB。

图3是信噪比取7dB时，本发明与忽略信号频谱相关性的算法、SBL-SML算法的均方根误差随帧数变化的比较图，帧数的取值范围是[2:14]，间隔为3。

从实验结果可以看出，针对准平稳信号，通过分帧考虑信号频谱之间的相关性，本发明获得了更高的方向估计精度。

以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合。

Claims

1.基于块稀疏贝叶斯学习的准平稳宽带阵列信号波达方向估计方法，其特征在于，包括下列步骤：

将空间角度平均划分成M个网格，每个网格的角度表示为θ_m，其中m＝1,2,...,M；

步骤2：分帧构建块稀疏贝叶斯学习模型：y_l＝Φx_l+n_l，其中x_l表示每帧待重构信号源向量，n_l表示每帧的噪声向量，Φ表示N×F的稀疏贝叶斯学习模型的感知矩阵，N为阵列数；定义每帧信号频谱的帧内相关性矩阵为B_l，控制各x_l的稀疏性的超参数向量均为Γ＝diag(γ₁,γ₂,...,γ_M)；

步骤3：联合每帧的阵列接收信号向量y_l，构建总稀疏贝叶斯学习模型：y＝Ψx+n，其中I_L表示L×L的单位矩阵，阵列接收信号向量y＝vec([y₁,y₂,...,y_L])，信号源向量x＝vec([x₁,x₂,...,x_L])，噪声向量n＝vec([n₁,n₂,...,n_L])，其中vec(·)为向量化算符；

并基于迭代更新公式对γ_m、B_l和噪声方差λ进行迭代更新，其中γ_m、B_l和λ的初始值为预设值：

<mrow> <msub> <mi>&gamma;</mi> <mi>m</mi> </msub> <mo>&LeftArrow;</mo> <mfrac> <mrow> <munderover> <mo>&Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <mi>T</mi> <mi>r</mi> <mo>&lsqb;</mo> <msubsup> <mi>B</mi> <mi>l</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mrow> <mo>(</mo> <msubsup> <mo>&Sigma;</mo> <mi>x</mi> <mi>J</mi> </msubsup> <mo>+</mo> <msubsup> <mi>&mu;</mi> <mi>x</mi> <mi>J</mi> </msubsup> <msup> <mrow> <mo>(</mo> <msubsup> <mi>&mu;</mi> <mi>x</mi> <mi>J</mi> </msubsup> <mo>)</mo> </mrow> <mi>H</mi> </msup> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> <mrow> <mi>F</mi> <mo>&times;</mo> <mi>L</mi> </mrow> </mfrac> <mo>,</mo> <mi>m</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>M</mi> <mo>,</mo> </mrow>

<mrow> <msub> <mi>B</mi> <mi>l</mi> </msub> <mo>&LeftArrow;</mo> <mfrac> <mn>1</mn> <mi>M</mi> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <mfrac> <mrow> <msubsup> <mo>&Sigma;</mo> <mi>x</mi> <mi>J</mi> </msubsup> <mo>+</mo> <msubsup> <mi>&mu;</mi> <mi>x</mi> <mi>J</mi> </msubsup> <msup> <mrow> <mo>(</mo> <msubsup> <mi>&mu;</mi> <mi>x</mi> <mi>J</mi> </msubsup> <mo>)</mo> </mrow> <mi>H</mi> </msup> </mrow> <msub> <mi>&gamma;</mi> <mi>m</mi> </msub> </mfrac> <mo>,</mo> <mi>l</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>L</mi> <mo>,</mo> </mrow>

<mrow> <mi>&lambda;</mi> <mo>&LeftArrow;</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <mi>y</mi> <mo>-</mo> <msub> <mi>&Psi;&mu;</mi> <mi>x</mi> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>+</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <mi>T</mi> <mi>r</mi> <mrow> <mo>(</mo> <msubsup> <mo>&Sigma;</mo> <mi>x</mi> <mi>J</mi> </msubsup> <msup> <mrow> <mo>(</mo> <msup> <mi>&Phi;</mi> <mi>m</mi> </msup> <mo>)</mo> </mrow> <mi>H</mi> </msup> <msup> <mi>&Phi;</mi> <mi>m</mi> </msup> <mo>)</mo> </mrow> </mrow> <mrow> <mi>F</mi> <mo>&times;</mo> <mi>N</mi> <mo>&times;</mo> <mi>L</mi> </mrow> </mfrac> <mo>,</mo> </mrow>

其中，符号Tr(·)表示求矩阵的迹运算，(·)^H表示矩阵的共轭转置操作，信号源向量的后验均值μ_x＝Σ₀Ψ^H(λI_NFL+ΨΣ₀Ψ^H)^-1y，Σ₀＝diag(Σ₁,...,Σ_L)，I_NFL表示NFL×NFL的单位矩阵，J＝M(l-1)+m；信号源向量的后验协方差矩阵Φ^m＝Φ(:,(m-1)F+1:mF)；

当满足迭代停止条件时，停止对γ_m、B_l和λ的迭代更新，并执行步骤4；

2.如权利要求1所述的方法，其特征在于，迭代停止条件为迭代次数达到指定次数。

3.如权利要求1所述的方法，其特征在于，迭代停止条件为连续两次迭代估计信号的误差小于指定门限，其中表示基于当前迭代更新后的参数γ_m、B_l和λ计算得到的信号源向量的后验均值，表示基于当前迭代的上一次迭代后的参数γ_m、B_l和λ计算得到的信号源向量的后验均值。

4.如权利要求1所述的方法，其特征在于，各γ_m的初始值为1。

5.如权利要求1所述的方法，其特征在于，各B_l的初始值单位矩阵。