一种基于强散射点目标探测的弱散射点目标跟踪逼近方法

技术领域

本发明涉及雷达探测与目标跟踪技术领域，具体地说，是一种利用单脉冲雷达测 角技术与无味卡尔曼滤波预测技术实现的，基于强散射点目标探测的，弱散射点目标跟踪 逼近方法。

背景技术

目前，单脉冲雷达因其较好的测角精度而被广泛地应用于目标跟踪过程中。它通 过对不同方向测角天线接收到的电磁波的振幅或相位的处理，从而得到测量目标的角度信 息，再经过数据关联和滤波更新等数据处理过程完成对目标的跟踪，这一过程也是传统的 先检测后跟踪(DetectBeforeTrack，DBT)技术的处理流程。这种方式实现简单，计算量 小，但在处理散射特性较弱的目标(弱目标)跟踪问题时存在检测困难。

近年来，又有学者提出检测前跟踪(TrackBeforeDetect，TBD)技术来解决弱目 标的跟踪问题，它的核心思想是考虑空间目标运动的连续性，在检测决策前，通过在目标轨 迹上的回波能量积累，以提高信噪比，再通过最终的门限判决，实现弱目标检测并输出其轨 迹。相较于DBT技术，TBD技术虽然适用于弱目标的检测与跟踪，但其要求时间的积累和较高 的存储量及计算量，实现难度大，实时性较差。

发明内容

本发明在DTB技术的思想上，提出基于强目标探测的弱目标跟踪逼近模型，即是通 过对相关强目标的量测而进一步完成对相关弱目标的跟踪逼近，以期在满足实时性的前提 下，完成对弱目标的跟踪。为了进一步描述该问题，考虑分布式面目标。对于该类目标，雷达 探测到的强散射点往往不是导弹的感兴趣跟踪点(弱散射点)，因此需要根据强散射点与弱 散射点的相对位置关系，借助强散射点的量测与预测信息，完成制导控制，以实现对弱散射 点的跟踪逼近。具体包括以下几个基本步骤：

步骤一、利用单脉冲雷达，实现对当前时刻强散射点目标的定位量测；

步骤二、利用无味卡尔曼滤波技术，实现对下一时刻强散射点目标的定位预测；

步骤三、利用强散射目标与弱散射目标的相对位置关系(先验知识)，实现对下一 时刻弱散射点目标的定位预测；

步骤四、根据弱散射点目标预测结果对跟踪器飞行方向的调整，完成单位时间跟 踪逼近；

返回步骤一进行迭代，可实现基于强散射点目标探测的弱散射点目标跟踪逼近。

在上述四个基本步骤的基础上，为了减小单脉冲雷达测角系统误差对跟踪性能的 影响，步骤四完成单位时间跟踪后，可根据强散射点目标预测结果完成对单脉冲雷达姿态 的调整。同时，步骤四中跟踪器飞行方向的调整(轨迹规划)可遵循追踪法或比例法两种方 案。

本发明的优点在于：

(1)实现弱点目标的实时跟踪；

(2)有效地避免了弱点目标的探测或提取过程，简化跟踪实现。

附图说明

图1是本发明的方法流程图；

图2是相邻时刻跟踪器坐标系相对位置示意图；

图3是相邻时刻雷达坐标系相对位置示意图；

图4是本发明的初始化几何关系示意图；

图5是跟踪轨迹仿真结果图；

图6是100次MonteCarlo仿真跟踪误差半径(均方根误差)曲线图；

图7是100次MonteCarlo仿真偏轴角均值变化对比图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种基于强散射点目标探测的弱散射点目标跟踪逼近方法，流程如图1 所示，包括以下几个步骤：

步骤一、利用单脉冲雷达，实现对k时刻强散射点目标的定位量测。

k时刻，单脉冲雷达测得强散射点在雷达坐标系下的测量值为sT_R_m(k)：

sT_R_m(k)＝[R_kα_kβ_k]^T(1)

其中，R_k为强散射点在雷达坐标系下的径向距离，α_k为强散射点在雷达坐标系下的 量测方位角，β_k为强散射点在雷达坐标系下的量测俯仰角。

步骤二、利用无味卡尔曼滤波技术，实现对k+1时刻强散射点目标的定位预测。

将sT_R_m(k)、sT_G_e(k-1)入下面步骤，完成强点的定位预测：

其中，sT_G_e(k-1)为地面坐标系下，k-1时刻对k时刻强散射点位置的预测值。

具体包括以下步骤：

(1)设置k时刻对称分布sigma点集

$\left(\begin{array}{cc}{\zeta }_k^{\left(0\right)}=sT\_G\_e(k-1)& \\ {\zeta }_k^{\left(i\right)}=sT\_G\_e(k-1)+{\left(\sqrt{(n+\lambda )P_{k|k-1}}\right)}_i,& i=1,2,...,n\\ {\zeta }_k^{\left(i\right)}=sT\_G\_e(k-1)-{\left(\sqrt{(n+\lambda )P_{k|k-1}}\right)}_{i-n},& i=n+1,n+2,...,2n\end{array}\right)---\left(2\right)$

并设置对应的两组权系数如下：

