基于IMM－SRCDKF的高速机动目标运动状态估计算法研究

付 磊，张 龙，张 丁

（上海机电工程研究所，上海 201109）

0 引 言

随着各国对临近空间高速飞行器研究的不断深入，该类飞行器具备长时间大范围的机动能力，能够打破现有的战略平衡，已经成为未来战争中的重要武器之一。各国飞行试验的成功也标志着高速飞行器技术正在逐步成熟，其突出的战略意义对临近空间高速武器的防御技术提出了更高的要求［1－3］。因此，在发展临近空间高速飞行器技术的同时，也应该加强其防御策略的研究。

在对临近空间高速飞行器的拦截制导过程中，必然需要利用外部信息支援预警系统对机动弹道、运动信息等进行估计，并且需要在拦截器中末制导阶段为拦截器提供与目标间相对运动的精确信息［4］。针对高速机动目标运动状态估计问题，主要从构造目标模型、机动检测等不同角度来提高估计精度。交互式多模型方法（interaction multiple model，IMM）是目前工程中广泛应用的方法之一，它主要通过设计多种模型来代表系统所有可能的行为方式，并对各模型的估计结果进行融合，从而综合不同模型的优点。目前，IMM方法在应用中的实现主要是基于扩展卡尔曼滤波（Extended Kalman filter，EKF）算法。随着 EKF算法研究和应用的深入，其缺点也慢慢显现。如有时需要计算复杂的Jacobi矩阵，EKF算法精度不高，甚至发散。

针对这些不足，有研究人员提出了粒子滤波（particle filter，PF）等算法，但此类方法计算量过大，很难在工程中应用。而近些年出现的一类非线性滤波方 法，如 Sigma－Point卡 尔 曼 滤 波 （Sigma－Point Kalman filter，SPKF），以其运算量小、计算稳定、精度高和不需要计算Jacobi矩阵等优点，渐渐被工程技术人员所关注，正成为替代EKF的较为理想的算法。SPKF方法利用加权统计线性化回归技术，通过一组确定性采样点（Sigma点）来捕获系统的相关统计参量，从而避免了EKF方法中的Jacobi矩阵计算，也比PF方法节约了大量运算过程。根据Sigma点选取的不同，SPKF方法主要分为无迹卡尔曼滤波（unscented Kalman filter，UKF）和中心差 分卡尔曼滤 波（central difference Kalman filter，CDKF）。与UKF相比，CDKF具有理论精度稍高和更易实现等特点。为克服舍入误差可能引起的滤波发散，有学者进一步提出了平方根（square root，SR）形式的CDKF方法，保证了协方差矩阵在传播过程中的半正定性，从而具有良好的数值稳定性。

本文针对临近空间高速目标的运动状态估计问题，结合交互多模型算法提出了IMM－SRCDKF目标信息估计算法，通过对比仿真验证了该算法对临近空间高速机动目标具有更好适应性的特点。

1 目标运动模型

对于临近空间高速目标，目标运动模型的准确性是影响目标跟踪性能的关键因素，建立符合实际目标飞行特性的目标运动模型以及便于数学计算的模型是实现目标信息提取的首要条件。由于目标机动状态具有较多的外界因素以及人为控制的影响，无法事先确切知道目标机动情况，因此建立和选取合适的目标运动模型是实现目标运动状态估计的难点之一。

“当前”统计模型的主要思想是当目标以某一加速度做机动时，在下一时刻其机动加速度只可能在该机动加速度的邻域范围内改变。因此，对于不同的目标机动情况，只需考虑目标在当前状态的机动加速度的概率密度，并不需要考虑目标机动加速度变化的所有可能性。相关研究和应用都表明“当前”统计模型能更准确地描述目标机动模型［5］。

以临近空间高速飞行器为例，“当前”统计模型假设飞行器机动加速度服从修正的瑞利分布，且均值为当前时刻的机动加速度，并利用机动加速度的均值与方差之间的关系来建立机动加速度的均值和方差自适应算法。

目标运动模型离散状态方程方程如式（1）所示。

式中：X（k｜k）为k时刻目标状态矩阵；

X（k＋1｜k）为k＋1时刻状态预测矩阵；

W（k）为系统噪声矩阵，其协方差矩阵为

2 观测模型

根据雷达的工作原理，雷达观测模型建立在雷达球坐标系，以距离R、方位角q1、高低角q2描述目标位置。该坐标系与雷达直角坐标系Orxryrzr重合，是一种常见的非惯性坐标。雷达直角坐标系和飞行器的地理坐标系相似，坐标原点Or为雷达站中心，坐标轴Orxr和Oryr是地球参考球体的切线，分别指向地理东向和北向，而Orzr轴垂直于水平面向上，因此该坐标系也被称为东北天地理坐标系。地心坐标系与雷达直角坐标系之间的关系如图1所示，雷达直角坐标系与雷达球坐标系的关系如图2所示。

