恽鹏 李星秀 吴盘龙 何山

摘要：为了解决在实际系统中因野值干扰带来的高超声速飞行器跟踪精度下降的问题，提出了一种交互多模型变分贝叶斯滤波算法（IMM-VB），该算法通过子模型权重与马尔可夫转移矩阵获取子模型的状态预测值。随后采用具有重尾特性的学生t分布取代高斯分布来描述量测模型，并利用VB算法实现子模型的量测协方差与状态的联合估计。最后在交互式多模型（IMM）框架下更新子模型权重与目标状态的融合输出。仿真结果表明，在野值观测条件下该算法比IMM算法具有更高的跟踪精度。

关键词：高超声速飞行器;野值;学生七分布;重尾分布;IMM-VB

中图分类号：TP24 文献标志码：A

高超声速飞行器具有飞行速度快、机动灵活、突防能力强、有效载荷大等优点，因此受到了世界各国的关注。目前而言，交互式多模型（IMM）算法是解决高超声速飞行器跟踪的有效途径之一。然而该算法是在广义伪贝叶斯理论基础上形成的，其前提是量测噪声与过程噪声均为高斯分布。在实际的高超声速飞行器跟踪过程中，由于大气、云层等因素的干扰，噪声的统计特性往往不符合高斯分布[1～3]，导致IMM算法的跟踪精度较差，无法满足实际系统的需求。

为了提高实际工程中的目标跟踪精度，许多研究学者对当前的跟踪算法进行了改进。为了解决重尾噪声线性系统的濾波问题，参考文献[4]提出了一种学生t滤波，该滤波基于噪声的重尾特性，将噪声建模为学生t分布，有效地克服了野值的干扰，提高了算法的鲁棒性。然而当过程噪声与量测噪声的自由度参数不同时，该滤波无法实现状态的更新。为了处理野值干扰导致的非高斯重尾噪声下的目标跟踪问题，参考文献[5]提出了一种最大熵卡尔曼滤波器，然而在误差协方差估计中缺少足够的理论基础，导致该滤波器的跟踪精度不佳。参考文献[6]利用变分贝叶斯法联合估计目标的状态分布与噪声分布，有效地提高了目标的跟踪精度，然而该算法是基于单模型设计的，因此不适合跟踪具有强机动性的目标，如高超声速飞行器等。

IMM因其计算量小、跟踪精度高等优点被广泛应用于各类飞行器的跟踪。参考文献[7]将常加速度模型（CA）/匀速模型（CV）/SINGER模型相结合，并基于径向速度估计补偿技术，有效地提高了临近空间高超声速飞行器的跟踪精度。参考文献[8]将IMM与容积卡尔曼滤波[9]（CKF）相结合，提出了一种IMM-CKF算法，用以提高高速滑翔飞行器的跟踪精度。参考文献[10]将标签伯努利[11]（ LMB）与无迹卡尔曼滤波[12]（UKF）相结合，提出了一种广义LMB-UKF算法，实现了多个高超声速飞行器的检测与跟踪。然而，上述算法将噪声均建模为高斯分布，当噪声为非高斯分布或者噪声协方差与建模时给定数值差距过大时[13～15]，都将严重降低算法的跟踪精度。因此不适用于野值干扰下高超声速目标的跟踪。

本文提出了一种IMM-VB算法用来解决野值观测条件下的高超声速目标跟踪问题。首先将CA模型与当前统计模型（CS）[16]进行交互，获取相应的预测值与权重，随后基于噪声的重尾特性，将量测噪声建模为学生t分布，利用VB算法，实现状态与量测噪声的联合估计。最后利用IMM-VB算法和IMM算法对高超声速飞行器进行跟踪。仿真结果表明，IMM-VB算法比IMM算法具有更高的跟踪精度。

1 基础算法

1.1 IMM算法

IMM算法是一种基于马尔可夫转移矩阵的广义伪贝叶斯算法[17，18]。该算法根据目标的先验信息，首先给予相应子模型的状态分量与权重，随后，对子模型进行滤波，同时根据量测信息与基于各个模型更新后的状态对不同模型赋予相应的权重，并通过加权求和更新目标状态。针对机动性较为复杂的目标，IMM算法相较于单模型算法具有更好的跟踪精度。下面介绍一下IMM算法的4个步骤：状态交互、状态滤波、模型概率更新和状态融合输出。

1.1.1 状态交互

假设k时刻各个子模型的目标状态、协方差、权重分别为X_j（k|k）、P_j（k|k）、μ（k|k），j为子模型的编号，对目标各个子模型进行交互运算，获取各个子模型的状态X_j⁰（k|k）与协方差P_j⁰（k|k）：式中：N为模型总数，X_j（k|k）为模型i到模型j的转移权重，π_ij为模型i到模型j的马尔科夫转移概率。

1.1.2 状态滤波

