简介：面向航天器姿态控制的研究者和工程师，这份资料围绕四元数与自适应滑模控制的结合展开。四元数可避免万向锁并稳定描述旋转，自适应滑模则在线调整参数以应对系统不确定性和外界干扰，二者互补用于高精度姿态跟踪。包内共7个文件，以MATLAB脚本（.m）、Simulink模型（.mdl）和图形（.fig）为主，其中SMC.mdl还提供了r2012a与r2014a两个版本，便于不同环境打开运行；plant.m、contr.m等脚本对应对象建模与控制律实现，37KB的小体积适合快速查阅。已有476人学习，特别适合课程设计、科研预研和工程入门。通过研究模型结构、脚本代码与仿真图形，可掌握基于四元数的自适应滑模控制器设计流程，包括滑模面选取、参数自适应律构造及稳定性分析，并学会处理陀螺效应、重力梯度等干扰因素，为实际航天器控制方案提供参考实现。

1. 为什么航天器姿态控制绕不开四元数与自适应滑模

把一颗对地观测卫星从初始姿态转到目标姿态，这看似简单，但燃料消耗导致惯量变化、太阳翼展开后的结构扰动、重力梯度力矩和光压带来的外干扰，都会让固定增益控制器在某个阶段突然“失灵”。用四元数做姿态表示，配合自适应滑模控制，能同时处理姿态表示的全局无奇异问题和对不确定性的在线补偿，这也是这个标题背后最常见的工程动机。

适合读这篇文章的人，是正在搭建航天器姿态仿真、被欧拉角奇异或滑模抖振折腾过一遍的控制方向工程师和研究生。读者手里不一定有现成代码包，但只要能复现模型和控制器，就能把闭环跑通。这篇文章按“原理 → 设计 → 仿真 → 排查 → 进阶验证”的顺序展开，每一步都给可运行的细节和参数。

2. 四元数姿态模型与自适应滑模：为什么误差四元数是控制器的眼睛

2.1 四元数为什么比欧拉角靠谱：从公式看奇异性

航天器姿态控制里，表示姿态的方式不外乎欧拉角、方向余弦矩阵和四元数。欧拉角直观，3-2-1转序在俯仰角接近90度时会出现万向锁，对应的姿态运动学方程分母趋近零，数值上直接爆掉。四元数用四个分量表示一次旋转，归一化后没有奇异性，这也是深空任务普遍采用四元数的原因。

四元数 (\mathbf{q} = [q_0, q_v^T]^T)，其中 (q_0) 是标量部分，(q_v) 是矢量部分，满足 (q_0^2 + q_v^T q_v = 1)。姿态运动学方程写成：

[ \dot{\mathbf{q}} = \frac{1}{2} \begin{bmatrix} -q_v^T \ q_0 I_3 + q_v^\times \end{bmatrix} \boldsymbol{\omega} ]

其中 (\boldsymbol{\omega}) 是本体系相对惯性系的角速度，(q_v^\times) 是叉乘反对称矩阵。这个方程的好处是全局连续，不需要处理欧拉角跳变。实际编程时，我不建议手推这个矩阵的每个元素，直接写一个四元数微分函数更不容易出错。

误差四元数定义为 (\mathbf{q}_e = \mathbf{q}_d^{-1} \otimes \mathbf{q})，其中 (\mathbf{q}_d) 是目标姿态。误差四元数的矢量部分反应了当前姿态与目标姿态的旋转轴和旋转角，直接可以作为滑模面的状态量，这是后面控制器设计的基础。

2.2 姿态动力学方程与不确定性来源

刚体航天器的姿态动力学方程为：

[ \mathbf{J} \dot{\boldsymbol{\omega}} + \boldsymbol{\omega} \times \mathbf{J} \boldsymbol{\omega} = \mathbf{u} + \mathbf{d} ]

(\mathbf{J}) 是惯量矩阵，(\mathbf{u}) 是控制力矩，(\mathbf{d}) 是外干扰力矩。工程里没有谁真能用标称惯量做出完美控制器。推进剂消耗会让惯量在任务周期内变化百分之几到百分之十几，太阳翼的振动模态虽然不直接进刚体模型，但等效为高频扰动。这些在模型里合并成惯量偏差 ( \Delta \mathbf{J} = \mathbf{J} - \mathbf{J}_0 ) 和扰动项 (\mathbf{d})。

