CN104977558A

CN104977558A - 一种基于贝叶斯压缩感知的分布源中心波达方向估计方法

Info

Publication number: CN104977558A
Application number: CN201510336392.3A
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: 2015-06-16
Filing date: 2015-06-16
Publication date: 2015-10-14
Anticipated expiration: 2035-06-16
Also published as: CN104977558B

Abstract

本发明提供一种基于贝叶斯压缩感知的分布源中心波达方向估计方法，属于无线移动通信技术领域。主要解决当信源中心波达角不在角度采样栅格上时，中心波达方向估计的固有误差问题。本发明通过设置一个平行的均匀线阵组成的天线阵列，建立分布源的近似阵列数据接收模型，并对空域角度进行采样，利用阵列导向矢量构造参数化的过完备冗余字典，使得分布源中心波达方向估计问题被转化成为稀疏矩阵方程求解问题，进而采用贝叶斯压缩感知方法对该方程组进行求解，并获得未知稀疏向量的最稀疏解，根据稀疏解与空域角度一一对应的关系获得中心波达方向的估计值。本发明方法计算复杂度低，且在小快拍数情况下具有分辨率和精确度高等特点。

Description

一种基于贝叶斯压缩感知的分布源中心波达方向估计方法

技术领域

本发明属于无线通信技术领域，具体涉及一种基于贝叶斯压缩感知的分布源波达方向估计方法。

背景技术

随着无线通信系统对通信能力需求的不断增加，利用无线收发器的自适应波束形成技术能够有效提升通信系统的效率。其中，无线信号源的波达方向(direction of arrival，简称DOA)是波束形成技术不可或缺的参数。因此，DOA估计是现代无线通信系统的一个基础并且非常重要的研究问题。另外，DOA参数的估计精度对无线通信系统性能有着至关重要的影响。

由于无线移动通信的多径传播现象，传统的DOA估计方法都是基于点源模型，该模型假设各多径分量在无线信道中是可分离的。但是，实际中信号的多径分量通常以簇的形式出现，不能被简单地看成一个点源。这时，采用空间分布源模型更能反应信号空间分布特性。因而，分布源参数估计方法具有重要的实际意义和广阔的应用前景。

目前，分布源中心DOA的估计技术大多是子空间类方法，如DISPARE方法[Y.Meng,P.Stoica,and K.M.Wong,“Estimation of the directions of arrival of spatially dispersed signals inarray processing,”IEE Proc.-Radar,Sonar Navig.,vol.143,no.1,pp.1-9,1996.]和TLS-ESPRIT方法[S.Shahbazpanahi,S.Valaee,and M.H.Bastani,“Distributed source localization using ESPRITalgorithm,”IEEE Trans.Signal Process.,vol.49,no.10,pp.2169-2178,2001.]。然而，这些方法都需要大量的快拍数和较高的信噪比才能获得精确的DOA估计。值得注意的是，最近几年，根据信号源的空域稀疏性，压缩感知方法被研究者们引入到点源DOA估计方法中，如最为成功的L1-SVD方法[D.Malioutov,M.Cetin,and A.Willsky,“A sparse signal reconstructionperspective for source localization with sensor arrays,”IEEE Trans.Signal Process.,vol.53,no.8,pp.3010–3022,2005.]。其中，L1-SVD方法的计算复杂度比较大，且工程上不容易实现。尽管如此，贝叶斯压缩感知方法结合了以上方法的优点，具有较低的计算复杂度，而且信号的恢复精度很高。另外，压缩感知方法应用于分布源中心DOA估计的方法还很少。

发明内容

本发明的目的在于针对上述已有技术的不足，提供一种基于贝叶斯压缩感知的分布源中心波达方向的估计方法，以达到在多径环境中有效提升小快拍数和低信噪比情况下波达方向的估计性能，应用于无线通信系统可提高无线通信系统的实际应用价值等目的。

本发明方法的流程如图1所示，具体包括以下步骤：

