做四旋翼控制仿真，PID调参调到吐是必经之路，但当你发现同样的参数换一组风扰就垮掉的时候，滑模控制就该上场了。这个仿真项目我前后复现了三轮，踩了不少坑，把最有价值的一条路线完整整理出来：基于Matlab/Simulink的双闭环滑模控制，位置外环加姿态内环，配合全局滑模面、指数趋近律和饱和函数抗抖振。这套架构做完之后，轨迹跟踪精度和抗扰能力比普通PID有质的提升，而且全局滑模消除了到达阶段，从初始时刻起就在滑模面上运动，鲁棒性覆盖全时间段。不管你是自动化专业做毕设，还是刚接触无人机控制想系统走一遍仿真，这篇文章都能给你一条可以直接跑通的路。

1. 双闭环滑模方案，为什么非这么拆不可

1.1 欠驱动系统与内外环分工的逻辑

四旋翼飞行器严格来说是一个欠驱动系统，四个电机输入，六个自由度输出，横滚、俯仰、偏航、三个位置，数学上无法做到对每个输出独立配一个输入。这是决定控制架构的底层约束。想要破解这个矛盾，思路就是分层：先把位置控制算出来，反解出需要的姿态角，再让姿态环去快速跟踪这个期望姿态。于是有了典型的内外环串级结构。

这套结构在工程上非常好理解。你想象开一辆车，你想要在地图上从A点到B点，这是位置层目标；具体怎么转方向盘、踩油门刹车，这是姿态层（身体的朝向与速度）要干的事。上层只负责规划"该往哪走"，下层负责"怎么执行到位"。外环控制器输出的是虚拟控制量，经过姿态解算变成期望横滚角和俯仰角；内环控制器跟踪这些期望角度，输出力矩控制量。

这样拆分之后，每个单环的模型规格都降下来了，控制器设计难度随之降低。位置外环面对的是一个近似的二阶积分器模型，姿态内环面对的也是一个带惯量矩阵的旋转动力学模型。两个环各自设计滑模面、各自用趋近律驱动，比设计一个全状态滑模控制器要直观得多，也更容易调试、更容易定位发散问题。

1.2 滑模控制与全局滑模的选型考量

滑模控制的核心思想是设计一个滑模面，让系统状态被"吸"到滑模面上，之后沿滑模面滑动到平衡点。这个滑动阶段对匹配扰动不敏感，所以滑模的鲁棒性来自这里——只要干扰满足匹配条件且幅值在界内，系统在滑模面上的行为基本不受干扰影响。

传统滑模有个绕不开的阶段划分：到达阶段和滑动阶段。系统从任意初始状态到滑模面之间那段路程，严格来说鲁棒性没有保证。全局滑模的思想就是把这部分也消掉，通过在滑模面中嵌入一个衰减项，让初始时刻系统状态就已经落在滑模面上。具体做法是给滑模面加一项关于初始误差的指数衰减函数，随着时间推移该项趋于零，滑模面退化成常规形式。好处是整个运动过程都处于滑模机制控制下，抗干扰能力全程生效，不信你试一下给位置环加一个常值风扰，全局滑模的轨迹偏移明显比普通滑模小。

另一个选型关键在趋近律。经典滑模常用符号函数构成等速趋近项，但符号函数的非连续性会带来高频抖振。我选的是指数趋近律，形式上是等速项加指数项的组合，指数项让系统远离滑模面时趋近速度更快，接近后等速项主导、切换幅度受控。配合饱和函数替换符号函数做准连续化处理，抖振问题能压到工程可接受的水平。这里要说明，抗抖振不是把抖振完全消除，而是用边界层内连续增益换取可接受的稳态误差，这个权衡后面第5节会详细讲。

2. 四旋翼动力学模型与仿真参数准备

2.1 简化动力学方程与状态定义

仿真模型不需要复杂到推导每片桨叶的气动力，做控制层验证时，把四旋翼当成刚体加上两个合力来源（总推力、三轴力矩）就够了。我用的是牛顿-欧拉形式，在地球系下写质心平动方程，在机体系下写转动方程。惯性系下位置方程是：

