一种基于稀疏贝叶斯学习的快速目标角度估计方法与流程

本发明属于阵列信号处理领域，具体涉及一种基于稀疏贝叶斯学习的快速目标角度估计方法。

背景技术：

目标角度估计(directionofarrivalestimation,doae)是阵列信号处理领域中非常重要的一个问题，相关的算法可广泛应用于雷达探测，声纳导航，多信道通信等领域。传统的目标角度估计算法大多是基于子空间信号流型分析或者属于最大似然估计类算法。但这两类算法往往分辨率有限，且受目标信号间相关性，采集样本数以及系统信噪比等因素影响较大，因此难以适应实际中各类复杂的应用环境。

随着稀疏恢复算法的引入，目标角度估计技术得到了很大的提升。稀疏恢复算法将传统的目标角度估计模型等效成字典学习模型，通过各类压缩感知算法得到更为精确的角度估计值。采用基于稀疏恢复算法的目标角度估计可以使得目标角度的获取对目标信号相关性，采集样本数以及信噪比的容忍度都大大提升，因此更有利于实际的应用。但是由于稀疏恢复算法需要进行较大规模的矩阵运算，算法的运算量随着所需要的估计精度成指数倍增长，这一点又限制了算法的应用。

稀疏贝叶斯学习(sparsebayesianlearning,sbl)是稀疏恢复算法中的一种。它将稀疏恢复模型中的参数估计纳入到贝叶斯推导的框架中，利用模型中的变量间的统计分布规律，通过反复迭代学习获得对参数的精确估计。

考虑一个由m个天线组成的阵列收到来自不同方向θ＝[θ1θ2…θk]^t的k个目标的回波信号，这里[…]^t表示矩阵的转置。sbl算法的信号模型可以写为y(t)＝ax(t)+n(t)(1)

其中y(t)表示t时刻的阵列接收信号，表示字典矩阵，表示阵列的导向矢量，阵列的导向矢量表示远场条件下的平面波入射到阵列所形成的接收信号矢量，它是一个复向量，维数等于阵列的阵元个数，对于一个由m个天线组成的均匀线阵，导向矢量可以表示成：

其中d表示阵元之间的间距，λ表示阵列发射的电磁波波长，π为圆周率常数，θ表示任一入射角度。字典矩阵a中的组成了一个关于空间入射角的一种网格划分，该网格划分越细，目标角度估计的分辨率往往就越高，一般的做法是将一个可能的空域区间做等间隔的划分，n的值会远大于阵列阵元个数m。(1)中的x(t)表示接收到的目标信号矢量，它的元素与字典矩阵中的每一列一一对应，当其所对应的入射角度存在目标时，x(t)的该元素等于信号的复数值；当其所对应的角度不存在目标时，x(t)的该元素值为0。n(t)表示系统的加性噪声。这里我们假设在不同时刻的目标信号统计独立，这一假设在雷达信号处理中十分常见。在实际工程应用中，n(t)一般假设为零均值的高斯白噪声，即n(t)～cn(0,σ0i)，并且不同时刻间的噪声统计独立，这里cn(0,σ0i)表示复数零均值的复高斯分布，其方差为σ0i，σ0表示噪声的功率，i表示单位矩阵；下文中用cn(α,β)表示均值为α，方差为β的复高斯分布函数。

根据以上信号模型，我们可以得到阵列的回波信号满足分布y(t)～cn(ax(t),σ0i)。为了进行贝叶斯推导，sbl算法模型中一般假设信号矢量x(t)的先验概率分布为x(t)～cn(0,γ)，且不同时刻的信号之间统计独立，其中

为一个对角矩阵，其对角线上的元素γi,i＝1,2…n为其对应的x(t)中元素的方差。当γi≈0时，表示其所对应的入射角度上不存在目标；否则存在目标。这里需要指出的是，尽管不同时刻的信号统计独立，但是由于来自相同发射信号源，因此对于某一个入射角度上的信号，它在不同时刻的先验概率分布是相同的。综上，稀疏贝叶斯算法待估计的参数包括γ，噪声功率σ0以及目标个数k。

由于稀疏贝叶斯学习算法的每次迭代过程中都要涉及较为复杂的大规模矩阵求逆，因此稀疏贝叶斯学习算法的运算量也是制约其实际应用的关键因素。因此迫切需要一种新的稀疏贝叶斯学习算法，能够突破传统稀疏贝叶斯学习算法运算量的瓶颈。

