CN103364762B - 任意阵列流形的单基地mimo雷达波达方向估计方法 - Google Patents

Abstract

本发明公开了一种基于任意阵列流形的单基地MIMO雷达波达方向估计方法。主要解决现有技术只适用于线性阵列和计算量大的问题。其实现步骤是：1）根据阵列流形写出单基地MIMO雷达的导向矢量；2）对MIMO雷达的导向矢量流形展开，获得范得蒙导向矢量，采用空间平滑思想进行导向矢量变换，求出从MIMO导向矢量到变换后导向矢量的转换矩阵；3）利用接收阵列与发射波形进行匹配滤波，形成自相关矩阵；4）对自相关矩阵进行特征分解，得到特征值和特征向量，选取特征向量形成噪声子空间；5）利用噪声子空间形成空间零谱函数，采用多项式求根方法得到方位角。本发明可实现对任意阵列流形的单基地MIMO雷达波达方向快速估计，且运算量小，可用于目标定位和雷达追踪。

Description

任意阵列流形的单基地MIMO雷达波达方向估计方法

技术领域

本发明属于雷达技术领域，具体涉及单基地多输入多输出MIMO雷达的波达方向估计，可应用于目标定向和雷达追踪。

背景技术

MIMO雷达可充分利用波形分集增益，增加系统的自由度，改善其空间分辨率，提高参数估计精度和增加系统最大可定位目标数量。近年来MIMO雷达发展非常迅速，其中波达方向估计DOA方法是MIMO雷达研究的重点问题。在各种阵列结构中，均匀线阵因其结构简单、实现容易，并能采用各种快速DOA算法而成为目前MIMO雷达DOA算法研究的基础。但一维均匀线阵只能提供180度的无模糊方位角信息，且有时特定系统的需要也使得阵列可能不是均匀线阵。对于MIMO雷达来说，发射时不形成空间波束，其阵元方向图覆盖整个空域，用来探测360度的整个空间信息，这时采用线性阵列在360度范围内会出现模糊，就需要二维的阵列，因此针对任意阵列流行的波达方向估计成为亟待解决的问题。

目前，最基本的角度超分辨估计方法是采用多信号分类MUSIC方法，这种方法虽然可在任意阵列流形下工作，但却需要对角度进行全空域搜索，运算量非常大。

为了降低运算量，可采用Richard H.Roy和Kailath等人于1986年提出的旋转不变技术估计信号参数ESPRIT方法和对MUSIC方法改进后的多项式求根多信号分类root-MUSIC方法，这些方法虽避免了角度谱峰搜索，降低了DOA的运算量，但这两种方法只能用于均匀线性阵列，对于任意阵列却无能为力。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于任意阵列流形的单基地MIMO雷达波达方向估计方法，以在减小运算量的同时，实现对任意阵列MIMO雷达的波达方向快速估计。

实现本发明目的的技术思路是：利用流形展开技术把MIMO雷达导向矢量展开并按一定的准则截断误差项，把导向矢量分解为矩阵与范得蒙形式的导向矢量相乘的形式，并借鉴空间平滑的思想，使得角度估计算法更加稳健；采用多项式求根方法获得目标角度信息，取代传统方法中采用谱峰搜索的方式，实现波达方向快速估计。其技术步骤包括如下：

1）根据阵列流形写出单基地MIMO雷达的导向矢量A(θ)，其中θ为方位角；

2）对MIMO雷达的导向矢量A(θ)流形展开，获得范得蒙导向矢量D(θ)，对D(θ)进行空间平滑处理，处理后的导向矢量记为D_S(θ)，利用最小二乘法得出从MIMO雷达导向矢量A(θ)到平滑处理后的导向矢量D_S(θ)的转换矩阵G；

3）对雷达的接收数据与发射波形进行匹配滤波，将匹配滤波结果记为X(t)，再利用滤波结果X(t)形成自相关矩阵R_x；

