圆型限制性三体问题（CR3BP）：状态转移矩阵、微分修正与打靶法

#0830 笔记

本文承接《圆型限制性三体问题》和《圆型限制性三体问题（二）》前 19 节。此前已经建立 CR3BP 的旋转坐标系、无量纲运动方程、伪势函数、五个拉格朗日点、雅可比常数与零速度面。

前两篇主要回答：

给定一个初始状态，航天器会怎样运动？哪些区域在给定雅可比常数下允许访问？

本篇进一步回答一个反问题：

如果希望航天器最终满足指定条件，应该怎样修改初始位置、初始速度或飞行时间？

二点边值问题 $L_4,L_5$ 附近的周期轨道族。

沿用前文的六维无量纲状态向量：

x = {[\begin{matrix} x & y & z & \dot{x} & \dot{y} & \dot{z} \end{matrix}]}^{T} .

CR3BP 的三条二阶运动方程已经改写为六条一阶方程：

\dot{x} = f (x, τ),

$\tau$ $\boldsymbol\phi(\boldsymbol x,\tau)$ $\boldsymbol f$ $\boldsymbol\Phi$ 混淆。

20. 从初值问题到二点边值问题

20.1 初值问题：给定起点，向未来传播

普通数值积分解决的是初值问题。给定

x (τ_{0}) = x_{0},

就可以用 Runge–Kutta 等方法计算

x (τ) = ψ (x_{0}, τ - τ_{0}),

$\boldsymbol\psi$ 表示数值传播映射。它把“初始状态和传播时间”映射为“终止状态”。

这里的逻辑方向是

已知初始状态 ⟶ 数值积分 ⟶ 得到终止状态

例如，已知航天器的初始位置和速度，就能计算 5 天后的状态。

20.2 二点边值问题：给定两端要求，反求未知量

在任务设计中，经常只知道：

起点必须位于某条地球停泊轨道；
终点必须到达月球附近的指定位置；
终点速度方向必须满足交会要求；
飞行时间可能固定，也可能允许调整。

这时，完整初始状态并没有被直接给出。需要反过来寻找某些初始位置、初始速度或飞行时间，使传播后的终止状态满足目标。

这种问题称为 two-point boundary value problem，缩写为 TPBVP，中文通常译作二点边值问题。

其逻辑方向是

给定两端的部分条件 ⟶ 寻找合适的初值与时间 ⟶ 形成满足条件的轨道

“二点”指轨迹的两个边界，不表示空间中只有两个点；“边值”指在初始端和终止端规定的条件。

20.3 为什么叫打靶法

打靶法可以类比为投飞镖：

先选一个初始方向和速度；
投出飞镖，观察落点；
计算落点与靶心之间的偏差；
根据“瞄准方向变化一点会让落点怎样变化”修正下一次投掷；
重复，直到命中靶心。

在轨道设计中：

飞镖问题	轨道设计问题
投掷方向和速度	初始状态、机动量或飞行时间
飞镖轨迹	数值积分得到的航天器轨迹
靶心	目标边界条件
脱靶距离	约束误差向量
调整瞄准	修改自由变量

打靶法不是没有方向的反复试错。它使用偏导数描述终点对初值的敏感度，再系统计算下一次修正量。

20.4 微分修正的核心思想

$\delta\boldsymbol x_0$ $\delta\boldsymbol x_f$ 。

在参考轨道附近，可以近似写成

δ x_{f} \approx \frac{\partial x_{f}}{\partial x_{0}} δ x_{0} .

中间的偏导数矩阵告诉我们：初值中的各分量分别会怎样影响终点。这就是后文的状态转移矩阵。

所谓微分修正，就是利用这种一阶敏感度关系，把终点误差反推成初值修正。

21. 变分方程：小扰动怎样沿轨道传播

21.1 参考轨道与邻近轨道

$\boldsymbol x_0(\tau)$ 是一条参考轨道上的状态。考虑附近另一条轨道：

x (τ) = x_{0} (τ) + δ x (τ) .

其中：

$\boldsymbol x_0(\tau)$ 是参考解；
$\delta\boldsymbol x(\tau)$ 是小状态偏差；
$\boldsymbol x(\tau)$ 是邻近解。

这里的“小”非常重要。后续一阶 Taylor 展开只有在两条轨道足够接近时才可靠。

例如

δ x (0) = {[\begin{matrix} 0 & 0 & 0 & 0 & 10^{- 6} & 0 \end{matrix}]}^{T}

$y$ 方向速度增加一个很小的无量纲量。

21.2 对动力学方程作一阶 Taylor 展开

参考轨道满足

{\dot{x}}_{0} = f (x_{0}, τ) .

邻近轨道满足

{\dot{x}}_{0} + δ \dot{x} = f (x_{0} + δ x, τ) .

$\boldsymbol x_0$ 附近作一阶 Taylor 展开：

f (x_{0} + δ x, τ) \approx f (x_{0}, τ) + {\frac{\partial f}{\partial x} |}_{x_{0}} δ x .

消去参考轨道本身满足的项，得到线性变分方程：

δ \dot{x} = A (τ) δ x

其中

A (τ) = {\frac{\partial f}{\partial x} |}_{x_{0} (τ)}

是动力学雅可比矩阵。

$\boldsymbol x_0(\tau)$ 上求值。

$\boldsymbol x_0$ $A$ 是常矩阵。若参考解是周期轨道，则

A (τ + T) = A (τ) .

21.3 CR3BP 中的动力学雅可比矩阵

$\Omega(x,y,z)$ 。六条一阶运动方程为

\dot{x} = v_{x}, \dot{y} = v_{y}, \dot{z} = v_{z},

{\dot{v}}_{x} = 2 v_{y} + Ω_{x},

{\dot{v}}_{y} = - 2 v_{x} + Ω_{y},