技术实现要素：

本发明要解决的技术问题是在低信噪比或者数据欠采样条件下，目标角度估计精度与鲁棒性较低，难以满足工程实际需求，且运算效率较低，难以满足探测系统对多目标进行实时角度估计的需求。

本发明的思路是针对现有sbl算法对多目标角度估计在运算效率方面的不足，提出一种基于稀疏贝叶斯学习的快速目标角度估计方法，该方法先基于近似信息传递(amp)算法的新型sbl算法，大幅提高角度估计的运算效率；再结合最大期望法最大化当前的联合概率密度函数来优化待估计参量，通过多次迭代求出最优的待估计参数值；最后根据对角矩阵γ中非零元素的位置得到相应的目标角度估计。

本发明解决其技术问题所采取的技术方案包括以下步骤：

s1待估计参量γj，j＝1，2，…，n以及σ0的初始化

s1.1根据所需的目标角度估计分辨率需求构造字典矩阵α，α的每一列所对应的角度构成了amp算法对空间角度的网格化划分，网格划分越密则角度估计的分辨率越高。对于本发明所述方法而言，采用对全角度空间的等间隔划分，记网格所对应的角度为

s1.2在这一步骤中，将对后续所需要的参数进行初始化。这里需要初始化的参数为γ＝[γ1γ2…γn]^t以及噪声功率σ0。一个好的初始化参数值可以大大加快下面算法的收敛速度，快速得到正确的结果。由于一般应用中不存在关于目标角度的先验信息，因此初始化γ时初始成γ0i的形式，即各个方向上的信号先验方差相等。根据采样得到的来自t个不同时刻的接收数据数据y＝[y(1)y(2)…y(t)]，γ0与σ0可由下面的公式得到：

上式中，m为天线组成的阵列阵元的个数，||…||2表示矩阵的二范数，snr表示预先估计得到的系统信噪比，tr(…)表示矩阵的迹，(…)^h表示矩阵的共轭转置；

s2利用amp算法快速获得各个时刻的信号后验概率密度函数

s2.1按照s1中的初始化结果，对每一个不同时刻的接收数据y(t),t＝1,2…t分别进行如下步骤的计算(即对于每一个t,t＝1,2…t，重复进行步骤s2.1.1—s2.1.6直至对于t时刻的数据，amp算法收敛。这样的重复步骤一共需要进行t次。为了叙述的方便，在下面的步骤描述中省略掉时间标号t，例如将t时刻接收到的目标信号矢量x(t)简写为x，t时刻的阵列接收信号y(t)简写为y)；

说明：步骤s2.1.1—s2.1.6中所出现的变量以及都是中间变量，不具备实际的物理含义；而代表通过如下步骤估计得到的x的估计量；下面所说的算法迭代是指步骤s2.1.1—s2.1.6的迭代；

s2.1.1amp参数初始化：对于x的每一个元素，设置初始的估计参数值如下

这里，表示的第j个元素，表示的第j个元素的初值，xj表示x真实的第j个元素，表示对估计得到的初值，表示对概率密度函数p(xj|γj)求期望，这里p(xj|γj)表示在已知γj值的条件下xj的概率密度函数，表示对估计得到的初值，k表示第k次算法迭代，k＝0表示初始化步骤。

由于我们一般假设概率密度函数p(xj|γj)为零均值的高斯分布，因此从(5)中我们可以得到

s2.1.2线性输出步骤：对于每一个i＝1,2…m，计算

上式中表示第k次算法迭代过程中的值，aij表示字典矩阵α的第i行第j列的元素，(…)i表示矢量的第i个元素，|…|表示复数的模，表示第k次算法迭代过程中的值，表示第k次算法迭代过程中的值，表示第k次算法迭代过程中的值，表示第k次算法迭代过程中的值，表示第k次算法迭代过程中的值。

s2.1.3非线性输出步骤：对于每一个i＝1,2…m，计算

yi表示接收数据y的第i个元素，表示第k次算法迭代过程中的值，表示在第第k次算法迭代过程中更新的的值，上面的函数

s2.1.4线性输入步骤：对于每一个j＝1,2…n，计算

表示第k次算法迭代过程中的值，表示第k次算法迭代过程中的值，这里(…)-¹表示矩阵求逆，(…)*表示复数的共轭。

s2.1.5非线性输入步骤：对于每一个j＝1,2…n，计算

这里表示第k+1次迭代的值，表示第k+1次迭代的值，上面的函数

