哈工大深空探测轨道设计课程作业.docx

地球—火星转移轨道设计
1.研究问题描述
在未来3年内寻找发射机会，设计地球－火星转移轨道，以总能量最小为指标。假设地球停泊轨道为高度200km、轨道倾角28.5º的圆轨道，目标轨道为高度500km、倾角90 º的火星绕飞轨道。轨道动力学模型只考虑太阳、地球、火星的引力作用及地球J2项摄动。
分别给出初步轨道设计参数和精确轨道设计参数。
2.解决方案
首先，基于pork-chop图寻找合适的发射窗口，然后采用圆锥曲线拼接法，进行轨道的初步设计，分段用二体模型设计各阶段轨道，最后在精确力学模型下，采用微分修正法，基于B平面参数进行精确轨道设计。
2.1地火转移轨道的初步设计
对于初步设计，是以圆锥曲线拼接法(Patched Conic Method)为基本原理，将探测轨道划分为几个分段过程，每个分段过程可以近似成二体问题处理，然后通过接口，将每一段拼接在一起。它是在二体模型假设下，通过求解Lambert问题确定发射窗口，利用圆锥曲线拼接法确定地心段、日心段、火星段的初始轨道参数。
2.1.1日心轨道设计及发射窗口的搜索
本文采用等高线图法搜索发射窗口，通过给出既定时间段内所有的发射和到达时间情况，获得初始和终了位置，进而解算Lambert问题，从而绘制“猪排”图，观察得到最优发射机会。
基于等高线图的最优发射机会搜索算法的主要步骤如下：
（1） 根据任务的需要确定出发射机会搜索的目标函数（性能指标）、发射时间的区间以及飞行时间的区间；
（2） 选取一组出发时刻t0和达到时刻tf，根据行星历表计算地球的位置RE(t0)、速度VE(t0)和火星的位置RM(tf)、速度VM(tf)；
（3） 利用转移时间（tf - t0）、RE(t0)、以及RM(tf)，通过求解Lambert问题，可以得到探测器在始末位置处的速度矢量V1(t0), V2(tf)；
（4） 确定发射机会的目标函数，并绘制出等高线图；
（5） 根据等高线图，找到目标函数取值较小的区域，确定出性能指标指标最优的发射时间；
以上步骤中涉及到的目标函数（性能指标）通常是指双曲线超速v∞，发射能量C3，速度增量Δv，它们的具体定义为：
双曲线在达到地球引力影响球边缘时速度有剩余，这个双曲线剩余速度v∞通常称为双曲线超速，计算公式为
(1-1)
式中v1是飞行器的速度矢量，vLE是发射时刻地球绕太阳公转的速度矢量。
发射能量是影响任务初始设计的关键参数，在飞行器质量一定的情况下，发射能量越大，所需运载火箭的运载能力越强，它的大小是发射时双曲线超速的平方
(1-2)
当假设飞行器从停泊轨道开始转移，这在停泊轨道上施加的速度增量Δv1为
(1-3)
类似地，可以定义飞行器到达目标天体的双曲线超速以及到达目标停泊轨道，交会时的速度增量
(1-4)
(1-5)
式中v2是飞行器到达目标行星时的速度矢量，vAT是发射时刻目标行星绕太阳公转的速度矢量，rTp为飞行器在目标天体俘获是近心点的高度。
在整个飞行过程中，所需的总的速度增量Δvtotal为
(1-6)

2.1.2地心段参数的确定
对地心段参数进行初步设计时，假设在地球影响球内探测器只受到地球的引力作用，从而按照二体轨道特性进行轨道参数的计算，逃逸速度增量假设为脉冲。假设探测器的发射过程为:探测器从地面发射后，首先进入圆形的停泊轨道;然后在经历一段时间的滑行后，在特定的时间，特定的位置，经过末级火箭加速后进入双曲线轨道进行逃逸，如图 1所示。为了尽量避免不必要的燃料消耗，假定停泊轨道和双曲线轨道是共面的，且双曲线的近地点半径与圆轨道半径相同，速度相切。因此，地心段参数初步设计的内容主要包括两个方面:停泊轨道参数的特性分析以及停泊轨道参数选择。

