CN106226729B - 基于四阶累量的互质阵列波达方向角估计方法 - Google Patents

Abstract

本发明公开了一种基于四阶累量的互质阵列波达方向角估计方法，主要解决现有技术中阵元利用率低，信号识别数量少的问题，其实现步骤是：1)分别构建一层互质阵列和二层互质阵列；2)获取一层互质阵列的输出信号和二层互质阵列的输出信号；3)根据一层互质阵列输出信号与二层互质阵列输出信号计算所有四阶累量；4)根据所有四阶累量构造四阶累量矩阵；5)计算四阶累量矩阵的噪声子空间；6)根据四阶累量的噪声子空间和阵列流型矩阵计算空间谱；7)根据空间谱绘制幅度谱图，得到波达方向角。本发明在阵元数量有限的情况下大大提高了阵列可识别的信源数目，适用于目标侦察和无源定位。

Description

基于四阶累量的互质阵列波达方向角估计方法

技术领域

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

背景技术

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

针对该问题的研究中，出现较早、应用较为广泛的是多重信号分类MUSIC子空间的模型，对于一个L阵元的典型线性均匀阵列，传统的MUSIC类计算方法可检测的信源数目是L-1个。之后的大部分算法都是利用该模型生成的，例如信号参数估计旋转不变技术ESPRIT。这些算法采用典型的线性均匀阵列，造成估计的信号数目低于阵元数目，目标个数很多时甚至无法识别，导致目标捕获失败。

为了在少的阵元条件下得到尽量大的角度自由度，检测更多的信源，一些新的阵列结构被提出，比较有代表性的是嵌套阵列以及互质阵列。Piya Pal等人在其发表的论文“Coprime Sampling and the MUSIC algorithm”(《Digital signal processingworkshop and IEEE signal processing education workshop》,pp.289-294,2011.)中公开了一种基于互质阵列的DOA估计方法，该方法能够使用N+2M-1个阵元，生成2MN+1个虚拟阵元，可检测MN个信号。该方法具有估计多于阵元数目的信号数的能力，但是，该方法仍然存在的不足之处是，在阵元数量一定的情况下，能估计的信号数量有时仍不能满足实际应用。

为解决上述问题，近年来，出现了使用四阶累量增加阵列虚拟阵元的方法，可使信号估计数量得到较大提升。但已有的基于四阶累量的波达角度估计方法多数是针对均匀阵列，而对于互质阵列的算法，由于没有合理设置阵元位置而不能使阵元达到最大利用率。在实际应用中，给定一定数量的阵元，如果不能合理利用这些阵元获得足够多的虚拟阵元，就不能估计足够多的信号，造成侦察和定位资源的浪费。

发明内容

本发明的目的在于针对上述现有技术存在的不足，提出一种基于四阶累量的互质阵列波达方向角估计方法，以在阵元数量一定的情况下，提高能够进行估计的信号数量，避免因不能合理利用阵元造成的资源浪费。

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

(1)用2M+N-1个天线接收机形成第一层互质阵列：

(1a)将每个天线接收机称为一个阵元，用N个天线接收机形成第一均匀线性阵列a，其阵元间距为Md；用2M-1个天线接收机形成第二均匀线性阵列b，其阵元间距为Nd；定义第一均匀线性阵列a的第一个阵元为第一层互质阵列的阵元0，其中，N＞M≥2且M与N互质，0＜d≤λ/2，λ为入射到阵列的窄带信号波长；

(1b)将第一均匀线性阵列a与第二均匀线性阵列b组合为第一层互质阵列：将第二均匀线性阵列b的第一个阵元放置于与一层互质阵列的阵元0相距为Nd的位置；将第二均匀线性阵列b的所有阵元依次插于第一均匀线性阵列a中，形成第一层互质阵列；

(1c)计算第一均匀线性阵列a与第二均匀线性阵列b形成的第一层互质阵列的虚拟阵元数A＝2MN+1；

(2)用2P+Q-2个天线接收机形成第二层互质阵列：

