CN114397620A - 改进和差非均匀阵列的高精度波达方向估计方法 - Google Patents

Abstract

本发明为了充分利用非均匀阵列的和差共阵信息，公开了一种改进和差非均匀阵列的高精度DOA估计方法；本发明能够在阵元个数有限的情况下估计更多的信号源，并且具有较高的估计精度；具体实施方式如下：首先用两个均匀子阵构成改进的一维非均匀阵列，对其做求和、求差运算求得虚拟阵元总的位置坐标集合；然后用改进的和差非均匀阵列接收K个非相干远场窄带信号，利用阵元接收数据生成和差伪接收数据矩阵，在对其进行矢量化去重处理后按照虚拟阵元位置大小顺序进行排列，得到和差连续虚拟阵元接收数据矩阵Z；最后对Z进行Toeplitz重构并利用ESPRIT算法得到DOA估计值；本发明提出的改进和差非均匀阵列，能够同时在阵元个数不变时扩充虚拟阵列孔径，提高估计性能。

Description

改进和差非均匀阵列的高精度波达方向估计方法

技术领域

本发明属于阵列信号处理领域，涉及一种改进和差非均匀阵列的高精度波达方向估计方法。

背景技术

到达角估计是阵列信号处理领域中一个重要的研究方向，被广泛的应用于雷达、通信、勘测等领域。在过去的几十年中，多重信号分类算法(multiple signalclassification，MUSIC)和旋转不变子空间算法(estimating signal via rotational，ESPRIT)被广泛地应用于该领域。在均匀线阵中，可精确估计的信号源数目必须小于阵元个数，并且当阵元间距为半波长时相邻阵元间产生耦合作用，极大地影响了到达角的估计精度。相关学者为了估计出比阵元个数更多的信号源，并且减少阵元间的耦合作用，提出了非均匀阵列的理论。常见的非均匀阵列有最小冗余阵列，互质阵列和嵌套阵列。最小冗余阵列的差分阵列没有孔洞，但是最小冗余阵列没有固定形式的表达式，因此，当给定传感器数量时，无法确定阵列的自由度大小。嵌套阵列和互质阵列有封闭的表达式，当阵元个数确定时就可以知道阵列孔径和阵列自由度，嵌套阵列没有孔洞，互质阵列却存在大量孔洞，极大地影响了互质阵列的性能，当两种阵列自由度相同时，互质阵列需要更多的阵元。

为了减少虚拟阵列孔洞，许多改进的互质阵列被陆续提出，例如扩展互质阵列，具有压缩因子的互质阵列，具有置换子阵列的互质阵列等。除了优化阵列结构外，许多新的增加稀疏阵列自由度的方法被提出，例如在无网格波达方向估计时，可以利用网格偏移量构造更大的差分阵列。此外还可以利用阵列插值算法进行孔洞填充来扩大阵列的自由度。近年来，相关学者提出了一种利用差值与和值进行联合波达方向估计的算法，该方法可以进一步增大阵列自由度，其中和差优化阵的位置集不仅包含阵元位置的差集，也包括阵元位置的和集，其虚拟阵元数远多于差分优化阵，因此其估计性能要优于差分优化阵，最大化和差优化阵虚拟阵元个数是增大阵列自由度提高波达方向(direction of arrival,DOA)估计精度的一个重要途径。本发明为了更好地利用稀疏阵列的和集和差集信息，提出了一种改进的和差非均匀阵列，该阵列与和差优化互质阵相比有更多的虚拟阵元，更大的自由度和阵列孔径，也有更高的估计精度。

发明内容

本发明的目的是提供一种改进和差非均匀阵列的高精度DOA估计方法。

一种改进和差非均匀阵列的高精度DOA估计方法步骤如下：

步骤一、用两个均匀子阵构成改进的一维非均匀阵列，子阵1由M个阵元构成，阵元间距为d，子阵2由N-1个阵元构成，阵元间距为(2M-1)d，将第2个子阵与第1个子阵排成一列并使第1个子阵的最后一个阵元和第2个子阵的第1个阵元重合，对合成阵列的阵元位置分别做求和、求差运算得虚拟阵元位置坐标的差集S_diff与和集S_sum，取集合S_diff与S_sum的并集即为虚拟阵元总的位置坐标集合S_total；

