CN109946643B - 基于music求解的非圆信号波达方向角估计方法 - Google Patents

Abstract

本发明提出了一种基于MUSIC求解的非圆信号波达方向角估计方法，用于解决现有技术中存在的非圆信号波达方向角估计精度较低的技术问题，实现步骤为：1)构建嵌套阵列c；2)获取嵌套阵列c的输出信号Y(t)；3)计算嵌套阵列c的虚拟均匀阵列协方差向量

和椭圆协方差向量

4)构造扩展协方差矩阵G；5)计算扩展协方差矩阵G的噪声子空间U_n；6)构造目标波达方向角θ的阵列流型矩阵Q(θ)；7)基于MUSIC求解的方法，计算目标波达方向角θ的空间谱P(θ)；8)绘制目标波达方向角θ和其空间谱P(θ)的幅度谱图。本发明可用于对飞机、舰船运动目标无源定位。

Description

基于MUSIC求解的非圆信号波达方向角估计方法

技术领域

本发明属于信号处理技术领域，涉及一种非圆信号波达方向角估计方法，特别涉及一种基于MUSIC求解的非圆信号波达方向角估计的方法，可用于对飞机、舰船运动目标无源定位。

背景技术

在现代通信、雷达和卫星系统中，常用的二相相移键控、M进制幅移键控以及幅度调制等调制信号因其非圆信号特性而受到广大专家学者的关注。非圆信号波达方向角是非圆信号的重要特征参数之一，非圆信号波达方向角估计是指利用天线阵列对电磁空间中具有非圆特性的声学信号、电磁信号进行感应接收，再运用现代信号处理方法快速准确的估计出信号源的方向，在雷达、声纳、无线通信等领域具有重要应用价值。非圆信号波达方向角估计精度是评价估计方法的一项重要指标。针对非圆信号波达方向角估计构建的模型中，应用较为广泛的是基于线性均匀阵列的模型，之前的大部分非圆信号波达方向角的估计都是利用该模型生成的。

目前关于非圆信号波达方向角的估计方法是利用多项式运算的方式估计非圆信号波达方向角，这种方法旨在构建波达角选择矩阵的多项式方程，通过求解多项式方程估计非圆信号波达方向角，与传统的非圆信号波达方向角的估计方法相比，大大提高了非圆信号波达方向角的估计精度。例如申请公布号为CN107300686A，名称为“基于多项式求解的非圆信号波达方向角的估计方法”的专利申请中，公开了一种基于多项式求解的非圆信号波达方向角的估计方法。该方法利用嵌套阵列的输出信号构造出两个波达角选择矩阵，经过计算得到四个子噪声矩阵，利用子噪声矩阵构造多项式方程并求解。该方法的估计精度取决于其构造的目标函数方程，利用方程的求解得到估计的角度，在一定程度上提高了估计精度，但是由于其在构造目标函数方程时，将目标函数的最小值直接设为零，而在实际应用中目标函数的最小值不一定为零，因此估计精度往往达不到预期。

发明内容

本发明的目的在于针对上述现有技术存在的不足，提出一种基于MUSIC求解的非圆信号波达方向角估计方法，用于解决现有技术中非圆信号波达方向角的估计精度较低的技术问题。

为实现上述目的，本发明采取的技术方案包括如下步骤：

(1)构建嵌套阵列c：

(1a)以天线接收机作为阵元，构建包括M个阵元的第一均匀线性阵列a和包括N个阵元且相邻阵元的间距为(M+1)d的第二均匀线性阵列b，M≥1，N≥1，a中相邻阵元的间距为d，b中相邻阵元的间距为(M+1)d，0＜d＜λ/2，λ为入射到阵元的窄带信号波长，并定义a中各阵元的位置依次为D(1)＝1，D(2)＝2，…，D(M)＝M，定义b中各阵元的位置依次为D(M+1)＝M+1，D(M+2)＝2(M+1)，…，D(M+N)＝N(M+1)；

(1b)根据第一均匀线性阵列a和第二均匀线性阵列b，构造嵌套阵列c；

(2)获取嵌套阵列c的输出信号Y(t)：

