临近空间高超声速目标跟踪动力学模型

李 凡，熊家军，张 凯，韩春耀，席秋实，张怀念

(1. 空军预警学院研究生大队，武汉 430019；2. 空军预警学院四系，武汉 430019)

0 引 言

临近空间高超声速飞行器(Near space hypersonic vehicle, NSHV)一般指飞行于20～100 km高空、飞行马赫数大于5的飞行器[1]，随着NSHV在军事领域发展及应用，其巨大的战略价值受到全世界的广泛关注，西方各军事大国投入大量资源进行NSHV关键技术的研究。2017年10月30日美国于夏威夷考艾岛导弹试验基地成功进行了潜射型高超声速助推滑翔导弹的首飞试验(试验代号CPS FE1)，此后，美国海军战略系统项目办公室(SSP)声称，已着手在上装载这一新型的垂发模块，“弗吉尼亚”级攻击型核潜艇将成为这一常规快速打击武器的发射平台，这一举措预示着美国NSHV已初步具备实战能力，有分析认为美国有望在未来几年启动NSHV装备采办程序。NSHV快速发展提高了进攻方的软杀伤及硬摧毁能力，但对于防御方而言，各国现有的预警监视体系对这类飞行器跟踪问题还没有很好的解决方案，因此，进行NSHV目标跟踪的关键技术研究具有重要意义。

机动目标跟踪问题实质上是对未知系统的描述，其单目标跟踪中主要包括滤波算法与机动模型两个方面[2]。滤波算法聚焦于局部量测噪声的误差修正，其本质是通过多个时刻残差样本提取量测噪声的分布特征参数，自适应修正当前时刻目标的状态估计值，滤波算法受噪声环境影响较大，对不同运动模式具有较强的适应性，此外，滤波算法发展依赖于基础数学理论的研究，难度相对较大；当前跟踪算法研究集中于机动模型改进或创新，针对不同机动样式目标跟踪的区别主要在于机动模型建立，机动模型致力于目标运动特性的描述，通过建立合理模型逼近目标真实运动模式，在观测数据有限的情况下，好的模型对于机动目标跟踪问题至关重要。当前高速强机动目标机动建模方法大体分为两条主线，一是从运动几何的角度研究目标状态随时间变化的规律[3-4]，直接以加速度(或加加速度)为基点，分析加速度之间的相关特性，并建立线性离散化的运动学模型跟踪状态方程，通过对加速度在时间维的积分求和得到目标状态估计，因此建模过程中不涉及坐标变换，该方法将目标看作一个等效质点，不考虑目标加速度产生的原因、目标与地球之间的引力关系及目标的姿态变化。常用的模型包括单模型及多模型，现有单模型主要有匀速(Constant velocity, CV)、匀加速(Constant acceleration, CA)、Singer、“当前”统计模型(Current statistical, CS)、Jerk[5-6]，这类模型结构简单，计算量小，但对复杂运动及强机动运动适应性较差；多模型主要包括固定结构多模型(Fixed structure multiple model, FSMM)和变结构多模型(Variable structure multiple model, VSMM)[7-8]，模型集覆盖的运动模式广，机动适应性较强，但噪声模型选取及模型交互参数设置较为困难。总体而言，运动学模型状态量之间具有很好的线性关系，通过状态转移矩阵递归实现状态预测，先验信息需求较弱，这一方法实质是对目标运动模式的趋近。二是从目标受力的角度出发[9]，对目标进行力学分析推导各方向的加速度特性，该方法分别在不同的坐标系下考虑重力、空气动力及推力(有动力的情况)的影响，由于目标受力分析分别在不同的坐标系下表达较为简洁，需要转换到同一坐标系下进行目标状态量估计及滤波，因此模型的非线性程度较高，此外，考虑力产生机动大多存在较为严重的耦合，该方法优点在于从弹道设计角度考虑目标机动的产生，模型建立较为合理时能较为准确地估计气动参数，但先验信息要求较强。现有NSHV目标跟踪动力学模型存在两个方面的问题，其一是动力学模型大多只考虑了空气动力因素即无动力情况，但NSHV存在有动力巡航方式，该方式飞行器控制性能、机动性能更佳，同时跟踪难度也更大。其二是动力学模型大多假设气动力参数服从为白噪声或一阶马尔科夫过程，需要建立更加合理的模型描述其气动加速度特性。

