基于有向图切换IMM-CKF高速滑翔目标跟踪算法

崔乃刚，蔡李根，荣思远

(1. 哈尔滨工业大学 航天学院，黑龙江 哈尔滨 150001； 2. 北京航天长征飞行器研究所，北京 100076)

0 引 言

近年来，探测、跟踪及拦截武器等各种反导技术得到了长足的发展，传统的弹道导弹面临突防效能的巨大挑战，战略威慑力逐步减弱。与此同时，飞行速度快和突防能力强的临近空间高速滑翔飞行器开始成为国内外学术和工程研究的热点。美国已经开展过多次飞行试验，积累了丰富的工程经验。在此背景下，临近空间高速滑翔飞行器的轨迹跟踪方法也开始成为热点，该方法不仅能够为拦截系统提供准确的目标运动信息，还可以为目标的轨迹预报提供依据，具有极大的理论研究价值。

临近空间高速滑翔飞行器的飞行轨迹一般包括助推段和滑翔段等多个阶段，其中滑翔段飞行模式可分为跳跃滑翔和平滑滑翔两种，具有机动能力强、机动样式多、机动时机不确定等特点[1-2]，是临近空间高速滑翔飞行器轨迹跟踪的关键阶段。李昌玺等[3]在“当前”统计模型基础上，引入一种利用位置估计值与加速度之间的函数关系自适应调整加速度方差的方法，对“当前”统计模型进行了改进，弥补了“当前”统计模型机动加速度的极值参数难以确定的缺陷；王国宏等[4]针对临近空间高速滑翔飞行器典型的跳跃滑翔机动模式，设计了一种SW模型，仿真证明SW模型对于跳跃滑翔目标的跟踪具有明显优势；张裕禄等[5]将气动加速度中的阻力和爬升力加速度建模为正弦自相关模型，转弯加速度建模为一阶马尔科夫模型，采用分离滤波模型对目标状态和气动加速度分别进行估计，仿真结果表明，该模型的跟踪精度比气动参数增广的经典动力学模型的跟踪精度更高。

由于机动形式复杂多样，对临近空间高速滑翔飞行器的运动建模十分困难，采用单个模型描述其运动会出现跟踪精度不理想甚至出现发散问题，因此有必要应用多模型算法对这类目标进行精确跟踪。苗伟等[6]提出了基于修正转弯模型的交互多模型跟踪算法，跟踪性能明显优于单模型算法；秦雷等[7]利用自适应网格交互式多模型(adaptive grid interacting multiple model，AGIMM)算法对进行转弯运动的高马赫数目标进行跟踪，具有比传统交互式多模型算法更高的精度；肖楚晗等[8]在AGIMM算法的基础上，提出了一种自适应变结构交互式多模型(adaptive variable structure interacting multiple model，AVSIMM)算法跟踪高马赫数再入滑翔目标，通过模型集参数与算法结构的双重自适应调整进一步提高了目标的跟踪精度。

在临近空间高速滑翔飞行器轨迹跟踪系统中，状态模型和测量模型均具有较强的非线性，需要利用非线性滤波方法。非线性滤波方法是探测传感器测量信息最优融合的关键技术，影响着高速滑翔飞行器轨迹跟踪精度。工程上应用最为广泛的扩展卡尔曼滤波方法由于其计算量小、技术成熟的显著优点也常常被应用于临近空间高速滑翔飞行器轨迹跟踪[9-10]，但是扩展卡尔曼滤波存在较大的截断误差，会导致滤波精度下降甚至发散。因此，学者们利用确定性的采样点和权值近似随机变量的高斯分布，提出了无迹卡尔曼滤波[11]、容积卡尔曼滤波[12]等多种非线性近似、精度更高的滤波算法，在临近空间高速滑翔飞行器跟踪领域得到了广泛应用。

综上所述，临近空间高速滑翔飞行器的轨迹跟踪存在机动目标机动形式复杂多样性与突变性、模型不确定性造成的运动建模困难、模型非线性造成的非线性估计误差等难点。本文首先设计了基于有向图切换的多模型结构，在交互式多模型算法基础上，引入有向图切换方法得到了基于有向图切换交互式多模型结构；然后介绍了容积卡尔曼滤波算法，并基于以SW模型为主的目标运动模型集和地基雷达测量模型完成了临近空间目标高精度跟踪方法的设计；最后，采用数学仿真手段对所提出的算法进行了验证，并进行相关分析。

1 基于有向图切换的多模型结构

