基于交互式多模型鲁棒滤波的目标机动估计算法*

彭一洋，奚 勇，陈 辉，孙 逊，仲科伟

0 引 言

目标加速度估计是机动目标跟踪领域的研究热点，对于提高导弹的制导精度和弹道品质具有重要意义[1,2].目标加速度估计的有效性和可靠性取决于数学模型的匹配度以及滤波算法的估计精度.交互式多模型算法(interacting multiple model，IMM)能够较好的解决单一模型难以匹配复杂目标运动问题.通过预先设定运动模型集，IMM算法面对复杂机动目标具有较强的鲁棒性[3].

滤波算法是影响IMM算法估计精度的重要因素之一，目标加速度估计是一个典型的非线性滤波问题.传统的扩展卡尔曼滤波[4](extended Kalman filter，EKF)在面对强非线性系统时会引入较大的截断误差，同时计算雅克比矩阵也会提高滤波器的设计难度.无迹卡尔曼滤波[5](unscented Kalman filter，UKF)、容积卡尔曼滤波[6](cubature Kalman filter，CKF)等基于确定性采样的非线性滤波方法虽然不需要计算雅克比矩阵，避免了线性截断误差.但UKF在面对高维系统时容易出现采样点权值为负的情况，CKF无法准确计算某些简单多项式的权值积分.高阶容积卡尔曼滤波[7](high-degree cubature Kalman filter，HCKF)是近几年提出的一种确定性采样滤波方法.相比于传统CKF，HCKF将球面-径向容积准则的阶次提高为五阶，避免了传统CKF无法计算某些简单多项式权值的缺点，滤波稳定性更好，精度更高.文献[8]将高阶容积卡尔曼滤波应用于极区惯性/重力梯度组合导航.文献[9]针对传统的容积卡尔曼滤波(CKF)估计精度有限的问题，提出了一种基于任意阶容积规则的高阶容积卡尔曼滤波方法并应用于机动目标跟踪问题.

以上滤波算法皆假设测量噪声具有高斯特性.但在实际应用中，由于目标不同部位对雷达波的反射能力不同，雷达导引头对同一目标的测量位置容易出现随机摆动，使导引头的测量噪声呈现出闪烁噪声的特性.在这种情况下，传统非线性卡尔曼滤波算法的估计精度会下降甚至发散.Huber提出的Huber-Based滤波能够在传统的高斯滤波理论框架内解决非高斯问题，对具有闪烁噪声的系统具有渐进最优鲁棒性[10].文献[11]提出了一种基于Huber理论的鲁棒广义HCKF，该方法在鲁棒性和估计精度上均具有优越的性能.文献[12]基于Huber理论提出了一种鲁棒差分滤波算法，数字仿真结果表明该滤波算法对混合高斯噪声的精度更高.

本文针对雷达导引头测量噪声为闪烁噪声的问题，在HCKF的测量更新部分采用Huber滤波理论，提出高阶容积鲁棒滤波，结合交互式多模型算法框架，提出多模型鲁棒滤波算法.数字仿真结果表明，采用所提算法设计的目标机动估计器具有较强的鲁棒性，在面对闪烁噪声时，具有更高的滤波精度.

1 高阶容积鲁棒滤波

针对测量噪声为闪烁噪声的非线性滤波问题，采用高阶容积卡尔曼滤波进行时间更新，Huber-Based滤波进行测量更新，提出高阶容积鲁棒滤波方法.

将目标机动估计问题考虑为如下非线性离散系统:

xk=f(xk-1)+wk-1

(1)

zk=h(xk)+vk

(2)

式中，f和h分别为系统的状态方程和测量方程，x∈Rn是状态矢量，z∈Rm是测量矢量，wk-1是系统噪声矢量，vk是测量噪声矢量.

1.1 状态更新

高阶容积卡尔曼滤波采用五阶容积-径向准则计算高斯权值积分公式来确定容积点和权值，通过容积点在非线性状态方程中的传播进行状态更新.高阶容积鲁棒滤波的状态更新方法如下.

(3)

式中，n为状态向量的维数，点集ξi表达式为

(4)

(5)

(6)

(2)使用系统状态方程传播容积点为

(7)

(3)式(8)、式(9)为状态一步预测及其误差协方差矩阵

(8)

(9)

式(9)中，Qk-1为系统噪声协方差矩阵，权值wi表达式

(10)

(4)使用状态一步预测及其误差协方差计算容积点预测