图1 地心坐标系与雷达直角坐标系Fig．1 Geocentric coordinate system and Radar rectangular coordinate system

图2 雷达直角坐标系与雷达球坐标系Fig．2 Radar rectangular coordinate system and radar spherical coordinate system

对于图2雷达球坐标系，雷达测量值Z＝［r，q1，q2］T，由几何关系可得

3 平方根中心差分卡尔曼滤波器

相比于EKF算法，平方根中心差分卡尔曼滤波（square root central difference Kalman filter，SRCDKF）算法与UKF算法都是基于Sigma点的卡尔曼滤波算法。SRCDKF借助Sterling插值获得后验均值和协方差，不需要对系统的状态方程以及观测方程进行线性化处理，避免了扩展卡尔曼滤波算法中状态转移矩阵的雅克比矩阵求解过程，从而避免了由线性化造成的精度损失。同时，工程应用中数值计算会引入舍入误差，可能导致UKF算法与EKF算法在滤波过程中出现滤波发散或状态协方差矩阵非对称或负定情况，从而使滤波器失效。SRCDKF算法直接利用协方差的平方根矩阵进行迭代计算，提高了实际应用中的数值稳定性［6－7］。

SRCDKF的计算流程如下：

1）计算所用的初值

式中：x0／0为系统状态初值；为系统状态协方差平方根矩阵初值。

2）计算用作时间更新的σ点

状态变量预测值的均值和均方误差分别为

3）时间更新

扩张状态变量xaw的维数；对于高斯分布，半步长h的最优值为；qr（·）表示对矩阵进行QR分解运算。4）计算观察点更新σ点

5）测量更新

式中：yk＋1／k为k＋1时刻观测更新值；Yi，k＋1／k为k＋1时刻观测采样点i的更新值；为k＋1时刻观测协方差矩阵；Pxk＋1，yk＋1为状态、观测量互协方差矩阵；cholupdate表示一阶Cholesky分解。

职业高中是培养社会所需人才的基地，所以教学方式一定要满足时代发展的最新特征，特别是会计电算化和数字化。一来是因为现代先进的教学手段，可以有效吸引学生更专注的课堂注意力，提高职业高中会计教学的效率和质量；二来让学生在学校就掌握现代会计的操作流程，更有利于学生尽早适应社会，为学生日后的就业前景奠定夯实的基础。

4 交互多模型SRCDKF算法（IMM－SRCDKF）

交互多模型算法是由BLOM和BAR－SHALOM在20世纪80年代基于广义贝叶斯算法提出的一种具有马尔可夫转移概率的结构自适应算法。交互多模型算法作为一种软切换算法，能够对不同滤波器的估计值进行加权融合［8］，有效避免了单模型不能准确描述不同类型机动的问题，并且不需要对目标的实时机动参数进行测量，在工程上具有较强的实用性。本文将SRCDKF滤波器与IMM算法相结合，得到了一种适用于临近空间高速目标跟踪滤波的多模型目标运动信息估计算法。

交互模型算法的融合过程如图3所示。针对目标跟踪滤波问题，根据目标的机动能力以及特性，选择多个符合目标不同状态的机动模型来建立一个有限数量的模型集。交互多模型算法利用上一时刻不同模型滤波器的估计残差计算得到当前时刻不同模型的匹配概率值；对该匹配概率值归一化后，再对当前时刻不同模型滤波器的滤波估计进行混合，从而得到当前时刻状态估计；同时，利用该值计算获得下一时刻的不同模型滤波器的滤波初值。

图3 交互多模型算法流程图Fig．3 Flow chart of IMM algorithm

各模型之间的相互转换概率可以用状态转移概率矩阵表示为

设vj（k）为滤波器j所对应的滤波器的估计残差，表示滤波估计所对应的观测值与目标实际观测值的差值，相应的协方差矩阵为Sj（k），假设其服从高斯分布，则滤波器j的当前状态的概率为

式中：vj（k）＝Z（k）－H（k）（k｜k－1），Sj（k）＝H（k）Pj（k｜k－1）H（k）T＋R（k）。为了更方便地利用模型概率，需要对Λj（k）进行归一化处理，即uj（k）＝

最后，根据上述计算流程得到的不同滤波器的概率、状态估计值以及协方差矩阵，对不同模型滤波器的输出值进行加权融合，即可得到混合状态滤波值（k｜k）和混合状态方差P（k｜k）为

