CN105353367A - 一种双基地mimo雷达空间机动目标跟踪方法 - Google Patents

Abstract

本发明公开了一种双基地MIMO雷达的空间机动目标跟踪方法，能够实现对空间机动目标的三维位置、速度等运动参数的准确估计与预测；该方法将双基地MIMO雷达均匀圆阵体制应用于空间机动目标跟踪，在双基地MIMO雷达均匀圆阵角度多普勒频率估计的基础上，基于目标状态空间，在发射均匀圆阵与接收均匀圆阵的测量直角坐标系下，根据测元与目标运动状态变量之间的函数关系，进而建立测量方程，并采用适用于机动目标跟踪的“当前”统计模型建立描述目标运动参数随时间演化的状态方程，利用适用于非线性状态滤波的无迹卡尔曼滤波算法实时估计目标的运动参数，并可对目标运动参数进行预测，实现对空间机动目标的跟踪。

Description

一种双基地MIMO雷达空间机动目标跟踪方法

技术领域

本发明属于雷达数据处理领域，具体涉及一种双基地MIMO雷达空间机动目标跟踪方法。

背景技术

双基地雷达是采用接收机和发射机系统分置结构实现的。这种结构的主要特点是：发射机置于后方，接收机置于前方，可以接近目标进行隐蔽侦查，同时避免了雷达电磁波双程传播带来的威力损失，提高了目标的信噪比。传统的双基地米波雷达利用目标相对于接收端的角度及距离和对空间目标进行定位，由于双基地雷达的接收机与发射机难以满足精确的时间同步，且其角度分辨率与测距精度均比较低，导致对目标定位的精度低。采用MIMO(多输入多输出)技术的双基地雷达，即双基地MIMO雷达，可以在接收端获得发射端的角度信息，在不需要时间同步和目标距离的情况下，能够对目标进行精确定位，获得目标的位置坐标。

西安电子科技大学申请的专利技术“多输入多输出雷达系统目标定位方法”(申请号200810150754.X，公开号CN101349748A)中公开了一种多输入多输出MIMO雷达目标定位的方法。该方法只能确定雷达目标的二维位置坐标，无法估计空间目标的三维位置坐标。

哈尔滨工程大学申请的专利技术“双基地多输入多输出雷达多目标定位方法”(申请号201110001351.0，公开号CN102135617A)中公开了一种双基地多输入多输出雷达多目标定位方法。该方法也仅适用于二维目标定位，无法估计空间目标的三维位置坐标。

西安电子科技大学申请的专利技术“双基地米波雷达目标三维精确定位方法”(专利申请号：201218001807.9)中公开了一种双基地多输入多输出雷达的多目标三维定位方法，将双基地雷达配置为发射阵为均匀圆阵、接收阵均匀线阵，估计出目标相对于均匀圆阵的方位角与俯仰角、目标相对于均匀线阵的接收角，利用这三个角度以及基线长度通过几何解算方法得到目标的三维位置坐标。该方法虽然可用于空间目标的三维定位，但仅利用了三个角度信息，且采用的是几何解算方法，当目标距离较远时，测角误差放大效应会使目标定位精度明显下降，无法对空间机动目标进行高精度跟踪。

西安电子科技大学申请的专利技术“双基地MIMO雷达均匀圆阵角度多普勒频率估计方法”(专利申请号：201510487104.4)中提出了一种双基地MIMO雷达均匀圆阵角度多普勒频率估计方法，将双基地雷达配置为发射阵、接收阵均为均匀圆阵，可同时估计出多个目标相对于接收阵与发射阵的方位角、俯仰角以及归一化多普勒频率。但该方法未涉及目标运动状态估计及预测跟踪问题。

发明内容

有鉴于此，本发明提供了一种双基地MIMO雷达空间机动目标的跟踪方法，其目的基于现有技术的不足，综合利用目标相对于双基地MIMO雷达的角度、多普勒频率估计信息，实现对空间机动目标的三维位置、速度等运动参数的准确估计与预测。

实施本发明的具体方案如下：

一种双基地MIMO雷达空间机动目标跟踪方法，包括以下步骤：

步骤1，分别将双基地MIMO雷达的发射机配置为M个阵元的均匀圆阵、接收机配置为N个阵元的均匀圆阵，并使发射机中M个阵元发射相互正交的波形信号；其中，M表示发射机阵元个数，N表示接收机阵元个数，且M、N均为自然数；

步骤2，利用所述发射机中M个阵元发射相互正交的波形信号，接收机中的N个阵元分别接收所述发射机中M个阵元发射的相互正交的波形信号，并进行匹配滤波，获得匹配滤波后的NM×1维雷达回波信号s和L次快拍积累得到的NM×L维雷达回波信号矩阵S，进而获得矩阵S中第n个接收阵元的M×L维切片矩阵形式S_n，其中，L表示快拍次数，n为自然数，n＝1，2，...，N；

步骤3，根据N个接收阵元的M×L维切片矩阵S₁～S_N，利用平行因子算法分别获得发射方向向量的估计接收方向向量的估计和归一化多普勒频率方向向量的估计

步骤4，根据所述发射方向向量的估计获得目标相对于发射机的方位角估计值和俯仰角估计值根据所述接收方向向量的估计获得目标相对于接收机的方位角估计值和俯仰角估计值

步骤5，根据所述归一化多普勒频率方向向量的估计利用最小二乘算法获得目标的归一化多普勒频率估计值

步骤6，利用地心直角坐标系建立目标的运动状态空间，基于目标的运动状态空间，采用CS模型建立目标状态方程其中，为目标运动状态变量，包括目标三维位置[x，y，z]、速度加速度k表示目标状态的采样时刻序号，k＝1，2，3，...，X(k-1)和X(k)分别表示第k-1个和第k个采样时刻的目标状态，Φ表示状态转移矩阵，U表示控制矩阵，表示第k个采样时刻的目标机动加速度均值，W(k)表示第k个采样时刻的目标状态噪声向量；

步骤7，基于目标状态空间，在发射均匀圆阵与接收均匀圆阵的测量直角坐标系下，根据测元与目标运动状态变量X之间的函数关系，建立测量方程Y(k)＝h(X(k))+v(k)，其中，v(k)表示各测元的测量噪声向量；Y(k)表示第k个采样时刻的目标测元，h(X(k))表示第k个采样时刻的目标测元与第k个采样时刻的目标运动状态之间的映射关系；

步骤8，基于测量方程，根据初始两个采样时刻获取的目标测元Y(1)，Y(2)，计算X(1)＝[x(1)，y(1)，z(1)]^T和X(2)＝[x(2)，y(2)，z(2)]^T，进而计算目标运动状态变量的初值和对应的协方差矩阵P(0)，得到目标运动状态的滤波初值；

步骤9，获得滤波初值后，采用UKF算法，利用第k个采样时刻的测元Y(k)、目标在第k-1个采样时刻的状态估计与协方差矩阵P(k-1)获得目标在第k个采样时刻的状态估计与协方差矩阵P(k)，依次获得目标在各采样时刻的运动状态，完成对空间机动目标的跟踪。

进一步地，在L次快拍积累中第n个接收阵元的M×L维切片矩阵S_n，其表达式为：

S_n＝A_R(n)A_TB^T+W_n，n＝1，2，...，N

其中，A_R(n)表示接收方向向量A_R的第n个元素，W_n表示第n个接收阵元L次快拍积累中的噪声矩阵；A_T表示发射方向向量，B表示归一化多普勒频率方向向量，其中，ρ为目标雷达波反射系数，f_d为目标的归一化多普勒频率。

进一步地，在步骤3中，所述发射方向向量的估计所述接收方向向量的估计和所述归一化多普勒频率方向向量的估计的具体子步骤为：

a)计算第m个发射阵元的第n个接收阵元的第l次快拍的平行因子三线性模型形式s_m，n，l，其具体表达式为

s_m，n，l＝A_T(m)A_R(n)B(l)+w_m，n，l，m∈{1，…，M}，n∈{1，…，N}，l∈{1，…，L}

其中，A_T(m)表示发射方向向量A_T的第m个元素，A_R(n)表示接收方向向量A_R的第n个元素，B(l)表示归一化多普勒频率方向向量B的第l个元素，w_m，n，l表示s_m，n，l的三维噪声数据集；

b)根据第m个发射阵元的第n接收阵元的第l次快拍的平行因子三线性模型形式s_m，n，l的对称性，对s_m，n，l的第二维和第三维分别进行切片，分别获得第m个发射阵元的L×N维切片矩阵形式F_m、第l次快拍的N×M维切片矩阵形式Z_l，其表达式分别为：

F_m＝BA_T(m)A_R ^T+W_m，m∈{1，2，...，M}

Z_l＝A_RB(l)A_T ^T+W_l，l∈{1，2，...，L}

其中，W_m表示第m个发射阵元噪声，W_l表示第l次快拍噪声；

c)将N个接收阵元的M×L维切片矩阵形式S₁～S_N按列形成NM×L维矩阵S，将M个发射阵元的L×N维切片矩阵形式F₁～F_M按列形成ML×N维矩阵F，将L次快拍的N×M维切片矩阵形式Z₁～Z_L按列形成NL×M维矩阵Z，具体地，

S＝[A_RοA_T]B^T+W_S＝[S₁，...，S_N]^T

F＝[A_TοB]A_R ^T+W_F＝[F₁，…，F_M]^T

Z＝[BοA_R]A_T ^T+W_Z＝[Z₁，…，Z_L]^T

其中，ο表示Khatri-Rao积，W_S表示NM×L维矩阵S的噪声矩阵，W_F表示ML×N维矩阵F的噪声矩阵，W_Z表示NL×M维矩阵Z的噪声矩阵；

d)根据NM×L维矩阵S、ML×N维矩阵F和NL×M维矩阵Z，利用平行因子算法得到发射方向向量的估计接收方向向量的估计和归一化多普勒频率方向向量的估计具体地，

其中，[·]⁺表示取伪逆。

进一步地，在步骤8中，完成滤波初值的计算，具体为：

a)利用Y(1)中的4个角度计算目标在第1个采样时刻的位置参数X(1)＝[x(1)，y(1)，z(1)]^T；利用Y(2)中的4个角度计算目标在第2个采样时刻的位置参数X(2)＝[x(2)，y(2)，z(2)]^T；

