作为阵列信号处理技术的重要研究方向之一，信号到达角(DOA)估计被广泛应用到雷达、声呐等领域。以MUSIC[和ESPRIT[为代表的传统空间谱估计算法，实现简单，且空间分辨率高，但要求高信噪比以及多快拍数。随着压缩感知理论的提出与发展，将压缩感知理论应用到阵列信号处理中成为DOA估计算法新的研究方向，学者们提出了大量性能优越的DOA估计算法。其中最为经典的是由D.Malioutov等人提出的

${l_1}{\rm{ - SVD}}$算法，其首先通过空间网格划分的方式，构造了基于阵列流型的过完备基矩阵，然后利用奇异值(singular value decomposition, SVD)分解对阵列接收信号矩阵降维，最后利用二阶锥规划求解

${l_1}$范数优化问题，并估计出信号的DOA。随后，在

${l_1}{\rm{ - SVD}}$算法的基础上，学者们又先后提出了许多基于稀疏表示的DOA估计算法，用以解决宽带[、二维角度[、相干信号源[等DOA估计问题。基于稀疏表示的DOA估计算法多是通过求解

${l_1}$范数优化问题估计出信号的DOA，但是在求解

${l_1}$范数优化问题的过程中，正则化参数多是通过人工设置的方式选取。鉴于稀疏问题的求解完全可以放在贝叶斯估计框架中分析和表示，而且这种统计优化方法更容易被理解和接受。因此近几年一些学者也在不断开展基于稀疏贝叶斯学习的DOA估计研究工作，主要研究二维DOA估计[、空间网格划分失配[以及算法收敛速度

针对此问题，首先假设空间中存在一个均匀线阵用于接收空间窄带目标信号，对空间进行网格划分使得空间信号稀疏化，从而完成信号稀疏表示，并得出可通过求解

${l_1}$范数优化问题估计信号DOA的结论，为应用稀疏贝叶斯理论估计信号的DOA，将阵列接收信号由复数数据转为实数数据。然后根据稀疏贝叶斯理论可知，阵列接收信号服从高斯分布。为估计出未知参数，对稀疏信号指定高斯分布，并为其方差以及噪声方差的倒数指定Gamma分布，得出可通过求解待估计参数的最大后验概率分布估计出未知参数。最后为了避免后验概率分布的直接求解，降低算法复杂度，通过变分贝叶斯学习算法寻找后验概率分布的近似分布，并通过监控下界值判断算法是否收敛，估计出未知参数，从而获得信号的DOA。

1 基于稀疏表示的DOA估计模型

假设

$K$个远场窄带信号入射到

$M$个各向同性的均匀线阵上，阵元间距为

$d$，各信号来波方向为

${\theta _{\rm{i}}}(i = 1,2, \cdots ,K)$，则阵列接收信号为

${{y}}(t) = {{As}}(t) + {{n}}(t)$

式中：${{y}}\left( t \right) = {\left[ {{y_1}(t),{y_2}(t), \cdots {y_M}(t)} \right]^{\rm{T}}}$为

$M \times 1$维的阵列接收信号矢量；${{s}}(t) = {\left[ {{s_1}(t),{s_2}(t), \cdots ,{s_N}(t)} \right]^{\rm{T}}}$为

$K \times 1$维的空间信号矢量；${{A}} =\left[ {{{a}}\left( {{\theta _1}} \right),{{a}}\left( {{\theta _2}} \right), \cdots ,{{a}}\left( {{\theta _K}} \right)} \right]$为

$M \times K$维阵列流形矩阵，${ a}\left( {{\theta _k}} \right) = {\rm{[}}1,\exp ( - {\rm{j}}2{\text π} f{\tau _k}), \cdots ,\exp $

$( - {\rm{j}}2{\text π} f(M - 1){\tau _k})]$T为

$M \times 1$维的方向矢量，${\tau _k} = d\sin {\theta _k}/c$，${\theta _k}$为第

$k$个信号的来波方向；${ n}(t) = [{n_1}(t),{n_2}(t), \cdots ,$

${n_M}(t){{\rm{]}}^{\rm{T}}}$为

$M \times 1$维的噪声矢量。

空间信号是可稀疏的，通过特定的网格划分可将空间信号稀疏化。将空间均匀划分为

$N$份，$\left\{ {{\theta _1},{\theta _2}, \cdots ,{\theta _N}} \right\}$，$N \gg K$。假设每一个

${\theta _n}(n = 1,2, \cdots ,N)$都对应一个空间信号

${x_n}(n = 1,2, \cdots ,N)$，这样便构造出了稀疏度为

$K$的

$N \times 1$维稀疏空间信号矢量：

${{x}}(t) = {[{x_1}(t),{x_2}(t), \cdots ,{x_N}(t)]^{\rm{T}}}$

${{x}}(t)$中只有

$K$个位置的元素是非零的，对应着空间中实际存在的信号

$s(t)$，其余

$N - K$个位置的元素为零。则稀疏化后的空间信号矢量对应的

$M \times N$维阵列流型矩阵

${{\varPhi}} $为

${{\varPhi}} {\rm{ = }}[{{a}}({\theta _1}),{{a}}({\theta _2}), \cdots ,{{a}}({\theta _N})]$

则稀疏表示模型下的DOA估计数学模型为

${{y}}(t) = {{\varPhi}} {{x}}(t) + {{n}}(t)$

(1)

通过求解如下

${l_1}$范数优化问题得出重构信号：

$\left\{ {\begin{array}{*{20}{l}}{\hat { x} = \arg \min {{\left\| { x} \right\|}_1}}\\{{\rm s.t.}{{\left\| {{ y} - { {Ax}}} \right\|}_2}\sigma }\end{array}} \right.$

式中

$\sigma $为噪声标准差。

稀疏贝叶斯理论适用于对实数数据的处理，而阵列接收数据是复数，为了将该理论应用到DOA估计中，需要将观测数据实数化以构建实数域的优化模型，因此式(1)需改写为

$\left( {\begin{array}{*{20}{c}} {{\rm Re}({{y}})} \\ {\rm{Im}} ({{y}})\end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{\rm Re}({{\varPhi}} )}&{ - \rm{Im} ({{\varPhi}} )} \\ {\rm{Im} ({{\varPhi}} )}&{R({{\varPhi}} )} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{\rm Re}({{x}})} \\ {\rm{Im}} ({{x}}) \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {{\rm Re}({{n}})} \\ {\rm{Im}} ({{n}}) \end{array}} \right)$

式中

${\rm Re}$和

$\operatorname{Im} $分别表示数据的实部与虚部。

2 基于稀疏变分贝叶斯学习的DOA估计算法

2.1 稀疏贝叶斯模型

对于压缩感知下的DOA估计模型：

${{y}} = {{\varPhi}} {{x}} + {{n}}$

式中：${{y}}$为观测数据，${{\varPhi}} $为观测矩阵，${{x}}$为待求稀疏系数矢量，${{n}}$为均值为0、方差为

${\sigma ^2}$的高斯噪声。

从而可以得到

${{y}}$服从高斯分布，即

$\begin{gathered} p({{y}}\left| {{x}} \right.,{\sigma ^2}) = \prod\limits_{n = 1}^N {p({{{y}}_n}\left| {{{{\varPhi}} _n}} \right.{{x}},{\sigma ^2})} {\rm{ = }} \hfill \\\quad\quad\quad\quad\quad{\rm{ (2{\text π}}}{\sigma ^2}{{\rm{)}}^{ - \frac{N}{2}}}\exp {( - \frac{1}{{2{\sigma ^2}}}\left\| {{{y}} - {{\varPhi}} {{x}}} \right\|)^2} \hfill \\ \end{gathered} $

根据贝叶斯估计理论，需要对参数

${{x}}$和

${\sigma ^2}$进行最大后验估计。为避免过匹配问题，需要对参数指定先验分布。根据相关向量机理论，为待估计参数

${{x}}$中的每一个元素指定均值为0、方差为

${\alpha _i}^{ - 1}$的高斯先验，即

$p({{x}}\left| {{\alpha}} \right.) = \prod\limits_{i = 1}^N {N({x_i}\left| 0 \right.,{\alpha _i}^{ - 1})} $

式中

${{\alpha}} = ({\alpha _1},{\alpha _2}, \cdot \cdot \cdot ,{\alpha _N})$为超参数，控制着待估计参数

${x_i}$的估计精度。由于高斯分布的方差的倒数的共轭分布为Gamma分布，因此为超参数

${\alpha _i}$以及噪声参数

${\alpha _0}$指定Gamma分布，即

$p({{\alpha}} ) = \prod\limits_{i = 1}^N {\varGamma ({\alpha _i}\left| a \right.,b)} $

$p({\alpha _0}) = \varGamma ({\alpha _0}\left| c \right.,d)$

式中

${\alpha _0} = {1 / {{\sigma ^2}}}$。Gamma分布为

$\varGamma (\xi \left| {a,b} \right.) = \frac{{{b^a}}}{{\varGamma (a)}}{\xi ^{a - 1}}\exp ( - b\xi )$

根据贝叶斯理论，可以得到待估计参数的最大后验概率分布为

$p({{x}},{{\alpha}} ,{\alpha _0}|{{y}}) = \frac{{p({{y}}|{{x}},{{\alpha}} ,{\alpha _0})p({{x}},{{\alpha}} ,{\alpha _0})}}{{p({{y}})}}$

$p({{y}}) = \iiint {p({{y}}|{{x}},{{\alpha}} ,{\alpha _0})}p({{x}},{{\alpha}} ,{\alpha _0}){\rm{d}}{{x}}{\rm{d}}{{\alpha}} {\rm{d}}{\alpha _0}$

$p(y)$的计算通常要通过高维、复杂的积分，是不易求解的。稀疏贝叶斯模型下通过分解方法求待估计参数最大后验概率分布，即

$p({{x}},{{\alpha}} ,{\alpha _0}|{{y}}) = p({{x}}|{{y}},{{\alpha}} ,{\alpha _0})p({{\alpha}} ,{\alpha _0}|{{y}})$

然后再利用积分求解

$p({{x}}|{{y}},{{\alpha}} ,{\alpha _0})$，算法实现复杂。为了简化

$p({{y}})$的求解，将变分贝叶斯理论引入稀疏贝叶斯估计中。

2.2 变分稀疏贝叶斯估计

设由未知待估计参数

${{x}}$，$\alpha $和

${{{\alpha}} _0}$组成的参数集

${{\vartheta}} $为

${{\vartheta}} = \left\{ {{{{\vartheta}} _1},{{{\vartheta}} _2},{\vartheta _3}} \right\} = \left\{ {{{x}},{{\alpha}} ,{\alpha _0}} \right\}$

观测量

${ y}$的边缘似然函数为

$p({{y}}){\rm{ = }}\frac{{p({{y}},{{\vartheta}} )}}{{p({{x}},{{\vartheta}} |{{y}})}}$

(2)

依据变分理论，在式(2)中引入一个关于参数集

${{\vartheta}} $的

$q$分布，即

$p({{y}}){\rm{ = }}{{\frac{{p({{y}},{{\vartheta}} )}}{{q({{\vartheta}} )}}} / {\frac{{p({{{{\vartheta}}}} |{{y}})}}{{q({{\vartheta}} )}}}}$

(3)

对式(3)两边取对数得

$\ln (p({{y}})) = \ln \frac{{p({{y}},{{\vartheta}} )}}{{q({{\vartheta}} )}} - \ln \frac{{p({{\vartheta}} |{{y}})}}{{q({{\vartheta}} )}}$

(4)

设

$\iiint {q({{\vartheta}} )d{{{\vartheta}} _1}d{{{\vartheta}} _2}d{\vartheta _3}} = 1$，式(4)可转换为

$\begin{gathered} \ln (p({{y}})) = \iiint {q({{\vartheta}} )\ln \frac{{p(y,{{\vartheta}} )}}{{q({{\vartheta}} )}}}{\rm{d}}{{{\vartheta}} _1}{\rm{d}}{{{\vartheta}} _2}{\rm{d}}{\vartheta _3} - \hfill \\ \quad\quad\quad\quad{\rm{ }}\iiint {q({{\vartheta}} )\ln \frac{{p({{\vartheta}} |{{y}})}}{{q({{\vartheta}} )}}}{\rm{d}}{{{\vartheta}} _1}{\rm{d}}{{{\vartheta}} _2}{\rm{d}}{\vartheta _3} \hfill \\ \end{gathered} $

(5)

式(5)简记为

$\ln (p({{y}})){\rm{ = }}L(q({{\vartheta}} )) + {\rm{KL}}(q({{\vartheta}} )||p({{\vartheta}} |{{y}}))$

$L(q({{\vartheta}} )){\rm{ = }}\iiint {q({{\vartheta}} )\ln \frac{{p({{y}},{{\vartheta}} )}}{{q({{\vartheta}} )}}}{\rm d}{{{\vartheta}} _1}{\rm d}{{{\vartheta}} _2}{\rm d}{\vartheta _3}$

${\rm{KL}}(q({{\vartheta}} )||p({{\vartheta}} |{{y}})){\rm{ = }} - \iiint {q({{\vartheta}} )\ln \frac{{p({{\vartheta}} |{{y}})}}{{q({{\vartheta}} )}}}{\rm{d}}{{{\vartheta}} _1}{\rm{d}}{{{\vartheta}} _2}{\rm{d}}{\vartheta _3}$

式中

${\rm{KL}}(q({{\vartheta}} )||p({{\vartheta}} |{{y}}))$称为KL散度，表示概率分布

$q({{\vartheta}} )$与后验分布

$p({{\vartheta}} |{{y}})$的近似程度，KL散度越小，近似程度越高。变分贝叶斯理论通过最小化KL散度寻找与后验分布

$p({{\vartheta}} |{{y}})$最近似的概率分布

$q({{\vartheta}} )$。$\ln (p({{y}}))$只与观测数据

$y$有关，与待估计参数

$x$、$\alpha $和

${\alpha _0}$无关；概率分布

$q({{\vartheta}} )$与后验分布

$p({{\vartheta}} |{{y}})$相等时，$L(q({{\vartheta}} ))$值最大。由于需要寻求与后验分布

$p({{\vartheta}} |{{y}})$最近似的概率分布

$q({{\vartheta}} )$代替后验分布

$p({{\vartheta}} |{{y}})$，而KL散度与后验分布

$p({{\vartheta}} |{{y}})$有关，所以通过最小化KL散度寻求概率分布

$q({{\vartheta}} )$是行不通的。而概率分布

$p({{y}},{{\vartheta}} )$是易求的，所以通过最大化

$L(q({{\vartheta}} ))$寻求关于变量

${{x}}$、${{\alpha}} $和

${\alpha _0}$的

$q({{\vartheta}} )$分布。

根据均值域理论，关于

${{x}}$、${{\alpha}} $和

${\alpha _0}$的联合分布为

$q({{\vartheta}} ) = q({{x}},{{\alpha}} ,{\alpha _0}){\rm{ = }}q({{x}})q({{\alpha}} )q({\alpha _0})$

各参数服从以下分布：

$q({{x}}){\rm{ = }}N({{x}}|\overline {{u}} ,\overline {{\varSigma}} )$

$q({{\alpha}} ){\rm{ = }}\sum\limits_{m = 1}^N {\varGamma ({\alpha _m}|{{\overline a }_m},{{\overline b }_m})} $

$q({\alpha _0}){\rm{ = }}\varGamma ({\alpha _0}|\overline c ,\overline d )$

式中：

$\overline {{u}} = {\alpha _0}\overline {{\varSigma}} {{{\varPhi}} ^{\rm{T}}}{{y}}$

(6)

$\overline {{\varSigma}} {\rm{ = (diag}}({\alpha _m}) + {\alpha _0}{{{\varPhi}} ^{\rm{T}}}{{\varPhi}} {)^{ - 1}}$

(7)

${\overline a_m} = a + 1/2$

(8)

${\overline b_m} = b + \frac{{|\overline {{u}} {|^2} + {{\bar {{\varSigma}} }_{mm}}}}{2}$

(9)

$\overline c = c + (N + 1)/2$

(10)

$\overline d = d + \frac{{||{{y}} - {{\varPhi}} \overline {{\mu}} |{|^2} + {\rm{tr}}(\overline {{\varSigma}} {{{\varPhi}} ^{\rm{T}}}{{\varPhi}} )}}{2}$

(11)

${\alpha _m} = {\overline a_m}/{\overline b_m}$

${\alpha _0} = \overline c/\overline d$

式中：${\alpha _{_m}}$表示

$q({{\alpha}} )$的均值，${\alpha _0}$表示

$q({\alpha _0})$的均值，${\overline {{\mu}} _m}$表示向量

$\overline {{u}} $的第

$m$个元素，${\overline {{\varSigma}} _{mm}}$表示矩阵

$\overline {{\varSigma}} $主对角线的第

$m$个元素，$a = b = c = d = {10^{ - 6}}$[。

对于变分贝叶斯学习方法，是通过监控下界[$L(q({{x}},{{\alpha}} ,{\alpha _0}))$来控制算法收敛：

$\begin{gathered} L(q({{x}},{{\alpha}} ,{\alpha _0})) = \left\langle {\ln p({{y}}|{{x}},{\alpha _0})} \right\rangle + \left\langle {\ln p({{x}}|{{\alpha}} )} \right\rangle + \hfill \\ \quad\quad\quad\quad\quad\quad\;\,\left\langle {\ln p({{\alpha}} )} \right\rangle + \left\langle {\ln p({\alpha _0})} \right\rangle - \left\langle {\ln q({{x}})} \right\rangle - \hfill \\ \quad\quad\quad\quad\quad\quad\;\,\left\langle {\ln q({{\alpha}} )} \right\rangle - \left\langle {\ln q({\alpha _0})} \right\rangle \hfill \\ \end{gathered} $

在计算

$L(q({{x}},{{\alpha}} ,{\alpha _0}))$时，每一项的计算量都很大，因此对

$L(q({{x}},{{\alpha}} ,{\alpha _0}))$进行简化，简化后的

$ L(q$

$({{x}},{{\alpha}} ,{\alpha _0}))$表达式如下：

$\begin{gathered} L(q({{x}},{{\alpha}} ,{\alpha _0})){\rm{ = }}\frac{1}{2}\ln |\overline {{\varSigma}} | - {\overline a _m}\sum\limits_{m = 1}^M {\ln {{\overline b }_m}} + \hfill \\ \quad\quad\quad\quad\quad\;\;\;\;{\rm{ }}\frac{1}{2}\sum\limits_{m = 0}^{N - 1} {\ln {{\overline {{\varSigma}} }_{mm}}} - \overline c \ln \overline d + {L_{{\rm{const}}}} \hfill \\ \end{gathered} $

(12)

$\begin{gathered} {L_{{\rm{const}}}} = - \frac{N}{2}\ln 2{\rm{{\text π}}} + \frac{{M + N}}{2} - \hfill \\ \quad\quad\;\;\;{\rm{ }}(M + N + 1)\ln \varGamma (a) + (M + N + 1)a\ln b + \hfill \\ \quad\quad\;\;\;{\rm{ }}(M + N)\ln \varGamma ({{\overline a}_m}) + \ln \varGamma (\overline c) \hfill \\ \end{gathered} $

式中

${L_{{\rm{const}}}}$是与变量

${{x}}$、${{\alpha}} $和

${\alpha _0}$无关的常数部分。

基于变分稀疏贝叶斯的DOA估计算法流程如下：

输入：观测向量

${{y}}$、测量矩阵

${{\varPhi}} $和收敛条件

$\sigma $。

输出：重构信号

${{x}}$。

1)初始化

$\overline {{u}} $、$\overline {{\varSigma}} $、${\overline a_m}$、${\overline b_m}$、$\overline c$、$\overline d$；

2)利用式(6)和(7)更新

$\overline {{u}} $、$\overline {{\varSigma}} $；

3)利用式(8)~(11)更新

${\overline a_m}$、${\overline b_m}$、$\overline c$、$\overline d$；

4)利用式(12)更新

${L_t}(q({{x}},{{\alpha}} ,{\alpha _0})$，$t$为当前更新次数；

5)如果：$\displaystyle\frac{{{L_t}(q({{x}},{{\alpha}} ,{\alpha _0})) - {L_{t - 1}}(q({{x}},{{\alpha}} ,{\alpha _0}))}}{{{L_{t - 1}}(q({{x}},{{\alpha}} ,{\alpha _0}))}} < \sigma $，则停止迭代，输出重构信号

${{x}}$，否则重返步骤2)，进行下一次迭代。

3 仿真结果

为评估基于稀疏变分贝叶斯学习的DOA估计算法的性能，本文通过MATLAB仿真比较在单快拍条件下基于VSBL的DOA估计算法和基于SBL的DOA估计算法分别在不同信噪比、不同阵元数下的DOA估计误差及其成功率以及算法运行时间。

3.1 不同信噪比下的DOA估计结果比较

假设空间中存在由20个阵元组成的均匀线阵，同时有2个信号分别以入射角10°和20°入射到该阵列中，信号的信噪比以2 dB为步进，从−10 dB到10 dB变化。用2种算法比较不同信噪比下的DOA估计结果，如果所估计出角度在误差允许范围内，则认为本次实验是成功的，每个条件下的实验进行500次，并统计各条件下的DOA估计误差及其成功率，如

图 1

图 1 不同信噪比下的DOA估计误差

图 2

图 2 不同信噪比下的DOA估计成功率

3.2 不同阵元数下的DOA估计结果比较

假设空间中存在一个均匀线阵，且同时有2个信号分别以入射角10°和20°入射到该阵列中，阵元数以2为步进从6阵元到20阵元变化，信噪比为15 dB，用2种算法比较不同阵元数下的DOA估计结果。如果所估计出角度在误差允许范围内，则认为本次实验是成功的，每个条件下的实验进行500次，并统计各条件下的DOA估计误差及其成功率，如

图 3

图 3 不同阵元数下的DOA估计误差

图 4

图 4 不同阵元数下的DOA估计成功率

3.3 算法运行时间比较

仿真条件：阵元数为20，信噪比比为15 dB，蒙特卡罗仿真次数为500次。仿真所用CPU为：Intel(R) Core(TM) i5-4570，运行内存为4 GB，MATLAB版本为R2014。

基于BCS的DOA估计算法运行时间为212 s。而基于VBCS的DOA估计算法的运行时间为92 s。

由算法运行时间结果可知，基于VBCS的DOA估计算法的运行时间要小于基于BCS的DOA估计算法的运行时间。因此基于VBCS的DOA估计算法的收敛速度要快于基于BCS的DOA估计算法收敛速度。

4 结论

本文提出了一种基于VSBL的DOA估计算法。首先建立了基于稀疏表示的DOA估计模型，并在此基础上建立了稀疏贝叶斯模型。然后通过变分贝叶斯学习算法的引入，简化了稀疏贝叶斯模型中最大后验概率的求解，并总结出该DOA估计算法的求解步骤。最后通过仿真比较了基于VSBL的DOA估计算法和基于SBL的DOA估计算法的性能，并得出以下结论：

1)本文算法在低信噪比下具有更高的DOA估计精度以及成功率，更有利于在复杂电磁环境下DOA估计的应用；

2)本文算法可通过更少的振元数高精度、高成功率估计出信号到达角，从而减轻硬件对大量数据存储、传输和处理的压力；

3)本文算法具有更低的算法复杂度，大幅度减少了算法的运行时间，更有利于DOA估计实时性要求的实现；

4)本文算法只适用于单快拍下的DOA估计问题，将该思想引用到多快拍下的DOA估计问题是下一步的研究方向。