基于广义近似消息传递的快速DOA估计方法

张 俊，张新禹，姜卫东，刘永祥，黎 湘

（国防科技大学电子科学学院，湖南 长沙410073）

0 引 言

波达方向（direction of arrival，DOA）估计是阵列信号处理研究的一类重点问题，其广泛应用于雷达、声纳、无线通信、地震及水下探测等众多领域。经过近几十年的发展，研究人员提出了大量高分辨DOA估计算法。但是这些传统的DOA估计算法在实际应用中受到诸多限制。例如，以多重信号分类方法（multiple signal classification，MUSIC）为代表的子空间分解类算法通常需要较高的信噪比（signal to noise ratio，SNR）、足够的快拍数量以及已知信源数量。子空间拟合类算法，如最大似然算法和加权子空间拟合算法等，虽然在低SNR和小快拍数的情况下表现比子空间分解类算法好，但因其有计算量大的缺点而难以实际应用。

近年来，稀疏恢复类DOA算法由于在低SNR、有限快拍数和相关信号等场景中的良好性能，成为研究热点。此类算法的基本思想是在满足一定约束条件下，利用过完备的基函数求解可以表示原始信号的稀疏系数向量。稀疏恢复类算法主要包括贪婪追踪、l （≥0）范数优化和稀疏贝叶斯学习（sparse Bayesian learning，SBL）三类算法。贪婪追踪类算法是在每次迭代中根据相关性选择与当前残余信号特征结构匹配的最优原子，实现对原始信号的逼近，典型算法有匹配追踪算法、正交匹配追踪算法等。虽然此类算法不需要信源数量的先验信息，也可用于处理相关信号，但是存在计算复杂度高、收敛速度慢等问题。l 范数优化算法是将l 范数惩罚项引入目标函数，保证解的稀疏性，防止出现过拟合，常用的算法有LASSO（least absolute shrinkage and selection operator）算法、-奇异值分解（singular value decomposition，SVD）算法等。此类算法通常也具有较高的计算复杂度，同时正则化参数的选择对性能的影响比较明显。SBL类方法是从贝叶斯观点出发，利用信号的稀疏先验信息对信号进行稀疏重构。相对于上述两类方法，SBL算法能够保证达到全局最优且计算效率也有显著提高。同时，SBL算法的在信号分量之间存在强相关性时仍然具有良好的重构性能。文献［6-7］提出了基于SBL的双基地无源雷达DOA估计方法，其思想主要是通过空间几何关系建立双站接收数据之间的关系，联合进行DOA估计，在一定程度上缓解了快拍数的限制，性能要优于只使用单站数据的情况，但是在计算后验概率密度时涉及到阶矩阵求逆问题，需要在计算效率和估计精度之间做出平衡。为避免大规模矩阵求逆问题，文献［8］提出了一种基于酉变换的近似消息传递（approximate message passing，AMP）的SBL-DOA估计算法。但是AMP算法只能解决线性问题，不适用于任意的信号分布，且要求字典矩阵满足独立同分布的零均值高斯分布。文献［12］将AMP算法进行了拓展，提出广义AMP（generalized AMP，GAMP）算法。GAMP算法对输入、输出信号的分布没有要求，且可以解决非线性问题，适用于具有精确渐近性能保证的非高斯估计问题。文献［13］进一步证明了在独立高斯先验和似然函数的条件下，当字典矩阵奇异值的峰—均值比较小时，可以保证GAMP算法收敛，并给出了一种阻尼GAMP算法，使得在上述条件不满足时也可以保证算法收敛，极大地拓展了算法的适用范围。

本文提出一种SBL框架下的基于GAMP的双基地联合DOA估计方法。首先，构建了双站快拍数据下的信号模型。通过数据间的独立性，利用GAMP网络将多维联合后验概率密度的计算简化为标量形式的边缘概率密度计算，并推导了参数更新的闭式表达式。然后，针对离网目标，提出了一种基于梯度的网格节点更新策略。仿真实验表明，在低SNR和有限快拍数的情况下，本算法相比于其他算法具有较高的DOA估计精度和较低的计算复杂度。

1 双基地无源雷达接收信号模型

1.1 系统模型

图1 双基地无源雷达系统空间几何关系Fig.1 Spatial geometry of the bi-static passive radar system

根据空间几何关系，容易得到第个目标在两个天线阵列上的入射角满足：

化简可得

因此，通过式（2）就可以建立两个雷达之间接收数据的关系。两个雷达站多快拍模型的接收数据可以写成

1.2 贝叶斯框架下的稀疏信号模型

