本发明属于阵列信号处理领域，具体涉及一种基于近似消息传递的不在网格(off-grid)目标角度估计方法。
背景技术：
在雷达、声纳和地震遥感等领域，目标角度估计(directionofarrivalestimation，doa)是一个基本问题。近些年，主流思想是引入稀疏恢复算法，根据接收信号的稀疏性特点，构建阵列接收信号的稀疏表示模型，完成信号的实时重建。相较于传统算法，其在提高对于噪声、采样数受限的信号及相关信号的鲁棒性中有着显著优势。
尽管解决这类稀疏恢复问题有困难，但是对于一个强的稀疏信号和设计良好的字典矩阵，精确恢复也是有可能的。稀疏贝叶斯学习(sparsebayesianlearning,sbl)就是其中一种经典方法。在这个方法中有两种类型的估计。第一类是基于稀疏诱导先验分布的最大后验估计。第二类估计是在隐变量空间中利用隐变量分布的变分表示进行操作的，该变分表示导致了超出模式之外的后验信息稀疏估计。针对不同应用场景提出了一系列基于这种方法的算法。这些算法或者通过期望最大化迭代或者靠证据函数的最大化迭代来求得贝叶斯分层模型中超参数的估计。但是sbl算法每次迭代都涉及到大规模的矩阵求逆，特别是在多目标求解中，该算法计算复杂度非常高，故限制了其在实际工程中的运用。为有效降低算法复杂度，在sbl的e-step中采用近似信息传递算法(approximatemessagepassingalgorithm,amp)(详情见中国专利：“一种基于稀疏贝叶斯学习的快速目标角度估计方法”，专利号：zl2019100130299)。amp算法是在底层因素图上使用中心极限定理的近似值进行环路置信度传播，通过将消息传递项引入迭代阈值方案，有效的对稀疏欠采样进行补偿，以低复杂度的运算得到信号的后验概率分布。
考虑一个由m个天线组成的阵列收到来自不同方向θ＝[θ1θ2…θk]t的k个目标的回波信号，这里[…]t表示矩阵的转置。sbl算法的信号模型可以写为
y(t)＝ax(t) n(t)(1)
其中y(t)表示t时刻的阵列接收信号，
其中d表示阵元之间的间距，λ表示阵列发射的电磁波波长，π为圆周率常数，θ表示任一入射角度。字典矩阵a中的
根据以上信号模型，我们可以得到阵列的回波信号满足分布
为一个对角矩阵，其对角线上的元素γj,j＝1,2…n为其对应的x(t)中元素的方差。当γj≈0时，表示其所对应的入射角度
由上述信号模型可知，稀疏贝叶斯算法的分辨率和角度估计精度由反映网格划分密度的字典矩阵a来决定，但是不管网格有多精细，总有处于网格之外的目标存在，这些目标的回波信号会造成字典矩阵a的失配，使系统的目标估计性能下降。因此需要一种新的稀疏贝叶斯算法来处理网格之外的目标，提高算法的角度估计精度。
技术实现要素：
本发明主要解决的技术问题是在目标不位于空间角度网格节点的情况下，传统的基于贝叶斯学习的目标角度估计算法只能认定该目标位于其临近的一个网格节点，难以对其真实角度进行精确的估计。
本发明的思路是针对现有的基于近似消息传递(amp)的贝叶斯学习算法(中国专利：“一种基于稀疏贝叶斯学习的快速目标角度估计方法”，专利号zl2019100131299)无法对离网(off-grids)目标角度进行更加精确估计的问题，提出一种基于近似消息传递的不在网格目标角度估计方法，该方法在采用amp算法对目标角度初步估计后，利用梯度衰减法更新目标点处的角度网格划分，即对字典矩阵a的部分列进行更新，使网格节点对准离网目标，避免失配情况出现，从而提高对目标角度估计的精度。
本发明采取的技术方案为：一种基于近似消息传递的不在网格目标角度估计方法，包括以下步骤：
s1待估计参量γj，j＝1，2，…，n、σ0和字典矩阵a的初始化
s1.1根据所需的目标角度估计分辨率需求构造字典矩阵α，α的每一列所对应的角度构成了amp算法对空间角度的网格化划分，网格划分越密则角度估计的分辨率越高。对于本发明所述方法而言，采用对全角度空间的等间隔划分对a进行初始化，记网格所对应的初始角度为
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中所出现的变量
s2.1.1amp参数初始化：对于x的每一个元素，设置初始的估计参数值如下
这里，
由于我们一般假设概率密度函数p(xj|γj)为零均值的高斯分布，因此从(5)中我们可以得到
s2.1.2线性输出步骤：对于每一个i＝1,2…m，计算
上式中
s2.1.3非线性输出步骤：对于每一个i＝1,2…m，计算
yi表示接收数据y的第i个元素，
s2.1.4线性输入步骤：对于每一个j＝1,2…n，计算
s2.1.5非线性输入步骤：对于每一个j＝1,2…n，计算
这里
s2.1.6判断amp算法是否收敛：计算
s2.2通过上述步骤可以得到t时刻的信号后验概率密度函数p(xj|y)的结果如下
p(xj)表示xj的概率密度函数；上面公式中的
s3利用em算法更新待估计参量γj,j＝1,2…n，求出目标数量最优解k*，更新噪声功率σ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；qi＞表示在已知接收数据y＝[y(1)y(2)…y(t)]和给定参数值qi的条件下求均值，上面的表达式中qi表示在第i次算法迭代过程中的q值，而qi 1表示在第i 1次算法迭代过程中的q值。本步骤包括以下步骤：
s3.1更新γj,j＝1,2…n：
由于各不同时刻的信号间统计独立，而待估计参量值相同。由于γj的更新只与
上面的表达式中
进一步对γj求偏导可以得到
从上式中可以看出，γj,j＝1,2…n的更新不需要经过矩阵运算，而是简单的标量运算，因此可以节省大量的运算时间。
s3.2估算目标数量最优解k*：不同于中国发明专利“一种基于稀疏贝叶斯学习的快速目标角度估计方法”，专利号zl2019100131299，本发明在估算噪声功率之前需要先确定目标的数量，便于在s3.3更新噪声功率σ0时可以利用本步骤计算的中间结果，降低计算资源消耗。
此步骤参考文献1、z.-m.liu,z.-t.huang,y.-y.zhou,anefficientmaximumlikelihoodmethodfordirection-of-arrivalestimationviasparsebayesianlearning,ieeetrans.wirelesscommunications11(10)(2012)3607–3617；2、p.stoica,y.selen,model-orderselection:areviewofinformationcriterionrules,ieeesignalprocess.mag.21(4)(2004)36–47；3、c.d.austin,r.l.moses,j.n.ash,e.ertin,ontherelationbetweensparsereconstructionandparameterestimationwithmodelorderselection,ieeej.sel.topicssignalprocess.4(3)(2010)560–570。
由上述参考文献可知，采用子空间分析法，目标数量最优解为
其中i是单位矩阵，
这里k是目标数量的估计值，k*为k的最优解，
s3.3更新噪声功率σ0：
此步骤参考文献p.gerstoft,c.f.mecklenbrauker,a.xenaki,s.nannuru,multisnapshotsparsebayesianlearningfordoa,ieeesignalprocess.letters23(10)(2016)1469–1473.
由上述参考文献可得
其中
根据s3.2定义的映射矩阵p，且p＝ph＝p2，可得
prph-σ0pph＝σy-σ0i(21)
结合tr(r-σy)≈0可以推出σ0的更新公式为：
对比公式(17)，可以看出求解k*和σ0可以同时进行，节省计算资源。
s4用梯度下降法更新字典矩阵a中的
通过上述步骤，已经粗略估算出目标的数量，及其所处的网格节点(即超过门限ε5的γj所对应的网格节点，也即目标的角度估计值)，但由于目标的真实位置可能处于两个网格节点之间，故本步骤采用梯度下降算法对网格节点进行调整，使网格节点更加接近目标真实位置，从而提高目标角度估计精度。由s3.2中计算出的目标数量最优解k*，可认为目标对应角度为γj,j＝1,2…n中前k*个最大的γj所对应的角度网格节点
其中
这里
通过(24)式对
其中，α是角度更新的步长，实际应用中步长取决于对角度估计的精度要求。(·)l表示本步骤内部迭代第l次的值。如果函数sign(·)的输入是正数则等于1，反之为-1。
当完成一次角度更新后，需要判断是否收敛。考虑到初始的字典网格可能已经满足精度需求，所以此步骤中迭代次数不会太多，当满足条件
s5判断算法迭代过程何时收敛，以及确定来波方向和数量
完成对γ＝[γ1γ2...γn]t以及σ0的值更新后，用如下表达式判断收敛：
ε2为判定门限，可根据系统实际进行设定，通常情况下取值为0.1到0.001之间。如果上式不成立，则返回s2继续进行迭代运算；如果上式成立，则可退出循环，来波数量即为最后一次计算出的k*，方向为γ＝[γ1γ2…γn]t中前k*个最大的γj所对应的方向
本发明取得的有益效果为：本发明可以较小的计算复杂度增加为代价，对位于角度网格节点以外的目标，自动调整网格划分，使之处于网格节点上，能够更加精确的估计目标的角度，运算效率较高，可应用于实时多目标高精度角度估计系统，具有重要的工程应用价值。
附图说明
图1处理流程图；
图2阵元数目为16时本发明方法与传统方法的空间功率谱图；
图3阵元数目为16时本发明方法与传统方法的性能随信噪比变化的比较；
图4阵元数目为16时本发明方法与传统方法的性能随采集样本数变化的比较；
图5阵元数目为16时本发明方法与传统方法的运算时间随样本数变化的比较图；
图6阵元数目为10时本发明方法与传统方法的运算时间随样本数变化的比较图；
图7阵元数目为16，角度间隔为0.5°时，本发明方法与传统方法的运算时间随样本数变化的比较图。
具体实施方式
下面结合附图对本发明进行进一步说明：
图1为本发明总处理流程。
本发明所述一种基于广义近似消息传递的贝叶斯学习快速目标角度估计算法，包括以下步骤：
s1待估计参量γj，j＝1，2，…，n、σ0和字典矩阵a的初始化；
s2利用amp算法快速获得各个时刻的信号后验概率密度函数；
s3利用em算法更新待估计参量γj,j＝1,2…n，求出目标数量最优解k*，更新噪声功率σ0；
s4用梯度下降法更新字典矩阵a中的
s5判断算法迭代过程何时收敛，以及确定来波方向和数量。
图2为本发明的方法与经典的lasso、rvm和atomicnorm算法的空间功率谱图。该仿真基于阵元数为16，间隔为半波长的均匀线阵。考虑两个不相干的目标分别于-5°和15.5°的位置处发射信号入射到阵列上，目标信号的信噪比均为-10db，阵列一共接收到50个采集样本。四种算法均采用相同的空间网格划分和字典矩阵，即对范围为-45°到45°的角度空间进行间隔1°的划分。由图可知，四种算法的空间普均在目标位置出现峰值，lasso算法峰值最窄，atomicnorm算法次之，rvm算法峰值最宽，本发明算法峰值宽度处于atomicnorm、rvm这两种算法中间。因此通过极点检测的方法可以得到准确的角度估计结果。
图3为atomicnorm算法、本算法以及不含s4步骤的本算法的估计精度随着信噪比变化的情况。估计精度由50次仿真的角度估计值的均方根误差来表示，其表达式为
图4为四种方法的估计精度随采集样本数的变化情况，此时的信噪比为-10db，(a)为目标角度为5°、30°，(b)为目标角度为-5°、15°。由图可知四种方法随样本数增加rmse都呈现下降趋势。rvm算法在小样本情况下估计精度最好，但随样本数增多rmse下降缓慢，本发明方法的表现与lasso算法相近。
图5为运算时间随样本数变化的比较图。通过运算时间来衡量各算法的计算复杂度，此时阵元数为16，snr为-5db，网格间距为1°，两目标入射角分别为5°和30.2°。由图可知，随着样本数目增加，运算时间都呈增加趋势，但是本发明方法的运算时间比lasso、rvm和atomicnorm算法要第几个量级。对比图3可知，加入梯度下降算法只是稍微增加了运算时间，但是doa估计精度却显著增加。
图6为在图5基础上讲阵元数设为10，所得到的运算时间随样本数变化的比较图。对比图5可以看出随着阵元数减少，各算法计算时间均有明显下降，本发明算法任然优于其它算法。
图7为在图5基础上将角度划分间隔变为0.5°，所得到的运算时间随样本数变化的比较图。对比图5可以看出，网格精度对运算时间有显著影响，各算法的运算时间都大幅增加，但是本发明算法运算时间依旧显著低于其它算法。
基于仿真的实验结果表明，本发明算法对噪声鲁棒性强，对小样本数据仍然有效，且运算效率和角度估计精度明显高于传统方法，满足实时目标角度估计要求。本发明可在雷达回波数据质量受限条件下，实现多目标的入射角度精确估计，尤其为强对抗条件下导弹防御、空间目标监视中的空间目标识别提供了技术支持，工程应用价值高。
1.一种基于近似消息传递的不在网格目标角度估计方法，其特征在于，该方法包括以下步骤：
s1待估计参量γj，j＝1，2，…，n、σ0和字典矩阵a的初始化
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次；
s2.1.1amp参数初始化：对于x的每一个元素，设置初始的估计参数值如下
这里，
假设概率密度函数p(xj|γj)为零均值的高斯分布，因此从(2)中我们可以得到
s2.1.2线性输出步骤：对于每一个i＝1,2…m，计算
上式中
s2.1.3非线性输出步骤：对于每一个i＝1,2…m，计算
yi表示接收数据y的第i个元素，
s2.1.4线性输入步骤：对于每一个j＝1,2…n，计算
s2.1.5非线性输入步骤：对于每一个j＝1,2…n，计算
这里
s2.1.6判断amp算法是否收敛：计算
s2.2通过上述步骤可以得到t时刻的信号后验概率密度函数p(xj|y)的结果如下
p(xj)表示xj的概率密度函数；上面公式中的
s3利用em算法更新待估计参量γj,j＝1,2…n，求出目标数量最优解k*，更新噪声功率σ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；qi＞表示在已知接收数据y＝[y(1)y(2)…y(t)]和给定参数值qi的条件下求均值，上面的表达式中qi表示在第i次算法迭代过程中的q值，而qi 1表示在第i 1次算法迭代过程中的q值；本步骤包括以下步骤：
s3.1更新γj,j＝1,2…n：
通过公式推导可以得到对γj,j＝1,2…n的更新表达式为
上面的表达式中
进一步对γj求偏导可以得到
s3.2估算目标数量最优解k*：
采用子空间分析法，目标数量最优解为
其中i是单位矩阵，
k是目标数量的估计值，k*为k的最优解，
s3.3更新噪声功率σ0：
根据公式：
其中
根据s3.2定义的映射矩阵p，且p＝ph＝p2，可得
prph-σ0pph＝σy-σ0i(18)
结合tr(r-σy)≈0可以推出σ0的更新公式为：
s4用梯度下降法更新字典矩阵a中的
由s3.2中计算出的目标数量最优解k*，可认为目标对应角度为γj,j＝1,2…n中前k*个最大的γj所对应的角度网格节点
其中
这里
通过(21)式对
其中，α是角度更新的步长，实际应用中步长取决于对角度估计的精度要求，(·)l表示本步骤内部迭代第l次的值，如果函数sign(·)的输入是正数则等于1，反之为-1；
当完成一次角度更新后，需要判断是否收敛；考虑到初始的字典网格可能已经满足精度需求，所以此步骤中迭代次数不会太多，当满足条件
s5判断算法迭代过程何时收敛，以及确定来波方向和数量
完成对γ＝[γ1γ2...γn]t以及σ0的值更新后，用如下表达式判断收敛：
ε2为判定门限，可根据系统实际进行设定；如果上式不成立，则返回s2继续进行迭代运算；如果上式成立，则可退出循环，来波数量即为最后一次计算出的k*，方向为γ＝[γ1γ2…γn]t中前k*个最大的γj所对应的方向
2.根据权利要求1所述基于稀疏贝叶斯学习的快速目标角度估计方法，其特征在于：s1.1中采用对全角度空间的等间隔划分对空间角度进行网格化划分。
3.一种根据权利要求1或2所述基于近似消息传递的不在网格目标角度估计方法，其特征在于：s2.1.6中，门限ε1的取值为0.1到0.001之间。
4.一种根据权利要求1或2所述基于近似消息传递的不在网格目标角度估计方法，其特征在于：s4中，判决阈值ε3取值为0.1°。
5.一种根据权利要求1或2所述基于近似消息传递的不在网格目标角度估计方法，其特征在于：s4中，判决阈值ε4取值为1000。
6.一种根据权利要求1或2所述基于近似消息传递的不在网格目标角度估计方法，其特征在于：s5中，判定门限ε2取值为0.1到0.001之间。
技术总结