$\left(\begin{array}{c}{\omega }_0^{\left(m\right)}=\lambda /(n+\lambda )\\ {\omega }_0^{\left(c\right)}=\lambda /(n+\lambda )+(1+{\alpha }^2+\beta )\\ {\omega }_i^{\left(m\right)}={\omega }_i^{\left(c\right)}=0.5/(n+\lambda ),i=1,2,...,2n\end{array}\right)---\left(3\right)$

式(2)中，P_k|k-1为sT_G_e(k-1)的协方差阵。式(2)、式(3)中，λ＝α²n-n，n为sT_G_e (k-1)的维数，α为决定sigma点的散布程度的参数，通常取一个小的正值(如0.01)，β为用来 描述x的分布信息的参数，高斯情况下β的最优取值为2。

(2)进行强散射点的定位一步更新准备

$\left(\begin{array}{c}{\zeta }_{k|k-1}^{\left(i\right)}=h_k\left({\zeta }_k^{\left(i\right)}\right),i=0,1...2n\\ {\hat{Z}}_{k|k-1}=\underset{i=0}{\overset{2n}{\Sigma }}{\omega }_i^{\left(m\right)}{\zeta }_{k|k-1}^{\left(i\right)}\\ P_{{\tilde{Z}}_k}=\underset{i=0}{\overset{2n}{\Sigma }}{\omega }_i^{\left(c\right)}({\zeta }_{k|k-1}^{\left(i\right)}-{\hat{Z}}_{k|k-1}){({\zeta }_{k|k-1}^{\left(i\right)}-{\hat{Z}}_{k|k-1})}^T+C_k\\ P_{{\tilde{X}}_k{\tilde{Z}}_k}=\underset{i=0}{\overset{2n}{\Sigma }}{\omega }_i^{\left(c\right)}({\zeta }_k^{\left(i\right)}-sT\_G\_e(k-1)){({\zeta }_{k|k-1}^{\left(i\right)}-{\hat{Z}}_{k|k-1})}^T\end{array}\right)---\left(4\right)$

其中，为sigma点集通过非线性量测函数h_k()的映射值，为对k时刻 量测的估计值，为对k时刻量测的预测误差协方差阵，为对k时刻量测和状态的预测 误差互协方差矩阵，C_k为k时刻量测噪声协方差矩阵。对于非线性映射h_k(·)，有 $h_k\left({\zeta }_k^{\left(i\right)}\right)=h_{k_2}\left(h_{k_1}\right({\zeta }_k^{\left(i\right)}\left)\right).$设${\zeta }_k^{\left(i\right)}={\left(\begin{array}{cccccc}x& y& z& v_x& v_y& v_z\end{array}\right)}_k^T,$(x为第i个sigma点的x坐标，y为第 i个sigma点的y坐标，z为第i个sigma点的z坐标，v_x为第i个sigma点的x方向速率，v_y为第i个 sigma点的y方向速率，v_z为第i个sigma点的z方向速率)，有

$h_{k_1}\left({\zeta }_k^{\left(i\right)}\right)=\left(\begin{array}{c}G2R\left(k\right)\times {\left(\begin{array}{c}x\\ y\\ z\end{array}\right)}_k+G_0\_R\left(k\right)\\ G2R\left(k\right)\times {\left(\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right)}_k\end{array}\right)={\left(\begin{array}{c}{x}_r\\ y_r\\ z_r\\ v_{xr}\\ v_{yr}\\ v_{zr}\end{array}\right)}_k---\left(5\right)$

$h_k\left({\zeta }_k^{\left(i\right)}\right)=h_{k_2}\left(h_{k_1}\right({\zeta }_k^{\left(i\right)}\left)\right)={\left(\begin{array}{c}\sqrt{{x_r}^2+{y_r}^2+{z_r}^2}\\ {\tan }^{-1}(y_r/z_r)\\ {\tan }^{-1}(x_r/z_r)\end{array}\right)}_k---\left(6\right)$

其中，${\left(\begin{array}{cccccc}{x}_r& y_r& z_r& v_{xr}& v_{yr}& v_{zr}\end{array}\right)}_k^T$为中间变量，G2R(k)表示k时刻地面坐标系与雷 达坐标系的旋转矩阵，G₀_R(k)表示地面坐标系原点在k时刻雷达坐标系下的值。且强散射 点的定位一步更新如下

$\left(\begin{array}{c}{K}_k=P_{{\tilde{X}}_k{\tilde{Z}}_k}{P}_{{\tilde{Z}}_k}^{-1}\\ {\hat{X}}_{k|k}={\hat{X}}_{k|k-1}+K_k(sT\_R\_m(k)-{\hat{Z}}_{k|k-1})\\ P_{k|k}=P_{k|k-1}-K_k{P}_{{\tilde{Z}}_k}{K}_k^T\end{array}\right)---\left(7\right)$

其中，K_k为k时刻的卡尔曼增益，为k时刻状态更新值，P_k|k为的协方差矩 阵，见式(4)；

最后，完成强散射点的定位一步预测，得到一步预测值和预测误差的协方差阵