(2a)用Q-1个天线接收机形成第三均匀线性阵列c，其阵元间距为PAd，用2P-1个天线接收机形成第四均匀线性阵列f，其阵元间距为QAd，其中，Q＞P≥2且P与Q互质；

(2b)将第三均匀线性阵列c与第四均匀线性阵列f组合为第二层互质阵列，即将第三均匀线性阵列c的第一个阵元放置于与一层互质阵列的阵元0相距为PAd的位置，将第四均匀线性阵列f的第一个阵元放置于与一层互质阵列的阵元0相距为QAd的位置，并将第四均匀线性阵列f的所有阵元依次插于第三均匀线性阵列c中，组合为第二层互质阵列；

(2c)计算第三均匀线性阵列c与第四均匀线性阵列f组合的第二层互质阵列的虚拟阵元数B＝2PQ+1；

(3)由第一层互质阵列和第二层互质阵列的天线接收机对空间目标信号进行采样，分别得到第一层互质阵列输出信号Y(t)和第二层互质阵列输出信号Z(t)，其中，Y(t)＝[y₀(t),…,y_i(t),…,y_2M+N-2(t)]，y_i(t)表示第一层互质阵列第i个阵元的输出信号，i的取值范围是0≤i≤2M+N-2，Z(t)＝[z₀(t),…,z_j(t),…,z_2P+Q-2(t)]，z_j(t)表示第二层互质阵列的第j个阵元的输出信号，j的取值范围是0≤j≤2P+Q-2；

(4)计算第一层互质阵列输出信号Y(t)与第二层互质阵列输出信号Z(t)的所有四阶累量c(k1,k2,k3,k4)＝cum(y_k1(t),y_k2(t)^*,z_k3(t),z_k4(t)^*)，其中，0≤k1,k2≤2M+N-2,1≤k3,k4≤2P+Q-2,cum表示求解四阶累量操作，(·)^*表示向量的共轭运算；

(5)定义中间变量D＝(AB-1)/2，从所有四阶累量c(k1,k2,k3,k4)中依次找出k1,k2,k3,k4满足k1+k3A-k2-k4A＝-D,...,-1,0,1,...,D条件的四阶累量，并依次定义为四阶累量元素g(-D),...,g(-1),g(0),g(1),...,g(D)，将这些四阶累量元素重新排列形成四阶累量矩阵G：

(6)计算四阶累量矩阵G的噪声子空间U_n；

(7)根据阵列流型矩阵α(θ)和四阶累量矩阵G的噪声子空间U_n，计算空间谱P(θ)；

(8)以波达方向角范围θ的值为x轴坐标，以空间谱P(θ)的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为目标的波达方向角度值，其中，K表示入射到一层互质阵列与二层互质阵列的空间目标信号个数，K≥1。

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

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

2)本发明将四阶累量应用到互质阵列的DOA估计中，通过使用四阶累量，设置阵元位置，可使用2M+N+2P+Q-3个阵元获得(2MN+1)(2PQ+1)个连续虚拟阵元，大大提高了阵列利用率，同时通过获得更多的虚拟阵元，进一步增加了阵列可识别的信源数目。

附图说明

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

图2是本发明中一层互质阵列与二层互质阵列的结构示意图。

具体实施方式

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

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

步骤1：用2M+N-1个天线接收机形成第一层互质阵列，得到第一层互质阵列虚拟阵元数A。

(1a)将每个天线接收机称为一个阵元，用N个天线接收机形成第一均匀线性阵列a，其阵元间距为Md；用2M-1个天线接收机形成第二均匀线性阵列b，其阵元间距为Nd；定义第一均匀线性阵列a的第一个阵元为第一层互质阵列的阵元0，其中，N＞M≥2且M与N互质，0＜d≤λ/2，λ为入射到阵列的窄带信号波长；

(1b)将第一均匀线性阵列a与第二均匀线性阵列b组合为一层互质阵列：将第二均匀线性阵列b的第一个阵元放置于与一层互质阵列的阵元0相距为Nd的位置；将第二均匀线性阵列b的所有阵元依次插于第一均匀线性阵列a中，形成第一层互质阵列；

(1c)计算第一均匀线性阵列a与第二均匀线性阵列b形成的第一层互质阵列的虚拟阵元数A＝2MN+1。