根据4个角度θ_t，φ_t，θ_r，φ_r解算目标位置x，y，z的具体过程如下：

θ_t表示目标相对于发射机的方位角，φ_t表示目标相对于发射机的俯仰角，θ_r表示目标相对于接收机的方位角，φ_r表示目标相对于接收机的俯仰角；

i)将4个角度从发射阵与接收阵测量直角坐标系转至测量地心直角系为：

$\left[\begin{array}{c}{l}_t\\ m_t\\ n_t\end{array}\right]={\mathrm{\Ω}}_t^T\left[\begin{array}{c}{\mathrm{sin\φ}}_t{\mathrm{cos\θ}}_t\\ {\mathrm{cos\φ}}_t\\ {\mathrm{sin\φ}}_t{\mathrm{sin\θ}}_t\end{array}\right]$

${\mathrm{\&theta;}}_t^{\&prime;}=\arctan \left(\frac{n_t}{l_t}\right)+\left\{\begin{array}{cc}0,& l_t>0,n_t>0\\ \mathrm{\&pi;},& l_t<0\\ 2\mathrm{\&pi;},& l_t>0,n_t<0\end{array}\right.$

${\mathrm{\&phi;}}_t^{\&prime;}=\frac{\mathrm{\&pi;}}{2}-\arctan \left(\frac{m_t}{\sqrt{{\left(l_t\right)}^2+{\left(n_t\right)}^2}}\right)$

$\left[\begin{array}{c}{l}_r\\ m_r\\ n_r\end{array}\right]={\mathrm{\&Omega;}}_r^T\left[\begin{array}{c}{\mathrm{sin\&phi;}}_r{\mathrm{cos\&theta;}}_r\\ {\mathrm{cos\&phi;}}_r\\ {\mathrm{sin\&phi;}}_r{\mathrm{sin\&theta;}}_r\end{array}\right]$

${\mathrm{\&theta;}}_r^{\&prime;}=\arctan \left(\frac{n_r}{l_r}\right)+\left\{\begin{array}{cc}0,& l_r>0,n_r>0\\ \mathrm{\&pi;},& l_r<0\\ 2\mathrm{\&pi;},& l_r>0,n_r<0\end{array}\right.$

${\mathrm{\&phi;}}_r^{\&prime;}=\frac{\mathrm{\&pi;}}{2}-\arctan \left(\frac{m_r}{\sqrt{{\left(l_r\right)}^2+{\left(n_r\right)}^2}}\right)$

其中，Ω_t ^T表示地心直角系到发射阵测量直角系的转换矩阵的转置，Ω_r ^T表示地心直角系到接收阵测量直角系的转换矩阵的转置，(l_t，m_t，n_t)为发射阵中心至目标连线的方向向量，θ′_t、φ′_t分别表示测量地心直角系下目标相对于发射机的方位角和俯仰角；(l_r，m_r，n_r)为接收阵中心至目标连线的方向向量，θ′_r、φ′_r分别表示测量地心直角系下目标相对于接收机的方位角和俯仰角；

ii)采用水平投影法解算目标在地心直角系下的位置

$\left\{\begin{array}{c}x\left(k\right)=X_{t0}+{\mathrm{\&Delta;X}}_1\\ y\left(k\right)=Y_{t0}+\frac{{\mathrm{\&Delta;X}}_1}{{\mathrm{cos\&theta;}}_t^{\&prime;}}{\mathrm{cot\&phi;}}_t^{\&prime;}\\ z\left(k\right)=Z_{t0}+{\mathrm{\&Delta;X}}_1{\mathrm{tan\&theta;}}_t^{\&prime;}\end{array}\right.$

其中

${\mathrm{\&Delta;X}}_1=\frac{(X_{t0}-X_{r0}){\mathrm{tan\&theta;}}_r^{\&prime;}-(Z_{t0}-Z_{r0})}{{\mathrm{tan\&theta;}}_t^{\&prime;}-{\mathrm{tan\&theta;}}_r^{\&prime;}}$

其中，(X_t0，Y_t0，Z_t0)为发射阵中心在地心直角坐标系下的位置；ΔX₁为发射阵中心至目标在X轴方向的距离，X_r0为接收阵中心在地心直角坐标系X轴方向的位置，Z_r0为接收阵中心在地心直角坐标系Z轴方向的位置；当k＝1时，获得采样时刻的位置参数X(1)＝[x(1)，y(1)，z(1)]^T，当k＝2时，获得采样时刻的位置参数X(2)＝[x(2)，y(2)，z(2)]^T；k取自然数；

b)对前两个采样时刻的位置参数进行差分计算目标在第2个采样时刻的近似速度

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}\left(2\right)=\frac{x\left(2\right)-x\left(1\right)}{T}\\ \stackrel{\&CenterDot;}{y}\left(2\right)=\frac{y\left(2\right)-y\left(2\right)}{T}\\ \stackrel{\&CenterDot;}{z}\left(2\right)=\frac{z\left(2\right)-z\left(1\right)}{T}\end{array}\right.$

其中，T表示采样周期；

c)根据上述计算结果构造滤波初值，

$\hat{X}\left(0\right)={\&lsqb;x\left(2\right),\stackrel{\&CenterDot;}{x}\left(2\right),0,y\left(2\right),\stackrel{\&CenterDot;}{y}\left(2\right),0,z\left(2\right),\stackrel{\&CenterDot;}{z}\left(2\right),0\&rsqb;}^T$

P(0)＝10⁶I₉

其中，I₉表示9维的单位矩阵。

进一步地，在步骤9中，所述完成对空间机动目标跟踪的具体步骤为：

a)状态向量采样

对状态估计向量进行sigma采样，k取自然数，χ_p(k-1)表示第p个状态采样点，p＝0，…，2n_x，n_x＝9；采样点为：

其中，κ为比例参数，表示矩阵的第i列，i＝1，2，…，n_x

b)时间更新

利用CS模型对目标状态采样点进行一步预测得到χ_p(k|k-1)，一步预测状态转移矩阵为Φ_a，根据目标测元与目标运动状态之间的映射函数h(·)计算观测一步预测的采样点Υ_p(k|k-1)，再对其加权求和，计算观测量的一步预测

χ_p(k|k-1)表示第p个状态采样点从k-1时刻到k时刻的预测值；Υ_p(k|k-1)表示第p个观测一步预测的采样点；表示第p个采样点的均值权重；为第k个采样时刻的观测量的一步预测；

χ_p(k/k-1)＝Φ_aχ_p(k-1)

Υ_p(k|k-1)＝h(χ_p(k|k-1))

c)状态更新

$\hat{X}\left(k\right)=\stackrel{\&OverBar;}{X}\left(k\right|k-1)+G\left(k\right)\left(Y\right(k)-\stackrel{\&OverBar;}{Y}(k|k-1\left)\right)$

P(k)＝P(k|k-1)-G(k)P_yy(k|k-1)G^T(k)

G(k)表示第k个采样时刻目标状态更新增益矩阵；P(k|k-1)表示目标状态协方差矩阵的一步预测；P_yy(k|k-1)表示观测协方差矩阵的一步预测；和分别表示第k-1到第k个采样时刻的目标状态一步预测和测元一步预测的均值；G^T(k)表示第k个采样时刻的目标状态更新增益矩阵的转置。

有益效果：

(1)本发明采用发射机、接收机均为均匀圆阵的配置，不但能够估计目标的四个角度和一个多普勒频率信息，还能够对空间目标进行三维定位。

(2)在进行目标运动状态参数解算时，传统的几何解算方法只能利用角度测元，无法利用多普勒频率测元，造成了测元信息浪费，而本发明采用了基于CS模型的UKF算法估计目标的运动状态参数，一方面提高了对测量噪声的抑制能力，另一方面能够融合利用所有测元信息，这两方面共同提高了目标运动状态参数的解算精度。

(3)本发明采用了适用于机动目标跟踪的CS模型，基于此可以对空间机动目标的运动状态进行一定时间的预测，实现对空间机动目标的有效稳定跟踪。

附图说明

图1为本发明的一种双基地MIMO雷达空间机动目标跟踪方法实现流程图

图2为本发明的双基地MIMO雷达配置示意图。

图3为用本发明方法仿真信噪比20dB情况下某临近空间运动目标相对于发射机的俯仰角和方位角的估计误差图。

图4为用本发明方法仿真信噪比20dB情况下某临近空间运动目标相对于发射机的俯仰角和方位角的估计误差图。

图5为用本发明方法仿真信噪比20dB情况下某临近空间运动目标多普勒频率的估计误差图。

图6为用本发明方法仿真信噪比为20dB时的目标位置滤波误差图。

图7为用本发明方法仿真信噪比为20dB时的目标速度滤波误差图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

图1为本发明的一种双基地MIMO雷达空间机动目标跟踪方法实现流程图，实现本发明的具体步骤如下：

步骤1，分别将双基地MIMO雷达的发射机配置为M个阵元的均匀圆阵、接收机配置为N个阵元的均匀圆阵，并使发射机中M个阵元发射相互正交的波形信号；其中，M表示发射机阵元个数，N表示接收机阵元个数，且M、N均为自然数；

图2为本发明的双基地MIMO雷达配置示意图。坐标系采用地心直角坐标系O_G-X_GY_GZ_G，原点O_G为地球中心，O_GX_G轴的正向为沿起始天文子午面与地球赤道平面的交线指向外方向，O_GZ_G轴沿地球自转轴指向北极，O_GY_G轴与O_GX_G、O_GZ_G轴构成右手直角坐标系。其中，图2中O_t为发射圆阵的中心点，O_r为接收圆阵的中心点。

发射圆阵与接收圆阵均布设在当地水平面内。设发射圆阵的圆点O_t的大地坐标(经度L_t，纬度B_t，高程h_t)为(L_t，B_t，h_t)，接收圆阵的圆点O_r的大地坐标(经度L_r，纬度B_r，高程h_r)为(L_r，B_r，h_r)。发射圆阵半径为r_t，其阵元个数M＝2floor(2πr_t/λ)+1，floor(·)表示向下取整运算，λ表示发射波的波长；接收圆阵半径为r_r，其阵元个数N＝2floor(2πr_r/λ)+1。

