CN108872930B - 扩展孔径二维联合对角化doa估计方法 - Google Patents

Abstract

本发明涉及雷达技术领域，为解决传统的二维DOA估计中，在俯仰角为70°～90°的实际移动通信俯仰角度范围内的角度估计失效问题、俯仰角和方位角的配对问题，以及奇异点问题，本发明，扩展孔径二维联合对角化DOA估计方法，利用非均匀L型天线阵列实现，具体步骤如下：(1)构造延时互相关矩阵；(2)构造选择矩阵，分别得到四个对角矩阵；(3)得到z轴低精度无模糊的方向余弦估计值以及高精度模糊的方向余弦估计值；(4)同理得到x轴上对应的两种方向余弦估计值；(5)通过解模糊方法得到高精度无模糊的方向余弦估计值

(6)得到俯仰角和方位角的估计值

和

本发明主要应用于估计接收到的信号的到达方向场合。

Description

扩展孔径二维联合对角化DOA估计方法

技术领域

本发明涉及采用阵列天线估计接收到的信号的到达方向的技术领域，尤其涉及采用扩展孔径的非均匀L型天线阵列的信号到达方向估计方法。

背景技术

空间信号到达方向(Direction of Arrival,DOA)估计是空间谱估计一个主要研究方向，被广泛应用在雷达、声呐、地震、通信等许多领域。DOA估计的基本问题就是确定各个信号到达阵列参考阵元的方向角，简称波达方向。经典的子空间分解类DOA估计算法有多重信号分类算法(MUSIC,Multiple Signal Classification)和基于旋转不变技术的信号参数估计算法(ESPRIT,Estimation of Signal Parameter via Rotational InvitationTechniques)。其中MUSIC算法是噪声子空间类算法，ESPRIT算法是信号子空间类算法，改进的MUSIC算法包括特征矢量法、求根MUSIC法、加权MUSIC算法等，改进的ESPRIT算法包括最小二乘ESPRIT、总体最小二乘ESPRIT、加权ESPRIT算法等。

传统的MUSIC算法和ESPRIT算法等高分辨率算法，虽然具有良好的估计性能，但是由于需要对接收信号协方差矩阵进行特征值分解，因此具有较大的计算量。传播算子算法使用线性运算使用线性运算代替了奇异值分解和特征值分解运算，显著地降低了计算复杂度。由于传播算子算法具有计算复杂度较低的优点，各国学者们对其进行了广泛研究，并提出大量基于传播算子的DOA估计算法。目前，存在大量基于传播算子的L型阵列、2-L型阵列、双平行线阵、三平行线阵等二维DOA估计算法。但是某些基于双平行线阵的传播算子算法在俯仰角为70°～90°的实际移动通信俯仰角度范围内存在角度估计失效问题，有些基于三平行线阵采用传播算子的二维DOA估计算法并没有充分利用所有的阵元信息。L型阵列形式简单且能够提供较好的角度估计性能，因此大量基于L型阵列的二维DOA估计算法被提出。扩展孔径可以有效地提高阵列的分辨率和角度估计精度，但会出现模糊的现象。有些算法提出了解模糊的算法，但是由于算法本身是基于ESPRIT算法，因此具有较大的计算复杂度。另一个值得注意的问题是阵元间距大于半波长的非均匀阵列用于二维波达方向估计时，即使无相同的方位角或俯仰角也存在具有相同的方向余弦的情况，即奇异点问题。

发明内容

为克服现有技术的不足，本发明旨在提出对角化二维DOA估计方法，解决传统的二维DOA估计中，在俯仰角为70°～90°的实际移动通信俯仰角度范围内的角度估计失效问题、俯仰角和方位角的配对问题，以及奇异点问题。为此，本发明采取的技术方案是，扩展孔径二维联合对角化DOA估计方法，利用非均匀L型天线阵列实现，其中在x轴和z轴上分别有两个阵元数目为M的均匀线阵，分别用X，Y，Z，W表示，各子阵中阵元间距为来波信号波长的一半；各坐标轴上的两个子阵阵元间距为d_s，d_s＝hλ/2,h为正整数；具体步骤如下：

(1)构造延时互相关矩阵；

(2)构造选择矩阵，分别得到四个对角矩阵；

(3)得到z轴低精度无模糊的方向余弦估计值以及高精度模糊的方向余弦估计值；

(4)同理得到x轴上对应的两种方向余弦估计值；

(5)通过解模糊方法得到高精度无模糊的方向余弦估计值

(6)得到俯仰角和方位角的估计值

和

