CN108020812B - 基于特殊三平行线阵结构的二维doa估计方法 - Google Patents

Abstract

本发明涉及阵列信号处理技术，为提出获得自动配对的二维角度值的方法，相比于传统的一些方法，本发明无需对接收数据进行奇异值或者特征值分解，无需进行谱峰搜索步骤，计算复杂度较低，并且可以获得良好的估计性能，具有一定的应用价值。为此，本发明采用的技术方案是，基于特殊三平行线阵结构的二维DOA估计方法，三平行线性阵列结构中的两条子阵列存在一些移位，并依据这移位的特点直接构造出包含角度信息的旋转矩阵，并引入欧拉方程，结合扩展孔径的传播算子算法和均匀线性阵列的旋转不变性，最终获得自动配对的二维角度值。本发明主要应用于阵列信号处理场合。

Description

基于特殊三平行线阵结构的二维DOA估计方法

技术领域

本发明涉及阵列信号处理技术，具体讲,涉及阵列信号处理领域中一种基于特殊三平行线阵的改进型二维DOA估计方法。

背景技术

阵列信号处理的一个基本问题是空间信号到达方向(Direction of Arrival,DOA)估计的问题，也是雷达、声呐、无线通信等众多领域的重要任务之一。DOA估计的基本问题就是确定同时处在空间中某一区域内多个感兴趣的信号源的空间位置(即各个信号到达阵列参考阵元的方向角，简称波达方向)。传统的二维DOA估计方法主要有MUSIC(MultipleSignal Classification,MUSIC)算法、ESPRIT(Estimation of signal parameters viarotational invariance techniques,ESPRIT)、高阶累积量的方法、传播算子的方法等。其中，由于基于传播算子的方法计算复杂度低，方法简单有效，并且可以取得良好的估计效果，受到广泛研究。

发明内容

为克服现有技术的不足，本发明旨在提出获得自动配对的二维角度值的方法，相比于传统的一些方法，本发明无需对接收数据进行奇异值或者特征值分解，无需进行谱峰搜索步骤，计算复杂度较低，并且可以获得良好的估计性能，具有一定的应用价值。为此，本发明采用的技术方案是，基于特殊三平行线阵结构的二维DOA估计方法，三平行线性阵列结构中的两条子阵列存在一些移位，并依据这移位的特点直接构造出包含角度信息的旋转矩阵，并引入欧拉方程，结合扩展孔径的传播算子算法和均匀线性阵列的旋转不变性，最终获得自动配对的二维角度值。

具体地：空间中存在K个远场窄带信号源入射到三维立体阵上，天线阵列被划分为三个子阵列X、Y和Z，其中，每个子阵中阵元之间的间距都是d_x，而三个子阵之间的距离为d_y，X、Y和Z子阵列分别具有M+1、M、M+1个阵元，提出的结构中有X和Z两个子阵列相对于坐标轴有一定的偏移量d_z，θ_k和

为第K个信号源的俯仰角和方位角，则X、Y和Z三个子阵列的输出数据矩阵表示为：

X(t)＝[x₁(t),x₂(t),x₃(t),…,x_M+1(t)]^T＝A_xΦ₁S(t)+n_x(t) (1)

Y(t)＝[y₁(t),y₂(t),y₃(t),…,y_M(t)]^T＝A_yS(t)+n_y(t) (2)

Z(t)＝[z₁(t),z₂(t),z₃(t),…,z_M+1(t)]^T＝A_zΦ₃S(t)+n_z(t) (3)

其中，X(t)、Y(t)、Z(t)指的是子阵列X、Y和Z的接收数据矩阵,t＝1,2,N,N是快拍数，S(t)＝[s₁(t),s₂(t),,s_k(t)]^T，s_k(t)是第K个信号源的矢量矩阵，n_x(t)、n_y(t)、n_z(t)是指噪声向量，三个子阵列的阵列流形矩阵分别为

其中每个信号源的向量表示形式为：

