嵌套阵中基于非高斯信源的求根谱峰搜索传播算子方法与流程

1.本发明涉及一种嵌套阵中基于非高斯信源的求根谱峰搜索传播算子方法，属于阵列信号处理领域。


背景技术：

2.信号来向及其空间分布通常称为空间谱，空间谱估计又称接收信号波达方向(direction of arrival,doa)估计是阵列信号处理一个重要的研究方向。其基本问题就确定同时处于空间中的多个目标信号的空间位置方向信息。传统的波达方向(doa)参数估计方法等大多利用了信号的二阶统计特性。但在实际应用中，我们遇到的通常是非高斯信号，其一阶、二阶统计量并不能完全描述信号的统计特性，这时采用高阶累积量不仅可以获得比二阶统计量更高的性能，而且可以解决二阶统计量不能解决的很多问题。此外，高阶累积量对高斯过程具有不敏感性，理论上可抑制噪声。高阶统计量作为一种强有力的信号处理工具,己经在通信、雷达、声纳及地球物理、生物医学等领域得到了极为广泛的应用。近年来，稀疏阵列由于其分辨目标比传感器数量多的特点，在doa估计中受到越来越多的关注。常见的稀疏阵列主要包括最小冗余阵列、互质阵列和嵌套阵列。嵌套阵列能够在得到阵元位置的固定闭式解的同时，还拥有比互质阵列更大的空间自由度。这些优点使得嵌套阵列在信号参数估计中有着较大优势。当信源为非高斯信源时，我们采用高阶累积量如四阶累积量进行doa估计可实现更大虚阵孔径。
3.谱峰搜索传播算子pm算法是一种快速的子空间算法，该算法利用阵列协方差矩阵的子矩阵得到噪声子空间，无需复杂的特征值分解，但该算法仅在信噪比比较好的情况下才有较好的波达方向估计性能，且谱峰搜索仍需要较大计算量。针对计算复杂度大的谱峰搜索部分，之前的学者研究有将该部分的内容转化为通过求解多项式的根来实现。


技术实现要素：

