CN109061555B - 嵌套阵列下混合相干doa估计方法 - Google Patents

Abstract

本发明公开了一种嵌套阵列下混合相干DOA估计方法，属于自适应阵列信号处理领域的DOA技术。本发明基于阵列非均匀结构的特殊性，在未知阵列信号是否存在相干或相关的情况下，重构出导向矩阵，然后使用基于压缩感知的稀疏重构的LASSO方法求得最终的DOA值。本发明不用提前知道信号是否存在相关或相干，且还能通过估计出的相关矩阵来判断出信号是否存在相关或相干性。本发明充分利用了嵌套阵列的高自由度的优势，使得估计信号的个数超出了物理阵列的孔径。

Description

嵌套阵列下混合相干DOA估计方法

技术领域

本发明涉及自适应阵列信号处理领域的波达方向角(DOA，direction-of-arrival)技术，具体是涉及在嵌套阵列条件下阵列接收信号存在相干的DOA估计方法。

背景技术

目前，绝大部分的国内外学者研究的DOA估计方法基本都是假设在均匀线阵(ULA,uniform linear array)下进行的。然而，在实际的复杂的通信环境中，由于多径反射等因素的存在，往往造成信号源之间相干或相关。此时阵列接收的源信号协方差矩阵将出现缺秩的现象，从而导致某些相干源的方向矢量将无法进行测量。而现有的基于ULA条件下的DOA方法却只能估计出相干源的角度，无法准确地得到是哪些信号源存在相干或相关的。

2010年，Piya Pal和P.P.Vaidyanathan等人提出的一种非均匀阵列结构，其方法是通过接收数据的相关矩阵向量化处理后得到的虚拟均匀阵列结构。由于多径反射等因素的存在，通过非均匀阵列结构得到的信号的相关矩阵R_ss本身是缺秩的，如果再将其转化到虚拟均匀阵列结构上，由于信号的相关矩阵R_ss存在互相关量，结果将会进一步的恶化DOA估计精度，甚至出现DOA值将无法估计的现象。

考虑一个N维远场窄带信号，阵元数为M。假设有K(K≥M)个来自不同方向的信号，其入射角度分别为θ_k,k＝1,2,...,K。理想情况下，n时刻的接收信号x(n)为：x(n)＝A(θ)s(n)+e(n),n＝1,2,…,N，其中A(θ)＝[a(θ₁),a(θ₂),…,a(θ_K)]为大小为M×K的阵列导向矢量矩阵，s(n)＝[s₁(n),s₂(n),…,s_K(n)]^T为n时刻信号的复包络，e(n)为零均值、方差为

的噪声向量，符号(·)^T表示转置。当噪声与源信号不相干，且各个信号源相互独立时，为了求得信号DOA值，一种解决思路是利用相关矩阵

进行向量化处理，其中，E{·}表示期望，A表示阵列导向矢量矩阵，I_M表示单位矩阵，

表示噪声功率，符号(·)^H表示共轭转置，即将非均匀阵列转化到虚拟阵列上，然后通过构造正交完备基矩阵，从而使用LASSO方法进行求解得到DOA值(具体可参考文献：TibshiraniR.Regression shrinkage and selection via the lasso[J].Journal of the RoyalStatistical Society.Series B(Methodological),1996:267–288)。该方案中，相关矩阵求解方法为

当信号存在相干时，即存在混合相干信号源时，此时在n时刻的接收信号x(n)可以重新表示为：

x(n)＝A(θ)s(n)+e(n)＝A₁(θ)s₁(n)+A₂(θ)s₂(n)+e(n),n＝1,2,…,N

其中，

为大小为M×K₁的导向矢量矩阵，

为大小为M×(K-K₁)的导向矢量矩阵，

为n时刻存在相干或相关时的信号的复包络，

为n时刻相互独立的信号的复包络。

此时，由于信号的相关矩阵R_ss＝E{s(n)s^H(n)}引入了互相关量，直接通过上述方法是无法估计出DOA值。目前国内外一些研究学者对于此问题展开了深度研究。SatoshiSHIRAI等人提出使用N-阶根方法进行DOA估计，即将相关矩阵R_X求N-阶根，然后利用求根后的

