基于嵌套阵列的非圆信号波达方向角估计方法与流程

本发明属于信号处理技术领域，特别涉及一种电磁信号的阵列信号波达方向角估计方法，可用于对飞机、舰船运动目标的侦察与无源定位。

背景技术：

信号的波达方向角doa估计是阵列信号处理领域的一个重要分支，它是指利用天线阵列对空间声学信号、电磁信号进行感应接收，再运用现代信号处理方法快速准确的估计出信号源的方向，在雷达、声纳、无线通信等领域具有重要应用价值。随着科技的不断进步，对阵列在进行信号波达方向估计时达到的自由度也有越来越高的要求。

在现代信号处理中，由于二相相移键控以及m进制幅移键控等非圆信号的应用越来越多，因此有关非圆信号的doa估计受到了越来越多的关注。pcharge等人在其发表的论文“anon-circularsourcesdirectionfindingmethodusingpolynomialrooting”(《signalprocessing》,vol81,pp.1765-17702001)中公开了一种利用多项式求解进行非圆信号doa估计的方法，但是，该方法仍然存在的不足是，该方法采用典型的线性均匀阵列，造成估计的信号数目低于阵元数目，目标个数很多时甚至无法识别，导致目标捕获失败。

为了解决上述问题，在较少的阵元条件下得到尽量大的角度自由度，检测更多的信源，一些新的阵列结构被提出，比较有代表性的是嵌套阵列以及互质阵列。ppiya等人在其发表的论文“nestedarrays:anovelapproachtoarrayprocessingwithenhanceddegreesoffreedom”(《ieeetransactionsonsignalprocessing》,vol58,no.8,august2010)中公开了一种基于嵌套阵列的doa估计方法，该方法能够使用m+n个阵元，生成2mn+2n-1个虚拟阵元，可检测mn+n-1个信号。该方法具有估计多于阵元数目的信号数的能力，但是，该方法中仍然存在的不足是，对于嵌套阵列的讨论都集中在接收信号为圆信号的条件下，对于如何利用嵌套阵列进行非圆信号的处理目前还没有研究。

在实际应用中，对于非圆信号环境，给定一定数量的阵元，如果不能合理利用这些阵元以及信号的非圆特性，就不能估计足够多的信号，造成侦察和定位资源的浪费。

技术实现要素：

本发明的目的在于针对上述现有技术存在的不足，提出一种基于嵌套阵列的非圆信号波达方向角估计方法，以在非圆信号环境下，利用嵌套阵列进行信号处理，提高能够进行估计的信号数量，避免因不能合理利用阵元和信号特性造成的资源浪费。

为实现上述目的，本发明技术方案包括如下：

(1)用m+n个天线接收机形成嵌套阵列，其中m、n分别表示两个天线接收阵列的阵元数，其取值范围为m≥1，n≥1；

(2)假设空间中有k个非圆目标信号，使用嵌套阵列天线接收机，对空间目标信号进行快拍采样和匹配滤波操作，得到嵌套阵列输出信号：y(t)＝[y1(t),…,yi(t),…,ym+n(t)]^t，其中，k的取值范围是k<mn+m+n-1，yi(t)表示嵌套阵列第i个阵元的输出信号，t的取值范围是1≤t≤l，l表示快拍数，i的取值范围是1≤i≤m+n，(·)^t表示矩阵转置运算；

(3)利用嵌套阵列输出信号y(t)，计算协方差矩阵rd和椭圆协方差矩阵rs：

其中，(·)^h表示矩阵共轭转置运算；

(4)根据协方差矩阵rd和椭圆协方差矩阵rs中的元素，分别构造等效协方差向量rd和等效椭圆协方差向量rs：

rd＝[rd(1,1),rd(2,1),…,rd(m+n,1),rd(1,2),…,rd(m+n,2),

…,rd(i,j),…,rd(1,m+n),…,rd(m+n,m+n)]^t

rs＝[rs(1,1),rs(2,1),…,rs(m+n,1),rs(1,2),…,rs(m+n,2),

…,rs(i,j),…,rs(1,m+n),…,rs(m+n,m+n)]^t

