簡(jiǎn)介這份資源聚焦航天工程中的蘭伯特轉(zhuǎn)移問(wèn)題面向天體力學(xué)、軌道設(shè)計(jì)與航天任務(wù)分析方向的學(xué)習(xí)者與工程師提供求解蘭伯特問(wèn)題的MATLAB實(shí)現(xiàn)思路。蘭伯特轉(zhuǎn)移以雙曲型軌道實(shí)現(xiàn)兩點(diǎn)間高效快速的軌道機(jī)動(dòng)其核心是在兩體問(wèn)題下確定初始速度、末端速度與轉(zhuǎn)移時(shí)間并區(qū)分順時(shí)針與逆時(shí)針兩種轉(zhuǎn)移情形廣泛用于近地軌道抬升、軌道面變更及地月、地火等星際轉(zhuǎn)移任務(wù)。壓縮包內(nèi)共1個(gè)文件為m格式的MATLAB腳本整體約2KB可直接用于輸入起止位置、轉(zhuǎn)移時(shí)間與航天器質(zhì)量等參數(shù)進(jìn)而計(jì)算升交點(diǎn)、降交點(diǎn)坐標(biāo)、飛行時(shí)間及所需總沖量幫助讀者理解數(shù)值與解析求解流程。目前已有2005人學(xué)習(xí)下載適合作為軌道轉(zhuǎn)移計(jì)算的入門參考與腳本模板便于在此基礎(chǔ)上開(kāi)展任務(wù)仿真與燃料優(yōu)化分析。1. 蘭伯特轉(zhuǎn)移到底在算什么從兩條軌道到一段飛行時(shí)間如果你手頭有兩組軌道根數(shù)或者兩個(gè)位置矢量再加上一個(gè)飛行時(shí)間想反推出中間那段轉(zhuǎn)移軌道長(zhǎng)什么樣、需要多大的速度增量那你碰上的就是蘭伯特問(wèn)題。它不關(guān)心你中途怎么飛只認(rèn)三個(gè)量起點(diǎn)、終點(diǎn)、時(shí)間。聽(tīng)起來(lái)簡(jiǎn)單但它幾乎是所有軌道轉(zhuǎn)移任務(wù)的總?cè)肟凇獜慕剀壍捞酵杰壍馈牡厍蛱右莸交鹦恰耐２窜壍狼腥氕h(huán)月軌道方案設(shè)計(jì)階段第一件事往往就是解一次蘭伯特。我最早接觸它是在做地月轉(zhuǎn)移窗口掃描的時(shí)候一開(kāi)始以為套個(gè)公式就行結(jié)果被多圈解、奇異區(qū)和收斂性折騰了好幾天。這篇筆記就按我實(shí)際做工程的順序來(lái)先把蘭伯特轉(zhuǎn)移的幾何和物理講清楚再落到可復(fù)現(xiàn)的求解流程、參數(shù)怎么設(shè)、代碼怎么寫最后把踩過(guò)的坑一條條擺出來(lái)。適合正在做軌道設(shè)計(jì)、任務(wù)分析、或者想自己寫一套轉(zhuǎn)移求解工具的從業(yè)者新手能照著跑通熟手能對(duì)著邊界條件摳細(xì)節(jié)。2. 蘭伯特轉(zhuǎn)移的幾何與時(shí)間方程為什么它是個(gè)邊值問(wèn)題2.1 從開(kāi)普勒軌道到蘭伯特定理蘭伯特定理說(shuō)的是一段開(kāi)普勒軌道上兩點(diǎn)之間的飛行時(shí)間只取決于這兩點(diǎn)的位置、軌道半長(zhǎng)軸以及兩點(diǎn)之間的弦長(zhǎng)跟軌道偏心率、近地點(diǎn)幅角這些形狀參數(shù)沒(méi)有直接關(guān)系。換句話說(shuō)只要給定起點(diǎn)位置矢量 r1、終點(diǎn)位置矢量 r2 和飛行時(shí)間 Δt轉(zhuǎn)移軌道的半長(zhǎng)軸就被唯一確定了在給定圈數(shù)下。這跟初值問(wèn)題正好相反。初值問(wèn)題是知道位置和速度往后積分蘭伯特是知道兩端位置和時(shí)間反推速度。所以它天然是個(gè)邊值問(wèn)題求解的核心就是找到一個(gè)半長(zhǎng)軸 a使得從 r1 沿軌道飛到 r2 恰好花 Δt。工程上我們真正要的是起點(diǎn)速度 v1 和終點(diǎn)速度 v2因?yàn)樗俣仍隽?Δv v1 - v_初始軌道、v2 - v_目標(biāo)軌道直接決定推進(jìn)劑預(yù)算。蘭伯特求解器輸出的就是這兩個(gè)速度矢量。2.2 轉(zhuǎn)移角 Δθ 與弦長(zhǎng) c 的幾何關(guān)系設(shè)起點(diǎn)位置矢量 r1、終點(diǎn) r2兩者夾角就是轉(zhuǎn)移角 Δθcos(Δθ) (r1 · r2) / (|r1| |r2|)弦長(zhǎng) c 由余弦定理給出c sqrt(|r1|^2 |r2|^2 - 2 |r1| |r2| cos(Δθ))這里有個(gè)必須注意的點(diǎn)Δθ 的取值不是 acos 直接給的那個(gè) [0, π]而是要根據(jù)飛行方向判斷。如果轉(zhuǎn)移是順行prograde且 Δθ 實(shí)際超過(guò) π就要取 2π - Δθ。判斷方法通常用 r1 × r2 的 z 分量符號(hào)結(jié)合任務(wù)規(guī)定的繞行方向。這一步搞錯(cuò)后面所有速度全錯(cuò)而且錯(cuò)得很隱蔽——因?yàn)楣秸諛邮諗恐皇墙獬鰜?lái)是另一條軌道。半周長(zhǎng) s 定義為s (|r1| |r2| c) / 2s 是后續(xù)時(shí)間方程里的關(guān)鍵中間量它把幾何信息壓縮成一個(gè)標(biāo)量。2.3 時(shí)間方程從拉格朗日形式到通用變量蘭伯特問(wèn)題的時(shí)間方程有幾種等價(jià)寫法我一般用拉格朗日形式的通用變量版本數(shù)值上比較穩(wěn)。核心是引入一個(gè)無(wú)量綱參數(shù) z Δθ 相關(guān)的變量或者用半長(zhǎng)軸 a 來(lái)表達(dá)。對(duì)橢圓軌道a 0時(shí)間方程可以寫成Δt sqrt(a^3 / μ) * [ (α - sin α) - (β - sin β) ]其中 α、β 由半周長(zhǎng)和半長(zhǎng)軸決定sin(α/2) sqrt(s / (2a)) sin(β/2) sqrt((s - c) / (2a))μ 是中心天體引力常數(shù)地球取 398600.4418 km3/s2月球取 4902.8 km3/s2這些值必須用對(duì)差一點(diǎn)在長(zhǎng)轉(zhuǎn)移時(shí)間里會(huì)放大成幾十公里的位置誤差。對(duì)雙曲軌道a 0用雙曲正弦形式Δt sqrt((-a)^3 / μ) * [ (sinh α - α) - (sinh β - β) ]拋物線情況a → ∞是奇異點(diǎn)實(shí)際工程里很少正好落在拋物線上但數(shù)值求解時(shí)如果迭代到 a 很大要小心溢出。常見(jiàn)做法是設(shè)一個(gè) a 的上限超過(guò)就按雙曲處理或直接報(bào)錯(cuò)。2.4 多圈解為什么同一個(gè) Δt 可能對(duì)應(yīng)多條軌道這是蘭伯特問(wèn)題最容易被忽略的地方。給定 r1、r2、Δt解可能不止一個(gè)。因?yàn)檗D(zhuǎn)移軌道可以繞中心天體轉(zhuǎn) 0 圈、1 圈、2 圈……每多轉(zhuǎn)一圈飛行時(shí)間就多一個(gè)軌道周期但起點(diǎn)終點(diǎn)位置不變。所以對(duì)同一個(gè) Δt可能存在多個(gè)半長(zhǎng)軸對(duì)應(yīng)不同的圈數(shù) N。工程上默認(rèn)取 N 0也就是最短的那條轉(zhuǎn)移角小于 2π 且不繞整圈。但在某些任務(wù)里比如長(zhǎng)時(shí)間滑行的轉(zhuǎn)移N 1 甚至 N 2 的解反而更省燃料。我一般會(huì)在求解器里把 N 作為輸入?yún)?shù)掃描 N 0, 1, 2把每個(gè)解的 Δv 都算出來(lái)對(duì)比。多圈解的存在性有前提Δt 必須大于該圈數(shù)對(duì)應(yīng)的最小時(shí)間。如果 Δt 太小N 1 無(wú)解求解器會(huì)不收斂。這時(shí)候不要硬迭代直接判斷并返回?zé)o解。3. 用 Python 實(shí)現(xiàn)蘭伯特求解從幾何輸入到速度輸出3.1 最小可運(yùn)行代碼牛頓迭代求半長(zhǎng)軸下面這段是我常用的核心求解器輸入 r1、r2、Δt、μ 和圈數(shù) N輸出 v1、v2。用的是牛頓法迭代半長(zhǎng)軸 a配合通用變量時(shí)間方程。import numpy as np def lambert_solver(r1, r2, dt, mu, N0, progradeTrue, tol1e-8, max_iter100): 蘭伯特轉(zhuǎn)移求解器 r1, r2: 起點(diǎn)/終點(diǎn)位置矢量 (km) dt: 飛行時(shí)間 (s) mu: 引力常數(shù) (km^3/s^2) N: 圈數(shù), 0 表示不繞整圈 prograde: 是否順行 返回: v1, v2 (km/s) r1 np.asarray(r1, dtypefloat) r2 np.asarray(r2, dtypefloat) r1_norm np.linalg.norm(r1) r2_norm np.linalg.norm(r2) # 轉(zhuǎn)移角 cos_dtheta np.dot(r1, r2) / (r1_norm * r2_norm) cos_dtheta np.clip(cos_dtheta, -1.0, 1.0) dtheta np.arccos(cos_dtheta) # 根據(jù)順行/逆行和叉乘方向修正轉(zhuǎn)移角 cross_z np.cross(r1, r2)[2] if prograde: if cross_z 0: dtheta 2 * np.pi - dtheta else: if cross_z 0: dtheta 2 * np.pi - dtheta # 弦長(zhǎng)和半周長(zhǎng) c np.sqrt(r1_norm**2 r2_norm**2 - 2 * r1_norm * r2_norm * cos_dtheta) s (r1_norm r2_norm c) / 2.0 # 初始猜測(cè)半長(zhǎng)軸 a s / 2.0 def time_of_flight(a): if a 0: alpha 2 * np.arcsin(np.sqrt(s / (2 * a))) beta 2 * np.arcsin(np.sqrt((s - c) / (2 * a))) if N 0: return np.sqrt(a**3 / mu) * ((alpha - np.sin(alpha)) - (beta - np.sin(beta))) else: return np.sqrt(a**3 / mu) * ((alpha - np.sin(alpha)) - (beta - np.sin(beta)) 2 * np.pi * N) else: alpha 2 * np.arcsinh(np.sqrt(s / (-2 * a))) beta 2 * np.arcsinh(np.sqrt((s - c) / (-2 * a))) return np.sqrt((-a)**3 / mu) * ((np.sinh(alpha) - alpha) - (np.sinh(beta) - beta)) # 牛頓迭代 for _ in range(max_iter): f time_of_flight(a) - dt da a * 1e-6 df (time_of_flight(a da) - time_of_flight(a - da)) / (2 * da) if abs(df) 1e-14: break a_new a - f / df if abs(a_new - a) tol: a a_new break a a_new # 由 a 反算 f 和 g 函數(shù), 再求速度 f 1 - (r2_norm / (np.sqrt(mu) * np.sqrt(a))) * np.sin( 2 * np.arcsin(np.sqrt(s / (2 * a))) - 2 * np.arcsin(np.sqrt(s / (2 * a))) ) if a 0 else None # 更穩(wěn)妥的做法: 用拉格朗日系數(shù)直接算 # 這里用標(biāo)準(zhǔn) f/g 表達(dá)式 if a 0: alpha 2 * np.arcsin(np.sqrt(s / (2 * a))) beta 2 * np.arcsin(np.sqrt((s - c) / (2 * a))) A np.sqrt(mu / (4 * a)) * (alpha - np.sin(alpha) - (beta - np.sin(beta))) else: alpha 2 * np.arcsinh(np.sqrt(s / (-2 * a))) beta 2 * np.arcsinh(np.sqrt((s - c) / (-2 * a))) A np.sqrt(mu / (-4 * a)) * (np.sinh(alpha) - alpha - (np.sinh(beta) - beta)) # 用 f/g 函數(shù)求 v1, v2 f_coef 1 - (r2_norm / (np.sqrt(mu) * np.sqrt(a))) * np.sin( (alpha - beta) / 2 ) if a 0 else 1 - (r2_norm / (np.sqrt(mu) * np.sqrt(-a))) * np.sinh( (alpha - beta) / 2 ) g_coef (r1_norm * r2_norm / np.sqrt(mu * a)) * np.sin( (alpha - beta) / 2 ) if a 0 else (r1_norm * r2_norm / np.sqrt(mu * (-a))) * np.sinh( (alpha - beta) / 2 ) v1 (r2 - f_coef * r1) / g_coef v2 (g_coef * r2 - r1) / g_coef # 注意: 這里需要 g_dot, 簡(jiǎn)化寫法 return v1, v2上面這段代碼里牛頓迭代部分是對(duì)的但 f/g 反算速度那段我故意留了個(gè)不完整的寫法因?yàn)閷?shí)際工程里更推薦用通用變量直接算 f、g、g_dot避免符號(hào)錯(cuò)誤。下面給一個(gè)更干凈的版本只算 v1 和 v2def lambert_velocity(r1, r2, dt, mu, N0, progradeTrue): r1 np.asarray(r1, dtypefloat) r2 np.asarray(r2, dtypefloat) r1n np.linalg.norm(r1) r2n np.linalg.norm(r2) cos_dtheta np.clip(np.dot(r1, r2) / (r1n * r2n), -1.0, 1.0) dtheta np.arccos(cos_dtheta) cross_z np.cross(r1, r2)[2] if prograde and cross_z 0: dtheta 2 * np.pi - dtheta if not prograde and cross_z 0: dtheta 2 * np.pi - dtheta c np.sqrt(r1n**2 r2n**2 - 2 * r1n * r2n * cos_dtheta) s (r1n r2n c) / 2.0 # 用二分法求 a, 比牛頓更穩(wěn) a_min s / 2.0 * 0.5 a_max s / 2.0 * 100.0 for _ in range(200): a 0.5 * (a_min a_max) if a 0: alpha 2 * np.arcsin(np.sqrt(s / (2 * a))) beta 2 * np.arcsin(np.sqrt((s - c) / (2 * a))) tof np.sqrt(a**3 / mu) * ((alpha - np.sin(alpha)) - (beta - np.sin(beta)) 2 * np.pi * N) else: alpha 2 * np.arcsinh(np.sqrt(s / (-2 * a))) beta 2 * np.arcsinh(np.sqrt((s - c) / (-2 * a))) tof np.sqrt((-a)**3 / mu) * ((np.sinh(alpha) - alpha) - (np.sinh(beta) - beta)) if tof dt: a_min a else: a_max a # 用 f/g 函數(shù) if a 0: alpha 2 * np.arcsin(np.sqrt(s / (2 * a))) beta 2 * np.arcsin(np.sqrt((s - c) / (2 * a))) f 1 - (a / r1n) * (1 - np.cos(alpha - beta)) g dt - np.sqrt(a**3 / mu) * ((alpha - beta) - (np.sin(alpha) - np.sin(beta))) g_dot 1 - (a / r2n) * (1 - np.cos(alpha - beta)) else: alpha 2 * np.arcsinh(np.sqrt(s / (-2 * a))) beta 2 * np.arcsinh(np.sqrt((s - c) / (-2 * a))) f 1 - ((-a) / r1n) * (1 - np.cosh(alpha - beta)) g dt - np.sqrt((-a)**3 / mu) * ((np.sinh(alpha) - np.sinh(beta)) - (alpha - beta)) g_dot 1 - ((-a) / r2n) * (1 - np.cosh(alpha - beta)) v1 (r2 - f * r1) / g v2 (g_dot * r2 - r1) / g return v1, v2這段代碼的邏輯說(shuō)明先用二分法把半長(zhǎng)軸 a 夾逼出來(lái)因?yàn)闀r(shí)間方程對(duì) a 是單調(diào)的在給定 N 下二分比牛頓更不容易發(fā)散。然后利用拉格朗日系數(shù) f、g、g_dot 直接由位置求速度避免顯式算 f_dot 帶來(lái)的符號(hào)混亂。參數(shù)說(shuō)明r1、r2 單位 kmdt 單位秒mu 單位 km3/s2。N 默認(rèn) 0prograde 默認(rèn) True。二分區(qū)間我取的是 [s/4, 50s]覆蓋了絕大多數(shù)近地和深空轉(zhuǎn)移。如果 dt 特別大比如幾個(gè)月的地火轉(zhuǎn)移a_max 要放大到 100s 以上否則會(huì)夾不到解。3.2 參數(shù)怎么設(shè)μ、圈數(shù)、順行逆行μ 的取值直接決定速度量級(jí)。地球 398600.4418月球 4902.8火星 42828.3太陽(yáng) 1.32712440018e11。這些值我一般寫成常量字典避免每次手敲。圈數(shù) N 的選擇近地軌道轉(zhuǎn)移通常 N 0。地月轉(zhuǎn)移 N 0 或 1 都可能取決于飛行時(shí)間。如果 Δt 超過(guò)一個(gè)軌道周期N 1 的解可能更省 Δv。我一般會(huì)掃 N 0, 1, 2把每個(gè)解的 Δv 列出來(lái)對(duì)比。順行逆行從地球出發(fā)去火星順行是常規(guī)選擇。但如果 r1 × r2 的 z 分量為負(fù)而任務(wù)要求順行就必須把 Δθ 修正到 2π - Δθ。這個(gè)判斷錯(cuò)了解出來(lái)的軌道會(huì)繞到另一側(cè)Δv 可能差好幾 km/s。3.3 驗(yàn)證解的正確性用二體積分回代解出 v1、v2 之后不要直接信。我一般會(huì)做一步回代驗(yàn)證用 r1、v1 作為初值用二體問(wèn)題積分到 Δt看終點(diǎn)位置跟 r2 差多少。如果差在幾米到幾十米量級(jí)說(shuō)明解是對(duì)的如果差了幾百公里說(shuō)明轉(zhuǎn)移角或圈數(shù)搞錯(cuò)了。from scipy.integrate import solve_ivp def propagate_two_body(r0, v0, dt, mu): def rhs(t, y): r y[:3] v y[3:] r_norm np.linalg.norm(r) a -mu * r / r_norm**3 return np.concatenate([v, a]) y0 np.concatenate([r0, v0]) sol solve_ivp(rhs, [0, dt], y0, rtol1e-10, atol1e-10) return sol.y[:3, -1], sol.y[3:, -1] # 驗(yàn)證 r1 np.array([7000.0, 0.0, 0.0]) r2 np.array([0.0, 8000.0, 0.0]) dt 3600.0 mu 398600.4418 v1, v2 lambert_velocity(r1, r2, dt, mu) r_check, v_check propagate_two_body(r1, v1, dt, mu) print(位置誤差 (km):, np.linalg.norm(r_check - r2))如果位置誤差在 1e-3 km 以內(nèi)基本可以放心用。這個(gè)回代步驟我強(qiáng)烈建議每次都做尤其是改了轉(zhuǎn)移角判斷邏輯之后。4. 蘭伯特轉(zhuǎn)移的避坑與排查那些讓 Δv 悄悄翻倍的細(xì)節(jié)4.1 轉(zhuǎn)移角判斷反了解出來(lái)是另一條軌道現(xiàn)象求解器收斂速度也正常但 Δv 比預(yù)期大很多或者軌道形狀明顯不對(duì)。原因Δθ 用了 acos 的默認(rèn)值 [0, π]沒(méi)有根據(jù)順行/逆行和叉乘方向修正。當(dāng)實(shí)際轉(zhuǎn)移角超過(guò) π 時(shí)解出來(lái)的是補(bǔ)角對(duì)應(yīng)的短程軌道方向完全反了。解決在算完 acos 之后強(qiáng)制判斷 cross_z 符號(hào)。順行且 cross_z 0 時(shí)取 2π - Δθ逆行且 cross_z 0 時(shí)取 2π - Δθ。這個(gè)邏輯我封裝成獨(dú)立函數(shù)每次調(diào)用前先確認(rèn)。4.2 多圈解漏掃錯(cuò)過(guò)更省燃料的窗口現(xiàn)象N 0 的解 Δv 很大任務(wù)看起來(lái)不可行但換一個(gè)飛行時(shí)間就突然可行了。原因只算了 N 0沒(méi)有掃 N 1、2。長(zhǎng)時(shí)間轉(zhuǎn)移里多繞一圈可能讓半長(zhǎng)軸更接近目標(biāo)軌道Δv 反而更小。解決把 N 作為循環(huán)變量對(duì)每個(gè) N 求解并記錄 Δv。如果某個(gè) N 無(wú)解Δt 小于該圈數(shù)最小時(shí)間直接跳過(guò)不要硬迭代。我一般會(huì)輸出一張表N、a、Δv1、Δv2、總 Δv人工挑最優(yōu)。4.3 二分區(qū)間設(shè)太窄深空轉(zhuǎn)移夾不到解現(xiàn)象二分法跑完 200 次a 停在邊界上回代誤差巨大。原因a_max 設(shè)成了 50s但地火轉(zhuǎn)移的 a 可能到幾個(gè) AU遠(yuǎn)超這個(gè)范圍。解決根據(jù)任務(wù)類型動(dòng)態(tài)設(shè) a_max。近地轉(zhuǎn)移 50s 夠用地月轉(zhuǎn)移設(shè)到 200s行星際轉(zhuǎn)移直接設(shè)到 1e4 s 量級(jí)。或者用自適應(yīng)擴(kuò)展先試一個(gè)區(qū)間如果解落在邊界就把區(qū)間翻倍再試。4.4 雙曲分支的 sinh 溢出現(xiàn)象迭代過(guò)程中報(bào) overflow或者 a 變成 NaN。原因a 接近 0 時(shí)sqrt(s / (-2a)) 變得很大sinh 直接溢出。解決在 a 0 的分支里加保護(hù)如果 sqrt(s / (-2a)) 50就認(rèn)為 a 太小直接返回?zé)o解或把 a 限制在一個(gè)下限。實(shí)際工程里 a 不會(huì)真的趨近 0因?yàn)槟菍?duì)應(yīng)拋物線能量無(wú)窮大。4.5 μ 用錯(cuò)速度整體偏移現(xiàn)象回代位置誤差不大但 Δv 跟別人對(duì)不上差一個(gè)固定比例。原因μ 用了 398600 而不是 398600.4418或者月球用了地球的 μ。解決把 μ 寫成常量字典調(diào)用時(shí)顯式傳參不要用全局變量。每次換中心天體先檢查 μ 值。5. 進(jìn)階技巧用 porkchop 圖快速鎖定發(fā)射窗口5.1 掃描出發(fā)和到達(dá)日期的 Δv 網(wǎng)格蘭伯特求解器最實(shí)用的進(jìn)階用法是畫 porkchop 圖。做法很簡(jiǎn)單固定起點(diǎn)軌道和終點(diǎn)軌道掃描出發(fā)日期 t1 和到達(dá)日期 t2對(duì)每個(gè) (t1, t2) 組合算一次蘭伯特轉(zhuǎn)移記錄總 Δv。把 Δv 畫成等高線圖低 Δv 的區(qū)域就是發(fā)射窗口。import numpy as np import matplotlib.pyplot as plt def porkchop(r1_func, r2_func, t1_range, t2_range, mu): dv_grid np.zeros((len(t1_range), len(t2_range))) for i, t1 in enumerate(t1_range): r1 r1_func(t1) for j, t2 in enumerate(t2_range): if t2 t1: dv_grid[i, j] np.nan continue r2 r2_func(t2) dt (t2 - t1) * 86400.0 try: v1, v2 lambert_velocity(r1, r2, dt, mu) dv1 np.linalg.norm(v1 - v1_initial(r1)) dv2 np.linalg.norm(v2 - v2_target(r2)) dv_grid[i, j] dv1 dv2 except Exception: dv_grid[i, j] np.nan return dv_grid這段代碼里 r1_func 和 r2_func 是起點(diǎn)和終點(diǎn)軌道在給定時(shí)刻的位置函數(shù)v1_initial 和 v2_target 是對(duì)應(yīng)軌道的速度。實(shí)際用時(shí)r1_func 可以用二體解析解或者數(shù)值積分得到。參數(shù)說(shuō)明t1_range 和 t2_range 單位是天dt 轉(zhuǎn)成秒。dv_grid 里 NaN 表示無(wú)解或 t2 t1。畫圖時(shí)用 contourf把 Δv 低于某個(gè)閾值的區(qū)域標(biāo)出來(lái)就是可行窗口。5.2 從 porkchop 圖讀窗口寬度和 Δv 裕度porkchop 圖上的低 Δv 區(qū)域通常是個(gè)斜橢圓長(zhǎng)軸方向?qū)?yīng)出發(fā)和到達(dá)日期的耦合關(guān)系。窗口寬度看的是這個(gè)橢圓在 t1 軸上的投影。如果投影只有幾天說(shuō)明窗口很窄發(fā)射機(jī)會(huì)稍縱即逝如果有幾周說(shuō)明容錯(cuò)空間大。我一般會(huì)在圖上疊加一條等 Δv 線比如 3.5 km/s然后看這條線包住的區(qū)域有多大。實(shí)際任務(wù)里還要留 5% 到 10% 的 Δv 裕度所以真正可用的窗口比圖上看到的還要窄一圈。5.3 用網(wǎng)格搜索代替手工調(diào)參早期我調(diào)蘭伯特參數(shù)是手工試改一個(gè)數(shù)跑一次效率極低。后來(lái)改成網(wǎng)格搜索把 N、prograde、a_max 這些參數(shù)做成組合批量跑自動(dòng)挑 Δv 最小的。這樣不僅快還能發(fā)現(xiàn)一些反直覺(jué)的解比如逆行軌道在某些窗口下反而更省。一個(gè)具體的習(xí)慣每次做新任務(wù)先跑一張粗網(wǎng)格 porkchop步長(zhǎng) 1 天看大趨勢(shì)再在低 Δv 區(qū)域跑細(xì)網(wǎng)格步長(zhǎng) 0.1 天精確定位。粗網(wǎng)格用 N 0細(xì)網(wǎng)格再掃 N 1、2。這樣既不會(huì)漏掉多圈解也不會(huì)在無(wú)解區(qū)域浪費(fèi)時(shí)間。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取