CN108398659B - 一种矩阵束与求根music结合的波达方向估计方法 - Google Patents

Abstract

本发明公开了一种矩阵束与求根MUSIC结合的波达方向估计方法，具体步骤包括：S1、构造分块Hankel矩阵，并进一步得到自相关数据矩阵，通过奇异值分解获取噪声子空间与信号子空间，构造矩阵束方程；S2、通过矩阵束方程求解第一维方向余弦，并通过实值映射判断重复特征值；S3、实值特征值映射回复值空间，对重复特征值参数取算术平均值，构造一维导向矢量；S4、通过正交方程构造多项式，通过多项式解出第二维方向余弦参数。本发明能够较好地解决重复特征值带来的多角度模糊问题，完成空间方位角与仰角的自动匹配，避免了重复的矩阵束特征值求解及额外的配对方法，一定程度上降低了方法的复杂度，同时可根据实际需求进行方法调整，改变空间参数的估计顺序。

Description

一种矩阵束与求根MUSIC结合的波达方向估计方法

技术领域

本发明属于雷达及阵列信号处理领域，特别涉及一种矩阵束与求根MUSIC结合的波达方向估计方法。

背景技术

波达方向(DOA)估计是指对空间信号到达天线阵元的方位角等参数的估计。在移动通信、无线电导航定位、雷达探测、声纳、医学诊断等许多场景中，空间信源的二维波达方向(2D-DOA)估计都有着十分重要的应用，特别在军事领域，如双基地雷达中，准确定位有用信息方位，抑制干扰信息等方面，是一个研究热点。

目前已有的DOA估计方法包括波束成形类算法和子空间算法，波束成形类算法的应用十分广泛，但是由于收到瑞利极限的限制，其空间角度分辨率较低，而自1979年Schmidt提出MUSIC算法后，基于子空间的DOA估计方法开始引起广泛的关注，该类方法以信号子空间和噪声子空间的正交性为理论基础，通过空间划分来进行空间信号参数估计，其中典型代表有多重信号分类(MUSIC)、旋转不变子空间(ESPRIT)与矩阵束(MP)等。

发明内容

本发明的主要目的在于克服现有技术的缺点与不足，提供一种矩阵束与求根MUSIC结合的波达方向估计方法。本发明通过矩阵束方法求解得到某维方位信息，并构造阵列导向矢量，利用导向矢量与噪声子空间的正交性构造多项式，通过对多项式求根得到另一维方位信息，完成空间方位角与仰角的自动配对。

建立天线阵元接收模型。

空间信源通过无线信道的传输情况较为复杂，为得到较为实用的参数模型，本发明简化了波形传输和信号类型。假设有P个波长为λ的远场窄带非相关信源入射到均匀分布在X-Y平面的均匀矩形阵列上，假设该矩形阵列在X与Y方向上的阵元数目分别为M和N，第P个信源的仰角与方位角分别为

和θ_p，则t时刻位于第(m,n)个阵元上的无噪接收信号为：

上式中s_p(t)为入射到阵列的第P个信源，x_p与y_p为X方向与Y方向上的阵列流形，其定义分别为

和

其中，d_x和d_y分别为X方向与Y方向上相邻两个阵元之间的间隔，u_x和u_y是两个一维的方向余弦，分别为

cosθ_p和

sinθ_p，由此可以定义二维空间频率

由上式可以看出，对于二维DOA估计，只要分别求出二维空间频率并正确配对，就能正确估计出空间信源的来波方向。

一种矩阵束与求根MUSIC结合的波达方向估计方法，具体步骤包括：

S1、由天线阵列每个阵元接收的叠加信号构造分块Hankel矩阵，并进一步得到自相关数据矩阵，通过奇异值分解获取噪声子空间与信号子空间，构造矩阵束方程；

S2、通过矩阵束方程求解第一维方向余弦，并通过实值映射判断重复特征值；

S3、对已分组的实值空间特征值进行映射得到复值空间特征值，对重复特征值参数取算术平均值，构造一维导向矢量；

S4、根据噪声子空间与阵列导向矢量的正交性，构造多项式，并通过多项式求根解出第二维方向余弦参数，由角度公式求解方位角与仰角。

进一步地，所述步骤S1，具体包括：

