CN105372652A - 基于接收线阵的mimo雷达空间机动目标跟踪方法 - Google Patents

Abstract

本发明公开了一种基于接收线阵的MIMO雷达空间机动目标跟踪方法，该方法去掉对目标缓慢运动或静止的约束，同时估计三个角度与目标多普勒频率；充分利用三个角度与多普勒频率信息以提高目标跟踪精度；该方法是将双基地MIMO雷达配置为发射阵为均匀圆阵、接收阵为均匀线阵，采用角度多普勒联合估计方法估计出雷达目标相对于发射阵的方位角和俯仰角、相对于接收阵的接收角以及雷达目标的归一化多普勒频率，然后建立方位角、俯仰角、接收角以及雷达目标的归一化多普勒频率测元与目标运动参数之间的测量方程，采用“当前”统计模型建立目标的状态方程，利用无迹卡尔曼滤波算法估计目标的运动参数并进行预测，实现目标的跟踪。

Description

基于接收线阵的MIMO雷达空间机动目标跟踪方法

技术领域

本发明属于雷达数据处理领域，具体涉及基于接收线阵的MIMO雷达空间机动目标跟踪方法。

背景技术

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

西安电子科技大学申请的专利技术“双基地米波雷达目标三维精确定位方法”(专利申请号：201218001807.9)中公开了一种双基地多输入多输出雷达的多目标三维定位方法，利用均匀圆阵同时估计信号的方位角和俯仰角的特性，将双基地雷达配置为发射阵为均匀圆阵，接收阵为均匀线阵，利用模式激励法估计目标相对于均匀圆阵的方位角与俯仰角，利用ESPRIT方法估计目标相对于均匀线阵的接收角，利用这三个角度以及基线的长度估计目标的三维位置坐标。

上述方法在估计三个角度时假定目标多普勒频率为零，但实际上要求目标缓慢运动或静止，因此上述方法不符合空间机动目标的实际；由于上述方法在确定目标位置坐标时，均采用的是多角度交会的几何解算方法，所以定位误差随目标距离增大而迅速放大。基于上述方法存在的缺点，本发明去掉对目标缓慢运动或静止的约束，采用平行因子三线性分解法估计三个角度与目标多普勒频率；采用UKF滤波算法而不是几何解算方法用于估计目标运动状态，充分利用三个角度与多普勒频率，增强对测量噪声的抑制能力，提高目标跟踪精度。

发明内容

有鉴于此，本发明提供基于接收线阵的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)中的3个角度计算目标在第1个采样时刻的位置参数X(1)＝[x(1)，y(1)，z(1)]^T；利用Y(2)中的3个角度计算目标在第2个采样时刻的位置参数X(2)＝[x(2)，y(2)，z(2)]^T；

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

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

将发射阵中心至目标连线的方向矢量从发射阵垂线测量坐标系转至接收阵垂线测量坐标系为：

$\left[\begin{array}{c}{l}_t^{\′}\\ m_t^{\′}\\ n_t^{\′}\end{array}\right]={\mathrm{\Ω}}_r{\mathrm{\Ω}}_t^T\left[\begin{array}{c}{l}_t\\ m_t\\ n_t\end{array}\right]={\mathrm{\Ω}}_r{\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]$

记的长度为ω，将 $\left[\begin{array}{c}{\mathrm{\Δx}}_r\\ {\mathrm{\Δy}}_r\\ {\mathrm{\Δz}}_r\end{array}\right]=\left[\begin{array}{c}{\mathrm{\ωl}}_t^{\′}+x_{O_t}\\ {\mathrm{\ωm}}_t^{\′}+y_{O_t}\\ {\mathrm{\ωn}}_t^{\′}+z_{O_t}\end{array}\right]$ 带入方程 $cos{\widehat{\mathrm{\α}}}_r=\frac{{\mathrm{\Δx}}_r{\mathrm{cosA}}_r+{\mathrm{\Δz}}_r\sin A_r}{\sqrt{{\mathrm{\Δx}}_r^2+{\mathrm{\Δy}}_r^2+{\mathrm{\Δz}}_r^2}}$ 得到关于ω的方程，求解该方程，记其正根为ω＝ω₁＞0，则目标的位置初值为

