CN106019213B - 一种部分稀疏l阵及其二维doa估计方法 - Google Patents

Abstract

本发明公开了一种部分稀疏L阵及其二维DOA估计方法。本发明包括两个子阵构成L阵，第一子阵的阵元间距为半倍波长，第二子阵的阵元间距为n倍波长；并在第二子阵上设置一个距参考阵元的距离为半倍波长的辅助阵元。在进行DOA估计处理时，将子阵分别置于x轴和z轴，利用互相关矩阵不受噪声影响的特性，求接收数据的互相关并提取信号子空间；利用ULA的平移不变性和信号子空间，求取z轴阵列流型矩阵的旋转矩阵，并通过对旋转矩阵进行特征值分解，得到可能的俯仰角估计值，再基于辅助阵元对估计结果去模糊处理；再基于其进行信源波形的估计及x轴的阵列流型矩阵，求取相应的方位角。本发明用于雷达，声呐等，其实现低成本，低运算量，测向精度高。

Description

一种部分稀疏L阵及其二维DOA估计方法

技术领域

本发明属于无线移动通信技术领域，特别是涉及一种利用线阵构造的部分稀疏L阵列及其二维波达方向(DOA)估计方法。

背景技术

空分多址技术是继时分复用、频分复用和码分复用之后又一个可提高网络容量的关键技术，空分多址的关键在于精确地测量信号源方位和形成高定向波束，因此，以DOA估计和数字波束形成为核心的阵列信号处理技术越来越得到广大研究者的研究和关注。

现有的DOA估计算法大多是基于传统满阵，即阵列中相邻阵元间距不大于入射信号半波长的天线阵列。但满阵由于阵元间距的限制，要增加阵列孔径，提升测向精度和分辨率就必须增加阵元数，造成系统过于复杂和系统成本增加。鉴于传统满阵存在的上述问题，人们提出了稀疏阵，即阵元间距大于半波长的均匀或非均匀阵列系统。和常规满阵相比，阵元数相同时它拥有更大的阵列孔径，减小了阵元间的互耦，改善了测向精度、测向分辨率和最大可处理信源数等性能；此外，阵列孔径相同时稀疏阵所需阵元数更少，这就意味着更小规模的接收系统和信号处理系统等，极大地降低了设备成本。

目前基于稀疏阵列的DOA估计主要分为两大类：一类是基于稀疏线阵的一维DOA估计；另一类是基于由几个稀疏线阵组成的简化面阵的二维DOA估计。其中，基于简化面阵的二维DOA估计中，L阵由于有更大的有效孔径、更小的运算量、更易实现、更强的方法适用性等优点得到了广泛的关注和应用。但现有的L阵大多是由常规ULA(均匀线阵)构成，利用稀疏线阵构成的L阵还比较少，已有的基于稀疏线阵的L阵的表达式不够简洁，且在计算二维DOA时，要么需要配对算法，要么就是需要谱峰搜索等，不能充分利用稀疏阵列的优势。

发明内容

本发明的发明目的在于：针对上述存在的问题，提供一种构造简单，基于线阵的部分稀疏L阵及其能获得较好性能且计算简单的DOA估计方法。

本发明的一种部分稀疏的L阵列，包括阵元数相同的线性第一子阵、第二子阵构成L字形的阵列，第一子阵和第二子阵的交点定义为参考阵元，其特征在于，所述第一子阵的阵元间距等于半倍波长，第二子阵的阵元间距等于n倍波长，其中整数n的取值范围为1～10；并在第二子阵上设置一个辅助阵元，所述辅助阵元距离参考阵元的距离为半倍波长。其中L阵列的第一子阵、第二子阵可分别放置于x轴和z轴，x轴对应方位角，z轴对应俯仰角；或者将L阵列的第一子阵、第二子阵分别放置于z轴和x轴，x轴对应方位角，z轴对应俯仰角。

同时，本发明还公开了一种用于所述的L阵列的二维波大方向的估计方法，包括下列步骤：

步骤1：将L阵列的第一子阵、第二子阵分别放置于x轴和z轴；

