）、（9)可得6与〇2,根 据式（10)、（11)可得夺，再将所得的P、〇2、v、E带入式化）、（7)得更新后的U、2,由此 完成第一次迭代过程；第二次迭代过程中，将上一次更新的均值U、协方差矩阵E及所得的 V带入式巧）、巧）、（12)，得更新后的U、2、g完成第二次迭代，后续迭代过程W此类推实现 交替迭代过程；也可给定均值U、协方差矩阵E及V的初始值开始第一次迭代过程，并按上 述方法进行交替迭代；
[0044] 设定最大迭代次数和容错口限T，所述容错口限T的一般不大于0.01，对U、Z、 P，t，2、不、￡按上述方法交替迭代，进行多次迭代过程直至迭代次数已达所设定的最大迭 代次数或最新一次迭代所得的P即护6"与上次迭代所得的P即护"满足如下关系时止；
[0045]
(15)
[004引步骤5.记步骤4-5最后一次迭代所得的护ew=[pr，...，pr'，...，pr]，其中 p^'，…，p^'，...，p^'分别与第l个、…，第n个、…、第N个采样位置相应的空域角度一一对 应，由此构建护6W关于空域角度的曲线，则所得曲线中K个最大峰值对应的K个空域角度即 为所述K个远场窄带信号源的中屯、D0A。
[0047] 本发明的有益效果是：
[0048] 1)本发明利用分布源近似模型，构造参数化的过完备冗余字典，首次将分布源中 屯、D0A估计问题转化为稀疏矩阵方程求解问题，采用贝叶斯压缩感知方法对稀疏矩阵进行 求解，有效地避免了传统子空间类方法依赖协方差矩阵估计的约束；
[0049] 2)本发明针对入射的中屯、D0A不在空域角度采样栅格上的情况，引入了空域角度 采样的量化误差，在小角度采样数和小快拍数的情况下，可W获得更好的分布源中屯、D0A 估计精度和分辨率；
[0050] 3)本发明通过引入角度量化误差，可W在空域角度采样数不大的条件下，有效地 解决由空域角度采样造成的分布源中屯、波达方向估计的固有误差问题；此外，本发明采用 贝叶斯压缩感知方法对稀疏矩阵方程进行求解，能够获得未知稀疏向量的最稀疏解；本发 明方法计算复杂度低，而且在小快拍数情况下，具有分辨率和精确度高等特点，有利于工程 实现。
【附图说明】
[0051] 图1为本发明提供的方法流程图；
[0052] 图2为本发明方法天线阵列设置及空域角度划分示意图；
[0053] 图3为实施例中四种方法的中屯、D0A的均方根误差随信噪比变化的曲线图；
[0054] 图4为实施例中四种方法的中屯、D0A的均方根误差随快拍数变化的曲线图；
[0055] 图5为实施例中四种方法的归一化运行时间随快拍数变化的曲线图；
[0056] 图6为实施例中S种方法的功率谱图。
【具体实施方式】
[0057] 下面参照附图和实施例说明本发明的方法实施过程。
[00则实施例
[0059] 根据图1，本发明的基于贝叶斯压缩感知的分布源中屯、D0A的估计方法步骤如下：
[0060] 步骤1.设置天线阵列，如图2所示，该阵列为M= 8元均匀线阵设于X轴，且每个 阵元间隔及阵列之间的距离相同均为山考虑K个远场窄带信号源入射到该阵列上，其中，K =2,d=A/2,A为信号载波的波长；
[0061] 步骤2.由于信号源在多径环境中传播，阵列接收到的数据矢量x(t)可W表示 为：
[0062]
[006引其中，马i=0k+免(t)，0k是入射信源的中屯、D0A，免阳的取值范围是0。~10。， 且4(0服从高斯分布，0kia)为服从高斯分布的随机复增益，Lk= 200是第k个分布源的 多径信号的总数，Sk(t)为第k个随机信号，w(t)为是加性高斯白噪声，且噪声与信号不相 关，Q= 20为数据采样的快拍数；
[0064] 根据一阶泰勒近似，有：
[0065]
[0066] 其中，将所述M元线阵左端的第一个阵元作为空间相位的参考点，得8X1维阵列 导向矢量a倘）= [l，パw?s6、"，…，e^/WM-ud?s4/叩，[?]嗦示矩阵转置运算，a' ( 0k)是 a( 0k)的一阶导数，和，可W得到近似模型的矩阵形式为：
[0067] X= 〇 (日，v)S+W= [A(日）+A'(日）V]S+W，
[0068] 其中，X= [x(l),x(2), . . . ,x(20)]，V= [Vi,V2]，S= [s(l),s(2), . . . ,s(20)]， s(t) = [Si(t)，S2(t)，...，SK(t)]T，〇>(0，、')=恤，(p?]，A(日）=[a(日i)，a(日2)]，A' (0) = [a' (0i)，a' (02)]，V=diag(v)，diag(〇为构造对角阵运算，W是未知的高斯白噪 声矩阵；
[0069] 步骤3.将上述中屯、D0A估计问题转化为稀疏矩阵方程求解问题；
[0070] 将空域角度空间进行栅格化，如图2所示，即将[0。，180° ]划分成N= 90份， 那么可W得到0=的，…，毎，...4。]，其中，&为第n个角度采样样本；
[00川利用阵列导向矢量a短）及其一阶导数a喊），n= 1，2,…，90,分别 构造8X90维的过完备冗余字典A和A'：A曲二[a的)，...，a(^)，...，a(4：)WP A'= [a'(各I)，...，a'(与,)，...，a'(每…)];
[0072] 当第k个入射信号源的中屯、DOA不在所划分栅格点上时，必然存在一个量化 误差Eg= 0k-0g，其中，马e.庇每，…，成,}是与空域角度0k最近的角度采样样本，gG{1，2,…，90};利用一阶泰勒展开，有9k 巧:) +a'峰)[Vk+￡g] +a"的)v'kV进而可W 得到阵列接收数据的稀疏表达式：
[0073]
[0074] 其中，又"=[3"(句)，...，3"诉)，...，3"成：)]，^ =内，...，巧。_.，[,,,]， ￥=化巧(;'）， E=[马，…，毎，…，马0]，8=阶1)，...，岭)，...，§(2〇11是稀疏矩阵，色=出江3倒；由于&仇巧包含了 两个未知稀疏向量豕和￡，因此，东代，引是一个参数化的过完备冗余字典；
[00巧]步骤4.采用贝叶斯压缩感知方法求解上述稀疏表达式；
[0076]初始化用户参数P=l(r2,K= 1〇4°，最大迭代次数为2000,T设定为10^4;
[0077] 步骤4-1.考虑实际中噪声服从均值为0,方差为0 2的高斯分布，则假设岭)中的 每一个元素都服从均值为0,方差为P。的先验高斯分布：
[0078]
[0079] 其中，p= [P|?...，p,,，...，pN]和Ndi￡lg巧）；令./'(戶）--^(p,,|l，ク)，n=l，2，…，90;将 噪声方差〇2初始化为矩阵XXH的最小特征值，将P始化为meaii(|&H快巧X|，2)，其中， mean(?，2)表示按行取平均，量化误差向量g和^都初始化为零向量；
[0080] 通过最大化§的后验概率估计，可W分别得到对应的均值U和协方差矩阵E:
[0083] 其中，（?r为矩阵求逆运算，[?]H表示矩阵共辆转置运算；
[0084] 步骤4-2.当U和S给定后，采用期望最大化（EM)方法可W分别得到P和0 2的 更新公式：
[0085]
[0086]
[0087] 其中，S= &(f，g皆&仇g)H，是协方差矩阵S的第n个对角元素，U。是均值U 的第n个行向量，II?II2和II?IIP分别为2范数和化obenius矩阵范数，tr( ?)为矩阵 的求迹运算，Im是M维的单位矩阵；
[0088] 步骤4-3.假设量化误差向量g给定，则稀疏向量乐可W通过下式得到：
[0091] 其中，?为Hardamard积运算，[? ]*表示矩阵共辆运算，A=(B,)H〇ii)'〇(UUH+i:)， B=A+A'丘和Bi二心+心丘，y(t)是U的第t个列向量；
[0092] 步骤4-4.假设f给定，量化误差向量g的估计值可W表示为：
[009引步骤4-5.对U，S和目，02，^，g交替迭代，直到
或达到 最大迭代次数；
[0097] 步骤5.记步骤4-5最后一次迭代所得的护=[pr"，...，pr，...，pr]，其中 pr，...，pr，...，pr分别与第1个、…，第n个、…、第N个采样位置相应的空域角度一一对 应，由此构建r"''关于空域角度的曲线，则所得曲线中K个最大峰值对应的K个空域角度即 为所述K个