近地轨道碰撞的最优脉冲机动

论文信息： Claudio Bombardelli 和 Javier Hernando-Ayuso, "Optimal Impulsive Collision Avoidance in Low Earth Orbit", Journal of Guidance, Control, and Dynamics, Vol. 38, No. 2, 2015. DOI: 10.2514/1.G000742

研究机构： 马德里理工大学 (Universidad Politécnica de Madrid)

1. 引言

随着地球轨道空间物体数量的持续增长，航天器碰撞规避机动（Collision Avoidance Maneuver, CAM）已成为在轨运营中不可或缺的环节。截至2025年，地球轨道上尺寸大于10厘米的编目碎片已超过40,000个，每年发布的碰撞预警数量以万计。每一次预测的近距离交会（conjunction）都可能需要航天器执行轨道机动以降低碰撞风险。

碰撞规避机动的核心问题可以表述为：在给定的脉冲速度增量（$\Delta v$）约束下，如何选择机动方向以最大化规避效果？ 这个问题的难点在于：

1. 轨道动力学高度非线性
2. 机动方向是一个三维连续空间
3. 碰撞概率计算涉及高维积分
4. 每次交会场景的轨道参数各不相同

传统方法通常依赖数值传播和参数扫描，计算耗时长（分钟到小时级），无法满足快速决策需求。Bombardelli 和 Hernando-Ayuso 在2015年发表的这篇论文，将最优规避机动问题解析化，求解时间从分钟级降至毫秒级，是空间碎片领域的里程碑式工作。

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

2.1 b平面坐标系

为描述两个空间物体的近距离交会，论文采用 b平面（b-plane）坐标系 $\langle \xi, \eta, \zeta \rangle$。这个坐标系以目标航天器 S2 为中心：

$$ \mathbf{u}_\xi = \frac{\mathbf{v}_2 \times \mathbf{v}_1}{\| \mathbf{v}_2 \times \mathbf{v}_1 \|} $$

$$ \mathbf{u}_\eta = \frac{\mathbf{v}_1 - \mathbf{v}_2}{\| \mathbf{v}_1 - \mathbf{v}_2 \|} $$

$$ \mathbf{u}_\zeta = \mathbf{u}_\xi \times \mathbf{u}_\eta $$

其中 $\mathbf{v}_1$ 和 $\mathbf{v}_2$ 分别是机动航天器 S1 和目标 S2 的速度矢量，$\mathbf{u}_\eta$ 沿相对速度方向，b平面（$\xi$-$\zeta$ 平面）则是垂直于相对速度的平面。交会时刻两物体的最接近点投影到b平面上的位置为 $(\xi_e, 0, \zeta_e)$。

2.2 碰撞概率的解析计算

在短时交会假设（short-term encounter hypothesis）下，两物体的相对运动可视为匀速直线运动。碰撞概率的计算简化为 b 平面上的二维积分：

$$ P = \iint_A \frac{1}{2\pi\sigma_\xi\sigma_\zeta\sqrt{1-\rho_{\xi\zeta}^2}} \exp\left[-\frac{1}{2(1-\rho_{\xi\zeta}^2)}\left(\frac{\xi-\xi_e}{\sigma_\xi}\right)^2 + \left(\frac{\zeta-\zeta_e}{\sigma_\zeta}\right)^2 - 2\rho_{\xi\zeta}\frac{(\zeta-\zeta_e)(\xi-\xi_e)}{\sigma_\zeta\sigma_\xi}\right] d\xi d\zeta $$

其中 $A$ 是以 S1 联合包络半径 $s_A$ 为半径的圆形区域，$\sigma_\xi, \sigma_\zeta, \rho_{\xi\zeta}$ 从 b 平面下的相对位置协方差矩阵 $C_{\xi\zeta}$ 中提取。

论文采用 Chan 方法 [Chan, 2008] 将上述积分简化为 Rician 积分：

