CN106443574B - 基于双层嵌套阵列的波达方向角估计方法 - Google Patents

Abstract

本发明公开了一种基于双层嵌套阵列的波达方向角估计方法，主要解决现有技术运算量大，识别信源数目少和估计误差大的问题，其实现步骤是：1)分别构建一层嵌套阵列和二层嵌套阵列；2)获取一层嵌套阵列的输出信号和二层嵌套阵列的输出信号；3)根据一层嵌套阵列输出信号与二层嵌套阵列输出信号计算所有四阶累量，形成四阶累量向量；4)根据四阶累量向量构造超完备基，并定义空域稀疏向量；5)将探测空域波达方向角估计转化为求解稀疏约束方程；6)利用凸优化方法，得到空域稀疏向量的最稀疏解；7)根据最稀疏解绘制幅度谱图，得到波达方向角。本发明在阵元数量有限的情况下大大提高了阵列可识别的信源数目，可用于目标侦察和无源定位。

Description

基于双层嵌套阵列的波达方向角估计方法

技术领域

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

背景技术

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

针对该问题的研究中，出现较早、应用较为广泛的是多重信号分类MUSIC子空间的模型，对于一个L阵元的典型线性均匀阵列，传统的MUSIC类计算方法可检测的信源数目是L-1个。这些算法采用典型的线性均匀阵列，造成估计的信号数目低于阵元数目，目标个数很多时甚至无法识别，导致目标捕获失败。为了能够在阵元数量较少情况下获得尽量大的角度自由度，检测更多的信源，出现了嵌套阵列以及互质阵列。嵌套阵列能够使用M+N个阵元，生成2MN+2N-1个虚拟阵元，可检测MN+N-1个信号。互质阵列能够使用2M+N-1个阵元，生成2MN+2M-1个虚拟阵元，可检测MN+M-1个信号。

基于这种信号检测效率较高的嵌套阵列，出现了进一步提高阵列虚拟阵元数的基于四阶累量方法，可使信号估计数量得到较大提升。Qing Shen等人在其发表的论文“Extension of nested arrays with the fourth-order difference co-arrayenhancement”(the 41st IEEE international conference on acoustics,speech andsignal processing,shanghai,China,March 2016,2991-2995)公开了一种基于四阶累量的嵌套阵列DOA估计方法，该方法使用M+N+P个阵元，生成2[(3P+2)N(M+1)-2P-1]-1个虚拟阵元，可检测(3P+2)N(M+1)-2P-2个信号。但该方法相当于用一个均匀阵列对嵌套阵列进行阵列扩展，没有合理设置阵元位置而不能使阵元达到最大利用率。在实际应用中，如果给定了一定数量的阵元，却不能合理利用这些阵元获得足够多的虚拟阵元，就不能估计足够多的信号，造成侦察和定位资源的浪费和因角度估计误差引起的目标侦察失误。

发明内容

本发明的目的在于针对上述现有技术存在的不足，提出一种基于双层嵌套阵列的波达方向角估计方法，用一个嵌套阵列对嵌套阵列进行阵列扩展，在阵元数量一定的情况下，提高了可估计信号的数量，避免了不合理利用阵元造成的资源浪费和因角度估计误差引起的目标侦察失误。

为实现上述目的，本发明的实现步骤包括如下：

(1)利用天线接收机形成第一层嵌套阵列，并获取虚拟阵元数：

(1a)定义每个天线接收机为一个阵元，用M个天线接收机形成第一均匀线性阵列a，其阵元间距为d；用N个天线接收机形成第二均匀线性阵列b，其阵元间距为(M+1)d；定义第一均匀线性阵列a的第一个阵元为起始阵元，其中，M≥1,N≥1，0<d≤λ/2，λ为入射到嵌套阵列的窄带信号波长；

(1b)将第一均匀线性阵列a与第二均匀线性阵列b组合为第一层嵌套阵列：

将第二均匀线性阵列b的第一个阵元放置于与起始阵元相距为(M+1)d的位置，再将第二均匀线性阵列b的所有阵元依次插于第一均匀线性阵列a中，得到由M+N个天线接收机形成的第一层嵌套阵列；

(1c)计算第一层嵌套阵列的虚拟阵元数：A＝2MN+2N-1；

(2)利用天线接收机形成第二层嵌套阵列：

(2a)用P个天线接收机形成第三均匀线性阵列f，其阵元间距为Ad，用Q个天线接收机形成第四均匀线性阵列h，其阵元间距为(P+1)Ad，其中，Q≥1，P≥1；