用两个均匀子阵构成改进的一维非均匀阵列，子阵1的阵元个数为M，阵元间距为d，阵元位置为{0,d,2d,…,(M-1)d}，子阵2的阵元个数为N-1,阵元间距为(2M-1)d，阵元位置为{(M-1)d,(3M-2)d,…,((2N-3)M-N+1)d}，d＝λ/2，λ为入射信号的波长，将第2个子阵与第1个子阵排成一列并使第1个子阵的最后一个阵元和第2个子阵的第1个阵元重合，那么改进非均匀阵列的阵元总数为P＝M+N-2(N≥3)，阵元位置的集合为D＝{d₁,d₂,d₃,…,d_P}＝{0,d,2d,…,(M-1)d,(3M-2)d,…,((2N-3)M-N+1)d}，分别对阵元位置做求和、求差运算，其中阵元位置的差集为

阵元位置的和集为

总的连续虚拟阵元位置坐标为

即S_total＝{-(2NM-2M-N)d，-(2NM-2M-N-1)d，…，0，…(2NM-2M-N-1)d，(2NM-2M-N)d}；虚拟阵元总数为2(2NM-2M-N)+1，当阵元总数为偶数且M和N取值满足M＝N时虚拟阵元总数最多，当阵元总数为奇数且M和N取值满足M＝N-1时虚拟阵元总数最多，以10阵元改进非均匀阵列为例，此时M+N-2＝10，M＝N＝6，阵元位置坐标为D＝{0,d,2d,3d,4d,5d,16d,27d,38d,49d}，

虚拟阵列位置坐标为S_total＝{-54d，-53d,-52d,…,-d,0,d,…,…,52d,53d,54d}，虚拟阵元总数为109，∪表示并集；

步骤二、将步骤一中的改进和差非均匀阵列作为接收阵列，接收K个非相干远场窄带信号，将第p个阵元接收信号x_p(t)与第1个阵元接收信号x₁(t)的相关函数

等价为d_p处虚拟阵元的接收数据，此时整个虚拟阵列接收数据为r_x(τ)，在多快拍下，r_x(τ)多次采样数据矩阵为Y₁，r_x ^*(-τ)多次采样数据矩阵为Y₂，将两个矩阵合成得到伪数据矩阵Y；

(2a)将改进和差非均匀阵列作为接收阵列得到接收数据矩阵X(t)；

K个非相干远场窄带信号入射到阵元个数为P的改进非均匀阵列上，信号入射角度分别为θ₁,…,θ_k,…,θ_K，在t时刻阵列接收信号为X(t)＝A(θ)s(t)+n(t)，其中，X(t)＝[x₁(t),…,x_p(t),…,x_P(t)]^T，s(t)＝[s₁(t),…,s_k(t),…,s_K(t)]^T表示入射信号矢量，

a_k表示第k个信号的幅度值，ω_k表示第k个信号的相位，n(t)＝[n₁(t),…,n_p(t),…,n_P(t)]^T表示与信号源不相关且相互独立的高斯白噪声，其均值为0，方差为

A(θ)＝[a(θ₁),…,a(θ_k),…,a(θ_K)]表示阵列流型矩阵，

为第k个信号的导向矢量，[·]^T为矩阵的转置；

(2b)利用接收信号生成伪数据矩阵Y；

在t时刻，将第1个阵元和第p个阵元的接收信号分别用x₁(t)，x_p(t)表示，取L次快拍，x_p(t)延时τ得到x_p(t+τ)，计算相关函数

[·]^*为共轭，将d₁位置的阵元设置为参考信源，即d₁＝0，

与入射信号的形式相同，可以将其等价为源信号，那么

可以等价为d_p处阵元的接收信号，即r_x(τ)＝Ar_s(τ)，其中

则r_x ^*(-τ)＝A^*r_s ^*(-τ)＝A^*r_s(τ)，令r(τ)＝[r_x(τ),r_x ^*(-τ)]^T＝[A,A^*]^Tr_s(τ)，设置伪采样周期为T_s，快拍数为L，此时得到伪数据矩阵Y，

设

