基于偏转距离近似模型的动能撞击小行星防御任务脉冲轨道优化研究1)

王艺睿 李明涛 周炳红

(中国科学院国家空间科学中心,北京100190)(中国科学院大学,北京 100049)

引言

太阳系内的小行星主要分布于火星轨道与木星轨道之间,也有一部分小行星会穿越地球轨道威胁到地球上的生命.诸多严重的小行星撞击地球事件(如6500 万年前的K-T 事件[1]、1908 年通古斯大爆炸事件[2]、2013 年车里雅宾斯克事件[3]等)已引发人类对于行星防御的关注.在能够对近地小行星提前预警的前提下,摧毁小行星结构或者偏转小行星轨道是避免其撞击地球的两种基本方式[4].目前已提出的防御手段主要有核爆[5]、动能撞击、激光烧蚀[6]、拖船[7]、引力拖船[8]、离子束牵引[9]等,其中动能撞击技术是目前最易实现且成熟度较高的手段.

动能撞击技术是指撞击器以一定的速度和角度撞击小行星,使小行星轨道发生改变.2005 年美国实施了深度撞击(deep impact)任务[10],任务释放一颗约370 kg 的撞击器,以约10.3 km/s 的相对速度撞击“Tempel 1”彗星的核,撞击后的3 年内该彗星位置变化了约10 km;双小行星重定向测试任务(double asteroid redirection test,DART)将于2021 年7 月执行,由一个约650 kg 的撞击器以约6.65 km/s 的相对速度撞击小行星Didymos 的卫星,预计产生0.8 ∼2 mm/s 的速度改变量(取决于动量传递因子的大小,动量传递因子用于描述溅射物对动量传递的影响[11]).Atchison 等[12]总结了DART 轨道设计需求,包括撞击日期、撞击角度以及为了满足光学需求的太阳相位角需求;Ozimek 和Atchison[13]研究了DART 轨道中逃逸段的月球借力方法;Atchison 等[14]介绍了DART的终端制导、载荷等信息.

早期对于动能撞击处置小行星的研究,多关注动能撞击对小行星的轨道偏转机理,一般不考虑运载能力约束.Ahrens 和Harris[15]的研究表明切向偏转速度增量对偏转小行星轨道效率最高,并估计了偏转小行星所需的速度增量,但如果小行星在撞击地球前会近距离飞掠大行星,实际所需的偏转速度增量会有所降低[16];Park 和Ross[17]推导了小行星与地球距离的一、二阶导数表达式,发现只有在撞击前某些特殊时刻最优偏转速度才平行于小行星速度方向;Ross 等[18]基于圆锥曲线拼接法考虑了地球引力对最优偏转速度增量的影响.Carusi 等[16]提出一种估算在多体动力学模型下将小行星偏转1 倍地球半径所需速度增量的方法,但仅研究了切向速度增量;Kahle 等[19]对其方法进行了改进,在多体问题下同时优化了切向、径向的速度增量,发现Carusi 等[16]的“切向速度增量”结论仅在二体阶段下是合理的,当小行星近距离飞掠行星时,最优速度增量方向并不是切向的了,近距离飞掠行星时受行星引力影响可能会降低偏转所需的速度增量.Vasile 和Colombo[20]发现当偏转时间小于小行星公转周期时,最优速度增量的径向分量占主导,当偏转时间更大时,切向分量占主导.Yamaguchi 和Yamakawa[21]提出了IGMs(impactgeometry maps),用于可视化偏转距离与IG(impact geometry)(小行星速度与航天器撞击速度的内积)之间的关系.Greenstreet 等[22]基于统计的方法,发现与行星发生近距离飞掠的小行星,所需的偏转速度增量可以非常小.

动能撞击任务设计的一种优化指标是最大化小行星被偏转前后的近地点距离改变量(偏转距离),换言之是使小行星被偏转后的近地点尽可能的离地球更远[23].采用多体型计算小行星被偏转前后的近地点距离,可以提供关于问题最完整的解,但轨道演化和近地点搜索计算过程中会产生冗余计算量,在任务设计初期的快速搜索研究(quick scoping studies)中,多体数值积分是不必要的[24].采用二体模型替代多体模型的方法可以显著降低轨道演化过程的计算量,但在搜索近地点的过程中仍存在冗余计算量.为同时降低轨道演化和近地点搜索的计算量,可采用解析偏转模型对偏转距离进行估算.