{\dot{v}}_{z} = Ω_{z} .

因此

A (τ) = {[\begin{matrix} 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ Ω_{x x} & Ω_{x y} & Ω_{x z} & 0 & 2 & 0 \\ Ω_{y x} & Ω_{y y} & Ω_{y z} & - 2 & 0 & 0 \\ Ω_{z x} & Ω_{z y} & Ω_{z z} & 0 & 0 & 0 \end{matrix}]}_{x_{0} (τ)}

其中

Ω_{i j} = \frac{\partial^{2} Ω}{\partial i \partial j}, i, j \in {x, y, z} .

矩阵的结构可以分成四块：

A = [\begin{matrix} 0_{3 \times 3} & I_{3 \times 3} \\ \nabla^{2} Ω & S \end{matrix}],

其中

S = [\begin{matrix} 0 & 2 & 0 \\ - 2 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}] .

各分块的含义是：

分块	含义
$0_{3\times3}$	位置导数不直接依赖位置本身
$I_{3\times3}$	位置的时间导数就是速度
$\nabla^2\Omega$	伪势的局部曲率，描述位置扰动怎样改变加速度
$S$	旋转坐标系中的科里奥利耦合

$\Upsilon$ $\Omega$ $\Omega$ $S$ 。

21.4 变分方程不是新的引力模型

变分方程

δ \dot{x} = A δ x

不是在原运动方程之外增加一种力。它只是描述：

同一动力学系统中的两条邻近轨道，其差值在一阶近似下如何演化。

原轨道仍由非线性的 CR3BP 方程积分；变分方程只在参考轨道旁边同步传播敏感度。

22. 状态转移矩阵

22.1 定义

线性变分方程的解可以写成

δ x (τ) = Φ (τ, τ_{0}) δ x (τ_{0})

$\boldsymbol\Phi(\tau,\tau_0)$ 称为状态转移矩阵，英文是 state transition matrix，缩写为 STM。

$\tau_0$ $\tau$ 的线性映射。

22.2 STM 的每个元素是什么

$6\times6$ 矩阵：

Φ (τ, τ_{0}) = \frac{\partial x (τ)}{\partial x (τ_{0})} .

例如

Φ_{12} = \frac{\partial x (τ)}{\partial y (τ_{0})}

$y$ $x$ 位置怎样变化。

又如

Φ_{45} = \frac{\partial \dot{x} (τ)}{\partial \dot{y} (τ_{0})}

$y$ $x$ 方向速度的影响。

所以 STM 可以理解成一张完整的“初始误差传播表”。

22.3 STM 的微分方程

从

δ x (τ) = Φ (τ, τ_{0}) δ x (τ_{0})

$\tau$ 求导：

δ \dot{x} (τ) = \dot{Φ} (τ, τ_{0}) δ x (τ_{0}) .

另一方面，变分方程给出

δ \dot{x} = A (τ) δ x = A (τ) Φ (τ, τ_{0}) δ x (τ_{0}) .

因为这对任意初始扰动都成立，所以

\dot{Φ} (τ, τ_{0}) = A (τ) Φ (τ, τ_{0})

初始条件为

Φ (τ_{0}, τ_{0}) = I_{6 \times 6}

初始时刻采用单位矩阵的原因是：

每个状态分量对自身的偏导是 1；
每个状态分量对其他分量的偏导是 0。

22.4 为什么一次积分有 42 个变量

原始状态向量有 6 个分量。STM 有

6 \times 6 = 36

个元素。

因此计算参考轨道和 STM 时，要同时积分

6 + 36 = 42

个一阶微分方程。

数值积分器实际接收的扩展状态可以写成

y = [\begin{matrix} x \\ vec (Φ) \end{matrix}] \in R^{42},

$\operatorname{vec}$ 表示按列或按行把矩阵排成一个长向量。采用哪一种排列方式都可以，但程序中的展开和还原规则必须保持一致。

22.5 STM 的组合性质

$\tau_0<\tau_1<\tau_2$ ，则

Φ (τ_{2}, τ_{0}) = Φ (τ_{2}, τ_{1}) Φ (τ_{1}, τ_{0})

$\tau_0$ $\tau_1$ $\tau_1$ $\tau_2$ 。

若 STM 可逆，则

Φ (τ_{0}, τ) = Φ (τ, τ_{0})^{- 1} .

22.6 STM 只是一阶局部近似

STM 给出

x_{new} (τ) \approx x_{ref} (τ) + Φ (τ, τ_{0}) δ x_{0} .

$\delta\boldsymbol x_0$ 太大，高阶 Taylor 项不能忽略，此近似就会变差。

因此微分修正必须迭代进行：

在当前参考轨道上计算 STM；
得到一个小修正；
使用修正后的真实初值重新积分非线性运动方程；
在新轨道上重新计算 STM；
重复直到约束满足。

23. 同时变分与非同时变分

23.1 同时变分

同一时刻 $\tau$ 的状态差：

δ x (τ) = x (τ) - x_{0} (τ) .

如果传播时间固定，这就是需要研究的状态差。

23.2 非同时变分

$\delta\tau$ ，则比较的是

x (τ + δ τ)

$\tau$ 的状态。

一阶展开给出

x (τ + δ τ) \approx x (τ) + \dot{x} (τ) δ τ .

因此非同时变分满足

δ q (τ) = δ x (τ) + \dot{x} (τ) δ τ

其中：

$\delta\boldsymbol x$ 是初始状态改变引起的同时状态差；
$\dot{\boldsymbol x}\delta\tau$ 是多传播或少传播一小段时间引起的变化；
$\delta\boldsymbol q$ 是同时考虑两种影响后的总差。

图 3.1：同时变分和非同时变分之间的关系

可以把它类比成沿公路行驶的两辆车：