S11、根据矩阵序列{D₀(t),D₁(t)...,D_M-1(t)}和无噪声信号序列{z_(m,0)(t),z_(m,1)(t),...,z_(m,N-1)(t)}进行处理，构造分块Hankel矩阵；

S12、根据阵元接收模型，对Hankel矩阵进行处理，得到自相关数据矩阵；

S13、对自相关矩阵进行奇异值分解，获取信号子空间和噪声子空间；

S14、构造矩阵束方程。

进一步地，所述分块Hankel矩阵，具体为：

其中，K表示x方向上的束参数，M表示x方向上阵元个数，D_m(t)表示一维Hankel矩阵，具体为：

其中，Q表示y方向上的束参数，N表示y方向上阵元个数，z_(m,n)(t)表示t时刻位于第(m,n)个阵元上的无噪接收信号，即阵元接收模型表达式，具体为：

进一步地，所述步骤S12中，根据阵元接收模型表达式，推导得到

其中，S(t)表示对角元素为P个信源表达式的对角矩阵，具体为：

S(t)＝diag(s₁(t),s₂(t),...,s_P(t))

X_d表示一维方向余弦，具体为：

X_d＝diag(x₁,x₂,...,x_P)

Q_L、Q_R为包含第二维参数信息的两个范德蒙矩阵，具体为：

其中，矩阵中元素表示不同阵元产生的相位延迟。

进一步地，根据D_m(t)的表达式以及D_e(t)与D_m(t)的关系，得到数据矩阵D_e(t)新的表示形式为：

D_e(t)＝E_LS(t)E_R

其中

为了获得更高的精度，根据上述数据矩阵D_e(t)，构造中心厄米特数据矩阵D_ce(t)，表示方式为：

根据中心厄米特数据矩阵，进一步得到自相关数据矩阵R_ss，表示方式为：

对R_ss进行奇异值分解，获取信号子空间和噪声子空间，具体分解方式为：

其中，min表示取值数据矩阵行列数中最小者，U_s与V_s分别为信号子空间的左右奇异矢量，同理U_n与V_n表示噪声子空间的左右奇异矢量。

表示对V_s进行共轭转置操作。Σ_s与Σ_n为对角矩阵，其对角线元素分别表示信号子空间与噪声子空间的奇异值；由于经过奇异值分解后，奇异值以从大到小方式进行排序，在无噪声情况下，前P个奇异值大于零，且对应的奇异矢量构成信号子空间，剩余奇异值均为零，且对应的奇异矢量构成噪声子空间。

进一步地，为了保证E_L、E_R及R_ss为满秩矩阵，则束参数K、Q选取时必须满足以下约束条件：

(M-K+1)(N-Q+1)≥P，K≥P且Q≥P

确保正确估计空间中的P个非相干信源。

进一步地，所述步骤S14中，U_s和阵列导向矢量E_L关系式如下：

U_s＝E_LT

其中，T为非奇异矩阵。

根据U_s的结构及其与E_L关系式，信号子空间分别去掉后Q行和前Q行得到U_s的两个子矩阵U_s1和U_s2。

由U_s＝E_LT可知，U_s1＝E_L1T，U_s2＝E_L2T,同样可知E_L1和E_L2分别为E_L去掉其后Q行和前Q行所得到，且有E_L2＝E_L1X_d。因此，构造矩阵束方程，具体为：

U_s2-λU_s1＝E_L1(X_d-λI)T

得到E_L2＝E_L1T^-1X_dT。定义矩阵Ψ_x＝T^-1X_dT，因此，Ψ_x与X_d为相似矩阵。

进一步地，所述步骤S2，具体包括：

S21、通过矩阵束方程求解空间中的P个信源，得到第一维的P个方向余弦，利用实值映射方法得到实值矢量空间；

S22、对实值矢量空间进行重复特征值判断，得到若干组重复特征值组。

进一步地，所述步骤S21中，通过矩阵束方程求解得到的第一维P个方向余弦，记为

为了判断P个方向余弦之间是否重复，而复值空间的特征值判定步骤较为繁琐，因此，将复值空间的特征值映射到实值空间进行判断。

