高超声速滑翔目标自适应跟踪方法

黄景帅，李永远，汤国建,*，包为民, 3

1. 国防科技大学 空天科学学院，长沙 410073 2. 中国运载火箭技术研究院，北京 100076 3. 中国航天科技集团有限公司，北京 100048

与传统弹道式目标不同，高超声速滑翔目标(Hypersonic Glide Target, HGT)具有大升阻比的气动外形，在临近空间以大于5的马赫数飞行，依靠气动力实现远距离滑翔和大范围机动[1]。随着HGT相关技术日趋成熟，对该类目标的跟踪引起了广泛的关注，主要的难点在于较低的飞行高度降低了雷达可探测的时间、所建立的跟踪模型难以囊括其复杂多变的机动模式[2-3]。

与跟踪其他机动目标相比，跟踪HGT存在一定的特殊性但也存在普遍性，相同的条件下量测传感器精度完全由所建立的目标运动模型和采用的滤波算法决定，因此如何实现高精度的HGT跟踪也是主要涉及这两方面。目标运动模型用于表征目标的真实运动，特别是机动目标，滤波算法基于目标运动模型从噪声量测中提取出目标的运动状态。常用的运动模型包括常速度(Constant Velovity, CV)模型(属于非机动模型)常加速度(Constant Acceleration, CA)模型、Singer模型、当前统计(Current Statistics, CS)模型、Jerk模型等，除了CV其他均属于机动模型[4]。由于HGT依靠气动力实施机动，因此气动力加速度对于跟踪系统来说是主要的未知输入项，对其建模的优劣一定程度上决定了目标运动模型的匹配性[5]。文献[2]将HGT气动力加速度中的未知气动系数项建模成维纳随机过程，并扩展至系统状态对其进行滤波估计，实现了稳定跟踪。文献[6]中采用与文献[2]中类似的建模方法，分别对HGT非跳跃式和跳跃式2种飞行轨迹进行了跟踪滤波，发现当非跳跃式轨迹中的倾侧角发生符号翻转时由于目标运动模型匹配不及时会引起较大的估计误差。文献[3,7]进一步考虑了各气动系数之间的耦合，给出了更加精确的气动力加速度模型，但是需要一定先验控制量的假设。上述模型均是从动力学的角度出发，对气动力加速度进行精细建模，称之为“动力学”模型。当然，从运动学的角度也可直接将HGT的加速度或气动力加速度建模成已有的机动模型，但是与“动力学”模型相比通常缺乏理论根据性和清晰的物理意义，跟踪精度对模型的参数较为敏感，需设置合适[6,8]。

由于HGT机动模式多样，利用单一固定模型的适应性较差，难以获得高精度的跟踪，甚至针对HGT的某种机动形式现有的机动模型均不能较好地匹配，例如HGT为了调整侧向运动会发生倾侧角符号的切换[9]。文献[10]基于非零均值正弦波相关模型对目标加速度进行均值补偿和方差自适应，提高了跟踪高超声速跳跃滑翔式目标时的模型适应性。文献[11]在认为控制变量服从一阶时滞过程的前提下，实现了气动参数模型与飞行状态的自适应匹配，仿真结果表明当HGT目标发生机动时所提模型的跟踪性能明显优于传统模型。文献[6]利用具有不同机动频率的CS模型构成的交互式多模型(Interacting Multiple Model, IMM)对HGT进行跟踪，并对目标气动力加速度方差进行自适应调整，估计精度优于传统的CS模型。文献[12]设计了一种自适应变结构的IMM算法用于跟踪HGT，提高了跟踪精度和效率。上述基于自适应模型的跟踪方法通常需要一定的假设条件，其普适性和准确性有待商榷，难以实用；而多模型方法中由于缺乏可靠先验信息使模型集个数和模型转移概率难以确定，模型集过少不能全面地表征目标的运动，模型集冗余会增加计算量、引起模型间的竞争，跟踪精度甚至得不到提升。此外，在跟踪HGT的仿真验证时，被跟踪轨迹的设计大都未结合HGT实际可能采用的轨迹设计方法且未考虑完成制导任务的需求，致使跟踪结果的可信度有所欠缺。

