企业官网建设流程全解析

简介：本资源是一套面向自动化、电子信息与应用数学等专业高年级本科生及飞行控制初学者的MATLAB教学仿真平台，聚焦小型旋翼无人机六自由度非线性动力学建模与动态响应分析。资源提供R2014a、R2019b、R2024b三版本兼容环境，涵盖113个文件，以34个参数配置.mat文件、22个核心算法.m脚本（如mavDynamics.m、FandM.m）、8个Simulink仿真模型.slx及17个.xml工程配置文件为主，总大小仅368KB，轻量易部署。已有75人学习下载，适用于控制系统设计、飞行器建模课程实验及运动规划算法验证场景。用户可直接运行runfg_small_UAV.bat启动仿真，通过修改质量特性、惯性张量与气动系数等配置参数，快速开展构型对比与动态特性扫描分析；所有关键模块均含中文注释，模型严格基于牛顿-欧拉方程构建，融合气动、推进与重力效应，具备工程级建模精度与教学级可读性。

基于Matlab的小型无人机6自由度非线性建模与仿真分析

以前总觉得无人机建模仿真是件"加分项"的事，真机飞得好就行。直到前段时间帮朋友调一台自制四旋翼，姿态控制参数在悬停时明明很稳，一做大角度翻转机动姿态就直接发散，光靠真机试飞反复试错烧掉好几组桨，我才意识到一个完整的6自由度非线性模型在Matlab里有多重要。后来我把整个飞行器的运动方程原原本本搬进Simulink，用非线性模型做故障复现和控制参数预调，问题才真正定位清楚。这篇东西就把我这次基于Matlab的小型无人机6自由度非线性建模与仿真分析过程完整梳理一遍，从坐标系定义讲到方程推导，再到Simulink实现和模型校核，适合正在做飞控算法验证、无人机课程设计或者想入门飞行器建模的读者参考。

1. 非线性建模为什么绕不开：从一次悬停实验说起

1.1 线性化模型的适用边界到底在哪

大多数飞控教材在讲控制器设计时，都会先做一步"小扰动线性化"：假设无人机悬停在某个平衡点附近，姿态角和角速度变化都很小，于是把含有三角函数、速度叉乘项的非线性方程直接用泰勒展开近似成线性方程。这个思路在理论分析阶段完全没问题，线性化模型最大的价值在于能直接用特征值、根轨迹、频域响应这些工具来分析系统的稳定性和控制器参数，设计思路清晰，运算成本低。

但问题藏在"小扰动"这三个字里面。所谓的"小"究竟是多小？从我的工程经验看，当滚转角或俯仰角的偏差超过15度、角速度超过0.5 rad/s、线速度超过5 m/s时，线性化模型和真实系统之间的误差就会迅速拉大。如果仿真任务只是验证悬停点的控制器稳定性，线性模型是够用的；可一旦涉及快速机动、抗风分析、轨迹跟踪这类大幅运动场景，线性模型的假设就撑不住了。更关键的是，很多控制器的参数在悬停点整定得再好，在大姿态角下也会失效，因为动力学本身的状态依赖关系已经变了。

我朋友那台四旋翼就是这个典型情况：悬停状态调试得很稳，PID参数放上去的姿态阶跃响应曲线看起来也漂亮，但真机做翻滚动作时，角速度一上去，机身突然剧烈抖动甚至翻转。用回放日志里的电机转速和姿态数据去驱动线性模型，发现模型输出和真机姿态完全是两回事。这时候才意识到，不是控制算法的问题，是模型根本没法描述大机动工况下的真实动力学。

1.2 非线性模型里到底装进了哪些"非线性"

那么所谓的非线性模型，到底比线性模型多了哪些东西？其实本质上就是没有省略运动方程里的完整物理项。对一个6自由度刚体飞行器来说，非线性来源主要集中在四个方面。

第一是姿态运动学里的三角函数。欧拉角速率和机体坐标系下的角速度之间有一个变换矩阵，里面含有正弦、余弦和正切项，这些函数天然是高度非线性的。姿态角一大，变换关系就不是简单的比例关系了。

第二是动力学方程里的惯性叉乘项。角动量在机体坐标系下求导时，会出现矢量叉乘形式的耦合项，比如角速度向量与惯量矩阵的乘积再叉乘角速度向量。这项描述的是旋转坐标系带来的"表观力"，飞行器在高速滚转时，这项的作用非常明显，也是大机动时线性化模型最容易丢掉的部分。