${\hat{X}}_{k+1|k}=F_k{\hat{X}}_{k|k}={\left(\begin{array}{cccccc}{\hat{x}}_{k+1}& {\hat{y}}_{k+1|k}& {\hat{z}}_{k+1|k}& {\left({\hat{v}}_x\right)}_{k+1|k}& {\left({\hat{v}}_y\right)}_{k+1|k}& {\left({\hat{v}}_z\right)}_{k+1|k}\end{array}\right)}^T---\left(8\right)$

$P_{k+1|k}=F_k{P}_{k|k}{F}_k^T+{\Gamma }_k{Q}_k{\Gamma }_k^T---\left(9\right)$

其中，F_k是系统状态转移矩阵，Q_k是过程演化噪声协方差矩阵，Γ_k是噪声矩阵， 为对k+1时刻x坐标的估计，为对k+1时刻y坐标的估计，为对k+1时刻z坐标的估 计，为对k+1时刻x方向速率的估计，为对k+1时刻y方向速率的估计，为 对k+1时刻z方向速率的估计。

得到k时刻对k+1时刻强散射点目标的定位预测sT_G_e(k)为：

$sT\_G\_e\left(k\right)=\left(\begin{array}{c}{\hat{x}}_{k+1|k}\\ {\hat{y}}_{k+1|k}\\ {\hat{z}}_{k+1|k}\end{array}\right)---\left(10\right)$

步骤三、利用强散射目标与弱散射目标的相对位置关系wT_sT(先验知识)，实现对 下一时刻弱散射点目标的定位预测。

将步骤二中得到的sT_G_e(k)转换到k时刻的雷达坐标系下，得到k时刻雷达坐标 系下，对k+1时刻强散射点位置的预测值sT_R_e(k)：

sT_R_e(k)＝G2R(k)×sT_G_e(k)+G₀_R(k)(11)

利用已知的wT_sT，得到k时刻雷达坐标系下对k+1时刻弱散射点位置的估计值：

wT_R_e(k)＝G2R(k)×wT_sT+sT_R_e(k)(12)

其中，wT_sT为弱点与强点的相对位置，是由先验知识求得的已知值。

转换坐标系，得到k时刻跟踪器坐标系下，k+1时刻弱点位置的估计值

wT_M_e(k)＝R2M(k)×wT_R_e(k)(13)

其中，R2M(k)表示k时刻雷达坐标系与跟踪器坐标系的旋转矩阵。

步骤四、根据弱散射点目标预测结果对跟踪器飞行方向的调整，完成单位时间跟 踪逼近。

根据步骤三中的估计值，求解k时刻预测偏转角

(1)追踪法$(\overrightarrow{oz}={\left(\begin{array}{ccc}0& 0& 1\end{array}\right)}^T)$

(2)比例法(factorK为比例系数)(此处需引入k-1时刻结果，因此在跟踪初始时刻 可先采用追踪法初始化轨迹，再采用比例法波进行跟踪)

其中，wT_M_k-1_e(k)为k时刻跟踪器坐标系原点在k-1时刻跟踪器坐标系中的坐标

wT_M_k-1_e(k)＝M2M^-1(k)×(wT_M_e(k)+S)(16)

其中，M2M(k)表示k-1时刻与k时刻跟踪器坐标系的旋转矩阵，S为k-1时刻跟踪器 在k时刻跟踪器坐标系下的坐标S＝[00-VrT]^T，Vr表示跟踪器飞行速率，设跟踪过程中 保持不变，T表示单脉冲测角时间间隔。

设ω_m为跟踪器最大角速度，ψ为T时间内跟踪器可偏转的最大角度，有ψ＝ω_mT。

比较与ψ。若在单位时间内，调整跟踪器飞行方向，使其在跟踪平面(由 当前预测wT_M_e(k)与跟踪器速度方向确定的平面)内转过ψ角度。若在单位时间内， 调整跟踪器飞行方向，使其在跟踪平面内转过角度。至此，完成单位时间跟踪逼近。

参看图2，图中所示为相邻时刻跟踪器坐标系相对位置示意图。图中O-XYZ_M(k)为k 时刻跟踪器坐标系，O-XYZ_M(k+1)为k+1时刻跟踪器坐标系，O-XYZ_M(k+1)由O-XYZ_M(k)在跟踪 平面旋转或ψ角得到，即O-XYZ_M(k)先绕Y轴逆时针旋转ξ₁角，再绕X轴顺时针旋转ξ₂角可得 到O-XYZ_M(k+1)坐标系。可得，k时刻跟踪器坐标系与k+1时刻跟踪器坐标系的转换矩阵更新 为

$M2M(k+1)=\left(\begin{array}{ccc}1& 0& 0\\ 0& cos(-{\xi }_2)& sin(-{\xi }_2)\\ 0& -sin(-{\xi }_2)& cos(-{\xi }_2)\end{array}\right)\left(\begin{array}{ccc}{\mathrm{cos\xi }}_1& 0& -{\mathrm{sin\xi }}_1\\ 0& 1& 0\\ {\mathrm{sin\xi }}_1& 0& {\mathrm{cos\xi }}_1\end{array}\right)---\left(17\right)$

