4. 轨道预报基础：数值积分方法（RK4、Adams）

各位同学，今天我们来聊聊轨道预报里最核心的“发动机”——数值积分方法。说白了，就是怎么用计算机去算卫星下一秒、下一分钟、下一天会飞到哪里去。

我刚开始做轨道预报那会儿，总觉得这玩意儿不就是解个微分方程嘛，有啥难的？结果第一次跑出来的轨道，卫星直接“飞”到火星上去了……嗯，从那以后我再也不敢小看数值积分了。

4.1 为什么需要数值积分？

轨道动力学方程，本质上是一个二阶常微分方程组：

r'' = -μ * r / r³ + a_perturb

这个方程，除了二体问题有解析解，其他情况基本无解。你想想看，地球非球形引力、大气阻力、太阳光压、第三体引力……这些摄动项加进去，方程就变得极其复杂。

所以，我们只能靠数值积分，一步一步往前推。就像走路一样，一步迈出去，算一步的位置和速度。

4.2 龙格-库塔法（RK4）—— 最常用的“万金油”

RK4，全称四阶龙格-库塔法。我个人习惯叫它“四阶龙哥”。为什么它这么流行？因为它在精度和计算量之间取得了很好的平衡。

它的思路其实很朴素：在每个积分步内，取四个点的斜率，然后加权平均，得到一个更准确的“平均斜率”。

4.2.1 RK4的数学形式

假设我们要解一阶微分方程：

dy/dt = f(t, y)

RK4的递推公式如下：

k1 = f(t_n, y_n)
k2 = f(t_n + h/2, y_n + h*k1/2)
k3 = f(t_n + h/2, y_n + h*k2/2)
k4 = f(t_n + h, y_n + h*k3)

y_{n+1} = y_n + (h/6) * (k1 + 2*k2 + 2*k3 + k4)

这里 h 是步长。你看，它计算了四个斜率，然后按 1:2:2:1 的比例加权。为什么是这个比例？这是从泰勒展开推导出来的，目的是让局部截断误差达到 O(h⁵)。

4.2.2 RK4的代码实现

下面是一个简单的 RK4 积分器，用于轨道预报：

def rk4_step(state, t, h, derivs):
    """
    state: [x, y, z, vx, vy, vz]
    t: 当前时间
    h: 步长
    derivs: 导数函数，返回 [vx, vy, vz, ax, ay, az]
    """
    k1 = derivs(t, state)
    
    s2 = [state[i] + 0.5 * h * k1[i] for i in range(6)]
    k2 = derivs(t + 0.5*h, s2)
    
    s3 = [state[i] + 0.5 * h * k2[i] for i in range(6)]
    k3 = derivs(t + 0.5*h, s3)
    
    s4 = [state[i] + h * k3[i] for i in range(6)]
    k4 = derivs(t + h, s4)
    
    new_state = [state[i] + (h/6.0) * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) 
                 for i in range(6)]
    
    return new_state

这段代码看着简单，但实际用的时候要注意：derivs 函数里要包含所有摄动力的计算。我见过有人只算了二体引力，结果预报误差一天就超过 100 公里……

4.3 Adams 方法 —— 多步法的代表

RK4 是单步法，只用当前步的信息。而 Adams 方法属于多步法，它会利用前面几步的信息来预测下一步。

为什么要用多步法？说白了，就是“吃老本”。前面几步的信息已经算出来了，不用白不用。这样可以在同样的计算量下，获得更高的精度。

4.3.1 Adams-Bashforth 预测公式

四阶 Adams-Bashforth 公式（显式）：

y_{n+1} = y_n + (h/24) * (55*f_n - 59*f_{n-1} + 37*f_{n-2} - 9*f_{n-3})

你看，它用了前面四步的导数信息。系数 55, -59, 37, -9 是从插值多项式推导出来的。

4.3.2 Adams-Moulton 校正公式

光有预测还不够，我们还要校正。Adams-Moulton 公式（隐式）：

y_{n+1} = y_n + (h/24) * (9*f_{n+1} + 19*f_n - 5*f_{n-1} + f_{n-2})

实际应用中，我们通常把预测和校正结合起来，形成 PECE 模式：

1. P（Predict）：用 Adams-Bashforth 预测 y_{n+1}
2. E（Evaluate）：计算 f_{n+1} = f(t_{n+1}, y_{n+1})
3. C（Correct）：用 Adams-Moulton 校正 y_{n+1}
4. E（Evaluate）：重新计算 f_{n+1}

4.4 RK4 vs Adams：怎么选？

这个问题，我经常被问到。其实没有绝对的答案，要看具体场景。

对比项	RK4	Adams（PECE）
每步计算量	4 次导数计算	2 次导数计算
精度	O(h⁵)	O(h⁵)
稳定性	较好	中等（需注意步长）
启动方式	自启动	需要 RK4 启动
变步长	容易实现	较复杂
适用场景	短弧段、高精度	长弧段、效率优先

我个人习惯：

• 短弧段（1-3天）：用 RK4，简单可靠
• 长弧段（1周以上）：用 Adams，效率更高
• 变步长需求：用 RK4，或者用更高阶的 Runge-Kutta-Fehlberg（RKF45）

4.5 步长选择 —— 一个容易被忽视的问题

步长选多大？这其实是个大学问。

我给大家一个经验公式：

h = T_orbit / N

其中 T_orbit 是轨道周期，N 是每圈积分步数。对于 LEO 卫星（周期约 90 分钟）：

• N = 1000（步长约 5.4 秒）：精度较好
• N = 500（步长约 10.8 秒）：精度一般
• N = 100（步长约 54 秒）：精度较差，可能发散

4.6 知识体系总览

下面这张图，是我总结的数值积分方法在轨道预报中的应用框架：

4.7 实际应用中的注意事项

最后，我给大家总结几条实战经验：

1. 初始化要谨慎：Adams 方法的前几步必须用 RK4 启动，否则误差会很大。
2. 步长要匹配：步长不能大于轨道周期的 1/100，否则积分会发散。
3. 摄动力要全：只算二体引力的预报，基本就是“睁眼瞎”。
4. 定期检查：每积分一圈，检查一下轨道能量是否守恒。如果能量变化超过 1%，说明步长太大了。

---

公众号:蓝海资料掘金营，微信deep3321