简介：兰伯特转移问题在航天工程中应用广泛，对应MATLAB求解脚本面向轨道设计、星际任务分析与天体力学教学场景，可帮助工程师和学生快速获得两体假设下的转移轨道参数。脚本基于牛顿第二定律与开普勒定律，计算双曲型转移轨道，确定航天器在两个位置间的初始速度、最终速度及转移时间，并支持顺行与逆行两种转移方式，便于方案对比与工况分析。压缩包内含一个m文件，整体仅2KB，代码精简、结构清晰，无需复杂配置即可直接运行或嵌入其他仿真流程。已有2005人浏览学习，具备良好参考价值。通过输入起止位置、转移时间、航天器质量等参数，脚本可输出升交点坐标、降交点坐标、飞行时间及所需总冲量等关键结果，可用于近地轨道变轨、行星际转移窗口分析以及燃料消耗估算，同时也有助于理解兰伯特问题的数值求解过程。

1. 兰伯特转移：给两个位置和一个时间，轨道怎么自己长出来？

做轨道交会或者深空转移设计时，最常碰到的需求是：追踪航天器当前在某个位置，目标航天器未来某个时刻会出现在另一个位置，发动机在中途点一次火，能不能刚好在指定时刻落到目标点上。这个问题抽象成数学表述，就是兰伯特转移——已知中心天体引力参数 μ、初始位置矢量 r1、目标位置矢量 r2 和转移时间 Δt，求解唯一连接这两点的开普勒轨道，并给出点火速度 v1 和到达速度 v2。兰伯特问题不是一道课后习题，它是交会对接、导弹拦截、行星际转移轨道设计的底层求解器。本文从兰伯特定理的几何条件开始，给出一套可复现的普适变量法解算代码，再讲清楚短程/长程、多圈解、归一化这些绕不开的参数，最后把 180 度奇异、初值漂移、单位混用这几个坑逐个排掉。

2. 从兰伯特定理到可运行算法：几何条件与短程、长程、多圈解怎么选

2.1 兰伯特定理的几何本质：三个几何量决定转移时间

兰伯特定理的表述很干净：在同一个中心引力场中，从 r1 到 r2 的椭圆转移时间只取决于三个几何量——半长轴 a、弦长 c = |r1 - r2| 以及两个位置矢量模的和 r1 + r2，与轨道的偏心率 e 无关。这意味着，轨道形状可以千变万化，但只要半长轴相同、两个端点位置确定，转移时间就相同。这个结论的反面才是工程上真正有用的形式：给定 r1、r2 和 Δt，反解出半长轴 a，再从 a 恢复出整条轨道和两端速度。

这带来一个直接推论：兰伯特问题本质上是一个单变量求根问题。弦长 c 和 r1 + r2 是输入条件直接算出来的几何量，唯一需要迭代求解的变量是半长轴 a。但直接用 a 做迭代变量并不明智，地球轨道 a 的量级是 7000 km，木星转移轨道的 a 可能是上亿公里，不同任务差四个量级以上，迭代初值和收敛域都会很难看。工程上常用普适变量 z 做中间变量，z 与 a 的关系在归一化单位下是 a = 1/z，z > 0 对应椭圆轨道，z = 0 是抛物线边界，z < 0 是双曲线轨道。这个变量把三种轨道类型统一在一个连续区间里，同时把 a 的巨大取值范围压缩到 z ∈ (-10, 10) 附近的经验收敛区，这是普适变量法能成为工程主流做法的根本原因。

2.2 解的族：短程、长程与多圈解，选错一个结果全错

给定 r1、r2 和一个转移时间 Δt，兰伯特问题的解并不唯一。几何上，从 r1 到 r2 的转移弧段可以沿小于 180 度的方向走，也可以沿大于 180 度的方向绕大圈，分别叫短程解和长程解。工程上默认取短程，因为转移时间短、燃料省，但交会任务里有时候目标星相位差大于 180 度，短程解的转移时间会非常短，速度增量反而大，这时候长程解反而是更优选择。这个选择在代码里通过一个 direction 标志位控制，short = -1，long = +1。