第三是执行机构的非线性。无刷电机加螺旋桨产生的推力，理想情况下近似和转速的平方成正比，转速响应本身也有延迟。这个平方关系在油门大范围变化时，绝不是一个小增量线性关系能够拟合的。

第四是气动力的非线性。低速飞行时可以近似简化，但速度稍微提高，气动阻力和速度的平方项就占了主导，线性阻力的误差可以大到几倍。

理解这些非线性来源，最大的好处是能帮你在仿真发散或结果异常时，快速判断到底是模型哪里没建对，还是算法本身有问题。

2. 建模前的物理基础：坐标系、欧拉角与刚体运动学

2.1 NED坐标系与机体坐标系的选择逻辑

在正式开始写方程之前，坐标系约定必须先定清楚，这一步如果含糊，后面所有推导都是一笔糊涂账。我用的惯性系是NED，也就是北东地坐标系：X轴指向正北，Y轴指向正东，Z轴垂直指向地面。机体坐标系的原点取在飞行器质心上，X轴指向机头方向，Y轴指向右翼，Z轴按右手定则指向机腹方向。

为什么选NED而不是大多数人习惯的ENU（东-北-天）坐标系？原因其实很实际：地面站、飞控、惯性测量单元输出的数据流方向都是按NED约定的，姿态解算出来的四元数和欧拉角也遵循这个约定。如果在仿真里用ENU，导航数据的符号和方向就得在代码里来回翻转，非常容易搞错。相信我，真实项目中因为坐标系方向写反导致仿真的水平位置方向和飞控日志反着的案例，比想象中多得多。

姿态角用Z-Y-X欧拉顺序定义：先绕惯性系Z轴偏航得到偏航角ψ，再绕新的Y轴俯仰得到俯仰角θ，最后绕新的X轴滚转得到滚转角φ。这个顺序是航空航天里的通用约定，三个旋转复合起来就构成了从机体坐标系到惯性坐标系的旋转矩阵。

旋转矩阵R的具体形式是教科书级的公式，但它本质上的作用，可以理解成一个"翻译机"：把机体坐标系下表示的向量（比如机体速度、推力方向）翻译到惯性坐标系下。翻译的规则取决于飞行器当前的头朝哪、机身倾了多少。如果姿态角很小，旋转矩阵可以近似成一个单位矩阵加一个小反对称阵，这就是线性化模型的坐标关系；而姿态角一大，三个三角函数互相耦合，只有完整的旋转矩阵才能准确描述。

2.2 运动学方程：位置和姿态的积分关系

运动学方程只关心"几何关系"，不关心力是怎么产生的。平动运动学很简单：惯性系下的位置变化率，等于机体速度向量通过旋转矩阵变换到惯性系下的结果。写成方程就是：

[ \begin{bmatrix} \dot{x} \ \dot{y} \ \dot{z} \end{bmatrix} = R \begin{bmatrix} u \ v \ w \end{bmatrix} ]

这里u、v、w是机体坐标系下的三个线速度分量。这个方程几乎没有争议，只是旋转矩阵的代入问题。

姿态运动学稍微复杂一点。很多人容易在这里犯迷糊：机体坐标系下的角速度p、q、r，和欧拉角速率ψ̇、θ̇、φ̇之间，并不是简单的相等关系。原因在于欧拉角本身是绕不同旋转轴定义的，而这套旋转轴是逐次转出来的，不是机体坐标系的固定轴。把它们投影到机体坐标系上，就得到了一个含三角函数的变换矩阵：

[ \begin{bmatrix} p \ q \ r \end{bmatrix} = \begin{bmatrix} 1 & 0 & -\sin\theta \ 0 & \cos\phi & \sin\phi\cos\theta \ 0 & -\sin\phi & \cos\phi\cos\theta \end{bmatrix} \begin{bmatrix} \dot{\phi} \ \dot{\theta} \ \dot{\psi} \end{bmatrix} ]