5 仿真对比

为验证所提出的临近空间高速飞行器运动状态估计方法的有效性，本文以Lockheed Martin公司的通用高超声速飞行器CAV－H的公开气动模型为研究对象，基于高斯伪谱法计算得到高速目标最优再入参考轨迹［9］，采用蒙特卡洛仿真进行检验。

对比仿真采用不同机动模型的EKF滤波算法、SRCDKF滤波算法以及多模型IMM－SRCDKF滤波算法的滤波结果，以目标的状态估计值和目标运动的实际值的均方误差（root mean square error，RMSE）数据作为评价指标分析估计结果。

仿真条件：目标初始位置为（－25，－1 500，80）km，以5 000m／s水平初始速度飞行，落点位置约束为（0，0，0）km，俯仰角约束为 －5°；假设雷达站位于地理坐标系原点，地面站雷达数据采样率为10 Hz，其测距精度为10 m（1σ），俯仰角及方位角的测角精度均为0．1°（1σ）。

5．1 不同滤波方法滤波效果对比

使用基于统计模型的EKF滤波算法和SRCDKF滤波算法对相同的目标轨迹、速度、加速度信息进行滤波估计，设模型机动频率参数为0．05，仿真结果如图4～12所示。

图4 观测目标轨迹Fig．4 The observed target trajectories

图5 滤波轨迹对比Fig．5 Contrast of filtered trajectories

图6 X轴向速度估计曲线Fig．6 Velocity estimate in Xdirection

图7 Y轴向速度估计曲线Fig．7 Velocity estimate in Ydirection

图8 Z轴方向速度估计曲线Fig．8 Velocity estimate in Zdirection

图9 X轴方向加速度估计曲线Fig．9 Acceleration estimate in Xdirection

图10 Y轴方向加速度估计曲线Fig．10 Acceleration estimate in Ydirection

图11 Z轴方向加速度估计曲线Fig．11 Acceleration estimate in Zdirection

图12 位置估计均方误差曲线Fig．12 Root－mean－square errors of position estimate

由仿真结果可以得出，目标在整个飞行过程中具有较大的机动特征不确定性，EKF和SRCDKF滤波算法均能对目标信息有效地做出估计。在100 s到200 s之间，就目标机动产生的估计误差而言，SRCDKF算法的估计误差比EKF算法的估计误差更小，且SRCDKF收敛速度更快。由图13～14可知，在600 s后，随着飞行机动特性减弱，飞行速度和加速度降低，SRCDKF算法的速度、加速度均方根误差降低，精度提高。

图13 速度估计均方误差曲线Fig．13 Root－mean－square errors of velocity estimate

5．2 多模型与单模型滤波效果对比

设定IMM－SRCDKF滤波算法分别采用3个不同机动频率的“当前”统计机动模型组合，其机动频率分别为0．01、0．07、0．2，单模型SRCDKF算法机动频率为0．07。仿真结果如图15～18所示。

图14 加速度均方误差曲线Fig．14 Root－mean－square errors of acceleration estimate

图15 位置估计均方根误差曲线Fig．15 Root－mean－square errors of position estimate

图16 速度估计均方根误差曲线Fig．16 Root－mean－square errors of velocity estimate

图17 加速度估计均方根误差曲线Fig．17 Root－mean－square errors of acceleration estimate

图18 IMM－SRCDKF模型概率分布Fig．18 Probability distributions of IMM－SRCDKF model

由图18可知，IMM－SRCDKF算法能够自动选择与当前实际机动状态相符的机动模型，与单模型SRCDKF滤波算法相比，在目标进行大机动时其精度更高。在0～400 s内，其位置、速度、加速度均方根误差明显较低，位置均方根误差峰值由单模型的2 889 m下降到2 415 m，误差均值由452 m 下降到385 m。在200 s、400s等几个机动结束时刻，IMM－SRCDKF算法估计得到的速度和加速度均方根误差均能以更快的速度收敛到零值附近，对不同大小的机动具有更强的鲁棒性。

5 结束语

本文针对临近空间高速目标给出了目标机动模型，并建立了雷达对目标飞行器的跟踪和观测模型，提出了基于SRCDKF方法的目标信息提取算法。通过仿真，验证了本文提出的IMM－SRCDKF目标信息估计算法的估计精度。对比EKF算法以及单模型SRCDKF算法，本文提出的算法具有收敛速度快、精度高特点。同时，基于SRCDKF的目标信息提取算法能够有效避免EKF算法所需的Jacobi矩阵计算，同时也避免了由数值计算舍入误差导致的状态协方差矩阵负定产生的滤波失效问题。