一种波达角的高精度估计方法与流程

1.本发明涉及阵列信号处理领域，尤其涉及一种阵列信号处理中的波达角的高精度估计方法。


背景技术：

2.利用天线阵列估计多个远场信号的波达角(doa)是阵列信号处理的基础问题，在雷达、声呐和地震预测等领域得到了广泛的应用。针对doa估计问题，国内外研究者提出了多种不同类型的方法，如子空间分解、支持向量机、矩阵特征空间分解等，每种方法都有各自的优缺点。随着压缩感知技术的出现，研究者们将doa估计归纳为稀疏恢复问题，并提出了多种不同的估计算法，如基于正交匹配追踪算法等。通常的稀疏估计方法中，将波达角范围划分为均匀网格，假定实际的波达角准确位于某个网格上。相比于子空间划分、特征空间分解等方法，压缩感知类方法通常能够获得更好的性能，特别是在快拍数量较少的情况下。稀疏贝叶斯学习是一类经典稀疏恢复算法，通过给变量增加稀疏引导先验，能够以较少的观测实现更优的估计性能。但是当真实来波方向并不精确位于网格上时，上述方案会导致较严重的估计偏差。更精细的网格划分一方面会导致较高的运算复杂度，另一方面会使得感知矩阵的列向量之间具有较高的相关度，导致估计算法收敛性变差。
3.近年来国内外多个团队提出了多种基于离格信号的稀疏重构算法，假定真实角度不必限于网格之上，能够计算出误差并对估计值进行修正。例如求根稀疏贝叶斯学习方法(rsbl)，在线性阵列场景下能够取得很好的估计性能。由于复杂度较低，经典的近似消息传递算法(amp)在稀疏贝叶斯学习框架下得到了广泛应用，但是普通的amp算法在遇到较非高斯分布的观测矩阵时会导致迭代过程中出现发散。


技术实现要素：

