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、升交点赤经Ω、近地点幅角ω、真近点角θ）适合描述轨道形状和朝向，但直接用于数值积分时，靠近i=0或e=0会出现奇异，所以计算时通常把六根数转换成惯性系下的位置矢量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，否则向量会落在错误的象限。反变换时则由位置速度计算角动量矢量，再分解出轨道根数。

注意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; % 中心引力加速度 end

onorbit.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, X+dt/2*k1); k3 = omgd(t_list(k)+dt/2, X+dt/2*k2); k4 = omgd(t_list(k)+dt, X+dt*k3); X = X + dt/6*(k1+2*k2+2*k3+k4); 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-X（3-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)); end

testangel.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); % 阻尼角速度，抑制章动 end

controlstar.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中的更新函数 end

visual.m里绘制姿态时，会调用angel.m得到欧拉角，然后画一个旋转过的三轴坐标系。注意动画帧率和实际物理时间的区别：仿真dt是10秒，动画每帧显示0.05秒，看到的速度不代表真实速度。我会在标题栏写“仿真时刻t=xxx秒”，这样看动画时不会误解运动速度。

如果需要把动画保存成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，只是四元数积分每步之后要重新归一化，避免模长漂移。把这套结构跑通后，再替换控制器或执行机构模型就非常方便了。

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