除了Ahrens 和Harris[15]提出的偏转距离线性解析模型外,目前使用较为广泛的的偏转距离解析模型还有:2001 年Conway[25]提出了基于拉格朗日系数解析递推偏转距离的方法(后文简称Conway 模型);2007 年Izzo[26]基于几何法推导了偏转距离的解析模型(后文简称Izzo 模型);2008 年Vasile 和Colombo[20]基于高斯摄动方程推导了偏转距离的解析模型(后文简称Vasile 模型),其相比于Conway 模型具有更高的计算效率.Vasile 模型是目前使用较为广泛的解析偏转模型,文献[27-28]均基于Vasile 模型,以偏转距离为指标优化了动能撞击任务.但是目前尚无研究对偏转距离解析模型的最优性进行分析.

小行星被偏转前后,虽然轨道根数发生了改变,但到达地球近地点的时刻可能变化不大.本文首先研究了近地点时刻受偏转影响的变化规律,提出在动能撞击任务的设计初期,可以假定小行星被偏转前后近地点时刻不变(近似偏转模型),由此可以省略掉近地点时刻搜索的计算量.相比于解析模型,这种近似模型可以在提升最优性的同时、降低求解复杂性.考虑运载约束,将化学推进变轨简化为脉冲推力变轨,建立了直接转移(两脉冲及三脉冲)和行星借力飞行转移(单次及两次借力)的动能撞击轨道优化模型.以偏转小行星Apophis 为例,设计了动能撞击任务,研究了利用两脉冲/三脉冲直接转移、单次/两次行星借力飞行技术是否可以提升偏转距离.通过对比二体动力学模型与高精度动力学模型对计算偏转距离的误差,进一步探讨了在动能撞击任务设计过程中采用二体动力学模型的合理性.

1 偏转距离的现有计算方法

本文中偏转距离的定义为:偏转前后小行星近地点(closest-approach,CA)距离的改变量(偏转距离).偏转距离δρCA的数学表达式为

其中,ρ(t)为t时刻小行星的地心位置矢量,a和b为时间搜索窗口的下界和上界.关于偏转距离的计算方法,本节将介绍数值模型、两种经典解析模型(Vasile 模型和Izzo 模型).

1.1 数值模型

对于小行星的高精度轨道演化动力学模型,本文考虑的摄动力有八大行星引力、月球引力、四大主带小行星引力(Ceres,Vesta,Pallas,Hygiea)、相对论效应(general relativity,GR)、地球非球形引力J2、太阳光压(solar radiation pressure,SRP)以及Yarkovsky效应(yarkovsky effect,YE).高精度动力学模型的数学表达式如下所示

其中,dpi表示第i个行星地日心矢量,ρpi表示第i个行星到小行星的矢量.行星和月球的星历调用JPL DE430 历表[29],小行星星历调用JPL Horizons On-Line Ephemeris System[30]公布的真实星历.动力学模型中考虑非引力项(太阳光压、Yarkovsky 效应)的原因为:太阳光压对一些尺寸较小的小行星影响较大,对于尺寸较大的小行星,Yarkovsky 效应不能忽略[31-33].利用Runge-Kutta-Fehlberg7(8)积分器[34-35]对轨道运动方程的一阶形式进行数值积分,并采用一种一维最小值算法[36]求解近地点距离,因为这种算法不需要地心距离的一阶导数信息.

为验证本文高精度轨道预报方法的准确性,本节对小行星Apophis 在2029 年的近距离飞掠地球事件进行了预报,并与JPL 公布的结果[30]进行了对比.小行星Apophis 是目前行星防御领域高度关注的小行星,根据JPL[30]公布的最新信息(表1),Apophis 将在2029 年4 月13 日与地球近距离飞掠,近地点距离为0.00025AU,其质量约为6.1×107t.将Apophis 在2019 年1 月1 日00:00:00(TDB)的位置速度(坐标系为ICRF)作为积分起点.

表1 小行星Apophis 初始位置速度[30]Table 1 Initial state of Apophis

表2 给出了JPL 公布的结果与本文预报的结果.仿真结果表明,利用本文方法计算的Apophis 近地事件,时间误差约为0.02 s,位置误差约为1 km.

表2 小行星Apophis 近地事件预报结果Table 2 Results of Apophis’closest approach

理论上,在动能撞击轨道优化问题中,需要利用高精度动力学模型进行递推以保证最优性,但是各种摄动力的引入会显著降低求解效率.采用二体模型(仅考虑太阳引力)替代多体模型,可以显著降低轨道演化过程的计算量,但仍需数值计算小行星的近地点.二体模型的数学表达式如下所示

