三体轨道动力学中的稳定与不稳定流形：低能量转移的免费通道

做轨道设计的朋友，第一眼看到那些低能量转移任务的轨迹图时，多半会皱眉：好好的霍曼转移不走，非得在太阳系里绕个大圈，时间多花好几年，图画的还密密麻麻。后来入了三体轨道动力学的门才明白，那些弯弯绕绕根本不是什么“炫技”，而是让探测器顺着稳定流形与不稳定流形的天然通道，在两个主天体之间的引力场里“搭便车”。

稳定流形与不稳定流形，听起来是纯数学概念，其实是三体轨道动力学里最贴近工程实用的工具之一。它回答了轨道设计中一个很实际的问题：在不烧燃料或者少烧燃料的前提下，怎么让一个航天器靠近拉格朗日点、绕着它转，或者从它身边溜走。这篇文章把概念、数学结构、计算流程和实操坑一次讲清楚。无论你是做任务初步设计的工程师，还是在读的轨道动力学研究生，甚至只是对三体问题好奇的自学者，都能从中找到可以直接落地的计算方法。

1. 先把概念放对位置：流形在轨道动力学里解决什么问题

1.1 低能量转移背后的“免费通道”

先说一个容易被忽略的事实：二体问题里，航天器的轨道是圆锥曲线，开普勒定律说了算，能量严格守恒，轨道类型不会无缘无故改变。但在三体问题里，事情不一样了。第三个物体受到的引力来自两个主天体，加上旋转坐标系带来的离心力和科里奥利力，系统不再是一个干净的“中心力场”。这时候，在拉格朗日点附近会出现一类特殊的动态结构：有些轨道方向会指数级地“吸引”物体靠近，有些方向会指数级地“推离”物体。这些方向分别对应稳定流形和不稳定流形。

低能量转移的本质，就是把这两种流形接起来。探测器先被稳定流形“吸”向某个平动点，然后又顺着不稳定流形“甩”出去，整个过程里主引擎几乎可以不点火。听上去很玄，但计算上完全可复现：只要你解出三体问题的运动方程，数值求一下特征向量，把一小段扰动沿特征方向放进去，再积分一段时间，流形曲线就出来了。

1.2 流形不是抽象的数学概念，而是可画出来的轨道曲面

很多人第一次接触“流形”这个词，脑子里浮现的是微分几何课本的抽象定义，总觉得这东西离工程十万八千里。实际上是反的。在CR3BP（圆型限制性三体问题）里，稳定流形和不稳定流形就是一些具体的、可数值生成的轨道集合。它们可以画成相空间里的曲线或曲面管，可以被初等线性代数和数值积分算出来，不需要抽象代数。

打个比方。如果一个平动点是旋转坐标系里的一个“漩涡中心”，稳定流形就是通向漩涡的滑梯，物体解开开始靠近它，最终慢慢盘旋进入；不稳定流形则是从漩涡甩出去的滑道，物体一旦进入这个方向就会快速离开。在拉格朗日点L1、L2这类共线平动点附近，滑梯和滑道都存在。它们的交叉、连接，构成了太阳系里四通八达的管道网络，业内常称“行星际运输网络”。

1.3 谁需要算流形

需要算流形的人，按接触深度大致分三类：

• 任务初步设计：要评估低能量转移方案、设计绕飞轨道、判断燃料预算，流形是必算项。
• 轨道动力学课程和科研项目：研究周期轨道族、计算转移轨道、做延拓分析，稳定/不稳定流形几乎是入门课。
• 自学三体问题的爱好者：如果只看平动点稳定性而忽略流形，相当于只学了地图上的几个地标，没学路网。

对这三类人来说，下面这套流程都通用：先搭CR3BP模型，再找平动点，然后做线性化分析，最后数值延拓出完整流形。接下来的内容，就是把这套流程一步步拆开。

2. CR3BP基础：方程、坐标和质量参数

2.1 无量纲化：让方程变得干净