假设空间有一目标P，在地心系下的位置坐标为(x，y，z)、速度坐标为发射圆阵的圆心与目标的连线O_tP与发射圆阵平面的夹角为俯仰角，记为φ_t，O_tP在发射圆阵平面的投影与当地正北方向的夹角为方位角，记为θ_t；接收圆阵的圆心至目标连线O_rP与接收圆阵平面的夹角为俯仰角，记为φ_r，O_rP在接收圆阵平面的投影与当地正北方向的夹角为方位角，记为θ_r；由目标运动引起的接收信号多普勒频率为f_d。

步骤2，利用所述发射机中M个阵元发射相互正交的波形信号，接收机中的N个阵元分别接收所述发射机中M个阵元发射的相互正交的波形信号，并进行匹配滤波，获得匹配滤波后的NM×1维雷达回波信号s和L次快拍积累得到的NM×L维雷达回波信号矩阵S，进而获得矩阵S中第n个接收阵元的M×L维切片矩阵形式S_n，其中，L表示快拍次数，n为自然数，n＝1，2，...，N；

具体地，利用发射机中M个阵元发射相互正交的波形信号，接收机中的N个阵元分别接收该发射机中M个阵元的发射信号，并进行匹配滤波，依次获得匹配滤波后的NM×1维雷达回波信号s，考察第l个脉冲匹配滤波后的回波信号s(l)，其表达式为

$s\left(l\right)={\mathrm{\&rho;e}}^{j2\mathrm{\&pi;}(l-1)f_d}\&lsqb;a_r({\mathrm{\&theta;}}_r,{\mathrm{\&phi;}}_r)\&CircleTimes;a_t({\mathrm{\&theta;}}_t,{\mathrm{\&phi;}}_t)\&rsqb;+\tilde{n}\left(l\right)$

其中，ρ为目标雷达波反射系数，假定目标对不同正交发射信号的反射系数在相同脉冲内是完全相关的，在不同脉冲间是相互独立的，a_r(θ_r，φ_r)为目标接收方向向量，a_t(θ_t，φ_t)为目标发射方向向量，表示Kronecker积，为第l个脉冲匹配滤波后的接收噪声，l＝1，2，…，L；

$a_r({\mathrm{\&theta;}}_r,{\mathrm{\&phi;}}_r)=\left[\begin{array}{c}{e}^{j\frac{2{\mathrm{\&pi;r}}_r}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_r\cos ({\mathrm{\&theta;}}_r-{\mathrm{\&beta;}}_{r1})}\\ .\\ .\\ .\\ e^{j\frac{2{\mathrm{\&pi;r}}_r}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_r\cos ({\mathrm{\&theta;}}_r-{\mathrm{\&beta;}}_{rn})}\\ .\\ .\\ .\\ e^{j\frac{2{\mathrm{\&pi;r}}_r}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_r\cos ({\mathrm{\&theta;}}_r-{\mathrm{\&beta;}}_{rN})}\end{array}\right]=\left[\begin{array}{c}{a}_{r1}({\mathrm{\&theta;}}_r,{\mathrm{\&phi;}}_r)\\ .\\ .\\ .\\ a_{rn}({\mathrm{\&theta;}}_r,{\mathrm{\&phi;}}_r)\\ .\\ .\\ .\\ a_{rN}({\mathrm{\&theta;}}_r,{\mathrm{\&phi;}}_r)\end{array}\right]\stackrel{\mathrm{\&Delta;}}{=}{A}_R$

其中，β_rn为第n个接收阵元的方位角，β_rn＝2π(n-1)/N，n＝1，2，…，N，符号表示“简记为”；a_rn(θ_r，φ_r)表示第n个接收阵元的接收方向因子。

$a_t({\mathrm{\&theta;}}_t,{\mathrm{\&phi;}}_t)=\left[\begin{array}{c}{e}^{j\frac{2{\mathrm{\&pi;r}}_t}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_t\cos ({\mathrm{\&theta;}}_t-{\mathrm{\&beta;}}_{t1})}\\ .\\ .\\ .\\ e^{j\frac{2{\mathrm{\&pi;r}}_t}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_t\cos ({\mathrm{\&theta;}}_t-{\mathrm{\&beta;}}_{tn})}\\ .\\ .\\ .\\ e^{j\frac{2{\mathrm{\&pi;r}}_t}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_t\cos ({\mathrm{\&theta;}}_t-{\mathrm{\&beta;}}_{tN})}\end{array}\right]=\left[\begin{array}{c}{a}_{t1}({\mathrm{\&theta;}}_t,{\mathrm{\&phi;}}_t)\\ .\\ .\\ .\\ a_{tm}({\mathrm{\&theta;}}_t,{\mathrm{\&phi;}}_t)\\ .\\ .\\ .\\ a_{tM}({\mathrm{\&theta;}}_t,{\mathrm{\&phi;}}_t)\end{array}\right]\stackrel{\mathrm{\&Delta;}}{=}{A}_T$

其中，β_tm为第m个发射阵元的方位角，β_tm＝2π(m-1)/M，m＝1，2，…，M；α_tm(θ_t，φ_t)表示第m个发射阵元的的发射方向因子。

经过L次快拍积累得到NM×L维矩阵S，在L次快拍内认为目标运动状态参数不变，S的表达式为

$S=\&lsqb;s\left(1\right),s\left(2\right),\mathrm{...},s\left(L\right)\&rsqb;=\&lsqb;a_r({\mathrm{\&theta;}}_r,{\mathrm{\&phi;}}_r)\&CircleTimes;a_t({\mathrm{\&theta;}}_t,{\mathrm{\&phi;}}_t)\&rsqb;B^T+\tilde{N}$

其中， $B={\&lsqb;\mathrm{\&rho;},{\mathrm{\&rho;e}}^{j2{\mathrm{\&pi;f}}_d},\mathrm{...},{\mathrm{\&rho;e}}^{j2\mathrm{\&pi;}(L-1)f_d}\&rsqb;}^T,\tilde{N}=\&lsqb;\tilde{n}\left(1\right),\tilde{n}\left(2\right),\mathrm{...},\tilde{n}\left(L\right)\&rsqb;,$ 上标T表示转置。表示由形成的矩阵，B表示归一化多普勒频率方向向量。

进而得到L次快拍积累中第n个接收阵元的M×L维切片矩阵S_n，其表达式为：

S_n＝A_R(n)A_TB^T+W_n，n＝1，2，...，N

其中，A_R(n)表示向量A_R的第n个元素，W_n表示第n个接收阵元L次快拍积累中的噪声矩阵；A_T表示发射方向向量，A_R表示接收方向向量。

步骤3，根据N个接收阵元的M×L维切片矩阵形式S₁～S_N，利用平行因子算法分别获得发射方向向量的估计接收方向向量的估计和多普勒频率方向向量的估计

步骤3的具体子步骤为：

3a)第n个接收阵元的M×L维切片矩阵S_n，根据N个接收阵元的M×L维切片矩阵S₁～S_N，获得M×N×L的三维数据集，进而得到第m个发射阵元的第n个接收阵元的第l次快拍的平行因子三线性模型形式s_m，n，l，其具体表达式为

s_m，n，l＝A_T(m)A_R(n)B(l)+w_m，n，l，m∈{1，…，M}，n∈{1，…，N}，l∈{1，…，L}，

其中，A_T(m)表示发射方向向量A_T的第m个元素，A_R(n)表示接收方向向量A_R的第n个元素，B(l)表示多普勒频率方向向量B的第l个元素，w_m，n，l表示三维噪声数据集；

3b)根据s_m，n，l，分别得到第m个发射阵元的L×N维切片矩阵形式F_m、第l次快拍的N×M维切片矩阵形式Z_l，进而分别获得M个发射阵元的L×N维切片矩阵形式F₁～F_M，L次快拍的N×M维切片矩阵形式Z₁～Z_L。

具体地，根据第m个发射阵元的第n接收阵元的第l次快拍的平行因子三线性模型形式s_m，n，l的对称性，对其第二维和第三维分别进行切片，分别获得第m个发射阵元的L×N维切片矩阵形式F_m、第l次快拍的N×M维切片矩阵形式Z_l，其表达式分别为：

F_m＝BA_T(m)A_R ^T+W_m，m∈{1，2，...，M}

Z_l＝A_RB(l)A_T ^T+W_l，l∈{1，2，...，L}

3c)将N个接收阵元的M×L维切片矩阵形式S₁～S_N按列形成NM×L维矩阵S，将M个发射阵元的L×N维切片矩阵形式F₁～F_M按列形成ML×N维矩阵F，将L次快拍的N×M维切片矩阵形式Z₁～Z_L按列形成NL×M维矩阵Z，具体地，

S＝[A_RοA_T]B^T+W_S＝[S₁，...，S_N]^T

F＝[A_TοB]A_R ^T+W_F＝[F₁，…，F_M]^T

Z＝[BοA_R]A_T ^T+W_P＝[Z₁，...，Z_L]^T

其中，ο表示Khatri-Rao积，W_S表示NM×L维矩阵S的噪声矩阵，W_F表示ML×N维矩阵F的噪声矩阵，W_P表示NL×M维矩阵P的噪声矩阵。

3d)根据NM×L维矩阵S、ML×N维矩阵F和NL×M维矩阵P，利用平行因子算法得到发射方向向量的估计接收方向向量的估计和多普勒频率方向向量的估计具体地，

其中，[·]⁺表示取伪逆。

步骤4，根据所述发射方向向量的估计获得目标相对于发射机的方位角估计值和俯仰角估计值根据所述接收方向向量的估计获得目标相对于接收机的方位角估计值和俯仰角估计值

具体地，发射方向向量的估计为接收方向向量的估计的为其表达式分别为：

${\hat{a}}_t({\mathrm{\&phi;}}_t,{\mathrm{\&theta;}}_t)={\&lsqb;e^{j\frac{2{\mathrm{\&pi;r}}_t}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_t\&CenterDot;\cos ({\mathrm{\&theta;}}_t-{\mathrm{\&beta;}}_{t1})},\mathrm{...},e^{j\frac{2{\mathrm{\&pi;r}}_t}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_t\&CenterDot;\cos ({\mathrm{\&theta;}}_t-{\mathrm{\&beta;}}_{tm})},\mathrm{...},e^{j\frac{2{\mathrm{\&pi;r}}_t}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_t\&CenterDot;\cos ({\mathrm{\&theta;}}_t-{\mathrm{\&beta;}}_{tM})}\&rsqb;}^T$