状态滤波过程与单模型滤波过程相同，依据卡尔曼滤波即可实现子模型下的状态更新。将X_j⁰（k|k），P_j⁰（k|k）作为子模型的输入进行状态预测与更新：式中：X_j（k+1|k）为状态预测值P_j（k+1|k）为预测协方差，K_j（k+1）为滤波增益，X_j（k+1|k+1）t状态更新值，P_j（k+1|k+1）为协方差更新值，F_j（k）为状态转移矩阵，H_j（k+1）为量测矩阵Q_j（k）为过程噪声协方差，R_j（k+1）为量测噪声协方差。

1.1.3 模型概率更新

通过滤波过程中获取的子模型信息与信息协方差，构建量测与模型匹配的似然函数∧_j（k+1）：式中：m为量测信息中的元素数目，子模型信息v_j（k+1）与信息协方差S_j（k+1）分别为：

通过状态交互后的子模型权重以及相应的似然函数更新模型概率μ_j（k+1|k+1）：

1.1.4 状态融合输出

对各子模型的模型概率与状态信息进行加权求和，获取融合后的状态X（k+1|k+1）与协方差P（k+1|k+1）：

1.2 VB算法

当实际的噪声不满足高斯分布时，贝叶斯公式很难推断出目标的后验分布。而VB算法通过迭代逼近的手段可以获取近似的后验分布[19]，有效地解决了非高斯条件下状态后验分布求解的问题。为了衡量近似分布的逼近程度，需要利用Kullback-Leibler散度（KLD）函数：式中：KL（q（Ψ））≥0且仅q（Ψ）=p（Ψ|X），p（Ψ|X）为参数集的后验分布，q（Ψ）为参数集的近似分布，且KL（q（Ψ））越小，表示q（Ψ）的逼近程度越高。令下界函数L（q（Ψ））=

P（Ψ，X）表示Ψ和X的联合概率密度函数，由式（16）可以看出，若想使KL（q（Ψ））最小化，则需要使L（q（Ψ））越大越好。基于平均场理论可以求得满足L（q（Ψ））最大化的参数近似分布：式中：E_Ψi-[1np（Ψ，X）]表示除参数Ψ_i外1np（Ψ，X）的数学期望，Ψ_i为参数集Ψ的第i个参数。在对式（17）进行求解的过程中，需要使用到上一次迭代步骤中的参数分布，因此，求解式（17）是一个不断迭代更新的过程，此时可以设置相应的迭代步数或者对相邻迭代的下界函数差值设置阈值来终止迭代过程。最后迭代过程获取的参数分布即为VB求解出的近似分布。

2 IMM-VB算法

由于量测信息存在野值的干扰，导致IMM算法的跟踪精度下降，本文将IMM与VB相结合提出一种IMM-VB算法，可以有效地克服跟踪过程中的野值干扰，提高高超声速飞行器的跟踪精度，其算法流程如图1所示。

该算法本质上还是基于IMM思想进行，其状态交互过程与IMM算法相同，这里就不再阐述。在滤波过程中由于量测信息存在野值干扰，导致其统计特性呈重尾分布，因此将量测噪声建模为学生t分布，可得子模型的似然函数式中： St（;μ，Σ，v）表示均值为μ、协方差为Σ、自由度参数为v的学生t分布，而学生t分布是高斯分布与伽马分布乘积的积分和，因此可以将式（18）改写为：式中：N（;μ，Σ）表示均值为μ，协方差为Σ的高斯分布，G（;α，β）表示形状参数α、尺度参数为β的伽马分布。引入辅助随机变量λ（k+1），并将其建模为伽马分布，则可以从式（19）获取概率密度函数p（λ_j（k+1））和p（z（k+1）|X_j（k+1），λ_j（k+1））：

利用式（20）、式（21）可以获取联合概率密度函数

结合式（4）、式（5）、式（20）、式（21）可以计算出联合概率密度函数p（λ_j（k+1），X_j（k+1），z_1：k+1）：

利用式（23）可求取1np（λ_j（k+1），X_j（k+1），z_1：k+1）式中：n为状态信息中的元素数目，利用式（17）、式（24）迭代获取相应参数的近似分布：式中：C为保证qⁱ⁺¹（λ_j（k+1》积分为1的常量，后面的C均指代这样的常量，但是对于每一个公式而言，其值可能不同。由式（27）可得：同理，可获取状态的近似分布：

由式（32）可推導出：

由高斯乘积定理可得：

进行完子模型滤波后，在模型概率更新过程中，由于辅助随机变量的存在，导致子模型的新息协方差发生改变，此时似然函数∧_j（k + 1）为：

最后，随后通过式（12）～式（14）、式（35）、式（36）、式（39）可以获取融合后的目标状态信息。

3 仿真试验

本文对高超声速飞行器[20]跟踪进行了仿真试验。以CA、CS为交互子模型，分别设计了相应的IMM算法和IMM-VB算法对该目标进行跟踪。其中目标的状态由二维坐标中的位置、速度、加速度组成：

马尔可夫转移矩阵为：

量测噪声满足：

式（42）表明野值干扰概率为0.1，且给出了有野值和无野值下的量测噪声分布。过程噪声满足：