$\left[\begin{array}{c}x\left(k\right)\\ y\left(k\right)\\ z\left(k\right)\end{array}\right]={\mathrm{\Ω}}_r^T\left[\begin{array}{c}{\mathrm{\Δx}}_r\\ {\mathrm{\Δy}}_r\\ {\mathrm{\Δz}}_r\end{array}\right]+\left[\begin{array}{c}{X}_{r0}\\ Y_{r0}\\ X_{r0}\end{array}\right]={\mathrm{\Ω}}_r^T\left[\begin{array}{c}{\mathrm{\ω}}_1{l}_t^{\′}+x_{O_t}\\ {\mathrm{\ω}}_1{m}_t^{\′}+y_{O_t}\\ {\mathrm{\ω}}_1{n}_t^{\′}+z_{O_t}\end{array}\right]+\left[\begin{array}{c}{X}_{r0}\\ Y_{r0}\\ X_{r0}\end{array}\right]$

其中，Ω_r表示地心直角系到接收阵垂线测量坐标系的转换矩阵，Ω_t表示地心直角系到发射阵垂线测量坐标系的转换矩阵，上标T表示矩阵的转置，(l_t，m_t，n_t)为发射阵中心至目标连线在发射阵垂线测量坐标系下的方向向量；A_r为接收线阵的安装方位角；(l′_t，m′_t，n′_t)为发射阵中心至目标连线在接收阵垂线测量坐标系下的方向向量；(Δx_r，Δy_r，Δz_r)表示目标在接收阵垂线测量坐标系下的位置坐标；表示发射阵中心O_t在接收阵垂线测量坐标系下的三维坐标；(X_r0，Y_r0，Z_r0)表示接收阵中心O_r在地心直角系下的三维坐标；当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{\·}{x}\left(2\right)=\frac{x\left(2\right)-x\left(1\right)}{T}\\ \stackrel{\·}{y}\left(2\right)=\frac{y\left(2\right)-y\left(1\right)}{T}\\ \stackrel{\·}{z}\left(2\right)=\frac{z\left(2\right)-z\left(1\right)}{T}\end{array}\right.$

其中，T表示采样周期；

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

$\hat{X}\left(0\right)={\[x\left(2\right),\stackrel{\·}{x}\left(2\right),0,y\left(2\right),\stackrel{\·}{y}\left(2\right),0,z\left(2\right),\stackrel{\·}{z}\left(2\right),0\]}^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))

$\stackrel{\‾}{Y}\left(k\right|k-1)=\underset{p=0}{\overset{2n_x}{\Σ}}{W}_p^{\left(m\right)}{\mathrm{\γ}}_p(k/k-1)$

c)状态更新

$\hat{X}\left(k\right)=\stackrel{\‾}{X}\left(k\right|k-1)+G\left(k\right)\left(Y\right(k)-\stackrel{\‾}{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)在进行目标运动状态参数解算时，传统的几何解算方法只能利用角度测元，无法利用多普勒频率测元，不但造成了测元信息浪费，而且定位误差随目标距离增加而放大，本发明采用UKF算法估计目标的运动状态参数，一方面能够融合利用所有测元信息，另一方面增强了对测量噪声的抑制能力，这两方面共同提高了目标运动状态解算精度。

(3)本发明采用了适用于机动目标跟踪的CS模型，可对空间机动目标的位置、速度进行一定时间的预测，实现对空间机动目标的有效稳定跟踪。

附图说明

图1为基于接收线阵的MIMO雷达空间机动目标跟踪方法实现流程图。

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

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

图4为用本发明方法仿真信噪比20dB情况下某临近空间运动目标相对于接收线阵的接收角的估计误差图；

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

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

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

具体实施方式

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

本发明提供了了基于接收线阵的MIMO雷达空间机动目标跟踪方法，如图1所示，实现本发明的具体步骤如下：

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

具体地，如图2所示，坐标系采用地心直角坐标系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)，且接收线阵的安装方位角为A_r(定义为与当地正北方向的夹角)。发射圆阵半径为r_t，其阵元个数M＝2floor(2πr_t/λ)+1，floor(·)表示向下取整运算，λ表示发射波的波长；接收线阵的阵元间距离为d_r，且取d_r＝λ/2。

