CN105824002A - 基于嵌套式子阵阵列的波达方向估计方法 - Google Patents

Abstract

本发明公开了一种基于嵌套式子阵阵列的波达方向估计方法，主要解决现有技术中存在自由度、阵列孔径和阵列密度不高的问题。其实现过程是：1.给定总阵元数确定子阵数目和子阵内阵元数目；2.根据子阵内阵元数目选取均匀线阵或最小冗余阵列或嵌套式线阵结构；3.根据子阵数目选取均匀线阵或最小冗余阵列或嵌套式线阵结构；4.根据选取的子阵内阵元结构和子阵结构构造嵌套式子阵阵列；5.根据嵌套式子阵阵列得到接收数据X(t)；6.根据X(t)得到差分合成阵列接收数据z_c，进而得到秩恢复数据协方差矩阵R_SS；7.对R_SS特征值分解，得到波达方向估计角度。本发明在同等条件下具有阵列配置灵活，波达方向估计测角性能好的优点，可用于雷达目标信号检测或功率估值。

Description

基于嵌套式子阵阵列的波达方向估计方法

本发明属于信号处理技术领域，特别涉及一种非等距线阵配置方法，用于估计目标个数大于阵元数情况下的波达方向，提高阵列的自由度、测角精度和密度。

背景技术

波达方向DOA估计是雷达、声纳信号处理中的一个重要研究方向。对于一个具有S个阵元的均匀线性阵列，采用传统的波达方向估计方法，如MUSIC，ESPRIT等所能分辨的最大目标个数为S-1个。而欠定波达方向估计问题，即目标个数大于阵元数的DOA估计经常出现并引起广泛的研究兴趣。解决该问题的一种有效方法是利用一个等效的虚拟阵列来提高波达方向估计的自由度DOF，该虚拟阵列是通过对一个特殊设计的阵列接收信号的协方差矩阵向量化来构造的。最近提出的嵌套式阵列NA和互质阵列CA就是利用稀疏阵的差分合成阵列来提高自由度的。而通过实际研究表明,最小冗余阵列MRA通过有效的阵元结构配置，可以获得最大的自由度DOF。在给定阵元数S及虚拟阵列为均匀线性阵列的情况下，最小冗余阵列MRA可以获得最大的阵列孔径。然而，最小冗余阵列的阵元位置及能够获得的自由度都没有明确的闭式解。虽然已有文献，如美国H.L.VanTrees教授的Optimumarrayprocessing:partIVofdetection,estimation,andmodulation通过穷搜的方法获得了阵元数小于等于17的最小冗余阵列，但却没有简单的方法来预测更大的最小冗余阵列MRA。

相对于最小冗余阵列MRA复杂的阵列结构获取过程，通过设计嵌套式阵列NA和互质阵列CA以获取更高的自由度DOF则很容易。

嵌套式阵列NA结构的阵元位置具有闭式解，并且通过利用接收数据的二阶统计信息，使用S个阵元就可以获得O(S²)的自由度。嵌套式阵列是由两个或多个具有不同阵元间隔的均匀线性子阵组成，它的缺点是第一个子阵的阵元间距通常较小，会引起阵元间的互耦问题。为此，P.P.Vaidyanathan等人后来提出的互质阵列CA可以通过增大阵元间距来减小阵元间的互耦问题。互质阵列是由两个阵元数目分别是M、N的均匀线性子阵构成，其中M、N为互质的整数，它能够获得多于阵元数目的自由度，但是由它构造的虚拟阵列并不是一个完全填充的均匀线性阵列，即在某些位置上没有虚拟阵元。通过利用虚拟阵列而不是原始阵列来进行DOA估计，这样由互质阵列获得的自由度就可以用来实现目标个数大于阵元数的DOA估计。但是在通过阵列接收信号的协方差矩阵构造虚拟阵列时，虚拟阵列中等效信源被实际信源的功率所替代，因此这些等效信源就表现为完全相干的信号。P.P.Vaidyanathan等提出采用一种基于空间平滑的MUSIC算法来解相干。然而由于基于空间平滑的估计算法一般只适用于均匀线阵，DOA估计时就只能利用互质阵列构造的虚拟阵列的一部分，这就会带来自由度和阵列孔径的降低。

然而，即便阵列结构拥有较高的自由度和阵列孔径，还要考虑到阵列密度问题，较低的阵列密度会引起发射功率受限，这在实际工程中是不希望出现的。

综上，现有的几种阵列虽然均能获得多于阵元数目的自由度，但都存在一定的局限，因此需要设计新的综合性能更优的阵列。

发明内容

本发明的目的在于针对上述已有阵列的不足，提出一种嵌套式子阵阵列的配置方法以减少阵元间的互耦，形成完全填充的均匀线性虚拟阵列，并获得较高的阵列密度、自由度和阵列孔径，提高波达方向DOA的估计精度。

本发明的技术思路是:根据总阵元数选择子阵数目及子阵内阵元数，利用现有文献结果得到相应的子阵结构和子阵内阵元结构；通过计算子阵间的最优单元间距，构造嵌套式子阵阵列；计算嵌套式子阵阵列的虚拟阵列以获得差分合成阵列，计算嵌套式子阵阵列的自由度；对差分合成阵列进行进行波达方向DOA估计。其实现方案包括如下：

1)构造嵌套式子阵阵列：

1a)给定总阵元数S，对S进行因式分解，得到子阵内阵元数目M及子阵数目N；

1b)根据子阵内阵元数目M，设计子阵内阵元结构，同时计算子阵内阵元的位置矢量u_M：

u_M＝[m₁,m₂,…,m_i,…,m_M]·d，

其中，m_i表示子阵内第i个阵元位置系数，1≤i≤M，且m₁＝0，d为入射信号的半波长；

1c)根据子阵数目N，设计子阵结构，同时计算N个子阵的阵位置矢量u_N：

u_N＝[n₁,n₂,…,n_j,…,n_N]·D，

其中，n_j表示第j个子阵位置系数，1≤j≤N，且n₁＝0，D＝L·d为子阵间的最优间隔单元，L表示子阵间最优间隔单元系数；

1d)根据上述设计参数，计算嵌套式子阵阵列位置矢量v：