比较同一时刻两车的位置，是同时变分；
$\delta t$ $\boldsymbol v\delta t$ 的位移，这是非同时变分。

23.3 时间偏导为什么等于状态导数

传播映射写成

x_{f} = ψ (x_{0}, t) .

改变终止时间而保持初始状态不变时：

\frac{\partial x_{f}}{\partial t} = {\dot{x}}_{f}

因此可变时间打靶法的敏感度矩阵，会在 STM 之外增加一列终点状态导数。

23.4 固定时间与可变时间

类型	时间是否作为自由变量	使用的变分
固定时间	$\delta t=0$	同时变分
可变时间	$\delta t$ 可调整	非同时变分

可变时间增加了设计自由度，通常能满足更多约束，也可能扩大收敛域；代价是问题维数和实现复杂度增加。

24. 自由变量、约束向量与雅可比矩阵

24.1 自由变量向量

把算法允许调整的量收集成自由变量向量：

X = {[\begin{matrix} X_{1} & X_{2} & \dots & X_{n} \end{matrix}]}^{T}

自由变量可以包括：

初始位置的某些分量；
初始速度的某些分量；
飞行时间；
多重打靶中的拼接点状态；
$\Delta\boldsymbol v$ ；
用于处理不等式的松弛变量。

“自由”只表示求解器可以调整，并不表示物理任务中毫无限制。

24.2 约束或误差向量

把希望归零的误差收集成

F (X) = {[\begin{matrix} F_{1} (X) & F_{2} (X) & \dots & F_{m} (X) \end{matrix}]}^{T}

$\boldsymbol X^*$ 满足

F (X^{*}) = 0

$\boldsymbol r_D$ ，可以定义

F (X) = r_{f} (X) - r_{D} .

$\boldsymbol v_D$ ，则可以扩展为

F (X) = [\begin{matrix} r_{f} - r_{D} \\ v_{f} - v_{D} \end{matrix}] .

24.3 约束雅可比矩阵

约束对自由变量的偏导数组成

D F (X) = \frac{\partial F}{\partial X}

$m\times n$ 矩阵：

D F = [\begin{matrix} \frac{\partial F_{1}}{\partial X_{1}} & \dots & \frac{\partial F_{1}}{\partial X_{n}} \\ ⋮ & ⋱ & ⋮ \\ \frac{\partial F_{m}}{\partial X_{1}} & \dots & \frac{\partial F_{m}}{\partial X_{n}} \end{matrix}] .

矩阵元素回答：

$j$ $i$ 个约束误差改变多少？

24.4 一阶 Taylor 展开与 Newton 修正

$i$ $\boldsymbol X^i$ 附近：

F (X^{i + 1}) \approx F (X^{i}) + D F (X^{i}) (X^{i + 1} - X^{i}) .

定义

Δ X^{i} = X^{i + 1} - X^{i} .

要求下一次误差的一阶预测为零：

D F (X^{i}) Δ X^{i} = - F (X^{i})

$\Delta\boldsymbol X^i$ 后更新

X^{i + 1} = X^{i} + Δ X^{i}

这就是微分修正最核心的两步。

24.5 方阵 Newton–Raphson 更新

若自由变量数与约束数相等，即

n = m,

$D\boldsymbol F$ 非奇异，则形式上有

X^{i + 1} = X^{i} - D F (X^{i})^{- 1} F (X^{i})

但实际计算通常不显式形成矩阵逆，而是直接求解线性方程

D F Δ X = - F .

这样通常更快、更稳定。

24.6 自由变量多于约束：最小范数修正

若

n > m,

问题是欠定的：满足一阶约束的修正可能有无穷多个。论文选择离当前猜测最近的最小范数修正：

Δ X = - D F^{T} {(D F D F^{T})}^{- 1} F

它可以理解为：在所有能消除一阶误差的修正中，选择长度最小的一项。

现代数值实现通常使用 QR 分解或 SVD 求伪逆，避免直接计算上式中的逆矩阵，并能更好地处理病态问题。

24.7 约束多于自由变量

若

n < m,

问题超定。一般无法让所有约束同时精确为零，只能寻找最小二乘近似：

min_{Δ X} {‖ D F Δ X + F ‖}_{2} .

若任务要求所有约束必须严格满足，就需要重新设计自由变量或删除互相冲突的约束。

24.8 一个一维直观例子

$1000\ \mathrm{s}$ $100\ \mathrm{km}$ 。当前约束误差为

F = + 100 km .

若终点位置对初速度的敏感度约为

\frac{\partial F}{\partial v_{0}} \approx 1000 s,

则一阶修正为

Δ v_{0} = - \frac{F}{\partial F / \partial v_{0}} = - \frac{100 km}{1000 s} = - 0.1 km / s .

$100\ \mathrm{m/s}$ 。

真实 CR3BP 打靶只是把“一个误差除以一个敏感度”扩展为多维矩阵方程。

25. 单重打靶法

25.1 基本结构

单重打靶只使用一整段轨道：

x_{0} \overset{传播 t}{\to} x_{f} .

$\boldsymbol X^i$ $\boldsymbol F(\boldsymbol X^i)$ 。随后使用雅可比矩阵计算新猜测。

图 3.2：单重打靶从初始弧逐步修正到目标弧

图中：

$\boldsymbol X^0$ $\boldsymbol F(\boldsymbol X^0)$ ；
$\boldsymbol X^i$ 是中间迭代，误差已经减小；
$\boldsymbol X^*$ $\boldsymbol F(\boldsymbol X^*)=0$ 。

25.2 固定时间单重打靶

若传播时间固定，自由变量只包含选定的状态分量。例如

X = {[\begin{matrix} {\dot{x}}_{0} & {\dot{y}}_{0} \end{matrix}]}^{T} .

若约束是终点的两个位置分量：