步骤2：L阵列接收K个不相关信源的入射信号，得到x轴、z轴各阵元的接收数据，其中K小于第一子阵的陈元数M；

步骤3：计算各入射信号的估计俯仰角；

步骤301：计算接收数据z′(t)和x′(t)在N(N的取值为大于等于10的整数)次采样下的互相关矩阵R_z′x′，并从R_z′x′中任意提取K列构造信号子空间U_z，其中接收数据x′(t)为x轴上除参考阵元外的所有陈元的接收数据，接收数据z′(t)为z轴上除辅助阵元外的所有阵元的接收数据；

将信号子空间U_z划分为上下两个(M-1)*K维的信号子空间U_z1和U_z2，基于U_z1和U_z2求取旋转矩阵Ω_z，令矩阵矩阵T_z为K×K的非奇异矩阵；

对矩阵F_z进行特征值分解，得到相应的K个特征值其中k＝1,2,…,K；

步骤302：计算每个入射信号的两个备选估计俯仰角：

根据计算K个第一备选估计俯仰角其中d_z＝nλ，λ表示信号波长，符号angle(·)表示取相位角；

根据公式计算估计空间相位差若则第二备选估计俯仰角否则

步骤4：对每个入射信号的备选估计俯仰角去模糊角处理，得到估计俯仰角

步骤401：计算接收数据z(t)在N次采样下的自相关矩阵R_zz，其中接收数据z(t)为z轴上所有阵元的接收数据；

步骤402：对R_zz进行特征值分解：其中U_s和U_n分别为(M+1)×K和(M+1)×(M+1-K)维的信号子空间和噪声子空间；Σ_s为从R_zz的M+1个特征值中，取前K个最大特征值组成对角矩阵Σ_s，剩余的M+1-K个特征值组成对角矩阵Σ_n，符号(·)^H表示矩阵共轭；

步骤403：根据公式计算每个入射信号的备选估计俯仰角在z轴的导向矢量其中包括第一备选估计俯仰角和第二备选估计俯仰角 为对应的估计空间相位差，其中e为自然常数，符号(·)^T表示矩阵转置；

将导向矢量与噪声子空间U_n正交的备选估计俯仰角作为入射信号k的估计俯仰角k＝1,2,…,K；

步骤5：基于K个估计俯仰角得到z轴的阵列流型矩阵

根据接收数据z(t)、阵列流型矩阵计算估计信号源即其中符号(·)⁺表示M-P广义逆；

步骤6：估计各入射信号的方位角

计算接收数据Z″(t)和X(t)在N次采样下的互相关矩阵R_xz″，其中接收数据Z″(t)为z轴上除参考阵元的所有阵元的接收数据；接收数据X(t)为x轴上所有阵元的接收数据；

根据公式计算x轴的阵列流型矩阵的估计值其中估计值 为估计信号源在N次采样下的自相关矩阵；

根据公式得到每个入射信号的方位角其中k＝1,2,…,K，表示矩阵的第k列的前M-1行元素的子向量，表示矩阵的第k列的后M-1行元素的子向量。

综上所述，由于采用了上述技术方案，本发明的有益效果是：可以实现低成本，低运算量，高测向精度的二维DOA估计，该方法可以应用于雷达，声呐及无线通信等领域，解决现有常规整理所不能解决的问题。

附图说明

图1是本发明提出的由稀疏ULA和常规ULA构造的部分稀疏L阵的阵列结构；

图2是本发明所提部分稀疏L阵的z轴阵元配置；

图3是利用本发明所提算法估计二维DOA的角度散布图；

图4是本发明所提算法和现有的算法的性能随信噪比变化的对比图；

图5是本发明所提算法和现有的算法的性能随采样快拍数变化的对比图。

具体实施方式

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

本发明为了利用稀疏阵有更大的阵列孔径和均匀线阵的平移不变性，提出一种基于稀疏ULA和常规ULA的L阵结构，然后利用互相关矩阵不受噪声影响的特性，ULA的平移不变性，解模糊算法和LS(最小二乘法)技术求解二维DOA。本发明可用于雷达，声呐，无线通信等领域，为了更好地理解本发明，以下结合附图对本申请作进一步描述。