(2b)将第三均匀线性阵列f与第四均匀线性阵列h组合为第二层嵌套阵列，即将第三均匀线性阵列f的第一个阵元放置于与起始阵元相距为(A+D)d的位置，将第四均匀线性阵列h的第一个阵元放置于与起始阵元相距为(P+1)Ad的位置，并将第四均匀线性阵列h的所有阵元依次插于第三均匀线性阵列f中，得到由P+Q个天线接收机形成的第二层嵌套阵列，其中D＝(AB-1)/2+MN+N，B＝2PQ+2Q-1；

(3)由第一层嵌套阵列和第二层嵌套阵列的天线接收机对空间目标信号进行采样，分别得到第一层嵌套阵列输出信号Y(t)和第二层嵌套阵列输出信号Z(t)，其中，

Y(t)＝[y₁(t),…,y_i(t),…,y_M+N(t)]，y_i(t)表示第一层嵌套阵列第i个阵元的输出信号，i的取值范围是1≤i≤M+N；

Z(t)＝[z₁(t),…,z_j(t),…,z_P+Q(t)]，z_j(t)表示第二层嵌套阵列的第j个阵元的输出信号，j的取值范围是1≤j≤P+Q；

(4)计算第一层嵌套阵列输出信号Y(t)与第二层嵌套阵列输出信号Z(t)的所有四阶累量R，形成四阶累量向量G；

(5)对探测空域波达方向角观测空间进行网格划分，构造出超完备基Φ(θ)，并定义一个Q×1维波达角空域稀疏向量s；

(6)根据四阶累量向量G和超完备基Φ(θ)建立约束关系，将探测空域波达方向角估计转化为求解如下稀疏约束方程：

min||s||₁

subject to||G-Φ(θ)s||₂≤η

s≥0

其中，||·₁表示表示求矩阵一阶范数操作，||·||₂表示求矩阵二阶范数操作，subject to表示约束关系，η为现实环境中的噪声标准差；

(7)利用凸优化方法求解上述稀疏约束方程，得到空域稀疏向量s的最稀疏解

(8)以波达方向角范围θ＝[θ₁,θ₂,...,θ_e,...,θ_E]的值为x轴坐标，以空域稀疏向量s的最稀疏解的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅度值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值，其中，K表示入射到第一层互质阵列与第二层嵌套阵列的空间目标信号个数，且假设空间目标信号在传播过程中加入了均值为零的复高斯白噪声，K≥1。

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

1)本发明将四阶累量应用到第一层嵌套阵列和第二层嵌套阵列形成的双层嵌套阵列模型的DOA估计中，通过使用四阶累量，合理设置阵元位置，可使用M+N+P+Q个阵元获得(2MN+2N-1)(2PQ+2Q-1)+2(P+1)(2MN+2N-1)+2(MN+N-1)个连续虚拟阵元，大大提高了阵列利用率，同时获得更多的虚拟阵元，增加了阵列可识别的信源数目。

2)本发明采用稀疏表示技术将波达方向角估计转化为稀疏信号的重构，是新理论技术与传统问题的结合；本发明利用信号源的空域稀疏性，得到的稀疏解中大系数对应的角度，突破了阵列分辨率的瑞利限，提高了目标侦察和无源定位在低信噪比、低快拍数、低先验知识条件下的DOA估计精度，避免了因角度估计误差引起的目标侦察失误。

3)本发明利用信号源的空域稀疏性得到的稀疏解中大系数对应的角度，为信号源的波达方向角，无目标的方向角对应的系数约等于0，因此无需预先知道目标的数目。

附图说明

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

图2是对本发明与现有基于四阶累量嵌套阵列的虚拟阵元数对比图。

具体实施方式

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

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

步骤1：利用天线接收机形成第一层嵌套阵列，并获取虚拟阵元数。

(1a)定义每个天线接收机为一个阵元，用M个天线接收机形成第一均匀线性阵列a，其阵元间距为d；用N个天线接收机形成第二均匀线性阵列b，其阵元间距为(M+1)d；定义第一均匀线性阵列a的第一个阵元为起始阵元，其中，M≥1,N≥1，0<d≤λ/2，λ为入射到嵌套阵列的窄带信号波长；

(1b)将第一均匀线性阵列a与第二均匀线性阵列b组合为第一层嵌套阵列：

将第二均匀线性阵列b的第一个阵元放置于与起始阵元相距为(M+1)d的位置，再将第二均匀线性阵列b的所有阵元依次插于第一均匀线性阵列a中，得到由M+N个天线接收机形成的第一层嵌套阵列；

