基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法

基于稀疏贝叶斯学习的机载三维异构阵doa估计方法
技术领域
1.本发明涉及雷达信号处理技术领域，具体涉及基于稀疏贝叶斯学习的机载三维异构阵doa(波达方向)估计方法。


背景技术：

2.机载三维异构阵是由多个不同的曲面阵构成的三维复杂阵列，它具有与载机外形相一致的空气动力外形，可以有效减小载机雷达反射截面积、降低载机负荷及增大有效发射孔径，能够满足复杂战场环境下雷达对目标探测的需求。相比于传统的二维平面阵，三维异构阵是一种更加广义的构型，可以看作是多个共形阵的组合阵。因此，三维异构阵的doa估计方法可以参考共形阵的研究方法，但是目前很多研究方法如music(多信号分类)和esprit(旋转不变性子空间)方法虽然具有一定效果，但是都忽略了载体曲率带来的对天线极化方向图的影响，并且受到算法需要大量快拍数据支撑的限制。
3.近年来，随着稀疏恢复类算法的兴起，如omp(正交匹配追踪)算法、iaa(迭代自适应)算法、sbl(稀疏贝叶斯学习)算法等被应用到阵列的doa估计中，通过构建压缩感知模型，使得少量快拍的doa估计成为可能。但是，由于机载三维异构阵的复杂性，线阵平面阵的稀疏doa估计方法无法适用，仍需要对相应模型进行改进。


技术实现要素：