当子模型为CA时：当子模型为CS时：式中：采样时间T=1s，a_max为目标加速度最大绝对值，a（k）为上一时刻目标加速度估计值，目标机动频率α取值1/60。辅助随机变量的自由度参数v取值50，初始时刻子模型的权重均为0.5，目标初始状态和协方差分别为：

根据所给的条件设计相应的算法，获取试验信息。图2展示了高超声速飞行器的真实轨迹和两种算法的估计轨迹。

图3、图4展示了两种算法在不同坐标轴上的位置跟踪均方根误差，图5、图6展示了两种算法不同坐标轴上的速度跟踪均方根误差。可以看出，IMM-VB相较于IMM算法具有更高位置和速度跟踪精度。

4 结束语

针对野值干扰下的高超声速目标跟踪问题，本文提出了一种IMM-VB算法，并通过仿真试验验证了该算法相较于IMM算法更适用于跟踪野值干扰下的高超声速目标。然而由于VB的存在，导致算法的计算时间过长，因此后续的研究过程中，应考虑如何降低算法的复杂度。另外，该算法只针对于单目标，因此如何设计针对于多目标跟踪的算法，也是后续工作的研究重点。

参考文献

[1]Agamennoni G，Nieto J I，Nebot E M.Approximate inferencein state-space models with heavy tailed noise[J].IEEE Transac-tions on Signal Processing，2012，60（10）：5024-5037.

[2]Huang Y L，Zhang Y G，Li N，et al.Robust Students t basednonlinear filter and smoother[J].IEEE Transactions on Aero-space and Electronic Systems，2016，52（5）：2586-2596.

[3]Karlgaard C D，Schaub H.Huber-based divided difference fil-tering[J].Journal of Guidance，Control and Dynamics，2007，30（3）：885-891.

[4]Roth M，Ozkan E，Gustafsson F.A students t filter for heavytailed process and measurement noise[C]//（ICASSP），2013IEEE International Conference on IEEE，2013：5770-5774.

[5]Chen B D，Wang J J，Zhao H Q，et al.Convergence of a fixed-point algorithm under maximum correntropy criterion[J].IEEESignal Processing Letters，2015，22（10）：1723-1727.

[6]Huang Y L，Zhang Y G，Li N，et al.A novel robust studentst-based Kalman filter[J].IEEE Transactions on Aerospace andElectronic Systems，2017，53（3）：1545-1554.

[7]Zhang X Y，Wang G H，Song Z Y，et al.Hypersonic sliding targettracking in near space[J].Defence Technology，2015，11（4）：370-381.

[8]Fan Y，Lu F，Zhu W X，et al.A hybrid model algorithm forhypersonic glide vehicle maneuver tracking based on the aerody-namic model[J].Applied Sciences，2017，7（2）：159.

[9]Arasaratnam I，Haykin S.Cubature kalman filters[J].IEEETransactions on Automatic Control，2009，54（6）：1254-1269.

[10]Cui H，Liu W，Zhu S，et al.A hypersonic vehicle tracking algo-rithm based on the UKF generalized labeled multi-Bemoulli fil-ter[C]//Control Conference（CCC），2016 35th Chinese，D”？E，2016：4911-4916.

[11]Reuter S，VoBT，Vo B N，et al.The labeled multi Bernoullifilter[J].IEEE Transactions on Signal Processing，2014，62（12）：3246-3260.

[12]Julier S J，Uhlmann J K，Durrant-Whyte H F.A new approachfor filtering nonlinear systems[C] //American Control Confer-ence，Proceedings ofthe 1995，IEEE，1995（3）：1628-1632.

[13]Wang C C，Wu P L，Bo Y M，et al.Maneuvering target trackingwith adaptive debiased converted measurement filter[J].Journalof Chinese Inertial Technology，2017，25（1）：136-140.

[14]Deng Y H，WuPL，LiXX，et al.Decoupled unbiased con-verted measurements adaptive gating CPHD filter for airbornemulti-target tracking[J].Journal of Aeronautics，Astronauticsand Aviation，2017，49（1）：123-131.

[15]Yu M J.INS/GPS integration system using adaptive filter forestimating measurement noise variance[J].IEEE Transac-tions on Aerospace and Electronic Systems，2012，48（2）：1786-1792.

[16]Jiang B，Sheng W，Zhang R，et al.Range tracking methodbased on adaptive“current”statistical model with velocity pre-diction[J].Signal Processing，2017，131：261-270.

[17]Cao H J，Zhang X Y，Wu P L.Combining theIMM with secondorder debiased converted measurement for detonation control[J].International Journal of Computational Science and Engineering，2017，14（1）：84-91.

[18]Wu P L，LiXX，Zhang L Z，et al.Tracking algorithm withradar and infrared sensors using a novel adaptive grid interactingmultiple model[J].IET Science，Measurement&Technology，2014，8（5）：270-276.

[19]Zhu H，Leung H， He Z S.A variational Bayesian approach torobust sensor fusion based on student-t distribution[J].Informa-tion Sciences，2013，221：201-214.

[20]何廣军，李槟槟，何其芳，等.基于IMM模型的周期性跳跃运动高超声速飞行器跟踪算法[J].弹箭与制导学报，2015，35（01）：37-40.