进一步地，所述将复值空间的特征值映射到实值空间，首先将数据矩阵D_ce(t)进行酉变换，具体为：

得到

其中，

与

为酉变换矩阵。

酉变换后得到的实值矩阵信号子空间为

由于U_s与二维阵列导向矢量E_L张成同一空间，因此，推导出实值阵列导向矢量与复值阵列导向矢量的关系为

进一步地，所述步骤S21中，根据实值阵列导向矢量与复值阵列导向矢量的关系式，对步骤S14中的矩阵束方程U_s2-λU_s1＝E_L1(X_d-λI)T进行改写，得到

D₁E_L＝D₂E_LX_d

其中，

因此，有

且

进一步地，根据改写后矩阵束方程的上述两个等式，分别构造矩阵W₁与W₂，具体为：

因为有

因此得到

(W₁+jW₂)E_L＝(W₁-jW₂)E_LX_d (1)

对公式(1)改写得到：

W₁E_L(X_d-I)＝jW₂E_L(X_d+I) (2)

进一步得到：

W₁E_Lj(I+X_d)(I-X_d)^-1＝W₂E_L (3)

定义矩阵Y＝j(I+X_d)(I-X_d)^-1，显然

进一步地，根据矩阵Y与X_d的关系式，将复值空间特征值一一映射到实值空间特征值来进行重复特征值判断，具体映射关系式为：

其中，r_p为对应于c_p的第p个实值空间一维特征值。

通过映射关系式可以将复值空间映射到实值空间，从而简化了重复特征值判定的步骤，对c中元素逐个进行实值空间映射，得到了实值矢量空间r＝{r₁,r₂,...,r_P}。

进一步地，在步骤S22中，由于环境噪声的存在，得到的实值矢量空间r中的每个最小二乘估计值并非实值，因此需要对r中的每个元素取实部，并对该实序列按升序排列，得到一组新的序列：

r'＝{Re(r_i),Re(r_j),...,Re(r_k)},(i,j,k＝1,2,...,P,i≠j≠k)

记序列中元素为r₁',r₂',...,r_P'，定义如下矩阵Q：

将矢量r'与Q相乘得到序列S,S＝r'Q＝[r₁'-r₂',r₂'-r₃',...,r_P-1'-r_P']。

进一步地，对得到的序列S进行重复特征值判断，具体为：

记S中元素为s₁,s₂,...,s_P,设定误差精度β＝0.05，对S中元素从首位开始判断，若元素S_j满足|S_j|＜β，则继续判断下一个元素，直到元素s_k不满足|S_k|＜β为止，则认为r中元素r_j',...,r_k'为重复特征值，并将该k-j+1个元素归为一组。

对S中所有元素进行上述重复特征值判断，直到所有特征值分组完毕，则可得到若干组重复特征值组：

{r_i',r_j',....,r_k'}，{r_q',r_r',....,r_s'},....(i,j,k,q,r,s＝1,2,...,P,i≠j≠p≠q≠r≠s)

进一步地，所述步骤S3中，对于第二维参数，由于利用了求根MUSIC方法构造多项式进行计算，因此，具体包括：

S31、对分组后的实值空间特征值映射得到复值空间特征值；

S32、对每组重复特征值取算术平均值，构造一维阵列导向矢量。

进一步地，所述步骤S31中，对步骤S22中公式(4)进行变形，得到：

其中，c_p'为对应于r_p的第p个实值空间一维特征值；将步骤S22中得到的若干组重复特征值组代入上述映射公式，得到分组后的复值空间特征值组为：

{c_i',c_j',....,c_k'},{c_q',c_r',....,c_s'},.....

分别记为K₁,K₂,...,且重数分别为J、L....。

进一步地，所述步骤S32中，具体步骤为：

(1)对每组参数求取算数平均：

(2)对参数K₁、K₂...分别构造一维阵列导向矢量：

a(μ_x1)＝(K₁ ⁰,K₁ ¹,...,K₁ ^K-1)，a(μ_x2)＝(K₂ ⁰,K₂ ¹,...,K₂ ^K-1)......

进一步地，所述步骤S4中，根据噪声子空间与阵列导向矢量的正交性，得到求根多项式：

其中

将一维阵列导向矢量a(μ_x1)、a(μ_x2)......分别代入上述求根多项式中求根。