其中，

γ_k＝cosθ_k，那么两个包含角度信息的对角矩阵Φ₁、Φ₃表示成如下：

噪声都是零均值的高斯白噪声，信号源之间是互不相关的，构造另外两个包含角度信息的旋转矩阵Φ₂和Φ₄：

然后，为了更好地进行处理数据，需要对接收数据矩阵做一些变换，利用基于扩展孔径的传播算子方法，将包含角度信息的旋转矩阵求解出来；

最后根据欧拉方程可以从两个包含角度信息的对角矩阵中解出自动配对的方位角和俯仰角的估计值。

一个实例中采用的具体步骤如下：

首先，对接收数据进行处理，考虑到A_x、A_y和A_z之间的关系，得到如下的关系式：

X₁(t)＝A_yΦ₁S(t)+n_x1(t) (11)

X₂(t)＝A_yΦ₂S(t)+n_x2(t) (12)

同理得到：

Z₁(t)＝A_yΦ₃S(t)+n_z1(t) (13)

Z₂(t)＝A_yΦ₄S(t)+n_z2(t) (14)

其中X₁、X₂、Z₁、Z₂为新构造的子阵列的部分数据接收矩阵；

接着构造如下的接收数据矩阵Q(t)：

其中A∈C^3M×K,取矩阵A的前K行构成矩阵A₁，后面的3M-K行构成矩阵A₂，那么接收数据矩阵的协方差R_Q为：

这里划分两个子矩阵为R₁∈C^3M×KR₁∈C^3M×K,R₂∈C^3M×(3M-^K)R₂∈C^3M×(3M-K),那么通过下式求出传播算子P∈C^K×(3M-K)P∈C^K×(3M-K)：

则满足下式：

P^HA₁＝A₂ (18)

定义一个新的矩阵P_e∈C^3M×K为：

其中P₁∈C^M×K、P₂∈C^M×K、P₃∈C^M×K。那么根据式(18)和(19)得到：

通过式(20)得到：

P₂A₁＝A_y(Φ₁+Φ₂)＝P₁A₁(Φ₁+Φ₂) (21)

P₃A₁＝A_y(Φ₃+Φ₄)＝P₁A₁(Φ₃+Φ₄) (22)

定义一个新矩阵Ψ：

式中

是广义逆矩阵。通过特征值分解可以得到一个包含角度信息的对角矩阵，得到

其中Ξ是置换矩阵，同理得到下式一个新矩阵ω：

引入欧拉公式：

通过式(23)、(24)、(25)，可以推导出下式：

这两个矩阵中对应元素的幅值理论上是一致的，但是实际中由于噪声的影响会产生一定的偏差，求其平均值来当做估计值

这样会有利于减小估计的误差，即:

其中，

和

分别是由式(26)和(27)中估计得到的第k个对角元素的幅值，令μ_k和ξ_k分别为(Φ₁+Φ₂)和(Φ₃+Φ₄)对角线上的第k项，那么得到

和θ_k的估计值

和

看出，方位角

和俯仰角θ_k由(29)和(30)式计算得到，其中

由

求得；

假设阵元的间距d＝λ/2，所有信号源有相同的功率

信噪比SNR定义为

表示噪声功率，定义DOA估计的均方根误差公式为：

其中V为蒙特卡洛仿真的次数，

和

分别表示第v次实验中第k个信号源的俯仰角测量值和方位角测量值，取d_x＝d_y＝d_z的条件，通过计算机仿真实验验证该方法的有效性。

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

本发明提供了一种特殊三平行线性阵列中基于改进的扩展孔径传播算子的二维DOA估计方法，可以比较精确的估计出空间中存在的若干个窄带信号发射源的方位角和俯仰角，计算复杂度低，估计误差较低，可以自动实现俯仰角和方位角的配对，可以适用于实时信号处理，在雷达、声呐、无线通信领域等领域的研究中具有非常重要的意义。