假定小行星被偏转前后日心位置r(td)不变、仅日心速度v(td)发生改变,即偏转使小行星状态由(r(td),v(td))改变为(r(td),v′(td)).将偏转前后小行星的状态利用二体模型由偏转时刻(td)解析递推(fkepler)至近地点时刻(tCA),可求得小行星近地点距离.小行星被偏转前近地点位置的日心矢量r(tCA)及地心矢量ρ(tCA)表示为

其中,近地点时刻tCA采用一维最小值算法[36]对式(2)进行求解,可得偏转距离为

1.2 解析模型

1.2.1 Vasile 模型

Vasile 模型可以解析计算小行星被偏转前后MOID (minimum orbit intersection distance)的改变量.其基本思想是将MOID 点处位置矢量的改变量分解为径向、切向及法向分量

其表示为轨道根数改变量的函数如下

其中,θMOID表示近地点时小行星的日心真近点角,.轨道根数改变量与偏转速度增量之间的关系可由高斯摄动方程计算,由于撞击改变了小行星的轨道周期,因此近地点处的平近点角改变量还需考虑由于轨道周期改变而产生的长期演化影响,公式如下

其中,θd表示撞击时小行星的日心真近点角,σ 表示撞击时刻距离近地点时刻的时间.

虽然MOID 概念并不等同于近地点距离[20,23],但通过观察公式可以发现,如果将MOID 点真近点角θMOID替换为近地点真近点角θCA,则该公式可以用于估算近地距离改变量.δα(td)表示偏转位置处小行星轨道根数的改变量,δv(td)表示偏转位置处小行星速度的改变量,δr(tCA)表示近地点位置的偏移量,ACA和Gd表示状态转移矩阵

近地点距离的改变量δρCA可表示为

其中,ρCA表示偏转前的近地点矢量.

1.2.2 Izzo 模型

Izzo 模型的基本思想是,通过计算偏转前后近地点时间的差异来计算小行星在B 平面上的偏移量∆ξ.其中,B 平面定义为垂直于小行星进入地球影响球双曲线速度U且经过地球质心的平面[37],öpik[38]在B 平面上定义了(ξ,η,ζ)坐标系:η 轴沿U的方向,ζ 轴为地球日心速度在B 平面上投影的反方向,ξ轴的指向由右手定则给出,在该坐标系下ξ 坐标可代表小行星轨道与地球轨道的MOID 值,B 平面上的偏移量∆ξ 可用来描述小行星的偏转距离.在日心视角下,小行星被偏转前后近地点位置的几何关系如图1所示.

图1 Izzo 模型原理示意图Fig.1 Schematic diagram of Izzo model

图1 中,∆σ 表示偏转前后近地点时刻的改变量,vast表示小行星的近地点日心速度,U表示小行星进入地球影响球时的相对速度,ψ 表示小行星速度vast与U之间的夹角,θaE表示地球速度vE与U夹角的补角.∆ξ 表示小行星被撞击后的偏转距离,按照图中的几何关系,可建立偏移量∆ξ 与交会时刻改变量∆σ 之间的等式如下

经过对∆σ 表达式的推导(具体步骤参考文献[26]),可将上式整理如下

其中,Usc表示撞击时撞击器相对于小行星的速度.

Izzo 模型可以比较直观地给出偏转距离的影响因素,比如小行星本身的物理属性(小行星质量mast、动量传递因子β 等)、小行星与地球的轨道关系(小行星半长轴aast、近地点时地球的速度vE等)以及轨道优化结果(撞击器质量msc、偏转时间σ、撞击速度Usc等).该公式的最后一项点乘表示在偏转时间较长的情况下,在小行星速度增量的切向分量是影响偏转距离的主要因素.另外需注意到,Izzo 模型只能计算小行星在B 平面内的偏移量,并不等于近地点距离的改变量,因此当采用该模型估算真实近地点距离改变量时会出现偏差.图2 展示了一种利用Izzo 模型计算产生偏差的情况.

在该情况下,即使偏转使小行星产生了较大的∆ξ,但实际上小行星被偏转后近地距离减小了.也就是说,Izzo 模型无法回答小行星被撞击后近地点距离会更近还是更远的问题.

图2 Izzo 模型产生偏差的情况Fig.2 Deviation of Izzo model

2 近地点时刻预估近似模型

如1.1 节所述,求解小行星近地点距离改变量的本质在于求出偏转前后的近地点时刻.如果能够提前估算出小行星被撞击后近地点的时刻,那么就可以利用轨道解析递推模型,快速计算出被偏转后的近地点矢量.下面通过公式推导分析近地点时刻的影响因素.

以Izzo 模型为基础,在已知小行星轨道根数的情况下,撞击时刻与近地点时刻之间的时间间隔σ可由下式计算

其中,td为偏转时刻,∆Mast为平近点角差值,nast为转速.对上式取微分并忽略掉短期影响项后得

式中,aast为小行星的轨道半长轴,vast为小行星被撞击时的日心速度,δvt为小行星被撞后产生的切向速度增量.对上式进行整理后,可得近地点时间改变量的解析表达式如下

经过推导可以发现,在小行星速度增量的各分量中,主要是切向分量影响近地点时刻的改变量.以偏转Apophis 为例,下面研究了当撞击产生不同的切向速度时,对近地点时刻的影响.以Apophis 在二体模型下求得的2029 年近地点时刻为参考(2029-04-13 21:40:09),当撞击使Apophis 产生不同的切向速度增量(δvt=0.2 ∼1 mm/s)时,Apophis 近地点时间改变量∆σ 随偏转时间σ (2 ∼10 a)的变化趋势,如图3所示.

图3 近地点时刻改变量的变化曲线Fig.3 Change of perigee epoch

图3 中坐标为(6,17.55)点的含义为:在Apophis到达近地点前6 年(2023/4/15 21:40:08)时进行撞击,如果撞击使Apophis 产生的切向速度增量为1 mm/s,那么将会使Apophis 在2029 年的近地点时刻晚17.55 s,也就是Apophis 被撞击后的近地点时刻为2029/4/13 21:40:26.从这幅图,可以获得两方面的信息:第一方面,在10 年预警时间的情况下,单颗动能撞击器能使Apophis 近地点时刻改变不到30 s;第二方面,由Izzo 模型可知,∆σ 与偏转距离成正比,因此随着偏转时间σ 的增大,近地点距离改变量(偏转距离)整体呈现波动上升的趋势.

下面仅着眼于偏转后的轨道,分析以不同时刻计算近地点时,近地近地距与真实近地距的差异.小行星进入地球影响球后,由于相对速度较高(数千米每秒),其运动轨迹近乎直线.小行星被偏转后,在地球影响球内的运动轨迹如图4 所示.其中,U表示小行星进入地球影响球时的相对速度,表示被偏转后的真实近地点时刻,表示近似近地点时刻,.近似后的近地距满足下面等式

图4 小行星真实近地距与近似近地距Fig.4 Real perigee distance and approximate perigee distance

根据1.1 节的数值模型,可求得Apophis 进入地球影响球的相对速度∥U∥≈5.9 km/s,假设被偏转后的真实近地距离为36 000 km,计算当近地点时间改变量∆σ ∈[0,30]s 时,近似近地距与真实近地距之间的相对误差,计算结果如图5 所示.

图5 近似近地距与真实近地距之间的相对误差Fig.5 Relative error between real perigee distance and approximate perigee distance

由图5 可知,在近地点改变时间小于30 s 时,近似近地距与真实近地距之间的相对误差不足0.01%.例如,当近地点时刻相差30 s 时,近似近地距为36 000.4 km,与真实近地距差别小于1 km.近地距越大,二者之间的误差越小.

因此,对于本文偏转Apophis 算例,在任务设计初期,可以忽略偏转前后小行星近地点时间的改变量,即.小行星的近地点距离改变量可由下面的近似模型进行估算

其中∥ρ′(tCA)∥表示小行星被偏转后,在偏转前近地点时刻tCA处的地心矢量.依据该模型,已知撞击后小行星的轨道根数,可以解析地快速求得tCA时刻轨道根数,省略了用于搜索近地点时刻的计算量.

3 动能撞击任务设计与优化建模

动能撞击任务流程图如图6 所示.本章将首先建立航天器(撞击器)轨道转移模型,然后建立动能撞击任务优化模型,最后介绍本文采用的优化方法.

图6 动能撞击任务流程Fig.6 Process of a kinetic impactor mission

3.1 转移模型

本文采用圆锥曲线拼接法求解行星际转移轨道.圆锥曲线拼接法的基本原理为将转移轨道按照天体引力影响球进行分段,航天器在影响球内只考虑中心天体的引力,最后通过航天器在影响球边界处的状态对多段圆锥曲线进行拼接,天体影响球边界的相对速度称为逃逸速度或双曲线超速.该简化算法在深空轨道设计中是通用且高效的[39].下面对两脉冲转移模型、三脉冲转移模型、借力飞行转移模型进行建模.

3.1.1 两脉冲转移模型

首先,航天器在t0时刻从地球逃逸,经过∆t1天后与撞击危地小行星.假定转移期间不施加任何速度增量,由运载火箭直接发射送入撞击轨道(图7).

图7 两脉冲直接转移模型Fig.7 Two-impulse direct transfer model

通过求解Lambert 问题[40],可以计算出在地球影响球边界处的逃逸速度v∞以及航天器撞击小行星的相对速度Usc,对应的数学表达式为

其中,r0和rf分别表示航天器在t0和t0+t1时的日心位置矢量,v0和vf为对应的日心速度矢量,vE和vast分别表示地球和小行星的日心速度矢量.v∞可由运载火箭提供,对应的C3为

假设转移过程不施加任何机动,则航天器在撞击时的质量等于发射入轨时的质量

其中,撞击时刻td=t0+t1,msc(t0)=f(C3).f(C3)为运载能力拟合公式,CZ-5 运载火箭的运载能力曲线如图8 所示.

图8 长征五号运载能力曲线Fig.8 CZ-5’s launch performance

当动能航天器撞击目标小行星后,超高速撞击会导致小行星表面产生溅射物,一部分溅射物最终会再次吸积在小行星表面,另一部分速度较高的溅射物将会从小行星的引力场中逃逸,这种现象类似于火箭产生动力的机理,将会放大动量传递效应,可采用动量传递因子β 描述溅射物对动量传递的影响[11].假定航天器与小行星的撞击过程为完全非弹性碰撞,根据动量守恒原理,小行星被撞击后产生的速度改变量可由下式估算[41-42]

其中,msc(td)和vsc(td)表示航天器在撞击小行星时的质量和速度,mast和vast(td)表示小行星被偏转时的质量和速度.本文不考虑溅射物的影响,因此取β=1.

综上所述,给定节点(t0,∆t1)后,可以确定航天器在撞击小行星前的速度vsc(td)及质量msc(td).由式(25)计算出小行星被偏转后的状态,从而可计算出小行星被偏转前后的近地点距离改变量.

3.1.2 三脉冲转移模型

航天器在t0时刻从地球逃逸,经过∆t1天后撞击危地小行星,期间在t0+α1∆t1时施加一次深空机动DS M1,α1的取值范围为0 到1(图9).

图9 三脉冲直接转移模型Fig.9 Three-impulse direct transfer model

对于三脉冲转移模型,需要给定航天器从地球逃逸时的速度大小及方向,对应3 个参数(C3,DLA,RLA).由于航天器在转移过程中施加了深空机动,所以航天器在撞击时的质量应减去转移过程的燃料消耗

其中,撞击时刻td=t0+t1,转移速度增量总和∆vtot=DS M1,msc(t0)=f(C3),Isp表示发动机比冲.

综上所述,给定入轨参数(C3,DLA,RLA)、时间节点(t0,∆t1,α1)后,可以确定航天器在撞击小行星前的速度及质量,由式(25)计算出小行星被偏转后的状态,从而可计算出小行星被偏转前后的近地点距离改变量.

3.1.3 借力飞行转移模型

为了提高动能撞击的偏转效率,本文引入借力飞行技术,探讨借力飞行技术是否对偏转效率有所提升.借力飞行的基本原理是利用自然天体的引力,对航天器的速度进行改变,如图10 所示.

图10 借力飞行示意图Fig.10 Schematic diagram of gravity assist

受自然天体引力作用,在自然天体引力影响球边界处航天器的进入与离开之间的夹角为

其中,µga和Rga分别表示借力天体的引力常数和半径,hp表示借力高度,vga表示借力天体的速度.定义一个新的坐标系o-ijk

因此,借力后矢量可描述如下

其中,φ 表示在j-k的投影与k轴正向之间的夹角,取值范围为[0,2π).

本文考虑的借力天体有金星、地球和火星.首先,航天器在t0时刻从地球逃逸,经过∆t1到达借力天体,再经过∆t2天后撞击危地小行星.转移期间,还需执行两次深空机动,DS M1施加于t0+α1∆t1时以及DS M2施加于t0+∆t1+α2∆t2时,α1与α2的取值范围为0 到1(图11).

图11 借力飞行转移模型Fig.11 Gravity assist transfer model

同样的,利用式(26)可以计算出航天器在撞击小行星时的质量,其中撞击时刻td=t0+t1+t2,转移速度增量总和∆vtot=DS M1+DS M2.

综上所述,给定入轨参数(C3,DLA,RLA)、时间节点(t0,∆t1,∆t2,α1,α2)、借力参数(hp1,φ1)后,可以确定航天器在撞击小行星前的速度及质量,由式(25)计算出小行星被偏转后的状态,从而可计算出小行星被偏转前后的近地点距离改变量.以此类推,每增加一个借力天体,需增加4 个新的优化变量(∆t,α,hp,φ).

3.2 优化模型

对于两脉冲直接转移模型,共有2 个决策变量

对于三脉冲直接转移模型,共有6 个决策变量

对于单次借力飞行转移模型,共有10 个决策变量

对于两次借力飞行转移模型,共有14 个决策变量

分别利用数值模型、Vasile 模型、Izzo 模型和近似模型构建目标函数(偏转距离),其表达式如下所示:

(1)利用数值模型构建目标函数

(2)利用Vasile 模型构建目标函数

(3)利用Izzo 模型构建目标函数

(4)利用近地点时刻预估模型构建目标函数

3.3 优化算法

通过下述优化算法,可计算出动能撞击任务中直接转移方案/借力飞行转移方案的最优转移轨迹.图12 给出了优化算法的流程图.

图12 优化算法流程图Fig.12 Optimization process

第一步:输入决策变量范围;

第二步:生成初始种群;

第三步:计算发射C3,并判断是否在拟合区间内(本文CZ-5 拟合曲线需C3<60 km2/s2).如果不在区间内,返回上一步;如果符合区间范围,进入下一步;

第四步:利用直接转移模型或行星借力借力飞行转移模型计算小行星被偏转后的状态;

第五步:按照指定的目标函数式(34)∼式(37),计算每个个体的适值函数;

第六步:判断是否达到终止条件.如果未达到条件,进入第七步;如果达到终止条件,直接进入第八步;

第七步:选择种群中的优质个体,进行遗传运算,更新种群,返回第三步;

第八步:输出适值度最高的个体,即输出最优转移轨迹.

4 仿真结果与讨论

本文以偏转Apophis 为例,在10 年预警时间的条件下,参考CZ-5 的运载能力,设计了直接转移和行星借力飞行转移的动能撞击偏转任务.

4.1 偏转距离计算方法比较

本节以直接转移方案为例,对比分析了利用不同模型进行优化的结果,对优化模型中的目标函数进行了选择.

航天器逃逸时间t0的范围为2020 年1 月1 日至2021 年1 月1 日,转移时间∆t1∈[100,1400]d.分别以数值模型、Vasile 模型、Izzo 模型、近似模型作为目标函数,利用遗传算法对直接转移方案进行了全局优化.该优化问题中,种群个体数设置为1000,种群代数设置为10.将优化出的发射窗口代入数值模型中计算数值偏转距离,表3 对比了采用不同目标函数时求得的最优解.

通过表3 优化结果的对比,可以发现在优化问题中,以Vasile 模型和近似模型作为目标函数时,可以较好地替代数值模型.尤其对于近似模型,利用其优化出的发射窗口与数值模型的基本一致,同时计算时间相比于数值模型降低了约90%;Vasile 模型也可以较好还原优化结果,但优化出的发射窗口与数值模型的存在偏差,会损失一定的最优性;利用Izzo 模型作为目标函数,即使其目标函数值达到398 km,但优化出的解不具备参考性,这是由于小行星被偏转前后出现了图2 所示的位置关系.这也验证了1.2.2 节中关于Izzo 模型初步分析的结论,不可将其作为优化的目标函数.从Izzo 模型的优化结果也可以得知,不管撞击器对小行星的偏转能力有多强,一旦选择了错误的偏转时机,则可能会给地球带来更大的威胁.

为了直观对比数值模型、解析模型、近似模型的发射窗口差异,下面在同一张图中绘制出利用不同模型求解出的偏转距离等于180 km 的等高线(图13).

表3 基于不同偏转距离模型的优化结果Table 3 Optimization results based on different deflection models

图13 不同模型的180 km 等高线图Fig.13 180 km contour maps of different models

图13 中,黄色粗线表示利用数值模型求解的180km 等高线,红色和蓝色细线分表表示用近似模型和Vasile 模型求解的结果,“*”标识出的位置表示利用数值模型优化出的发射窗口.红线与黄线以及蓝线与黄线的交叉点,表示与数值模型误差为0%的点.可以看出,红线与红线的拟合性较好,说明近似模型可以较好的替代数值模型;而蓝线与黄线交叉的点较少,说明Vasile 模型的解与数值解存在一定偏差.同时,可以明显的看到“*”标识落在蓝线上,这也说明如果以Vasile 模型的解作为优化目标,优化出的窗口会与实际最优窗口有一定偏差.

4.2 直接转移方案

4.2.1 两脉冲转移首先利用Pork-Chop 图分析发射窗口,具体计算的方法为:将航天器逃逸时间t0的范围选为2020年1 月1 日至2021 年1 月1 日,转移时间∆t1∈[600,1400]d,以1 d 为步长,利用近似模型计算所有情况产生的偏转距离,最终绘制得偏转距离等高线如图14 所示.

