为充分发挥新一代长征运载火箭^[1]氢氧入轨级动力系统的高比冲优势, 中国首次探测火星工程(简称探火工程)^[2]中, 采用了长征运载火箭在近地停泊轨道滑行至预定位置后, 二次点火工作将探测器直接送入地火转移轨道的发射方式, 显著提高了探测器有效仪器设备的装载量, 但同时也明显增加了运载火箭的发射难度, 如需要大型低温长征运载火箭在发射日实现相应发射轨道的零窗口准时发射等。而由于液氢约-253℃的严酷低温工作条件限制, 使得其点火工作条件保障困难, 主要包括地面起飞时预冷好的点火条件, 以及在停泊轨道上较长时间滑行后满足二次点火工作的条件等^[3-4]。国际上, 虽然美国已经实现了低温入轨级具备近地滑行一圈后再点火的能力, 但即便其成熟火箭Delta IV(采用氢氧推进剂), 在近年来也多有媒体报道其地面点火发射一再推迟。而为了满足低温氢氧入轨级的严酷二次点火工作要求, 长征运载火箭在当前采用了较为保守的设计处理措施^[3-4], 因此需要严格限制最长滑行时间不能超过设计门限。为了能够在每26个月才能有一组较合适奔火任务的发射机会中实现及时发射, 国际上普遍采用了持续2~3周的发射日期窗口与每天至少30 min的发射窗口, 以及相应的发射轨道设计以保证深空探测任务的顺利实施^[5-6]。新研制的大型长征运载火箭在滑行时间严格受限下, 再兼顾测控布站限制等, 完成持续2~3周发射日期窗口的探火发射轨道设计尤为艰难, 本文即研究如何解决这一问题。

不考虑滑行时间严格受限和发射场航落区严格约束下, 地火转移轨道的设计已有较成熟的理论^[7], 通常根据Lambert求解逃逸速度并采用圆锥曲线拼接法完成初值设计, 再借助B平面矢量法^[8]与微分修正完成行星际转移轨道设计。B平面矢量法由于能够使计算初值快速收敛而得以普遍应用。国内外此方面的文献较多, 如文献[9]结合B平面矢量法给出了地球到火星转移轨道的近地发射段双曲线轨道的几何约束分析, 但没有长征运载火箭需要关注的停泊轨道最长滑行时间等具体约束限制下的深入分析研究等。文献[10]分析了2009—2024年地球到火星间往返发射窗口机会及发射能量C₃与逃逸赤纬等参数的数值列表及扼要分析等, 但同样不考虑运载火箭停泊轨道的滑行时间等约束的限制。文献[11]主要以运载火箭发射重量及C₃为约束, 分析了探测器轨道设计, 指出火星探测器全部的双曲线出发轨道近地点包络为一沿地球逃逸速度矢量V_E∞为轴、以r₀sin β为半径的圆轨迹(r₀为近地点高度, β为近地点拱线与双曲渐近线夹角)。文献[12]根据探测器火星目标轨道高度、倾角及天顶距约束分析了地球出发轨道速度、近地点幅角与赤经等参数设计需求, 然而发射场位置、射向及运载火箭滑行时间等约束对发射轨道的限制关系并未涉及。文献[13]将B平面矢量法推广至制导策略设计中。文献[14]分析了火星探测转移轨道的偏差传播, 指出从近地点到近火点的位置误差放大约100万倍、速度误差放大约10万倍, 实际飞行中途修正不可缺省。此外, 一般的近地轨道卫星发射任务, 只要求运载火箭的起飞时刻、星箭分离点的半长轴、偏心率、轨道倾角和近地点幅角等5个参数满足要求即可以完成发射轨道设计, 而火星等深空探测发射任务的最终目标点是近火点目标轨道, 因此除要求运载火箭以上参数严格满足要求外, 还需要星箭分离点的升交点经度与真近点角等2个参数也必须满足探测器要求, 而这2个参数与实际飞行路径直接关联, 因此如果按照传统的轨道界面分工只能采用双方反复拼接的设计方法, 计算效率与精度等显然难以满足当前设计需求。实际上, 复杂、多约束的运载火箭发射轨道设计对探测任务能否顺利实施影响显著。例如, 2003年欧洲航天局采用联盟号及Fregat上面级火箭发射火星快车探测器时, 即遭遇到上述问题^[15]:起初仅有10天的预期发射窗口及1 120 kg的运载能力, 发射风险较高, 后来经过联盟号火箭优化停泊轨道倾角及二次点火条件, 以及优化调整了整流罩与火箭三子级落区, 最终将运载能力增加到1 200 kg, 发射日期窗口拓展到了30天, 并且优化后简化了控制系统装订参数——采用9套飞行程序满足了持续30天的发射日期窗口, 仅需要探测器在飞行第10天时深空机动10 m/s, 这些措施有力保障了探测工程的顺利实施。美国DeltaⅡ 7925型火箭发射探测火星等深空任务时, 通过固定93°和99°两个射向保障了航落区安全, 简化了控制系统装订参数, 优化确定了连续多日的发射日期窗口^[5-6]。新一代长征运载火箭从海南发射火星探测器时, 由于低温入轨级2次起动间最长滑行时间严格受限^[3-4], 使入轨点位置(即轨道近地点幅角)的调节范围明显减小, 为满足相应发射年份下2~3周的有效发射日期窗口, 分析发现不排除结合探测器深空速度脉冲机动以完成发射轨道总体优化^[15-16], 故当前文献资料中的设计算法难以直接满足长征运载火箭发射轨道深入优化的需要。依此背景, 本文重点研究了长征运载火箭在当前滑行时间、发射场射向及落区等多约束限制下, 发射探测器行星际转移轨道的优化设计算法。该算法直接采用探测器近火点目标轨道参数, 并结合近地出发轨道诸多设计约束开展设计, 在高精度力学计算模型基础上, 提出双向微分修正算法, 从而也一并解决了在转移轨道的深空对接点处增加速度脉冲机动, 以进一步完成发射轨道优化的问题。