${\hat{a}}_r({\mathrm{\&phi;}}_r,{\mathrm{\&theta;}}_r)={\&lsqb;e^{j\frac{2{\mathrm{\&pi;r}}_r}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_r\&CenterDot;\cos ({\mathrm{\&theta;}}_r-{\mathrm{\&beta;}}_{r1})},\mathrm{...},e^{j\frac{2{\mathrm{\&pi;r}}_r}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_r\&CenterDot;\cos ({\mathrm{\&theta;}}_r-{\mathrm{\&beta;}}_{rn})},\mathrm{...},e^{j\frac{2{\mathrm{\&pi;r}}_r}{\mathrm{\&lambda;}}{\mathrm{sin\&phi;}}_r\&CenterDot;\cos ({\mathrm{\&theta;}}_r-{\mathrm{\&beta;}}_{rN})}\&rsqb;}^T$

一般，β_r1＝0，β_t1＝0，中的每一项都除以第一项然后去掉其第一项，获得发射新向量a₁，再取a₁对数的虚部得到a₁＝Im(lna₁)；中的每一项都除以第一项然后去掉其第一项，得到接收新向量a₂，再取a₂对数的虚部得到a′₂＝Im(lna₂)，a′₁和a′₂的表达式分别为：

$a_1^{\&prime;}=\left[\begin{array}{c}{\mathrm{\&xi;}}_t{\mathrm{sin\&phi;}}_t{\mathrm{cos\&theta;}}_t({\mathrm{cos\&beta;}}_{t2}-1)+{\mathrm{\&xi;}}_t{\mathrm{sin\&phi;}}_t{\mathrm{sin\&theta;}}_t{\mathrm{sin\&beta;}}_{t2}\\ {\mathrm{\&xi;}}_t{\mathrm{sin\&phi;}}_t{\mathrm{cos\&theta;}}_t({\mathrm{cos\&beta;}}_{{\mathrm{tm}}^{\&prime;}}-1)+{\mathrm{\&xi;}}_t{\mathrm{sin\&phi;}}_t{\mathrm{sin\&theta;}}_t{\mathrm{sin\&beta;}}_{{\mathrm{tm}}^{\&prime;}}\\ {\mathrm{\&xi;}}_t{\mathrm{sin\&phi;}}_t{\mathrm{cos\&theta;}}_t({\mathrm{cos\&beta;}}_{tM}-1)+{\mathrm{\&xi;}}_t{\mathrm{sin\&phi;}}_t{\mathrm{sin\&theta;}}_t{\mathrm{sin\&beta;}}_{tM}\end{array}\right]$

$a_2^{\&prime;}=\left[\begin{array}{c}{\mathrm{\&xi;}}_r{\mathrm{sin\&phi;}}_r{\mathrm{cos\&theta;}}_r({\mathrm{cos\&beta;}}_{r2}-1)+{\mathrm{\&xi;}}_r{\mathrm{sin\&phi;}}_r{\mathrm{sin\&theta;}}_r{\mathrm{sin\&beta;}}_{r2}\\ .\\ .\\ .\\ {\mathrm{\&xi;}}_r{\mathrm{sin\&phi;}}_r{\mathrm{cos\&theta;}}_r({\mathrm{cos\&beta;}}_{{\mathrm{tn}}^{\&prime;}}-1)+{\mathrm{\&xi;}}_r{\mathrm{sin\&phi;}}_r{\mathrm{sin\&theta;}}_r{\mathrm{sin\&beta;}}_{{\mathrm{tn}}^{\&prime;}}\\ .\\ .\\ .\\ {\mathrm{\&xi;}}_r{\mathrm{sin\&phi;}}_r{\mathrm{cos\&theta;}}_r({\mathrm{cos\&beta;}}_{tN}-1)+{\mathrm{\&xi;}}_r{\mathrm{sin\&phi;}}_r{\mathrm{sin\&theta;}}_r{\mathrm{sin\&beta;}}_{tN}\end{array}\right]$

其中，ξ_t＝2πr_t/λ，ξ_r＝2πr_r/λ；

a′₁中的第m项除以(cosβ_tm′-1)，m′∈{2，3，...，M}，得到β₁，

a′₂中的第n项除以(cosβ_rn′-1)，n′∈{2，3，...，N}，得到β₂，则β₁和β₂的表达式分别为：

${\mathrm{\&beta;}}_1=\left[\begin{array}{c}{c}_{t0}+c_{t1}{\mathrm{sin\&beta;}}_{t2}/({\mathrm{cos\&beta;}}_{t2}-1)\\ .\\ .\\ .\\ c_{t0}+c_{t1}{\mathrm{sin\&beta;}}_{{\mathrm{tm}}^{\&prime;}}/({\mathrm{cos\&beta;}}_{{\mathrm{tm}}^{\&prime;}}-1)\\ .\\ .\\ .\\ c_{t0}+c_{t1}{\mathrm{sin\&beta;}}_{tM}/({\mathrm{cos\&beta;}}_{tM}-1)\end{array}\right]$

${\mathrm{\&beta;}}_2=\left[\begin{array}{c}{c}_{r0}+c_{r1}{\mathrm{sin\&beta;}}_{t2}/({\mathrm{cos\&beta;}}_{r2}-1)\\ .\\ .\\ .\\ c_{r0}+c_{r1}{\mathrm{sin\&beta;}}_{{\mathrm{rn}}^{\&prime;}}/({\mathrm{cos\&beta;}}_{{\mathrm{rn}}^{\&prime;}}-1)\\ .\\ .\\ .\\ c_{r0}+c_{r1}{\mathrm{sin\&beta;}}_{\mathrm{rN}}/({\mathrm{cos\&beta;}}_{\mathrm{tN}}-1)\end{array}\right]$

其中，c_t0＝ξ_tsinφ_tcosθ_t，c_t1＝ξ_tsinφ_tsinθ_t，c_r0＝ξ_rsinφ_rcosθ_r，c_r1＝ξ_rsinφ_rsinθ_r；

根据 $U_t\left[\begin{array}{c}{c}_{t0}\\ c_{t1}\end{array}\right]={\mathrm{\&beta;}}_1,$ 可知 $\left[\begin{array}{c}{c}_{t0}\\ c_{t1}\end{array}\right]$ 的求解是一个标准的参数估计问题，可以用最小二乘估计得到 $\left[\begin{array}{c}{c}_{t0}\\ c_{t1}\end{array}\right]$ 的估计值

$\left[\begin{array}{c}{\hat{c}}_{t0}\\ {\hat{c}}_{t1}\end{array}\right]={\left({U_t}^T{U}_t\right)}^{-1}{U_t}^T{\mathrm{\&beta;}}_1,$

其中， $U_t=\left[\begin{array}{cc}1& \sin {\mathrm{\&beta;}}_{t2}/(\cos {\mathrm{\&beta;}}_{t2}-1)\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ 1& \sin {\mathrm{\&beta;}}_{{\mathrm{tm}}^{\&prime;}}/(\cos {\mathrm{\&beta;}}_{{\mathrm{tm}}^{\&prime;}}-1)\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;\\ 1& \sin {\mathrm{\&beta;}}_{\mathrm{tM}}/(\cos {\mathrm{\&beta;}}_{\mathrm{tM}}-1)\end{array}\right].$

根据 $U_r\left[\begin{array}{c}{c}_{r0}\\ c_{r1}\end{array}\right]={\mathrm{\&beta;}}_2,$ 可知 $\left[\begin{array}{c}{c}_{r0}\\ c_{r1}\end{array}\right]$ 的求解是一个标准的参数估计问题，可以用最小二乘估计，得到 $\left[\begin{array}{c}{c}_{r0}\\ c_{r1}\end{array}\right]$ 的估计值

$\left[\begin{array}{c}{\hat{c}}_{r0}\\ {\hat{c}}_{r1}\end{array}\right]={\left({U_r}^T{U}_r\right)}^{-1}{U_r}^T{\mathrm{\&beta;}}_2,$

其中， $U_r=\left[\begin{array}{cc}1& {\mathrm{sin\&beta;}}_{r2}/({\mathrm{cos\&beta;}}_{r2}-1)\\ .& .\\ .& .\\ .& .\\ 1& {\mathrm{sin\&beta;}}_{{\mathrm{rn}}^{\&prime;}}/({\mathrm{cos\&beta;}}_{{\mathrm{rn}}^{\&prime;}}-1)\\ .& .\\ .& .\\ .& .\\ 1& {\mathrm{sin\&beta;}}_{rN}/({\mathrm{cos\&beta;}}_{rN}-1)\end{array}\right].$

进而获得雷达目标的发射方位角估计值发射俯仰角估计值接收方位角估计值和接收俯仰角估计值分别为 ${\widehat{\mathrm{\&phi;}}}_t={\sin }^{-1}(\sqrt{{{\hat{c}}_{t0}}^2+{{\hat{c}}_{t1}}^2}/{\mathrm{\&xi;}}_t),{\widehat{\mathrm{\&theta;}}}_t={\tan }^{-1}({\hat{c}}_{r1}/{\hat{c}}_{r0}),{\widehat{\mathrm{\&phi;}}}_r={\sin }^{-1}(\sqrt{{{\hat{c}}_{r0}}^2+{{\hat{c}}_{r1}}^2}/{\mathrm{\&xi;}}_r).$ 发射方位角估计值发射俯仰角估计值的估计误差如图3所示，接收方位角估计值和接收俯仰角估计值的估计误差如图4所示。

步骤5，根据所述归一化多普勒频率方向向量的估计利用最小二乘算法获得目标的归一化多普勒频率估计值

具体地，根据归一化多普勒频率方向向量令B的每一项都除以第一项ρ，得到再取对数的虚部，得到其相位 $\hat{h}=angle\left({\hat{b}}^{\&prime;}\right)={\&lsqb;0,2{\mathrm{\&pi;f}}_d,\mathrm{...},2\mathrm{\&pi;}(L-1)f_d\&rsqb;}^T,$

其中，ρ表示目标的雷达波发射系数，angle(·)表示取相位。