$v=u_N\⊕u_M=\{n_j\·D+m_i\·d\}=\[p_1,p_2,...,p_k,...,p_S\]\·d,$

其中p_k表示嵌套式子阵阵列第k个阵元的位置系数，k＝1,2,…,S，符号表示交叉求和；

2)根据嵌套式子阵阵列，得到嵌套式子阵阵列的接收数据X(t)；

3)根据阵列接收数据X(t)估计嵌套式子阵阵列的协方差矩阵R_XX，向量化该协方差矩阵得到向量化接收数据z_ct，去除z_ct中重复的元素，得到虚拟的差分合成阵列接收数据z_c：

4)将差分合成阵列接收数据z_c划分为l_c个子阵接收数据，其中z_cg是第g个子阵接收数据，g＝1,2,…,l_c；计算第g个子阵的协方差矩阵R_g，并对该R_g进行前后向平滑平均，得到秩恢复的数据协方差矩阵R_SS，其中l_c表示划分的子阵个数，f_v表示嵌套式子阵阵列的自由度；

5)根据秩恢复的数据协方差矩阵R_SS估计出空间谱函数S(θ)，计算S(θ)中Q个极大值S(θ_q)，该S(θ_q)对应的角度θ_q即为波达方向，其中，θ_q表示第q个入射信号的入射角，q＝1,2,…,Q，Q表示入射信号的总个数。

本发明与现有阵列结构相比具有以下优点：

1)能获得孔径更大的阵列

最小冗余阵列MRA的阵元位置及可用自由度都没有闭式解，其设计往往需要复杂的计算机穷搜，而本发明阵列的阵元位置及可用自由度可利用已知的阵列的参数预测获得，实现简单，易于获得孔径更大的阵列。

2)测角性能好

现有互质阵列CA构造的虚拟差分合成阵列由于不是完全填充的均匀线阵，会降低DOA估计可用的自由度及测角精度；而本发明阵列构造的差分合成阵列为完全填充的均匀线阵，不仅在同样物理阵元数情况下可获得比CA更大的阵列孔径和更高的自由度，而且能估计更多的目标，具有更好的测角性能。

3)易于生产、更换和扩展

本发明的所有子阵结构相同，便于批量生产，对出现问题的子阵更换简单，可方便安装和扩展成更大的阵列。

4)成本低

本发明阵列是稀疏阵，在保持阵列面积不变的情况下，能够有效减少阵元数目，降低成本，由于各个子阵完全相同，批量生产更利于降低成本。

5)机动性高，角度分辨能力强

本发明的所有子阵都可当成一个子雷达系统，单个子雷达体积较小，机动性好，便于变换阵地；同时由于阵列的稀疏配置，可获得更好的角度分辨能力。

附图说明

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

图2是本发明第2实施例中的16阵元的阵列配置几何结构示意图；

图3是本发明第6实施例中的16阵元的阵列配置几何结构示意图；

图4是本发明第8实施例中的16阵元的阵列配置几何结构示意图；

图5是对用本发明九种实施例阵列配置的差分合成阵列自由度与总阵元数关系比较图；

图6是对用本发明九种实施例阵列配置的阵列孔径与总阵元数关系比较图；

图7是对用本发明九种实施例的阵列配置的阵列密度与总阵元数关系比较图；

图8是对用本发明五种实施例的阵列配置的波达方向DOA估计均方根误差与信噪比关系比较图；

图9是对用本发明五种实施例的阵列配置的波达方向DOA估计均方根误差与快拍数关系比较图。

具体实施方式

以下参照附图对本发明的实施和技术效果作进一步说明：

参照图1，本发明给出如下9种实施例：

实施例1，子阵结构和子阵内阵元结构均采用均匀线阵ULA结构的嵌套式子阵阵列波达方向估计。

步骤1，给定总阵元数S构造子阵结构及子阵内阵元结构。

现有技术中的非均匀阵列结构包括三：一种是最小冗余阵列MRA，其构造是根据总阵元数S得到最小冗余序列和最小冗余阵列MRA的阵列结构；第二种是嵌套式阵列NA，其构造是根据总阵元数S得到阵列的嵌套级数，然后根据嵌套级数确定每一级上的阵元数目，从而得到嵌套式阵列NA的阵列结构；第三种是互质阵列CA，其结构是根据总阵元数S得到两个子阵的阵元数目，从而得到互质阵列CA的阵列结构。

本发明的子阵结构及子阵内阵元结构均采用均匀线阵ULA结构，其设计步骤如下：

1a)给定总阵元数S，对S进行因式分解，得到子阵内阵元数M及子阵数目N：

S＝M·N，

其中M和N的取值最接近，且M≥2,N≥2；

1b)根据子阵内阵元数目M，将子阵内阵元结构设计为与现有的均匀线阵ULA结构相同，同时计算子阵内阵元的位置矢量u_1M：

u_1M＝[m₁₁,m₁₂,…,m_1i,…m_1M,]·d，

其中，m_1i表示子阵内阵元结构为均匀线阵ULA结构的第i个阵元位置系数，1≤i≤M，且m₁₁＝0，d为入射信号的半波长；

1c)根据子阵数目N，将子阵结构设计为与现有的均匀线阵ULA结构相同，同时计算N个子阵的阵位置矢量u_1N：

u_1N＝[n₁₁,n₁₂,…,n_1j,…n_1M,]·D，

其中，n_1j表示子阵结构为均匀线阵ULA结构的第j个子阵位置系数，1≤j≤N，且n₁₁＝0，D＝L·d为子阵间的最优间隔单元，L表示子阵间最优间隔单元系数；

1d)根据上述设计参数，计算嵌套式子阵阵列位置矢量v₁、嵌套式子阵阵列的阵列孔径a_1v和嵌套式子阵阵列的自由度f_1v：

$v_1=u_{1N}\⊕u_{1M}=\{n_{1j}\·D+m_{1i}\·d\}=\[p_{11},p_{12},...,p_{1k},...,p_{1S}\]\·d,$

a_1v＝(N-1)(2M-1)+M-1，

f_1v＝(2N-1)(2M-1)，

其中符号表示交叉求和；p_1k表示嵌套式子阵阵列第k个阵元的位置系数，k＝1,2,…,S；