步骤1.设置天线阵列，如图2所示，天线阵列的几何中心为坐标原点，该阵列为M元均匀线阵设于x轴，且每个阵元间隔及阵列之间的距离相同均为d；假定有K个远场窄带信号源入射到该天线阵列上，其中，M≥4，K≥1，d＝λ/2，λ为信号载波的波长；

如图2所示，所述中心波达方向是x轴与信号源中线之间的夹角θ，其取值范围为0°～180°；相对于中心波达方向的角偏差的取值范围为0°～10°；

步骤2.由于信号源在多径环境中传播，天线阵列接收到的数据矢量x(t)可以表示为：

其中，θ_k是第k个入射信号源的中心DOA，且服从对称分布如高斯分布、均匀分布、柯西分布和指数分布等，β_kl(t)为服从高斯分布的随机复增益，L_k为第k个远场窄带信号源的多径信号的总数，s_k(t)为第k个随机信号，w(t)为是加性高斯白噪声，且噪声与信号不相关，Q为数据采样的快拍数；

根据一阶泰勒近似，有：

其中，将所述M元线阵左端的第一个阵元作为空间相位的参考点，则得M×1维阵列导向矢量[·]^T表示矩阵转置运算，a′(θ_k)是a(θ_k)的一阶导数；令由此得到近似模型的矩阵形式为：

X＝Φ(θ,v)S+W＝[A(θ)+A′(θ)V]S+W，X＝[x(1),x(2),...,x(Q)] (3)

其中，v＝[v₁,v₂,…,v_K]，V＝diag(v)，S＝[s(1),s(2),...,s(Q)]，s(t)＝[s₁(t),s₂(t),...,s_K(t)]^T，A(θ)＝[a(θ₁),a(θ₂),…,a(θ_K)]，A′(θ)＝[a′(θ₁),a′(θ₂),…,a′(θ_K)]，diag(·)为构造对角阵运算，W是未知的高斯白噪声矩阵；

步骤3.将上述中心DOA估计问题转化为稀疏矩阵方程求解问题，即构建天线阵列接收数据的系数表达式；

将空域角度空间进行栅格化，如图2所示，即将所述天线阵列上半空间的空域角度(0°,180°]划分成N份并由其中一端的第一个栅格点开始由1编号至N，其中第N个栅格点位于x轴上，每一个栅格点的位置为其相应空域角度的采样位置，其中N不小于K的9倍，则可得到空域角度向量n∈[1，N]，其中为第n个角度采样样本；

利用阵列导向矢量及其一阶导数分别构造M×N维的过完备冗余字典和

\tilde{A} = [a ({\tilde{θ}}_{1}), ..., a ({\tilde{θ}}_{n}), ..., a ({\tilde{θ}}_{N})], {\tilde{A}}^{'} = [a^{'} ({\tilde{θ}}_{1}), ..., a^{'} ({\tilde{θ}}_{n}), ..., a^{'} ({\tilde{θ}}_{N})];

当第k个入射信号源的中心DOA不在所划分的栅格点上时，必然存在一个量化误差ε_g＝θ_k-θ_g，其中，是与空域角度θ_k最近的角度采样样本，g∈{1,2,…,N}；利用一阶泰勒展开，有进而可以得到阵列接收数据的稀疏表达式：