当K个窄带非相关信号入射到阵列上，其中第k个信号的二维波达方向为(θ_k,φ_k)，k＝1,2,…K，θ_k和

分别为来波信号的方位角和俯仰角；将位于坐标原点的阵元作为参考阵元，则t时刻的接收数据矢量ρ_ε(t)表示为：

ρ_ε(t)＝A_εs(t)+n_ε(t) (1)

式中，ε＝z,w,x,y,n_ε(t)∈C^M×1是均值为0，方差为σ²的加性高斯白噪声，且与信号s(t)相互独立，

表示阵列流型矩阵；

A_z，A_x分别与子阵Z和子阵X相对应，相应的

具体形式如下式所示：

此外，其他两个子阵的阵列流型矩阵如下

A_y＝A_xΨ(θ)(3)

式中，

步骤(1)、(2)、(3)、(4)具体为：

根据接收数据矢量ρ_x(t)和ρ_w(t)构造互相关矩阵

如下

因此，根据KR运算得到的延时互相关矩阵如下

r_xw(l)＝vec{A_xR_sA_w ^H}＝(A_w ^*⊙A_x)r_s(l)(5)

式中r_xw(l)＝[r_1,1(l),r_2,1(l),...,r_M,1(l),...,r_1,M(l),...,r_M,M(l)]^T，R_s＝diag{r₁(l),r₂(l),...,r_k(l)}，r_s(l)＝(r₁(l),r₂(l),...,r_k(l))^T；

为了充分利用阵列接收信号的空时二维特性，对接收数据矢量ρ_x(t)和ρ_w(t)依据时域最大重叠原则分别划分为L帧数据，第l帧数据表示为：

ρ_x(t)＝[ρ_x(l),ρ_x(l+1),...,ρ_x(l+N-L)]

ρ_w(t)＝[ρ_w(l),ρ_w(l+1),...,ρ_w(l+N-L)] (6)

l＝1,2,...,L，因此，构造延时互相关矩阵R_xw如下

式中

表示延时自相关矩阵R_ss的第k行第l列；

按照同样地方式分别构建延时互相关矩阵R_yw、R_xz和R_yz，在此基础上，定义一个新的矩阵如下

方向余弦的估计具体步骤如下：

通过对R进行奇异值分解SVD(Singular Value Decomposition)，得到信号子空间U_s以及具有K个较大奇异值的对角矩阵Λ_s

从式(8)易知，U_s包含高精度模糊的方向余弦信息以及低精度无模糊的方向余弦信息，构造选择矩阵G₁＝[G₀₁,G₀₀,G₀₂,G₀₀]，G₂＝circshift(G₁，M²)，其中，

因此，包含x轴低精度无模糊的方向余弦对角矩阵表示如下

式中

是一个酉矩阵；

构造选择矩阵G₃＝[G₀₁,G₀₂,G₀₀,G₀₀]，G₄＝circshift(G₃,2M²)，包含z轴低精度无模糊的方向余弦对应对角矩阵表示如下：

阵列中均包含x轴高精度模糊的方向余弦信息，对应对角矩阵表示如下

式中

G₆＝circshift(G₅,1)，

为了得到z轴高精度模糊的方向余弦信息，需要调整U_s的顺序，

G₇＝blkdiag{H₀₁,H₀₁,H₀₁,H₀₁}，H₀₁＝[(circshift(H₀₀,0))^T,...,circshift(H₀₀,M-1))^T]^T，H₀₀＝blkdiag{[1,0,...,0]_1×M,...,[1,0,...,0]_1×M}，H₀₀∈C^M×M；阵列中均包含z轴高精度模糊的方向余弦信息，对应对角矩阵如下

式中

然后，通过联合对角化方法，得到自动配对的

和

因为d＝λ/2，则z轴低精度无模糊的方向余弦估计为

因为d_s>λ/2，方向余弦-1≤υ≤1，则z轴高精度模糊的方向余弦估计值为

式中

表示取不小于l的最小整数,

表示取不大于l的最大整数；

同理，x轴低精度无模糊的方向余弦估计值

以及对应的高精度模糊的方向余弦估计值

符号:(·)^T,(·)^*,(·)^H和

分别表示转置,共轭,共轭转置和伪逆运算，⊙和

分别表示Khatri-Rao积和Kronecker积，E[·]表示统计期望，arg(·)表示相位，I_M是一个维数M×M单位矩阵，diag{·}是由列向量元素组成的对角矩阵，blkdiag{·}表示块对角化，circshift(·,m)是沿着行向右循环移动m个单位。

步骤(5)、(6)具体为：