来进行DOA估计；Y.Hu等人提出基于互质阵列条件下使用四阶累计量来对相干信号进行估计。对于N-阶根方法，可以明显的发现，当信号完全相干时，其估计误差精度非常大；对于基于互质阵列条件下使用四阶累计量的方案，可以看出，如果不引入扩展的互质阵列，误差也很大，并且会出现虚假的峰，计算复杂度也相当的高。并且，两者提出的方法都没有很好的解决此问题：如果信号是相干的，哪些信号之间是相干的。

此外，现有的在ULA条件下，当信号存在相关时，比如：空间平滑方法(SSP，SpatialSmoothing Processing)，特征向量重构方法(Eigenvector Method)，利用相关矩阵重构算法。如果直接将ULA方法使用到嵌套阵列下，由于相关或相干信号的互相关项不为零，嵌套阵的虚拟阵列信号形式发生了变化，上述方法都不在适用。

发明内容

本发明的发明目的在于：针对现有的嵌套阵列DOA方法无法解决信号数超过阵列孔径且信号存在相干的情况下，本发明公开了一种信号存在有相关或相干源的同时存在有独立源时的DOA估计方法。

本发明的嵌套阵列下混合相干DOA估计方法，包括下列步骤：

步骤1：通过嵌套阵列获取接收数据x(n)，并计算接收数据x(n)的相关矩阵R_X：

其中N表示快拍数，即估计的次数，符号(·)^H表示共轭转置；

步骤2：将空间角度均匀划分为D份，并构造稀疏信号的导向矩阵

将感兴趣的信号的角度区域划均匀分为D份，得到D个搜索角度

构造导向矩阵

为：

其中，

表示搜索角度

的导向矢量，i＝1,…,D，符号(·)*表示共轭，符号

表示Kronecker积；

步骤3：对相关矩阵R_X进行向量化处理，得到矩阵

其中

I_M表示单位矩阵，

表示噪声功率，vec(·)表示向量化；

步骤4：根据DOA估计模型

求解向量化的相关矩阵

其中符号||·||₁表示L1范数，符号||·||₂表示L2范数，权值λ_t的取值范围为(0,1)；

噪声功率

的获取方式为：

将待估信号源划分为G组，且G的取值小于所述嵌套阵列的物理阵元数；

对相关矩阵R_X进行奇异值分解，取分解后的前K-G个最小特征值的取值作为噪声功率

其中K表示待估信号源数，即接收信号数；

步骤5：获取待估信号源的DOA估计值，以及获取待估信号源的相关或相干信息：

将步骤4求解得到的相关矩阵

转化为方阵形式，得到矩阵

提取矩阵

的主对角元素并进行谱峰搜索，取峰值对应的搜索角度作为待估信号源的DOA估计值；

提取矩阵

的下三角矩阵p_trial，再将矩阵p_trial的主对角元素置零，得到相干信号的矩阵

分别将矩阵

的每一列上的所有元素相加，得到参数

以及将

的每一行上的所有元素相加，得到参数

再分别对

进行谱峰搜索，若存在谱峰，则对应位置的待估信号源为相干或相关的。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

在嵌套阵列存在相关或相干的情况下,既能解决信号相关或相干问题，也能判断出是哪些信号存在相关或相干的。是一种新的DOA估计方法。本发明充分利用了嵌套阵列的高自由度的优势，使得估计信号的个数超出了物理阵列的孔径，并且，本发明不仅不需要对信号的互相关量进行消除，还利用了互相关量来判断哪些信号是存在相关或相干的。

附图说明

图1为阵元数为6的物理阵列排列示意图，其中，图1-a表示实际物理阵列的排列方式；图1-b表示将非均匀阵列填充为均匀阵列的方式；图1-c表示嵌套阵列条件下的虚拟阵元位置。

图2为相关或相干信号源的位置示意图；

图3为嵌套阵列信号完全相干的条件下阵列的波达方向图。

图4为嵌套阵列信号完全相干的条件下哪些信号是存在完全相干的波达方向图。