$$ P(u, v) = e^{-v/2} \sum_{m=0}^{\infty} \frac{v^m}{2^m m!} \left[1 - e^{-u/2} \sum_{k=0}^{m} \frac{u^k}{2^k k!}\right] $$

其中：

$$ u = \frac{s_A^2}{\sigma_\xi\sigma_\zeta \sqrt{1-\rho_{\xi\zeta}^2}} $$

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

$u$ 是碰撞截面面积与 1-$\sigma$ 协方差椭圆面积之比，$v$ 是"入侵深度"（depth of intrusion）的平方 [Lázaro & Righetti, 2012]。

3. 核心突破：线性模型

3.1 b平面位移与脉冲的线性关系

论文的核心创新之一是基于 Bombardelli [2013] 前期工作的 线性化模型。在机动点处施加脉冲 $\Delta v = (\Delta v_r, \Delta v_\theta, \Delta v_h)^T$（径向、横向、法向分量）后，b平面上的相对位置变化为：

$$ \mathbf{r} = \mathbf{R} \mathbf{K} \mathbf{D} \Delta v \equiv \mathbf{M} \Delta v $$

其中三个矩阵分别描述几何变换、轨道运动学变换和机动位置变换。

3.2 各矩阵的物理意义

矩阵 $\mathbf{R}$（几何旋转矩阵）将机动脉冲从 S1 的轨道坐标系变换到 b 平面：

$$ \mathbf{R} = \begin{bmatrix} 0 & 0 & -1 \\ \cos\beta & -\sin\beta & 0 \\ -\sin\beta & -\cos\beta & 0 \end{bmatrix} $$

其中 $\beta$ 是 $\mathbf{v}_1$ 与 $\mathbf{v}_1 - \mathbf{v}_2$ 之间的夹角。

矩阵 $\mathbf{K}$ 包含轨道参数（半长轴 $a_0$、偏心率 $e_0$、真近点角 $\theta_c$ 等）和交会几何（倾角差 $\phi$、轨道面夹角 $\psi$），是连接脉冲方向与b平面位移的关键。

矩阵 $\mathbf{D}$ 包含从机动时刻 $\theta_m$ 到交会时刻 $\theta_c$ 的轨道传播，其元素 $d_{tr}, d_{t\theta}, d_{rr}, d_{r\theta}, d_{wh}$ 是 $e_0, \theta_c, \theta_m$ 的非线性函数（见论文 Appendix A）。

4. 最大碰撞距离机动

4.1 问题建模

对于 正碰撞（direct impact，即 $r_e = 0$）的情况，最大碰撞距离优化问题可表述为：

$$ \begin{aligned} \text{最大化} \quad & J_r(\Delta v) = \xi^2 + \zeta^2 \\ \text{约束} \quad & f(\Delta v) = \Delta v^T \Delta v - \Delta v_0^2 \leq 0 \end{aligned} $$

改写为矩阵形式：

$$ J_r = \mathbf{r}^T \mathbf{Q} \mathbf{r} = \Delta v^T \mathbf{A} \Delta v $$

其中 $\mathbf{Q} = \text{diag}(1, 0, 1)$ 只选取 $\xi$ 和 $\zeta$ 分量（b平面），$\mathbf{A} = \mathbf{M}^T \mathbf{Q} \mathbf{M}$。

4.2 本征值问题的导出

使用拉格朗日乘子法：

$$ \mathcal{L}(\Delta v, \lambda) = J_r - \lambda f $$

令 $\partial\mathcal{L}/\partial\Delta v = 0$，得到：

$$ \mathbf{A} \Delta v = \lambda \Delta v $$

这就是一个标准的本征值问题。 最优机动方向就是矩阵 $\mathbf{A}$ 的最大本征值 $\lambda_1$ 对应的本征向量 $\mathbf{s}_1$。

最优脉冲为：

$$ \Delta v_{\text{opt}} = \Delta v_0 \mathbf{s}_1 $$

对应的最大碰撞距离为：

$$ \Delta r_{\max} = \sqrt{\lambda_1} \Delta v_0 $$

