圆形轨道的小推力碰撞规避

论文信息： Javier Hernando-Ayuso 和 Claudio Bombardelli, "Low-Thrust Collision Avoidance in Circular Orbits", Journal of Guidance, Control, and Dynamics, Vol. 44, No. 5, 2021, pp. 984–995. DOI: 10.2514/1.G005547

研究机构： ispace-inc（东京）/ 马德里理工大学 (Universidad Politécnica de Madrid)

1. 引言

本系列的第一篇解读了 Bombardelli & Hernando-Ayuso (2015) 的里程碑式工作——脉冲机动的最优碰撞规避：假设航天器在某一时刻瞬时施加 $\Delta v$，通过巧妙的数学变换将优化问题简化为 $3\times3$ 矩阵的本征值问题。但现实中，越来越多的航天器正采用电推进系统——尤其是即将部署的大型低轨星座（如 Starlink、OneWeb 等）——这些系统产生的是连续、小推力而非瞬时脉冲。

脉冲与小推力的本质区别：

本文将解读 Hernando-Ayuso & Bombardelli (2021) 发表于 JGCD 的论文，该论文首次系统地解决了圆形轨道下小推力碰撞规避的最优控制问题，并将前作 (2015) 的脉冲框架扩展到了连续推力领域。

2. 理论基础：b平面与碰撞概率

2.1 b平面坐标系

与前一篇文章一致，论文采用 b平面（b-plane）坐标系 $\langle \xi, \eta, \zeta \rangle$ 描述近距离交会。其中 $\mathbf{u}_\eta$ 沿相对速度方向，b平面（$\xi$-$\zeta$ 平面）垂直于相对速度，$\xi$ 轴指向两轨道的最小轨道交会距离（MOID）方向。

图1 展示了 b 平面在最近距离处的几何关系。 $\xi$ 轴指向纸面内，$\zeta$ 轴与 $\eta$ 轴构成右手系。机动航天器 S1 相对于目标 S2 的 b 平面位置向量为 $\mathbf{b} = (\xi, \zeta)^T$。图中还标注了两个物体的速度向量 $\mathbf{v}_1$、$\mathbf{v}_2$ 以及它们的相对速度 $\mathbf{v}_1 - \mathbf{v}_2$——这是整个碰撞规避理论的基础几何框架。

2.2 碰撞概率的代理目标函数

在短时交会假设下，论文采用 Chan 方法[Rician 积分]计算碰撞概率。关键洞察是：碰撞概率 $P$ 是 $v$ 的单调递减函数，其中：

$$v = \frac{\left(\frac{\xi}{\sigma_\xi}\right)^2 + \left(\frac{\zeta}{\sigma_\zeta}\right)^2 - 2\rho_{\xi\zeta}\frac{\xi}{\sigma_\xi}\frac{\zeta}{\sigma_\zeta}}{1-\rho_{\xi\zeta}^2}$$

因此，最小化碰撞概率等价于最大化以下代理目标函数：

$$J_P = \left(\frac{\xi}{\sigma_\xi}\right)^2 + \left(\frac{\zeta}{\sigma_\zeta}\right)^2 - 2\rho_{\xi\zeta}\frac{\xi\zeta}{\sigma_\xi\sigma_\zeta}$$

另一种可选的目标是最大化碰撞距离：

$$J_d = \xi^2 + \zeta^2$$

两个目标函数都可以写作 b 平面坐标的二次型 $J = \frac{1}{2}\mathbf{b}^T \mathbf{Q} \mathbf{b}$，区别仅在于 $\mathbf{Q}$ 矩阵的形式。

3. 圆形轨道的 b 平面运动学

3.1 与脉冲情形的简化差异

在 Bombardelli & Hernando-Ayuso (2015) 的脉冲模型中，机动航天器 S1 可以在任意偏心率的椭圆轨道上运行，需要处理完整的 $\mathbf{R}, \mathbf{K}, \mathbf{D}$ 三个变换矩阵。而在本文中，作者专注于两个物体均在圆形轨道上运行的情形——这对低轨大型星座具有重要实际意义——交会几何大大简化。

对于圆形轨道，两个物体的轨道速度大小近似相等（$\chi \simeq 1$），交会几何简化为轨道面夹角 $\kappa$ 的单一参数描述：