X = \overset{&OverBar;}{Φ} (\tilde{v}, \tilde{ϵ}) \tilde{S} + W = [\tilde{A} + {\tilde{A}}^{'} (\tilde{V} + \tilde{E}) + {\tilde{A}}^{''} \tilde{V} \tilde{E}] \tilde{S} + W - - - (4)

其中，

{\tilde{A}}^{''} = [a^{''} ({\tilde{θ}}_{1}), ..., a^{''} ({\tilde{θ}}_{n}), ..., a^{''} ({\tilde{θ}}_{N})], \tilde{v} = [{\tilde{v}}_{1}, ..., {\tilde{v}}_{n}, ..., {\tilde{v}}_{N}], \tilde{V} = d i a g (\tilde{v}),

\tilde{ϵ} = [{\tilde{ϵ}}_{1}, ..., {\tilde{ϵ}}_{n}, ..., {\tilde{ϵ}}_{N}], \tilde{S} = [\tilde{s} (1), ..., \tilde{s} (t), ..., \tilde{s} (Q)]

是稀疏矩阵，

\tilde{E} = d i a g (\tilde{ϵ}),

由于包含了两个未知稀疏向量和因此，是一个参数化的过完备冗余字典；

步骤4.采用贝叶斯压缩感知方法求解上述稀疏表达式(4)；

步骤4-1.考虑到实际中噪声服从均值为0、方差为σ²的高斯分布，则可假设中的每一个元素都服从均值为0、方差为p_n的先验高斯分布：

f (\tilde{S} | \tilde{p}) ~ Π_{t = 1}^{Q} C N (\tilde{s} (t) | 0, \tilde{P}) - - - (5)

其中，

\tilde{p} = [p_{1}, . . ., p_{n}, . . ., p_{N}]

和

\tilde{P} = diag (\tilde{p});

令n＝1,2,…,N，其中0<ρ<1；

通过最大化的后验概率估计，可以分别得到对应的均值U和协方差矩阵∑：

U = σ^{- 2} Σ \overset{&OverBar;}{Φ} {(\tilde{v}, \tilde{ϵ})}^{H} X - - - (6)

Σ = {[σ^{- 2} \overset{&OverBar;}{Φ} {(\tilde{v}, \tilde{ϵ})}^{H} \overset{&OverBar;}{Φ} (\tilde{v}, \tilde{ϵ}) + {\tilde{P}}^{- 1}]}^{- 1} - - - (7)

其中，(·)^-1为矩阵求逆运算，[·]^H表示矩阵共轭转置运算；

步骤4-2.当U和Σ给定后，采用期望最大化(EM)方法可以分别得到和σ²的更新公式：

p_{n}^{n e w} = \sqrt{Q^{2} + 4 ρ t r (Σ_{n n} + | | U_{n} | |_{2}^{2})} - Q / 2 ρ - - - (8)

{(σ^{2})}^{n e w} \approx {| | X - \overset{&OverBar;}{Φ} (\tilde{v}, \tilde{ϵ}) U | |_{F}^{2} + {(σ^{2})}^{o l d} t r (\overset{&OverBar;}{Ξ} {[{(σ^{2})}^{o l d} I_{M} + \overset{&OverBar;}{Ξ}]}^{- 1})} / M Q - - - (9)

其中，Σ_nn是协方差矩阵Σ的第n个对角元素，U_n是均值U的第n个行向量，||·||₂和||·||_F分别为2范数和Frobenius矩阵范数，tr(·)为矩阵的求迹运算，I_M是M维的单位矩阵；

步骤4-3.假设量化误差向量给定，则稀疏向量可以通过下式得到：

{\tilde{v}}^{n e w} = {(\overset{&OverBar;}{Λ} + {κI}_{N})}^{- 1} \overset{&OverBar;}{c} - - - (10)

\overset{&OverBar;}{c} = 2 Re {\frac{1}{Q} Σ_{t = 1}^{Q} [\tilde{U} * {(B_{1})}^{H} (x (t) - B μ (t))] - d i a g ({(B_{1})}^{H} B Σ)} - - - (11)

其中，⊙为Hardamard积运算，[·]^*表示矩阵共轭运算， μ(t)是U的第t个列向量，κ是一个很小的正实数，10^-12≤κ≤10^-8；

步骤4-4.假设给定，量化误差向量的估计值可以表示为：

{\tilde{ϵ}}^{n e w} = {(\hat{Λ} + {κI}_{N})}^{- 1} \hat{c} - - - (12)

\hat{c} = 2 Re {1 / Q Σ_{t = 1}^{Q} [d i a g (μ {(t)}^{*}) {\tilde{B}}_{1}^{H} (x (t) - \tilde{B} μ (t))] - d i a g ({\tilde{B}}_{1}^{H} \tilde{B} Σ)} - - - (13)

步骤4-5.给定均值U、协方差矩阵∑及的初始值，按式(8)、(9)可得与σ²，根据式(10)、(11)可得再将所得的σ²、带入式(6)、(7)得更新后的U、Σ，由此完成第一次迭代过程；第二次迭代过程中，将上一次更新的均值U、协方差矩阵∑及所得的带入式(8)、(9)、(12)，得更新后的U、Σ、完成第二次迭代，后续迭代过程以此类推实现交替迭代过程；也可给定均值U、协方差矩阵∑及的初始值开始第一次迭代过程，并按上述方法进行交替迭代；

设定最大迭代次数和容错门限τ，所述容错门限τ的一般不大于0.01，对U、Σ、σ²、按上述方法交替迭代，进行多次迭代过程直至迭代次数已达所设定的最大迭代次数或最新一次迭代所得的即与上次迭代所得的即满足如下关系时止：

| | {\tilde{p}}^{n e w} - {\tilde{p}}^{o l d} | |_{2} / | | {\tilde{p}}^{o l d} | |_{2} < τ - - - (15)

步骤5.记步骤4-5最后一次迭代所得的其中分别与第1个、…，第n个、…、第N个采样位置相应的空域角度一一对应，由此构建关于空域角度的曲线，则所得曲线中K个最大峰值对应的K个空域角度即为所述K个远场窄带信号源的中心DOA。

本发明的有益效果是：

1)本发明利用分布源近似模型，构造参数化的过完备冗余字典，首次将分布源中心DOA估计问题转化为稀疏矩阵方程求解问题，采用贝叶斯压缩感知方法对稀疏矩阵进行求解，有效地避免了传统子空间类方法依赖协方差矩阵估计的约束；

2)本发明针对入射的中心DOA不在空域角度采样栅格上的情况，引入了空域角度采样的量化误差，在小角度采样数和小快拍数的情况下，可以获得更好的分布源中心DOA估计精度和分辨率；