假设空间有一目标P，在地心系下的位置坐标为(x，y，z)、速度坐标为发射圆阵的圆心至目标连线O_tP与发射圆阵平面法线的夹角为俯仰角，记为φ_t，O_tP在发射圆阵平面的投影与当地正北方向的夹角为方位角，记为为θ_t；接收线阵的中点至目标连线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{\ρe}}^{j2\mathrm{\π}(l-1)f_d}\[a_r\left({\mathrm{\α}}_r\right)\⊗a_t({\mathrm{\θ}}_t,{\mathrm{\φ}}_t)\]+\tilde{n}\left(l\right)$

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

$a_r\left({\mathrm{\α}}_r\right)=\left[\begin{array}{c}1\\ .\\ .\\ .\\ e^{j\frac{2{\mathrm{\πd}}_r}{\mathrm{\λ}}(n-1){\mathrm{cos\α}}_r}\\ .\\ .\\ .\\ e^{j\frac{2{\mathrm{\πd}}_r}{\mathrm{\λ}}(N-1){\mathrm{cos\α}}_r}\end{array}\right]=\left[\begin{array}{c}{a}_{r1}\left({\mathrm{\α}}_r\right)\\ .\\ .\\ .\\ a_{rn}\left({\mathrm{\α}}_r\right)\\ .\\ .\\ .\\ a_{rN}\left({\mathrm{\α}}_r\right)\end{array}\right]\stackrel{\mathrm{\Δ}}{=}{A}_R$

这里，n＝1，2，…，N，符号表示“简记为”；a_rn(α_r)表示第n个接收阵元的接收方向因子。

$a_t({\mathrm{\θ}}_t,{\mathrm{\φ}}_t)=\left[\begin{array}{c}{e}^{j\frac{2{\mathrm{\πr}}_t}{\mathrm{\λ}}{\mathrm{sin\φ}}_t\cos ({\mathrm{\θ}}_t-{\mathrm{\β}}_{t1})}\\ .\\ .\\ .\\ e^{j\frac{2{\mathrm{\πr}}_t}{\mathrm{\λ}}{\mathrm{sin\φ}}_t\cos ({\mathrm{\θ}}_t-{\mathrm{\β}}_{tm})}\\ .\\ .\\ .\\ e^{j\frac{2{\mathrm{\πr}}_t}{\mathrm{\λ}}{\mathrm{sin\φ}}_t\cos ({\mathrm{\θ}}_t-{\mathrm{\β}}_{tM})}\end{array}\right]=\left[\begin{array}{c}{a}_{t1}({\mathrm{\θ}}_t,{\mathrm{\φ}}_t)\\ .\\ .\\ .\\ a_{tm}({\mathrm{\θ}}_t,{\mathrm{\φ}}_t)\\ .\\ .\\ .\\ a_{tM}({\mathrm{\θ}}_t,{\mathrm{\φ}}_t)\end{array}\right]\stackrel{\mathrm{\Δ}}{=}{A}_T$

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

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

$S=\[s\left(1\right),s\left(2\right),\mathrm{...},s\left(L\right)\]=\[a_r\left({\mathrm{\α}}_r\right)\⊗a_t({\mathrm{\θ}}_t,{\mathrm{\φ}}_t)\]B^T+\tilde{N}$

其中， $B={\[\mathrm{\ρ},{\mathrm{\ρe}}^{j2{\mathrm{\πf}}_d},\mathrm{...},{\mathrm{\ρe}}^{j2\mathrm{\π}(L-1)f_d}\]}^T,\tilde{N}=\[\tilde{n}\left(1\right),\tilde{n}\left(2\right),\mathrm{...},\tilde{n}\left(L\right)\],$ 上标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}

