基于四阶累量的嵌套阵列波达方向角估计方法
【专利摘要】本发明公开了一种基于四阶累量的嵌套阵列波达方向角估计方法，主要解决现有技术中阵元利用率低，信号识别数量少的问题，其实现步骤是：1)分别构建一层嵌套阵列和二层嵌套阵列；2)获取一层嵌套阵列的输出信号和二层嵌套阵列的输出信号；3)根据一层嵌套阵列输出信号与二层嵌套阵列输出信号计算所有四阶累量；4)根据所有四阶累量构造四阶累量矩阵；5)计算四阶累量矩阵的噪声子空间；6)根据四阶累量的噪声子空间和阵列流型矩阵计算空间谱；7)根据空间谱绘制幅度谱图，得到波达方向角。本发明在阵元数量有限的情况下大大提高了阵列可识别的信源数目，适用于目标侦察和无源定位。
【专利说明】
基于四阶累量的嵌套阵列波达方向角估计方法
技术领域
[0001] 本发明属于信号处理技术领域，特别设及一种电磁信号的阵列信号波达方向角估 计方法，可用于对飞机、舰船运动目标的侦察与无源定位。
【背景技术】
[0002] 信号的波达方向角DOA估计是阵列信号处理领域的一个重要分支，它是指利用天 线阵列对空间声学信号、电磁信号进行感应接收，再运用现代信号处理方法快速准确的估 计出信号源的方向，在雷达、声纳、无线通信等领域具有重要应用价值。随着科技的不断进 步，对阵列在进行信号波达方向估计时达到的自由度也有越来越高的要求。
[0003] 针对该问题的研究中，出现较早、应用较为广泛的是多重信号分类MUSIC子空间的 模型，对于一个L阵元的典型线性均匀阵列，传统的MUSIC类计算方法可检测的信源数目是 L-I个。之后的大部分算法都是利用该模型生成的，例如信号参数估计旋转不变技术 ESPRIT。运些算法采用典型的线性均匀阵列，造成估计的信号数目低于阵元数目，目标个数 很多时甚至无法识别，导致目标捕获失败。
[0004] 为了在少的阵元条件下得到尽量大的角度自由度，检测更多的信源，一些新的阵 列结构被提出，比较有代表性的是嵌套阵列W及互质阵列。Piya Pal等人在其发表的论文 "Nested Arrays=A Novel Approach to Array Processing With Enhanced Degrees of Freedom"(《IEEE transactions on signal processing》，VOL 58,NO.8,August 2010)中 公开了 一种基于嵌套阵列的DOA估计方法，该方法能够使用M+N个阵元，生成2MN+2N-1个虚 拟阵元，可检测MN+N-1个信号。该方法具有估计多于阵元数目的信号数的能力，但是，该方 法仍然存在的不足之处是，在阵元数量一定的情况下，能估计的信号数量有时仍不能满足 实际应用。
[0005] 为解决上述问题，近年来，出现了使用四阶累量增加阵列虚拟阵元的方法，可使信 号估计数量得到较大提升。但已有的基于四阶累量的波达角度估计方法多数是针对均匀阵 列，而对于嵌套阵列的算法，由于没有合理设置阵元位置而不能使阵元达到最大利用率。在 实际应用中，给定一定数量的阵元，如果不能合理利用运些阵元获得足够多的虚拟阵元，就 不能估计足够多的信号，造成侦察和定位资源的浪费。
【发明内容】
[0006] 本发明的目的在于针对上述现有技术存在的不足，提出一种基于四阶累量的嵌套 阵列波达方向角估计方法，W在阵元数量一定的情况下，提高能够进行估计的信号数量，避 免因不能合理利用阵元造成的资源浪费。
[0007] 为实现上述目的，本发明技术方案包括如下：
[000引（1)用M+N个天线接收机形成一层嵌套阵列：
[0009] (Ia)将每个天线接收机称为一个阵元，用M个天线接收机形成第一均匀线性阵列 曰，其阵元间距为d;用N个天线接收机形成第二均匀线性阵列b，其阵元间距为(M+l)d;定义 第一均匀线性阵列a的第一个阵元为起始阵元，其中，1>1，1^>1，0<(1《^/2，^为入射到阵列 的窄带信号波长；
[0010] (化)将第一均匀线性阵列a与第二均匀线性阵列b组合为一层嵌套阵列:将第二均 匀线性阵列b的第一个阵元放置于与起始阵元相距为Md的位置;将第二均匀线性阵列b的所 有阵元依次插于第一均匀线性阵列a中，形成一层嵌套阵列；
[0011] (Ic)计算一层嵌套阵列的虚拟阵元数A = 2(M+1)N-1;
[0012] (2)用P+Q个天线接收机形成二层嵌套阵列：
[0013] (2a)用P个天线接收机形成第S均匀线性阵列C，其阵元间距为AcU用Q个天线接收 机形成第四均匀线性阵列f，其阵元间距为(P+l)Ad，其中，Q>1，P>1;
[0014] (2b)将第S均匀线性阵列C与第四均匀线性阵列f组合为二层嵌套阵列，即将第S 均匀线性阵列C的第一个阵元放置于与起始阵元相距为(A-l)d的位置，将第四均匀线性阵 列f的第一个阵元放置于与起始阵元相距为[(P+l)A-l]d的位置，并将第四均匀线性阵列f 的所有阵元依次插于第=均匀线性阵列C中，组合为二层嵌套阵列；
[0015] (2c)计算二层嵌套阵列的虚拟阵元数B = 2(P+1)Q-1;
[0016] (3)由一层嵌套阵列和二层嵌套阵列的天线接收机对空间目标信号进行采样，分 别得到一层嵌套阵列输出信号Y(t)和二层嵌套阵列输出信号Z(t)，其中，Y(t) = [yo (t)，…，yi(t)，???，yM+^t)]，yi(t)表示一层嵌套阵列第i个阵元的输出信号，i的取值范围是 0《i《M+N，Z(t) = [zo(t)，???，zj(t)，…，zp+Q(t)]，zj(t)表示二层嵌套阵列的第j个阵元的 输出信号J的取值范围是0刮《P+Q;
[0017] (4)计算一层嵌套阵列输出信号Y(t)与二层嵌套阵列输出信号Z(t)的所有四阶累 量c(kl ,k2,k3,k4) = cum(yki(t) ,yk2(t)*,zk3(t) ,zk4(t)*)，
[001引其中，0《kl ,k2《M+N, l《k3 ,k4《P+Q，cum表示求解四阶累量操作，（?）*表示向 量的共辆运算；
[0019] (5)定义中间变量D=(AB+2A-l)/2,从所有四阶累量C化1，42,43,44)中依次找出 kl，k2，k3，k4满足kl+k3A-k2-k4A = -D，…，-1，0，1，…，D条件的四阶累量，并依次定义为四 阶累量元素g(-D)，…，旨(-1)，旨(0)，旨(1)，-'，旨化），将运些四阶累量元素重新排列形成四阶 累量矩阵G:
[0020]
;
[0021] (6)计算四阶累量矩阵G的噪声子空间Un;
[0022] (7)根据阵列流型矩阵a(0)和四阶累量矩阵G的噪声子空间Un,计算空间谱P(0);
[0023] (8) W波达方向角范围0的值为X轴坐标，W空间谱P(0)的幅度值为y轴坐标，绘制 幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，运些谱峰的峰 值点所对应的X轴坐标即为目标的波达方向角度值，其中，K表示入射到一层嵌套阵列与二 层嵌套阵列的空间目标信号个数，K>1。
[0024] 本发明与现有技术相比具有W下优点：
[0025] 1)本发明采用了双层嵌套阵列模型进行波达方向角度估计，克服了现有技术中采 用典型的线性均匀阵列，造成估计的信号数目低于阵元数目的缺点，提高了在阵元数目相 同的条件下的阵列可识别信源数目。
[0026] 2)本发明将四阶累量应用到嵌套阵列的DOA估计中，通过使用四阶累量，设置阵元 位置，可使用M+N+P+Q+1个阵元获得(2M+2N+1) (2PQ+2Q+1)个连续虚拟阵元，大大提高了阵 列利用率，同时通过获得更多的虚拟阵元，进一步增加了阵列可识别的信源数目。
[0027] 3)本发明采用了双层嵌套阵列模型进行波达方向角度估计，相比于使用其他新型 阵列模型对阵元的数目要求更低，提高了阵元数目使用的灵活性。
【附图说明】
[0028] 图1是本发明的实现流程图；
[0029] 图2是本发明中一层嵌套阵列与二层嵌套阵列的结构示意图。
【具体实施方式】
[0030] W下参照附图，对本发明的技术方案和效果作进一步的详细说明。
[0031 ]参附图1，本发明的具体步骤如下：
[0032] 步骤1:用M+N个天线接收机形成一层嵌套阵列，得到其一层嵌套阵列虚拟阵元数 A。
[0033] (Ia)将每个天线接收机称为一个阵元，用M个天线接收机形成第一均匀线性阵列 曰，其阵元间距为d;用N个天线接收机形成第二均匀线性阵列b，其阵元间距为(M+l)d;定义 第一均匀线性阵列a的第一个阵元为起始阵元，其中，1>1，1^>1，0<(1《^/2，^为入射到阵列 的窄带信号波长；
[0034] (化)将第一均匀线性阵列a与第二均匀线性阵列b组合为一层嵌套阵列:将第二均 匀线性阵列b的第一个阵元放置在与起始阵元相距为Md的位置;将第二均匀线性阵列b的所 有阵元依次插于第一均匀线性阵列a中，形成一层嵌套阵列；
[0035] (Ic)计算一层嵌套阵列的虚拟阵元数A = 2(M+1)N-1。
[0036] 步骤2:用P+Q个天线接收机形成二层嵌套阵列，得到二层嵌套阵列的虚拟阵元数 B。
[0037] (2a)用P个天线接收机形成第S均匀线性阵列C，其阵元间距为AcU用Q个天线接收 机形成第四均匀线性阵列f，其阵元间距为(P+l)Ad，其中，Q>1，P>1;
[0038] (2b)将第S均匀线性阵列C与第四均匀线性阵列f组合为二层嵌套阵列，即将第S 均匀线性阵列C的第一个阵元放置于与起始阵元相距为(A-l)d的位置，将第四均匀线性阵 列f的第一个阵元放置于与起始阵元相距为[(P+l)A-l]d的位置，并将第四均匀线性阵列f 的所有阵元依次插于第=均匀线性阵列C中，组合为二层嵌套阵列；
[0039] (2c)计算二层嵌套阵列的虚拟阵元数B = 2(P+1)Q-1。
[0040] 步骤3:获得一层嵌套阵列输出信号和二层嵌套阵列输出信号。
[0041] 由一层嵌套阵列和二层嵌套阵列的天线接收机对空间目标信号进行采样，分别得 到一层嵌套阵列输出信号Y(t)和二层嵌套阵列输出信号Z(t)，其中，Y(t) = [yo(t)，…，yi (t)，???，yM+N(t)]，yi(t)表示一层嵌套阵列第i个阵元的输出信号，i的取值范围是0《i《M+ N，Z(t) = [Z0(t)，…，Zj(t)，…，ZP+Q(t)]，Zj(t)表示二层嵌套阵列的第j个阵元的输出信号， j的取值范围是〇《j《P+Q。
[0042] 步骤4:计算输出信号的所有四阶累量。
[0043] 计算一层嵌套阵列输出信号Y(t)与二层嵌套阵列输出信号Z(t)的所有四阶累量
[0044] c(kl ,k2,k3,k4) = cum(yki(t) ,yk2(t)*,zk3(t) ,zk4(t)*)，
[0045] 其中，0《kl ,k2《M+N, l《k3 ,k4《P+Q,cum表示求解四阶累量操作，（?）*表示向 量的共辆运算。
[0046] 步骤5:构造四阶累量矩阵G。
[0047] 定义中间变量D= (AB+2A-1 )/2，从所有四阶累量C化1，k2，k3，k4)中依次找出kl， k2，k3，k4满足kl+k3A-k2-k4A =-D，…，-I，0，I，…，D条件的四阶累量，并依次定义为四阶累 量元素g(-D)，…，g(-l)，g(0)，g(l)，…，g(D)，将运些四阶累量元素重新排列，形成四阶累 量矩阵G:
[004引
[0049] 步骤6:计算四阶累量矩阵的噪声子空间Un。
[0050] (6a)对四阶累量矩阵G进行如下特征分解：
[0051] G = U-A-IjH,
[0052] 其中，A为四阶累量矩阵（.）H的特征值矩阵，U为四阶累量矩阵G的特征值所对应 的特征向量矩阵，（?）H表示矩阵的共辆转置远算；
[0053] (6b)将特征值矩阵A中的特征值按从大到小排序，取后D-K个较小特征值对应的 特征向量矩阵作为噪声子空间Un。
[0054] 步骤7:计算空间谱P(目）。
[0055] (7a)将观测空域[-90°，90° ]按等间隔L划分成E个角度，定义波达方向角范围为0 =[目1，目2，'。，目6,…目E]，其中，目e表示第e个目标波达方向角，e = l，2,…，E，E〉〉max{M，N，P, QhL的取值根据期望达到的角度估计精度进行设定，网格划分间隔越小则最终得到的角度 估计值精度越高；
[0056] (7b)计算第e个目标波达方向角06对应的导向矢量a(0e):
[0化7]
[005引其中，（?）嗦示矩阵转置运算，j为虚数单位；
[0059] (7c)计算所有的目标波达方向角对应的导向矢量，得到阵列流型矩阵a(0):
[0060] a(目）= [a(目 1),a(目e),a(目E)];
[0061] (7d)根据阵列流型矩阵a(0)和四阶累量矩阵G的噪声子空间Un，计算空间谱P(0):
[0062]
[0063] 其中，（?）H表示矩阵的共辆转置远算。
[0064] 步骤8:绘制幅度谱图。
[0065] (8a似波达方向角范围0的值为X轴坐标，W空间谱P(0)的幅度值为y轴坐标，绘制 幅度谱图；
[0066] (8b)从绘制的幅度谱图中，按照从高到低的顺序将所有谱峰幅度值进行排序，将 前K个谱峰的峰值点所对应的X轴坐标做为目标的波达方向角度值，其中，K表示入射到一层 嵌套阵列与二层嵌套阵列的空间目标信号个数，且假设空间目标信号在传播过程中加入了 均值为零的复高斯白噪声，K^l。
[0067] 下面结合仿真实例对本发明的效果做进一步的描述。
[0068] 仿真实例1.构建双层嵌套阵列结构。
[0069] 1.1)利用第一均匀线性阵列a与第二均匀线性阵列b形成一层嵌套阵列，设N = 2，M =2,第一均匀线性阵列a含有3个阵元，阵元间距为d，第二均匀线性阵列b含有3个阵元，阵 元间距为3d，第一均匀线性阵列a的第一个阵元放置在Id位置；
[0070] 1.2)计算第一均匀线性阵列a与第二均匀线性阵列b形成的一层嵌套阵列虚拟阵 元数 A = 2MN+2N-1 = 11;
[0071] 1.3)利用第S均匀线性阵列C与第四均匀线性阵列f形成二层嵌套阵列，设P = 2，Q =2,第S均匀线性阵列C含有2个阵元，阵元间距为lld，第四均匀线性阵列f含有2个阵元， 阵元间距为33d，第S均匀线性阵列C的第一个阵元在Ild的位置，计算虚拟阵元数B = 2PQ+ 2Q-1 = 11。
[0072] 其中第一均匀线性阵列a、第二均匀线性阵列b和第S均匀线性阵列c、第四均匀线 性阵列f形成了双层嵌套阵列，其阵列结构图如图2所示。
[0073] 由图2可见，嵌套阵列结构在阵元数量一定的情况下，可W得到更多的阵元位置信 息，从而增加阵列可识别信源数目，同时，嵌套阵列结构相比于使用其他新型阵列模型对阵 元的数目要求更低，提高了阵元数目使用的灵活性。
[0074] 仿真实例2,根据四阶累量计算双层嵌套阵列所有虚拟阵元。
[0075] 2.1)设第一均匀线性阵列a与第二均匀线性阵列b形成的一层嵌套阵列阵元位置 为[l，2,3,6]d，；
[0076] 2.2)设第S均匀线性阵列C与第四均匀线性阵列f形成的二层嵌套阵列阵元位置 为[ll，22,33,66]d;
[0077] 2.3)计算得到一层嵌套阵列每个阵元分别与二层嵌套阵列中每个阵元形成的虚 拟阵列为[12,13,15，17]d、[23,24,26，28]d、[34,35,36，39]d、[67,68,70，72]d;
[0078] 2.4)根据四阶累量计算的方法和原则，用上述运些虚拟阵列生成最终的所有虚拟 阵元，结果如表1。
[0079] 表1双层嵌套阵列所有虚拟阵元
[00801
[0081 ] 从表1可见，该双层嵌套阵列用M+N+P+Q = 8个原始阵元生成了一个从-60d到60d共 含有AB = 121个虚拟连续阵元的线性阵列。
[0082]综上，本发明解决了现有技术阵元利用率低，识别信源数目少，无源定位估计误差 大的问题，降低了对阵元数目的要求，保证了阵元数目使用的灵活性，提高了一定阵元数情 况下阵列可识别的信源数目W及低信噪比下对信号方向角的估计性能。
【主权项】
1. 一种基于四阶累量的嵌套阵列的波达方向角估计方法，其特征在于，包括以下步骤： (1) 用M+N个天线接收机形成一层嵌套阵列： (la) 将每个天线接收机称为一个阵元，用M个天线接收机形成第一均匀线性阵列a，其 阵元间距为d;用N个天线接收机形成第二均匀线性阵列b，其阵元间距为(M+l)d;定义第一 均匀线性阵列a的第一个阵元为起始阵元，其中，M多1，N多1，0〈(1<λ/2，λ为入射到阵列的窄 带信号波长； (lb) 将第一均匀线性阵列a与第二均匀线性阵列b组合为一层嵌套阵列:将第二均匀线 性阵列b的第一个阵元放置于与起始阵元相距为Md的位置;将第二均匀线性阵列b的所有阵 元依次插于第一均匀线性阵列a中，形成一层嵌套阵列； (lc) 计算一层嵌套阵列的虚拟阵元数A = 2 (M+1)N-1; (2) 用P+Q个天线接收机形成二层嵌套阵列： (2a)用P个天线接收机形成第三均匀线性阵列c，其阵元间距为AcU用Q个天线接收机形 成第四均匀线性阵列f，其阵元间距为(P+l)Ad，其中，Q多1，P彡1; (2b)将第三均匀线性阵列c与第四均匀线性阵列f组合为二层嵌套阵列，即将第三均匀 线性阵列c的第一个阵元放置于与起始阵元相距为(A-l)d的位置，将第四均匀线性阵列f的 第一个阵元放置于与起始阵元相距为[(P+l)A-l]d的位置，并将第四均匀线性阵列f的所有 阵元依次插于第三均匀线性阵列c中，组合为二层嵌套阵列； (2c)计算二层嵌套阵列的虚拟阵元数B = 2(P+1)Q-1; (3) 由一层嵌套阵列和二层嵌套阵列的天线接收机对空间目标信号进行采样，分别得 到一层嵌套阵列输出信号Y(t)和二层嵌套阵列输出信号Z(t)，其中，Y(t) = [y〇(t)，H_，yi (t)，···，yM+N(t)]，yi(t)表示一层嵌套阵列第i个阵元的输出信号，i的取值范围是(Xi彡M+ N，Z(t) = [ZO(t)，···，Zj(t)，…，ZP+Q(t)]，Zj(t)表示二层嵌套阵列的第j个阵元的输出信号， 的取值范围是〇彡j彡P+Q; (4) 计算一层嵌套阵列输出信号Y(t)与二层嵌套阵列输出信号Z(t)的所有四阶累量c (kl，k2，k3，k4) = cum(yki(t)，yk2(t)*，zk3(t)，zk4(t)*)，其中，0<kl，k2<M+N，Kk3，k4<P +Q，cum表示求解四阶累量操作，（· Γ表示向量的共辄运算； (5) 定义中间变量D = (AB+2A-1)/2，从所有四阶累量c (kl，k2，k3，k4)中依次找出kl， k2，k3，k4满足kl+k3A-k2-k4A =-D，…，-I，O，I，…，D条件的四阶累量，并依次定义为四阶累 量元素 g(_D)，…，g(-l)，g(0)，g(l)，~，g(D)，将这些四阶累量元素重新排列形成四阶累量 矩(6) 计算四阶累量矩阵G的噪声子空间Un; (7) 根据阵列流型矩阵α(θ)和四阶累量矩阵G的噪声子空间Un，计算空间谱Ρ(θ); (8) 以波达方向角范围Θ的值为X轴坐标，以空间谱Ρ(θ)的幅度值为y轴坐标，绘制幅度 谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点 所对应的X轴坐标即为目标的波达方向角度值，其中，K表示入射到一层嵌套阵列与二层嵌 套阵列的空间目标信号个数，K多1。2. 根据权利要求1所述的基于四阶累量的嵌套阵列波达方向角估计方法，其中步骤(6) 中计算四阶累量矩阵G的噪声子空间Un，按如下步骤进行： (6a)对四阶累量矩阵G进行如下特征分解： 其中，Λ为四阶累量矩阵（·）H的特征值矩阵，U为四阶累量矩阵G的特征值所对应的特 征向量矩阵，（·）H表示矩阵的共辄转置远算； (6b)将特征值矩阵Λ中的特征值按从大到小排序，取后D-K个较小特征值对应的特征 向量矩阵作为噪声子空间Un。3. 根据权利要求1所述的基于四阶累量的嵌套阵列波达方向角估计方法，其中步骤(7) 中的阵列流型矩阵α (Θ )，按如下步骤构造： (7a)将观测空域[-90°，90°]等间隔划分成E个角度，定义为波达方向角范围Q = M1, Θ2,…，9e，…ΘΕ]，其中，0e表示第e个目标波达方向角， e = l，2,…，E，E?max{M，N，P，Q}; (7b)计算第e个目标波达方向角对应的导向矢量Ci(Qe3):其中，（·）T表示矩阵转置运算，j为虚数单位； (7c)计算所有的目标波达方向角对应的导向矢量，得到阵列流型矩阵α(θ): α(θ) = [α(θι) ,---,Q(Oe)4. 根据权利要求1所述的基于四阶累量的嵌套阵列波达方向角估计方法，其中步骤(7) 中计算空间谱Ρ(θ)，按如下形式进行：其中，α(θ)表示阵列流型矩阵，（·）Η表示矩阵的共辄转置运算。
【文档编号】G01S3/14GK106019215SQ201610560482
【公开日】2016年10月12日
【申请日】2016年7月15日
【发明人】蔡晶晶, 宗汝, 赵晗希, 武斌, 李鹏, 苏瑶
【申请人】西安电子科技大学