(1)在无噪声的情况下，对于J重的特征值a(μ_x1)，有J个根分布在单位圆上，且这J个根表示感兴趣的J个信源的二维方向余弦信息。

(2)在有噪声的情况下，对于J重的特征值a(μ_x1)，则有2J个分布在单位圆附近的根，包含了该J个信源的第二维参数信息，可选取单位圆内最靠近单位的J个根，并根据方向余弦与角度关系公式：

将J重一维方向余弦分别与所求出的J个根代入上述角度计算公式中，求得J个信源的方位角与仰角信息。

将余下的导向矢量a(μ_x2)，a(μ_x3).....通过多项式求得的根代入上述角度计算公式中，直至最后一个参数计算完毕。

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明相较于传统矩阵束方法以及二维root-music方法，一定程度上降低了复杂度，并解决了传统的配对方法由于重复特征值导致的角度模糊问题；

2、在数据接收矩阵结构上，本发明利用了接收信源的共轭信息，使得低信噪比下角度估计信息更加精确；

3、在方法设计上，本发明结合了二维矩阵束和求根MUSIC方法，避免了重复特征值计算和配对方法。

附图说明

图1为本发明的步骤流程示意图；

图2为本发明中分块Hankel矩阵的构造示意图；

图3为本发明中天线阵列示意图；

图4为本发明中方位角与仰角为(30°,60°,60°)的空间信源经由多项式求根得到的单根分布图；

图5为本发明中方位角与仰角为(30°,60°,60°)的空间信源经由多项式求根得到的重根分布图；

图6为本发明中方位角与仰角为(30°,60°,30°)的空间信源经由多项式求根得到的重根分布图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例

如图3所示，本实施例中天线阵列是由M×N个等距均匀分布阵元组成，基于该天线阵列，本发明提供了一种矩阵束与求根MUSIC结合的波达方向估计方法，如图1所示，该波达方向估计方法包括下述步骤：

S1、由天线阵列每个阵元接收的叠加信号构造分块Hankel矩阵(如图2所示)，并进一步得到自相关数据矩阵，通过奇异值分解获取噪声子空间与信号子空间，构造矩阵束方程；

进一步地，所述步骤S1，具体包括：

S11、根据矩阵序列{D₀(t),D₁(t)...,D_M-1(t)}和无噪声信号序列{z_(m,0)(t),z_(m,1)(t),...,z_(m,N-1)(t)}进行处理，构造分块Hankel矩阵；如图2所示为本发明中分块Hankel矩阵构造示意图。

所述分块Hankel矩阵，具体为：

其中，K表示x方向上的束参数，M表示x方向上的阵元个数，D_m(t)表示一维Hankel数据矩阵，具体为：

其中，Q表示y方向上的束参数，N表示y方向上阵元个数，z_(m,n)(t)表示t时刻位于第(m,n)个阵元上的无噪接收信号，即阵元接收模型表达式，具体为：

其中，上式中s_p(t)为入射到阵列的第p个信源，x_p与y_p为X方向与Y方向上的阵列流型，其定义分别为

和

其中，d_x和d_y分别为X方向与Y方向上相邻两个阵元之间的间隔，u_x和u_y是两个一维的方向余弦，分别为

cosθ_p和

sinθ_p。

S12、根据阵元接收模型，对Hankel矩阵进行处理，得到自相关数据矩阵；

所述步骤S12中，根据阵元接收模型表达式，推导得到

其中，S(t)表示对角元素为P个信源表达式的对角矩阵，具体为：

S(t)＝diag(s₁(t),s₂(t),...,s_P(t))

X_d表示一维方向余弦，具体为：

X_d＝diag(x₁,x₂,...,x_P)

Q_L、Q_R为两个包含第二维参数信息的范德蒙矩阵，具体为：

根据D_m(t)的表达式以及D_e(t)与D_m(t)的关系，得到数据矩阵D_e(t)新的表示形式为：

D_e(t)＝E_LS(t)E_R

其中

为了获得更高的精度，根据上述数据矩阵D_e(t)，构造中心厄米特数据矩阵D_ce(t)，表示方式为：

进一步得到自相关数据矩阵R_ss，表示方式为：