1.   长征运载火箭海南发射场约束

工程中从飞行安全角度考虑, 一般允许的射向A₀范围约为90°~115°。

注意到, 地球自转(赤道上转速为465.1 m/s)对实现长征运载火箭目标入轨速度的贡献V_launch为

式中:ϕ_L为发射场纬度。

依据齐奥尔科夫斯基公式^[17], 动力比冲I_sp, 火箭起始质量m₀, 燃烧结束质量m_k产生的速度增量ΔV为

对纬度约20°的发射场, 如果射向从90°增加到115°, 按照氢氧发动机比冲约4 300 m/s, 结合探测器与火箭末级入轨质量的大致相对关系^[18], 易知发射火星探测器等深空任务的运载能力损失约3%, 故可能情况下需要优先选择东射向发射。

2.   火星探测任务发射轨道设计

2.1   发射约束下初始轨道搜索

2.1.1   Lambert问题求解

忽略地球和火星的引力及实际尺寸, 认为转移轨道直接从地球(在日心坐标系下的位置)出发, 到达火星(在日心坐标系下的位置), 如图 1(a)所示, 该问题为典型的日心Lambert问题:

图  1  从地球转移至火星的飞行轨道示意图

Figure  1.  Schematic diagram of Earth-to-Mars flight trajectory

下载: 全尺寸图片 幻灯片

式中:t_M和t_E分别为抵达和出发时刻; μ为日心天体引力常数; r_M和r_E分别为抵达和出发时刻地球与火星距日心距离; c为地球与火星的直线距离; a为转移轨道半长轴。

关于Lambert问题的求解已有非常成熟的算法, 本文不再赘述, 推荐采用Izzo的计算方法^[19], 如果轨道近似为抛物线, 推荐采用Gooding的算法以避免Izzo算法中的奇点问题^[20-21]。图 1(a)中, V_d和V_a分别为Lambert转移轨道影响球处出发和抵达速度矢量, V_E和V_M分别为地球和火星的速度矢量, V₀和V₁分别为Lambert转移轨道相对于地球与火星影响球的速度矢量。通过V₀和V₁可以进一步设计地球影响球(SOI)内的双曲线飞行轨迹及火星影响球内的双曲线飞行轨迹, 从而将整个转移轨道由地球影响球内的双曲线弧段+日心Lambert弧段+火星影响球内双曲线弧段构成(见图 1(b)), 此即经典圆锥曲线拼接模型。

2.1.2   出入影响球的双曲线轨道初值设计

