赵建磊, 李海阳,*

(1. 国防科技大学空天科学学院, 湖南 长沙 410073;2. 空天任务智能规划与仿真湖南省重点实验室, 湖南 长沙 410073)

0 引 言

随着人类探索地球外层空间活动的不断深入,空间飞行器尤其是非合作空间飞行器的数量日益增长。出于对有限空间资源的争夺,恶意抵近服务如干扰、侦照、抓捕等行为逐渐增加,飞行器在轨安全受到严重威胁。通常情况下,飞行器根据实际任务需要对自身进行的轨道机动、轨道维持、抵近观测、姿态控制等操作的信息是不公开的。在此背景下,非合作目标的轨道机动信息的感知、检测和识别对及时发现非合作飞行器的恶意抵近服务意图,维护在轨飞行器的安全显得尤为重要。

国内外学者从不同角度对轨道机动检测和机动信息识别等问题进行研究。Kelecy等人对脉冲小推力作用下的定轨与机动检测进行研究,对比分析了最小二乘估计方法与卡尔曼滤波方法的适用性。Lemmens等人分别提出了基于两行轨道要素(two line elements, TLE)一致性检验的低轨卫星机动检测方法和支持近实时机动检测的TLE时间序列分析方法。Pittelkau和Acquesta等人分别提出了机动检测中避免历史数据中噪声干扰的新方法,其中前者采用中值滤波器,而后者则采用被称为Box具有噪声的基于密度的聚类方法(density-based spatial clustering of applications with noise, DBSCAN)。Bai等人利用卫星的历史轨道信息采用聚类的方法进行了机动检测,分析了均值聚类、层次聚类和模糊聚类在机动检测中的应用。Yu等人通过分析机动前后的摄动,利用相对动力学方程推导出两种机动模型,提出了基于轨道动力学模型的机动检测方法。Daniel等人提出了一种新的基于马氏距离的二元假设检测方法,用于表征和检测航天器异常。Bao等人提出了一种基于贝叶斯理论和多模型的滤波器,使用最小均方误差方法来适应机动目标可能进行的多种运动,以此进行联合检测和弱目标跟踪。Clark等人通过将TLE分析和状态传播这两种众所周知的检测方法融合成混合检测方法来提高检测的计算速度和精度提高,提高了空间态势感知(space situational awareness, SSA)算法的性能。Justin等人利用光声特征检测方案实现了航天器近实时机动检测和参数估计,增强了对小机动事件的识别能力。Jiang等人探讨了跟踪执行未知连续机动的卫星问题,提出一种估计卫星状态和机动加速度的鲁棒扩展卡尔曼滤波器,与无偏最小方差输入和状态估计方法相比,该方法在具有相同的位置估计精度的同时机动期间速度估计误差降低了5倍左右。Christopher从目标历史TLE记录中获得的数据训练神经网络,用于发现新的TLE记录中的异常轨道机动、错误标记或其他原因引起的意外偏差,该方法能够避免依赖于个别目标轨道状态或特定的传播模型,提高了检测的准确性。Li等人利用目标的历史TLE数据,分别用统计的方法和自组织映射聚类的方法对LEO轨道卫星的机动检测问题进行了研究。崔红正等人提出针对不同推力非合作目标的机动检测方法,结合地基与天基观测数据,可以满足多数非合作目标轨道机动检测需求。王庆瑞等人将距离变化率作为特征量,提出一种基于概率判决模型的自适应轨道机动检测方法,该方法基于Neyman-Pearson准则求解序列数据的变轨判决门限。

以上方法中基于历史轨道信息的机动检测方法只能判断是否发生机动,不能对机动参数进行有效的定量识别,而对于机动信息的识别,大多是基于天基或地基平台对目标的连续测角或测距追踪,结合动力学模型和滤波方法进行机动轨道定轨及机动识别。但对于多数非合作目标,难以对其进行连续追踪观测并且动力学模型是未知的,只能通过稀疏短弧段观测数据在无机动假设下定轨得到若干个离散的轨道状态。而当卫星发生轨道机动时,通过轨道测量得到的卫星轨道根数相比于不发生轨道机动时的轨道根数,通常含有较大的误差,且在短时间内有限条轨道根数无法通过滤波等手段处理误差。