根据 $P\left[\begin{array}{c}{b}_0\\ f_d\end{array}\right]=\hat{h},$ 可知 $\left[\begin{array}{c}{b}_0\\ f_d\end{array}\right]$ 均求解是一个标准的参数估计问题，可以用最小二乘估计，得到其估计

$\left[\begin{array}{c}{\hat{b}}_0\\ {\hat{f}}_d\end{array}\right]={\left(P^{T}P\right)}^{-1}{P}^T\hat{h},$

其中， $P=\left[\begin{array}{cc}1& 0\\ 1& 2\mathrm{\&pi;}\\ .& .\\ .& .\\ .& .\\ 1& 2\mathrm{\&pi;}(L-1)\end{array}\right],b_0=0.$

归一化多普勒频率估计值的估计误差如图5所示。

步骤6，利用地心直角坐标系建立目标的运动状态空间，基于目标的运动状态空间，采用CS模型建立目标状态方程其中，为目标运动状态变量，包括目标在地心直角坐标系下的位置[x，y，z]、速度加速度k表示目标状态的采样时刻序号，k＝1，2，3，...，X(k-1)和X(k)分别表示第k-1和第k个采样时刻的目标状态，Φ表示状态转移矩阵，U表示控制矩阵，表示第k个采样时刻的目标机动加速度均值，W(k)表示第k个采样时刻的目标状态噪声向量；

具体地，

$\mathrm{\&Phi;}=\left[\begin{array}{ccc}{\mathrm{\&Phi;}}_0& 0& 0\\ 0& {\mathrm{\&Phi;}}_0& 0\\ 0& 0& {\mathrm{\&Phi;}}_0\end{array}\right],{\mathrm{\&Phi;}}_0=\left[\begin{array}{ccc}1& T& \frac{1}{\mathrm{\&alpha;}}(-1+\mathrm{\&alpha;}T+e^{-\mathrm{\&alpha;}T})\\ 0& 1& \frac{1}{\mathrm{\&alpha;}}(1-e^{-\mathrm{\&alpha;}T})\\ 0& 0& e^{-\mathrm{\&alpha;}T}\end{array}\right],$ 其中，α为目标机动常数，T为采样周期；

$U=\left[\begin{array}{ccc}{U}_0& 0& 0\\ 0& U_0& 0\\ 0& 0& U_0\end{array}\right],U_0=\left[\begin{array}{c}\frac{1}{\mathrm{\&alpha;}}(-T+\frac{{\mathrm{\&alpha;T}}^2}{2}+\frac{1-e^{-\mathrm{\&alpha;}T}}{\mathrm{\&alpha;}})\\ T-\frac{1-e^{-\mathrm{\&alpha;}T}}{\mathrm{\&alpha;}}\\ 1-e^{-\mathrm{\&alpha;}T}\end{array}\right],\stackrel{\&OverBar;}{a}\left(k\right)\left[\begin{array}{c}\stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{x}}\left(k\right)\\ \stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{y}}\left(k\right)\\ \stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{z}}\left(k\right)\end{array}\right],\stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{x}}\left(k\right),\stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{y}}\left(k\right),\stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{z}}\left(k\right)$ 分别表示第k个采样时刻x，y，z方向上的机动加速度均值；

$W\left(k\right)=\left[\begin{array}{c}{W}_x\\ W_y\\ W_z\end{array}\right],W_{p_1}={\&Integral;}_{kT}^{(k+1)T}\left[\begin{array}{c}\{-1+\mathrm{\&alpha;}\&lsqb;(k+1)T-\mathrm{\&tau;}\&rsqb;+e^{-\mathrm{\&alpha;}\left(\right(k+1)T-\mathrm{\&tau;})}\}/{\mathrm{\&alpha;}}^2\\ \{1-e^{-\mathrm{\&alpha;}\left(\right(k+1)T-\mathrm{\&tau;})}\}/\mathrm{\&alpha;}\\ e^{-\mathrm{\&alpha;}\left(\right(k+1)T-\mathrm{\&tau;}}\end{array}\right]w_{p_1}\left(\mathrm{\&tau;}\right)d\mathrm{\&tau;},p_1=x,y,z,$ w_x(t)，w_y(t)，w_z(t)表示各方向上的加速度驱动噪声；W(k)为状态噪声向量，满足Q(k)＝E[W(k)W^T(k)]，Q(·)表示状态噪声协方差矩阵，E[·]表示期望函数，kT表示第k个采样点对应的采样时刻，(k+1)T表示第k+1个采样点对应的采样时刻，τ表示从时刻kT到(k+1)T的时间变量。

$Q\left(k\right)=\left[\begin{array}{ccc}{Q}_x\left(k\right)& 0& 0\\ 0& Q_y\left(k\right)& 0\\ 0& 0& Q_z\left(k\right)\end{array}\right],$

$Q_{p_1}\left(k\right)=2{\mathrm{\&alpha;\&sigma;}}_{p_1}^2\left(k\right)Q_0$

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_{p_1}^2\left(k\right)=\frac{4-\mathrm{\&pi;}}{4}{[a_{\max }^{p_1}-{\stackrel{\stackrel{\&OverBar;}{\&CenterDot;\&CenterDot;}}{p}}_1\left(k\right)]}^2,{\stackrel{\stackrel{\&OverBar;}{\&CenterDot;\&CenterDot;}}{p}}_1\left(k\right)\&GreaterEqual;0\\ {\mathrm{\&sigma;}}_{p_1}^2\left(k\right)=\frac{4-\mathrm{\&pi;}}{4}{[a_{-\max }^{p_1}+{\stackrel{\stackrel{\&OverBar;}{\&CenterDot;\&CenterDot;}}{p}}_1\left(k\right)]}^2,{\stackrel{\stackrel{\&OverBar;}{\&CenterDot;\&CenterDot;}}{p}}_1\left(k\right)<0\end{array}\right.,$

其中，表示p₁方向上的正向极限加速度，表示p₁方向上的负向极限加速度，p₁＝x，y，z，表示p₁方向上的加速度均值。

$Q_0=\left[\begin{array}{ccc}{q}_{11}& q_{12}& q_{13}\\ q_{21}& q_{22}& q_{23}\\ q_{31}& q_{32}& q_{33}\end{array}\right]$

$q_{11}=0.5{\mathrm{\&alpha;}}^{-5}(1-e^{-2\mathrm{\&alpha;}T}+2\mathrm{\&alpha;}T+\frac{2{\mathrm{\&alpha;}}^3{T}^3}{3}-2{\mathrm{\&alpha;}}^2{T}^2-4{\mathrm{\&alpha;Te}}^{-\mathrm{\&alpha;}T})$

q₁₂＝0.5α^-4(e^-2αT+1-2e^-αT+2αTe^-αT-2αT+α²T²)

q₁₃＝0.5α^-3(1-e^-2αT-2αTe^-αT)

q₂₂＝0.5α^-3(4e^-αT-3-e^-2αT+2αT)

q₂₃＝0.5α^-2(e^-2αT+1-2e^-αT)

q₃₃＝0.5α^-1(1-e^-2αT)

步骤7，基于目标状态空间，在发射圆阵与接收圆阵的测量坐标系下，根据测元与目标运动状态变量X之间的函数关系建立测量方程Y(k)＝h(X(k))+v(k)，其中，v(k)表示各测元测量噪声构成的向量；

${\widehat{\mathrm{\&theta;}}}_t=\arctan \left(\frac{{\mathrm{\&Delta;z}}_t}{{\mathrm{\&Delta;x}}_t}\right)+\left\{\begin{array}{cc}0,& {\mathrm{\&Delta;x}}_t>0,{\mathrm{\&Delta;z}}_t>0\\ \mathrm{\&pi;},& {\mathrm{\&Delta;x}}_t<0\\ 2\mathrm{\&pi;},& {\mathrm{\&Delta;x}}_t>0,{\mathrm{\&Delta;z}}_t<0\end{array}\right.+{\mathrm{\&epsiv;}}_1$

${\widehat{\mathrm{\&phi;}}}_t=\frac{\mathrm{\&pi;}}{2}-\arctan \left(\frac{{\mathrm{\&Delta;y}}_t}{\sqrt{{\left({\mathrm{\&Delta;x}}_t\right)}^2+{\left({\mathrm{\&Delta;z}}_t\right)}^2}}\right)+{\mathrm{\&epsiv;}}_2$

${\widehat{\mathrm{\&theta;}}}_r=\arctan \left(\frac{{\mathrm{\&Delta;z}}_r}{{\mathrm{\&Delta;x}}_r}\right)+\left\{\begin{array}{cc}0,& {\mathrm{\&Delta;x}}_r>0,{\mathrm{\&Delta;z}}_r>0\\ \mathrm{\&pi;},& {\mathrm{\&Delta;x}}_r<0\\ 2\mathrm{\&pi;},& {\mathrm{\&Delta;x}}_r>0,{\mathrm{\&Delta;z}}_r<0\end{array}\right.+{\mathrm{\&epsiv;}}_3$

${\widehat{\mathrm{\&phi;}}}_r=\frac{\mathrm{\&pi;}}{2}-\arctan \left(\frac{{\mathrm{\&Delta;y}}_r}{\sqrt{{\left({\mathrm{\&Delta;x}}_r\right)}^2+{\left({\mathrm{\&Delta;z}}_r\right)}^2}}\right)+{\mathrm{\&epsiv;}}_4$

${\hat{f}}_d=\frac{v_t+v_r}{c-v_r}\&CenterDot;\frac{f_0}{PRF}+{\mathrm{\&epsiv;}}_5$

其中，v(k)＝[ε₁，ε₂，ε₃，ε₄，ε₅]^T为测量噪声向量，ε₁，ε₂，ε₃，ε₄，ε₅分别为对应的噪声，满足测量噪声协方差矩阵 $R=E\&lsqb;v\left(k\right)v^T\left(k\right)\&rsqb;=diag({\mathrm{\&sigma;}}_{{\mathrm{\&theta;}}_t}^2,{\mathrm{\&sigma;}}_{{\mathrm{\&phi;}}_t}^2,{\mathrm{\&sigma;}}_{{\mathrm{\&theta;}}_r}^2,{\mathrm{\&sigma;}}_{{\mathrm{\&phi;}}_r}^2,{\mathrm{\&sigma;}}_{f_d}^2),$ diag(·)表示取向量构成的对角矩阵，分别表示各种测元的测量噪声方差；Δx_t，Δy_t，Δz_t，v_t分别表示目标在发射圆阵测量直角坐标系下的三维位置坐标和径向速率，Δx_r，Δy_r，Δz_r，v_r分别表示目标在接收圆阵测量直角坐标系下的三维位置坐标和径向速率，c表示空气介质中的光速，f₀表示雷达载波频率，PRF表示雷达脉冲重复频率。

$\left[\begin{array}{c}{\mathrm{\&Delta;x}}_t\\ {\mathrm{\&Delta;y}}_t\\ {\mathrm{\&Delta;z}}_t\end{array}\right]={\mathrm{\&Omega;}}_t\left(\left[\begin{array}{c}x\\ y\\ z\end{array}\right]-\left[\begin{array}{c}{X}_{t0}\\ Y_{t0}\\ Z_{t0}\end{array}\right]\right)$

$\left[\begin{array}{c}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_t\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{y}}_t\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{z}}_t\end{array}\right]={\mathrm{\&Omega;}}_t\left[\begin{array}{c}\stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;}{y}\\ \stackrel{\&CenterDot;}{z}\end{array}\right]$