图5为嵌套阵列混合相干条件下阵列的波达方向图。

图6为嵌套阵列混合相干条件下哪些信号是存在完全相干的波达方向图。

图7为嵌套阵列存在相干源和独立源的条件下本发明与空间平滑方法的性能对比。

图8为嵌套阵列存在相干源和独立源的条件下N阶根的方法估计。

图9为嵌套阵列存在相干条件下最小均方误差(RMSE)随快拍数变化图。

图10为嵌套阵列存在相干条件下最小均方误差(RMSE)随SNR变化图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图，对本发明作进一步地详细描述。

本发明针对现有的嵌套阵列DOA方法无法解决信号数超过阵列孔径且信号存在相干的情况下，公开了一种信号存在有相关或相干源的同时存在有独立源时的DOA估计方法，该方法可测信号数超过了实际物理阵列个数。

本发明基于嵌套阵列结构的特殊性，在信号存在相干或相关时，重构出向量化后的接收数据相关矩阵，然后利用基于压缩感知的LASSO方法来求得最终的DOA值。本发明不用提前知道信号是否存在相干或者哪些信号是存在相干的，直接利用阵列的接收数据矩阵来估计出信号的DOA值以及哪些信号是存在相干的。在信噪比(SNR)较低时，也依然能够对信号进行准确的估计。

无论信号s(n)是否存在相干或相关性，信号间的相关矩阵R_ss均可表示为R_ss＝E{s(n)s^H(n)}。对于R_ss，其主对角线上的元素为信号的自相关函数，非对角线上的元素为互相关函数。如果信号之间相互独立，则非对角线上的元素均为0；如果信号之间存在相关或相干性，则非对角线上的元素存在不为零的值。由此可知，信号独立为信号相干或相关情况下的特殊情况。

首先，只考虑信号相干时的情况。在嵌套阵列下，阵元数为M。假设有K个来自不同方向的信号，其入射角度分别为θ_k,k＝1,2,...,K。理想情况下，n时刻阵列的接收数据矩阵为x(n)＝A(θ)s(n)+e(n)，其中，噪声与信号统计独立，其相关矩阵R_X为：

如果信号之间存在相关或相干性，则上述等式可转化为：

对于上式，如果信号s_i(n)和s_j(n)之间存在相干或相关性，则ρ_ij不等于0；如果信号s_i(n)和s_j(n)之间相互统计独立，则ρ_ij等于0。

将上述嵌套阵列转化到虚拟均匀阵列上，即将上述相关矩阵向量化后，可以得到：

其中，p₁＝[ρ₁₁ ρ₂₂ … ρ_KK]^T，p₂＝[ρ₁₂ … ρ_1K … ρ_K1 … ρ_K,K-1]^T。vec(·)表示向量化；

为了能够直观的使用基于压缩感知的LASSO方法来进行求解，将上述等式转化为如下形式：

其中，

为了能够估计出相干信号，采用基于压缩感知的稀疏重构思想，将感兴趣的信号的角度区域均匀划分为D份，得到搜索角度集合

目标假设在划分的角度区域内。故此，虚拟阵列下的接收数据矩阵为

其中，

则，本发明的加权l₁-范数模型的DOA估计求解模型即为下面的形式：

其中β表示无限接近于零的参数。

将加权l₁-范数模型的DOA估计求解模型转化为LASSO方法进行上述估计，有：

通过CVX工具箱(Matlab中凸优化工具包)进行求解后，可以得到向量化后的

然后，将其转化为方阵的形式。即：对于上述

得知，其矩阵的维度为D²×1，故此，可以将矩阵

每隔D个组成一个行向量，一共有D个行向量，然后按照顺序进行组合，便可以得到方阵

提取

主对角元素，通过谱峰搜索(即频谱画图)，取峰值对应的搜索角度作为待估信号源的DOA估计值，从而得到需估计的DOA值。