针对以上问题,本文基于稀疏短弧段观测数据定轨得到的有限条轨道测量状态,对非合作飞行器的机动参数进行识别。首先建立基于稀疏轨道根数的飞行器轨道反演模型,然后基于轨道反演对飞行器机动参数进行识别,最后通过算例对本文的方法进行校验。

1 非合作飞行器轨道反演模型

1.1 基本假设

在本文的研究中,做以下假设。

(1) 当变轨推力作用时间远远小于变轨前后的轨道周期时,可以假设推力随时间变化的函数是脉冲函数,使得脉冲矢量与原推力产生的冲量一样,因此本文中涉及的变轨冲量均考虑为脉冲冲量。

(2) 轨道测量仅存在随机误差,不存在系统误差。

(3) 以美国战略司令部发布的TLE中轨道机动后恢复时间为例,空间飞行器在发生机动后,通过测定轨得到稳定的TLE轨道通常需要2～7天,本文假设飞行器在得到稳定的轨道之前不会施加下一次轨道机动。

(4) 在非合作飞行器的一个轨道周期内,至少能得到一段短弧段观测数据,结合假设(3),即对于轨道周期小于2天的非合作飞行器,相邻短弧段观测数据之间至多考虑存在一次轨道机动。

1.2 冲量作用下飞行器轨道预报模型

在二体假设下,天体为质点,空间飞行器只受中心引力作用。而实际飞行任务中,中心天体一般为非理想球体,因此飞行器还会受中心天体非球形引力、大气阻力、太阳光压、第三体引力等各种摄动因素的影响。

在J2000地心惯性坐标系中的飞行器动力学模型如下所示:

(1)

当飞行器的初始运动状态可知,在脉冲机动作用下,飞行器在任意时刻的轨道状态表示如下:

(2)

1.3 飞行器轨道反演

(3)

定义轨道反演的状态偏差为在所有测量点时刻(=0,1,…,),反演轨道状态与测量状态的加权偏差和,即

(4)

(5)

=DU/VU

(6)

(7)

式中: ,(=,,)为对应于沿轨道不同方向的测量误差的比例系数。

状态偏差反映了通过轨道测量数据反演的轨道与测量轨道的接近程度,状态偏差越小,说明反演的轨道与测量轨道更接近,同时与测量轨道的离散状态相比,通过测量状态反演的轨道是一条连续的轨道,减小了测量数据中的随机误差的影响,反演得到的轨道更加接近真实轨道。因此,可以通过最优化方法使状态偏差最小,以此反演非合作飞行器的轨道。

2 基于轨道反演的机动参数识别方法

2.1 机动参数识别优化模型

(1) 优化变量

由式(3)和式(4)可知,状态偏差与反演轨道的初始轨道状态和中间机动参数有关,即

(8)

因此,将优化变量选为反演轨道初始状态和中间机动,即

(9)

(2) 优化变量边界约束

(10)

式中:,,,,,分别为飞行器轨道的半长轴、偏心率、轨道倾角、升交点赤经、近地点角距和真近点角;和分别为根据飞行器轨道类型设定的半长轴边界。

(3) 目标函数

目标函数取为状态偏差最小,即

(11)

值得注意的是,由于在设计变量中引入的中间机动通常是多于实际机动的,多余的机动在优化结果中为小值,因此在优化流程的最后进行一次小机动值的修正,即若存在设计变量|Δ|(=1,2,…,)的数值小于一定门限,认为在该时刻未发生机动。

(12)

2.2 “全局寻优+局部修正”的双重优化求解

上述问题可以归结为非线性规划问题,关于这种非线性规划问题的求解,有很多算法可以使用,如差分进化(differential evolution, DE)算法、序列二次规划(sequential quadratic programming, SQP)算法等。此处选取DE算法和SQP算法,其中DE算法全局搜索能力较好,但效率不高,尤其当轨道外推使用高精度模型进行时,会花费大量时间;序列二次规划算法局部寻优效果好,但对初值比较敏感。因此,在使用高精度轨道外推模型求解时,为了提高搜索效率以及避免初值陷入局部最优,提出一种“全局寻优+局部修正”的双重优化算法进行计算,计算流程如图1所示。

图1 双重优化求解流程

在二体条件下,通过DE算法寻找满足约束的全局最优解;

以步骤1中最优解为初值,利用SQP算在高精度模型条件下,寻找满足约束的最优解;

剔除小机动值,重新进行一遍SQP优化流程。