本文着重介绍地球影响球内双曲线弧段设计方法, 火星影响球内双曲线弧段设计方法类似。有别于传统的借助B平面矢量法, 本文算法直接根据近地端发射场的约束及火星端目标轨道的约束设计双曲线弧段。通过日心Lambert问题求解, 得探测器相对于地球影响球边界速度矢量V_o(地球与火星端分别对应图 1中V₀和V₁), 地球段采用下述计算单位:

式中:a_E为地球赤道半径; m_E为地球质量; []表示括号内变量在不同区域飞行段是不同参量, 此即指无量纲化, 如地球段、日心飞行段、火星飞行段等。

记地球影响球半径为r_o, v_o=|V_o|, 从式(5)得双曲线半长轴绝对值a_o:

设定近地点高度h_pE, 从式(6)得偏心率e_o:

双曲线轨道倾角i_o与发射场射向A_L及发射场纬度ϕ_L关系为

从式(8)中, 得双曲线轨道近地点幅角ω_o和升交点赤经Ω_o为

式中:; ; μ_E为地球引力常数; 为速度矢量V_o的分量; E_o为地球影响球半径r_o处偏近点角, 即

结合式(10)得对应的真近点角f_o:

式(8b)一般同时求得2个近地点幅角ω_o, 代入式(8a)分别得对应的升交点赤经Ω_o。根据式(11)计算双曲线轨道纬度幅角u_o:

为便于区别, 工程上通常以u_o∈(－π/2, π/2)作为升交出发, 以u_o∈(π/2, 3π/2)作为降交出发。

由式(5)~式(11)可以给出地球至火星转移轨道在脱离地球影响球处的Kepler轨道根数a_o、e_o、i_o、Ω_o、ω_o、f_o, 而星箭分离点的地心轨道根数与该轨道只是真近点角f_I不同, f_I由火箭发射段确定。研究表明, 不同射向及滑行时间对应的f_I变化不大, 可近似认为是常值。

综上, 在地球影响球内, 以地球和探测器组成二体问题动力学模型, 根据日心Lambert弧段给出的地球影响球边界处的速度矢量V_o, 可得J2000系下双曲线出发轨道星箭分离点轨道根数a_o、e_o、i_o、Ω_o、ω_o、f_I, 并由之得位置矢量r_J2000和速度矢量_J2000, 再通过式(12)坐标转换^[22]得星箭分离时刻与地固CGCS2000系瞬时重合(之后保持惯性空间不动)的地心赤道惯性系(以下标G表示)下位置r_G和速度_G:

式中:Q(t)为岁差章动矩阵; R(t)为地球转角矩阵; W(t)为极移矩阵。由r_G、_G再得相应的轨道根数需求:a_o、e_o、i_oE、Ω_oE、ω_oE、f_I。

抵达近火点双曲线轨道计算在火星J2000平赤道坐标系(用下标ME表示)下描述, 该坐标系与J2000火心坐标系(用下标MJ表示)转换关系^[22]见式(13), IAU年度工作报告对火星北极指向α₀、δ₀的精确数据每三年互联网发布一次。

式中:矩阵M下标x, y, z分别表示绕相应坐标轴的转换矩阵^{[17, 22]}。火星影响球内的双曲线弧段设计方法与地球影响球的双曲弧段设计方法类似但不是本文关注的重点, 所以火星端的双曲线轨道设计这里不再赘述。

2.1.3   长征运载火箭滑行时间及发射场设计约束

长征运载火箭发射轨道星下点轨迹几何示意见图 2, 采用圆停泊轨道, 忽略各摄动因素, 得

图  2  探测火星地球出发轨道的几何图示

Figure  2.  Earth-to-Mars launch trajectory geometry illustration

下载: 全尺寸图片 幻灯片

设θ₁为从长征运载火箭起飞到进入停泊轨道前的工作弧段, θ₂为停泊轨道滑行结束至探测器分离的弧段(分析表明θ₁、θ₂可以近似为常值), 其他符号含义同前, 据各参量物理含义, 易得所需长征运载火箭入轨级2次起动间滑行时间T_h为

