CN112014792B - 一种改进的二维互质阵列doa估计方法 - Google Patents

Abstract

本发明公开了一种基于L型改进互质阵列的二维DOA估计方法，利用互质阵列的结构优势实现高性能DOA估计。L型均匀阵列的阵元间距受到半波长的限制，阵元间距很小，阵元间存在严重的互耦效应，并且阵列孔径较小从而制约了角度分辨力，本发明提出的改进的压缩平移互质阵列DOA估计方法，先计算了x轴和z轴子阵接收数据的协方差矩阵，并对其进行矢量化和去冗余处理，利用处理后的数据构造Toeplize矩阵，最终求得俯仰角和方位角的估计值。本发明给出的改进的压缩平移互质L型阵列，通过对阵元接收数据进行自差和互差运算可以得到虚拟阵元位置连续的均匀完整满布阵列，该阵列相比于现有阵列有更大的阵列孔径，更大的阵列自由度，更高的角度分辨率，应用价值较高。

Description

一种改进的二维互质阵列DOA估计方法

技术领域

本发明属于阵列信号处理技术领域，尤其涉及一种新型二维互质阵列DOA估计方法。

背景技术

一维DOA估计假定信号是从某一个平面入射的，这种假设限定了一维DOA估计算法的应用范围，二维DOA可以打破上述约束条件，对从三维空间的任意方向入射的信号进行准确刻画，因此研究二维DOA估计具有重要意义。常用来进行二维DOA估计的天线阵列结构包括L型均匀阵列、均匀面阵、平行阵列等。目前基于均匀阵列的二维到达角估计算法已经日渐成熟，而均匀阵列具有一定的局限性和不足。相对于均匀阵列，非均匀阵列结构则更加灵活，且具有更好的估计性能，因此基于非均匀阵列的二维DOA估计研究越来越多，本发明研究了基于L型改进互质阵列的二维DOA估计，利用互质阵列的结构优势实现高性能的DOA估计。L型均匀阵列的阵元间距受到半波长的限制，用于高频信号DOA估计的均匀阵列的阵元间距很小，存在严重的互耦效应，且阵列孔径小从而制约了角度分辨力，因此阵元间距稀疏的非均匀二维阵列受到了越来越多的关注。L型互质阵列即是互质阵列与传统L型阵列相结合。互质阵列结构有常规互质阵列、扩展互质阵列、压缩互质阵列、移位互质阵列等等，常规互质阵列由两个均匀阵列组成，两个子阵的阵元相互交叉放置，压缩平移互质阵列也是由两个均匀阵列组成，但是两个子阵的阵元并不是相互交叉放置，而是将一个子阵整体转移到另一个子阵之后，本发明给出了一种改进的压缩平移互质L型阵列，对阵元接收数据进行自差和互差运算可以得到虚拟阵元位置连续的均匀完整满布阵列，该阵列相比于现有阵列有更大的阵列孔径，更大的阵列自由度，更高的角度分辨率。

发明内容

本发明的目的是提供一种新型二维互质阵列DOA估计方法。

为了实现上述目的，本发明采取如下的技术解决方案：

利用设计的新型互质阵列作为接收阵列接收来自K个不同方向的窄带、非高斯的平稳远场入射信号，目标信号的二维到达角可表示为

其中θ_k∈[0,π]表示第k个来波信号的俯仰角，/>

表示第k个来波信号的方位角。

一种改进的二维互质阵列DOA估计方法步骤如下：

步骤一、布置改进的L型压缩互质阵列。

该L型阵列由结构完全相同的x轴和z轴子阵构成，下面以x轴子阵为例来说明，x轴子阵由两个子阵构成，其中子阵1是阵元间距为d的M个阵元构成的均匀阵列，子阵2是由N个阵元构成的一个非均匀阵列，子阵2最后两个阵元间距为Md，其它阵元的间距为Md+λ/2，两个子阵之间存在(M+1)d的位移，x轴坐标上的阵元可用集合D_x表示为：{0,d,…,(M-1)d，2Md,3Md+λ/2,…,2Md+(Md+λ/2)(N-3),3Md+(Md+λ/2)(N-3)}；

步骤二、利用L型互质阵列接收信号计算x轴和z轴子阵数据协方差矩阵R_X和R_Z。

x轴上的M+N-1个阵元接收的信号的T次采样数据为X，可以得到x轴子阵数据协方差矩阵

z轴上的M+N-1个阵元接收的信号的T次采样数据为Z，可以得到z轴子阵数据协方差矩阵/>

