2. 二体问题与开普勒轨道

轨道预报的数学起点是二体问题。把地球视为质量集中在中心的质点,航天器视为另一个质点,忽略大气阻力、非球形引力、日月摄动等次要力,则航天器相对地心的运动由牛顿万有引力定律与牛顿第二定律共同决定。设航天器位置矢量为 r,地心引力常数为 μ=GM(地球取 μ≈3.986004418×10¹⁴ m³/s²),则二体运动方程为:

d²r/dt² = -μ·r / |r|³

该方程是一个二阶非线性矢量微分方程。它的解析解表明:轨道是圆锥曲线,地心位于其中一个焦点上。这一结论正是开普勒第一定律的动力学来源。同时,由角动量守恒可推出面积速度不变(第二定律),由能量守恒与半长轴关系可推出周期平方正比于半长轴立方(第三定律)。因此,二体问题不是经验拟合,而是有严格解析解的动力学系统,这也是它能作为轨道预报基准的原因。

六个轨道根数:把三维轨迹压缩成六个数字

运动方程的解包含六个积分常数,工程上习惯用一组具有明确几何意义的参数——轨道六根数——来描述:

  • 半长轴 a:决定轨道大小与周期,T=2π√(a³/μ)。
  • 偏心率 e:决定轨道形状,e=0 为圆,0<e<1 为椭圆。
  • 轨道倾角 i:轨道面与赤道面的夹角,0°≤i≤180°。
  • 升交点赤经 Ω:从春分点方向量到升交点方向的赤经角。
  • 近地点幅角 ω:在轨道面内,从升交点量到近地点方向的角度。
  • 真近点角 f(或平近点角 M、偏近点角 E):描述航天器当前在轨道上的位置。

前五个根数确定轨道在惯性空间中的几何形状与朝向,第六个根数确定航天器在轨道上的瞬时位置。需要特别强调参照约定:本文所有角度均在地心赤道惯性系(ECI)中定义,春分点方向为 X 轴,赤道面为 XY 平面,Z 轴指向北极。方向与正负号均以此为准,后续图示也保持一致。

位置速度互算:预报的核心运算

轨道预报在计算层面就是两类互算:由六根数求某时刻的位置速度(正向预报),以及由位置速度反求六根数(定轨结果转换)。

由根数到位置速度的典型流程是:先由平近点角 M 通过开普勒方程 M=E-e·sinE 解出偏近点角 E(迭代法,如牛顿迭代),再由 E 求轨道面内的位置分量:

x_op = a(cosE - e)
y_op = a√(1-e²)·sinE
r    = a(1 - e·cosE)

速度分量可由对时间求导得到,或利用活力公式与角动量关系计算。随后通过三次旋转(近地点幅角 ω、倾角 i、升交点赤经 Ω)把轨道面坐标转换到 ECI 系,得到 r 与 v。反向计算则先由 r、v 求比角动量 h=r×v 与偏心率矢量 e,再依次解出 i、Ω、ω、a、e 和真近点角,最后换算为平近点角。

开普勒方程 M=E-e·sinE 是超越方程,没有初等闭式解。工程上常用牛顿迭代:Eₖ₊₁=Eₖ-(Eₖ-e·sinEₖ-M)/(1-e·cosEₖ),初值可取 E₀=M,通常 3~5 次迭代即可收敛到 1e-10 rad 量级。

近地圆轨道与椭圆轨道的解析预报练习

下面给出两个可直接手算或编程验证的示例。参数均为教学假设值,用于说明方法,不代表任何实际任务数据。

练习一:近地圆轨道。设 a=7000 km,e=0,i=98°,Ω=30°,ω=0°,取 t=0 时 M=0°。由于 e=0,开普勒方程退化为 E=M,真近点角等于平近点角。轨道周期 T=2π√(a³/μ)≈2π√(3.43×10¹⁴/3.986×10¹⁴)≈5828 s≈97.1 min。在 t=0 时,航天器位于近地点(对圆轨道即升交点方向),位置矢量在 ECI 中可由旋转矩阵直接算出。若取 t=1457 s(约 T/4),则 M=90°,位置到达轨道面内与初始方向成 90° 处。通过这一练习可以验证:圆轨道预报只需处理平近点角到真近点角的线性关系,重点在于旋转矩阵的正确性。

练习二:椭圆轨道。设 a=12000 km,e=0.3,i=45°,Ω=60°,ω=90°,t=0 时 M=0°。周期 T≈2π√(1.728×10¹⁵/3.986×10¹⁴)≈13080 s≈3.63 h。近地点半径 r_p=a(1-e)=8400 km,远地点半径 r_a=a(1+e)=15600 km。取 t=2000 s,先求平近点角 M=n·t,其中平均角速度 n=√(μ/a³)≈4.80×10⁻⁴ rad/s,得 M≈0.96 rad≈55°。用牛顿迭代解 E:取 E₀=M,迭代后得 E≈1.13 rad。代入位置公式得 r≈a(1-e·cosE)≈12000×(1-0.3×0.427)≈10463 km,轨道面内 x_op≈12000×(0.427-0.3)≈1524 km,y_op≈12000×√(1-0.09)×0.904≈10360 km。再经 ω、i、Ω 三次旋转得到 ECI 坐标。此练习的关键观察点是:偏近点角 E 与平近点角 M 的差异随 e 增大而显著,e=0.3 时两者相差约 10° 以上,直接令 E=M 会带来明显位置误差。

二体问题轨道预报:根数 → 位置速度 与 位置速度 → 根数 正向预报:六根数 → r, v 输入:a, e, i, Ω, ω, M₀ M = M₀ + n·(t-t₀),n = √(μ/a³) 开普勒方程迭代:M = E - e·sinE 轨道面坐标:x_op=a(cosE-e), y_op=a√(1-e²)sinE 三次旋转:Rz(-Ω)·Rx(-i)·Rz(-ω) 输出:ECI 系 r, v 反向定轨:r, v → 六根数 输入:ECI 系 r, v h = r × v,e_vec = (v×h)/μ - r/|r| i = arccos(h_z/|h|),Ω = arctan2(h_x, -h_y) ω = arccos(n·e_vec/(|n||e_vec|)),a = 1/(2/|r| - v²/μ) f = arccos(e_vec·r/(|e_vec||r|)),E → M 输出:a, e, i, Ω, ω, M 参照系:地心赤道惯性系(ECI),春分点方向为 X 轴,Z 轴指向北极

结果判读与常见误差来源

完成上述练习后,应从三个层面判读结果。第一,量级检查:近地轨道周期应在 90~100 min 量级,椭圆轨道近地点半径必须大于地球半径(约 6378 km)加大气层高度,否则轨道不可持续。第二,守恒量检查:由算出的 r、v 反算比机械能 ε=v²/2-μ/r 与比角动量 h=r×v,在同一轨道上应保持常数,偏差过大说明旋转矩阵或角度单位有误。第三,互逆检查:把正向算出的 r、v 再输入反向流程,应能还原出原始六根数,这是验证代码正确性最直接的手段。

常见误差来源包括:角度制与弧度制混用;开普勒方程迭代未收敛就停止;旋转顺序写反(正确顺序是先绕 Z 轴转 -ω,再绕 X 轴转 -i,最后绕 Z 轴转 -Ω);以及把真近点角 f 与平近点角 M 直接等同。对于 e<0.01 的近圆轨道,这种等同带来的位置误差通常可接受;但当 e>0.1 时,必须严格求解开普勒方程。掌握这两类轨道的解析预报流程,就为后续引入 J2 摄动、大气阻力等更复杂模型打下了可对比的基准。