把惯量偏差和扰动打包，定义集总扰动：

[ \mathbf{L} = -\mathbf{J}_0^{-1}(\Delta \mathbf{J} \dot{\boldsymbol{\omega}} + \boldsymbol{\omega} \times \Delta \mathbf{J} \boldsymbol{\omega} - \mathbf{d}) ]

这个集总扰动有界，但界是多少，设计阶段说不准。自适应滑模的核心思路就是：不假设一个固定的扰动界，而是用一个自适应律在线估计出这个界，再把它用在切换项的增益里。

2.3 自适应滑模控制的设计骨架

滑模面取为：

[ \mathbf{s} = \boldsymbol{\omega}e + \lambda \mathbf{q}{ev} ]

其中 (\boldsymbol{\omega}_e = \boldsymbol{\omega} - \boldsymbol{\omega}_d)，姿态机动任务里 (\boldsymbol{\omega}d = 0)。(\lambda > 0) 是滑模面系数，它决定姿态误差和角速度误差在滑模面里的权重。(\mathbf{q}{ev}) 是误差四元数的矢量部分。

有了滑模面，控制力矩写成三项叠加：第一项是科氏力矩前馈 (\boldsymbol{\omega} \times \mathbf{J}_0 \boldsymbol{\omega})，用来抵消标称模型里的非线性项；第二项是线性趋近项 (-\mathbf{K} \mathbf{s})，驱动滑模面向零点收敛；第三项是自适应切换项 (-\hat{\rho} \tanh(\mathbf{s}/\epsilon))，用来压制集总扰动。

自适应律设计成：

[ \dot{\hat{\rho}} = \gamma (|\mathbf{s}| - \sigma \hat{\rho}) ]

(\gamma) 是自适应增益，(\sigma) 是泄漏系数。很多人对自适应律感到像黑匣子，其实它只是在做一个估计动作：滑模面远离零点时，(|\mathbf{s}|) 大，(\hat{\rho}) 增大，切换力度加强；滑模面接近零点时，泄漏项防止估计值无限涨上去。

3. 自适应滑模控制器设计：从滑模面到自适应律的完整递推

3.1 把控制器公式逐行落定

这一节不给理论推导，只给可以直接落到代码的控制器骨架。控制目标是对每个仿真步根据当前四元数、角速度和目标姿态，算出一组控制力矩。控制器内部顺序是：

第一步，根据当前姿态与目标姿态计算误差四元数；第二步，用误差四元数的矢量部分和角速度差拼出滑模面；第三步，由滑模面计算线性反馈项和切换项；第四步，加上科氏力矩前馈，得到最终控制力矩。

误差四元数计算注意顺序：(\mathbf{q}_e = \mathbf{q}_d^{-1} \otimes \mathbf{q})，不是反过来的 (\mathbf{q}_d \otimes \mathbf{q}^{-1})。方向错了，控制器会往反方向推，姿态越控越偏。

这里有一个关键的工程修正：( \mathbf{q} ) 与 ( -\mathbf{q} ) 表示同一个姿态，数值上如果四元数随时取反，误差四元数会在两个符号版本之间跳变，滑模面的符号也跳，控制力矩来回翻转。我一般在每个仿真步固定四元数标量部分为正，即 ( q_0 < 0 ) 时整体取负，这个操作很小但能避免大量怪问题。

切换项不用 (\text{sign}(\mathbf{s}))，用 (\tanh(\mathbf{s}/\epsilon))。符号函数在零点附近理想切换，数值仿真里就是高频抖振的来源。双曲正切函数在边界层 (\epsilon) 内是连续过渡，能同时抑制抖振，代价是加大一点点稳态误差。

3.2 用 MATLAB 函数封装控制器，接口对准仿真循环

MATLAB 是航天器姿态仿真最常用的环境，控制器写成独立函数比较干净，方便替换成 PID 或固定增益滑模做对比。下面这段是控制器核心函数：