• 当 $0 < \kappa < \pi/2$：$\phi = 0, \psi = \mp\kappa$
• 当 $\pi/2 < \kappa < \pi$：$\phi = \pi, \psi = \pi \pm \kappa$

其中 $\phi$ 为轨道面内转角，$\psi$ 为面外转角。

3.2 简化的变换矩阵

b 平面位移与推力之间的关系简化为一个 $2\times3$ 的矩阵 $\mathbf{R}_0$：

$$\mathbf{R}_0 = \begin{bmatrix} 0 & 0 & -1 \\ -\cos\frac{\kappa}{2} & -\sin\frac{\kappa}{2} & 0 \end{bmatrix}$$

这比脉冲情形下的完整 $3\times3$ 矩阵 $\mathbf{RKD}$ 简洁得多，为后续的小推力优化奠定了基础。

4. 小推力碰撞规避动力学

4.1 动力学方程

令 $\mathbf{u} = (u_r, u_\theta, u_h)^T$ 为推力加速度矢量（在轨道坐标系的径向、横向、法向分量），其幅值恒定且远小于当地引力加速度：

$$\sqrt{u_r^2 + u_\theta^2 + u_h^2} = a_0 \ll \frac{\mu}{r_1^2}$$

在推力弧段内，考虑一个无限小的时间间隔 $\delta t$，推力产生的速度变化 $\delta\mathbf{v} = \mathbf{u}\delta t$ 引起 b 平面位移 $\delta\mathbf{b} = (\delta\xi, \delta\zeta)^T$。利用线性关系 [Bombardelli & Hernando-Ayuso, 2015]：

$$\frac{d\mathbf{b}}{dt} = \mathbf{M}_0(t) \mathbf{u}$$

其中 $\mathbf{M}_0(t) = \mathbf{R}_0(t) \mathbf{K}(t) \mathbf{D}(t)$ 是一个 $2\times3$ 的时变矩阵，将推力方向映射到 b 平面运动。

这是从脉冲到小推力的关键推广： 在脉冲情形下，$\Delta\mathbf{b} = \mathbf{M}\Delta\mathbf{v}$ 是一次性的代数关系；而在小推力情形下，这变成了一个微分方程——控制变量 $\mathbf{u}(t)$ 是一个随时间变化的函数，优化的复杂度从有限维提升到了无限维。

4.2 机动方案

论文考虑的 CAM 方案为：一个持续的推力弧段 $\Delta\theta_t > 0$（推力幅值恒定），后跟一个滑行弧段 $\Delta\theta_c \geq 0$。推力弧段的长度受限于电推进系统的功率约束（通常只能在日照段工作）。本文不涉及多弧段机动（可作为未来工作扩展）。

5. 切向推力的解析解

5.1 纯切向机动

在深入最优控制之前，作者首先给出了一个极具工程价值的简化解：推力方向始终保持切向（沿速度方向或反方向）。基于 Bombardelli 等人 [2011] 关于常切向推力下二体问题解析解的前期工作，推导出推力弧段结束后的 b 平面坐标：

$$\xi \simeq \xi_0 + 2a_0\frac{r_1^3}{\mu}\left[\Delta\theta_t + \sin(\Delta\theta_c) - \sin(\Delta\theta_t + \Delta\theta_c)\right]$$

$$\zeta \simeq \zeta_0 + a_0\frac{r_1^3}{\mu}\cos\frac{\kappa}{2}\left[3\Delta\theta_t + \Delta\theta_t\left(2 + \frac{\Delta\theta_c}{\Delta\theta_t}\right) - 8\sin\frac{\Delta\theta_t}{2}\sin\left(\frac{\Delta\theta_t}{2} + \Delta\theta_c\right)\right]$$

这两个解析公式的美妙之处在于：不需要求解任何微分方程或优化问题——给定机动弧长 $\Delta\theta_t$ 和滑行弧长 $\Delta\theta_c$，可以直接计算碰撞概率。这为任务初步设计提供了一个极快的评估工具。

5.2 切向解的特征

注意 $\xi$ 坐标中含有一个长期项（secular term）$2a_0(r_1^3/\mu)\Delta\theta_t$，这意味着随着推力弧长增加，沿 $\xi$（MOID 方向）的位移线性增长，这是提升碰撞距离的最有效方式。而 $\zeta$ 坐标同时包含长期分量和振荡分量。

6. 最优控制公式

6.1 Pontryagin 极大值原理的应用

将碰撞规避表述为 Bolza 型最优控制问题。目标函数和动力学分别为：

$$J = \frac{1}{2}\mathbf{b}^T\mathbf{Q}\mathbf{b}\bigg|_{t=t_f} + \int_{t_0}^{t_f} \boldsymbol{\lambda}^T(t)\mathbf{M}_0(t)\mathbf{u}\,dt$$

构建哈密顿函数 $H(\mathbf{u}, \boldsymbol{\lambda}, t) = \boldsymbol{\lambda}^T\mathbf{M}_0\mathbf{u}$，其中 $\boldsymbol{\lambda} = (\lambda_\xi, \lambda_\zeta)^T$ 为协态向量。

应用 Pontryagin 极大值原理：

1. 推力幅值取最大值 $a_0$（bang-bang 控制的连续推力版本）
2. 最优推力方向为：

$$\mathbf{u} = a_0\frac{\mathbf{M}_0^T\boldsymbol{\lambda}}{\|\mathbf{M}_0^T\boldsymbol{\lambda}\|}$$

1. 协态方程（由于哈密顿函数不显含 $\mathbf{b}$）：

$$\frac{d\boldsymbol{\lambda}}{dt} = -\frac{\partial H}{\partial \mathbf{b}} = 0 \quad\Rightarrow\quad \boldsymbol{\lambda}(t) = \text{常数}$$

1. 横截条件：$\boldsymbol{\lambda}(t_f) = \mathbf{Q}\mathbf{b}(t_f)$

6.2 降维为二维 TPBVP

这是本文最优控制推导的精华。 由于协态向量为常数，且与终端 b 向量通过 $\mathbf{Q}$ 矩阵线性关联，整个优化问题被简化为一个二维两点边值问题（TPBVP）：

$$\frac{d\mathbf{b}}{dt} = a_0\frac{\mathbf{M}_0\mathbf{M}_0^T\mathbf{Q}\mathbf{b}(t_f)}{\|\mathbf{M}_0^T\mathbf{Q}\mathbf{b}(t_f)\|}, \quad \mathbf{b}(t_0) = (\xi_0, \zeta_0)^T, \quad \mathbf{b}(t_f) = (\xi_f, \zeta_f)^T$$

与脉冲情形相比，这里的根本区别在于：
- 脉冲：$3\times3$ 矩阵的本征值问题（有限维）
- 小推力：二维微分方程的两点边值问题（无限维 → 数值求解）

但通过 b 平面动力学的线性特性，TPBVP 的维度从完整轨道动力学的 6 维（位置 + 速度）降至仅 2 维（b 平面坐标），使得数值求解非常高效。

6.3 最优推力方向的参数化

与脉冲情形类似，最优推力方向可以用面内角 $\gamma_{\text{in}}$ 和面外角 $\gamma_{\text{out}}$ 参数化：

$$\mathbf{u} = \begin{pmatrix} u_r \\ u_\theta \\ u_h \end{pmatrix} = \begin{pmatrix} a_0\cos\gamma_{\text{out}}\sin\gamma_{\text{in}} \\ a_0\cos\gamma_{\text{out}}\cos\gamma_{\text{in}} \\ a_0\sin\gamma_{\text{out}} \end{pmatrix}$$

这与 Bombardelli & Hernando-Ayuso (2015) 中的参数化完全一致，体现了两个工作的连续性和统一性。

7. 数值结果与讨论

论文基于当前 Starlink 星座的两行根数（TLE）信息，选取了一个典型场景：300 kg 航天器，标称推力 10 mN。以下按协方差椭圆朝向 $\Theta$ 分类呈现数值结果。$\Theta = 5^\circ$ 是论文的主要基准情形。

7.1 $\Theta = 5^\circ$ 基准情形：最优 vs 切向的大图景

论文首先以 $\Theta = 5^\circ$ 为基准，在推力弧中点对应的轨道圈数 $n_{\text{rev}}$ 和推力弧长 $\Delta\theta_t$ 的二维参数空间中，对比了最优解与切向解。

图2 展示了最优机动对应的碰撞概率（对数尺度），作为 $n_{\text{rev}}$ 和 $\Delta\theta_t$ 的函数。碰撞概率随 $n_{\text{rev}}$ 呈振荡变化，局部极小值出现在整圈和半圈附近（$n_{\text{rev}} = N + 0.5$ 处），这是由 $\zeta$ 坐标中的振荡分量所导致。图中用红色虚线标出了 $\Delta\theta_t = 100^\circ, 200^\circ, 300^\circ$ 三个截面，后续的线图即沿这些截面展开。右侧颜色条显示的是不同推力水平（2–30 mN）下碰撞概率的对数值，范围从约 $10^{-8.5}$ 到 $10^{-5}$。

图3 展示了纯切向机动对应的碰撞概率（对数尺度）。 与图2 采用相同的坐标系和颜色尺度。定性来看，切向解的等高线结构与最优解非常相似——振荡模式、局部极小值的位置基本一致。这初步表明在大部分参数空间，切向解是对最优解的良好近似。

图4 定量地展示了切向解相对于最优解的相对误差。 这是一张关键图：绿色区域表示切向解与最优解性能接近，而蓝色/深色区域表示切向解明显劣于最优解。可以看到，在大部分参数空间（绿色和黄色区域），相对误差在 $10^{-1}$ 以内（即切向解的碰撞概率最多比最优解高几十个百分点）。但在 $n_{\text{rev}} \in [0.5, 3]$ 且 $\Delta\theta_t > 100^\circ$ 的区域（蓝色），相对误差显著增大，切向解可能比最优解差一个数量级以上。

7.2 $\Theta = 5^\circ$, $\Delta\theta_t = 200^\circ$：深远机动的最优 vs 切向对比

图5（$\Theta = 5^\circ, \Delta\theta_t = 200^\circ$ 截面） 是论文的核心对比图。蓝色虚线和绿色点划线分别对应顺行和逆行的切向解，橙色实线为最优解，黑色虚线为无机动的基准碰撞概率（约 $4 \times 10^{-5}$）。图中标注了三个特征点 a、b、c——在这些位置，最优解与切向解的性能差距最为显著。

一个重要发现是：在大部分 $n_{\text{rev}}$ 范围内，最优解与逆行切向解几乎重合——这意味着对于 $\Theta = 5^\circ$，最优推力方向在大部分提前量下接近于逆行切向。但在点 a（$n_{\text{rev}} \approx 0.5-1.5$，即临交会机动区域）和点 c（$n_{\text{rev}} \approx 5.5$），最优解显著优于任何切向策略，碰撞概率可降低近一个数量级。

图6 展示了点 a 处最优推力的面内角 $\gamma_{\text{in}}$（蓝线）和面外角 $\gamma_{\text{out}}$（红线）在推力弧段内的演化。 横轴为推力弧段内的纬度辐角 $\theta$（从约 220° 到 460°），纵轴为角度（度）。可见在点 a（临交会机动），最优推力方向在整个弧段内剧烈变化：面内角从约 $-60^\circ$ 平滑过渡到 $+40^\circ$，面外角也在 $-50^\circ$ 到 $-10^\circ$ 之间波动。这解释了为什么切向解（固定方向）在此区域表现不佳——最优解需要连续大幅度调整推力方向。

图7 展示了点 b 处最优推力方向的演化。 点 b 位于 $n_{\text{rev}} \approx 2.5$——这是一个最优解与逆行切向解几乎完全重合的区域。如图所示，面内角（蓝线）在整个弧段内保持在约 $0^\circ$ 附近，面外角（红线）也维持在约 $0^\circ$——这意味着最优推力方向几乎就是纯逆行切向。这是切向解表现最好的参数区域。

图8 展示了点 c 处最优推力方向的演化。 点 c 位于 $n_{\text{rev}} \approx 5.5$——提前量很大。与点 a 类似，面内角（蓝线）在弧段内显著变化（约 $-150^\circ$ 到 $+30^\circ$），面外角（红线）也有约 $20^\circ$ 的波动。这种大幅度的方向变化再次使切向近似失效。

7.3 推力弧长的影响：$\Delta\theta_t = 120^\circ, 60^\circ, 30^\circ$

为了展示推力弧长对最优解与切向解差异程度的影响，论文给出了 $\Theta = 5^\circ$ 下不同 $\Delta\theta_t$ 的截面对比。

图9（$\Delta\theta_t = 120^\circ$）显示， 当推力弧长缩短至 120° 时，最优解（橙色实线）与逆行切向解（绿色点划线）的差异已经比图5（$\Delta\theta_t = 200^\circ$）显著缩小。在大部分 $n_{\text{rev}}$ 范围内，两者几乎无法区分。碰撞概率的整体水平也有所上升（约 $10^{-8}$ 到 $10^{-3}$），因为较短的推力弧段产生的 b 平面位移更小。

图10（$\Delta\theta_t = 60^\circ$）进一步展示了这一趋势。 标注了三个新特征点 d、e、f。最优解与逆行切向解在绝大多数 $n_{\text{rev}}$ 下几乎重合——只有在点 d（$n_{\text{rev}} \approx 0.5$，临交会）附近，最优解略优于切向解。这说明对于较短的推力弧段，切向解几乎等同于最优解。

图11（$\Delta\theta_t = 30^\circ$）是这一序列的极限。 在 30° 的短弧段下，最优解、顺行切向解和逆行切向解三条曲线几乎完全重合——此时小推力解已经平滑收敛到 Bombardelli & Hernando-Ayuso (2015) 的脉冲解。碰撞概率随 $n_{\text{rev}}$ 的变化也更为平缓，振荡幅度减小。

图12 展示了点 d 处最优推力的方向演化。 与点 a（200° 弧段）相比，面内角（蓝线）和面外角（红线）的变化幅度显著减小——面内角在约 $-140^\circ$ 到 $-110^\circ$ 的窄范围内变动，面外角几乎恒定在约 $20^\circ$。这反映了短弧段机动下，最优推力方向更接近常数。

图13 展示了点 e 处的方向演化。 点 e 位于 $n_{\text{rev}} \approx 2$ 附近。面内角（蓝线）在约 $0^\circ-120^\circ$ 之间变化，面外角（红线）保持在约 $4^\circ$ 附近——推力主要保持在轨道面内，接近切向。

图14 展示了点 f 处的方向演化。 点 f 位于 $n_{\text{rev}} \approx 3.5$ 附近。与点 e 类似，推力主要保持在轨道面内（面外角接近 0°），面内角在约 $-6^\circ$ 到 $+12^\circ$ 之间微小波动。这再次印证了切向近似在此参数区域的有效性。

7.4 $\Theta = 0^\circ$ 情形：协方差椭圆与 b 平面对齐

当协方差椭圆的主轴与 b 平面坐标轴完全对齐时（$\Theta = 0^\circ$），交会几何具有特殊性质。

图15 展示了 $\Theta = 0^\circ$ 时最优机动对应的碰撞概率。 与 $\Theta = 5^\circ$（图2）相比，一个显著区别是：碰撞概率在 $n_{\text{rev}}$ 上的振荡模式更规则——局部极小值整齐地排列在 $n_{\text{rev}} = N + 0.5$ 处，且整体碰撞概率水平更低（大部分区域在 $10^{-6}$ 到 $10^{-5}$ 之间），说明 $\Theta = 0^\circ$ 的协方差构型更有利于碰撞规避。

图16 展示了 $\Theta = 0^\circ$ 时切向机动的碰撞概率。 与图15 对比，切向解的等高线结构与最优解几乎完全一致——这一点将在图17 的相对误差图中得到定量确认。

图17 是理解 $\Theta = 0^\circ$ 情形最重要的图。 相对误差的对数值在整个参数空间几乎全部在 $-0.5$ 到 $0.5$ 之间（即相对误差在约 $0.3$ 到 $3$ 倍之间），大部分区域甚至为 0（白色/浅色区域）。这意味着当 $\Theta = 0^\circ$ 时，切向解与最优解几乎不可区分——在工程实践中，直接使用切向解析公式即可，无需运行数值优化。

图18（$\Theta = 0^\circ, \Delta\theta_t = 60^\circ$ 截面） 提供了线图视角的确认。最优解（橙色实线）与逆行切向解（绿色点划线）在全部 $n_{\text{rev}}$ 范围内几乎完美重合。标注了点 g，其方向演化见下图。

图19 展示了点 g 处的方向演化。 面内角（蓝线）在 $-150^\circ$ 到 $+180^\circ$ 之间变化——基本上覆盖了 360° 的全范围——但面外角（红线）始终保持在 $0^\circ$ 附近。这说明在 $\Theta = 0^\circ$ 时，最优推力完全在轨道面内，不需要面外分量。推力方向虽然随时间变化，但其效果等价于切向推力。

7.5 $\Theta = -5^\circ$ 情形：协方差椭圆反向倾斜

图20 展示了 $\Theta = -5^\circ$ 时最优机动对应的碰撞概率。 与 $\Theta = +5^\circ$（图2）相比，碰撞概率整体水平大幅降低——大部分区域在 $10^{-11}$ 到 $10^{-6}$ 之间，比 $\Theta = +5^\circ$ 的情形低了 2–3 个数量级。同时，振荡模式也更规则，局部极小值整齐排列。这揭示了一个重要现象：当协方差椭圆朝负方向倾斜时，切向推力的长期分量（沿 $\xi$ 方向增大碰撞距离）恰好沿着协方差椭圆的低概率方向，从而大幅提高规避效率。

图21 展示了 $\Theta = -5^\circ$ 时切向机动的碰撞概率。 与图20 对比，切向解仍然很好地捕捉了最优解的整体结构——碰撞概率水平、振荡模式、局部极小值位置都高度一致。

图22 定量展示了 $\Theta = -5^\circ$ 下切向解的相对误差。 与 $\Theta = 0^\circ$（图17）相比，误差略大——在 $n_{\text{rev}} \in [0.5, 2]$ 且 $\Delta\theta_t > 150^\circ$ 的区域，相对误差可达约 1–3。但整体而言，切向解仍提供了很好的近似，尤其在大 $n_{\text{rev}}$ 和小 $\Delta\theta_t$ 区域。

图23（$\Theta = -5^\circ, \Delta\theta_t = 60^\circ$ 截面） 显示，在短弧段下，最优解（橙色实线）与逆行切向解（绿色点划线）几乎不可区分。与 $\Theta = +5^\circ$ 的同参数图（图10）相比，碰撞概率整体低了约 4 个数量级——这再次凸显了协方差椭圆朝向对规避效率的决定性影响。标注了点 h。

图24 展示了点 h 处的方向演化。 面内角（蓝线）在约 $-150^\circ$ 到 $+180^\circ$ 之间变化，但面外角（红线）始终接近 $0^\circ$——与 $\Theta = 0^\circ$ 的情形类似，最优推力几乎完全在轨道面内。

7.6 扩展提前量分析：$\Theta = 0^\circ, 1^\circ, 3^\circ, 5^\circ$

为了研究协方差椭圆朝向在更大 $n_{\text{rev}}$ 范围内的影响，论文将分析扩展到了 $n_{\text{rev}} \in [0, 14]$。

图25（$\Theta = 0^\circ$，扩展 $n_{\text{rev}}$ 到 14） 显示碰撞概率在 $n_{\text{rev}} > 7$ 后持续振荡下降，最低可达约 $10^{-8}$。最优解、顺行切向解和逆行切向解三条曲线在整个范围内几乎完美重合，进一步印证了 $\Theta = 0^\circ$ 时切向解等同于最优解的结论。

图26（$\Theta = 1^\circ$） 是论文中最具启发性的图之一。仅仅 $1^\circ$ 的协方差椭圆倾斜，就使得最优解（橙色实线）与切向解之间出现了可辨识的差异：在 $n_{\text{rev}} \in [5, 9]$ 区间，切向解的局部极小值消失或移位，而最优解仍然保持较低的碰撞概率。在大 $n_{\text{rev}}$（>10）下，切向解的碰撞概率甚至可能比最优解高出两个数量级。这说明即使是最微小的协方差椭圆倾斜，也会对长提前量机动的策略选择产生显著影响。

图27（$\Theta = 3^\circ$） 进一步放大了 $\Theta = 1^\circ$ 中观察到的趋势。最优解与切向解的差异更加显著：在大 $n_{\text{rev}}$ 下，切向解几乎完全失效，而最优解仍能将碰撞概率维持在 $10^{-6}$ 以下。

图28（$\Theta = 5^\circ$，扩展 $n_{\text{rev}}$ 到 14） 展示了最大倾斜角的情形。在大 $n_{\text{rev}}$（>10）区域，切向解的碰撞概率不再下降，甚至可能出现回升，而最优解则持续下降。这清晰地表明：协方差椭圆朝向 $\Theta$ 是决定长提前量机动策略有效性的关键参数，而 $1^\circ$ 的微小倾斜就足以改变最优策略的性质。

7.7 b 平面轨迹可视化：理解切向 vs 最优的本质差异

图29 是理解最优解与切向解本质差异的关键可视化。 它在 $\xi$-$\zeta$ 平面上同时绘制了三条轨迹：
- 蓝色虚线：顺行切向解的 b 平面轨迹
- 绿色虚线：逆行切向解的 b 平面轨迹
- 彩色散点：最优解的 b 平面轨迹（颜色代表 $n_{\text{rev}}$，从黄色（0 圈）到深紫色（约 4.5 圈））

背景中的斜线是恒定碰撞概率的等高线（$P = 10^{-4}, 10^{-5}, 10^{-6}$），它们的方向由协方差椭圆的 $\Theta = 5^\circ$ 倾斜决定。

核心洞察：
- 逆行机动（绿色）在提前量较小（$n_{\text{rev}}$ 低）时更优——它的轨迹在 $\xi$ 增大方向上有较大位移，沿低概率等高线方向移动
- 顺行机动（蓝色）在提前量较大（$n_{\text{rev}}$ 高）时更优——它的 $\zeta$ 变化顺应协方差椭圆的倾斜方向
- 最优解（彩色散点）在两者之间切换：在低 $n_{\text{rev}}$ 时偏向逆行，在高 $n_{\text{rev}}$ 时偏向顺行

这张图直观地解释了为什么 $\Theta > 0^\circ$ 时，最优机动策略不是简单地"提前越早越好"——而是在不同的 $n_{\text{rev}}$ 区间，最优推力方向需要顺应协方差椭圆等高线的倾斜方向。

8. 精度验证：解析/半解析方法 vs 高精度数值传播

论文通过高精度数值传播（EGM96 重力场模型截断至 50 阶，包含日月和行星引力）对解析/半解析方法进行了严格验证。

图30 展示了切向解析公式的绝对误差。 四条曲线分别对应 $\xi$ 和 $\zeta$ 坐标在 Keplerian 和高精度力模型下的误差。Keplerian 情形下，绝对误差在 $10^{-4}$ m 量级——几乎可以忽略。高精度力模型下，绝对误差增大至约 $10^{-1}$ 到 $10^0$ m（即分米到米级），但在 LEO 交会距离通常为百米到公里级的背景下，这一误差完全可以接受。$\xi$ 的误差（蓝/绿色）略大于 $\zeta$ 的误差（红/橙色），这与 $\xi$ 坐标包含长期项（累积效应更大）的事实一致。

图31 展示了切向解析公式的相对误差。 Keplerian 情形下，相对误差在 $10^{-6}$ 到 $10^{-2}$ 之间。高精度力模型下，相对误差约在 $10^{-2}$ 到 $10^0$（即 1% 到 100%）之间——在 $\zeta \to 0$ 的奇点附近，相对误差较大，但此时绝对误差本身很小，因此在碰撞概率计算中不产生实际影响。

图32 展示了最优控制解的绝对误差。 与图30（切向解）相比，最优解在高精度力模型下的绝对误差略大——约 $10^0$ 到 $10^1$ m（米到十米级）。这是因为最优控制解涉及对 $\mathbf{M}_0$ 矩阵的多次运算，数值误差会因 TPBVP 的迭代求解而略有放大。但即便如此，10 m 量级的误差在实际操作中仍然可以接受。