(11)

(5)使用系统测量方程传播容积点

(12)

(6)计算测量一步预测

(13)

(7)计算协方差矩阵

(14)

(15)

1.2 测量更新

高阶容积鲁棒滤波采用Huber-Based滤波进行测量更新的方法，该方法将测量更新转化为线性回归问题，利用广义极大似然估计进行求解.下面给出高阶容积鲁棒滤波的测量更新方法.

令xk为k时刻目标状态的真实值，δk为状态估计误差.则有：

(16)

线性化测量方程可得：

(17)

式中，测量矩阵Hk=[(Pk/k-1)-1Pxz].

对式进行移项可得如下线性回归问题：

(18)

yk=Mkxk+ξk

(19)

式中，各变量表达式如式(20)～(21)所示.

(20)

(21)

(22)

(23)

测量更新即为求解如下代价函数极小值：

(24)

式中，ζi为ζ的第i个元素，ζ为残差.

ρ为目标函数，Huber证明了当ρ满足式(25)时，Huber-Based滤波对噪声为闪烁噪声的系统具备渐进最优鲁棒性.

(25)

式中，γ为可调因子.

使代价函数J最小的解应满足式(26).

(26)

定义权值函数ψ(ζi)和矩阵Ψ分别如式(27)～(28)所示

ψ(ζi)=ρ′(ζi)/ζi

(27)

Ψ(ζi)=diag{ψ(ζi)}

(28)

则式(26)可写成如式(29)所示的矩阵形式.

(29)

使用加权迭代法得到式(29)的解如式(30)所示

(30)

迭代初值如式(31)所示

(31)

状态估计误差协方差矩阵的计算方法如式(32)所示

(32)

2 IMM-HCHF算法

IMM算法需要预先设定不同的目标跟踪模型组成目标机动模型集.选取惯性系下目标的位置矢量、速度矢量、加速度矢量作为状态量，本文所采用的目标机动运动模型的离散形式如2.1节所示.

2.1 目标机动模型

本文选取常加速度模型、Singer模型[13]、“当前”统计模型[14]3种模型组成目标机动运动模型集，下面将给出一维条件下3种模型的离散形式，三维条件下可直接进行矩阵扩展.

(1)常加速度模型

常加速度模型认为目标进行匀加速运动，同时考虑实际系统中的加速度摄动，加入了过程噪声.常加速度模型的优点是形式简单，易于使用，缺点是假设条件过于简化，估计精度较低.

常加速度模型的离散形式状态方程如式(33)所示

(33)

(34)

(35)

式中，T为采样步长，q为加速度方差.

(2)Singer模型

Singer模型认为目标加速度服从零均值假设，即认为目标加速度在一定范围内为均匀分布.Singer模型面对大多数机动形式均能保证较高的精度，但由于该模型基于零均值假设，当目标的机动形式为非零均值时，模型的估计精度则会下降.

Singer模型的离散形式状态方程如式(36)所示

(36)

(37)

(38)

(39)

(3)“当前”统计模型

“当前”统计模型认为目标加速度服从修正的瑞利分布，当目标长时间存在较大机动加速度时，“当前”统计模型的估计精度明显高于Singer模型，但由于“当前”统计模型服从的是非零均值假设，在跟踪匀速目标时，该模型往往不能有效反映目标的机动情况.“当前”统计模型和Singer模型能够在一定程度上实现互补.

“当前”统计模型的离散形式状态方程如式(37)所示

(40)

(41)

式中，T为采样步长，α为机动频率，σ2为目标加速度方差，则“当前”统计模型的目标加速度方差σ2的表达式如式(42)所示

(42)

2.2 测量模型

(43)

(44)

qε=arcsin(Δy/r)

(45)

qβ=arctan2(-Δz,Δx)

(46)

(47)

(48)

式中，Δ表示目标和导弹之间的位置或速度差值.

2.3 IMM-HCHF算法

IMM算法的基本思想是选取不同的目标机动模型组成模型集，用于匹配不同的目标机动形式.在每个滤波周期内，各个模型对同一量测数据进行并行滤波处理，模型与模型之间通过时齐马尔科夫链进行转换，最后采用贝叶斯公式对各并行滤波器的估计结果进行数据融合作为滤波输出.

IMM算法中并行滤波器的估计精度将直接影响整个算法的精度，本文提出的IMM-HCHF算法将前文介绍的HCHF算法应用于传统IMM算法的并行滤波器中，其他流程不变.