因为方向余弦估计值一一对应，因此分别估计n_z和n_x，利用解模糊的方法，z轴高精度无模糊的方向余弦估计值为：

其中，

用下式进行估计：

同理，x轴精度无模糊的方向余弦估计值为

第k个信号的方位角和俯仰角估计表达式如下

本发明的特点及有益效果是：

扩展孔径以及KR运算提高了角度估计性能；通过构造联合对角矩阵，能够实现方位角和俯仰角的自动配对，且有效地解决了奇异点问题；在俯仰角为70°～90°的实际移动通信的俯仰角度范围内不会出现角度模糊。

附图说明：

图1非均匀L型天线阵列结构示意图。

图2方位角估计直方图。

图3俯仰角估计直方图。

图4不同角度组合估计联合均方误差。

图5本发明流程图。

具体实施方式

针对已有DOA估计算法存在的问题，本发明提出了一种基于非均匀L型阵列的扩展二维DOA估计算法，该天线阵列为非均匀L型阵列，其中在x轴和z轴上分别有两个阵元数目为M的均匀线阵，分别用X，Y，Z，W表示。各子阵中阵元间距为来波信号波长的一半；各坐标轴上的两个子阵阵元间距为d_s(d_s＝hλ/2,h为正整数)。

本发明采用的技术方案：扩展孔径二维联合对角化DOA估计算法，包括以下步骤：

(1)构造延时互相关矩阵。

(2)构造选择矩阵，分别得到四个对角矩阵。

(3)得到z轴低精度无模糊的方向余弦估计值以及高精度模糊的方向余弦估计值。

(4)同理得到x轴上对应的两种方向余弦估计值。

(5)通过解模糊方法得到高精度无模糊的方向余弦估计值

(6)得到俯仰角和方位角的估计值

和

下面将结合附图和实施例对本发明做进一步的描述：

构造如图1所示的非均匀L型阵列。假设空间中有K个窄带非相关信号入射到阵列上，其中第k个信号的二维波达方向为(θ_k,φ_k)(k＝1,2,…K)，θ_k和

分别为来波信号的方位角和俯仰角。

1延时互相关矩阵的构造

将位于坐标原点的阵元作为参考阵元，则t时刻的接收数据矢量ρ_ε(t)可表示为

ρ_ε(t)＝A_εs(t)+n_ε(t)(1)

式中，ε＝z,w,x,y,n_ε(t)∈C^M×1是均值为0，方差为σ²的加性高斯白噪声，且与信号s(t)相互独立。

表示阵列流型矩阵。

A_z，A_x分别与子阵Z和子阵X相对应。相应的

具体形式如下式所示：

此外，其他两个子阵的阵列流型矩阵如下

A_y＝A_xΨ(θ) (3)

式中，

为了消除高斯白噪声的影响，根据接收数据矢量ρ_x(t)和ρ_w(t)构造互相关矩阵

如下

因此，根据KR运算得到的延时互相关矩阵如下

r_xw(l)＝vec{A_xR_sA_w ^H}＝(A_w ^*⊙A_x)r_s(l) (5)

式中r_xw(l)＝[r_1,1(l),r_2,1(l),...,r_M,1(l),...,r_1,M(l),...,r_M,M(l)]^T，R_s＝diag{r₁(l),r₂(l),...,r_k(l)}，r_s(l)＝(r₁(l),r₂(l),...,r_k(l))^T。

为了充分利用阵列接收信号的空时二维特性，对接收数据矢量ρ_x(t)和ρ_w(t)依据时域最大重叠原则分别划分为L帧数据，第l(l＝1,2,...,L)帧数据可以表示为：

ρ_x(t)＝[ρ_x(l),ρ_x(l+1),...,ρ_x(l+N-L)]

ρ_w(t)＝[ρ_w(l),ρ_w(l+1),...,ρ_w(l+N-L)] (6)

因此，我们可以构造延时互相关矩阵R_xw如下

式中

表示延时自相关矩阵R_ss的第k行第l列。

按照同样地方式分别构建延时互相关矩阵R_yw、R_xz和R_yz。在此基础上，定义一个新的矩阵如下

2方向余弦的估计

通过对R进行奇异值分解(Singular Value Decomposition，SVD)，我们可以得到信号子空间U_s以及具有K个较大奇异值的对角矩阵Λ_s

从式(8)易知，U_s包含高精度模糊的方向余弦信息以及低精度无模糊的方向余弦信息。

构造选择矩阵G₁＝[G₀₁,G₀₀,G₀₂,G₀₀]，G₂＝circshift(G₁，M²)。其中，