将K个非圆目标信号入射到嵌套阵列c中的每一个阵元，K≥1，并对c的输出信号进行采样，然后对所有采样信号滤波，得到嵌套阵列c的输出信号Y(t)：

Y(t)＝[y₁(t)，…，y_i(t)，…，y_M+N(t)]^T

其中，y_i(t)表示c中第i个阵元的输出信号，1≤i≤M+N，t表示离散时间，1≤t≤L，L表示信号在时域上采样点数，[·]^T表示矩阵[·]的转置运算；

(3)计算嵌套阵列c的虚拟均匀阵列协方差向量

和椭圆协方差向量

(3a)计算等效协方差向量r_d和等效椭圆协方差向量r_s；

(3b)计算虚拟均匀阵列协方差向量

和虚拟均匀阵列椭圆协方差向量

(4)构造扩展协方差矩阵G；

(5)计算扩展协方差矩阵G的噪声子空间U_n：

对扩展协方差矩阵G特征分解，得到G的特征值矩阵Λ和G的特征值所对应的特征向量矩阵U，并对Λ中的特征值按照从大到小的顺序进行排列，然后将(L₁+L₂-K)个较小特征值在U中对应的特征向量作为G的噪声子空间U_n；

(6)构造目标波达方向角θ的阵列流型矩阵Q(θ)：

(6a)将嵌套阵列c[0°，180°]的观测方向范围按0°到180°的递增方向均匀划分成E+1个子观测方向区域，E≥M+N，并将第e个子观测方向区域和第e+1个子观测方向区域之间的边界射线与0°射线之间的夹角作为第e个目标波达方向角θ_e，其中e＝1，2，…，E，定义E个目标波达方向角θ：

θ＝[θ₁，θ₂，…，θ_e，…，θ_E]

(6b)构造目标波达方向角θ的阵列流型向量α₁(θ)和阵列流型向量α₂(θ)，其中θ_e的阵列流型向量α₁(θ_e)和阵列流型向量α₂(θ_e)的表达式分别为：

其中，j为虚数单位；

(6c)构造目标波达方向角θ的阵列流型矩阵Q(θ)，其中θ_e的阵列流型矩阵Q(θ_e)的表达式为：

(7)基于MUSIC求解的方法，计算目标波达方向角θ的空间谱P(θ)，其中θ_e的空间谱P(θ_e)计算公式为：

其中，det[·]表示矩阵[·]的求行列式运算；

(8)绘制目标波达方向角θ与空间谱P(θ)的幅度谱图：

以目标波达方向角θ为x轴坐标，以空间谱P(θ)的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照幅度值从高到低的顺序寻找幅度值较大的前K个谱峰，K个谱峰的峰值点所对应的x轴坐标即为非圆信号波达方向角的估计值。

本发明与现有技术相比，具有以下优点：

本发明对嵌套阵列输出的扩展协方差矩阵进行特征值分解，由于分解得到的噪声子空间和阵列流型矩阵正交的性质，信号入射方向会出现极小值，这样有效地利用矩阵计算将噪声干扰和目标信号分离出来。通过构造空间谱函数确定方向角的估计精度，并绘制空间谱幅度谱图，然后对空间谱幅度谱图进行谱峰搜索，确认波达方向角的估计值。避免了现有技术中，构造目标函数方程时最小值选取对精度的影响，有效地提高了非圆信号波达方向角的估计精度。

附图说明

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

图2是本发明中嵌套阵列的结构示意图；

图3是本发明中E个波达方向角的角度定义图。

具体实施方式

下面结合附图和具体实施例，对本发明作进一步详细描述：

参照图1、一种基于MUSIC求解的非圆信号波达方向角估计方法，包括如下步骤：

步骤1)构建嵌套阵列c：

步骤1a)以天线接收机作为阵元，构建包括M个阵元的第一均匀线性阵列a和包括N个阵元且相邻阵元的间距为(M+1)d的第二均匀线性阵列b，M≥1，N≥1，a中相邻阵元的间距为d，b中相邻阵元的间距为(M+1)d，0＜d＜λ/2，λ为入射到阵元的窄带信号波长，并定义a中各阵元的位置依次为D(1)＝1，D(2)＝2，…，D(M)＝M，定义b中各阵元的位置依次为D(M+1)＝M+1，D(M+2)＝2(M+1)，…，D(M+N)＝N(M+1)；