其中Y₁＝A[r_s(T_s),…,r_s(lT_s),…,r_s(LT_s)]为对应阵元位置为D＝{d₁,…,d_p,…,d_P}时的伪接收数据矩阵，Y₂＝A^*[r_s(T_s),…,r_s(lT_s),…,r_s(LT_s)]为对应阵元位置是-D＝{-d₁,…,-d_p,…,-d_P}时的接收数据矩阵；

步骤三、分别求解矩阵Y₁的自协方差数据矩阵R_r1，以及Y₁和Y₂的互协方差矩阵R_r2，R_r3，分别对R_r1，R_r2，R_r3进行矢量化处理得到矢量化后的差伪数据矩阵Z₁及矢量化后的和伪数据矩阵Z₂，Z₃；

首先求解Y₁的自协方差数据矩阵R_r1，

即差伪数据协方差矩阵，然后求解Y₁，Y₂的互协方差矩阵R_r2，R_r3，

即和伪数据协方差矩阵，其中

[·]^H为复共轭转置，分别对R_r1，R_r2，R_r3对进行矢量化处理，得到矢量化后的和差伪数据矩阵Z₁，Z₂，Z₃，Z₁＝vec(R_r1)＝C₁Γ＝(A^*⊙A)Γ，C₁＝A^*⊙A中的任意元素可以表示为

Z₂＝vec(R_r2)＝C₂Γ＝(A⊙A)Γ，C₂＝A⊙A中的任意元素可以表示为

Z₃＝vec(R_r3)＝C₃Γ＝(A^*⊙A^*)Γ，C₃＝A^*⊙A^*中的任意元素可以表示为

其中

vec(·)表示对矩阵进行矢量化操作,即将矩阵元素按列排成一个列矢量，⊙表示求解Khatri-Rao积；

步骤四、将矢量化后的和差伪数据矩阵Z₁，Z₂，Z₃中的重复元素仅保留一个，并且将其按照连续虚拟阵元从小到大的顺序进行排列，得到和差连续虚拟阵元接收数据Z，对Z进行Toeplitz矩阵重构得到重构矩阵T；

根据阵元位置差集S_diff中的虚拟阵元位置找到Z₁中对应位置的一个数据，并将数据按照虚拟阵元位置从小到大排列，根据阵元位置和集子集

中的虚拟阵元位置找到Z₂中对应位置的一个数据，并将数据按照虚拟阵元位置从小到大排列，根据阵元位置和集子集

中的虚拟阵元位置找到Z₃中对应位置的一个数据，并将数据按照虚拟阵元位置从小到大排列，最后将所有找到的数据按照总的虚拟阵元位置集合S_total中对应位置元素的大小顺序进行排列，取连续部分得到和差连续虚拟阵元接收数据Z，设总的连续虚拟阵元个数为

和差连续虚拟阵元接收数据Z可以表示为

由Z构造Toeplitz矩阵T：

步骤五、对Toeplitz矩阵T进行特征分解，利用ESPRIT算法得到DOA的估计

对Toeplitz矩阵T进行特征分解，得到对应的信号子空间U_s，取U_s的前

行元素将其看做子阵1的信号子空间U_s1，取U_s的后

行元素将其看做子阵2的信号子空间U_s2，Ψ_z＝(U_S1)⁺U_S2，对Ψ_z进行特征分解，得到特征值

目标的到达角估计值可以通过

得到，即

angle(·)表示取相位。

前述步骤中，K表示信号源数目，k＝1,2,...,K表示信号源的标号，p＝1,2,...,P表示阵元位置的标号，l＝1,2,…,L表示快拍数，

表示虚拟阵元位置标号。

1.本发明提出了一种新型非均匀和差阵列结构，进一步扩大了阵列孔径，利用此阵列进行DOA估计，达到了预期的效果。

2.本发明不仅利用差协同阵列得到虚拟阵元接收数据，还利用和协同阵列得到虚拟阵元接收数据，极大程度的增加了虚拟阵元个数，具有很高的估计精度。

3.本发明的模型使用较少阵元数估计多于阵元个数的信号数，降低了成本，具有实际应用价值。

附图说明

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

图1为本发明的流程图；

图2是本发明非均匀阵列的结构图；

图3为本发明阵列与均匀阵列谱峰图；

图4为本发明方法空域目标DOA估计结果图；

图5为本发明方法与其他方法DOA随信噪比变化的均方根误差图；

具体实施方式

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

