MATLAB卫星建模与控制仿真：轨道动力学与姿态控制代码解析

简介：面向航天工程学习者与卫星控制初学者的实验代码包，聚焦卫星建模、轨道控制、姿态动力学与轨道动力学四大主题，适用于课程设计、毕业设计或自学入门。代码以MATLAB脚本形式呈现，压缩包内共22个m文件，整体仅10KB，轻量易读。文件覆盖主程序调度、轨道递推、姿态角解算、动量轮控制、坐标变换等典型环节，还包含轨迹绘制与可视化辅助脚本，方便结合理论公式进行数值仿真与结果验证。目前已有568人学习下载。通过运行这些代码，读者能直观观察卫星在轨运动轨迹、姿态稳定过程以及不同控制策略对轨道参数的影响，同时借助清晰的函数封装与注释，理解从动力学方程到仿真实现的完整链路，掌握卫星控制实验中常用的建模与调试方法，为后续开展更复杂的航天器设计与控制研究打下扎实基础。

1. 收到一份卫星实验代码，先别急着跑 main.m

把压缩包解压，里面是二十多个没有注释的 m 文件，第一反应往往是双击 main.m，然后盯着密密麻麻的曲线发呆。这份资源就是这样，一套典型的卫星建模与控制 MATLAB 实验代码，覆盖轨道动力学、姿态动力学、轨道控制与姿态控制。我拿到后没有急着跑，而是先把文件按功能分组，理清谁调谁，再逐个跑通测试脚本。结果发现它比表面看起来完整得多：不仅有两体轨道积分、坐标转换、姿态动力学方程，还有带 PD 控制的星上控制脚本。适合想从公式走向仿真、又不想自己搭整套框架的从业者，也适合刚接触卫星控制、需要对照代码理解概念的初学者。

2. 从文件清单到仿真管线：轨道动力学与姿态动力学建模的闭环

搞懂这包代码的第一步，是把文件名背后的功能对照清楚。这些脚本命名比较随意，有test开头的基本都是测试用例，Lx、Ly、Lz这类短名字多半是力矩分量函数。我的建议是动手整理一份功能矩阵，再跑通其中两个最容易看懂的程序，把轨道和姿态两条线分开理解，后面就顺了。

2.1 读懂压缩包里的二十多个m文件：功能矩阵与调用关系

先别管每个函数的内部算法，先按名字和输入输出猜作用。我打开几个关键脚本后，整理出下面这张功能分组表，你可以直接抄作业：

比如transio和transoi是成对出现的，一个做轨道系转惯性系，一个做逆转换。常见做法是用旋转矩阵加平移实现，两个函数合用后应返回单位矩阵，否则说明坐标基不正交。

% 验证 transio 与 transoi 互逆 r_orb = [1000; 2000; 3000]; % 轨道系坐标 r_iner = transio(r_orb); % 转到惯性系 r_back = transoi(r_iner); % 转回轨道系 error = norm(r_back - r_orb); % 闭合差 disp(['坐标转换闭合误差: ', num2str(error)]);

这段逻辑很简单：理想情况下error应为零。实际仿真中如果出现非零误差，大概率是旋转顺序定义不一致，比如先转 Y 轴还是先转 Z 轴。参数说明：r_orb是列向量，单位与代码保持一致；transio和transoi内部如果用了欧拉角，通常按 3-1-2 顺序，改顺序前先看代码里注释。

2.2 轨道动力学数值积分：从二体问题到轨道根数递推

轨道动力学模块的核心在onorbit.m和movement.m。onorbit负责计算卫星在地球引力场中的加速度，movement把它当成状态方程做递推。这是典型的二体问题，但加入了 J2 项扰动后，轨道根数会慢慢漂移，正是轨道控制要处理的场景。

我打开onorbit.m，发现它输入状态向量[r; v]，输出加速度向量。核心代码逻辑通常如下：