m * p_ddot = R * F_t + m * g * e3

其中R是机体坐标系到惯性系的旋转矩阵，F_t是机体z轴方向的总推力。转动方程在机体系下表示为：

J * omega_dot = -omega × (J * omega) + tau

J是惯量矩阵，omega是机体角速度，tau是控制力矩向量。这个模型忽略气动阻力、桨叶挥舞效应和陀螺力矩，对于控制算法对比来说足够了。想加入干扰，直接在位置方程或力矩方程后面加一个扰动项即可，方便测试鲁棒性。

状态变量取12个：p=[x,y,z]以及速度v=[vx,vy,vz]，姿态角phi=[φ,θ,ψ]以及角速度omega=[p,q,r]。Simulink里用积分器串联起来，输入是总推力U1和力矩向量[τφ,τθ,τψ]，输出12个状态。姿态角从角速度积分得到，旋转矩阵R用ZYX欧拉角形式计算。小角度假设下，位置方程可以线性化为x方向加速度约等于g乘以俯仰角，y方向约等于负的g乘以横滚角，这也是后面从虚拟控制量反解期望姿态角的数学基础。

2.2 模型参数的物理参考

仿真参数的合理性直接影响控制器的执行效果。我参考了一类常见小型四旋翼平台的物理参数：质量m取1.2kg，重力加速度g取9.8m/s²，惯量矩阵近似对角，Ixx=Iyy=0.0123kg·m²，Izz=0.0243kg·m²。推力系数和力矩系数这种与执行器相关的参数，在控制器直接输出力和力矩的仿真架构里不需要显式建模，如果后面要加电机动态，可以把控制量映射到PWM或转速，再通过系数转换成力和力矩。

期望轨迹我常用两类：一类是空间螺旋线，x和y随时间正弦变化，z线性上升，这样轨迹曲率变化平滑，容易观察跟踪效果；另一类是方波阶跃目标，用来检验系统的响应速度和超调。仿真时间一般设置20秒，位置初始在原点附近，姿态初始为零。为了让全局滑模的初始状态满足设计条件，控制器里需要读取当前实际状态来计算初始误差项，这个在Simulink里通过反馈信号实时计算即可。

3. 双闭环全局滑模控制器设计

3.1 位置外环滑模面与控制律推导

位置外环负责x、y、z三个通道。以x通道为例，定义跟踪误差e_x = x_d - x，其中x_d是期望轨迹。全局滑模面不直接用一阶形式，我会在滑模面中加入初始误差的衰减项：

s_x = e_x_dot + c_x * e_x - f_x(t)

其中f_x(t) = (e_x_dot(0) + c_x * e_x(0)) * exp(-a_x * t)，a_x是衰减速率。这样f_x(0)恰好等于括弧里的初始值，保证s_x(0)=0，系统从一开始就贴着滑模面运动。随着时间推移，f_x(t)趋于零，滑模面退化为常规形式。

对位置方程做小角度线性化，x通道的二阶模型近似为x_ddot = g * θ_d + d_x，其中θ_d是期望俯仰角，d_x是扰动。设计虚拟控制量u_x = x_ddot，令其满足滑模到达条件。取趋近律为指数形式：

s_dot_x = -ε_x * sat(s_x / Φ_x) - k_x * s_x

令v_x = x_ddot_d + c_x * e_x_dot + f_x_dot(t)，组合得到位置控制律：

u_x = v_x + ε_x * sat(s_x / Φ_x) + k_x * s_x

这里不要忘了，滑模面里的c_x和趋近律里的k_x、ε_x是三个作用完全不同的参数。c_x决定滑模面的收敛几何形态，k_x决定远离滑模面时的趋近速度，ε_x决定边界层内的切换强度。y通道完全同理，z通道则要考虑重力补偿，控制律里额外加上重力加速度项，形式为u_z = m * (g + v_z + ε_z * sat(s_z/Φ_z) + k_z * s_z) / (cosφ_d * cosθ_d)。