反解欧拉角速率时，这个矩阵的逆矩阵里会出现(\tan\theta)项。当俯仰角θ接近±90度时，正切项趋向无穷大，数值上就会爆炸——这就是欧拉角表示的奇异性问题。做悬停和普通飞行仿真时θ一般不会接近90度，所以用欧拉角没问题；但如果你要对倒飞、筋斗这类机动做仿真，就一定要换成四元数表示姿态，否则仿真会在接近奇异点时直接崩溃。

2.3 动力学方程：力、力矩与加速度的关系

动力学方程处理的是"为什么物体会动"的问题。用的是牛顿-欧拉方程，但要注意，方程必须写在机体坐标系下才最方便，因为转动惯量、推力、力矩都是相对于机体轴给出的。

平动动力学方程在机体坐标系下写为：

[ m \dot{\vec{v}} + \vec{\omega} \times (m \vec{v}) = \vec{F}_{ext} + m R^T \vec{g}_0 ]

其中(\vec{\omega})是机体角速度向量，(\vec{F}_{ext})是作用在机身上的合力（推力、阻力等），(\vec{g}_0)是惯性系下的重力加速度向量（在NED系里是([0,0,g]^T)），乘上旋转矩阵转置把它变换到机体坐标系里。

这里那个叉乘项(\vec{\omega} \times (m \vec{v}))就是惯性坐标系求导和机体坐标系求导之间的转换项。原理可以这样理解：惯性系下动量对时间求导，等于机体系下的变化率加上坐标旋转带来的"表观变化"。如果飞行器在高速旋转，哪怕机体系下的速度不变，惯性系下动量方向其实也在变，这一项正是描述这个现象的。很多简化模型会把这一项丢掉，低速飞行时可能感觉不到，但大角速度机动时，它是不可忽略的。

角运动方程也一样：

[ I \dot{\vec{\omega}} + \vec{\omega} \times (I \vec{\omega}) = \vec{M} ]

这里I是3x3转动惯量矩阵，在机体坐标系下可以近似为对角阵。(\vec{\omega} \times (I \vec{\omega}))是角动量转移项，描述的是转动坐标系下的角动量变化。一般四旋翼的Ixx和Iyy比较接近，Izz略大，这一项的耦合效应在滚转俯仰同时发生时特别明显。

把运动学方程和动力学方程合起来，六个自由度对应12个状态变量：位置x、y、z，姿态角φ、θ、ψ，线速度u、v、w，角速度p、q、r。后续所有仿真实现都围绕这12个状态展开。

3. 完整的六自由度非线性方程组与参数来源

3.1 作用在机身上的力与力矩清单

要把方程写完整，关键是把所有作用在机身上的力和力矩列全。很多初版模型要么少算了一项，要么把方向搞反，结果仿真出来的行为非常诡异。完整的受力清单如下：

重力。大小是质量乘以重力加速度，方向在NED系下恒为Z轴正方向（向下）。在机体坐标系下表示时，需要把重力向量乘上旋转矩阵的转置，所以姿态角会直接影响重力在机体轴上的分量。这个投影关系看似简单，实际上是大姿态角模型的核心特征之一。

旋翼推力。每个旋翼产生的推力沿机体Z轴负方向（也就是向上方向），大小与电机转速的平方近似成正比：(F_i = k_F \omega_i^2)，其中kF是推力系数。四个旋翼合力之和就是：

[ F_{thrust} = k_F (\omega_1^2 + \omega_2^2 + \omega_3^2 + \omega_4^2) ]

反扭矩。每个旋翼转动时，空气会给机身一个反方向力矩，大小为(M_i = k_M \omega_i^2)，方向与旋翼转向相反。四旋翼相邻电机转向相反，正是为了在悬停时让反扭矩相互抵消。这个力矩差是偏航控制的核心。

陀螺力矩。螺旋桨本身是一个高速旋转的转子，当机身姿态变化时，转子角动量方向跟着改变，会产生陀螺效应力矩。数学上是(M_{gyro} = J_p \cdot \vec{\omega}_b \times \vec{\Omega})，其中Jp是单个旋翼转子的转动惯量，(\vec{\Omega})是电机转子的转速矢量（沿机体Z轴）。这个力矩在快速滚转俯仰时会体现为姿态交叉耦合，做精准控制仿真时必须包含。

