一种基于改进的离网格稀疏贝叶斯学习的DOA估计方法

一种基于改进的离网格稀疏贝叶斯学习的doa估计方法
技术领域
1.本发明涉及雷达定位技术领域，对接收信号进行波达方向估计。

背景技术：

2.接近现有技术：
3.1.通过对搜索范围均匀划分网格，从贝叶斯的角度出发，yang z基于离网格模型提出了离网格稀疏贝叶斯推理(off-grid sparse bayesian inference,ogsbi)，并开发了一种迭代算法。这种新的方法假设所有快照的信号为拉普拉斯先验，以来利用不同快照之间的联合稀疏性，适用于单快照和多快照情况的回波方向定位问题。现有的方法仍然是一个主要问题，它们的性能依赖于精度和计算工作量之间的权衡。当采用粗网格时，建模误差较大，而当采用密集采样网格时，涉及的计算量较大，不利于实际应用。
4.2.为了扩大阵列的阵元孔径，dai j提出了一个嵌套阵列的接收模型，该模型将噪声方差作为一个感兴趣的未知信号的一部分，然后通过稀疏贝叶斯学习(sparse bayesian learning,sbl)的方法对网格点进行迭代，消除了由离格间隙引起的模型误差。此外，yang j用二阶泰勒展开式来代替一阶泰勒展开式，以减轻来自网格的间隙误差，但它会增加计算量。虽然引入嵌套阵列结构大大提高了阵列可探测源数，且突破了传统天线阵阵元间距半波长的限制，使得天线孔径得到极大的扩展，能够获得角度估计性能的提升，但是其缺点是估计精度有限。
5.3.非圆(non-circular，nc)信号，如二进制相位键控(binary phase shift keying，bpsk)、脉冲振幅调制(pulse amplitude modulation，pam)和幅度移位键控(amplitude shift keying，ask)等具有非零椭圆协方差矩阵属性，它可以扩大阵列孔径，增加自由度和估计精度，广泛应用于现代数字调制方案，如通信或卫星系统。qiu t将具有数控信号的接收模型转换为实值形式，并通过重构实值稀疏模型，提出了一种基于sbl的脉冲噪声实值信号处理方法。此外，zheng r提出了利用sbl方法解决圆信号和非圆信号共存时的离网格doa估计问题，该方法可以自动识别候选角度网格中的源，并匹配协方差向量和伪协方差向量中感兴趣的方向信息。这些方法采用非圆信号的椭圆协方差旋转不变性，扩大了接收数据信息，但是由于采用的接收阵列为均匀线阵，精度有待提高，尚需改进才能运用到实际场景中。

技术实现要素：