(1c)计算第一层嵌套阵列的虚拟阵元数：A＝2MN+2N-1。

步骤2：利用天线接收机形成第二层嵌套阵列。

(2a)用P个天线接收机形成第三均匀线性阵列f，其阵元间距为Ad，用Q个天线接收机形成第四均匀线性阵列h，其阵元间距为(P+1)Ad，其中，Q≥1，P≥1；

(2b)将第三均匀线性阵列f与第四均匀线性阵列h组合为第二层嵌套阵列，即将第三均匀线性阵列f的第一个阵元放置于与起始阵元相距为(A+D)d的位置，将第四均匀线性阵列h的第一个阵元放置于与起始阵元相距为(P+1)Ad的位置，并将第四均匀线性阵列h的所有阵元依次插于第三均匀线性阵列f中，得到由P+Q个天线接收机形成的第二层嵌套阵列，其中D＝(AB-1)/2+MN+N，B＝2PQ+2Q-1。

步骤3：获得第一层嵌套阵列输出信号和第二层嵌套阵列输出信号。

由第一层嵌套阵列和第二层嵌套阵列的天线接收机对空间目标信号进行采样，分别得到第一层嵌套阵列输出信号Y(t)和第二层嵌套阵列输出信号Z(t)：

Y(t)＝[y₁(t),…,y_i(t),…,y_M+N(t)]，

Z(t)＝[z₁(t),…,z_j(t),…,z_P+Q(t)]，

其中，y_i(t)表示第一层嵌套阵列第i个阵元的输出信号，i的取值范围是1≤i≤M+N；z_j(t)表示第二层嵌套阵列的第j个阵元的输出信号，j的取值范围是1≤j≤P+Q。

步骤4：计算第一层嵌套阵列输出信号Y(t)与第二层嵌套阵列输出信号Z(t)的所有四阶累量R，形成四阶累量向量G。

(4a)根据第一层嵌套阵列输出信号Y(t)与第二层嵌套阵列输出信号Z(t)，得到双层嵌套阵列输出信号X(t)，即将第一层嵌套阵列输出信号Y(t)的序列与第二层嵌套阵列输出信号Z(t)的序列依次排列，再对应赋值给双层嵌套阵列输出信号X(t)的序列：

(4b)计算第一层嵌套阵列输出信号Y(t)与第二层嵌套阵列输出信号Z(t)的所有四阶累量R:

R＝c(k1,k2,k3,k4)＝cum(x_k1(t),x_k2 ^*(t),x_k3(t),x_k4 ^*(t))，

其中，c(k1,k2,k3,k4)表示四阶累量，cum表示求解四阶累量操作；k1，k2，k3，k4分别表示双层嵌套阵列第k1，k2，k3，k4个阵元，0≤k1,k2,k₃,k₄≤M+N+P+Q；x_k1(t),x_k2(t)^*,x_k3(t),x_k4(t)^*分别表示双层嵌套阵列第k1，k2，k3，k4个阵元的输出信号，(·)^*表示向量的共轭运算；

(4c)定义中间变量L＝(AB-1)/2+(P+1)(2MN+2N-1)+(MN+N-1),定义四阶累量c(k1,k2,k3,k4)相应的阵元k1,k2,k3,k4分别到起始阵元的距离为d1,d2,d3,d4，依次找出满足d1+d3-d2-d4＝0,1,…,L条件的四阶累量，并将其依次定义为四阶累量元素g(0),g(1),…,g(L)，形成四阶累量向量G：

G＝[g(0),g(1),...,g(L)]。

步骤5：对探测空域波达方向角观测空间进行网格划分，构造出超完备基Φ(θ)，并定义一个Q×1维波达角空域稀疏向量s。

(5a)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成E个角度，定义为波达方向角范围θ＝[θ₁,θ₂,...,θ_e,...,θ_E]，其中θ_e为目标信号的来波方向角，e＝1,2,...,E，E≥M,N,P,Q；

(5b)根据波达方向角范围θ＝[θ₁,θ₂,...,θ_e,...,θ_E]，得到一个空域稀疏化后的超完备基：Φ(θ)＝[w(θ₁),...,w(θ_e),...,w(θ_E)]，其中，w(θ_e)表示方向角θ_e对应的导向矢量：

其中,(·)^T表示矩阵转置运算，j为虚数单位。

步骤6：建立约束关系，将探测空域波达方向角估计转化为求解稀疏约束方程。

