文章快速检索  
  高级检索
基于HCKS-EM的战斗机蛇形机动目标跟踪算法
卢春光, 周中良, 刘宏强, 寇添, 杨远志     
空军工程大学 航空工程学院, 西安 710038
摘要: 针对有色量测噪声背景下战斗机蛇形机动模式转弯角速度辨识问题,考虑到目标状态与转弯角速度之间相互耦合的特性,基于期望最大化(EM)算法框架,提出了一种带有色量测噪声的联合估计与辨识算法。通过采用量测差分法实现了有色噪声白化,从而将有色量测噪声背景下的转弯角速度辨识问题转换成具有一步状态延迟的转弯角速度辨识问题。基于EM算法实现了战斗机蛇形机动目标状态与转弯角速度的联合估计与辨识:在E-step,通过利用有色量测噪声背景下的高阶容积卡尔曼平滑(HCKS)算法,获得了目标状态的后验估计;在M-step,通过极大化条件似然函数,进而获得转弯角速度的解析解。通过仿真验证了本文算法的目标状态估计与角速度辨识的精度均优越于传统的扩维法以及交互多模型法。而且又从窗口长度以及最大迭代次数2个方面评估分析了算法的性能,仿真结果表明,窗口长度以及最大迭代次数越大,精度越高。
关键词: 蛇形机动     有色量测噪声     量测差分法     期望最大化(EM)     联合估计与辨识    
Fighter zigzag maneuver target tracking algorithm using HCKS-EM
LU Chunguang, ZHOU Zhongliang, LIU Hongqiang, KOU Tian, YANG Yuanzhi     
Aeronautics Engineering College, Air Force Engineering University, Xi'an 710038, China
Received: 2018-01-18; Accepted: 2018-04-13; Published online: 2018-05-18 10:48
Foundation item: National Natural Science Foundation of China (61472443)
Corresponding author. ZHOU Zhongliang.E-mail:zzl_panda@163.com
Abstract: Motivated by identifying the turn rate of fighter zigzag maneuver under the background of colored measurement noise, the joint estimation and identification algorithm with colored measurement noise is proposed based on expectation maximization (EM) algorithm by considering the characteristics of the coupling between the target state and the turn rate. The colored noise whitening is realized by using the measurement difference scheme, and thus, the turn rate identification problem with colored measurement noise is transformed into the turn rate identification problem with one-step delayed state. The joint estimation and identification of both fighter zigzag maneuver target states and turn rate are achieved by EM algorithm:in the E-step, the target state posteriori estimation is achieved accurately using the high-degree cubature Kalman smoothers (HCKS) algorithm with colored measurement noise; in the M-step, the analytical identification result of turn rate is obtained by maximizing the conditional likelihood function. It is verified in the final simulation that the proposed algorithm performs better in terms of target state estimation and turn rate identification accuracy than the traditional augmentation method and interacting multi-model algorithm. Furthermore, the performance of the proposed algorithm is evaluated and analyzed from two aspects of window length and maximum number of iterations. The simulation results show that the larger the window length and the maximum number of iterations are, the higher the precision is.
Keywords: zigzag maneuver     colored measurement noise     measurement difference scheme     expectation maximization (EM)     joint estimation and identification    

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

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

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

1 问题描述

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

(1)

式中:T为采样间隔; Ωk为待辨识的转弯角速度; , xkRn为状态变量; zkRm为系统量测向量; h(·)为非线性量测函数; wk为零均值高斯白噪声, 协方差为Q; vk为高斯有色噪声, 满足一阶自回归(AR)方程[9]:

(2)

其中:ψk, k-1为已知的常数矩阵; ζkRm为零均值高斯白噪声, 协方差为Rwkζk互不相关。初始状态x0wkζk无关, 并且满足x0~N(x0; , P0)。

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

(3)

式中:

则通过辨识出参数θ, 进而根据可以求得蛇形机动目标的转弯角速度Ωk

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

(4)

式中:zk+1*为重新构造的白色量测。

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

(5)

定理1[11]  如果pθ(xs|Z1*k)服从高斯分布, pθ(xs|Z1k)也服从高斯分布, 则在最小均方差意义下满足:

也就是说在最小均方差误差意义下, 基于白化后的量测Z1*k的状态估计与基于有色量测Z1k的状态估计是等价的。

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

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

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

在系统状态X1k={x1, x2, …, xk}缺失, 量测序列Z1*k={z1*, z2*, …, zk*}已知的前提下, k为当前时刻量测数据样本数, 则可知完备数据似然函数为[14-15]

(6)

给定第t次迭代后得到的参数θ的估计值, 则完备数据似然函数的关于概率密度函数 (xk-l-1|Z1*k)的期望为

(7)

式中:

(8)
(9)
(10)

