航天器面内机动规避最优脉冲分析

苏 飞,刘 静,张 耀,杨 旭,程昊文

(中国科学院国家天文台,北京 100012)

0 引 言

航天事业的发展为地球外层空间带来大量空间碎片,根据最新资料显示,目前在轨10 cm以上空间碎片数量约2万个,厘米级空间碎片接近50万个[1-4]。恶劣的空间环境对航天器的安全运行造成巨大威胁,机动规避是航天器躲避空间碎片的一种主动防护措施,受到航天机构和卫星运行商的广泛关注[5-8]。

国内外学者对航天器机动规避问题进行了大量研究。基于碰撞概率的防撞规避机动方法精度高,文献[9-11]研究了碰撞概率和推力的解析表达,结合航天器轨道控制分析规避机动问题;在此基础上,文献[12-14]推导了机动脉冲与碰撞概率的解析关系,通过拉格朗日乘子法得到固定脉冲的最优方向;文献[15-17]基于概率计算模型,通过一维碰撞概率积分方法推导速度增量的最优解,并利用两体轨道积分代替高精度轨道积分研究碰撞规避机动策略;文献[18]将机动脉冲解耦,分步求解脉冲方向和脉冲大小,在碰撞概率降低到安全值的前提下得到最优的机动冲量。许多学者通过分析规避机动动力学模型,利用寻优方法获取避碰策略,文献[19]分析了基于遗传算法的燃料消耗和机动时长,采用多目标优化求解碰撞规避机动策略;文献[20]利用非线性规划算法求解最小速度增量,得到最优的碰撞规避机动策略;文献[21]对近圆轨道提出了可以快速得到推力三维解的解析方法,利用梯度搜索算法找到最优解。一些学者从实际应用出发,建立量化的分析模型,不通过寻优直接分析次优解得到规避机动策略,文献[22]研究了卫星规避相关的轨道机动方式,结合碰撞风险建立了卫星规避方案量化分析方法;文献[23]利用根数长期项分析了初始轨道根数与危险交会点的关系,研究了卫星冲量对交会点的影响;文献[24]分析了高度规避、时间规避、正常轨控结合的空间目标碰撞规避策略;文献[25-27]研究了执行机动导致的额外风险和燃料消耗对航天器任务的影响;文献[28]在单次脉冲的基础上,分析了多次脉冲规避空间碎片能使卫星回到原始的运行轨道,同时不改变卫星原始状态。

然而多数已有成果在研究航天器机动规避时,通常将轨道摄动方程外推,并在交会点处采用数值法穷举寻优,虽然一些研究分析了碰撞概率与机动脉冲的解析关系,但在建立联系时通常利用高斯方程直接求解,计算量大且只能得到次优解。为了解决最优解问题,首先通过小偏差线性化等处理将高斯方程解析化,得到解析的轨道规避机动动力学模型,直接分析机动脉冲与航天器相对运动的线性关系。通过航天器位置协方差信息及交会平面的相对位置信息得到航天器危险交会的碰撞概率和交会距离,并建立与机动脉冲的解析关系,然后采用机动方向和机动大小分步求解的策略计算机动冲量,高效准确地求解不同时机下的最优面内脉冲。

1 机动规避动力学建模

1.1 系统描述及建模假设

利用航天器运动方程对初始状态和状态误差分布的协方差矩阵外推,就能预判航天器在运行过程中是否遭遇危险交会,如果危险交会的瞬时碰撞概率或交会距离大于安全阈值,则航天器需要进行机动来规避风险,保证运行安全。

在航天器危险交会规避任务中,往往通过改变轨道形状来躲避碎片,即机动力或机动脉冲约束在轨道面内,这主要是因为面内机动控制简单,同时改变轨道方位对航天器任务影响较大。航天器危险交会示意图如图1所示。

图1 航天器危险交会示意图Fig.1 Illustration of dangerous rendezvous trajectory of spacecraft

建模过程做如下假设:

假设1作用于航天器上的脉冲为小量,航天器在机动脉冲的作用下,到达危险交会点的速度大小及方向不变;

假设2危险交会过程为短期交会,碰撞时刻两航天器的速度为匀速直线运动;

假设3两空间目标等效为已知半径的球体,位置误差椭球在相遇期间保持不变;