4）对自相关矩阵R_x进行特征分解，得到一系列特征值λ_m及其对应的特征向量e_m，m=1,2,…,N²,其中N为阵元数；将特征值λ_m进行从小到大排序，依次选取N²-P个小特征值，P为同一距离门内的目标个数，取这些特征值对应的特征向量形成噪声子空间 ${\hat{E}}_n=[e_{P+1},\·\·\·e_m\·\·\·,e_{N^2}];$

5）利用噪声子空间形成MUSIC空间零谱函数： $f\left(\mathrm{\θ}\right)={\left|\right|{{\hat{E}}_n}^{H}A\left(\mathrm{\θ}\right)\left|\right|}^2=A^H\left(\mathrm{\θ}\right){\hat{E}}_n{{\hat{E}}_n}^{H}A\left(\mathrm{\θ}\right),$ 其中|| ||²表示2-范数，表示矩阵的共轭转置；

6）把步骤5）中的A(θ)用步骤2）中的GD_s(θ)代替，并令z=e^jθ，则f(θ)可转化为：令f(z)=0，采用多项式求根的方法求得P个方位角的复指数形式z_k，其中k=1,2…P；

7）利用接收角的复指数形式z_k求得第k个接收方位角θ_k：θ_k=angle(z_k)，其中angle()表示求相位角。

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

（1）现有的不依赖流形的DOA估计方法是MUSIC谱估计方法，该方法需要全空域搜索。尤其对于MIMO雷达，导向矢量的维数成平方性增加，当搜索间隔很小时，传统的MUSIC算法的运算量达到O(10⁸)，而本发明把空域搜索变为多项式求根得到波达方位角，运算量可降为O(10⁵)；

（2）目前可减小运算量的快速算法为旋转不变技术估计信号参数ESPRIT方法和多项式求根多信号分类root-MUSIC方法，但这些方法局限于对线性均匀阵列的处理，对于任意阵列无能为力。本发明采用Jacobi–Anger展开技术，可把MIMO雷达导向矢量变为范得蒙形式，适用于任意阵列。

附图说明

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

图2是利用本发明对目标20次独立实验所得的结果；

图3是利用本发明对目标1定位时的目标角度均方误差随信噪比SNR变化图；

图4是利用本发明对目标2定位时的目标角度均方误差随信噪比SNR变化图。

具体实施方式

参照图1，本实施例的具体实现步骤如下：

步骤1，根据单基地MIMO雷达的阵元分布情况，获得该矩阵导向矢量A(θ)。

根据单基地MIMO雷达的阵元分布，获得各阵元的极坐标为(r_n,β_n)，n=1,2,…N，N为阵元数；

根据各阵元的极坐标(r_n,β_n)，得到MIMO雷达的发射导向矢量a_t(θ)和接收导向矢量a_r(θ)分别为：

$a_t\left(\mathrm{\θ}\right)={[e^{\mathrm{jk}{r}_1\cos ({\mathrm{\β}}_1-\mathrm{\θ})}\·\·\·e^{\mathrm{jk}{r}_n\cos ({\mathrm{\β}}_n-\mathrm{\θ})}\·\·\·e^{\mathrm{jk}{r}_N\cos ({\mathrm{\β}}_N-\mathrm{\θ})}]}^T,$

$a_r\left(\mathrm{\θ}\right)={[e^{\mathrm{jk}{r}_1\cos ({\mathrm{\β}}_1-\mathrm{\θ})}\·\·\·e^{\mathrm{jk}{r}_n\cos ({\mathrm{\β}}_n-\mathrm{\θ})}\·\·\·e^{\mathrm{jk}{r}_N\cos ({\mathrm{\β}}_N-\mathrm{\θ})}]}^T,$

其中，k=2π/λ，λ为波长，r_n和β_n分别为第n个阵元距坐标原点的距离和在极坐标中相对x轴的角度位置，θ为方位角，[]^T表示矩阵转置；

根据发射导向矢量a_t(θ)和接收导向矢量a_r(θ)获得MIMO雷达的导向矢量A(θ)为：

$A\left(\mathrm{\θ}\right)=a_t\left(\mathrm{\θ}\right)\⊗a_r\left(\mathrm{\θ}\right),$

