简介：兰伯特转移用于求解航天器在两固定位置间的最优转移轨道，是轨道设计与任务规划中的基础算法。该压缩包内含1个MATLAB脚本lambert.m，面向天体力学与航天轨道计算学习者，以及需要快速验算转移参数的工程师。脚本通过输入起始/目标位置、转移时间等参数，解算顺/逆时针转移的初速与终速，并给出飞行时间、升交点/降交点坐标及总冲量等结果，可辅助评估任务可行性并降低燃料消耗。资源包以单个m文件为主，整包约2KB，无需安装依赖，直接在MATLAB中运行即可，适合作为理解兰伯特问题的入门工具或嵌入个人计算流程。目前已有1986人浏览/学习，脚本轻量、结构清晰，便于对照算法公式梳理输入输出逻辑，并在此基础上扩展多圈转移或摄动修正，可显著缩短轨道转移算法的上手时间。 搞轨道计算的人，几乎都会遇到兰伯特问题。简单说：你手里有个飞行器，在 t1 时刻位于 r1，我希望它在 t2 时刻到达另一个位置 r2，那么它应该飞一条什么样的轨道？这就是兰伯特问题（Lambert's problem），求解出来的轨道叫兰伯特转移轨道（Lambert transfer orbit），很多时候也直接叫兰伯特轨道。它几乎是所有轨道转移、交会对接、行星际飞行任务设计里绕不开的基础模块，也是很多刚入门的小伙伴被劝退的第一个硬骨头。这篇文章我想用尽量不绕弯的方式，把兰伯特问题背后的原理、求解套路和实测踩坑经验整理出来，希望对正在啃轨道力学的你有帮助。
1. 兰伯特问题到底在解决什么
1.1 教科书里的经典定义
兰伯特问题的标准说法是：在中心引力场中，给定引力常数 μ，给定两个位置矢量 r1 和 r2，再给定从 r1 飞到 r2 所花的飞行时间 Δt，要求解出一条满足这些条件的圆锥曲线轨道，并给出轨道在两端的速度矢量 v1 和 v2。
这个定义看起来很短，但信息量很大。它不关心你出发前在一条什么轨道上，也不关心你到达后要进入哪条轨道，只关心“从 r1 到 r2，花 Δt，走哪条弧段”。这正好是很多任务设计的核心问题：你打算让航天器在某个时刻从当前位置出发，在指定时刻出现在目标位置，中间这段轨道怎么走。
1.2 为什么航天工程离不开兰伯特转移
最简单的答案就是：真实任务里很少有两个目标刚好共面、共圆心、还要求转 180 度的理想条件。卫星交会对接时，追踪星和目标星可能不在同一轨道面；深空探测器进行行星际转移时，地球和火星的位置关系一直在变；甚至每一次中段轨道修正，都要重新计算一条从当前位置到目标点的转移轨道。这些情况都能化解成“已知两个位置和飞行时间，求轨道”的兰伯特问题。
所以兰伯特问题不是某一类特殊问题的名字，而是一大类轨道机动问题的公共底座。无论是近地轨道的相位调整、远距离拦截追踪，还是行星际探测器的发射窗口设计，底层都要反复调用兰伯特求解器。
1.3 它和霍曼转移是什么关系
很多教材先讲霍曼转移，再讲兰伯特问题，容易让人以为两者是并列关系。实际上霍曼转移是兰伯特问题的一个特例：当两个轨道都是圆轨道、共面、且转移角固定为 180 度时，飞行时间被半长轴唯一确定，于是可以直接写出解析解。但真实任务中很难满足这些条件，比如你可能要从椭圆轨道移到另一个非共面的椭圆轨道，或者要求在一个精确的时间点完成转移，这时候霍曼公式就不够用了。
| 对比项 | 霍曼转移 | 兰伯特转移 |
|---|---|---|
| 轨道关系 | 两圆轨道共面 | 任意椭圆、任意倾角 |
| 转移角 | 固定 180° | 任意角度，0~360° 都可 |
| 飞行时间 | 由半长轴唯一决定 | 由用户指定 |
| 求解方式 | 解析公式 | 数值迭代 |
| 典型应用 | 理论教学、简单轨道设计 | 交会、深空探测、变轨规划 |
2. 兰伯特定理：核心一把钥匙
2.1 兰伯特定理在说什么
18 世纪兰伯特提出了一个关键定理：在中心引力作用下，航天器从 r1 飞到 r2 的飞行时间，只取决于轨道的半长轴 a、两个位置矢量的模之和 r1 + r2，以及两位置之间的弦长 c，与轨道的偏心率无关。
这句话初看很反直觉。两条完全不同的椭圆轨道，偏心率差很多，但只要 a、r1 + r2、c 三个参数一致，飞行时间就完全一致。正是这个定理，把“给定两点和飞行时间求轨道”的问题变成了一个可以数值迭代的问题：我们只需要找到合适的半长轴，使对应的飞行时间等于给定值，剩下的轨道形状会自动确定。
用普适变量法表达时，飞行时间方程写为：
sqrt(μ) * Δt = χ³ * S(z) + A * sqrt(y)
其中 χ 是普适变量，z = χ² / a，S(z) 是 Stumpff 函数，y 是中间变量。这个公式把椭圆、抛物线、双曲线统一在一个方程里，避免了解题时不停分类讨论。
2.2 从几何参数理解“转移角”
转移角 Δν 是从 r1 矢量转到 r2 矢量沿轨道运动方向扫过的角度。数学上通过夹角公式：
cos(Δν) = (r1 · r2) / (r1 * r2)
可以直接算出角度，但这里有个坑：arccos 的结果范围是 0 ~ 180°，它只告诉你 r1 和 r2 之间的最小夹角，无法区分“从 r1 顺时针转过去”还是“逆时针转过去”，也无法表达超过 180° 的长弧转移。
打个比方，从北京飞上海，可以一路向东南飞，也可以绕地球转个大半圈再从另一个方向到达。两个方案都满足起点和终点，但飞行时间和燃料消耗完全不同。兰伯特问题的“转移方向”必须明确指定，否则解出来的轨道可能完全不是你想要的那条。
2.3 短程、长程与多圈
根据转移角 Δν 的大小，兰伯特转移通常分为三种常见情况：
- 短程（short-way）：Δν 在 0 ~ 180° 之间，转移时间通常较短，是工程中最常见的选择。
- 长程（long-way）：Δν 在 180° ~ 360° 之间，相当于绕个大弯，时间更长，适合某些能耗约束场景。
- 多圈（multi-revolution）：在短程或长程基础上，再额外绕中心天体一整圈或多圈，总转移角为 Δν + 2πN。多圈问题比单圈复杂得多，因为同样的飞行时间可能对应多条不同轨道，设计时需要额外处理根的选择。
在后面的普适变量法中，短程和长程通常用一个符号参数 dm 区分：dm = +1 表示短程，dm = -1 表示长程。很多新手在这里搞反，导致算出来的速度矢量方向完全错误。
3. 用普适变量法实现兰伯特求解
3.1 为什么不用解析法，而用迭代
兰伯特问题很难写出封闭解析解，因为它本质上是开普勒方程的一种推广，最终都会落在一个超越方程上。既然无法直接反解，就退而求其次：先假设一个轨道参数，算出飞行时间，再和目标时间比较，不断修正参数直到收敛。
传统做法是按椭圆、抛物线、双曲线分别推导公式，但这会导致程序里充满 if-else，维护起来很痛苦。普适变量法的优势在于用 Stumpff 函数统一了三种圆锥曲线，一套代码走天下。工程上更稳健的还有 Gooding 算法、p-迭代法等，但我建议初学者先啃懂普适变量法，因为它能把计算逻辑讲得很清楚。
3.2 普适变量法的四个关键公式
第一，定义 Stumpff 函数：
C(z) = (1 - cos(√z)) / z（z > 0）
S(z) = (√z - sin(√z)) / (√z)³（z > 0）
当 z < 0 时，用双曲函数替换；z 接近 0 时取级数展开的前几项。这两个函数的作用是把椭圆和双曲线的情况统一成一套表达式。
第二，计算辅助量 A：
A = dm * sqrt(r1 * r2 * (1 + cos(Δν)))
这里 dm 控制短程或长程，A 的符号决定了转移方向。
第三，计算中间变量 y：
y = r1 + r2 + A * (z * S(z) - 1) / sqrt(C(z))
y 必须大于 0，否则对应的轨道几何不存在，迭代就要退出或调整搜索区间。
第四，可用的飞行时间方程：
sqrt(μ) * Δt = χ³ * S(z) + A * sqrt(y)
其中 χ = sqrt(y / C(z))。对给定的 Δt，我们只需要找 z 使这个方程成立，然后用 Lagrange 系数 f、g 回代出 v1、v2。
3.3 可直接运行的 Python 演示代码
下面是教学演示用的兰伯特求解器，用普适变量法加二分求根实现。为了简洁，我只确保短程、单圈、椭圆轨道下可用；工程项目请改用更稳健的 Gooding 算法。
import numpy as np from scipy.optimize import brentq def stumpff_c(z): if z > 1e-8: return (1.0 - np.cos(np.sqrt(z))) / z if z < -1e-8: return (np.cosh(np.sqrt(-z)) - 1.0) / (-z) return 0.5 def stumpff_s(z): if z > 1e-8: sz = np.sqrt(z) return (sz - np.sin(sz)) / (sz**3) if z < -1e-8: sz = np.sqrt(-z) return (np.sinh(sz) - sz) / (sz**3) return 1.0 / 6.0 def lambert_uni(r1_vec, r2_vec, dt, mu, dm=1): """ r1_vec, r2_vec: 两个位置矢量，长度单位 km dt: 飞行时间，单位 s mu: 中心天体引力常数，单位 km^3/s^2 dm: +1 短程，-1 长程（本演示仅保证短程） """ r1 = np.linalg.norm(r1_vec) r2 = np.linalg.norm(r2_vec) cos_nu = np.dot(r1_vec, r2_vec) / (r1 * r2) cos_nu = np.clip(cos_nu, -1.0, 1.0) dnu = np.arccos(cos_nu) A = dm * np.sqrt(r1 * r2 * (1.0 + cos_nu)) def tof(z): c = stumpff_c(z) s = stumpff_s(z) y = r1 + r2 + A * (z * s - 1.0) / np.sqrt(c) if y < 0.0: return np.nan chi = np.sqrt(y / c) return (chi**3 * s + A * np.sqrt(y)) / np.sqrt(mu) # 教学演示：只在 z∈(1e-6, 10) 内找根，覆盖大多数短程单圈椭圆转移 fa = tof(1e-6) - dt fb = tof(10.0) - dt if np.isnan(fa) or np.isnan(fb): raise ValueError("函数在端点无定义") if fa * fb > 0: raise ValueError("区间两端符号相同，无法求根，请扩大搜索范围") z_root = brentq(lambda z: tof(z) - dt, 1e-6, 10.0, xtol=1e-12) c = stumpff_c(z_root) s = stumpff_s(z_root) y = r1 + r2 + A * (z_root * s - 1.0) / np.sqrt(c) chi = np.sqrt(y / c) f_lag = 1.0 - y / r1 g_lag = A * np.sqrt(y / mu) gdot = 1.0 - y / r2 v1 = (r2_vec - f_lag * r1_vec) / g_lag v2 = (gdot * r2_vec - r1_vec) / g_lag return v1, v2这段代码的逻辑不复杂：先计算辅助量 A，然后定义飞行时间函数tof(z)，用brentq找 z 使飞行时间等于给定 Δt，最后用 Lagrange 系数算出两个端点的速度。注意我这里刻意把搜索区间限制在 z ∈ (1e-6, 10)，为什么可以这样做？因为对大多数短程椭圆转移，z 不会太大；如果飞行时间特别长，或者涉及双曲线轨道，就需要重新规划搜索区间。
4. 实战：设计一条地球到火星的兰伯特转移轨道
4.1 输入条件与单位处理
假设地球和火星都在同一个平面内的近圆轨道上，忽略真实星历和相位关系，只为了演示计算流程。设出发时刻地球位置为 r1 = [1 AU, 0, 0]，到达时刻火星位置为 r2 = [0, 1.524 AU, 0]，也就是两个位置之间的夹角是 90 度。设转移时间 Δt = 200 天。
单位必须统一。轨道力学里常见的长度单位是 km，时间单位是 s。所以先把 AU 转成 km：
AU = 1.495978707e8 # 1 AU = 1.495978707亿 km r1 = np.array([AU, 0.0, 0.0]) r2 = np.array([0.0, 1.524 * AU, 0.0]) mu_sun = 1.32712440018e11 # 太阳引力常数，km^3/s^2 dt = 200 * 86400 # 200 天，转换为秒然后用前面写的函数求速度：
v1, v2 = lambert_uni(r1, r2, dt, mu_sun, dm=1) print("出发速度 v1 =", v1) print("到达速度 v2 =", v2)跑通之后，你会得到两个三维速度矢量，单位是 km/s。这里的 v1 是航天器在 r1 处为了进入兰伯特转移轨道需要具备的速度，v2 是它到达 r2 时的速度。如果后续要对接火星轨道，还需要拿 v2 和火星轨道速度做差，算机动速度增量。
4.2 从速度结果反推轨道半长轴
得到 v1 后，可以立刻用能量方程检查这条转移轨道是椭圆还是双曲线。能量方程：
ε = v²/2 - μ/r = -μ / (2a)
所以半长轴：
a = -μ / (2ε)
对应代码如下：
def semi_major_axis(r_vec, v_vec, mu): r = np.linalg.norm(r_vec) v = np.linalg.norm(v_vec) epsilon = v**2 / 2.0 - mu / r return -mu / (2.0 * epsilon) a_transfer = semi_major_axis(r1, v1, mu_sun) print("转移轨道半长轴 a =", a_transfer)200 天转移比霍曼转移（约 258.5 天）更短，因此需要更大能量，半长轴会小于霍曼转移的 1.262 AU，偏心率也会明显偏离 0。你如果打印出来，会发现 a 是一个小于 1.262 AU 的正值，这正是“快转移”的典型特征。
4.3 用数值积分验证转移时间
兰伯特求解器算完就完事了吗？建议务必验证一遍。最直接的方法，是用初始状态 (r1, v1) 在太阳引力场里做二体数值积分，看积分到 Δt 时刻时航天器是否真的到达 r2。
from scipy.integrate import solve_ivp def two_body_rhs(t, state, mu): r_vec = state[:3] v_vec = state[3:] r = np.linalg.norm(r_vec) a = -mu * r_vec / r**3 return np.concatenate([v_vec, a]) state0 = np.concatenate([r1, v1]) sol = solve_ivp(two_body_rhs, [0, dt], state0, args=(mu_sun,), rtol=1e-9) diff = sol.y[:3, -1] - r2 print("终点位置误差 =", diff)如果误差在公里量级甚至更小，说明求解器和积分器没问题。如果误差大到离谱，多半是兰伯特求解代码里 z 的区间、符号或者单位出了问题。这一步是检验算法正确性的黄金标准，强烈建议每个人在做完兰伯特计算后都跑一遍。
5. 写给新手的避坑指南
5.1 转移角一定要先明确“走哪边”
我见过不少初次接触兰伯特问题的人，包括我自己早期，都会在转移角上翻车。arccos 只会给出 0 ~ 180 度的小角，但实际任务中你可能需要走大于 180 度的大弧。如果不先画图确认短程还是长程，直接用默认的 dm=1，算出来的轨道很可能把航天器送到目标相反方向。
正确做法：拿到 r1 和 r2 后，先画一张几何示意图，标出中心天体位置、两个位置矢量、以及你计划沿着哪个方向飞行。确认转移角小于 180 度用 dm=1，大于 180 度用 dm=-1。一句话，算法可以不知道你的意图，但你必须知道自己要往哪边飞。
5.2 迭代不收敛时先查什么
兰伯特迭代不收敛，十有八九是下面几个原因。
第一，搜索区间没选对。演示代码里我把 z 限制在 (1e-6, 10)，但真实问题如果飞行时间极长，或者轨道能量很高，z 的真实根可能落在区间外。遇到这种情况，先扩大区间，或者观察 tof(z) 随 z 的变化曲线，手动确认根的大致范围。
第二，y < 0。它表示在当前 z 下对应的几何构型不存在。很多实现会直接让函数返回异常值，但如果异常值混入求根算法，很容易导致诡异行为。稳妥的办法是在迭代循环里显式排除 y < 0 的情况。
第三，单位不一致。μ 用 km³/s²，位置矢量却用了 AU，时间用了小时，算出来的结果自然是一堆天文数字。写代码时强制规定一套单位体系，所有输入都先转换好。
5.3 结果合不合理的快速检查方法
兰伯特求解完成后，有三个快速检查手段。
一是能量检查：用 v²/2 - μ/r 判断轨道类型，如果半长轴为负就是双曲线，为正就是椭圆。对比任务预期是否一致。
二是方向检查：v1 和目标点方向应该大致顺向，不会是反向飞。你可以画图确认 v1 确实指向 r2 方向。
三是数值积分验证：用 solve_ivp 或者任何靠谱积分器，把 (r1, v1) 推演到 Δt，看终点是否落在 r2 附近。这个检查最可信，也是我每次做完兰伯特计算后必做的动作。
最后说点我自己的感受。刚接触兰伯特问题时，我觉得它就是一道数学题，直到有一次做任务仿真，因为转移角方向判断错误，生成了一条与目标卫星背道而驰的轨道，我才意识到算法细节背后的几何直觉有多重要。如果你也卡在代码里，别硬啃公式，先画一张“从 r1 到 r2 怎么走”的示意图，再回来调参，很多坑会好踩很多。
本文还有配套的精品资源，点击获取