气动阻力。完整的多旋翼气动阻力建模非常复杂，工程上常采用简化模型：低速时用线性阻尼(F_d = -k_d \vec{v})，中高速时用二次阻力(F_d = -\frac{1}{2}\rho C_d A |v| \vec{v})。在小型无人机仿真里，线性阻尼项通常足够描述低速悬停附近的特性，但如果做高速前飞或抗风仿真，就得用二次项才比较贴合实际。

3.2 参数测定与估算：一个典型小型四旋翼的数值

模型参数对不对，直接决定仿真结果可不可信。参数不能拍脑袋随意估，至少要有一个来源依据。以我这次搭建的轴距450mm小型四旋翼为例，关键参数如下表：

这里转动惯量我强烈建议实测，不要只靠SolidWorks估算。空机CAD模型算出来的转动惯量往往偏低，因为电调、电池扎带、减震球这类装配件在3D模型里不够精确，而且电池的实际布局位置对惯量影响很大。三线摆的测量方法其实不复杂：做一个圆形托盘，用三根等长细线悬挂起来，把无人机放在托盘上让系统做微小扭转振动，测周期就可以算出总惯量，再减去托盘自身的惯量，就得到无人机惯量。这个方法精度足够工程使用，成本也低。

推力系数需要在静拉力台上标定。把电机和桨正装在一个力传感器上，给不同油门信号，测量对应推力。转速可以用转速计或者电调遥测数据读出来，然后拟合(F = k_F \omega^2)。反扭矩系数更难直接测，通常的做法是：让四旋翼在真机悬停（此时偏航力矩平衡），通过飞控日志读出四个电机的转速差，反推出kM的大致范围。还有一个更省事的办法，用厂家提供的螺旋桨推力扭矩系数数据表来查。没有条件实测时，至少也要保证参数的量级合理，否则仿真悬停需要的转速和真机差太多，那模型就没有意义了。

3.3 完整状态方程的组装

所有力和力矩都齐了之后，就可以组装成12维非线性状态方程了。为了后续在Matlab里写代码方便，我习惯用一个结构体数组来组织状态向量：

[ x = [x,y,z,\phi,\theta,\psi,u,v,w,p,q,r]^T ]

然后每个状态的变化率写成显式方程。这里全部写出来：

位置方程： [ \dot{x} = (\cos\theta\cos\psi)u + (\sin\phi\sin\theta\cos\psi - \cos\phi\sin\psi)v + (\cos\phi\sin\theta\cos\psi + \sin\phi\sin\psi)w ] [ \dot{y} = (\cos\theta\sin\psi)u + (\sin\phi\sin\theta\sin\psi + \cos\phi\cos\psi)v + (\cos\phi\sin\theta\sin\psi - \sin\phi\cos\psi)w ] [ \dot{z} = (-\sin\theta)u + (\sin\phi\cos\theta)v + (\cos\phi\cos\theta)w ]

欧拉角速率方程（把姿态运动学矩阵求逆）： [ \dot{\phi} = p + q\sin\phi\tan\theta + r\cos\phi\tan\theta ] [ \dot{\theta} = q\cos\phi - r\sin\phi ] [ \dot{\psi} = q\frac{\sin\phi}{\cos\theta} + r\frac{\cos\phi}{\cos\theta} ]

线速度方程： [ \dot{u} = rv - qw - g\sin\theta ] [ \dot{v} = pw - ru + g\sin\phi\cos\theta ] [ \dot{w} = qu - pv + g\cos\phi\cos\theta - \frac{k_F}{m}(\omega_1^2+\omega_2^2+\omega_3^2+\omega_4^2) ]

角速度方程： [ \dot{p} = \frac{I_{yy}-I_{zz}}{I_{xx}}qr - \frac{J_p}{I_{xx}}q\Omega_{sum} + \frac{l k_F}{I_{xx}}(\omega_2^2 - \omega_4^2) ] [ \dot{q} = \frac{I_{zz}-I_{xx}}{I_{yy}}pr + \frac{J_p}{I_{yy}}p\Omega_{sum} + \frac{l k_F}{I_{yy}}(\omega_1^2 - \omega_3^2) ] [ \dot{r} = \frac{I_{xx}-I_{yy}}{I_{zz}}pq + \frac{k_M}{I_{zz}}(\omega_1^2 - \omega_2^2 + \omega_3^2 - \omega_4^2) ]