步骤1b)将b的头部与α的尾部进行拼接，且b中位置为D(M+1)＝M+1的阵元与a中位置为D(M)＝M的阵元之间的距离为d，得到嵌套阵列c，c的结构如图2所示，使用这样的阵列结构模型，解决了一般的线性均匀阵列模型中估计的信号数目低于阵元数目的问题，增加了阵列可识别信源数目，具体实施例中M＝10，N＝10，a中相邻阵元的间距d＝2米。

步骤2)获取嵌套阵列c的输出信号Y(t)：

将K个非圆目标信号入射到嵌套阵列c中的每一个阵元，K≥1，并对c的输出信号进行采样，然后对所有采样信号滤波，得到嵌套阵列c的输出信号Y(t)：

Y(t)＝[y₁(t)，…，y_i(t)，…，y_M+N(t)]^T

其中，y_i(t)表示c中第i个阵元的输出信号，1≤i≤M+N，t表示离散时间，1≤t≤L，L表示信号在时域上采样点数，[·]^T表示矩阵[·]的转置运算；

信号在时域上采样点数L在阵列信号处理中称之为快拍数，理论上来说快拍数越大，目标波达方向角的估计精度就会越高。但考虑到大快拍数会增加阵列系统的软硬件复杂度，提高了系统成本，因此合理选取适应系统的快拍数才能有效地提高目标波达方向角的估计精度，具体实施例中上＝500。在利用MUSIC求解估计目标波达方向角时，需要对入射到嵌套阵列c中的非圆目标信号数量已知，具体实施例中，K＝3；

步骤3)计算嵌套阵列c的虚拟均匀阵列协方差向量

和椭圆协方差向量

步骤3a)根据嵌套阵列c的输出信号Y(t)，计算等效协方差矩阵R_d和等效椭圆协方差矩阵R_s并将R_d中的各列元素按照由小到大的列号进行拼接，得到等效协方差向量r_d，同时将R_s中的各列元素按照由小到大的列号进行拼接，得到椭圆协方差向量r_s，

R_d和R_s计算公式分别为：

其中，[·]^H表示矩阵[·]的共轭转置运算，[·]^T表示矩阵[·]的共轭转置运算；

r_d和r_s表达式分别为：

r_d＝[R_d(1，1)，R_d(2，1)，…，R_d(M+N，1)，R_d(1，2)，…，R_d(M+N，2)，…，R_d(1，M+N)，…，R_d(M+N，M+N)]^T

r_s＝[R_s(1，1)，R_s(2，1)，…，R_s(M+N，1)，R_s(1，2)，…，R_s(M+N，2)，…，R_s(1，M+N)，…，R_s(M+N，M+N)]^T

其中，R_d(i，j)表示协方差矩阵R_d中位于第i行，第j列的元素，1≤i≤M+N，1≤j≤M+N，M+N为嵌套阵列c的阵元数，R_s(i，j)表示椭圆协方差矩阵R_s中位于第i行，第j列的元素；

步骤3b)计算等效协方差向量r_d中所有元素的维数E_i，j和等效椭圆协方差向量r_s中所有元素的维数F_i，j，并删除r_d中维数相同的元素和维数不连续的元素后，将剩余元素按维数从小到大进行排列，得到元素个数为C_d的虚拟均匀阵列协方差向量

同时删除r_s中维数相同的元素和维数不连续的元素后，将剩余元素按维数从小到大进行排列，得到元素个数为C_s的虚拟均匀阵列椭圆协方差向量

其中，C_d＝2MN+2N-1，C_s＝MN+M+N；

E_i，j和F_i，j的计算公式分别为：

E_i，j＝D(j)-D(i)

F_i，j＝D(j)+D(i)