假设模型集所含模型数目为n，则模型间的转移概率矩阵为马尔科夫矩阵Π.

(49)

IMM算法流程如下：

(1)输入交互

(50)

(51)

交互后各并行滤波器的状态估计和协方差矩阵输入如式(52)和式(53)所示.

(52)

(53)

(2)并行滤波器滤波

(3)模型概率估计

依据第②步中各滤波器的估计结果对模型更新模型概率.

(54)

(55)

(4)模型估计混合

利用更新后的模型概率对并行滤波器的输出进行加权融合得到最终的状态估计结果，状态估计如式(56)所示，状态估计协方差矩阵如式(57)所示

(56)

(57)

3 数字仿真与结果分析

3.1 仿真场景

本文选择战斗机为目标，采用2.2节建立的雷达导引头测量模型对目标进行跟踪，认为目标在雷达导引头开机后开始进行机动，分别在导引头的噪声为闪烁噪声和高斯白噪声时进行仿真.

闪烁噪声的概率密度函数和参数如式(58)所示.目标初始状态如表1所示，导引头测量误差如表2所示.

表1 目标初始状态Tab.1 Initial states of target

表2 导引头测量误差Tab.2 Noise characteristics of seeker

p(x)=(1-ε)p1+εp2

(58)

3.2 仿真结果与分析

在导引头开机后目标首先进行转弯机动逃离，然后进行加速机动，最后进行滚筒机动.目标的飞行轨迹如图1所示.进行100次蒙特卡洛数字仿真.当导引头测量噪声为闪烁噪声时，仿真结果如图2～图10所示.

图1 目标飞行轨迹Fig.1 The trajectory of target

图2 x向位置估计结果Fig.2 The estimate results of the target position in x direction

图3 y向位置估计结果Fig.3 The estimate results of the target position in y direction

图5 x向速度估计结果Fig.5 The estimate results of the target velocity in x direction

图6 y向速度估计结果Fig.6 The estimate results of the target velocity in y direction

图7 z向速度估计结果Fig.7 The estimate results of the target velocity in z direction

图8 x向加速度估计结果Fig.8 The estimate results of the target acceleration in x direction

图9 y向加速度估计结果Fig.9 The estimate results of the target acceleration in y direction

图10 z向加速度估计结果Fig.10 The estimate results of the target acceleration in z direction

从以上仿真结果中可以看出，在测量信息带有闪烁噪声时，IMM-EKF和IMM-HCKF的精度基本一致，IMM-HCHF的目标加速度估计精度明显高于IMM-EKF和IMM-HCKF，位置估计精度提高了约52%，速度估计精度提高了约45%，加速度估计精度提高了约41%.

Y向的速度和加速度估计均方误差在4.2 s～4.5 s突然增大，主要原因是在该时间段内发生了加速度值的突变，由于系统的惯性，滤波器的输出无法瞬间跟上目标状态的变化，导致误差瞬间增大，经过约0.4s的收敛后均方误差回到了合理水平.Z向发生的加速度估计均方误差突然增大现象的原因与Y向相同，但由于加速度幅值变化相对较小，对目标速度变化的影响较小，因此速度估计的均方误差并不明显.

当导引头测量噪声为高斯白噪声时，加速度的均方误差结果如图11所示.

图11 加速度估计结果Fig.11 The estimate results of the target acceleration

从仿真结果中可以看出，IMM-HCHF能够在量测噪声为高斯白噪声和闪烁噪声两种典型状态下均保持较高的精度，具有较强的鲁棒性和适应性，当导引头测量噪声为高斯白噪声时，IMM-HCHF的精度与IMM-EKF和IMM-HCKF的精度相当.具有较强的鲁棒性和适应性.

4 结 论

本文以导弹的雷达导引头对复杂机动目标的轨迹跟踪问题研究为背景，在IMM目标状态估计框架下，针对滤波过程中传统卡尔曼滤波无法应对闪烁噪声的问题，提出了基于Huber算法的高阶容积鲁棒滤波算法.选取做复杂机动的战斗机为状态估计对象，进行了三维数字仿真验证.仿真结果表明，该方法可以有效解决测量噪声为闪烁噪声时的目标加速度估计问题，估计精度高于IMM-EKF和IMM-HCKF，在面对高斯白噪声时具有与传统高斯滤波相当的滤波精度，具有较强的鲁棒性和适应性.