步骤三、对x轴和z轴子阵数据协方差矩阵R_X和R_Z分别进行矢量化和去冗余处理得到无冗余数据

和/>

首先矢量化协方差矩阵R_X得到一个列矢量r_x＝vec(R_X)，vec(·)表示将矩阵中的元素按行排成一个列矢量，例如设矩阵Β大小为3×2维，则其内部经过矢量化操作vec(B)后，将原尺寸3×2维拉伸至6×1维，矢量化后具体的表现形式为vec(B)＝[b₁₁,b₁₂,b₁₃,b₂₁,b₂₂,b₂₃]^T，D_x中的元素两两作差得到的结果存在相同元素时，协方差矩阵中存在冗余元素，记阵元位置为u的虚拟阵元对应元素为r_u，冗余度为p(r_u)，虚拟阵元接收数据可由所有冗余元素相加再平均得到，即

其中p为阵元位置为u的冗余元素数，r_ui表示位置为u的阵元对应的第i个元素。则虚拟接收信号可由得到的结果按阵元位置进行重新排序得到。最终得到一个无冗余且有序的接收数据，记作/>

例如：阵列结构阵元位置为{0,d,2d,6d,10d,14d,17d}，得到的虚拟阵元个数为35个，其阵元位置坐标可记为：U＝{-17,-35+1,…,-1,0,1,…,17-1,17}

按照同样的方法，矢量化协方差矩阵r_z＝vec(R_z)，再对r_z中元素进行上述冗余处理，则虚拟接收信号可由得到的结果按阵元位置进行重新排序得到。最终得到一个无冗余且有序的接收数据，记作

步骤四、分别利用接收数据

和/>

构造Toeplize矩阵T_x和T_z。

首先利用接收数据

构造一个Toeplize矩阵：

其中

表示/>

中对应阵元位置为p的元素，p∈U表示虚拟阵元位置坐标。

利用接收数据

构造一个Toeplize矩阵：

其中

表示/>

中对应阵元位置为p的元素，p∈U表示虚拟阵元位置坐标。/>

步骤五、俯仰角的估计

以及特征矢量矩阵T₁。

对T_z进行特征值分解，得到信号子空间U_S，并取U_S的前m-1行元素构成信号子空间U_S1，取U_S的后m-1行元素构成信号子空间U_S2；由U_S1和U_S2计算得到

对Ψ_z进行特征值分解，得到的特征值对应于Φ_z矩阵对角线上元素的值，特征值矢量构成的特征矢量矩阵T₁；利用Φ_z矩阵对角线上的元素即可求解各个目标信号的俯仰角的估计值/>

其中

第k个对角元素记为γ_zk，则入射角度的估计值/>

可表示为：/>

步骤六、方位角的估计

以及特征矢量矩阵T₂。

对T_x进行特征值分解，得到信号子空间U_xS，并取U_xS的前m-1行元素构成信号子空间U_xS1，取U_xS的后m-1行元素构成信号子空间U_xS2；由U_xS1和U_xS2计算得到

对Ψ_x进行特征值分解，得到的特征值对应于Φ_x矩阵对角线上元素的值，特征值矢量构成的特征矢量矩阵T₂；利用Φ_z矩阵对角线上的元素即可求解各个目标信号的方位角的估计值/>

其中/>

第k个对角线元素记为γ_xk，则方位角的估计值/>

可表示为/>

步骤七、俯仰角的估计值

和方位角的估计值/>

进行配对运算。

利用步骤五和六中的两个矩阵特征矢量矩阵T₁和T₂,T₁由Ψ_z特征分解的特征向量组成，T₂由Ψ_x特征分解的特征向量组成。由其计算得到一个矩阵

找出矩阵G每一行里绝对值最大的那一个元素，利用其对应的坐标来对Φ_z和Φ_x中元素的顺序进行调整，形成一一匹配。

前述步骤中，K表示信号源数目，k＝1,2,...,K表示信号源的标号，i＝1,2,...P表示重复阵元个数，T表示采样次数。

本发明提出了一种新型的压缩平移互质阵列下的二维DOA估计算法，本发明所提阵列结构通过协方差矩阵的矢量化，可以得到一个满布阵列，该阵列无阵列空洞，且具有更大的有效孔径，因此改进的压缩平移L型互质阵列具有更大的阵列自由度和更高的角度分辨能力。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中需要使用的附图做简单介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为常规L型互质阵列结构；

图2为压缩平移互质阵列结构图；

图3为最优配置阵列结构；

图4为本发明阵列结构图；

图5为本发明的流程图；

图6为本发明方法6个入射信号谱峰图；