对于信号的先验概率，根据目标的独立性假设可以认为、中的元素是相互独立的，且每一行的元素是独立同分布的，服从均值为零、方差为γ≥0（＝1，2，…，）的复高斯分布。这里假设入射信号在两个阵列处的方差相同，这样设定主要是为了简化模型复杂度。实际应用中，可以根据雷达方程求解出两个阵列接收到同一目标信号的功率比值，故可以通过补偿使两处的接收信号功率一致。故信号幅值的先验概率可写为

式中：x 表示X （＝1，2）中第行列的元素；（·）为冲击函数。那么接收信号的先验概率密度可表示为

在SBL框架下，本文假设这些超参数服从伽马分布：

2 基于SBL的快速DOA估计算法

2.1 利用GAMP算法求解信号幅值的后验概率

根据模型中的独立性假设及传输信道的湍流效应，不同阵列之间、不同信源之间信号幅值都是相互独立的，故入射信号的后验概率可以分解为

式中：、（＝1，2，…，）分别表示、的第行。由于、的先验分布为复高斯分布，根据共轭先验的性质，其后验分布也应为复高斯分布。故只需求出均值和协方差矩阵即可获得后验概率密度。

根据各样本之间的独立性，有

式中：、 （＝1，2，…，）分别是、中的第个元素。

通过式（11）可以看出，矢量形式的后验概率密度转化成了一系列标量后验概率密度的乘积。这里可以采用文献［12］提出的GAMP算法计算、的后验概率密度，具体步骤如算法1所示。

2.2 估计γ

按照标准EM算法的操作，的最大后验估计为

对于式（19）的求解一般是构建ln（（，，））的紧凑下界（），根据Jensen不等式可得

式中：［·］表示在概率密度（）上求期望。

将式（20）代入式（19）并写成标量形式，忽略（）中与γ（＝1，2，…，）无关的项，可以得到

式中：目标函数结合式可以写为

在无信息先验的情况中，是一个很小的值，忽略并令式（23）等于零，可以得到γ的迭代更新表达式：

经实验验证，式（25）比式（24）收敛更快。

2.3 估计目标数量K和噪声功率 和

一般情况下，目标数量可通过搜索的峰值来估计，但是在样本数较少、SNR较低的条件下往往会出现多个虚假尖峰。参考文献［15］的方法，本文采用如下表达式对进行估计：

准确的噪声功率估计可以加快算法的收敛速度，然而通过类似式（19）的方式求解也会涉及阶矩阵求逆的运算，不但计算复杂度较高，而且估计结果比实际值要小。因此，本算法采用下列实证方法来估计噪声功率：

式中：和可能在每次迭代中会改变，但是在角度估计正确的条件下，上式是对噪声功率的一致无偏渐近估计。当→∞时，该估计的方差趋近于克拉美罗下限。

2.4 算法初始值的选择

虽然本算法对初始值不敏感，但是一个好的初始值可以加快算法的收敛速度。在已知SNR的情况下，噪声功率的初始值可设为

如果缺少SNR的信息，噪声功率也可初始化为一个较小的值。因为当角度估计不准时，通过式（29）对噪声功率的估计值会偏高。

信号功率可通过数据协方差矩阵来进行初始化。初始时假设所有方向的信号功率均为，即＝1，其中

2.5 对离网目标的处理

从模型的构建可以看出，角度网格的密度决定着本算法DOA估计的精度。对于离网目标，文献［19］给出了一个简单的网格更新策略，但是会增加网格的数量和密度，从而使计算复杂度提高；文献［10］根据阵列流行矩阵的范德蒙德结构提出了Root-SBL方法，但此方法只适用于均匀线阵，不具有普遍性。因此，本算法在噪声功率估计完成后，考虑使用梯度方法求解目标所处网格的偏移量来实现对网格位置的动态调整。

最优的网格分布应满足如下条件：

式中：上标为网格更新步骤中的迭代次数；sign（·）为符号函数；为更新步长。使用单站数据分别计算对同一目标的网格更新，然后换算成一个站的角度，取平均值作为最终的网格更新。

2.6 算法复杂度分析

最后，对本文所提算法的具体实现步骤总结如下。

3 仿真分析

本节通过数值仿真的方式对所提出的基于GAMP的快速SBL算法（GAMP-SBL）进行有效性验证，并从空间谱、不同SNR和快拍数下的均方根误差以及计算时间等4个方面与经典的Lasso算法、BM-SBL算法、M-SBL算法进行对比分析。其中，LASSO算法采用CVX工具箱求解，M-SBL算法默认已知目标数量，故该算法没有估计这一步骤。

（1）仿真1：算法空间谱分析