F = [\begin{matrix} x_{f} - x_{D} \\ y_{f} - y_{D} \end{matrix}],

则雅可比矩阵中的偏导可以直接从 STM 对应行列提取：

D F = [\begin{matrix} \frac{\partial x_{f}}{\partial {\dot{x}}_{0}} & \frac{\partial x_{f}}{\partial {\dot{y}}_{0}} \\ \frac{\partial y_{f}}{\partial {\dot{x}}_{0}} & \frac{\partial y_{f}}{\partial {\dot{y}}_{0}} \end{matrix}] .

25.3 可变时间单重打靶

$t$ 也是自由变量：

X = [\begin{matrix} 初始状态中的可调分量 \\ t \end{matrix}] .

雅可比矩阵相应增加时间列：

D F = [\begin{matrix} 由 STM 提取的状态敏感度 & \frac{\partial F}{\partial t} \end{matrix}] .

若约束直接是终点状态，则时间列通常包含

{\dot{x}}_{f} .

25.4 单重打靶的优点和局限

优点：

变量较少；
结构简单；
初始猜测良好、轨迹较短时计算效率高。

局限：

长时间传播会放大初值误差；
靠近主天体或不稳定区域时，STM 元素可能非常大；
一次初值修正可能导致终点发生巨大变化；
收敛结果强烈依赖初始猜测。

单重打靶不是“不准确”，而是在高度敏感的长弧问题中数值条件可能很差。

26. 多重打靶法

26.1 把长轨道切成多个短段

多重打靶将一条长轨道切成若干短弧。每个短弧的起点称为 patch point，可译为拼接点或匹配点。

$k$ 个拼接点：

x_{1}, x_{2}, \dots, x_{k} .

$\boldsymbol x_i$ $t_i$ 得到

x_{i}^{+} = ψ (x_{i}, t_{i}) .

理想情况下，这个终点必须等于下一段的起点：

x_{i}^{+} = x_{i + 1} .

26.2 连续性约束

$i$ 个连续性误差：

c_{i} = ψ (x_{i}, t_{i}) - x_{i + 1} = 0

把所有连续性误差堆叠起来：

F_{k} (X) = [\begin{matrix} c_{1} \\ c_{2} \\ ⋮ \\ c_{k - 1} \end{matrix}] .

$\boldsymbol c_i$ $k$ 个拼接点之间共有

6 (k - 1)

个内部连续性约束。

图 3.3：多重打靶中的分段轨道与连续性误差

图中的虚线不是航天器真实执行的跳跃，而是当前初始猜测中各分段没有接上的误差。微分修正要同时移动拼接点，使所有虚线逐渐缩到零。

26.3 固定时间多重打靶的自由变量

若每段传播时间固定，自由变量通常包含全部拼接点以及其他设计量：

X = [\begin{matrix} x_{1} \\ x_{2} \\ ⋮ \\ x_{k} \\ X_{l} \end{matrix}] .

$\boldsymbol X_l$ $l$ ，则

n = 6 k + l .

约束向量由连续性约束和任务特有约束组成：

F (X) = [\begin{matrix} F_{k} (X) \\ F_{c} (X) \end{matrix}] .

26.4 连续性雅可比矩阵的带状结构

对

c_{i} = ψ (x_{i}, t_{i}) - x_{i + 1}

求偏导：

\frac{\partial c_{i}}{\partial x_{i}} = Φ_{i} (t_{i}, 0),

\frac{\partial c_{i}}{\partial x_{i + 1}} = - I_{6 \times 6} .

它对其他不相邻拼接点的偏导为零。因此连续性雅可比矩阵具有形式

\frac{\partial F_{k}}{\partial X_{k}} = [\begin{matrix} Φ_{1} & - I & 0 & \dots & 0 \\ 0 & Φ_{2} & - I & \dots & 0 \\ ⋮ & ⋱ & ⋱ & ⋮ \\ 0 & \dots & 0 & Φ_{k - 1} & - I \end{matrix}]

大量零块使其成为稀疏带状矩阵。多重打靶虽然变量更多，但可以使用稀疏线性代数高效求解。

26.5 可变时间多重打靶

如果每段时间也允许修改，则

T = {[\begin{matrix} t_{1} & t_{2} & \dots & t_{k - 1} \end{matrix}]}^{T}

加入自由变量向量：

X = [\begin{matrix} X_{k} \\ T \\ X_{l} \end{matrix}] .

连续性误差对本段时间的偏导为

\frac{\partial c_{i}}{\partial t_{i}} = {\dot{x}}_{i}^{+}

因此时间敏感度部分通常是一个块对角结构：

\frac{\partial F_{k}}{\partial T} = [\begin{matrix} {\dot{x}}_{1}^{+} & 0 & \dots & 0 \\ 0 & {\dot{x}}_{2}^{+} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & {\dot{x}}_{k - 1}^{+} \end{matrix}] .

26.6 只使用总时间作为一个自由变量

$T$ 作为自由变量，并让各段时间相等：

t_{i} = \frac{T}{k} .

根据链式法则：

\frac{\partial c_{i}}{\partial T} = \frac{1}{k} {\dot{x}}_{i}^{+} .

因而时间敏感度减少为一列：

\frac{\partial F_{k}}{\partial T} = \frac{1}{k} [\begin{matrix} {\dot{x}}_{1}^{+} \\ {\dot{x}}_{2}^{+} \\ ⋮ \\ {\dot{x}}_{k - 1}^{+} \end{matrix}] .

这种做法自由变量较少，但不能独立调整每一段的传播时间。

26.7 多重打靶为什么更稳健

设一条长轨道经过高度敏感区域。单重打靶中的初值误差需要传播整段时间，可能被放大许多数量级。

多重打靶把轨道切短后：

每个分段的误差只传播较短时间；
每个拼接点都可以被调整；
很大的全局误差被拆成多个较小的局部连续性误差；
数值条件通常得到改善。