根据四阶累量向量G和超完备基Φ(θ)建立约束关系，将探测空域波达方向角估计转化为求解如下稀疏约束方程：

min||s||₁

subject to||G-Φ(θ)s||₂≤η

s≥0

其中，||·||₁表示表示求矩阵一阶范数操作，||·||₂表示求矩阵二阶范数操作，subjectto表示约束关系，η为现实环境中的噪声标准差。

步骤7：利用凸优化方法求解上述稀疏约束方程，得到空域稀疏向量s的最稀疏解

步骤8：绘制幅度谱图。

(8a)以波达方向角范围θ＝[θ₁,θ₂,...,θ_e,...,θ_E]的值作为x轴坐标，以空域稀疏向量s的最稀疏解的幅度值为y轴坐标，绘制幅度谱图；

(8b)设入射到第一层嵌套阵列与第二层嵌套阵列的空间目标信号个数为K，且假设空间目标信号在传播过程中加入了均值为零的复高斯白噪声，K≥1；

(8c)在(8a)绘制出的幅度谱图中，按照从高到低的顺序寻找幅度值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

下面结合仿真实例对本发明的效果做进一步的描述。

1.仿真条件：

用给定总数量的阵元形成二层嵌套阵列，每个子阵阵元数目的不同将影响最终四阶累量虚拟阵元个数。在给定阵元总数的情况下，子阵阵元数目参数均按照孔径扩展能力最优进行选取。阵元数最少的第一层嵌套阵列为M＝1，N＝1，阵元数为2，阵元数最少的第二层嵌套阵列为P＝1，Q＝1，阵元数为2。因此，双层嵌套阵列的阵元数最少为4。

用本发明与现有基于四阶累量嵌套阵列最多所能形成的虚拟阵元数进行比较，结果如图2所示。从图2中可以看出，当阵元总数等于4时，现有基于四阶累量嵌套阵列生成的四阶累量虚拟阵元略高于本发明；但当阵元总数大于等于5时，本发明都能够生成更多的四阶累量虚拟阵元，且随着阵元数的增大，本发明的优势越明显，因此，在阵元数量一定且大于等于5的情况下，本发明可以得到更多的虚拟阵元，从而增加了阵列可识别信源数目。

综上，本发明解决了现有技术阵元利用率低，识别信源数目少，无源定位估计误差大的问题，在有限阵元检测多信号中更加有利，提高了角度估计数量的同时，保证了目标侦察和无源定位的高效性和准确性，可用于目标侦察和无源定位。

Claims (3)

1.一种基于双层嵌套阵列的波达方向角估计方法，包括：

(1)利用天线接收机形成第一层嵌套阵列，并获取虚拟阵元数：

(1a)定义每个天线接收机为一个阵元，用M个天线接收机形成第一均匀线性阵列a，其阵元间距为d；用N个天线接收机形成第二均匀线性阵列b，其阵元间距为(M+1)d；定义第一均匀线性阵列a的第一个阵元为起始阵元，其中，M≥1,N≥1，0＜d≤λ/2，λ为入射到嵌套阵列的窄带信号波长；

(1b)将第一均匀线性阵列a与第二均匀线性阵列b组合为第一层嵌套阵列：

将第二均匀线性阵列b的第一个阵元放置于与起始阵元相距为(M+1)d的位置，再将第二均匀线性阵列b的所有阵元依次插于第一均匀线性阵列a中，得到由M+N个天线接收机形成的第一层嵌套阵列；

(1c)计算第一层嵌套阵列的虚拟阵元数：A＝2MN+2N-1；

(2)利用天线接收机形成第二层嵌套阵列：

(2a)用P个天线接收机形成第三均匀线性阵列f，其阵元间距为Ad，用Q个天线接收机形成第四均匀线性阵列h，其阵元间距为(P+1)Ad，其中，Q≥1，P≥1；

(2b)将第三均匀线性阵列f与第四均匀线性阵列h组合为第二层嵌套阵列，即将第三均匀线性阵列f的第一个阵元放置于与起始阵元相距为(A+D)d的位置，将第四均匀线性阵列h的第一个阵元放置于与起始阵元相距为(P+1)Ad的位置，并将第四均匀线性阵列h的所有阵元依次插于第三均匀线性阵列f中，得到由P+Q个天线接收机形成的第二层嵌套阵列，其中D＝(AB-1)/2+MN+N，B＝2PQ+2Q-1；

(3)由第一层嵌套阵列和第二层嵌套阵列的天线接收机对空间目标信号进行采样，分别得到第一层嵌套阵列输出信号Y(t)和第二层嵌套阵列输出信号Z(t)，其中，