图 2中, 表示春分点, α_o和δ_o分别为地球至火星转移轨道在地球影响球处的逃逸赤经和赤纬, 其他符号含义均在上下文中有明确说明。根据长征运载火箭发射探测器入轨时刻升交点赤经Ω及起飞时刻发射场恒星时角α_L定义, 结合图 2球面三角关系易得

式中:λ_L为发射场经度; Ω_e为探测器分离点时刻发射轨道计算的升交点经度; ω_e为地球自转角速度; t_f0为火箭滑行时间T_h=0 s时从起飞到探测器分离的时间。对选定的射向A_L, 预先选择很小的阈值ε_i> 0(i=1, 2), 调整发射当日出发时刻t_E使:

易得满足长征运载火箭最长滑行时间限制的飞行轨道参数及起飞时刻初解。

2.1.4   射向对飞行轨道及航程角的影响

长征运载火箭关注的发射日出发轨道及航程角与射向的初步关系分析如下。结合图 2, 得飞行航程角θ(θ=ω_oE+f_o－d)与地球影响球处逃逸赤纬δ_o、射向A_L及发射场纬度ϕ_L关系:

设ϕ_L=20°, 可得图 3。最长滑行时间限定了近地点幅角ω_oE与飞行航程角θ的变化范围, 每一组航程角对应一个起飞时刻。

图  3  射向A_L、飞行航程角θ与δ_o关系

Figure  3.  Contour map of δ_o with respect to launch azimuth A_L and flight range angle θ

下载: 全尺寸图片 幻灯片

结合图 3, 在射向范围0°~180°情况下, 结合发射轨道基础理论, 分析易知:

1) 当|δ_o|＜ϕ_L(如δ_o=0°)时, 同一射向下有航程角θ相差约180°的2组解, 2组解的滑行时间需求和起飞时刻差异较明显。

2) 当δ_o=ϕ_L, 仅有一组射向随航程角θ减小而增大的解; 当δ_o=－ϕ_L, 除了一组射向随航程角θ减小而增大的解外, 还有航程角恰好为180°、射向任意的另外一组解。

3) 当|δ_o|>ϕ_L, 同一射向下有航程角θ比较接近的2组解, 射向越大两航程角越接近, 2组解的滑行时间需求和起飞时刻较接近。

2.1.5   结合Pork-chop图完成初始发射轨道优选

采用上述加入滑行时间限制再精确双向微分修正的设计算法可以较快速完成初步发射轨道搜索, 通过对近地点与近火点速度增量辅以Pork-chop图^[7]分析, 即可以初步完成多设计约束下, 长征运载火箭发射火星探测轨道的初步运载能力确定、连续发射日期窗口初选和初始发射轨道设计。限于篇幅, 此部分内容从略。

2.2   精确设计发射轨道

2.2.1   精确动力学方程

在起飞时刻地心惯性系下, 长征运载火箭飞行的质心动力学矢量方程^{[17, 23]}为

式中:m为火箭质量; r_G0为惯性系原点到火箭质心矢径; P为火箭推力; R_f为作用在火箭上的气动力; g为作用在火箭上的地球引力加速度。通过式(20b)转换得入轨时刻地心惯性系下位置和速度:

进而得相应的发射轨道根数:a_o、e_o、i_oE、Ω_e、ω_oE、f_I。

从探测器分离开始, 设定整个地火转移轨道都是在日心J2000坐标系下完成积分, 动力学方程为

式中:μ_S为太阳引力系数; j∈⊙表示太阳系八大行星及月球等; μ_j为太阳系下扰动天体引力系数, 对于地火转移轨道设计, 需要考虑八大行星及月球的影响; F^E_ns和F^M_ns分别为地球和火星的非球形摄动力^[24-25]; δ_E、δ_M为开关函数, 当探测器在地球(火星)影响球界内取1, 其余取0;F_others指地球辐射压、高层空间大气、火星稀薄大气、小行星引力等; F_SRP为太阳辐射压, 设定探测器面质比为S/m, S为表面积, m_s为质量, 忽略行星阴影区影响, 则

其中:κ为探测器反射系数; c₀为光速; ρ_⊙为距离r_⊙处的太阳辐射流。