$v_t=\frac{{\mathrm{\&Delta;x}}_t}{R_t}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_t+\frac{{\mathrm{\&Delta;y}}_t}{R_t}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{y}}_t+\frac{{\mathrm{\&Delta;z}}_t}{R_t}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{z}}_t$

$\left[\begin{array}{c}{\mathrm{\&Delta;x}}_r\\ {\mathrm{\&Delta;y}}_r\\ {\mathrm{\&Delta;z}}_r\end{array}\right]={\mathrm{\&Omega;}}_r\left(\left[\begin{array}{c}x\\ y\\ z\end{array}\right]-\left[\begin{array}{c}{X}_{r0}\\ Y_{r0}\\ Z_{r0}\end{array}\right]\right)$

$\left[\begin{array}{c}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_r\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{y}}_r\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{z}}_r\end{array}\right]={\mathrm{\&Omega;}}_r\left[\begin{array}{c}\stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;}{y}\\ \stackrel{\&CenterDot;}{z}\end{array}\right]$

$v_r=\frac{{\mathrm{\&Delta;x}}_r}{R_r}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_r+\frac{{\mathrm{\&Delta;y}}_r}{R_r}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{y}}_r+\frac{{\mathrm{\&Delta;z}}_r}{R_r}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{z}}_r$

${\mathrm{\&Omega;}}_t=R_y(-\frac{\mathrm{\&pi;}}{2})R_x\left(B_t\right)R_z(-\frac{\mathrm{\&pi;}}{2}+L_t)$

${\mathrm{\&Omega;}}_r=R_y(-\frac{\mathrm{\&pi;}}{2})R_x\left(B_r\right)R_z(-\frac{\mathrm{\&pi;}}{2}+L_r)$

其中，Ω_t表示地心直角系到发射阵测量直角系的转换矩阵，Ω_r表示地心直角系到接收阵测量直角系的转换矩阵，R_x(α)，R_y(α)，R_z(α)分别表示绕x，y，z轴逆时针旋转α角的旋转矩阵，具体为

$R_x\left(\mathrm{\&alpha;}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\&alpha;}& \sin \mathrm{\&alpha;}\\ 0& -\sin \mathrm{\&alpha;}& \cos \mathrm{\&alpha;}\end{array}\right],R_y\left(\mathrm{\&alpha;}\right)=\left[\begin{array}{ccc}\cos \mathrm{\&alpha;}& 0& -\sin \mathrm{\&alpha;}\\ 0& 1& 0\\ \sin \mathrm{\&alpha;}& 0& \cos \mathrm{\&alpha;}\end{array}\right],R_z\left(\mathrm{\&alpha;}\right)=\left[\begin{array}{ccc}\cos \mathrm{\&alpha;}& \sin \mathrm{\&alpha;}& 0\\ -\sin \mathrm{\&alpha;}& \cos \mathrm{\&alpha;}& 0\\ 0& 0& 1\end{array}\right]$

$(x_{t0},y_{t0},z_{t0})=f_{DXZJ}^{DD}(L_t,B_t,h_t)$

$(x_{r0},y_{r0},z_{r0})=f_{DXZJ}^{DD}(L_r,B_r,h_r)$

其中，表示大地坐标到地心直角坐标的转换函数，的具体形式为

$\left\{\begin{array}{c}x=(N+h)\cos B\cos L\\ y=(N+h)\cos B\sin L\\ z=\&lsqb;N(1-e^2)+h\&rsqb;\sin B\end{array}\right.$

其中a为地球参考椭球半长轴，e²为参考椭球第一偏心率的平方。

步骤8，基于测量方程，根据初始两个采样时刻获取的目标测元Y(1)，Y(2)，计算X(1)＝[x(1)，y(1)，z(1)]^T和X(2)＝[x(2)，y(2)，z(2)]^T，进而计算目标运动状态变量的初值和对应的协方差矩阵P(0)，得到滤波初值；具体为：

8a)利用Y(1)中的4个角度计算目标在第1个采样时刻的位置参数x(1)，y(1)，，z(1)；利用Y(2)中的4个角度计算目标在第2个采样时刻的位置参数x(2)，y(2)，z(2)；

根据4个角度θ_t，φ_t，θ_r，φ_r解算目标位置x，y，z的具体过程如下：

i)将4个角度从发射阵与接收阵的测量直角坐标系转至测量地心直角系下(其坐标原点与测量直角坐标系的原点相同，各坐标轴与地心直角坐标系平行)

$\left[\begin{array}{c}{l}_t\\ m_t\\ n_t\end{array}\right]={\mathrm{\&Omega;}}_t^T\left[\begin{array}{c}{\mathrm{sin\&phi;}}_t{\mathrm{cos\&theta;}}_t\\ {\mathrm{cos\&phi;}}_t\\ {\mathrm{sin\&phi;}}_t{\mathrm{sin\&theta;}}_t\end{array}\right]$

${\mathrm{\&theta;}}_t^{\&prime;}=\arctan \left(\frac{n_t}{l_t}\right)+\left\{\begin{array}{cc}0,& l_t>0,n_t>0\\ \mathrm{\&pi;},& l_t<0\\ 2\mathrm{\&pi;},& l_t>0,n_t<0\end{array}\right.$

${\mathrm{\&phi;}}_t^{\&prime;}=\frac{\mathrm{\&pi;}}{2}-\arctan \left(\frac{m_t}{\sqrt{{\left(l_t\right)}^2+{\left(n_t\right)}^2}}\right)$

$\left[\begin{array}{c}{l}_r\\ m_r\\ n_r\end{array}\right]={\mathrm{\&Omega;}}_r^T\left[\begin{array}{c}{\mathrm{sin\&phi;}}_r{\mathrm{cos\&theta;}}_r\\ {\mathrm{cos\&phi;}}_r\\ {\mathrm{sin\&phi;}}_r{\mathrm{sin\&theta;}}_r\end{array}\right]$

${\mathrm{\&theta;}}_r^{\&prime;}=\arctan \left(\frac{n_r}{l_r}\right)+\left\{\begin{array}{cc}0,& l_r>0,n_r>0\\ \mathrm{\&pi;},& l_r<0\\ 2\mathrm{\&pi;},& l_r>0,n_r<0\end{array}\right.$

${\mathrm{\&phi;}}_r^{\&prime;}=\frac{\mathrm{\&pi;}}{2}-\arctan \left(\frac{m_r}{\sqrt{{\left(l_r\right)}^2+{\left(n_r\right)}^2}}\right)$

其中，Ω_t ^T表示地心直角系到发射阵测量直角系的转换矩阵的转置，Ω_r ^T表示地心直角系到接收阵测量直角系的转换矩阵的转置，(l_t，m_t，n_t)为发射阵中心至目标连线的方向向量，θ′_t、φ′_t分别表示测量地心直角系下目标相对于发射机的方位角和俯仰角；(l_r，m_r，n_r)为接收阵中心至目标连线的方向向量，θ′_r、φ′_r分别表示测量地心直角系下目标相对于接收机的方位角和俯仰角；

ii)采用水平投影法解算目标在地心直角系下的位置

$\left\{\begin{array}{c}x=X_{t0}+{\mathrm{\&Delta;X}}_1\\ y=Y_{t0}+\frac{{\mathrm{\&Delta;X}}_1}{{\mathrm{cos\&theta;}}_t^{\&prime;}}{\mathrm{cot\&phi;}}_t^{\&prime;}\\ z=Z_{t0}+{\mathrm{\&Delta;X}}_1{\mathrm{tan\&theta;}}_t^{\&prime;}\end{array}\right.$

其中

${\mathrm{\&Delta;X}}_1=\frac{(X_{t0}-X_{r0}){\mathrm{tan\&theta;}}_r^{\&prime;}-(Z_{t0}-Z_{r0})}{{\mathrm{tan\&theta;}}_t^{\&prime;}-{\mathrm{tan\&theta;}}_r^{\&prime;}}$

其中，(X_t0，Y_t0，Z_t0)为发射阵中心在地心直角坐标系下的位置；ΔX₁为发射阵中心至目标在X轴方向的距离，X_r0为接收阵中心在地心直角坐标系X轴方向的位置，Z_r0为接收阵中心在地心直角坐标系Z轴方向的位置；当k＝1时，获得采样时刻的位置参数X(1)＝[x(1)，y(1)，z(1)]^T，当k＝2时，获得采样时刻的位置参数X(2)＝[x(2)，y(2)，z(2)]^T；k取自然数；

8b)对前两个采样时刻的位置参数进行差分计算目标在第2个采样时刻的近似速度

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}\left(2\right)=\frac{x\left(2\right)-x\left(1\right)}{T}\\ \stackrel{\&CenterDot;}{y}\left(2\right)=\frac{y\left(2\right)-y\left(2\right)}{T}\\ \stackrel{\&CenterDot;}{z}\left(2\right)=\frac{z\left(2\right)-z\left(1\right)}{T}\end{array}\right.$

其中，T表示采样周期；

8c)根据计算结果构造滤波初值，