其中，W_m表示第m个发射阵元噪声，W_l表示第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_Z＝[Z₁，…，Z_L]^T

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

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

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

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

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

${\hat{a}}_t({\mathrm{\φ}}_t,{\mathrm{\θ}}_t)={\[e^{j\frac{2{\mathrm{\πr}}_t}{\mathrm{\λ}}{\mathrm{sin\φ}}_t\·\cos ({\mathrm{\θ}}_t-{\mathrm{\β}}_{t1})},\mathrm{...},e^{j\frac{2{\mathrm{\πr}}_t}{\mathrm{\λ}}{\mathrm{sin\φ}}_t\·\cos ({\mathrm{\θ}}_t-{\mathrm{\β}}_{tm})},\mathrm{...},e^{j\frac{2{\mathrm{\πr}}_t}{\mathrm{\λ}}{\mathrm{sin\φ}}_t\·\cos ({\mathrm{\θ}}_t-{\mathrm{\β}}_{tM})}\]}^T$

一般地，β_t1＝0，中的每一项都除以第一项然后去掉其第一项，得到新向量a₁，再取a₁对数的虚部，即a′₁＝Im(lna₁)，a′₁的表达式为：

$a_1^{\′}=\left[\begin{array}{c}{\mathrm{\ξ}}_t{\mathrm{sin\φ}}_t{\mathrm{cos\θ}}_t({\mathrm{cos\β}}_{t2}-1)+{\mathrm{\ξ}}_t{\mathrm{sin\φ}}_t{\mathrm{sin\θ}}_t{\mathrm{sin\β}}_{t2}\\ .\\ .\\ .\\ {\mathrm{\ξ}}_t{\mathrm{sin\φ}}_t{\mathrm{cos\θ}}_t({\mathrm{cos\β}}_{{\mathrm{tm}}^{\′}}-1)+{\mathrm{\ξ}}_t{\mathrm{sin\φ}}_t{\mathrm{sin\θ}}_t{\mathrm{sin\β}}_{{\mathrm{tm}}^{\′}}\\ .\\ .\\ .\\ {\mathrm{\ξ}}_t{\mathrm{sin\φ}}_t{\mathrm{cos\θ}}_t({\mathrm{cos\β}}_{tM}-1)+{\mathrm{\ξ}}_t{\mathrm{sin\φ}}_t{\mathrm{sin\θ}}_t{\mathrm{sin\β}}_{tM}\end{array}\right]$

其中，ξ_t＝2πr_t/λ，ξ_r＝2πr_r/λ，m′∈{2，3，...，M}；

a′₁中的第m′项除以(cosβ_tm′-1)，得到β₁，β₁的表达式为：

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

其中，c_t0＝ξ_tsinφ_tcosθ_t，c_t1＝ξ_tsinφ_tsinθ_t；

根据 $U_t\left[\begin{array}{c}{c}_{t0}\\ c_{t1}\end{array}\right]={\mathrm{\β}}_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{\β}}_1,$

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

进而得到雷达目标的发射方位角估计值发射俯仰角估计值

${\widehat{\mathrm{\θ}}}_t=\arctan ({\hat{c}}_{t1}/{\hat{c}}_{t0}),{\widehat{\mathrm{\φ}}}_t=\arcsin (\sqrt{{{\hat{c}}_{t0}}^2+{{\hat{c}}_{t1}}^2}/{\mathrm{\ξ}}_t),$

接收方向向量的估计为其表达式为：

${\hat{a}}_r\left({\mathrm{\α}}_r\right)={\[1,e^{j\frac{2{\mathrm{\πd}}_r}{\mathrm{\λ}}{\mathrm{cos\α}}_r},\mathrm{...},e^{j\frac{2{\mathrm{\πd}}_r}{\mathrm{\λ}}(n-1){\mathrm{cos\α}}_r},\mathrm{...},e^{j\frac{2{\mathrm{\πd}}_r}{\mathrm{\λ}}(N-1){\mathrm{cos\α}}_r}\]}^T$

