基于普适变量法的火星探测器轨道初步设计及仿真

周 杰，张树瑜，刘付成

（1.上海航天控制工程研究所，上海 200233；2.上海航天技术研究院，上海 201109）

0 引言

火星探测是深空探测中的一个热点领域，历来受到世界航天大国的重视。嫦娥探月工程拉开了我国深空探测的序幕，火星探测成为下一个目标。火星探测器轨道初步设计方法中，重点是求解兰伯特问题。GAUSS最早获得了巨大成功，但其算法远不能满足当代航天动力学的要求［1］。BATTIN，VAUGHAN分别改进了Guass算法，使它对所有可能实现的轨道均收敛，同时满足较高的数值精度，但计算方法较复杂［2－3］。BATE，BOND 提出并改进的普适变量法，其优点是使用同一组公式就可计算各种圆锥曲线的轨道，且计算过程较简单，便于计算机编程，同时其精度和收敛性也较好［4－5］。本文基于该方法对基于普适变量法的火星探测器轨道初步设计及仿真进行了研究。

1 火星探测器轨道初步设计

根据能量获取方式，火星探测器轨道可分为大推力变轨、小推力变轨和利用天体引力辅助变轨三种基本类型。本文讨论典型大推力变轨方式。

探测器在从地球逃逸轨道向环绕火星的目标轨道飞行过程中，始终受到地球、太阳、火星和其他行星引力的影响。这是一个N体问题，无法获得解析解。本文引入影响球概念，影响球模型是将N体问题转为一系列的二体问题，可描述为：每个大质量的天体都有一定的引力支配区域，探测器在一个天体影响球内只受该天体单独的引力作用，当探测器远离该天体时，飞行到该天体引力影响球边界，则开始进入另一个天体的引力支配区域。

根据影响球概念，完整的火星探测器轨道需经过地球发射和逃逸段的地球影响球、地火转移段的太阳影响球和被火星捕获段的火星影响球三个不同的区域。采用的求解影响球模型方法为圆锥曲线拼接法，它基于假设：

a）在探测器飞行的任意时刻，探测器仅受到一个天体引力作用；

b）假设地球和火星的运行轨道为同面的圆轨道。

在上述假设下，地火转移能量最省轨道是霍曼转移轨道，其近日点位于地球的日心轨道轨道上，远日点位于火星的日心轨道上，如图1所示。

图1 火星探测器轨道Fig.1 Trajectory of Mars probe

1.1 地火转移轨道段

地火转移轨道的初步设计步骤如下。

步骤1：确定从地球出发日期和到达火星日期。

步骤2：查星历表，确定出发时在日心惯性坐标系中地球的位置R1和速度v1，以及到达时火星的位置R2和速度v2。

步骤3：求解兰伯特问题，确定出发和到达时探测器的双曲剩余速度。

步骤1中，理论上出发日期可任选，但实际受能量的限制，只能在特定的一段时间内选取，称之为发射窗口。对不同的出发和到达日期，探测器实现地火转移所需的能量均不同。实际上，能量最省发射机会出现的频率取决于地球和火星的会合周期。

设地球和火星环绕太阳运动的角速度分别为ωE，ωM，则角速度差

地球和火星相对位置重复出现两次的时间间隔

若令PE，PM分别为地球和火星环绕太阳运动的轨道周期，则地球和火星的会合周期还可表示为

由式（2）或式（3）可得，地球和火星的会合周期为777.9d，即能量最省的发射机会约每26个月会出现1次。发射日期可选在能量最省发射日期附近，此时地球、火星与太阳连线间的夹角约44°。

步骤2中，查星历表，得到出发和到达日期地球和火星在日心黄道惯性坐标系（X轴指向J2000春分点方向，Z轴沿黄道面法向，Y轴由右手准则确定）中的状态向量。