图 1地心逃逸双曲线轨道面内参数的几何关系
设停泊轨道半径(即双曲线近地点半径)为今1，根据双曲线轨道参数的相关公式容易算得双曲线轨道半长轴、偏心率、动量矩以及近地点速度:
(1-7)
其中a1为轨道半长轴，e1为偏心率，为地心引力常量，v∞1为地心双曲线
剩余速度大小，hl为双曲线轨道的动量矩大小，vp1为双曲线近地点速度大小。
而根据剩余速度矢量v∞1可以计算出转移双曲线渐近线对应的赤经α∞1（RLA: Right ascension）和赤纬δ∞1（DLA: Declination）
(1-8)
(1-9)
进一步可以确定B平面中的渐近线的单位方向矢量
(1-10)
而B平面的其他两个坐标系可以表示为
(1-11)
(1-12)
式中，B平面的夹角可以由渐近线的赤纬和轨道倾角确定
(1-13)
双曲线的单位角动量矢量，可以表示为
(1-14)
当速度趋于无穷时，真近点角的正弦和余弦可以确定为
(1-15)
进一步可以确定转移双曲线轨道的近地点矢量
(1-16)
从而可以确定航天器在双曲线转移轨道的近地点处的位置矢量为，而速度的单位矢量为
(1-17)
而航天器在近地点的速度矢量可以表示为。根据近地点的位置矢量和速度矢量，就可以计算出双曲线轨道的轨道六根数。
考虑到行星际距离比地球影响球大很多，故地心段双曲线轨道与地球影响球的交点位置对轨道几乎没有影响，而v∞1矢量的大小和方向完全决定了探测器在脱离地球作用范围之外的运动，因此控制v∞1矢量大小和方向的精度非常重要。实际中存在无数条满足剩余速度大小方向约束的逃逸双曲线存在，这些双曲线位于将图 1中双曲线轨道绕过地心的v∞1矢量(即渐近线的方向)旋转一周组成的曲面上。并且所有双曲线的近地点组成的圆形轨迹称为轨道发射圆，因此停泊轨道必然经过OL点，等探测器运行到轨道发射圆的位置，沿着速度方向施加脉冲，进入双曲线轨道，如图 2所示。

图 2轨道发射圆示意图
设v∞1矢量在地心惯性下可以表示为:
(1-18)
则地心到OL点的矢量的赤经为（π+α∞1），赤纬为（-δ∞1），表示为：
(1-19)

2.1.3火心段参数的确定
对火心段参数进行初步设计时，可以认为在火星影响球内探测器只受到火星的中心引力的作用，从而按照二体轨道特性进行轨道参数的计算，捕获速度增量假设为脉冲。火心段双曲线轨道的特性与地心段双曲线的性质完全类似。在地心段双曲线轨道的剩余速度约束是设计的目标，而在火星段双曲线剩余速度的约束是初始条件，其大小和方向是固定的。
火星段轨道参数的选择主要内容是确定双曲线的近火点半径和轨道倾角。假设进入双曲线的轨道和探测器的工作轨道是共面的。为了实现火星探测器对于火星的全面覆盖，火星探测的轨道选取为极地轨道，倾角为i2=90o，设双曲线的近火点高度为rp2。
在进行轨道设计时，通常将目标的双曲线的参数用B平面参数来描述。探测器进入火星影响球时，可以根据火心双曲线进入速度v∞2与目标火心轨道的倾角和近心距rp2计算标称B平面参数BT和BR。
(1-20)

其中为火星引力常量，φ为B矢量与T矢量的夹角，如图 3所示，i2为目标火心轨道的轨道倾角，速度v∞2为火心双曲线进入速度大小，δ∞2为该矢量的火星赤纬。

图 3火心进入双曲线轨道与B平面示意图

2.2基于B平面参数的精确轨道设计
初步轨道设计中，只考虑二体模型，存在很大误差，必须进行修正，本文采用微分修正法，基于B平面参数进行修正，从而得到精确的转移轨道。
2.2.1问题描述
本文考虑地球、火星和太阳的引力和地球J2项摄动力的影响下，其动力学模型如下
(1-21)
其中是地球J2项摄动力加速度。
在探测任务中，火星探测器从地球逃逸轨道出发飞到火星目标停泊轨道，需要对初始设计参数进行修正，从而使火星探测器在精确力学模型下，准确到达火星目标停泊轨道。火星探测器的终端参数选择为目标轨道的倾角和B平面参数。选择好终端参数，记初始时刻为的控制参数为，火星探测器抵达目标区间的终端参数记做，则火星探测器初始状态和终了状态之间可以用某个函数来表示，即
(1-22)
微分修正问题就是采用一定的迭代方法使探测器抵达目标区间的实际状态与期望状态的误差小于规定值
(1-23)
2.2.2 制导方法
我们将实际轨道在标称轨道附近进行泰勒展开后只保留线性项，得
(1-24)
其中，是被控制量，敏感矩阵。
求解式(1-27)的方法主要考虑和的维数，其分别为和
1）
此种情况较容易，只需对敏感矩阵求逆便可求得
(1-25)
2）
此种情况控制量的数量多，在式的约束下，求得幅值最小的修正量即可。其性能指标为
(1-26)
其中是Lagrange乘子。可利用变分法求解式，即
(1-27)
3）
此种情况控制量较少，可采用最小二乘法求解该问题。即求式的最小值
(1-28)
其解是，但此解不能保证小于规定值。

2.2.3 轨道精确设计求解
本文采用微分修正的方法，给出一种数值求解式中的偏导数矩阵的方法。下面是以B平面参数为终端参数并采用微分修正进行解精确轨道设计，具体步骤如下：
步骤一、通过地火转移轨道的初步设计，得到地球停泊轨道参数、停泊轨道上的速度增量、日心转移轨道参数和进入火星影响球时的速度增量的初值，选取地球停泊轨道的升焦点赤经、近地点幅角、停泊轨道上的速度增量和进入火星影响球时的速度增量为控制参数。
步骤二、根据控制参数为初值在精确动力学模型下进行轨道数值积分运算，求得终端参数值，终端参数选取B平面参数（、）和目标轨道倾角。
步骤三、计算得到的参数值与标准参数进行比较，获得参数偏差量
，从而求得新的控制参数。
步骤四、利用新的控制参数重新对动力学模型进行轨道积分运算，得到新的终端参数值偏差，这个偏差量逐渐减小。
步骤五、重复上述过程，直到终端参数满足精度要求。
其程序流程图如图 4所示：

图 4基于微分修正的精确轨道设计流程图

2.3仿真分析
2.3.1初步轨道参数设计结果
地球和火星的会合周期为约为780天，为了寻找未来三年内的发射机会，以2016年5月10日起始，搜索之后3*365天的发射机会，可以得到从地球发射飞行器所需的能量C3的变化情况如图 5所示，而轨道转移所需的总速度增量变化情况如图 6所示：

图 5 未来3年内地球到火星任务所需发射能量的变化等高线图

图 6 未来3年内地球到火星任务的总速度增量变化等高线图
通过对等高线图进行分析，以最小速度增量为优化目标时，卫星的发射时间为15-May-2018，总的飞行时间为208天，此时速度增量取最小为5.7653 km/s。
对地球—火星探测轨道进行初步设计，可以确定出转移时刻停泊轨道以及双曲线轨道对应的轨道六根数如表 1所示
表 1 初步设计的轨道参数
轨道类型
半长轴a（km）
偏心率e
轨道倾角i(rad)
近地点幅角ω(rad)
升交点赤经Ω(rad)
真近点角θ(rad)
停泊轨道
6578.13700
0
0.4972188
0
0.8971506
2.1701513
转移轨道
-54029.004
1.121751
0.4974188
2.1701513
0.8971506
0

可以计算出在停泊轨道需要施加的速度脉冲大小为=3.9981km/s，而速度脉冲在地球惯性坐标系下的矢量可以表示为[-1.78419414174206，
-2.70261001844018，-2.16699941349961]km/s。

2.3.2 精确轨道参数设计结果
基于B平面参数对轨道进行精确设计，可以确定发射时刻为12-May-2018，00:15:33.137，到达时刻为10-Dec-2018，00:12:47.156，飞行器的转移时间为211.9981天，整个转移轨道如图 7所示：

图 7 J2000惯性坐标系（以太阳质心为坐标原点）下探测器的转移轨道图形
通过计算可以得到在地球停泊轨道以及双曲线轨道的轨道参数，如表 2所示：
表 2 精确设计出的轨道参数
轨道类型
半长轴a（km）
偏心率e
轨道倾角i(rad)
近地点幅角ω(rad)
升交点赤经Ω(rad)
真近点角θ(rad)
停泊轨道
6578.13700
0
0.4974188
0
0.8971770
2.1701669
转移轨道
-54028.155
1.121753
0.4974188
2.1701669
0.8971770
0

发射时，在地球停泊轨道上，在探测器所施加的脉冲为dv1=
[-1.78420583339426，-2.70264143209095，-2.16701775917320]km/s，脉冲大小为=3.8966km/s；从发射时刻起，当t = 18057600s时，探测器到达火星影响球，此时施加的脉冲为dv2=[0.0901724915339645，3.86006626112751，0.0383287816881011]km/s，脉冲大小为=3.8613 km/s；
当t = 18316634.0185162 s时，探测器到达火星停泊轨道，此时施加的脉冲为dv3= [0.898193335636322，-0.561658888501794，1.62920970970374]km/s，脉冲大小为= 1.9433 km/s。dv1、dv2和dv3均为J200惯性系下的矢量。
最终到达火星目标停泊轨道，目标轨道参数如表3所示：
表 3 火星目轨道参数
轨道类型
半长轴a（km）
偏心率e
轨道倾角i(rad)
近地点幅角ω(rad)
升交点赤经Ω(rad)
真近点角θ(rad)
停泊轨道
3896
0
1.5707963
0
2.9540761
2.7754371