1f)根据嵌套式子阵阵列位置矢量v₁，得到所有阵元距离参考阵元的间距，从而确定虚拟阵列为完全填充型均匀线阵的嵌套式子阵阵列。

步骤2，根据所设计的嵌套式子阵阵列位置矢量，得到嵌套式子阵阵列的接收数据X(t)。

现有技术中对阵列接收数据X(t)的表达形式有窄带信号接收数据，宽带信号接收数据，本发明采用窄带信号接收数据表达形式，其具体实现如下：

2a)根据嵌套式子阵阵列位置矢量v，计算第q个入射信号的导向矢量a(θ_q)：

$a\left({\mathrm{\θ}}_q\right)={\[e^{{\mathrm{j\πp}}_1{\mathrm{sin\θ}}_q},e^{{\mathrm{j\πp}}_2{\mathrm{sin\θ}}_q},...,e^{{\mathrm{j\πp}}_k{\mathrm{sin\θ}}_q},...,e^{{\mathrm{j\πp}}_S{\mathrm{sin\θ}}_q}\]}^T,$

其中，表示第q个入射信号在第k个阵元处的旋转因子，p_k表示嵌套式子阵阵列第k个阵元的位置系数，q＝1,2,…,Q，Q表示入射信号的总个数；

2b)根据入射信号的导向矢量a(θ_q)，构造导向矢量矩阵A：

A＝[a(θ₁),a(θ₂),…,a(θ_q),…,a(θ_Q)]，

2c)根据导向矢量矩阵A，入射信号S(t)，复高斯白噪声n(t)，得到嵌套式子阵阵列的接收数据X(t)：

X(t)＝AS(t)+n(t)，

其中，S(t)＝[s₁(t),s₂(t),…,s_q(t),…,s_Q(t)]^T，s_q(t)表示第q个入射信号，()^T表示向量的转置，n(t)表示均值为0，方差为的复高斯白噪声，并且与入射信号不相关，t表示采样时间t＝1,2,…,T，T表示快拍数。

步骤3，根据嵌套式子阵阵列的接收数据X(t)估计阵列协方差矩阵R_XX，得到虚拟的差分合成阵列接收数据z_c。

3a)根据嵌套式子阵阵列的接收数据X(t)估计阵列协方差矩R_XX：

$\begin{array}{c}{R}_{XX}=E\[X\left(t\right)X^H\left(t\right)\]={\mathrm{AR}}_{SS}{A}^H\\ =\underset{q=1}{\overset{Q}{\Σ}}{\mathrm{\σ}}_q^{2}a\left({\mathrm{\θ}}_q\right)a^H\left({\mathrm{\θ}}_q\right)+{\mathrm{\σ}}_n^2{I}_s\end{array},$

其中，其中E表示求数学期望；(·)^H表示共轭转置；I_s表示S×S的单位矩阵，表示第q个入射信号的功率；

3b)根据阵列协方差矩R_XX，计算向量化接收数据z_ct

$z_{ct}=vec\left(R_{XX}\right)=Bc+{\mathrm{\σ}}_n^2{1}_n,$

其中，vec表示矩阵的向量化操作；B表示虚拟的差分合成阵列的方位矢量矩阵， 表示Kronecker积，()^*表示取共轭，e_i表示长度为S,第i个元素为1，其他元素全为0的行向量，()^T表示转置，表示入射信号的功率矢量；

3c)根据向量化接收数据z_ct，得到虚拟的差分合成阵列接收数据z_c：

z_c＝Dis(z_ct)，

其中，Dis表示去除向量中重复元素的操作。

步骤4，根据虚拟的差分合成阵列接收数据z_c，计算秩恢复的数据协方差矩阵R_SS。

4a)根据虚拟的差分合成阵列接收数据z_c，得到子阵接收数据z_cg，即取差分合成阵列接收数据z_c中位于(g+n-1)处的l_c个数据构成子阵接收数据z_cg，n＝1,2,…,l_c,g＝1,2,…,l_c，

其中，l_c表示划分的子阵个数，f_v表示嵌套式子阵阵列的自由度；

4b)根据子阵接收数据z_cg，计算第g个子阵的协方差矩阵R_g：

$R_g=z_{cg}{z}_{cg}^H;$

4c)根据第g子阵的协方差矩阵R_g，计算前向平滑的协方差矩阵R^f：

$R^f=\frac{1}{l_c}\underset{g=1}{\overset{l_c}{\Σ}}{R}_g;$

4d)根据子阵接收数据z_cg计算第g个子阵的后向接收数据y_cg；

y_cg＝Trans(z_cg)，

其中，Trans表示左右翻转并取共轭；

4e)根据第g个子阵的后向接收数据y_cg；计算第g个子阵的后向协方差矩阵R_g ^b：

${R_g}^b=y_{cg}{y}_{cg}^H;$

4f)根据第g个子阵的后向协方差矩阵R_g ^b，计算后向平滑的协方差矩阵R^b：

$R^b=\frac{1}{l_c}\underset{g=1}{\overset{l}{\Σ}}{R}_g^b;$

4g)根据前向平滑的协方差矩阵R^f和后向平滑的协方差矩阵R^b，计算秩恢复的数据协方差矩阵R_SS：

$R_{SS}=\frac{R^f+R^b}{2}.$

以上步骤3-步骤4内容皆为现有常用方法。

步骤5，根据R_SS估计出空间谱函数S(θ)，计算S(θ)中Q个极大值S(θ_q)，S(θ_q)对应的角度θ_q即为波达方向。

现有技术中利用R_SS估计出空间谱函数S(θ)有最大似然算法ML，加权子空间拟合法WSF，波束成形法DBF，基于特征空间的多重信号分类法MUSIC，旋转因子不变法ESPRIT，本发明中采用基于特征空间的多重信号分类法MUSIC，其具体实现如下：

5a)对R_SS进行特征值分解：

$R_{SS}=\underset{t=1}{\overset{T}{\Σ}}{\mathrm{\λ}}_t{V}_t{V_t}^H,$

其中，λ_t表示矩阵R_SS的特征值，V_t表示矩阵R_SS对应于λ_t的特征矢量，t表示采样时间，t＝1,2,…，T，T表示快拍数，Q表示入射信号的总个数；