3)本发明通过引入角度量化误差，可以在空域角度采样数不大的条件下，有效地解决由空域角度采样造成的分布源中心波达方向估计的固有误差问题；此外，本发明采用贝叶斯压缩感知方法对稀疏矩阵方程进行求解，能够获得未知稀疏向量的最稀疏解；本发明方法计算复杂度低，而且在小快拍数情况下，具有分辨率和精确度高等特点，有利于工程实现。

附图说明

图1为本发明提供的方法流程图；

图2为本发明方法天线阵列设置及空域角度划分示意图；

图3为实施例中四种方法的中心DOA的均方根误差随信噪比变化的曲线图；

图4为实施例中四种方法的中心DOA的均方根误差随快拍数变化的曲线图；

图5为实施例中四种方法的归一化运行时间随快拍数变化的曲线图；

图6为实施例中三种方法的功率谱图。

具体实施方式

下面参照附图和实施例说明本发明的方法实施过程。

实施例

根据图1，本发明的基于贝叶斯压缩感知的分布源中心DOA的估计方法步骤如下：

步骤1.设置天线阵列，如图2所示，该阵列为M＝8元均匀线阵设于x轴，且每个阵元间隔及阵列之间的距离相同均为d，考虑K个远场窄带信号源入射到该阵列上，其中，K＝2，d＝λ/2，λ为信号载波的波长；

步骤2.由于信号源在多径环境中传播，阵列接收到的数据矢量x(t)可以表示为：

其中，θ_k是入射信源的中心DOA，的取值范围是0°～10°，且服从高斯分布，β_kl(t)为服从高斯分布的随机复增益，L_k＝200是第k个分布源的多径信号的总数，s_k(t)为第k个随机信号，w(t)为是加性高斯白噪声，且噪声与信号不相关，Q＝20为数据采样的快拍数；

根据一阶泰勒近似，有：

其中，将所述M元线阵左端的第一个阵元作为空间相位的参考点，得8×1维阵列导向矢量[·]^T表示矩阵转置运算，a′(θ_k)是a(θ_k)的一阶导数，和，可以得到近似模型的矩阵形式为：

X＝Φ(θ,v)S+W＝[A(θ)+A′(θ)V]S+W，