其中，rd(i,j)表示协方差矩阵rd中位于第i行，第j列的元素，i的取值范围为1≤i≤m+n，j的取值范围为1≤j≤m+n，rs(i,j)表示椭圆协方差矩阵rs中位于第i行，第j列的元素；

(5)计算等效协方差向量rd中所有元素的维数ei,j和等效椭圆协方差向量rs中所有元素的维数fi,j：

ei,j＝d(j)-d(i)

fi,j＝d(j)+d(i)

其中，d(i)表示嵌套阵列中第i个阵元的位置，d(j)表示嵌套阵列中第j个阵元的位置；

(6)删除等效协方差向量rd中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列协方差向量删除等效椭圆协方差向量rs中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列椭圆协方差向量

(7)根据虚拟均匀阵列协方差向量和虚拟均匀阵列椭圆协方差向量得到波达角选择矩阵g；

(8)利用矩阵特征值分解的方法，计算波达角选择矩阵g的噪声子空间un；

(9)提取噪声子空间un的前l1行和前(l1+l2-k)列的所有元素构成第一子矩阵，将该第一子矩阵作为第一噪声矩阵un1；提取噪声子空间un的后l2行和后(l1+l2-k)列的所有元素构成第二子矩阵，将该第二子矩阵作为第二噪声矩阵un2；

(10)根据第一噪声矩阵un1和第二噪声矩阵un2，构造如下多项式方程：

其中，p14(x)表示根据第一噪声矩阵un1和第二噪声矩阵un2构造的第一复合向量p14中第x个元素，p23(x)表示根据第一噪声矩阵un1和第二噪声矩阵un2构造的第二复合向量p23中第x个元素，z表示多项式方程的根，x的取值范围是1≤x≤2(l1+l2)-3；

(11)计算多项式方程的所有根，由多项式方程的每一个根的辐角与目标波达方向角度值的关系，得到目标波达方向角度值θ。

本发明与现有技术相比具有以下优点：

1)本发明采用了嵌套阵列模型进行波达方向角度估计，克服了现有技术中采用典型的线性均匀阵列，造成估计的信号数目低于阵元数目的缺点，增加了在阵元数目相同的条件下的阵列可识别信源数目，大大提高了阵列利用率。

2)本发明利用了非圆信号的信号特性，该类信号不仅具有协方差矩阵的特性，还具有椭圆协方差矩阵的特性，同时采用这两个矩阵进行信号估计相比于单使用协方差矩阵进行信号估计，使得可估计的信源数更多。

附图说明

图1是本发明的实现流程图；

图2是本发明中嵌套阵列的结构示意图。

具体实施方式

以下参照附图，对本发明的技术方案和效果作进一步的详细说明。

参附图1，本发明的具体步骤如下：

步骤1：用m+n个天线接收机形成嵌套阵列。

(1a)将每个天线接收机称为一个阵元，用m个天线接收机形成第一均匀线性阵列a，其阵元间距为d，定义第一均匀线性阵列a的第一个阵元为起始阵元，定义起始阵元位置d(1)＝1，第一均匀线性阵列a的其它阵元位置依次为d(2)＝2，d(3)＝3，d(4)＝4，…，d(m)＝m；其中，m的取值范围为m≥1，d的取值范围为0<d≤λ/2，λ为入射到阵列的窄带信号波长；

(1b)用n个天线接收机形成第二均匀线性阵列b，其阵元间距为(m+1)d，第二均匀线性阵列b的阵元位置依次设置为d(m+1)＝m+1，d(m+2)＝2(m+1)，d(m+2)＝3(m+1)，…，d(m+n)＝n(m+1)，其中，n的取值范围为n≥1；

(1c)将第二均匀线性阵列b的第一个阵元放置于与起始阵元相距为md的位置，将第二均匀线性阵列b的所有阵元依次插于第一均匀线性阵列a中，形成嵌套阵列。