步骤1：构建稀疏L阵：

本实施方式中，设第一、二线阵的阵元个数均为M，将第一、二线阵的交点阵元定义为参考阵元，其中，第一子阵的陈元间距d₁＝λ/2，即常规的ULA；第二子阵的陈元间距d₂＝nλ，其中λ表示信号波长，本实施例中取n＝1，即稀疏ULA。参见图1，将第一子阵放置于x轴、第二子阵放置于z轴(疏ULA和常规ULA在坐标轴上的位置可以互换)，参考阵元所在位置定义为原点。且在第二子阵设置辅助阵元，其与原点处的参考阵元间距为d₁。其中z轴的系数ULA用来估计俯仰角，x轴的常规ULA用来估计方位角。

步骤2：构建系统模型，从而得到本发明的DOA估计的理论支持：

设K个(K<M)个远场窄带不相关信号S(t)＝[s₁(t),...,s_K(t)]^T从K个方向(θ_i,φ_i),i＝1,...,K入射到本发明的发明L阵上，参见图2，其x轴和z轴的阵列流型矩阵A_z和A_x分别为：A_x＝[a_x(φ₁),...,a_x(φ_K)]、A_z＝[a_z(θ₁),...,a_z(θ_K)]，其中， 分别为第k个信号在x轴和z轴的导向矢量，其中d_x＝d₁，d_z＝d₂。考虑系统噪声n_x(t)＝[n_x,0(t),...,n_x,M-1(t)]^T和n_z(t)＝[n_z,0(t),...,n_z,M(t)]^T分别为x轴和z轴的M×1维和(M+1)×1维的独立同分布的加性高斯白噪声。从而，可将x轴和z轴的阵元在t时刻的接收数据矢量分别表示为：

其中，L阵在原点处的参考阵元是x轴和z轴所共有的，因此有x₀(t)＝z₀(t)，n_x,0(t)＝n_z,0(t)。

步骤3：求俯仰角：

在得到本发明的L阵得到接收数据X(t)和Z(t)之后，为了充分利用ULA的平移不变性和互相关矩阵不受加性噪声影响的特性，选取图2所示的x轴第1到M-1个阵元的接收数据x′(t)＝[x₁(t),...,x_M-1(t)]^T和z轴除辅助阵元外其余阵元的接收数据z′(t)＝[z₀(t),z₂(t),...,z_M(t)]^T，求N(取值范围为20～1000)次采样下，x′(t)和z′(t)的互相关矩阵R_z′x′，即：

由上式可知：互相关矩阵R_z′x′的每一列都是A_z的行的线性组合，符号“A(a:b,:)”表示对应矩阵A的第a至第b行。本发明中选取R_z′x′的前K列来构造信号子空间U_z，则存在一个K×K的非奇异矩阵T_z，满足：U_z＝A_zT_z。

在得到信号子空间U_z后，将U_z划分为上下两个(M-1)×K维的信号子空间U_z1和U_z2，即：U_z1＝U_z(1:M-1,:)＝A_z1T_z＝A_z(1:M-1,:)T_z、U_z2＝U_z(2:M,:)＝A_z2T_z＝A_z(2:M,:)T_z。

根据U_z1和U_z2，利用z轴稀疏ULA的平移不变性，有其中，旋转矩阵Ω_z表示包含K个入射信号俯仰角的所有信息的对角阵。令则可得U_z2＝U_z1F_z。从而，由于矩阵F_z和旋转矩阵Ω_z相似，因此，F_z的特征值是Ω_z的对角线元素的排列。通过对F_z进行特征值分解，得到相应的K个特征值λ₁,λ₂,…,λ_K，再计算每个入射信号的两个备选估计俯仰角：

根据计算K个第一备选估计俯仰角其中d_z＝nλ，λ表示信号波长，符号angle(·)表示取相位角；

根据公式计算估计空间相位差若则第二备选估计俯仰角否则