反解期望姿态角这一步最容易出错。虚拟控制量u_x、u_y是惯性系下的期望加速度，需要转换成姿态角。在小角度假设下，u_x ≈ g * θ_d，u_y ≈ -g * φ_d，所以θ_d = u_x / g，φ_d = -u_y / g。如果是大机动场景必须用完整旋转矩阵反解，但仿真里常规轨迹用小角度公式就够，简单且稳定。

3.2 姿态内环滑模面与力矩控制律

姿态内环的期望输入来自位置环解算：φ_d、θ_d，偏航角ψ_d通常直接给定常数，比如保持0度。姿态环的任务是让实际姿态角跟踪这些期望值，输出三轴力矩。

以横滚通道φ为例，定义姿态误差e_φ = φ_d - φ，全局滑模面为：

s_φ = e_φ_dot + c_φ * e_φ - f_φ(t)

f_φ(t)同样取初始误差指数衰减项。姿态动力学中φ通道的近似方程为φ_ddot = τ_φ / Ixx + d_φ，其中τ_φ是横滚力矩，d_φ是扰动。设计虚拟角加速度v_φ = φ_ddot_d + c_φ * e_φ_dot + f_φ_dot(t)，力矩控制律为：

τ_φ = Ixx * (v_φ + ε_φ * sat(s_φ / Φ_φ) + k_φ * s_φ)

俯仰通道τ_θ用Iyy同理，偏航通道τ_ψ用Izz同理。姿态环的带宽要明显高于位置环，这体现在c_φ的取值远大于c_x。我给一组常用起始值：位置环c_x=2，姿态环c_φ=12，这样内外环时间尺度分开，仿真稳定且轨迹平滑。总推力U1由z通道控制律直接给出，标量方程，不需要额外解算。

3.3 指数趋近律、饱和函数与参数关系

趋近律的设计直接决定控制信号品质。s_dot = -εsat(s/Φ) - ks可以拆开看：-ks是一个比例型趋近项，s大时该项贡献大，系统快速靠近滑模面，这就是"指数"的由来，不过实际是指数收敛特性，不是指数函数项；-εsat(s/Φ)是等速型趋近项，保证s进入边界层后仍然有切换力。两者组合的意义是：远快近稳，没有纯符号函数那种一路满幅切换的暴力感。

饱和函数取代符号函数是抗抖振的关键。sat(s/Φ)的表达式为：当|s|>Φ时等于sign(s)，当|s|≤Φ时等于s/Φ。数学上这是一个连续函数，边界层Φ越厚，控制信号越平滑，但边界层内滑模面不再精确为零，会带来稳态误差。我仿真里先用Φ=0.01跑，抖振明显；放大到Φ=0.05后，控制力矩曲线平滑很多，位置跟踪误差只增加了一点点。所以Φ不是越大越好，需要结合精度和抖动容忍度折中。

全局滑模的衰减速率a也值得留意。a太小，f(t)衰减太慢，滑模面长期像个"拖尾"，响应变慢；a太大，f(t)瞬间消失，全局滑模退化回普通滑模。我建议取a=2~5，保证过渡过程在1秒内完成，同时不影响稳态行为。

4. Matlab仿真搭建与代码实现

4.1 Simulink顶层模块划分

整个仿真模型可以分成六块：信号源（期望轨迹生成）、位置环控制器、姿态解算、姿态环控制器、四旋翼动力学模型、数据显示。我习惯用Simulink里的MATLAB Function模块写控制器逻辑，动力学模型单独用一个S函数或者另一个MATLAB Function，这样逻辑清晰，改参数也方便。

模块间的信号流向是：信号源输出12维期望状态，位置环读取期望[x,y,z]和实际[x,y,z]及速度，输出虚拟控制量[u_x,u_y,u_z]，姿态解算模块把u_x,u_y转成期望姿态角φ_d,θ_d，姿态环读取期望姿态角和实际姿态角及角速度，输出力矩[τφ,τθ,τψ]，动力学模型接收U1和力矩，积分输出12个状态，反馈回控制器。

模块清单我整理一下：

4.2 核心代码：饱和函数、位置环与姿态环

饱和函数是整个控制器里被调用最频繁的子函数，我单独写在MATLAB Function里：