在具体应用中，虽然存在信号的相干，但是实际应用中有噪声的存在，对R_X进行奇异值分解后，其特征值是不会存在等于0的情况。如果将K个接收信号划分为G个组，组与组之间相互独立，组内是相干的。当划分的组数G小于物理阵元数时，对R_X进行奇异值分解后，其前(K-G)个最小特征值对应的为信号的噪声功率，可以将其取平均得到信号的噪声功率

当划分的组数G大于或等于物理阵元数时，则不可以通过奇异值分解得到信号的噪声。对于使用LASSO方法，λ_t的优选值取为0.25。

通过上述方法，便可以求出信号的DOA值。然而，嵌套阵列的自由度虽然很高，但是真正能估计出的信号的DOA值的个数却没有那么高。对阵列所能处理的信号最大数以及哪些信号是相干的分析如下：

1.考虑信号完全相干的情况。

对于上述的方法，可分辨的目标最大的数目由物理阵列和信号的相干数决定。这是由于对于每一对相干信号，通过R_ss可以明显的发现，R_ss非对角线上的元素都会有两个互相干或相关量，并且每一个信号在主对角元素上也有一个自相关量。由于R_ss的共轭对称性，可以只考虑R_ss的下三角元素。也就是说，估计的R_ss的元素的个数从K²变为

故此，可以得出，分辨出的目标最大的数目满足的关系为：L_PA≥2L_VA，其中，L_PA在本发明中表示为物理阵列的个数，L_VA表示R_ss中最少所需要的相关或相干的系数，也就是说，在信号完全相干或相关的情况下，利用其中一个信号来得到与其他信号相关或相干的对数。即L_VA＝M+P，P表示最少所需要的相关或相干的系数的个数，M表示主对角元素的个数。

例如：对于存在三个完全相干的信号，则有

(只考虑下三角元素)，L_VA＝M+P＝3+2＝5，故L_PA≥2L_VA＝2×5＝10。

2.考虑信号混合相干或相关的情况。

嵌套阵列是利用相关矩阵来进行扩展信号的自由度的，下表1给出了嵌套阵列的阵元数M分别为奇数和偶数时的最大自由度。实际上，对于本发明，嵌套阵列所能处理的信号数并没有达到嵌套阵列的自由度所给出的那么高。这是由于虚拟阵列也是一个均匀阵列，也需要满足上述相干情况下分辨出的目标最大的数目满足的关系。对于嵌套阵列，是在原来的物理阵列上，以第一级嵌套阵列的阵元间距

将第一个物理阵元到最后一个阵元之间补齐，然后将此物理阵列关于参考阵元0点位置做左右对称，然后再加上参考阵元0。这样便可以得到虚拟阵列。

表1阵元数M分别为奇数和偶数时的最大自由度(DOF)

M	最优M<sub>1</sub>,M<sub>2</sub>	DOF
			偶数	M<sub>1</sub>＝M<sub>2</sub>＝M/2	(M<sup>2</sup>-2)/2+M
奇数	M<sub>1</sub>＝(M-1)/2,M<sub>2</sub>＝(M+1)/2	(M<sup>2</sup>-1)/2+M

例如：对于阵元数为6的物理阵列，其排列方式如图1所示，对于此时的混合相干或相关信号，可以将存在相干或相关的信号划分为G个组，组与组之间是满足统计独立的，组内是相干或相关的。虽然此时是混合相干或相关，但是各组内相干或相关信号数也需要满足物理阵列下的条件，即L_PA≥2L_VA。由于将信号分为了G组，但是，组与组之间是相互独立的，故此，通过求解各组中的相干或相关条件数，即各组中的相关或相干为：2L_VAg，然后将各组的相关或相干条件数相加得到最终的值2L_VAG。由于是通过嵌套阵列来进行求解的，所以，在转化为虚拟阵列后，其阵列的自由度增加了，如果阵列数为奇数，那么此时的虚拟阵列的最大自由度为L_VIR＝(M²-1)/2+M；如果阵列数为偶数，那么此时的虚拟阵列的最大自由度为L_VIR＝(M²-2)/2+M。所以最终的混合相干或相关信号满足的条件为：L_VIR≥2L_VAG。如果某一组只存在统计独立信号的话，可以将其看成是自己与自己相干或相关的信号。