$\hat{X}\left(0\right)={\&lsqb;x\left(2\right),\stackrel{\&CenterDot;}{x}\left(2\right),0,y\left(2\right),\stackrel{\&CenterDot;}{y}\left(2\right),0,z\left(2\right),\stackrel{\&CenterDot;}{z}\left(2\right),0\&rsqb;}^T$

P(0)＝10⁶I₉

这里I₉表示9维的单位阵。

步骤9，在获得滤波初值后，采用UKF算法，利用第k个采样时刻的测元Y(k)、目标在第k-1个采样时刻的状态估计与协方差矩阵P(k-1)获得目标在第k个采样时刻的状态估计与协方差矩阵P(k)，依次获得目标在各采样时刻的运动状态，完成对空间机动目标的跟踪。

9a)状态向量采样

对状态向量进行sigma点采样，k取自然数，χ_p(k-1)表示第p个状态采样点，p＝0，…，2n_x，n_x＝9，公别表示第p个采样点的均值权重与方差权重，κ为设计参数：

$\begin{array}{cccc}{\mathrm{\&chi;}}_p(k-1)=\hat{X}(k-1)& W_p^{\left(m\right)}=\mathrm{\&kappa;}/(n_x+\mathrm{\&kappa;})& W_p^{\left(c\right)}=\mathrm{\&kappa;}/(n_x/\mathrm{\&kappa;})+(1-{\mathrm{\&alpha;}}^2+\mathrm{\&beta;})& p=0\end{array}$

其中，表示矩阵的第i列，i＝1，2，…，n_x

采样后的向量：

$\mathrm{\&chi;}(k-1)=\&lsqb;\hat{X}(k-1),\hat{X}(k-1)\&PlusMinus;\sqrt{(n_x+\mathrm{\&kappa;})P(k-1)}\&rsqb;$

9b)时间更新

利用CS模型对目标状态采样点进行一步预测，采用加速度均值自适应算法，有

χ_p(k|k-1)＝Φ_aχ_p(k-1)

式中χ_p(k|k-1)表示第p个采样点χ_p(k-1)从k-1时刻到k时刻的预测值。

${\mathrm{\&Phi;}}_a=\left[\begin{array}{ccc}{\mathrm{\&Phi;}}_{a0}& 0& 0\\ 0& {\mathrm{\&Phi;}}_{a0}& 0\\ 0& 0& {\mathrm{\&Phi;}}_{a0}\end{array}\right],{\mathrm{\&Phi;}}_{a0}=\left[\begin{array}{ccc}1& T& \frac{T^2}{2}\\ 0& 1& T\\ 0& 0& 1\end{array}\right]$

$\stackrel{\&OverBar;}{Y}\left(k\right|k-1)=\underset{p=0}{\overset{2n_x}{\&Sigma;}}{W}_p^{\left(m\right)}{\mathrm{\&chi;}}_p(k/k-1)$

$P\left(k\right|k-1)=\underset{p=0}{\overset{2n_x}{\&Sigma;}}{W}_p^{\left(c\right)}\left({\mathrm{\&chi;}}_p\right(k|k-1)-\stackrel{\&OverBar;}{X}(k|k-1\left)\right){\left({\mathrm{\&chi;}}_p\right(k|k-1)-\stackrel{\&OverBar;}{X}(k|k-1\left)\right)}^T+Q\left(k\right|k-1)$

$Q\left(k\right|k-1)=\left[\begin{array}{ccc}{Q}_x\left(k\right|k-1)& 0& 0\\ 0& Q_y\left(k\right|k-1)& 0\\ 0& 0& Q_z\left(k\right|k-1)\end{array}\right],$

$Q_{p_1}\left(k\right|k-1)=2{\mathrm{\&alpha;\&sigma;}}_{p_1}^2\left(k\right|k-1)Q_0$

$\left\{\begin{array}{c}{\mathrm{\&sigma;}}_{p_1}^2\left(k\right|k-1)=\frac{4-\mathrm{\&pi;}}{4}{[a_{\max }^{p_1}-{\stackrel{\stackrel{\&OverBar;}{\&CenterDot;\&CenterDot;}}{p}}_1\left(k\right|k-1)]}^2,{\stackrel{\stackrel{\&OverBar;}{\&CenterDot;\&CenterDot;}}{p}}_1\left(k\right|k-1)\&GreaterEqual;0\\ {\mathrm{\&sigma;}}_{p_1}^2\left(k\right|k-1)=\frac{4-\mathrm{\&pi;}}{4}{[a_{-\max }^{p_1}+{\stackrel{\stackrel{\&OverBar;}{\&CenterDot;\&CenterDot;}}{p}}_1\left(k\right|k-1)]}^2,{\stackrel{\stackrel{\&OverBar;}{\&CenterDot;\&CenterDot;}}{p}}_1\left(k\right|k-1)<0\end{array}\right.,$

表示p₁方向上机动加速度的一步预测，p₁＝x，y，z。

根据状态Sigma采样点的一步预测，利用测元与目标状态变量之间的映射函数h(·)计算观测一步预测的Sigma点，Υ_p(k|k-1)表示观测一步预测的第p个Sigma点，再对其加权求和，计算观测量的一步预测

Υ_p(k|k-1)＝h(χ_p(k|k-1))

9c)状态更新

$G\left(k\right)=P_{xy}\left(k\right|k-1)P_{yy}^{-1}\left(k\right|k-1)$

$\hat{X}\left(k\right)=\stackrel{\&OverBar;}{X}\left(k\right|k-1)+G\left(k\right)\left(Y\right(k)-\stackrel{\&OverBar;}{Y}(k|k-1\left)\right)$

P(k)＝P(k|k-1)-G(k)P_yy(k|k-1)G^T(k)

G(k)表示第k个采样时刻目标状态更新增益矩阵；P_yy(k|k-1)表示观测协方差矩阵的一步预测，P_yy(k|k-1)表示状态观测互协方差矩阵的一步预测。和分别表示第k-1到第k个采样时刻的目标状态一步预测和测元一步预测的均值；G^T(k)表示第k个采样时刻的目标状态更新增益矩阵的转置。

上述式中Q为状态噪声阵，R为测量噪声阵，κ为比例参数，用于调节sigma点和的距离，仅影响二阶之后的高阶矩带来的偏差，一般取κ＝0；α为正值的比例缩放因子，控制sigma点分布的范围，通常取(0，1)的一个较小的值；β为引入f(·)高阶项信息的参数，通常取β＝2。目标位置跟踪误差如图6所示和目标速度跟踪误差如图7所示。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种双基地MIMO雷达空间机动目标跟踪方法，其特征在于，包括以下步骤：

步骤1，分别将双基地MIMO雷达的发射机配置为M个阵元的均匀圆阵、接收机配置为N个阵元的均匀圆阵，并使发射机中M个阵元发射相互正交的波形信号；其中，M表示发射机阵元个数，N表示接收机阵元个数，且M、N均为自然数；

步骤2，利用所述发射机中M个阵元发射相互正交的波形信号，接收机中的N个阵元分别接收所述发射机中M个阵元发射的相互正交的波形信号，并进行匹配滤波，获得匹配滤波后的NM×1维雷达回波信号s和L次快拍积累得到的NM×L维雷达回波信号矩阵S，进而获得矩阵S中第n个接收阵元的M×L维切片矩阵形式S_n，其中，L表示快拍次数，n为自然数，n＝1，2，...,，N；

步骤3，根据N个接收阵元的M×L维切片矩阵S₁～S_N，利用平行因子算法分别获得发射方向向量的估计接收方向向量的估计和归一化多普勒频率方向向量的估计

步骤4，根据所述发射方向向量的估计获得目标相对于发射机的方位角估计值和俯仰角估计值根据所述接收方向向量的估计获得目标相对于接收机的方位角估计值和俯仰角估计值

步骤5，根据所述归一化多普勒频率方向向量的估计利用最小二乘算法获得目标的归一化多普勒频率估计值

步骤6，利用地心直角坐标系建立目标的运动状态空间，基于目标的运动状态空间，采用CS模型建立目标状态方程其中，为目标运动状态变量，包括目标三维位置[x，y，z]、速度加速度k表示目标状态的采样时刻序号，k＝1，2，3，…，X(k-1)和X(k)分别表示第k-1个和第k个采样时刻的目标状态，Φ表示状态转移矩阵，U表示控制矩阵，表示第k个采样时刻的目标机动加速度均值，W(k)表示第k个采样时刻的目标状态噪声向量；

步骤7，基于目标状态空间，在发射均匀圆阵与接收均匀圆阵的测量直角坐标系下，根据测元与目标运动状态变量X之间的函数关系，建立测量方程Y(k)＝h(X(k))+v(k)，其中，v(k)表示各测元的测量噪声向量；Y(k)表示第k个采样时刻的目标测元，h(X(k))表示第k个采样时刻的目标测元与第k个采样时刻的目标运动状态之间的映射关系；

步骤8，基于测量方程，根据初始两个采样时刻获取的目标测元Y(1)，Y(2)，计算X(1)＝[x(1)，y(1)，z(1)]^T和X(2)＝[x(2)，y(2)，z(2)]^T，进而计算目标运动状态变量的初值和对应的协方差矩阵P(0)，得到目标运动状态的滤波初值；

步骤9，获得滤波初值后，采用UKF算法，利用第k个采样时刻的测元Y(k)、目标在第k-1个采样时刻的状态估计与协方差矩阵P(k一1)获得目标在第k个采样时刻的状态估计与协方差矩阵P(k)，依次获得目标在各采样时刻的运动状态，完成对空间机动目标的跟踪。

2.如权利要求1所述一种双基地MIMO雷达空间机动目标跟踪方法，其特征在于，在L次快拍积累中第n个接收阵元的M×L维切片矩阵S_n，其表达式为：

S_n＝A_R(n)A_TB^T+W_n，n＝1，2，...，N

其中，A_R(n)表示接收方向向量A_R的第n个元素，W_n表示第n个接收阵元L次快拍积累中的噪声矩阵；A_T表示发射方向向量，B表示归一化多普勒频率方向向量，其中，ρ为目标雷达波反射系数，f_d为目标的归一化多普勒频率。