本发明的目的是提供一种改进和差非均匀阵列的DOA估计方法。

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

步骤一、用两个均匀子阵构成改进的一维非均匀阵列，子阵1由M个阵元构成，阵元间距为d，子阵2由N-1个阵元构成，阵元间距为(2M-1)d，将第2个子阵与第1个子阵排成一列并使第1个子阵的最后一个阵元和第2个子阵的第1个阵元重合，对合成阵列的阵元位置分别做求和、求差运算得虚拟阵元位置坐标的差集S_diff与和集S_sum，取集合S_diff与S_sum的并集即为虚拟阵元总的位置坐标集合S_total；

用两个均匀子阵构成改进的一维非均匀阵列，子阵1的阵元个数为M，阵元间距为d，阵元位置为{0,d,2d,…,(M-1)d}，子阵2的阵元个数为N-1,阵元间距为(2M-1)d，阵元位置为{(M-1)d,(3M-2)d,…,((2N-3)M-N+1)d}，d＝λ/2，λ为入射信号的波长，将第2个子阵与第1个子阵排成一列并使第1个子阵的最后一个阵元和第2个子阵的第1个阵元重合，那么改进非均匀阵列的阵元总数为P＝M+N-2(N≥3)，阵元位置的集合为D＝{d₁,d₂,d₃,…,d_P}＝{0,d,2d,…,(M-1)d,(3M-2)d,…,((2N-3)M-N+1)d}，分别对阵元位置做求和、求差运算，其中阵元位置的差集为

阵元位置的和集为

总的连续虚拟阵元位置坐标为

即S_total＝{-(2NM-2M-N)d，-(2NM-2M-N-1)d，…，0，…(2NM-2M-N-1)d，(2NM-2M-N)d}；虚拟阵元总数为2(2NM-2M-N)+1，当阵元总数为偶数且M和N取值满足M＝N时虚拟阵元总数最多，当阵元总数为奇数且M和N取值满足M＝N-1时虚拟阵元总数最多，以10阵元改进非均匀阵列为例，此时M+N-2＝10，M＝N＝6，阵元位置坐标为D＝{0,d,2d,3d,4d,5d,16d,27d,38d,49d}，

虚拟阵列位置坐标为S_total＝{-54d，-53d,-52d,…,-d,0,d,…,…,52d,53d,54d}，虚拟阵元总数为109，∪表示并集；

步骤二、将步骤一中的改进和差非均匀阵列作为接收阵列，接收K个非相干远场窄带信号，将第p个阵元接收信号x_p(t)与第1个阵元接收信号x₁(t)的相关函数

等价为d_p处虚拟阵元的接收数据，此时整个虚拟阵列接收数据为r_x(τ)，在多快拍下，r_x(τ)多次采样数据矩阵为Y₁，r_x ^*(-τ)多次采样数据矩阵为Y₂，将两个矩阵合成得到伪数据矩阵Y；

(2a)将改进和差非均匀阵列作为接收阵列得到接收数据矩阵X(t)；

K个非相干远场窄带信号入射到阵元个数为P的改进非均匀阵列上，信号入射角度分别为θ₁,…,θ_k,…,θ_K，在t时刻阵列接收信号为X(t)＝A(θ)s(t)+n(t)，其中，X(t)＝[x₁(t),…,x_p(t),…,x_P(t)]^T，s(t)＝[s₁(t),…,s_k(t),…,s_K(t)]^T表示入射信号矢量，

a_k表示第k个信号的幅度值，ω_k表示第k个信号的相位，n(t)＝[n₁(t),…,n_p(t),…,n_P(t)]^T表示与信号源不相关且相互独立的高斯白噪声，其均值为0，方差为

A(θ)＝[a(θ₁),…,a(θ_k),…,a(θ_K)]表示阵列流型矩阵，

为第k个信号的导向矢量，[·]^T为矩阵的转置；

(2b)利用接收信号生成伪数据矩阵Y；

在t时刻，将第1个阵元和第p个阵元的接收信号分别用x₁(t)，x_p(t)表示，取L次快拍，x_p(t)延时τ得到x_p(t+τ)，计算相关函数

[·]^*为共轭，将d₁位置的阵元设置为参考信源，即d₁＝0，

与入射信号的形式相同，可以将其等价为源信号，那么