代价是：

自由变量显著增加；
连续性约束显著增加；
需要更复杂的稀疏矩阵组织。

可以把它类比成长距离修路：单重打靶只调整整条道路的起点方向；多重打靶允许在沿途多个控制点同时调整方向，最后再把所有路段平滑接起来。

27. 打靶法的完整计算流程

论文将轨道设计过程归纳为九个步骤。结合前文，可以整理为：

27.1 选择动力学模型

根据任务精度选择：

二体模型；
CR3BP；
双圆限制性四体模型；
椭圆限制性三体模型；
$N$ 体模型。

模型越复杂，不一定越适合直接生成初始猜测。工程上常先用简单模型得到近似解，再逐步提高模型精度。

27.2 选择打靶形式

确定：

单重还是多重；
固定时间还是可变时间。

选择取决于轨迹长度、敏感性、可用初始猜测和任务约束。

27.3 定义自由变量

明确求解器允许修改哪些量，并记录其单位和排列顺序。

27.4 定义约束向量

把所有任务条件写成“希望等于零”的误差函数。

27.5 构造雅可比矩阵

可以使用：

解析偏导；
STM；
有限差分；
自动微分。

有限差分容易实现，但步长太大会产生截断误差，太小又会受浮点舍入误差影响。STM 通常更适合精确、重复的轨道修正。

27.6 产生初始轨道猜测

初始猜测可以来自：

线性化解；
二体 Lambert 解；
已知周期轨道；
延拓得到的相邻解；
更简单动力学模型中的轨迹；
人工选取的拼接点。

27.7 初始化自由变量向量

$\boldsymbol X^0$ 。

27.8 反复微分修正

每一轮迭代：

用当前自由变量积分轨迹及 STM；
$\boldsymbol F(\boldsymbol X^i)$ ；
$\|\boldsymbol F\|$ $\varepsilon$ ；
$D\boldsymbol F(\boldsymbol X^i)$ ；
$\Delta\boldsymbol X^i$ ；
$\boldsymbol X^{i+1}=\boldsymbol X^i+\Delta\boldsymbol X^i$ 。

常用停止条件为

‖ F (X^{i}) ‖_{2} < ε

还应设置最大迭代次数，并监测修正量是否异常增大。

27.9 重新生成并验证完整轨迹

$\boldsymbol X^*$ 通常只包含初始状态或拼接点。仍需再次进行高精度数值积分，生成完整连续轨迹，并验证：

约束是否确实满足；
拼接处是否连续；
雅可比常数是否按预期守恒；
是否撞击主天体；
时间和速度是否符合物理要求；
对小误差是否过度敏感。

27.10 为什么初始猜测十分关键

Newton 法只在当前猜测附近使用一阶近似，因此属于局部收敛方法。

初始猜测太差时可能出现：

一阶 Taylor 近似失效；
约束误差不降反升；
雅可比矩阵接近奇异；
修正量变得非常大；
数值积分进入碰撞奇点或完全不同的动力学区域。

因此“算法没有收敛”不一定说明目标轨道不存在，也可能说明自由变量、约束、缩放、初始猜测或打靶形式选择不合适。

28. 用可变时间单重打靶生成平面 Lyapunov 轨道

28.1 平面 Lyapunov 轨道

$xy$ $L_1,L_2,L_3$ 附近的周期轨道。

它不是航天器静止在拉格朗日点，而是在旋转坐标系中围绕该区域形成闭合运动。

周期条件是

x (T) = x (0) .

28.2 在共线点附近线性化

$L_i$ 的平衡状态为

x_{L_{i}} = {[\begin{matrix} x_{L_{i}} & 0 & 0 & 0 & 0 & 0 \end{matrix}]}^{T} .

定义相对状态

ξ = x - x_{L_{i}} .

在平衡点附近线性化：

\dot{ξ} = A_{L_{i}} ξ .

共线点附近通常存在：

一对实特征值，对应指数增长和衰减方向；
一对平面纯虚特征值，对应平面振荡；
一对垂直纯虚特征值，对应面外振荡。

为了构造平面周期初值，抑制实特征值对应的不稳定模式和面外模式，只保留平面纯虚特征值对应的正弦运动。

线性近似可写成

ξ (t) = ξ_{0} \cos ω t + \frac{η_{0}}{β_{3}} \sin ω t,

η (t) = η_{0} \cos ω t - ξ_{0} β_{3} \sin ω t,

其中

β_{3} = \frac{ω^{2} + Ω_{x x, L_{i}}}{2 ω} .

线性解只负责产生合理初始猜测，不是非线性 CR3BP 中的精确周期轨道。

$x$ 轴的镜像对称性

$y=0$ $y=0$ ，则另一半可以由对称性恢复。

因此不必直接施加完整六维首尾周期约束，只需要求解半圈。

图 3.4：利用 $y=0$ 对称面修正平面 Lyapunov 轨道

$\Sigma$ 表示

Σ : y = 0.

初始状态可选为

x_{0} = {[\begin{matrix} x_{0} & 0 & 0 & 0 & {\dot{y}}_{0} & 0 \end{matrix}]}^{T} .

$x_0$ $\dot y_0$ 由线性近似给出初始猜测。

28.4 自由变量和约束

论文使用可变时间单重打靶，自由变量为

X = [\begin{matrix} {\dot{y}}_{0} \\ t \end{matrix}]

$y=0$ ，并垂直穿越该平面：

F (X) = [\begin{matrix} {\dot{x}}_{f} \\ y_{f} \end{matrix}] = 0

这里：

$y_f=0$ 保证到达对称面；
$\dot x_f=0$ $x$ 轴上的对称面。

28.5 对应雅可比矩阵

两个约束、两个自由变量，所以