步骤2：获得嵌套阵列输出信号y(t)。

假设空间中有k个非圆目标信号，使用嵌套阵列天线接收机，对空间目标信号进行快拍采样和匹配滤波操作，得到嵌套阵列输出信号：y(t)＝[y1(t),…,yi(t),…,ym+n(t)]^t，其中，k的取值范围是k<mn+m+n-1，yi(t)表示嵌套阵列第i个阵元的输出信号，t的取值范围是1≤t≤l，l表示快拍数，i的取值范围是1≤i≤m+n，(·)^t表示矩阵转置运算。

步骤3：计算协方差矩阵rd和椭圆协方差矩阵rs。

利用嵌套阵列输出信号y(t)，计算协方差矩阵rd和椭圆协方差矩阵rs：

其中，(·)^h表示矩阵共轭转置运算。

步骤4：构造等效协方差向量rd和等效椭圆协方差向量rs。

根据协方差矩阵rd和椭圆协方差矩阵rs中的元素，分别构造等效协方差向量rd和等效椭圆协方差向量rs：

rd＝[rd(1,1),rd(2,1),…,rd(m+n,1),rd(1,2),…,rd(m+n,2),

…,rd(i,j),…,rd(1,m+n),…,rd(m+n,m+n)]^t

rs＝[rs(1,1),rs(2,1),…,rs(m+n,1),rs(1,2),…,rs(m+n,2),

…,rs(i,j),…,rs(1,m+n),…,rs(m+n,m+n)]^t

其中，rd(i,j)表示协方差矩阵rd中位于第i行，第j列的元素，i的取值范围为1≤i≤m+n，j的取值范围为1≤j≤m+n，rs(i,j)表示椭圆协方差矩阵rs中位于第i行，第j列的元素。

步骤5：计算等效协方差向量和等效椭圆协方差向量中所有元素的维数。

根据等效协方差向量rd和等效椭圆协方差向量rs中每一个元素所在的行和列在嵌套阵列中对应的阵元位置，计算等效协方差向量rd中所有元素的维数ei,j和等效椭圆协方差向量rs中所有元素的维数fi,j：

ei,j＝d(j)-d(i)

fi,j＝d(j)+d(i)

其中，d(i)表示嵌套阵列中第i个阵元的位置，d(j)表示嵌套阵列中第j个阵元的位置。

步骤6：获得虚拟均匀阵列协方差向量和虚拟均匀阵列椭圆协方差向量

根据等效协方差向量rd中所有元素的维数，删除等效协方差向量rd中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列协方差向量

根据等效椭圆协方差向量rs中所有元素的维数，删除等效椭圆协方差向量rs中维数相同的元素和维数不连续的元素，并将剩余元素按维数从小到大排列，得到虚拟均匀阵列椭圆协方差向量

步骤7：构造波达角选择矩阵g。

根据虚拟均匀阵列协方差向量和虚拟均匀阵列椭圆协方差向量中所有元素进行行列排列，得到波达角选择矩阵g：

其中，l1＝(cd+1)/2，l2＝cs+1-(cd+1)/2，cd表示虚拟均匀阵列协方差向量中元素的个数，cd的取值为2mn+2n-1，cs表示虚拟均匀阵列椭圆协方差向量中元素的个数，cs的取值为mn+m+n，(·)*表示向量的共轭运算。

步骤8：计算波达角选择矩阵的噪声子空间un。

(8a)将波达角选择矩阵g进行特征分解，得到特征值矩阵和特征向量矩阵：

g＝u·∧·u^h

其中，λ为波达角选择矩阵g的特征值矩阵，u为矩阵g的特征值所对应的特征向量矩阵，(·)^h表示矩阵的共轭转置运算；

(8b)将特征值矩阵λ中的特征值按从大到小排序，取其后(l1+l2-k)个较小特征值对应的特征向量矩阵作为噪声子空间un。

步骤9：根据噪声子空间un获得第一噪声矩阵un1和第二噪声矩阵un2。