多圈解更隐蔽。当 Δt 足够长，转移轨道可以绕中心天体一圈甚至多圈后再到达 r2。转移时间越长，可行解的数量越多——每增加一圈，短程和长程各多一组解。实际交会任务为了等目标星相位或者节省燃料，经常会用 N=1 或 N=2 的多圈解。多圈解在数学上表现为 z 的迭代区间整体平移，初始圈数 N 需要作为输入参数给定，且 z 的收敛区间限制在 (2Nπ)² 附近。很多初写兰伯特解算器的人在这里翻车：没给 N 却跑出多圈解，或者给了 N=2 却收敛回 N=0 的短程解，这就是 2.3 节要讲的算法选型和初值策略问题。

2.3 Gooding 算法与普适变量法：不同场景下的选型理由

兰伯特问题的数值解法主要有三支：Gauss/Battin 经典迭代法、Gooding 算法和普适变量法（Universal Variables）。Gauss 法是最早的解析-数值混合方法，数学优雅但对初值极度敏感，r1 与 r2 夹角较大或偏心率接近 1 时经常发散，现代工程代码里已经很少直接使用。Gooding 算法是英国皇家飞机研究院的 R. H. Gooding 在 1988 年和 1990 年两篇论文里提出的，核心是用一个对初值不敏感的分数函数 plus 一套精心构造的初值估计公式，鲁棒性极强，是飞行软件和商业软件（如 GMAT）的常见选择。缺点是实现复杂，初值估计公式有一堆经验系数，抄错一个就全盘崩溃。

普适变量法是我个人最推荐先写一遍的算法。它的逻辑链条最短——一个 z 变量统一三种轨道类型，Stumpff 函数用级数展开，牛顿迭代求根，全部代码在一百行以内。对做任务设计验证来说够用，对理解兰伯特问题数学结构来说也是最清晰的路径。普适变量法对初值敏感的问题可以通过一个非常简单的策略缓解：z 的初值从抛物线边界 z0 = 0 起步，配合迭代过程中的括号法（bracketing）约束 z 的取值区间，就能覆盖绝大多数工程场景。下面的代码就是按这个思路实现的。

3. 用普适变量法跑通最小兰伯特解算器：代码框架与迭代参数

3.1 输入归一化：先把距离、时间、引力参数变成无量纲量

写兰伯特解算器的第一个动作不是写求解函数，而是处理输入单位。直接拿 km、s、km³/s² 混算，迭代变量 z 的取值范围会变得非常奇怪，收敛半径和容差都不好设。归一化的做法是选一个参考距离 r0（习惯取 r1 的模），定义距离单位 DU = r0，时间单位 TU = sqrt(r0³/μ)，速度单位 VU = DU / TU。归一化之后，μ = 1，所有长度和时间都变成无量纲量，z 的迭代区间稳定在经验值上。这样同一套代码换到月球、火星或者日心转移，只需要改 r0 和 μ 两个输入，不用改任何迭代参数。

import numpy as np def normalize(r1_vec, r2_vec, dt, mu): r0 = np.linalg.norm(r1_vec) # 参考距离，取初始位置矢量模 v0 = np.sqrt(mu / r0) # 参考速度，对应圆轨道速度 t0 = r0 / v0 # 参考时间，对应特征时间 r1_n = r1_vec / r0 # 归一化位置矢量 r2_n = r2_vec / r0 dt_n = dt / t0 # 归一化转移时间 return r1_n, r2_n, dt_n, r0, v0, t0

这个函数做四件事：把两组位置矢量按 r0 缩放到无量纲量，把转移时间按特征时间缩放，同时把 r0、v0、t0 三个换算因子保存下来供最后一步恢复到物理单位。注意 v0 的定义是 sqrt(μ / r0)，它刚好是 r0 处的圆轨道速度，不是随便选的，这样归一化后速度的量级落在 0.5~3 之间，后面恢复 f、g 系数时数值条件最好。如果你只是做单次计算，不恢复物理单位，那这个函数可以省掉一半，但工程上最后一定是要拿到真实速度增量去算燃料的，所以换算因子从一开始就要留好。

3.2 Stumpff 函数与核心迭代：用普适变量 z 统一椭圆与双曲线

兰伯特解算器的核心是一个关于 z 的非线性方程。标准形式是：

F(z) = ( y(z) / C(z) )^(3/2) · S(z) + A · sqrt( y(z) ) − Δt = 0