D F = [\begin{matrix} \frac{\partial {\dot{x}}_{f}}{\partial {\dot{y}}_{0}} & \frac{\partial {\dot{x}}_{f}}{\partial t} \\ \frac{\partial y_{f}}{\partial {\dot{y}}_{0}} & \frac{\partial y_{f}}{\partial t} \end{matrix}] .

从 STM 和终点状态导数提取：

D F = [\begin{matrix} Φ_{45} (t, 0) & {\ddot{x}}_{f} \\ Φ_{25} (t, 0) & {\dot{y}}_{f} \end{matrix}]

其中状态顺序采用

[x, y, z, \dot{x}, \dot{y}, \dot{z}]^{T} .

$\dot y_0$ $t$ ，直到两项约束都足够接近零。

28.6 从一条轨道生成整个轨道族

$x_0$ ，并把上一条轨道的解作为下一条的初始猜测。这称为参数延拓。

基本过程是：

已知轨道 ⟶ 参数稍微变化 ⟶ 预测新初值 ⟶ 微分修正 ⟶ 下一条轨道 .

当简单参数延拓在折返点附近失效时，可以使用伪弧长延拓。它沿轨道族在解空间中的切向方向前进，而不是强迫某一个物理参数始终单调变化。

图 3.5：地月系统中 $L_1,L_2,L_3$ 的平面 Lyapunov 轨道族

图中的每条闭合曲线都是不同初始条件对应的一条独立周期轨道，而不是同一航天器多圈运动留下的所有轨迹。

29. 通用周期轨道多重打靶法

29.1 为什么需要更通用的方法

$y=0$ 镜像对称性。它不适合：

三维周期轨道；
不关于坐标轴对称的轨道；
$L_4,L_5$ 附近的非对称轨道；
复杂星历模型中的周期或近周期解。

因此论文构造了可变时间多重打靶周期轨道算法。

29.2 周期闭合约束

$k$ 个拼接点

x_{0}, x_{1}, \dots, x_{k - 1} .

$T$ ，每段传播时间取

Δ t = \frac{T}{k} .

内部连续性约束为

ψ (x_{i}, T / k) - x_{i + 1} = 0, i = 0, \dots, k - 2.

周期轨道还必须把最后一段接回第一段：

ψ (x_{k - 1}, T / k) - x_{0} = 0

这条首尾连续性就是周期性约束。

29.3 为什么不能约束全部六个周期分量再任意增加周期

自治系统中的周期轨道具有相位自由度：沿同一条周期轨道移动起始点，仍然代表同一条几何轨道。

因此“所有拼接点状态加周期”通常会产生一个自由度冗余。需要额外的相位条件或固定一个几何条件，避免雅可比矩阵因时间平移对称性而退化。

$y$ $H$ $y_H$ ，并规定初始运动方向，消除这种不确定性。

29.4 用松弛变量表示运动方向

若希望初始运动为顺时针，例如

{\dot{y}}_{0} < 0,

$\beta$ ，把不等式写成

{\dot{y}}_{0} + β^{2} = 0

$\beta^2\ge0$ ，所以自动得到

{\dot{y}}_{0} = - β^{2} \leq 0.

$\beta^2$ 项的符号。

29.5 自由变量与约束维数

$k$ $T$ $\beta$ ：

X = [\begin{matrix} x_{0} \\ x_{1} \\ ⋮ \\ x_{k - 1} \\ T \\ β \end{matrix}]

因此

n = 6 k + 2.

论文的约束向量维数为

m = 6 k + 1,

所以

n > m .

该问题使用最小范数修正，而不是方阵 Newton 逆。

30. 三维 Halo 轨道

30.1 什么是 Halo 轨道

$z$ 方向振幅，在旋转坐标系中形成闭合的三维环状结构。

$L_i$ ”是对空间形态的描述。拉格朗日点不是一个具有实体表面的中心天体，航天器也不是像近地卫星一样由单一中心引力维持圆轨道。

30.2 单周期矩阵

$T$ 的周期轨道，经过一整周期的 STM 为

M = Φ (T, 0)

它称为 monodromy matrix，可译为单周期矩阵。

它将初始扰动映射到一整周期后的扰动：

δ x (T) = M δ x (0) .

$M$ 的某个特征值模长大于 1，对应方向的扰动每经过一圈都会放大；若模长小于 1，则会衰减；位于单位圆上的特征值对应中性振荡行为。

30.3 分岔与新轨道族

沿 Lyapunov 轨道族变化时，单周期矩阵的特征值也会变化。若特征值结构在某一点发生特定改变，可能出现分岔：一个新的周期轨道族从原轨道族中产生。

Halo 轨道族与平面 Lyapunov 轨道族在一个分岔轨道处相交。

生成 Halo 初值的基本过程是：

沿 Lyapunov 轨道族寻找分岔附近的轨道；
在该平面轨道上等时间选取若干拼接点；
$z$ 方向扰动；
用通用周期轨道多重打靶法恢复连续性和周期性；
得到第一条三维 Halo 轨道；
使用延拓生成整个 Halo 轨道族。

30.4 北部和南部 Halo 轨道族

$xy$ 平面对称。因此面外扰动的方向决定得到哪一支轨道族：

$z>0$ 一侧，称为北部 Halo 轨道；
$z<0$ 一侧，称为南部 Halo 轨道。

$xy$ 平面的镜像。

图 3.9：地月系统中共线点附近的北部 Halo 轨道族

理解三维轨道图时，应同时观察：

透视图；
$xy$ 投影；
$xz$ 投影；
$yz$ 投影。

$xy$ 投影，可能误以为轨道仍然是平面的。

30.5 Halo 轨道与雅可比常数

$z$ 振幅时，通常需要更大的速度幅值。由于

C = 2 Ω - v^{2},

速度增大通常对应雅可比常数下降。