3.如权利要求1所述一种双基地MIMO雷达空间机动目标跟踪方法，其特征在于，在步骤3中，所述发射方向向量的估计所述接收方向向量的估计和所述归一化多普勒频率方向向量的估计的具体子步骤为：

a)计算第m个发射阵元的第n个接收阵元的第l次快拍的平行因子三线性模型形式s_m，n，l，其具体表达式为

s_m，n，l＝A_T(m)A_R(n)B(l)+w_m，n，l，m∈{1，…，M}，n∈{1，…，N}，l∈{1，…，L}

其中，A_T(m)表示发射方向向量A_T的第m个元素，A_R(n)表示接收方向向量A_R的第n个元素，B(l)表示归一化多普勒频率方向向量B的第l个元素，w_m，n，l表示s_m，n，l的三维噪声数据集；

b)根据第m个发射阵元的第n接收阵元的第l次快拍的平行因子三线性模型形式s_m，n，l的对称性，对s_m，n，l的第二维和第三维分别进行切片，分别获得第m个发射阵元的L×N维切片矩阵形式F_m、第l次快拍的N×M维切片矩阵形式Z_l，其表达式分别为：

F_m＝BA_T(m)A_R ^T+W_m，m∈{1，2，…，M}

Z_l＝A_RB(l)A_T ^T+W_l，l∈{1，2，...，L}

其中，W_m表示第m个发射阵元噪声，W_l表示第l次快拍噪声；

c)将N个接收阵元的M×L维切片矩阵形式S₁～S_N按列形成NM×L维矩阵S，将M个发射阵元的L×N维切片矩阵形式F₁～F_M按列形成ML×N维矩阵F，将L次快拍的N×M维切片矩阵形式Z₁～Z_L按列形成NL×M维矩阵Z，具体地，

S＝[A_RοA_T]B^T+W_S＝[S₁，…，S_N]^T

F＝[A_TοB]A_R ^T+W_F＝[F₁，…，F_M]^T

Z＝[BοA_R]A_T ^T+W_Z＝[Z₁，…，Z_L]^T

其中，ο表示Khatri-Rao积，W_S表示NM×L维矩阵S的噪声矩阵，W_F表示ML×N维矩阵F的噪声矩阵，W_Z表示NL×M维矩阵Z的噪声矩阵；

d)根据NM×L维矩阵S、ML×N维矩阵F和NL×M维矩阵Z，利用平行因子算法得到发射方向向量的估计接收方向向量的估计和归一化多普勒频率方向向量的估计具体地，

其中，[·]⁺表示取伪逆。

4.如权利要求1所述一种双基地MIMO雷达空间机动目标跟踪方法，其特征在于，在步骤8中，完成滤波初值的计算，具体为：

a)利用Y(1)中的4个角度计算目标在第1个采样时刻的位置参数X(1)＝[x(1)，y(1)，z(1)]^T：利用Y(2)中的4个角度计算目标在第2个采样时刻的位置参数X(2)＝[x(2)，y(2)，z(2)]^T；

根据4个角度θ_t，φ_t，θ_r，φ_r解算目标位置x，y，z的具体过程如下：

θ_t表示目标相对于发射机的方位角，φ_t表示目标相对于发射机的俯仰角，θ_r表示目标相对于接收机的方位角，φ_r表示目标相对于接收机的俯仰角；

i)将4个角度从发射阵与接收阵测量直角坐标系转至测量地心直角系为：

$\left[\begin{array}{c}{l}_t\\ m_t\\ n_t\end{array}\right]={\mathrm{\&Omega;}}_t^T\left[\begin{array}{c}sin{\mathrm{\&phi;}}_{t}cos{\mathrm{\&theta;}}_t\\ {\mathrm{cos\&phi;}}_t\\ \sin {\mathrm{\&phi;}}_t\sin {\mathrm{\&theta;}}_t\end{array}\right]$

${\mathrm{\&theta;}}_t^{\&prime;}=arctan\left(\frac{n_t}{l_t}\right)+\left\{\begin{array}{cc}0,& l_t>0,n_t>0\\ \mathrm{\&pi;},& l_t<0\\ 2\mathrm{\&pi;},& l_t>0,n_t<0\end{array}\right.$

${\mathrm{\&phi;}}_t^{\&prime;}=\frac{\mathrm{\&pi;}}{2}-arctan\left(\frac{m_t}{\sqrt{{\left(l_t\right)}^2+{\left(n_t\right)}^2}}\right)$

$\left[\begin{array}{c}{l}_r\\ m_r\\ n_r\end{array}\right]={\mathrm{\&Omega;}}_r^T\left[\begin{array}{c}sin{\mathrm{\&phi;}}_{r}cos{\mathrm{\&theta;}}_r\\ {\mathrm{cos\&phi;}}_r\\ \sin {\mathrm{\&phi;}}_r\sin {\mathrm{\&theta;}}_r\end{array}\right]$

${\mathrm{\&theta;}}_r^{\&prime;}=arctan\left(\frac{n_r}{l_r}\right)+\left\{\begin{array}{cc}0,& l_r>0,n_r>0\\ \mathrm{\&pi;},& l_r<0\\ 2\mathrm{\&pi;},& l_r>0,n_r<0\end{array}\right.$

${\mathrm{\&phi;}}_r^{\&prime;}=\frac{\mathrm{\&pi;}}{2}-arctan\left(\frac{m_r}{\sqrt{{\left(l_r\right)}^2+{\left(n_r\right)}^2}}\right)$

其中，Ω_t ^T表示地心直角系到发射阵测量直角系的转换矩阵的转置，Ω_r ^T表示地心直角系到接收阵测量直角系的转换矩阵的转置，(l_t，m_t，n_t)为发射阵中心至目标连线的方向向量，θ′_t、φ′_t分别表示测量地心直角系下目标相对于发射机的方位角和俯仰角；(l_r，m_r，n_r)为接收阵中心至目标连线的方向向量，θ′_r、φ′_r分别表示测量地心直角系下目标相对于接收机的方位角和俯仰角；

ii)采用水平投影法解算目标在地心直角系下的位置

$\left\{\begin{array}{c}x\left(k\right)=X_{t0}+\mathrm{\&Delta;}{X}_1\\ y\left(k\right)=Y_{t0}+\frac{{\mathrm{\&Delta;X}}_1}{{\mathrm{cos\&theta;}}_t^{\&prime;}}{\mathrm{cot\&phi;}}_t^{\&prime;}\\ z\left(k\right)=Z_{t0}+{\mathrm{\&Delta;X}}_1{\mathrm{tan\&theta;}}_t^{\&prime;}\end{array}\right.$

其中

${\mathrm{\&Delta;X}}_1=\frac{(X_{t0}-X_{r0}){\mathrm{tan\&theta;}}_r^{\&prime;}-(Z_{t0}-Z_{r0})}{{\mathrm{tan\&theta;}}_t^{\&prime;}-{\mathrm{tan\&theta;}}_r^{\&prime;}}$

其中，(X_t0，Y_t0，Z_t0)为发射阵中心在地心直角坐标系下的位置；ΔX₁为发射阵中心至目标在X轴方向的距离，X_r0为接收阵中心在地心直角坐标系X轴方向的位置，Z_r0为接收阵中心在地心直角坐标系Z轴方向的位置；当k＝1时，获得采样时刻的位置参数X(1)＝[x(1)，y(1)，z(1)]^T，当k＝2时，获得采样时刻的位置参数X(2)＝[x(2)，y(2)，z(2)]^T；k取自然数；

b)对前两个采样时刻的位置参数进行差分计算目标在第2个采样时刻的近似速度

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}\left(2\right)=\frac{x\left(2\right)-x\left(1\right)}{T}\\ \stackrel{\&CenterDot;}{y}\left(2\right)=\frac{y\left(2\right)-y\left(1\right)}{T}\\ \stackrel{\&CenterDot;}{z}\left(2\right)=\frac{z\left(2\right)-z\left(1\right)}{T}\end{array}\right.$

其中，T表示采样周期；

c)根据上述计算结果构造滤波初值，

$\hat{X}\left(0\right)={\&lsqb;x\left(2\right),\stackrel{\&CenterDot;}{x}\left(2\right),0,y\left(2\right),\stackrel{\&CenterDot;}{y}\left(2\right),0,z\left(2\right),\stackrel{\&CenterDot;}{z}\left(2\right),0\&rsqb;}^T$

P(0)＝10⁶I₉

其中，I₉表示9维的单位矩阵。

5.如权利要求1所述一种双基地MIMO雷达空间机动目标跟踪方法，其特征在于，在步骤9中，所述完成对空间机动目标跟踪的具体步骤为：

a)状态向量采样

对状态估计向量进行sigma采样，k取自然数，χ_p(k-1)表示第p个状态采样点，p＝0，…，2n_x，n_x＝9；采样点为：

其中，κ为比例参数，表示矩阵的第i列，i＝1，2，…，n_x

b)时间更新

利用CS模型对目标状态采样点进行一步预测得到χ_p(k|k-1)，一步预测状态转移矩阵为Φ_a，根据目标测元与目标运动状态之间的映射函数h(·)计算观测一步预测的采样点γ_p(k|k-1)，再对其加权求和，计算观测量的一步预测

χ_p(k|k-1)表示第p个状态采样点从k-1时刻到k时刻的预测值；γ_p(k|k-1)表示第p个观测一步预测的采样点；表示第p个采样点的均值权重；为第k个采样时刻的观测量的一步预测；

χ_p(k/k-1)＝Φ_aχ_p(k-1)

γ_p(k|k-1)＝h(χ_p(k|k-1))

c)状态更新

$\hat{X}\left(k\right)=\stackrel{\&OverBar;}{X}\left(k\right|k-1)+G\left(k\right)\left(Y\right(k)-\stackrel{\&OverBar;}{Y}(k|k-1\left)\right)$

P(k)＝P(k|k-1)-G(k)P_yy(k|k-1)G^T(k)

G(k)表示第k个采样时刻目标状态更新增益矩阵；P(k|k-1)表示目标状态协方差矩阵的一步预测；P_yy(k|k-1)表示观测协方差矩阵的一步预测；和分别表示第k-1到第k个采样时刻的目标状态一步预测和测元一步预测的均值；G^T(k)表示第k个采样时刻的目标状态更新增益矩阵的转置。