为了减小单脉冲雷达测角系统误差对跟踪性能的影响，完成单位时间跟踪后，单 脉冲雷达姿态的调整有两种选择：

(a)根据强散射点目标预测结果进行调整

参看图3，图3为相邻时刻雷达坐标系相对位置示意图,k时刻雷达坐标系可由k+1 时刻雷达坐标系先绕Y轴逆时针旋转θ₁角，再绕X轴顺时针旋转θ₂角得到。k时刻雷达坐标系 与k+1时刻雷达坐标系的转换矩阵更新为

R2R(k+1)＝B_X(-θ₂)B_Y(θ₁)(18)

(b)不进行单独调整

R2R(k+1)＝M2M(k+1)(19)

此外，k+1时刻地面坐标系到雷达坐标系的旋转矩阵更新为：

G2R(k+1)＝R2R(k+1)G2R(k)(20)

同理，k+1时刻地面坐标系到跟踪器坐标系的旋转矩阵更新为：

G2M(k+1)＝M2M(k+1)G2M(k)(21)

通过矩阵计算，得k+1时刻雷达坐标系到跟踪器坐标系的旋转矩阵为

R2M(k+1)＝G2M(k+1)×G2R(k+1)^-1(22)

地面坐标系原点在k+1时刻雷达坐标系下的坐标可由以下迭代式求得：

G₀_R(k+1)＝R2R(k+1)×G₀_R(k)+S(23)

至此，返回步骤一进行迭代，当R_k<R_max(最小跟踪距)时，停止迭代，如此即可实现 基于强散射点目标探测的弱散射点目标跟踪逼近。实施例：

k＝1时：

步骤一、单脉冲雷达测得强散射点在雷达坐标系下的测量值(角坐标系表示)为：

sT_R_m(1)＝[R₁α₁β₁]^T(24)

步骤二、将Z₁＝sT_R_m(1)、${\hat{X}}_{1|0}=sT\_G\_e\left(0\right)={\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\end{array}\right)}^{\prime }$代入下面步骤， 完成强点的定位预测。

首先，设置对称分布sigma点集

$\left(\begin{array}{cc}{\zeta }_1^{\left(0\right)}={\hat{X}}_{1|0}& \\ {\zeta }_1^{\left(i\right)}={\hat{X}}_{1|0}+{\left(\sqrt{(n+\lambda )P_{1|0}}\right)}_i,& i=1,2,...,n\\ {\zeta }_1^{\left(i\right)}={\hat{X}}_{1|0}-{\left(\sqrt{(n+\lambda )P_{1|0}}\right)}_{i-n},& i=n+1,n+2,...,2n\end{array}\right)---\left(25\right)$

设置权系数取值如下：

$\left(\begin{array}{c}{\omega }_0^{\left(m\right)}=\lambda /(n+\lambda )\\ {\omega }_0^{\left(c\right)}=\lambda /(n+\lambda )+(1-{\alpha }^2+\beta )\\ {\omega }_i^{\left(m\right)}={\omega }_i^{\left(c\right)}=0.5/(n+\lambda ),i=1,2,...,2n\end{array}\right)---\left(26\right)$

其中，P_1|0＝10⁶I₆，I₆为六维单位阵，λ＝α²(n+κ)-n，n＝6，α＝0.01，κ＝0，β＝2。

然后，完成强散射点的定位一步更新准备

$\left(\begin{array}{c}{\zeta }_{1|0}^{\left(i\right)}=h_1\left({\zeta }_1^{\left(i\right)}\right),i=0,1...2n\\ {\hat{Z}}_{1|0}=\underset{i=0}{\overset{2n}{\Sigma }}{\omega }_i^{\left(m\right)}{\zeta }_{1|0}^{\left(i\right)}\\ P_{{\tilde{Z}}_1}=\underset{i=0}{\overset{2n}{\Sigma }}{\omega }_i^{\left(c\right)}({\zeta }_{1|0}^{\left(i\right)}-{\hat{Z}}_{1|0}){({\zeta }_{1|0}^{\left(i\right)}-{\hat{Z}}_{1|0})}^T+R_k\\ P_{{\tilde{X}}_1{\tilde{Z}}_1}=\underset{i=0}{\overset{2n}{\Sigma }}{\omega }_i^{\left(c\right)}({\zeta }_1^{\left(i\right)}-{\hat{X}}_{1|0}){({\zeta }_{1|0}^{\left(i\right)}-{\hat{Z}}_{1|0})}^T\end{array}\right)---\left(27\right)$

其中，$R_k=\left(\begin{array}{ccc}{\sigma }_{21}^2& 0& 0\\ 0& {\sigma }_{22}^2& 0\\ 0& 0& {\sigma }_{22}^2\end{array}\right),{\sigma }_{21}^2=1,{\sigma }_{22}^2=0.01,h_1\left({\zeta }_1^{\left(i\right)}\right)=h_{1_2}\left(h_{1_1}\right({\zeta }_1^{\left(i\right)}\left)\right)$为量测映 射函数，且${\zeta }_1^{\left(i\right)}={\left(\begin{array}{cccccc}x& y& z& v_x& v_y& v_z\end{array}\right)}_1^{T\left(i\right)},$有