5b)根据V_t，计算空间谱函数S(θ)：

$S\left(\mathrm{\θ}\right)=\frac{1}{\left|\underset{t=Q+1}{\overset{T}{\Σ}}a\left(\mathrm{\θ}\right)V_t{V_t}^{H}a\left(\mathrm{\θ}\right)\right|},$

其中，a(θ)表示入射角度为θ的入射信号的导向矢量，

表示第q个入射信号在第k个阵元处的旋转因子，p_k表示嵌套式子阵阵列第k个阵元的位置系数；

5c)根据空间谱函数S(θ)，计算S(θ)中Q个极大值S(θ_q)，S(θ_q)对应的角度θ_q即为波达方向。

实施例2，子阵结构采用均匀线阵ULA结构、子阵内阵元结构采用最小冗余阵列MRA结构的嵌套式子阵阵列波达方向估计。

该实施例2与实施例1的区别在于步骤1构造子阵结构及子阵内阵元结构不同，其他步骤均与实施例1相同，其构造子阵结构及子阵内阵元结构步骤如下：

(1a)给定总阵元数S，对S进行因式分解，得到子阵内阵元数M及子阵数目N：

S＝M·N，

其中M和N的取值最接近，且M≥2,N≥2；

(1b)根据子阵内阵元数目M，将子阵内阵元结构设计为与现有的最小冗余阵列MRA结构相同，同时计算子阵内阵元的位置矢量u_2M：

u_2M＝[m₂₁,m₂₂,…,m_2i,…m_2M,]·d，

其中，m_2i表示子阵内阵元结构为最小冗余阵列MRA结构的第i个阵元位置系数，1≤i≤M，且m₂₁＝0，d为入射信号的半波长；

(1c)根据子阵数目N，将子阵结构设计为与现有的均匀线阵ULA结构相同，同时计算N个子阵的阵位置矢量u_2N：

u_2N＝[n₂₁,n₂₂,…,n_2j,…n_2M,]·D，

其中，n_2j表示子阵结构为均匀线阵ULA结构的第j个子阵位置系数，1≤j≤N，且n₂₁＝0，D＝L·d为子阵间的最优间隔单元，L表示子阵间最优间隔单元系数；

(1d)根据上述设计参数，计算嵌套式子阵阵列位置矢量v₂、嵌套式子阵阵列的阵列孔径a_2v和嵌套式子阵阵列的自由度f_2v：

$v_2=u_{2N}\⊕u_{2M}=\{n_{2j}\·D+m_{2i}\·d\}=\[p_{21},p_{22},...,p_{2k},...,p_{2S}\]\·d,$

a_2v＝(N-1)[M(M-1)+1-2f_M]+M(M-1)/2-f_M，

f_2v＝(2N-1)[M(M-1)+1-2f_M]，

其中符号表示交叉求和；p_2k表示嵌套式子阵阵列第k个阵元的位置系数，k＝1,2,…,S，f_M表示子阵内阵元结构为最小冗余阵列MRA结构的差分合成阵列的冗余度；

(1e)根据嵌套式子阵阵列位置矢量v₂，得到所有阵元距离参考阵元的间距，从而确定虚拟阵列为完全填充型均匀线阵的嵌套式子阵阵列，其结构如图2所示。

实施例3，子阵结构采用均匀线阵ULA结构、子阵内阵元结构采用嵌套式线阵NA结构的嵌套式子阵阵列波达方向估计。

该实施例3与实施例1的区别在于步骤1构造子阵结构及子阵内阵元结构不同，其他步骤均与实施例1相同，其构造子阵结构及子阵内阵元结构步骤如下：

1-a)给定总阵元数S，对S进行因式分解，得到子阵内阵元数M及子阵数目N：

S＝M·N，

其中M和N的取值最接近，且M≥2,N≥2；

1-b)根据子阵内阵元数目M，将子阵内阵元结构设计为与现有的嵌套式线阵NA结构相同，同时计算子阵内阵元的位置矢量u_3M：

u_3M＝[m₃₁,m₃₂,…,m_3i,…m_3M,]·d，

其中，m_3i表示子阵内阵元结构为嵌套式线阵NA结构的第i个阵元位置系数，1≤i≤M，且m₃₁＝0，d为入射信号的半波长；

1-c)根据子阵数目N，将子阵结构设计为与现有的均匀线阵ULA结构相同，同时计算N个子阵的阵位置矢量u_3N：

u_3N＝[n₃₁,n₃₂,…,n_3j,…n_3M,]·D，

其中，n_3j表示子阵结构为均匀线阵ULA结构的第j个子阵位置系数，1≤j≤N，且n₃₁＝0，D＝L·d为子阵间的最优间隔单元，L表示子阵间最优间隔单元系数；

1-d)根据上述设计参数，计算嵌套式子阵阵列位置矢量v₃、嵌套式子阵阵列的阵列孔径a_3v和嵌套式子阵阵列的自由度f_3v：

$v_3=u_{3N}\⊕u_{3M}=\{n_{3j}\·D+m_{3i}\·d\}=\[p_{31},p_{32},...,p_{3k},...,p_{3S}\]\·d,$

其中符号表示交叉求和；p_3k表示嵌套式子阵阵列第k个阵元的位置系数，k＝1,2,…,S；

1-e)根据嵌套式子阵阵列位置矢量v₃，得到所有阵元距离参考阵元的间距，从而确定虚拟阵列为完全填充型均匀线阵的嵌套式子阵阵列。

实施例4，子阵结构和子阵内阵元结构均采用最小冗余阵列MRA结构的嵌套式子阵阵列波达方向估计。

该实施例4与实施例1的区别在于步骤1构造子阵结构及子阵内阵元结构不同，其他步骤均与实施例1相同，其构造子阵结构及子阵内阵元结构步骤如下：

(1-a)给定总阵元数S，对S进行因式分解，得到子阵内阵元数M及子阵数目N：

S＝M·N，

其中M和N的取值最接近，且M≥2,N≥2；

(1-b)根据子阵内阵元数目M，将子阵内阵元结构设计为与现有的最小冗余阵列MRA结构相同，同时计算子阵内阵元的位置矢量u_4M：