本文针对跳跃式NSHV目标提出了一种新型的动力学跟踪模型，从目标运动特性出发，对气动参数建立了新模型，使之同时体现周期性与衰减性，同时，为增强算法对机动突变点的适应性，引入强跟踪滤波进行加速度突变检测及滤波修正。

1 NSHV运动特性分析

平衡滑翔及跳跃式滑翔是NSHV最为常用的两种飞行方式，其中跳跃式滑翔能够在纵向平面进行多次跳跃，机动性能较强，是当前机动目标跟踪领域的一大难题，本文仅讨论跳跃式NSHV的跟踪。而对目标运动特性进行分析并抽象出其运动模式的一般规律，是构建合理跟踪模型的前提条件，因此，本节主要分析了NSHV在不同维度的运动特点以及可能的机动样式。

1.1 NSHV动力学分析

飞行器动力学分析有利于理解NSHV跳跃式飞行原理、分析其运动特性，当飞行器在空间内飞行时一般都有三个力作用于该物体上，即重力、推力、空气动力，即飞行器受力表达式为[10]：

F=FT+FA+FG=FT+FD+FL+FG

(1)

其中,FT为发动机推力(无动力时推力为0)，FA为空气动力主要分为空气阻力FD及升力FL，FG为重力。则NSHV在助推段的离散化动力学方程如下：

式中：k为时间刻度，ν(k)为目标速度，α(k)为攻角，M(k)为机体质量，g为重力加速度，θ(k)为航迹倾角，σ(k)为航向角，r(k)为目标质心到地心距离，φ(k)为目标纬度。

在助推段分析推力及空气动力对加速度dν(k)/dk及纵向平面角速度dθ(k)/dk的影响程度，分别假定式(2a)及式(2b)中除推力及空气动力外参数为常数，易知在FT(k)cosα(k)与空气阻力FD(k)在加速度dν(k)/dk中的灵敏性相同，同理，FT(k)sinα(k)与空气升力FL(k)在dθ(k)/dk中的灵敏度一致，此外，推力在加速度及纵向平面角速度中分别为正弦和余弦对偶分量，因此，推力FT(k)至少有一个分量与空气动力分量在助推段的灵敏性相同。

NSHV目标助推拉升后进入无动力滑翔阶段，这一方式能有效规避雷达探测，降低视距。其中无动力滑翔阶段动力学方程如下：

(3)

式中：φ(k)为倾侧角，当目标无动力飞行时受力较为简单，只存在空气动力与重力，飞行器通过自身调整实现横向机动。目标航向角σ(k)较小[11](σ(k)∈[0°,10°])，即目标在横向机动强度较小；根据文献[11]可知航迹倾角变化率dθ(k)/dk较大，目标在高度层机动强度较大。

当有动力与无动力互相转换时相当于加入一个阶跃的推力，根据分析可知至少有一个推力分量与空气动力对运动影响程度处于同一量级，仅考虑无动力情况时，动力学模型难以描述转换阶段的加速度突变，转换阶段容易出现跟踪误差升高。

1.2 横纵向机动特性分析

纵向平面内“打水漂”跳跃飞行是这类目标最鲜明的特点，其轨迹跳跃幅度呈衰减形式，类似于振荡运动，可以推测相邻整数周期的目标状态量呈现类周期性变化[12-15]。

从飞行器控制参量方面看，攻角与倾侧角控制飞行器的主控量，其中攻角控制机体纵轴与速度方向间的俯仰幅度，通过不同的升力/阻力系数决定目标在纵向上的飞行高度、飞行速度、航迹倾角等，攻角的设定需要满足飞行走廊的限制条件，不可能非常复杂，攻角经常采用常值、分段或线性等简单函数来描述。倾侧角同时影响目标在横纵向的飞行轨迹，决定了飞行器的横向机动，为了有效规避、突防禁飞区，倾侧角设定通常是在飞行走廊内幅度不变符号反转。

攻角方面。文献[16]指出在满足飞行走廊条件限制前提下无论如何设定攻角控制量，只要不是一直满足平衡滑翔条件，其弹道都呈现出跳跃滑翔样式。实际上平衡滑翔比跳跃滑翔更难实现，平衡滑翔方式对控制参量设定要求更高，一般来说，跳跃滑翔攻角取值可以是满足飞行走廊条件下的一个区间，而平衡滑翔攻角取值只能是满足式(4)条件的一个定值。

(4)

平衡滑翔状态下，航迹倾角很小，一般可以认为cosθ(k)=1。其中

FL(k)=CL(k)q(k)S

(5)