其中:l为滑动窗口的长度, 其长度大小可依据实际需求自行设定。当l=k-1时, I1是一个常数, 与待辨识量θ无关; 当0≤l < k-1时, I1θ相关性较弱, 而且为了满足实时性的要求, 辨识算法应尽量关注当前时刻的数据或者邻近时刻的数据, 所以本文仅仅考虑I2I3对条件期望函数的影响。

假设先验概率密度函数pθ(xk-i|xk-i-1)和pθ(zk-i*|xk-i, xk-i-1)均服从高斯分布, 并且满足

(11)
(12)

则将式(11)和式(12)分别代入I2I3当中, 可知I3与未知参数θ无关, 所以, 在辨识参数θ时, 仅仅考虑I2对条件期望函数的影响, 即条件期望函数的大小仅取决于I2的大小。

EM算法通过不断迭代求得的极大似然估计, 每次迭代包含2个部分:E-step和M-step。

2.1 E-step

E-step的主要任务是根据先验概率密度、后验平滑概率密度以及xk-ixk-i-1的联合概率密度进行确定性采样计算条件期望函数

通过上述的分析可知, 条件期望函数大小仅取决于I2的大小, 故在此仅仅考虑I2的计算。将式(11)代入到式(9)中, 可求解得到I2关于θ的解析表达式为

(13)

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

其中:平滑量Pk-i-1sPk-is可通过带有色量测噪声的高阶容积平滑获得[10]; Ak-i-1为固定区间平滑增益。

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

(14)
(15)

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

式中:Λj为常数矩阵; xij表示xi的第j个分量; Ψj表示n×n单位矩阵的第j列; nxi的维数。

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

(16)
(17)

根据状态平滑的误差协方差Pk-i-1s(i=0, 1, …, l)相关定义可进一步推导出:

(18)
(19)

(20)
(21)
2.2 M-step

M-step的主要任务是解决I2的极大化问题, 即求解使I2满足极大值时所对应的值, 用于EM算法的下一次迭代更新, 表达式为

(22)

I2取得极大值时满足

(23)

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

(24)

通过不断迭代最终求得k时刻参数θ的估计值, 通过反解参数θ与参数Ωk之间的函数关系, 可以求得k时刻目标的转弯角速度的辨识结果。

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

由于vk为高斯有色噪声, 通过白化有色噪声vk而获得的量测zk+1*xk+1xk相关, 为方便后续的计算, 在此将目标状态进行扩维, 表示为

在最小均方误差意义下, 扩维状态的一步预测及误差协方差Pk+1|ka

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

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

步骤1  状态预测。

(25)

步骤2  状态修正。

(26)
(27)

式中:hk+1*(xk+1|ka)=hk+1(xk+1)-ψk+1, khk(xk)。

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

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

步骤1  前向一步平滑。

(28)

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

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

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

(30)

式中:Px=SxSxT; ξiHCKF=ξiSx+x

(31)

式中:ein维单位向量, 且其第i个元素为1。

(32)

容积点所对应的权重ωi

(33)

可将本文算法总结如下:

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

初始化:给定量测集合Z1*k, 目标初始状态x0和协方差P0, 角速度初始值为Ω0, 最大迭代次数为5。

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

步骤2  M-step。通过式(23)和式(24)可求解得到参数θ的估计值用于下一次迭代。

步骤3  迭代终止。给定一个阈值ε, 如果或者t≤5则tt+1, 否则, 立即终止迭代。

递归:t=0, kk+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

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

式中:vk为高斯有色噪声, 与其对应的ζk是零均值高斯白噪声, 协方差为Rk=[σr2 σφ2], σr=10 m, , 量测噪声相关系数为ψk+1, k=diag[0.2 0.1]。

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

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

为了评估分析本文联合估计和辨识算法的性能, 设置了如下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
表 1 k=50 s时的角速度和状态所耗费的时间 Table 1 Turn rate and state calculation time when k=50 s
滑动窗口长度 时间/s
5 0.046 5
10 0.104 8
15 0.122 2
20 0.164 9

仿真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
图 7 不同迭代次数下位置和速度估计均方根误差 Fig. 7 RMSE of position and velocity estimation with different numbers of iterations
表 2 k=120 s时的角速度和状态所耗费的时间 Table 2 Turn rate and state calculation time when k=120 s
最大迭代次数 时间/s
3 0.049 5
5 0.082 7
10 0.159 9
15 0.257 7

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

图 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算法用于消除这种影响, 提高了联合辨识与估计的精度。