u_4M＝[m₄₁,m₄₂,…,m_4i,…m_4M,]·d，

其中，m_4i表示子阵内阵元结构为最小冗余阵列MRA结构的第i个阵元位置系数，1≤i≤M，且m₄₁＝0，d为入射信号的半波长；

(1-c)根据子阵数目N，将子阵结构设计为与现有的最小冗余阵列MRA结构相同，同时计算N个子阵的阵位置矢量u_4N：

u_4N＝[n₄₁,n₄₂,…,n_4j,…n_4M,]·D，

其中，n_4j表示子阵结构为均匀线阵ULA结构的第j个子阵位置系数，1≤j≤N，且n₄₁＝0，D＝L·d为子阵间的最优间隔单元，L表示子阵间最优间隔单元系数；

(1-d)根据上述设计参数，计算嵌套式子阵阵列位置矢量v₄、嵌套式子阵阵列的阵列孔径a_4v和嵌套式子阵阵列的自由度f_4v：

$v_4=u_{4N}\⊕u_{4M}=\{n_{4j}\·D+m_{4i}\·d\}=\[p_{41},p_{42},...,p_{4k},...,p_{4S}\]\·d,$

a_4v＝[N(N-1)/2-f_N][M(M-1)+1-2f_M]+M(M-1)/2-f_M，

f_4v＝[N(N-1)+1-2f_N][M(M-1)/2-f_M]，

其中符号表示交叉求和；p_4k表示嵌套式子阵阵列第k个阵元的位置系数，k＝1,2,…,S，f_N表示子阵结构为最小冗余阵列MRA结构的差分合成阵列的冗余度；

(1-e)根据嵌套式子阵阵列位置矢量v₄，得到所有阵元距离参考阵元的间距，从而确定虚拟阵列为完全填充型均匀线阵的嵌套式子阵阵列。

实施例5，子阵结构采用最小冗余阵列MRA结构、子阵内阵元结构采用均匀线阵ULA结构的嵌套式子阵阵列波达方向估计。

该实施例5与实施例1的区别在于步骤1构造子阵结构及子阵内阵元结构不同，其他步骤均与实施例1相同，其构造子阵结构及子阵内阵元结构步骤如下：

1A)给定总阵元数S，对S进行因式分解，得到子阵内阵元数M及子阵数目N：

S＝M·N，

其中M和N的取值最接近，且M≥2,N≥2；

1B)根据子阵内阵元数目M，将子阵内阵元结构设计为与现有的均匀线阵ULA结构相同，同时计算子阵内阵元的位置矢量u_5M：

u_5M＝[m₅₁,m₅₂,…,m_5i,…m_5M,]·d，

其中，m_5i表示子阵内阵元结构为均匀线阵ULA结构的第i个阵元位置系数，1≤i≤M，且m₅₁＝0，d为入射信号的半波长；

1C)根据子阵数目N，将子阵结构设计为与现有的最小冗余阵列MRA结构相同，同时计算N个子阵的阵位置矢量u_5N：

u_5N＝[n₅₁,n₅₂,…,n_5j,…n_5M,]·D，

其中，n_5j表示子阵结构为最小冗余阵列MRA结构的第j个子阵位置系数，1≤j≤N，且n₅₁＝0，D＝L·d为子阵间的最优间隔单元，L表示子阵间最优间隔单元系数；

1D)根据上述设计参数，计算嵌套式子阵阵列位置矢量v₅、嵌套式子阵阵列的阵列孔径a_5v和嵌套式子阵阵列的自由度f_5v：

$v_5=u_{5N}\⊕u_{5M}=\{n_{5j}\·D+m_{5i}\·d\}=\[p_{51},p_{52},...,p_{5k},...,p_{5S}\]\·d,$

a_5v＝[N(N-1)/2-f_N](2M-1)+M-1，

f_5v＝[N(N-1)+1-2f_N](2M-1)，

其中符号表示交叉求和；p_5k表示嵌套式子阵阵列第k个阵元的位置系数，k＝1,2,…,S；

1E)根据嵌套式子阵阵列位置矢量v₅，得到所有阵元距离参考阵元的间距，从而确定虚拟阵列为完全填充型均匀线阵的嵌套式子阵阵列。

实施例6，子阵结构采用最小冗余阵列MRA结构、子阵内阵元结构采用嵌套式线阵NA结构的嵌套式子阵阵列波达方向估计。

该实施例6与实施例1的区别在于步骤1构造子阵结构及子阵内阵元结构不同，其他步骤均与实施例1相同，其构造子阵结构及子阵内阵元结构步骤如下：

(1A)给定总阵元数S，对S进行因式分解，得到子阵内阵元数M及子阵数目N：

S＝M·N，

其中M和N的取值最接近，且M≥2,N≥2；

(1B)根据子阵内阵元数目M，将子阵内阵元结构设计为与现有的嵌套式线阵NA结构相同，同时计算子阵内阵元的位置矢量u_6M：

u_6M＝[m₆₁,m₆₂,…,m_6i,…m_6M,]·d，

其中，m_6i表示子阵内阵元结构为嵌套式线阵NA结构的第i个阵元位置系数，1≤i≤M，且m₆₁＝0，d为入射信号的半波长；

(1C)根据子阵数目N，将子阵结构设计为与现有的最小冗余阵列MRA结构相同，同时计算N个子阵的阵位置矢量u_6N：

u_6N＝[n₆₁,n₆₂,…,n_6j,…n_6M,]·D，

其中，n_6j表示子阵结构为最小冗余阵列MRA结构的第j个子阵位置系数，1≤j≤N，且n₆₁＝0，D＝L·d为子阵间的最优间隔单元，L表示子阵间最优间隔单元系数；

(1D)根据上述设计参数，计算嵌套式子阵阵列位置矢量v₆、嵌套式子阵阵列的阵列孔径a_6v和嵌套式子阵阵列的自由度f_6v：

$v_6=u_{6N}\⊕u_{6M}=\{n_{6j}\·D+m_{6i}\·d\}=\[p_{61},p_{62},...,p_{6k},...,p_{6S}\]\·d,$

其中符号表示交叉求和；p_6k表示嵌套式子阵阵列第k个阵元的位置系数，k＝1,2,…,S；