图33 展示了最优控制解的相对误差。 与图31（切向解）的模式相似：Keplerian 情形下误差极小，高精度力模型下相对误差约在 $10^{-2}$ 到 $10^0$ 之间。$\zeta$ 坐标在奇点附近的相对误差较大，但绝对误差仍处于可接受范围。

精度验证的总体结论：
- Keplerian 情形： b 平面坐标的绝对误差不超过 0.01 m，相对误差可忽略
- 含摄动情形（LEO）： 绝对误差小于 10 m，相对误差约 1%——与脉冲论文 [Bombardelli & Hernando-Ayuso, 2015] 的结论一致
- J2 和其他环境摄动对 b 平面预测的影响在实际操作中可以忽略

9. 结论与启示

Hernando-Ayuso & Bombardelli (2021) 成功地将碰撞规避的最优控制理论从脉冲域推广到小推力域，其核心贡献可总结为：

1. 方法的连续性与统一性 — 论文使用与脉冲工作相同的 b 平面框架和推力方向参数化，保证了两个领域在理论和工程实践上的无缝衔接

2. 解析与数值的优雅结合 — 通过 Pontryagin 极大值原理将无限维最优控制问题降维为二维 TPBVP，同时提供切向推力的完全解析解，实现了亚毫秒级的快速评估

3. 切向解的工程实用性 — 在大部分参数空间（尤其是 $\Theta = 0^\circ$ 和最优提前量附近），切向解与最优解性能接近，为任务初步设计提供了一个无需优化的快速评估工具

4. 协方差朝向的关键作用 — 论文揭示了一个反直觉的现象：$1^\circ$ 的协方差椭圆倾斜就可能导致"提前机动不如临交会机动"的结论，这对星座运营策略有直接影响

5. 精度满足工程需求 — 高精度力模型下的误差不超过 10 m（绝对）和 1%（相对），验证了 b 平面线性理论的工程适用性

对工程实践的启示：

五年前，Bombardelli & Hernando-Ayuso (2015) 将脉冲碰撞规避从"数值扫描"变为"解析求解"；六年后，Hernando-Ayuso & Bombardelli (2021) 又将这一范式延伸到了小推力时代。随着数万颗电推进星座卫星的部署，这篇论文提供的方法论将成为碰撞规避操作的核心工具。

如果说第一篇论文回答了"往哪个方向打一个脉冲最有效"，那么这篇论文回答的是"如何在有限推力弧段内连续调节推力方向，使碰撞概率最小"——前者是一个代数问题，后者是一个控制问题。两者的结合，构成了现代碰撞规避理论的完整图景。

参考文献

1. Bombardelli, C., and Hernando-Ayuso, J., "Optimal Impulsive Collision Avoidance in Low Earth Orbit", J. Guidance, Control, and Dynamics, Vol. 38, No. 2, 2015, pp. 217–225. DOI: 10.2514/1.G000742

2. Hernando-Ayuso, J., and Bombardelli, C., "Low-Thrust Collision Avoidance in Circular Orbits", J. Guidance, Control, and Dynamics, Vol. 44, No. 5, 2021, pp. 984–995. DOI: 10.2514/1.G005547

3. Bombardelli, C., Baù, G., and Peláez, J., "Asymptotic Solution for the Two-Body Problem with Constant Tangential Thrust Acceleration", Celestial Mechanics and Dynamical Astronomy, Vol. 110, No. 3, 2011, pp. 239–256.

4. Bombardelli, C., "Analytical Formulation of Impulsive Collision Avoidance Dynamics", Celestial Mechanics and Dynamical Astronomy, Vol. 118, No. 2, 2014, pp. 99–114.

5. Chan, F. K., Spacecraft Collision Probability, Aerospace Press, El Segundo, CA, 2008.

6. Salemme, G., Armellin, R., and Di Lizia, P., "Continuous-Thrust Collision Avoidance Manoeuvres Optimization", AIAA SciTech 2020 Forum, AIAA Paper 2020-0231.

7. Bryson, A. E., and Ho, Y.-C., Applied Optimal Control, Taylor & Francis, 1975.

8. García-Pelayo, R., and Hernando-Ayuso, J., "Series for Collision Probability in Short-Encounter Model", J. Guidance, Control, and Dynamics, Vol. 39, No. 8, 2016, pp. 1908–1916.