取对数的虚部，得到a′₂＝Im(lna₂)，a′₂的表达式为：

$a_2^{\′}={\[0,\frac{2{\mathrm{\πd}}_r}{\mathrm{\λ}}\cos \mathrm{\α},\mathrm{...},\frac{2{\mathrm{\πd}}_r}{\mathrm{\λ}}(n-1){\mathrm{cos\α}}_r,\mathrm{...},\frac{2{\mathrm{\πd}}_r}{\mathrm{\λ}}(N-1){\mathrm{cos\α}}_r\]}^T,$

根据 $\mathrm{\Γ}\left[\begin{array}{c}{c}_{r0}\\ c_{r1}\end{array}\right]=a_2^{\′},$ 可知 $\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({\mathrm{\Γ}}^T\mathrm{\Γ}\right)}^{-1}{\mathrm{\Γ}}^T{a}_2^{\′}$

其中， $\mathrm{\Γ}=\left[\begin{array}{cc}1& 0\\ .& .\\ .& .\\ .& .\\ 1& 2\mathrm{\π}\frac{d_r}{\mathrm{\λ}}(n-1)\\ .& .\\ .& .\\ .& .\\ 1& 2\mathrm{\π}\frac{d_r}{\mathrm{\λ}}(N-1)\end{array}\right],$ c_r0＝0，c_r1＝cosα_r

进而得到雷达目标相对接收线阵的接收角估计值为：

${\widehat{\mathrm{\α}}}_r={\cos }^{-1}\left({\hat{c}}_{r1}\right)$

发射方位角估计值发射俯仰角估计值的估计误差如图3所示；接收线阵的接收角估计值如图4为所示。

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

具体地，根据归一化多普勒频率方向向量令B的每一项都除以第一项ρ，得到再取对数的虚部，得到其相位 $\hat{h}=angle\left({\hat{b}}^{\′}\right)={\[0,2{\mathrm{\πf}}_d,\mathrm{...},2\mathrm{\π}(L-1)f_d\]}^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{\π}\\ .& .\\ .& .\\ .& .\\ 1& 2\mathrm{\π}(L-1)\end{array}\right],$ b₀＝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{\Φ}=\left[\begin{array}{ccc}{\mathrm{\Φ}}_0& 0& 0\\ 0& {\mathrm{\Φ}}_0& 0\\ 0& 0& {\mathrm{\Φ}}_0\end{array}\right],{\mathrm{\Φ}}_0=\left[\begin{array}{ccc}1& T& \frac{1}{\mathrm{\α}}(-1+\mathrm{\α}T+e^{-\mathrm{\α}T})\\ 0& 1& \frac{1}{\mathrm{\α}}(1-e^{-\mathrm{\α}T})\\ 0& 0& e^{-\mathrm{\α}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{\α}}(-T+\frac{{\mathrm{\αT}}^2}{2}+\frac{1-e^{-\mathrm{\α}T}}{\mathrm{\α}})\\ T-\frac{1-e^{-\mathrm{\α}T}}{\mathrm{\α}}\\ 1-e^{-\mathrm{\α}T}\end{array}\right],\stackrel{\‾}{a}\left(k\right)=\left[\begin{array}{c}\stackrel{\‾}{\stackrel{\·\·}{x}}\left(k\right)\\ \stackrel{\‾}{\stackrel{\·\·}{y}}\left(k\right)\\ \stackrel{\‾}{\stackrel{\·\·}{z}}\left(k\right)\end{array}\right],\stackrel{\‾}{\stackrel{\·\·}{x}}\left(k\right),\stackrel{\‾}{\stackrel{\·\·}{y}}\left(k\right),\stackrel{\‾}{\stackrel{\·\·}{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}={\∫}_{kT}^{(k+1)T}\left[\begin{array}{c}\{-1+\mathrm{\α}\[(k+1)T-\mathrm{\τ}\]+e^{-\mathrm{\α}\left(\right(k+1)T-\mathrm{\τ})}\}/{\mathrm{\α}}^2\\ \{1-e^{-\mathrm{\α}\left(\right(k+1)T-\mathrm{\τ})}\}/\mathrm{\α}\\ e^{-\mathrm{\α}\left(\right(k+1)T-\mathrm{\τ})}\end{array}\right]w_{p_1}\left(\mathrm{\τ}\right)d\mathrm{\τ},p_1=x,y,z,$ w_x(t)，w_y(t)，w_z(t)表示x，y，z方向上的加速度驱动噪声；W(k)为状态噪声向量，满足 $E\[W\left(k\right)W^T(k+\mathrm{\ξ})\]=0,(\∀\mathrm{\ξ}\≠0),Q\left(k\right)=E\[W\left(k\right)W^T\left(k\right)\],$ 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{\α\σ}}_{p_1}^2\left(k\right)Q_0$