如前所述，目标跟踪精度和滤波算法也紧密相关，相同目标运动模型采用不同的滤波算法通常会产生不同的估计精度，在目标跟踪中自适应或鲁棒的滤波方法已经广泛地用于处理模型误差。针对系统模型中未知的噪声统计特性，提出了多种噪声估计器，其中基于极大后验准则的Sage-Husa估计器由于结构简单、原理清晰得到了大量关注[13-15]。但是，文献[16]指出Sage-Husa只有当过程噪声或量测噪声已知其中之一才是有效的。而且，引入估计器后易引起滤波的不收敛[17]。对于状态模型误差，基于正交性原理的强跟踪滤波(Strong Tracking Filtering, STF)表现出较强的鲁棒性，其通过渐消因子放大预测状态误差的协方差，强迫模型误差出现时正交性原理依然能够尽可能满足，进而调整增益矩阵来增加新观测数据的比重，以降低状态估计误差[18]。文献[19]将STF用于跟踪脉冲机动卫星，可快速地降低脉冲机动引起的估计误差，实现了高精度稳定跟踪。文献[20]将STF应用至高阶容积卡尔曼滤波(Kalman Filtering, KF)，增强了其跟踪目标时的鲁棒性。文献[21]提出了一种归一化的STF用于跟踪小脉冲机动的航天器。由此，针对HGT运动模型误差几乎无法避免的情况，也可配合鲁棒的滤波算法来提高跟踪精度。

众所周知，系统噪声的统计特性对卡尔曼滤波结果的影响较大。上述所提跟踪方法大多通过不同方式对系统过程噪声的协方差实施自适应调整，最终提高跟踪滤波精度。文献[22]在跟踪机动飞机时通过CS模型中的目标加速度递归方程和均值输入估计方法提出了一种自适应CS模型，实现了对过程噪声协方差的自适应调整。为解决跟踪机动目标时Singer模型参数的失配问题，文献[23]采用多模型方法中的似然函数并根据其变化自适应地调整最大加速度参数。对机动模型中的关键参数实施自适应调整的方法结构简单、原理清晰，对跟踪HGT具有一定借鉴意义。

为了提高HGT跟踪方法的适应性和精度，本文在现有研究的基础上提出了一种机动频率自适应的跟踪方法。首先将气动力加速度建模成形式简洁的Singer模型，基于正交性原理和无迹卡尔曼滤波算法利用滤波新息计算得到可反映状态模型误差大小的调整因子，然后用于放大Singer模型中的机动频率参数，实现对当前模型误差的自适应。仿真结果表明所提方法的适应性好、计算量小，能够对HGT 2种典型的机动飞行轨迹实现高精度跟踪。

1 目标跟踪模型

目标跟踪模型通常由状态方程和量测方程组成，所建模型与目标真实运动模式的匹配程度直接影响着跟踪精度。

1.1 状态方程

图1给出了HGT相对于地基雷达再入运动的示意图，oxyz表示雷达坐标系，也即东北天坐标系。考虑地球旋转，则HGT的相对运动方程可表示为

(1)

式中：δr/δt和δ2r/δt2分别为相对速度v和相对加速度；g和a分别为目标所受的引力加速度和气动力加速度；ω和r分别为地球自转角速度和目标地心矢径；2ω×(δr/δt)和ω×(ω×r)分别为科式加速度项和牵连加速度项。

图1 目标相对地基雷达的再入运动Fig.1 Target reentry motion relative to ground-based radar

为了获得系统的状态方程，需将式(1)投影至雷达坐标系下。设v=[vx,vy,vz]T，则相对加速度的各分量为

(2)

考虑椭圆地球，则引力加速度的形式为[24]