其中，X＝[x(1),x(2),...,x(20)]，v＝[v₁,v₂]，S＝[s(1),s(2),...,s(20)]，s(t)＝[s₁(t),s₂(t),...,s_K(t)]^T，A(θ)＝[a(θ₁),a(θ₂)]，A′(θ)＝[a′(θ₁),a′(θ₂)]，V＝diag(v)，diag(·)为构造对角阵运算，W是未知的高斯白噪声矩阵；

步骤3.将上述中心DOA估计问题转化为稀疏矩阵方程求解问题；

将空域角度空间进行栅格化，如图2所示，即将[0°,180°]划分成N＝90份，那么可以得到其中，为第n个角度采样样本；

利用阵列导向矢量及其一阶导数n＝1,2,…,90，分别构造8×90维的过完备冗余字典和

A (\tilde{θ}) = [a ({\tilde{θ}}_{1}), ..., a ({\tilde{θ}}_{n}), ..., a ({\tilde{θ}}_{90})]

和

{\tilde{A}}^{'} = [a^{'} ({\tilde{θ}}_{1}), ..., a^{'} ({\tilde{θ}}_{n}), ..., a^{'} ({\tilde{θ}}_{90})];

当第k个入射信号源的中心DOA不在所划分栅格点上时，必然存在一个量化误差ε_g＝θ_k-θ_g，其中，是与空域角度θ_k最近的角度采样样本，g∈{1,2,…,90}；利用一阶泰勒展开，有进而可以得到阵列接收数据的稀疏表达式：

