★ FEATURED ARTICLE 简介资源包围绕MATLAB实现Halo轨道设计展开面向航天动力学、轨道设计方向的科研人员、工程技术人员及高年级本科生/研究生也可作为计算轨道力学课程的实践案例。内容系统涵盖三体问题动力学建模、Richardson三阶展开法生成初始轨道猜测、微分修正法优化轨道周期性与稳定性并通过三维轨迹绘制和相平面分析展示地月L2点Halo轨道的空间形态与动态特性。包内共15个文件容量约5.55MB以MATLAB脚本.m和Python脚本.py为主配合PNG结果图、CSV轨道数据、Markdown说明、文本报告以及依赖清单文件结构清晰既便于在MATLAB中复现主流程也可用Python快速交叉验证。目前已有105人学习/下载。配套代码覆盖参数定义、解析初值计算、轨道优化、可视化和数值对比全流程还给出轨道振幅、周期、能量耗散与逃逸速度等关键参数的分析结果通过解析解与数值解的对比验证了模型和算法的有效性读者可直接运行和修改。1. Halo轨道这个代码包到底能帮你解决什么Halo轨道设计在工程里是个典型的看起来有标准答案实际每一步都在踩坑的问题。标题里这套MATLAB项目核心不是画一条好看的螺旋线而是解决如何在圆限制性三体模型下用数值方法找到绕拉格朗日点L2稳定的三维周期轨道这件事。你可能是做深空探测器任务设计、卫星编队或者课程设计拿到代码包后最想知道的是它凭什么能收敛参数怎么改为什么我改了一个数就发散。下面我会按实际调通这类代码的顺序把建模、打靶、参数和坑位讲清楚保证你照着装完就能跑出第一条Halo轨道。2. 圆限制性三体问题与拉格朗日点Halo轨道建模前必须弄懂的底层逻辑2.1 圆限制性三体模型为什么它够用又方便传统上设计Halo轨道最常用的模型是圆限制性三体问题CR3BP。它把两个大天体比如地球和月球假设为围绕共同质心做匀速圆周运动第三个物体质量忽略不计。这个模型忽略太阳摄动、非球形引力和太阳光压但它保留了拉格朗日点附近的非线性动力学结构用于初步设计完全够用。我见过不少同行一上来就直接用全引力模型结果Halo轨道怎么都收敛不了最后发现问题不在算法而在模型太复杂。CR3BP的好处是方程是自治的、无量纲的而且只用一个参数mu就能描述系统方便对比和验证。你先在这个模型里把初值调出来再考虑加摄动这个顺序千万别反。真实工程里后续加太阳摄动等微扰是另一层工作如果基础轨道都没算对加再多修正项都是白搭。无量纲化是CR3BP最容易卡住的地方。长度单位取两个主天体之间的距离质量单位取两个主天体总质量时间单位取系统角速度的倒数。这样运动方程里只剩下一个参数mu也就是小质量天体占系统总质量的比值。地月系统mu约为0.01215日地系统mu约为3e-6这个数字会出现在所有代码参数里。我第一次用的时候就因为忘了把距离和时间归一化画出来的轨道长度和实际差了两个量级还以为程序写错了。2.2 拉格朗日点求根从无量纲方程到L1/L2坐标CR3BP的拉格朗日点是三力平衡点。对于L2点它位于两个主天体连线上、远离较大天体的外侧。求解拉格朗日点位置本质上是在x轴上解一个零点方程。用MATLAB的fzero非常方便前提是给一个靠近解的初值。对于地月系统L2大约在月球外侧1.1557无量纲距离处所以初值给1.1没问题。以下代码求地月L2坐标mu 0.012150584672391; xL2 fzero((x) x - (1-mu)*(xmu)/abs(xmu)^3 - mu*(x-1mu)/abs(x-1mu)^3, 1.1);这个方程来自x方向平衡条件引力、离心力和科氏力在平衡时抵消。注意abs的使用避免了在对负数取三次方时出现复数。如果你是日地系统mu约3e-6L2的x坐标约1.01fzero的初值要相应改到1.01附近否则可能收敛到L1甚至L3。我一般会先画出这个函数曲线看一眼再选初值这一步能避免很多玄学问题。2.3 线性化近似给微分改正准备一套像样的初值Halo轨道没有解析解只能靠数值迭代。但任何非线性迭代都需要一个靠谱的初值否则会发散。线性化模型能给出第一版猜测在L2点附近轨道在x-z平面内做小幅振荡在y方向存在中心流形。对于小幅振幅Halo初始z向速度vz取0初始y向速度vy0则需要在迭代中修正。常见的初值做法是固定一个x方向振幅Ax令初始状态为 [xL2 - Ax, 0, z0, 0, vy0, 0]。其中z0和vy0是待猜的量。黑匣子就在这里不同振幅下z0和vy0的关系只能靠迭代自己去逼近。我一般先给z0 Ax * 0.5vy0 0.2 * z0 或 -0.2 * z0然后交给微分改正。符号取决于你想要的轨道是绕L2上方还是下方这个在后文会讲。线性化理论还能帮你判断Halo轨道是否存在。对于经典Halo第I类x方向振幅必须大于某个临界值否则真正的三维周期轨道不存在只会收敛到平面Lissajous轨道。地月L2附近这个临界值大约在无量纲0.01量级。如果你给的Ax只有0.001微分改正怎么迭代都得不到vx0的闭合条件所以别在小振幅上浪费时间。这也是为什么很多初学脚本在给定小振幅时翻车。注意线性化只负责提供初值真正的轨道频率和振幅是由非线性迭代决定的。这也是微分改正算法成为Halo轨道设计核心的原因。下一章我会把它拆开讲。3. 微分改正法怎么在MATLAB里落地打靶、约束与迭代3.1 穿越xz平面的约束条件把周期轨道变成两个方程Halo轨道虽然复杂但有一个几何特点它关于xz平面对称。轨道在穿越xz平面时速度只有y方向分量不为零x和z方向分量为0。基于这个性质可以把寻找周期轨道转换为一个打靶问题。具体来说从一个位于xz平面的初始点 [x0, 0, z0, 0, vy0, 0] 出发积分等到轨道再次穿越xz平面时如果终点的速度 vx0 且 vz0那么这个轨道就是周期轨道。这是一个两点边值问题。由于轨道对称这半圈轨道映射到另一半就形成完整的闭合三维环线。这样我们有两个未知数z0和vy0和两个约束方程终点vx0vz0用牛顿法迭代就可以求解。这里的打靶体现在初值像子弹一样射出去看它在目标平面的落点偏差然后调整枪口。我第一次写的时候总觉得这个target是终点状态但实际上是两个速度分量约束想通了才不会写错。3.2 打靶法的迭代框架约束函数与雅可比矩阵迭代的第一步是定义约束函数。我会用一个独立的函数来算给定初始猜测下的终点速度约束。代码框架如下function g halo_constraints(guess, mu, xL2, Ax) % guess [z0; vy0] z0 guess(1); vy0 guess(2); state0 [xL2 - Ax; 0; z0; 0; vy0; 0]; opts odeset(Events, (t,s) event_xz(t,s), RelTol, 1e-12, AbsTol, 1e-12); [t, state, te, ye, ie] ode45((t,s) cr3bp_rhs(t,s,mu), [0 10], state0, opts); if isempty(ie) error(未穿越xz平面请检查初值或积分时长); end vx_end state(end,4); vz_end state(end,6); g [vx_end; vz_end]; end事件函数负责让积分停在y0的穿越瞬间而不是某个拍脑袋的固定时间。如果积分时间不精确终点并不是真实的穿越平面约束判断就没意义了。这里RelTol和AbsTol都设为1e-12是为了保证穿越时刻的状态精度。积分上限10是无量纲时间足够覆盖半个Halo周期地月L2大约是1.3个无量纲时间单位上限只是保护实际由事件终止。得到约束函数后我用有限差分计算雅可比矩阵。对这个2x2问题只需要扰动z0和vy0各一次再调用两次约束函数加上原始一次一共五次积分。然后更新猜测J finite_diff_jacobian(guess, mu, xL2, Ax); delta - J \ g; guess guess delta;与解析状态转移矩阵相比有限差分省去了推导和实现变分方程的麻烦代价是每次迭代多算几次积分。对于初学或验证算法这是最稳的路径。我最早写程序时也尝试过直接推导状态转移矩阵翻车几次后还是先用差分把流程跑通再用解析方法加速。3.3 初值生成与ode45积分别一上来就盲调ode45是比较适合这个问题的积分器因为它有自适应步长。但需要设置足够高的精度否则对轨道周期属性影响很大。建议至少RelTol1e-12AbsTol1e-12。不要用默认的1e-3那会让收敛解看起来闭合但稍微拉长时间就飞了。另一个常见坑是半个周期的位置选取。初始点在xz平面上但并不是所有穿越都能用。你需要让事件函数避免初始点本身触发。我用direction1只捕获y从负到正的穿越而初始点y0不满足这个方向所以不会在t0触发。但如果你选的初值让轨道先向上运动这个方向可能就不合适了。具体是北Halo还是南Halo会导致穿越方向不同。如果事件没有触发把direction改成-1或0试试。调这个迭代时如果g的范数在5次迭代后还不低于1e-8回头看初值猜测。一个常用技巧是先固定z0只迭代vy0等vx满足后再迭代z0。但2x2联合迭代通常也能直接工作前提是初值别太离谱。我曾经因为vy0符号反了让迭代在z00附近来回震荡最后靠打印中间状态才发现是符号问题。4. MATLAB代码实现从运动方程到Halo轨道绘制4.1 主脚本参数定义、求拉格朗日点、调用打靶下面是完整可运行的核心流程。我用的是函数式写法方便你改参数。主脚本负责设定参数、求拉格朗日点、迭代初值% main_halo_design.m mu 0.012150584672391; xL2 fzero((x) x - (1-mu)*(xmu)/abs(xmu)^3 - mu*(x-1mu)/abs(x-1mu)^3, 1.1); Ax 0.02; % x方向振幅 z0_guess 0.03; % z方向初值猜测 vy0_guess -0.04; % y方向初速度猜测 guess [z0_guess; vy0_guess]; for it 1:15 g halo_constraints(guess, mu, xL2, Ax); J finite_diff_jacobian(guess, mu, xL2, Ax); delta - J \ g; guess guess delta; fprintf(iter %d: |g|%.6e, z0%.10f, vy0%.10f\n, it, norm(g), guess(1), guess(2)); if norm(g) 1e-10, break; end end主脚本里Ax决定了Halo轨道在x方向的大小通常根据任务需求给定。z0和vy0的初值需要试凑不同Halo族方向不同。迭代循环里用了有限差分雅可比没有直接使用状态转移矩阵是为了降低理解门槛。运行后你会看到约束残差快速下降。如果前两步的|g|反而增大多半是符号或方向问题停下的那一刻就要检查初值。4.2 运动方程函数CR3BP的标准写法运动方程是动力学核心必须写对。下面是最常用的CR3BP无量纲速度-加速度形式function dstate cr3bp_rhs(t, state, mu) % state [x; y; z; vx; vy; vz] x state(1); y state(2); z state(3); vx state(4); vy state(5); vz state(6); r1 sqrt((x mu)^2 y^2 z^2); r2 sqrt((x - 1 mu)^2 y^2 z^2); ax 2*vy x - (1-mu)*(xmu)/r1^3 - mu*(x-1mu)/r2^3; ay -2*vx y - (1-mu)*y/r1^3 - mu*y/r2^3; az - (1-mu)*z/r1^3 - mu*z/r2^3; dstate [vx; vy; vz; ax; ay; az]; end注意ax和ay中包含了科氏力项2vy和-2vx。如果你把运动方程写成二阶加速度形式漏掉这两项是最大翻车点之一。r1和r2分别是航天器到两个主天体的距离这里用的是质心坐标系。mu是质量比无量纲后两个主天体位于 (-mu,0,0) 和 (1-mu,0,0)。这个函数会被ode45反复调用所以不要在里面写任何全局变量或IO操作保持纯粹。4.3 打靶迭代函数修正z与vy的细节接下来是事件函数和打靶约束函数。事件函数让积分停在y0的穿越面。这里我加了时间判断避免初始时刻的y0被当成事件触发function [value, isterminal, direction] event_xz(t, state) value state(2); if t 1e-10 value 1; % 忽略初始时刻的y0 end isterminal 1; % 触发后停止 direction 1; % 只捕获y从负到正的穿越 end function g halo_constraints(guess, mu, xL2, Ax) z0 guess(1); vy0 guess(2); state0 [xL2 - Ax; 0; z0; 0; vy0; 0]; opts odeset(Events, (t,s) event_xz(t,s), RelTol, 1e-12, AbsTol, 1e-12); [t, state, te, ye, ie] ode45((t,s) cr3bp_rhs(t,s,mu), [0 10], state0, opts); if isempty(ie) error(未穿越xz平面请检查初值或积分时长); end vx_end state(end,4); vz_end state(end,6); g [vx_end; vz_end]; enddirection1的选择对应从xz平面出发后轨道先进入y0一侧然后从负到正返回。如果程序报错未穿越xz平面优先尝试direction-1。有限差分雅可比函数如下function J finite_diff_jacobian(guess, mu, xL2, Ax) delta 1e-8; g0 halo_constraints(guess, mu, xL2, Ax); gp halo_constraints([guess(1)delta; guess(2)], mu, xL2, Ax); gm halo_constraints([guess(1)-delta; guess(2)], mu, xL2, Ax); gq halo_constraints([guess(1); guess(2)delta], mu, xL2, Ax); gqm halo_constraints([guess(1); guess(2)-delta], mu, xL2, Ax); J(:,1) (gp - gm)/(2*delta); J(:,2) (gq - gqm)/(2*delta); end这里用了中心差分比单边差分误差小。每一步迭代会调用5次约束函数含原始一次每次都要积分半圈整体计算量不大。如果嫌慢可以把delta取1e-7但精度会略降。不要取1e-10因为ode45的数值噪声会淹没差分信号。4.4 绘制三维轨道与关键检查收敛后用同样的初始状态积分若干个完整轨道周期然后画图。Halo轨道的周期可以从事件间隔估计。我一般积分5到6圈看轨迹是否闭合state0_final [xL2 - Ax; 0; guess(1); 0; guess(2); 0]; opts odeset(RelTol,1e-12, AbsTol,1e-12); T_guess 2.6; % 地月L2附近Halo无量纲周期约为2.6左右可以先初估 [t, state_traj] ode45((t,s) cr3bp_rhs(t,s,mu), [0 5*T_guess], state0_final, opts); plot3(state_traj(:,1), state_traj(:,2), state_traj(:,3), LineWidth, 1.2); xlabel(x); ylabel(y); zlabel(z); axis equal; grid on;注意T_guess只是初估严格周期需要从事件时间计算。你也可以直接积分[0 30]看y轨迹是否始终有规律。如果漂移明显说明精度不够或收敛精度不够。画图前最好先检查一下state_traj有没有NaN我遇到过精度设置太低导致积分步长崩溃的情况。5. 设计Halo轨道最常见的5个翻车现场与排查方案5.1 现象初值看着合理但微分改正发散迭代几步后|g|不降反升甚至出现NaN。原因通常是z0和vy0的初始组合落在了另一个解流形上或者事件函数的direction与轨道实际穿越方向不匹配。我遇到过一次因为z0取了负值而事件direction1导致事件根本没触发。解决方法是先检查事件是否触发、观察穿越时刻的y和速度符号。如果direction设反了直接把direction1改为direction-1。另一种可能是积分时间太短轨道没来得及穿越平面可以把积分上限从10改到20同时用事件控制积分上限只起保护作用。注意打印te事件时间如果te接近0说明初始点也被当成事件了检查事件函数的时间保护。5.2 现象轨道画出来是闭合的但YZ振幅不符合设计值收敛到闭合轨道没错但z方向的振幅不是你想要的值。这是因为你固定的是x方向振幅Ax而yz平面的振幅是由系统动力学确定的不能单独指定。如果任务要求特定z振幅就得把固定量改成Az把Ax作为未知数参与迭代。常见做法是重写约束函数未知数改为 [Ax; vy0] 或 [Ax; z0]约束仍然是终点vx和vz等于0。改起来不复杂但你要明白Halo轨道的振幅参数是耦合的。这是新手最容易把设计两个字理解错的地方。5.3 现象积分几十圈后轨迹开始螺旋漂移短期几圈没问题积分长了就漂。这几乎都是积分器精度不够。ode45默认的RelTol1e-3对短期可视化可能看不出问题但对周期轨道是致命伤。另外收敛残差也要严格到1e-10以上。否则轨道状态存在微小切向误差几十圈后会沿流形方向发散。把RelTol和AbsTol都设到1e-12同时把微分改正的收敛门限设为1e-10漂移现象会基本消失。如果还漂就换用ode113或显式RK8但这通常不是刚需。5.4 现象打靶收敛到平面Lissajous轨道而不是三维Halo这太常见了。你给了一个非零的z0初值但迭代结束后z0收敛到很小的值画出来是一条扁平的类似8字的轨道本质是平面Lissajous。原因是x方向振幅Ax太小没有达到Halo轨道的临界振幅。地月L2的临界Ax大约在0.01量级如果Ax小于这个值非线性耦合不足以把z方向运动锁定成垂直方向的大幅振荡。解决方法是把Ax调大一倍再试同时注意收敛后的z0是否明显非零。如果z0收敛到0基本就是Lissajous了。也可以用Jacobi常数来检查但最直观的还是看三维图里z方向的包络。5.5 现象ode45积分太慢长时间积分卡死Halo轨道设计里打靶迭代要反复积分如果精度要求高ode45在穿越平面附近会步长变得很小速度慢得让人想砸电脑。我遇到过一次单次积分耗时几秒整个迭代要一分钟。解决方法是先用低精度RelTol1e-8做粗收敛等约束残差降到1e-6后再把精度提高到1e-12做精修正。另一个技巧是给ode45一个合适的初值步长或者改用ode113在光滑问题上更快。不要一开始就上1e-12那会让你的迭代在早期白白浪费大量时间。6. 从设计到验证用周期连续性检验你的Halo轨道6.1 多圈积分与状态偏差统计确认轨道可用的第一步是把收敛后的初值积分至少10个周期然后对比每个周期结束时的状态和初始状态。我习惯的做法是计算最后状态与初始状态的六维向量差并记录任意两组周期穿越点之间的最大偏差。% 用事件函数捕捉每次穿越xz平面的时刻和状态 opts odeset(Events, (t,s) event_xz(t,s), RelTol, 1e-12, AbsTol, 1e-12); [t, state, te, ye, ie] ode45((t,s) cr3bp_rhs(t,s,mu), [0 30], state0_final, opts); % 对比最后一次穿越与前一次同方向穿越的偏差 dev norm(ye(end, :) - ye(end-1, :));dev在1e-9量级说明轨道具备很好的周期性。注意这里事件函数的方向是固定的所以相邻同方向事件间隔正好是一个完整轨道周期。如果dev偏大回到微分改正把收敛门限收紧到1e-11并用ode113重新积分验证。6.2 验证通过后还能做什么轨道延拓、太阳帆与转移段Halo轨道设计不是终点它往往是更大任务设计的第一块拼图。你可以在收敛的Halo轨道上做参数延拓连续改变mu或振幅得到一族Halo轨道这对工程选轨特别有用。我一般会写一个简单的延拓循环以上一个收敛解作为下一个初值能够稳定地扫出一条轨道族。另外如果你需要从地球附近进入这条Halo轨道可以借助不变流形理论设计转移轨道。CR3BP的不变流形实际上是对Halo轨道初值做小扰动后向前或向后积分的轨迹簇在MATLAB中实现并不复杂。把这里的运动方程和事件检测保留下来直接在扰动后的状态上再积分就能得到流形结构。这也是很多项目里从设计Halo到任务设计之间的最短路径。我在做某跨平台系统项目时曾经因为没有仔细检查事件方向连续发了三个晚上才定位到是direction的问题。后来我养成了一个习惯每写一个打靶迭代第一步不是看残差而是打印事件时刻和事件状态确认积分确实停在我想要的位置上。这个习惯帮我省了无数时间。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
