1. 引言
随着世界航空航天技术的飞速发展以及临近空间环境认知水平的不断提高,各大军事强国逐渐发现了这片空域的潜在军事价值,加速了对于临近空间的开发利用 [1]。高超声速飞行器是指以吸气式发动机或其组合发动机为主要动力的在临近空间中以大于5马赫(Ma)速度飞行的空天飞机、高超声速轰炸机、高超声速巡航导弹以及高超声速再入滑翔导弹,具有超高的飞行速度,良好的机动性能,较强的突防能力,广阔的作战空间,具备快速精确打击、远程快速投送等能力,将成为未来完成反介入作战、精确打击时间敏感目标和重大战略节点的新式武器。
高超声速武器的发展,给各国现有防空预警系统带来了很大的挑战,相信在可预见的未来,此类飞行器的潜在威胁将会逐渐转为现实,出于保护国家国防安全的需要,必须发展并加强相应的防御手段。因而,分析低轨预警卫星系统对于临近空间高超声速飞行器的定位能力,为临近空间高超声速飞行器的预警体系研究、论证、建设提供参考具有十分重要的意义。另外,高超声速飞行器飞行速度快,严重压缩预警系统反应时间,对中/近程防空武器的拦截能力提出了较高的要求 [2]。作为防御方,研究临近空间高超声速飞行器的轨迹预测问题,一方面有助于预警探测网络高效地管理探测资源,保持连续跟踪目标;另一方面有助于制定拦截作战方案,对于防御该类目标具有重要意义 [3]。本文运用理论分析与仿真实验相结合的方式以及先定位跟踪后轨迹预测的次序对于临近空间高超声速飞行器展开研究。
2. 高超声速飞行器运动方程
忽略地球自转与扁率,在地球惯性球坐标系下,建立目标飞行器的运动方程。本文主要针对的是周期性高超声速巡航弹道,该飞行器以类似幅度逐渐减小的正弦曲线的运动方式进行跳跃滑翔,是飞行器中一种典型的运动特性,具有重要的研究价值以及广阔的运用前景。
高超声速飞行器跳跃滑翔运动方程:
(1)
其中:V为目标飞行器的飞行速率;分别为目标所在的经纬度;r为飞行器与地球的质心距;表示高超声速滑翔飞行器的空间三维位置;s为飞行器射程;分别为弹道倾角和弹道偏角。分别为大气密度、重力加速度、阻力、升力以及飞行器质量,可计算为:
(2)
式中:是海平面标准大气压为1.225 kg/m3,h为飞行器高度,为引力3.978655 × 105 km3/s2,地球半径为6371 km,为升力系数,为阻力系数,S为机翼参考面积。
以美国波音公司设计的通用再入飞行器(CAV-L)的气动数据为例计算气动系数,当攻角固定,马赫数大于8时,气动系数随马赫数变化缓慢,因此气动系数可简化为攻角的函数。可采用如下经验公式对高超声速滑翔飞行器进行气动系数的计算。
(3)
其中,
飞行器质量、机翼参考面积以及攻角借鉴CAV总体参数表,取质量m为2200 kg,参考面积S为12.66 m2,攻角范围为−4˚~20˚。
3. 双星定位算法以及仿真实验
3.1. 构建双星观测系统
双星观测系统由两颗观测卫星组成,每颗卫星上都载有目标观测平台,通过求得每颗卫星在空间中的位置以及卫星与目标间连线的视线矢量,从而获取目标的空间位置信息。定位过程如图1所示。
Figure 1. Double star positioning diagram
图1. 双星定位示意图
首先建立双星对于目标的观测坐标系,如图2所示。
Figure 2. The definition of satellite instantaneous orbit coordinate system in geocentric inertial coordinate system
图2. 地心惯性坐标系中定义卫星瞬时轨道坐标系
其中将双星上的瞬时轨道坐标系设定为对于目标的观测坐标系,其中卫星质心是该坐标系的原点,是由卫星质心指向地球质心(又称为当地垂线),轴指向卫星运行轨道的负法向,与其他两轴组成右手坐标系,卫星角度测量如图3所示。
定义各轴测角信息:
1) 方位角e,卫星的视线矢量在面内的投影与坐标轴的夹角,由指向为正,其取值范围为。
2) 俯仰角a,卫星视线矢量在面内的投影与视线的夹角,由指向为正,其取值范围为。
Figure 3. Schematic diagram of satellite angle measurement
图3. 卫星角度测量示意图
3.2. 目标几何定位算法
目前常用的无源定位算法有3种,依次是几何定位算法、基于奇异值分解的最小二乘算法、总体最小二乘算法。对于本文的双星观测体系来说,几何定位算法不仅可以得到理论上的最优解,并且计算速度最快,定位精度最高。故本文采用几何定位算法得出目标的定位模型。
在卫星观测平台获得目标方位角e、俯仰角a信息,目标在瞬时轨道坐标系下的视线单位矢量:
(4)
观测矢量:,,为瞬时轨道坐标系到坐标系的转换矩阵。则观测矢量有下式:
(5)
式中,为绕三轴的旋转矩阵,轨道参数为地心惯性坐标系到卫星轨道坐标系的轴转动角。
(6)
对式(8)取逆过程得到观测矢量模型为:
(7)
(8)
双星定位目标过程如图4所示:图中的C与D分别表示两颗观测卫星在地惯系下的空间位置,CA与DB分别表示两颗卫星的观测视线,其单位矢量即为观测矢量,AB为两视线的公垂线。本文中利用双星系统实现对于临近空间高超声速飞行器的定位,构建目标定位模型,基于双星交汇定位的几何算法进行计算,理论上两条直线交点就是目标点,但在实际观测过程中,由于测量误差的存在,观测系统中两颗卫星的观测视线可能存在一定偏差,得到的两条直线并非同面相交,导致两条观测视线在空间中并不能交于一点。这时,可以采用两条视线的唯一公垂线与两条视线交点的中点作为目标的实际空间位置,使得定位误差最小。
设,可得:
(9)
(10)
(11)
根据空间异面直线定理得:
(12)
(13)
则空间目标在地惯系下的位置为:
(14)
3.3. 仿真实验
3.3.1. 飞行器弹道仿真
本文基于飞行器在大气层内进行无动力跳跃滑翔飞行过程中速度倾角、高度与速度之间的关系,如下(15)式,对高超声速飞行器的初始状态进行设定,并利用四阶龙格–库塔法对目标弹道进行仿真分析。设置目标的初始高度为60 km,初始经纬度均为0˚,目标初始速度为6.2 km/s,同时目标初始弹道倾角和弹道偏角分别为0.05 rad、0.1 rad,选取飞行器运动规律特征明显的跳跃滑翔段主弹道,即设定1200 s的飞行时间。经过四阶龙格–库塔迭代后的仿真弹道如下图5~8所示。
(15)
Figure 5. Hypersonic vehicle altitude-time trajectory
图5. 飞行器高度–时间轨迹图
Figure 6. Hypersonic vehicle range-altitude trajectory
图6. 飞行器航程–高度轨迹图
Figure 7. Three dimensional trajectory of hypersonic vehicle
图7. 飞行器三维运动轨迹图
Figure 8. Hypersonic vehicle speed-time diagram
图8. 飞行器速度–时间图
从以上仿真结果可以看出(图5~8),目标飞行器跳跃滑翔轨道比较稳定,高速飞行器速度均大于5马赫数即大于1.7 km/s,滑翔高度稳定在20~100 km之间,满足超高声速分析器的技术指标要求。
3.3.2. STK目标观测场景的搭建
高超声速飞行器的经纬度高度三维轨迹如图7所示。通过经纬度和立体三维坐标之间的换算关系,将目标飞行器的空间三维信息转换成xyz三维坐标信息,并做成星历文件导入STK中,如图9。
利用STK可以方便地建立Walker星座,设置轨道平面数量为3,每个轨道平面中卫星数量为8,相邻轨道卫星相位为1,采用极地轨道,将星座高度设置为1000 km。基于此观测场景和STK中给出的相应数据进行目标的定位分析,最后点击应用后生成Walker星座。
Figure 9. Target ephemeris file import
图9. 目标星历文件导入
从STK仿真双星定位目标飞行器的过程中,选择观测时间段为4:00:15~04:14:37,即飞行器在空间飞行的时间共计862 s。取S111和S123两颗卫星共同可见的7分零5秒,即4:00:15~4:07:20的观测时长对目标进行定位性能分析。双星定位过程如图10所示。
3.3.3. 仿真分析
根据STK观测系统给定的误差参数范围以及观测星上跟踪设备自身的测角误差,设置观测卫星自身的位置误差为18 m (3σ),速度误差为0.3 m/s (3σ),观测角度误差为0.5e−3rad (3σ)。为了更加精确地得到并分析目标的定位误差,这里采用误差计算中的均方根误差(RMSE)来进行计算,其中RMSE的计算公式如下:
(16)
其中N是蒙特卡洛打靶仿真次数,是目标的真实位置,是定位算法得到的目标位置,经过Matlab仿真,可以得到目标定位精度的仿真图,如图11所示。
Figure 11. Root mean square error of target positioning
图11. 目标定位的均方根误差
从图11看出,双星对于目标的定位精度在800 m左右,其中x轴定位精度最高,误差在200 m左右,y轴定位误差较大,误差在900 m左右,z轴定位误差在400 m左右。分析可得,双星定位算法对于目标的定位精度较好,误差基本上控制在公里量级以下。影响观测系统最终定位精度的误差源种类繁多,双星观测体系中的误差源主要包括: 观测系统静态指向误差、动态稳定误差、观测系统与卫星间的对准误差、卫星的测量误差以及定位算法误差等。后续研究可针对减小以上误差,进一步提高卫星定位精度。
4. 高超声速飞行器的轨迹跟踪与预测
针对高超声速飞行器机动能力强、运动特征复杂以及目标先验知识掌握情况不同等问题,目前一般有基于动力学模型和运动学模型两种轨迹预测方法。前者通过分析高超声速飞行器的受力情况,对目标进行动力学建模或进行特征参数拟合,利用飞行约束条件进行弹道反设计,获得运动轨迹的解析形式,但该类方法对没有足够目标先验信息的防御方而言实现难度较大。而后者利用运动学模型拟合高超声速飞行器的运动特征,对目标进行持续跟踪观测,该方法对目标的先验信息需求较小。所以本文采用第二种方法,针对非弹道式高超声速飞行器周期跳跃的运动的特点,应用成熟的运动学模型描述目标的运动模型,结合改进的滤波方法对目标运动状态进行估计,通过函数逼近的方法对目标航迹拟合得到目标模型,最后将滤波估计状态值代入双正弦和函数拟合的速度曲线中,进一步递推目标运动轨迹 [4]。
4.1. Singer运动模型
Singer模型又称为时间相关模型,是一种机动模型,其将目标的机动认为是与时间相关的有色噪声序列,并不是如匀速运动(CV)模型和匀加速运动(CA)模型那样将机动认为是高斯白噪声序列。Singer模型认为目标的加速度是零均值的一阶时间相关的稳态马尔科夫过程,则加速度满足时间相关函数:
(17)
式中,分别为区间内决定目标机动特性的待定参数,为机动加速度,为机动加速度方差,为机动时间常数的倒数,即机动频率,通常其经验取值范围为:转弯机动;逃避机动;大气扰动。
对时间相关函数应用Wiener-Kolmogorov白程序化处理后,机动加速度可用输入为白噪声的一阶时间相关模型来表示为:
(18)
式中,是零均值高斯白噪声,方差为。
则目标机动模型可以描述为如下的一阶时间相关模型:
(19)
将上式离散化,其状态—空间模型可以表示为:
(20)
当时,时,此时目标做匀加速运动,时间相关模型趋近于CA模型。当时,,目标做匀速直线运动,时间相关模型可简化为CV模型。当时,时间相关模型介于匀速和匀加速之间。因此,时间相关模型有较大的目标机动覆盖范围,是更有效的目标机动模型 [5]。
一般情况下一阶时间相关模型就可以很好地描述目标的运动状态,但当目标做持续长时间的剧烈机动时,该模型具有较大误差,必须建立如下的二阶时间相关模型:
(21)
其中为单位强度白噪声,参数。
4.2. 无迹卡尔曼算法
无迹卡尔曼滤波(UKF)是基于线性回归的思想,通过UT变换技术构造Sigma点,基于Kalman的基本框架来进行状态和协方差的更新。UKF是针对状态概率密度函数进行近似,与扩展卡尔曼(EKF)对非线性模型进行近似不同,即便模型比较复杂,也不增加算法的实现难度,而且其不受限于系统的形式,不需要计算雅可比矩阵,可以处理不可导非线性函数。其中UT变换所得到的非线性函数的统计量的准确性可以达到三阶Taylor级数相当或更高,UKF的采用形式为确定性采样,Sigma点与状态的维数有关。
4.2.1. UKF非线性状态方程
UKF离散非线性状态方程如下 [6]:
(22)
其中状态向量为维向量,观测量为维向量,模型的过程噪声是均值为零的高斯白噪声向量,其方差阵为对角矩阵,观测噪声也是均值为零的高斯白噪声向量,其方差阵也是对角阵矩阵。
4.2.2. 量测方程
系统的量测信息是基于双星观测系统探测得到的,目标的位置信息为,通过转换可得质心距、方位角和俯仰角的量测信息。
建立高超声速目标量测方程如下:
(23)
其中,为状态向量中位置信息部分,也是均值为零的高斯白噪声。
4.2.3. UKF滤波算法步骤(图12)
下面将从4个方面介绍UKF滤波方法的核心思想:
1) 初始化
状态向量初始估计值:
原状态向量协方差初始估计值:
2) 计算采样点
UT变换构造2n + 1个对称的Sigma采样点如下:
(24)
其中和分别是当前时刻的状态估计值和协方差矩阵,n为状态量的维数,,决定采样点距离均值的远近程度,通常取值比较小;保证方差阵的半正定性。
3) 时间更新方程
系统状态更新方程和一步预测协方差更新方程为:
(25)
计算权值为:
(26)
4) 量测更新方程
系统测量预测协方差为:
(27)
权值:
(28)
系统增益矩阵为:
(29)
系统状态估计及其协方差矩阵分别为:
(30)
Figure 12. Sketch map of Unscented Kalman filter algorithm based on Singer model
图12. 基于Singer模型的无迹卡尔曼滤波算法示意图
4.3. 基于双因子加强的自适应Singer-UKF算法
一方面本文将基于双星观测系统而得到的目标定位信息作为滤波的状态初值与先验信息,由于双星定位算法本身存在误差,所以随着误差的引入,必然会影响后面跟踪滤波算法的精度,并且量测模型中将三维坐标间接计算转换为量测量,使得量测噪声不准确。另一方面很难存在某一种模型能将目标的整个运动状态描述的非常完美,若运动模型与实际目标运动状态失配较严重就会导致滤波效果发散,进而增大跟踪滤波算法误差。基于以上问题,本文引入双调节因子,通过选择合适的调节因子自适应调节状态预测值与观测值在滤波算法中的所占的比重,进而减小因滤波初值不准确以及状态方程偏差对滤波的影响,同时降低因异常扰动引起的滤波跟踪误差。而引入合适的调节因子能增强滤波器的跟踪能力,使输出残差矩阵一直保持相互正交,并且能实时调节状态协方差等项。
4.3.1. 构造调节因子
令
(31)
从上式中可以看出,,将调节因子代入UKF滤波算法中,
(32)
由上式可以看出,当状态滤波初值以及状态方程存在误差时,改进后的UKF算法通过加入调节因子,将、以及自适应膨胀,减小状态方程对滤波结果的影响,实时平衡初值和状态方程权重,从而减小滤波初始误差 [7]。
4.3.2. 构造调节因子
正交性原理是使得滤波器加强跟踪滤波性能的理论依据,因为增益是实时变化的,通过引入调节因子不断调整系统增益,使得输出残差呈正交性。
系统状态估计为:
(33)
令残差,则上式可写为。
通过调整,并由正交性原理可得:
(34)
构造调节因子如下:
(35)
其中:
(36)
在上式中,分别为系统噪声协方差阵和测量噪声方差矩阵,为弱化因子,,为状态转移矩阵,为系统协方差矩阵,为残差协方差矩阵,如下式:
(37)
其中,为遗忘因子,且。
将调节因子代入UKF滤波算法中,
(38)
4.3.3. 实时调整系统噪声
在实际应用中,系统噪声协方差阵在非线性系统往往是时变且不确定的,可根据极大后验 (Maximum a posterior, MAP)估计原理 [8] 以及指数加权的方法,设计出次优无偏MAP常值噪声统计估计器,并推导出时变噪声统计估计器的递推公式。通过在滤波器中加入调节因子,让系统噪声能实时的调整。
在渐消记忆时变噪声统计估计器中,有如下递推公式 [9]:
(39)
其中,b为遗忘因子,取值范围为,通过多次仿真可知,b取0.89时能得到
较好的跟踪滤波效果。
将选择好的b值代入递推公式中,可以得到如下:
(40)
将加入到UKF算法中就可以使其进行自适应滤波,提高了滤波精度。使滤波数值和真实数值更加接近,减小了跟踪误差。
4.4. 仿真实验与轨迹预测
4.4.1. 仿真实验
Singer模型有较大的目标机动覆盖范围,是有效的目标机动模型,所以Singer模型相对机动目标具有较好的跟踪性能。无迹卡尔曼滤波是基于统计线性回归的思想,与扩展卡尔曼(EKF)相比,不存在高阶截断误差,具有更大的发展优势。由于高超声速飞行器跳跃滑翔阶段机动性能好,运动特征复杂,故本文采用基于改进后的Singer-UKF对目标飞行器进行定位跟踪,为外推预测目标轨迹提供持续的跟踪信息(图13)。
取由双星定位的得到的目标信息作为目标状态向量初始值,其中初始位置(6433.53 km, 59.12 km, 5.93 km)、速度(−0.15 km/s, 6.43 km/s, 0.39 km/s)以及加速度(−0.0067 km/s2, −0.0016 km/s2, 0.0018km/s2),系统误差方阵为(单位:(km/s2)2,(km/s2)2,(km/s2)2)。目标量测初始值为(6433.8 km, 0.52˚, 0.05˚),量测噪声方差阵(单位:km2,(˚)2,(˚)2)。因为本节是基于Singer模型的双因子自适应无迹卡尔曼滤波算法对目标飞行器进行定位跟踪,经过多次仿真实验,在实时调整系统噪声中对弱化因子取0.8,遗忘因子为0.95,对时间相关模型取机动频率为0.6具有较好的跟踪效果。UKF卡尔曼滤波器取,系统仿真时间200 s,周期为1 s。结果如图14和表1。
Table 1. Statistical table of target error characteristics
表1. 目标误差特性统计表
从上面仿真图可以看出,与传统卡尔曼滤波算法对比,基于双因子加强的Singer-UKF算法在滤波初值存在误差的同时明显的减小了跟踪滤波的误差,并保持强跟踪性,跟踪性能稳定,改善了滤波发散,没有像卡尔曼滤波算法那样偶尔出现比较大的跟踪波动,同时保证在整个滤波过程中有较高的滤波精度。
4.4.2. 轨迹预测
对目标持续跟踪200 s后,开始预测目标飞行器轨迹。对目标三维速度曲线分析可知,其满足幅值衰减的周期性运动,可利用如下双正弦和函数拟合目标位置曲线,短期内这样既能达到理想精度,又能简化计算难度。
(41)
上式中,为振幅,为频率,为每一个基波的初相。
利用经过滤波后的目标状态估计值得到拟合参数,然后计算速度函数的解析导数即得到目标加速度随时间变化的表达式,以跟踪最后数据作为初值代入拟合公式中,结合常加速度模型即可递推得到目标的预测轨迹。轨迹预测误差曲线如图15。
由于目标实际运动特性十分复杂,利用跟踪数据结合拟合函数进行飞行器的轨迹预测,轨迹预测误差后面随时间后移而不断增大,这是高超声速飞行器轨迹预测现在普遍面临的问题 [10],但早期轨迹预测精度较高,误差较小,可为短期轨迹预测提供思路参考。
5. 结论
本文首先分析高超声速飞行器的运动特性得到目标运动方程,然后基于双星观测系统的目标几何定位算法得到了目标定位信息,将定位信息作为状态初值代入改进的Singer-UKF自适应算法中,从仿真曲线看出,相较于传统滤波算法,该算法具有跟踪高超声速目标更高的稳定性和优越性。最后结合滤波跟踪值,对目标飞行器进行轨迹预测,并得到较好的预测精度。