s2.1.6判断amp算法是否收敛：计算的值，其中||…||1表示矩阵的1范数，表示第k+1次迭代的值，同理，表示第k次迭代的值。如果该值大于某一设定门限ε1，则回到s2.1.2重新迭代；否则跳出循环进入s2.2得到p(xj|y)的结果。门限ε1取决于系统的信噪比等因素，需要根据实际情况进行调整。通常情况下，门限ε1的取值为0.1到0.001之间。

s2.2通过上述步骤可以得到不同时刻的信号后验概率密度函数p(xj|y)的结果如下

p(xj)表示xj的概率密度函数；上面公式中的均由s2.1中获得，为amp算法被判断收敛后，最后一次迭代过程中的的值；这里γj的值由s1或者s3获得，在第一次em算法迭代中，γj的值由s1中的初始值确定，在其他情况下，γj的值由上一次em算法循环中s3计算出来的γj确定。

s3利用em算法更新待估计参量γj,j＝1,2…n以及σ0的值

在s2中已经获得了信号的后验概率密度函数p(xj|y)，根据em算法，本步骤就是利用下面的表达式来逐个更新待估计参量的值

上式中q＝[γ1…γnσ0]^t，x＝[x(1)x(2)…x(t)]，n＝[n(1)n(2)…n(t)]，<…|y；qⁱ>表示在已知接收数据y＝[y(1)y(2)…y(t)]和给定参数值qⁱ的条件下求均值，上面的表达式中qⁱ表示在第i次算法迭代过程中的q值，而qⁱ⁺¹表示在第i+1次算法迭代过程中的q值。本步骤包括以下步骤：

s3.1由于各不同时刻的信号间统计独立，而待估计参量值相同。由于γj的更新只与有关，因此取期望时的概率密度函数可以变成p(xj(t)|y(t)；qⁱ)，即已知接收数据y(t)和给定参数值qⁱ的条件下的xj(t)的概率密度函数。通过公式推导可以得到对γj,j＝1,2…n的更新表达式为

上面的表达式中表示第i次算法迭代过程中的γj，表示第i+1次算法迭代过程中的γj。

进一步对γj求偏导可以得到

从上式中可以看出，γj,j＝1,2…n的更新不需要经过矩阵运算，而是简单的标量运算，因此可以节省大量的运算时间。

s3.2更新各个时刻的噪声功率σ0，σ0的更新首先需要计算出p(x|y；q)的表达式，这里p(x|y；q)表示在已知接收数据y和给定参数值q的条件下x的概率密度函数。结合s3.1中所计算得到的γj,j＝1,2…n，p(x|y；q)服从高斯分布，其方差σ以及均值μ可以通过下面的表达式获得

根据(17)的结果，σ0的更新可以采用如下的表达式：

这里||…||f表示矩阵的frobenius范数。

s4判断参数更新的迭代过程何时收敛，以及确定目标角度和数量

当完成了对γj,j＝1,2…n以及σ0的更新后，需要进行收敛判断。一般可以采用如下的表达式作为判断准则

这里γ＝[γ1γ2…γn]^t，γⁱ⁺¹表示第i+1次算法迭代中的γ值，γⁱ表示第i次算法迭代中的γ值；ε2是用户根据系统实际情况自定的门限值，可通过实际调试的方法进行确定。如果(19)左边的表达式值不小于ε2，则应返回s2；如果(19)左边表达式的值小于或等于ε2，则可退出循环，进行目标角度和数量的确定。通常情况下，门限ε2的取值为0.1到0.001之间。

目标角度和数量的确定可以通过设定门限ε3来进行判定得到。如果γj的值超过了门限ε3，则可判定γj所对应的来波方向存在目标；如果γj的值低于门限ε3，则不存在目标。得到了目标角度也就确定了目标的数量。通常情况下，门限ε3可以取值为除去当前被检测角度上的γj的值之后，所有γi,i＝1,2…n,i≠j的平均值的两倍。

当信噪比较低的情况下，实际目标来波方向周围的角度上对应γj的值也可能较高，为了获得准确的角度估计，则可以对γ＝[γ1γ2…γn]^t进行峰值检测，认为有且只有在γ＝[γ1γ2…γn]^t的峰值点处存在目标。

本发明取得的有益效果为：通过本发明可提升低信噪比，小样本条件下多目标角度估计精度，并且本发明估计目标角度的迭代收敛速度快，运算效率高，可应用于实时多目标角度估计系统，具有重要的工程应用价值。