$h_{1_1}\left({\zeta }_1^{\left(i\right)}\right)=\left(\begin{array}{c}G2R\left(1\right)\times {\left(\begin{array}{c}x\\ y\\ z\end{array}\right)}_1^{\left(i\right)}+G_0\_R\left(1\right)\\ G2R\left(1\right)\times {\left(\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right)}_1^{\left(i\right)}\end{array}\right)={\left(\begin{array}{c}{x}_r\\ y_r\\ z_r\\ v_{xr}\\ v_{yr}\\ v_{zr}\end{array}\right)}_1^{\left(i\right)}---\left(28\right)$

$h_k\left({\zeta }_1^{\left(i\right)}\right)=h_{k_2}\left(h_{k_1}\right({\zeta }_1^{\left(i\right)}\left)\right)={\left(\begin{array}{c}\sqrt{{x_r}^2+{y_r}^2+{z_r}^2}\\ {\tan }^{-1}(y_r/z_r)\\ {\tan }^{-1}(x_r/z_r)\end{array}\right)}_1^{\left(i\right)}---\left(29\right)$

G2R(1)表示k＝1时刻，地面坐标系与雷达坐标系的旋转矩阵，此时位置关系如图4 所示，图中O_R/MZ_R/M轴、O_R/MX_R/M轴、OZ轴共纸面，此时跟踪器速度方向为O_R/MZ_R/M轴方向，跟踪器 速度方向与跟踪器-目标连线的夹角(Δ＝5°)、跟踪器-目标连线与OZ轴的夹角(θ＝87°)、 O_R/MY_R/M轴与目标初始速度方向OY轴的夹角(γ＝60°)、跟踪器-目标初始距(R₀＝1000m)为 已知量。地面坐标系与雷达(跟踪器)坐标系的转换矩阵初始化如下

$G2R\left(1\right)=\left(\begin{array}{ccc}\cos (\theta -\Delta )& 0& -\sin (\theta -\Delta )\\ 0& 1& 0\\ \sin (\theta -\Delta )& 0& \cos (\theta -\Delta )\end{array}\right)\left(\begin{array}{ccc}\cos (-\gamma )& \sin (-\gamma )& 0\\ -\sin (-\gamma )& \cos (-\gamma )& 0\\ 0& 0& 1\end{array}\right)---\left(30\right)$

G₀_R(1)表示k＝1时刻，地面坐标系原点在k时刻雷达坐标系下的值，且有

$G_0\_R\left(1\right)=\left(\begin{array}{c}tan{\beta }_1\times R_1/\sqrt{1+{\tan }^2{\alpha }_1+{\tan }^2{\beta }_1}\\ {\mathrm{tan\alpha }}_1\times R_1/\sqrt{1+{\tan }^2{\alpha }_1+{\tan }^2{\beta }_1}\\ R_1/\sqrt{1+{\tan }^2{\alpha }_1+{\tan }^2{\beta }_1}\end{array}\right)---\left(31\right)$

强散射点的定位一步更新如下

$\left(\begin{array}{c}{K}_k=P_{{\tilde{X}}_k{\tilde{Z}}_k}{P}_{{\tilde{Z}}_k}^{-1}\\ {\hat{X}}_{k|k}={\hat{X}}_{k|k-1}+K_k(Z_k-{\hat{Z}}_{k|k-1})\\ P_{k|k}=P_{k|k-1}-K_k{P}_{{\tilde{Z}}_k}{K}_k^T\end{array}\right)---\left(32\right)$

最后，完成强散射点的定位一步预测，得到一步预测值和预测误差的协方差阵

${\hat{X}}_{k+1|k}=F_k{\hat{X}}_{k|k}={\left(\begin{array}{cccccc}{\hat{x}}_{k+1|k}& {\hat{y}}_{k+1|k}& {\hat{z}}_{k+1|k}& {\left({\hat{v}}_x\right)}_{k+1|k}& {\left({\hat{v}}_y\right)}_{k+1|k}& {\left({\hat{v}}_z\right)}_{k+1|k}\end{array}\right)}^T---\left(33\right)$

$P_{k+1|k}=F_k{P}_{k|k}{F}_k^T+{\Gamma }_k{Q}_k{\Gamma }_k^T---\left(34\right)$

其中，$F_k=\left(\begin{array}{cccccc}1& 0& 0& T& 0& 0\\ 0& 1& 0& 0& T& 0\\ 0& 0& 1& 0& 0& T\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),{\Gamma }_k=\left(\begin{array}{ccc}{T}^2/2& 0& 0\\ 0& T^2/2& 0\\ 0& 0& T^2/2\\ T& 0& 0\\ 0& T& 0\\ 0& 0& T\end{array}\right),T=0.05,$$Q_k=\left(\begin{array}{ccc}{\sigma }_{11}^2& 0& 0\\ 0& {\sigma }_{11}^2& 0\\ 0& 0& {\sigma }_{12}^2\end{array}\right),{\sigma }_{11}^2=1,{\sigma }_{12}^2=\mathrm{0.01.}$

至此，得到k＝1时刻对k＝2时刻强散射点目标的定位预测