4.本发明所要解决的技术问题是：提供一种嵌套阵中基于非高斯信源的求根谱峰搜索传播算子方法，将空间平滑算法(spatial smoothing，ss)、求根算法(root)和谱峰搜索传播算子算法(propagator method，pm)结合并应用于基于非高斯信源的嵌套阵列，该发明算法名称为基于空间平滑的求根谱峰搜索传播算子算法(ss-rootpm)。在基于非高斯信源的嵌套阵列结构中，该发明算法既能保证一定的角度估计性能，也大大降低了计算复杂度。
5.本发明为解决上述技术问题采用以下技术方案：
6.嵌套阵中基于非高斯信源的求根谱峰搜索传播算子(ss-rootpm)方法，包括以下步骤：
7.步骤1、建立嵌套阵接收信号的数学模型；
8.步骤2、构建接收信号的四阶累积量矩阵，对其进行矢量化处理，对完成矢量化后的向量进行排序、删除冗余；
9.步骤3、利用空间平滑算法，将步骤2中矢量化处理得到的向量的方向矩阵对应的
二级虚拟阵列划分为若干个相互重叠的子阵列，对全部子阵列的协方差矩阵求和并计算均值，得到空间平滑协方差矩阵；
10.步骤4、根据空间平滑协方差矩阵分块，估计传播算子并构造q矩阵；
11.步骤5、构造多项式并进行求根处理得到波达方向估计角度。
12.进一步，步骤1中所述嵌套阵由两个均匀线性子阵串联而成，其中第一级线性子阵的阵元数为m1，阵元之间的间距为d＝λ/2；第二级线性子阵的阵元数是m2，阵元之间的间距为d2＝(m1 1)d；两个线性子阵之间的间距为d＝λ/2，λ表示波长。
13.进一步，步骤1中对于整个嵌套阵，接收信号表示为：
14.x＝a(θ)s n
15.其中，表示方向矩阵，m表示嵌套阵中线性子阵的阵元数，k表示信源数，θk为第k个信源的入射角，k＝1，2，...，k，d表示嵌套阵中两个线性子阵之间的间距，λ表示波长表示均值为0的非高斯信号信源，且sk(l)表示第ll个快拍数的信源，l＝1，2，...，l，l表示快拍数，表示为均值为0、方差为σ2的加性高斯白噪声。
16.进一步，步骤3具体为：
17.利用空间平滑算法，将二级虚拟阵列划分为mb 1个相互重叠的子阵列，每一个子阵列均包含mb 1个阵元，其中，第i个子阵列的阵元位置为：
18.{(-i 1 n)d，n＝0，1，
…
，mb}；
19.式中，d表示嵌套阵中两个线性子阵之间的间距；
20.第i个子阵列的接收信号矩阵为z1的第mb 2-i行至第2mb 2-i行，并构造第i个子阵列的协方差矩阵：
[0021][0022]
式中，z1表示对矢量化处理得到的向量进行排序、删除冗余后的结果，z
1i
表示z1的第mb 2-i行至第2mb 2-i行；
[0023]
将全部mb 1个子阵列的协方差矩阵求和并计算均值，得到空间平滑协方差矩阵：
[0024][0025]
进一步，步骤4具体为：
[0026]
将空间平滑协方差矩阵进行分块：
[0027]
r＝[g，h]
[0028]
其中，g和h分别为(mb 1)
×
k和(mb 1)
×
(mb 1-k)矩阵；
[0029]
传播算子为p＝(ghg)-1gh
h，构造q矩阵：
[0030][0031]
进一步，步骤5中构造求根多项式的表达式为：
[0032][0033]
式中，且z＝exp(jw)(当z＝exp(jw)时，即多项式的根正好在单位圆上，w表示空间频率)，w表示空间频率，p(exp(jw))是一个以空间频率为w的导向矢量，且与噪声子空间是正交的；
[0034]
对上式进行求根处理，取单位圆内具有最大幅值的k个根z1，z2，...，zk的相位给出波达方向估计，即：
[0035][0036]
其中，λ为波长，d为阵元间距。
[0037]
本发明采用以上技术方案与现有技术相比，具有以下技术效果：
[0038]
1、本发明算法(ss-rootpm)提出结合空间平滑算法、求根算法和谱峰搜索传播算子算法，将计算复杂度大的谱峰角度搜索转变为求解多项式的根，进而降低了部分计算复杂度；
[0039]
2、本发明算法在计算复杂度方面远低于基于空间平滑算法的谱峰搜索传播算子(ss-pm)算法，低于基于空间平滑算法的求根多重信号分类(multiple signal classification，ss-rootmusic)算法；
[0040]
3、本发明算法在参数估计性能方面略优于ss-pm算法，低于ss-rootmusic算法。
附图说明
[0041]
图1是嵌套阵示意图；
[0042]
图2是本发明的算法流程图；
[0043]
图3是一级虚拟阵元位置及权重示意图；
[0044]
图4是二级虚拟阵元位置及权重示意图；
[0045]
图5是几种不同阵元位置对比图；
[0046]
图6是采用空间平滑算法后的虚拟阵列图；
[0047]
图7是本发明算法多次参数估计点图；
[0048]
图8是本发明算法在不同快拍数下算法doa估计性能比较图；
[0049]
图9是本发明算法在不同阵元数下算法doa估计性能比较图；
[0050]
图10是本发明算法和ss-rootmusic算法、ss-pm算法的角度估计性能在相同阵列结构和相同快拍数条件下不同信噪比的对比图；
[0051]
图11是本发明算法和ss-rootmusic算法、ss-pm算法在不同快拍数下的计算复杂度图；
[0052]
图12是本发明算法和ss-rootmusic算法、ss-pm算法在不同阵元数下的计算复杂度图。
具体实施方式
[0053]
下面结合附图以及具体实施例对本发明的技术方案做进一步的详细说明。
[0054]
假设空间有k个窄带远场不相干信源入射到此嵌套阵上，其一维波达方向为θk(k
＝1，2，
…
，k)。嵌套阵可以表示为两个均匀线阵串联而成，第一级线性子阵的阵元数为m1，其阵元之间的间距为d＝λ/2，第二级线性子阵的阵元数是m2，其阵元之间的间距为d2＝(m1 1)d，两个子阵之间的间距为d＝λ/2，其中λ表示波长。本发明中涉及的嵌套阵如图1所示，不失一般性，假设二级嵌套阵的两个子阵的阵元数均为m1＝m2＝m。首先根据阵列信号数学模型得到累积量矩阵，对累积量矩阵进行矢量化并排序、去除冗余。
[0055]
本发明算法的流程图概要如图2所示，具体实现如下：
[0056]
步骤1：建立嵌套阵接收信号的数学模型。
[0057]
对于整个嵌套阵，接收信号可以表示为：
[0058]
x＝a(θ)s n(1)
[0059]
其中，表示方向矩阵，且
[0060]
表示均值为0的非高斯信号信源，且l表示快拍数。表示为均值为0、方差为σ2的加性高斯白噪声。其中加粗大写字母表示矩阵，加粗小写字母表示矢量，(
·
)
t
表示为转置运算。
[0061]
步骤2：求出四阶累积量矩阵r4，对其进行向量化处理得到向量z。
[0062]
四阶累积量方法利用接收到的数据的四阶累积量特征来构建高阶虚拟阵列模型，扩展虚拟阵列的孔径，并实现高dof的doa估计。我们首先求解接收信号的四阶累积量矩阵
[0063]
其中
[0064][0065][0066]
b(θ)为利用四阶累积量阵列扩展后的阵列流形其中：a
*
(θ)表示共轭运算，(
·
)h表示为共轭转置运算，表示kronecker积。
[0067]
为进一步提高doa估计精度，我们采用虚拟化的方法对四阶累积量矩阵r4进行处理：
[0068]
z＝vec(r4)＝vec[b(θ)c
sbh
(θ)]＝λ(θ)p
ꢀꢀꢀ
(5)
[0069]
其中
[0070]
[0071]
且表示第k个信源的功率。vec(
·
)表示矩阵矢量化，
⊙
表示khatri-rao积。
[0072]
矢量化后的r4对应的方向矢量矩阵a(θ)可考虑由a(θ)两级差分阵列扩展而来b(θ)对应的一级连续虚拟阵元位置如图3所示，其中横坐标为location(位置)，纵坐标weight value(权重)为对应着每个横坐标的位置的阵元所占据的权重比。位置集合数学表达式可以表示为：
[0073]
l1＝{-mad，-(m
a-1)d，...，0，...，(m
a-1)d，mad}，ma＝m(m 1)-1
ꢀꢀꢀꢀ
(7)
[0074]
a(θ)对应的二级虚拟阵元位置如图4所示，其闭式解可以表示为：
[0075]
l2＝{-mbd，-(m
b-1)d，...，0，...，(m
b-1)d，mbd}，mb＝2ma＝2(m(m 1)-1)
ꢀꢀꢀꢀ
(8)
[0076]
显然，二级虚拟阵元的位置也是完全连续的。同时我们给出了几种不同阵元的位置对比图，分别为phisical array(物理阵元)、coarray(共阵)、foc coarray(四阶累积量的共阵)，如图5所示，其中横坐标location(位置)表示阵元所处的位置。
[0077]
根据l2中二级虚拟阵元的位置对z进行排序并去除冗余阵元构造z1，
[0078]
z1＝λ1p
ꢀꢀꢀꢀ
(9)
[0079]
其中，是λ(θ)删除冗余并排序的矩阵。根据以上可知，差分阵列的方向矩阵λ1为
[0080][0081]
其中：θk(k＝1，2，
…
，k)为k个窄带远场不相干信源的一维波达方向。
[0082]
步骤3：利用空间平滑算法，得到空间平滑协方差矩阵。
[0083]
空间平滑算法的基本思想是将等距线阵划分成若干相互重叠的子阵列，若各子阵列结构相同，则它们的协方差矩阵可以相加后平均取代原来意义上的接收信号协方差矩阵。
[0084]
式(8)已经表明了z的方向矩阵a(θ)对应的是完全连续的二级虚拟阵，且为范围为[-mbd，mbd]经阵元间距为d的长均匀线阵，阵元数为t＝2mb 1＝4(m(m 1)-1) 1。显然，对z按照二级虚拟阵元位置排序并去除冗余阵元获得z1并未改变虚拟阵元的范围和连续性。利用空间平滑算法，如图6所示，将二级虚拟阵列划分为mb 1个相互重叠的子阵列，每一个子阵列均包含mb 1个阵元。其中，第i个子阵的阵元位置为：
[0085]
{(-i 1 n)d，n＝0，1，
…
，mb}
ꢀꢀꢀꢀ
(11)
[0086]
其接收信号矩阵为z1的第mb 2-i行至第2mb 2-i行，记为z
1i
，并构造协方差矩阵：
[0087][0088]
将全部mb 1个子阵的协方差矩阵求和并计算均值，得到空间平滑协方差矩阵：
[0089][0090]
由此我们得到了空间平滑协方差矩阵r，该矩阵与基于经典子空间类算法的信号协方差矩阵具有相同的形式，接下来的步骤可直接对该矩阵进行操作。
[0091]
步骤4：根据空间平滑协方差矩阵分块估计传播算子并构造q矩阵。
[0092]
将空间平滑协方差矩阵进行分块如下：
[0093]
r＝[g，h]
ꢀꢀꢀꢀ
(14)
[0094]
其中，g和h分别为(mb 1)
×
k和(mb 1)
×
(mb 1-k)矩阵。
[0095]
在无噪声情况下，
[0096]
h＝gp
ꢀꢀꢀꢀꢀ
(15)当存在噪声时，虽然式(14)仍然可以分块，但是式(15)已经不再满足条件了，p的估计值可以通过对代价函数最小化获得，即
[0097]
j(p)＝||h-gp||2ꢀꢀꢀꢀ
(16)
[0098]
其中，||.||表示frobenius范数。代价函数j为关于p的二次凸函数，其最优解表示为
[0099]
p＝(ghg)-1ghhꢀꢀꢀꢀ
(17)
[0100]
得到了传播算子p后，接下来构造q矩阵：
[0101][0102]
步骤5、构造多项式并进行求根处理得到波达方向估计角度；
[0103]
pm算法中构造出的q矩阵的列span{q}组成的子空间是包含在span{en}(en为噪声子空间)中的。因为q包含它的mb 1-k列是线性无关的。因此有：
[0104]
span{q}＝span{en}
ꢀꢀꢀꢀ
(19)
[0105]
传播算子定义了与r的最小特征值对应的特征向量矩阵en构成噪声子空间。pm算法中构造的估计函数为
[0106]fpm
(θ)＝ah(θ)qqha(θ)
ꢀꢀꢀ
(20)
[0107]
求根pm算法是pm算法的一种多项式求根形式，其基本思想是pisarenko分解。定义多项式
[0108][0109]
式中，u
l
是矩阵r的第l个特征向量，且z＝exp(jw)(当z＝exp(jw)时，即多项式的根正好在单位圆上，w表示空间频率)，p(exp(jw))是一个以空间频率为w的导向矢量，且与噪声子空间是正交的。
[0110]
从噪声特征向量矩阵en中提取信息可转化为从q矩阵中提取信息，我们希望求pm
[0111]
ph(z)qqhp(z)
ꢀꢀꢀꢀꢀꢀ
(22)
[0112]
的零点。由于只对单位圆上的z值感兴趣，所以可以用p
t
(z-1
)代替ph(z)，这就给出了求根pm多项式，即
[0113]
[0114]
注意，p(z)是2mb次多项式，它的根相对于单位圆为镜像对。其中，取单位圆内具有最大幅值的k个根z1，z2，...，zk的相位给出波达方向估计，即
[0115][0116]
其中，λ为波长，d为阵元间距。
[0117]
本发明方法计算复杂度分析如下：
[0118]
对本发明算法的运算复杂度进行分析，具体如下：其中，嵌套阵子阵1和子阵2的大小为均为m。信源数为k，快拍数为l。以复乘次数作为依据，本节算法的主要复杂度包括：计算四阶累积量矩阵需要o{(2m)4l}，t表示二级虚拟阵元总数，且t＝4(m(m 1)-1) 1，二级虚拟阵列划分为个相互重叠的子阵列，每个子阵列均包含个w阵元，且w＝(t 1)/2。空间平滑后计算每个虚拟子阵的协方差矩阵需要o{w3}，对空间平滑协方差矩阵r进行特征值分解需要o{w3}，估计传播算子需要o{k3 kw(k w)}，多项式进行求根需要o{(2w-1)3}，谱峰搜索的复杂度为o{nw(2(w-k) 1)}其中搜索次数n＝μ/0.01(μ是角度搜索范围且0.01为搜索精度)。
[0119]
因此本发明提出算法总复杂度为o{(2m)4l w3 k3 kw(k w) (2w-1)3}，ss-rootmusic算法的总复杂度为o{(2m)4l 2w3 (2w-1))3}，ss-pm算法总复杂度为o{(2m)4l w3 k3 kw(k w) nw(2(w-k) 1)}，其中角度搜索范围μ取90。
[0120]
各算法及其相应的算法复杂度如表1所示。
[0121]
表1不同算法的复杂度
[0122]
算法计算复杂度ss-rootmusico{(2m)4l 2w3 (2w-1))3}ss-pmo{(2m)4l w3 k3 kw(k w) nw(2(w-k) ss-rootpmo{(2m)4l w3 k3 kw(k w) (2w-1)3}
[0123]
图7是本发明算法角度估计点阵图。考虑有k＝13个不相关窄带信号入射至m1＝m2＝m＝6的二级嵌套线阵中，信号的方位角在0
°
到65
°
间均匀分布，l＝1000，snr＝0db。其中横坐标trial number(测试次数)，纵坐标doa estimates(到达角估计)，表示通过改变测试次数观测角度估计的结果。图8和9展示了所提算法的估计结果。横坐标snr(信噪比)表示在以信噪比为变量的条件下观测纵坐标rmse(均方根误差)的变化，从而展示算法的估计性能。从图中可以看出，该算法能有效估计出信源角度，且可估计的信源数大于实际阵元数。
[0124]
图8是本发明算法角度估计性能在不同快拍下的曲线图。快拍数增加，即采样数据增多。由图可以得出，算法的角度估计性能随着快拍数增加变得更好。其中，入射信号的角度参数(θ1，θ2)＝(10
°
，40
°
)，嵌套阵的大小为m1＝m2＝m＝6。
[0125]
图9是本发明算法角度估计性能在不同阵元数下的曲线图。阵元数增加，即接收天线获得的分集增益增加。由图可以得出，算法的角度估计性能随着阵元数增加变得更好。其中，入射信号的角度参数(θ1，θ2)＝(10
°
，40
°
)，快拍数l＝1000。
[0126]
图10是所提ss-rootpm算法和ss-rootmusic算法、ss-pm算法在不同信噪比的仿真对比结果。由图9表明，所提算法的角度估计性能高于ss-pm算法，低于ss-rootmusic算法。其中，入射信号的角度参数(θ1，θ2)＝(10
°
，40
°
)，嵌套阵的大小为m1＝m2＝m＝7。
[0127]
图11是本发明算法和ss-rootmusic算法、ss-pm算法在不同快拍数下的计算复杂
度。其中，横坐标number of snapshots(快拍数)，纵坐标complex multiplication times(复乘次数)，表示复乘次数随着快拍数的变化而变化。由图可以看出在计算复杂度上，本发明算法远远低于ss-pm算法，略低于ss-rootmusic算法；且随着快拍数的增加，三者的计算复杂度也变得越大。
[0128]
图12是本发明算法和ss-rootmusic算法、ss-pm算法在不同阵元数下的计算复杂度。其中，横坐标number of elements(阵元数)，纵坐标complex multiplication times(复乘次数)，表示复乘次数随着阵列阵元元素个数的变化而变化。由图可以看出在计算复杂度上，本发明算法远远低于ss-pm算法，略低于ss-rootmusic算法；且随着阵元数目的增加，三者的计算复杂度也变得越大。
[0129]
以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。