一种改进SRCKF强跟踪滤波的机动目标跟踪方法与流程

本发明涉及一种针对光学卫星观测改进平方根容积卡尔曼滤波器(square-rootcubaturekalmanfilter，srckf)强跟踪滤波的机动目标跟踪方法，属于数据关联与跟踪技术领域。

背景技术：

近年来，光学图像卫星由于其观测地域广，不易受国界限制，不易受电磁干扰等优势，一直以来都是各国研究的重点，而在光学卫星观测中，对空间目标进行状态估计以及态势感知也是研究的一大重点。

但是，在现有条件下，空间目标的运动模式日趋复杂，且会产生抛洒诱饵、变轨以及分导等行为，故会偏离原有的运动模式，产生一定的机动现象，进而对目标状态估计以及态势感知造成了较大困难。

与此同时，现有的机动目标跟踪方法一般分为两类，第一类是首先针对目标运动进行建模，尔后利用交互式多模方法(interactivemulti-models,imm)对目标进行跟踪，但是此类方法需要建立较为准确的目标运动模型，而且在模型转换过程中往往需要消耗大量时间，空间目标抛洒诱饵、变轨、分导等往往是短时间的行为，且其运动模型相对难以建立，故此类方法并不适用于空间目标的跟踪。第二类则是利用自适应滤波的思路，基于平方根容积卡尔曼滤波器的强跟踪滤波方法即是其中的一种，这种方法不改变运动模型，而是通过改变滤波，使之能够对目标的运动产生适应性，但是此类方法在目标观测初期容易导致较大的误差，且对观测噪声的鲁棒性不好，并不适用于空间目标跟踪的全流程，本发明主要是针对此类方法的缺点，对其进行改进，从而使之能够适应目标不发生机动的情况且使其跟踪结果更具有鲁棒性。

技术实现要素：

本发明所要解决的技术问题是，针对现有技术对空间机动目标跟踪存在的缺点，提出一种改进srckf强跟踪滤波的机动目标跟踪方法。

为解决上述技术问题，本发明所采用的技术方案是：一种改进srckf强跟踪滤波的机动目标跟踪方法，包括以下步骤：

s1根据目标运动情况和卫星的观测情况建立运动模型和观测模型；

s1.1根据目标运动情况建立运动模型，当目标为自由段运动的导弹目标时，在不机动的条件下，可以视作仅受重力影响，故考虑使用j2摄动模型对目标运动进行建模；

目标在eci坐标系下的状态向量为x(k)＝[x(k),y(k),z(k),vx(k),vy(k),vz(k)]^t,其中x(k),y(k),z(k)、vx(k),vy(k),vz(k)分别代表目标在k时刻的位置和速度，有目标j2摄动模型如下：

其中μe＝3.98604418m³/sec²代表地球引力常数，re＝6378137m代表地球半径，j2＝0.001082626075代表摄动系数；

由于上述方程难以直接表述为关于x(k)的方程形式，故在已知k时刻目标状态为x(k)的条件下，考虑使用四阶runge-kutta数值积分方法来得到k+1时刻的值，此时有：

x(k+1)＝f(x(k))+q(k)

其中q(k)为零均值高斯矢量，t代表目标状态离散化时的时间间隔；

s1.2根据卫星的观测情况建立观测模型，卫星为被动式观测，仅仅能够测得目标的方位角和俯仰角，在eci坐标系下，假设k时刻卫星的状态向量为：

sa(k)＝[xs(k),ys(k),zs(k),vsx(k),vsy(k),vsz(k)]^t

其中xs(k),ys(k),zs(k)、vsx(k),vsy(k),vsz(k)分别代表卫星在k时刻的位置和速度；

以卫星速度方向为xb轴，以卫星到地球球心连线为zb轴，按右手螺旋定则给出yb轴，定义xb轴，yb轴，zb轴的轴矢量分别为sxb，syb，szb，给出星体坐标系，αs代表目标方位角，βs代表目标俯仰角，sx代表卫星视线，此时可由如下方式计算αs、βs：

sxb＝[vx(k),vy(k),vz(k)]^t

szb＝-[x(k),y(k),z(k)]^t

sx＝[xs(k)-x(k),ys(k)-y(k),zs(k)-z(k)]^t

为了与滤波模型切合，一般将上式写为：

其中z(k)代表卫星k时刻的量测，dβs和dαs代表传感器测量噪声，r(k)为零均值高斯矢量；

s2：滤波初始化，由于使用了srckf方法，故在滤波过程中，迭代量并非协方差矩阵，而是协方差矩阵的平方根矩阵，故在初值时需要将初值协方差矩阵进行平方根分解，即有滤波初始化过程为：