步骤3中，确定出发和到达时探测器的双曲剩余速度v∞。该速度是指探测器在地球影响球边界处相对地球的速度。可将步骤2中求出的出发时地球的日心位置矢量和到达时火星的日心位置矢量近似视作飞行器在离开地球影响球和到达火星影响球时的位置矢量，通过求解兰伯特问题可得探测器在地球影响球边界处（或火星影响球边界处）的日心位置矢量和日心速度矢量，这样就可解得探测器的逃逸双曲剩余速度（或到达双曲剩余速度）。

求解兰伯特问题是地火转移轨道设计中的关键。兰伯特问题可描述为：椭圆弧上两点间的飞行时间只取决于椭圆的半长轴、弧上两点至焦点的距离之和，以及连接弧上两点的弧长。其数学表达式为

式中：t为椭圆弧上两点间的飞行时间；a为椭圆轨道的半长轴；r1，r2为弧上两点到焦点的距离；c为连接弧上两点的弧长。对兰伯特问题，一般解法为：给定r1，r2以及从r1至r2的飞行时间t，求航天器在r1，r2处的速度矢量v1，v2［6］。

由共面矢量的基本定理：若A，B，C为共面矢量，且A，B不共线，则C可由A，B线性组合表示。因开普勒运动限制在一平面内，有

式中：f，g为拉格朗日系数。可推导出拉格朗日系数与真近点角变化值Δθ间的关系以及v1，v2的表达式为

式中：μ为中心天体引力常数；L为角动量，且L＝｜r1×v2｜［6］。

拉格朗日系数确定后，兰伯特问题即可获解。本文用BATE的普适变量法求解兰伯特问题。

设普适变量χ与近点角间的关系为（近地点处t0＝0）

式中：E，F分别为椭圆和双曲线轨道的偏近点角。

设t0为χ等于零的时刻，则时刻t0＋Δt的χ值可由普适开普勒方程

迭代解出。此处：r0，vr0分别为时刻t＝t0的半径和径向速度；α为长半轴a的倒数，α＜0，α＝0，α＞0分别对应为双曲线、抛物线和椭圆；C（z），S（z）为Stumpff函数，且

其中：z＝αχ2，z＜0，z＝0，z＞0分别对应为双曲线、抛物线和椭圆。

用χ，C（z），S（z）表示的拉格朗日系数为

上述各式中未知量为L，χ，z，已知量为Δθ，Δt，r1，r2。用普适变量法求解兰伯特问题时需用逐次逼近法，求解过程如下：

a）计算r1，r2，Δθ。

c）令

e）计算拉格朗日系数

f）计算v1，v2。

至此，采用普适变量法求解兰伯特问题过程结束。

在上述第二步中已求得地球和火星在日心惯性系中的状态向量Rearth，vearth，RMars，vMars。在第三步中，假设探测器在地球影响球边界处的日心位置矢量即为地球的位置矢量R1＝Rearth，在火星影响球边界处的日心位置矢量即为火星的位置矢量R2＝RMars，求解兰伯特问题得到探测器在地球影响球边界处的日心速度vdeparture和在火星影响球边界处的日心速度varrive。则逃逸双曲剩余速度v∞，departure（探测器在地球影响球边界处相对地球的速度）和到达双曲剩余速度v∞，arrive（探测器在火星影响球边界处相对火星的速度）分别为

引入深空探测器轨道设计中的重要参数C3，定义C3＝（v∞）2。设探测器在半径为r的地球停泊圆轨道上加速进入一逃逸双曲轨道，其速度增量可表示为

可见，发射C3越大，探测器逃逸地球引力场所需的速度增量就越大，其大小表征了转移轨道对发射能量需求。同样，到达C3与探测器到达火星近火点时制动所需的Δv直接相关，表征了转移轨道对探测器机动能力需求。地火转移轨道设计中，节省能量是考虑的首要问题，故希望的发射C3和到达C3较小。

不同出发日期和到达日期对应的发射C3和到达C3不同。本文选定一段范围的出发日期和到达日期（如出发日期为2013年10月1日至2014年3月1日，到达日期为2014年5月20日至2015年1月20日），用上述解法求得每一对出发日期和到达日期相对应的发射C3和到达C3，结果分别如图2、3所示。