式中，表示Kronecker积运算；

步骤2，对MIMO雷达的导向矢量A(θ)流形展开，获得范得蒙导向矢量D(θ)，对D(θ)进行空间平滑处理，处理后的导向矢量记为D_S(θ)，利用最小二乘法得出从MIMO雷达导向矢量A(θ)到平滑处理后的导向矢量D_S(θ)的转换矩阵G。

2a）对MIMO雷达的导向矢量A(θ)流形展开，获得范得蒙导向矢量D(θ)：

令R=max(r_n)，r_n为第n个阵元距坐标原点的距离，n=1,2,…N；选择范得蒙导向矢量D(θ)的维数M>8kR，k=2π/λ，λ为波长；

忽略截断误差，根据所述维数M，得到范得蒙导向矢量D(θ)：

$D\left(\mathrm{\θ}\right)={[e^{\frac{-j(M-1)\mathrm{\θ}}{2}},\·\·\·,e^{\frac{j(M-1)\mathrm{\θ}}{2}}]}^T,$

其中，[]^T表示矩阵转置，j表示虚数；

2b）对D(θ)进行空间平滑处理，得到变换后的导向矢量D_S(θ)：

将范得蒙导向矢量D(θ)划分为相互重叠的L个子阵，对应每个子阵的阵元数为Q=M-L+1，其中，M为范得蒙导向矢量D(θ)的维数；

根据所述阵元数Q构造平滑矩阵：F_k=[0_Q×(k-1)|I_Q|0_Q×(L-k)]，其中，k=1,2,…L，0_Q×(k-1)表示Q×(k-1)阶零阵，I_Q表示Q阶单位矩阵，0_Q×(L-k)表示Q×(L-k)阶零阵；

将范得蒙形式的导向矢量D(θ)与平滑矩阵F_k相乘，得到空间平滑后的导向矢量D_S(θ)为：

$D_S\left(\mathrm{\θ}\right)=\frac{1}{L}\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{F}_{k}D\left(\mathrm{\θ}\right);$

2c）利用最小二乘法得出从MIMO雷达导向矢量A(θ)到平滑处理后的导向矢量D_S(θ)的转换矩阵G：

在0-360度范围内等间隔的取N个角度，即：θ₁,θ₂,...,θ_N，N为阵元数；

根据θ₁,θ₂,...,θ_N，分别形成最小二乘导向矩阵A^c和最小二乘平滑矩阵D_c：

A^c=[A(θ₁),A(θ₂),…,A(θ_N)]，

D_c=[D_s(θ₁),D_s(θ₂),…,D_s(θ_N)]；

根据最小二乘导向矩阵A^c和最小二乘平滑矩阵D_c，利用最小二乘法得到转换矩阵G：

G=A^cD_c ^H(D_cD_c ^H)^-1，

其中，D_c ^H表示矩阵D_c的共轭转置，()^-1表示矩阵求逆。

步骤3，对雷达的接收数据与发射波形进行匹配滤波，并将滤波结果记为X(t)，再利用X(t)形成自相关矩阵R_x。

3a）对雷达的接收数据与发射波形进行匹配滤波：

假设各目标的雷达截面积RCS相同，则接收信号表示为：

$X_r=\underset{k=1}{\overset{P}{\mathrm{\Σ}}}A\left({\mathrm{\θ}}_k\right)S{e}^{j2\mathrm{\π}{f}_{\mathrm{dk}}t}+V,$

其中，P为同一距离门内的目标个数，θ_k为第k个目标的方向角，f_dk是第k个目标的多普勒频率，k=1,2,…P，V是接收阵列的噪声矩阵;

对接收阵列分别用发射波形进行匹配滤波，得到匹配滤波后的结果X(t)：

X(t)=E(SX_r)，

其中，E()表示数学期望,S=[s₁ ^T,…s_i ^T…,s_N ^T]^T表示由N个发射信号组成的矩阵，s_i表示第i个发射信号，表示s_i的转置，i=1…N；

3b）利用匹配滤波后的结果形成自相关矩阵R_x：