滑翔段是临近空间高速滑翔飞行器跟踪过程中最关键的阶段，具有机动样式突变性和机动能力不确定性的特点。跳跃滑翔作为高速滑翔飞行器的显著特点，其机动周期会随时间呈近似指数衰减，即周期机动角速度会随时间增大。考虑到目标在纵向跳跃滑翔机动情况下存在的周期衰减现象，需要采用多个不同周期机动角速度对应的SW模型集应对这一情况。所以，需要利用多模型方法设计滑翔段高速滑翔飞行器的跟踪算法。

目前应用最为广泛的多模型算法是交互式多模型(interacting multiple model，IMM)算法，模型集的设计很大程度上决定了IMM算法的性能。为了提高IMM算法的性能，需要尽可能多地覆盖目标所有可能的运动形式，那就必须增加模型集的模型数。但是模型数的增加会使计算量成倍增加，同时还会导致模型之间出现竞争，从而降低算法的性能。因此，需要在模型数和计算量之间寻求一个较好的平衡，为此本文引入有向图切换方法，设计了基于有向图切换的多模型结构。

1.1 交互式多模型结构

IMM算法用一个覆盖目标可能运动形式的固定模型集来描述目标的运动，每个模型分别设计相应子滤波器并行运行，各子滤波器间存在数据交互，在每一时刻基于最大似然方法计算出各个模型的后验概率之后，通过对各模型的状态估计加权求和得到最终的目标估计结果。IMM算法的主要步骤如下：

1) 数据输入交互

计算k-1时刻混合概率

i,j=1,2,...,r

(1)

式中：pij是模型i到模型j的转移概率；r为模型集模型数目；归一化常数μj(k|k-1)为

(2)

计算输入混合状态估计值和协方差阵

(3)

(4)

2) 状态滤波

3) 模型概率更新

似然函数Λj(k)根据极大似然函数N(·)定义，如式(5)所示。

Sj[k;P0j(k-1)]}

(5)

其具体形式为

(6)

模型概率更新

(7)

4) 滤波综合输出

(8)

1.2 有向图切换交互式多模型结构

变结构多模型算法使用实时可变的模型集，模型集包含不同模型的组合，在每一个滤波时刻根据测量输入选取最适合的模型集，基于IMM算法的框架进行状态估计。变结构多模型有一个完备的模型总集，包含各种可能的运动模型，实时的模型集通过实时的量测从完备模型集中选取固定数量的模型集。为了更好地表示变结构多模型算法的模型切换，用M表示模型，M0表示模型总集，C表示模型集，C0表示模型集总集。

递推的变结构多模型算法一般是在IMM算法的框架上，加入自适应调整模型集的步骤。本文设计基于有向图切换的多模型结构，以模型集为单位进行切换，预先设定一系列同样模型数的模型集C0={C1,C2,C3,…,CN}。对于比较对称的拓扑结构和具有较明确物理意义的有向图，有向图切换法比较实用。有向图切换法将模型总集分为若干个模型集，一个模型可以属于不同的模型集，但是不存在完全一样的模型集。一般某个模型集的边缘模型是其相交的模型集的核心模型，所有模型集的总集C0可以完成对模型总集M0的完全覆盖。

如图1所示，有向图切换法根据实时模型集中边缘模型的模型概率大小来判断是否发生转移，从而进行模型集间的切换，即如果当前模型集中某一边缘模型的模型概率最大且大于设定阈值，则下一时刻模型集切换为以该边缘模型为核心模型的模型集。

图1 有向图切换示意图Fig. 1 Schematic diagram of directed graph switching

2 基于容积卡尔曼滤波的非线性滤波方法

容积卡尔曼滤波方法是一种以确定性采样逼近非线性分布的方法，经严格数学推导可知容积卡尔曼滤波方法能达到三阶近似精度，比扩展卡尔曼滤波具有更高的非线性近似精度，比无迹卡尔曼滤波方法具有更强的高维稳定性。因此，本文选择容积卡尔曼滤波方法，将其应用于临近空间高速滑翔飞行器的跟踪。

容积卡尔曼滤波算法首先计算加权函数为标准正态分布密度的积分的基本容积点和对应的权值，如式(9)所示。

(9)

式中：n是系统状态维数；m是容积点数，满足m=2n；ξj为第j个基本容积点，其对应的权值均为1/m；[1]j为第j个基本容积点对应的n维列向量，[1]具体形式为

(10)

容积卡尔曼滤波的具体步骤如下：

1) 时间更新

Step 1.计算容积点：

(11)

Step 2. 计算通过非线性状态方程传播的容积点：

(12)

式中：f(·)为状态方程函数。

Step 3. 计算状态和误差协方差矩阵一步预测：

(13)