图2 地球发射C3等高图Fig.2 Lnuch C3contour map

图3 到达火星C3等高图Fig.3 Arrival C3contour map

由图2、3可知：发射能量呈山脊分布，共分为上下两个区域，每个区域有一局部极小点，其中下半区对应的是Ⅰ类转移轨道，其日心转角小于180°，上半区对应的是Ⅱ类转移轨道，其日心转角大于180°。相对而言，Ⅱ类转移轨道的最小能量略小于Ⅰ类转移轨道，但所花时间更长，两者各有优缺点。

同理可得到达火星C3等高图。对某些探测任务，常希望用较小的速度脉冲制动就能形成环绕火星的目标轨道，这样就可通过到达火星C3等高图确定合适的发射机会。

另外，从发射和到达C3等高图还可知：在2013年12月1日至2013年12月20日期间发射所需的C3能量较小，且到达火星制动所需的能量也较小，相对应的到达火星的日期为2014年10月1日至2014年10月30日。对于发射质量约1 000kg的探测器，我国CZ－3A运载火箭所能提供的最大C3能量大小约9.5km2／s2，即在发射C3等高图中C3值小于9.5的区域均可行。如以2013年12月15日为出发日期，2014年10月20日为到达日期，由兰伯特问题求解得到发射C3≈8.955 49km2／s2，到达C3≈12.052 8km2／s2，飞行时间约309d。

1.2 地球逃逸轨道段

作过地心与v∞平行的直线，绕此直线旋转逃逸双曲线可得一旋转曲面，其近地点在空间形成一圆形轨迹，该曲面上任一双曲线均可作为逃逸双曲线。

为保证发射面与地球停泊轨道面及逃逸双曲线在同一平面内，设A为过地心和v∞反向延长线上停泊轨道高度上的点（如图4所示），则发射场、地心和点A确定的平面即为发射轨道面，探测器应在点A从发射轨道进入停泊轨道，滑行至前述圆形轨迹上任一点，即在逃逸双曲线近地点上加速进入逃逸轨道。

图4 发射方位角和发射窗口Fig.4 Lunch azimuth and Lunch window

对某一发射点，因安全考虑对发射方位角有限制。1d中，待发射点随地球自转至点1时，对应发射的最小方位角，此时发射窗口打开；转至点2时，对应最大方位角，此时发射窗口关闭。

1.3 中途修正原理

根据前述圆锥曲线拼接法设计了初步的地火转移轨道，但由于发射产生的误差及探测器在飞行过程中受到的各种摄动力的影响，实际轨道会产生偏差，偏离预定的轨道，若不进行轨道修正，则其到达火星时会由于非线性效应产生巨大偏差。计算表明：在300km的地球停泊轨道，近地点发动机关机时速度变化0.01%（1.1m／s）时，日心霍曼转移椭圆目标轨道半径将变化153 000km；发动机关机时近地点半径变化0.01%（0.67km）时也将产生偏差70 000km。因此，为到达预定的火星探测目标轨道，须在地火巡航段完成对探测器进行多次轨道修正，一般地火转移段需作2～3次轨道修正，首次修正一般在探测器脱离地球影响球后初期进行，以后根据实际再进行数次微小修正。

轨道修正一般是改变探测器的速度。如当探测器飞行到火星的影响球范围时，其相对火星的速度可认为是v∞，则通常对其进行调整，设其改变量为Δv∞，以使调整后的探测器按预定的要求到达火星附近（如轨道倾角为90°，近火点高度为250km），但该关系为非线性，利用微分修正法方法寻找合适的Δv∞较难，有时甚至不能收敛。

研究发现用迭代算法修正B平面上的误差时收敛性非常好，易计算获得修正误差所需的速度增量，因此被广泛用于行星际轨道的制导和精确设计。B平面是通过目标星中心且垂直于双曲轨道入射渐近线矢量S＾的平面（如图5所示），在B平面上建立一平面坐标系，其原点O为目标星质心；N＾为行星赤道面法向或黄道面发向，则B平面上两坐标轴单位矢量为