其中(\Omega_{sum} = \omega_1 - \omega_2 + \omega_3 - \omega_4)，注意这里是带符号的转速和，和电机转向有关。这个方程组的推导过程虽然长了点，但每一步都有明确的物理来源，后面在Matlab里实现时就变成了一个几乎可以"翻译"成代码的公式清单。

4. Simulink模型搭建：从方程到仿真模块

4.1 总体架构设计：为什么我要拆成三个子系统

在Simulink里搭模型，我强烈建议不要把所有方程塞进一个巨大的框图里，而是按物理逻辑拆成三个子系统：力与力矩计算模块、刚体动力学模块、运动学更新模块。

之所以这样拆，最直接的好处是复用性强。力与力矩模块的输入是四个电机转速命令，输出是合力和合力矩；刚体动力学模块的输入是合力和合力矩，输出是机体线加速度和角加速度；运动学更新模块的输入是线速度和角速度，输出是位置和姿态角。这样拆完之后，每个模块都可以单独测试。更实际的好处是，后面要接入控制器时，你只需要在力与力矩模块前加一个控制信号接口，完全不用改动动力学核心；要做故障注入（比如某个电机失效），也只需要在力与力矩模块里单独处理对应电机的出力就行。

顶层接口设计如下：模型最外层输入是四个电机转速命令(\omega_1)到(\omega_4)，输出是12维状态向量，也可以直接把位置和姿态角单独引出来给Scope或者Record模块。

第二个好处是方便调试。仿真结果不对时，可以先把每个子系统的中间输出拉出来看。比如机体角速度是否符合物理直觉，推力计算有没有在某个转速区间出现明显突变，这些中间层信号的可见性，能帮你省下大量排查时间。

4.2 MATLAB Function实现核心方程

实现方式上，我选择用MATLAB Function模块写状态方程，而不是用一堆积分器、增益、求和模块去搭框图。原因很简单：对12维非线性方程来说，框图的连线复杂度太高，容易连错且不易检查；而MATLAB Function里的代码几乎可以直接对照第三节的公式逐行翻译，可读性和可维护性都强得多。

下面这段代码就是核心状态方程的完整实现。我把它写在MATLAB Function模块里，输入是时间t（实际未用）、状态向量x和参数结构体params，输出是状态导数x_dot。

function x_dot = quadrotor_dynamics(x, params) % 状态向量解包 % x = [x,y,z,phi,theta,psi,u,v,w,p,q,r] phi = x(4); theta = x(5); psi = x(6); u = x(7); v = x(8); w = x(9); p = x(10); q = x(11); r = x(12); % 电机转速输入（由控制器或参数配置给出） w1 = params.w1; w2 = params.w2; w3 = params.w3; w4 = params.w4; % 三角函数预计算，避免重复计算 sp = sin(phi); cp = cos(phi); st = sin(theta); ct = cos(theta); ss = sin(psi); cs = cos(psi); % 旋转矩阵R（机体坐标系->NED惯性系） R11 = ct*cs; R12 = sp*st*cs - cp*ss; R13 = cp*st*cs + sp*ss; R21 = ct*ss; R22 = sp*st*ss + cp*cs; R23 = cp*st*ss - sp*cs; R31 = -st; R32 = sp*ct; R33 = cp*ct; % 位置导数 x_dot = zeros(12,1); x_dot(1) = R11*u + R12*v + R13*w; x_dot(2) = R21*u + R22*v + R23*w; x_dot(3) = R31*u + R32*v + R33*w; % 姿态角导数（欧拉角速率方程） if abs(ct) < 1e-4 % 接近俯仰±90度，欧拉角奇异，这里做保护 x_dot(4) = 0; x_dot(5) = 0; x_dot(6) = 0; else x_dot(4) = p + q*sp*tan(theta) + r*cp*tan(theta); x_dot(5) = q*cp - r*sp; x_dot(6) = (q*sp + r*cp)/ct; end % 合力与合力矩 Fz = -params.kF * (w1^2 + w2^2 + w3^2 + w4^2); % 机体Z轴负方向推力 tau_phi = params.l * params.kF * (w2^2 - w4^2); tau_theta = params.l * params.kF * (w1^2 - w3^2); tau_psi = params.kM * (w1^2 - w2^2 + w3^2 - w4^2); % 陀螺力矩（Omega_sum为带符号转速和） Omega_sum = w1 - w2 + w3 - w4; gyro_p = -params.Jp * q * Omega_sum / params.Ixx; gyro_q = params.Jp * p * Omega_sum / params.Iyy; % 线速度导数 x_dot(7) = r*v - q*w - params.g * st; x_dot(8) = p*w - r*u + params.g * sp * ct; x_dot(9) = q*u - p*v + params.g * cp * ct + Fz/params.m; % 角速度导数 x_dot(10) = ((params.Iyy - params.Izz)/params.Ixx)*q*r + gyro_p + tau_phi/params.Ixx; x_dot(11) = ((params.Izz - params.Ixx)/params.Iyy)*p*r + gyro_q + tau_theta/params.Iyy; x_dot(12) = ((params.Ixx - params.Iyy)/params.Izz)*p*q + tau_psi/params.Izz; end