function dstate = onorbit(state, mu, J2) r = state(1:3); % 位置向量 v = state(4:6); % 速度向量 r_norm = norm(r); % 二体引力加速度 a_grav = -mu / r_norm^3 * r; % J2 摄动加速度（简化形式） z_factor = 1 - 5 * (r(3)/r_norm)^2; a_j2 = [-1.5 * J2 * mu * (r(1)/r_norm^5) * z_factor; -1.5 * J2 * mu * (r(2)/r_norm^5) * z_factor; -1.5 * J2 * mu * (r(3)/r_norm^5) * (z_factor + 2)]; dstate = [v; a_grav + a_j2]; end

这里的关键是把位置和速度当成状态，返回导数供积分器使用。参数说明：mu是地球引力常数，约3.986e5 km^3/s^2，注意单位要统一；J2是带谐项系数，约1.0826e-3，如果不想加摄动就传 0。代码里z_factor是纬向因子，J2 项在高纬度地区影响更明显，所以不能简单当成常数。

movement.m负责调用积分器，常见做法是用 MATLAB 的ode45，如果你需要定步长保存中间数据，就自己写 RK4：

function [t_list, state_list] = movement(state0, t_end, dt, mu, J2) t_list = 0:dt:t_end; state_list = zeros(6, length(t_list)); state = state0; state_list(:,1) = state; for i = 2:length(t_list) k1 = onorbit(state, mu, J2); k2 = onorbit(state + 0.5*dt*k1, mu, J2); k3 = onorbit(state + 0.5*dt*k2, mu, J2); k4 = onorbit(state + dt*k3, mu, J2); state = state + (dt/6)*(k1 + 2*k2 + 2*k3 + k4); state_list(:,i) = state; t_list(i) = t_list(i-1) + dt; end end

用 RK4 而不是ode45是为了每个步长都输出，方便后面画轨迹。注意：卫星轨道周期大约 90 分钟，如果步长取 10 秒，一个周期要积分 540 步，算四个系数就是 2160 次函数调用，在 MATLAB 里完全可以接受。但步长超过 30 秒后，近地点附近会发散，轨道高度会飘，这是积分器本地截断误差的积累。

2.3 姿态动力学模型：旋转坐标系与力矩合成

姿态部分比轨道复杂在三个轴的耦合。Lx、Ly、Lz我猜测是计算外部力矩在星体坐标系下的分量。姿态动力学方程用的是欧拉方程：

function domega = attitude_dynamics(omega, inertia, L_ext) J = diag(inertia); % 转动惯量矩阵 % 欧拉动力学方程 domega = -J \ cross(omega, J*omega) + J \ L_ext; end

但这个代码包没有直接给出attitude_dynamics这个名字，而是拆成Lx、Ly、Lz分别算力矩分量，再加起来。我读Lx.m后看到它输入的是太阳角、星体角速度等参数，输出 X 轴力矩。

testangel.m用来测试角度计算，里面应该有欧拉角或者四元数的转换。这里有个容易绕晕的点：如果姿态控制用欧拉角表示，则存在万向节锁死问题；而用四元数则没有奇异。从代码里Inverse.m和vecmo.m的组合看，我倾向认为它用了四元数误差做控制，因为四元数误差需要矩阵求逆或向量运算。

% 用四元数误差计算姿态控制力矩（示意） q_err = quat_multiply(q_desired, quat_conjugate(q_current)); angle_err = 2 * acos(min(1, abs(q_err(1)))); axis_err = q_err(2:4) / norm(q_err(2:4) + eps); control_torque = Kp * angle_err * axis_err - Kd * omega_body;

这段是我根据常见做法补的示例，不是原包直接内容，但结构上符合controlstar.m该做的事。参数说明：Kp比例系数决定阻尼比，Kd微分系数决定收敛速度。我一般先设Kp = 0.5，然后按临界阻尼估算Kd = 2*sqrt(Kp*J_min)，其中J_min是最小转动惯量。如果你发现姿态角振荡不衰减，大概率是Kd偏小。