$\left\{\begin{array}{cc}{\mathrm{\&sigma;}}_{p_1}^2\left(k\right)=\frac{4-\mathrm{\&pi;}}{4}{\&lsqb;a_{\max }^{p_1}-{\stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{p}}}_1\left(k\right)\&rsqb;}^2,& {\stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{p}}}_1\left(k\right)\&GreaterEqual;0\\ {\mathrm{\&sigma;}}_{p_1}^2\left(k\right)=\frac{4-\mathrm{\&pi;}}{4}{\&lsqb;a_{-\max }^{p_1}+{\stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{p}}}_1\left(k\right)\&rsqb;}^2,& {\stackrel{\&OverBar;}{\stackrel{\&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)表示各测元测量噪声构成的向量；

具体地，垂线测量坐标系的定义为：原点o为发射圆阵的圆心或接收线阵的中点，ox轴在原点o所在水平面内且指向天文北，oy轴沿过原点o的铅垂线方向指向上，oz轴、ox轴和oy轴构成右手直角坐标系。

目标相对发射圆阵的方位角估计值与俯仰角估计值为

${\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$

接收线阵的安装方位角为A_r，线阵在其测量坐标系内的坐标为[cosA_r，0，sinA_r]，目标相对接收线阵的接收角估计值为

${\widehat{\mathrm{\&alpha;}}}_r=\mathrm{arc}cos\frac{{\mathrm{\&Delta;x}}_r\cos A_r+{\mathrm{\&Delta;z}}_r\sin A_r}{\sqrt{{\mathrm{\&Delta;x}}_r^2+{\mathrm{\&Delta;y}}_r^2+{\mathrm{\&Delta;z}}_r^2}}+{\mathrm{\&epsiv;}}_3$

目标运动引起的多普勒频率为

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

其中，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{\&alpha;}}_r}^2,{\mathrm{\&sigma;}}_{f_d}^2),$ 分别表示各种测元的测量噪声方差；Δ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]=\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}}_t\\ \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)中的3个角度计算目标在第1个采样时刻的位置参数X(1)＝[x(1)，y(1)，z(1)]^T，利用Y(2)中的3个角度计算目标在第2个采样时刻的位置参数X(2)＝[x(2)，y(2)，z(2)]^T；

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

发射阵圆心至目标连线在发射阵垂线测量坐标系下的方向向量为

$\left[\begin{array}{c}{l}_t\\ m_t\\ n_t\end{array}\right]=\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]$

将其转至接收阵垂线测量坐标系下，的方向向量为

$\left[\begin{array}{c}{l}_t^{\&prime;}\\ m_t^{\&prime;}\\ n_t^{\&prime;}\end{array}\right]={\mathrm{\&Omega;}}_r{\mathrm{\&Omega;}}_t^T\left[\begin{array}{c}{l}_t\\ m_t\\ n_t\end{array}\right]$

发射阵圆心O_t在接收阵垂线测量坐标系下的坐标为

$\left[\begin{array}{c}{x}_{O_t}\\ y_{O_t}\\ z_{O_t}\end{array}\right]={\mathrm{\&Omega;}}_r\left[\begin{array}{c}{X}_{t0}-X_{r0}\\ Y_{t0}-Y_{r0}\\ Z_{t0}-Z_{r0}\end{array}\right]$

