CN104345306B - 基于Khatri‑Rao子空间的目标波达角估计方法 - Google Patents

Abstract

本发明属于目标波达角估计技术领域，特别涉及基于Khatri‑Rao子空间的目标波达角估计方法。其具体步骤为：首先建立信号模型，然后对接收信号求协方差，并对协方差应用Khatri‑Rao子空间的知识改进，然后消去噪声协方差部分，再对其进行奇异值分解，最后构造求根系数并对其求根，最终得出DOA估计。实验证明，在非均匀线阵的情况下，本发明可以比较精确的估计信源空间方向，而且可以增加最大估计信源个数，且可以减去计算负担。

Description

基于Khatri-Rao子空间的目标波达角估计方法

技术领域

本发明属于目标波达角估计技术领域，特别涉及基于Khatri-Rao子空间的目标波达角估计方法，尤其在空间有多个目标天线阵列为非均匀线性阵列的情况下，本发明不仅可以精确的估计多个目标的角度，而且还能够减少计算量。

背景技术

20世纪40年代自适应天线组合技术发展起来的阵列天线技术，充分显示了它在信号源定位、信源分离、信道参数估计等方面的优势和巨大潜力。阵列信号处理也就成为雷达探测技术的重要领域。在阵列信号处理的主要问题中，空间谱估计用于空间信号方位的超分辨估计，在干扰情况下的现代雷达多目标分辨中有着重要的意义。而且在复杂电磁场环境下对待特定目标的方位检测在现代电子战中有着重要的意义。方位估计的基本问题就是确定同时出现在空间某个区域内的多个信号的空间位置，这就是我们通常说的DOA(波达角)估计。1968年Schmidit提出的MUSIC算法对于空间谱估计的贡献是不可限量的。对于非均匀线性阵列(NLA)，也可以看做均匀线阵中间缺少了部分传感器元件或者部分传感器元件不能工作，这个时候我们就要考虑怎样处理接收的信号来进行参数的估计。一大批基于MUSIC的新方法也出现了，Piya Pal和P.P.Vaidyanathan于2010年发表的signalprocessing的《Nested Arrays:A Novel Approach to Array Processing With EnhancedDegrees of Freedom》中，提出了一种空域滤波处理的方法，想要应用空域滤波来处理虚拟传感器，作者尝试为到达信号的协方差矩阵构造一个Hermitian矩阵。这种方法可以估计到N²/4+N/2-1个阵源，N为非均匀线性阵列的阵元数这种方法提供了一个很好的DOA估计，不管是对于似稳健信号还是稳健信号。然而，它要求虚拟阵列必须是均匀线阵。

发明内容

本发明的目的在于提出基于Khatri-Rao子空间的目标波达角估计方法。本发明在雷达信号接收阵列为非均匀线性阵列的情况下，不仅可以精确估计更多目标，而且可以减少计算量。

为实现上述技术目的，本发明采用如下技术方案予以实现。

基于Khatri-Rao子空间的目标波达角估计方法包括以下步骤：

步骤1，利用雷达发射信号，利用雷达信号接收阵列接收回波信号；雷达信号接收阵列为由N个阵元组成的非均匀线性阵列，用r_j表示所述非均匀线性阵列第j+1阵元与第1个阵元之间的距离，j＝0,1…,N-1；r_j为λ/2的整数倍，其中λ为雷达发射信号的波长；

步骤2，得出所述非均匀线性阵列接收的第m帧回波信号的估计协方差矩阵m＝1,2,...M，M表示所述非均匀线性阵列接收的回波信号的帧数；得出经矢量运算处理后的回波信号的协方差矩阵Y， $Y=[\mathrm{vec}\left({\stackrel{\‾}{R}}_1\right),\mathrm{vec}\left({\stackrel{\‾}{R}}_2\right),...,\mathrm{vec}\left({\stackrel{\‾}{R}}_m\right),...,\mathrm{vec}\left({\stackrel{\‾}{R}}_M\right)],$ 其中，vec(·)为列矢量化运算符；

步骤3，设置投影矩阵I_M为M维的单位矩阵，上标T表示矩阵或向量的转置；1_M为M×1维的列向量，1_M中的每个元素为1；

步骤4，对矩阵进行奇异值分解，其中，U为N²×N²维的酉矩阵，Σ为N²×M维的矩阵，V^*为M×M维的酉矩阵；利用矩阵U的前p列组成信号子空间U_s，利用矩阵U的后N²-p列组成噪声子空间U_n；