$L_3$ $L_1,L_2$ 轨道族。这与前篇笔记中各拉格朗日点临界 Jacobi 常数的大小关系一致。

$L_4,L_5$ 附近的非对称短周期轨道

31.1 等边三角形点附近的线性运动

$L_4$ $L_5$ 附近定义偏差状态

ξ = x - x_{L_{i}}, i = 4, 5.

线性化后：

\dot{ξ} = A_{L_{i}} ξ .

$s_1,s_2$ ：

s_{1} = \frac{\sqrt{2}}{2} {(1 - \sqrt{1 - 27 μ (1 - μ)})}^{1 / 2},

s_{2} = \frac{\sqrt{2}}{2} {(1 + \sqrt{1 - 27 μ (1 - μ)})}^{1 / 2} .

由于

s_{1} < s_{2},

所以：

$s_1$ 对应长周期模式；
$s_2$ 对应短周期模式。

31.2 构造短周期初值

把长周期模式的系数设为零，只保留短周期模式，则线性近似可写成

ξ (t) = ξ_{0} \cos s_{2} t + \frac{{\dot{ξ}}_{0}}{s_{2}} \sin s_{2} t,

η (t) = η_{0} \cos s_{2} t + \frac{{\dot{η}}_{0}}{s_{2}} \sin s_{2} t .

线性周期初值为

T_{0} = \frac{2 π}{s_{2}} .

$k$ 个拼接点，再把它们交给通用周期轨道多重打靶法修正。

31.3 为什么使用通用多重打靶

$L_4,L_5$ $x$ 轴镜像结构。因此不能只传播半圈并施加两个对称性约束。

通用多重打靶直接施加：

所有分段连续；
最后一段接回第一段；
周期和运动方向满足要求。

所以它不依赖某一条特殊对称轴。

31.4 轨道族的演化

$L_4,L_5$ 的第一条小幅短周期轨道后，使用伪弧长延拓逐渐生成更大的轨道。

图 3.10：由 $L_4,L_5$ 附近线性短周期解延拓得到的周期轨道族

靠近平衡点的轨道近似为倾斜椭圆。随着轨道族向外延伸：

非线性效应越来越明显；
轨道形状不再接近简单椭圆；
$L_3$ 附近；
对应 Jacobi 常数发生系统变化。

这也说明：线性化解只在平衡点附近有效，但它可以作为非线性轨道族的种子。

32. 面向图形界面的通用多重打靶算法

32.1 为什么需要通用求解框架

如果每一种轨道任务都重新编写自由变量、约束和雅可比矩阵，设计效率会很低。论文因此把多种打靶形式统一到同一个图形界面框架中。

用户可以选择：

固定时间或可变时间；
只要求各段连续；
约束初始位置；
约束终止位置；
同时约束初始和终止位置。

32.2 四种位置约束选项

设第一和最后一个拼接点的位置分别为

r_{1} = {[\begin{matrix} x_{1} & y_{1} & z_{1} \end{matrix}]}^{T},

r_{k} = {[\begin{matrix} x_{k} & y_{k} & z_{k} \end{matrix}]}^{T} .

$\boldsymbol r_0$ $\boldsymbol r_f$ 。

四种约束选项是：

只约束内部连续性；
$\boldsymbol r_1-\boldsymbol r_0=0$ ；
$\boldsymbol r_k-\boldsymbol r_f=0$ ；
连续性加初始和终止位置约束。

再分别配合固定时间和可变时间，总共得到

4 \times 2 = 8

种打靶配置。

32.3 最完整的自由变量和约束

可变时间多重打靶的完整自由变量可写成

X = [\begin{matrix} x_{1} \\ x_{2} \\ ⋮ \\ x_{k} \\ t_{1} \\ t_{2} \\ ⋮ \\ t_{k} \end{matrix}] .

约束向量由三部分组成：

F (X) = [\begin{matrix} F_{k} (X) \\ F_{c 0} (X) \\ F_{c f} (X) \end{matrix}]

其中：

F_{c 0} = r_{1} - r_{0},

F_{c f} = r_{k} - r_{f} .

其他七种形式都可以从最完整形式中删除相应行列得到：

删除终点约束及对应雅可比矩阵行，就取消终点位置约束；
删除起点约束，就取消起点位置约束；
删除时间自由变量及对应雅可比矩阵列，就切换为固定时间；
只保留连续性部分，就得到纯拼接修正。

这说明打靶算法的通用性来自统一的

自由变量 + 约束 + 雅可比矩阵

表示方式。

33. 数值实现中的实际问题

33.1 不要显式计算矩阵逆

虽然公式常写成

Δ X = - D F^{- 1} F,

程序中更合适的做法是直接解

D F Δ X = - F .

对于方阵可以使用 LU 或 QR；对于欠定、超定或接近奇异的问题，可以使用 SVD。

33.2 变量和约束的尺度

自由变量中可能同时包含：

$1$ ；
很小的速度修正；
很长的传播时间；
不同物理意义的任务参数。

若数量级差异太大，雅可比矩阵会病态。可以引入缩放矩阵：

\tilde{X} = S_{X}^{- 1} X,

\tilde{F} = S_{F}^{- 1} F .

让各变量和约束在数值上具有相近数量级，通常能改善求解稳定性。

33.3 阻尼修正

若完整 Newton 步过大，可以采用

X^{i + 1} = X^{i} + α Δ X^{i}, 0 < α \leq 1.

$\alpha<1$ 时称为阻尼 Newton 步。它牺牲单步速度，换取更稳定的误差下降。

33.4 有限差分检查 STM 和解析雅可比

可以用中心差分检查某一列敏感度：

\frac{\partial x_{f}}{\partial X_{j}} \approx \frac{x_{f} (X + h e_{j}) - x_{f} (X - h e_{j})}{2 h} .

