基于HCKS-EM的战斗机蛇形机动目标跟踪算法

卢春光， 周中良， 刘宏强， 寇添， 杨远志

(空军工程大学航空工程学院, 西安 710038)

蛇形机动是战斗机飞行员在隐蔽接敌、机动规避及协同探测等战术动作中经常采用的战术机动类型。在一定的空战态势下，飞行员进行蛇形机动反映了其战术意图，因此在空战中如何快速而又准确地识别出目标的蛇形机动模式对于明确其战术意图以及评估当前的战场态势具有十分重要的意义。目前，机动模式的识别主要从机动模式的几何特征(宏观)以及运动参量特征(微观)2个方面着手：文献[1]通过将战术机动进行分割并按时间序列进行编码，采用隐马尔可夫模型进行训练、识别；文献[2]通过提取目标航向角变化率、高度变化率等运动参量特征，设计了一个2级识别方案，采用模糊推理和时间自动机方法实现了机动模式的识别。对于蛇形机动目标而言，转弯角速度是其最为关键的运动参量，因此通过载机雷达的量测数据精确辨识出转弯角速度对于提高战斗机蛇形机动模式的识别率具有重要意义。

目前，机动目标转弯角速度辨识的问题已经得到了许多学者的关注。文献[3-5]将转弯角速度作为目标状态变量，通过重新构造状态方程，分别使用扩展卡尔曼滤波(EKF)/无迹卡尔曼滤波(UKF)、基于5阶球径容积规则的高阶容积卡尔曼滤波(HCKF)等非线性滤波算法估计角速度和其他状态变量，其缺点是角速度估计的精度依赖于目标初始状态及协方差、过程噪声等因素。文献[6]根据角速度与目标速度方向角之间的物理关系，通过采用滤波算法对目标速度方向角进行估计，间接求出了角速度的估计值，该方法的缺点是采用方向角等中间变量间接地估计角速度，将引入新的估计误差。文献[7]将角速度当作系统的未知参数，采用期望最大化(EM)算法将角速度与目标状态进行联合估计与辨识，该方法的缺点是使用的是全量测数据，算法运行时间长，不能满足空战中实时识别的要求。

在实际空战中，目标的角速度是未知且时变的，并且与目标状态相互耦合，因此在进行联合状态估计与角速度辨识时，为了获得角速度辨识的解析解，需要解除这种耦合关系。对在传统意义上的状态估计与参数辨识中，通常认为量测噪声之间是相互独立的，但是在实际跟踪过程中，雷达的采样频率过高[8]，使得相邻量测噪声之间的相关性往往不能够被忽略，因此，在对角速度辨识和目标状态估计的过程中必须考虑这种相关性。针对上述问题，本文基于EM算法框架提出一种带有色量测噪声的联合估计与辨识算法,并引入滑窗思想，以提高辨识的实时性：E-step利用当前估计的角速度以及经过带有色量测噪声的高阶容积卡尔曼平滑(HCKS)获得的后验平滑概率密度和联合概率密度，计算完备数据似然函数的条件期望；M-step通过使期望最大化进而更新获得下一次迭代的角速度估计量。

1 问题描述

蛇形机动可以分解成若干个具有不同转弯角速度的圆弧转弯机动，因此蛇形机动目标转弯角速度的辨识问题转化成圆弧转弯机动的角速度辨识问题。在有色量测噪声背景下，假设目标在二维平面中运动，转弯机动模型状态方程及量测方程为

(1)

f(xk,Ωk)=

vk=ψk,k-1vk-1+ζk-1

(2)

由式(1)可知，蛇形机动目标的转弯角速度Ωk这一未知参数耦合在状态转移矩阵之中，为了获取转弯角速度的解析解，需要解除转弯角速度与状态转移矩阵之间的非线性耦合关系，从而将状态方程等价转换成如式(3)形式：

xk+1=F(xk)θ+xk+wk

(3)

式中：

θ1=sin(ΩkT)/Ωk

θ2=(1-cos(ΩkT))/Ωk

θ3=1-cos(ΩkT)

θ4=sin(ΩkT)