可以等价为d_p处阵元的接收信号，即r_x(τ)＝Ar_s(τ)，其中

则r_x ^*(-τ)＝A^*r_s ^*(-τ)＝A^*r_s(τ)，令r(τ)＝[r_x(τ),r_x ^*(-τ)]^T＝[A,A^*]^Tr_s(τ)，设置伪采样周期为T_s，快拍数为L，此时得到伪数据矩阵Y，

设

其中Y₁＝A[r_s(T_s),…,r_s(lT_s),…,r_s(LT_s)]为对应阵元位置为D＝{d₁,…,d_p,…,d_P}时的伪接收数据矩阵，Y₂＝A^*[r_s(T_s),…,r_s(lT_s),…,r_s(LT_s)]为对应阵元位置是-D＝{-d₁,…,-d_p,…,-d_P}时的接收数据矩阵；

步骤三、分别求解矩阵Y₁的自协方差数据矩阵R_r1，以及Y₁和Y₂的互协方差矩阵R_r2，R_r3，分别对R_r1，R_r2，R_r3进行矢量化处理得到矢量化后的差伪数据矩阵Z₁及矢量化后的和伪数据矩阵Z₂，Z₃；

首先求解Y₁的自协方差数据矩阵R_r1，

即差伪数据协方差矩阵，然后求解Y₁，Y₂的互协方差矩阵R_r2，R_r3，

即和伪数据协方差矩阵，其中

[·]^H为复共轭转置，分别对R_r1，R_r2，R_r3对进行矢量化处理，得到矢量化后的和差伪数据矩阵Z₁，Z₂，Z₃，Z₁＝vec(R_r1)＝C₁Γ＝(A^*⊙A)Γ，C₁＝A^*⊙A中的任意元素可以表示为

Z₂＝vec(R_r2)＝C₂Γ＝(A⊙A)Γ，C₂＝A⊙A中的任意元素可以表示为

Z₃＝vec(R_r3)＝C₃Γ＝(A^*⊙A^*)Γ，C₃＝A^*⊙A^*中的任意元素可以表示为

其中

vec(·)表示对矩阵进行矢量化操作,即将矩阵元素按列排成一个列矢量，⊙表示求解Khatri-Rao积；

步骤四、将矢量化后的和差伪数据矩阵Z₁，Z₂，Z₃中的重复元素仅保留一个，并且将其按照连续虚拟阵元从小到大的顺序进行排列，得到和差连续虚拟阵元接收数据Z，对Z进行Toeplitz矩阵重构得到重构矩阵T；

根据阵元位置差集S_diff中的虚拟阵元位置找到Z₁中对应位置的一个数据，并将数据按照虚拟阵元位置从小到大排列，根据阵元位置和集子集

中的虚拟阵元位置找到Z₂中对应位置的一个数据，并将数据按照虚拟阵元位置从小到大排列，根据阵元位置和集子集

中的虚拟阵元位置找到Z₃中对应位置的一个数据，并将数据按照虚拟阵元位置从小到大排列，最后将所有找到的数据按照总的虚拟阵元位置集合S_total中对应位置元素的大小顺序进行排列，取连续部分得到和差连续虚拟阵元接收数据Z，设总的连续虚拟阵元个数为

和差连续虚拟阵元接收数据Z可以表示为

由Z构造Toeplitz矩阵T：

步骤五、对Toeplitz矩阵T进行特征分解，利用ESPRIT算法得到DOA的估计

对Toeplitz矩阵T进行特征分解，得到对应的信号子空间U_s，取U_s的前

行元素将其看做子阵1的信号子空间U_s1，取U_s的后

行元素将其看做子阵2的信号子空间U_s2，Ψ_z＝(U_s2)⁺U_s2，对Ψ_z进行特征分解，得到特征值

目标的到达角估计值可以通过

得到，即

angle(·)表示取相位。

前述步骤中，K表示信号源数目，k＝1,2,...,K表示信号源的标号，p＝1,2,...,P表示阵元位置的标号，l＝1,2,…,L表示快拍数，

表示虚拟阵元位置标号。