式中：b＝‖B‖；θ为双曲轨道入射渐近线与火星近焦点坐标系Xp间的夹角；iM，ωM分别为探测器轨道相对于火星赤道的轨道倾角和近火点角距。

图5 B平面Fig.5 B－plane

则由式（25）、（26），B平面参数可表示为

当iM＝90°（目标轨道为火星极地轨道）时，可得特殊形式

由上建立了B平面参数和轨道状态参数间的关系，在地火转移轨道的终端约束条件（一般是轨道倾角和近火点高度）给定后，可将其转换为B平面参数，再用不同搜索算法进行计算，直至B平面参数得到满足，此时的转移轨道便为期望轨道。

2 仿真

基于本文方法求得的转移轨道参数，用STK软件的轨道机动模块astrogator工具设计一完整的火星探测器轨道：探测器在任意时刻的受力模型为真实力模型，即考虑各种星体的引力摄动，并在转移轨道段增加两次中途修正。采用的星历为JPL的DE405星表，中途修正方法为B平面法，目标轨道为近火点高度250km，轨道周期50h的大椭圆极地轨道。

选择发射日期为2013年12月15日，到达日期为2014年10月20日，用兰伯特定理可得地火转移轨道 的 初 始 参 数 为：C3＝8.955 49km2／s2，α＝203.84°，δ＝25.006 2°。假设探测器从西昌卫星发射中心（东经102.0°，北纬28.3°）发射，先发射至高度300km的圆停泊轨道，然后滑行至合适的地点点火加速进入地球逃逸双曲轨道。双曲轨道的参数满足地火转移轨道的初始参数，发射方位角98°。地火转移阶段进行2次轨道修正。第一次在探测器刚离开地球影响球时进行，第二次在探测器刚进入火星影响球时进行。在火星轨道捕获阶段，当探测器运行到近火点为250km时，给其一速度增量进行制动，使探测器进入一周期为50h的大椭圆极火轨道。仿真所得参数见表1，地球逃逸轨道、地火转移轨道和火星捕获轨道如图6所示。

3 结束语

本文对火星探测器轨道的设计进行了研究。先求解了基于二体假设的兰伯特问题，得到了能量较优的发射窗口，再讨论了基于B平面法的中途轨道修正，并用STK软件模拟了火星探测器从我国西昌卫星发射中心发射、中途轨道修正并被火星捕获的全过程。仿真结果可知：圆锥曲线拼接法是行星际轨道初步设计的基本方法，该法获得的参数满足行星际轨道初步设计的精度要求，可作为行星际轨道精确设计的良好初值；用B平面法作为中途轨道修正方法，可使搜索算法有良好的收敛性；组合使用基于普适变量法的兰伯特问题数值求解和STK仿真软件，可有效获得真实力模型下的行星际轨道设计。本文的火星探测器轨道设计未考虑中途轨道修正策略，尚有待深入研究。

表1 仿真所得参数Tab.1 Parameter gained by simulation

图6 地球逃逸轨道、地火转移轨道和火星捕获轨道仿真结果Fig.6 Simulation results of earth departure orbit，earth－Mars transfer trajectory and mars capture orbit

［1］ GAUSS D F.Theory of the motion of the heavenly bodies moving about the sun in conic sections［M］.New York：Dover Publications，1963.

［2］ BATTIN R H.An introduction to the mathematics and method of astrodynamics［M］.New York：AIAA Education Serious，1987.

［3］ BATTIN R H，VAUGHAN R M.An elegant lambert algorithm［J］.Journal of Guidance，Control and Dynamics，1984，7：662－670.

［4］ BATE R R，MUELLER D，WHITE J E.Fundamentals of astrodynamics［M］.New York：Dover Publications，1971.

［5］ BOND V，MODERN A.Astrodynamics：fundmentals and perturbation methods［M］.Princeton University Press，1996.

［6］ CURTIS H D.Orbital mechanics for engineering students［M］. Elsevier： Butterworth－Heinemann，2005.