假设4重力场为理想中心力场,忽略大气阻力、太阳光压等摄动因素的影响,忽略航天器姿态误差对脉冲方向的影响。

1.2 坐标系定义

相关动力学建模参照坐标系如图2所示,其中OXYZ为偏心率矢量坐标系,其原点O位于地球质心,X轴与航天器偏心率矢量eo重合,Z轴垂直于航天器轨道平面,正方向与航天器角速度方向相同,Y轴由右手螺旋法则确定;sxyz为航天器轨道坐标系,其原点s位于航天器质心,x轴由地心指向航天器质心,z轴与Z轴平行且指向相同,y轴由右手螺旋法则确定,三轴单位矢量分别为ex、ey、ez。

图2 建模参照坐标系Fig.2 Reference frame for dynamics modeling

1.3 机动动力学模型

如图2所示,航天器在偏心率矢量坐标系的状态矢量表示为

(1)

式中,eX和eY为偏心率矢量坐标系对应坐标轴的单位矢量;a为航天器轨道半长轴;E为航天器对应真近点角处的偏近点角;e为航天器轨道偏心率;r为航天器地心距;μ为万有引力常数。

令E=Em,取惯性系内航天器机动点即真近点角θm处的状态rm、vm代入式(1),求解方程式得

(2)

式中,‖‖符号表示对矢量求模,将式(2)代入式(1)得偏近点角为E处航天器的位置矢量r，表示为

(3)

结合高斯方程,冲量使航天器轨道产生的瞬时变化为

(4)

式中,Δvx、Δvy、Δvz分别为冲量在航天器轨道坐标系三轴上的分量;h、p分别为航天器变轨前的轨道角动量和轨道半通径,表示为

(5)

将航天器的状态量代入式(3),则航天器在机动点θm处机动前和机动后,运行到危险交会点的位置矢量r(Ec)、r′(Ec)分别表示为

(6)

(7)

(8)

式中,Δrx、Δry分别为航天器在轨道系中相应坐标轴的位置变化。由于航天器机动后,其偏心率矢量坐标系绕Z轴转动Δω,将式(6)、式(7)表示到机动前的偏心率矢量坐标系,则

Δr=A′r′(Ec)-Ar(Ec)

(9)

其中

(11)

(12)

(13)

(14)

为了推导方便同时减小因量级的巨大差异引起的仿真失真,进行归一化处理,定义

(15)

(16)

(17)

沿迹方向的距离变化考虑机动后,航天器运行相同相位条件下的时间差,基于文献[22]对摄动参数的定义得

(18)

其中

(19)

(20)

(21)

将式(19)～式(21)代入式(18),分别在无摄动力存在摄动的情况积分,然后做差并泰勒展开,推导过程不再赘述,得沿迹方向的距离变化为

fx2(sinEc-sinEm)+fx3(cosEc-cosEm)+

fx4(sin 2Ec-sin 2Em)+fx5(cos 2Ec-cos 2Em)]+

fy2(sinEc-sinEm)+fy3(cosEc-cosEm)+

fy4(sin 2Ec-sin 2Em)+fy5(cos 2Ec-cos 2Em)]}

(22)

其中

(23)

从动力学模型可以看出,在机动脉冲为小量时,航天器机动导致的距离变化可以简化为航天器轨道参数的解析表达,与外推模型相比可以大大节省计算时间,并能得到机动大小与距离变化的等式关系,有利于模型的物理分析和控制参数选取。

2 航天器危险交会分析

航天器危险交会过程中通常采用交会距离判定法或碰撞概率判定法来判断发生碰撞的可能性。碰撞概率判定法适用于两星的定位误差较大且有累积效应时,交会距离判定法则适用于能实时测量两星相对状态。

2.1 交会距离

航天器运行时,设s(t)为t时刻航天器s-1与s-2的相对位置矢量,则任意时刻t+Δt两个目标的相对位置矢量为

(24)

式中,vr=v1-v2,若t时刻两目标间的距离最小,则有

(25)

根据式(25)得两目标距离最小时,其相对位置矢量和相对速度矢量互相垂直。定义交会坐标系,其三轴单位矢量分别为

(26)

航天器相对运动方向平行于x轴,交会平面(b平面)定义为y-z平面,则航天器危险交会关系描述为两航天器在b平面的交会距离,航天器轨道坐标系到交会坐标系的转换矩阵为

(27)