S13、对自相关矩阵进行奇异值分解，获取信号子空间和噪声子空间；

对R_ss进行奇异值分解，获取信号子空间和噪声子空间，具体分解方式为：

其中，min表示取值数据矩阵行列数中最小者，U_s与V_s分别为信号子空间的左右奇异矢量，同理U_n与V_n表示噪声子空间的左右奇异矢量。

表示对V_s进行共轭转置操作。Σ_s与Σ_n为对角矩阵，其对角线元素分别表示信号子空间与噪声子空间的奇异值。由于经过奇异值分解后，奇异值以从大到小方式进行排序，在无噪声情况下，前P个奇异值大于零，且对应的奇异矢量构成信号子空间，剩余奇异值均为零，且对应的奇异矢量构成噪声子空间。

进一步地，为了保证E_L、E_R及R_ss为满秩矩阵，则束参数K、Q选取时必须满足以下约束条件：

(M-K+1)(N-Q+1)≥P，K≥P且Q≥P

确保正确估计空间中的P个非相干信源。

S14、构造矩阵束方程；

所述步骤S14中，U_s和阵列导向矢量E_L关系式如下：

U_s＝E_LT

其中，T为非奇异矩阵。

根据U_s的结构及其与E_L关系式，信号子空间分别去掉后Q行和前Q行得到U_s的两个子矩阵U_s1和U_s2。

由U_s＝E_LT可知，U_s1＝E_L1T，U_s2＝E_L2T,同样可知E_L1和E_L2分别为E_L去掉其后Q行和前Q行所得到，且有E_L2＝E_L1X_d。因此，构造矩阵束方程，具体为：

U_s2-λU_s1＝E_L1(X_d-λI)T

得到E_L2＝E_L1T^-1X_dT。定义矩阵Ψ_x＝T^-1X_dT，因此，Ψ_x与X_d为相似矩阵。

S2、通过矩阵束方程求解第一维方向余弦，并通过实值映射判断重复特征值；

进一步地，所述步骤S2，具体包括：

S21、通过矩阵束方程求解空间中的P个信源，得到第一维的P个方向余弦，利用实值映射方法得到实值矢量空间；

所述步骤S21中，通过矩阵束方程求解得到的第一维P个方向余弦，记为c＝{c₁,c₂,...,c_P}。

为了判断P个方向余弦之间是否重复，而复值空间的特征值判定步骤较为繁琐，因此，将复值空间的特征值映射到实值空间进行判断。

进一步地，所述将复值空间的特征值映射到实值空间，首先将数据矩阵D_ce(t)进行酉变换，具体为：

得到

其中，

与

为酉变换矩阵。

酉变换后得到的实值矩阵信号子空间为

由于U_s与二维阵列导向矢量E_L张成同一空间，因此，推导出实值阵列导向矢量与复值阵列导向矢量的关系为

所述步骤S21中，根据实值阵列导向矢量与复值阵列导向矢量的关系式，对步骤S14中的矩阵束方程U_s2-λU_s1＝E_L1(X_d-λI)T进行改写，得到

D₁E_L＝D₂E_LX_d

其中，

因此，有

且

根据改写后矩阵束方程的上述两个等式，分别构造矩阵W₁与W₂，具体为：

因为有

因此得到

(W₁+jW₂)E_L＝(W₁-jW₂)E_LX_d

对上式改写得到：

W₁E_L(X_d-I)＝jW₂E_L(X_d+I)

进一步得到：

W₁E_Lj(I+X_d)(I-X_d)^-1＝W₂E_L

定义矩阵Y＝j(I+X_d)(I-X_d)^-1，显然

根据矩阵Y与X_d的这一关系，将复值空间特征值一一映射到实值空间特征值来进行重复特征值判断，具体映射关系式为：

其中，r_p为对应于c_p的第p个实值空间一维特征值。

通过映射关系式可以将复值空间映射到实值空间，从而简化了重复特征值判定的步骤，对c中元素逐个进行实值空间映射，得到了实值矢量空间r＝{r₁,r₂,...,r_P}。

S22、对实值矢量空间进行重复特征值判断，得到若干组重复特征值组；