提取噪声子空间un的前l1行和前(l1+l2-k)列的所有元素构成第一子矩阵，将该第一子矩阵作为第一噪声矩阵un1；

提取噪声子空间un的后l2行和后(l1+l2-k)列的所有元素构成第二子矩阵，将该第二子矩阵作为第二噪声矩阵un2。

步骤10：构造多项式方程。

(10a)根据第一噪声矩阵un1，计算第一噪声向量c1：

c1＝[c1(1),c1(2),…,c1(u),…,c1(2l1-1)]

其中，c1(u)表示第一噪声向量c1中的第u个元素，u和为两个不同的中间变量，1≤u≤2l1-1，

(10b)根据第一噪声矩阵un1和第二噪声矩阵un2，计算第二噪声向量c2：

c2＝[c2(1),c2(2),…,c2(v),…,c2(l1+l2-1)]

其中，c2(v)表示第二噪声向量c2中第v个元素，v和为两个不同的中间变量，1≤v≤l1+l2-1，

(10c)根据第一噪声矩阵un1和第二噪声矩阵un2，计算第三噪声向量c3：

c3＝[c3(1),c3(2),…,c3(w),…,c3(l1+l2-1)]

其中，c3(w)表示第三噪声向量c3中的第w个元素，w和为两个不同的中间变量，1≤w≤l1+l2-1，

(10d)根据第二噪声矩阵un2，计算第四噪声向量c4：

c4＝[c4(1),c4(2),…,c4(z),…,c4(2l2-1)]

其中，c4(z)表示第四噪声向量c4中的第z个元素，z和为两个不同的中间变量，1≤z≤2l2-1，

(10e)根据第一噪声向量c1和第四噪声向量c4，计算第一复合向量p14：

p14＝[p14(1),p14(2),…,p14(x),…,p14(2l1+2l2-3)]

其中，p14(x)表示第一复合向量p14中的第x个元素，x和为两个不同的中间变量，1≤x≤2l1+2l2-3，

(10f)根据第二噪声向量c2和第三噪声向量c3，计算第二复合向量p23：

p23＝[p23(1),p23(2),…,p23(x),…,p23(2l1+2l2-3)]

其中，p23(x)表示第二复合向量p23中第x个元素，g为中间变量，其取值范围是max[1,x-l1-l2+2]≤g≤min[l1+l2-1,x]；

(10g)根据第一复合向量p14和第二复合向量p23中的元素，得到构造多项式方程：

其中，z表示多项式方程的根，z＝[z1,…,zh,…,zk]，zh表示多项式方程的第h个根，h的取值范围是1≤h≤k。

步骤11：获得目标波达方向角度值θ。

(11a)根据多项式方程，计算多项式方程的所有根z：

该多项式该多有2q＝2l1+2l2-4个根，其中每个根都有一个与其近似的根，每对近似根中根保留其中一个，就得到了该多有q个根z1，…zn,….zq，如果信号数k<q，则此处得到根的数量应该是k个，分别为z1,…,zh,…,zk，将其表示为：

z＝[z1,…,zh,…,zk]，

其中，z表示多项式方程的根，zh表示多项式方程的第h个根，h的取值范围是1≤h≤k。

(11b)由多项式方程的每一个根的辐角与相应的目标波达方向角度值的关系，得到相应的目标波达方向角度值：

θh＝arcsin(λ/(2πd)arg(zh))，

其中，θh表示第h个目标信号波达方向角度值；

(11c)由每一个的目标波达方向角度值，得到目标波达方向角度值θ：

θ＝[θ1,θ2,…,θh,…,θk]。

实施例：假设空间中有4个bpsk入射信号，其波长为λ，获取目标波达方向角度值θ。

第一步，根据虚拟均匀阵列协方差向量和虚拟均匀阵列椭圆协方差向量构造的原则及维数计算的方法，可得以下两个结果：

一是虚拟均匀阵列协方差向量中元素的个数为cd＝2mn+2n-1＝7，元素的维数范围为-3～3；

