3. 二体轨道计算实战
本节要打通二体问题中最常用的一条计算闭环：给定某时刻的位置矢量 r 与速度矢量 v，先求出六个轨道根数；再由轨道根数和目标时刻，反算出该时刻的位置与速度。这条闭环是轨道确定、轨道预报、交会分析与任务规划的基础。为避免歧义，本节统一采用地心惯性系（ECI）：+X 指向春分点，+Z 指向北极，+Y 按右手系补齐；角度单位在公式中用弧度，在界面展示时换算为度；引力参数取 μ = 398600.4418 km³/s²。
一、从 r、v 求轨道根数的计算链
六个根数中，a（半长轴）、e（偏心率）、i（轨道倾角）、Ω（升交点赤经）、ω（近地点幅角）、ν（真近点角）各有明确的几何来源，计算顺序不能随意打乱。核心中间量是角动量矢量 h = r × v、节点矢量 n = ẑ × h 和偏心率矢量 e_vec = ((v² − μ/r)·r − (r·v)·v)/μ。
具体计算步骤如下：
- 能量与半长轴：比机械能
ε = v²/2 − μ/r，则a = −μ/(2ε)。若ε ≥ 0，说明是抛物线或双曲线轨道，a为负或无穷，需改用p = h²/μ描述。 - 偏心率：
e = |e_vec|。当e接近 0 时，e_vec方向数值不稳定，近地点幅角与真近点角应改用其他方式处理。 - 倾角：
i = arccos(h_z/|h|)，取值[0, π]。 - 升交点赤经：
Ω = arccos(n_x/|n|)，若n_y < 0则取2π − Ω；当|n| ≈ 0（赤道轨道）时Ω无定义，可约定为 0。 - 近地点幅角：
ω = arccos((n·e_vec)/(|n||e_vec|))，若e_vec,z < 0则取2π − ω；赤道轨道时改用经度代替。 - 真近点角：
ν = arccos((e_vec·r)/(|e_vec||r|))，若r·v < 0（正在向近地点运动）则取2π − ν。
下面这段 Python 代码把上述步骤封装为函数，并处理了赤道轨道与圆轨道的退化情形：
import numpy as np
MU = 398600.4418 # km^3/s^2
def rv_to_elements(r, v):
r = np.asarray(r, float); v = np.asarray(v, float)
rn = np.linalg.norm(r); vn = np.linalg.norm(v)
h = np.cross(r, v); hn = np.linalg.norm(h)
n = np.cross([0, 0, 1.0], h); nn = np.linalg.norm(n)
e_vec = ((vn**2 - MU/rn)*r - np.dot(r, v)*v) / MU
e = np.linalg.norm(e_vec)
energy = vn**2/2 - MU/rn
a = -MU/(2*energy) if abs(energy) > 1e-12 else np.inf
i = np.arccos(np.clip(h[2]/hn, -1, 1))
if nn < 1e-10: # 赤道轨道
Omega = 0.0
argp = np.arctan2(e_vec[1], e_vec[0]) if e > 1e-10 else 0.0
else:
Omega = np.arccos(np.clip(n[0]/nn, -1, 1))
if n[1] < 0: Omega = 2*np.pi - Omega
argp = np.arccos(np.clip(np.dot(n, e_vec)/(nn*e), -1, 1))
if e_vec[2] < 0: argp = 2*np.pi - argp
if e > 1e-10:
nu = np.arccos(np.clip(np.dot(e_vec, r)/(e*rn), -1, 1))
if np.dot(r, v) < 0: nu = 2*np.pi - nu
else:
nu = np.arctan2(np.dot(np.cross(n, r), h)/(nn*hn*rn),
np.dot(n, r)/(nn*rn))
return dict(a=a, e=e, i=i, Omega=Omega, argp=argp, nu=nu)二、由根数生成任意时刻的位置与速度
反算的关键是先把真近点角推进到目标时刻。对椭圆轨道，用偏近点角 E 作为桥梁：tan(E/2) = sqrt((1−e)/(1+e))·tan(ν/2)，平近点角 M = E − e·sinE，平均角速度 n = sqrt(μ/a³)，则 M(t) = M₀ + n·(t − t₀)。由 M 解 E 用牛顿迭代：E ← E − (E − e·sinE − M)/(1 − e·cosE)，初值取 E₀ = M（小偏心率）或 E₀ = π（大偏心率），通常 5～8 次迭代即可收敛到 1e-12。
得到 E 后，在近焦点坐标系（PQW）中写出位置与速度，再经三次旋转 Rz(Ω)·Rx(i)·Rz(ω) 变换到 ECI：
def elements_to_rv(a, e, i, Omega, argp, nu, MU=MU):
p = a*(1 - e**2)
r_pqw = np.array([p*np.cos(nu)/(1+e*np.cos(nu)),
p*np.sin(nu)/(1+e*np.cos(nu)), 0.0])
v_pqw = np.sqrt(MU/p)*np.array([-np.sin(nu), e+np.cos(nu), 0.0])
cO, sO = np.cos(Omega), np.sin(Omega)
ci, si = np.cos(i), np.sin(i)
cw, sw = np.cos(argp), np.sin(argp)
R = np.array([
[cO*cw - sO*sw*ci, -cO*sw - sO*cw*ci, sO*si],
[sO*cw + cO*sw*ci, -sO*sw + cO*cw*ci, -cO*si],
[sw*si, cw*si, ci]])
return R @ r_pqw, R @ v_pqw
def propagate(a, e, i, Omega, argp, nu0, dt, MU=MU):
E0 = 2*np.arctan2(np.sqrt(1-e)*np.sin(nu0/2),
np.sqrt(1+e)*np.cos(nu0/2))
M0 = E0 - e*np.sin(E0)
M = M0 + np.sqrt(MU/a**3)*dt
E = M if e < 0.8 else np.pi
for _ in range(50):
dE = (E - e*np.sin(E) - M)/(1 - e*np.cos(E))
E -= dE
if abs(dE) < 1e-12: break
nu = 2*np.arctan2(np.sqrt(1+e)*np.sin(E/2),
np.sqrt(1-e)*np.cos(E/2))
return elements_to_rv(a, e, i, Omega, argp, nu)三、算例与闭环验证
取一组教学示例初值（合理假设，非实测）：r₀ = [7000, 0, 0] km，v₀ = [0, 7.5, 1.0] km/s。代入上述函数可得近似结果：a ≈ 6800 km 量级、e ≈ 0.05 量级、i ≈ 7.6°、Ω ≈ 0°、ω 与 ν 由几何关系确定。验证闭环时，最直接的做法是“往返一致性检验”：用 rv_to_elements 得到根数后，立即用 elements_to_rv 在 dt = 0 处重建状态，比较重建的 r、v 与原始输入之差，位置误差应在 1e-9 km 量级、速度误差在 1e-12 km/s 量级；若误差明显偏大，多半是角度象限判断或旋转矩阵顺序写错。
第二步验证是“物理量守恒”：把状态推进一个轨道周期 T = 2π·sqrt(a³/μ)，位置速度应回到初值附近；推进半个周期，位置应接近远地点。第三步验证是“能量与角动量不变”：任意时刻计算 ε 与 |h|，与初值相比相对偏差应小于 1e-10。这三类检验分别覆盖了几何一致性、周期性正确性和动力学守恒性，是工程实现中最常用的自检手段。
四、实现中最容易出错的几个点
第一是象限判断。反余弦只能给出 [0, π]，Ω、ω、ν 都必须结合矢量分量的符号或 r·v 的符号补到 [0, 2π)，否则轨道会被“镜像”到错误象限。第二是退化情形：赤道轨道 n ≈ 0 时 Ω 无定义，圆轨道 e ≈ 0 时 ω 无定义，代码中必须显式分支，不能直接除零。第三是单位与角度制：μ 用 km³/s² 时长度必须是 km、时间必须是 s，混用 m 或分钟会导致 a 差几个数量级。第四是牛顿迭代的初值：e > 0.8 时若仍取 E₀ = M，在高偏心率近地点附近可能不收敛，应改用 E₀ = π 或先做一次粗略修正。
把上述流程串起来，就得到一条可复用的计算管线：输入状态矢量 → 求 h、n、e_vec → 输出六个根数 → 由目标时刻推进 M、解 E、得 ν → 旋转回 ECI 输出 r(t)、v(t) → 用往返一致性、周期性与守恒量三重检验确认结果可信。掌握这条闭环，后续的轨道机动、交会窗口与摄动修正才有可靠的计算底座。