CN111257845A - 一种基于近似消息传递的不在网格目标角度估计方法 - Google Patents

一种基于近似消息传递的不在网格目标角度估计方法 Download PDF

Info

Publication number
CN111257845A
CN111257845A CN202010086824.0A CN202010086824A CN111257845A CN 111257845 A CN111257845 A CN 111257845A CN 202010086824 A CN202010086824 A CN 202010086824A CN 111257845 A CN111257845 A CN 111257845A
Authority
CN
China
Prior art keywords
algorithm
value
angle
representing
target
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010086824.0A
Other languages
English (en)
Other versions
CN111257845B (zh
Inventor
张新禹
张俊
霍凯
张双辉
刘永祥
姜卫东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202010086824.0A priority Critical patent/CN111257845B/zh
Publication of CN111257845A publication Critical patent/CN111257845A/zh
Application granted granted Critical
Publication of CN111257845B publication Critical patent/CN111257845B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G01S7/418—Theoretical aspects

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明属于阵列信号处理领域，具体涉及一种基于近似消息传递的不在网格目标角度估计方法，包括以下步骤：S1待估计参量γj，j＝1，2，…，N、σ0和字典矩阵A的初始化；S2利用AMP算法快速获得各个时刻的信号后验概率密度函数；S3利用EM算法更新待估计参量γj,j＝1,2…N，求出目标数量最优解K*，更新噪声功率σ0；S4用梯度下降法更新字典矩阵A中的
Figure DDA0002382349660000011
S5判断算法迭代过程何时收敛，以及确定来波方向和数量。本发明取得的有益效果为：本发明可以较小的计算复杂度增加为代价，对位于角度网格节点以外的目标，自动调整网格划分，使之处于网格节点上，能够更加精确的估计目标的角度，运算效率较高，可应用于实时多目标高精度角度估计***，具有重要的工程应用价值。

Description

一种基于近似消息传递的不在网格目标角度估计方法
技术领域
本发明属于阵列信号处理领域，具体涉及一种基于近似消息传递的不在网格(off-grid) 目标角度估计方法。
背景技术
在雷达、声纳和地震遥感等领域，目标角度估计(Direction of ArrivalEstimation，DOA) 是一个基本问题。近些年，主流思想是引入稀疏恢复算法，根据接收信号的稀疏性特点，构 建阵列接收信号的稀疏表示模型，完成信号的实时重建。相较于传统算法，其在提高对于噪 声、采样数受限的信号及相关信号的鲁棒性中有着显著优势。
尽管解决这类稀疏恢复问题有困难，但是对于一个强的稀疏信号和设计良好的字典矩阵， 精确恢复也是有可能的。稀疏贝叶斯学习(Sparse Bayesian Learning,SBL)就是其中一种经典 方法。在这个方法中有两种类型的估计。第一类是基于稀疏诱导先验分布的最大后验估计。 第二类估计是在隐变量空间中利用隐变量分布的变分表示进行操作的，该变分表示导致了超 出模式之外的后验信息稀疏估计。针对不同应用场景提出了一系列基于这种方法的算法。这 些算法或者通过期望最大化迭代或者靠证据函数的最大化迭代来求得贝叶斯分层模型中超参 数的估计。但是SBL算法每次迭代都涉及到大规模的矩阵求逆，特别是在多目标求解中，该 算法计算复杂度非常高，故限制了其在实际工程中的运用。为有效降低算法复杂度，在SBL 的E-Step中采用近似信息传递算法(Approximate MessagePassing algorithm,AMP)(详情见 中国专利：“一种基于稀疏贝叶斯学习的快速目标角度估计方法”，专利号：ZL2019100130299)。 AMP算法是在底层因素图上使用中心极限定理的近似值进行环路置信度传播，通过将消息传 递项引入迭代阈值方案，有效的对稀疏欠采样进行补偿，以低复杂度的运算得到信号的后验 概率分布。
考虑一个由m个天线组成的阵列收到来自不同方向θ＝[θ1 θ2 … θK]T的K个目标的 回波信号，这里[…]T表示矩阵的转置。SBL算法的信号模型可以写为
y(t)＝Ax(t)+n(t) (1)
其中y(t)表示t时刻的阵列接收信号，
Figure BDA0002382349640000011
表示字典矩阵，
Figure BDA0002382349640000012
表示阵列的导向矢量，阵列的导向矢量表示远场条件下的平面波入射到阵 列所形成的接收信号矢量，它是一个复向量，维数等于阵列的阵元个数，对于一个由m个天 线组成的均匀线阵，导向矢量可以表示成：
Figure BDA0002382349640000013
其中d表示阵元之间的间距，λ表示阵列发射的电磁波波长，π为圆周率常数，θ表示任 一入射角度。字典矩阵A中的
Figure RE-GDA0002433713180000014
组成了一个关于空间入射角的一种网格划分， 该网格划分越细，目标角度估计的分辨率往往就越高，一般的做法是将一个可能的空域区间 做等间隔的划分，N的值会远大于阵列阵元个数m。(1)中的x(t)表示接收到的目标信号矢 量，它的元素与字典矩阵中的每一列一一对应，当其所对应的入射角度存在目标时，x(t)的 该元素等于信号的复数值；当其所对应的角度不存在目标时，x(t)的该元素值为0。n(t)表示 ***的加性噪声。这里我们假设在不同时刻的目标信号统计独立，这一假设在雷达信号处理 中十分常见。在实际工程应用中，n(t)一般假设为零均值的高斯白噪声，即
Figure RE-GDA0002433713180000021
并且不同时刻间的噪声统计独立，这里
Figure RE-GDA0002433713180000022
表示复数零均值的复高斯分布，其方差为 σ0I，σ0表示噪声的功率，I表示单位矩阵；下文中用
Figure RE-GDA0002433713180000023
表示均值为α，方差为β的 复高斯分布函数。
根据以上信号模型，我们可以得到阵列的回波信号满足分布
Figure BDA0002382349640000024
为了 进行贝叶斯推导，SBL算法模型中一般假设信号矢量x(t)的先验概率分布为
Figure BDA0002382349640000025
且不同时刻的信号之间统计独立，其中
Figure BDA0002382349640000026
为一个对角矩阵，其对角线上的元素γj,j＝1,2…N为其对应的x(t)中元素的方差。当 γj≈0时，表示其所对应的入射角度
Figure BDA0002382349640000027
上不存在目标；否则存在目标。这里需要指出的是，尽 管不同时刻的信号统计独立，但是由于来自相同发射信号源，因此对于某一个入射角度上的 信号，它在不同时刻的先验概率分布是相同的。综上，稀疏贝叶斯算法待估计的参数包括Γ， 噪声功率σ0以及目标个数K。
由上述信号模型可知，稀疏贝叶斯算法的分辨率和角度估计精度由反映网格划分密度的 字典矩阵A来决定，但是不管网格有多精细，总有处于网格之外的目标存在，这些目标的回 波信号会造成字典矩阵A的失配，使***的目标估计性能下降。因此需要一种新的稀疏贝叶 斯算法来处理网格之外的目标，提高算法的角度估计精度。
发明内容
本发明主要解决的技术问题是在目标不位于空间角度网格节点的情况下，传统的基于贝 叶斯学习的目标角度估计算法只能认定该目标位于其临近的一个网格节点，难以对其真实角 度进行精确的估计。
本发明的思路是针对现有的基于近似消息传递(AMP)的贝叶斯学习算法(中国专利： “一种基于稀疏贝叶斯学习的快速目标角度估计方法”，专利号ZL2019100131299)无法对离 网(off-grids)目标角度进行更加精确估计的问题，提出一种基于近似消息传递的不在网格目 标角度估计方法，该方法在采用AMP算法对目标角度初步估计后，利用梯度衰减法更新目标 点处的角度网格划分，即对字典矩阵A的部分列进行更新，使网格节点对准离网目标，避免 失配情况出现，从而提高对目标角度估计的精度。
本发明采取的技术方案为：一种基于近似消息传递的不在网格目标角度估计方法，包括 以下步骤：
S1待估计参量γj，j＝1，2，…，N、σ0和字典矩阵A的初始化
S1.1根据所需的目标角度估计分辨率需求构造字典矩阵Α，Α的每一列所对应的角度构 成了AMP算法对空间角度的网格化划分，网格划分越密则角度估计的分辨率越高。对于本发 明所述方法而言，采用对全角度空间的等间隔划分对A进行初始化，记网格所对应的初始角 度为
Figure BDA0002382349640000031
S1.2在这一步骤中，将对后续所需要的参数进行初始化。这里需要初始化的参数为 γ＝[γ1 γ2 … γN]T以及噪声功率σ0。一个好的初始化参数值可以大大加快下面算法的收敛 速度，快速得到正确的结果。由于一般应用中不存在关于目标角度的先验信息，因此初始化 Γ时初始化成γ0I的形式，即各个方向上的信号先验方差相等。根据采样得到的来自T个不同 时刻的接收数据Y＝[y(1) y(2) … y(T)]，γ0与σ0可由下面的公式得到：
Figure BDA0002382349640000032
上式中，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中所出现的变量
Figure BDA0002382349640000033
以及
Figure BDA0002382349640000034
都是中间变量，不具备实际的物理含义；而
Figure BDA0002382349640000035
代表通过如下步骤估计得到 的x的估计量；下面所说的算法迭代是指步骤S2.1.1—S2.1.6的迭代；
S2.1.1AMP参数初始化：对于x的每一个元素，设置初始的估计参数值如下
Figure BDA0002382349640000036
这里，
Figure BDA0002382349640000037
表示
Figure BDA0002382349640000038
的第j个元素，
Figure BDA0002382349640000039
表示
Figure BDA00023823496400000310
的第j个元素的初值，xj表示x真实的第j个 元素，
Figure BDA00023823496400000311
表示对
Figure BDA00023823496400000312
估计得到的初值，
Figure BDA00023823496400000313
表示对概率密度函数p(xj|γj)求期望，这 里p(xj|γj)表示在已知γj值的条件下xj的概率密度函数，
Figure BDA00023823496400000314
表示对
Figure BDA00023823496400000315
估计得到的初值，k 表示第k次算法迭代，k＝0表示初始化步骤。
由于我们一般假设概率密度函数p(xj|γj)为零均值的高斯分布，因此从(5)中我们可以得 到
Figure BDA0002382349640000041
S2.1.2线性输出步骤：对于每一个i＝1,2…m，计算
Figure BDA0002382349640000042
上式中
Figure BDA0002382349640000043
表示第k次算法迭代过程中的
Figure BDA0002382349640000044
值，aij表示字典矩阵Α的第i行第j列的元 素，(…)i表示矢量的第i个元素，|…|表示复数的模，
Figure BDA0002382349640000045
表示第k次算法迭代过程中的
Figure BDA0002382349640000046
值，
Figure BDA0002382349640000047
表示第k次算法迭代过程中的
Figure BDA0002382349640000048
值，
Figure BDA0002382349640000049
表示第k次算法迭代过程中的
Figure BDA00023823496400000410
值，
Figure BDA00023823496400000411
表 示第k次算法迭代过程中的
Figure BDA00023823496400000412
值，
Figure BDA00023823496400000413
表示第k次算法迭代过程中的
Figure BDA00023823496400000414
值。
S2.1.3非线性输出步骤：对于每一个i＝1,2…m，计算
Figure BDA00023823496400000415
yi表示接收数据y的第i个元素，
Figure BDA00023823496400000416
表示第k次算法迭代过程中的
Figure BDA00023823496400000417
值，
Figure BDA00023823496400000418
表示 在第第k次算法迭代过程中更新的
Figure BDA00023823496400000419
的值，上面的函数
Figure BDA00023823496400000420
S2.1.4线性输入步骤：对于每一个j＝1,2…N，计算
Figure BDA00023823496400000421
Figure BDA00023823496400000422
表示第k次算法迭代过程中的
Figure BDA00023823496400000423
值，
Figure BDA00023823496400000424
表示第k次算法迭代过程中的
Figure BDA00023823496400000425
值，这里 (…)-1表示矩阵求逆，(…)*表示复数的共轭。
S2.1.5非线性输入步骤：对于每一个j＝1,2…N，计算
Figure BDA00023823496400000426
这里
Figure BDA00023823496400000427
表示第k+1次迭代的
Figure BDA00023823496400000428
值，
Figure BDA00023823496400000429
表示第k+1次迭代的
Figure BDA00023823496400000430
值，上面的函数
Figure BDA0002382349640000051
S2.1.6判断AMP算法是否收敛：计算
Figure BDA0002382349640000052
的值，其中||…||1表示矩阵的1范数，
Figure BDA0002382349640000053
表示第k+1次迭代的
Figure BDA0002382349640000054
值，同理，
Figure BDA0002382349640000055
表示第k次 迭代的
Figure BDA0002382349640000056
值。如果该值大于某一设定门限ε1，则回到S2.1.2重新迭代；否则跳出循环进入S2.2 得到p(xj|y)的结果。门限ε1取决于***的信噪比等因素，需要根据实际情况进行调整。通常 情况下，门限ε1的取值为0.1到0.001之间。
S2.2通过上述步骤可以得到t时刻的信号后验概率密度函数p(xj|y)的结果如下
Figure BDA0002382349640000057
p(xj)表示xj的概率密度函数；上面公式中的
Figure BDA0002382349640000058
为S2.1步骤中根据t时刻的样本 值y(t)进行计算获得，为AMP算法被判断收敛后，最后一次迭代过程中的
Figure BDA0002382349640000059
的值；这里 γj的值由S1或者S3获得，在第一次EM算法迭代中，γj的值由S1中的初始值确定，在其他情况下，γj的值由上一次EM算法循环中S3计算出来的γj确定。
S3利用EM算法更新待估计参量γj,j＝1,2…N，求出目标数量最优解K*，更新噪声功 率σ0
在S2中已经获得了信号的后验概率密度函数p(xj|y)，根据EM算法，本步骤利用下面 的表达式来逐个更新待估计参量的值
Figure BDA00023823496400000511
上式中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的更新只与
Figure BDA00023823496400000510
有关，因此取期望时的概率密度函数可以变成p(xj(t)|y(t)；qi)，即已知接收数据y(t)和给定参 数值qi的条件下的xj(t)的概率密度函数。通过公式推导可以得到对γj,j＝1,2…N的更新表达 式为
Figure BDA0002382349640000061
上面的表达式中
Figure BDA0002382349640000062
表示第i次算法迭代过程中的γj，
Figure BDA0002382349640000063
表示第i+1次算法迭代过程中的 γj。
进一步对γj求偏导可以得到
Figure BDA0002382349640000064
从上式中可以看出，γj,j＝1,2…N的更新不需要经过矩阵运算，而是简单的标量运算， 因此可以节省大量的运算时间。
S3.2估算目标数量最优解K*：不同于中国发明专利“一种基于稀疏贝叶斯学习的快速目 标角度估计方法”，专利号ZL2019100131299，本发明在估算噪声功率之前需要先确定目标的 数量，便于在S3.3更新噪声功率σ0时可以利用本步骤计算的中间结果，降低计算资源消耗。
此步骤参考文献1、Z.-M.Liu,Z.-T.Huang,Y.-Y.Zhou,An efficient maximumlikelihood method for direction-of-arrival estimation via sparse bayesianlearning,IEEE Trans.Wireless Communications 11(10)(2012)3607–3617；2、P.Stoica,Y.Selen,Model-order selection:a review of information criterion rules,IEEESignal Process.Mag.21(4)(2004)36–47；3、C.D.Austin,R. L.Moses,J.N.Ash,E.Ertin,On the relation between sparse reconstruction and parameter estimation withmodel order selection,IEEE J.Sel.Topics Signal Process.4(3)(2010)560–570。
由上述参考文献可知，采用子空间分析法，目标数量最优解为
Figure BDA0002382349640000065
其中I是单位矩阵，
Figure BDA0002382349640000066
Figure BDA0002382349640000067
这里K是目标数量的估计值，K*为K的最优解，
Figure BDA0002382349640000068
由A中对应假定目标角度
Figure BDA0002382349640000069
的列组 成的矩阵，假定目标角度
Figure BDA00023823496400000610
指的是γj,j＝1,2…N中最大的K个值所对应的角度。
S3.3更新噪声功率σ0：
此步骤参考文献P.Gerstoft,C.F.Mecklenbrauker,A.Xenaki,S.Nannuru,Multisnapshot sparse bayesian learning for doa,IEEE Signal Process.Letters 23(10)(2016)1469–1473.
由上述参考文献可得
Figure BDA00023823496400000611
其中
Figure BDA0002382349640000071
这里
Figure BDA0002382349640000072
是
Figure BDA0002382349640000073
中各列的协方差矩阵，从上式可得
Figure BDA0002382349640000074
根据S3.2定义的映射矩阵P，且P＝PH＝P2，可得
PRPH-σ0PPH＝ΣY-σ0I (21)
结合tr(R-ΣY)≈0可以推出σ0的更新公式为：
Figure BDA0002382349640000075
对比公式(17)，可以看出求解K*和σ0可以同时进行，节省计算资源。
S4用梯度下降法更新字典矩阵A中的
Figure BDA0002382349640000076
通过上述步骤，已经粗略估算出目标的数量，及其所处的网格节点(即超过门限ε5的γj所对应的网格节点，也即目标的角度估计值)，但由于目标的真实位置可能处于两个网格节点 之间，故本步骤采用梯度下降算法对网格节点进行调整，使网格节点更加接近目标真实位置， 从而提高目标角度估计精度。由S3.2中计算出的目标数量最优解K*，可认为目标对应角度 为γj,j＝1,2…N中前K*个最大的γj所对应的角度网格节点
Figure BDA0002382349640000077
使用梯度下降法对这些目标 粗略对应的角度网格节点进行调整，需要注意的是这里只需要更新目标当前对应的角度网格 节点
Figure BDA0002382349640000078
不必更新整个角度网格，以免提高运算复杂度(即对于S3.2中确定的每一个目标所 对应的角度
Figure BDA0002382349640000079
进行本步骤迭代计算，直到算法收敛。S4步骤一共需要执行K*次)。可由下式 推导更新规则
Figure BDA00023823496400000710
其中
Figure BDA00023823496400000711
A是带角度变量的字典矩阵。根据上式可求出关于
Figure BDA00023823496400000712
的梯 度
Figure BDA00023823496400000713
这里
Figure BDA00023823496400000714
表示矩阵的实部，X是由S2步骤求出的对X的估计量，
Figure BDA00023823496400000715
是A中对应于除
Figure BDA00023823496400000716
以外角度
Figure BDA00023823496400000717
的列所组成的矩阵，
Figure BDA00023823496400000718
表示X中对应于除
Figure BDA00023823496400000719
以外角度
Figure BDA00023823496400000720
的行所组成的矩阵，
Figure BDA00023823496400000721
表示X中对应于角度
Figure BDA00023823496400000722
的行矩阵，
Figure BDA00023823496400000723
是A中的列向量，d(·)是a(·)的导数。
通过(24)式对
Figure BDA00023823496400000724
进行优化，有
Figure BDA00023823496400000725
其中，α是角度更新的步长，实际应用中步长取决于对角度估计的精度要求。(·)l表示 本步骤内部迭代第l次的值。如果函数sign(·)的输入是正数则等于1，反之为-1。
当完成一次角度更新后，需要判断是否收敛。考虑到初始的字典网格可能已经满足精度 需求，所以此步骤中迭代次数不会太多，当满足条件
Figure BDA0002382349640000081
或者l≥ε4时，认为对
Figure BDA0002382349640000082
的估计已经满足要求，退出循环，ε3、ε4为判决阈值。其中，ε3可以根据***对角度估计的 精度要求确定，一般取值为0.1°；ε4根据***的实时性要求确定，值越大本算法的实时性就 越低，通常情况可取值为1000。
S5判断算法迭代过程何时收敛，以及确定来波方向和数量
完成对γ＝[γ1 γ2 ... γN]T以及σ0的值更新后，用如下表达式判断收敛：
Figure BDA0002382349640000083
ε2为判定门限，可根据***实际进行设定，通常情况下取值为0.1到0.001之间。如果上 式不成立，则返回S2继续进行迭代运算；如果上式成立，则可退出循环，来波数量即为最后 一次计算出的K*，方向为γ＝[γ1 γ2 … γN]T中前K*个最大的γj所对应的方向
Figure BDA0002382349640000084
由于在 S4中对
Figure BDA0002382349640000085
进行了更精确的计算，所以本算法对离网目标的角度估计精度大大提高。
本发明取得的有益效果为：本发明可以较小的计算复杂度增加为代价，对位于角度网格 节点以外的目标，自动调整网格划分，使之处于网格节点上，能够更加精确的估计目标的角 度，运算效率较高，可应用于实时多目标高精度角度估计***，具有重要的工程应用价值。
附图说明
图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中的
Figure BDA0002382349640000091
S5判断算法迭代过程何时收敛，以及确定来波方向和数量。
图2为本发明的方法与经典的LASSO、RVM和Atomic Norm算法的空间功率谱图。该仿真基于阵元数为16，间隔为半波长的均匀线阵。考虑两个不相干的目标分别于-5°和15.5° 的位置处发射信号入射到阵列上，目标信号的信噪比均为-10dB，阵列一共接收到50个采集 样本。四种算法均采用相同的空间网格划分和字典矩阵，即对范围为-45°到45°的角度空间 进行间隔1°的划分。由图可知，四种算法的空间普均在目标位置出现峰值，LASSO算法峰 值最窄，Atomic Norm算法次之，RVM算法峰值最宽，本发明算法峰值宽度处于Atomic Norm、 RVM这两种算法中间。因此通过极点检测的方法可以得到准确的角度估计结果。
图3为Atomic Norm算法、本算法以及不含S4步骤的本算法的估计精度随着信噪比变 化的情况。估计精度由50次仿真的角度估计值的均方根误差来表示，其表达式为
Figure BDA0002382349640000092
这里
Figure BDA0002382349640000093
表示第i次仿真得到的角度估计值。由图可知，本发明算法在两组目标(角度分别为-5°、15.5°和5°、30.2°)测试中性能优于Atomic Norm算法和 不含S4步骤的算法，尤其在信噪比较低的情况下性能更好。
图4为四种方法的估计精度随采集样本数的变化情况，此时的信噪比为-10dB，(a)为目 标角度为5°、30°，(b)为目标角度为-5°、15°。由图可知四种方法随样本数增加RMSE都呈现下降趋势。RVM算法在小样本情况下估计精度最好，但随样本数增多RMSE下降缓 慢，本发明方法的表现与LASSO算法相近。
图5为运算时间随样本数变化的比较图。通过运算时间来衡量各算法的计算复杂度，此 时阵元数为16，SNR为-5dB，网格间距为1°，两目标入射角分别为5°和30.2°。由图可知，随着样本数目增加，运算时间都呈增加趋势，但是本发明方法的运算时间比LASSO、RVM和Atomic Norm算法要第几个量级。对比图3可知，加入梯度下降算法只是稍微增加了运算时间，但是DOA估计精度却显著增加。
图6为在图5基础上讲阵元数设为10，所得到的运算时间随样本数变化的比较图。对比 图5可以看出随着阵元数减少，各算法计算时间均有明显下降，本发明算法任然优于其它算 法。
图7为在图5基础上将角度划分间隔变为0.5°，所得到的运算时间随样本数变化的比较 图。对比图5可以看出，网格精度对运算时间有显著影响，各算法的运算时间都大幅增加， 但是本发明算法运算时间依旧显著低于其它算法。
基于仿真的实验结果表明，本发明算法对噪声鲁棒性强，对小样本数据仍然有效，且运 算效率和角度估计精度明显高于传统方法，满足实时目标角度估计要求。本发明可在雷达回 波数据质量受限条件下，实现多目标的入射角度精确估计，尤其为强对抗条件下导弹防御、 空间目标监视中的空间目标识别提供了技术支持，工程应用价值高。

Claims (6)

1.一种基于近似消息传递的不在网格目标角度估计方法，其特征在于，该方法包括以下步骤：
S1待估计参量γj，j＝1，2，…，N、σ0和字典矩阵A的初始化
S1.1根据所需的目标角度估计分辨率需求构造字典矩阵Α，Α的每一列所对应的角度构成了AMP算法对空间角度的网格化划分，网格划分越密则角度估计的分辨率越高，记网格所对应的初始角度为
Figure FDA0002382349630000011
S1.2对γ＝[γ1 γ2…γN]T以及噪声功率σ0进行初始化；初始化Γ时初始化成γ0I的形式，即各个方向上的信号先验方差相等，根据采样得到的来自T个不同时刻的接收数据Y＝[y(1) y(2)…y(T)]，γ0与σ0可由下面的公式得到：
Figure FDA0002382349630000012
上式中，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的每一个元素，设置初始的估计参数值如下
Figure FDA0002382349630000013
这里，
Figure FDA0002382349630000014
表示
Figure FDA0002382349630000015
的第j个元素，
Figure FDA0002382349630000016
表示
Figure FDA0002382349630000017
的第j个元素的初值，xj表示x真实的第j个元素，
Figure FDA0002382349630000018
表示对
Figure FDA0002382349630000019
估计得到的初值，
Figure FDA00023823496300000110
表示对概率密度函数p(xj|γj)求期望，p(xj|γj)表示在已知γj值的条件下xj的概率密度函数，
Figure FDA00023823496300000111
表示对
Figure FDA00023823496300000112
估计得到的初值，k表示第k次算法迭代，k＝0表示初始化步骤；
假设概率密度函数p(xj|γj)为零均值的高斯分布，因此从(2)中我们可以得到
Figure FDA00023823496300000113
S2.1.2线性输出步骤：对于每一个i＝1,2…m，计算
Figure FDA00023823496300000114
上式中
Figure FDA0002382349630000021
表示第k次算法迭代过程中的
Figure FDA0002382349630000022
值，aij表示字典矩阵Α的第i行第j列的元素，(…)i表示矢量的第i个元素，|…|表示复数的模，
Figure FDA0002382349630000023
表示第k次算法迭代过程中的
Figure FDA0002382349630000024
值，
Figure FDA0002382349630000025
表示第k次算法迭代过程中的
Figure FDA0002382349630000026
值，
Figure FDA0002382349630000027
表示第k次算法迭代过程中的
Figure FDA0002382349630000028
值，
Figure FDA0002382349630000029
表示第k次算法迭代过程中的
Figure FDA00023823496300000210
值，
Figure FDA00023823496300000211
表示第k次算法迭代过程中的
Figure FDA00023823496300000212
值；
S2.1.3非线性输出步骤：对于每一个i＝1,2…m，计算
Figure FDA00023823496300000213
yi表示接收数据y的第i个元素，
Figure FDA00023823496300000214
表示第k次算法迭代过程中的
Figure FDA00023823496300000215
值，
Figure FDA00023823496300000216
表示在第第k次算法迭代过程中更新的
Figure FDA00023823496300000217
的值，上面的函数
Figure FDA00023823496300000218
S2.1.4线性输入步骤：对于每一个j＝1,2…N，计算
Figure FDA00023823496300000219
Figure FDA00023823496300000220
表示第k次算法迭代过程中的
Figure FDA00023823496300000221
值，
Figure FDA00023823496300000222
表示第k次算法迭代过程中的
Figure FDA00023823496300000223
值，这里(…)-1表示矩阵求逆，(…)*表示复数的共轭；
S2.1.5非线性输入步骤：对于每一个j＝1,2…N，计算
Figure FDA00023823496300000224
这里
Figure FDA00023823496300000225
表示第k+1次迭代的
Figure FDA00023823496300000226
值，
Figure FDA00023823496300000227
表示第k+1次迭代的
Figure FDA00023823496300000228
值，上面的函数
Figure FDA00023823496300000229
S2.1.6判断AMP算法是否收敛：计算
Figure FDA00023823496300000230
的值，其中||…||1表示矩阵的1范数，
Figure FDA00023823496300000231
表示第k+1次迭代的
Figure FDA00023823496300000232
值，同理，
Figure FDA00023823496300000233
表示第k次迭代的
Figure FDA00023823496300000234
值；如果该值大于某一设定门限ε1，则回到S2.1.2重新迭代；否则跳出循环进入S2.2得到p(xj|y)的结果；门限ε1取决于***的信噪比等因素，需要根据实际情况进行调整；
S2.2通过上述步骤可以得到t时刻的信号后验概率密度函数p(xj|y)的结果如下
Figure FDA0002382349630000031
p(xj)表示xj的概率密度函数；上面公式中的
Figure FDA0002382349630000032
为S2.1步骤中根据t时刻的样本值y(t)进行计算获得，为AMP算法被判断收敛后，最后一次迭代过程中的
Figure FDA0002382349630000033
的值；这里γj的值由S1或者S3获得，在第一次EM算法迭代中，γj的值由S1中的初始值确定，在其他情况下，γj的值由上一次EM算法循环中S3计算出来的γj确定；
S3利用EM算法更新待估计参量γj,j＝1,2…N，求出目标数量最优解K*，更新噪声功率σ0
在S2中已经获得了信号的后验概率密度函数p(xj|y)，根据EM算法，本步骤利用下面的表达式来逐个更新待估计参量的值
Figure FDA0002382349630000034
上式中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的更新表达式为
Figure FDA0002382349630000035
上面的表达式中
Figure FDA0002382349630000036
表示第i次算法迭代过程中的γj，
Figure FDA0002382349630000037
表示第i+1次算法迭代过程中的γj；
进一步对γj求偏导可以得到
Figure FDA0002382349630000038
S3.2估算目标数量最优解K*：
采用子空间分析法，目标数量最优解为
Figure FDA0002382349630000039
其中I是单位矩阵，
Figure FDA0002382349630000041
Figure FDA0002382349630000042
K是目标数量的估计值，K*为K的最优解，
Figure FDA0002382349630000043
由A中对应假定目标角度
Figure FDA0002382349630000044
的列组成的矩阵，假定目标角度
Figure FDA0002382349630000045
指的是γj,j＝1,2…N中最大的K个值所对应的角度；
S3.3更新噪声功率σ0：
根据公式：
Figure FDA0002382349630000046
其中
Figure FDA0002382349630000047
这里
Figure FDA0002382349630000048
是
Figure FDA0002382349630000049
中各列的协方差矩阵，从上式可得
Figure FDA00023823496300000410
根据S3.2定义的映射矩阵P，且P＝PH＝P2，可得
PRPH-σ0PPH＝ΣY-σ0I (18)
结合tr(R-ΣY)≈0可以推出σ0的更新公式为：
Figure FDA00023823496300000411
S4用梯度下降法更新字典矩阵A中的
Figure FDA00023823496300000412
由S3.2中计算出的目标数量最优解K*，可认为目标对应角度为γj,j＝1,2…N中前K*个最大的γj所对应的角度网格节点
Figure FDA00023823496300000413
使用梯度下降法对这些目标粗略对应的角度网格节点进行调整，需要注意的是这里只需要更新目标当前对应的角度网格节点
Figure FDA00023823496300000414
不必更新整个角度网格，以免提高运算复杂度，即对于S3.2中确定的每一个目标所对应的角度
Figure FDA00023823496300000415
进行本步骤迭代计算，直到算法收敛，S4步骤一共需要执行K*次；可由下式推导更新规则
Figure FDA00023823496300000416
其中
Figure FDA00023823496300000417
A是带角度变量的字典矩阵；根据上式可求出关于
Figure FDA00023823496300000418
的梯度
Figure FDA00023823496300000419
这里
Figure FDA00023823496300000420
表示矩阵的实部，X是由S2步骤求出的对X的估计量，
Figure FDA00023823496300000421
是A中对应于除
Figure FDA00023823496300000422
以外角度
Figure FDA00023823496300000423
的列所组成的矩阵，
Figure FDA00023823496300000424
表示X中对应于除
Figure FDA00023823496300000425
以外角度
Figure FDA00023823496300000426
的行所组成的矩阵，
Figure FDA00023823496300000427
表示X中对应于角度
Figure FDA00023823496300000428
的行矩阵，
Figure FDA00023823496300000429
是A中的列向量，d(·)是a(·)的导数；
通过(21)式对
Figure FDA00023823496300000430
进行优化，有
Figure FDA0002382349630000051
其中，α是角度更新的步长，实际应用中步长取决于对角度估计的精度要求，(·)l表示本步骤内部迭代第l次的值，如果函数sign(·)的输入是正数则等于1，反之为-1；
当完成一次角度更新后，需要判断是否收敛；考虑到初始的字典网格可能已经满足精度需求，所以此步骤中迭代次数不会太多，当满足条件
Figure FDA0002382349630000052
或者l≥ε4时，认为对
Figure FDA0002382349630000053
的估计已经满足要求，退出循环，ε3、ε4为判决阈值；其中，ε3可以根据***对角度估计的精度要求确定；ε4根据***的实时性要求确定，值越大本算法的实时性就越低；
S5判断算法迭代过程何时收敛，以及确定来波方向和数量
完成对γ＝[γ1 γ2...γN]T以及σ0的值更新后，用如下表达式判断收敛：
Figure FDA0002382349630000054
ε2为判定门限，可根据***实际进行设定；如果上式不成立，则返回S2继续进行迭代运算；如果上式成立，则可退出循环，来波数量即为最后一次计算出的K*，方向为γ＝[γ1γ2…γN]T中前K*个最大的γj所对应的方向
Figure FDA0002382349630000055
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之间。
CN202010086824.0A 2020-02-11 2020-02-11 一种基于近似消息传递的不在网格目标角度估计方法 Active CN111257845B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010086824.0A CN111257845B (zh) 2020-02-11 2020-02-11 一种基于近似消息传递的不在网格目标角度估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010086824.0A CN111257845B (zh) 2020-02-11 2020-02-11 一种基于近似消息传递的不在网格目标角度估计方法

Publications (2)

Publication Number Publication Date
CN111257845A true CN111257845A (zh) 2020-06-09
CN111257845B CN111257845B (zh) 2020-09-22

Family

ID=70945657

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010086824.0A Active CN111257845B (zh) 2020-02-11 2020-02-11 一种基于近似消息传递的不在网格目标角度估计方法

Country Status (1)

Country Link
CN (1) CN111257845B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112346123A (zh) * 2020-11-06 2021-02-09 中国地震灾害防御中心 一种地震数据处理via双参数分析方法
CN113281702A (zh) * 2021-04-30 2021-08-20 中国人民解放军战略支援部队信息工程大学 协同短波多站角度与卫星时频的超视距目标直接定位方法
CN116879862A (zh) * 2023-09-08 2023-10-13 西安电子科技大学 基于分层稀疏迭代的单快拍稀疏阵空间角度超分辨方法
CN117388791A (zh) * 2023-09-13 2024-01-12 桂林电子科技大学 一种6gisca***宽带信号doa估计算法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104375133A (zh) * 2014-11-11 2015-02-25 西北大学 一种空间二维doa的估算方法
KR101783777B1 (ko) * 2016-05-17 2017-10-10 국방과학연구소 공분산 Fitting을 통한 도래각 추정 방법
CN109116293A (zh) * 2018-08-22 2019-01-01 上海师范大学 一种基于离格稀疏贝叶斯的波达方向估计方法
CN109490819A (zh) * 2018-11-16 2019-03-19 南京邮电大学 一种基于稀疏贝叶斯学习的离格波达方向估计方法
CN109658467A (zh) * 2018-12-12 2019-04-19 浙江工业大学 一种基于多字典改进型压缩感知框架的内窥镜图像感知重构方法
CN109752710A (zh) * 2019-01-07 2019-05-14 中国人民解放军国防科技大学 一种基于稀疏贝叶斯学习的快速目标角度估计方法
CN110208735A (zh) * 2019-06-12 2019-09-06 西北工业大学 一种基于稀疏贝叶斯学习的相干信号doa估计方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104375133A (zh) * 2014-11-11 2015-02-25 西北大学 一种空间二维doa的估算方法
KR101783777B1 (ko) * 2016-05-17 2017-10-10 국방과학연구소 공분산 Fitting을 통한 도래각 추정 방법
CN109116293A (zh) * 2018-08-22 2019-01-01 上海师范大学 一种基于离格稀疏贝叶斯的波达方向估计方法
CN109490819A (zh) * 2018-11-16 2019-03-19 南京邮电大学 一种基于稀疏贝叶斯学习的离格波达方向估计方法
CN109658467A (zh) * 2018-12-12 2019-04-19 浙江工业大学 一种基于多字典改进型压缩感知框架的内窥镜图像感知重构方法
CN109752710A (zh) * 2019-01-07 2019-05-14 中国人民解放军国防科技大学 一种基于稀疏贝叶斯学习的快速目标角度估计方法
CN110208735A (zh) * 2019-06-12 2019-09-06 西北工业大学 一种基于稀疏贝叶斯学习的相干信号doa估计方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Z.M.LIU等: "An effeicient maximum likelihood method for direction of arrival estimation via sparse bayesian learning", 《IEEE TRANS. WIRELESS COMMUNICATIONS》 *
刘本源: "快速块稀疏贝叶斯学习算法的理论与应用", 《中国博士学位论文全文数据库 信息科技辑》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112346123A (zh) * 2020-11-06 2021-02-09 中国地震灾害防御中心 一种地震数据处理via双参数分析方法
CN113281702A (zh) * 2021-04-30 2021-08-20 中国人民解放军战略支援部队信息工程大学 协同短波多站角度与卫星时频的超视距目标直接定位方法
CN113281702B (zh) * 2021-04-30 2024-02-09 中国人民解放军战略支援部队信息工程大学 协同短波多站角度与卫星时频的超视距目标直接定位方法
CN116879862A (zh) * 2023-09-08 2023-10-13 西安电子科技大学 基于分层稀疏迭代的单快拍稀疏阵空间角度超分辨方法
CN116879862B (zh) * 2023-09-08 2023-12-01 西安电子科技大学 基于分层稀疏迭代的单快拍稀疏阵空间角度超分辨方法
CN117388791A (zh) * 2023-09-13 2024-01-12 桂林电子科技大学 一种6gisca***宽带信号doa估计算法

Also Published As

Publication number Publication date
CN111257845B (zh) 2020-09-22

Similar Documents

Publication Publication Date Title
CN111257845B (zh) 一种基于近似消息传递的不在网格目标角度估计方法
WO2018094565A1 (zh) 脉冲噪声下的波束成形方法及装置
CN109490819B (zh) 一种基于稀疏贝叶斯学习的离格波达方向估计方法
CN106646344B (zh) 一种利用互质阵的波达方向估计方法
CN108802683B (zh) 一种基于稀疏贝叶斯学习的源定位方法
CN107450045B (zh) 基于focuss二次加权算法的doa估计方法
CN108872926B (zh) 一种基于凸优化的幅相误差校正及doa估计方法
CN109239649B (zh) 一种阵列误差条件下的互质阵列doa估计新方法
CN110244272B (zh) 基于秩一去噪模型的波达方向估计方法
CN111337873A (zh) 一种基于稀疏阵的doa估计方法
CN114720938A (zh) 基于深度展开的大规模天线阵列单比特采样doa估计方法
CN117092585B (zh) 单比特量化DoA估计方法、***和智能终端
CN109783960B (zh) 一种基于网格部分细化的波达方向估计方法
CN115236584A (zh) 基于深度学习的米波雷达低仰角估计方法
CN116224219A (zh) 一种阵列误差自校正原子范数最小化doa估计方法
CN112731273B (zh) 一种基于稀疏贝叶斯的低复杂度信号波达方向估计方法
CN113759303B (zh) 一种基于粒子群算法的无网格波达角估计方法
CN111880143B (zh) 改进稀疏贝叶斯学习的高精度定位方法、存储介质及设备
Yang et al. A correlation-aware sparse Bayesian perspective for DOA estimation with off-grid sources
CN116299193A (zh) 一种mimo雷达智能doa估计方法
CN115015832A (zh) 一种非均匀噪声下大规模阵列幅相误差及目标方位联合估计方法
CN114415106A (zh) 基于改进lamp网络的互耦阵列doa估计方法
Tan et al. An iterative adaptive dictionary learning approach for multiple snapshot DOA estimation
CN110471026B (zh) 一种相位增强的米波雷达目标低仰角doa估计方法
CN109298384B (zh) 一种基于变分贝叶斯推断的非均匀直线阵波达方向角估计方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant