小行星 Bennu 最小控制代价软着陆

背景

小天体着陆并不是按比例缩小的行星着陆。这里的引力极弱，离心加速度可能与引力处于同一量级；在惯性空间中看似平缓的推力过程，换到天体固连坐标系后还会因科里奥利项产生明显横向漂移。因此，用天体固连系表达模型，或把惯性系模型严格转换到固连系，有助于直接处理相对于表面的约束。

本例规划探测器在小行星 (101955) Bennu 赤道平面内的下降轨迹。探测器从一个相对 Bennu 静止的近域位置出发，在 \(2.5\) 小时后以规定的 \(2\ \mathrm{mm/s}\) 向内法向速度首次接触表面。有限而很小的接触速度使事件定义清楚且可以直接验证：接触前轨迹必须位于天体外部，终端速度则必须具有给定的向内分量。因此，这里求解的是软着陆，而不是零相对速度的悬停会合。

它与 Earth-to-Venus 示例的物理问题并不重复。本例研究的是数百米尺度、天体固连旋转系中的表面接近过程，还要满足局部着陆走廊；另一个示例则是行星之间的日心轨道转移。

Bennu 数据与任务假设

球形模型采用 OSIRIS-REx 研究在 2019 年任务早期给出的 Bennu 整体参数估计，而不是把它们作为当前的最佳拟合值。参数经过适当取整以便复现；这些数值不表示球形引力能够替代 Bennu 的实测不规则形状，也不包含完整的参数协方差。

Bennu 参数	符号	数值
平均半径	\(R\)	\(245.03\ \mathrm{m}\)
引力参数	\(\mu\)	\(4.892\ \mathrm{m^3/s^2}\)
恒星自转周期	\(P\)	\(4.296057\ \mathrm{h}\)
自转角速度	\(\omega=2\pi/P\)	\(4.063\times10^{-4}\ \mathrm{rad/s}\)

下表参数是为了构造示例而明确给出的任务设计条件，不是 Bennu 的固有属性，也不是对某一实际飞行器性能的描述。

任务参数	数值
飞行时间 \(T\)	\(9{,}000\ \mathrm{s}=2.5\ \mathrm{h}\)
最大指令加速度 \(a_{\max}\)	\(2.0\times10^{-4}\ \mathrm{m/s^2}\)
最小下滑角 \(\gamma\)	\(20^\circ\)
向内法向接触速度 \(v_c\)	\(2.0\times10^{-3}\ \mathrm{m/s}\)

旋转坐标系动力学

坐标原点位于 Bennu 质心，\(x\) 轴穿过目标着陆点，\(y\) 轴位于赤道平面内；坐标系以恒定角速度 \(\boldsymbol\Omega=\omega\boldsymbol e_z\) 随 Bennu 转动。平面状态与推力加速度指令定义为

\[ \boldsymbol r=(x,y)^{\mathsf T},\qquad \boldsymbol v=(v_x,v_y)^{\mathsf T},\qquad \boldsymbol a_T=(a_x,a_y)^{\mathsf T}. \]

采用球对称引力势后，天体固连系中的方程为

\[ \dot{\boldsymbol r}=\boldsymbol v, \]

\[ \dot{\boldsymbol v} =-\frac{\mu}{\lVert\boldsymbol r\rVert^3}\boldsymbol r -2\boldsymbol\Omega\times\boldsymbol v -\boldsymbol\Omega\times (\boldsymbol\Omega\times\boldsymbol r) +\boldsymbol a_T. \]

展开为分量形式：

\[ \begin{aligned} \dot x &= v_x, & \dot y &= v_y,\\ \dot v_x &= -\frac{\mu x}{(x^2+y^2)^{3/2}} +2\omega v_y+\omega^2x+a_x,\\ \dot v_y &= -\frac{\mu y}{(x^2+y^2)^{3/2}} -2\omega v_x+\omega^2y+a_y. \end{aligned} \]

各项符号可由旋转坐标系中的矢量二阶求导直接得到。其中 \(+\omega^2(x,y)\) 是离心加速度，不能误当作附加引力项。

边界条件、目标与约束

初始和终端条件为

\[ \begin{aligned} \boldsymbol r(0)&=(539.066,-122.515)^{\mathsf T}\ \mathrm{m} =R(2.20,-0.50)^{\mathsf T},\\ \boldsymbol v(0)&=(0,0)^{\mathsf T}\ \mathrm{m/s},\\ \boldsymbol r(T)&=(R,0)^{\mathsf T},\\ \boldsymbol v(T)&=(-0.002,0)^{\mathsf T}\ \mathrm{m/s}. \end{aligned} \]

着陆点处的 \(x\) 轴正向指向天体外部，因此终端 \(v_x<0\) 表示规定的向内接触速度。固定时域目标函数为

\[ \min_{\boldsymbol a_T} J =\int_0^T\lVert\boldsymbol a_T(t)\rVert_2^2\,\mathrm dt. \]

这一 \(L^2\) 控制代价有利于得到平滑指令，适合用于制导研究，但它不等于推进剂消耗。脚本还单独报告

\[ \Delta v_{\mathrm{proxy}} =\int_0^T\lVert\boldsymbol a_T(t)\rVert_2\,\mathrm dt, \]

作为便于解释的积分加速度指标。若要求真正的最小推进剂解，还必须加入飞行器质量、推力、比冲和质量消耗方程。

推力加速度圆盘与球面避碰约束为

\[ \lVert\boldsymbol a_T\rVert_2\le a_{\max}, \qquad \lVert\boldsymbol r\rVert_2\ge R. \]