(1E)根据嵌套式子阵阵列位置矢量v₆，得到所有阵元距离参考阵元的间距，从而确定虚拟阵列为完全填充型均匀线阵的嵌套式子阵阵列，其结构如图3所示。

实施例7，子阵结构和子阵内阵元结构均采用嵌套式线阵NA结构的嵌套式子阵阵列波达方向估计。

该实施例7与实施例1的区别在于步骤1构造子阵结构及子阵内阵元结构不同，其他步骤均与实施例1相同，其构造子阵结构及子阵内阵元结构步骤如下：

1-1)给定总阵元数S，对S进行因式分解，得到子阵内阵元数M及子阵数目N：

S＝M·N，

其中M和N的取值最接近，且M≥2,N≥2；

1-2)根据子阵内阵元数目M，将子阵内阵元结构设计为与现有的嵌套式线阵NA结构相同，同时计算子阵内阵元的位置矢量u_7M：

u_7M＝[m₇₁,m₇₂,…,m_7i,…m_7M,]·d，

其中，m_7i表示子阵内阵元结构为嵌套式线阵NA结构的第i个阵元位置系数，1≤i≤M，且m₇₁＝0，d为入射信号的半波长；

1-3)根据子阵数目N，将子阵结构设计为与现有的嵌套式线阵NA结构相同，同时计算N个子阵的阵位置矢量u_7N：

u_7N＝[n₇₁,n₇₂,…,n_7j,…n_7M,]·D，

其中，n_7j表示子阵结构为嵌套式线阵NA结构的第j个子阵位置系数，1≤j≤N，且n₇₁＝0，D＝L·d为子阵间的最优间隔单元，L表示子阵间最优间隔单元系数；

1-4)根据上述设计参数，计算嵌套式子阵阵列位置矢量v₇、嵌套式子阵阵列的阵列孔径a_7v和嵌套式子阵阵列的自由度f_7v：

$v_7=u_{7N}\⊕u_{7M}=\{n_{7j}\·D+m_{7i}\·d\}=\[p_{71},p_{72},...,p_{7k},...,p_{7S}\]\·d,$

其中符号表示交叉求和；p_7k表示嵌套式子阵阵列第k个阵元的位置系数，k＝1,2,…,S；

1-5)根据嵌套式子阵阵列位置矢量v₇，得到所有阵元距离参考阵元的间距，从而确定虚拟阵列为完全填充型均匀线阵的嵌套式子阵阵列。

实施例8，子阵结构采用嵌套式线阵NA结构、子阵内阵元结构采用均匀线阵ULA结构的嵌套式子阵阵列波达方向估计。

该实施例8与实施例1的区别在于步骤1构造子阵结构及子阵内阵元结构不同，其他步骤均与实施例1相同，其构造子阵结构及子阵内阵元结构步骤如下：

(1-1)给定总阵元数S，对S进行因式分解，得到子阵内阵元数M及子阵数目N：

S＝M·N，

其中M和N的取值最接近，且M≥2,N≥2；

(1-2)根据子阵内阵元数目M，将子阵内阵元结构设计为与现有的均匀线阵ULA结构相同，同时计算子阵内阵元的位置矢量u_8M：

u_8M＝[m₈₁,m₈₂,…,m_8i,…m_8M,]·d，

其中，m_8i表示子阵内阵元结构为均匀线阵ULA结构的第i个阵元位置系数，1≤i≤M，且m₈₁＝0，d为入射信号的半波长；

(1-3)根据子阵数目N，将子阵结构设计为与现有的嵌套式线阵NA结构相同，同时计算N个子阵的阵位置矢量u_8N：

u_8N＝[n₈₁,n₈₂,…,n_8j,…n_8M,]·D，

其中，n_8j表示子阵结构为嵌套式线阵NA结构的第j个子阵位置系数，1≤j≤N，且n₈₁＝0，D＝L·d为子阵间的最优间隔单元，L表示子阵间最优间隔单元系数；

(1-4)根据上述设计参数，计算嵌套式子阵阵列位置矢量v₈、嵌套式子阵阵列的阵列孔径a_8v和嵌套式子阵阵列的自由度f_8v：

$v_8=u_{8N}\⊕u_{8M}=\{n_{8j}\·D+m_{8i}\·d\}=\[p_{81},p_{82},...,p_{8k},...,p_{8S}\]\·d,$

其中符号表示交叉求和；p_8k表示嵌套式子阵阵列第k个阵元的位置系数，k＝1,2,…,S；

(1-5)根据嵌套式子阵阵列位置矢量v₈，得到所有阵元距离参考阵元的间距，从而确定虚拟阵列为完全填充型均匀线阵的嵌套式子阵阵列，其结构如图4所示。

实施例9，子阵结均采用嵌套式线阵NA结构、子阵内阵元结构采用最小冗余阵列MRA结构的嵌套式子阵阵列波达方向估计。

该实施例9与实施例1的区别在于步骤1构造子阵结构及子阵内阵元结构不同，其他步骤均与实施例1相同，其构造子阵结构及子阵内阵元结构步骤如下：

(1-A)给定总阵元数S，对S进行因式分解，得到子阵内阵元数M及子阵数目N：

S＝M·N，

其中M和N的取值最接近，且M≥2,N≥2；

(1-B)根据子阵内阵元数目M，将子阵内阵元结构设计为与现有的最小冗余阵列MRA结构相同，同时计算子阵内阵元的位置矢量u_9M：

u_9M＝[m₉₁,m₉₂,…,m_9i,…m_9M,]·d，

其中，m_9i表示子阵内阵元结构为最小冗余阵列MRA结构的第i个阵元位置系数，1≤i≤M，且m₉₁＝0，d为入射信号的半波长；

(1-C)根据子阵数目N，将子阵结构设计为与现有的嵌套式线阵NA结构相同，同时计算N个子阵的阵位置矢量u_9N：

u_9N＝[n₉₁,n₉₂,…,n_9j,…n_9M,]·D，

其中，n_9j表示子阵结构为嵌套式线阵NA结构的第j个子阵位置系数，1≤j≤N，且n₉₁＝0，D＝L·d为子阵间的最优间隔单元，L表示子阵间最优间隔单元系数；