其中 S(z) 和 C(z) 是 Stumpff 函数，y(z) 是中间变量，A 由 r1、r2 的几何关系和转移方向决定。Stumpff 函数本质上是把开普勒方程里的三角函数展开成对 z 的幂级数，这样 z 取负值（双曲线）时仍然有定义，不需要分支判断。z 接近 0 时级数收敛慢，需要单独处理，标准做法是 z 的绝对值小于 1e-7 时用 z/6 − z²/120 这样的截断展开，防止对 0 取整导致精度崩塌。

def stumpff(z): # Stumpff 函数 C(z) 和 S(z)，z > 0 对应椭圆，z < 0 对应双曲线 if z > 1e-8: sqrt_z = np.sqrt(z) c = (1.0 - np.cos(sqrt_z)) / z s = (sqrt_z - np.sin(sqrt_z)) / (z * sqrt_z) elif z < -1e-8: sqrt_neg = np.sqrt(-z) c = (np.cosh(sqrt_neg) - 1.0) / (-z) s = (np.sinh(sqrt_neg) - sqrt_neg) / (-z * sqrt_neg) else: # z 接近于零时用级数截断，避免除零误差 c = 1.0 / 6.0 - z / 120.0 + z**2 / 5040.0 s = 1.0 / 6.0 - z / 120.0 + z**2 / 5040.0 return c, s

这段代码要特别注意 z 接近 0 的分支。很多人直接写 cos(z)/z，z 一掉到 1e-10 量级结果就变成 NaN，迭代直接崩掉。工程上 z 的迭代经常会在抛物线边界附近来回震荡，这个分支是必须写的。Stumpff 函数的级数展开在 z 绝对值大于 50 之后高次项会开始放大误差，所以后面迭代时要对 z 做区间约束，不允许它跑到 ±50 之外——真跑到那个量级，多半是初值给错了，不是物理上有解。

3.3 z 的初值选择与收敛判据：怎样让牛顿迭代不漂移

迭代主循环用牛顿法，但裸的牛顿法在兰伯特问题上容易漂移。我的做法是初值 z0 = 0（抛物线边界），然后用一个括号保护：每次迭代后检查新的 z 是否落在区间 [z_min, z_max] 内，超出去就换用二分法回退一步。这个策略牺牲一点收敛速度，换来的是即使给了很离谱的转移时间，程序也能稳定返回一个解而不是 NaN。收敛判据用相对误差，F(z) 的绝对值小于 1e-10 就算收敛，同时校验 y(z) > 0——y(z) 必须为正，否则轨道不存在。

def solve_lambert_z(r1_n, r2_n, dt_n, direction): r1_mag = np.linalg.norm(r1_n) r2_mag = np.linalg.norm(r2_n) cos_dnu = np.dot(r1_n, r2_n) / (r1_mag * r2_mag) cos_dnu = np.clip(cos_dnu, -1.0, 1.0) if direction == -1: # 短程解：转移角小于 180 度 dnu = np.arccos(cos_dnu) else: # 长程解：转移角大于 180 度 dnu = 2.0 * np.pi - np.arccos(cos_dnu) # A 是兰伯特几何常数，由 r1、r2 的模和转移角唯一确定 A = direction * np.sqrt(r1_mag * r2_mag) * np.sin(dnu) / np.sqrt(1.0 - cos_dnu) def F(z): c, s = stumpff(z) y = r1_mag + r2_mag + A * (z * s - 1.0) / np.sqrt(c) if y < 0: return -100.0, y # 非物理域，返回大负值让迭代回退 F_val = (y / c)**1.5 * s + A * np.sqrt(y) - dt_n return F_val, y z = 0.0 # 从抛物线边界起步 for _ in range(50): c, s = stumpff(z) F_val, y = F(z) if abs(F_val) < 1e-10 and y > 0: return z, y # 数值微分求导，比解析求导少一个推导出错的可能 dz = 1e-7 F_prime = (F(z + dz)[0] - F_val) / dz z_new = z - F_val / F_prime # 括号保护：新值超出经验区间就退回二分 if z_new < -8.0 or z_new > 8.0: z_new = 0.5 * z z = z_new raise RuntimeError("lambert iteration did not converge")