4.本发明的目的是提供一种波达角的高精度估计方法，能够解决离格波达角估计中复杂度和鲁棒性的矛盾，有效提高波达角的计算精度，在机载雷达等快拍数较少的场景下具有明显的性能优势。
5.本发明采用下述技术方案：
6.一种波达角的高精度估计方法，依次包括以下步骤：
7.a：首先建立离格doa双线性系统模型y＝(a+∑nβ
nen
)s+w＝φ{β}s+w；
8.其中，φ{β}＝a+ediag(β)＝a+∑nβ
nen
，en＝[0,
…
,e(θn),
…
,0]∈cm×n，离格doa双线性系统模型为接收信号向量y＝(a(θ)+e(θ)diag(β))s+w的简记；y＝(a(θ)+e(θ)diag(β))s+w为y
t
＝(a(θ)+e(θ)diag(β))s
t
+w
t
的多重观测形式；设配置有m个无方向阵元所构成的均匀直线阵列，存在k个窄带远场来波信号，y
t
表示第t个快拍下的接收信号向量，a(θ)＝[a(θ1),...,a(θn)]∈cm×n，cm×k表示维度为m
×
k的复值矩阵；设来波信号的角度空间θ被划分为n个网格，表示为向量形式θ＝[θ1,
…
,θn]，每个网格大小为δ
θ
＝π/n；矩阵e(θ)＝[e(θ1),
…
,e(θn)]表示误差矩阵，向量β＝[β1,
…
βn]
t
表示离格误差参数向量，离格误差参数
n＝1,
…
,n，s
t
∈cn×1表示长度为n的待求稀疏来波信号向量；l为快拍数量，w
t
表示均值为0且方差为1/λ的高斯白噪声向量；矩阵y＝[y1,y2,
…
,y
l
]、s＝[s1,s2,
…
,s
l
]和w＝[w1,w2,
…
,w
l
]分别表示接收信号向量、来波信号向量和高斯白噪声向量的集合矩阵，矩阵s具有相同的稀疏特征；a和e分别为阵列流型矩阵a(θ)和误差矩阵e(θ)的简记；θ为来波信号的角度空间；
[0009]
然后对离格doa双线性系统模型利用贝叶斯公式进行因式分解；
[0010][0011]
其中，将稀疏来波信号向量s
t
定义为两个相同的等效向量和和均表示稀疏来波信号向量，数学表达为δ(
·
)为delta函数；表示在第t个快拍中，接收信号向量y
t
的似然函数，似然函数表示为复高斯形式其中表示均值为协方差矩阵为λ-1
i的复高斯分布，i表示对角阵；概率分布代表向量的先验分布，γ表示超先验参数向量；所有来波信号向量t＝1:l都具有共同的稀疏特性，共享同一个超先验参数向量γ；离格误差参数向量β服从均匀分布p(β)＝u[-δ
θ
/2,δ
θ
/2]，噪声方差的先验分布假设为p(λ)＝1/λ；
[0012]
b：对稀疏来波信号s进行估计；
[0013]
c：对离格误差参数向量β进行估计；
[0014]
d：根据得到的稀疏来波信号s和离格误差参数向量β，通过步骤a得到的离格doa双线性系统模型，进行波达角的高精度计算。
[0015]
所述的步骤b中，应用近似消息传递算法对稀疏来波信号s进行估计。
[0016]
所述的步骤c中，应用采用期望最大化方法对离格误差参数向量β进行估计
[0017]
所述的步骤a中，离格doa双线性系统模型的建立步骤如下，
[0018]
a1：设配置有m个无方向阵元所构成的均匀直线阵列，存在k个窄带远场来波信号，则第t个快拍下的接收向量表示为
[0019][0020]
其中，y
t
∈cm×1为第t个快拍下全部m个阵元的接收信号向量，符号cm×1表示维度为m
×
1的复值向量，为阵列流型矩阵，符号cm×k表示维度为m
×
k的复值矩阵，表示第t个快拍下的真实来波信号向量，ck×1表示维度为k
×
1的复值向量，w
t
表示均值为0且方差为1/λ的高斯白噪声向量，λ表示噪声精度；
[0021]
a2：将阵列流型矩阵分解为向量形式
[0022]
其中，向量表示第k个来波信号的导向矢量，分解表示为
j表示虚数单位，表示第k个来波信号的真实波达角度，k＝1,
…
,k；来波信号的角度空间θ被划分为n个网格，表示为向量形式θ＝[θ1,
…
,θn]，每个网格大小为δ
θ
＝π/n；每个角度θn均对应于一个潜在的来波信号源n＝1,
…
,n；来波信号方向在角度域是稀疏分布的，即n＞＞k，构成稀疏来波信号向量
[0023]
a3：将阵列流型矩阵重构为a(θ)＝[a(θ1),...,a(θn)]∈cm×n；
[0024]
a4：引入离格模型，假定信号源的真实波达角度与任何一个网格都不完全重合，即假定存在离格误差参数根据离格模型，将导向矢量一阶泰勒展开为：
[0025][0026]
其中，离格误差参数βn定义为向量表示误差向量，其中表示对导向矢量依波达角度求导数；下标nk表示距离真实波达角度最近的网格，表示网格nk对应的波达角度；
[0027]
a5：根据一阶泰勒展开，将第t个快拍下的接收信号向量重写为：
[0028]yt
＝(a(θ)+e(θ)diag(β))s
t
+w
t
；
[0029]
其中，矩阵e(θ)＝[e(θ1),
…
,e(θn)]表示误差矩阵，向量β＝[β1,
…
βn]
t
表示离格误差参数向量，s
t
∈cn×1表示长度为n的待求稀疏来波信号向量；doa估计模型中有l个快拍时，
[0030]
a6：当doa估计模型中有l个快拍时，将步骤a5中的第t个快拍下的接收信号向量写为多重观测形式：
[0031]
y＝(a(θ)+e(θ)diag(β))s+w；
[0032]
其中，矩阵y＝[y1,y2,
…
,y
l
]、s＝[s1,s2,
…
,s
l
]和w＝[w1,w2,
…
,w
l
]分别表示接收信号向量、来波信号向量和高斯白噪声向量的集合矩阵，矩阵s具有相同的稀疏特征，即每列非零元素位置相同，波信号的角度空间θ为固定值；
[0033]
a7：将阵列流型矩阵a(θ)和误差矩阵e(θ)分别简记为a和e，将步骤a6中的多重观测形式简记为
[0034]
y＝(a+ediag(β))s+w；
ꢀꢀꢀꢀꢀꢀꢀ
(1)
[0035]
a8：由于(1)式中的表达式(a+ediag(β))只有离格误差参数向量β为未知，则将其重新定义为
[0036]
φ{β}＝a+ediag(β)＝a+∑nβ
nen
；
ꢀꢀꢀꢀꢀ
(2)
[0037]
其中，βn表示离格误差参数，en定义为en＝[0,
…
,e(θn),
…
,0]∈cm×n，将(1)式中的doa估计模型转化为双线性形式
[0038]
y＝(a+∑nβ
nen
)s+w＝φ{β}s+w
ꢀꢀꢀꢀꢀ
(3)
[0039]
其中，s＝[s1,...,s
l
]表示l个待求稀疏来波信号向量；只有当时，
βn取非零值，向量β和矩阵s具有完全相同的稀疏特性；
[0040]
a9：根据公式(3)所示的离格doa双线性系统模型，利用数据之间的约束关系，对系统中所有已知和未知变量的联合分布函数，利用贝叶斯公式进行因式分解，
[0041][0042]
所述的步骤b包括以下具体步骤：
[0043]
b1：初始化离格误差参数向量β和中间变量为长度为n的全零向量，中间变量初始化为1，噪声精度初始化为λ＝1，超先验参数向量γ初始化为长度n的全1向量；然后进入步骤b2；
[0044]
b2:计算来稀疏波信号向量的后验概率为高斯形式其中均值和方差分别计算为
[0045][0046]
然后进入步骤b3；
[0047]
式中，运算符号《
·
》表示对向量求平均，i表示对角阵，y表示接收信号向量；
[0048]
b3:利用步骤b2得到的均值和方差以及步骤b1中初始化的中间变量和计算中间变量和分别为
[0049][0050]
然后进入步骤b4；
[0051]
b4:利用步骤b1中初始化的超先验参数γ以及步骤b3得到的中间变量和计算中间变量的后验分布为高斯形式其中均值和方差分别计算为
[0052][0053]
然后进入步骤b5；
[0054]
b5:利用步骤b4得到的中间变量的均值和方差计算超先验参数γ为
[0055][0056]
然后进入步骤b6；其中l表示doa估计中的快拍个数；
[0057]
b6:利用步骤b4得到的中间变量的均值和方差以及步骤b3得到的中间变量和计算中间变量和
[0058]
[0059]
然后进入步骤b7；
[0060]
b7:判断是否到达设定的迭代次数，若到达则结束迭代；若未到达则返回步骤b2；最终确定稀疏来波信号向量
[0061]
所述的步骤b2至b6为迭代过程，迭代次数设定为15次。
[0062]
所述的步骤c包括以下具体步骤：
[0063]
c1:初始化离格误差参数β为长度为n的全0向量，将β代入公式(2)构建矩阵φ{β}；然后进入步骤c2；
[0064]
c2:利用步骤b2中求得的方差t＝1:l计算中间变量矩阵d为
[0065][0066]
然后进入步骤c3；
[0067]
c3:利用步骤b7中求得的稀疏来波信号向量和公式(2)中定义的en,n＝1,2,...,n分别构造矩阵h∈cn×n和向量α∈cn×1为
[0068]
α＝(α1…
αi…
αn)
t
；
[0069]
其中，矩阵和向量中的元素h
ij
,i,j＝1,2,...,n，αi,i＝1,2,...,n，
[0070][0071][0072]
i,j表示矩阵h和向量α的元素标号，矩阵ei,ej与公式(2)中的矩阵en具有相同的定义；然后进入步骤c4；
[0073]
c4:求离格误差参数向量β的估计值，β＝h-1
α；然后进入步骤c5；
[0074]
c5:判断是否到达设定的迭代次数，若到达则结束迭代；若未到达则返回步骤c2；最终得到离格误差参数向量β。
[0075]
所述的，步骤c2至c4为迭代过程，迭代次数设定为5次。
[0076]
本发明将doa估计模型转换为双线性问题，利用双线性vamp算法解决该问题。本发明首先将离格doa估计模型转换为双线性问题，利用双线性vamp算法对双线性模型进行推导和运算；然后建立doa估计仿真模型对所提算法与文献中已有方法进行对比，结果表明本发明所提算法在机载雷达等快拍数较少的场景下具有明显的性能优势。
附图说明
[0077]
图1为本发明的流程示意图；
[0078]
图2为本发明在一次蒙特卡洛仿真中估计值和真实值的对比图；
[0079]
图3为本发明与归一化均方误差随快拍次数的变化曲线图；
[0080]
图4为本发明与现有算法在快拍个数l为2时估计性能随信噪比的变化曲线图；
[0081]
图5为本发明与现有算法在快拍个数l为4时估计性能随信噪比的变化曲线图；
[0082]
图6为本发明与现有算法在快拍个数l为6时估计性能随信噪比的变化曲线图；
[0083]
图7为本发明与现有算法在快拍个数l为8时估计性能随信噪比的变化曲线图。
具体实施方式
[0084]
以下结合附图和实施例对本发明作以详细的描述：
[0085]
如图1所示，本发明所述的波达角的高精度估计方法，包括以下步骤：
[0086]
a：建立离格doa双线性系统模型y＝(a+σnβ
nen
)s+w＝φ{β}s+w，并利用贝叶斯公式进行因式分解；
[0087]
a1：设配置有m个无方向阵元所构成的均匀直线阵列，存在k个窄带远场来波信号，则第t个快拍下的接收向量表示为
[0088][0089]
其中，y
t
∈cm×1为第t个快拍下全部m个阵元的接收信号向量，符号cm×1表示维度为m
×
1的复值向量，为阵列流型矩阵，符号cm×k表示维度为m
×
k的复值矩阵，表示第t个快拍下的真实来波信号向量，ck×1表示维度为k
×
1的复值向量，w
t
表示均值为0且方差为1/λ的高斯白噪声向量，其中λ表示噪声精度。a2：将阵列流型矩阵分解为向量形式：
[0090][0091]
其中，向量表示第k个来波信号的导向矢量，可以分解表示为j表示虚数单位，表示第k个来波信号的真实波达角度，k＝1,
…
,k；由于doa估计的目标是利用接收向量y
t
估计所有k个来波信号的波达角度和幅度，因此在本发明中，将来波信号的角度空间θ划分为n个网格，表示为向量形式θ＝[θ1,
…
,θn]，每个网格大小为δ
θ
＝π/n。假设每个角度θn均对应于一个潜在的来波信号源n＝1,
…
,n；由于来波信号方向在角度域是稀疏分布的，即n＞＞k，从而构成稀疏来波信号向量
[0092]
a3：将阵列流型矩阵重构为a(θ)＝[a(θ1),...,a(θn)]∈cm×n；
[0093]
a4：由于要提高doa的精度，最直接的方法是将来波信号的角度空间θ划分为更小的网格，即加大网格数量n。但网格数量n的加大会导致复杂度提升，且使得矩阵a(θ)的相关性变大，而相关性高的矩阵a(θ)将会导致估计性能变差，引起消息传递等迭代类算法的发散，造成估计失败。为了避免上述问题，本发明中引入离格模型，即假定信号源的真实波达角度与任何一个网格都不完全重合，即假定存在离格误差参数根据离格模型，将导向矢量一阶泰勒展开为：
[0094]
[0095]
其中，离格误差参数βn定义为向量表示误差向量，其中表示对导向矢量依波达角度求导数；下标nk表示距离真实波达角度最近的网格，表示网格nk对应的波达角度。
[0096]
a5：根据一阶泰勒展开，将第t个快拍下的接收信号向量重写为：
[0097]yt
＝(a(θ)+e(θ)diag(β))s
t
+w
t
；
[0098]
其中，矩阵e(θ)＝[e(θ1),
…
,e(θn)]表示误差矩阵，向量β＝[β1,
…
βn]
t
表示离格误差参数向量，s
t
∈cn×1表示长度为n的待求稀疏来波信号向量。
[0099]
a6：当doa估计模型中有l个快拍时，将上式写为多重观测形式
[0100]
y＝(a(θ)+e(θ)diag(β))s+w；
[0101]
其中，矩阵y＝[y1,y2,
…
,y
l
]、s＝[s1,s2,
…
,s
l
]和w＝[w1,w2,
…
,w
l
]分别表示接收信号向量、来波信号向量和高斯白噪声向量的集合矩阵，且矩阵s具有相同的稀疏特征，即每列非零元素位置相同。上式来波信号的角度空间θ为固定值。
[0102]
a7：将阵列流型矩阵a(θ)和误差矩阵e(θ)分别简记为a和e，则doa估计模型简记为
[0103]
y＝(a+ediag(β))s+w；
ꢀꢀꢀꢀꢀꢀꢀ
(1)
[0104]
观察可知，上述doa估计模型中未知参数有s和β，在信号处理领域，上述doa估计模型属于利用一组观测同时估计两组独立参数的双线性模型。
[0105]
a8：由于(1)式中的表达式(a+ediag(β))只有离格误差参数向量β为未知，则可以将其重新定义为
[0106]
φ{β}＝a+ediag(β)＝a+∑nβ
nen
；
ꢀꢀꢀꢀꢀ
(2)
[0107]
其中，βn表示离格误差参数，en定义为en＝[0,
…
,e(θn),
…
,0]∈cm×n，则将(1)式中的doa估计模型转化为双线性形式
[0108]
y＝(a+σnβ
nen
)s+w＝φ{β}s+w
ꢀꢀꢀꢀꢀ
(3)
[0109]
其中，s＝[s1,...,s
l
]表示l个待求稀疏来波信号向量。根据(1)式中向量β的物理含义可知，只有当时βn取非零值，从而可确定向量β和矩阵s具有完全相同的稀疏特性，即非零元素位置一样。
[0110]
a9：根据公式(3)所示的离格doa双线性系统模型，利用数据之间的约束关系，对系统中所有已知和未知变量的联合分布函数，利用贝叶斯公式进行因式分解，
[0111][0112]
上式中，将稀疏来波信号向量s
t
定义为两个相同的等效向量和三者完全相等，即均表示稀疏来波信号向量，数学表达为此处δ(
·
)为delta函数；表示在第t个快拍中，接收信号向量y
t
的似然函数，由于噪声服从方差为1/λ的复高斯分布，则似然函数可以表示为复高斯形式其中表示均值为
协方差矩阵为λ-1
i的复高斯分布，i表示对角阵；概率分布代表向量的先验分布，γ表示超先验参数向量。由于本发明假定对于所有来波信号向量t＝1:l都具有共同的稀疏特性，则共享同一个超先验参数向量γ。离格误差参数向量β服从均匀分布p(β)＝u[-δ
θ
/2,δ
θ
/2]，噪声方差的先验分布假设为p(λ)＝1/λ。根据(4)式所示双线性模型的因式分解，可以利用双线性近似消息传递和期望最大化算法对双线性模型进行算法推导。具体步骤分为稀疏来波信号s和离格误差参数β估计两部分。
[0113]
b：对稀疏来波信号s进行估计；
[0114]
b1：设离格误差参数向量β已知，则矩阵φ{β}也为已知，使得公式(3)所示双线性模型退化为单线性模型。在解决单线性问题时，应用近似消息传递算法对稀疏来波信号向量进行估计，具体估计过程如下：
[0115]
初始化离格误差参数向量β和中间变量为长度为n的全零向量，中间变量初始化为1，噪声精度初始化为λ＝1，超先验参数向量γ初始化为长度n的全1向量；然后进入步骤b2；
[0116]
b2:计算来稀疏波信号向量的后验概率为高斯形式其中均值和方差分别计算为
[0117][0118]
然后进入步骤b3；
[0119]
式中，运算符号《
·
》表示对向量求平均，i表示对角阵，y表示接收信号向量。本发明中使用方差向量的平均，能够提升算法的鲁棒性。
[0120]
b3:利用步骤b2得到的均值和方差以及步骤b1中初始化的中间变量和计算中间变量和分别为
[0121][0122]
然后进入步骤b4；
[0123]
b4:利用步骤b1中初始化的超先验参数γ以及步骤b3得到的中间变量和计算中间变量的后验分布为高斯形式其中均值和方差分别计算为
[0124][0125]
然后进入步骤b5；
[0126]
b5:利用步骤b4得到的中间变量的均值和方差计算超先验参数γ为
[0127][0128]
然后进入步骤b6；其中l表示doa估计中的快拍个数。
[0129]
b6:利用步骤b4得到的中间变量的均值和方差以及步骤b3得到的中间变量和计算中间变量和
[0130][0131]
然后进入步骤b7；
[0132]
b7:判断是否到达设定的迭代次数，若到达则结束迭代；若未到达则返回步骤b2；最终确定稀疏来波信号向量
[0133]
本实施例中，步骤b2至b6为迭代过程，迭代次数设定为15次。
[0134]
c：对离格误差参数向量β进行估计；
[0135]
本发明对离格误差参数向量β的估计采用期望最大化方法，计算步骤如下：
[0136]
c1:初始化离格误差参数β为长度为n的全0向量，将β代入公式(2)构建矩阵φ{β}；然后进入步骤c2；
[0137]
c2:利用步骤b2中求得的方差t＝1:l计算中间变量矩阵d为
[0138][0139]
然后进入步骤c3；
[0140]
c3:利用步骤b7中求得的稀疏来波信号向量和公式(2)中定义的en,n＝1,2,...,n分别构造矩阵h∈cn×n和向量α∈cn×1为
[0141]
α＝(α1…
αi…
αn)
t
；
[0142]
其中，矩阵和向量中的元素h
ij
,i,j＝1,2,...,n，αi,i＝1,2,...,n，
[0143][0144][0145]
其中i,j表示矩阵h和向量α的元素标号，此处矩阵ei,ej与公式(2)中的矩阵en具有相同的定义；然后进入步骤c4；
[0146]
c4：求离格误差参数向量β的估计值，β＝h-1
α；然后进入步骤c5；
[0147]
c5：判断是否到达设定的迭代次数，若到达则结束迭代；若未到达则返回步骤c2；最终得到离格误差参数向量β。
[0148]
本实施例中，步骤c2至c4为迭代过程，迭代次数设定为5次。
[0149]
d：利用步骤b中求得的稀疏来波信号向量s，和步骤c中求得的离格误差参数β，结合步骤a中建立的doa估计模型y＝(a+∑nβ
nen
)s+w＝φ{β}s+w，最终得到波达角的高精度估计值。
[0150]
复杂度方面，本发明可以划分为b和c两个步骤，其中b2步需要矩阵乘法φ
t
{β}φ{β}，具有复杂度b步的其余运算均为标量形式，复杂度为其中符号表示复杂度的阶数。c步骤需要计算矩阵求逆，复杂度为所以本发明的整体复杂度可写为
[0151]
为了说明本发明的性能，可以建立仿真环境进行数值仿真，并且与国内外现有算法进行对比。仿真环境选择均匀线性阵列，阵元个数为m＝30，则虚拟网格个数为n＝30，即网格宽度为6
°
，设置快拍数为l＝2-12，设置阈值γ
max
＝103。为衡量本发明的有效性，仿真中选择已有的网格化稀疏贝叶斯学习(ogsbl)、求根sbl(rtsbl)和高斯-赛尔德根(gsroot)作为对比，此外还选择克拉美罗界(crb)作为参考。
[0152]
图2为一次蒙特卡洛仿真中估计值estimated和真实值true angle的对比图，仿真中选择真实角度为[-28.6,-18.6,3.5,15.6,31.7]。图2中横坐标为角度，纵坐标为估计向量s2的绝对值。本仿真中设置信噪比snr＝30db。从图2中可以看出，本发明能够精确捕捉到5个来波信号的离格角度和准确功率，在无来波信号的角度估计结果趋近于0。，
[0153]
图3为本发明与已有算法归一化均方误差(nmse)随快拍次数的变化曲线图。数值仿真结果为500次蒙特卡洛仿真的平均值。每次仿真中波达方向以[-30,-18,6,18,30]为基础，叠加(-3,3)的均匀随机角度偏移产生。图3中设置信噪比snr＝20db，来波信号个数为k＝5。从仿真结果中可以看出，本发明即图中bad-vamp在快拍数较多时，性能与现有的rtsbl和gsroot较为接近。但是随快拍数的减少，rtsbl和gsroot方法性能明显恶化，使得本发明表现出较大的性能增益。从图3中还能看出，由于不考虑网格偏差，ogsbl算法性能较差，从而证明了离格估计算法的必要性。
[0154]
图4至图7分别为快拍个数l为2、4、6和8时估计性能随信噪比的变化曲线图。从图4至图7中可以得出与图3相同的结论，即在快拍个数较少时本发明具有明显的性能优势。要达到相同的估计性能，本发明需要的快拍个数较少，更适合于机载雷达等快时变场景。