简介兰伯特转移用于求解航天器在两固定位置间的最优转移轨道是轨道设计与任务规划中的基础算法。该压缩包内含1个MATLAB脚本lambert.m面向天体力学与航天轨道计算学习者以及需要快速验算转移参数的工程师。脚本通过输入起始/目标位置、转移时间等参数解算顺/逆时针转移的初速与终速并给出飞行时间、升交点/降交点坐标及总冲量等结果可辅助评估任务可行性并降低燃料消耗。资源包以单个m文件为主整包约2KB无需安装依赖直接在MATLAB中运行即可适合作为理解兰伯特问题的入门工具或嵌入个人计算流程。目前已有1986人浏览/学习脚本轻量、结构清晰便于对照算法公式梳理输入输出逻辑并在此基础上扩展多圈转移或摄动修正可显著缩短轨道转移算法的上手时间。 搞轨道计算的人几乎都会遇到兰伯特问题。简单说你手里有个飞行器在 t1 时刻位于 r1我希望它在 t2 时刻到达另一个位置 r2那么它应该飞一条什么样的轨道这就是兰伯特问题Lamberts 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 χ² / aS(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)) / zz 0S(z) (√z - sin(√z)) / (√z)³z 0当 z 0 时用双曲函数替换z 接近 0 时取级数展开的前几项。这两个函数的作用是把椭圆和双曲线的情况统一成一套表达式。第二计算辅助量 AA dm * sqrt(r1 * r2 * (1 cos(Δν)))这里 dm 控制短程或长程A 的符号决定了转移方向。第三计算中间变量 yy 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, dm1): 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, xtol1e-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 转成 kmAU 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, dm1) 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,), rtol1e-9) diff sol.y[:3, -1] - r2 print(终点位置误差 , diff)如果误差在公里量级甚至更小说明求解器和积分器没问题。如果误差大到离谱多半是兰伯特求解代码里 z 的区间、符号或者单位出了问题。这一步是检验算法正确性的黄金标准强烈建议每个人在做完兰伯特计算后都跑一遍。5. 写给新手的避坑指南5.1 转移角一定要先明确“走哪边”我见过不少初次接触兰伯特问题的人包括我自己早期都会在转移角上翻车。arccos 只会给出 0 ~ 180 度的小角但实际任务中你可能需要走大于 180 度的大弧。如果不先画图确认短程还是长程直接用默认的 dm1算出来的轨道很可能把航天器送到目标相反方向。正确做法拿到 r1 和 r2 后先画一张几何示意图标出中心天体位置、两个位置矢量、以及你计划沿着哪个方向飞行。确认转移角小于 180 度用 dm1大于 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 怎么走”的示意图再回来调参很多坑会好踩很多。本文还有配套的精品资源点击获取 兰伯特问题详解：从原理到普适变量法求解轨道转移
简介兰伯特转移用于求解航天器在两固定位置间的最优转移轨道是轨道设计与任务规划中的基础算法。该压缩包内含1个MATLAB脚本lambert.m面向天体力学与航天轨道计算学习者以及需要快速验算转移参数的工程师。脚本通过输入起始/目标位置、转移时间等参数解算顺/逆时针转移的初速与终速并给出飞行时间、升交点/降交点坐标及总冲量等结果可辅助评估任务可行性并降低燃料消耗。资源包以单个m文件为主整包约2KB无需安装依赖直接在MATLAB中运行即可适合作为理解兰伯特问题的入门工具或嵌入个人计算流程。目前已有1986人浏览/学习脚本轻量、结构清晰便于对照算法公式梳理输入输出逻辑并在此基础上扩展多圈转移或摄动修正可显著缩短轨道转移算法的上手时间。 搞轨道计算的人几乎都会遇到兰伯特问题。简单说你手里有个飞行器在 t1 时刻位于 r1我希望它在 t2 时刻到达另一个位置 r2那么它应该飞一条什么样的轨道这就是兰伯特问题Lamberts 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 χ² / aS(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)) / zz 0S(z) (√z - sin(√z)) / (√z)³z 0当 z 0 时用双曲函数替换z 接近 0 时取级数展开的前几项。这两个函数的作用是把椭圆和双曲线的情况统一成一套表达式。第二计算辅助量 AA dm * sqrt(r1 * r2 * (1 cos(Δν)))这里 dm 控制短程或长程A 的符号决定了转移方向。第三计算中间变量 yy 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, dm1): 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, xtol1e-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 转成 kmAU 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, dm1) 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,), rtol1e-9) diff sol.y[:3, -1] - r2 print(终点位置误差 , diff)如果误差在公里量级甚至更小说明求解器和积分器没问题。如果误差大到离谱多半是兰伯特求解代码里 z 的区间、符号或者单位出了问题。这一步是检验算法正确性的黄金标准强烈建议每个人在做完兰伯特计算后都跑一遍。5. 写给新手的避坑指南5.1 转移角一定要先明确“走哪边”我见过不少初次接触兰伯特问题的人包括我自己早期都会在转移角上翻车。arccos 只会给出 0 ~ 180 度的小角但实际任务中你可能需要走大于 180 度的大弧。如果不先画图确认短程还是长程直接用默认的 dm1算出来的轨道很可能把航天器送到目标相反方向。正确做法拿到 r1 和 r2 后先画一张几何示意图标出中心天体位置、两个位置矢量、以及你计划沿着哪个方向飞行。确认转移角小于 180 度用 dm1大于 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 怎么走”的示意图再回来调参很多坑会好踩很多。本文还有配套的精品资源点击获取