步骤2：用2P+Q-2个天线接收机形成第二层互质阵列，得到第二层互质阵列的虚拟阵元数B。

(2a)用Q-1个天线接收机形成第三均匀线性阵列c，其阵元间距为PAd，用2P-1个天线接收机形成第四均匀线性阵列f，其阵元间距为QAd，其中，Q＞P≥2且P与Q互质；

(2b)将第三均匀线性阵列c与第四均匀线性阵列f组合为第二层互质阵列，即将第三均匀线性阵列c的第一个阵元放置于与第一层互质阵列的阵元0相距为PAd的位置，将第四均匀线性阵列f的第一个阵元放置于与第一层互质阵列的阵元0相距为QAd的位置，并将第四均匀线性阵列f的所有阵元依次插于第三均匀线性阵列c中，组合为第二层互质阵列；

(2c)计算第三均匀线性阵列c与第四均匀线性阵列f组合的第二层互质阵列的虚拟阵元数B＝2PQ+1。

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

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

其中，Y(t)＝[y₀(t),…,y_i(t),…,y_2M+N-2(t)]，y_i(t)表示第一层互质阵列第i个阵元的输出信号，i的取值范围是0≤i≤2M+N-2，Z(t)＝[z₀(t),…,z_j(t),…,z_2P+Q-2(t)]，z_j(t)表示第二层互质阵列的第j个阵元的输出信号，j的取值范围是0≤j≤2P+Q-2。

步骤4：计算输出信号的所有四阶累量。

根据第一层互质阵列输出信号Y(t)和第二层互质阵列输出信号Z(t)计算所有四阶累量：c(k1,k2,k3,k4)＝cum(y_k1(t),y_k2(t)^*,z_k3(t),z_k4(t)^*)，

其中，0≤k1,k2≤2M+N-2,1≤k3,k4≤2P+Q-2,cum表示求解四阶累量操作，(·)^*表示向量的共轭运算。

步骤5：构造四阶累量矩阵G。

定义中间变量D＝(AB-1)/2，从所有四阶累量c(k1,k2,k3,k4)中依次找出k1,k2,k3,k4满足k1+k3A-k2-k4A＝-D,...,-1,0,1,...,D条件的四阶累量，并依次定义为四阶累量元素g(-D),...,g(-1),g(0),g(1),...,g(D)，将这些四阶累量元素重新排列形成四阶累量矩阵G：

步骤6：计算四阶累量矩阵的噪声子空间U_n。

(6a)对四阶累量矩阵G进行如下特征分解：

G＝U·Λ·U^H，

其中，Λ为四阶累量矩阵G的特征值矩阵，U为四阶累量矩阵G的特征值所对应的特征向量矩阵，(·)^H表示矩阵的共轭转置远算；

(6b)将特征值矩阵Λ中的特征值按从大到小排序，取后D-K个较小特征值对应的特征向量矩阵作为噪声子空间U_n。

步骤7：计算空间谱P(θ)。

(7a)将观测空域[-90°,90°]按等间隔L划分成E个角度，定义波达方向角范围为θ＝[θ₁,θ₂,…,θ_e,…θ_E]，

其中，θ_e表示第e个目标波达方向角，e＝1,2,...,E，E＞＞max{M,N,P,Q}；L的取值根据期望达到的角度估计精度进行设定，网格划分间隔越小则最终得到的角度估计值精度越高；

(7b)计算第e个目标波达方向角θ_e对应的导向矢量α(θ_e)：

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

(7c)计算所有的目标波达方向角对应的导向矢量，得到阵列流型矩阵α(θ)：

α(θ)＝[α(θ₁),…,α(θ_e),…,α(θ_E)]；

(7d)根据阵列流型矩阵α(θ)和四阶累量矩阵G的噪声子空间U_n，计算空间谱P(θ)：

其中，(·)^H表示矩阵的共轭转置远算。

步骤8：绘制幅度谱图。

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

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

仿真实例1.构建双层互质阵列结构。