3 仿真校验

本节对上述的飞行器轨道反演模型及机动参数识别方法进行仿真校验。仿真在计算机处理器为Intel i7-7700 CPU 3.60 GHz,运行内存为16 GB的C++环境下进行。

3.1 仿真参数设定

地球同步静止轨道由于具有独特的高轨和地球同步特性,成为了通信、气象、导弹预警、广播电视和数据中继等重要卫星的首选轨道,但同样也将面临严峻的恶意服务考验。因此,对地球同步静止轨道卫星机动信息的感知、检测和识别显得尤为重要。本节以地球同步静止轨道为例,对本文提出的模型和方法进行校验。

表1 仿真初始轨道根数

表2 标称轨道中间机动参数

对上述标称轨道每隔7 h选取该时刻的轨道状态作为测量样本,并施加不同的均值为0的高斯随机定轨测量误差,以此模拟实际情况下稀疏短弧段定轨得到的有限条离散轨道测量数据。不同工况下的定轨误差如表3所示。

表3 定轨误差标准差

3.2 算例分析

使用本文描述的机动参数识别方法,不同测量点时刻对应状态偏差的权重因子均取为1,优化目标函数取为式(7),目标函数中沿轨道不同方向的测量误差的比例系数, (=,,)取值如下所示：

(13)

其中，取为015 m/s,DE算法种群规模为240,代数为300,SQP算法最大迭代代数为50。

每种工况分别进行50次蒙特卡罗仿真,机动参数识别结果与标称机动参数的偏差如表4所示。从结果来看,定轨误差越大,则机动参数的识别误差越大。

表4 机动参数识别结果

为进一步说明双重优化算法的收敛情况及机动参数识别效果,选取工况2中的一次仿真,初始轨道测量数据如表5所示。双重优化算法的迭代收敛情况如图2所示,仿真显示双重优化算法中DE算法在100代左右基本收敛,SQP算法在20代左右基本收敛。

图2 双重优化算法收敛情况

表5 离散轨道测量数据

在表5中，定轨结果以经典轨道根数的形式给出,括号中数据依次代表半长轴、偏心率、轨道倾角、升交点赤经、近地点角距、真近点角,其中半长轴的单位为m,角度的单位为(°)。

同时,为比较说明本文优化算法的优势,分别用DE算法和SQP算法对算例进行求解。表6给出了不同算法求解的最大机动识别误差和计算时间,比较可知,SQP算法机动识别误差较大,原因是SQP算法对初值较为敏感,随机给定初值容易使算法陷入局部收敛;DE算法计算时间较长,原因在于直接使用高精度轨道外推模型进行全局搜索时,由于计算量骤增,导致计算效率极低。相对而言,本文提出的“全局寻优+局部修正”的双重优化算法对本文的模型适应性很好。

表6 不同算法机动识别误差及计算时间

经SQP高精度修正后的飞行器机动参数反演结果如表7和表8所示,分别给出了剔除小机动前后的机动参数识别情况。可以看出,剔除小机动后的识别结果与不剔除小机动相比结果更加接近标称轨道施加的机动,机动识别的速度误差在0.08 m/s以内,时间误差在3 min以内。

表7 飞行器机动参数识别结果

表8 飞行器机动参数识别结果(剔除小机动)

为进一步说明反演轨道与实际轨道的偏差,图3给出了在仿真时段反演得到的轨道与标称轨道的位置偏差及轨道迹向、法向和径向偏差。可以看出,在整个仿真时段,反演轨道与标称轨道的最大位置偏差不超过1 200 m。

图3 反演轨道与标称轨道的偏差

4 结 论

本文在稀疏轨道信息条件下,针对非合作飞行器轨道机动参数的识别问题,探讨了一种可行的思路。建立了飞行器轨道反演模型及基于轨道反演的机动参数识别方法,并提出相应的“全局寻优+局部修正”的双重优化求解算法。仿真结果显示,该方法能够在稀疏轨道信息的条件下,对非合作飞行器的机动参数进行较为有效的识别;空间飞行器定轨误差越大,则机动参数的识别误差越大;工况2中算例机动参数识别结果显示,剔除小机动后的识别结果与不剔除小机动相比结果偏差更小,反演轨道与标称轨道的位置偏差不超过1 200 m,机动识别的速度误差在0.08 m/s以内,时间误差在3 min以内。