X = \overset{&OverBar;}{Φ} (\tilde{v}, \tilde{ϵ}) \tilde{S} + W = [\tilde{A} + {\tilde{A}}^{'} (\tilde{V} + \tilde{E}) + {\tilde{A}}^{''} \tilde{V} \tilde{E}] \tilde{S} + W,

其中，

{\tilde{A}}^{''} = [a^{''} ({\tilde{θ}}_{1}), ..., a^{''} ({\tilde{θ}}_{n}), ..., a^{''} ({\tilde{θ}}_{N})], \tilde{v} = [{\tilde{v}}_{1}, ..., {\tilde{v}}_{n}, ..., {\tilde{v}}_{90}], \tilde{V} = d i a g (\tilde{v}),

\tilde{ϵ} = [{\tilde{ϵ}}_{1}, ..., {\tilde{ϵ}}_{n}, ..., {\tilde{ϵ}}_{90}], \tilde{S} = [\tilde{s} (1), ..., \tilde{s} (t), ..., \tilde{s} (20)]

是稀疏矩阵，

\tilde{E} = d i a g (\tilde{ϵ});

步骤4.采用贝叶斯压缩感知方法求解上述稀疏表达式；

初始化用户参数ρ＝10^-2，κ＝10^-10，最大迭代次数为2000，τ设定为10^-4；

步骤4-1.考虑实际中噪声服从均值为0，方差为σ²的高斯分布，则假设中的每一个元素都服从均值为0，方差为p_n的先验高斯分布：

f (\tilde{S} | \tilde{p}) ~ Π_{t = 1}^{Q} C N (\tilde{s} (t) | 0, \tilde{P}),

其中，

\tilde{p} = [p_{1}, . . ., p_{n}, . . ., p_{N}]

和

\tilde{P} = diag (\tilde{p});

令n＝1,2,…,90；将噪声方差σ²初始化为矩阵XX^H的最小特征值，将始化为其中，mean(·，2)表示按行取平均，量化误差向量和都初始化为零向量；

U = σ^{- 2} Σ \overset{&OverBar;}{Φ} {(\tilde{v}, \tilde{ϵ})}^{H} X

Σ = {[σ^{- 2} \overset{&OverBar;}{Φ} {(\tilde{v}, \tilde{ϵ})}^{H} \overset{&OverBar;}{Φ} (\tilde{v}, \tilde{ϵ}) + {\tilde{P}}^{- 1}]}^{- 1}

其中，(·)^-1为矩阵求逆运算，[·]^H表示矩阵共轭转置运算；

p_{n}^{n e w} = \sqrt{Q^{2} + 4 ρ t r (Σ_{n n} + | | U_{n} | |_{2}^{2})} - Q / 2 ρ

{(σ^{2})}^{n e w} \approx {| | X - \overset{&OverBar;}{Φ} (\tilde{v}, \tilde{ϵ}) U | |_{F}^{2} + {(σ^{2})}^{o l d} t r (\overset{&OverBar;}{Ξ} {[{(σ^{2})}^{o l d} I_{M} + \overset{&OverBar;}{Ξ}]}^{- 1})} / M Q

{\tilde{v}}^{n e w} = {(\overset{&OverBar;}{Λ} + {κI}_{N})}^{- 1} \overset{&OverBar;}{c}

\overset{&OverBar;}{c} = 2 Re {\frac{1}{Q} Σ_{t = 1}^{Q} [{\tilde{U}}^{*} {(B_{1})}^{H} (x (t) - B μ (t))] - d i a g ({(B_{1})}^{H} B Σ)}

其中，⊙为Hardamard积运算，[·]^*表示矩阵共轭运算，

B = \tilde{A} + {\tilde{A}}^{'} \tilde{E}

和

B_{1} = {\tilde{A}}^{'} + {\tilde{A}}^{''} \tilde{E},

μ(t)是U的第t个列向量；

步骤4-4.假设给定，量化误差向量的估计值可以表示为：

{\tilde{ϵ}}^{n e w} = {(\hat{Λ} + {κI}_{N})}^{- 1} \hat{c}

\hat{c} = 2 Re {1 / Q Σ_{t = 1}^{Q} [d i a g (μ {(t)}^{*}) {\tilde{B}}_{1}^{H} (x (t) - \tilde{B} μ (t))] - d i a g ({\tilde{B}}_{1}^{H} \tilde{B} Σ)}

步骤4-5.对U，Σ和σ²，交替迭代，直到或达到最大迭代次数；

步骤5.记步骤4-5最后一次迭代所得的其中分别与第1个、…，第n个、…、第N个采样位置相应的空域角度一一对应，由此构建关于空域角度的曲线，则所得曲线中K个最大峰值对应的K个空域角度即为所述K个远场窄带信号源的中心DOA，如：

{\hat{θ}}_{k} = {\tilde{θ}}_{j_{k}} + {\tilde{ϵ}}_{j_{k}}, k = 1, 2,

其中，j_k∈{1,2,…,90}，k＝1,2是的两个最大峰值所在的坐标。

为了验证本发明的性能，进行如下三个仿真试验。

试验一：考察两个分布源的中心DOA估计的均方根误差随信噪比和快拍数变化的性能。

两个分布源的中心DOA分别为61.8°和92.2°，角度扩展都设置为1°，阵元数为8的情况下，进行400次独立的Monte Carlo仿真。均方根误差的计算式定义为：

{RMSE}_{θ} = \sqrt{\frac{1}{400 K} Σ_{r = 1}^{400} Σ_{k = 1}^{K} {({\hat{θ}}_{r, k} - θ_{k})}^{2}}

其中，和θ_k分别表示估计得到的中心DOA和入射中心DOA，r＝1,2,...,400和k＝1,2,...,K。

图3表示采样快拍数Q＝20时，中心DOA估计的均方根误差随信噪比变化的情况。从图中可以看出，在整个信噪比范围内，本发明方法所获得的中心DOA估计的均方根误差小于其它相比较的方法。

图4表示信噪比为10dB时，中心DOA估计的均方根误差随快拍数变化的情况。从图中可以看出，在小快拍数的情况下，本发明方法所获得的中心DOA估计的均方根误差明显小于其它方法。

试验二：考察三个分布源入射到阵列的计算复杂度。

图5表示信源数为3，阵元数为10时，归一化运行时间随快拍数变化的情况。从图5中可以看出，随着快拍数的增加，本发明方法的运行时间小于DISPARE方法和L1-SVD方法。

试验三：考察两个靠得很近的中心DOA的分布源功率谱。

两个分布源的中心DOA分别为62.6°和69.4°，角度扩展都设置为1°，阵元数为8，信噪比为10dB，快拍数Q＝100。

图6表示三种方法的两个分布源的功率谱情况。

从图6中可以看出，本发明方法所提供的功率谱有两个峰值，能够对两个靠得很近的分布源中心DOA进行估计，而DISPARE方法和L1-SVD方法失效，只能给出一个峰值。因此，本发明方法的分辨率高于DISPARE方法和L1-SVD方法。

综上，本发明方法不仅降低了分布源中心DOA估计的运行时间，还提高了中心DOA估计的精度和分辨率，更适用于工程实现。

Claims

1.一种基于贝叶斯压缩感知的分布源中心波达方向估计方法，具体包括以下步骤：

步骤1.设置天线阵列，天线阵列的几何中心为坐标原点，该阵列为M元均匀线阵设于x轴，且每个阵元间隔及阵列之间的距离相同均为d；假定有K个远场窄带信号源入射到该天线阵列上，其中，M≥4，K≥1，d＝λ/2，λ为信号载波的波长；

所述中心波达方向是x轴与信号源中线之间的夹角θ，其取值范围为0°～180°；相对于中心波达方向的角偏差的取值范围为0°～10°；

其中，θ_k是第k个入射信号源的中心DOA，服从对称分布，β_kl(t)为服从高斯分布的随机复增益，L_k为第k个远场窄带信号源的多径信号的总数，s_k(t)为第k个随机信号，w(t)为是加性高斯白噪声，且噪声与信号不相关，Q为数据采样的快拍数；

根据一阶泰勒近似，有：

其中，v＝[v₁，v₂，...，v_K]，V＝diag(v)， A(θ)＝[a(θ₁),a(θ₂),...,a(θ_K)]，A′(θ)＝[a′(θ₁),a′(θ₂),...,a′(θ_K)]，diag(·)为构造对角阵运算，W是未知的高斯白噪声矩阵；

将空域角度空间进行栅格化：将所述天线阵列上半空间的空域角度(0°,180°]划分成N份并由其中一端的第一个栅格点开始由1编号至N，其中第N个栅格点位于x轴上，每一个栅格点的位置为其相应空域角度的采样位置，其中N不小于K的9倍，则可得到空域角度向量n∈[1，N]，其中为第n个角度采样样本；

\tilde{A} = [a ({\tilde{θ}}_{1}), ..., a ({\tilde{θ}}_{n}), ..., a ({\tilde{θ}}_{N})], {\tilde{A}}^{'} = [a^{'} ({\tilde{θ}}_{1}), ..., a^{'} ({\tilde{θ}}_{n}), ..., a^{'} ({\tilde{θ}}_{N})];

当第k个入射信号源的中心DOA不在所划分的栅格点上时，必然存在一个量化误差ε_g＝θ_k-θ_g，其中，是与空域角度θ_k最近的角度采样样本，g∈{1,2,...,N}；利用一阶泰勒展开，有进而可以得到阵列接收数据的稀疏表达式：

X = \overset{&OverBar;}{Φ} (\tilde{v}, \tilde{ϵ}) \tilde{S} + W = [\tilde{A} + {\tilde{A}}^{'} (\tilde{V} + \tilde{E}) + {\tilde{A}}^{''} \tilde{V} \tilde{E}] \tilde{S} + W - - - (4)

其中，

{\tilde{A}}^{''} = [a^{''} ({\tilde{θ}}_{1}), ..., a^{''} ({\tilde{θ}}_{n}), ..., a^{''} ({\tilde{θ}}_{N})], \tilde{v} = [{\tilde{v}}_{1}, ..., {\tilde{v}}_{n}, ..., {\tilde{v}}_{N}], \tilde{V} = d i a g (\tilde{v}),

\tilde{ϵ} = [{\tilde{ϵ}}_{1}, ..., {\tilde{ϵ}}_{n}, ..., {\tilde{ϵ}}_{N}],

是稀疏矩阵，

\tilde{E} = d i a g (\tilde{ϵ}),

步骤4.采用贝叶斯压缩感知方法求解上述稀疏表达式(4)；

f (\tilde{S} | \tilde{p}) Π_{t = 1}^{Q} C N (\tilde{s} (t) | 0, \tilde{P}) - - - (5)

其中，和令n＝1，2，...，N，其中0<ρ<1；

U = σ^{- 2} Σ \overset{&OverBar;}{Φ} {(\tilde{v}, \tilde{ϵ})}^{H} X - - - (6)

Σ = {[σ^{- 2} \overset{&OverBar;}{Φ} {(\tilde{v}, \tilde{ϵ})}^{H} \overset{&OverBar;}{Φ} (\tilde{v}, \tilde{ϵ}) + {\tilde{P}}^{- 1}]}^{- 1} - - - (7)

其中，(·)^-1为矩阵求逆运算，[·]^H表示矩阵共轭转置运算；

步骤4-2.当U和Σ给定后，采用期望最大化方法可以分别得到和σ²的更新公式：

p_{n}^{n e w} = \sqrt{Q^{2} + 4 ρ t r (Σ_{n n} + | | U_{n} | |_{2}^{2})} - Q / 2 ρ - - - (8)

{(σ^{2})}^{n e w} \approx {| | X - \overset{&OverBar;}{Φ} (\tilde{v}, \tilde{ϵ}) U | |_{F}^{2} + {(σ^{2})}^{o l d} t r (\overset{&OverBar;}{Ξ} {[{(σ^{2})}^{o l d} I_{M} + \overset{&OverBar;}{Ξ}]}^{- 1})} / M Q - - - (9)

{\tilde{v}}^{n e w} = {(\overset{&OverBar;}{Λ} + {κI}_{N})}^{- 1} \overset{&OverBar;}{c} - - - (10)

\overset{&OverBar;}{c} = 2 Re {\frac{1}{Q} Σ_{t = 1}^{Q} [{\tilde{U}}^{*} {(B_{1})}^{H} (x (t) - B μ (t))] - d i a g ({(B_{1})}^{H} B Σ)} - - - (11)

其中，为Hardamard积运算，[·]^*表示矩阵共轭运算，和μ(t)是U的第t个列向量，κ为正实数且10^-12≤κ≤10^-8；

步骤4-4.假设给定，量化误差向量的估计值可以表示为：

{\tilde{ϵ}}^{n e w} = {(\hat{Λ} + {κI}_{N})}^{- 1} \hat{c} - - - (12)

\hat{c} = 2 Re {1 / Q Σ_{t = 1}^{Q} [d i a g (μ {(t)}^{*}) {\tilde{B}}_{1}^{H} (x (t) - \tilde{B} μ (t))] - d i a g ({\tilde{B}}_{1}^{H} \tilde{B} Σ)} - - - (13)

步骤4-5.设定最大迭代次数和容错门限τ，所述容错门限τ的一般不大于0.0，1对U、Σ、σ²、进行多次交替迭代过程直至迭代次数已达所设定的最大迭代次数或最新一次迭代所得的即与上次迭代所得的即满足如下关系时止：

| | {\tilde{p}}^{n e w} - {\tilde{p}}^{o l d} | |_{2} / | | {\tilde{p}}^{o l d} | |_{2} < τ - - - (15)

2.根据权利要求1所述的基于贝叶斯压缩感知的分布源中心波达方向估计方法，其特征在于，所述的相对于中心波达方向的角偏差的分布形式是高斯分布、均匀分布、柯西分布、指数分布中的一种。

3.根据权利要求1所述的基于贝叶斯压缩感知的分布源中心波达方向估计方法，其特征在于，所述的接收数据的快拍数Q的范围是10～200次。

4.根据权利要求1所述的基于贝叶斯压缩感知的分布源中心波达方向估计方法，其特征在于，步骤4-5所述交替迭代过程具体如下：

给定均值U、协方差矩阵∑及的初始值，按式(8)、(9)可得与σ²，根据式(10)、(11)可得再将所得的σ²、带入式(6)、(7)得更新后的U、Σ，由此完成第一次迭代过程；第二次迭代过程中，将上一次更新的均值U、协方差矩阵∑及所得的带入式(8)、(9)、(12)，得更新后的U、Σ、完成第二次迭代，后续迭代过程以此类推实现交替迭代过程；或者，可给定均值U、协方差矩阵∑及的初始值开始第一次迭代过程，并按上述方法进行交替迭代。