首先，利用第一均匀线性阵列a与第二均匀线性阵列b形成第一层互质阵列，设N＝2,M＝3，第一均匀线性阵列a含有3个阵元，阵元间距为2d，第二均匀线性阵列b含有3个阵元，阵元间距为3d，第一均匀线性阵列a的第一个阵元放置在0d位置；

其次，计算第一均匀线性阵列a与第二均匀线性阵列b形成的第一层互质阵列虚拟阵元数A＝2MN+1＝13；

然后，利用第三均匀线性阵列c与第四均匀线性阵列f形成第二层互质阵列，设P＝2,Q＝3，第三均匀线性阵列c含有2个阵元，阵元间距为PA＝26d，第四均匀线性阵列f含有3个阵元，阵元间距为QA＝39d，第三均匀线性阵列c的第一个阵元在26d的位置，计算虚拟阵元数B＝2PQ+1＝13。

其中第一均匀线性阵列a、第二均匀线性阵列b和第三均匀线性阵列c、第四均匀线性阵列f形成了双层互质阵列，其阵列结构图如图2所示。

由图2可见，互质阵列结构在阵元数量一定的情况下，可以得到更多的阵元位置信息，从而增加阵列可识别信源数目。

仿真实例2，根据四阶累量计算双层互质阵列所有虚拟阵元。

设第一均匀线性阵列a与第二均匀线性阵列b形成的第一层互质阵列阵元位置为[0，2，3，4，6，9]d；

设第三均匀线性阵列c与第四均匀线性阵列f形成的第二层互质阵列阵元位置为[26，39，52，78，117]d；

计算得到第一层互质阵列每个阵元分别与第二层互质阵列中每个阵元形成的虚拟阵列为[26，28，29，30，32，35]d、[39，41，42，43，45，48]d、[52，54，55，56，58，61]d、[78，80，81，82，84，87]d、[117，119，120，121，123，126]d；

根据四阶累量计算的方法和原则，用上述这些虚拟阵列生成最终的所有虚拟阵元，结果如表1。

表1 双层互质阵列所有虚拟阵元

从表1可见，该双层互质阵列用2M+N+2P+Q-3＝11个原始阵元生成了一个从-84d到84d共含有AB＝169个虚拟连续阵元的线性阵列。

综上，本发明主要解决了现有技术阵元利用率低，识别信源数目少，无源定位估计误差大的问题，提高了一定阵元数情况下阵列可识别的信源数目以及低信噪比下对信号方向角的估计性能。

本发明在有限阵元检测多信号中更加有利，在提高角度估计数量的同时，保证了目标侦察和无源定位的高效性和准确性，在多目标识别上表现出突出的性能。

Claims (4)

1.一种基于四阶累量的互质阵列的波达方向角估计方法，包括：

(1)用2M+N-1个天线接收机形成第一层互质阵列：

(1a)将每个天线接收机称为一个阵元，用N个天线接收机形成第一均匀线性阵列a，其阵元间距为Md；用2M-1个天线接收机形成第二均匀线性阵列b，其阵元间距为Nd；定义第一均匀线性阵列a的第一个阵元为第一层互质阵列的阵元0，其中，N＞M≥2且M与N互质，0＜d≤λ/2，λ为入射到阵列的窄带信号波长；

(1b)将第一均匀线性阵列a与第二均匀线性阵列b组合为第一层互质阵列：将第二均匀线性阵列b的第一个阵元放置于与第一层互质阵列的阵元0相距为Nd的位置；将第二均匀线性阵列b的所有阵元依次插于第一均匀线性阵列a中，形成第一层互质阵列；

(1c)计算第一均匀线性阵列a与第二均匀线性阵列b形成的第一层互质阵列的虚拟阵元数A＝2MN+1；

(2)用2P+Q-2个天线接收机形成第二层互质阵列：

(2a)用Q-1个天线接收机形成第三均匀线性阵列c，其阵元间距为PAd，用2P-1个天线接收机形成第四均匀线性阵列f，其阵元间距为QAd，其中，Q＞P≥2且P与Q互质；