g=grr0+gωω0

(3)

式中：r0和ω0分别为r和ω方向上的单位矢量，其各分量分别为

(4)

其中：r为目标地心矢径r的模;px、py和pz为p的3个分量;Ro为Ro的模；Bo为地理纬度;μo=Bo-φo，φo为地心纬度。gr和gω的表达式分别为

(5)

式中：μ为地球引力常数;J2为带谐修正系数;aE为地球的长半轴。式(5)中，为了计算简便近似地用地基雷达纬度φo替代目标纬度φ。

科式和牵连加速度项的表达式分别为

(6)

(7)

式中：r=[rx,ry,rz]T;ω为ω的模。

气动力加速度作为HGT跟踪系统中主要的未知输入项，需进行重点考虑。这里直接从运动学的角度，将a建模成现有的机动模型。考虑HGT一直受到气动力的作用，且机动强度有所起伏，选用介于CV模型和CA模型之间的Singer模型来匹配a的真实变化，因此有[4,25]

(8)

式中：λx、λy和λz分别为3个方向上的机动频率；wx、wy和wz为互不相关的零均值高斯白噪声，功率谱密度分别为2λxσ2、 2λyσ2和2λzσ2，σ2为当前时刻目标加速度的方差，表达式为

(9)

式中：amax为目标加速度的最大幅值，目标以amax或-amax实施机动的概率均为Pmax，P0为目标加速度为零的概率。

针对上述运动学建模，系统状态方程可写

(10)

式中：X为系统状态向量，具体表达式为

X=[px，py，pz，vx，vy，vz，ax，ay，az]T

(11)

F(X)的表达式见附录A;w表示过程噪声，为互不相关的零均值白噪声向量，其自相关函数为

E[w(t1)wT(t2)]=Wδ(t1-t2)

(12)

式中：W为对角矩阵;δ(·)为狄拉克δ函数。

为了后续的滤波估计，需将式(10)进行离散化。设采样周期为T，则式(10)转化为

(13)

(14)

1.2 量测方程

地基雷达提供的量测量包括相对距离R、方位角A和高低角E，如图2所示。

图2 地基雷达量测量Fig.2 Ground-based rabar measurements

由图2，易得量测方程为

(15)

(16)

式中：p=[px，py，pz]T。

由于相对距离和角度之间的量级相差较大，为了降低滤波过程中出现矩阵条件数病态的可能性，将量测量转换至直角坐标系下，可表示为

(17)

式中:Hk为系数矩阵。

(18)

(19)

(20)

2 基于UKF算法的机动频率自适应

鉴于HGT机动模式多样，仅用单一固定的Singer模型难以表征气动力加速度的变化。在卡尔曼滤波过程中，新息反映着系统模型所存在的误差。假设量测方程准确的前提下，新息代表着状态模型误差，下面将其用于自适应Singer模型中的机动频率参数。起始将式(8)中的机动频率取至很小，此时Singer模型无限接近于CA模型，认为机动加速度变化平缓。当HGT机动加速度变化幅度增加时，放大机动频率参数调整过程噪声协方差。

由于状态方程是非线性的，因此只能选用非线性的滤波算法。相比于扩展卡尔曼滤波(Extended KF, EKF)，基于无迹变换的UKF省去了复杂雅克比矩阵的运算，且具有更高的理论精度[26]。因此，选用UKF作为滤波估计的框架。

2.1 机动频率自适应

对于线性KF，若系统模型能够完全匹配真实模型，则滤波器输出的新息序列满足正交性原理，即[26]

(21)

式中：ξk为新息向量，即量测量与其预测值的差。当式(21)满足时，意味着量测量中的有用信息被完全提取出来，因此新息序列的正交性可作为衡量滤波性能的标志。

对于次优的非线性滤波器，新息序列不可能是完全正交的。此时若是弱自相关的，即可认为其具有较好的滤波性能。以EKF为例，考虑如下形式的非线性系统：