代码里的 A 有一个重要边界：当 dnu 接近 π（即 r1 和 r2 几乎反向共线）时，sin(dnu) 和 sqrt(1 - cos_dnu) 同时趋向 0，但比值趋于定值，数值上会非常不稳定。这是兰伯特问题著名的 180 度奇异点，后面第 5 章专门讲怎么处理。迭代过程中的经验区间 [-8, 8] 不是随便拍的，z = 8 对应椭圆半长轴 a = 0.125，已经是极高偏心率轨道，再大就是接近抛物线的双曲线，工程任务里几乎不会出现，这个区间足够覆盖 99% 的转移场景。

4. 速度恢复与轨道验证：从 z 到点火速度，再验证解不是数学幻觉

4.1 用拉格朗日系数 f、g 恢复两端速度

迭代收敛到 z 之后，下一步是把 z 转成实际的速度矢量 v1 和 v2。做法是先算出 y(z)，再通过拉格朗日系数 f、g 建立 r1、v1 到 r2、v2 的转移关系。f 和 g 的定义是 r2 = f · r1 + g · v1，v2 = f_dot · r1 + g_dot · v1，其中 f_dot 和 g_dot 是 f、g 对时间的导数。用 z 表示的公式高度集中：f = 1 − y / r1_mag，g = A · sqrt(y)，f_dot 和 g_dot 由 f、g 与 r2 的模推出来。恢复物理单位时乘回 v0 即可。

def lambert_solve(r1_vec, r2_vec, dt, mu): r1_n, r2_n, dt_n, r0, v0, t0 = normalize(r1_vec, r2_vec, dt, mu) z, y = solve_lambert_z(r1_n, r2_n, dt_n, direction=-1) r1_mag = np.linalg.norm(r1_n) r2_mag = np.linalg.norm(r2_n) # 拉格朗日系数，全部为无量纲量 f = 1.0 - y / r1_mag g = A_geometric(r1_n, r2_n, direction=-1) * np.sqrt(y) v1_n = (r2_n - f * r1_n) / g v2_n = 1.0 / g * (r2_n - r1_n) # 备用公式，与 f、g_dot 等价 return v1_n * v0, v2_n * v0, f, g

这段代码里 v2_n 用了简化的备用公式，省去 f_dot、g_dot 的显式计算。速度恢复的精度直接取决于 z 的收敛精度，这也是为什么前面迭代要卡到 1e-10——如果你只卡到 1e-6，恢复出来的速度误差会在线性传播中被放大，最终位置误差可能是几公里量级。第 4.2 节的验证方法能把这个误差直接暴露出来。

4.2 验证方法：用二体传播把解出来的轨道推回去

解算器输出一个速度矢量，但验证是硬碰硬的——把 v1 当作初始速度做二体传播 Δt，看传播结果是否精确等于 r2。这一步必须做，因为兰伯特解算器本质上是迭代求根，迭代收敛只说明 F(z) 接近零，不直接说明轨道在物理上闭合。二体传播用 f、g 系数法：给定 r1、v1、Δt，计算对应的偏近点角变化量，再合成新的位置。写一个最小传播函数即可。

def propagate(r1_vec, v1_vec, dt, mu): # 二体传播：从初始位置和速度推 dt 秒后的位置 r_mag = np.linalg.norm(r1_vec) v_mag = np.linalg.norm(v1_vec) energy = v_mag**2 / 2.0 - mu / r_mag a = -mu / (2.0 * energy) # 计算初始偏近点角 h_vec = np.cross(r1_vec, v1_vec) h_mag = np.linalg.norm(h_vec) e_vec = np.cross(v1_vec, h_vec) / mu - r1_vec / r_mag e_mag = np.linalg.norm(e_vec) cos_E0 = 1.0 / e_mag * (1.0 - r_mag / a) sin_E0 = np.dot(r1_vec, v1_vec) / (np.sqrt(a * mu) * e_mag) E0 = np.arctan2(sin_E0, cos_E0) # 解开普勒方程 E - e sinE = M M0 = E0 - e_mag * np.sin(E0) M_target = M0 + np.sqrt(mu / a**3) * dt E = M_target for _ in range(30): E = E - (E - e_mag * np.sin(E) - M_target) / (1.0 - e_mag * np.cos(E)) # 用 E 合成传播后的位置 f = 1.0 - a / r_mag * (1.0 - np.cos(E - E0)) g = dt + np.sqrt(a**3 / mu) * (np.sin(E - E0) - (E - E0)) r_new = f * r1_vec + g * v1_vec return r_new