式中：Qk-1为过程噪声协方差。

2) 量测更新

Step 1.计算用于量测更新的容积点：

(14)

Step 2. 计算通过非线性量测方程传播的容积点：

(15)

式中：h(·)为量测方程函数。

(16)

式中：Rk为测量噪声矩阵。

Step4. 计算增益矩阵：

(17)

(18)

3 临近空间目标高精度跟踪方法

本章首先建立临近空间目标运动模型集合地基雷达测量模型，结合前文提出的基于有向图切换的多模型结构和容积卡尔曼滤波方法得到临近空间目标高精度跟踪方法。

3.1 临近空间目标运动模型集

临近空间高速滑翔飞行器纵向的跳跃滑翔是一种类似“打水漂”的跳跃飞行，是临近空间高速滑翔飞行器最鲜明的特点，其轨迹跳跃幅度呈衰减形式，类似于振荡运动，可以采用文献[4]所述的SW模型进行目标运动建模。SW模型将目标加速度a(t)的自相关函数建模为正弦波自相关时间函数Rτ，如式(19)所示。

(19)

选取状态量为目标的位置、速度、加速度与加加速度，即

(20)

SW模型离散化的状态方程为

Xk=Φk|k-1Xk-1+Wk-1

(21)

式中：Xk-1表示k-1时刻的状态量；W(k)为离散时间过程噪声序列，过程噪声协方差Qk-1的取值详见文献[4]。Φk|k-1为状态转移矩阵，满足

(22)

其中，T为计算周期。

周期机动的角速度ω是SW模型的主要模型参数，考虑到临近空间高速滑翔飞行器在纵向跳跃滑翔机动情况下存在的周期衰减现象，需要设计多个不同周期机动角速度的SW模型构成模型集，以应对这一情况。首先确定跳跃滑翔机动情况下周期机动角速度ω变化的可能范围[ωmin,ωmax]，并在[ωmin,ωmax]区间里设计一系列不同周期机动角速度对应的SW模型构成模型集，以覆盖跳跃滑翔情况下周期机动角速度ω的变化。

假设在[ωmin,ωmax]区间里设计了5个周期机动角速度对应的SW模型，且对应的周期机动角速度ω取值分别为ωa、ωb、ωc、ωd和ωe，ωa<ωb<ωc<ωd<ωe。设实时模型集数量为3，则可以设计3个可能的模型集，分别是C1⊃{ωa,ωb,ωc}、C2⊃{ωb,ωc,ωd}和C3⊃{ωc,ωd,ωe}，则模型切换逻辑如下：

1) 当前模型集为C1。如果ωc对应的SW模型的实时模型概率最大且大于阈值t1，下一时刻模型集切换为C2；否则，下一时刻模型集保持不变；

2) 当前模型集为C2。如果ωb对应的SW模型的实时模型概率最大且大于阈值t1，下一时刻模型集切换为C1；如果ωd对应的SW模型的实时模型概率最大且大于阈值t1，下一时刻模型集切换为C3；否则，下一时刻模型集保持不变；

3) 当前模型集为C3。如果ωc对应的SW模型的实时模型概率最大且大于阈值t1，下一时刻模型集切换为C2；否则，下一时刻模型集保持不变。

3.2 地基雷达测量模型

考虑到临近空间高速滑翔飞行器在滑翔段为无动力飞行，无明显红外等辐射特征，且具有较大的机动能力，机动范围较大，一般利用地基雷达通过主动发射电磁波对目标进行探测。地基雷达具有功率大、作用距离远、探测精度高等优势。如图 2所示，利用地基雷达对处于滑翔段飞行的临近空间高速滑翔飞行器进行探测。

图2 地基雷达测量示意图Fig. 2 Schematic diagram of ground-based radar measurement

地基雷达的测量量包括雷达站距目标的斜距ρ、高低角γ和方位角η，测量量可定义为

Z=[γ,η,ρ]T

(23)

相应的测量方程为

(24)

式中：xe,ye和ze为临近空间高速滑翔飞行器在雷达体坐标系中的位置分量；vγ,vη和vρ分别为高低角、方位角和斜距的测量误差。

3.3 算法流程

基于临近空间目标运动模型集、地基雷达测量方程、有向图切换方法、容积卡尔曼滤波方法，可以得到临近空间目标高精度跟踪方法基本流程如下：

1) 模型集自适应：利用有向图切换方法从初始时刻到k-1时刻的模型集序列Ck-1和k时刻的观测Zk来确定k时刻的模型集C0(k)。