4.针对现有技术存在的问题，本发明的目的在于提供一种基于稀疏贝叶斯学习的机载三维异构阵doa估计方法，将欧拉旋转变换引入到机载三维异构阵的方向图求解中，得到更为通用的机载三维异构阵阵列流形模型，从而构建出适用于机载三维异构阵的doa估计模型。
5.为了达到上述目的，本发明采用以下技术方案予以实现。
6.基于稀疏贝叶斯学习的机载三维异构阵doa估计方法，包括以下步骤：
7.步骤1，建立机载三维异构阵的阵列模型；
8.步骤2，利用欧拉旋转变换，对机载三维异构阵的每个阵元进行局部直角坐标系到全局直角坐标系的转换，得到机载三维异构阵的全局天线方向图；
9.步骤3，获取机载三维异构阵的接收信号，并建立对应的过完备字典矩阵；通过所述过完备字典矩阵将目标的角度估计问题转化为入射信号矩阵的稀疏恢复问题，并建立对应的压缩感知模型；
10.步骤4，利用稀疏贝叶斯学习算法求解所述压缩感知模型，得到机载三维异构阵的doa估计结果。
11.与现有技术相比，本发明的有益效果为：
12.(1)本发明方法考虑了机载三维异构阵的载体曲率对方向图综合的影响，doa估计模型更加符合实际情况。
13.(2)本发明方法可以实现少量快拍情况下机载三维异构阵的doa估计；且本发明方
法无需目标数量的先验信息。
14.(3)本发明方法首次将压缩感知模型以及稀疏贝叶斯方法引进到机载三维异构阵的doa估计上面，对后续学者的相关研究具有一定的指引作用。
附图说明
15.下面结合附图和具体实施例对本发明做进一步详细说明。
16.图1是本发明提供一种基于稀疏贝叶斯学习的机载三维异构阵的doa估计方法的流程示意图；
17.图2是本发明实施例提供的机载三维异构阵的几何特性示意图
18.图3是本发明实施例提供的机载三维异构阵的阵元布置示意图；
19.图4是本发明实施例提供的机载三维异构阵的阵元坐标图；
20.图5是本发明实施例提供的异构阵三轴欧拉旋转变换过程示意图；
21.图6本发明实施例提供的圆柱阵和圆锥阵的欧拉旋转变换过程示意图；其中，6(a)对应圆柱阵，6(b)对应圆锥阵；
22.图7本发明实施例的单次快拍的doa估计结果图；其中，7(a)对应目标1的doa估计结果，7(b)对应目标2的doa估计结果，7(c)对应目标3的doa估计结果。
具体实施方式
23.下面将结合实施例对本发明的实施方案进行详细描述，但是本领域的技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限制本发明的范围。
24.参考图1，本发明提供的一种基于稀疏贝叶斯学习的机载三维异构阵doa估计方法，包括以下步骤：
25.步骤1，建立机载三维异构阵的阵列模型；
26.参考图2，机载三维异构阵由四个部分的共形阵组成，分别为机头的圆锥阵、机身圆柱阵、机身圆台阵和机翼阵。下面以其中的圆柱共形阵和圆锥共形阵组成的三维异构阵为例，进行本发明方法的介绍。载机沿y轴正方向飞行，速度矢量为v,飞行高度为h，异构阵的几何关系图如图3所示，阵元坐标图如图4所示。
27.参考图3，机载圆柱阵由m1个圆环阵组成，每个圆环阵由关于圆心对称的n1个天线阵元组成。m1个圆环阵彼此平行，阵元间距和圆环阵间距都为d,且共同垂直于y轴。每个圆环阵半径为r，第1个圆环阵位于xoz平面。以圆环阵上顶部阵元为第1个阵元，阵元序号为逆时针排序，则圆柱阵第m1个阵列圆环上第n1个阵元的位置矢量为：
[0028][0029]
圆柱阵对应的阵元位置坐标为：
[0030]
[0031]
圆柱阵整体的位置矢量：
[0032][0033]
与圆柱阵相比，圆锥阵由不同半径、不同阵元个数的m2个圆环阵组成。将第m2(m2＝1,2,
…
,m2)个圆环阵的半径记为阵元个数记为底层半径与阵元个数与圆柱阵一致，圆柱阵半径和阵元数的表达式为：
[0034][0035][0036]
第1个圆环阵离原点距离为d2，同样以圆环阵上顶部阵元为第1个阵元，阵元序号为逆时针排序，则圆锥阵第m2个阵列圆环上第n2个阵元的位置矢量为：
[0037][0038]
其中，n2的最大取值为随着圆环阵序号m2变换。
[0039]
圆锥阵对应的阵元位置坐标为：
[0040][0041]
则圆锥阵整体的位置矢量：
[0042][0043]
则整个异构阵的坐标为：
[0044]
p＝[p1，p2]
[0045]
其中，p∈c3×
n
，n为几个共形阵的阵元数总和，对应的第n个阵元的位置矢量为：
[0046][0047]
其中，n∈1,2,
…
,n。
[0048]
如果再加入其他共形阵，以列向量形式将坐标进行合并即可。
[0049]
步骤2，利用欧拉旋转变换，对机载三维异构阵的每个阵元进行局部直角坐标系到全局直角坐标系的转换，得到机载三维异构阵的全局天线方向图；
[0050]
由于三维异构阵面的非线性特性，即使每个阵元的辐射方向图具有相同极化，在阵列的全局坐标系下，由于载体曲面的弯曲，载体表面各阵元辐射方向图具有不同极化，亦即产生了交叉极化。根据交叉极化定义的极化基和θ来计算极化阵列的辐射方向图。因为考虑极化因素，阵列天线在方向的场强方向图可以表示为：
[0051][0052]
其中，中，和分别为在和θ方向的场分量，和分别表示表达式为
[0053][0054][0055]
其中，δφ
n
表示第n个阵元的相位差，上标h表示共轭转置操作；和w
nθ
分别代表沿方向极化分量加权和沿θ方向极化分量加权，它们的比值由第n个阵元姿态决定，不能独立调节。
[0056]
接下来通过欧拉坐标旋转变换将局部坐标系下的天线方向图变换到全局直角坐标系下的天线方向图。
[0057]
欧拉旋转变换是一种实现局部坐标和全局坐标之间变换的有效方向。不同坐标系、不同共形阵对应的欧拉旋转矩阵不同，本发明提到的异构阵要实现坐标变换需要进行三次欧拉旋转变换。下面说明一种变换情况，如图5所示，其旋转轴依次为z轴，y轴，x轴，欧拉旋转角分别为d
n
、e
n
和f
n
，三次欧拉旋转矩阵表达式分别为：
[0058]
沿z轴顺时针旋转：
[0059][0060]
沿y轴顺时针旋转：
[0061][0062]
沿x轴顺时针旋转：
[0063][0064]
则最终的变换矩阵为：r＝r
x
r
y
r
z
[0065]
局部坐标(x,y,z)和全局坐标之间的变换关系为：
[0066][0067]
不同的坐标系，不同共形阵的旋转情况都不一样。
[0068]
如图6所示，以圆柱阵和圆锥阵为例，圆柱阵只需要进行一个沿y轴的旋转，如图6
(a)所示，则第m个阵元对应的旋转角分别为：
[0069]
d
n
＝0
[0070]
e
n
＝π
‑
θ0+(m
‑
1)β
[0071]
f
n
＝0；
[0072]
如图6(b)所示，圆锥阵需要先进行一个y轴的旋转，然后沿x轴旋转，则第m个阵元对应的旋转角分别为：
[0073]
d
n
＝0
[0074]
e
n
＝π
‑
θ0+(m
‑
1)β
[0075]
f
n
＝θ0；
[0076]
其中，θ0为圆锥阵的锥角，m表示阵元序号，β表示相邻阵元之间的夹角；
[0077]
首先需要知道阵元的局部方向图，然后进行旋转变换得到全局的方向图。阵元局部球坐标系下的方向图跟具体天线有关，一般选用圆形微带天线作为共形阵的阵元天线，其局部球坐标系下的第n个阵元的天线方向图表示为：
[0078][0079][0080][0081]
其中，和分别为局部球坐标下和θ
l
方向上的单位向量；和分别为第n个阵元局部方向图在和θ
l
方向的场分量；j2()和j0()分别为2阶和0阶bessel函数；为入射信号方向在局部坐标系下的坐标，其与全局球坐标系下的坐标的转换关系为实现步骤为：
[0082]
(a1)将对应的全局球坐标转换为全局直角坐标系下(x,y,z)的形式，转换关系即：
[0083][0084]
(b1)利用欧拉旋转变换将全局直角坐标系下的入射方向(x,y,z)变换为局部直角坐标系下的入射方向(x
l
,y
l
,z
l
)，转换关系为：
[0085][0086]
(c1)由于阵元的辐射方向图一般由方位角和俯仰角所确定，因此需要将局部直角坐标(x
l
,y
l
,z
l
)转换为局部球坐标转换关系：
[0087][0088]
θ
l
＝arccosz
l
[0089]
将局部球坐标代入之前的微带天线的方向图函数就可得到局部坐标系下的阵元辐射方向图，然后再用上述欧拉变换的逆变换转换成全局球坐标系方向图，具体步骤如下：
[0090]
(a2)利用球坐标系和直角坐标系的转换关系，得到天线方向图在局部直角坐标系下的分量，转换后的阵元方向图为：
[0091][0092]
转换关系：
[0093][0094]
其中，分别为局部直角坐标下各正交方向上的单位向量。
[0095]
(b2)将在局部直角坐标系下的方向图转换为全局直角坐标系下的方向图，表示成使用欧拉旋转变换，即：
[0096][0097]
(c2)利用球坐标系和直角坐标系的转换关系将全局直角坐标系下天线方向图转换为全局球坐标系下的方向图可表示为：
[0098][0099]
其中，和是全局球坐标系下方位向和俯仰向上的单位方向向量；是全局球坐标系下方位向和俯仰向上的单位方向向量；分别为局部直角坐标下各正交方向上的单位向量。
[0100]
步骤3，获取机载三维异构阵的接收信号，并建立对应的过完备字典矩阵；通过所述过完备字典矩阵将目标的角度估计问题转化为入射信号矩阵的稀疏恢复问题，并建立对应的压缩感知模型；
[0101]
(3.1)接受信号获取：
[0102]
假设有k个远场信号s
i
以入射角到达阵列，其中θ
i
和分别表示第i个入射信号的俯仰角和方位角。入射角为的信号的单位方向矢量为
[0103]
[0104]
在得到阵元的位置矢量后，机载三维异构阵的阵元的空间相位差为：
[0105][0106]
对应的导向矢量为：
[0107][0108][0109]
其中，γ
n
为第n个阵元对单位信号的响应；g
n
为阵元方向图g
nθ
和分别为和k
θ
，分别表示为信源来波极化状态，θ
ngk
为g
n
与的夹角。
[0110]
接收信号快拍数据模型可表示为：
[0111][0112]
各向量分别定义为：
[0113]
x(t)＝[x1(t),x2(t),
…
,x
n
(t)]
t
[0114]
s(t)＝[s1(t),s2(t),
…
,s
k
(t)]
t
[0115]
n(t)＝[n1(t),n2(t),
…
,n
n
(t)]
t
[0116][0117]
其中，n为机载三维异构阵的阵元总数，t为采样时间，长度即快拍数为n
t
，为机载三维异构阵的阵列流形矩阵，为阵列入射信号，为加性高斯白噪声。
[0118]
(3.2)考虑到上述信号模型的doa估计是一个二维非线性问题，本发明假设信号快拍数为n
t
＝1，将阵元角度问题转变为稀疏表达问题，通过构建过完备字典矩阵建立压缩感知模型，通过稀疏恢复算法估计方位角和俯仰角。
[0119]
将方位角度
‑
俯仰角度组成的二维空间域沿方位角轴和俯仰角轴划分为和n
θ
个均匀格子，则对应的所有网格点的导向矢量组成的完备矩阵为：
[0120][0121]
均匀量化后的快拍数据模型可表示为：
[0122][0123]
其中，γ
i',j'
为信号方向为的接收信号幅度；为高斯噪声；入射信号矩阵为
[0124]
(3.3)在上述量化后的模型中，通过过完备字典矩阵φ使得角度估计问题转化为γ的稀疏恢复问题。由于信号源个数远小于空间域的量化的网格点数，γ在空间域是稀疏的，可以构建成l1范数最小化的求解问题，其目标函数为：
[0125][0126]
步骤4，利用稀疏贝叶斯学习算法求解所述压缩感知模型，得到机载三维异构阵的doa估计结果。
[0127]
稀疏贝叶斯学习(sbl，sparse bayesian learning)是在贝叶斯原理的基础上，加入最大后验概率估计和最大似然估计的一种稀疏恢复算法。稀疏贝叶斯在无需稀疏度先验信息和噪声干扰阈值信息的条件下，可以较为准确得到稀疏信号的估计值，本发明采用该算法来进行l1范数的求解。直接给出sbl步骤：
[0128]
(4.1)为了便于迭代计算开始，需要初始化参数，该初始化过程并不会影响最后的计算结果。初始化噪声水平σ2＝0.001，初始化超参数α
i
＝0.001，设置收敛条件δ＝0.001；
[0129]
(4.2)根据下式计算稀疏信号后验概率的方差σ和均值μ，并令μ
old
＝μ；
[0130]
∑＝(σ
‑2φ
t
φ+λ)
[0131]
μ＝σ
‑2σφ
t
x
[0132]
其中，
[0133]
(4.3)根据最大似然估计方法计算超参和噪声水平(σ2)
new
的估计值：
[0134][0135][0136]
其中，γ
i
＝1
‑
α
i
∑
ii
，∑
ii
表示∑的对角线元素组成的向量，
[0137]
(4.4)根据步骤(4.2)计算的和(σ2)
new
重新计算方差和均值，并记作σ
new
和μ
new
；
[0138]
(4.5)计算dμ＝abs(max(μ
new
‑
μ
old
))，并判断dμ是否小于等于δ，若是，则停止迭代，μ
new
即为重构信号γ；否则令μ
old
＝μ
new
，返回步骤(4.3)。
[0139]
仿真实验
[0140]
以下通过仿真实验验证本发明的有效性，以圆柱
‑
圆锥阵的组成的异构阵为例，仿真参数如表所示：
[0141]
表1系统仿真参数
[0142][0143][0144]
目标1、2和3对应的估计结果如图7(a)、7(b)和7(c)所示，图7(a)估计的目标1的角度值(60
°
,100
°
)，(70
°
,110
°
)，幅度值为0.7db和0.1db；7(b)估计的角度值为(64.76
°
,
105.2
°
)，(70
°
,110
°
)，幅度值为0.9999db，1db；7(c)估计的角度值为(64.76
°
,105.2
°
)，(74.76
°
,115.2
°
),幅度值为0.75db，0.9db，可以看出，三个目标幅度值的估计基本正确，目标2和目标3角度估计有一个小于网格精度0.476误差，这是由于网格划分的精度不够而产生的误差，在误差允许范围。综上，本发明方法能很好的进行机载三维异构阵的doa估计。
[0145]
虽然，本说明书中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。