附图说明：

图1为提出的特殊三平行天线阵列的结构示意图。

图2为信噪比为15dB时，天线阵列对两个入射信号的分辨图。

图3为在不同快拍数情况下，方位角角度估计误差随信噪比变化趋势图。

图4为在不同快拍数情况下，俯仰角角度估计误差随信噪比变化趋势图。

图5为在不同阵元总数下，方位角角度估计误差随信噪比变化趋势图。

图6为在不同阵元总数下，俯仰角角度估计误差随信噪比变化趋势图。

图7为二维联合角度估计误差随信噪比变化趋势图。

图8为一个入射信号时在不同方位角和俯仰角组合时的估计误差mesh图。

具体实施方式

提出一种特殊的三平行线性阵列，与传统的三平行线性阵列不同，此阵列结构中的两条子阵列存在一些移位，并依据这移位的特点直接构造出包含角度信息的旋转矩阵。并引入欧拉方程，结合改进的扩展孔径的传播算子算法和均匀线性阵列的旋转不变性，最终可以获得自动配对的二维角度值。相比于传统的一些方法，该方法无需对接收数据进行奇异值或者特征值分解，无需进行谱峰搜索步骤，计算复杂度较低，并且可以获得良好的估计性能，具有一定的应用价值。

本发明利用三平行线性天线阵列上不同阵元接收的数据信息估计空间中同时存在的K个远场窄带信号源，获得每个发射信号源相对于原点坐标的方位角和俯仰角信息。

本发明的具体实现方案如下：首先建立系统的数学分析模型，如图1所示，假设空间中存在K个远场窄带信号源入射到三维立体阵上，天线阵列被划分为三个子阵列X、Y和Z。其中，每个子阵中阵元之间的间距都是d_x，而三个子阵之间的距离为d_y。不同于传统的平行阵列结构，提出的结构中有X和Z两个子阵列相对于坐标轴有一定的偏移量d_z。假定θ_k和

为第K个信号源的俯仰角和方位角，则阵列的输出数据可以表示为：

X(t)＝[x₁(t),x₂(t),x₃(t),…,x_M+1(t)]^T＝A_xΦ₁S(t)+n_x(t) (1)

Y(t)＝[y₁(t),y₂(t),y₃(t),…,y_M(t)]^T＝A_yS(t)+n_y(t) (2)

Z(t)＝[z₁(t),z₂(t),z₃(t),…,z_M+1(t)]^T＝A_zΦ₃S(t)+n_z(t) (3)

其中，X(t)、Y(t)、Z(t)指的是子阵列X、Y和Z的接收数据矩阵，t＝1,2,N,N是快拍数，S(t)＝[s₁(t),s₂(t),,s_k(t)]^T，n_x(t)、n_y(t)、n_z(t)是指噪声向量，s_k(t)是第K个信号源。阵列流形

其中每个信号源的向量表示形式为：

其中，

γ_k＝cosθ_k，包含角度信息的对角矩阵如下：

假设噪声都是零均值的高斯白噪声，信号源之间是互不相关的。构造另外两个包含角度信息的旋转矩阵：

然后，为了更好地进行处理数据，需要对接收数据矩阵做一些变换，利用基于扩展孔径的传播算子方法，可以将包含角度信息的旋转矩阵求解出来。

最后根据欧拉方程可以从两个包含角度信息的对角矩阵中解出自动配对的方位角和俯仰角的估计值。

下面是本发明的一个实例，采用的具体步骤如下：

首先，对接收数据进行处理，考虑到A_x、A_y和A_z之间的关系，可以得到如下的关系式：

X₁(t)＝A_yΦ₁S(t)+n_x1(t) (11)

X₂(t)＝A_yΦ₂S(t)+n_x2(t) (12)

同理可以得到：

Z₁(t)＝A_yΦ₃S(t)+n_z1(t) (13)

Z₂(t)＝A_yΦ₄S(t)+n_z2(t) (14)