6.针对传统的均匀线阵限制了阵列孔径和波达方向估计精度的问题，本发明提出了一种嵌套阵列下基于改进的离网格稀疏贝叶斯学习的doa估计方法。为了增加阵列孔径，利用非圆形信号的协方差矩阵不具有旋转不变性，构造了阵列接收数据的扩展矩阵，充分利用了差阵阵列与和阵阵列。同时，将噪声作为信号的一部分，采用稀疏贝叶斯学习方法对网格点进行迭代，每次迭代时对网格点进行更新，以消除离网格间隙造成的模型误差。
7.技术方案
8.一种嵌套阵列下非圆信号的基于离网格稀疏贝叶斯推理的doa估计方法，具体包括如下步骤：
9.s1：通过嵌套阵列结构的阵列天线接收窄带远场波非圆信号，得到接收信号信息为y0。根据所述接收信号信息，利用非圆信号的特性，计算得到协方差矩阵
10.s2：对进行向量化处理，得到虚拟阵列接收信号引入一个行交换矩阵对向量化后的接收信号进行行交换得到量测数据z。
11.s3：收到的量测数据z可以用稀疏近似方法求解，采用均匀分布的网格划分空间，量测数据z的表达式可以改写成过完备基的形式。
12.s4：对z中两个和阵和两个相同的差分阵的虚拟阵元的位置从小到大进行排序，得到两个差分阵以及两个和阵的分量信息，可以把量测数据z分成四部分。
13.s5：计算权值矩阵w，用权值矩阵对向量化协方差矩阵进行归一化，再经过去冗余矩阵f得到标准的sbl形式。
14.s6：计算s7中所需的参数，先对去冗余后的量测数据假设隐藏变量d，并计算其均值μ和协方差σ，且对构成sbl条件的参数δ和β进行更新。如果满足最终的收敛准则或达到最大的迭代次数，则截断迭代过程。
15.s7：更新网格最后，根据网格峰值估计目标信源。
16.技术效果
17.本发明将嵌套阵列与非圆信号doa估计问题相结合，利用了非圆信号的非圆特性，得到了更高的自由度。本发明采用了嵌套阵列来接收信源发出的非圆信号，利用了信号的非圆特性，构造了阵列接收数据的扩展矩阵，充分利用了差阵阵列和和阵阵列，提高了doa估计精度。同时，将噪声作为信号的一部分，采用稀疏贝叶斯学习方法对网格点进行迭代，迭代后对网格点进行更新，以消除离网格间隙造成的模型误差。同传统均匀线阵阵列接收信源发出的圆信号相比，本发明方法提高了估计的精确度。
附图说明
18.图1本发明方法流程示意图
19.图2本发明的嵌套阵列结构示意图
20.图3本发明的嵌套线阵虚拟阵列结构示意图
21.图4实施例信号个数为2时不同方法的振幅对比图
22.图5实施例信号个数为2时在不同信噪比下误差对比图
23.图6实施例信号个数为2时在不同快照下误差对比图
24.图7实施例信号个数为2时在不同阵列元组个数下误差对比图
25.图8实施例信号个数为13时不同方法的振幅对比图
26.图9实施例信号个数为13时在不同信噪比下误差对比图
27.图10实施例信号个数为13时在不同快照下误差对比图
具体实施方式
28.以下结合附图和实施例进一步介绍本发明技术方案。
29.本发明方法流程示意图如图1所示，其中过程可以按照如下方式介绍：
30.s1：利用非圆信号的特性，计算接收信号信息y0和协方差矩阵
31.考虑k个窄带远场源θk,k＝1,2,
…
k撞击一个嵌套阵列，如图2所示，该嵌套阵列是由内部和外部分别有m和n个固定间距的均匀线性阵列传感器组成，其中d＝λ/2是阵元间距，λ是信号载波的波长。传感器的位置表示为
[0032][0033]
li为i第个阵元的位置，接收信号信息可以表示为
[0034][0035]
所述为t时刻接收到的数据，引入非圆相位后，扩展的导向矢量可以表示为矢量可以表示为是复数空间，第k个扩展的导向向量为方向矢量可以表示为θk为第k个信号的doa，为符合高斯分布的增广脉冲噪声项，s(t)＝[s1(t),
…
,sk(t),
…
sk(t)]
t
为k个信号组成的矢量，其中(
·
)
*
为共轭，(
·
)
t
为共轭转置，j为虚数单位，sk(t)表示第k个信号，t表示时域快拍。利用非圆信号的特性，计算得到协方差矩阵为
[0036][0037]
s2：计算虚拟阵列接收信号引入一个行交换矩阵后得到量测数据z
[0038]
接收数据求完协方差矩阵后进行向量化处理，考虑其渐进误差t表示总的快拍数，则向量化后的接收信号为
[0039][0040]
其中其中是第k个虚拟阵列导向矢量，是一个单快拍的多个信号向量，是除第l个位置1以外的所有零的列向量,其中是除第l个位置1以外的所有零的列向量,其中是实数空间，
⊙
和分别表示khatri-rao积和kronecker积。引入的一个行交换矩阵满足和可以表示为
[0041][0042]
其中0表示零矩阵，i表示单位矩阵，令e0＝e
*
⊙
e，对向量化后的接收信号进行行交换得到量测数据z，即
[0043][0044]
其中表示虚拟差分阵列和和阵阵列的方向矩阵，且满足和
[0045]
s3：将量测数据z的表达式可以改写成过完备基的形式
[0046]
接收数据z可以采用稀疏逼近的方式表示，将平面空间均匀划分为i个采样网格数，满足i＞＞k，则采样空间可以表示为经行交换矩阵变化后的接收数据可重新表述为
[0047][0048]
其中表示新的导向矢量，是p的零扩展矩阵，它的非零元素与真实的角度位置一致。理想情况下，估计的角度被认为在网格上，但是在实际场景中，这种现象几乎不存在。本发明方法采用一种基于一阶泰勒展式线性逼近的方法解决这个问题，假设是真实角度临近的网格点，导向向量可以线性化表示为
[0049][0050]
其中是对于的一阶求导，的一阶求导，而β＝[β1,β2,
…
βi]
t
是除了中第ik个元素的零向量，可以得到接收数据z可以采用稀疏逼近的方式表示
[0051][0052]
其中是一个虚拟阵列方向矩阵，且1
′n＝j1n和
[0053]
s4：利用差分阵以及和阵的分量信息，把量测数据z分成四部分
[0054]
差阵和阵信息等价于对阵列转向向量中的位置信息的和差运算，差分阵列为对应的和阵阵列中存在孔洞，所以本发明方法只选择连续的虚拟数组元素，即和例如图3所示为一个嵌套线阵虚拟阵
列结构示意图，其中m＝3,n＝3，则差分数组为《-11,11》d，连续负的数组为《-14,0》d，连续正的数组为《0,14》d，ula可以表示为《-14,0》d∪《-11,11》d∪《0,14》d。通过引入行交换矩阵，接收数据z可以被分成4部分，分别于差阵、负和阵、正和阵、差阵相对应，即
[0055][0056]
其中其中为零扩展矩阵和为渐进误差。
[0057]
s5：用权值矩阵对向量化协方差矩阵进行归一化，再去冗余操作
[0058]
先用权值矩阵对向量化协方差矩阵进行归一化。由于渐进误差符合一个高斯分布其中j行交换矩阵，t表示总的快拍数，为了使误差服从一个标准的稀疏贝叶斯格式，可以推导出
[0059][0060]
接收数据z由差分和阵被分成4部分中包含冗余信息，因此需要对接收数据进行去冗余操作，若去除z1,z2,z3,z4上的冗余项，这取决于选择矩阵f1,f2,f3,f4。根据虚拟数组的位置位置和本发明方法定义的权值函数如下
[0061][0062]
其中card(μ)表示μ的个数。令dv＝mn+m-1,sv＝mn+m+n-1,m和n为固定间距的均匀线性阵列的传感器个数，则对应的选择矩阵分别为
[0063]
[0064]
通过去除接收数据z上的冗余项，可以获得
[0065][0066]
所有有效的接收信息都可以表示为去冗余后用稀疏逼近的方式表示为
[0067][0068]
这里d表示稀疏信号，是一个隐藏变量，同时，去除冗余后的符合一个高斯分布，即
[0069][0070]
s6：计算隐藏变量的均值和协方差，并对构成sbl条件的参数进行更新
[0071]
计算隐藏变量d的均值μ和协方差σ，一个的典型sbl处理是d为每个元素分配一个具有方差的非平稳高斯先验分布δi，假设该超参数λ＝diag(δ)，且δ＝[δ1,
…
,δi,
…
,δ
i+1
]
t
，i表示稀疏网格，δ可以建模为一个伽马分布这里υ一个小的正常数，满足其中q＝fwf
t
，f为选择矩阵，w为权值矩阵。d是一个隐藏变量，它的后验概率密度为其中均值和方差为
[0072][0073]
σ＝(ψh(β)f
t
q-1
fψ(β)+λ-1
)-1
ꢀꢀꢀ
(18)
[0074]
进一步的，可以应用期望最大化(expectation-maximization，em)算法进行递归计算，直到达到规定的精度，在此过程中，需要达到的下界，即
[0075][0076]
δi的超参数更新可以简化为
[0077][0078]
[
·
]
ii
表示矩阵的第(i,i)个元素。此外，有其中p和v有特定的表达式，β的更新可以简化为
[0079]
β＝p-1vꢀꢀꢀꢀ
(21)
[0080]
如果p是可逆的，上式成立；否则，
[0081]
s7：更新网格，估计目标信源
[0082]
为了减少网格间隔的误差，需要更新网格本发明方法定义要更新的网格如下
[0083][0084]
其中是真实角度临近的网格点，为对应网格中构成sbl条件的参数。该定义可以提供一个更好的逼近真实的转向矩阵ψ(β)的方法，我们用新的网格更新矩阵ψ(β)，使下一次迭代更接近真实的值，即最后，根据网格峰值估计目标信源。
[0085]
实施例
[0086]
分别对信号个数为2和13的非圆信号进行仿真，仿真中采样频率均设置为300khz，两种仿真的具体参数设置如下：
[0087]
信号个数为2时，考虑10个物理阵元，即m＝5,n＝5，非圆相位设定为假设两个等功率独立的窄带非圆形信号分别从方向-50
°
+η和方向30
°
+η撞击嵌套阵列，其中η为在的范围10*randn(-1,1)内选择的随机值，以消除预定义的doa集中可能包含的先验信息。性能估计标准为联合均方根误差(root mean square error,rmse)定义为：
[0088][0089]
其中为第i次蒙特卡罗过程第k个信号doa的精确估计值，θk表示第k个信号doa的真实值，mc表示蒙特卡罗试验次数。
[0090]
假设两个信噪比为0db的信号位于一个随机生成的doa中，快照数为t＝200，角分辨率为2度，所有算法的归一化振幅如图4所示，可以看出，l
1-svd和ss-music算法的源信号振幅最小，与真实的doa有偏差。虽然ogsbi-ula可以估计源的位置，但精度还需要提高。显然，ogsbi-na和提出的方法都成功地到了两个源；然而，ogsbi-na信号峰值不如提出的算法。这表明，本发明专利的方法在检测方面优于现有的算法。
[0091]
假设每个信号的信噪比从-10db到10db，快照数量固定在200，mc＝200，角度分辨率为2度，给出五种方法的rmse，包括l
1-svd、ss-music、ogsbi-ula、ogsbi-na和提出的算法的对比结果如图5所示。随着信噪比的提高，本发明专利的方法性能逐渐提高。在接收信号是嵌套阵列的情况下，由于使用了由非循环信号形成的共轭信息，该算法比ogsbi-na方法具有更高的性能。
[0092]
假设信噪比保持在0db，并将快照的数量从50增加到800，运行蒙特卡罗200次。五种方法的rmse与快照数之间的关系如图6所示。结果表明，ogsbi-na和该算法只需300个快照就能获得超分辨率，且该方法的rmse值小于ogsbi-na。随着快照数量的增加，本发明专利方法的rmse值最低。证明了非圆信号的共轭特性可以改进估计精度。
[0093]
为了验证传感器的数量与该算法的估计性能之间的关系，本发明选择了不同数量的传感器进行比较，结果在图7中给出。当mc＝200和t＝200时，本发明专利方法的精确度随着阵列元素数的增加而逐渐增加。同时，随着信噪比的增加，本发明专利方法的rmse值的增加呈下降趋势。
[0094]
信号个数为13时，考虑10个物理阵元，即m＝5,n＝5，13个非圆相位设定为假设13个等功率独立的窄带非圆信号分别从方向《-60
°
:10:60
°
》+η撞击嵌套阵列，其中η
为在的范围10*randn(-1,1)内选择的随机值，以消除预定义的doa集中可能包含的先验信息。
[0095]
假设13个信噪比为20db的信号位于一个随机生成的doa中，快照数为t＝200，角分辨率为2度，所有算法的归一化振幅如图8所示，可以看出，本发明专利方法和ogsbi-na算法都成功地定位了所有13个源。然而，ogsbi-na方法的一些信号峰的振幅较低，这表明本发明专利方法优于现有的算法。
[0096]
假设每个信号的信噪比从10db到20db，快照数量固定在200，mc＝200，网格分辨率为2度，给出三种方法的rmse，包括ss-music、ogsbi-na和本发明专利方法的对比结果如图9所示。由于使用了非循环信号形成的共轭信息，随着信噪比的提高，本发明专利的方法性能逐渐提高。
[0097]
假设信噪比保持在10db，并将快照的数量从100增加到1200，运行蒙特卡罗200次。三种方法的rmse与快照数量的关系如图10所示。随着快照数量的增加，ss-music对doa估计的精度略有提高，但由于估计源的数量较多，误差仍然相对较大。当快照t＝200时，本发明专利方法的rmse与ogsbi-na的rmse接近。随着快照数量的增加，本发明专利方法的rmse值最低。结果表明，非圆信号的共轭性质可以提高估计的精度。
[0098]
通过对欠定和非欠定的信号估计仿真中可以看出，本发明提出的一种基于非循环信号嵌套阵列的非网格稀疏贝叶斯推理的doa估计算法可以有效提高doa估计精度，该方法利用非圆形信号的统计特性具有旋转不变性，对所接收到的数据进行扩展，基于此优点，构造差分共阵列和和共阵列，大大提高了阵列的孔径。