R_x=E{X(t)X^H(t)}，

其中E{}表示数学期望，X^H(t)表示X(t)的共轭转置。

步骤4，对自相关矩阵R_x进行特征分解，得到噪声子空间。

对自相关矩阵R_x进行特征分解，得到一系列特征值λ_m及其对应的特征向量e_m，m=1,2,…,N²；

将特征值λ_m进行从小到大排序，依次选取N²-P个小特征值，P为同一距离门内的目标个数，取这些特征值对应的特征向量形成噪声子空间

${\hat{E}}_n=[e_{P+1},\·\·\·e_i\·\·\·,e_{N^2}],$

其中，i=P+1,…N²。

步骤5，利用噪声子空间形成MUSIC空间谱函数f(θ)：

$f\left(\mathrm{\θ}\right)={\left|\right|{{\hat{E}}_n}^{H}A\left(\mathrm{\θ}\right)\left|\right|}^2=A^H\left(\mathrm{\θ}\right){\hat{E}}_n{{\hat{E}}_n}^{H}A\left(\mathrm{\θ}\right),$

其中，表示的共轭转置，|| ||²表示2-范数；

步骤6，把步骤2）中的A(θ)用GD_s(θ)代替，并令z=e^jθ，则：

$f\left(z\right)={D_s}^H\left(z\right)G^H{\hat{E}}_n{{\hat{E}}_n}^H{\mathrm{GD}}_s\left(z\right),$

令f(z)=0，采用多项式求根的方法求得P个方位角的复指数形式z_p，其中p=1…P；

步骤7，利用接收角的复指数形式z_p求得第p个接收方位角θ_p：

θ_p=angle(z_p)，

其中angle()表示求相位角。

至此求得信号源的方位角，完成对雷达波达方向的估计。

本发明的效果通过以下仿真试验进一步说明：

1．仿真条件：

不失一般性，令阵列均为均匀圆形阵列，发射阵列数为8，阵列半径为0.6倍波长，发射信号为同载频正交的相位编码信号，设存在两个信号源，信号源1的方位角为20°，信号源2的方位角为40°，接收脉冲重复周期数为100，信噪比SNR＝10dB，导向矢量傅立叶级数截断长度M=23，空间平滑子阵长度为Q=21，进行20次独立的蒙特卡罗实验，信号源p的均方误差计算采用公式其中为信号源p方位角的估计值，θ_p为信号源p的方位角。

2．仿真内容：

仿真1，采用本发明对目标方位角度进行目标定位仿真，仿真结果如图2所示。从图2可以看出，采用本发明可以实现对任意阵列单基地MIMO雷达快速波达方向的估计。

仿真2，采用本发明对信号源1定位时，仿真目标角度均方误差随信噪比SNR变化，仿真结果如图3所示；

仿真3，采用本发明对信号源2定位时，仿真目标角度均方误差随信噪比SNR变化，仿真结果如图4所示。

从图3和图4可以看出，目标角度估计的均方误差随信噪比SNR增加而降低，且信噪比为5dB时目标角度估计的均方误差就可达到0.1°以下，定位精度高，表明本发明是实际可行的。

综上，本发明可实现对任意阵列单基地MIMO雷达波达方向的估计，且定位精度高。

Claims (4)

1.一种任意阵列流形的单基地MIMO雷达波达方向估计方法，包括如下步骤：

1)根据阵列流形写出单基地MIMO雷达的导向矢量A(θ)，其中θ为方位角；

2)对任意形状的MIMO雷达的导向矢量A(θ)采用Jacobi-Anger展开，并根据阵列的最大孔径进行维数截断，获得关于e^jθ的范得蒙导向矢量D(θ)，对D(θ)进行空间平滑处理，处理后的导向矢量记为D_S(θ)，利用最小二乘法得出从MIMO雷达导向矢量A(θ)到平滑处理后的导向矢量D_S(θ)的转换矩阵G；

3)对雷达的接收数据与发射波形进行匹配滤波，将匹配滤波结果记为X(t)，再利用滤波结果X(t)形成自相关矩阵R_x；