$sT\_G\_e\left(1\right)=\left(\begin{array}{c}{\hat{x}}_{2|1}\\ {\hat{y}}_{2|1}\\ {\hat{z}}_{2|1}\end{array}\right)---\left(\mathrm{35}\right)$

步骤三、实现对下一时刻弱散射点目标的定位预测。

将步骤二中得到的sT_G_e(1)转换到k＝1时刻的雷达坐标系下，得到k＝1时刻雷 达坐标系下，对k＝2时刻强散射点位置的预测值：

sT_R_e(1)＝G2R(1)×sT_G_e(1)+G₀_R(1)(36)

利用已知的弱点与强点的相对位置wT_sT＝[LL0]^T,L＝10，得到k＝1时刻雷达 坐标系下对k＝2时刻弱散射点位置的估计值

wT_R_e(1)＝G2R(1)×wT_sT+sT_R_e(1)(37)

转换坐标系，得到k＝1时刻跟踪器坐标系下，k＝2时刻弱点位置的估计值

wT_M_e(1)＝R2M(1)×wT_R_e(1)(38)

其中，$R2M\left(1\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right).$

步骤四、根据弱散射点目标预测结果对跟踪器飞行方向的调整，完成单位时间跟 踪逼近。

根据步骤三中的估计值，求解k＝1时刻预测偏转角

k＝1时，采用追踪法$(\overrightarrow{oz}={\left(\begin{array}{ccc}0& 0& 1\end{array}\right)}^T)$

至此，比较与ψ(ψ＝ω_mT，ω_m＝5°/s，T＝0.05s)。在单位时间内，调整跟踪 器飞行方向，使其在跟踪平面内转过ψ角度。

参看图2，k＝1时刻跟踪器坐标系与k＝2时刻跟踪器坐标系的转换矩阵更新为

$M2M\left(2\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& cos(-{\xi }_2)& sin(-{\xi }_2)\\ 0& -sin(-{\xi }_2)& cos(-{\xi }_2)\end{array}\right)\left(\begin{array}{ccc}{\mathrm{cos\xi }}_1& 0& -{\mathrm{sin\xi }}_1\\ 0& 1& 0\\ {\mathrm{sin\xi }}_1& 0& {\mathrm{cos\xi }}_1\end{array}\right)---\left(40\right)$

根据强散射点目标预测结果进行调整，参看图3，k＝1时刻雷达坐标系与k＝2时刻 雷达坐标系的转换矩阵更新为

R2R(2)＝B_X(-θ₂)B_Y(θ₁)(41)

此外，k＝2时刻地面坐标系到雷达坐标系的旋转矩阵更新为：

G2R(2)＝R2R(2)G2R(1)(42)

同理，k+1时刻地面坐标系到跟踪器坐标系的旋转矩阵更新为：

G2M(2)＝M2M(2)G2M(1)＝M2M(2)G2R(1)(43)

通过矩阵计算，得k+1时刻雷达坐标系到跟踪器坐标系的旋转矩阵为

R2M(2)＝G2M(2)G2R(2)^-1(44)

地面坐标系原点在k+1时刻雷达坐标系下的坐标可由以下迭代式求得：

G₀_R(2)＝R2R(2)×G₀_R(1)+S(45)

S＝[00-VrT]^T，Vr＝150m/s，，T＝0.05s。至此，完成单位时间跟踪逼近，返回步 骤一迭代。

k>1，且R_k<R_max时(R_max＝100m)：

步骤一、单脉冲雷达测得强散射点在雷达坐标系下的测量值(角坐标系表示)为：

sT_R_m(k)＝[R_kα_kβ_k]^T(46)

步骤二、将Z_k＝sT_R_m(k)、入下面步骤，完成强点的定位预 测：

首先，设置对称分布sigma点集

$\left(\begin{array}{cc}{\zeta }_k^{\left(0\right)}={\hat{X}}_{k|k-1}& \\ {\zeta }_k^{\left(i\right)}={\hat{X}}_{k|k-1}+{\left(\sqrt{(n+\lambda )P_{k|k-1}}\right)}_i,& i=1,2,...,n\\ {\zeta }_k^{\left(i\right)}={\hat{X}}_{k|k-1}-{\left(\sqrt{(n+\lambda )P_{k|k-1}}\right)}_{i-n},& i=n+1,n+2,...,2n\end{array}\right)---\left(47\right)$

设置权系数取值如下：

$\left(\begin{array}{c}{\omega }_0^{\left(m\right)}=\lambda /(n+\lambda )\\ {\omega }_0^{\left(c\right)}=\lambda /(n+\lambda )+(1-{\alpha }^2+\beta )\\ {\omega }_i^{\left(m\right)}={\omega }_i^{\left(c\right)}=0.5/(n+\lambda ),i=1,2,...,2n\end{array}\right)---\left(48\right)$

其中，λ＝α²(n+κ)-n，n＝6，α＝0.01，κ＝0，β＝2。

然后，完成强散射点的定位一步更新准备