技术特征：

1.一种基于改进的离网格稀疏贝叶斯学习的doa估计方法，具体包括如下步骤：s1：通过嵌套阵列结构的阵列天线接收窄带远场波非圆信号，得到接收信号信息为y0；根据所述接收信号信息，利用非圆信号的特性，计算得到协方差矩阵s2：对进行向量化处理，得到虚拟阵列接收信号引入一个行交换矩阵对向量化后的接收信号进行行交换得到量测数据z；s3：收到的量测数据z用稀疏近似方法求解，采用均匀分布的网格划分空间，量测数据z的表达式改写成过完备基的形式；s4：对z中两个和阵和两个相同的差分阵的虚拟阵元的位置从小到大进行排序，得到两个差分阵以及两个和阵的分量信息，把量测数据z分成四部分；s5：计算权值矩阵w，用权值矩阵对向量化协方差矩阵进行归一化，再经过去冗余矩阵f得到标准的sbl形式；s6：计算s7中所需的参数，先对去冗余后的量测数据假设隐藏变量d，并计算其均值μ和协方差σ，且对构成sbl条件的参数δ和β进行更新；如果满足最终的收敛准则或达到最大的迭代次数，则截断迭代过程；s7：更新网格最后，根据网格峰值估计目标信源。2.如权利要求1所述的方法，其特征在于：所述s1，具体为：考虑k个窄带远场源θ
k
,k＝1,2,
…
,k撞击一个嵌套阵列，该嵌套阵列是由内部和外部分别有m和n个固定间距的均匀线性阵列传感器组成，其中d＝λ/2是阵元间距，λ是信号载波的波长；传感器的位置表示为l
i
为i第个阵元的位置，接收信号信息表示为所述为t时刻接收到的数据，引入非圆相位后，扩展的导向矢量表示为表示为是复数空间，第k个扩展的导向向量为方向矢量表示为θ
k
为第k个信号的doa，为符合高斯分布的增广脉冲噪声项，s(t)＝[s1(t),
…
,s
k
(t),
…
s
k
(t)]
t
为k个信号组成的矢量，其中(
·
)
*
为共轭，(
·
)
t
为共轭转置，j为虚数单位，s
k
(t)表示第k个信号，t表示时域快拍；利用非圆信号的特性，计算得到协方差矩阵为3.如权利要求1所述的方法，其特征在于：
所述s2，具体为：接收数据求完协方差矩阵后进行向量化处理，考虑其渐进误差t表示总的快拍数，表示复高斯分布，则向量化后的接收信号为其中其中是第k个虚拟阵列导向矢量，是一个单快拍的多个信号向量，是除第l个位置1以外的所有零的列向量,其中是除第l个位置1以外的所有零的列向量,其中是实数空间，
⊙
和分别表示khatri-rao积和kronecker积；引入的一个行交换矩阵满足和表示为其中0表示零矩阵，i表示单位矩阵，令e0＝e
*
⊙
e，对向量化后的接收信号进行行交换得到量测数据z，即其中表示虚拟差分阵列和和阵阵列的方向矩阵，且满足和4.如权利要求1所述的方法，其特征在于：所述s3，具体为：接收数据z采用稀疏逼近的方式表示，将平面空间均匀划分为i个采样网格数，满足i＞＞k，则采样空间表示为经行交换矩阵变化后的接收数据可重新表述为其中表示新的导向矢量，是p的零扩展矩阵，它的非零元素与真实的角度位置一致；采用一种基于一阶泰勒展式线性逼近的方法，假设是真实角度临近的网格点，导向向量线性化表示为其中是对于的一阶求导，
而β＝[β1,β2,
…
β
i
]
t
是k＝1,2,
…
,k除了中第i
k
个元素的零向量，得到接收数据z采用稀疏逼近的方式表示其中是一个虚拟阵列方向矩阵，且1
′
n
＝j1
n
和5.如权利要求1所述的方法，其特征在于：所述s4，具体为：差阵和阵信息等价于对阵列转向向量中的位置信息的和差运算，差分阵列为对应的和阵阵列中存在孔洞，选择连续的虚拟数组元素，即和则新形成的连续均匀线阵表示为其中d是阵元间距；通过引入行交换矩阵，接收数据z被分成4部分，分别于差阵、负和阵、正和阵、差阵相对应，即其中为零扩展矩阵和为渐进误差。6.如权利要求1所述的方法，其特征在于：所述s5，具体为：先用权值矩阵对向量化协方差矩阵进行归一化；渐进误差符合一个高斯分布其中j行交换矩阵，t表示总的快拍数，为了使误差服从标准的稀疏贝叶斯格式，推导出接收数据z由差分和阵被分成4部分中包含冗余信息，对接收数据进行去冗余操作，若去除z1,z2,z3,z4上的冗余项，这取决于选择矩阵f1,f2,f3,f4；根据虚拟数组的位置和定义的权值函数如下
其中card(μ)表示μ的个数；令d
v
＝mn+m-1,s
v
＝mn+m+n-1,m和n为固定间距的均匀线性阵列的传感器个数，则对应的选择矩阵分别为通过去除接收数据z上的冗余项，获得所有有效的接收信息都表示为去冗余后用稀疏逼近的方式表示为这里d表示稀疏信号，是一个隐藏变量，同时，去除冗余后的符合一个高斯分布，即7.如权利要求1所述的方法，其特征在于：所述s6，具体为：计算隐藏变量d的均值μ和协方差σ，d为每个元素分配一个具有方差的非平稳高斯先验分布δ
i
，假设该超参数λ＝diag(δ)，且δ＝[δ1,
…
,δ
i
,
…
,δ
i+1
]
t
，i表示稀疏网格，δ建模为一个伽马分布这里υ一个小的正常数，满足其中q＝fwf
t
，f为选择矩阵，w为权值矩阵；d是隐藏变量，它的后验概率密度为其中均值和方差为
σ＝(ψ
h
(β)f
t
q-1
fψ(β)+λ-1
)-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)进一步的，应用期望最大化(expectation-maximization，em)算法进行递归计算，直到达到规定的精度，在此过程中，需要达到的下界，即δ
i
的超参数更新简化为[
·
]
ii
表示矩阵的第(i,i)个元素；此外，有其中p和v有特定的表达式，β的更新简化为β＝p-1
v(21)如果p是可逆的，上式成立；否则，8.如权利要求1所述的方法，其特征在于：所述s7，具体为：定义要更新的网格如下其中是真实角度临近的网格点，为对应网格中构成sbl条件的参数；该定义提供逼近真实的转向矩阵ψ(β)，用新的网格更新矩阵ψ(β)，使下一次迭代更接近真实的值，即最后，根据网格峰值估计目标信源。

技术总结

针对传统的均匀线阵限制了阵列孔径和波达方向估计精度的问题，本发明提出了一种基于改进的离网格稀疏贝叶斯学习的DOA估计方法。为了增加阵列孔径，利用非圆形信号的协方差矩阵不具有旋转不变性，构造了阵列接收数据的扩展矩阵，充分利用了差阵阵列与和阵阵列。同时，将噪声作为信号的一部分，采用稀疏贝叶斯学习方法对网格点进行迭代，每次迭代时对网格点进行更新，以消除离网格间隙造成的模型误差。以消除离网格间隙造成的模型误差。以消除离网格间隙造成的模型误差。