因此，包含x轴低精度无模糊的方向余弦对角矩阵可表示如下

式中

是一个酉矩阵。

构造选择矩阵G₃＝[G₀₁,G₀₂,G₀₀,G₀₀]，G₄＝circshift(G₃,2M²)。包含z轴低精度无模糊的方向余弦对应对角矩阵表示如下

阵列中均包含x轴高精度模糊的方向余弦信息，对应对角矩阵表示如下

式中

G₆＝circshift(G₅,1)，

为了得到z轴高精度模糊的方向余弦信息，我们需要调整U_s的顺序，

G₇＝blkdiag{H₀₁,H₀₁,H₀₁,H₀₁}，H₀₁＝[(circshift(H₀₀,0))^T,...,circshift(H₀₀,M-1))^T]^T，H₀₀＝blkdiag{[1,0,...,0]_1×M,...,[1,0,...,0]_1×M}，H₀₀∈C^M×M。

阵列中均包含z轴高精度模糊的方向余弦信息，对应对角矩阵如下

式中

然后，通过联合对角化方法，可得到自动配对的

和

因为d＝λ/2，则z轴低精度无模糊的方向余弦估计为

因为d_s>λ/2，方向余弦-1≤υ≤1，则z轴高精度模糊的方向余弦估计值为

式中

表示取不小于l的最小整数,

表示取不大于l的最大整数。

同理，x轴低精度无模糊的方向余弦估计值

以及对应的高精度模糊的方向余弦估计值

3二维DOA估计的实现

因为方向余弦估计值一一对应，因此分别估计n_z和n_x即可。利用解模糊的方法，z轴高精度无模糊的方向余弦估计值为

其中，

用下式进行估计

同理，x轴精度无模糊的方向余弦估计值为

根据以上的分析，第k个信号的方位角和俯仰角估计表达式如下

符号:(·)^T,(·)^*,(·)^H和

分别表示转置,共轭,共轭转置和伪逆运算。⊙和

分别表示Khatri-Rao(KR)积和Kronecker积。E[·]表示统计期望，arg(·)表示相位。I_M是一个维数M×M单位矩阵。diag{·}是由列向量元素组成的对角矩阵。blkdiag{·}表示块对角化。circshift(·,m)是沿着行向右循环移动m个单位。(特别注意的是，由于MathType软件内置广义逆符号

和“*”有冲突，在一个公式内无法同时显示二者。所以部分矩阵的广义逆用“+”表示。)

结合上述步骤中的实施方式，对本发明的有效性进行仿真验证如下：

仿真中取M＝3，即L型阵列共有11个阵元，阵列间距d＝0.5λ，其中λ为信号波长，进行M＝500次蒙特卡洛仿真。

仿真实验1：假设有K＝2个等功率非相关信号入射到天线阵列，其中SNR＝10dB，快拍数为200数据帧数L为10。信号的方位角和俯仰角为(θ₁,φ₁)＝(45°,65°),(θ₂,φ₂)＝(70°,85°)。图2和图3显示了方位角估计直方图和俯仰角估计直方图。从图中可以看出，本文提出的算法能够清晰的分辨这两个来波信号。

仿真实验2：假设有K＝2个信号入射到天线阵列，两个信号分别为(90°,60°)，(120°,90°)或(65°,33°)，(85°,60°)。其中快拍数N、数据帧数L和信噪比SNR分别为1000，500，30dB。图4为角度估计值分布散点图。

Claims (2)

1.一种扩展孔径二维联合对角化DOA估计方法，其特征是，利用非均匀L型天线阵列实现，其中在x轴和z轴上分别有两个阵元数目为M的均匀线阵，分别用X、Y、Z、W表示，各子阵中阵元间距与来波信号波长成正比，即x和z坐标轴上的两个子阵阵元间距为d_s，d_s＝hλ/2，h为正整数；具体步骤如下：

(1)构造延时互相关矩阵；

(2)构造选择矩阵，分别得到四个对角矩阵；

(3)得到z轴低精度无模糊的方向余弦估计值以及高精度模糊的方向余弦估计值；

(4)同理得到x轴上对应的两种方向余弦估计值；

(5)通过解模糊方法得到z、x轴高精度无模糊的方向余弦估计值

(6)由

得到俯仰角和方位角的估计值

和

当K个窄带非相关信号入射到阵列上，其中第k个信号的二维波达方向为(θ_k，φ_k)，k＝1，2，…K，θ_k和

分别为来波信号的方位角和俯仰角；将位于坐标原点的阵元作为参考阵元，则t时刻的接收数据矢量ρ_ε(t)表示为：