其中，代表对目标初始时刻的状态估计，e(·)代表对矩阵求期望，p(0)代表目标初始时刻状态估计值的协方差矩阵，s(0)代表目标初始时刻协方差矩阵的平方根，chol(·)代表对矩阵进行cholesky分解得到下三角矩阵；

同时由于s5中需要统计渐消因子变化次数，初始时刻，定义变化次数flag＝0；

s3：进行滤波预测，也即在已知k时刻的系统状态时，对k+1时刻的系统状态进行外推的过程，仍按照srckf方法的框架，有：

对于k∈{0,…,∞}：

计算容积点：

其中代表k时刻容积点矩阵υk的第i列，n代表状态向量的维数，代表第i个容积点，代表k时刻的系统状态估计值，s(k)代表k时刻系统状态估计值协方差矩阵的平方根，[1]i代表矩阵[1]中第i个元素，矩阵[1]如下所示：

预测过程如下：

υk+1|k＝f(υk)

s^e(k+1|k)＝qr([vk+1|ksq]^t)

其中代表利用k时刻容积点值预测出k+1时刻值容积点矩阵υk+1|k的第i列，vk+1|k代表容积点矩阵υk+1|k与k+1时刻的预测值的偏差，qr(·)即表示对矩阵进行qr分解得到的上三角矩阵求转置的函数，sq＝chol(q(k))，代表预测过程中运动模型噪声的协方差矩阵的平方根，q(k)＝q(k)q(k)^t代表k时刻预测过程中运动模型噪声的协方差矩阵；

即得到了k+1时刻未带有渐消因子的预测值及其协方差的平方根s^e(k+1|k)；

s4：进行第一次量测更新，本发明是基于强跟踪滤波方法改进而来，由于强跟踪滤波方法是在拓展卡尔曼滤波的框架下提出的，在计算渐消因子过程中，需要给出运动方程以及观测方程的雅克比矩阵，本发明为了在不给出雅可比矩阵的条件下，将其拓展到srckf方法的框架下，需要先增加一次量测更新过程，具体步骤如下：

传播容积点：

其中代表利用k+1时刻的预测值和协方差平方根得到的容积点矩阵的第i列，表示对容积点矩阵按照s1.2的观测方程进行观测，代表对进行观测的容积点矩阵ψk+1|k的第i列；

未加渐消因子时量测更新过程如下：

szz^e(k+1)＝qr([zk+1sr]^t)

pxz^e(k+1)＝xk+1zk+1^t

其中zk+1代表k+1时刻容积点矩阵ψk+1|k与预测量测值的偏差，xk+1代表k+1时刻容积点矩阵与k+1时刻的预测值的偏差，sr＝chol(r(k+1))，代表k+1时刻观测方程噪声所产生的协方差矩阵的平方根，r(k+1)＝r(k+1)r(k+1)^t即是k+1时刻观测方程噪声所产生的协方差矩阵，z(k+1)代表k+1时刻的真实量测值，γ(k+1)代表对新息增量协方差的估计值，ρ为遗忘因子，0＜ρ≤1，一般取值为ρ＝0.95。

即得到了未加渐消因子时的预测量测值及其协方差的平方根szz^e(k+1)以及预测值和预测量测值的互协方差矩阵pxz^e(k+1)，以及新息增量ε(k+1)；

s5：计算带有时间因子的渐消因子并进行门限判决，此时需要利用s3和s4得出的s^e(k+1|k)、szz^e(k+1)和pxz^e(k+1)，此时有：

hx(k+1)＝pxz^e(k+1)^t/s^e(k+1|k)/s^e(k+1|k)^t

其中norm(·)代表求矩阵各个元素的平方和，hx(k+1)代表k+1时刻对观测过程中对预测过程中运动模型噪声的协方差矩阵平方根的增益矩阵，l(k+1)即代表k+1时刻的时间因子项，式中tb代表观测时的采样间隔，ts代表卫星观测时的收敛时间，λ0代表渐消因子的计算初值，λ(k+1)代表k+1时刻的渐消因子终值；

此时若λ(k+1)＞1，则对变化次数flag加1，并将flag与预先设定的门限值nλ比较，nλ设为5，若flag大于等于nλ，则令k＝k-nλ，即是向前回溯nλ个时刻，将变化次数flag重置为0，将滤波过程中所有含有k的量都变为新的时刻所代表的值，并令若flag小于nλ，则令s(k+1|k)＝s^e(k+1|k)，继续正常滤波过程；若λ(k+1)＝1，则直接将变化次数flag置为0，同时令s(k+1|k)＝s^e(k+1|k)，继续正常滤波过程；

s6：根据s5的门限判决结果再次进行量测更新；

根据s5得出的s(k+1|k)的值，对k+1时刻得到的真实量测进行更新：