q(k)=ρ(k)v2(k)/2=ρ0e-βh(k)v2(k)/2

(6)

式中:CL(k)为升力系数，是与控制量攻角α(k)和马赫数Ma相关的函数，q(k)为动压，S为飞行器的参考面积，ρ(k)为当地大气密度，ρ0为海平面大气密度，β为高度常数，h(k)为目标高度。

文献[16]指出当其他参数一致仅改变初始高度时，若飞行器滑翔初始高度不满足平衡滑翔状态时，NSHV轨迹都会产生跳跃。当初始高度h2(k)小于平衡滑翔高度h1(k)时，当地大气密度ρ0e-βh2(k)大于平衡滑翔的ρ0e-βh1(k)，导致飞行器升力增加，无法满足平衡滑翔条件，其合力作用向上，NSHV向上做加速运动，当地大气密度降低，升力减小直到合力为0，向上的速度达到最大，然后向上做减速运动，爬升到最高点后，在重力的作用下向下做加速运动，高度降低，升力增加并再次达到平衡时，向下的速度最大，高度继续降低直至向下速度为0，飞行器到达最低点，此后继续向上飞行，完成一个轨迹周期。(初始高度大于平衡滑翔高度也会形成跳跃滑翔轨迹)。此外，仅改变初始速度或初始速度倾角时，其弹道同样为跳跃滑翔样式。因此，NSHV升力纵向的升力加速度相关性也具有类似的周期特性。

倾侧角方面。倾侧角的设定主要是规避敌方探测拦截，美国最初将NSHV目标定义为全球快速打击武器，旨在高速对敌方大纵深目标进行软杀伤或硬摧毁，可以通过大C或S形横向弹道规避敌方重点防御区域后依靠高速打击纵深目标，横向高频跳跃或其他方式复杂机动会导致在敌方区域的无效路径增长，被发现的概率增大，攻击的时效性降低。此外，NSHV飞行速度较快，横向机动半径较大，飞行

器动能有限，无法如纵向进行多次跳跃。

2 动力学跟踪模型构建

假定雷达站坐标为东-北-天坐标系(East-north-up, ENU)，在半速度(Velocity-turn-climb, VTC)中分析气动力模型及推力模型，并在地心惯性坐标系(Earth central inertial, ECI)中构建重力及表视力模型[10]。

(7)

2.1 重力及表视力模型

假设目标在ECI及ENU坐标系下位置状态分别为XECI,XENU，假定雷达站的经纬高为L,B,H，可知ECI到ENU的转换为[5]：

(8)

其中转换矩阵为

Re为地球等效半径，ωet为从参考时刻到t时刻ECI和地心坐标系(Earth central, EC)的夹角(一般假定雷达探测到目标开始EC和ECI坐标系重合)[17]，则ENU坐标系下目标加速度为：

(9)

2.2 气动力模型

推力产生加速度是通过对空气作用而实现的，因此，将推力模型认为是一种叠加的特殊空气动力，则aP=aT+aA可表示为：

(10)

根据VTC坐标系与ENU坐标系的转换规则，ENU坐标系下新的气动加速度为[17]：

(11)

2.3 气动加速度建模

目标在飞行过程中重力加速度变化较小，其运动模式变化主要取决于气动加速度的改变，第2.2节中分析了各项加速度的物理含义，结合第1.2节中NSHV机动特性分析，对阻力加速度，横向转弯加速度及爬升加速度进行分布假设建模。

1)对于爬升加速度，由纵向机动特性可知，其自相关存在周期性与衰减性，因此，采用阻尼函数表示爬升加速度的相关特性。其自相关函数为：

(12)

1/[α+(β-w)j]+1/[α-(β+w)j]}

(13)

其中，F(·)为傅里叶变换，w为角速率，令G(w)=[α-(β-w)j][α+(β+w)j]，则有

W(wj)H(wj)H(-wj)

(14)

其中，W(wj)为白噪声的傅里叶变换。则滤波器的传递函数为：

(15)

其中，s为拉普拉斯算子，则机动加速度a(t)的连续时间微分方程为：

(16)

基于动力学模型跟踪NSHV，需要将VTC坐标系下的气动加速度进行转换到ENU坐标系下，这一过程存在高度的转换耦合，因此难以得到过程噪声协方差的显性表达式，但状态协方差的迭代更新可以一定程度适应过程噪声协方差取值变化，所以更为关注状态转移矩阵的离散化，式(16)可简化为：

(17)