步骤4：利用辅助阵元解模糊，即对每个入射信号的备选估计俯仰角去模糊角处理，得到估计俯仰角

由于信号子空间和噪声子空间正交，因此，本发明中，先利用z轴所有的接收数据，求其自相关矩阵然后通过特征值分解得到：其中，U_s和U_n分别为(M+1)×K和(M+1)×(M+1-K)维的信号子空间和噪声子空间；Σ_s为从R_zz的M+1个特征值中，取前K个最大特征值组成对角矩阵Σ_s，剩余的M+1-K个特征值组成对角矩阵Σ_n；

根据公式计算每个入射信号的备选估计俯仰角在z轴的导向矢量其中包括和 为对应的估计空间相位差，其中e为自然常数；

将导向矢量与噪声子空间U_n正交的备选估计俯仰角作为入射信号k的估计俯仰角k＝1,2,…,K，即将每个信号的备选估计俯仰角中的不正交删除。

步骤5：利用已得的俯仰角与LS技术估计信源波形。

本发明中通过求解信源波形建立俯仰角和方位角之间一一对应的关系，可以实现俯仰角和方位角的自动配对，不需要额外的配对过程。在求解信源波形时，首先基于K个估计俯仰角得到z轴的阵列流型矩阵即其中然后将求解的问题转化为最小二乘问题：从而可得：

步骤6：估计方位角。

先求解X(t)和Z″(t)＝[z₁(t),z₂(t),...,z_M(t)]^T的互相关矩阵R_xz″，即：

然后利用最小二乘法求A_xR_s的估计值，即：

由于方位角的所有信息都包含在A_x中，因此，求方位角的问题就转变为求解A_x的估计值的问题，因而本发明在得到后，先求解其自相关矩阵即然后利用和求x轴阵列流型矩阵A_x的估计值，即：通过对的每一列应用平移不变性，可以求得每一列的旋转因子，从而可以得到相应的方位角φ_i，即：i＝1,…,K。

另外，本发明的第一子阵和第二子阵所在的位置可以互换，即将L阵列的第一子阵、第二子阵分别放置于z轴和x轴，x轴对应方位角，z轴对应俯仰角。则在DOA估计的处理过程中，采用上述相同的原理，先求取第二子阵所在的x轴的信号子空间U_x，然后再估计每个入射信号的两个备选估计方位角，进行去模糊处理，基于每个信号唯一对应的方位角得到其在x轴的阵列流型矩阵，再基于信号波形的估计求取z轴的阵列流型矩阵，进而基于已估计的方位角得到估计俯仰角，实现DOA估计。

图3是取K＝2，M＝5，N＝200，λ＝0.8m，SNR＝5dB，(θ₁,φ₁)＝(55°,60°)，(θ₂,φ₂)＝(65°,70°)时的仿真结果图，由该角度散布图可知，本发明的俯仰角和方位角比较集中的分布于真实值附近。

图4、5是本发明和现有方式(其中ULA-CCM是文献“DOA Estimation using cross-correlation matrix”中所采用的方案，JSVD是文献“Joint svd of twocrosscorrelation matrices to achieve automatic pairing in 2-D angleestimation problems”中所采用的方案)的性能随信噪比变化的对比图、以及随采样快拍数变化的对比图。图4分别是第一个俯仰角和方位角的估计误差随信噪比的变化曲线，其中，K＝2，M＝5，N＝200，λ＝0.8m，SNR＝(0～40)dB，(θ₁,φ₁)＝(50°,60°)，(θ₂,φ₂)＝(70°,50°)，由图4可知，采用本发明提出的稀疏L阵和相应的二维DOA估计算法之后，俯仰角的性能可以比对比算法提升大约15dB，而方位角可以提升大约10dB。图5分别是第一个俯仰角和方位角的估计误差随采样快拍数N的变化曲线，其中，K＝2，M＝5，N＝200，λ＝0.8m，SNR＝5dB，(θ₁,φ₁)＝(50°,60°)，(θ₂,φ₂)＝(70°,50°)，由图5可知，采用本发明提出的部分稀疏L阵和相应的二维DOA估计方法之后，俯仰角和方位角的性能均可以比对比算法提升至少12dB。