function u = ads_smc_controller(q, omega, qd, omega_d, J0, K, lambda, epsilon, rho_hat) % 自适应滑模控制器 % 输入: % q - 当前姿态四元数 [q0; q1; q2; q3]，单位四元数 % omega - 当前角速度 [wx; wy; wz]，rad/s % qd - 目标姿态四元数 % omega_d - 目标角速度，姿态机动通常为 [0;0;0] % J0 - 惯量矩阵标称值，3x3 % K - 线性反馈增益，3x3 对角阵 % lambda - 滑模面系数标量，通常 0.5~2 % epsilon - 边界层厚度，通常 0.01~0.1 % rho_hat - 自适应扰动上界估计值，标量 % 输出: % u - 控制力矩 [Tx; Ty; Tz]，N*m % 误差四元数: qe = qd_inv ⊗ q qe = quat_multiply(quat_inverse(qd), q); qe = quat_normalize(qe); % 角速度误差 omega_e = omega - omega_d; % 滑模面 s = omega_e + lambda * qe_v s = omega_e + lambda * qe(2:4); % 科氏力矩前馈项 u_cross = cross(omega, J0 * omega); % 线性反馈 + 自适应切换项 u_feedback = -K * s; u_switch = -rho_hat * tanh(s / epsilon); % 最终控制力矩 u = u_cross + u_feedback + u_switch; end

这个函数里依赖两个四元数操作函数：quat_multiply 计算四元数乘法，quat_inverse 计算逆四元数。四元数乘法顺序很重要，注释里已经标了。增益 K 我习惯取对角阵，三个轴的反馈力度独立调，后面的参数表给一组能直接收敛的初值。

3.3 自适应律的更新代码与参数初值

自适应律 (\dot{\hat{\rho}} = \gamma(\lVert\mathbf{s}\rVert - \sigma\hat{\rho})) 本身是微分方程，仿真里用欧拉法离散即可。注意 (\hat{\rho}) 是标量，在控制器函数外面更新，否则每个仿真步看不到滑模面的完整轨迹。

function rho_hat = update_rho(rho_hat, s, gamma, sigma, Ts) % 自适应律离散更新 % s - 当前滑模面 3x1 % gamma - 自适应增益，控制估计速度 % sigma - 泄漏系数，防止估计值发散 % Ts - 仿真步长，秒 rho_hat = rho_hat + (gamma * norm(s) - sigma * rho_hat) * Ts; rho_hat = max(rho_hat, 0); % 扰动上界不可能是负数 end

gamma 取 0.5 到 5，越小估计越稳但收敛越慢；sigma 取 0.01 到 0.1，主要作用是给估计值加阻尼。(\hat{\rho}) 初始值取 0 就行，自适应会让它自己长起来。如果仿真后期发现力矩饱和，检查 (\hat{\rho}) 是否涨到太大，优先调大 sigma，而不是调小 gamma。

这一节的公式和代码合在一起，就是一个完整的自适应滑模控制器闭环框架。很多做固定增益滑模的工程抱怨抖振无法接受，其实问题不在方法，而在切换项和自适应参数没有按这种方式耦合调参。下面一节把整个闭环仿真搭起来。

4. 用 MATLAB 搭闭环仿真：四元数更新脚本与参数表

4.1 姿态运动学和动力学更新：欧拉法与四元数归一化

闭环仿真需要一个主循环，每个步长内完成四步：控制器算力矩、动力学更新角速度、运动学更新四元数、自适应律更新扰动界。数值积分我用欧拉法演示，步长取 0.01 秒以下，精度足够用于验证控制算法。步长放大到 0.05 秒时，四元数更新容易漂移，需要在每步末尾强制归一化。

四元数运动学更新的核心是构造前面第 2 章的矩阵。代码直接放在主循环里：

for k = 1:N_steps % 第1步: 控制器输出力矩 u = ads_smc_controller(q, omega, qd, omega_d, J0, K, lambda, epsilon, rho_hat); % 第2步: 动力学更新，J0只用于标称模型，实际被控对象带惯量偏差和扰动 L = -J0_inv * (cross(omega, J0*omega) - u - disturbance(t(k))); domega = L; omega = omega + domega * Ts; % 第3步: 运动学更新，四元数微分 dq = 0.5 * Omega(omega) * q Omega = [0, -omega(1), -omega(2), -omega(3); omega(1), 0, omega(3), -omega(2); omega(2), -omega(3), 0, omega(1); omega(3), omega(2), -omega(1), 0]; q = q + 0.5 * Omega * q * Ts; q = q / norm(q); % 强制归一化 % 第4步: 自适应律更新 qe = quat_multiply(quat_inverse(qd), q); s = omega - omega_d + lambda * qe(2:4); % 注意与控制器内部保持一致 rho_hat = update_rho(rho_hat, s, gamma, sigma, Ts); end