(22)

(23)

式中：Kk为滤波增益；Φk-1和Γk的表达式为

(24)

因此，对于EKF使新息正交性满足的充分条件为

(25)

(26)

(27)

那么，量测量的预测值为

(28)

因此，ξk表示为

(29)

将式(29)和式(26)代入式(25)得

(30)

(31)

(32)

当状态模型存在误差时，会导致式(32)不再满足，左右两端不相等的程度反映着误差的大小，并采用左右两端主对线上元素的比值对其进行量化，具体形式为

(33)

式中：ζd称为调整因子，反映了第d个量测通道上误差大小。当目标进行大幅度机动引起状态模型误差时，可利用ζd对机动模型中的机动频率实施相应调整，以降低模型误差，提高滤波精度。当存在模型误差时，Vk的真实值在实际计算中是未知的，可以通过式(34)对其进行估算：

(34)

式中：κ为遗忘因子，通常设置为0.95。

针对上述所建立的HGT跟踪模型，量测量转换后三通道上的调整因子可分别用于调整x、y和z方向上的机动频率，x方向的调整策略为

(35)

2.2 基于UKF算法的自适应实现

由于式(33)是基于EKF推导得出的，不能直接嵌入至UKF框架下，因此采用如式(36)所示的等效计算式：

(36)

式中：Py为依赖式(22)所描述的模型计算得到的新息协方差矩阵的标称值，表达式为

(37)

对于式(22)所描述的非线性系统，首先对滤波器进行如式(38)的初始化：

(38)

自适应UKF算法可分为如下几个步骤[26]：

1) 计算simga点和权重系数

(39)

(40)

式中：β为与系统状态先验分布有关的常数，通常设为2。

2) 时间更新

(41)

(42)

3) 机动频率自适应

(43)

(44)

(45)

(46)

(47)

4) 量测更新

(48)

(49)

(50)

(51)

(52)

3 数值仿真

为了验证所提方法跟踪滤波的性能，首先设计2条典型HGT的飞行轨迹用于跟踪，并与其他几种跟踪模型进行对比。

3.1 HGT轨迹设计

以美国洛马公司发布的通用飞行器CAV-H模型为研究对象，分别基于阻力加速度-速度剖面跟踪和高斯伪谱优化方法生成2条再入飞行轨迹。考虑飞行过程中受多种约束的影响，主要设置最大过载为3g、最大动压为100 kPa和最大驻点热流密度为1 800 kW/m2，控制量攻角和倾侧角的最大变化率分别为20 (°)/s和85 (°)/s，初终端的状态约束如表1所示[12]。

表1 CAV-H再入初终端条件[12]Table 1 Initial and terminal reentry conditions for CAV-H[12]

跟踪阻力加速度-速度剖面是HGT再入段轨迹规划与制导中常用的一种方法，其将再入段分为初始下降段和滑翔段，以是否满足平衡滑翔条件为交班点，攻角剖面事先给定，初始下降段采用零倾侧角飞行，在滑翔段跟踪设计的阻力加速度-速度剖面用于控制纵向运动，由此确定倾侧角的幅值，侧向运动由视线方位角误差走廊进行控制，用于确定倾侧角符号的翻转时机。以纵向跳跃幅度最小为性能指标，基于高斯伪谱方法进行轨迹优化。具体两轨迹的形态和状态量的变化如图3～图5所示，“轨迹1”和“轨迹2”分别为基于跟踪剖面和高斯伪谱方法获得的轨迹。可看出，“轨迹1”存在倾侧角符号的翻转，会引起HGT的阶跃机动，当倾侧角不翻转时机动较平缓；“轨迹2”的倾侧角连续变化，机动幅度强于“轨迹1”中倾侧角符号不翻转时的幅度。上述2种典型轨迹可充分验证跟踪模型对不同机动模式的适应性。