二是虚拟均匀阵列椭圆协方差向量中元素的个数为cs＝mn+m+n＝5，元素的维数范围为2～6；

第二步，根据第一步的结果，计算两个中间变量l1和l2：

l1＝(cd+1)/2＝4，l2＝cs+1-(cd+1)/2＝2

第三步，根据第二步计算的中间变量l1＝4和l2＝2，可得2l1-1＝5，l1+l2-1＝5，l1-l2+1＝3，将虚拟均匀阵列协方差向量和虚拟均匀阵列椭圆协方差向量中的元素带入上式，得到波达角选择矩阵g的计算公式中，得到本实例的波达角选择矩阵g：

第四步，根据波达角选择矩阵g，计算波达角选择矩阵g的噪声子空间un，并提取噪声子空间un的前4行和前3列所有元素构成的子矩阵，生成4×3维的第一噪声矩阵un1；提取噪声子空间un的后2行和后3列所有元素构成的子矩阵，生成2×3维的第二噪声矩阵un2；

第五步，根据第一噪声矩阵un1和第二噪声矩阵un2，构造多项式方程：

5.1)根据第一噪声矩阵un1，计算第一噪声向量c1中各元素为：

5.2)根据第一噪声矩阵un1和第二噪声矩阵un2，计算第二噪声向量c2中各元素为：

5.3)根据第一噪声矩阵un1和第二噪声矩阵un2，计算第三噪声向量c3中各元素为：

5.4)根据第二噪声矩阵un2，计算第四噪声向量c4中各元素为：

5.5)根据第一噪声向量c1和第四噪声向量c4，计算第一复合向量p14中各元素为：

5.6)根据第二噪声向量c2和第三噪声向量c3，计算第二复合向量p23中各元素为：

5.7)利用第一复合向量p14和第二复合向量p23中元素为系数，生成多项式方程：

[p14(1)-p23(1)]z^-4+[p14(2)-p23(2)]z^-3+[p14(3)-p23(3)]z^-2+[p14(4)-p23(4)]z^-1+[p14(5)-p23(5)]

+[p14(6)-p23(6)]z+[p14(7)-p23(7)]z²+[p14(8)-p23(8)]z³+[p14(9)-p23(9)]z⁴＝0其中z为多项式方程的根；

第六步，根据第五步构造的多项式方程，求解该多项式方程的根，可得到4对近似根，将每对近似根中去掉1个根，该后得到4个根z1、z2、z3、z4；

第七步，根据多项式方程的根，计算目标信号波达方向角度值：θh＝arcsin(λ/(2πd)arg(zh))，

其中，zh表示多项式方程的第h个根，θh表示第h个目标信号波达方向角度值，h的取值范围是1≤h≤4；

第八步，将第六步得到的4个根z1、z2、z3、z4，带入第七步的关系式，计算出每一个目标波达方向角度值θ1、θ2、θ3、θ4，得到所有目标的波达方向角度值θ＝[θ1,θ2,θ3,θ4]。

本发明的效果通过以下仿真实验进一步描述。

利用第一均匀线性阵列a与第二均匀线性阵列b形成嵌套阵列，设m＝1，n＝2，第一均匀线性阵列a含有1个阵元，阵元间距为d，第二均匀线性阵列b含有2个阵元，阵元间距为2d，第一均匀线性阵列a的第一个阵元放置在1d位置，形成的嵌套阵列阵元位置为[1,2,4]d，其阵列结构图如2所示。

由图2可见，嵌套阵列结构在阵元数量一定的情况下，可以得到更多的阵元位置信息，从而增加阵列可识别的信源数目，同时，嵌套阵列结构相比于使用均匀线性阵列结构，对阵元的数目要求更低，提高了阵元数目使用的灵活性。

综上，本发明解决了现有技术中阵元利用率低，识别信源数目少，非圆信号特性未能充分利用的问题，降低了对阵元数目的要求，根证了阵元数目使用的高效性，提高了一定阵元数情况下阵列可识别的信源数目以及低信噪比下对非圆信号方向角的估计性能。