步骤5，构造N个逻辑矩阵G₀,G₁,...,G_j,...,G_N-1，j＝0,1,2...N-1；其中，G₀＝[g₀,g₁,…g_j,…g_N-1]^T，上标T表示矩阵或向量的转置，g_j为(2J+1)行的列矢量，J为雷达信号接收阵列的孔径，g_j中第J+1+2r_j/λ个元素为1，其余元素为0；当j≠N-2时，G_j+1＝＜＜(G_j,d_j)，＜＜(G_j,d_j)表示矩阵G_j循环向左移动d_j列得出的新矩阵，d_j＝2(r_j+1-r_j)/λ；

将所述N个逻辑矩阵组合成矩阵G，

$G=\left[\begin{array}{c}{G}_0\\ G_1\\ ...\\ G_j\\ ...\\ G_{N-1}\end{array}\right];$

步骤6，构造矩阵W，上标H表示矩阵的共轭转置；根据矩阵W构造噪声子空间系数矢量c_n，c_n为(4J+1)×1维的矢量；根据所述非均匀线性阵列的阵元位置构造阵元方位系数矢量c_p，c_p为(4J+1)×1维的矢量；得出4J+1个求根系数，令l＝1,2,...,4J+1，则第l个求根系数c(l)为：c(l)＝c_n(l)c_p(l)，c_n(l)表示噪声子空间系数矢量c_n的第l个元素，c_p(l)表示阵元方位系数矢量c_p的第l个元素；

步骤7，构造多项式f(z)：令f(z)＝0，得出关于z的方程，求解所述关于z的方程；将求解出的所有z的解中绝对值超过设定阈值的解保留，对于保留的任一个z的解，根据以下方程得出对应的目标信号的波达角θ：

$z=e^{-j\frac{2\mathrm{\πd}}{\mathrm{\λ}}\sin \left(\mathrm{\θ}\right)}.$

d＝λ/2

本发明的有益效果为：1)采用对ROOT-MUSIC的方法来代替谱搜索，极大的降低了运算负担。2)用噪声子空间结合阵元位置信息获得求根系数矢量，较直接地对谱函数求根，大大的降低了运算难度。

附图说明

图1为本发明的基于Khatri-Rao子空间的目标波达角估计方法的流程图；

图2为仿真实验1中采用本发明和4D-MUSIC方法得出目标信号波达角的最小均方误差和信噪比的关系示意图；

图3为仿真实验2采用本发明和4D-MUSIC方法得出目标信号波达角的最小均方误差和角分布间隔的关系示意图；

图4，为仿真实验3中，采用本发明得出的目标信号波达角和频谱的关系示意图。

具体实施方式

下面结合附图对本发明作进一步说明：

参照图1，为本发明的基于Khatri-Rao子空间的目标波达角估计方法的流程图。该基于Khatri-Rao子空间的目标波达角估计方法包括以下步骤：

步骤1，利用雷达发射信号，利用雷达信号接收阵列接收回波信号；雷达信号接收阵列为由N个阵元组成的非均匀线性阵列，用r_j表示所述非均匀线性阵列第j+1阵元与第1个阵元之间的距离，j＝0,1…,N-1；r_j为λ/2的整数倍，其中λ为雷达发射信号的波长。

其具体步骤为：

利用雷达发射信号(窄带信号)，利用雷达信号接收阵列接收回波信号(远场窄带回波信号)。本发明实施例中，雷达信号接收阵列为由N个阵元组成的非均匀线性阵列(NLA)，对于非均匀线性阵列，用r_j表示第j+1阵元与第1个阵元之间的距离，j＝0,1…,N-1；r_j为λ/2的整数倍，其中λ为雷达发射信号的波长。显然，r₀＝0，则得出阵元位置矢量r＝[r₀,r₁,...,r_N-1]。

本发明实施例中，雷达信号接收阵列接收的第m帧回波信号X_m为：

X_m＝AS_m+N_m

其中，m＝1,2,...M，M表示雷达信号接收阵列接收的回波信号的帧数。

X_m＝[x(1)_m,x(2)_m,....x(k)_m]，k表示每帧回波信号的快拍数。x(1)_m至x(k)_m分别表示第m帧回波信号的第1次快拍数据至第m帧回波信号的第k次快拍数据。S_m表示第m帧回波信号的包络，N_m表示第m帧回波信号的空间噪声。