其中X₁、X₂、Z₁、Z₂为新构造的子阵列的部分数据接收矩阵。

接着构造如下的接收数据矩阵Q(t)：

其中A∈C^3M×K,取矩阵A的前K行构成矩阵A₁，后面的3M-K行构成矩阵A₂。那么接收数据矩阵的协方差R_Q为：

这里划分两个子矩阵为R₁∈C^3M×K,R₂∈C^3M×(3M-K),那么通过下式即可求出传播算子P∈C^K×(3M-K)：

则满足下式：

P^HA₁＝A₂(18)

定义一个新的矩阵P_e∈C^3M×K为：

其中P₁∈C^M×K、P₂∈C^M×K、P₃∈C^M×K。那么根据式(18)和(19)可以得到：

通过式(20)我们可以得到：

P₂A₁＝A_y(Φ₁+Φ₂)＝P₁A₁(Φ₁+Φ₂) (21)

P₃A₁＝A_y(Φ₃+Φ₄)＝P₁A₁(Φ₃+Φ₄) (22)

定义一个新的矩阵Ψ：

式中

是广义逆矩阵。通过特征值分解可以得到一个包含角度信息的对角矩阵，我们可以得到

其中Ξ是置换矩阵，同理下式一个新矩阵ω：

引入欧拉公式：

通过式(23)、(24)、(25)，可以推导出下式：

可以看到这两个对角矩阵中包含了方位角和俯仰角的角度信息，由于(24)式的计算方式，所以式(26)和(27)中的矩阵中同一位置上的元素是配对的。这两个矩阵中对应元素的幅值理论上是一致的，但是实际中由于噪声的影响会产生一定的偏差，可以求其平均值来当做估计值

这样会有利于减小估计的误差，即:

其中，

和

分别是由式(26)和(27)中估计得到的第k个对角元素的幅值。令μ_k和ξ_k分别为(Φ₁+Φ₂)和(Φ₃+Φ₄)对角线上的第k项，那么可以得到

和θ_k的估计值

和

可以看出，方位角

和俯仰角

可以由(29)和(30)式计算得到，其中

可以由

求得。

不失一般性，我们假设阵元的间距d＝λ/2，所有信号源有相同的功率

信噪比SNR定义为

表示噪声功率。定义DOA估计的均方根误差公式为：

其中V为蒙特卡洛仿真的次数，

和

分别表示第v次实验中第k个信号源的俯仰角测量值和方位角测量值，实验中我们取d_x＝d_y＝d_z的条件，通过计算机仿真实验可以验证该方法的有效性。

Claims (1)

1.一种基于特殊三平行线阵结构的二维DOA估计方法，其特征是，三平行线性阵列结构中的两条子阵列存在一些移位，并依据这移位的特点直接构造出包含角度信息的旋转矩阵，并引入欧拉方程，结合扩展孔径的传播算子算法和均匀线性阵列的旋转不变性，最终获得自动配对的二维角度值；具体步骤如下：

空间中存在K个远场窄带信号源入射到三维立体阵上，天线阵列被划分为三个子阵列X、Y和Z，其中，每个子阵中阵元之间的间距都是d_x，而三个子阵之间的距离为d_y，X、Y和Z子阵列分别具有M+1、M、M+1个阵元，提出的结构中有X和Z两个子阵列相对于坐标轴有一定的偏移量d_z，θ_k和

为第K个信号源的俯仰角和方位角，则X、Y和Z三个子阵列的输出数据矩阵表示为：

X(t)＝[x₁(t),x₂(t),x₃(t),…,x_M+1(t)]^T＝A_xΦ₁S(t)+n_x(t) (1)

Y(t)＝[y₁(t),y₂(t),y₃(t),…,y_M(t)]^T＝A_yS(t)+n_y(t) (2)

Z(t)＝[z₁(t),z₂(t),z₃(t),…,z_M+1(t)]^T＝A_zΦ₃S(t)+n_z(t) (3)