(1-D)根据上述设计参数，计算嵌套式子阵阵列位置矢量v₉、嵌套式子阵阵列的阵列孔径a_9v和嵌套式子阵阵列的自由度f_9v：

$v_9=u_{9N}\⊕u_{9M}=\{n_{9j}\·D+m_{9i}\·d\}=\[p_{91},p_{92},...,p_{9k},...,p_{9S}\]\·d,$

其中符号表示交叉求和；p_9k表示嵌套式子阵阵列第k个阵元的位置系数，k＝1,2,…,S；

(1-E)根据嵌套式子阵阵列位置矢量v₉，得到所有阵元距离参考阵元的间距，从而确定虚拟阵列为完全填充型均匀线阵的嵌套式子阵阵列。

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

仿真1：对本发明的九种实施例的阵列配置的差分合成阵列的自由度、有效孔径和阵列密度与总阵元数的关系进行仿真。

1.1)仿真条件：总阵元数为9～36间的一些整数。

1.2)仿真内容与结果

1.2a)在上述1.1)仿真条件下，对本发明九种实施例阵列配置的差分合成阵列自由度与总阵元数的关系进行仿真，结果如图5所示；

1.2b)在上述1.1)仿真条件下，对本发明九种实施例阵列配置的有效孔径与总阵元数关系进行仿真，结果如图6所示；

1.2c)在上述1.1)仿真条件下，对本发明九种实施例阵列配置的阵列密度与总阵元数关系进行仿真，结果如图7所示。

从图5可以看出，差分合成阵列的自由度随总阵元数的增加而增大；当阵元数目相同时，实施例4具有最高的自由度，实施例1的自由度最低。

从图6可以看出，实施例4具有最高的阵列孔径。随着总阵元数的增大，各实施例的阵列孔径差距也愈加明显。

从图7可以看出，阵列密度随总阵元数的增加而减小；实施例1的阵列密度最高。

仿真2：对本发明五种实施例的阵列配置方法的波达方向DOA估计角度均方根误差与信噪比关系进行仿真。

2.1)仿真条件：阵元数S＝12，信源数Q＝16，探测角度θ，-70^°≤θ≤70^°快拍数T＝500,500次Monte-Carlo实验。

2.2)仿真内容与结果：

在上述2.1)仿真条件下对本发明五种实施例的阵列配置的波达方向DOA估计角度均方根误差与信噪比关系进行仿真，结果如图8所示。

从图8可以看出，DOA估计性能随着信噪比的提高而提高，实施例4具有最小的DOA估计均方根误差。这里只选择五种实施例作比较的原因是总阵元数S＝12时，子阵内阵元数为3的最小冗余阵列MRA和嵌套式阵列NA阵列结构相同，因此实施例2和实施例3结构相同，实施例4和实施例6结构相同，实施例7和实施例9结构相同。

仿真3：对本发明五种实施例的阵列配置的波达方向DOA估计角度均方根误差与快拍数关系进行仿真。

3.1)仿真条件：设阵元数S＝12，信源数Q＝16，探测角度θ，-70°≤θ≤70°，信噪比SNR＝0dB，500次Monte-Carlo实验。

3.2)仿真内容与结果：

在上述3.1)仿真条件下，对本发明五种实施例的阵列配置的波达方向DOA估计角度均方根误差与快拍数关系进行仿真，结果如图9所示。

从图9可以看出，DOA估计性能随着快拍数的增加而提高，实施例4具有最小的DOA估计均方根误差。这里只选择五种实施例作比较的原因是总阵元数S＝12时，子阵内阵元数为3的最小冗余阵列MRA和嵌套式阵列NA阵列结构相同，因此实施例2和实施例3结构相同，实施例4和实施例6结构相同，实施例7和实施例9结构相同。

综上，本发明阵列配置灵活，具有较高的自由度、阵列孔径和阵列密度，较好的波达方向DOA估计测角性能。

Claims (9)

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

1)构造嵌套式子阵阵列：

1a)给定总阵元数S，对S进行因式分解，得到子阵内阵元数目M及子阵数目N；

1b)根据子阵内阵元数目M，设计子阵内阵元结构，同时计算子阵内阵元的位置矢量u_M：

u_M＝[m₁,m₂,…,m_i,…,m_M]·d，

其中，m_i表示子阵内第i个阵元位置系数，1≤i≤M，且m₁＝0，d为入射信号的半波长；

1c)根据子阵数目N，设计子阵结构，同时计算N个子阵的阵位置矢量u_N：

u_N＝[n₁,n₂,…,n_j,…,n_N]·D，

其中，n_j表示第j个子阵位置系数，1≤j≤N，且n₁＝0，D＝L·d为子阵间的最优间隔单元，L表示子阵间最优间隔单元系数；

1d)根据上述设计参数，计算嵌套式子阵阵列位置矢量v：

$v=u_N\⊕u_M=\{n_j\·D+m_i\·d\}=\[p_1,p_2,...,p_k,...,p_S\]\·d,$

其中p_k表示嵌套式子阵阵列第k个阵元的位置系数，k＝1,2,…,S，符号表示交叉求和；

2)根据嵌套式子阵阵列，得到嵌套式子阵阵列的接收数据X(t)；

3)根据阵列接收数据X(t)估计嵌套式子阵阵列的协方差矩阵R_XX，向量化该协方差矩阵得到向量化接收数据z_ct，去除z_ct中重复的元素，得到虚拟的差分合成阵列接收数据z_c：

4)将差分合成阵列接收数据z_c划分为l_c个子阵接收数据，其中z_cg是第g个子阵接收数据，g＝1,2,…,l_c；计算第g个子阵的协方差矩阵R_g，并对该R_g进行前后向平滑平均，得到秩恢复的数据协方差矩阵R_SS，其中l_c表示划分的子阵个数，f_v表示嵌套式子阵阵列的自由度；

5)根据秩恢复的数据协方差矩阵R_SS估计出空间谱函数S(θ)，计算S(θ)中Q个极大值S(θ_q)，该S(θ_q)对应的角度θ_q即为波达方向，其中，θ_q表示第q个入射信号的入射角，q＝1,2,…,Q，Q表示入射信号的总个数。