A＝[a(θ₁),a(θ₂),...a(θ_i),...a(θ_p)]，A表示p个目标的阵列流型，p为大于1的自然数；a(θ_i)表示第i个目标的导向矢量，i＝1,2,...,p；a(θ_i)为：

$a\left({\mathrm{\θ}}_i\right)={[1,e^{-\frac{j2\mathrm{\π}{r}_1}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_i\right)},...e^{-\frac{j2\mathrm{\π}{r}_{N-1}}{\mathrm{\λ}}\sin \left({\mathrm{\θ}}_i\right)}]}^T$

其中，θ_i表示第i个目标的波达角，上标T表示矩阵或向量的转置。

步骤2，得出所述非均匀线性阵列接收的第m帧回波信号的估计协方差矩阵m＝1,2,...M，M表示所述非均匀线性阵列接收的回波信号的帧数；得出经矢量运算处理后的回波信号的协方差矩阵Y， $Y=[\mathrm{vec}\left({\stackrel{\‾}{R}}_1\right),\mathrm{vec}\left({\stackrel{\‾}{R}}_2\right),...,\mathrm{vec}\left({\stackrel{\‾}{R}}_m\right),...,\mathrm{vec}\left({\stackrel{\‾}{R}}_M\right)],$ 其中，vec(·)为列矢量化运算符。

其具体子步骤为：

(2.1)第m帧回波信号的协方差矩阵R_m的表达式为：

$R_m=E\left\{X_m{X}_m^H\right\}=\mathrm{AE}\left\{S_m{S_m}^H\right\}{A}^H+E\left\{N_m{N_m}^H\right\}$

其中，上标H表示矩阵的共轭转置，E{·}表示期望，令E{S_mS_m ^H}＝D_m，令E{N_mN_m ^H}＝C，D_m的展开形式为：D_m＝diag(d_m1,d_m2,...d_mp)，diag(d_m1,d_m2,...d_mp)表示以d_m1至d_mp作为主对角线元素而构成的对角矩阵，p表示目标数；C表示每帧回波信号的空间噪声的协方差矩阵。因此，R_m＝AD_mA^H+C，AD_mA^H为第m帧回波信号的空间信号的协方差矩阵。

本发明实施例中，在步骤1之后，得出第m帧回波信号的估计协方差矩阵

${\stackrel{\‾}{R}}_m=\frac{1}{k}\underset{i^{\′}=1}{\overset{k}{\mathrm{\Σ}}}x{\left(i^{\′}\right)}_{m}x{{\left(i^{\′}\right)}_m}^H$

其中，k表示每帧回波信号的快拍数，x(i)_m表示第m帧回波信号的第i'次快拍数据，上标H表示矩阵的共轭转置，i'取1至k。此时，

(2.2)然后，根据Khatri-Rao子空间的概念，对每一帧的回波信号的协方差矩阵进行矢量化处理。具体地说，定义即：

$y_m=\mathrm{vec}(A{D}_m{A}^H+C)=\mathrm{vec}\left(A{D}_m{A}^H\right)+\mathrm{vec}\left(C\right)=(A^{*}\⊗A)d_m+\mathrm{vec}\left(C\right)\≈\mathrm{vec}\left({\stackrel{\‾}{R}}_m\right)$