由于在有色噪声条件下，传统的高斯近似滤波器和平滑器估计性能不佳，进而影响EM算法的辨识效果，所以为了实现目标状态的精确估计和转弯角速度的准确辨识，需要解除相邻量测之间相互耦合的关系，实现有色量测噪声的白化[10]。本文采取量测差分的方法实现有色量测噪声的白化，重新构造量测方程如下：

(4)

依据式(1)和式(2)可以得到具有一步状态延迟的量测方程：

(5)

2 基于HCKS-EM的联合估计与辨识算法

基于极大似然估计准则，本文提出了一种带有色量测噪声的联合估计与辨识算法对未知角速度进行辨识，具体的算法框架如图1所示。在第t次迭代时，通过将带有色噪声的量测序列经过HCKS，获得目标状态的平滑估计量，并采用EM算法进行参数θ的辨识，辨识后的参数θ用于第t+1次的状态估计，不断迭代直到满足设定的要求为止[12-13]。

图1 联合估计与辨识算法Fig.1 Joint estimation and identification algorithm

(6)

(7)

式中：

(8)

dxk-idxk-i-1

(9)

(10)

其中：l为滑动窗口的长度，其长度大小可依据实际需求自行设定。当l=k-1时，I1是一个常数，与待辨识量θ无关；当0≤l

pθ(xk-i|xk-i-1)～N(xk-i;F(xk-i-1)θ+xk-i-1,Q)

(11)

ψh(xk-i-1),R)

(12)

2.1 E-step

xk-i-1)TQ-1F(xk-i-1)·

(13)

式中：后验平滑概率密度和基于滑窗内数据的xk-i和xk-i-1的联合概率密度满足如下分布：

为了求解I2，首先定义如下2种积分函数：

(14)

(15)

然后分别将F(xi)和xi按xi的分量进行分解，可得

Γ(xk-i-1)和Δ(xk-i,xk-i-1)可分别表示为

(16)

(17)

(18)

(19)

则

(20)

(21)

2.2 M-step

(22)

当I2取得极大值时满足

(23)

则可求得参数θ的迭代表达式为

(24)

2.3 带有色量测噪声的HCKS算法的设计

众所周知，平滑是在滤波的基础之上进行的，因此首先设计一种带有色量测噪声的HCKF算法，算法如下：

算法1带有色量测噪声的高斯滤波算法框架。

步骤1状态预测。

(25)

步骤2状态修正。

(26)

(27)

带有色量测噪声的HCKS算法是一类前向后向平滑器，其主要包括2个部分：前向一步平滑和后向固定区间平滑，具体算法如下：

算法2带有色量测噪声的高斯平滑算法框架。

步骤1前向一步平滑。

(28)

步骤2后向固定区间平滑。

(29)

2.4 带有色量测噪声的HCKS算法的实现

本文中采用了Jia等[5]提出的5阶球径容积规则计算式(25)、式(27)和式(28)中的非线性高斯积分：

(30)

ξi=

(31)

式中：ei为n维单位向量，且其第i个元素为1。

(32)

容积点所对应的权重ωi为

(33)

可将本文算法总结如下：

算法3基于HCKS-EM的联合估计与辨识算法

步骤1E-step。已知第t次迭代后参数θ的估计值为θt，由式(28)和式(29)可以求得后验平滑概率密度和基于滑窗内数据的xk-i和xk-i-1的联合概率密度，进而可以求解得到I2的值。

递归：t=0，k←k+1，并将参数θ的初值设为上次迭代结束时的辨识值，进入下一轮迭代。

3 仿真分析

本文利用水平方向上的转弯机动非线性动态模型，仿真出一条蛇形机动轨迹。假设机动目标在1～80 s以Ω1=-2(°)/s作转弯运动，在k=81 s时转弯角速度突变为Ω2=2(°)/s并持续到300 s。设置转弯角速度初始值为Ω0=-1(°)/s，采样周期T=1 s，q1=0.1 m2/s3，过程噪声wk的协方差为Q。

Q=diag(q1M,q1M)

通过载机雷达可以获得目标与载机之间的相对距离r、方向角φ的信息，则可获得系统的非线性量测方程为

基于HCKS-EM联合估计与辨识算法和扩维法的初始状态以及协方差分别为

扩维法初始状态和协方差的分别设置为

100 mrad/s]