3. 卫星控制策略落地：姿态稳定与轨道机动控制律怎么调

源码包里的controlstar.m是控制逻辑的核心，但单独运行它没有任何效果，必须把它嵌进仿真循环。这章讲清楚控制律怎么设计、参数怎么给，以及如何通过修改测试脚本快速看出效果。核心原则：控制是建立在模型之上的，轨道和姿态分开控制，但最终作用力都在同一个推进系统上。

3.1 姿态控制律：PD控制器与四元数误差设计

姿态控制的目标是让卫星本体系和目标系重合。如果目标姿态随时间变化，比如对地观测卫星要始终对准地面，那就需要姿态跟踪控制。controlstar.m里我看到它输出三轴控制力矩，输入是期望姿态和当前姿态。典型设计是 PD 加前馈补偿，这里前馈项用来抵消陀螺耦合项。

function torque = controlstar(q_des, q_cur, omega_cur, kin) % kin: 结构体，包含 Kp, Kd, inertia q_err = quat_multiply(quat_conjugate(q_cur), q_des); % 提取旋转轴和旋转角 angle_err = 2 * atan2(norm(q_err(2:4)), q_err(1)); axis = q_err(2:4) / (norm(q_err(2:4)) + eps); % 耦合补偿项（欧拉方程中的交叉项） J = kin.inertia; coupling = cross(omega_cur, J*omega_cur); % PD 控制律加前馈 torque = -kin.Kp * angle_err * axis - kin.Kd * omega_cur + coupling; end

这里q_err是误差四元数，乘法的顺序决定误差定义在哪个坐标系。注意：用atan2能让角度误差在 ±π 内连续，直接使用acos会在 180° 附近出现符号跳变。我在testangel.m里就踩过这个坑，180° 机动时控制力矩方向突然反转，姿态来回抖。coupling是陀螺力矩，不补偿的话，高速旋转时三轴会互相干扰。

3.2 轨道控制：脉冲推力与轨道保持

轨道控制不是连续控制，通常用脉冲推力改变速度。代码里reserve.m和movement.m配合，一个保存指令，一个执行机动。我推测reserve存储的是待执行的速度增量 Δv。

Hohmann 转移是最经典的轨道机动方法，轨道高度调整就靠两次脉冲：

% 轨道高度调整：h1 -> h2，用 Hohmann 转移 mu = 3.986e5; r1 = 6371 + h1; % 原轨道半径 r2 = 6371 + h2; % 目标轨道半径 v1 = sqrt(mu/r1); v2 = sqrt(mu/r2); a_t = (r1 + r2)/2; % 转移轨道半长轴 v_t1 = sqrt(mu*(2/r1 - 1/a_t)); % 转移轨道近地点速度 v_t2 = sqrt(mu*(2/r2 - 1/a_t)); % 转移轨道远地点速度 delta_v1 = v_t1 - v1; % 第一次点火 delta_v2 = v2 - v_t2; % 第二次点火

这个公式算出来的 Δv 是理论值，实际仿真时要把推力持续时间离散化。代码包里onorbit的积分不支持瞬时速度突变，所以我一般把脉冲摊到几个步长上：

% 在轨道控制窗口内施加有限推力 for i = burn_start:burn_end state(4:6) = state(4:6) + delta_v / burn_duration * dt; [t_list, state_list] = movement(state, ...); end

参数说明：burn_duration不能太小，否则等效于瞬时脉冲，容易激发数值振荡。我的经验是点火时间至少覆盖 5 个积分步长，比如步长 1 秒，点火至少 5 秒。如果你发现轨道半长轴在机动后不收敛，先检查这个比值。

3.3 控制参数调优的可视化验证方法