采用RKF7(8)阶龙格库塔算法, 分析计算表明, 探测器在靠近地球或火星飞行段均采用10×10阶非球形引力项, 比两者均采用J₂项时设计出的转移轨道在近地端的速度差异小于0.1 m/s、位置差异小于30 m, 在近火端的速度差异约0.3 m/s、位置差异约140 m; 光压对探测器的摄动影响量一般可能超过近地端高阶摄动影响而小于火星端高阶项影响; 注意到光压影响与探测器具体形态关系较密切等因素, 故综合后下面以考虑大行星及月球等中心天体引力、地球与火星飞行段均采用10×10阶次非球形引力项给出分析算例。

2.2.2   双向微分修正算法实现精确对接

记星箭分离点时刻(记为t₀)实际飞行轨道与参考轨道的状态量(即位置和速度)偏差为ΔX(t₀), 小偏差条件下, t时刻状态量偏差ΔX(t)近似满足:

式中:Φ(t, t₀)为参考轨道t₀时刻至t时刻的状态转移矩阵^[7]。由于从近地点到近火点转移轨道的误差传递特性急剧放大, 因此直接采用近火星目标轨道根数求解近地点出发轨道极其困难, 一般需要转化处理。

若分别从近地点、近火点向地球至火星转移轨道中间的某个时间点t_m积分, 记地球出发端的状态量偏差为ΔX(t₀), 火星到达端的状态量偏差为ΔX(t_f), 类似式(21), 有

式中:ΔX_E(t_m)为从地球端正向积分至t_m时刻轨道的状态量偏差; ΔX_M(t_m)为从火星端逆向积分至t_m时刻轨道的状态量偏差。对同一条转移轨道, 要求这2个偏差相同, 即

再将式(21)中的t取为t_f, 并结合式(22)易得

即时刻t₀~t_f的状态转移矩阵Φ(t_f, t₀)可分解为2个分弧段的状态转移矩阵Φ_E(t_m, t₀)和Φ_M(t_m, t_f); 由于将状态转移矩阵Φ(t, t₀)一分为二后, 两积分段的误差传递特性急剧降低, 即可以较顺利地实现问题求解。实际上B平面矢量法的本质原理亦是如此, 只是固定为目标天体影响球边界处的B平面参数。分析表明, 由于日心转移弧段的良好解析性, 使中间积分结合点可以自由移动而不影响计算结果, 所以就解决了后续在结合点处增加深空机动以进一步优化长征火箭发射轨道的问题。下面对算法进一步介绍。

以r_J2000和_J2000代表在地心J2000坐标系下地球出发时探测器与长征运载火箭分离点的位置和速度, 以r_MJ和_MJ代表探测器火心J2000坐标系下抵达近火点的位置和速度, 组成一个12维的状态矢量Y, 即

式中:符号R表示相对日心系的位置, 下标d表示地球影响球界出发点, a表示火星影响球界抵达点。记下标E表示探测器分离点时刻地球的日心系位置和速度, 下标M表示抵达近火点时刻火星的日心位置和速度, 则

在地球出发点, 存在3个约束:

式中: h_pE、i_E和f_E分别为长征运载火箭发射探测火星轨道分离时近地点高度、倾角和真近点角约束, 此3个约束结合式(25)和式(12)易得。

在近火双曲线轨道段, 也存在3个类似约束:

式中: h_pM、i_M和f_M分别为探测火星目标轨道近火点高度、倾角和真近点角约束, 此3个约束结合式(26)和式(13)易得。

当将R_d、_d正向积分, 并且反向积分R_a、_a, 使得2条积分轨道空间参数匹配, 对应6个约束, 记出发历元为T₀, 抵达历元为T_f, 正反双向积分的相遇历元时刻为T_m, 则

式中:上标“+”表示正向积分轨道, 上标“-”表示反向积分轨道。基于长征运载火箭发射轨道优化中, 后续关注的深空速度脉冲机动优化问题, 分析表明, T_m可以从地球影响球边界至火星影响球边界之间自由选择, 相应根据深空对接点处速度脉冲增量为g_10~12赋值^[16]; 对于无需深空脉冲速度机动的情况, 一般可取T_m=(T₁+T₂)/2。

将式(27)~式(29)统一记为

此时迭代方程为

式中:

可将雅可比矩阵分解为如下形式:

式中:

12项约束对应12个变量, 通过迭代求解即完成火星目标轨道的地球出发轨道设计。

2.3   计算流程