再与 STM 或解析雅可比对应列比较。若差异很大，可能存在：

状态排列错误；
STM 展开顺序错误；
二阶偏导公式错误；
时间尺度混淆；
有限差分步长不合适。

33.5 收敛误差小不等于任务一定合理

即使

‖ F ‖ < ε,

也只表示所定义的数学约束满足。仍要检查：

是否穿过主天体内部；
是否需要不可实现的速度变化；
飞行时间是否合理；
是否违反零速度面的可达性；
周期轨道是否极度不稳定；
高保真动力学模型下是否仍然有效。

34. 概念梳理

变分分析不是重新求一条轨道

它研究参考轨道附近的小扰动如何传播。真实的新轨道仍需重新积分非线性运动方程。

$A$ $D\boldsymbol F$ 不同

A = \frac{\partial f}{\partial x}

描述瞬时动力学对状态的敏感度；

D F = \frac{\partial F}{\partial X}

描述任务约束对设计自由变量的敏感度。STM 把二者联系起来。

STM 是局部线性映射，不是全局精确映射

δ x_{f} \approx Φ δ x_{0}

只在扰动足够小时可靠。

上游和下游指时间传播顺序，不是空间上下方向

上游状态是较早时刻的状态，下游状态是由数值积分得到的较晚时刻状态。

打靶不是盲目试错

$D\boldsymbol F$ 计算误差最敏感的修正方向。

单重打靶不是只迭代一次

“单重”指轨道只由一个传播段组成，不是只作一次 Newton 修正。

多重打靶不是多条不同任务轨道

它把同一条期望轨道分成多个短段，并用连续性约束把它们连接起来。

可变时间多出的雅可比列来自终点状态导数

\frac{\partial x_{f}}{\partial t} = {\dot{x}}_{f} .

周期轨道族不是同一航天器的多圈轨迹

图中的每条闭合曲线对应一组不同的初始条件和周期。

平衡点、周期轨道和稳定性是三个不同概念

平衡点：在旋转系中状态不变；
周期轨道：经过一个周期后状态重复；
稳定性：小扰动经过传播后是增长、衰减还是保持有界。

允许区域不保证打靶法一定收敛

零速度面只给出能量型必要条件。打靶法能否收敛还取决于动力学结构、初始猜测、约束选择和数值条件。

35. 本篇涉及的符号术语表

符号	含义
$\boldsymbol x$	$[x,y,z,\dot x,\dot y,\dot z]^T$
$\boldsymbol f(\boldsymbol x,\tau)$	六维一阶动力学函数
$\boldsymbol x_0(\tau)$	参考轨道上的状态
$\delta\boldsymbol x$	同一时刻参考轨道与邻近轨道的状态差
$A(\tau)$	$\partial\boldsymbol f/\partial\boldsymbol x$
$\nabla^2\Omega$	伪势函数的 Hessian 矩阵
$\boldsymbol\Phi(\tau,\tau_0)$	状态转移矩阵 STM
$\delta\boldsymbol q$	同时考虑状态变化和时间变化的非同时变分
$\delta\tau$	传播时间的小变化
$\boldsymbol\psi(\boldsymbol x_0,t)$	$t$ 的数值映射
TPBVP	two-point boundary value problem，二点边值问题
$\boldsymbol X$	自由变量或设计变量向量
$\boldsymbol F(\boldsymbol X)$	约束或误差向量
$D\boldsymbol F$	约束对自由变量的雅可比矩阵
$\Delta\boldsymbol X$	一次微分修正量
$\varepsilon$	约束收敛容差
patch point	多重打靶中的拼接点或匹配点
$\boldsymbol c_i$	$i$ 段与下一拼接点之间的连续性误差
$T$	周期轨道的完整周期或任务总传播时间
$M=\Phi(T,0)$	周期轨道的单周期矩阵 monodromy matrix
continuation	延拓法，由已知解逐步生成相邻解
pseudo-arclength continuation	伪弧长延拓
bifurcation	分岔，旧轨道族上产生新轨道族的结构变化
$\beta$	把运动方向不等式转为等式的松弛变量
Lyapunov orbit	共线点附近的平面周期轨道
Halo orbit	共线点附近的三维周期轨道

36. 本篇总结

给定初始状态后，CR3BP 运动方程可以产生一条轨道；但轨道设计需要反过来寻找能满足终点条件的初始状态和传播时间。

在参考轨道附近对动力学方程线性化，得到变分方程：

δ \dot{x} = A (τ) δ x

状态转移矩阵将初始扰动映射到终点扰动：

δ x (τ) = Φ (τ, τ_{0}) δ x (τ_{0})

它满足

\dot{Φ} = A Φ, Φ (τ_{0}, τ_{0}) = I

若传播时间也变化，则

δ q = δ x + \dot{x} δ τ

$\boldsymbol X$ $\boldsymbol F(\boldsymbol X)=0$ ，微分修正归结为线性方程：

D F Δ X = - F

单重打靶只传播一整段轨道，结构简单但对长时间敏感传播较脆弱；多重打靶把轨道切成短段，通过

ψ (x_{i}, t_{i}) - x_{i + 1} = 0

强制各段连续，变量更多但通常更稳健。

$y=0$ $L_4,L_5$ 附近的非对称轨道则使用更通用的周期性多重打靶和延拓法生成。

本篇的核心逻辑可以压缩为：

CR3BP 运动方程 ⟶ 变分方程 ⟶ 状态转移矩阵 ⟶ 约束雅可比矩阵 ⟶ 微分修正 ⟶ 单重／多重打靶 ⟶ 周期轨道与转移轨道

前两篇笔记给出了 CR3BP 中“轨道为什么这样运动”的动力学骨架；本篇进一步给出“怎样主动寻找一条满足任务条件的轨道”的数值方法。