2.根据权利要求1所述的基于嵌套式子阵阵列的波达方向估计方法，其中步骤1a)中的子阵内阵元数M及子阵数目N，通过按如下公式计算：

S＝M·N，

其中M和N的取值最接近，且M≥2,N≥2。

3.根据权利要求1所述的基于嵌套式子阵阵列的波达方向估计方法，其中步骤1b)中子阵内阵元结构的设计，采用现有均匀线阵ULA结构或最小冗余阵列MRA结构或嵌套式线阵NLA的结构。

4.根据权利要求1所述的基于嵌套式子阵阵列的波达方向估计方法，其中步骤1c)中子阵结构的设计，采用现有均匀线阵ULA结构或最小冗余阵列MRA结构嵌套式线阵NLA的结构。

5.根据权利要求1所述的基于嵌套式子阵阵列的波达方向估计方法，其中步骤1c)中子阵间最优间隔单元系数L的选取，由子阵内阵元结构决定：

当子阵内阵元结构设计为均匀线阵ULA结构时，子阵间最优间隔单元系数为L＝2M-1；

当子阵内阵元结构设计为最小冗余阵列MRA时，子阵间最优间隔单元系数为：

L＝M(M-1)+1-2f_M；

当子阵内阵元结构设计为嵌套式线阵NLA时，根据子阵内阵元数目M的奇偶性确定：

若子阵内阵元数目M为偶数，则子阵间最优间隔单元系数L＝M²/2+M-1；

若子阵内阵元数目M为奇数，则子阵间最优间隔单元系数L＝(M+1)²/2-1。

6.根据权利要求1所述的基于嵌套式子阵阵列的波达方向估计方法，其中步骤2)中的嵌套式子阵阵列的接收数据X(t)，其表示如下：

X(t)＝AS(t)+n(t)，

其中A为导向矢量矩阵，A＝[a(θ₁),a(θ₂),…,a(θ_q),…,a(θ_Q)]，a(θ_q)表示第q个入射信号的导向矢量， 表示第q个入射信号在第k个阵元处的旋转因子，p_k表示嵌套式子阵阵列第k个阵元的位置系数，q＝1,2,…,Q，Q表示入射信号的总个数，S(t)表示入射信号，S(t)＝[s₁(t),s₂(t),…，s_q(t),…，s_Q(t)]^T，s_q(t)表示第q个入射信号，(·)^T:表示向量的转置，n(t)表示均值为0，方差为的复高斯白噪声，且与入射信号不相关。

7.根据权利要求1所述的基于嵌套式子阵阵列的波达方向估计方法，其中步骤3)中虚拟的差分合成阵列接收数据z_c，按如下步骤计算：

7.1)根据嵌套式子阵阵列的接收数据X(t)估计阵列协方差矩R_XX：

R_XX＝E[X(t)X^H(t)]，

其中，其中E表示求数学期望；(·)^H表示共轭转置；

7.2)根据阵列协方差矩R_XX，计算向量化接收数据z_ct：

z_ct＝vec(R_XX)，

其中，vec表示矩阵的向量化操作；

7.3)根据向量化接收数据z_ct，得到虚拟的差分合成阵列接收数据z_c：

z_c＝Dis(z_ct)，

其中，Dis表示去除向量中重复元素的操作。

8.根据权利要求1所述的基于嵌套式子阵阵列的波达方向估计方法，其中步骤4)中秩恢复的数据协方差矩阵R_SS，按如下步骤计算：

8.1)根据虚拟的差分合成阵列接收数据z_c，得到的子阵接收数据z_cg，其表示如下：

z_c＝[z_c1,z_c2,…,z_cg,…,z_clc]，

其中z_cg是第g个子阵接收数据，g＝1,2,…,l_c；

8.2)根据子阵接收数据z_cg，计算第g个子阵的协方差矩阵R_g：

$R_g=z_{cg}{z}_{cg}^H;$

8.3)根据第g子阵的协方差矩阵R_g，计算前向平滑的协方差矩阵R^f：

$R^f=\frac{1}{l_c}\underset{g=1}{\overset{l_c}{\Σ}}{R}_g;$

8.4)根据子阵接收数据z_cg计算第g个子阵的后向接收数据y_cg；

y_cg＝Trans(z_cg)，

其中，Trans表示左右翻转并取共轭；

8.5)根据第g个子阵的后向接收数据y_cg，计算第g个子阵的后向协方差矩阵R_g ^b：

${R_g}^b=y_{cg}{y}_{cg}^H;$

8.6)根据第g个子阵的后向协方差矩阵R_g ^b，计算后向平滑的协方差矩阵R^b：

$R^b=\frac{1}{l_c}\underset{g=1}{\overset{l_c}{\Σ}}{R}_g^b;$

8.7)根据前向平滑的协方差矩阵R^f和后向平滑的协方差矩阵R^b，计算秩恢复的数据协方差矩阵R_SS：

$R_{SS}=\frac{R^f+R^b}{2}.$

9.根据权利要求1所述的基于嵌套式子阵阵列的波达方向估计方法，其中步骤5)中空间谱函数S(θ)，按如下步骤计算：

9.1)对秩恢复的数据协方差矩阵R_SS进行如下特征值分解：

$R_{SS}=\underset{t=1}{\overset{T}{\Σ}}{\mathrm{\λ}}_t{V}_t{V}_t^H,$

其中，λ_t表示矩阵R_SS的特征值，V_t表示矩阵R_SS对应于λ_t的特征矢量，t表示采样时间，t＝1,2,…,T，T表示快拍数，Q表示入射信号的总个数；

9.2)根据特征矢量V_t，计算空间谱函数S(θ)：

$S\left(\mathrm{\θ}\right)=\frac{1}{\left|\underset{t=Q+1}{\overset{T}{\Σ}}a\left(\mathrm{\θ}\right)V_t{V}_t^{H}a\left(\mathrm{\θ}\right)\right|},$

其中，a(θ)表示入射角度为θ的入射信号的导向矢量，p_k表示嵌套式子阵阵列第k个阵元的位置系数；

9.3)根据空间谱函数S(θ)，计算S(θ)中Q个极大值S(θ_q)，S(θ_q)对应的角度θ_q即为波达方向。