图7为二维俯仰角估计均方根误差随信噪比变化；

图8为二维方位角估计均方根误差随信噪比变化；

具体实施方式

为了让本发明的上述和其它目的、特征及优点能更明显，下文特举本发明实施例，并配合所附图示，做详细说明如下。

本发明的目的是提供一种新型二维互质阵列DOA估计方法。

为了实现上述目的，本发明采取如下的技术解决方案：

利用设计的新型互质阵列作为接收阵列接收来自K个不同方向的窄带、非高斯的平稳远场入射信号，目标信号的二维到达角可表示为

其中θ_k∈[0,π]表示第k个来波信号的俯仰角，/>

表示第k个来波信号的方位角。

一种改进的二维互质阵列DOA估计方法步骤如下：

步骤一、布置改进的L型压缩互质阵列。

该L型阵列由结构完全相同的x轴和z轴子阵构成，下面以x轴子阵为例来说明，x轴子阵由两个子阵构成，其中子阵1是阵元间距为d的M个阵元构成的均匀阵列，子阵2是由N个阵元构成的一个非均匀阵列，子阵2最后两个阵元的间距为Md，其它阵元的间距为Md+λ/2，两个子阵之间存在(M+1)d的位移，x轴坐标上的阵元位置坐标可用集合D_x表示为：{0,d,…,(M-1)d，2Md,3Md+λ/2,…,2Md+(Md+λ/2)(N-3),3Md+(Md+λ/2)(N-3)}；

步骤二、利用L型互质阵列接收信号计算x轴和z轴子阵数据协方差矩阵R_X和R_Z。

x轴上的M+N-1个阵元接收的信号的T次采样数据为X，可以得到x轴子阵数据协方差矩阵

z轴上的M+N-1个阵元接收的信号的T次采样数据为Z，可以得到z轴子阵数据协方差矩阵/>

步骤三、对x轴和z轴子阵数据协方差矩阵R_X和R_Z分别进行矢量化和去冗余处理得到无冗余数据

和/>

首先矢量化协方差矩阵R_X得到一个列矢量r_x＝vec(R_X)，vec(·)表示将矩阵中的元素按行排成一个列矢量，例如设矩阵Β大小为3×2维，则其内部经过矢量化操作vec(B)后，将原尺寸3×2维拉伸至6×1维，矢量化后具体的表现形式为vec(B)＝[b₁₁,b₁₂,b₁₃,b₂₁,b₂₂,b₂₃]^T，D_x中的元素两两作差得到的结果相同时，列矢量r_x中存在冗余元素，记阵元位置为u的虚拟阵元对应元素为r_u，冗余度为p(r_u)，虚拟阵元接收数据可由所有冗余元素相加再平均得到，即

其中r_ui表示位置为u的阵元对应的第i个元素，P为位置为u的阵元对应的重复冗余元素个数。则虚拟接收信号可由得到的结果按阵元位置进行重新排序得到。最终得到一个无冗余且有序的接收数据，记作/>

例如：阵列结构阵元位置为{0,d,2d,6d,10d,14d,17d}，得到的虚拟阵元个数为35个，其阵元位置坐标可记为：U＝{-17,-35+1,…,-1,0,1,…,17-1,17}

按照同样的方法，矢量化协方差矩阵r_z＝vec(R_z)，再对r_z中元素进行上述冗余处理，则虚拟接收信号可由得到的结果按阵元位置进行重新排序得到。最终得到一个无冗余且有序的接收数据，记作

步骤四、分别利用接收数据

和/>

构造Toeplize矩阵T_x和T_z。

首先利用接收数据

构造一个Toeplize矩阵：

其中

表示/>

中对应阵元位置为p的元素，p∈U表示虚拟阵元位置坐标。

利用接收数据

构造一个Toeplize矩阵：

其中

表示/>

中对应阵元位置为p的元素，p∈U表示虚拟阵元位置坐标。

步骤五、俯仰角的估计

以及特征矢量矩阵T₁。

对T_z进行特征值分解，得到信号子空间U_S，并取U_S的前m-1行元素构成信号子空间U_S1，取U_S的后m-1行元素构成信号子空间U_S2；由U_S1和U_S2计算得到

对Ψ_z进行特征值分解，得到的特征值对应于Φ_z矩阵对角线上元素的值，特征值矢量构成的特征矢量矩阵T₁；利用Φ_z矩阵对角线上的元素即可求解各个目标信号的俯仰角的估计值/>

其中

第k个对角元素记为γ_zk，则入射角度的估计值

可表示为：/>

步骤六、方位角的估计

以及特征矢量矩阵T₂。