关键发现：$\text{rank}(\mathbf{A})$ 要么为1（当 $\theta_c - \theta_m = 2\pi n$ 时），要么为2（一般情况下）。这意味着总存在至少一个脉冲方向使碰撞距离保持不变（零本征值方向）。

5. 最小碰撞概率机动

对于最小化碰撞概率，目标函数变为：

$$ \tilde{J}_P = v = \frac{\xi^2}{\sigma_\xi^2} + \frac{\zeta^2}{\sigma_\zeta^2} - 2\rho_{\xi\zeta} \frac{\xi\zeta}{\sigma_\xi\sigma_\zeta} \bigg/ (1-\rho_{\xi\zeta}^2) $$

将其写作二次型：

$$ \tilde{J}_P = \mathbf{r}^T \mathbf{Q}' \mathbf{r} $$

其中：

$$ \mathbf{Q}' = \frac{1}{1-\rho_{\xi\zeta}^2} \begin{bmatrix} 1/\sigma_\xi^2 & 0 & -\rho_{\xi\zeta}/\sigma_\xi\sigma_\zeta \\ 0 & 0 & 0 \\ -\rho_{\xi\zeta}/\sigma_\xi\sigma_\zeta & 0 & 1/\sigma_\zeta^2 \end{bmatrix} $$

代入 $\mathbf{r} = \mathbf{M}\Delta v$ 得：

$$ \tilde{J}_P = \Delta v^T \mathbf{A}' \Delta v $$

同样转化为本征值问题 $\mathbf{A}'\Delta v = \lambda \Delta v$，最优解对应最大本征值的本征向量。

最大碰撞距离与最小碰撞概率两种准则在形式上完全等价，只是二次型矩阵不同。

6. 非正碰撞的一般情况

6.1 问题推广

当两物体并非正碰撞（即预期最接近距离 $r_e \neq 0$）时，b平面相对位置为：

$$ \mathbf{r} = \mathbf{r}_e + \mathbf{M} \Delta v $$

碰撞距离平方变为：

$$ J_r = J_{r0} + \Delta v^T \mathbf{A} \Delta v + 2 \mathbf{r}_e^T \mathbf{Q} \mathbf{M} \Delta v $$

这是一个 非凸二次优化问题。论文巧妙地将其转化为凸问题。

6.2 凸化处理

令 $\mathbf{u} = \Delta v / \Delta v_0$，$\mathbf{b}^T = \mathbf{r}_e^T \mathbf{Q} \mathbf{M} / \Delta v_0$，原问题变为：

$$ \begin{aligned} \text{最大化} \quad & \tilde{J}_r(\mathbf{u}) = \mathbf{u}^T \mathbf{A} \mathbf{u} + 2\mathbf{b}^T \mathbf{u} \\ \text{约束} \quad & \mathbf{u}^T \mathbf{u} - 1 \leq 0 \end{aligned} $$

根据 Boyd [2004] 的凸优化理论，可进一步转化为：

$$ \text{最小化} \quad \frac{(\mathbf{s}_1^T \mathbf{b})^2}{\lambda - \lambda_1} + \frac{(\mathbf{s}_2^T \mathbf{b})^2}{\lambda - \lambda_2} + \lambda \quad \text{约束: } \lambda \geq \lambda_1 $$

其中 $\lambda_1 \geq \lambda_2$ 是 $\mathbf{A}$ 的两个非零本征值，$\mathbf{s}_1, \mathbf{s}_2$ 是对应的本征向量。

6.3 求解步骤

步骤 1： 求解非线性方程（使用牛顿法）：

$$ \left(\frac{\mathbf{s}_1^T \mathbf{b}}{\lambda - \lambda_1}\right)^2 + \left(\frac{\mathbf{s}_2^T \mathbf{b}}{\lambda - \lambda_2}\right)^2 - 1 = 0, \quad \lambda \geq \lambda_1 $$

步骤 2： 计算最优脉冲：