(2b)将第三均匀线性阵列c与第四均匀线性阵列f组合为第二层互质阵列，即将第三均匀线性阵列c的第一个阵元放置于与第一层互质阵列的阵元0相距为PAd的位置，将第四均匀线性阵列f的第一个阵元放置于与第一层互质阵列的阵元0相距为QAd的位置，并将第四均匀线性阵列f的所有阵元依次插于第三均匀线性阵列c中，组合为第二层互质阵列；

(2c)计算第三均匀线性阵列c与第四均匀线性阵列f组合的第二层互质阵列的虚拟阵元数B＝2PQ+1；

(3)由第一层互质阵列和第二层互质阵列的天线接收机对空间目标信号进行采样，分别得到第一层互质阵列输出信号Y(t)和第二层互质阵列输出信号Z(t)，其中，Y(t)＝[y₀(t),…,y_i(t),…,y_2M+N-2(t)]，y_i(t)表示第一层互质阵列第i个阵元的输出信号，i的取值范围是0≤i≤2M+N-2，Z(t)＝[z₀(t),…,z_j(t),…,z_2P+Q-2(t)]，z_j(t)表示第二层互质阵列的第j个阵元的输出信号，j的取值范围是0≤j≤2P+Q-2；

(4)计算第一层互质阵列输出信号Y(t)与第二层互质阵列输出信号Z(t)的所有四阶累量c(k1,k2,k3,k4)＝cum(y_k1(t),y_k2(t)^*,z_k3(t),z_k4(t)^*)，其中，0≤k1,k2≤2M+N-2,1≤k3,k4≤2P+Q-2,cum表示求解四阶累量操作，(·)^*表示向量的共轭运算；

(5)定义中间变量D＝(AB-1)/2，从所有四阶累量c(k1,k2,k3,k4)中依次找出k1,k2,k3,k4满足k1+k3A-k2-k4A＝-D,...,-1,0,1,...,D条件的四阶累量，并依次定义为四阶累量元素g(-D),...,g(-1),g(0),g(1),...,g(D)，将这些四阶累量元素重新排列形成四阶累量矩阵G：

(6)计算四阶累量矩阵G的噪声子空间U_n；

(7)根据阵列流型矩阵α(θ)和四阶累量矩阵G的噪声子空间U_n，计算空间谱P(θ)；

(8)以波达方向角范围θ的值为x轴坐标，以空间谱P(θ)的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为目标的波达方向角度值，其中，K表示入射到第一层互质阵列与第二层互质阵列的空间目标信号个数，K≥1。

2.根据权利要求1所述的一种基于四阶累量的互质阵列的波达方向角估计方法，其中步骤(6)中计算四阶累量矩阵G的噪声子空间U_n，按如下步骤进行：

(6a)对四阶累量矩阵G进行如下特征分解：

G＝U·Λ·U^H，

其中，Λ为四阶累量矩阵G的特征值矩阵，U为四阶累量矩阵G的特征值所对应的特征向量矩阵，(·)^H表示矩阵的共轭转置运算；

(6b)将特征值矩阵Λ中的特征值按从大到小排序，取后D-K个较小特征值对应的特征向量矩阵作为噪声子空间U_n。

3.根据权利要求1所述的一种基于四阶累量的互质阵列的波达方向角估计方法，其中步骤(7)中的阵列流型矩阵α(θ)，按如下步骤构造：

(7a)将观测空域[-90°,90°]等间隔划分成E个角度,定义为波达方向角范围θ＝[θ₁,θ₂,…,θ_e,…θ_E]，其中，θ_e表示第e个目标波达方向角，e＝1,2,...,E，E＞＞max{M,N,P,Q}；

(7b)计算第e个目标波达方向角θ_e对应的导向矢量α(θ_e)：

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

(7c)计算所有的目标波达方向角对应的导向矢量，得到阵列流型矩阵α(θ)：

α(θ)＝[α(θ₁),…,α(θ_e),…,α(θ_E)]。

4.根据权利要求1所述的一种基于四阶累量的互质阵列的波达方向角估计方法，其中步骤(7)中计算空间谱P(θ)，按如下形式进行：

其中，α(θ)表示阵列流型矩阵，(·)^H表示矩阵的共轭转置运算。