function y = sat(x, phi) % 饱和函数，代替符号函数抑制抖振 if abs(x) <= phi y = x / phi; elseif x > phi y = 1; else y = -1; end end

位置外环控制器核心代码，输入期望位置xd、yd、zd和实际位置x、y、z及速度vx、vy、vz。全局滑模的初始误差项通过两个持续变量保存：

function [ux, uy, uz] = pos_control(xd, yd, zd, x, y, z, vx, vy, vz, t) persistent ex0 ey0 ez0 dex0 dey0 dez0 init if isempty(init) ex0 = xd - x; ey0 = yd - y; ez0 = zd - z; dex0 = 0; dey0 = 0; dez0 = 0; init = true; end % 位置环参数 cx = 2; cy = 2; cz = 3; kx = 3; ky = 3; kz = 5; ex = 0.5; ey = 0.5; ez = 0.8; phi_x = 0.05; phi_y = 0.05; phi_z = 0.05; a = 3; % 全局滑模面（带衰减项） fx = (dex0 + cx*ex0) * exp(-a*t); fy = (dey0 + cy*ey0) * exp(-a*t); fz = (dez0 + cz*ez0) * exp(-a*t); sx = (vx - (xd-x)*0 + (xd-x)*0 ... % 修正为正确的误差导数

这里要特别说明，误差导数的计算在仿真中不能直接对误差做数值微分，噪声会很大。我一般把期望轨迹设计成光滑可导的函数，直接用解析导数，或者借助状态变量里的速度反馈。期望速度已知，实际速度从状态直接反馈，所以e_dot = vx_d - vx，代码里直接用这个表达式。

姿态内环控制器接收期望角φd、θd、ψd和实际角度φ、θ、ψ及角速度p、q、r：

function [tau_phi, tau_theta, tau_psi] = att_control(phid, thetad, psid, phi, theta, psi, p, q, r, t) persistent ephi0 etheta0 epsi0 ... cphi = 12; ctheta = 12; cpsi = 10; kphi = 8; ktheta = 8; kpsi = 6; ephi = 1.0; etheta = 1.0; epsi = 0.8; phi_b = 0.05; theta_b = 0.05; psi_b = 0.05; Ixx = 0.0123; Iyy = 0.0123; Izz = 0.0243; % 全局滑模面 fphi = (d_ephi0 + cphi*ephi0) * exp(-3*t); sphi = d_ephi + cphi*ephi - fphi; % 控制力矩 tau_phi = Ixx * (dd_phid + cphi*d_ephi + dfphi + ephi*sat(sphi/phi_b) + kphi*sphi); ... end

实际写代码时我需要提醒你，MATLAB Function里的persistent变量在每次仿真开始时要确保清零，否则多次运行仿真的初始误差状态会残留，造成结果不一致。我通常加一个手动复位开关，或者用Simulink的Reset端口，调试时特别管用。

4.3 仿真参数配置与轨迹设置

仿真求解器我选ODE45变步长，相对误差设为1e-6，最大步长限制在0.01秒。固定步长也能跑，但滑模控制中的切换动作对求解精度敏感，变步长能更好捕捉边界层内的动态。如果后期要部署到实时硬件，再改用固定步长0.001秒并加离散化处理，仿真阶段没这个必要。

期望轨迹生成函数示意：

function [xd, yd, zd, psid] = trajectory(t) % 螺旋上升轨迹 xd = 2 * sin(0.2 * t); yd = 2 * cos(0.2 * t); zd = 0.5 + 0.1 * t; psid = 0; end

跑完之后我一般把状态量导到工作区，画三张图：位置跟踪对比图、姿态角响应图、控制输入（力矩）曲线图。位置对比图能直观看出跟踪误差；控制输入曲线最容易暴露抖振问题——如果力矩曲线像锯齿一样密密麻麻高频翻转，就是抖振没压住；如果曲线平滑但位置误差偏大，就要检查是边界层太厚还是趋近律增益不够。

5. 调试记录与避坑清单

5.1 高频抖振问题

抖振是滑模控制第一次跑仿真最常撞上的问题。现象描述：控制力矩输出呈现高频振荡，频率可以到几百赫兹级别，轨迹倒是还能跟踪，但控制信号物理上不可执行，仿真步长稍大甚至会出现状态发散。仿真里看起来还能"跑"，实际上这个结果离硬件应用差十万八千里。

抖振的直接来源就是符号函数sign(s)在滑模面两侧的来回切换。如果滑模面附近有噪声或数值误差，切换信号会变得极其频繁。我处理抖振的优先级依次是：先用饱和函数替换符号函数，边界层厚度从0.01开始试；如果还抖，检查是不是趋近律的ε取太大，ε是等速趋近项的幅值，相当于单位阶跃切换力，ε越大越容易抖；最后再检查仿真步长，变步长最大步长收紧到0.001能明显改善。

这里有个误区要避开：很多人为了消抖把边界层调得特别厚，比如Φ=0.5，抖振确实没了，但位置稳态误差也上去了。滑模控制的精髓在鲁棒性，边界层把鲁棒性削弱了，就失去意义了。我个人的经验是，位置环Φ在0.03~0.08之间，姿态环Φ在0.02~0.05之间，是一个比较合理的区间。

5.2 发散与超调排查

仿真发散是调试期的家常便饭。我复现这个项目时遇到过两次典型的发散，根源完全不同。

第一次是姿态解算符号反了。位置环解算出θ_d = u_x/g没问题，但在姿态内环里误差定义写成e_θ = θ - θ_d，控制律方向全反了，形成正反馈，姿态角几秒内就冲到几十度然后发散。这种问题排查最快的方法是看姿态响应曲线：如果姿态角往期望角反方向冲，基本就是误差符号或者控制增益符号错了。

第二次是内外环带宽不匹配。我把位置环的c_x调成和姿态环的c_φ差不多大，结果位置环快速输出大幅度姿态期望，姿态环跟不上，整个系统进入竞争状态，轨迹出现明显的蛇形振荡。后来把姿态环带宽保持在外环的4~6倍，这个振荡立刻消失。整定原则：先把内环调稳，再把外环参数慢慢加上去，不要一上来就追求外环的收敛速度。

单位问题也值得说。我一开始把期望轨迹用角度制写，控制律全是弧度制，位置环直接炸。检查办法是在控制器入口打印一次信号值，看看数量级是否合理，横滚角应该是0.1弧度级别而不是几度级别。

5.3 参数整定的实操顺序

滑模控制器参数看着多，其实整定顺序有规律可循。我建议按以下步骤操作：

第一步，先整定内环。把姿态环独立出来，期望姿态角给阶跃信号，调c_φ和k_φ让姿态响应无超调、快速稳定。c_φ决定滑模面的收敛速度，k_φ决定趋近速度，这两个参数优先从偏小的值往上加。

第二步，整定外环。接上位置环和姿态解算，期望轨迹先给慢速直线，把c_x、c_y、c_z从1~2开始调，让位置误差收敛。外环的输出是期望姿态角，如果期望姿态角频繁突变，说明c_x偏大，要降。

第三步，加趋近律增益。ε和k一起调，观察控制力矩的幅值和抖振。k偏大控制响应快但易振荡，ε偏大容易抖振，两者搭配时优先保证k适中、ε偏小。

第四步，加干扰验证。在位置模型里加一个幅值1~2的常值扰动，看系统能否保持跟踪精度。这一步是全局滑模和普通滑模差异最明显的地方，你对比一下就能直接看出全局滑模在初始阶段几乎不产生位置偏移。

这个项目做完之后，你手里就有了一套可以反复扩展的基础平台。后面想加自适应律估参数、加干扰观测器补偿匹配扰动、甚至切换成四元数姿态表示做全姿态机动，都是在这个框架上做增量开发。我用这套仿真对比过理想PID、普通滑模和全局滑模三种方案，在同样风扰条件下的位置误差峰值差距非常直观，谁抗扰强、谁切换平滑，一条曲线就说明问题。这也是我建议你把原理吃透而不是只跑通Demo的原因——控制方案的取舍判断力，比跑通一个仿真值钱得多。