2) 数据输入交互：基于C0(k)进行类似交互式多模型算法的数据输入交互，如式(1)～(4)所示，模型的初始状态估计和模型概率均由k-1时刻的模型集C0(k-1)转移或者加权得到。

4) 模型概率更新：结合k时刻的观测Zk与每个子滤波器估计的量测值，利用极大似然函数计算模型集C0(k)中各模型的后验概率，即式(5)～(7)，进行下一时刻模型概率的更新。

5) 综合输出：利用更新所得的模型概率对各子滤波器状态滤波结果进行加权输出，如式(8)。

4 仿真分析

4.1 仿真场景

为了验证本文所提临近空间目标高精度跟踪方法的有效性，采用数学仿真方法进行验证，具体仿真场景如下。

临近空间高速滑翔飞行器的初始位置为[160°E,28°N,30 km]，做跳跃滑翔机动，初始速度v=5 600 m/s，初始航向角ψ0=80°，航迹倾角γ0=0°，临近空间高速滑翔飞行器的初始跳跃周期为ΔT0=209.4s，初始周期机动角速度ω0=0.03 rad/s，跳跃周期呈指数衰减规律。获得仿真场景如图 3所示，目标的航迹倾角变化曲线如图 4所示，目标的速度变化曲线如图 5所示，目标的跳跃周期和周期机动角速度的变化曲线如图 6所示。地基雷达的位置为[122°E,33°N,0 m]，采样频率为5Hz，测角精度为1 mrad(1σ)，测距精度为10 m(1σ)。

图3 仿真场景Fig.3 Simulation scene

图4 航迹倾角变化曲线Fig.4 Time history of track inclined angle

图5 目标速度变化曲线Fig.5 Time history of target velocity

图6 跳跃周期和周期机动角速度变化曲线Fig.6 Time histories of jump period and maneuver angular velocity

4.2 仿真结果

在上述仿真场景下，采用容积卡尔曼滤波方法分别结合文献[4]所提的SW模型单模型、文献[6]所提的交互式多模型和本文的有向图切换交互式多模型(digraph switching IMM,DSIMM)算法进行临近空间高速滑翔飞行器的轨迹跟踪。设置单个SW模型的周期机动角速度ω=0.03 rad/s；设置交互式多模型算法的模型集由3个SW模型构成，其周期机动角速度分别为{ω=0.04 rad/s，ω=0.05 rad/s，ω=0.06 rad/s}；有向图切换交互式多模型算法的完备模型集有5个SW模型构成，其周期机动角速度分别为ωa=0.03 rad/s、ωb=0.04 rad/s、ωc=0.05 rad/s、ωd=0.06 rad/s和ωe=0.07 rad/s；可能模型集分别为C1⊃{ωa,ωb,ωc}、C2⊃{ωb,ωc,ωd}和C3⊃{ωc,ωd,ωe}；有向图切换阈值t1=0.75。进行100次相互独立的Monte Carlo仿真，仿真结果如图7～8所示。

图7 位置RMSE对比图Fig.7 Position RMSE comparison

图8 速度RMSE对比图Fig.8 Velocity RMSE comparison

图9 周期机动角速度参数变化曲线Fig.9 Time histories of maneuver angular velocity parameter

定义平均均方根误差(root mean squared error,RMSE)为

将3种算法对目标跟踪结果的平均RMSE汇总至表1。

表1 各种算法的平均RMSETab.1 Average RMSE of different algorithms

由以上结果可知，文献[4]所提单个SW模型的跟踪算法在初始周期机动角速度ω接近0.03 rad/s的阶段内跟踪效果较好，随着周期机动角速度ω的不断增大，由于模型参数越来越不符合实际，跟踪效果越来越差；基于文献[6]所提的交互式多模型的跟踪算法相比于单个SW模型的跟踪算法而言，跟踪性能有明显提升，但是由于模型集固定不变，在周期机动角速度ω超出模型集覆盖范围时，跟踪性能明显下降；本文提出的有向图切换交互式多模型的跟踪算法跟踪性能最高，由图 9可以看到，有向图切换交互式多模型算法实时的模型集可以随周期机动角速度ω的变化进行实时调整，跟踪性能显著提升，从而验证了有向图切换交互式多模型跟踪算法的优越性和模型集切换的有效性。

5 结束语

本文研究了临近空间高速滑翔飞行器做典型的跳跃滑翔机动时进行轨迹跟踪的方法，将容积卡尔曼滤波方法与有向图切换交互式多模型算法相结合，有效地进行了模型集的切换。与交互式多模型算法相比，在周期机动角速度变化范围较广时，有向图切换交互式多模型算法具有一定的模型集自适应性。