参考文献
[1] 宋华, 章新华, 许林周. 基于离散隐马尔科夫模型的空中目标战术机动识别[J]. 仪器仪表学报, 2007, 28 (4): 588–592.
SONG H, ZHANG X H, XU L Z. Aerial combat maneuver identification based on discrete hidden Markov model[J]. Chinese Journal of Scientific Instrument, 2007, 28 (4): 588–592. (in Chinese)
[2] 钟友武, 柳嘉润, 申功璋. 自主近距空战中敌机的战术动作识别方法[J]. 北京航空航天大学学报, 2007, 33 (9): 1056–1059.
ZHONG Y W, LIU J R, SHEN G Z. Recognition method for tactica maneuver of target in autonomous close-in air combat[J]. Journal of Beijing University of Aeronautics and Astronautics, 2007, 33 (9): 1056–1059. DOI:10.3969/j.issn.1001-5965.2007.09.013 (in Chinese)
[3] ROTH M, HENDEBY G, GUSTAFSSON F.EKF/UKF maneuvering target tracking using coordinatd turn models with polar/Cartesian velocity[C]//17th International Conference on Information Fusion (FUSION).Piscataway, NJ: IEEE Press, 2014: 1-8.
[4] ARASARATNAM I, HAYKIN S. Cubature Kalman smoother[J]. Automatica, 2011, 47 (10): 2245–2250. DOI:10.1016/j.automatica.2011.08.005
[5] JIA B, XIN M, CHENG Y. High-degree cubature Kalman filter[J]. Automatica, 2013, 49 (2): 510–518. DOI:10.1016/j.automatica.2012.11.014
[6] 黄伟平, 徐毓, 王杰. 基于改进"当前"统计模型的转弯机动跟踪算法[J]. 控制与决策, 2011, 26 (9): 1412–1416.
HUANG W P, XU Y, WANG J. Algorithm based on modified current statistic mode fo turn maneuver[J]. Control and Decision, 2011, 26 (9): 1412–1416. (in Chinese)
[7] SONG B, WANG X X, LIANG Y, et al.Analytical identification of system parameter nonlinearly coupled in dynamic transition matrix[C]//American Control Conference(ACC).2016: 1832-1837.
[8] 王小旭, 梁彦, 潘泉, 等. 带有色量测噪声的非线性系统Unscented卡尔曼滤波器[J]. 自动化学报, 2012, 38 (6): 986–998.
WANG X X, LIANG Y, PAN Q, et al. Unscented Kalman filter for nonlinear systems with colored measurement noise[J]. Acta Automatica Sinica, 2012, 38 (6): 986–998. (in Chinese)
[9] WANG X X, LIANG Y, PAN Q, et al. Nonlinear Gaussian smoothers with colored measurement noise[J]. IEEE Trasactions on Automatic Control, 2015, 60 (3): 870–876. DOI:10.1109/TAC.2014.2337991
[10] 黄玉龙, 张勇刚, 李宁, 等. 一种带有色量测噪声的非线性系统辨识方法[J]. 自动化学报, 2015, 41 (11): 1877–1892.
HUANG Y L, ZHANG Y G, LI N, et al. An identification method for nonlinear systems with colored measurement noise[J]. Acta Automatica Sinica, 2015, 41 (11): 1877–1892. (in Chinese)
[11] WANG X X, LIANG Y, PAN Q, et al. Nonlinear Gaussian smoothers with colored measurements noise[J]. IEEE Transactions on Automatic Control, 2015, 60 (3): 870–876. DOI:10.1109/TAC.2014.2337991
[12] NOBUHIRO Y. Parameter estimation of aircraft dynamics via unscented smoother with expectation-maximization algorithm[J]. Journal of Guidance, Control, and Dynamics, 2011, 34 (2): 426–436. DOI:10.2514/1.51696
[13] LAN H, LIANG Y, YANG F, et al. Joint estimation and identification for stochastic systems with unknown inputs[J]. IET Control Theory & Appications, 2013, 7 (10): 1377–1386.
[14] SCHON T, WILLS A, NINNESS B. System identification of nonlinear state-space models[J]. Automatica, 2011, 47 (1): 39–49. DOI:10.1016/j.automatica.2010.10.013
[15] HUANG Y, ZHANG Y G, LI N, et al. Latency probability estimation of non-linear systems with one-step randomly delayed measurements[J]. IET Control Theory & Appications, 2016, 10 (7): 843–852.
http://dx.doi.org/10.13700/j.bh.1001-5965.2018.0047
北京航空航天大学主办。
0

文章信息

卢春光, 周中良, 刘宏强, 寇添, 杨远志
LU Chunguang, ZHOU Zhongliang, LIU Hongqiang, KOU Tian, YANG Yuanzhi
基于HCKS-EM的战斗机蛇形机动目标跟踪算法
Fighter zigzag maneuver target tracking algorithm using HCKS-EM
北京航空航天大学学报, 2018, 44(9): 2004-2012
Journal of Beijing University of Aeronautics and Astronsutics, 2018, 44(9): 2004-2012
http://dx.doi.org/10.13700/j.bh.1001-5965.2018.0047

文章历史

收稿日期: 2018-01-18
录用日期: 2018-04-13
网络出版时间: 2018-05-18 10:48

相关文章

工作空间