动力学更新里我故意用了 (-\mathbf{J}_0^{-1} (\dots)) 并让实际被控对象叠加扰动，是为了让仿真更贴近真实。实际任务中控制器用的惯量矩阵是地面装订值，卫星上天后惯量已经变了，所以仿真里在动力学方程中混入惯量偏差是检验自适应滑模能力的重要手段。

4.2 参数表：一组能直接收敛的初值

(\mathbf{K}) 与 (\lambda) 有关系但不算耦合：(\lambda) 决定滑模面里姿态误差的分量，(\mathbf{K}) 决定滑模面整体趋近零的速度。先把 (\lambda) 固定为 1，(\mathbf{K}) 从 3 加到 10，观察姿态误差曲线。如果你发现角速度先剧烈震荡后收敛，多半是 (\lambda) 太大，先降回 0.5 再重新调 (\mathbf{K})。

4.3 结果怎么判：收敛时间、稳态精度和抖振幅度

一套闭环跑完，主要看三个指标。第一个是收敛时间，从初始姿态到姿态误差小于 0.1 度的时间，机动任务通常要求在几十秒内完成；第二个是稳态精度，最终姿态误差和角速度误差是否落在合理范围；第三个是控制力矩的抖振幅度，切换项没有边界层时力矩曲线像锯齿，加了 (\tanh) 后应当平滑很多。

我一般会把三个图放在一起看：姿态误差四元数矢量部分的曲线、角速度曲线、控制力矩曲线。姿态误差曲线如果有高频分量，先怀疑四元数符号问题；角速度曲线如果在零点附近来回穿越，检查边界层厚度；力矩曲线如果有明显毛刺，基本是自适应增益 (\gamma) 太大。

5. 自适应滑模仿真常见问题排查：现象、原因与解决

5.1 姿态不收敛，角速度来回振荡

现象：姿态误差长时间不归零，角速度曲线像“开关”一样来回摆动，控制力矩正负交替，仿佛控制器在跟被控对象斗气。

原因：最常见的是四元数符号跳变。四元数 (\mathbf{q}) 和 (-\mathbf{q}) 表示同一姿态，但误差四元数会在两个版本之间跳，滑模面符号跟着跳，控制器方向就乱了。另一个常见原因是 (\lambda) 取值过大，滑模面里姿态误差项占了主导，导致过阻尼一样的振荡。

解决：每个仿真步末尾加上符号修正：if q(1) < 0, q = -q; end，保证标量部分为正。(\lambda) 先取 0.5 试跑，确认无误再往上调。这个坑我见到的概率非常高，很多情况下不是控制器设计的问题，而是四元数符号管理不到位。

5.2 控制力矩高频抖振

现象：控制力矩曲线出现密集的锯齿状毛刺，幅值不大但频率很高，放到幅频特性里能看到高频分量。更严重的情况下，抖振会激励未建模的柔性模态，太阳翼之类的大柔性结构随之振动。

原因：切换项用了符号函数或者等效的硬切换。sign 函数在滑模面过零时产生理想切换，这在连续系统里没问题，但在离散仿真里就是高频抖振源。边界层厚度 (\epsilon) 取太小，也会让 tanh 逼近 sign 的行为。

解决：切换项全部改为 (\tanh(\mathbf{s}/\epsilon))，(\epsilon) 从 0.01 到 0.1 试。边界层加厚，抖振减小，但稳态误差变大；边界层变薄，稳态精度好但抖振加重。这是一个显式的权衡，没有一劳永逸的参数，我一般先取 0.05 跑通闭环再微调。

5.3 自适应参数 (\hat{\rho}) 越涨越大

现象：仿真后半段，(\hat{\rho}) 从个位数涨到几十上百，控制力矩持续增大直到进入饱和限幅，姿态反而失控。

原因：自适应律只有正增益没有泄漏项，或者泄漏系数 (\sigma) 太小。滑模面只要有一点偏差，(\hat{\rho}) 就持续上升，而扰动实际没有那么大的上界，估计值过度放大。