图14 Apophis 偏转距离等高线图(2020 年)Fig.14 Apophis deflection distance contour map(2020)

由图14 中可知,在2020 年中,利用两脉冲转移方案偏转Apophis 共有两个发射窗口,一个位于2020年5 月2 日前后,另一个位于2020 年5 月13 日前后.这2020 年中两个发射窗口的具体信息总结如表4.

可以看到,2 号窗口虽然拥有更大的撞击器质量和Apophis ∆v,但最终与1 号窗口的偏转距离相差不大.这是由于2 号窗口所需的转移时间更长,偏转相位次于1 号窗口,从而导致Apophis 被偏转后到达2029 年近地点的偏转时间降低.因此,不能直接用撞击体质量、速度增量等因素单方面来评估偏转距离,而是需要综合考虑多项影响因素,如Izzo 模型所示,包括小行星本身的物理属性(小行星质量mast,β 等)、小行星与地球的轨道关系(小行星半长轴aast、近地点时地球的速度vE等)以及轨道优化结果(撞击器质量msc、偏转时间σ、撞击速度Usc等).

表4 2020 年内两个发射窗口的优化结果Table 4 Optimization results of two launch windows in 2020

以总任务时间较少的1 号窗口为例,直接转移方案的转移轨迹如图15 所示.航天器自2020 年5 月2日从地球逃逸,经过670 d 后,在2022 年3 月4 日与Apophis 相撞.经过撞击,Apophis 在2029 年的近地距可增加约192 km.

图15 2020 年1 号发射窗口的两脉冲直接转移轨迹Fig.15 Two-impulse direct transfer trajectory of the 1st launch window in 2020

4.2.2 三脉冲转移

将航天器逃逸时间t0的范围选为2020 年1 月1日至2021 年1 月1 日,转移时间∆t1∈[100,1400]d.优化结果如表5 所示,转移轨迹如图16 所示.

表5 三脉冲直接转移方案优化结果Table 5 Optimization results of three-impulse direct transfer trajectory

图16 三脉冲直接转移轨迹Fig.16 Three-impulse direct transfer trajectory

仿真结果表明,三脉冲直接转移解与两脉冲直接转移的1 号窗口非常接近,三脉冲解虽然具有更高的发射入轨质量,但是由于施加了中途机动,同时又会损失一部分撞击器质量,利用三脉冲直接转移方案对偏转距离不会有明显提升.

4.3 借力飞行转移方案

航天器逃逸时间t0范围为2020 年1 月1 日至2021 年1 月1 日,转移时间∆t∈[100,1000]d,hp∈[100,3000]km,φ ∈[0,2π]rad,C3∈[0,50]km2/s2,DLA∈[0,π]rad,RLA∈[0,2π]rad.分别设计了如下的借力飞行偏转方案:金星、地球、火星、金星−金星、金星−地球、金星−火星、地球−金星、地球−地球、地球−火星、火星−金星、火星−地球、火星−火星.遗传算法的参数设置:单次借力的种群个体数设置为5000,种群代数设置为30;两次借力的种群个体数设置为8000,种群代数设置为30,设计结果见表6.

表6 借力飞行转移方案优化结果Table 6 Gravity assist transfer trajectory optimization results

仿真结果表明,只有单次地球借力的转移方案对偏转距离有所提升,相比直接转移方案提升了约20 km.关于单次地球借力转移方案的细节如表7所示.

表7 地球借力转移方案优化结果Table 7 Optimization results of Earth gravity assist transfer trajectory

地球借力转移方案的转移轨迹如图17 所示.航天器自2020 年2 月19 日从地球逃逸,2021 年5 月17 日进行地球借力,2023 年1 月25 日与Apophis 相撞,航天器转移时间共计1071 d.经过撞击,Apophis在2029 年的近地距可增加约210 km.

图17 地球借力转移轨迹Fig.17 Earth gravity assist transfer trajectory

由借力飞行转移方案的设计结果可知,利用借力飞行未必能够提升偏转距离.这是因为借力飞行概念的初衷是降低深空探测过程的燃料消耗,而对于动能撞击这种特定的深空探测任务来说,其性能指标是偏转距离,这是一个需要综合考虑燃料消耗、转移时间、撞击相位等因素的指标.如果在航天器转移的过程中执行不当的借力飞行,航天器为了与借力天体相位匹配,可能会错过最佳发射窗口,从而降低偏转距离.