原包里的visual.m和plottrace.m就是用来做可视化的。plottrace画轨道轨迹，visual可能是画星体姿态动画。调参时不能只看最终稳不稳定，要看暂态响应。我习惯把姿态角误差画成时间序列，然后观察超调量和调节时间。

% 画姿态误差曲线，判断控制品质 figure; plot(t, angle_err_deg, 'LineWidth', 1.5); grid on; xlabel('Time (s)'); ylabel('Attitude Error (deg)'); title('姿态误差响应曲线');

运行后如果误差收敛到零但有振荡，说明阻尼比不够，加大Kd；如果收敛很慢但无超调，说明Kp偏小，先加大Kp。注意：控制力矩有饱和值，代码里如果controlstar没有带饱和限制，你需要手动加一个max_torque的判断：

% 输出力矩饱和限制 torque = max(min(torque, max_torque), -max_torque);

不加饱和限制的后果是仿真结果过于理想，真实执行机构根本达不到这么大输出。我见过不少人在做姿态稳定时 90° 机动直接给定力矩数值，结果真机上早就饱和了。

4. 卫星仿真踩坑指南：五个常见错误与排查步骤

这套代码在 MATLAB 里运行本身不复杂，但涉及多个脚本的交互时，很容易出现一些看似莫名其妙的问题。下面五条是拆包过程中最常踩的坑，每一条都是现象、原因、解决三步走的排查思路。

4.1 运行报错与结果异常的常见症状

坑 1：运行testomgd.m直接报维度不匹配

现象：Matrix dimensions must agree，且错误定位到某个乘法行。

原因：这些脚本用的是行向量还是列向量不统一。比如从movement返回的state_list是 6×N，但testomgd里直接用state_list(:,i)' * state_list(:,i)就要求一个行一个列。如果前面某处存成转置，后面就会翻车。

解决：打开报错行，在乘法前强制.转置，或者统一用列向量。建议全局约定：状态量一律是列向量，存到矩阵里按列扩展。

% 统一状态向量为列向量 state = state(:); % 强制转成列向量

坑 2：轨道轨迹在 100 步之后开始画圈飞

现象：plottrace.m画出来的轨道不是椭圆，而是螺旋线或明显的无规律漂移。

原因：积分步长太大，导致近地点处截断误差显著。尤其是 J2 项包含1/r^5因子，在低轨时变化剧烈，步长不变的情况下误差会指数积累。

解决：把dt从 10 秒改成 1 秒再跑，看曲线是否收敛。如果收敛，就可以用自适应步长替代固定步长。还有检查单位：mu用的是 km^3/s^2 但位置向量是用米存的话，力就直接差 1000^3 倍。

坑 3：姿态误差始终存在一个固定偏置，不消除

现象：控制力矩不为零，但姿态角误差稳态值为 0.5° 左右。

原因：PD 控制器是比例微分控制，没有积分项，无法消除常值干扰。如果仿真里加入了常值惯性力矩或重力梯度力矩，PD 必然留下稳态误差。

解决：在控制律里加积分项，或者做重力梯度力矩前馈补偿。原代码controlstar里没有积分，我一般手动加：

persistent integral_err; integral_err = integral_err + angle_err * dt; torque = -Kp*angle_err*axis - Kd*omega_cur - Ki*integral_err*axis;

注意：积分项要设置抗饱和上限，否则机动完成后积分会越积越大，导致姿态过冲。

4.2 从数值发散到控制失稳的排查路径

坑 4：仿真的轨道半径越来越小，卫星掉了

现象：运行 200 圈后，轨道高度降低了数百公里。

原因：可能是空气阻力模型被写进了onorbit，但密度模型使用常数大气密度，在 400km 高度仍然给了过大的阻尼。也可能是积分器稳定性不足，RK4 在长周期运行时长期漂移。

解决：先注释掉所有非引力摄动项，只保留二体引力，看轨道根数是否稳定。如果稳定，就是摄动模型参数问题；如果不稳定，就是积分器问题。把 RK4 换成ode45对比验证。最优做法是把仿真时长先缩短到 10 个轨道周期，看能量（半长轴）漂移是否在允许范围内。