这段代码里有两个细节值得特别说明。第一，三角函数项我全部预先算好再复用，避免在每个方程里重复计算，虽然对单个仿真影响不大，但做蒙特卡洛参数扫描时能明显省时间。第二，姿态角导数里加了俯仰角接近90度时的保护，实际项目中欧拉角奇异是仿真崩溃最常见的原因之一，宁可在这里做保护，也不能让积分器直接算出一个NaN。

线速度方程里的Fz是四个推力求和，方向沿机体Z轴负方向，也就是向上。很多新手容易在这里把符号搞反，结果仿真一开始无人机就往下掉，怎么调参数都救不回来。

4.3 初始化脚本、配平初值与求解器设置

Simulink模型里的参数不应该硬编码在模块里，而是通过一个初始化脚本统一管理。我习惯用一个.m文件定义所有的参数并打包成结构体：

% quad_initialize.m params.m = 1.5; % kg params.g = 9.81; % m/s^2 params.l = 0.25; % m params.Ixx = 0.021; % kg*m^2 params.Iyy = 0.024; % kg*m^2 params.Izz = 0.037; % kg*m^2 params.kF = 1.63e-6; % N/(rad/s)^2 params.kM = 3.6e-8; % N*m/(rad/s)^2 params.Jp = 1.2e-5; % kg*m^2 params.kd = 0.03; % N*s/m % 悬停配平转速：总推力=重力 % m*g = kF * 4 * w_hover^2 w_hover = sqrt(params.m * params.g / (4 * params.kF)); params.w1 = w_hover; params.w2 = w_hover; params.w3 = w_hover; params.w4 = w_hover;

配平初值这一步非常关键。直接给一个所有状态都为0的初值，电机转速也给一个随意值，那仿真一开始的加速度就会非常大，飞机会迅速坠落，曲线很难看，也看不出真实动力学特性。正确的做法是先把悬停转速算出来，让初始时刻合力矩为零、推力等于重力，然后在这个配平点附近做扰动实验。这样模型在时间原点附近的行为才可控，便于分析。

求解器设置方面，如果只是做一般性的开环动力学仿真，固定步长ode4（四阶龙格库塔）加1毫秒步长就够了。四旋翼的刚体动力学带宽一般在10Hz到50Hz这个量级，1ms的采样足够覆盖。但如果后面把电机动态响应也建模进去，电机转速环的响应时间可能只有几十毫秒，这时候建议把步长缩到0.1ms或0.2ms，否则高频动态会失真。做硬件在环仿真时，步长还必须要和飞控的控制频率对齐，一般250Hz或500Hz，这个要单独确认。

5. 仿真分析与模型校核：模型建得对不对

5.1 开环悬停响应：什么样的结果算"模型正确"

模型搭完，第一步要做的不是急着接控制器，而是先跑开环悬停扰动响应。具体做法是：用配平初值，给一个微小初始扰动，比如让滚转角初始值设为2度，然后观察姿态角的响应曲线。

这里有个反直觉的经验：一个没接控制器的四旋翼开环模型，姿态响应应该是发散的或者等幅振荡的，而不是稳定的。原因很简单，四旋翼本身是一个静不稳定的系统，没有姿态反馈，任何扰动都会让它继续偏离。如果你发现开环悬停仿真里姿态角自己慢慢回到了0，那大概率模型里有问题，可能是陀螺力矩符号搞反了，或者某个阻尼项加得过大，掩盖了真实动力学特性。