式中,α为危险交会点处s-2速度矢量v2在s-1轨道面投影与s-1速度矢量v1的夹角(-π<α<π),且v2×v1·eZ>0时为正;β为v2在航天器s-1轨道面投影与v2的夹角(-π/2<β<π/2),且v2·eZ>0时为正;B为v1与相对速度v1-v2的夹角(0

(28)

(29)

交会距离为

(30)

2.2 碰撞概率

航天器在运行过程中位置误差服从3维正态分布,并可以通过分布中心和位置误差协方差矩阵描述。此时航天器间交会的碰撞概率均可表示为

P=∭Vf(x,y,z)dxdydz

(31)

式中,V为以一个航天器为圆心,两个航天器包络半径之和ssum为半径的球体;f(x,y,z)为高斯概率密度函数(probability density function，PDF),表示为

(32)

式中,se=(xbm,0,zbm)为b平面内两航天器交会距离最小时的相对位置矢量;|C|为相对位置s误差协方差矩阵C的行列式值,假设航天器交会时为匀速直线运动,则积分球变为圆柱,在误差协方差矩阵中选取与相对位置信息有关的部分得

(33)

此时碰撞概率表示为圆域内的积分

(34)

文献[6]提出一种不等方差PDF在圆域内积分问题的解决方法。通过将不等方差的等概率密度椭圆用与其面积相等的等概率密度圆代替,并用无穷级数表示二重积分,对式(34)简化得到碰撞概率的解析表示，即

(35)

式中,u和v均为无量纲变量,表示为

(36)

(37)

3 最优规避分析

轨道维持条件下的交会距离及碰撞概率最优的规避机动计算采用两步策略实现[16]:

步骤1通过求碰撞概率/交会距离对于机动方向的梯度确定机动速度方向。由于规避机动速度通常很小,可以将规避机动引起的交会点的航天器移动看成机动速度大小的线性函数,这样碰撞概率/交会距离的梯度与机动大小就不相关,因此,最优的机动方向与机动大小可以独立求解。

步骤2机动方向确定后,碰撞概率及交会距离即为机动大小的函数。求解满足安全碰撞值和轨道维持要求的机动速度大小就是求解一个非线性方程。

(1)交会距离最优问题描述为

(38)

机动方向可以通过梯度矢量D得到

(39)

其中,D表示为

(40)

(41)

ρ2=fun(Ec,Em,e,φ)

(42)

当确定交会目标的Ec、Em、e,最优交会距离值所需要的机动脉冲速度方向即为式(43)的解。

(43)

(2) 同理,碰撞概率最优问题描述为

(44)

机动方向可以通过梯度矢量D′得到

(45)

其中,D′表示为

(46)

此时碰撞概率相关系数描述为

v=Fun(Ec,Em,e,φ)

(47)

当确定交会目标的Ec、Em、e,最优碰撞概率所需要的机动脉冲速度方向即为式(48)的解。

(48)

4 仿真分析及校验

根据理论建模与分析的结果,建立航天器规避机动数值仿真模型和航天器危险交会分析模型。为了验证所设计的规避机动策略的有效性,本小节利用美俄卫星碰撞事件为例进行仿真。Cosmos-2251和Iridium-33中,前者为俄罗斯的废弃卫星没有任何机动能力,后者为美国的在用通信卫星可以执行规避机动,选择两颗卫星碰撞前的各自最新的一组2行轨道根数(two line mean element，TLE)根数,利用简化摄动预报模型(simplified general perturbations,SGP4)轨道预报模型,可知两颗卫星交会时间(time of the closest approach，TCA)的参数如表1所示。

表1 美俄卫星碰撞交会几何参数Table 1 Encounter geometry parameters of US & Russian satellite collision

其中,a、e、θ分别为Iridium-33的轨道半长轴、偏心率和真近点角,RI、RC分别为Iridium-33和Cosmos-2251的等效半径,α为危险交会点处Cosmos-2251速度矢量在Iridium-33轨道面投影与Iridium-33速度矢量的夹角。β为投影与Cosmos-2251速度矢量夹角。

根据Cosmos-2251和Iridium-33的预报数据,相对位置s误差协方差矩阵假设为对角协方差矩阵,其标准差在航天器相对运动切线方向为1 km,在b平面内的两个正交方向上为100 m,此时在误差协方差矩阵中与相对位置信息有关的部分[9]为

(49)

图3是航天器机动规避时不同机动提前量下的最优交会距离,从图3可得在提前量不是整数周期或不靠近整数周期时,提前机动时间越大,得到的最优交会距离越大;提前量为整数周期或在其附近时,最优交会距离回落,稍有减小。

图3 不同机动提前量的最大交会距离(Δv=0.01 m/s)Fig.3 Maximum achievable miss distance under different maneuver opportunity (Δv=0.01 m/s)

航天器执行机动规避时,通常需要考虑轨道维持,即需要在避碰的同时完成轨道抬高,所以在仿真过程中,取脉冲方向与轨道系x轴夹角为0≤φ≤π,图4是最优交会距离约束下航天器机动规避时不同机动提前量下的最优脉冲方向,从图4得航天器规避提前量在0～0.2周期时,最优脉冲方向收敛至180°(0°),提前量超过约0.2周期后,最优脉冲方向收敛至90°附近。

图4 不同机动提前量时的最优脉冲方向(Δv=0.01 m/s)Fig.4 Optimal direction of pulse under different maneuver opportunity (Δv=0.01 m/s)

图5是航天器机动规避时不同机动提前量下的最优碰撞概率,与最优交会距离分析类似,从图5可得在提前量不是整数周期或不靠近整数周期时,提前机动时间越大,得到的最优碰撞概率越小;提前量为整数周期或在其附近时,最优碰撞概率振荡且稍有增大。分析式(44)可知,最优碰撞概率是交会平面内两个距离量的二元二次函数,对两式求偏导，得

(50)

从式(50)可知,最优碰撞概率的梯度是两个距离量耦合,与误差协方差矩阵相对位置信息式(49)对应,误差协方差是表征两个航天器位置误差的累积,通过碰撞概率表示两个航天器的危险交会具有实际意义。

图5 不同机动提前量时的最优碰撞概率(Δv=0.01 m/s)Fig.5 Optimal achievable collision probability under different maneuver opportunity (Δv=0.01 m/s)

取脉冲方向与轨道系x轴夹角为0≤φ≤π,图6是最优碰撞概率约束下航天器机动规避时不同机动提前量下的最优脉冲方向,从图6得航天器规避提前量在0～0.25周期时,最优脉冲方向收敛至180°(0°),提前量超过约0.25周期后,最优脉冲方向收敛至90°附近。

图6 不同机动提前量时的最优脉冲方向(Δv=0.01 m/s)Fig.6 Optimal direction of pulse under different maneuver opportunity (Δv=0.01 m/s)

图7是航天器机动规避中,采用不同脉冲时不同机动提前量下的最优交会距离,从图7可得相同提前机动时机下机动脉冲越大,最优的交会距离越大。

图7 不同机动提前量时的最大交会距离 (不同脉冲)Fig.7 Maximum achievable miss distance under different maneuver opportunity (different impulsive)

图8是最优交会距离约束下,航天器机动规避中,采用不同脉冲时不同机动提前量下的最优脉冲方向,从图8得脉冲变化不影响最优机动方向。

图8 不同机动提前量时的最优脉冲方向 (不同脉冲)Fig.8 Optimal direction of pulse under different maneuver opportunity (different impulsive)

图9是航天器机动规避中,采用不同脉冲时不同机动提前量下的最优碰撞概率,从图9可得相同提前机动时机下，机动脉冲越大,最优的碰撞概率越小。

图9 不同机动提前量时的最优碰撞概率 (不同脉冲)Fig.9 Optimal achievable collision probability under different maneuver opportunity (different impulsive)

图10是最优碰撞概率约束下,航天器机动规避中,采用不同脉冲时不同机动提前量下的最优脉冲方向,从图10得脉冲变化不影响最优机动方向。

图10 不同机动提前量时的最优脉冲方向 (不同脉冲)Fig.10 Optimal direction of pulse under different maneuver opportunity (different impulsive)

5 结 论

基于数理变化推导出解析的轨道规避机动动力学模型,建立了脉冲大小与机动位移的量化关系;利用航天器交会过程的协方差信息等得到碰撞概率和交会距离,通过代入机动位移分析了规避模型的优化目标函数;最后搭建仿真模型研究了航天器面内最优机动规避问题,为机动规避的工程问题提供支持。