ρ_ε(t)＝A_εs(t)+n_ε(t) (1)

式中，ε＝z、w、x或y，n_ε(t)∈C^M×1是均值为0、方差为σ²的加性高斯白噪声，且与信号s(t)相互独立，

表示阵列流型矩阵；

A_z、A_x分别与子阵Z和子阵X相对应，相应的

具体形式如下式所示：

此外，其他两个子阵的阵列流型矩阵如下：

A_y＝A_xΨ(θ) (3)

式中，

步骤(1)、(2)、(3)、(4)具体为：

根据接收数据矢量ρ_x(t)和ρ_w(t)构造互相关矩阵

如下：

因此，根据KR运算得到的延时互相关矩阵R_xw如下：

式中r_xw(l)＝[r_1，1(l)，r_2，1(l)，...，r_M，1(l)，...，r_1，M(l)，...，r_M，M(l)]^T，R_s＝diag{r₁(l)，r₂(l)，...，r_k(l)}，r_s(l)＝(r₁(l)，r₂(l)，...，r_k(l))^T；

为了充分利用阵列接收信号的空时二维特性，对接收数据矢量ρ_x(t)和ρ_w(t)依据时域最大重叠原则分别划分为L帧数据，第l帧数据表示为：

ρ_x(t)＝[ρ_x(l)，ρ_x(l+1)，...，ρ_x(l+N-L)]

ρ_w(t)＝[ρ_w(l)，ρ_w(l+1)，...，ρ_w(l+N-L)] (6)

l＝1，2，...，L，因此，构造延时互相关矩阵R_xw如下：

式中

表示延时自相关矩阵R_ss的第k行第l列；

按照同样地方式分别构建延时互相关矩阵R_yw、R_xz和R_yz，在此基础上，定义一个新的矩阵如下：

方向余弦的估计具体步骤如下：

通过对R进行奇异值分解SVD，得到信号子空间U_s以及具有K个较大奇异值的对角矩阵Λ_s

式(9)中，U_s包含高精度模糊的方向余弦信息以及低精度无模糊的方向余弦信息，构造选择矩阵G₁＝[G₀₁，G₀₀，G₀₂，G₀₀]，G₂＝circshift(G₁，M²)，其中，

因此，包含x轴低精度无模糊的方向余弦对角矩阵表示如下：

式中

是一个酉矩阵；

构造选择矩阵G₃＝[G₀₁，G₀₂，G₀₀，G₀₀]，G₄＝circshift(G₃，2M²)，包含z轴低精度无模糊的方向余弦对应对角矩阵表示如下：

阵列中均包含x轴高精度模糊的方向余弦信息，对应对角矩阵表示如下：

式中

G₆＝circshift(G₅，1)，

为了得到z轴高精度模糊的方向余弦信息，需要调整U_s的顺序，

G₇＝blkdiag{H₀₁，H₀₁，H₀₁，H₀₁}，H₀₁＝[(circshift(H₀₀，0))^T，...，(circshift(H₀₀，M-1))^T]^T，H₀₀＝blkdiag{[1，0，...，0]_1×M，...，[1，0，...，0]_1×M}，H₀₀∈C^M×M；阵列中均包含z轴高精度模糊的方向余弦信息，对应对角矩阵如下：

式中

然后，通过联合对角化方法，得到自动配对的

和

因为d＝λ/2，则z轴低精度无模糊的方向余弦估计为

因为d_s＞λ/2，方向余弦-1≤υ≤1，则z轴高精度模糊的方向余弦估计值为

式中

表示取不小于l的最小整数，

表示取不大于l的最大整数；

同理得到x轴低精度无模糊的方向余弦估计值

以及对应的高精度模糊的方向余弦估计值

符号：(·)^T、(·)^*、(·)^H和

分别表示转置、共轭、共轭转置和伪逆运算，⊙和

分别表示Khatri-Rao积和Kronecker积，E[·]表示统计期望，arg(·)表示相位，I_M是一个维数M×M单位矩阵，diag{·}是由列向量元素组成的对角矩阵，blkdiag{·}表示块对角化，circshift(·m)是沿着行向右循环移动m个单位。

2.如权利要求1所述的扩展孔径二维联合对角化DOA估计方法，其特征是，步骤(5)、(6)具体为：

因为方向余弦估计值一一对应，因此分别估计n_z和n_x，利用解模糊的方法，z轴高精度无模糊的方向余弦估计值为：

其中，

用下式进行估计：

同理，x轴高精度无模糊的方向余弦估计值为

第k个信号的方位角和俯仰角估计表达式如下：

Images