图3 HGT跟踪轨迹Fig.3 HGT trajectories to be tracked

图4 控制量随时间变化Fig.4 Control variables versus time

图5 升阻加速度随时间变化Fig.5 Lift and drag accelerations versus time

3.2 跟踪滤波

为了充分地验证所提机动频率自适应方法的有效性，将其分别与基于Singer模型的方法、文献[6]中提出的基于CS模型的方法、文献[7]中提出的基于动力学模型的方法和文献[6]中提出的利用IMM实现机动频率自适应的方法进行性能对比。基于Singer模型的方法将气动力加速度建模成Singer模型，如式(8)所示，基于动力学模型的方法见附录B，基于CS模型的方法将气动力加速度建模成CS模型，形式为[6]

(53)

为了后续书写方便，将基于Singer模型的方法简记为“Singer”，在“Singer”基础上引入自适应机动频率的方法记为“Adaptive Singer”，基于动力学模型的方法记为“Dynamics”，基于CS模型的方法记为“Current Statistics”，图6所示的基于IMM的方法记为“IMM-Singer”。为了保持仿真验证的一致性，除“Adaptive Singer”中基于UKF算法对机动频率进行自适应调整外，其他方法也均采用UKF算法。每种方法进行300次的蒙特卡罗仿真，采用均方根误差(Root Mean Square Error, RMSE)来评估滤波结果的估计精度，位置RMSE的计算公式为

图6 交互式多模型算法结构Fig.6 Structure of interacting multiple model algorithm

(54)

3.2.1 轨迹1

将雷达部署在经纬高为(35°, 0.5°, 0)的位置，设其最大量测距离为800 km。考虑地球曲率影响，易得雷达可量测的轨迹弧段及相应状态量的变化，如图7～图9所示。由于图4中t≈600 s时，倾侧角发生了符号翻转，相应地在图9中引起了雷达坐标系下气动力加速度的突变。

图7 可观测弧段(轨迹1)Fig.7 Observable segment (trajectory 1)

图8 量测量真值(轨迹1)Fig.8 Noise-free measurements (trajectory 1)

图9 雷达坐标系下气动力加速度真值(轨迹1)Fig.9 True aerodynamic accelerations in radar frame (trajectory 1)

(55)

式中：Λij为从模型i转移到模型j的概率。

基于上述方法和仿真设置，分别进行300次的蒙特卡罗仿真，图10～图12给出了位置、速度和气动力加速度的均方根误差随时间的变化。可看出，当气动力加速度突变引起状态模型误差时，“Adaptive Singer”能够快速地降低估计误差，显著优于“Singer”，说明对自适应机动频率的有效性。在700 s以后，目标机动幅值有变化但较平缓，该方法表现出轻微的自适应过度，导致精度略低于其他几种方法，此时“Singer”和“Current Statistics”的精度稍高。“IMM-Singer”中由于采用IMM方法对机动频率进行自适应，因此其估计精度与“Adaptive Singer”最为接近。“Current Statistics”理论上也具有噪声方差的自适应能力，但是当加速度均值估计不准确时，可能适得其反，其结果劣于没有自适应能力的“Singer”。“Dynamics”对HGT阶跃机动的适应性最差。

图10 位置RMSE(轨迹1)Fig.10 RMSE of position (trajectory 1)

图11 速度RMSE(轨迹1)Fig.11 RMSE of velocity (trajectory 1)

图12 气动力加速度RMSE(轨迹1)Fig.12 RMSE of aerodynamic acceleration (trajectory 1)

3.2.2 轨迹2

将雷达部署在经纬高为(23°, 5°, 0)的位置，与“轨迹1”相同，易得雷达可量测的轨迹弧段及相应状态量的变化，如图13～图15所示。相比于轨迹1，轨迹2的机动是连续的，但是机动幅度的变化强于轨迹1无倾侧角翻转的时刻。

图13 可观测弧段(轨迹2)Fig.13 Observable segment (trajectory 2)

