一种面向临近空间高超声速再入滑翔目标跟踪算法

何 山，吴盘龙，恽 鹏，李星秀

(1. 南京理工大学自动化学院， 南京 210094； 2. 南京理工大学理学院， 南京 210094)

0 引 言

近年来，随着航空航天技术的不断发展，临近空间高超声速飞行器以其飞行跨域大、速度快、机动强等复杂特性受到了世界各国的广泛关注。尤其是美俄等军事强国通过不断投入大量资源，先后开展了高超声速目标的多个关键技术验证，以期望作为精确打击武器使其成为未来空天作战中的新型战略威胁[1]。因此，对该类目标的跟踪预报与状态估计技术进行研究具有重要意义。

目前，针对临近空间高超声速目标跟踪均围绕跟踪模型的构建展开，并设计相应的滤波算法，归纳起来分为两类：1)动力学模型状态估计，即通过对目标动力学建模，将气动参数扩维到状态矢量中进行联合估计，其本质是机动的辨识，具有较高的估计精度[2-4]。但是在实际跟踪过程中，目标大部分参数是未知，需要通过辨识得到，当气动力变化较快或动力学模型不匹配时，跟踪能力会大幅下降，从而很难满足实际工程应用需求；2)运动学模型估计，即利用匀速(Constant velocity, CV)、匀加速(Constant acceleration, CA)、匀速转弯、Jerk、Singer以及“当前”统计等模型对目标的未知机动建模[5-7]，其算法设计简单，易于工程实现，但加速度估计精度有限。对此，通过构造输入估计器和回顾成本函数，有效地提高加速度的估计精度。

通常基于运动学模型估计的滤波算法大致可分为具有机动检测的跟踪算法和自适应跟踪算法两类[8]。其中，具有机动检测跟踪算法的基本思想是通过量测信息来观测目标运动的残差变化，并利用数理统计理论进行机动检测[9-12]。该类方法特点是不需要对目标的机动特性做任何先验假设，计算量小，但由于机动检测机制的存在而产生不可避免的时间延迟；自适应跟踪算法是通过对目标进行估计的同时对滤波器增益进行修正，有效地提高了跟踪过程的稳定性。

然而，临近空间高超声速飞行器机动的复杂性使得常规的运动模型很难与目标运动方式相匹配，从而影响跟踪精度甚至滤波发散。文献[13]对临近空间高超声速飞行器的飞行特性进行分析，提出一种基于CV + CA + Singer的交互多模型(Interacting multiple model, IMM)跟踪算法，通过不同模型的组合匹配出最优的目标运动模型，提高算法的跟踪精度。同时，为解决IMM算法实时性不高问题，文献[14]针对临近空间目标4种典型的非弹道式机动模式，设计了一种修正变结构IMM算法，相比传统变结构多模型算法，提高了模型切换速度。尽管多模型方法以其良好的跟踪性能备受青睐，但其设计过程复杂且实现时需要极大的计算资源，在许多资源受限的场合，基于自适应的单模算法往往也可以获得更好的跟踪性能。文献[15]针对临空飞行器机动特性强和周期滑跃的运动特点，提出了临空高速飞行器滑跃机动跟踪的Sine模型，更加高效地描述了高超声速目标滑跃式运动特性，但当目标不发生机动或突发的大机动时，该算法的跟踪精度会下降。

鉴于此，本文将输入估计的思想引入到跟踪模型中，将机动加速度看成是未知的确定输入；然后通过构建输入估计器使残差收敛到零，并利用回顾成本对未知加速度进行重构，采用递推最小二乘法估计出系统的参数；最后将估计的加速度引入到卡尔曼滤波的预测中，实现目标状态的有效估计。

1 跟踪建模

基于临近空间高超声速再入滑翔目标的质心运动方程可知[16]，气动力不仅随飞行高度、速度等变化，还与目标飞行攻角、倾角等姿态量密切相关，而临近空间高超声速飞行器再入滑翔段的强机动主要是受到气动力作用而使得加速度变化所致。因此，为了准确地描述目标状态随着时间变化的过程，本文将机动加速度看成是未知的确定输入，构建目标的运动方程

Xk=Fk-1Xk-1+Ck-1uk-1+Wk-1

(1)

式中：Xk∈Rs为k时刻的状态向量，Fk-1∈Rs×s为状态转移矩阵，Ck-1∈Rs×d为机动输入矩阵，uk-1∈Rd为机动大小，Wk-1是均值为零且方差为Qk-1的高斯白噪声。

在雷达探测系统中，雷达的量测信息通常是建立在三维球坐标中，目标与雷达之间的相对位置如图1所示。

图1 雷达测量坐标系Fig.1 Radar measurement coordinate

雷达的量测值Zk∈Rp主要包括径向距离rk、高低角θk和方位角ηk，其测量方程为

(2)

式中：xk,yk和zk是目标在测量直角坐标系下相对各方向的坐标位置，Vr,k,Vθ,k和Vη,k是相互独立、均值为零且恒定方差为σr,k,σθ,k和ση,k的高斯白噪声。

球坐标可通过下式转化到笛卡尔坐标系下的量测

(3)

因此，根据式(2)的量测值Zk，对真实均值和协方差矩阵求数学期望得到无偏量测偏差μk和协方差Rk

(4)

式中:各变量取值参照文献[17]。