$\left(\begin{array}{c}{\zeta }_{k|k-1}^{\left(i\right)}=h_k\left({\zeta }_k^{\left(i\right)}\right),i=0,1...2n\\ {\hat{Z}}_{k|k-1}=\underset{i=0}{\overset{2n}{\Sigma }}{\omega }_i^{\left(m\right)}{\zeta }_{k|k-1}^{\left(i\right)}\\ P_{{\tilde{Z}}_k}=\underset{i=0}{\overset{2n}{\Sigma }}{\omega }_i^{\left(c\right)}({\zeta }_{k|k-1}^{\left(i\right)}-{\hat{Z}}_{k|k-1}){({\zeta }_{k|k-1}^{\left(i\right)}-{\hat{Z}}_{k|k-1})}^T+R\\ P_{{\tilde{X}}_k{\tilde{Z}}_k}=\underset{i=0}{\overset{2n}{\Sigma }}{\omega }_i^{\left(c\right)}({\zeta }_k^{\left(i\right)}-{\hat{X}}_{k|k-1}){({\zeta }_{k|k-1}^{\left(i\right)}-{\hat{Z}}_{k|k-1})}^T\end{array}\right)---\left(49\right)$

其中，$R_k=\left(\begin{array}{ccc}{\sigma }_{21}^2& 0& 0\\ 0& {\sigma }_{22}^2& 0\\ 0& 0& {\sigma }_{22}^2\end{array}\right),{\sigma }_{21}^2=1,{\sigma }_{22}^2=0.01,h_k\left({\zeta }_k^{\left(i\right)}\right)=h_{k_2}\left(h_{k_1}\right({\zeta }_k^{\left(i\right)}\left)\right)$为量测映 射函数，且${\zeta }_k^{\left(i\right)}={\left(\begin{array}{cccccc}x& y& z& v_x& v_y& v_z\end{array}\right)}_k^T,$有

$h_{k_1}\left(X_k\right)=\left(\begin{array}{c}G2R\left(k\right)\times {\left(\begin{array}{c}x\\ y\\ z\end{array}\right)}_k+G_0\_R\left(k\right)\\ G2R\left(k\right)\times {\left(\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right)}_k\end{array}\right)={\left(\begin{array}{c}{x}_r\\ y_r\\ z_r\\ v_{xr}\\ v_{yr}\\ v_{zr}\end{array}\right)}_k---\left(50\right)$

$h_k\left(X_k\right)=h_{k_2}\left(h_{k_1}\right(X_k\left)\right)={\left(\begin{array}{c}\sqrt{{x_r}^2+{y_r}^2+{z_r}^2}\\ {\tan }^{-1}(y_r/z_r)\\ {\tan }^{-1}(x_r/z_r)\end{array}\right)}_k---\left(51\right)$

G2R(k)表示k时刻地面坐标系与雷达坐标系的旋转矩阵，G₀_R(k)表示地面坐标系 原点在k时刻雷达坐标系下的值。且强散射点的定位一步更新如下

$\left(\begin{array}{c}{K}_k=P_{{\tilde{X}}_k{\tilde{Z}}_k}{P}_{{\tilde{Z}}_k}^{-1}\\ {\hat{X}}_{k|k}={\hat{X}}_{k|k-1}+K_k(Z_k-{\hat{Z}}_{k|k-1})\\ P_{k|k}=P_{k|k-1}-K_k{P}_{{\tilde{Z}}_k}{K}_k^T\end{array}\right)---\left(52\right)$

最后，完成强散射点的定位一步预测，得到一步预测值和预测误差的协方差阵

${\hat{X}}_{k+1|k}=F_k{\hat{X}}_{k|k}={\left(\begin{array}{cccccc}{\hat{x}}_{k+1}& {\hat{y}}_{k+1|k}& {\hat{z}}_{k+1|k}& {\left({\hat{v}}_x\right)}_{k+1|k}& {\left({\hat{v}}_y\right)}_{k+1|k}& {\left({\hat{v}}_z\right)}_{k+1|k}\end{array}\right)}^T---\left(53\right)$

$P_{k+1|k}=F_k{P}_{k|k}{F}_k^T+{\Gamma }_k{Q}_k{\Gamma }_k^T---\left(54\right)$

其中，$F_k=\left(\begin{array}{cccccc}1& 0& 0& T& 0& 0\\ 0& 1& 0& 0& T& 0\\ 0& 0& 1& 0& 0& T\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),{\Gamma }_k=\left(\begin{array}{ccc}{T}^2/2& 0& 0\\ 0& T^2/2& 0\\ 0& 0& T^2/2\\ T& 0& 0\\ 0& T& 0\\ 0& 0& T\end{array}\right),T=0.05,$$Q_k=\left(\begin{array}{ccc}{\sigma }_{11}^2& 0& 0\\ 0& {\sigma }_{11}^2& 0\\ 0& 0& {\sigma }_{12}^2\end{array}\right),{\sigma }_{11}^2=1,{\sigma }_{12}^2=\mathrm{0.01.}$

至此，得到k时刻对k+1时刻强散射点目标的定位预测

$sT\_G\_e\left(k\right)=\left(\begin{array}{c}{\hat{x}}_{k+1|k}\\ {\hat{y}}_{k+1|k}\\ {\hat{z}}_{k+1|k}\end{array}\right)---\left(55\right)$

步骤三、实现对下一时刻弱散射点目标的定位预测。