在步骤S22中，由于环境噪声的存在，得到的实值矢量空间r中的每个最小二乘估计值并非实值，因此需要对r中的每个元素取实部，并对该实序列按升序排列，得到一组新的序列：

r'＝{Re(r_i),Re(r_j),...,Re(r_k)},(i,j,k＝1,2,...,P,i≠j≠k)

记序列中元素为r₁',r₂',...,r_P'，定义如下矩阵Q：

将矢量r'与Q相乘得到序列S,S＝r'Q＝[r₁'-r₂',r₂'-r₃',...,r_P-1'-r_P']。

对得到的序列S进行重复特征值判断，具体为：

记S中元素为s₁,s₂,...,s_P,设定误差精度β＝0.05，对S中元素从首位开始判断，若元素S_j满足|S_j|＜β，则继续判断下一个元素，直到元素s_k不满足|S_k|＜β为止，则认为r中元素r_j',...,r_k'为重复特征值，并将该k-j+1个元素归为一组。

重复上述步骤直到所有特征值分组完毕，则可得到若干组重复特征值组：

{r_i',r_j',....,r_k'}，{r_q',r_r',....,r_s'},....(i,j,k,q,r,s＝1,2,...,P,i≠j≠p≠q≠r≠s)

S3、对已分组的实值空间特征值进行映射得到复值空间特征值，对重复特征值参数取算术平均值，构造一维导向矢量；

进一步地，所述步骤S3中，对于第二维参数，由于利用了求根MUSIC方法构造多项式进行计算，因此，具体包括：

S31、对分组后的实值空间特征值映射得到复值空间特征值；

所述步骤S31中，对步骤S21中映射关系式进行变形，得到：

其中，c_p'为对应于r_p的第p个实值空间一维特征值；将步骤S22中得到的若干组重复特征值组代入上述映射公式，得到分组后的复值空间特征值组为：

{c_i',c_j',....,c_k'},{c_q',c_r',....,c_s'},.....

分别记为K₁,K₂,...,且重数分别为J、L....。

S32、对每组重复特征值取算术平均值，构造一维阵列导向矢量；

所述步骤S32中，具体步骤为：

(1)对每组参数求取算数平均：

(2)对参数K₁、K₂...分别构造一维阵列导向矢量：

a(μ_x1)＝(K₁ ⁰,K₁ ¹,...,K₁ ^K-1)，a(μ_x2)＝(K₂ ⁰,K₂ ¹,...,K₂ ^K-1)......

S4、根据噪声子空间与阵列导向矢量的正交性，构造多项式，并通过多项式求根解出第二维方向余弦参数，由角度公式求解方位角与仰角；

进一步地，所述步骤S4中，根据噪声子空间与阵列导向矢量的正交性，得到求根多项式：

其中

将一维阵列导向矢量a(μ_x1)、a(μ_x2)......分别代入上述求根多项式中求根。

(1)在无噪声的情况下，对于J重的特征值a(μ_x1)，有J个根分布在单位圆上，且这J个根表示感兴趣的J个信源的二维方向余弦信息。

(2)在有噪声的情况下，对于J重的特征值a(μ_x1)，则有2J个分布在单位圆附近的根，包含了该J个信源的第二维参数信息，可选取单位圆内最靠近单位的J个根，并根据方向余弦与角度关系公式：

将J重一维方向余弦分别与所求出的J个根代入上述角度计算公式中，求得J个信源的方位角与仰角信息。

将余下的导向矢量a(μ_x2)，a(μ_x3).....通过多项式求得的根代入上述角度计算公式中，直至最后一个参数计算完毕。图4、图5分别显示了方位角与仰角为(30°,60°,60°)的空间信源经由多项式求根得到的单根、重根分布图。如图6所示为方位角与仰角为(30°,60°,30°)的空间信源经由多项式求根得到的重根分布图。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (8)

1.一种矩阵束与求根MUSIC结合的波达方向估计方法，其特征在于，具体步骤包括：

S1、由天线阵列每个阵元接收的叠加信号构造分块Hankel矩阵，并进一步得到自相关数据矩阵，通过奇异值分解获取噪声子空间与信号子空间，构造矩阵束方程，包括：