图14 量测量真值(轨迹2)Fig.14 Noise-free measurements (trajectory 2)

图15 雷达坐标系下气动力加速度真值(轨迹2)Fig.15 True aerodynamic accelerations in radar frame (trajectory 2)

为验证各方法的适应性，仿真设置与轨迹1中保持一致，图16～图18给出了300次蒙特卡罗仿真后位置、速度和气动力加速度的均方根误差随时间的变化。可看出，“Adaptive Singer”方法的滤波精度最高，对HGT机动幅度变化的适应性最好，特别是在目标机动幅度变化较大时，而“IMM-Singer”次之。由于机动幅度变化要强于轨迹一非阶跃机动的时刻，导致其他几种方法的估计精度较差，难以适应HGT多样的机动模式。

图16 位置RMSE(轨迹2)Fig.16 RMSE of position (trajectory 2)

图17 速度RMSE(轨迹2)Fig.17 RMSE of velocity (trajectory 2)

图18 气动力加速度RMSE(轨迹2)Fig.18 RMSE of aerodynamic acceleration (trajectory 2)

通过对不同轨迹的跟踪滤波表明，所设计的“Adaptive Singer”对HGT的阶跃机动和连续幅值变化的机动均具有较好的适应性。图19给出了各方法单次跟踪滤波的计算时间，可看出“IMM-Singer”约为其他方法的2倍，计算量大。相比于其他几种方法，“Adaptive Singer”不仅滤波精度高，而且计算量几乎没有增加，实时性好。

图19 单次跟踪滤波计算时间Fig.19 Computational time for single tracking filtering

4 结 论

为了精确跟踪机动模式多样的高超声速滑翔目标，基于正交性原理和UKF算法提出了一种机动频率自适应的跟踪方法。通过仿真验证，得到如下3点结论：

1) 由于HGT机动形式多样，几乎无法用固定的跟踪模型将其囊括，因此状态模型误差固然会存在。

2) 介于CV和CA模型之间的Singer模型可作为表征HGT机动的基础模型，再结合其他方法合理地调整其机动频率，降低模型误差。

3) 所提的“Adaptive Singer”方法能够较好地适应HGT的多种机动形式，与利用交互式多模型进行机动频率自适应的方法相比，跟踪精度高，计算量小。

(A1)

式中:F1=vx，F2=vy，F3=vz，F4、F5和F6的表达式分别为

(A2)

Rosinμo)sin2Bo-(pz+Rocosμo)·

sinBocosBo]-2ωvxsinBo+ay

(A3)

Rocosμo)cos2Bo-(py-Rosinμo)·

sinBocosBo]+2ωvxcosBo+az

(A4)

F7、F8和F9的表达式为

(A5)

附录B

若从动力学角度对HGT气动加速度进行建模，由于其采用倾侧转弯控制方式，则气动力加速度的形式为[1]

(B1)

式中：ρ为大气密度，近似为高度的指数函数;v为相对速度v的模;CVtoR为目标速度坐标系至雷达坐标系的转换矩阵;η为目标的面质比。将式(B1)中的右边部分分为可用目标位置和速度表示的量和不可表示的量，具体形式为[7,24]

(B2)

式中：M的表达式为

(B3)

(B4)

式中：γ为雷达坐标系至目标速度坐标系转换过程中的滚转角。通常情况下，HGV的面质比、气动系数和制导控制规律等先验信息均难以获知，因此广义气动系数是待估计的未知量。由于HGT在飞行过程中受到热流密度、动压、过载和控制量等因素的制约，其控制变量—攻角和倾侧角的变化应相对平缓，可用维纳随机过程来表征广义气动系数的变化，形式为[6]

(B5)

式中：wD、wLS和wLC为互不相关的零均值高斯白噪声。将式(11)中ax、ay和az分别替换为ηD、ηLS和ηLC即为动力学模型的状态向量，其状态方程的形式与式(10)一致。