4)对自相关矩阵R_x进行特征分解，得到一系列特征值λ_m及其对应的特征向量e_m，m＝1,2,…,N²,其中N为阵元数；将特征值λ_m进行从小到大排序，依次选取N²-P个小特征值，P为同一距离门内的目标个数，取这些特征值对应的特征向量形成噪声子空间 ${\hat{E}}_n=[e_{P+1},\·\·\·e_m\·\·\·,e_{N^2}];$

5)利用噪声子空间形成MUSIC空间零谱函数： $f\left(\mathrm{\θ}\right)={\left|\right|{{\hat{E}}_n}^{H}A\left(\mathrm{\θ}\right)\left|\right|}^2=A^H\left(\mathrm{\θ}\right){\hat{E}}_n{{\hat{E}}_n}^{H}A\left(\mathrm{\θ}\right),$ 其中||||²表示2-范数，表示矩阵的共轭转置；

6)把步骤5)中的A(θ)用步骤2)中的GD_s(θ)代替，并令z＝e^jθ，则f(θ)可转化为：令f(z)＝0，采用多项式求根的方法求得P个方位角的复指数形式z_k，其中k＝1,2…P；

7)利用接收角的复指数形式z_k求得第k个接收方位角θ_k：θ_k＝angle(z_k)，其中angle()表示求相位角。

2.根据权利要求1所述的方法，其中，步骤2)所述的对MIMO雷达的导向矢量A(θ)流形展开，获得范得蒙导向矢量D(θ)，按如下步骤进行：

令R＝max(r_n)，r_n为第n个阵元距坐标原点的距离，n＝1,2,…N；选择范得蒙导向矢量D(θ)的维数M＞8kR，k＝2π/λ，λ为波长；

忽略截断误差，根据所述维数M，得到范得蒙导向矢量D(θ)：

$D\left(\mathrm{\θ}\right)={[e^{\frac{-j(M-1)\mathrm{\θ}}{2}},\·\·\·,e^{\frac{j(M-1)\mathrm{\θ}}{2}}]}^T,$

其中，[]^T表示矩阵转置，j表示虚数。

3.根据权利要求1所述的方法，其中，步骤2)所述的对D(θ)进行空间平滑处理，按如下步骤进行：

将范得蒙导向矢量D(θ)划分为相互重叠的L个子阵，对应每个子阵的阵元数为Q＝M-L+1，其中，M为范得蒙导向矢量D(θ)的维数；

根据所述阵元数Q构造平滑矩阵：F_k＝[0_Q×(k-1)|I_Q|0_Q×(L-k)]，其中，k＝1,2,…L，0_Q×(k-1)表示Q×(k-1)阶零阵，I_Q表示Q阶单位矩阵，0_Q×(L-k)表示Q×(L-k)阶零阵；

将范得蒙形式的导向矢量D(θ)与平滑矩阵F_k相乘，得到空间平滑后的导向矢量D_S(θ)为：

$D_S\left(\mathrm{\θ}\right)=\frac{1}{L}\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{F}_{k}D\left(\mathrm{\θ}\right).$

4.根据权利要求1所述的方法，其中，步骤2)所述的利用最小二乘法得出从MIMO雷达导向矢量A(θ)到平滑处理后的导向矢量D_S(θ)的转换矩阵G，按如下步骤进行：

在0-360度范围内等间隔的取N个角度，即：θ₁,θ₂,...,θ_N，N为阵元数；

根据θ₁,θ₂,...,θ_N，分别形成最小二乘导向矩阵A^c和最小二乘平滑矩阵D_c：

A^c＝[A(θ₁),A(θ₂),…,A(θ_N)]，

D_c＝[D_s(θ₁),D_s(θ₂),…,D_s(θ_N)]；

根据最小二乘导向矩阵A^c和最小二乘平滑矩阵D_c，利用最小二乘法得到转换矩阵G：

G＝A^cD_c ^H(D_cD_c ^H)^-1，

其中，D_c ^H表示矩阵D_c的共轭转置，()^-1表示矩阵求逆。