为了评估分析本文联合估计和辨识算法的性能，设置了如下3组仿真：

仿真1窗口长度l分别设置为5、10、15、20，最大迭代次数均为5次，各执行100次蒙特卡罗仿真。图2为不同窗口长度下角速度辨识结果，图3为不同窗口长度下角速度辨识均方根误差(RMSE)。从图2和图3可以看出，随着窗口长度的增大，该算法在第一阶段收敛于真实值的时刻就越早，精度越高，并且越稳定，但是当角速度发生突变时，对于突变的角速度反应的就越慢。并且从图4和表1可以看出，窗口长度越大，该算法估计的目标状态整体精度就越高，但是消耗的时间越长，这显然是时间与精度之间的“博弈”问题，当窗口大于10时，由窗口长度即量测数据带来的精度收益优势不太明显，相反运算时间的增加确实比较可观的，也就是说算法的执行效率降低了。并且在角速度发生突变时刻及其附近，窗口长度为10时的估计效果明显好于其他3个窗口。

图2 不同窗口长度下转弯角速度辨识结果Fig.2 Turn rate identification results with different window lengths

图3 不同窗口长度下转弯角速度辨识的均方根误差Fig.3 RMSE of turn rate identification with different window lengths

图4 不同窗口长度下位置和速度估计均方根误差Fig.4 RMSE of position and velocity estimation with different window lengths

滑动窗口长度时间/s50.0465100.1048150.1222200.1649

仿真2窗口长度设置为10，最大迭代次数r分别设置为3、5、10、15，各执行100次蒙特卡罗仿真。从图5和图6可以看出，随着迭代次数的增加，该算法在初始时刻收敛于真实值的时刻就越早，并且对于角速度突变反应的也比较灵敏。从图7目标状态4个分量的RMSE可以看出，迭代次数越大，该算法估计的目标状态整体精度就越高，但是根据表2中的运算时间可知花费的时间也会相应地增大。尤其是当最大迭代次数大于5时，最大迭代次数增加所带来的精度收益越发的不明显。

图5 不同最大迭代次数下转弯角速度辨识结果Fig.5 Turn rate identification result with different numbers of maximum iterations

图6 不同迭代次数下转弯角速度辨识的均方根误差Fig.6 RMSE of turn rate identification with different numbers of iterations

仿真3将本文联合估计与辨识算法窗口长度设置为10，最大迭代次数设置为5，将该算法与带有色的UKF算法和带有色的HCKF算法各执行100次蒙特卡罗仿真。从图8和图9可以看出，本文EM算法收敛于角速度真实值的时刻比带有色的UKF算法和带有色的HCKF算法都要早，并且稳定、精度高，但是当角速度发生突变时，对于突变的角速度反应的速度比较慢。从图10目标状态4个分量的RMSE可以看出，除了在角速度突变时刻及其邻近时刻之外，本文EM算法估计的精度都要明显高于带有色的UKF算法和带有色的HCKF算法估计的精度，带有色的HCKF算法比带有色的UKF算法估计的精度要高。

图7 不同迭代次数下位置和速度估计均方根误差Fig.7 RMSE of position and velocity estimation with different numbers of iterations

最大迭代次数时间/s30.049550.0827100.1599150.2577

图8 EM算法与扩维法之间的对比Fig.8 Comparison between EM algorithm and augmentation method

图9 转弯角速度辨识均方根误差Fig.9 RMSE of turn rate identification

图10 位置和速度估计均方根误差Fig.10 RMSE of position and velocity estimation

4 结 论

针对蛇形机动目标角速度辨识与目标状态联合估计的问题，基于EM算法框架提出一种带有色量测噪声的联合估计与辨识算法。主要结论有：

1) EM算法设计：在E-Step，通过将完备似然函数进行分解，从而将其转换成带参数θ的解析表达式，提高了辨识与估计的精度，降低了由状态扩维法扩维所带来的复杂性；在M-Step，通过极大化完备似然函数求得高精度的解析解。

2) 带有色量测噪声的HCKS算法的设计：考虑到在实际空战中，相邻噪声序列之间的相关性对于目标状态估计与转弯角速度辨识的影响，设计了带有色量测噪声的HCKS算法用于消除这种影响，提高了联合辨识与估计的精度。