$$ \Delta v_{\text{opt}} = -\Delta v_0 (\mathbf{A} - \lambda_{\text{opt}} \mathbf{I})^\dagger \mathbf{b} $$

其中 $\dagger$ 表示伪逆矩阵。

步骤 3： 代入式 (22) 得到最大碰撞距离。

6.4 最优解的几何特性

论文证明，最优解必须满足：

$$ \mathbf{r}_e^T \mathbf{Q} (\mathbf{r}_{\text{opt}} - \mathbf{r}_e) \geq 0 $$

这一条件将可行解限制在由初始b平面相对位置 $\mathbf{r}_e$ 的垂线界定的半空间中。当机动角距 $\Delta \theta$ 变化时，最优解可能出现不连续性（discontinuity），这对应于可达域（reachable domain）椭圆与约束边界相切的情形。

7. 算例验证 I：2009年铱星-宇宙号碰撞

7.1 交会几何

2009年2月10日，活跃的 铱星33（Iridium 33）与失效的 宇宙2251（Cosmos 2251）卫星在西伯利亚上空约 788.6 km 高度相撞，产生了超过 2000 块可追踪碎片。这是历史上首次两颗完整卫星的意外碰撞。

论文以 Iridium 作为机动航天器 S1，主要参数：

协方差矩阵（b平面）：

$$ C_{\xi\zeta} = \begin{bmatrix} 0.02 & 0 \\ 0 & 0.80 \end{bmatrix} \text{km}^2 $$

非正碰撞假设：$\xi_e = -70$ m, $\zeta_e = -70$ m。

7.2 最优机动方向

图1 展示了 Iridium–Cosmos 交会的最优机动方向。 图1(a) 为面内转角 $\sigma$，图1(b) 为面外转角 $\gamma$。两条曲线对比了最小化碰撞概率（$P_{\min}$，虚线）与最大化碰撞距离（$r_{\max}$，实线）两种准则。

关键发现：
- 当 $\Delta\theta$ 较小时（即临近交会时才机动），$\sigma$ 和 $\gamma$ 均远离零，说明最优脉冲既不是纯切向也不是纯法向，而是一个组合方向
- 两种准则在小 $\Delta\theta$ 下差异显著，切不可混为一谈
- 当 $\Delta\theta \to \infty$（提前多个轨道周期机动），两种准则趋于一致

7.3 碰撞距离与概率对比

图2 展示了两种准则下的机动效果对比。 图2(a) 显示最大碰撞距离（$r_{\max}$）准则在 $\Delta\theta > 2\pi$ 后可达 3000 m 以上，而最小概率（$P_{\min}$）准则略低。图2(b) 显示碰撞概率随 $\Delta\theta$ 增大呈周期性振荡下降。

7.4 b平面可达域

图3 的 b 平面地图展示了最优解在 b 平面的轨迹。 恒定概率等高线呈椭圆形，两种准则在 $\xi$-$\zeta$ 平面上遵循不同的路径。

尤为重要的是图4，它揭示了最优解的不连续性：

图4 展示了当 $\Delta\theta$ 从 $100^\circ$ 增加到 $170^\circ$ 时，可达域（灰色椭圆区域）的演变。 当 $\Delta\theta \approx 136^\circ$ 时，可达域椭圆的长轴恰好与约束边界相切，最优解发生不连续跳跃。这不是数值误差，而是问题的内在几何特性。

8. 算例验证 II：RapidEye4–UoSat2 交会

8.1 交会几何

2013年5月，RapidEye4 卫星执行了一次真实的碰撞规避机动，与失效的 UoSat2 卫星交会。主要参数：

8.2 给定碰撞概率下的最优机动

图5 展示了为达到指定碰撞概率（$10^{-6}, 10^{-8}, 10^{-10}$）所需的最优机动方向。 图5(a) 的面内角 $\sigma$ 在 $\Delta\theta$ 接近 1 轨道周期时出现不连续跳跃（从约 $-90^\circ$ 突跳到约 $180^\circ$），这对应着从顺行机动到逆行机动的切换。图5(b) 的面外角 $\gamma$ 同样在相同位置出现特征性尖峰。