计算流形之前，先把三体问题的模型固定下来。完整的三体问题，即使牛顿自己来了也头疼，所以工程分析几乎都用它的简化版本：圆型限制性三体问题。

这个模型的前提有三个：两个主天体（大天体和小天体）质量远大于第三个物体，所以第三个物体不影响前两者的运动；两个主天体做圆轨道运动；第三个物体只在两个主天体的引力场里运动。这比完全三体问题简单得多，但已经足够反映出拉格朗日点流形的大部分基本性质。

为了让数学形式干净，我们会做无量纲化。三个基本单位是：

• 长度单位取两个主天体之间的平均距离；
• 质量单位取两个主天体质量之和；
• 时间单位取两个主天体相对运动的角速度为1所需要的时长。

做这个无量纲化之后，方程里不再出现引力常数G，不再出现真实轨道半径，只剩下一个关键参数：

μ = m2 / (m1 + m2)

其中m2是较小主天体的质量。μ的值决定了引力场的“不对称程度”，是CR3BP里最重要的一个数。地月系统差不多是0.012左右，太阳系里很多小行星-主星系统也在这个量级附近。

2.2 运动方程与雅可比常数

在质心旋转坐标系里，令两个主天体分别位于(-μ, 0, 0)和(1-μ, 0, 0)，第三个物体的位置是(x, y, z)，速度是(x', y', z')。运动方程可以写成：

x'' - 2y' = U_x

y'' + 2x' = U_y

z'' = U_z

这里U是一个“有效势”，包含了引力势和旋转坐标系引入的离心势：

U = (1/2)(x² + y²) + (1-μ)/r1 + μ/r2

其中 r1 = sqrt((x+μ)² + y² + z²)，r2 = sqrt((x-1+μ)² + y² + z²)。U的下标表示对相应坐标求偏导。

注意方程中那两个速度耦合项-2y'和+2x'，它们来自旋转坐标系里的科里奥利效应，正是这个耦合，让三体系统比二体系统多出一大堆复杂动力学。

这个系统存在一个沿轨道守恒的标量，叫雅可比常数：

C = x² + y² + 2(1-μ)/r1 + 2μ/r2 - (x'² + y'² + z'²)

雅可比常数不是一个能量，但它扮演了类似能量的角色，用来约束航天器可达的区域。计算过程中，它也是一个非常实用的“体检指标”：数值积分做得好不好，看一个轨迹上的C漂移多少就心里有数了。

2.3 共线平动点的数值求解

平动点是相对旋转坐标系静止的点，满足速度为零且加速度为零。令运动方程右端为零，解出的点就是拉格朗日点，共五个。其中三个在x轴上，叫共线平动点L1、L2、L3；另外两个与两个主天体构成等边三角形，叫三角平动点L4、L5。

共线平动点的y、z坐标都是零，问题退化成一维求根。把y=z=0代入U_x=0，得到一个关于x的非线性方程：

x - (1-μ)(x+μ)/|x+μ|³ - μ(x-1+μ)/|x-1+μ|³ = 0

这个方程直接用牛顿迭代求根，初值给得好，几步就收敛。实用的初值经验是：

• 求L1，取x = 1-μ - (μ/3)^(1/3)附近；
• 求L2，取x = 1-μ + (μ/3)^(1/3)附近；
• 求L3，取x = -1 + 5μ/8附近。

这些表达式只是迭代初值，不是精确解，别直接拿来当最终坐标。真正算流形的时候，平动点坐标必须精确到足够位，否则后续特征向量方向都会偏。

3. 稳定/不稳定流形的数学定义与几何结构

3.1 从线性稳定方向说起

理解流形之前，先回到线性系统。假设有一个平衡点x*，把运动方程在它附近做泰勒展开，取一阶项，得到线性系统：

δx' = A δx

其中A是状态雅可比矩阵。A的特征值决定这个平衡点的“性格”。如果某个特征值有负实部，对应方向上的扰动会随时间衰减，这是稳定方向；如果某个特征值有正实部，对应方向上的扰动会被放大，这是不稳定方向。特征值为纯虚数或者模为1，对应中性方向，扰动既不衰减也不放大，围绕平衡点做小幅振荡。

对于共线平动点L1、L2，A的特征值结构很有规律：一个正的实特征值，一个负的实特征值，还有两对纯虚特征值。所以它们属于“鞍点×中心×中心”的组合。正实特征值对应的方向是不稳定方向，负实特征值对应的方向是稳定方向。在平衡点附近，沿稳定方向给一个微小扰动，物体会慢慢滑回平衡点；沿不稳定方向给一个微小扰动，物体会立刻加速离开。

3.2 稳定流形定理：局部理论到全局延拓

线性分析只告诉我们平衡点附近一小块区域的行为。非线性系统的完整答案，由稳定流形定理给出。它的意思大体是：在双曲平衡点附近，所有最终趋于平衡点的状态，构成一条局部的稳定流形W^s；所有最终远离平衡点的状态，构成一条局部的不稳定流形W^u。这些流形是光滑的曲线或曲面，并且在平衡点处与对应的稳定/不稳定特征空间相切。

关键点在于“局部”这个词。稳定流形定理保证的只是平衡点附近一个小邻域里存在这些流形。但实际数值计算中，我们会把这个局部方向向外延拓，沿着数值积分一直走，得到大范围的稳定/不稳定流形。这个延拓过程是否严格成立，有深入的理论讨论，但工程上早就默认它是可靠的：沿着特征方向给一个很小的初始偏移，然后用高精度积分器积分，得到的曲线就是全局流形的数值近似。

3.3 平衡点流形与周期轨道流形的区别

又一个容易混淆的地方：不同对象对应的流形，维度不同。

• 对于平衡点L1/L2本身，稳定流形是一条一维曲线，不稳定流形也是一条一维曲线。每条曲线在平衡点两侧各有一个分支，所以一共两条稳定分支、两条不稳定分支。
• 对于周期轨道（比如晕轨道、李萨如轨道），稳定流形和不稳定流形是二维的管状曲面。因为轨道上的每个点都有它自己的稳定方向，所有点的稳定方向拼在一起，就围成了一个管道。

所以，如果任务目标是“在L1点附近的晕轨道上运行”，计算流形时要处理的是二维流形管。如果只是“从L1点附近飞过”，用平衡点的流形曲线就够用了。两种需求对应两种算法，下文都会介绍。

还有一个细节值得提前提一下：共线平动点虽然不稳定，但这种不稳定是慢速法向失稳。换句话说，正实特征值通常不大，物体沿不稳定方向离开，不会像爆炸一样瞬间飞走。正是这种“温和的不稳定”，让流形管道可以被航天器利用，也给了轨道设计数天到数周的机动窗口。

4. 流形计算的完整流程

4.1 第一步：变分方程与状态转移矩阵

计算流形，核心数学工具是状态转移矩阵。它是“附近状态偏差如何传递”的线性映射。设参考轨迹为X(t)，相邻轨迹为X(t)+δ(t)。把δ(t)的运动方程线性化，得到变分方程：

δ' = A(t) δ

A(t)是CR3BP右侧状态对状态的6×6雅可比矩阵，结构是：

A = [0, I; K, C]

这里K是有效势U的Hessian矩阵，C是科里奥利反对称矩阵：

C = [[0, 2, 0], [-2, 0, 0], [0, 0, 0]]

对位置速度都做小量展开，变分方程可以随主方程一起数值积分。具体做法是：把状态从6维扩到42维，后36维存一个6×6状态转移矩阵Φ，初值为单位阵。积分一个周期T后，Φ(t0, t0+T)就是单值矩阵，也叫Monodromy矩阵。

K矩阵的具体元素不复杂，都是三个空间坐标的二阶偏导。例如：

U_xx = 1 - (1-μ)/r1³ + 3(1-μ)(x+μ)²/r1⁵ - μ/r2³ + 3μ(x-1+μ)²/r2⁵

其他交叉项类似，都是标准的势函数二阶导数。实际写代码时，可以直接用解析表达式，不建议用数值差分代替。因为流形方向对特征向量的精度很敏感，数值差分引入的截断误差会污染特征分解结果。

4.2 第二步：平衡点或周期轨道的单值矩阵

如果研究对象是平衡点本身，不需要积分变分方程。直接把平衡点坐标带进雅可比矩阵A，做特征分解，就能得到稳定和不稳定方向。

如果研究对象是周期轨道，那就要先解出轨道。常见套路是微分修正打靶法：从平动点附近的小振幅周期轨道出发，用预测-修正迭代逐步增大振幅，延拓出晕轨道族。这个过程本身是个小工程，但假设你已经有一条周期轨道X*(t)，周期为T，接下来就简单了：在轨道上取一点X*(t0)，从它出发积分一个周期的变分方程，得到单值矩阵M。

单值矩阵的特征值分布是判断轨道类型的核心依据。对于一条Halo轨道，M的特征值长这样：

• 1和1：对应沿轨道切向的平凡方向；
• λ和1/λ，λ>1：对应不稳定和稳定方向；
• 两个模为1的共轭复特征值对：对应中心方向。

如果λ比1大得多，说明轨道的不稳定程度高，流形管展开很快，航天器进入管道后离开也更迅速。如果λ很接近1，说明轨道只是弱不稳定，流形方向不明显，数值计算时要格外小心。

4.3 第三步：提取局部稳定/不稳定方向

特征分解之后，要对特征值和特征向量做筛选。具体做法是：

• 在所有特征值里找到那个大于1的实特征值λ，它对应的特征向量就是不稳定方向v_u；
• 找到它的倒数1/λ，对应特征向量就是稳定方向v_s。

注意，实特征值在数值计算里可能附带有很小的虚部，这是数值误差，不必惊慌。遇到这种情况，直接取实部归一化即可。真正需要警惕的是复数中心方向的共轭特征对，它们模为1，不产生双曲流形，不要拿来当流形方向。

把特征向量按某种范数归一化后，初始扰动可以写成：

X_u0 = X* + ε · v_u

X_s0 = X* + ε · v_s

其中ε是一个很小的扰动幅度，经验区间见下一节。至于取+号还是-号，流形有两个分支，分开算即可。如果你想看两条分支，从周期轨道上每个采样点出发做±两组扰动。

4.4 第四步：数值延拓与流形管采样

初始扰动准备好了，倒数第二步也是决定性的一步：选对时间积分方向。

这条规则必须刻在脑子里：

• 不稳定流形：沿不稳定方向给扰动后，正向积分时间t，让扰动被指数放大，看到的就是从轨道向外逃逸的曲线。
• 稳定流形：沿稳定方向给扰动后，反向积分时间t，反向时间会让原本衰减的扰动变成放大，看到的就是从远处渐渐汇入轨道的曲线。

很多新手卡在稳定流形上，反复问“我沿稳定方向积分了，怎么看到的只是贴着轨道的一小段？”答案就是时间方向搞反了。正向积分稳定方向，扰动会被压回轨道附近，什么有趣的几何都看不到。必须反向积分，才能让全局流形“展开”。

积分时长一般取若干个轨道周期。停止条件不唯一，常见两种：一是固定无量纲时间，比如积分5～10个周期；二是按事件截断，比如到达离某个主天体指定距离的球面就停。具体用哪种，取决于你关心的任务事件发生在哪里。

如果对象是周期轨道，需要做流形管采样。做法是：在轨道上取N个采样点，对每个采样点都重复一次特征分解和初始扰动生成，然后积分得到N条曲线。把这些曲线拼在一起，就围成流形管的表面。采样密度不是越高越好，可视化和初步分析时取20～40个点就非常清楚了，太多反而糊成一团。

4.5 一个可以照抄的最小实现流程

下面给一段最小可运行的思路骨架，不提供完整函数实现，但流程是真实可用的。这里以L1平衡点的稳定/不稳定流形为例。

import numpy as np from scipy.integrate import solve_ivp mu = 0.012 # 数值求解共线平动点L1坐标，自备函数 xL = solve_lagrange(mu, 'L1') XL = np.array([xL, 0, 0, 0, 0, 0]) # CR3BP的6x6雅可比矩阵，自备函数 A = jacobian_cr3bp(XL, mu) # 特征分解 w, V = np.linalg.eig(A) idx_unstable = np.argmax(w.real) # 正实部对应不稳定方向 idx_stable = np.argmin(w.real) # 负实部对应稳定方向 vu = V[:, idx_unstable].real vs = V[:, idx_stable].real vu /= np.linalg.norm(vu) vs /= np.linalg.norm(vs) eps = 1e-6 # 不稳定流形：沿不稳定方向扰动后正向积分 X_u0 = XL + eps * vu sol_u = solve_ivp(cr3bp_rhs, [0, 10], X_u0, rtol=1e-13, atol=1e-13, method='DOP853') # 稳定流形：沿稳定方向扰动后反向积分 X_s0 = XL + eps * vs sol_s = solve_ivp(cr3bp_rhs, [0, -10], X_s0, rtol=1e-13, atol=1e-13, method='DOP853')

cr3bp_rhs是CR3BP右侧函数，输入状态，输出六个导数。这个骨架换到周期轨道也基本成立，差异只在第二步：特征分解的对象从平衡点处的A矩阵，换成周期轨道上某一点的单值矩阵M。

5. 实操细节：参数选择、调试与避坑

5.1 初始扰动量ε到底该取多大

ε的选择，是流形计算里最典型的“细节决定成败”。

ε太大，初始状态离真实流形太远，后面的曲线根本不是流形上的轨迹，而是一条从流形附近漂出去的杂散轨道。ε太小，初始扰动被积分器舍入误差淹没，数值上“看不到”不稳定方向，尤其对弱不稳定轨道，这个问题更严重。我自己的经验区间是：

• 平衡点流形：ε取1×10⁻⁶到1×10⁻⁵比较稳；
• 周期轨道流形：ε取1×10⁻⁷到1×10⁻⁶，原因是你已经选在轨道附近，对误差更敏感。

怎么确定ε合不合适？做一个收敛性检查：分别用ε、ε/10、ε×10三档跑同一条流形。如果三条曲线在中段以后基本重合，说明这个ε在合理区间。如果三档曲线差异很大，要么ε太大，初始点不在流形上；要么积分方向或特征向量算错了。

5.2 积分器、容差与雅可比常数校验

CR3BP本身是光滑的常微分方程组，没有刚性问题的典型特征，所以普通高精度积分器就够。不过流形计算往往要长时间积分，尤其是反向积分稳定流形，轨迹可能要绕很多圈，对精度要求高。我建议用DOP853这类7/8阶显式Runge-Kutta，相对容差和绝对容差都设到1×10⁻¹³以下。

一个有效的误差监控手段是雅可比常数。理论上游离的CR3BP轨迹雅可比常数严格守恒。数值积分过程中，计算每个输出点的C值，如果相对变化在1×10⁻¹²量级，说明积分质量很好；如果漂到了1×10⁻⁶，说明轨迹已经明显偏离真实解，流形位置大概率也偏了。

不要只看积分器自带的误差估计就觉得万事大吉。误差估计是局部的，长时间积分后全局误差会累积，雅可比常数是少数便宜的全局诊断量，每跑一次流形都应该顺手算一下。

5.3 复数特征向量的处理

共线平动点的雅可比矩阵，特征值里混着一对纯虚数和一对纯虚共轭，对应中心方向。做特征分解时，这些特征值对应的特征向量也是复的。如果取错向量，把中心方向的复数特征向量直接塞进积分器，数值上会立刻出问题。

处理原则很简单：找双曲方向时，只认实数部分为正或负的那个实特征值。数值程序里可以用w.real.argmax()和w.real.argmin()去选。如果因为数值误差导致“实”特征值带了0.001量级的虚部，就把对应特征向量取实部，再归一化使用。

更稳妥的做法是，尽量避免在轨道上直接做特征分解，而是回到单值矩阵M。M的特征值结构比瞬时雅可比矩阵更清晰：1、1、λ、1/λ、两个单位模共轭对。从中挑出λ>1和1/λ<1的两个方向，对应流形方向，符号信息干净得多。

5.4 常见问题速查表

这张表里前两个问题，基本覆盖了绝大多数入门阶段的失败案例。我自己带过的项目里，新人跑流形卡住，十有八九是时间方向或者±分支的问题，而不是数学推导错。

6. 进阶：从流形到转移轨道设计

6.1 稳定流形是“免费刹车”

如果你想让航天器从远处自然减速，最终进入平动点附近的周期轨道，稳定流形就是那条“免费刹车”路径。航天器一旦进入稳定流形管道，即使在关发动机状态下，也会被引力场逐渐拉向目标轨道附近。任务轨道设计做的事，就是调整入轨速度和位置，让探测器落在稳定管道的“入口”。

当然，现实任务极少能只靠稳定流形“零误差”入轨。但利用稳定流形可以显著减小末端速度修正量，有时能省下几十到上百m/s的速度增量，这在实际工程里可能就是节省数百公斤推进剂。

6.2 不稳定流形是“免费推力”

从周期轨道离开时，不稳定流形又成了“免费推力”。探测器在晕轨道上，沿不稳定方向给一个很小的扰动，它就会自然沿着不稳定管道离开，速度增量可能只有几m/s。很多深空任务的轨迹设计，离开时就是刻意踩上不稳定流形，让引力场把人“推”出去，而不是靠发动机强行变轨。

低能量转移的秘诀，就是把一条稳定流形和一条不稳定流形在合适的位置“接起来”。接点上通常需要一个小的速度脉冲来修正两段管道的指向差。这也就是传统霍曼转移轨道的燃料消耗很大，而低能量转移却能让燃料花销大幅下降的原因。

6.3 一段小经验：怎么判断计算结果可不可信

最后分享一个我自己长期使用的判断技巧。算完一条轨道或一小段流形后，别急着画图、别急着对接，先做三件事的“体检”。

第一，看端点状态是否满足任务约束。流形曲线上的每一点，都有明确的坐标和速度，直接看它是否穿过你关心的目标区域。

第二，看关键物理量是否守恒。雅可比常数的漂移只允许在输出精度附近。如果曲线末端的C值比起点变了1×10⁻⁴，那这条曲线作为真实流形的可信度已经很低，重新积分，不要硬用。

第三，看相邻采样点的流形是否交叉杂乱。真正的流形管是光滑的，采样点之间的流形段不会剧烈交叉。如果画出来的管子像一团乱麻，大概率是特征向量方向算错，或者初始扰动偏大。

这三步都不需要额外工具，几分钟就能跑一遍，但能避免后续任务对接时浪费大量时间排查数值问题。我自己刚学流形计算那阵，就是少了这个习惯，结果花了很长时间在一条早已数值失真的“伪流形”上做转移方案，教训相当深刻。

算流形这件事，说难并不难，本质就是“求特征方向、给初始扰动、选对时间方向、积分”。但说容易也不容易，每一个细节，从平动点的精度、特征向量的归一化，到积分器的容差、ε的选取，都可能让最终结果差之千里。把基础概念和这套流程吃透，再去分析复杂的管道连接，你会觉得那些蜿蜒的流形曲线，其实到处都是明明白白的物理逻辑。