S11、根据矩阵序列{D₀(t),D₁(t)...,D_M-1(t)}和无噪声信号序列{z_(m,0)(t),z_(m,1)(t),...,z_(m,N-1)(t)}进行处理，构造分块Hankel矩阵；

S12、根据阵元接收模型，对Hankel矩阵进行处理，得到自相关数据矩阵；

S13、对自相关矩阵进行奇异值分解，获取信号子空间和噪声子空间；

S14、构造矩阵束方程；

S2、通过矩阵束方程求解第一维方向余弦，并通过实值映射判断重复特征值；

S3、对已分组的实值空间特征值进行映射得到复值空间特征值，对重复特征值参数取算术平均值，构造一维导向矢量；

S4、根据噪声子空间与阵列导向矢量的正交性，构造多项式，并通过多项式求根解出第二维方向余弦参数，由角度公式求解方位角与仰角；

所述步骤S4中根据正交性方程构造求根多项式，具体为：

根据噪声子空间与阵列导向矢量的正交性，得到求根多项式：

其中

2.根据权利要求1所述的一种矩阵束与求根MUSIC结合的波达方向估计方法，其特征在于，所述步骤S1中构造的分块Hankel矩阵，具体表示为：

其中，K表示x方向上的束参数，M表示x方向上阵元个数，D_m(t)表示x方向上构造的一维Hankel数据矩阵；

其中，D_m(t)具体表示为：

其中，Q表示y方向上的束参数，N表示y方向上阵元个数，z_(m,n)(t)表示t时刻位于第(m,n)个阵元上的无噪接收信号。

3.根据权利要求2所述的一种矩阵束与求根MUSIC结合的波达方向估计方法，其特征在于，所述步骤S1中，通过奇异值分解获取信号子空间与噪声子空间，具体方法为：

所述自相关数据矩阵R_ss，表示方式为：

其中，M、N分别表示x与y方向上阵元数，K、Q分别表示x与y方向上的束参数，D_ce(t)表示中心厄米特数据矩阵，表示方式为：

对R_ss进行奇异值分解，获取信号子空间和噪声子空间，具体分解方式为：

其中，min表示取值数据矩阵行列数中最小者，U_s与V_s分别为信号子空间的左右奇异矢量，U_n与V_n分别表示噪声子空间的左右奇异矢量；

表示对V_s进行共轭转置操作；Σ_s与Σ_n为对角矩阵，其对角线元素分别表示信号子空间与噪声子空间的奇异值。

4.根据权利要求3所述的一种矩阵束与求根MUSIC结合的波达方向估计方法，其特征在于，所述步骤S2中，构造矩阵束方程求取一维方向信息，具体为：

U_s和阵列导向矢量E_L关系式如下：

U_s＝E_LT

其中，T为非奇异矩阵，E_L表示阵列导向矢量，具体为：

其中，X_d表示一维方向余弦，Q_L表示y方向上的一维导向矢量，具体为：

根据U_s的结构及其与E_L关系式，信号子空间分别去掉后Q行和前Q行得到U_s的两个子矩阵U_s1和U_s2；

由U_s＝E_LT可知，U_s1＝E_L1T，U_s2＝E_L2T,同样可知E_L1和E_L2分别为E_L去掉其后Q行和前Q行所得到，且有E_L2＝E_L1X_d；因此，构造矩阵束方程，具体为：

U_s2-λU_s1＝E_L1(X_d-λI)T

得到E_L2＝E_L1T^-1X_dT；定义矩阵Ψ_x＝T^-1X_dT，因此，Ψ_x与X_d为相似矩阵。

5.根据权利要求3所述的一种矩阵束与求根MUSIC结合的波达方向估计方法，其特征在于，所述步骤S2中，实值映射过程具体为：

首先将数据矩阵D_ce(t)进行酉变换，得到

其中，

与

为酉变换矩阵；

酉变换后得到的实值矩阵信号子空间为

推导出实值阵列导向矢量与复值阵列导向矢量的关系为

根据实值阵列导向矢量与复值阵列导向矢量的关系式，对矩阵束方程U_s2-λU_s1＝E_L1(X_d-λI)T进行改写，得到

D₁E_L＝D₂E_LX_d

其中，

因此，有

且

根据改写后矩阵束方程的上述两个等式，分别构造矩阵W₁与W₂，具体为：