以我搭建的模型为例，给2度滚转角初始偏差，仿真10秒，滚转角会呈振荡发散趋势，同时偏航角也会因为反扭矩耦合出现漂移。这才是健康模型的输出。如果曲线一开始就冲出去变成NaN，那通常是数值问题，优先检查姿态角导数里的奇异保护有没有生效，以及积分器步长是否合适。

5.2 非线性模型与线性化模型的对比实验

模型通过开环合理性检查之后，接下来做线性化对比实验。这一步主要是为了验证非线性模型的正确性，同时也能量化线性化模型的误差边界。

在Simulink里可以把模型在悬停配平点做线性化分析，得到状态矩阵A，再画出线性模型的阶跃响应，与非线性模型在同样输入下的响应做对比。小扰动情况下，比如滚转阶跃指令只有5度，两个模型的响应曲线几乎重合；把阶跃指令加大到30度，非线性模型和线性模型的峰值角速度、超调量就开始有肉眼可见的差异；指令到60度时，两个模型的响应行为就完全不同了，线性化模型可能预测出一个稳定收敛的响应，而非线性模型则表现出明显的耦合与发散。

这个对比实验的价值在于：第一，它验证了非线性模型在平衡点附近和线性模型的一致性，说明方程推导没有大问题；第二，它给出了线性化模型的适用范围，为后续控制器设计提供参考边界。做论文或者技术报告时，这组对比数据也是非常有力的素材。

5.3 与真实飞行日志校核：参数修正的闭环

纯仿真层面的验证还不够，模型最终要能用，必须和真机数据校核。方法不复杂：把一次真实试飞的飞控日志导出来，取出电机转速指令、姿态角、角速度这些信号。然后把电机转速作为模型输入，在Matlab里重新跑一遍仿真，对比模型输出的姿态角和真机记录的姿态角。

第一次对比通常会有偏差，原因就在气动阻尼系数、推力系数这些参数和真实值有差距。这里有一个很实用的调参思路：先观察悬停段的误差，如果模型姿态振荡频率比真机高，说明转动惯量偏小了；如果模型响应比真机迟钝，说明惯量偏大或者阻尼系数给大了。然后做一段机动段对比，重点看大姿态角下误差走向，如果误差在大姿态时快速增大，优先检查陀螺力矩项和气动阻力模型是否需要从线性升级成二次阻尼。

这个过程本质上是一个以模型输出和真机输出之间的误差最小化为目标的参数辨识，不需要什么复杂算法，手动根据试飞数据逐项修正几轮，精度就能达到工程可用的水平。以我的经验，经过三四轮修正后，模型姿态输出和真机日志的平均误差能控制在5度以内，这足够支撑控制器参数预调了。

5.4 仿真发散原因排查清单

以我踩过的坑和帮别人排查的经验，把最常见的问题列成一张表，方便仿真异常时快速定位：

这里面最容易被忽略的是单位问题。我见过一个同学在Simulink里把角度初值设置成了30，运行几秒姿态就完全乱掉，排查到最后发现其实他把30度直接当成30弧度用了。这在仿真里是灾难性的错误，Matlab本身不会提示，因为数值上积分器照样能跑，只是结果完全失真。

排查的方法论其实就一句话：不要先怀疑求解器和积分器，先相信你的模型方程和参数。大多数发散问题都出在建模阶段——符号方向、单位、参数量级这三类错误占据了九成以上的排查时间。

搞完模型校核之后，这个非线性6自由度模型就可以干很多事了：接PID控制器做姿态轨迹仿真、做电机失效故障注入、做参数敏感性分析、甚至配合硬件在环平台做半实物仿真。我个人的体会是，模型的价值不在于"复杂"本身，而在于你遇到问题时能快速定位到底是模型假设错误还是控制器设计错误。从那次帮朋友调试四旋翼之后，我基本养成了习惯：任何飞行器项目起步，先在Matlab里搭完整非线性模型，哪怕真机还没装配完成，模型也能把很多问题提前暴露出来。这套流程对一个刚开始做无人机建模的人而言，一开始可能会觉得方程推导枯燥、参数标定繁琐，但咬牙把模型搭完并且和真机数据校核通过之后，后续的控制算法开发会很轻松。

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