这个 propagation 函数的逻辑完整覆盖了从轨道要素到开普勒方程求解的全过程。注意到它先用能量算半长轴 a，用角动量和离心率矢量算偏近点角 E0，再通过开普勒方程推进到目标时刻的 E，最后用拉格朗日系数合成新位置。对比传播结果和 r2 的误差，如果位置差超过 1 米，基本可以断定解算器或输入有问题。我一般会做一个快速自检：随机生成 100 组近地轨道的 r1、Δt，解算后用 propagate 回推，统计位置误差的中位数和最大发散点，误差中位数超过 1 米就要回头查迭代容差。

4.3 与成熟工具交叉验证：用 GMAT、SPICE 的习惯

写到自己实现的解算器与商用或科研工具交叉验证，是轨道设计里最容易被省略但最重要的习惯。常见做法是用 GMAT 或 JPL 的 SPICE 工具做一个同样的转移算例：给定相同的 r1、r2、Δt，运行它们的兰伯特解算器，对比 v1、v2 的差值和传播后的终端位置误差。两套独立实现的结果如果速度差超过毫米每秒量级，说明至少一方有 bug。GMAT 使用随遇配置的 DSN 或 Gooding 算法，SPICE 则不直接提供兰伯特解算器，需要用 lmb 库或 SPICE 工具箱拼装。值不值得做这件事？如果你在做的是真正的飞行任务，性能验证一步都省不得；如果只是课题预研，至少挑一个最危险的算例——大转移角加长程解、多圈解、近 180 度奇异点这三类——跑一遍交叉验证，能把隐藏的数值问题暴露大半。

5. 兰伯特转移避坑与排查：5 个让结果翻车的细节

5.1 现象：r1 和 r2 夹角接近 180 度时，解算器直接输出 NaN

这是兰伯特问题最经典的数值奇异点。当两个位置矢量几乎反向共线时，转移弧段对应的角度接近 π，3.3 节代码里 A 的表达式分母 sqrt(1 - cos_dnu) 趋向 0，分子 sin(dnu) 也趋向 0，理论上比值趋向定值，但实际计算中两个接近 0 的数相除，误差会放大到不可接受，迭代完全发散。原因本质上是几何退化——两个共线反向的点无法确定唯一的轨道平面，任何经过这两个点的平面都能构成一条弧段。解决思路分两层：如果夹角在 179~181 度之间，直接放弃当前转移方向，改用另一个方向（长程替代短程）；如果必须处理 180 度附近的场景，给 A 的计算加一个基于 r1、r2 叉积方向的约束，强行指定轨道平面，再用退化公式 A = 未知量的极限值。

5.2 现象：给了多圈解 N=2，结果收敛到 N=0 的短程解

多圈解的 z 值必须落在 (2Nπ)² 附近的区间，不同 N 的解在 z 轴上相互分隔。直接让牛顿迭代从 z=0 起步，大概率收敛到 N=0 的解，因为 F(z) 的零点对初值有捕获域。解决方法是把 z 的初值直接设为 z0 = (2Nπ)²，同时把迭代区间限制在 (2Nπ)² ± 2π 的范围内，超出就强制回退。调试时的一个技巧是把 F(z) 画出来看零点分布——你会看到 z 轴上多个零点交错，只有限制区间才能稳定命中目标。还有一个细节：N 取太大时，F(z) 的曲线会变得非常陡峭，牛顿法容易在相邻零点间跳跃，所以实际工程里 N 超过 3 圈时我会改用二分法而不是牛顿法。

5.3 现象：转移时间太短，迭代不收敛或速度值数量级离谱