因为有

得到

(W₁+jW₂)E_L＝(W₁-jW₂)E_LX_d (1)

对公式(1)改写得到：

W₁E_L(X_d-I)＝jW₂E_L(X_d+I)

进一步得到：

W₁E_Lj(I+X_d)(I-X_d)^-1＝W₂E_L

定义矩阵Y＝j(I+X_d)(I-X_d)^-1，显然

根据矩阵Y与X_d的关系式，将复值空间特征值一一映射到实值空间特征值来进行重复特征值判断，具体映射关系式为：

其中，r_p为对应于c_p的第p个实值空间一维特征值；

根据映射关系式将复值空间映射到实值空间，简化了重复特征值判定的步骤，对c中元素逐个进行实值空间映射，得到了实值矢量空间r＝{r₁,r₂,...,r_P}。

6.根据权利要求5所述的一种矩阵束与求根MUSIC结合的波达方向估计方法，其特征在于，所述步骤S2中重复特征值判断，具体方法为：

在步骤S2中，对r中的每个元素取实部，并对该实序列按升序排列，得到一组新的序列：

r'＝{Re(r_i),Re(r_j),...,Re(r_k)},(i,j,k＝1,2,...,P,i≠j≠k)

记序列中元素为r₁',r₂',...,r_P'，定义如下矩阵Q：

将矢量r'与Q相乘得到序列S,S＝r'Q＝[r₁'-r₂',r₂'-r₃',...,r_P-1'-r_P']；

对得到的序列S进行重复特征值判断，具体为：

记S中元素为s₁,s₂,...,s_P,设定误差精度β＝0.05，对S中元素从首位开始判断，若元素S_j满足|S_j|＜β，则继续判断下一个元素，直到元素s_k不满足|S_k|＜β为止，则认为r中元素r_j',...,r_k'为重复特征值，并将该k-j+1个元素归为一组；

对S中所有元素进行上述重复特征值判断，直到所有特征值分组完毕，则可得到若干组重复特征值组：

{r_i',r_j',....,r_k'}，{r_q',r_r',....,r_s'},....(i,j,k,q,r,s＝1,2,...,P,i≠j≠p≠q≠r≠s)。

7.根据权利要求1所述的一种矩阵束与求根MUSIC结合的波达方向估计方法，其特征在于，所述步骤S3中复值映射，具体方法为：

对步骤S2中的复值空间特征值一一映射到实值空间特征值的具体映射关系式进行变形，得到：

其中，c_p'为对应于r_p的第p个实值空间一维特征值；将步骤S2中得到的若干组重复特征值组代入上述映射公式，得到分组后的复值空间特征值组为：

{c_i',c_j',....,c_k'},{c_q',c_r',....,c_s'},.....

分别记为K₁,K₂,...,且重数分别为J、L....；

对得到的每组赋值空间特征值参数求取算数平均：

再对得到的算术平均值K₁、K₂...分别构造一维阵列导向矢量：

a(μ_x1)＝(K₁ ⁰,K₁ ¹,...,K₁ ^K-1)，a(μ_x2)＝(K₂ ⁰,K₂ ¹,...,K₂ ^K-1)......。

8.根据权利要求1所述的一种矩阵束与求根MUSIC结合的波达方向估计方法，其特征在于，所述步骤S4中根据角度公式求解DOA，具体方法为：

将一维阵列导向矢量a(μ_x1)、a(μ_x2)......分别代入上述多项式中求根；

(1)在无噪声的情况下，对于J重的特征值a(μ_x1)，有J个根分布在单位圆上，且这J个根表示感兴趣的J个信源的二维方向余弦信息；

(2)在有噪声的情况下，对于J重的特征值a(μ_x1)，则有2J个分布在单位圆附近的根，包含了该J个信源的第二维参数信息，选取单位圆内最靠近单位的J个根，并根据方向余弦与角度关系公式：

其中，

和θ_p分别表示仰角与方位角；

将J重一维方向余弦分别与所求出的J个根代入上述角度计算公式中，求得J个信源的方位角与仰角信息；

将余下的导向矢量a(μ_x2)，a(μ_x3).....通过多项式求得的根代入上述角度计算公式中，直至最后一个参数计算完毕。

Images