目标P在接收阵垂线测量坐标系下的坐标为(Δx_r，Δy_r，Δz_r)，于是的方向向量为

这里表示向量的模，即目标相对于发射阵圆心的距离，令则有

$\left[\begin{array}{c}{\mathrm{\&Delta;x}}_r\\ {\mathrm{\&Delta;y}}_r\\ {\mathrm{\&Delta;z}}_r\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&omega;l}}_t^{\&prime;}+x_{O_t}\\ {\mathrm{\&omega;m}}_t^{\&prime;}+y_{O_t}\\ {\mathrm{\&omega;n}}_t^{\&prime;}+z_{O_t}\end{array}\right]$

将其带入方程整理后获得关于ω的一元二次方程，该方程的正根即为ω的取值，令求解出的ω＝ω₁＞0，则目标的位置初值为

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]={\mathrm{\&Omega;}}_r^T\left[\begin{array}{c}{\mathrm{\&Delta;x}}_r\\ {\mathrm{\&Delta;y}}_r\\ {\mathrm{\&Delta;z}}_r\end{array}\right]+\left[\begin{array}{c}{X}_{r0}\\ Y_{r0}\\ X_{r0}\end{array}\right]={\mathrm{\&Omega;}}_r^T\left[\begin{array}{c}{\mathrm{\&omega;}}_1{l}_t^{\&prime;}+x_{O_t}\\ {\mathrm{\&omega;}}_1{m}_t^{\&prime;}+y_{O_t}\\ {\mathrm{\&omega;}}_1{n}_t^{\&prime;}+z_{O_t}\end{array}\right]+\left[\begin{array}{c}{X}_{r0}\\ Y_{r0}\\ X_{r0}\end{array}\right]$

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(1\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-1)，k取自然数，χ_p(k-1)表示第p个状态采样点，p＝0，…，2n_x，这里n_x＝9，分别表示第p个采样点的均值权重与方差权重，κ为设计参数：

采样策略可根据应用要求选择，通常选择2n_x+1个sigma点对称采样：

$\begin{array}{cccc}\mathrm{\&chi;}(k-1)=\hat{X}(k-1)& W_0^{\left(m\right)}=\mathrm{\&kappa;}/(n_x+\mathrm{\&kappa;})& W_0^{\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)=\left[\begin{array}{cc}\hat{X}(k-1)& \hat{X}(k-1)\&PlusMinus;\sqrt{(n_x+\mathrm{\&kappa;})P(k-1)}\end{array}\right]$

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;}{X}\left(k\right|k-1)\underset{p=0}{\overset{2n_x}{\&Sigma;}}{W}_p^{\left(m\right)}{\mathrm{\&chi;}}_p\left(k\right|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}{cc}{\mathrm{\&sigma;}}_{p_1}^2\left(k\right|k-1)=\frac{4-\mathrm{\&pi;}}{4}{\&lsqb;a_{\max }^{p_1}-{\stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{p}}}_1\left(k\right|k-1)\&rsqb;}^2,& {\stackrel{\&OverBar;}{\stackrel{\&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}{\&lsqb;a_{-\max }^{p_1}+{\stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{p}}}_1\left(k\right|k-1)\&rsqb;}^2,& {\stackrel{\&OverBar;}{\stackrel{\&CenterDot;\&CenterDot;}{p}}}_1\left(k\right|k-1)<0\end{array}\right.,$

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

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

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

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

9c)状态更新

$P_{yy}\left(k\right|k-1)=\underset{p=0}{\overset{2n_x}{\&Sigma;}}{W}_p^{\left(c\right)}\left({\mathrm{\&gamma;}}_p\right(k|k-1)-\stackrel{\&OverBar;}{Y}(k|k-1\left)\right){\left({\mathrm{\&gamma;}}_p\right(k|k-1)-\stackrel{\&OverBar;}{Y}(k|k-1\left)\right)}^T+R\left(k\right)$