(5)

2 基于RCIE-UCMKF的目标状态估计

2.1 输入估计器

为了准确估计出当前时刻的状态，可构造如下的输入估计器[10-12,18-19]

(6)

(7)

式中：Mi,k∈Rd×d，Ni,k∈Rd×p，n是输入估计器的阶数。将式(7)改写成矩阵的形式

(8)

为计算出参数向量Φk，选取最小二乘法的指标函数

(9)

根据最小二乘法的递推可得到

(10)

(11)

(12)

2.2 RCIE算法

为了便于计算，假设Fk=F,Ck=C,Hk=H，考虑到系统矩阵F,C和H是已知的，且(F,H)是能观的，定义系统的马尔可夫参数为

Li=HFi-1C∈Rp×di≥1

(13)

取任意正整数r，对于所有k≥r，则有

(14)

其中

(15)

(16)

(17)

(18)

其中，

(19)

然后，将式(17)～(19)中k用k-kj(j=1,…,m)替代，其中0≤k1≤…≤km，则有

(20)

(21)

(22)

定义扩展性能变量：

(23)

并将式(21)代入到(23)中有

(24)

式中：

(25)

(26)

(27)

(28)

(29)

用式(29)减去式(24)可得

(30)

最后，定义回顾成本函数为[18-19]

(31)

(32)

(33)

2.3 RCIE-UCMKF算法

为了准确的估计出系统的状态量，在无偏转换Kalman滤波的框架下，利用回顾成本输入估计和递推最小二乘计算出未知机动大小，整个算法的框图如图2所示。

图2 RCIE-UCMKF算法框图Fig.2 The structure of RCIE-UCMKF algorithm

根据前述部分的推导，所提RCIE-UCMKF算法的完整步骤为：

(34)

Pk+1|k=FPkFT+Qk

(35)

(36)

(37)

4)状态估计和误差协方差估计为

(38)

Pk+1=Pk+1|k-Kk+1HPk+1|k

(39)

3 仿真校验

3.1 仿真场景

本文以HTV-2为参考对象，对临近空间高超声速飞行器再入滑翔段进行仿真试验。根据文献[16]，假设目标为质点，忽略地球自转及非球形摄动因素等影响。飞行器质量907.2 kg，特征参考面积0.4837 m2，攻角10°，倾斜角0.5°，阻力系数和升力系数选取参考文献[12]，目标的起始位置为(6450.2450 km, 139.270 °, 35.270 °)，初始速度为6000 m/s，航迹倾角2.0 °，偏航角为155.00 °，则目标的飞行轨迹如图3所示。

图3 目标运动轨迹Fig.3 Target trajectory

为便于计算，假设地球是球体，且半径为6371.3930 km，测量雷达位置(6372.3930 km, 90.50°, -40.0°)，雷达的采样间隔T为1 s，距离误差为100 m，方位和俯仰角误差为0.5 mrad，则HTV-2在雷达界面的显示如图4所示。

在构建的三维仿真场景中，目标状态和输入向量为

(40)

(41)

其余的仿真参数设置如表1所示。

表1 仿真参数设置Table 1 The simulation parameters

3.2 仿真分析

根据上述的仿真场景，RCIE-UCMKF算法对目标位置估计和目标真实轨迹的对比如图5所示，可以看出所提算法能够对目标进行有效地跟踪。

图5 目标位置估计Fig.5 Estimation of target position

为了测试该算法的跟踪精度，本文分别用CKF[20]、UCMKF[17]、IE-UCMKF[9]、IMM-UCMKF[13]和RCIE-UCMKF算法在Matlab2014环境下进行200次蒙特卡罗实验。图6、图7展示了5种算法的位置均方根误差和速度均方根误差，可以看出，RCIE-UCMKF算法相比于其他4种算法，无论是位置还是速度都具有较高精确度和良好的性能。

图6 位置均方根误差Fig.6 RMSEs of the position

图7 速度均方根误差Fig.7 RMSEs of the velocity

CKF、UCMKF、IE-UCMKF、IMM-UCMKF和RCIE-UCMKF算法的各项算法性能如表2所示。由于CKF通过一组具有权重的采样点集近似计算所需的一、二阶矩，避免了非线性量测模型的线性化处理，而UCMKF是将量测通过坐标变换转换成直角坐标系中的伪线性形式，然后估计转换量测误差的前两阶矩，所以CKF和UCMKF算法在状态估计精度上相差不大，但是在计算量上，由于CKF需要求容积点，计算量远远高于UCMKF算法。RCIE-UCMKF算法的运行时间虽高于UCMKF，但却低于CKF和IMM-UCMKF。这是因为比起UCMKF算法，本文所提算法需要利用回顾成本和最小二乘计算每一时刻的加速度，因此时间会长于UCMKF算法。

表2 各个算法的性能Table 2 The performance of each algorithm

4 结 论

针对临近空间高超声速再入滑翔目标跟踪问题，本文提出了一种RCIE-UCMKF算法，蒙特卡罗试验验证了该算法的有效性和可行性。然而，实际跟踪过程中，由于该类目标的飞行特性使得量测数据出现一定的延迟和丢失的现象，因此后续的工作是考虑在该类情况下，提高算法的稳定性，并在一定的基础上降低算法的复杂度，使得该算法可应用于预警探测系统等工程领域中。