Y(t)＝[y₁(t),…,y_i(t),…,y_M+N(t)]，y_i(t)表示第一层嵌套阵列第i个阵元的输出信号，i的取值范围是1≤i≤M+N；

Z(t)＝[z₁(t),…,z_j(t),…,z_P+Q(t)]，z_j(t)表示第二层嵌套阵列的第j个阵元的输出信号，j的取值范围是1≤j≤P+Q；

(4)计算第一层嵌套阵列输出信号Y(t)与第二层嵌套阵列输出信号Z(t)的所有四阶累量R，形成四阶累量向量G；

(5)对探测空域波达方向角观测空间进行网格划分，构造出超完备基Φ(θ)，并定义一个Q×1维波达角空域稀疏向量s；

(6)根据四阶累量向量G和超完备基Φ(θ)建立约束关系，将探测空域波达方向角估计转化为求解如下稀疏约束方程：

min||s||₁

subject to||G-Φ(θ)s||₂≤η

s≥0

其中，||·||₁表示表示求矩阵一阶范数操作，||·||₂表示求矩阵二阶范数操作，subject to表示约束关系，η为现实环境中的噪声标准差；

(7)利用凸优化方法求解上述稀疏约束方程，得到空域稀疏向量s的最稀疏解

(8)以波达方向角范围θ＝[θ₁,θ₂,...,θ_e,...,θ_E]的值为x轴坐标，以空域稀疏向量s的最稀疏解的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅度值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值，其中，K表示入射到第一层互质阵列与第二层嵌套阵列的空间目标信号个数，且假设空间目标信号在传播过程中加入了均值为零的复高斯白噪声，K≥1；θ_e为目标信号的来波方向角，e＝1,2,...,E,E≥M,N,P,Q。

2.根据权利要求1所述的方法，其中步骤(4)计算所有四阶累量R，形成四阶累量向量G，按如下步骤进行：

(4a)根据第一层嵌套阵列输出信号Y(t)与第二层嵌套阵列输出信号Z(t)，得到双层嵌套阵列输出信号X(t)，即将第一层嵌套阵列输出信号Y(t)的序列与第二层嵌套阵列输出信号Z(t)的序列依次排列，再对应赋值给双层嵌套阵列输出信号X(t)的序列：

(4b)计算第一层嵌套阵列输出信号Y(t)与第二层嵌套阵列输出信号Z(t)的所有四阶累量R:

R＝c(k1,k2,k3,k4)＝cum(x_k1(t),x_k2 ^*(t),x_k3(t),x_k4 ^*(t))，

其中，c(k1,k2,k3,k4)表示四阶累量，cum表示求解四阶累量操作；k1，k2，k3，k4分别表示双层嵌套阵列第k1，k2，k3，k4个阵元，0≤k1,k2,k₃,k₄≤M+N+P+Q；x_k1(t),x_k2(t)^*,x_k3(t),x_k4(t)^*分别表示双层嵌套阵列第k1，k2，k3，k4个阵元的输出信号，(·)^*表示向量的共轭运算；

(4c)定义中间变量L＝(AB-1)/2+(P+1)(2MN+2N-1)+(MN+N-1),定义四阶累量c(k1,k2,k3,k4)相应的阵元k1,k2,k3,k4分别到起始阵元的距离为d1,d2,d3,d4，依次找出满足d1+d3-d2-d4＝0,1,...,L条件的四阶累量，并将其依次定义为四阶累量元素g(0),g(1),...,g(L)，形成四阶累量向量G：

G＝[g(0),g(1),...,g(L)]。

3.根据权利要求1所述的方法，其中步骤(5)对探测空域波达方向角观测空间进行网格划分，构造超完备基Φ(θ)，按如下步骤进行：

(5a)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成E个角度，定义为波达方向角范围θ＝[θ₁,θ₂,...,θ_e,...,θ_E]，θ_e为目标信号的来波方向角，e＝1,2,...,E，E≥M,N,P,Q；

(5b)根据波达方向角范围θ＝[θ₁,θ₂,...,θ_e,...,θ_E]，得到一个空域稀疏化后的超完备基：Φ(θ)＝[w(θ₁),...,w(θ_e),...,w(θ_E)]，其中，w(θ_e)表示方向角θ_e对应的导向矢量：

其中,(·)^T表示矩阵转置运算，j为虚数单位，L为中间变量，定义为L＝(AB-1)/2+(P+1)(2MN+2M-1)+(2MN-N)。