4. 姿态动力学与轨道动力学
姿轨联合控制系统的“被控对象”由两部分动力学共同描述：姿态动力学决定航天器绕质心的转动，轨道动力学决定质心在惯性空间中的平动。二者通过推力矢量方向、质心位置变化、太阳翼转动、燃料消耗等途径耦合。本节先建立刚体姿态动力学方程与动量矩、陀螺效应的物理图像，再给出轨道动力学方程与推力、摄动力的建模方法，为后续控制器设计提供模型基础。
一、刚体姿态动力学方程
设航天器本体坐标系为 Obxbybzb，原点取在质心，三轴与惯量主轴重合。刚体绕质心转动的角速度在本体系下记为 ω = [ωx, ωy, ωz]T，惯量矩阵为 J = diag(Jx, Jy, Jz)。由角动量定理，在惯性系中角动量 H 的时间导数等于外力矩 M；转换到本体系后得到欧拉方程：
Jx·ω̇x - (Jy - Jz)·ωy·ωz = Mx
Jy·ω̇y - (Jz - Jx)·ωz·ωx = My
Jz·ω̇z - (Jx - Jy)·ωx·ωy = Mz式中 M = [Mx, My, Mz]T 为本体系下的合外力矩，包括控制力矩、气动力矩、重力梯度力矩、太阳光压力矩和磁力矩等。方程左端的 (Jy−Jz)ωyωz 等项称为陀螺力矩项，它体现了角速度分量之间的交叉耦合：即使某一轴没有外力矩，只要另外两轴存在角速度，该轴也会产生角加速度。这是姿态动力学区别于一般线性系统的关键非线性特征。
姿态运动学用于描述姿态参数的变化。若采用四元数 q = [q0, q1, q2, q3]T，则运动学方程为：
q̇0 = -0.5·(q1·ωx + q2·ωy + q3·ωz)
q̇1 = 0.5·(q0·ωx - q3·ωy + q2·ωz)
q̇2 = 0.5·(q3·ωx + q0·ωy - q1·ωz)
q̇3 = 0.5·(-q2·ωx + q1·ωy + q0·ωz)动力学方程给出 ω 的演化，运动学方程把 ω 积分成姿态参数，二者串联构成完整的姿态模型。工程实现中通常以四元数积分、再归一化，避免欧拉角在 90° 附近的奇异。
二、动量矩与陀螺效应
刚体角动量 H = Jω。当星上装有高速旋转部件（如动量轮、控制力矩陀螺）时，系统总角动量为 H = Jω + hw，其中 hw 为轮系相对本体的角动量。此时欧拉方程需改写为：
J·ω̇ + ω × (J·ω + h_w) = M - ḣ_w其中 ω × hw 就是陀螺力矩。它的物理含义是：当本体以角速度 ω 转动时，旋转的轮子会产生一个垂直于 ω 与 hw 的力矩，迫使航天器绕第三轴转动。这正是控制力矩陀螺（CMG）的工作原理——通过改变框架角使 hw 方向变化，从而输出控制力矩；同时也是动量轮系统产生轴间耦合、需要进行解耦控制的原因。
三、轨道动力学方程
把航天器视为质点，在惯性系（如 J2000）中位置为 r，速度为 v，则轨道动力学方程为：
ṙ = v
v̇ = -μ·r/|r|³ + a_pert + a_thrust其中 μ = GMe 为地球引力常数（取 3.986×105 km³/s²），第一项为二体中心引力加速度，apert 为摄动加速度，athrust 为推力加速度。若用轨道六根数 (a, e, i, Ω, ω, M) 描述，则二体运动下除平近点角 M 外其余根数均为常数，M 随时间线性变化。摄动力和推力则使各根数缓慢漂移，这正是轨道保持与机动要处理的对象。
四、推力与摄动力建模
推力建模的关键是方向与作用点。设推力器在本体系下的安装方向为 ut，大小为 T，则惯性系推力加速度为：
a_thrust = (T/m) · R_bi · u_t其中 Rbi 为本体系到惯性系的姿态旋转矩阵，m 为质量。可见推力对轨道的作用依赖当前姿态，姿态误差会直接造成推力方向偏差；同时推力作用点若偏离质心，还会产生干扰力矩 Mt = rt × (RibutT)，反过来影响姿态。这就是姿轨耦合最直接的通道。
摄动力建模通常包括以下几项：
| 摄动源 | 典型量级（LEO） | 建模要点 |
|---|---|---|
| 地球非球形 J2 项 | 约 10-3 m/s² | 引起 Ω、ω 长期漂移，是轨道保持主因 |
| 大气阻力 | 10-6～10-5 m/s² | 与面质比、大气密度模型相关，降低半长轴 |
| 太阳光压 | 10-7 m/s² | 与受照面积、反射系数有关，高轨更显著 |
| 三体引力 | 随轨道高度变化 | 高轨需考虑日月引力 |
工程上常把摄动分为保守力与非保守力：J2、三体引力属保守力，可用势函数梯度表示；大气阻力、光压属非保守力，需按面积和姿态实时计算。对于姿轨联合控制，摄动力不仅改变轨道，还会通过重力梯度、气动压心偏移产生姿态干扰力矩，因此建模时应在同一仿真框架内统一处理。
五、建模实现要点
把上述方程落到仿真或星上算法时，需注意以下几点。第一，坐标系定义必须全文一致：本体系取质心为原点、三轴为惯量主轴；惯性系取 J2000；旋转矩阵 Rbi 表示本体系到惯性系的转换，方向不能混用。第二，积分步长要与动力学带宽匹配：姿态角速度可达 0.1～5 °/s，轨道周期约 90 min，姿轨联合仿真常取 0.01～0.1 s 步长，用四阶龙格库塔积分，四元数每步归一化。第三，惯量矩阵随太阳翼转动和燃料消耗变化，可在参数表中按工况给出 Jx、Jy、Jz 的标称值与变化范围。第四，摄动力按轨道高度选择：LEO 重点建 J2 与大气阻力，GEO 重点建光压与日月三体引力。
下面给出一组教学示例参数，用于验证模型实现是否正确（数值为合理假设，非实测）：Jx=1200、Jy=1500、Jz=1000 kg·m²，初始 ω=[0.02, −0.01, 0.03]T rad/s，无外力矩。按欧拉方程积分 100 s，可观察到 ωx、ωy、ωz 出现周期性交换，角动量模长 |Jω| 保持不变——这是无外力矩刚体的自由运动特征，可作为代码正确性的第一道检查。若 |Jω| 随时间明显漂移，说明陀螺项符号或积分器有误。
轨道部分可用同一组参数检查：取 a=7000 km、e=0.001、i=98°，仅保留 J2 项积分 1 天，应看到 Ω 与 ω 的线性漂移，而 a、e、i 基本不变。若 a 出现明显下降，说明阻力项被误加；若 i 快速变化，说明 J2 项公式中的 sin i 因子写错。通过这些守恒量与漂移趋势的对照，可以在进入控制器设计前确认动力学模型可信。