其中，D(i)表示嵌套阵列c中第i个阵元的位置，D(j)表示嵌套阵列c中第j个阵元的位置，1≤i≤M+N，1≤j≤M+N；

步骤4)构造扩展协方差矩阵G：

根据虚拟均匀阵列协方差向量

和椭圆协方差向量

构造扩展协方差矩阵G；

其中，L₁和L₂表示中间变量，L₁＝(C_d+1)/2，L₂＝C_s+1-(C_d+1)/2，且有L₁＞L₂，(·)^*表示向量的共轭运算；

步骤5)计算扩展协方差矩阵G的噪声子空间U_n：

对扩展协方差矩阵G特征分解，得到G的特征值矩阵Λ和G的特征值所对应的特征向量矩阵U，并对Λ中的特征值按照从大到小的顺序进行排列，排列后的特征值矩阵中前K个较大特征值对应的特征向量被称为信号子空间，Λ中剩余的(L₁+L₂-K)个较小特征值对应的特征向量被称为噪声子空间U_n；

扩展协方差矩阵G特征分解，其公式为：

G＝U.Λ.U^H

其中，人为扩展协方差矩阵G的特征值矩阵，U为矩阵G的特征值所对应的特征向量矩阵；

步骤6)构造目标波达方向角θ的阵列流型矩阵Q(θ)：

步骤6a)将嵌套阵列c[0°，180°]的观测方向范围按0°到180°的递增方向均匀划分成E+1个子观测方向区域，E≥M+N，并将第e个子观测方向区域和第e+1个子观测方向区域之间的边界射线与0°射线之间的夹角作为第e个目标波达方向角θ_e，其中e＝1，2，…，E，定义E个目标波达方向角θ：

θ＝[θ₁，θ₂，…，θ_e，…，θ_E]

E个目标波达方向角的角度定义方法如图3所示，在对嵌套阵列c的观测方向范围的进行均匀划分时，子观测方向区域的数量E+1的选择会影响对目标方向角估计的精度，在确定E+1时通常会根据期望达到的角度估计精度进行设定，E+1取值越大对目标方向角估计的精度就越高，具体实施例中E＝180；

步骤6b)构造目标波达方向角θ的阵列流型向量α₁(θ)和阵列流型向量α₂(θ)，其中θ_e的阵列流型向量α₁(θ_e)和阵列流型向量α₂(θ_e)的表达式分别为：

其中，j为虚数单位；

步骤6c)构造目标波达方向角θ的阵列流型矩阵Q(θ)，其中θ_e的阵列流型矩阵Q(θ_e)的表达式为：

步骤7)基于MUSIC求解的方法，计算目标波达方向角θ的空间谱P(θ)，其中θ_e的空间谱P(θ_e)计算公式为：

其中，det[·]表示矩阵[·]的求行列式运算；

MUSIC求解的基本思想是对任意阵列输出信号的协方差矩阵进行特征分解，得到与信号分类相对应的信号子空间和信号分量相正交的噪声子空间后，利用信号子空间和噪声子空间的正交性构造空间谱函数，通过谱峰搜索，对目标信号波达角方向进行估计。本步骤利用步骤5)中得到的扩展协方差矩阵G的噪声子空间U_n和步骤6)中的得到的目标波达方向角θ的阵列流型矩阵Q(θ)，计算出目标波达方向角θ的空间谱P(θ)。为后续的谱峰搜索以及目标信号波达角方向估计值的确定提供了数据基础；

步骤8)绘制目标波达方向角θ与空间谱P(θ)的幅度谱图：

以目标波达方向角θ为横轴坐标，以空间谱P(θ)的幅度值为纵轴坐标，绘制幅度谱图，从该幅度谱图中按照幅度值从高到低的顺序寻找幅度值较大的前K个谱峰，K个谱峰的峰值点所对应的横轴坐标即为非圆信号波达方向角的估计值；