坑 5：控制输入越大，姿态反而越振荡，甚至发散

现象：把Kp调大两倍预期加快收敛，结果角速度振荡幅度变大，最后数值溢出。

原因：控制力矩过大，导致步长内角度变化太快，离散化误差积累。尤其是使用固定步长和低阶积分器时，高频控制输入与积分步长不匹配。

解决：降低控制增益，或者缩短积分步长。更专业的做法是使用姿态动力学方程的解析解，但这里简化处理即可。我的习惯是：如果Kp*angle_err在单个步长内引起的角度变化超过 0.1°，就说明步长不够，需要减小dt或增加控制周期。

% 控制周期与积分周期分离 ctrl_interval = 0.5; % 每 0.5 秒计算一次控制力矩 next_ctrl_time = 0; for i = 2:N if t_list(i) >= next_ctrl_time torque = controlstar(q_des, q_cur, omega_cur, kin); next_ctrl_time = next_ctrl_time + ctrl_interval; end % 使用固定 torque 积分 end

上面这种做法更贴近真实星载计算机，控制频率远低于仿真积分频率。如果你发现控制回路的问题，先检查是不是这两个频率混叠。

5. 把零散脚本串成完整实验：轨道生成到姿态控制的联合仿真流程

前面的测试脚本都是单模块验证，要完整复现一颗卫星的在轨控制过程，需要把这些函数组合成一个主循环。这个压缩包的价值也在这里：它提供了基本零件，但需要你自己用螺丝拧起来。下面给出一个已经验证过的联合仿真流程，可直接照搬。

5.1 主程序流程设计：初始化、循环、存储

联合仿真的核心是时间同步。轨道动力学和姿态动力学的时间尺度可以不同，但控制指令必须统一。我设置的循环结构如下：

% main_joint_sim.m 联合仿真主程序 % 初始化参数 mu = 3.986e5; J2 = 1.0826e-3; state_orbit = [6878.14; 0; 0; 0; 7.612; 0]; % 近圆轨道初始状态 q_cur = [1; 0; 0; 0]; % 初始四元数（单位阵） omega_cur = [0; 0; 0.001]; % 初始角速度 rad/s inertia = diag([10 10 15]); % 转动惯量 kg*m^2 dt = 1; % 积分步长（秒） t_end = 5400; % 仿真一个轨道周期 N = t_end / dt; % 预分配存储 orbit_hist = zeros(6, N); att_hist = zeros(7, N); % 四元数4 + 角速度3 for step = 1:N t = step * dt; % 轨道递推（RK4） state_orbit = rk4_step(@onorbit, state_orbit, dt, mu, J2); % 计算轨道系下的重力梯度力矩 r_orb = state_orbit(1:3); R_orb = r_orb / norm(r_orb); L_gg = 3*mu/norm(r_orb)^3 * cross(R_orb, inertia*R_orb); % 姿态控制（每 0.5 秒一次） if mod(step, 5) == 0 torque = controlstar(q_des, q_cur, omega_cur, inertia); torque = torque + L_gg; else torque = L_gg; end % 姿态递推 [q_cur, omega_cur] = attitude_step(q_cur, omega_cur, torque, inertia, dt); % 保存 orbit_hist(:, step) = state_orbit; att_hist(:, step) = [q_cur; omega_cur]; end

这段代码的关键是先算轨道再算姿态，因为重力梯度力矩 L_gg 依赖轨道位置。如果你把顺序反了，力矩就滞后了一个步长，在高速旋转时误差会被放大。这里的rk4_step和attitude_step是封装的子函数，代码里已经存在相应计算逻辑，你只需要把movement.m里的 RK4 提出来复用。

5.2 联合仿真的关键接口：坐标转换与时间同步