附图说明

图1处理流程图；

图2新方法与传统方法的空间功率谱图；

图3空间两目标相隔较远时新方法与传统方法的性能随信噪比变化的比较；

图4空间两目标相隔较近时新方法与传统方法的性能随信噪比变化的比较；

图5空间两目标相隔较近时新方法与传统方法的性能随采集样本数变化的比较；

图6空间两目标相隔较远时新方法与传统方法的性能随采集样本数变化的比较；

图7新方法与传统方法的运算时间随不同阵元个数变化的比较图；

图8新方法与传统方法的运算时间随不同采集样本数变化的比较图。

具体实施方式

下面结合附图对本发明进行进一步说明：

图1为本发明处理流程图。

本发明所述一种基于稀疏贝叶斯学习的快速目标角度估计算法，包括以下步骤：

s1进行待估计参量γj,j＝1,2…n以及σ0的初始化，这些待估计参量的初值设置为后续em算法迭代提供基础；

s2利用amp算法快速获得各个时刻的信号后验概率密度函数，该步骤需要进行多次amp算法迭代，直至amp算法收敛；

s3利用em算法更新待估计参量γj,j＝1,2…n以及σ0的值，该步骤结合s2需要进行多次迭代，即每当通过s2计算得到一次p(xj|y),j＝1,2…n就进行一次s3计算，重复该过程直到收敛；

s4判断待估计参数更新的迭代过程是否收敛，如果不收敛，则回到s2重新迭代；如果收敛，则跳出循环并确定目标来波方向和数量。

图2为本发明的方法(下面称为amp算法)与经典的lasso算法的空间功率谱图。该仿真基于一个32阵元的均匀线阵，考虑两个不相干的目标分别从20°和30°的位置处发射信号入射到阵列上。目标信号的信噪比均为0db，阵列一共接收到10个采集样本。lasso算法的计算采用了matlab的cvx凸优化工具包。由图可知，本发明算法与lasso法相似，两者的空间谱都在目标入射的方向上出现了尖峰。因此通过极点检测的方法可以得到准确的角度估计结果。

图3为两种算法的估计精度随着信噪比变化的情况。估计精度由角度估计的均方根误差来表示。其表达式为这里表示第i次仿真得到的角度估计值。由图可知，本发明算法在信噪比低至-14db条件下仍然能够有效的得到目标角度的精确估计。并且相比于lasso算法，本发明方法在信噪比较低时能得到更为精确的目标角度估计，估计精度提高较为显著。

图4为两目标分别从25°和30°入射时，两种方法的估计精度随信噪比的变化情况。由图可知，尽管目标间隔变小，本发明方法的估计精度相较于传统的lasso算法仍然有较为显著的提升。

图5为两目标分别从25°和30°入射时，两种方法的估计精度随采集样本数变化的比较图。此时目标的信噪比为-8db。由图可知，尽管两种方法在样本数足够多的时候都能得到准确的目标角度估计，它们在采集样本数不足时的性能差异较为明显。新方法比传统方法在采集样本数不足时能得到更好的估计性能。

图6为两目标改为从20°和30°入射情况下，新方法与传统方法性能随采集样本数的变化情况，此时的信噪比仍然为-8db。此时，新方法相较于传统方法在采集样本数不足时，估计精度有明显提升。

图7为两种方法的运算时间随不同阵元个数变化的比较图。此时信噪比为-8db，两目标分别从20°和30°入射，采集样本数为50。由图可知，lasso算法的运算时间随阵元个数增加而显著增加；而本发明的新方法则变化不明显。此外，可以看出新方法大大降低了运算时间，相较于传统算法，大约降低了两个数量级。

图8为两种方法的运算时间随不同采集样本个数变化的比较图。此时信噪比为-8db，两目标分别从20°和30°入射，阵元个数为32。由图可知，两种方法的运算时间都随着采集样本数的增多而有所提升，但新方法的提升并不明显。此外，新方法相较于传统方法的运算时间降低了大约两个数量级。

基于仿真与实测数据的实验结果表明，本发明对噪声鲁棒性强，对小样本数据仍然有效，且运算效率明显高于传统方法，满足实时目标角度估计要求。本发明可在雷达回波数据质量受限条件下，实现多目标的入射角度精确估计，尤其为强对抗条件下导弹防御、空间目标监视中的空间目标识别提供了技术支持，工程应用价值高。