在绘制θ与P(θ)的幅度谱图的过程中，将每个目标波达方向角θ_e的角度值作为横坐标并将它对应的空间谱P(θ_e)的幅度值作为纵坐标，在幅度谱图中以点的形式标出。按这样的方式，绘制E个目标波达方向角θ与空间谱P(θ)的幅度谱图。绘制好的幅度谱图中，总共应含有E个幅度谱点。然后对所有点的进行搜索，寻找出幅度值较大的前K个谱峰，这些谱峰的峰值点所对应的横轴坐标就是本方法MUSIC求解对非圆信号波达方向角的估计值。

Claims (5)

1.一种基于MUSIC求解的非圆信号波达方向角估计方法，其特征在于，包括如下步骤：

(1)构建嵌套阵列c：

(1a)以天线接收机作为阵元，构建包括M个阵元的第一均匀线性阵列a和包括N个阵元且相邻阵元的间距为(M+1)d的第二均匀线性阵列b，M≥1，N≥1，a中相邻阵元的间距为d，b中相邻阵元的间距为(M+1)d，0＜d＜λ/2，λ为入射到阵元的窄带信号波长，并定义a中各阵元的位置依次为D(1)＝1，D(2)＝2，…，D(M)＝M，定义b中各阵元的位置依次为D(M+1)＝M+1，D(M+2)＝2(M+1)，…，D(M+N)＝N(M+1)；

(1b)根据第一均匀线性阵列a和第二均匀线性阵列b，构造嵌套阵列c；

(2)获取嵌套阵列c的输出信号Y(t)：

将K个非圆目标信号入射到嵌套阵列c中的每一个阵元，K≥1，并对c的输出信号进行采样，然后对所有采样信号滤波，得到嵌套阵列c的输出信号Y(t)：

Y(t)＝[y₁(t),…,y_i(t),…,y_M+N(t)]^T

其中，y_i(t)表示c中第i个阵元的输出信号，1≤i≤M+N，t表示离散时间，1≤t≤L，L表示信号在时域上采样点数，[·]^T表示矩阵[·]的转置运算；

(3)计算嵌套阵列c的虚拟均匀阵列协方差向量

和椭圆协方差向量

(3a)计算等效协方差向量r_d和等效椭圆协方差向量r_s；

(3b)计算虚拟均匀阵列协方差向量

和虚拟均匀阵列椭圆协方差向量

(4)构造扩展协方差矩阵G：

根据虚拟均匀阵列协方差向量

和椭圆协方差向量

构造扩展协方差矩阵G：

其中，L₁和L₂表示中间变量，L₁＝(C_d+1)/2,L₂＝C_s+1-(C_d+1)/2，且有L₁＞L₂，(·)^*表示向量的共轭运算；

(5)计算扩展协方差矩阵G的噪声子空间U_n：

对扩展协方差矩阵G特征分解，得到G的特征值矩阵Λ和G的特征值所对应的特征向量矩阵U，并对Λ中的特征值按照从大到小的顺序进行排列，然后将(L₁+L₂-K)个较小特征值在U中对应的特征向量作为G的噪声子空间U_n；

(6)构造目标波达方向角θ的阵列流型矩阵Q(θ)：

(6a)将嵌套阵列c[0°,180°]的观测方向范围按0°到180°的递增方向均匀划分成E+1个子观测方向区域,E≥M+N，并将第e个子观测方向区域和第e+1个子观测方向区域之间的边界射线与0°射线之间的夹角作为第e个目标波达方向角θ_e，其中e＝1,2,...,E，定义E个目标波达方向角θ：

θ＝[θ₁,θ₂,…,θ_e,…θ_E]

(6b)构造目标波达方向角θ的阵列流型向量α₁(θ)和阵列流型向量α₂(θ)，其中θ_e的阵列流型向量α₁(θ_e)和阵列流型向量α₂(θ_e)的表达式分别为：

其中，j为虚数单位；

(6c)构造目标波达方向角θ的阵列流型矩阵Q(θ)，其中θ_e的阵列流型矩阵Q(θ_e)的表达式为：

(7)基于MUSIC求解的方法，计算目标波达方向角θ的空间谱P(θ)，其中θ_e的空间谱P(θ_e)计算公式为：

其中，det[·]表示矩阵[·]的求行列式运算；

(8)绘制目标波达方向角θ与空间谱P(θ)的幅度谱图：