解决：加入泄漏项 (\sigma\hat{\rho})，并给 (\hat{\rho}) 加硬限幅。泄漏系数 0.01 到 0.1 是合理范围，硬限幅的值根据执行机构力矩上限决定，比如反应轮单轴最大力矩 0.2 N·m，(\hat{\rho}) 限幅到 5 就足够。加限幅后自适应收敛变慢，但系统稳定性明显提升。

5.4 单位混用导致结果差几个数量级

现象：同样的控制器参数，换了一台电脑或一份数据，仿真结果完全不对，控制力矩要么小到没作用，要么大到直接饱和。

原因：航天器模型里惯量的单位是公斤·平方米，角速度是弧度每秒，力矩是牛顿·米。工程脚本里常见的错误是把角速度用了度每秒，或者惯量用了克·平方厘米。四元数运动学方程里的角速度必须是弧度制，否则微分方程本身的系数就不对。

解决：在仿真脚本头部单独写一个单位常量区，把所有量纲统一声明。角速度在输入控制器前核对单位，是弧度每秒而不是度每秒。这个坑属于最没技术含量但翻车率最高的类型。

5.5 自适应增益 (\gamma) 调大后引入了低频振荡

现象：(\gamma) 从 1 调到 5，扰动估计收敛快了，但姿态误差曲线出现明显的周期性起伏，频率在 0.1 Hz 量级。

原因：自适应回路的时间尺度与姿态控制回路耦合。(\gamma) 太大时，扰动上界估计本身变成一个振荡源，它估计的不是真实扰动上界，而是去追滑模面的高频误差，导致控制器输出里混入低频振荡分量。

解决：(\gamma) 不要一开始就调大，先用固定增益滑模跑通闭环，记录滑模面残差水平，再逐步加大 (\gamma)。理想状态是 (\hat{\rho}) 曲线平滑上升并稳定，如果有明显波纹，减小 (\gamma) 或者增大 (\sigma)。这一条是经验性质的调参规律，数学模型也能解释，但仿真里看曲线比算特征值来得快。

6. 验证自适应能力的三个实验：惯量突变、初值撒点与周期扰动

控制器闭环跑通只是第一步，真正判断自适应滑模值不值得用，要看它面对不确定性时的表现。我建议动手加三个实验，难度递增但每个改动都很小。

第一个实验是惯量突变。仿真进行到 20 秒时，把实际惯量矩阵换成加了 15% 偏差的版本，模拟推进剂消耗或太阳翼展开后的惯量变化。如果控制器在突变后几秒内重新收敛，说明集总扰动上界的自适应估计确实在工作；如果姿态误差明显增大但不恢复，检查自适应增益 (\gamma) 是否足够大。

第二个实验是初值撒点。随机生成 30 组初始姿态，覆盖从 0 度到 150 度范围内的不同旋转轴和旋转角，跑批量仿真。统计每组收敛时间和稳态误差，画出散点分布。这个实验能暴露控制器在某个特定姿态附近的不稳定行为，比如滑模面系数 (\lambda) 在小角度机动时表现良好、在大角度机动时振荡的问题。

第三个实验是周期扰动。在动力学方程里叠加一个幅值 0.01 N·m、频率 0.05 Hz 的正弦扰动，模拟轨道周期性外扰。观察 (\hat{\rho}) 是否跟着扰动周期波动，以及切换项是否随之调整力度。自适应滑模的核心优势恰恰在这里：固定增益滑模必须用一个保守的大切换增益去覆盖扰动，而自适应滑模只需要估计出实际扰动上界，省下的力矩余量可以用在更大的姿态机动上。

有一次我做惯量突变实验时，加了泄漏项的 (\hat{\rho}) 在突变后无法快速升到新扰动上界，系统在过渡期出现了明显超调，最后靠增大 (\gamma) 和降低 (\sigma) 组合调整才压下来。这件事让我养成了一个习惯：每次改完参数，先跑惯量突变场景看过渡过程，再回去调稳态指标。顺序反了，往往会被稳态的良好表现误导。自适应滑模不是一个可以永远信任固定参数的控制器，它的参数本身就是自适应回路的一部分，希望这个方法能让你的姿态仿真实验少走一段弯路，希望帮到你。

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