将步骤二中得到的sT_G_e(k)转换到k时刻的雷达坐标系下，得到k时刻雷达坐标 系下，对k+1时刻强散射点位置的预测值：

sT_R_e(k)＝G2R(k)×sT_G_e(k)+G₀_R(k)(56)

利用已知的弱点与强点的相对位置wT_sT＝[LL0]^T,L＝10，得到k时刻雷达坐标 系下对k+1时刻弱散射点位置的估计值

wT_R_e(k)＝G2R(k)×wT_sT+sT_R_e(k)(57)

转换坐标系，得到k时刻跟踪器坐标系下，k+1时刻弱点位置的估计值

wT_M_e(k)＝R2M(k)×wT_R_e(k)(58)

步骤四、根据弱散射点目标预测结果对跟踪器飞行方向的调整，完成单位时间跟 踪逼近。

根据步骤三中的估计值，采用比例法(factorK＝3)求解k时刻预测偏转角

其中，wT_M_k-1_e(k)＝M2M^-1(k)×(wT_M_e(k)+S)，S＝[00-VrT]^T，Vr＝150m/s，T ＝0.05s。

至此，比较与ψ(ψ＝ω_mT，ω_m＝5°/s，T＝0.05s)。若在单位时间内，调整跟 踪器飞行方向，使其在跟踪平面内转过ψ角度；若在单位时间内，调整跟踪器飞行方 向，使其在跟踪平面内转过角度。

参看图2，k时刻跟踪器坐标系与k+1时刻跟踪器坐标系的转换矩阵更新为

$M2M\left(2\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& cos(-{\xi }_2)& sin(-{\xi }_2)\\ 0& -sin(-{\xi }_2)& cos(-{\xi }_2)\end{array}\right)\left(\begin{array}{ccc}{\mathrm{cos\xi }}_1& 0& -{\mathrm{sin\xi }}_1\\ 0& 1& 0\\ {\mathrm{sin\xi }}_1& 0& {\mathrm{cos\xi }}_1\end{array}\right)---\left(60\right)$

根据强散射点目标预测结果进行调整，参看图3，k时刻雷达坐标系与k+1时刻雷达 坐标系的转换矩阵更新为

R2R(k+1)＝B_X(-θ₂)B_Y(θ₁)(61)

此外，k+1时刻地面坐标系到雷达坐标系的旋转矩阵更新为：

G2R(k+1)＝R2R(k+1)G2R(k)(62)

同理，k+1时刻地面坐标系到跟踪器坐标系的旋转矩阵更新为：

G2M(k+1)＝M2M(k+1)G2M(k)(63)

通过矩阵计算，得k+1时刻雷达坐标系到跟踪器坐标系的旋转矩阵为

R2M(k+1)＝G2M(k+1)×G2R(k+1)^-1(64)

地面坐标系原点在k+1时刻雷达坐标系下的坐标可由以下迭代式求得：

G₀_R(k+1)＝R2R(k+1)×G₀_R(k)+S(65)

至此，返回步骤一进行迭代k＝k+1，当R_k<R_max(最小跟踪距)时，停止迭代，如此即 可实现基于强散射点目标探测的弱散射点目标跟踪逼近。

仿真验证：

将本实施案例注入实际仿真系统中，步骤四中采用比例法，且根据强散射点目标 预测结果调整雷达姿态，得到的跟踪轨迹如图5所示，进行100次MonteCarlo仿真的跟踪误 差曲线如图6所示，平均最终跟踪误差(均方根误差)半径为：

Δ＝5m

采用同样参数，对比步骤四中：

采用比例法，根据强散射点目标预测结果调整雷达姿态

采用跟踪法，根据强散射点目标预测结果调整雷达姿态

两种方案的跟踪结果，进行100次MonteCarlo仿真得到的偏轴角(跟踪器速度方 向与跟踪器-目标连线的夹角)变化曲线如图7。可见，跟踪法对跟踪器的机动特性要求更 高，比例法的跟踪机动调整更为平缓。

通过改变弱点与强点的相对位置取值，对比步骤四中：

采用比例法，根据强散射点目标预测结果调整雷达姿态

采用比例法，不根据强散射点目标预测结果调整雷达姿态

两种方案的跟踪结果，平均最终跟踪误差(均方根误差)半径对比见表3

表3平均最终跟踪误差(均方根误差)半径

可见，在该组参数下，当弱点与强点相隔的时候，根据强散射点目标预测 结果调整雷达姿态的方案，可以减小单脉冲雷达测角系统误差对跟踪性能的影响。相较之 下，雷达姿态不调整的方案跟踪性能恶化严重。

本发明主要针对弱散射点目标的跟踪问题，利用单脉冲雷达测角测距技术以及基 于坐标系转换的无味卡尔曼滤波技术，建立了基于强散射点目标探测的弱散射点目标跟踪 逼近模型。通过实例分析，验证了模型的正确性，对比了轨迹规划过程中采用比例法和跟踪 法两种方案对跟踪器机动性能的要求，并分析了强弱点距离对调整/不调整雷达姿态两种 方案的跟踪性能的影响，给出了模型具体方案的选择依据。