综上所述，对于混合相干或相关的信号满足的条件为：

(1)各组内相干或相关信号数需要满足物理阵列下的条件，即L_PA≥2L_VA；

(2)所有组的混合相干或相关满足虚拟阵列下的条件，即：L_VIR≥2L_VAG。

由于是通过阵列孔径来进行分析的，有时候由于估计的角度不合适，就存在这样的情况：在刚好达到上述条件下，存在测不出来的现象。故此，在具体应用中，应将阵元数调节的高一点。

3.信号之间相干或相关分析。

通过上述方法，可以得到主对角线上的元素，然后利用主对角线上的元素，通过谱峰搜索，得到所有信号源所在的位置。然而，没有利用非对角线上的元素，如果利用非对角线上的元素，便可以得到各个信号之间哪些信号是相干或相关的。

现在，考虑对于存在K个信号源相干或相关的情况。如果信号存在相关或相干，那么，

对角线以下的元素存在非零值，并且，此时的非零值为信号源之间的相关或相干系数。一个相关或相干系数的行与列指向两个信号源，并且，这两个信号源是存在相关或相干的。故此，可以通过基于压缩感知的稀疏重构的LASSO方法求得的

来提取

对角线以下的元素，得到新的相关矩阵

然后，将相关矩阵

的每一列进行相加，便可以判断出哪些列中所对应的信号源存在与其他信号源相干或相关；将相关矩阵

的每一行进行相加，便可以判断出哪些行中所对应的信号源存在与其他信号源相干或相关。最后通过相关矩阵

的行与列来进行分析，从而可以得出哪些信号源是相干或相关的。

例如，对于图2，如果将感兴趣的信号的角度区域划分为D份，其中第2份、第3份和第4份是存在相干或相关的。则可以得知，图2中的灰色区域是存在互相关系数的，通过上述方法便可以得到哪些信号源是存在相干或相关的。

具体执行步骤如下：通过提取

的下三角矩阵p_trial，然后将p_trial的主对角元素置零，便可以得到相干信号的矩阵

通过

分别将

的每一列上的所有元素相加，得到

通过谱峰搜索，可以判断出哪些列中所对应的信号源存在与其他信号源相干或相关；通过

分别将

的每一行上的所有元素相加，得到

通过谱峰搜索，可以判断出哪些行中所对应的信号源存在与其他信号源相干或相关。即，通过相关矩阵

的行与列来进行分析谱峰分析，谱峰位置所对应的信号源为相干或相关信号源。

为了表明本发明在阵元存在相关或相干的条件下本发明DOA估计性能优势，通过仿真来证明此方法的性能优势。

仿真参数1：本次实验考虑2个完全相干的信号的情况。采用本发明(proposed)的嵌套阵列思想进行估计，快拍数为1000，信号角度在[-60°，60°]上均匀分布，共计2个信号，阵元数为8个，物理阵列的位置为{1,2,3,4，5，10，15，20}d。信号的信噪比(SNR)均为SNR＝10dB，信号的相干系数为1(即信号是完全相干的)，实验仿真结果如图3，4所示；

为了仿真混合相干的情况，选择信号角度在[-60°，60°]上均匀分布，共计12个信号，其中第i(i＝1,2,3,4,5,6)个信号与第(i+6)个信号是完全相干的，其余的条件不变，本发明与空间平滑方法(SSP，Spatial Smoothing Processing)进行了仿真对比，仿真结果如图5，6所示；

为了仿真同时存在相干源和独立源的情况，选择信号角度在[-60°，60°]上均匀分布，共计10个信号。相干源分为两组，第一组，即第1,2,3个信号是完全相干，第二组，即第4，5个信号是完全相干的。其余的信号为统计独立信号，其余的条件不变，本发明与N-阶根进行了仿真，仿真结果如图7，8所示。

从图3可以发现，选择使用阵元数为8的嵌套阵列，便可以进行估计信号的DOA值。故此，此次实验仿真参数的设置直接证明了信号完全相干的条件下信号数的设置；从图4，6可以得到，当信号完全相干的情况下，可以通过相关矩阵