对T_x进行特征值分解，得到信号子空间U_xS，并取U_xS的前m-1行元素构成信号子空间U_xS1，取U_xS的后m-1行元素构成信号子空间U_xS2；由U_xS1和U_xS2计算得到

对Ψ_x进行特征值分解，得到的特征值对应于Φ_x矩阵对角线上元素的值，特征值矢量构成的特征矢量矩阵T₂；利用Φ_z矩阵对角线上的元素即可求解各个目标信号的方位角的估计值/>

其中/>

第k个对角线元素记为γ_xk，则方位角的估计值/>

可表示为/>

步骤七、俯仰角的估计值

和方位角的估计值/>

进行配对运算。

利用步骤五和六中的两个矩阵特征矢量矩阵T₁和T₂,T₁由Ψ_z特征分解的特征向量组成，T₂由Ψ_x特征分解的特征向量组成。由其计算得到一个矩阵

找出矩阵G每一行里绝对值最大的那一个元素，利用其对应的坐标来对Φ_z和Φ_x中元素的顺序进行调整，形成一一匹配。

前述步骤中，K表示信号源数目，k＝1,2,...,K表示信号源的标号，i＝1,2,...P表示重复阵元个数，T表示采样次数。

本发明所提阵列结构通过协方差矩阵的矢量化，可以得到一个满布阵列，该阵列无阵列空洞，且具有更大的有效孔径，因此改进的压缩平移L型互质阵列具有更大的阵列自由度和更高的角度分辨能力。

本发明的效果可以通过以下的仿真结果进一步说明：

仿真实验一：该仿真实验主要是对不同阵列结构DOA估计的自由度进行仿真实验。分别以常规互质阵列、压缩互质阵列和本发明研究的阵列作为接收数据的阵列，参数M＝3、N＝4,一共6个物理阵元。入射信号的俯仰角在30°～150°之间，信噪比为20dB,快拍数为1000，图6为不同阵列结构基于子空间类算法得到的俯仰角估计谱峰图。

图6是入射信号为6个时不同阵列结构的DOA估计谱峰图，从图中可以看出当入射信号个数为6个时，本发明所提出的阵列结构可以有效估计出角度，性能优于其它算法。

仿真实验二：该仿真实验主要是研究不同信噪比下算法的均方根误差变化情况，同样将本发明方法与其它算法做比较。L型阵列结构参数取M＝3,N＝4，分别为本发明提出的改进L型互质阵列和常规L型互质阵列。接收来自两个不同方向的目标信号，两个信号的二维到达角分别为(23°,50°)和(70°，111°)，快拍数T＝1000,信噪比的值在-5dB～25dB区间内变化，每5dB进行100次实验。仿真结果如图7、图8所示。

图7、图8为本发明方法与其它方法角度估计均方根误差对比图，其中图7为连续信噪比下俯仰角估计均方根误差变化，图8为连续信噪比下方位角估计均方根误差变化。可以看出随着信噪比的增大，各个方法的均方根误差都在减小。而本发明方法在信噪比较低时均方根误差最小，说明估计性能较好。这是因为本发明方法是基于改进的阵列结构进行的DOA估计，改进后的阵列相比于一般互质阵列孔径增大了许多，因此有较好的信号分辨力，而且能够形成一个完整的虚拟阵列，不会有孔洞造成信息损失，形成的虚拟阵元个数也比常规互质阵列多。而其它方法都是基于常规互质阵列进行的DOA估计，阵列自身孔径较小而且形成的虚拟阵列存在孔洞，MR算法虽然能够弥补孔洞，但是由于常规互质阵列本身孔径小，因此形成的虚拟阵列最大孔径也是有限的，所以较本发明方法的性能还是差一些。

以上所述，仅是本发明的较佳实施例而已，并非对本发明做任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (1)

1.一种改进的二维互质阵列DOA估计方法，其特征在于：

一种改进的二维互质阵列DOA估计方法，利用设计的新型互质阵列作为接收阵列接收来自K个不同方向的窄带、非高斯的平稳远场入射信号，目标信号的二维到达角可表示为

其中θ_k∈[0,π]表示第k个来波信号的俯仰角，/>

表示第k个来波信号的方位角；

一种改进的二维互质阵列DOA估计方法，步骤如下：

步骤一、布置改进的L型压缩互质阵列；