传播容积点：

量测更新过程如下：

szz(k+1)＝qr([zk+1sr]^t)

pxz(k+1)＝xk+1zk+1^t

w(k+1)＝pxz(k+1)/szz(k+1)/szz(k+1)^t

s(k+1)＝qr([xk+1-w(k+1)zk+1w(k+1)sr]^t)

其中代表k+1时刻的预测量测值，szz(k+1)代表协方差的平方根，pxz(k+1)代表k+1时刻的预测值和预测量测值的互协方差矩阵，w(k+1)表示k+1时刻的滤波增益矩阵。

即得到了k+1时刻的状态估计向量及其协方差矩阵的平方根s(k+1)，从而可以继续迭代，推导出所有时刻的状态估计值，滤波过程的算法流程图如图3所示。

与现有技术相比，本发明所具有的有益效果为：首先，无需建立较为准确的目标运动模型，仅使用目标正常运动时的运动模型也可以达到对机动目标的跟踪效果；其次，相较于一般的适应性滤波方法，本发明在目标不机动时的跟踪效果与普通不考虑机动的滤波方法相当，且受噪声干扰较少，跟踪结果更具有鲁棒性。

附图说明

图1为卫星观测模型；

图2为改进srckf强跟踪滤波的机动目标跟踪方法实施流程图；

图3为改进srckf强跟踪滤波的机动目标跟踪方法滤波过程实施流程图；

图4为目标不发生机动的位置速度rmse对比图：(a)位置误差(b)速度误差；

图5为目标发生长时间小机动的位置速度rmse对比图：(a)位置误差(b)速度误差；

图6为目标发生长时间小机动的位置速度局部放大的rmse对比图：(a)位置误差(b)速度误差；

图7为目标发生短时间大机动的位置速度rmse对比图：(a)位置误差(b)速度误差；

图8为目标发生短时间大机动的位置速度局部放大的rmse对比图：(a)位置误差(b)速度误差；

具体实施方式

下面结合附图对本发明具体实施方式作进一步说明。

本发明的技术方案是：一种改进srckf强跟踪滤波的机动目标跟踪方法，具体包括下述步骤：

s1：根据目标运动情况和卫星的观测情况建立运动模型和观测模型；

s2：滤波初始化；

s3：进行滤波预测；

s4：进行第一次量测更新；

s5：计算带有时间因子的渐消因子并进行门限判决；

s6：根据s5的门限判决结果再次进行量测更新。

为了验证本发明方法的有效性，本发明选取了三种场景，分别是目标不发生机动，发生长时间小机动以及短时间大机动，对比了本发明方法和不考虑机动的ukf滤波方法以及srckf强跟踪滤波方法，蒙特卡洛仿真次数为100次，以均方根误差(rootmeansquarederror，rmse)为效果衡量标准，给出了不同算法的位置速度rmse对比图。

图4为目标不发生机动的位置速度rmse对比图，从图中可以看出，本发明方法在目标不发生机动的条件下，速度位置收敛结果都和ukf滤波方法相同，而srckf强跟踪滤波方法由于初始时刻估计值不准确，从而误将初始时刻判定为运动状态，导致一开始收敛速度很慢，且最终收敛效果也不如另外两者，同时srckf强跟踪滤波方法收敛过程中容易受噪声影响，在100次蒙特卡洛的条件下，其结果仍有较大波动。

图5、图6分别为目标发生长时间小机动的位置速度rmse对比图以及其局部放大图，默认机动时刻为400s，机动时长为100s，机动加速度为从图中可以看出，与不考虑机动的ukf滤波方法相比，本发明方法和srckf强跟踪滤波方法都能在短时间内对目标的机动产生响应，并能够在目标机动停止后收敛到比较小的值；而与srckf强跟踪滤波相比，本发明方法在目标机动后发散速度较慢，对目标机动的跟踪精度更好，且srckf强跟踪滤波在跟踪过程中会产生波动性，本文方法更加鲁棒。

图7、图8分别为目标发生短时间大机动的位置速度rmse对比图以及其局部放大图，默认机动时刻为400s，机动时长为20s，机动加速度为从图中可以看出，与不考虑机动的ukf滤波方法相比，本发明方法和srckf强跟踪滤波方法都能在短时间内对目标的机动产生响应，并能够在目标机动停止后收敛到比较小的值；而与srckf强跟踪滤波相比，本发明方法在目标机动后发散速度较慢，对目标机动的跟踪精度更好，且srckf强跟踪滤波在跟踪过程中会产生波动性，本发明方法更加鲁棒。

从实验结果可以看出本发明方法不仅适用于机动目标跟踪，对目标不机动的状态估计也较好，表明本发明方法有效。