以目标波达方向角θ为x轴坐标，以空间谱P(θ)的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照幅度值从高到低的顺序寻找幅度值较大的前K个谱峰，K个谱峰的峰值点所对应的x轴坐标即为非圆信号波达方向角的估计值。

2.根据权利要求1所述的基于MUSIC求解的非圆信号波达方向角估计方法，其特征在于，步骤(1b)中所述的构造嵌套阵列c，构造方法为：

将b的头部与a的尾部进行拼接，且b中位置为D(M+1)＝M+1的阵元与a中位置为D(M)＝M的阵元之间的距离为d，得到嵌套阵列c；

其中a为第一均匀线性阵列，b为第二均匀线性阵列，M为a中阵元数个阵元，M≥1，D(M+1)＝M+1为b中第一个阵元的位置，D(M)＝M为a中第M个阵元的位置。

3.根据权利要求1所述的基于MUSIC求解的非圆信号波达方向角估计方法，其特征在于，步骤(3a)中所述的计算等效协方差向量r_d和等效椭圆协方差向量r_s,计算方法如下：

根据嵌套阵列c的输出信号Y(t)，计算等效协方差矩阵R_d和等效椭圆协方差矩阵R_s，计算公式分别为：

其中，Y(t)表示嵌套阵列c的输出信号，t表示离散时间，1≤t≤L，L表示信号在时域上采样点数，[·]^H表示矩阵[·]的共轭转置运算，[·]^T表示矩阵[·]的共轭转置运算；

将R_d中的各列元素按照由小到大的列号进行拼接，得到等效协方差向量r_d，同时将R_s中的各列元素按照由小到大的列号进行拼接，得到等效椭圆协方差向量r_s，r_d和r_s的表达式分别为：

r_d＝[R_d(1,1),R_d(2,1),…,R_d(M+N,1),R_d(1,2),…,R_d(M+N,2),…,R_d(1,M+N),…,R_d(M+N,M+N)]^T，

r_s＝[R_s(1,1),R_s(2,1),…,R_s(M+N,1),R_s(1,2),…,R_s(M+N,2),…,R_s(1,M+N),…,R_s(M+N,M+N)]^T，

其中，R_d(i,j)表示协方差矩阵R_d中位于第i行，第j列的元素，1≤i≤M+N，1≤j≤M+N，M+N为嵌套阵列c的阵元数，R_s(i,j)表示椭圆协方差矩阵R_s中位于第i行，第j列的元素，[·]^T表示矩阵[·]的共轭转置运算。

4.根据权利要求1所述的基于MUSIC求解的非圆信号波达方向角估计方法，其特征在于，步骤(3b)中所述的计算虚拟均匀阵列协方差向量

和虚拟均匀阵列椭圆协方差向量

其计算方法如下：

计算等效协方差向量r_d中所有元素的维数E_i,j和等效椭圆协方差向量r_s中所有元素的维数F_i,j，计算公式分别为：

E_i,j＝D(j)-D(i)

F_i,j＝D(j)+D(i)

其中，D(i)表示嵌套阵列c中第i个阵元的位置，D(j)表示嵌套阵列c中第j个阵元的位置，1≤i≤M+N，1≤j≤M+N，M+N为嵌套阵列c的阵元数；

删除r_d中维数相同的元素和维数不连续的元素后，将剩余元素按维数从小到大进行排列，得到元素个数为C_d的虚拟均匀阵列协方差向量

同时删除r_s中维数相同的元素和维数不连续的元素后，将剩余元素按维数从小到大进行排列，得到元素个数为C_s的虚拟均匀阵列椭圆协方差向量

其中，C_d＝2MN+2N-1，C_s＝MN+M+N。

5.根据权利要求1所述的基于MUSIC求解的非圆信号波达方向角估计方法，其特征在于，步骤(5)中所述的扩展协方差矩阵G特征分解，其公式为：

G＝U·Λ·U^H

其中，Λ为扩展协方差矩阵G的特征值矩阵，U为矩阵G的特征值所对应的特征向量矩阵，[·]^H表示矩阵[·]的共轭转置运算。

Images