图2为各算法在阵列1上的空间谱。实验条件为SNR＝-7 dB，SNR＝-10 dB，图2（a）快拍数为40，图2（b）快拍数为10。从图2（a）中可以看出，各算法均能正确识别出两个目标，本算法与BM-SBL算法的波束宽度基本一致，M-SBL算法波束稍宽。当快拍数进一步降低，如图2（b）所示，BM-SBL算法副瓣起伏较大，M-SBL算法只估计出了10°方向的目标，LASSO算法性能最好，但这是调整正则化参数后的结果，本算法虽然也出现了峰值较低的副瓣，但不影响目标方位的估计，说明了本算法在小快拍数条件下的稳健性。

图2 算法空间谱分析Fig.2 Analysis of algorithm spatial spectrum

（2）仿真2：算法DOA估计精度分析

图3为不同SNR下各算法的估计精度。各算法所用快拍数为40，站1的SNR从-10 dB变化到2 dB。从图3中可以看出，各算法的估计误差均随SNR增大而减小，本算法与M-SBL算法的误差基本一致，均明显小于BM-SBL和LASSO算法。

图3 SNR对不同算法估计精度的影响Fig.3 Influence of SNR on the estimation accuracy of different algorithms

图4为快拍数对算法估计精度的影响。SNR为-10 dB，快拍数为10到150，其中M-SBL算法在10个快拍下不能正确估计，故只统计了快拍数大于20的情况。由图4可知，随快拍数的增多，各算法的估计误差呈下降趋势，BM-SBL的误差下降最快，本算法在小快拍时的估计误差最小，当快拍数较大时，估计误差与BM-SBL算法接近。

图4 快拍数对不同算法估计精度的影响Fig.4 Influence of snapshot number on the estimation accuracy of different algorithms

（3）仿真3：算法计算时间分析

平均计算时间可以体现算法的计算复杂度，本实验中各算法均不采用角度网格更新策略。

图5为不同SNR下各算法的平均计算时间对比，实验条件同图3实验。可以看出，除BM-SBL算法外，其他算法的平均计算时间基本不随SNR变化。另外，本算法的平均计算时间远小于其他算法，一定程度上说明本算的计算复杂度要小于其他算法。

图5 不同SNR下各算法的平均计算时间Fig.5 Average calculation time of each method with different SNR

图6为不同快拍下各算法的平均计算时间对比，实验条件同图5实验。从图6中可以看出，各算法的平均计算时间均随快拍数的增加而变大，但是M-SBL、BM-SBL和本算法上升的比较平缓，LASSO算法的增加速度最快，受快拍数的影响最大。本算法平均计算时间仍然远小于其他算法。对比图5和图6可以看出，快拍数对计算时间的影响较大。

图6 不同快拍下各算法的平均计算时间Fig.6 Average calculation time of each method with different number of snapshots

（4）仿真4：角度空间网格更新分析

为更好地衡量角度更新策略的有效性，实验设置空间目标源在阵列1上的入射角度为-5°和5.4°，其他实验条件同图3实验，此时无网格更新算法的理论均方根误差为0.28°。

图7为不同SNR下算法采用网格更新策略前后的估计误差。从图7中可以看出，在SNR较小时，使用更新策略后会导致误差增大，但是当SNR较大时，更新策略会明显降低误差，突破0.28°的理论误差，说明本算法提出的网格更新策略是有效的。对比图7与图3，虽然在小SNR时误差有所上升，但依然比其他算法要低。

图7 不同SNR下网格更新前后算法的估计误差Fig.7 RMSE of the proposed algorithm versus SNR before and after the refinement

图8为在上述实验条件下，算法采用网格更新策略前后的平均计算时间。从图8中可知，算法的平均计算时间对SNR的变化不敏感，增加网格更新策略后，算法计算时间有所提高，但与图5对比可知仍然显著低于其他算法。

图8 不同SNR下网格更新前后算法的平均计算时间Fig.8 Average calculation time of the proposed algorithm versus different SNR before and after the refinement

4 结 论

本文提出一种基于GAMP算法改进的双基地无源雷达SBL的DOA估计方法。该算法根据双站雷达接收数据之间的相互独立性，构建了适用于GAMP算法的信号模型，将SBL中后验概率密度的计算转变为标量边缘概率密度的乘积，避免了大规模矩阵求逆的计算。同时针对离网目标，推导了基于梯度下降的角度网格空间更新方法。仿真结果表明，与经典的LASSO、M-SBL以及BM-SBL算法相比，本算法在低SNR和有限快拍数条件下估计误差低、稳健性好，尤其是计算复杂度上有显著的优势。但本算法主要针对的是相互独立的目标信号，后续工作重点将是相关性目标信号下的算法改进及应用研究。