2)对于转弯力加速度，根据横向机动特性分析可知，其机动强度相对纵向较弱，常用的弹道形状为直线、C形或S形，其横向加速度相关性呈较强的衰减性，因此假设横向加速度服从一阶马尔科夫模型，转弯加速度的微分方程为：

(18)

式中：γ为机动时间常数。

3)对于阻力加速度，当目标质量、等效面积等恒定时，目标所受阻力仅与速度、空气密度有关。

(19)

根据此前分析可知，跳跃滑翔式目标在纵向上的高度呈类周期的跳跃形式，而目标总体速度呈阶梯下降的趋势，因此，其阻力加速度也呈现振荡下降的形式。因此，将其建模为衰减振荡自相关随机过程。

3 动力学跟踪模型状态方程

为方便讨论，在此区分运动状态量与气动状态量的概念，运动状态量为ENU坐标系下状态量(包括位置、速度、加速度)，气动状态量为VTC坐标系下状态量(包括加速度及加加速度)。

3.1 运动状态量更新方程

k时刻运动状态量更新方程fI(X(k))为[17]：

(20)

3.2 气动状态量更新方程

式(18)的离散化方程为：

at(k+1)=e-γTat(k)+wt(k+1)

(21)

其中，T为采样时间间隔。式(17)的离散状态方程可通过式(22)求得。

(22)

(23)

式中:

则完整的跟踪状态离散方程为：

4 加速度突变点的滤波修正

根据式(2)～(3)可知，有无动力段转换时目标产生较大的机动加速度，将导致跟踪误差增大，在跟踪初始模型经过一段时间的调整与目标运动模式匹配，之后较为稳定地描述目标运动，并通过增益进行预测量的修正，模型需要保证跟踪算法整体的鲁棒性，其加速度不会出现突变，当目标出现机动突变时模型无法进行瞬态变化，同时，滤波算法修正误差需要一定的量测样本，因此会出现机动跳变位置的误差升高。针对这一问题引入强跟踪滤波对机动突变点进行补偿修正。考虑到有动力转换主要用于推动飞行器从弹道低点拉升，选用升力加速度ac(k)作为机动突变检测参数。判决量：

(24)

当相邻加速度相差达到两个数量级以上即ξ(k)≥2时，认为发生机动突变，采用强跟踪滤波STF对当前时刻量测进行滤波修正[18]。

5 仿真分析

图1 仿真轨迹

5.1 仿真试验1

航迹I为有动力巡航轨迹，航迹I初始位于助推段从雷达视距中出现，达到助推顶点后经过一个完整的滑翔—助推，完成跳跃巡航。三种算法的跟踪仿真结果如图2所示。

图2 位置及速度误差均方根比较

如图2所示，算法I直接通过运动加速度描述目标的位置及速度变化，在滤波初始阶段,目标运动加速度变化的量级远超过气动加速度，算法I需要一定的时间适应加速度变化速度，其位置误差及速度误差收敛速度较慢；在200～300 s内由于运动加速度出现突变，导致位置及速度滤波误差都有一定程度跃起，其机动适应性相对较差。其余时刻目标机动较弱，运动模式趋于匀速及匀加速之间，三种算法位置滤波精度相近，但该算法速度滤波精度较高，略优于算法II，与算法III相当。

算法II基于气动加速度扩维实现目标状态预测及滤波，气动加速度变化的量级小，但本身的量级也较小(参见图3中气动加速度量级)，因此，需要一段时间才能稳定描述气动加速度变化，这是限制动力学模型收敛时间的关键性因素之一；在200～300 s内强机动时其位置及速度误差相对算法I较小，动力学模型对运动加速度的突变具有一定的适应性；其余时刻起误差水平较为平缓。

算法III基于衰减振荡函数构建新型动力学跟踪模型，并通过强跟踪滤波检测气动加速度突变点实时修正滤波状态，其收敛速度较快；在200～300 s内速度滤波误差起伏较小，但位置误差保持平稳趋势。整体而言本文算法整体水平较为稳定，滤波效果最好。为更好地比较算法性能，给出各算法气动加速度滤波估计仿真图。(其中算法I的气动加速度估计是通过气动加速度到运动加速度逆运算实现)。

图3 气动加速度参数估计

在气动加速度方面，算法I通过滤波状态解算得到气动加速度量，如图3所示，算法I的气动加速度估计最大特点为过调整幅度较大，因运动加速度变化量较大导致调整速度较快、过调整较为严重，此外，在200～300 s内，由于运动加速度在这一时间段内变化较大，该算法解算的气动加速度在此时间段内出现局部的大幅度跳跃，但该时间段内目标实际气动加速度并没有剧烈变化。