因此，本发明所提的新的L阵列及其相应的二维DOA估计方法能够很好的提高二维DOA估计的测向精度，和已有的方式相比，本发明以较少的成本和计算量实现了比已有方式更好的性能。

Claims (2)

1.一种部分稀疏L阵列的二维波达方向的估计方法，其特征在于，包括下列步骤：

步骤1：设置部分稀疏L阵列：阵元数相同的线性第一子阵、第二子阵构成L字形的阵列，第一子阵和第二子阵的交点定义为参考阵元，所述第一子阵的阵元间距等于半倍波长，第二子阵的阵元间距等于n倍波长，其中整数n的取值范围为1～10；并在第二子阵上设置一个辅助阵元，所述辅助阵元距离参考阵元的距离为半倍波长；

将L阵列的第一子阵、第二子阵分别放置于x轴和z轴；

步骤2：L阵列接收K个不相关信源的入射信号，得到x轴、z轴各阵元的接收数据，其中K小于第一子阵的阵元数M；

步骤3：计算各入射信号的估计俯仰角；

步骤301：计算接收数据z′(t)和x′(t)在N次采样下的互相关矩阵R_z′x′，并从R_z′x′中任意提取K列构造信号子空间U_z，其中接收数据x′(t)为x轴上除参考阵元外的所有阵元的接收数据，接收数据z′(t)为z轴上除辅助阵元外的所有阵元的接收数据；

将信号子空间U_z划分为上下两个(M-1)*K维的信号子空间U_z1和U_z2，基于U_z1和U_z2求取旋转矩阵Ω_z，令矩阵F_z＝T_z ^-1Ω_zT_z，矩阵T_z为K×K的非奇异矩阵；

对矩阵F_z进行特征值分解，得到相应的K个特征值其中k＝1,2,…,K；

步骤302：计算每个入射信号的两个备选估计俯仰角：

根据计算K个第一备选估计俯仰角其中d_z＝nλ，λ表示信号波长，符号angle(·)表示取相位角；

根据公式计算估计空间相位差若则第二备选估计俯仰角否则

步骤4：对每个入射信号的备选估计俯仰角去模糊角处理，得到估计俯仰角

步骤401：计算接收数据z(t)在N次采样下的自相关矩阵R_zz，其中接收数据z(t)为z轴上所有阵元的接收数据；

步骤402：对R_zz进行特征值分解：其中U_s和U_n分别为(M+1)×K和(M+1)×(M+1-K)维的信号子空间和噪声子空间；Σ_s为从R_zz的M+1个特征值中，取前K个最大特征值组成对角矩阵Σ_s，剩余的M+1-K个特征值组成对角矩阵Σ_n；

步骤403：根据公式计算每个入射信号的备选估计俯仰角在z轴的导向矢量其中包括第一备选估计俯仰角和第二备选估计俯仰角 为对应的估计空间相位差，其中e为自然常数；

将导向矢量与噪声子空间U_n正交的备选估计俯仰角作为入射信号k的估计俯仰角

步骤5：基于K个估计俯仰角得到z轴的阵列流型矩阵

根据接收数据z(t)、阵列流型矩阵计算估计信号源即其中符号(·)⁺表示M-P广义逆；

步骤6：估计各入射信号的方位角

计算接收数据Z″(t)和X(t)在N次采样下的互相关矩阵R_xz″，其中接收数据Z″(t)为z轴上除参考阵元的所有阵元的接收数据；接收数据X(t)为x轴上所有阵元的接收数据；

根据公式计算x轴的阵列流型矩阵的估计值其中估计值 为估计信号源在N次采样下的自相关矩阵；

根据公式得到每个入射信号的方位角其中k＝1,2,…,K，表示矩阵的第k列的前M-1行元素的子向量，表示矩阵的第k列的后M-1行元素的子向量。

2.如权利要求1所述的方法，其特征在于，采样数N的取值为大于等于10的整数。