本发明提出的新型非均匀和差阵列结构，进一步扩大了阵列孔径，在进行数据处理的过程中，不仅利用差协同阵列得到虚拟阵元接收数据，还利用和协同阵列得到虚拟阵元接收数据，同时利用时间和空间信息，极大程度地增加了虚拟阵元个数，具有很高的估计精度，此外使用本发明阵列可以利用较少阵元估计多于阵元个数的信号数，降低了成本，具有实际应用价值。

仿真1：采用图1所示的改进和差非均匀阵列和均匀线阵作为接收阵列，改进和差非均匀阵列阵元总数为10，此时M+N-2＝10，M＝N＝6，阵元位置坐标为D＝{0,d,2d,3d,4d,5d,16d,27d,38d,49d}，均匀线阵阵元总数为50，设有17个非相干信号入射到阵列中，信号入射角度分别为[-80°,-70°,-60°,-50°,-40°,-30°,-20°,-10°,0°,10°,20°,30°,40°,50°,60°,70°,80°]，信噪比为SNR＝10dB，快拍数为L＝1000，阵列谱峰图如图3所示。

从图3可以看出在采用10阵元改进和差非均匀阵列和50阵元均匀线阵对17个入射信源的到达角方向进行估计时，改进和差非均匀阵列的谱峰图非常接近均匀线阵的谱峰图，能够精准地估计出入射信号的到达角方向，证明了改进非均匀阵列的有效性。

仿真2：采用图2所示的改进和差非均匀阵列作为接收阵列，阵元总数为10，此时M+N-2＝10，M＝N＝6，阵元位置坐标为D＝{0,d,2d,3d,4d,5d,16d,27d,38d,49d}，对空间K＝17个目标源进行方向识别，其到达阵列的角度分别为[-80°,-70°,-60°,-50°,-40°,-30°,-20°,-10°,0°,10°,20°,30°,40°,50°,60°,70°,80°]，在信噪比SNR＝8dB，快拍数为L＝200的环境下进行30次仿真实验，仿真结果如图4所示。

从图4中我们可以看到本发明所提阵列的估计结果与真实角度非常接近，角度偏离小，具有很好的估计精度。

仿真3：采用图2所示的改进和差非均匀阵列和传统和差互质阵列作为接收阵列，改进和差非均匀阵列和传统和差互质阵列均有6个阵元，改进和差非均匀阵列阵元位置为{0,d,2d,3d,4d,5d,16d,27d,38d,49d}，传统和差互质阵列阵元位置为{0,5d,6d,10d,12d,15d,18d,20d,24d,25d}，利用本发明和差优化算法对空间K＝2个目标源进行方向识别，其到达阵列的角度分别为(58°,60°)，快拍数为L＝120，信噪比从4dB到20dB的环境下，进行30次实验统计均方根误差，仿真结果如图5所示。

从图5可以清楚看到在阵元个数相同且和差优化处理方法相同时，本发明所提阵列进行DOA估计的均方根误差小于传统和差互质阵列，且本发明所提阵列的虚拟阵元个数明显多于互质阵列，在低信噪比下本发明方法同样有较高的估计精度。

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

Claims (5)

1.改进和差非均匀阵列的高精度波达方向估计方法，其特征在于，包括以下步骤：K个非相干远场窄带信号入射到阵元个数为P的改进非均匀阵列上，所述阵列是用两个均匀子阵构成改进的一维非均匀阵列，

步骤一、，子阵1由M个阵元构成，阵元间距为d，子阵2由N-1个阵元构成，阵元间距为(2M-1)d，将第2个子阵与第1个子阵排成一列并使第1个子阵的最后一个阵元和第2个子阵的第1个阵元重合，对合成阵列的阵元位置分别做求和、求差运算得虚拟阵元位置坐标的差集S_diff与和集S_sum，取集合S_diff与S_sum的并集即为虚拟阵元总的位置坐标集合S_total；

用两个均匀子阵构成改进的一维非均匀阵列，子阵1的阵元个数为M，阵元间距为d，阵元位置为{0,d,2d,…,(M-1)d}，子阵2的阵元个数为N-1,阵元间距为(2M-1)d，阵元位置为{(M-1)d,(3M-2)d,…,((2N-3)M-N+1)d}，d＝λ/2，λ为入射信号的波长，将第2个子阵与第1个子阵排成一列并使第1个子阵的最后一个阵元和第2个子阵的第1个阵元重合，那么改进非均匀阵列的阵元总数为P＝M+N-2(N≥3)，阵元位置的集合为D＝{d₁,d₂,d₃,…,d_P}＝{0,d,2d,…,(M-1)d,(3M-2)d,…,((2N-3)M-N+1)d}，分别对阵元位置做求和、求差运算，其中阵元位置的差集为