算法II气动加速度量的估计明显优于算法I，过调整量相对较小，但调整时延较长，在200～300 s内，能较为准确反应真实气动加速度的变化趋势。

算法III采用了新型模型，对气动加速度量估计效果最好，其包络与真实值匹配程度最高，此外，采用气动加速度检测及强跟踪滤波补偿，其时延相对较小。

图4 真实运动加速度

如图3所示，气动加速度方面在0～100 s内快速变化，其余时刻曲线较为平缓；而运动加速度方面如图4所示，出现两次变化较大，造成两种加速度不对应的原因分析如下：

1)气动加速度是运动加速度的分量，由式(20)可知，气动加速度为运动加速度分量，运动加速度还包括重力及表视力，特别是重力不可忽略。

2)气动加速度的耦合，由式(11)转换矩阵可知，气动加速度分量对运动加速度的影响与当前各方向的速度有关，气动加速度与运动加速度的转换存在严重的耦合，从图3及图4也可以看出,气动加速度与运动加速度值不在一个数量级，其对应的变化规律难以定性分析。

3)机动的坐标属性，目标运动都是具有坐标属性的，在不同的坐标系下运动模式会呈现不同的形式，本文仿真的雷达观测位置也一定程度造成了气动加速度与运动加速度的不对应关系。

因此，三种算法位置及速度误差在200～300 s段出现一定的误差起伏。无论动力学或运动学跟踪模型最终都必须过渡到运动加速度、速度进行状态量估计，运动加速度的突变一定会引起速度及位置急剧变化，这也是气动加速度平缓而出现滤波误差突变的原因。

5.2 仿真试验2

航迹II为无动力滑翔轨迹，跟踪算法的参数设置不变，得到算法I、算法II及算法III的滤波结果如下图(篇幅所限，仅给出气动加速度估计结果)。

图5 气动加速度估计

无动力滑翔时，气动加速度调整量较小，三种算法都能大致反应气动加速度变化，并呈现局部小范围的波动，在阻力加速度估计中，算法I及算法II都出现较大波动，算法I容易过调整，估计精度一直处于波动状态；算法II前期误差较大，但后期估计精度较为稳定；算法III整体水平最为稳定，精度较高。在转弯力加速度估计中，算法I同样存在过调整现象；算法II调整速度较慢，试验较长；算法III效果较好。对于爬升力加速度估计，前期三种算法估计误差较大，后期误差趋于稳定，且三种算法精度相近。

为定量比较算法的滤波性能，位置、速度、气动加速度误差统计平均如表1所示。

综上所述两个部分的仿真结果，算法I适应匀速到匀加速之间的运动，直接通过反应运动加速度变化，导致对气动加速度估计容易出现过调整。算法II整体滤波误差较为稳定，气动加速度估计略优，但由于气动加速度量级较小，导致位置误差收敛速度较慢。将本文算法算法III与算法I及算法II进行比较，可知本文算法所构建模型的优越性。本文主要将气动加速度建模为衰减振荡函数形式，对气动加速度的变化规律匹配程度更高，在位置及速度跟踪误差方面精度更高；此外，引入强跟踪滤波的机动检测与修正，降低估计时延提高跟踪精度。

表1 位置、速度及气动加速度误差统计平均Tabble 1 Position, velocity and aerodynamic acceleration error statistical average

本文的研究重点在于将气动加速度自相关假设为具有衰减振荡的二阶马尔科夫模型，因此，在仿真中仅对比了两种经典的指数衰减一阶马尔科夫模型，在本文的仿真场景下本文模型跟踪性能略优，但并不能说明其效果一定优于改进后的经典模型，所幸本文模型从底层的建模假设进行研究，适用于经典模型的改进方法同样可以运用于本文模型。此外，当前IMM依然是NSHV跟踪的首选方案，本文模型为后续跟踪NSHV目标IMM中模型选择，提供了一个更具潜力的单模型。

6 结 论

本文基于滑跃式NSHV运动方程，结合NSHV控制参量分析了目标的横纵向机动特性，阐述了滑跃式轨迹产生的原因，在此基础上构建了一种新型NSHV跟踪动力学模型，仿真表明本文模型相对现有经典模型有一定优越性，但还存在衰减系数、机动频率设置，以及在IMM中的应用等其他方面问题有待进一步研究。