$P_{xy}\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{\&gamma;}}_p\right(k|k-1)-\stackrel{\&OverBar;}{Y}(k|k-1\left)\right)}^T$

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

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_xy(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，nl＝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)中的3个角度计算目标在第1个采样时刻的位置参数X(1)＝[x(1)，y(1)，z(1)]^T；利用Y(2)中的3个角度计算目标在第2个采样时刻的位置参数X(2)＝[x(2)，y(2)，z(2)]^T；

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

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

将发射阵中心至目标连线的方向矢量从发射阵垂线测量坐标系转至接收阵垂线测量坐标系为：

$\left[\begin{array}{c}{l}_t^{\&prime;}\\ m_t^{\&prime;}\\ n_t^{\&prime;}\end{array}\right]={\mathrm{\&Omega;}}_r{\mathrm{\&Omega;}}_t^T\left[\begin{array}{c}{l}_t\\ m_t\\ n_t\end{array}\right]={\mathrm{\&Omega;}}_r{\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]$

记的长度为ω，将 $\left[\begin{array}{c}{\mathrm{\&Delta;x}}_r\\ {\mathrm{\&Delta;y}}_r\\ {\mathrm{\&Delta;z}}_r\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&omega;l}}_t^{\&prime;}+x_{O_t}\\ {\mathrm{\&omega;m}}_t^{\&prime;}+y_{O_t}\\ {\mathrm{\&omega;n}}_t^{\&prime;}+z_{O_t}\end{array}\right]$ 带入方程 $cos{\widehat{\mathrm{\&alpha;}}}_r=\frac{{\mathrm{\&Delta;x}}_r{\mathrm{cosA}}_r+{\mathrm{\&Delta;z}}_r{\mathrm{sinA}}_r}{\sqrt{{\mathrm{\&Delta;x}}_r^2+{\mathrm{\&Delta;y}}_r^2+{\mathrm{\&Delta;z}}_r^2}}$ 得到关于ω的方程，求解该方程，记其正根为ω＝ω₁＞0，则目标的位置初值为

$\left[\begin{array}{c}x\left(k\right)\\ y\left(k\right)\\ z\left(k\right)\end{array}\right]={\mathrm{\&Omega;}}_r^T\left[\begin{array}{c}{\mathrm{\&Delta;x}}_r\\ {\mathrm{\&Delta;y}}_r\\ {\mathrm{\&Delta;z}}_r\end{array}\right]+\left[\begin{array}{c}{X}_{r0}\\ Y_{r0}\\ Z_{r0}\end{array}\right]={\mathrm{\&Omega;}}_r^T\left[\begin{array}{c}{\mathrm{\&omega;}}_1{l}_t^{\&prime;}+x_{O_t}\\ {\mathrm{\&omega;}}_1{m}_t^{\&prime;}+y_{O_t}\\ {\mathrm{\&omega;}}_1{n}_t^{\&prime;}+z_{O_t}\end{array}\right]+\left[\begin{array}{c}{X}_{r0}\\ Y_{r0}\\ Z_{r0}\end{array}\right]$

其中，Ω_r表示地心直角系到接收阵垂线测量坐标系的转换矩阵，Ω_t表示地心直角系到发射阵垂线测量坐标系的转换矩阵，上标T表示矩阵的转置，(l_t，m_t，n_t)为发射阵中心至目标连线在发射阵垂线测量坐标系下的方向向量；A_r为接收线阵的安装方位角；(l′_t，m′_t，n′_t)为发射阵中心至目标连线在接收阵垂线测量坐标系下的方向向量；(Δx_r，Δy_r，Δz_r)表示目标在接收阵垂线测量坐标系下的位置坐标；表示发射阵中心O_t在接收阵垂线测量坐标系下的三维坐标；(X_r0，Y_r0，Z_r0)表示接收阵中心O_r在地心直角系下的三维坐标；当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)-x\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_z，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-¹))

$\stackrel{\&OverBar;}{Y}\left(k\right|k-1)=\underset{p=0}{\overset{2n_x}{\mathrm{\&Sigma;}}}{W}_p^{\left(m\right)}{\mathrm{\&gamma;}}_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个采样时刻目标状态更新增益矩阵的转置。