的行与列来进行分析，从而判断出是哪些信号存在相关或相干的。

从图5，7，8，发现，当信号混合相干或者同时存在相干源和独立源时，本发明也可以进行估计。通过图5，可以看到使用SSP方法无法估计出；通过图8的N-阶根方法，虽然可以估计信号的独立源，但是对于相干源却无法进行估计。本发明也直接证明了在信号混合相干的情况下，信号数的设置。此外，本次实验仿真也充分利用了嵌套阵列的高自由度的优势。

仿真参数2：本次实验考虑2个完全相干的信号。采用嵌套阵列思想进行估计，信号角度为[-10°，20°]，共计2个信号，阵元数为10个，物理阵列的位置为{1,2,3,4，5,6,12,18,24,30}d。当考虑最小均方误差(RMSE，Root Mean Square Error)随快拍数变化时，SNR选择为SNR＝10dB；当考虑RMSE随SNR变化时，快拍数选择为1000。划分的网格精度为0.2。为消除实验的随机性对实验的影响，本次实验的蒙特卡罗次数为30次。实验仿真结果如图9，10所示。

从图9，图10的仿真实验分析得出，随着SNR的由小变大或者快拍数的由少变多，本发明可以很好的估计出信号的DOA值。此外，通过图9，10可以得知，当SNR、快拍数很低的条件下，RMSE的值很小，也就是说，本发明依然可以很好的估计出信号的DOA。

以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合。

Claims (3)

1.嵌套阵列下混合相干DOA估计方法，其特征在于，包括下列步骤：

步骤1：通过嵌套阵列获取接收数据x(n)，并计算接收数据x(n)的相关矩阵R_X：

其中N表示快拍数，符号(·)^H表示共轭转置；

步骤2：将空间角度均匀划分为D份，并构造稀疏信号的导向矩阵

将感兴趣的信号的角度区域划均匀分为D份，得到D个搜索角度

构造导向矩阵

为：

其中，

表示搜索角度

的导向矢量，i＝1,…,D，符号(·)^*表示共轭，符号

表示Kronecker积；

步骤3：对相关矩阵R_X进行向量化处理，得到矩阵

其中

θ₁，θ₂，…，θ_K表示K个来自不同方向的信号的入射角度，a(θ₁)，a(θ₂)，…，a(θ_K)表示各入射角度的导向矢量；

ρ₁₁，ρ₁₂，…，ρ_KK表示信号间的相关矩阵的元素；I_M表示单位矩阵，

表示噪声功率，vec(·)表示向量化；

步骤4：根据DOA估计模型

求解向量化的相关矩阵

其中，相关矩阵

的维度为D²×1，符号||·||₁表示L1范数，符号||·||₂表示L2范数，权值λ_t的取值范围为(0,1)；

噪声功率

的获取方式为：将待估信号源划分为G组，且G的取值小于所述嵌套阵列的物理阵元数；

对相关矩阵R_X进行奇异值分解，取分解后的前K-G个最小特征值的取值作为噪声功率

其中K表示待估信号源数，即接收信号数；

步骤5：获取待估信号源的DOA估计值，以及获取待估信号源的相关或相干信息：

将步骤4求解得到的相关矩阵

转化为方阵形式，得到矩阵

提取矩阵

的主对角元素并进行谱峰搜索，取峰值对应的搜索角度作为待估信号源的DOA估计值；

提取矩阵

的下三角矩阵p_trial，再将矩阵p_trial的主对角元素置零，得到相干信号的矩阵

分别将矩阵

的每一列上的所有元素相加，得到参数

以及将

的每一行上的所有元素相加，得到参数

再分别对

进行谱峰搜索，若存在谱峰，则对应位置的待估信号源为相干或相关的。

2.如权利要求1所述的方法，其特征在于，相关矩阵

转化为方阵具体为：

将相关矩阵

每隔D个组成一个行向量，得到D个行向量，然后按照顺序进行组合，得到相关矩阵

的方阵形式的矩阵

3.如权利要求1所述的方法，其特征在于，权值λ_t的优选值取为0.25。

Images