MATLAB卫星仿真：从轨道动力学到姿态控制的完整闭环

简介卫星建模、卫星控制、轨道动力学与姿态动力学方向的 MATLAB 实验代码包聚焦卫星系统仿真与控制中的轨迹递推、姿态稳定与机动策略内容从基础轨道力学延伸到复杂姿态控制策略。压缩包共 22 个 .m 文件整体大小约 10KB全部为 MATLAB 脚本包含主程序、轨迹绘制、角动量计算、轨道机动与姿态控制等模块脚本间通过参数传递互相调用便于按需修改和验证。代码覆盖轨道根数计算、姿态运动学与动力学建模、喷气/飞轮控制等典型环节运行后可直观观察卫星位置、速度及姿态的演化通过调整初始参数或控制增益还能对比不同控制策略对轨道保持与姿态定向的影响适合开展课程设计或毕业设计专题。资源目前已有 568 人学习航天、飞行器设计相关专业学生或工程师可据此进行仿真验证与二次开发。1. 卫星建模与控制实验代码从main.m到轨道、姿态耦合仿真拿到“实验代码.zip”这份压缩包我建议先别急着运行main.m。里面plottrace.m、controlstar.m、omgd.m、onorbit.m这些文件名基本能看出来是一套用MATLAB写的卫星仿真框架轨道部分负责回答“卫星往哪飞”姿态部分负责回答“卫星朝哪看”最后用visual.m和plottrace.m把结果画出来。对做卫星控制算法验证的人来说这套代码的价值不是某个单独函数而是它把轨道动力学、姿态动力学、控制器串成了一个完整闭环。它适合刚接触航天仿真的研究生快速建立整体概念也适合需要搭六自由度仿真平台的工程师做模块替换。这一章不展开具体算法先把文件结构和数据流理清后面几章再一个模块一个模块地拆。2. 轨道动力学与轨道控制解析omgd.m、onorbit.m与轨道递推实现2.1 状态量与坐标系transio.m和transoi.m的坐标变换卫星轨道仿真第一步是坐标变换。轨道六根数半长轴a、偏心率e、轨道倾角i、升交点赤经Ω、近地点幅角ω、真近点角θ适合描述轨道形状和朝向但直接用于数值积分时靠近i0或e0会出现奇异所以计算时通常把六根数转换成惯性系下的位置矢量r和速度矢量v。transio.m和transoi.m就是这一对正反变换。我在工程里见过的典型接口是这样% 轨道六根数 - 惯性系位置、速度 [r_eci, v_eci] transio(a, e, i, OMEGA, omega, theta); % 惯性系位置、速度 - 轨道六根数 [a, e, i, OMEGA, omega, theta] transoi(r_eci, v_eci);这段代码的逻辑是先在近焦点坐标系中构造位置速度再绕近地点幅角、轨道倾角、升交点赤经做三次旋转得到J2000惯性系下的状态。正变换最关键的一步是用真近点角求偏近点角时要用atan2而不是acos否则向量会落在错误的象限。反变换时则由位置速度计算角动量矢量再分解出轨道根数。参数物理含义常用单位备注a半长轴km圆轨道时为半径e偏心率无量纲0为圆0~1为椭圆i轨道倾角rad0为赤道面OMEGA升交点赤经rad指春分点方向omega近地点幅角rad升交点到近地点夹角theta真近点角rad瞬时角度随时间变化注意MATLAB的三角函数默认弧度。很多人在这里用角度制导致transio算出来的轨道形状完全不对半长轴看起来正常但远地点方向偏掉几十度。我一般会在所有脚本开头统一用deg2rad转换并在函数入口处检查输入范围超过2π就直接报错。reserve.m从名字看是“保留/预留”在这类实验代码里通常承担保存初始化参数的职责。把轨道根数、仿真时长、步长集中放在一个脚本里比在main.m里散着定义更容易排查。实测中我只改reserve.m的参数不碰其他文件就能在testorbit.m里观察到不同初始轨道下的覆盖变化。2.2 轨道递推omgd.m与onorbit.m的分工omgd.m从命名看像是“orbit motion gravity dynamics”的缩写我把它理解为轨道动力学右函数也就是给定当前时刻t和状态向量X返回状态导数。二体模型下加速度只包含中心引力项function dX omgd(t, X) mu 398600.4418; % 地球引力常数km^3/s^2 r X(1:3); v X(4:6); r_norm norm(r); dX zeros(6,1); dX(1:3) v; dX(4:6) -mu / r_norm^3 * r; % 中心引力加速度 endonorbit.m则是在这个右函数外面套一层积分器。常见做法是直接用ode45但对轨道长期预报ode45的变步长控制在高度偏心轨道上会把步长压得很小效率不高。我通常换成定步长RK4逻辑比ode45更可控function [t_list, X_list] onorbit(X0, tspan, dt) t_list tspan(1):dt:tspan(2); X_list zeros(6, length(t_list)); X X0(:); for k 1:length(t_list) X_list(:,k) X; k1 omgd(t_list(k), X); k2 omgd(t_list(k)dt/2, Xdt/2*k1); k3 omgd(t_list(k)dt/2, Xdt/2*k2); k4 omgd(t_list(k)dt, Xdt*k3); X X dt/6*(k12*k22*k3k4); end end这里dt的选择要跟轨道周期匹配。近地轨道周期约90分钟dt取1秒到10秒都能得到稳定结果。如果dt超过30秒数值粘性会把轨道能量慢慢耗散直观表现是远地点高度逐渐下降。要验证积分器是否可靠可以跟踪轨道比机械能E v^2/2 - mu/r在纯二体模型下这个值应该保持不变。testomgd.m就是干这个的。它先调用transio生成初始状态再用onorbit递推若干周期最后打印能量偏差。如果能量偏差超过万分之一先检查mu和半径单位是否一致再检查积分步长。还有一个常见问题是把G和地球质量分开写导致mu算错轨道周期会成倍数漂移。如果需要在近地轨道精度更高可以在omgd.m中叠加J2摄动项。J2加速度的量级是中心引力的千分之一但对长期轨道预报影响很大尤其是升交点赤经漂移。加入J2后倾角i和升交点赤经Ω的长周期变化就能被模拟出来这对设计太阳同步轨道非常关键。2.3 轨道控制策略testorbit.m里为什么先看Δv轨道控制首先要算清速度增量Δv而不是直接给推力。霍曼转移是最基础的两圆轨道间最经济转移方式从低轨h1到高轨h2两次冲量Δv1和Δv2的关系如下% 霍曼转移计算 R_EARTH 6378.137; % km r1 R_EARTH h1; r2 R_EARTH h2; a_trans (r1 r2) / 2; v1 sqrt(mu / r1); v2 sqrt(mu / r2); vt1 sqrt(2*mu/r1 - mu/a_trans); % 转移椭圆近地点速度 vt2 sqrt(2*mu/r2 - mu/a_trans); % 转移椭圆远地点速度 delta_v1 vt1 - v1; delta_v2 v2 - vt2; delta_v_total abs(delta_v1) abs(delta_v2);testorbit.m如果把这段代码放在循环里不断改变h1、h2并绘制曲线能看到一个反直觉的结论从200km升到35786km所需的Δv大约3.9km/s而从35786km再往高处升单位高度需要的Δv反而下降。这说明低轨附近的轨道机动代价最高所以在实际任务中低轨卫星很少大幅调轨。实际工程里推力器点火不是瞬时的有限推力会让转移轨道偏离理想霍曼椭圆所以仿真需要把控制加速度加进omgd.m的加速度项中而不是直接改状态量。比如切向推力建模为thrust_acc F / (m0 - mdot * t); % 单位m/s^2需转换到km/s^2 dX(4:6) dX(4:6) thrust_acc * v_hat;其中v_hat是速度方向单位向量。注意这里的单位一致性omgd.m中位置单位是km那么推力加速度也要从m/s^2除以1000后加入。很多实验代码在这一步漏掉单位转换导致轨道半长轴出现完全不合理的跳动。控制任务常用策略执行机构备注轨道高度调整霍曼转移/连续小推力化学推进/电推进需要两个点火点轨道倾角修正在升交点或降交点垂直点火化学推进倾角变化大时代价很高相位调整先降轨再升轨利用周期差推进器用于同轨道面编队编队保持相对轨道要素控制冷气/电推进需要高精度轨道预报碰撞规避沿速度方向微小冲量推进器需结合轨道预报与误差分析3. 姿态动力学与稳定控制从angel.m到controlstar.m的实现3.1 姿态描述与测量angel.m和testangel.m的欧拉角计算姿态仿真的第一步是定义“姿态角”。angel.m命名上应该是angle但代码里常见这种笔误计算的是卫星本体坐标系相对轨道坐标系的欧拉角。轨道坐标系z轴指向地心x轴沿速度方向y轴垂直轨道面卫星的滚动角φ、俯仰角θ、偏航角ψ就是本体轴相对这个参考系的转角。默认旋转顺序为Z-Y-X3-2-1时由姿态旋转矩阵C反解欧拉角的代码通常写成function [phi, theta, psi] angel(C) % C为3x3姿态旋转矩阵从参考系到本体系 theta asin(-C(1,3)); phi atan2(C(2,3), C(3,3)); psi atan2(C(1,2), C(1,1)); endtestangel.m的作用是把已知欧拉角正算得到旋转矩阵再反算回去检查闭环误差。注意theta接近90度时会出现万向节锁此时phi和psi无法唯一确定atan2给出的值可能不连续。对地定向卫星的正常姿态theta不会到90度但做全姿态机动或失效模式仿真时一定要处理这个边界。工程上我更推荐用四元数做状态量和积分欧拉角只作为显示和遥测数据。四元数没有奇异点用四元数反解欧拉角的代码在visual.m中做可视化显示这样既能避免奇异性又能让人看懂姿态。3.2 姿态动力学方程Lx.m、Ly.m、Lz.m的力矩模型姿态动力学用欧拉方程描述I·ω_dot ω×(I·ω) M_ext M_ctrl。其中ω是本体角速度I是惯量张量M_ext是环境力矩。Lx.m、Ly.m、Lz.m这三个文件我倾向于认为它们分别计算x、y、z轴的力矩分量。这里的关键是环境力矩模型要分清楚类别重力梯度力矩近地轨道最大的持续干扰之一由惯量差和重力场梯度共同产生。地磁力矩磁力矩器与地磁场相互作用产生控制力矩同时也是一个干扰源。气动阻力400-600km低轨不可忽略产生与质心压心差相关的力矩。太阳光压力矩高轨或大帆板卫星明显与表面积和反射系数有关。重力梯度力矩的简化实现可以写成function [Mx, My, Mz] gravity_gradient_torque(I, R_eci, C_body2eci) mu 398600.4418; r_norm norm(R_eci); n2 mu / r_norm^3; % 将地心矢量转到本体系 r_body C_body2eci * (R_eci / r_norm); % 重力梯度力矩矢量式3*n^2 * (r_hat × (I * r_hat)) M_gg 3*n2 * cross(r_body, I * r_body); Mx M_gg(1); My M_gg(2); Mz M_gg(3); end实际项目里Lx.m、Ly.m、Lz.m应该接收姿态矩阵和轨道位置把重力梯度力矩从轨道系转到本体系再叠加。如果直接写成常值会忽略姿态变化对力矩方向的调制导致控制律仿真结果过于乐观。我见过有人用固定姿态角跑完整个轨道周期结果重力梯度力矩的正负号变化都没体现出来。验证姿态动力学模型最直接的方法是检查角动量守恒。在没有外力矩的纯二体姿态仿真中总角动量H I·ω 在惯性系下应保持不变。如果H的模长在变化说明力矩计算里混入了数值误差或坐标系旋转错误。3.3 姿态控制律vecmo.m与controlstar.m的执行机构配合有了力矩模型还要有控制律。vecmo.m从命名看更像是“vector momentum”或“magnetic torque”相关函数。用磁力矩器做速率阻尼时最常用的是B-dot控制律function M_cmd vecmo(B_body, B_body_pre, dt, k_gain) % 输入地磁场矢量在三个时刻的值输出磁力矩指令 B_dot (B_body - B_body_pre) / dt; M_cmd -k_gain * cross(B_body, B_dot); % 阻尼角速度抑制章动 endcontrolstar.m则更像是整体姿态控制的入口它把期望姿态与实际姿态的误差转换为控制力矩。工程上最常见的是PD控制% 期望姿态角和角速度可以由星敏或飞轮目标生成 attitude_error angel_desired - angel_current; omega_error omega_desired - omega_current; M_ctrl -Kp * attitude_error - Kd * omega_error; % 加上执行机构饱和限制 M_ctrl max(min(M_ctrl, M_max), -M_max);Kp和Kd不能随便给要参考转动惯量。惯量大的轴用小的Kp否则容易激发挠性振动。控制力矩限制也要加进去如果是磁力矩器力矩只能在与磁场垂直的平面内产生必须把期望力矩投影到可用方向如果是飞轮则要额外考虑动量管理和饱和卸载。testmo.m可以放这组控制器的闭环测试。把初始姿态设置成10度目标为0度看姿态角和角速度在饱和之后的收敛曲线。稳定时间、超调量、稳态误差三个指标一起看只看一个会漏问题。比如某些参数组合下稳态误差很小但超调量超过30度会导致实际任务中天线指向完全失锁。4. 可视化与结果验证plottrace.m、trace.m与visual.m的工程用法4.1 计算与绘图分离trace.m计算轨迹plottrace.m只管画很多仿真脚本把计算和画图混在一起导致改一个坐标单位全部重来。这套代码里面trace.m和plottrace.m分开设计上是正确的。trace.m返回轨迹数组plottrace.m接收轨迹并绘制。% main.m中的典型调用 X0 [r0; v0]; [t_list, X_list] onorbit(X0, [0, 5400], 10); % 1.5小时轨道 plottrace(X_list, t_list);plottrace.m里我建议至少画两个图三维轨迹图和轨道高度随时间变化图。三维图用plot3地面轨迹投影到经纬度平面则用geoplot或自己画经纬度网格。轨道高度变化曲线最能暴露模型错误如果高度线性衰减说明有人造阻尼项如果是周期波动说明摄动力模型未归一化如果出现阶跃跳变基本可以确定是坐标系旋转出错。地面轨迹投影时需要先算卫星的经纬度。给定ECI位置通过格林尼治恒星时转换到ECEF再反算经纬度。这个过程在代码里要写清楚当前时刻对应的格林尼治恒星时否则绘图时卫星会“跑偏”经度。4.2 动态可视化movement.m与visual.m的联合使用姿态控制结果适合用动画展示。movement.m负责逐帧更新卫星位置visual.m负责把姿态画成旋转的三轴坐标系。用MATLAB的animatedline做轨迹动画比每次都set(XData)更高效figure; h animatedline(Color, b, LineWidth, 1.5); axis equal; grid on; for k 1:size(X_list,2) addpoints(h, X_list(1,k), X_list(2,k), X_list(3,k)); drawnow limitrate; % 同时更新姿态三轴调用visual.m中的更新函数 endvisual.m里绘制姿态时会调用angel.m得到欧拉角然后画一个旋转过的三轴坐标系。注意动画帧率和实际物理时间的区别仿真dt是10秒动画每帧显示0.05秒看到的速度不代表真实速度。我会在标题栏写“仿真时刻txxx秒”这样看动画时不会误解运动速度。如果需要把动画保存成GIF或视频不建议直接用animatedline的drawnow截屏。更好的做法是循环中先定位当前数据点执行movement函数更新位置再用exportgraphics或VideoWriter写入。GIF文件要注意尺寸和帧率一般1024x768、10fps就够了太大文件会超过100MB。4.3 验证模块间的数据流testangel.m、testomgd.m、reserve.m的测试链这套代码里有不少test前缀文件它们的作用不是摆设。testomgd.m验证轨道积分器testangel.m验证姿态变换testmo.m验证力矩控制器。我按依赖关系把它们排成一条链先通过testomgd确认轨道递推稳定再用testangel确认姿态测量无偏差最后跑testmo确认控制环路收敛。任何一层不过单独调试对应模块不要带着错误数据跑完整仿真。reserve.m如果内容是保存工作区建议改成在main.m末尾自动保存关键变量避免每次手动敲save命令。保存时带上时间戳save(sprintf(results_%s.mat, datestr(now,yyyymmdd_HHMMSS)), X_list, t_list, X0, a, e, i);文件命名带上日期后再回看仿真结果就能知道是哪一版代码跑的。重视结果可追溯性对写论文和工程评审都很有帮助。另外验证时要把二体模型下的轨道周期和理论值对比。圆轨道理论周期T 2π*sqrt(a^3/mu)如果仿真得到的星下点重复周期与理论值差超过1%就要回头查积分器。5. 进阶把实验代码改造成批量参数扫描与姿态-轨道联合仿真5.1 用parfor扫描控制参数寻找最优增益如果要用这套代码做控制参数整定最直接的做法是把main.m里的Kp、Kd改成输入参数外面再套一个批量循环。MATLAB的parfor可以让多组参数组合并行跑再把结果汇总成表格Kp_list [0.01, 0.05, 0.1]; Kd_list [0.1, 0.5, 1.0]; results table(); for i 1:length(Kp_list) for j 1:length(Kd_list) [settle, overshoot, steady] run_single_sim(Kp_list(i), Kd_list(j)); results [results; table(Kp_list(i), Kd_list(j), settle, overshoot, steady)]; end end扫描前先确保单次仿真能在1秒内跑完再考虑并行。用profile查看哪个函数耗时最长轨道积分通常是瓶颈可以把omgd.m编译成mex或减少输出点数。输出点数只保留每个轨道周期的极值而不是每一积分步都存内存占用能降一个数量级。5.2 让姿态控制力矩参与轨道递推形成六自由度闭环大多数情况下轨道和姿态是分开仿真的轨道用onorbit.m姿态用controlstar.m。但如果要模拟推力器点火对姿态的干扰或者验证推力矢量控制就需要联合仿真。做法是每一仿真步先由轨道计算期望姿态姿态控制律输出控制力矩再把控制力矩对卫星质心的作用换算成合力与合力矩喷气推力则要把推力加速度加到omgd.m的加速度项中。此时Lx.m、Ly.m、Lz.m不只是环境力矩模型还要叠加控制力矩接口。把联合逻辑封装成一个大右函数sixdof_dynamics(t, X)再统一调用onorbit.m比在两个循环里同步时钟要省心得多。状态向量扩展为12维前6维是位置速度后6维是四元数或欧拉角和角速度。统一的积分器对刚体动力学方程也能用RK4只是四元数积分每步之后要重新归一化避免模长漂移。把这套结构跑通后再替换控制器或执行机构模型就非常方便了。本文还有配套的精品资源点击获取