图6 显示达到指定碰撞概率所需的 $\Delta v$ 大小。 关键结果：
- 要求越严格（$P$ 越小），所需 $\Delta v$ 越大
- 当机动时机接近交会（$\Delta\theta \to 0$），所需 $\Delta v$ 急剧增加
- 对于 $P = 10^{-6}$，最优机动 $\Delta v < 0.02$ m/s 即可，极其高效

9. 方法精度评估

9.1 J2 摄动影响

论文采用 Keplerian 轨道模型，但评估了最重要的 J2 地球非球形摄动的影响。

图7(a) 显示对于 Iridium-Cosmos 实际场景，J2 摄动引起的碰撞距离误差不超过 0.25%， 对于工程应用完全可以忽略。图7(b) 展示了对应的碰撞距离曲线。

图8 针对更极端的近正碰撞场景（$\psi = 1^\circ$，倾角 $45^\circ$ 以最大化 J2 效应）。 图8(a) 显示误差峰值可达 9%，但这些峰值恰好对应 规划不良的机动（图8b 中碰撞距离很小的区域）。在实际应用中，操作者不会选择这些低效的机动方案，因此 J2 误差在实际操作中可以忽略。

论文进一步评估了大气阻力影响。即使是在极端情况（350 km 高度近正碰撞），只要机动提前量不太大（数圈以内），大气阻力误差仍然很小。

9.2 与其他方法的比较

与完全数值算法相比，该解析方法的优势：

10. 结论与启示

Bombardelli & Hernando-Ayuso (2015) 的这篇论文通过巧妙的数学变换，将复杂的碰撞规避机动优化问题简化为：

1. 本征值问题 — 求解 $3\times 3$ 矩阵的本征值和本征向量
2. 单变量非线性方程 — 用牛顿法轻松求解

这一简化使得最优碰撞规避机动的计算时间从 分钟级降至亚毫秒级，可以轻松嵌入实时碰撞规避决策系统中。

对工程实践的启示：

1. 临近交会机动方向不直观 — 最优脉冲不是简单的"沿速度方向加速"，而是一个包含面外分量的组合方向，需要通过本征值分析精确确定

2. 最小碰撞概率 ≠ 最大碰撞距离 — 两种准则在小 $\Delta\theta$ 下存在显著差异，应根据任务风险偏好选择合适的优化目标

3. 解析模型的精度足够工程使用 — J2 和大气阻力等主要摄动的引入误差不超过 1%（合理规划场景下），完全可以替代耗时的数值优化

4. 已有开源工具 — 论文作者开发的 OCCAM（Optimal Computation of Collision Avoidance Maneuver）工具已提供在线演示，可供实际操作和教学使用

这篇论文代表了空间碎片碰撞规避领域从"数值试错"到"解析求解"的范式转变，是每一个从事空间碎片、轨道力学和航天器运营的研究者和工程师的必读经典。

参考文献

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. Bombardelli, C., "Analytical Formulation of Impulsive Collision Avoidance Dynamics", Celestial Mechanics and Dynamical Astronomy, Vol. 118, No. 2, 2013, pp. 99-114.

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

4. Klinkrad, H., Space Debris: Models and Risk Analysis, Springer–Verlag, New York, 2006.

5. Sanchez-Ortiz, N., et al., "Collision Risk Assessment and Avoidance Manoeuvres — The New CORAM Tool for ESA", 64th IAC, 2013.

6. Boyd, S. P., and Vandenberghe, L., Convex Optimization, Cambridge Univ. Press, 2004.

7. Lázaro, D., and Righetti, P., "Evolution of EUMETSAT LEO Conjunctions Events Handling Operations", 12th SpaceOps, 2012.

8. Valsecchi, G., et al., "Resonant Returns to Close Approaches: Analytical Theory", Astronomy and Astrophysics, Vol. 408, No. 3, 2003, pp. 1179-1196.