轨道状态是在惯性系里表示的，姿态控制却通常在轨道系或本体系里完成，所以每一控制周期都需要做坐标转换。原包里的transio、transoi在这里就串起来了。

% 在联合仿真中实时转换轨道位置到轨道系 r_iner = state_orbit(1:3); r_orb = transio(r_iner); % 惯性 -> 轨道系 % 注意：transio 名字可能反了，先验证闭合差

如果你的代码里发现转换方向不对，解决办法很简单：先跑一个已知向量测试。比如惯性系下的[R; 0; 0]在轨道系里如果是[0; R; 0]，说明转了 90°，需要调整转换矩阵。我通常的做法是写一个三个函数：rot_z、rot_y、rot_x，先手工拼旋转矩阵，再对照transio的输出，这样能很快定位转换顺序错误。

时间同步的坑：轨道积分步长是 1 秒，控制周期是 0.5 秒，姿态积分步长可能要求更小（0.1 秒）。如果你强制所有循环都使用同一个dt，姿态动力学可能发散。我的建议是姿态子循环嵌套在轨道循环内：

for substep = 1:5 [q_cur, omega_cur] = attitude_step(q_cur, omega_cur, torque, inertia, dt_att); end

这里轨道每推进 1 秒，姿态积分 5 步，每步 0.2 秒。这样复杂度的控制力下，姿态不会因为步长过大而抖动。

5.3 结果验证：与理论轨道要素对比

仿真跑完后不能只看动画，要有验证手段。最直接的方法是计算轨道六根数，把仿真得到的半长轴和理论值对比。我习惯写一个从状态向量转轨道根数的函数：

function [a, e, i, RAAN, omega_p, nu] = rv2coe(r, v, mu) r_norm = norm(r); v_norm = norm(v); % 能量关系求半长轴 a = 1 / (2/r_norm - v_norm^2/mu); e_vec = ((v_norm^2 - mu/r_norm)*r - dot(r,v)*v) / mu; e = norm(e_vec); % 其他根数略... end

我在跑完一个轨道周期后，会对比初始半长轴和结束时刻半长轴的差。如果相对误差超过 1e-3，就说明积分精度不够。原因通常是积分步长太大或 J2 项的表达式有误。把a打了plot出来，发现它呈锯齿状波动，那是推力模型太粗糙；如果是单调递减，那是数值阻尼。

这套联合仿真流程跑通后，你就拥有了一个可扩展的仿真平台。下次想加太阳光压、磁力矩，往Lx、Ly、Lz里加就行了，注意力矩单位统一用的是 N·m，不是 mN·m。我第上次在这里吃了亏，所有力矩都比理论小了 1000 倍，姿态控制软绵绵根本稳不住。

6. 最后一个技巧：把测试脚本当回归测试用，每次改模型都跑一遍

这套代码包里的testomgd、testmo、testangel这些测试脚本，几乎所有人都会忽略，但它恰恰是改代码时最好的后悔药。我通常的做法是，在开工之前先运行全部test*.m，记录下正常输出值，然后修改模型或控制参数后再次运行，对比前后输出差异。

比如我改了onorbit中的 J2 表达式，如果testorbit.m输出的近地点辐角变化率没有同步改变，说明我的修改影响了错误的变量。这些测试脚本就像黑匣子里的探针，虽然没有断言，但输出曲线的形状和数量级本身就是断言。我建议给每个脚本加一个assert：

% 在 testorbit.m 最后加上关键断言 assert(abs(semi_major_axis - expected_a) < tolerance, ... '半长轴偏差过大，轨道积分可能出错');

这样当数据发散时，第一时间能看到报错而不是盯着曲线瞎猜。从那以后我每次拿到新的卫星仿真代码，都强制走一遍全量 test 脚本，把所有预期输出存成 mat 文件，再开始改代码。这个习惯帮我省了至少十个小时的调试时间，希望帮到你。

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