阵元位置的和集为

总的连续虚拟阵元位置坐标为

即S_total＝{-(2NM-2M-N)d，-(2NM-2M-N-1)d，…，0，…(2NM-2M-N-1)d，(2NM-2M-N)d}；虚拟阵元总数为2(2NM-2M-N)+1，当阵元总数为偶数且M和N取值满足M＝N时虚拟阵元总数最多，当阵元总数为奇数且M和N取值满足M＝N-1时虚拟阵元总数最多，以10阵元改进非均匀阵列为例，此时M+N-2＝10，M＝N＝6，阵元位置坐标为D＝{0,d,2d,3d,4d,5d,16d,27d,38d,49d}，

虚拟阵列位置坐标为S_total＝{-54d，-53d,-52d,…,-d,0,d,…,…,52d,53d,54d}，虚拟阵元总数为109，∪表示并集；

步骤二、将步骤一中的改进和差非均匀阵列作为接收阵列，接收K个非相干远场窄带信号，将第p个阵元接收信号x_p(t)与第1个阵元接收信号x₁(t)的相关函数

等价为d_p处虚拟阵元的接收数据，此时整个虚拟阵列接收数据为r_x(τ)，在多快拍下，r_x(τ)多次采样数据矩阵为Y₁，r_x ^*(-τ)多次采样数据矩阵为Y₂，将两个矩阵合成得到伪数据矩阵Y；

步骤三、分别求解矩阵Y₁的自协方差数据矩阵R_r1，以及Y₁和Y₂的互协方差矩阵R_r2，R_r3，分别对R_r1，R_r2，R_r3进行矢量化处理得到矢量化后的差伪数据矩阵Z₁及矢量化后的和伪数据矩阵Z₂，Z₃；

步骤四、将矢量化后的和差伪数据矩阵Z₁，Z₂，Z₃中的重复元素仅保留一个，并且将其按照连续虚拟阵元从小到大的顺序进行排列，得到和差连续虚拟阵元接收数据Z，对Z进行Toeplitz矩阵重构得到重构矩阵T；

步骤五、对Toeplitz矩阵T进行特征分解，利用ESPRIT算法得到DOA的估计

2.根据权利要求1所述的一种改进和差非均匀阵列的高精度DOA估计方法，将步骤一中的改进和差非均匀阵列作为接收阵列，接收K个非相干远场窄带信号，将第p个阵元接收信号x_p(t)与第1个阵元接收信号x₁(t)的相关函数

等价为d_p处虚拟阵元的接收数据，此时整个虚拟阵列接收数据为r_x(τ)，在多快拍下，r_x(τ)多次采样数据矩阵为Y₁，r_x ^*(-τ)多次采样数据矩阵为Y₂，将两个矩阵合成得到伪数据矩阵Y；

(2.1)将改进和差非均匀阵列作为接收阵列得到接收数据矩阵X(t)；

K个非相干远场窄带信号入射到阵元个数为P的改进非均匀阵列上，信号入射角度分别为θ₁,…,θ_k,…,θ_K，在t时刻阵列接收信号为X(t)＝A(θ)s(t)+n(t)，其中，X(t)＝[x₁(t),…,x_p(t),…,x_P(t)]^T，s(t)＝[s₁(t),…,s_k(t),…,s_K(t)]^T表示入射信号矢量，

a_k表示第k个信号的幅度值，ω_k表示第k个信号的相位，n(t)＝[n₁(t),…,n_p(t),…,n_P(t)]^T表示与信号源不相关且相互独立的高斯白噪声，其均值为0，方差为

A(θ)＝[a(θ₁),…,a(θ_k),…,a(θ_K)]表示阵列流型矩阵，

为第k个信号的导向矢量，[·]^T为矩阵的转置；

(2.2)利用接收信号生成伪数据矩阵Y；

在t时刻，将第1个阵元和第p个阵元的接收信号分别用x₁(t)，x_p(t)表示，取L次快拍，x_p(t)延时τ得到x_p(t+τ)，计算相关函数

