由于很难获取目标运动的先验信息，仅采用单模型的跟踪算法在不同机动情形下的跟踪效果往往较差。为了解决这个问题，多模型(Multiple Model, MM)算法[1]得到了深入研究并被广泛应用于目标跟踪[2-3]和轨迹预测[4]等领域。MM算法的核心在于假设模型集包含目标所有可能运动模式的基础上，利用可能模型的概率加权匹配目标运动。MM算法的状态观测一般是模式依赖的，而且量测序列嵌入模式信息，即系统模式序列是Morkov过程。根据模型集设计方法，可以将MM算法分为固定结构和可变结构MM算法[5]。其中，固定结构MM算法的代表是n阶广义伪贝叶斯(n-order Generalized Pseudo-Bayesian, GPBn)算法[6]和交互多模型(Interacting Multiple Model, IMM)算法[7]。IMM算法通过应用更好的假设管理技术，具有GPB2的性能和GPB1计算上的优势, 逐渐成为研究热点。针对固定结构MM算法高度依赖模型集和模型竞争导致的性能下降等问题[8]，可变结构MM算法[9-10]通过模型集自适应算法的开发具有较好的跟踪性能和费效比。
IMM算法是在GPB算法的基础上提出的一种具有系统模式Morkov切换的算法。IMM算法进行目标状态估计的计算时，考虑每个模型滤波器都有可能成为当前有效的系统模型滤波器，每个滤波器的初始条件都是基于前一时刻各条件模型滤波器结果的合成，即模型条件重初始化。由于合成过程中存在非等维状态向量的混合，例如当匀速(Constant Velocity, CV)模型和匀加速(Constant Acceleration, CA)模型进行混合估计时，加速度分量即为“额外”分量，而标准IMM算法并未给出关于此问题的描述。IMM混合估计需要匹配状态向量的前两阶矩(均值和协方差)，因此低维状态向量需要扩维以匹配高维状态向量，这也是本文讨论的主题。然而，目前此问题只有少量的研究[11-14]。一种简单的方法[11]是对低维状态采用零均值和方差进行扩维，但会导致有偏估计。文献[12, 14]提出了一种无偏方法，即利用高维状态中的均值和协方差对低维状态进行扩维。该方法通过弹道导弹数据进行了检验并展现了良好的性能。文献[13]认为被扩维的分布应该尽可能反映系统潜在情况(模式切换)，采用均匀分布描述“额外”分量，因而提高了模式切换时的跟踪精度。
本文所提方法基于以下思想：当前时刻的最优混合方法应根据时变的目标状态进行动态调整。本文方法在分析现有方法的基础上，将“切换”状态的概念引入目标状态，并利用模型概率和滤波新息进行联合判决，从而选择最优的混合方法，提高跟踪精度。
1 问题描述对于一个跳变线性系统，其数学描述为
|
(1) |
|
(2) |
式中:xk和zk分别为k时刻的状态向量和量测；Fj、Hj和Gj分别为模型j的状态转移矩阵、量测矩阵和噪声分布矩阵；wkj~
|
(3) |
式中：Φ表示模型集；mkj表示k时刻模型j与系统真实模式匹配。为描述简便起见，假设模型集M只包含2个非等维状态的模型。事实上，本文算法容易推广到包含任意多个模型的模型集。
假设k时刻状态估计的后验概率是高斯分布的，即
|
(4) |
式中:
|
(5) |
|
(6) |
式中：μki, j为在给定k+1时刻模式为j时，k时刻模型为i的概率；μki为k时刻模型i概率；cj为归一化因子。混合估计的状态估计及协方差为
|
(7) |
|
(8) |
假设模型1状态向量的维度高于模型2，即nx1>nx2，且模型2的状态分量包含于模型1，即
|
(9) |
|
(10) |
式中:ck为模型1和2的公共分量；
|
(11) |
|
(12) |
然而，对于模型1的重初始化，低维的状态向量必须按照某种方法进行扩维，这也正是本文所考虑的问题。
IMM算法是递推的，包括模型条件重初始化、模型条件滤波、模型概率更新和估计融合。其中，式(5)~式(8)给出了模型条件重初始化。
模型条件滤波是在给定重初始化的状态和协方差条件后进行状态估计更新，包含以下步骤：
1) 状态预测及协方差，计算式分别为
|
(13) |
|
(14) |
2) 新息(滤波残差)及其协方差，计算式分别为
|
(15) |
|
(16) |
3) 滤波更新，包括滤波增益、状态估计更新及其协方差，计算式分别为
|
(17) |
|
(18) |
|
(19) |
模型概率更新，计算式分别为
|
(20) |
|
(21) |
式中:Λk+1i为与mk+1i匹配的似然函数，服从高斯分布。
估计融合给出k+1时刻的总体估计及其协方差，计算式分别为
|
(22) |
|
(23) |
如第1节所述，本文所研究的焦点是如何将“额外”分量体现到扩维后的状态向量上。假设模型2扩维后的状态xk2, a分布为
|
(24) |
式中：q(ek2)为模型2中“额外”分量ek的分布。显然，不同混合方法的区别在于q(ek2)的选择。本文总结了关于此问题以往的研究[11-14]，下面介绍现有的3种方法，分别记为A1，A2和A3。
1) A1采取了一种简单的策略，认为由于ek是模型1的分量并不存在于模型2，故可以用零均值和协方差进行扩维，则相应的分布q(ek2)和扩维后的状态xk2, a为
|
(25) |
式中：δ为冲激函数。
|
(26) |
2) 文献[12, 14]提出了一种无偏方法A2，认为A1会导致模型1状态估计有偏，而ek2利用模型1的“额外”分量ek信息进行扩维是无偏的，相应的分布q(ek2)和扩维后的状态xk2, a为
|
(27) |
|
(28) |
3) A3考虑到系统模式切换的可能性，认为ek2在其取值区间内服从均匀分布，则相应的分布q(ek2)和扩维后的状态xk2, a为
|
(29) |
|
(30) |
式中:ne=nx1-nx2为ek的维度；参数el取值区间为[al, bl]；Cov(ek2)表示ek2的协方差。
事实上，上述3种方法是基于不同假设下设计的。A1假设系统模式为mk2，因此将“额外”分量ek置零；A2假设系统模式为mk1，因此ek2直接利用ek1信息作为其分布；A3假设系统模式可能切换，因此认为ek2为均匀分布。显然，在假设成立的条件下，这些算法能够有效提高估计精度，但跳变线形系统的模式是互相转换的，也就是说目标的运动模式是时变的，在没有先验的情况下假设模式固定是不合理的。因此，本文提出方法的核心在于认为更为合理的混合估计方法应当基于目标当前状态而决定。为了实现对目标状态的准确估计，本文引入“切换”状态的概念。对于一个两模型交互的系统，目标状态包括S1、S1→2、S2和S2→1。显然，如果目标真实状态是已知的，则可以选择最佳的混合估计方法。假设目标真实状态为S1，则q(ek2)的最优选择为
|
(31) |
同样地，当目标真实模式为S2时，q(ek2)的最优选择为
|
(32) |
对于“切换”状态S1→2，由于在此过程中ek应置零，因此q(ek2)的最优选择同式(32)。当目标状态为S2→1时，意味着ek可能在其取值区间跳变，q(ek2)的分布应体现出此特性，故采用A3的均匀分布假设，即
|
(33) |
基于上述分析，可以发现A1、A2和A3分别是不同目标状态下的最优方法。因此，在已知目标状态时，3种方法可以被纳入一个统一的框架，根据当前状态选择最优的一个。然而，跟踪系统很难准确知道目标的真实状态。尽管如此，依然可以利用一些滤波信息，例如模型概率和新息，尽可能地确定目标状态。
相比于状态S1和S2可以通过模型概率μki体现，“切换”状态由于模式切换导致的机动，往往体现在新息变化上。因此，为了检测“切换”状态S2→1，本文采用一种常用于目标跟踪的机动检测器[15-17]来确定目标当前状态。显然，当目标从模式2切换到模式1时，由于模型失配将导致模型2的新息增大。因此，令检验统计量为新息的衰减记忆平均值，计算式为
|
(34) |
式中:ρ(s)=(1-λ)/(1-λs+1)，0＜λ < 1为衰减因子; vk2和Sk2分别为模型2的新息及其协方差; s为滑窗长度; [k-s+1, k]为滑窗检测区间。在高斯假设下，由于vk2和Sk2都服从零均值的高斯分布，故(vk2)T(Sk2)-1vk2服从自由度为nz(量测向量维度)的卡方分布。作为一个独立同分布的高斯变量加权和，检验统计量εk虽然不是标准的卡方分布，但可以通过矩匹配被视作一个经过缩放的卡方变量[16-17]，有
|
(35) |
因此，检验统计量εk服从自由度为nz的卡方分布。令假设H0代表目标状态切换未发生，假设H1则代表模型切换发生，“切换”状态S2→1由式(36)检测：
|
(36) |
由式(36)可知，当εk超过检测阈值χ2nz(α)时拒绝假设H0，否则接受该假设。由概率论可知，虚警概率Pfa=α，即在目标状态未改变的情况下错误判定目标状态发生变化的概率为α。因此，置信水平1-α通常设的相当高(例如95%或99%)以降低虚警概率。当检测门限χ2nz(α)过高时，必然导致漏警概率Pla增大，即在目标状态已改变的情况下未成功检测。很显然，检测器往往需要在虚警概率和检测概率间做性能平衡，尽量降低误判造成的影响。
当检测器触发阈值，则状态S2→1被激活，利用A3重新估计滑窗区间[k-s+1, k]内的目标状态估计。当检验统计量εk低于式(36)阈值时，则直接利用模型概率确定目标状态，即
|
(37) |
当μk1>μk2时，认为目标状态为S1，故采用A2进行混合估计；否则认为目标状态为S2，采取A1。图 1给出了本文方法(记为A4)的具体流程。
|
| 图 1 本文方法原理框图 Fig. 1 Functional block diagram of proposed method |
如前所述，由于门限检测导致的一定的虚警及漏警。下面简要分析各种情况下A4的跟踪性能差异。当目标真实状态未发生变化(H0)，则正确接受H0时的跟踪精度趋近于A1或A2，否则虚警导致A4的跟踪精度趋近于A3。当目标真实状态发生变化(H1)，若正确接纳H1，A4则采用A3重新估计[k-s+1, k]区间的状态，其跟踪性能趋近于A3。相比于A1和A2，A4能够在成功检测到状态变化后再启用A3，这意味着A4可以利用量测数据的后验信息来校正目标状态，从而抑制因模型切换造成的跟踪误差剧增，极大地降低峰值误差，其作用类似于IMM平滑器[18]，区别在于前者只在检测器触发时作用。若错误地拒绝接纳H1，由于漏警发生会导致跟踪误差增大，跟踪性能接近A1或A2。此外，滑窗长度s的选取同样影响跟踪性能。若滑窗过长会导致检测延迟，不利于系统快速响应模式切换；滑窗过短则导致检测的可靠性降低。因此，在虚警率较低的条件下尽可能地选择相对较短的滑窗长度以快速响应系统模式切换，具体设置方法可参考文献[15-16]。
基于上述分析可以看出，本文方法并不固定采用某种混合方法，而是通过新息检验和模型概率确定目标状态以选取最佳方法，从而将3种混合方法的优势结合起来。
3 仿真与分析在本节中，对由不同运动模型组成的运动场景进行了仿真实验，用于比较本文所提方法和现有3种方法的性能。在仿真中，本文考虑了过程噪声的影响，并给出2种评价指标，最后还讨论了仿真结果。
3.1 轨迹与模型场景1：由2个常见的线形运动模型组成。其中，模型1和模型2分别为CA和CV模型，用于产生一个2D运动轨迹(即nz=2)。2种运动模型的状态向量分别为
|
(38) |
|
(39) |
式中:xk和yk为位置分量；
|
(40) |
|
(41) |
|
(42) |
式中:I2表示2阶单位矩阵; T=1 s为采样间隔。
场景2：由线形的CV模型和非线性的协同转弯(Coordinated Turn, CT)模型组成。CT模型的相应参数设置为
|
(43) |
|
(44) |
|
(45) |
|
(46) |
式中:ωk为转弯率。
目标初始状态为[1 000 m 0 m 100 m/s0 m/s]T。在0~50 s和101~150 s内做匀速运动，在51~100 s内，场景1目标做加速度为30 m/s2的匀加速运动，场景2做转弯率3.6 (°)/s的匀速转弯运动。目标的真实轨迹如图 2所示。
|
| 图 2 场景1和2的真实轨迹 Fig. 2 True trajectories of Scenarios 1 and 2 |
对于上述2个场景，系统量测噪声方差为σm2=502 m2，其协方差为Rki=σm2I2；过程噪声协方差为Qki=σa2I2，考虑到过程噪声的大小对跟踪性能的影响，故设为
|
(47) |
除了混合方法之外，4种方法的IMM估计器参数设置是相同的，滤波算法采用卡尔曼滤波[19]。A1和A2的原理在第2节已介绍过。由于A3认为“额外”分量ek服从均匀分布，式(29)中场景1的加速度的取值范围为[-100, 100]m/s2，场景2中转弯率取值范围为[-3.6, 3.6](°)/s。对于A4，所需参数设置如下：衰减因子λ=0.95，置信度α=0.01；滑窗长度s=2，阈值χ22(α)=9.21。系统模式的Markov转移矩阵和初始模型概率分别为
|
(48) |
|
(49) |
本文的性能评价指标如下:
|
(50) |
|
(51) |
式中:M=500为蒙特卡罗仿真次数；N=150为目标运动时间；τk(σa)是过程噪声方差为σa2时的均方根误差(Root Mean Square Error, RMSE)；τσa为不同σa下的RMSE均值。
3.3 结果与分析限于篇幅，本文只选择展示仿真结果中的一部分内容，但这并不影响4种算法间的性能比较。对于场景1来说，图 3~图 6分别展示了过程噪声方差σa2为10-2和10-1时的位置RMSE以及速度RMSE。显然，A1、A2和A4在目标进行CV和CA运动时性能良好，但系统模式从CV切换到CA时，前2种算法的跟踪误差会急剧增大。A3和A4则能在模型切换时保持较好的跟踪精度。但A3的缺陷在于，当系统模式不发生切换时，它的跟踪精度是4种方法中最差的。随着过程噪声方差的增大，4种算法间的性能差别逐渐变小，这种趋势如图 7所示。同样地，对于场景2，图 8和图 9给出了过程噪声方差σa2为100时的位置RMSE以及速度RMSE。图 10给出了不同过程噪声方差时的位置RMSE均值，容易看到场景2的仿真结果类似于场景1。
|
| 图 3 场景1位置均方根误差(σa2=10-2) Fig. 3 Position RMSE of Scenario 1 (σa2=10-2) |
|
| 图 4 场景1速度均方根误差(σa2=10-2) Fig. 4 Velocity RMSE of Scenario 1 (σa2=10-2) |
|
| 图 5 场景1位置均方根误差(σa2=101) Fig. 5 Position RMSE of Scenario 1 (σa2=101) |
|
| 图 6 场景1速度均方根误差(σa2=101) Fig. 6 Velocity RMSE of Scenario 1 (σa2=101) |
|
| 图 7 不同过程噪声下场景1位置均方根误差均值 Fig. 7 Mean value of position RMSE of Scenario 1 under different process noises |
|
| 图 8 场景2位置均方根误差(σa2=100) Fig. 8 Position RMSE of Scenario 2 (σa2=100) |
|
| 图 9 场景2速度均方根误差(σa2=100) Fig. 9 Velocity RMSE of Scenario 2 (σa2=100) |
|
| 图 10 不同过程噪声下场景2位置均方根误差均值 Fig. 10 Mean value of position RMSE of Scenario 2 under different process noises |
由仿真结果可见，无论目标处于何种状态或过程噪声，A4都能保持较好的性能，因此本文所提算法A4的性能要优于其他3种算法。基于上述仿真结果，本文总结如下：
1) 过程噪声对不同混合方法有显著影响。当过程噪声较小时，系统模式切换导致的机动要求增大“额外”分量的协方差，这可以通过q(ek2)的设计实现。
2) 实际上，A3采用的均匀分布始终保持一个较大的协方差，因此对模型切换时的跟踪误差抑制起到了较好的作用，但相比于A1和A2所采取的策略，A3牺牲了目标状态处于稳态时的跟踪精度。
3) A4性能优越的原因在于，通过对目标状态的重新划分，并利用检测器对目标状态进行确认可以灵活地调整当前混合方法，使其具有了3种方法的优势。
4) 当过程噪声较大时，其协方差已足够大导致各方法间的差异逐渐消失。
4 结论本文提出了一种基于非等维状态的IMM混合估计方法，通过仿真结果验证了该方法具有良好的性能，得出以下结论：
1) 过程噪声的大小会影响IMM混合方法的性能，较小的过程方差会导致较大的模型切换误差，而较大的过程噪声方差会削弱这种影响。
2) 通过引入“切换”状态实现目标状态的重新划分，结合滤波信息可以实现对目标状态的确定，从而灵活地选取当前最佳的混合方法。
为使本文提出的算法能具有更好的适应性，仍需要对目标状态确认逻辑做更为深入理论分析。
| [1] |
LI X R, JILKOV V P. Survey of maneuvering target tracking.Part V.Multiple-model methods[J]. IEEE Transactions on Aerospace and Electronic Systems, 2005, 41(4): 1255-1321. DOI:10.1109/TAES.2005.1561886 |
| [2] |
MAZOR E, AVERBUCH A, BAR-SHALOM Y, et al. Interacting multiple model methods in target tracking:A survey[J]. IEEE Transactions on Aerospace and Electronic Systems, 1998, 34(1): 103-123. DOI:10.1109/7.640267 |
| [3] |
FOO P H, NG G W. Combining the interacting multiple model method with particle filters for manoeuvring target tracking[J]. IET Radar, Sonar & Navigation, 2011, 5(3): 234-255. |
| [4] |
翟岱亮, 雷虎民, 李炯, 等. 基于自适应IMM的高超声速飞行器轨迹预测[J]. 航空学报, 2016, 37(11): 3466-3475. ZHAI D L, LEI H M, LI J, et al. Trajectory prediction of hypersonic vehicle based on adaptive IMM[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(11): 3466-3475. (in Chinese) |
| [5] |
LI X R. Multiple-model estimation with variable structure.II.Model-set adaptation[J]. IEEE Transactions on Automatic Control, 2000, 45(11): 2047-2060. DOI:10.1109/9.887626 |
| [6] |
JAFFER A G, GUPTA S C. On estimation of discrete processes under multiplicative and additive noise conditions[J]. Information Sciences, 1971, 3(3): 267-276. DOI:10.1016/S0020-0255(71)80010-5 |
| [7] |
BLOM H A P, BAR-SHALOM Y. The interacting multiple model algorithm for systems with Markovian switching coefficients[J]. IEEE Transactions on Automatic Control, 1988, 33(8): 780-783. DOI:10.1109/9.1299 |
| [8] |
LI X R, BAR-SHALOM Y. Multiple-model estimation with variable structure[J]. IEEE Transactions on Automatic Control, 1996, 41(4): 478-493. DOI:10.1109/9.489270 |
| [9] |
LAN J, LI X R. Equivalent-model augmentation for variable-structure multiple-model estimation[J]. IEEE Transactions on Aerospace and Electronic Systems, 2013, 49(4): 2615-2630. DOI:10.1109/TAES.2013.6621840 |
| [10] |
XU L, LI X R, DUAN Z. Hybrid grid multiple-model estimation with application to maneuvering target tracking[J]. IEEE Transactions on Aerospace and Electronic Systems, 2016, 52(1): 122-136. DOI:10.1109/TAES.2015.140423 |
| [11] |
BAR-SHALOM Y, WILLETT P K, TIAN X. Tracking and data fusion[M]. Storrs: YBS Publishing, 2011.
|
| [12] |
YUAN T, BAR-SHALOM Y, WILLETT P, et al. A multiple IMM estimation approach with unbiased mixing for thrusting projectiles[J]. IEEE Transactions on Aerospace and Electronic Systems, 2012, 48(4): 3250-3267. DOI:10.1109/TAES.2012.6324701 |
| [13] |
GRANSTRÖM K, WILLETT P, BAR-SHALOM Y. Systematic approach to IMM mixing for unequal dimension states[J]. IEEE Transactions on Aerospace and Electronic Systems, 2015, 51(4): 2975-2986. DOI:10.1109/TAES.2015.150015 |
| [14] |
LOPEZ R, DANÈS P, ROYER F.Extending the IMM filter to heterogeneous-order state space models[C]//49th IEEE Conference on Decision and Control (CDC).Piscataway, NJ: IEEE Press, 2010: 7369-7374.
|
| [15] |
BAR-SHALOM Y, BIRMIWAL K. Variable dimension filter for maneuvering target tracking[J]. IEEE Transactions on Aerospace and Electronic Systems, 1982, 27(5): 621-629. |
| [16] |
JIN B, JIU B, SU T, et al. Switched Kalman filter-interacting multiple model algorithm based on optimal autoregressive model for manoeuvring target tracking[J]. IET Radar, Sonar & Navigation, 2014, 9(2): 199-209. |
| [17] |
LI X R, JILKOV V P.Survey of maneuvering target tracking.Part IV: Decision-based methods[C]//SPIE Proceedings Series.Society of Photo-Optical Instrumentation Engineers.Piscataway, NJ: IEEE Press, 2002: 511-534.
|
| [18] |
LOPEZ R, DANES P. Low-complexity IMM smoothing for jump Markov nonlinear systems[J]. IEEE Transactions on Aerospace and Electronic Systems, 2017, 53(3): 1261-1272. DOI:10.1109/TAES.2017.2669698 |
| [19] |
KALMAN R E. A new approach to linear filtering and prediction problems[J]. Journal of basic Engineering, 1960, 82(1): 35-45. DOI:10.1115/1.3662552 |