5 二体模型对防御窗口的影响

以上仿真均是在二体模型的假设下进行的,本章将主要讨论基于二体模型与基于高精度模型优化结果的差别.

实际上,在高精度模型中某些摄动力对于偏转距离的影响非常微弱,在优化过程中是可以忽略的因素.因此,首先对小行星Apophis 在运动过程中所受的摄动力进行分析,筛选出对偏转距离影响最大的摄动力,对高精度模型进行简化,从而减少优化过程中不必要的计算量.以小行星Apophis 为例,下图给出了Apophis 在2019 年1 月1 日至2030 年1 月1 日之间各项摄动加速度的变化趋势.其中,图18为内太阳系4 颗行星、广义相对论、太阳光压以及Yarkovsky 效应的摄动加速度变化曲线;图19 为外太阳系5 颗行星、3 颗矮行星、4 颗主带小行星、地球非球形引力J2 的摄动加速度变化曲线.

图18 Apophis 摄动力量级(1)Fig.18 Orders of magnitude of Apophis perturbation(1)

由图18 可知,Apophis 在2029 年到达近地点之前,金星、地球、火星和月球的摄动加速度会发生量级上的显著变化.不同摄动源按照摄动加速度量级的分布区间总结于表8 中.

图19 Apophis 摄动力量级(2)Fig.19 Orders of magnitude of Apophis perturbation(2)

表8 Apophis 摄动力量级Table 8 Orders of magnitude of Apophis perturbation

基于3.1.1 节所述的两脉冲转移模型,以第4 章两脉冲转移仿真结果中的1 号窗口(2022 年5 月2日)为例,在2022 年3 月2 日对小行星Apophis 施加0.38 mm/s 的速度增量,分别利用高精度模型和二体模型对偏转后的小行星进行递推,计算二体模型与高精度模型对计算偏转距离的误差.仿真发现,利用二体模型计算偏转距离为192.03 km,利用高精度模型计算的偏转距离176.27 km.进一步,在二体模型上每次增加一项摄动力,用不同力模型递推求得的偏转距离,结果如表9 所示中,发现对偏转距离影响最大的摄动力为地球和金星引力.在考虑金星+地球摄动模型基础上,表10 每次增加一项摄动力,分析除金星和地球外,对偏转距离影响最大的因素.由表9和表10 可知,对小行星Apophis 偏转距离影响最大的因素为金星、地球、木星.在二体模型基础上引入地球、金星、金星的引力摄动后,计算出的偏转距离与高精度解的误差在1%以内.

在保证计算精度的前提下,为了提高计算速度,本文引入金星、地球、木星引力摄动对二体模型修正,利用修正模型替代高精度模型.修正模型的数学表达式如下

表9 在二体模型基础上每次增加一项摄动力Table 9 Add one perturbation source at a time based on the two-body model

表10 在金星+地球摄动模型基础上每次增加一项摄动Table 10 Add one perturbation source at a time based on the Venus+Earth model

以直接转移方案为例,分别基于二体模型和修正模型对偏转Apophis 的动能撞击任务进行优化.航天器的逃逸时间t0的范围选为2020 年1 月1 日至2021 年1 月1 日,转移时间为∆t1∈[600,1400]d.种群个体数设置为200,种群代数设置为10.利用二体模型及修正模型优化的结果如表11 所示.提取出二者优化出的防御窗口,代入高精度模型计算出对应的高精度结果.

表11 二体模型与修正模型优化解对比Table 11 Comparison of optimization results between two-body model and modified two-body model

仿真结果表明,利用二体模型进行优化,在不影响防御窗口的情况下,计算时间会比精度更高的修正模型更少.在动能撞击任务前期设计过程中,可以利用二体模型快速评估防御效果.

6 结论

本文对动能撞击小行星防御任务的脉冲轨道优化方法开展了研究.为提升动能撞击任务的求解效率,将两种经典的偏转距离解析模型引入优化模型,同时提出一种基于近地点时刻预估的偏转距离近似模型;考虑运载约束,建立了直接转移和行星借力飞行转移的动能撞击优化模型.

仿真结果表明,相比于经典的偏转距离解析模型,本文提出的基于近地点时刻预估的近似模型更加简洁直观,在提升计算效率的同时,也能保证优化问题的最优性.三脉冲直接转移方案与两脉冲直接转移方案的最优偏转效果基本一致,借力飞行转移方案相比于直接转移方案对偏转距离的提升效果并不明显.在动能撞击任务的前期设计中,可以基于二体模型进行防御效果的快速评估,虽然对偏转距离有一定误差,但对防御窗口的优化结果影响不大.