改进的L型压缩互质阵列由结构完全相同的x轴和z轴子阵构成，下面以x轴子阵为例来说明，x轴子阵由两个子阵构成，其中子阵1是阵元间距为d的M个阵元构成的均匀阵列，子阵2是由N个阵元构成的一个非均匀阵列，子阵2最后两个阵元的间距为Md，其它阵元的间距为Md+λ/2，两个子阵之间存在(M+1)d的位移，x轴坐标上的阵元位置可用集合D_x表示为：{0,d,…,(M-1)d，2Md,3Md+λ/2,…,2Md+(Md+λ/2)(N-3),3Md+(Md+λ/2)(N-3)}；

步骤二、利用L型互质阵列接收信号计算x轴和z轴子阵数据协方差矩阵R_X和R_Z；

x轴上的M+N-1个阵元接收的信号的T次采样数据为X，可以得到x轴子阵数据协方差矩阵

z轴上的M+N-1个阵元接收的信号的T次采样数据为Z，可以得到z轴子阵数据协方差矩阵/>

步骤三、对x轴和z轴子阵数据协方差矩阵R_X和R_Z分别进行矢量化和去冗余处理得到无冗余数据

和/>

首先矢量化协方差矩阵R_X得到一个列矢量r_x＝vec(R_X)，vec(·)表示将矩阵中的元素按行排成一个列矢量；D_x两两作差得到的结果相同，协方差矩阵中存在冗余元素，记阵元位置为u的虚拟阵元对应元素为r_u，冗余度为p(r_u)，虚拟阵元接收数据可由所有冗余元素相加再平均得到，即

p为位置为u的重复冗余元素的个数，r_ui表示位置为u的阵元对应的第i个元素，则虚拟接收信号可由得到的结果按阵元位置进行重新排序得到，最终得到一个无冗余且有序的接收数据，记作/>

阵元位置为{0,d,2d,6d,10d,14d,17d}，得到的虚拟阵元个数为35个，其阵元位置坐标可记为：U＝{-17,-35+1,...,-1,0,1,...,17-1,17}；

按照同样的方法，矢量化协方差矩阵r_z＝vec(R_z)，再对r_z中元素进行上述冗余处理，则虚拟接收信号可由得到的结果按阵元位置进行重新排序得到，最终得到一个无冗余且有序的接收数据，记作

步骤四、分别利用接收数据

和/>

构造Toeplize矩阵T_x和T_z；

首先利用接收数据

构造一个Toeplize矩阵：

/>

其中

表示/>

中对应阵元位置为p的元素，p∈U表示虚拟阵元位置坐标；

利用接收数据

构造一个Toeplize矩阵：

其中

表示/>

中对应阵元位置为p的元素，p∈U表示虚拟阵元位置坐标；

步骤五、俯仰角的估计

以及特征矢量矩阵T₁；

对T_z进行特征值分解，得到信号子空间U_S，并取U_S的前m-1行元素构成信号子空间U_S1，取U_S的后m-1行元素构成信号子空间U_S2，由U_S1和U_S2计算得到

对Ψ_z进行特征值分解，得到的特征值对应于Φ_z矩阵对角线上元素的值，特征值矢量构成的特征矢量矩阵T₁，利用Φ_z矩阵对角线上的元素即可求解各个目标信号的俯仰角的估计值/>

其中

diag[·]表示以括号中元素为对角元素的对角矩阵，Φ_z的第k个对角元素记为γ_zk，则入射角度的估计值/>

可表示为：

λ表示入射信号的波长；

步骤六、方位角的估计

以及特征矢量矩阵T₂；

对T_x进行特征值分解，得到信号子空间U_xS，并取U_xS的前m-1行元素构成信号子空间U_xS1，取U_xS的后m-1行元素构成信号子空间U_xS2；由U_xS1和U_xS2计算得到

对Ψ_x进行特征值分解，得到的特征值对应于Φ_x矩阵对角线上元素的值，特征值矢量构成的特征矢量矩阵T₂；利用Φ_x矩阵对角线上的元素即可求解各个目标信号的方位角的估计值/>

其中

第k个对角线元素记为γ_xk，则方位角的估计值/>

可表示为/>

步骤七、俯仰角的估计值

和方位角的估计值/>

进行配对运算；

利用步骤五和六中的两个矩阵特征矢量矩阵T₁和T₂,T₁由Ψ_z特征分解的特征向量组成，T₂由Ψ_x特征分解的特征向量组成，由其计算得到一个矩阵

找出矩阵G每一行里绝对值最大的那一个元素，利用其对应的坐标来对Φ_z和Φ_x中元素的顺序进行调整，形成一一匹配；

前述步骤中，K表示信号源数目，k＝1,2,...,K表示信号源的标号，i＝1,2,...P表示重复阵元个数，T表示采样次数。

Images