考虑发射场约束, 分析与精确设计长征运载火箭发射火星探测器转移轨道的设计流程见图 4。

图  4  长征运载火箭设计探测火星轨道示意图

Figure  4.  Schematic diagram of Long March launch vehicle for Earth-to-Mars trajectory design

下载: 全尺寸图片 幻灯片

3.   分析及算例

计算采用的火星探测器目标轨道约束为:近火点高度500 km, 倾角93°; 结合初步选定的射向, 近地约束为:射向A_L=107°, 出发近地点高度200 km, 以UTC时间2020-06-19出发, 转移时间198天, 2021-01-02T23:59升交抵达火星为算例, 计算得地球影响球边界处逃逸赤纬约24°, 由图 3知有滑行时间T_h相近的2组初解, 采用双向微分修正得表 1、表 2精确解, 表中EJ2000表示地心J2000赤道系, ME表示J2000火心IAU平赤道系^[22], SJ2000代表距近火点时刻t_m=90天时日心J2000系参数, ΔV指对应地球或火星端各自圆停泊轨道的速度增量。

表  1  发射轨道T_h=1 849 s, 探测器05 : 25释放, 地球升交出发

Table  1.  Trajectory T_h=1 849 s, probe released on 05 : 25, ascending from the Earth

轨道根数	EJ2000	SJ2000	ME
a	-21 424.437 km	1.321 527 5 AU	-3 816.831 km
e	1.307 039 0	0.241 839 5	2.020 742 1
i/(°)	25.405	23.424	93.000
Ω/(°)	316.214	356.094	156.275
ω/(°)	293.745	293.291	118.136
M/(°)	2.781	57.042	0
Hp	200.0 km	1.001 930 0 AU	500.0 km
ΔV/(m·s-1)	4 039.207	0	2 446.974
注:Hp为轨道的近拱点高度(分别相对近地、椭圆近日、近火), M为平近点角。

下载: 导出CSV 

| 显示表格

表  2  发射轨道T_h=2 312 s, 探测器03 : 15释放, 地球降交出发

Table  2.  Trajectory T_h=2 312 s, probe released on 03 : 15, descending from the Earth

轨道根数	EJ2000	SJ2000	ME
a	-21 373.589 km	1.321 389 5 AU	-3 815.267 km
e	1.307 769 4	0.241 822 3	2.021 160 4
i/(°)	25.372	23.420	93.000
Ω/(°)	281.242	356.067	156.217
ω/(°)	325.518	293.283	118.191
M/(°)	2.790	57.077	0
Hp	200.0 km	1.001 848 1 AU	500.0 km
ΔV/(m·s-1)	4 041.078	0	2 447.373

下载: 导出CSV 

| 显示表格

分析表明:①t_m可以大幅移动而不影响优化设计结果; ②在结合点处速度增量可以为0 m/s或增加到超过500 m/s, 都能找到良好求解。

若限定最长滑行时间不能超过2 000 s^[3-4], 则只有表 1中解满足要求; 而如果限定最长滑行时间不超过1 100 s或更短, 则需要采用上述加入滑行时间限制再精确双向微分修正的设计算法, 通过加大射向或结合探测器深空机动等完成发射轨道优化; 具体的分析是繁琐复杂的工程问题, 已不在本文讨论范围。

将表 1、表 2中计算结果在同等条件下, 与STK仿真分析软件积分对算表明, 二者计算精度与所需时长相当, 算例中从近地点到深空对接点、从深空对接点到近火点, 位置差异约为1 km。

4.   结论

本文给出了长征运载火箭允许滑行时间受限和发射场可用射向范围约束下, 长征运载火箭发射火星探测任务从地面起飞至近火点轨道设计的完整模型算法。

1) 解决了适用于工程运算所需精度的分析研究模型问题。

2) 解决了火星探测任务的发射轨道优化设计问题, 确定了运载能力、优选了发射日期窗口, 完成了航落区优化选择等。

双向微分修正算法具有结合点可以大幅移动, 以便于后续分析中探测器采用深空机动后, 深入开展长征运载火箭发射轨道再优化分析的优势。工程应用表明此方法稳定、可靠好用。除了可以用于地火转移发射轨道设计外, 本文算法还可用于长征运载火箭发射其他深空探测器的发射轨道设计。