在目标点附近定义相对于局部切平面的高度 \(h=x-R\)。下滑走廊要求轨迹始终位于切平面上方，并随着接触时刻临近将横向误差收紧至零：

\[ h\ge\tan\gamma\,|y|, \qquad \gamma=20^\circ. \]

非线性规划内部实际采用 \(21^\circ\)，预留一度数值裕量；随后在 4,001 个时刻按规定的 \(20^\circ\) 约束检查重构轨迹。这个裕量用于抵消多项式插值在节点之间的微小偏移，并没有改变最终报告的工程要求。

尺度化与离散

小天体弱引力会使 SI 制下的方程系数相差很大，因此优化模型采用 Bennu 半径和自然引力时间进行尺度化：

\[ t_* = \sqrt{\frac{R^3}{\mu}}=1734.146\ \mathrm{s},\qquad v_* = \frac{R}{t_*}=0.141297\ \mathrm{m/s},\qquad a_* = \frac{\mu}{R^2}=8.14794\times10^{-5}\ \mathrm{m/s^2}. \]

定义

\[ \bar{\boldsymbol r}=\frac{\boldsymbol r}{R},\quad \bar t=\frac{t}{t_*},\quad \bar{\boldsymbol v}=\frac{\boldsymbol v}{v_*},\quad \bar{\boldsymbol a}=\frac{\boldsymbol a_T}{a_*}, \]

则引力项的量级为一。相应的无量纲参数为 \(\bar\omega=0.704519\)、\(\bar T=5.189874\) 和 \(\bar a_{\max}=2.454608\)。

默认转录采用 48 个区间、每个区间 3 个 Lobatto 点；--quick 使用 24 个区间。初始猜测是一条同时满足两端位置和速度的三次 Hermite 曲线，再把该曲线代入旋转系动力学反算控制，为优化器提供具有动力学依据的初值。

独立验证

优化器返回成功并不足以证明连续轨迹可用。脚本在求解后执行两层检查：

1. 在 4,001 个均匀时刻重构状态和控制，检查球面、\(20^\circ\) 下滑走廊、加速度边界以及终端条件。
2. 独立重构每一段二次控制多项式，并使用 SciPy 的 DOP853 方法重新积分原始旋转系常微分方程，相对容差设为 \(2\times10^{-11}\)。积分轨迹会与配点多项式比较，也会重新检查表面和着陆走廊。

在独立积分之前，局部分段二次求值器还会在每个网格区间的四分之一点和四分之三点与 Pockit 插值矩阵对比。这样可以避免某个未经抽查的区间在无提示的情况下积分了不同的控制多项式。

运行示例

在仓库根目录运行经过完整验证的默认算例：

python -m examples.asteroid_soft_landing

使用较小网格进行快速冒烟测试：

python -m examples.asteroid_soft_landing --quick

无界面运行并保存图片：

python -m examples.asteroid_soft_landing --save asteroid-soft-landing.png --no-show

程序保持与其他示例一致的调用结构：

system, phase = build_problem()
guess = initial_guess(phase)
solution = solve_problem(system, guess)
plot_solution(solution)

验证结果

默认 \(48\times3\) 网格在验证机器上得到：

指标	验证值
物理平方加速度目标 \(J\)	\(2.220023066\times10^{-5}\ \mathrm{m^2/s^3}\)
积分加速度指标	\(0.385324\ \mathrm{m/s}\)
最大指令加速度	\(0.115282\ \mathrm{mm/s^2}\)
加速度上限	\(0.200000\ \mathrm{mm/s^2}\)
密集重构的表面约束违反量	\(0\ \mathrm{m}\)
密集重构的 \(20^\circ\) 走廊违反量	\(0\ \mathrm{m}\)
独立积分终端位置误差	\(2.49\times10^{-4}\ \mathrm{m}\)
独立积分终端速度误差	\(5.30\times10^{-8}\ \mathrm{m/s}\)
独立积分最小离表高度	\(2.3\times10^{-5}\ \mathrm{m}\)
独立积分最小下滑走廊裕量	\(-6.7\times10^{-5}\ \mathrm{m}\)
配点与独立积分最大尺度化差值	\(3.282\times10^{-5}\)

独立积分轨迹对下滑走廊的最大偏移仅为 \(0.067\ \mathrm{mm}\)，小于脚本设置的 \(0.5\ \mathrm{mm}\) 积分容差，而且轨迹仍位于球形表面之外。峰值指令约使用加速度上限的 58%。脚本直接检查完整终端速度向量，而不是从图中推断接触速度。

适用范围与局限

该模型保留了真实小天体着陆问题中的关键最优控制结构，但不能直接用于 Bennu 的工程着陆设计。它假设球对称引力、匀速自转、赤道平面运动、质点飞行器和理想的加速度指令跟踪；没有考虑 Bennu 的多面体形状与不规则引力、姿态和质量动力学、太阳光压、导航误差、地形风险、羽流与表面的相互作用、执行器带宽及接触力学。实际任务必须替换为经过验证的三维形状与引力模型，并加入鲁棒制导和具体飞行器约束。

数据来源

• Barnouin, O. S. et al., "Shape of (101955) Bennu indicative of a rubble pile with internal stiffness," Nature Geoscience 12, 247-252 (2019), doi:10.1038/s41561-019-0330-x。
• Scheeres, D. J. et al., "The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements," Nature Astronomy 3, 352-361 (2019), doi:10.1038/s41550-019-0721-3。

源代码

完整可运行示例见： examples/asteroid_soft_landing.py。