[·]^*为共轭，将d₁位置的阵元设置为参考信源，即d₁＝0，

与入射信号的形式相同，可以将其等价为源信号，那么

可以等价为d_p处阵元的接收信号，即r_x(τ)＝Ar_s(τ)，其中

则r_x ^*(-τ)＝A^*r_s ^*(-τ)＝A^*r_s(τ)，令r(τ)＝[r_x(τ),r_x ^*(-τ)]^T＝[A,A^*]^Tr_s(τ)，设置伪采样周期为T_s，快拍数为L，此时得到伪数据矩阵Y，

设

其中Y₁＝A[r_s(T_s),…,r_s(lT_s),…,r_s(LT_s)]为对应阵元位置为D＝{d₁,…,d_p,…,d_P}时的伪接收数据矩阵，Y₂＝A^*[r_s(T_s),…,r_s(lT_s),…,r_s(LT_s)]为对应阵元位置是-D＝{-d₁,…,-d_p,…,-d_P}时的接收数据矩阵。

3.根据权利要求1所述的一种改进和差非均匀阵列的高精度DOA估计方法，分别求解矩阵Y₁的自协方差数据矩阵R_r1，以及Y₁和Y₂的互协方差矩阵R_r2，R_r3，分别对R_r1，R_r2，R_r3进行矢量化处理得到矢量化后的差伪数据矩阵Z₁及矢量化后的和伪数据矩阵Z₂，Z₃；

首先求解Y₁的自协方差数据矩阵R_r1，

即差伪数据协方差矩阵，然后求解Y₁，Y₂的互协方差矩阵R_r2，R_r3，

即和伪数据协方差矩阵，其中

[·]^H为复共轭转置，分别对R_r1，R_r2，R_r3对进行矢量化处理，得到矢量化后的和差伪数据矩阵Z₁，Z₂，Z₃，Z₁＝vec(R_r1)＝C₁Γ＝(A^*⊙A)Γ，C₁＝A^*⊙A中的任意元素可以表示为

Z₂＝vec(R_r2)＝C₂Γ＝(A⊙A)Γ，C₂＝A⊙A中的任意元素可以表示为

Z₃＝vec(R_r3)＝C₃Γ＝(A^*⊙A^*)Γ，C₃＝A^*⊙A^*中的任意元素可以表示为

其中

vec(·)表示对矩阵进行矢量化操作,即将矩阵元素按列排成一个列矢量，⊙表示求解Khatri-Rao积。

4.根据权利要求1所述的一种改进和差非均匀阵列的高精度DOA估计方法，将矢量化后的和差伪数据矩阵Z₁，Z₂，Z₃中的重复元素仅保留一个，并且将其按照连续虚拟阵元从小到大的顺序进行排列，得到和差连续虚拟阵元接收数据Z，对Z进行Toeplitz矩阵重构得到重构矩阵T；

根据阵元位置差集S_diff中的虚拟阵元位置找到Z₁中对应位置的一个数据，并将数据按照虚拟阵元位置从小到大排列，根据阵元位置和集子集

中的虚拟阵元位置找到Z₂中对应位置的一个数据，并将数据按照虚拟阵元位置从小到大排列，根据阵元位置和集子集

中的虚拟阵元位置找到Z₃中对应位置的一个数据，并将数据按照虚拟阵元位置从小到大排列，最后将所有找到的数据按照总的虚拟阵元位置集合S_total中对应位置元素的大小顺序进行排列，取连续部分得到和差连续虚拟阵元接收数据Z，设总的连续虚拟阵元个数为

和差连续虚拟阵元接收数据Z可以表示为

由Z构造Toeplitz矩阵T：

5.根据权利要求1所述的一种改进和差非均匀阵列的高精度DOA估计方法，对Toeplitz矩阵T进行特征分解，利用ESPRIT算法得到DOA的估计

对Toeplitz矩阵T进行特征分解，得到对应的信号子空间U_s，取U_s的前

行元素将其看做子阵1的信号子空间U_s1，取U_s的后

行元素将其看做子阵2的信号子空间U_s2，Ψ_z＝(U_s2)⁺U_s2，对Ψ_z进行特征分解，得到特征值

目标的到达角估计值可以通过

得到，即

angle(·)表示取相位。

Images