其中，vec(·)为列矢量化运算符(假设任意矩阵B＝[a₁,a₂,...a_n]，则上标*表示取共轭，是Khatri-Rao乘积，d_m＝[d_m1,d_m2,...d_mp]^T，m＝1,2…M，M表示雷达信号接收阵列接收的回波信号的帧数。

这样，经矢量运算处理后的回波信号的协方差矩阵Y表示为：

$Y=(A^{*}\⊗A){\mathrm{\Ψ}}^T+\mathrm{vec}\left(C\right)1_M^T\≈[\mathrm{vec}\left({\stackrel{\‾}{R}}_1\right),\mathrm{vec}\left({\stackrel{\‾}{R}}_2\right),...,\mathrm{vec}\left({\stackrel{\‾}{R}}_m\right),...,\mathrm{vec}\left({\stackrel{\‾}{R}}_M\right)]$

其中，上标*表示取共轭，表示Khatri-Rao乘积，上标T表示矩阵或向量的转置。1_M为M×1维的列向量，1_M中的元素全为1，M表示雷达信号接收阵列接收的回波信号的帧数，Ψ＝[d₁,d₂...d_m...d_M]^T，Ψ表示信号包络自相关矩阵的矢量化处理。C＝E{N_mN_m ^H}，上标H表示矩阵的共轭转置，E{·}表示期望；N_m表示第m帧回波信号的空间噪声。

因此，根据至即可得出经矢量运算处理后的回波信号的协方差矩阵Y。

步骤3，设置投影矩阵I_M为M维的单位矩阵，上标T表示矩阵或向量的转置；1_M为M×1维的列向量，1_M中的每个元素为1。

其具体步骤为：

本发明实施例中，通过一个投影矩阵来消除未知的每帧回波信号的空间噪声的协方差矩阵C。具体地说，设置投影矩阵使其满足：

$\begin{array}{c}{\mathrm{YP}}_{1M}^{\⊥}=((A^{*}\⊗A){\mathrm{\Ψ}}^T+\mathrm{vec}\left(C\right)1_M^T)P_{1M}^{\⊥}\\ =(A^{*}\⊗A){\left(P_{1M}^{\⊥}\mathrm{\Ψ}\right)}^T\end{array}$

其中，上标*表示取共轭，表示Khatri-Rao乘积，上标T表示矩阵或向量的转置。1_M为M×1维的列向量，1_M中的元素全为1，M表示雷达信号接收阵列接收的回波信号的帧数，Ψ＝[d₁,d₂...d_m...d_M]^T，vec(·)表示列矢量化处理。可知，I_M为M维的单位矩阵，因此，Ψ和的秩相等，当和都是列满秩时，有：

其中，(·)表示·的值域(也称·的列向量张成的子空间集合)。

步骤4，对矩阵进行奇异值分解，其中，U为N²×N²维的酉矩阵，Σ为N²×M维的矩阵，V^*为M×M维的酉矩阵；利用矩阵U的前p列组成信号子空间U_s，利用矩阵U的后N²-p列组成噪声子空间U_n。

具体地说，对矩阵按照以下公式进行奇异值分解：

${\mathrm{YP}}_{1M}^{\⊥}={\mathrm{U\ΣV}}^{*}=U\left[\begin{array}{cc}{\mathrm{\Σ}}_s& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{V}_s^H\\ V_n^H\end{array}\right]=[U_s,U_n]\left[\begin{array}{cc}{\mathrm{\Σ}}_s& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{V}_s^H\\ V_n^H\end{array}\right]$

其中，U为N²×N²维的酉矩阵，Σ为N²×M维的矩阵，V^*为M×M维的酉矩阵；矩阵U的前p列组成信号子空间U_s，U_s为N²×p维的矩阵；矩阵U的后N²-p列组成噪声子空间U_n，U_n为N²×(N²-p)维的矩阵，U＝[U_s，_Un]； $\mathrm{\Σ}=\left[\begin{array}{cc}{\mathrm{\Σ}}_s& 0\\ 0& 0\end{array}\right],V^{*}=\left[\begin{array}{c}{V}_s^H\\ V_n^H\end{array}\right],$ ∑_s表示信号的奇异值对角矩阵，上标H表示矩阵的共轭转置。

步骤5，构造N个与天线阵元位置有关的逻辑矩阵，N个逻辑矩阵依次表示为G₀,G₁,...,G_j,...,G_N-1，j＝0,1,2...N-1；其中，G₀＝[g₀,g₁,…g_j,…g_N-1]^T，上标T表示矩阵或向量的转置，g_j为(2J+1)行的列矢量，J为雷达信号接收阵列的孔径(阵列天线孔径)，g_j中第J+1+2r_j/λ个元素为1，其余元素为0，r_j表示第j+1阵元与第1个阵元之间的距离，G₀为N×(2J+1)维的矩阵；当j≠N-2时，G_j+1＝＜＜(G_j,d_j)，＜＜(G_j,d_j)表示矩阵G_j循环向左移动d_j列得出的新矩阵，d_j＝2(r_j+1-r_j)/λ。

将所述N个逻辑矩阵组合成矩阵G，

$G=\left[\begin{array}{c}{G}_0\\ G_1\\ ...\\ G_j\\ ...\\ G_{N-1}\end{array}\right]$

则G为N²×(2J+1)维的矩阵。

步骤6，构造矩阵W，上标H表示矩阵的共轭转置；根据矩阵W构造噪声子空间系数矢量c_n，c_n为(4J+1)×1维的矢量，J为雷达信号接收阵列的孔径；根据所述非均匀线性阵列的阵元位置构造阵元方位系数矢量c_p，c_p为(4J+1)×1维的矢量；得出4J+1个求根系数，令l＝1,2,...,4J+1，则第l个求根系数c(l)为：c(l)＝c_n(l)c_p(l)，c_n(l)表示噪声子空间系数矢量c_n的第l个元素，c_p(l)表示阵元方位系数矢量c_p的第l个元素。

其具体子步骤为：

(6.1)构造矩阵W，上标H表示矩阵的共轭转置，矩阵W为(2J+1)×(2J+1)维的矩阵。

求解噪声子空间系数矢量c_n：

c_n＝[c_n(1),c_n(2)...,c_n(l),...,c_n(4J+1)]^T

其中，上标T表示矩阵或向量的转置，c_n(l)为矩阵W中的第l-2J-1个对角线元素的和，l＝1,2,...,4J+1；矩阵W的第0个对角线元素依次为a_1，1,a_2，2,…a_2J+1,2J+1，a_x,y表示矩阵W第x行第y列的元素，x取1至2J+1，y取1至2J+1；当l-2J-1>0时，矩阵W的第l-2J-1个对角线元素依次为a_1,l-2J,a_2,l-2J+1,…a_4J+2-l,2J+1；当l-2J-1<0时，矩阵W的第l-2J-1个对角线元素依次为a_2J+2-l,1,a_2J+3-l,2,…a_2J+1,l。

(6.2)求解阵元方位系数矢量c_p。

求解阵元方位系数矢量c_p的过程为：定义阵元存在的逻辑矢量为c_t，如果雷达信号接收阵列为均匀线阵，则阵元存在的逻辑矢量c_t中的每个元素为1；如果雷达信号接收阵列为非均匀线阵，则阵元存在的逻辑矢量c_t中会有0元素。

当雷达信号接收阵列为非均匀线阵时，阵元存在的逻辑矢量c_t中每个元素按如下方式取值：从雷达信号接收阵列的第一个阵元开始，沿着雷达信号接收阵列的排列方向每隔λ/2获取一个阵元位置，直至获取到雷达信号接收阵列的最后第一个阵元所在位置；将获取的第l个阵元位置表示为p_l，l＝1,2,...，获取的第1个阵元位置p₁为雷达信号接收阵列的第一个阵元所在位置；如果获取的第l个阵元位置p_l有阵元存在，则阵元存在的逻辑矢量c_t中第l个元素取1，反之，如果获取的第l个阵元位置p_l没有阵元存在，则阵元存在的逻辑矢量c_t中第l个元素取0。

c_p为逻辑虚拟位置的自卷积，而虚拟位置又是表示阵元存在的逻辑矢量c_t的自卷积，即c_p＝(c_t*c_t)*(c_t*c_t)，其中*表示卷积符号，c_p为(4J+1)×1维的矢量，其元素由0和1组成。矢量c_p的展开式为：

c_p＝[c_p(1),c_p(2)...,c_p(l),...,c_p(4J+1)]^T

其中，l＝1,2,...,4J+1。

(6.3)噪声子空间系数矢量c_n非0元素的个数要比c_p中的多，所以，利用方位系数矢量c_p把c_n中对应位置的非0元素置0，进而得出求根系数矢量c，

c＝[c(1),c(2)...,c(l),...,c(4J+1)]

其中，c(l)＝c_n(l)c_p(l)，l＝1,2,...,4J+1。

步骤7，构造多项式f(z)：令f(z)＝0，得出关于z的方程，求解所述关于z的方程；将求解出的所有z的解中绝对值超过设定阈值的解保留，然后针对保留的每个z的解，得出对应的目标信号的波达角。

对于保留的任一个z的解，根据以下方程得出对应的目标信号的波达角θ：

$z=e^{-j\frac{2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin \left(\mathrm{\&theta;}\right)}$

d＝λ/2

其中，λ为雷达发射信号的波长。

本发明的效果可以通过以下仿真实验进一步说明。

在仿真实验中，我们设定每一帧信号长度k＝512，帧次数设定为M＝50。

仿真实验1：考虑最小均方误差RMSE(也就是估计误差)和信噪比SNR的关系。在仿真实验1中，雷达信号接收阵列阵元数为4，阵元位置矢量r＝[0,2,3,8]，共有4个目标信号，它们的波达角分别为-40°、-10°、10°、30°，信噪比变化范围为-3dB到15dB。在仿真实验1中，分别采用本发明和4D-MUSIC方法得出目标信号波达角的最小均方误差。参照图2，为仿真实验1中采用本发明和4D-MUSIC方法得出目标信号波达角的最小均方误差和信噪比的关系示意图。图2中，横轴表示信噪比，单位为dB，纵轴表示目标信号波达角的最小均方误差，单位为度。从图2中可以看出，本发明的目标信号波达角的最小均方误差大约仅仅是4D-MUSIC方法的一半。

仿真实验2：考虑均方误差与角度分离的关系。在仿真实验2中，雷达信号接收阵列阵元数为4，阵元位置矢量r＝[0,2,3,8]，目标信号的波达角分别为-15°、Δ、15°，信噪比为10dB，Δ表示一个变化的目标信号的波达角(角分布间隔)，其变化范围为[-14.7°,14.7°]。在仿真实验2中，分别采用本发明和4D-MUSIC方法得出目标信号波达角的最小均方误差。参照图3，为仿真实验2采用本发明和4D-MUSIC方法得出目标信号波达角的最小均方误差和角分布间隔的关系示意图。图3中，横轴表示角分布间隔，单位为度，纵轴表示目标信号波达角的最小均方误差，单位为度。从图3可以看出，当波达方向相互之间不是太近时，均方误差是相当的稳定的；其次，当目标的角分布间隔靠近1°时，本发明的目标信号波达角的最小均方误差小于4D-MUSIC方法，仍然能够很好的工作，而4D-MUSIC方法需要角分布间隔大于4°才能确保它正常工作。

仿真实验3：我们观察非均匀线阵下，阵元数一定情况下，估计的目标数和精确度。在仿真实验3中，雷达信号接收阵列阵元数为4，阵元位置矢量r＝[0,1,3,7]，目标的波达角分别为65°、-50°、-40°、-30°、-20°、-10°、0°、10°、25°、40°、50°和70°。在仿真实验3中，采用本发明得出的目标信号波达角和归一化频率的关系示意图。参照图4，为仿真实验3中，采用本发明得出的目标信号波达角和频谱的关系示意图，图4中，横轴表示角度，单位为度，纵轴表示归一化频率。图4中，12条竖直平行的线与12个波达角对应，从图4可以看出，频谱的峰值出现在12个波达角的位置，说明本发明可以很精确的区分12个目标。

综上所述，在处理非均匀线阵的情况下，本发明不仅可以更加精确得进行DOA估计，而且可以增加最大信源估计数目，且减少了计算负担，为相应的信号处理带来不小的方便。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (4)

1.基于Khatri-Rao子空间的目标波达角估计方法，其特征在于，包括以下步骤：

步骤1，利用雷达发射信号，利用雷达信号接收阵列接收回波信号；雷达信号接收阵列为由N个阵元组成的非均匀线性阵列，用r_j表示所述非均匀线性阵列第j+1阵元与第1个阵元之间的距离，j＝0，1…，N-1；r_j为λ/2的整数倍，其中λ为雷达发射信号的波长；

步骤2，得出所述非均匀线性阵列接收的第m帧回波信号的估计协方差矩阵m＝1，2，...M，M表示所述非均匀线性阵列接收的回波信号的帧数；得出经矢量运算处理后的回波信号的协方差矩阵Y，其中，vec(·)为列矢量化运算符；

步骤3，设置投影矩阵 I_M为M维的单位矩阵，上标T表示矩阵或向量的转置；1_M为M×1维的列向量，1_M中的每个元素为1；

步骤4，对矩阵进行奇异值分解，其中，U为N²×N²维的酉矩阵，∑为N²×M维的矩阵，V^*为M×M维的酉矩阵；利用矩阵U的前p列组成信号子空间U_s，利用矩阵U的后N²-p列组成噪声子空间U_n；p为雷达信号接收阵列接收的回波信号中目标信号的个数；

步骤5，构造N个逻辑矩阵G₀，G₁，...，G_j，...，G_N+1，j＝0，1，2...N-1；其中，G₀＝[g₀，g₁，…g_j，…g_N-1]^T，上标T表示矩阵或向量的转置，g_j为(2J+1)行的列矢量，J为雷达信号接收阵列的孔径，g_j中第J+1+2r_j/λ个元素为1，其余元素为0；当j≠N-2时，G_j+1＝＜＜(G_j，d_j)，＜＜(G_j，d_j)表示矩阵G_j循环向左移动d_j列得出的新矩阵，d_j＝2(r_j+1-r_j)/λ；

将所述N个逻辑矩阵组合成矩阵G，

$G=\left[\begin{array}{c}{G}_0\\ G_1\\ ...\\ G_j\\ ...\\ G_{N-1}\end{array}\right];$

步骤6，构造矩阵W，上标H表示矩阵的共轭转置；根据矩阵W构造噪声子空间系数矢量c_n，c_n为(4J+1)×1维的矢量；根据所述非均匀线性阵列的阵元位置构造阵元方位系数矢量c_p，c_p为(4J+1)×1维的矢量；得出4J+1个求根系数，令l＝1，2，...，4J+1，则第l个求根系数c(l)为：c(l)＝c_n(l)c_p(l)，c_n(l)表示噪声子空间系数矢量c_n的第l个元素，c_p(l)表示阵元方位系数矢量c_p的第l个元素；

步骤7，构造多项式f(z)：令f(z)＝0，得出关于z的方程，求解所述关于z的方程；将求解出的所有z的解中绝对值超过设定阈值的解保留，对于保留的任一个z的解，根据以下方程得出对应的目标信号的波达角θ：

$z=e^{-j\frac{2\mathrm{\&pi;}d}{\mathrm{\&lambda;}}sin\left(\mathrm{\&theta;}\right)}$

d＝λ/2。

2.如权利要求1所述的基于Khatri-Rao子空间的目标波达角估计方法，其特征在于，在步骤2中，所述非均匀线性阵列接收的第m帧回波信号的估计协方差矩阵为：

${\stackrel{\&OverBar;}{R}}_m=\frac{1}{k}\underset{i=1}{\overset{k}{\&Sigma;}}x{\left(i^{\&prime;}\right)}_{m}x{{\left(i^{\&prime;}\right)}_m}^H$

其中，k表示每帧回波信号的快拍数，x(i′)_m表示第m帧回波信号的第i′次快拍数据，上标H表示矩阵的共轭转置，i′取1至k。

3.如权利要求1所述的基于Khatri-Rao子空间的目标波达角估计方法，其特征在于，在步骤6中，噪声子空间系数矢量c_n为：

c_n＝[c_n(1)，c_n(2)...，c_n(l)，...，c_n(4J+1)]^T

其中，上标T表示矩阵或向量的转置，c_n(l)为矩阵W中的第l-2J-1个对角线元素的和，l＝1，2，...，4J+1；矩阵W的第0个对角线元素依次为a_1，1，a_2，2，…a_{2J+1，2J+1，}a_x，y表示矩阵W第x行第y列的元素，x取1至2J+1，y取1至2J+1；当l-2J-1＞0时，矩阵W的第l-2J-1个对角线元素依次为a_1，l-2J，a_2，l-2J+1，…a_{4J+2-l，2J+1}；当l-2J-1＜0时，矩阵W的第l-2J-1个对角线元素依次为a_2J+2-l，1，a_2J+3-l，2，…a_2J+1，l。

4.如权利要求1所述的基于Khatri-Rao子空间的目标波达角估计方法，其特征在于，在步骤6中，求解阵元方位系数矢量c_p的过程为：定义阵元存在的逻辑矢量c_t，从雷达信号接收阵列的第一个阵元开始，沿着雷达信号接收阵列的排列方向每隔λ/2获取一个阵元位置，直至获取到雷达信号接收阵列的最后第一个阵元所在位置；将获取的第l个阵元位置表示为p_l，l＝1，2，...，获取的第1个阵元位置p₁为雷达信号接收阵列的第一个阵元所在位置；如果获取的第l个阵元位置p_l有阵元存在，则阵元存在的逻辑矢量c_t中第l个元素为1，反之，如果获取的第l个阵元位置p_l没有阵元存在，则阵元存在的逻辑矢量c_t中第l个元素为0；然后得出阵元方位系数矢量c_p，c_p＝(c_t*c_t)*(c_t*c_t)，其中，*表示卷积符号。