当 Δt 趋近于 0，兰伯特问题的物理极限趋近于一条直线碰撞轨道——航天器沿 r1 到 r2 的直线落向中心天体，速度趋向无穷大。这不是数值问题，是物理边界。现象表现为两种：迭代过程中 z 一直往负方向跑（双曲线轨道），或者恢复出来的 v1 模长超过逃逸速度的几十倍。解决方法是先算当前解域的极限转移时间 Δt_min，在外部做一个判断：如果用户传入的 Δt 小于 Δt_min，直接报错而不是让解算器硬算。Δt_min 的近似估算可以这样：把转移轨道视作通过 r1、r2 的能量最小椭圆，用椭圆转移时间公式做一个下限估计，误差在 10% 以内足够用。

5.4 现象：迭代变量 z 在 ±50 之外反复震荡，程序卡死

z 跑到 ±50 之外基本宣告输入有问题或者单位混用。最常见的原因是 3.1 节的归一化没有做，直接把 km 和秒交给迭代器，z 的量级被原始数据的量级放大。另一种原因是最初我给过 cr3bp 场景里传递兰伯特参数给其他模块，用户把 μ 写成 398600.4418 但距离单位用的还是 km，归一化后的 μ=1 不成立，方程右端整体失效。解决办法只有一条：任何输入都要走统一入口做归一化，在解算器内部检查 r1_mag、r2_mag、dt_n 三个量是否落在合理范围内（0.5~2、0.5~2、0.1~20），超出就抛异常。血泪经验是把这个检查放在归一化函数里做，不要放在调用端，因为你永远不知道下游调用者会不会直接传原始值。

5.5 现象：算出来的速度在地球系下方向正常，但传播后终端位置偏了半个地球

这个问题多数不是解算器的问题，而是输入 r2 与 r1 不在同一个惯性系里。轨道设计里 r1 常用瞬时轨道要素转换得到，r2 从星历查询，星历可能用的 J2000 系，而轨道要素转位置矢量用的是 TEME 系，两个坐标系之间差一个岁差章动旋转，位置偏差几百公里。解决方法是统一坐标系：在调用兰伯特解算器之前，把所有输入转到同一个惯性系（推荐 J2000 或 EME2000），并且在输入结构体里坐上标的标签，后续任何模块只是用数据都带着 coords 字段。排查时先看 r1、r2 夹角是否接近轨道倾角差，再看两个位置矢量与速度和轨道面的垂直关系，基本上能锁定坐标系旋转量。

6. 工程落地：从单次解算到交会搜索的三个进阶技巧

6.1 技巧一：把兰伯特解算器封装成时间搜索函数

单次兰伯特求解只能解决给定时间的转移，交会任务设计真正需要的是一个搜索器：给定期望交会时刻的窗口，扫描多个 Δt 候选值，每个 Δt 算一次 lambert_solve，画出转移所需的 v1 模长曲线，曲线的谷底就是燃料最优转移方案。做法是外层循环每 30 秒一个步长扫一小时窗口，内层调用解算器时要让 direction 参数跟着夹角变化——夹角小于 90 度用短程，大于 90 度就要同时尝试两个方向取速度增量小者。这个搜索器写在代码里不到五十行，但能把兰伯特解算器的价值真正放大到工程决策层面。

6.2 技巧二：多圈解的参数约束先于数学求解

多圈解不是无条件的。给定转移时间 Δt 和圈数 N，轨道还要满足实际的能源约束——点火速度增量不能超过某个阈值。实际操作中我会先用一个粗估公式：Δv ≈ sqrt(μ/r1) 乘一个与转移角相关的系数，初步判断 Δv 是否在可接受范围，再调用多圈解算器。这个顺序很重要，因为在搜索循环里调多圈解算器的代价远高于零圈解，先把不可行区域过滤掉，能省下大量计算时间。

6.3 技巧三：把验证脚本变成解算器的一部分

自己的解算器每一版改动后，验证不能只靠一次手跑。我习惯把第 4.2 节的传播验证写成一个单元测试文件，每次改完代码就随机生成一百组算例，自动统计终端位置误差分布。误差超过 1 米就报警，超过 1 公里的会直接拒绝发布。这个习惯帮我挡掉过好几次回归 bug——最严重的一次是改了 Stumpff 函数的截断常数，本来只有 1e-12 的精度损失，在验证脚本里被放大到位置误差 20 公里，没有验证脚本根本发现不了。希望这个验证习惯同样帮到你，把你自己的兰伯特工具链搭稳固。

本文还有配套的精品资源，点击获取