其中，X(t)、Y(t)、Z(t)指的是子阵列X、Y和Z的接收数据矩阵,t＝1,2,…N,N是快拍数，S(t)＝[s₁(t),s₂(t),…,s_k(t)]^T，s_k(t)是第K个信号源的矢量矩阵，n_x(t)、n_y(t)、n_z(t)是指噪声向量，三个子阵列的阵列流形矩阵分别为

其中每个信号源的向量表示形式为：

其中，

γ_k＝cosθ_k，那么两个包含角度信息的对角矩阵Φ₁、Φ₃表示成如下：

噪声都是零均值的高斯白噪声，信号源之间是互不相关的，构造另外两个包含角度信息的旋转矩阵Φ₂和Φ₄：

然后，为了更好地进行处理数据，需要对接收数据矩阵做一些变换，利用基于扩展孔径的传播算子方法，将包含角度信息的旋转矩阵求解出来；

最后根据欧拉方程可以从两个包含角度信息的对角矩阵中解出自动配对的方位角和俯仰角的估计值，具体步骤如下：

首先，对接收数据进行处理，考虑到A_x、A_y和A_z之间的关系，得到如下的关系式：

X₁(t)＝A_yΦ₁S(t)+n_x1(t) (11)

X₂(t)＝A_yΦ₂S(t)+n_x2(t) (12)

同理得到：

Z₁(t)＝A_yΦ₃S(t)+n_z1(t) (13)

Z₂(t)＝A_yΦ₄S(t)+n_z2(t) (14)

其中X₁、X₂、Z₁、Z₂为新构造的子阵列的部分数据接收矩阵；

接着构造如下的接收数据矩阵Q(t)：

其中A∈C^3M×K,取矩阵A的前K行构成矩阵A₁，后面的3M-K行构成矩阵A₂，那么接收数据矩阵的协方差R_Q为：

这里划分两个子矩阵为R₁∈C^3M×KR₁∈C^3M×K,R₂∈C^3M×(3M-K)R₂∈C^3M×(3M-K),那么通过下式求出传播算子P∈C^K×(3M-K)P∈C^K×(3M-K)：

则满足下式：

P^HA₁＝A₂ (18)

定义一个新的矩阵P_e∈C^3M×K为：

其中P₁∈C^M×K、P₂∈C^M×K、P₃∈C^M×K，那么根据式(18)和(19)得到：

通过式(20)得到：

P₂A₁＝A_y(Φ₁+Φ₂)＝P₁A₁(Φ₁+Φ₂) (21)

P₃A₁＝A_y(Φ₃+Φ₄)＝P₁A₁(Φ₃+Φ₄) (22)

定义一个新矩阵Ψ：

式中

是广义逆矩阵，通过特征值分解可以得到一个包含角度信息的对角矩阵，得到A₁＝A₁Ξ,其中Ξ是置换矩阵，同理得到下式一个新矩阵ω：

引入欧拉公式：

通过式(23)、(24)、(25)，可以推导出下式：

这两个矩阵中对应元素的幅值理论上是一致的，但是实际中由于噪声的影响会产生一定的偏差，求其平均值来当做估计值w_k，这样会有利于减小估计的误差，即:

其中，w_k1和w_k2分别是由式(26)和(27)中估计得到的第k个对角元素的幅值，令μ_k和ξ_k分别为(Φ₁+Φ₂)和(Φ₃+Φ₄)对角线上的第k项，那么得到

和θ_k的估计值

和θ_k：

看出，方位角

和俯仰角θ_k由(29)和(30)式计算得到，其中

由

求得；

假设阵元的间距d＝λ/2，所有信号源有相同的功率

信噪比SNR定义为

表示噪声功率，定义DOA估计的均方根误差公式为：

其中V为蒙特卡洛仿真的次数，

和

分别表示第v次实验中第k个信号源的俯仰角测量值和方位角测量值。

Images