简介：一套面向固定翼无人机系统建模与仿真的MATLAB程序包，适合航空航天工程教学、无人机系统研发及自动控制理论研究人员使用。包内整合六自由度非线性动力学模型、自主导航控制系统、三维路径规划算法、传感器数据融合模块与飞行姿态控制器等核心模块，采用模块化架构，各子系统可独立运行，并配有空气动力学系数数据库、多控制策略对比分析工具、实时仿真与离线数据处理能力，以及可视化分析界面。通过修改参数配置，可针对不同机型特性开展仿真验证，为固定翼无人机控制系统开发提供完整理论验证平台。资源共54个文件，以m格式脚本为主，辅以zbak备份、asv自动保存、txt说明及zip备份包，整体仅52KB，结构紧凑。已有124人学习下载，便于快速上手与二次开发。

1. 固定翼无人机建模与仿真MATLAB实现：实机试飞前，先把飞机在电脑里摔一遍

固定翼无人机建模与仿真MATLAB实现，说直白点就是：真机还没上跑道，先用数学方程和一组气动参数，在电脑里把飞机的每一次抬头、偏航、滚转都算出来。我接触这个方向是因为飞控算法调试，参数在实机上乱调一次就得提心吊胆半天，后来改成先在仿真里把问题逼出来。这套方案能解决三件事：验证飞行包线内飞机是否稳定、给飞控算法提供一个能反复“炸机”的测试床、以及为后续硬件在环仿真准备一个可信的飞机模型。适合正在做固定翼飞控、飞行器总体设计或者毕业设计里涉及飞行动力学的学生和工程师。而且真正做过之后你会发现，难点根本不在六自由度方程，而在：气动数据怎么给、配平初值怎么猜、线性化结果怎么验证。

2. 从机体坐标系入手：固定翼六自由度运动方程怎么落进MATLAB

2.1 先定坐标系约定：这一步定了，后面所有代码都不用返工

写固定翼模型前，我吃过一次大亏：坐标轴方向没统一，导致气动力符号反了，仿真出来的飞机一给升降舵就往下栽，排查了两天才发现是滚转角方向约定不一致。所以现在无论项目多急，我都会先花十分钟把坐标系约定写死在代码注释里。

最常见的是北美惯用约定：地面坐标系用NED（北-东-地），机体坐标系原点在重心，x轴指向机头，y轴指向右翼，z轴指向机腹。欧拉角顺序为yaw-pitch-roll，也就是先绕z轴偏航，再绕y轴俯仰，最后绕x轴滚转。这套约定和大多数飞控代码、以及MATLAB航空航天工具箱里的方向余弦矩阵是一致的。

对应的状态量我一般选12维：线速度u、v、w（机体轴），角速度p、q、r（滚转、俯仰、偏航），位置pN、pE、高度h（或高度pd），以及欧拉角phi、theta、psi。注意一个细节：高度方向，如果你喜欢用负的“地”位置，那整个传感器模型和对地速度都要跟着变号，容易乱。我更推荐直接把高度写成正的h，单独作为一个状态，重力势能在方程里单独处理。

% 固定翼六自由度状态向量定义 % x(1) u 机体轴x向速度 (m/s) % x(2) v 机体轴y向速度 (m/s) % x(3) w 机体轴z向速度 (m/s) % x(4) p 滚转角速率 (rad/s) % x(5) q 俯仰角速率 (rad/s) % x(6) r 偏航角速率 (rad/s) % x(7) pN 北向位置 (m) % x(8) pE 东向位置 (m) % x(9) h 高度 (m)，正值向上 % x(10) phi 滚转角 (rad) % x(11) theta 俯仰角 (rad) % x(12) psi 偏航角 (rad)

这段注释不是摆设，它决定了你写动力学方程时每一项的正负号。特别是w朝下为正还是朝上为正，直接决定攻角alpha = atan2(w, u)还是atan2(-w, u)——写反了，升力方向就反了，整个模型就是一台“反着飞的飞机”。我在团队里见过不止一个人在这里栽跟头，最后全是靠对照单位阶跃响应才查出来。

2.2 六自由度运动方程的最小实现：一个函数文件就是一个飞机

方程本身是牛顿欧拉方程，没什么花活，关键是用MATLAB写得体。我不推荐在Simulink里从头画积分器连线的第一个版本，调试起来太痛苦；先用一个m函数把飞机动力学写出来，跑通后再决定要不要搬到Simulink里。

function dx = fixedwing_dynamics(t, x, param) % 六自由度刚体动力学 + 重力 + 气动力/力矩 % param是结构体，包含质量、转动惯量、气动系数等 % 拆状态 u = x(1); v = x(2); w = x(3); p = x(4); q = x(5); r = x(6); phi = x(10); theta = x(11); psi = x(12); % 体轴系下重力分量 Fg = param.m * 9.80665 * ... [-sin(theta); cos(theta)*sin(phi); cos(theta)*cos(phi)]; % 计算空速、攻角、侧滑角 V = sqrt(u^2 + v^2 + w^2); alpha = atan2(w, u) * (V > 1e-3); % 低速时避免除零 beta = asin(v / max(V, 1e-3)); % 气动力/力矩，单独函数封装 [Fa, Ma] = aero_force_moment(V, alpha, beta, p, q, r, param); % 线运动方程（体轴系） F_total = Fg + Fa + param.F_engine; du = F_total(1)/param.m - q*w + r*v; dv = F_total(2)/param.m - r*u + p*w; dw = F_total(3)/param.m - p*v + q*u; % 角运动方程（体轴系，假设惯量积为0） dphi = p + tan(theta)*(q*sin(phi) + r*cos(phi)); dtheta = q*cos(phi) - r*sin(phi); dpsi = (q*sin(phi) + r*cos(phi)) / max(cos(theta), 1e-3); % 位置与高度 dpN = (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; dpE = (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; dh = u*sin(theta) - v*sin(phi)*cos(theta) - w*cos(phi)*cos(theta); dx = [du; dv; dw; dp; dq; dr; dpN; dpE; dh; dphi; dtheta; dpsi]; end

代码逻辑上，这个函数把飞机当成一个带质量、转动惯量的刚体，外部只输入一个结构体param。这样每个模块独立，后续换飞机型号只改param，不改方程。参数说明：param.m是质量，F_engine按体轴系三个分量给；惯量积我默认忽略，这对常规固定翼够用，但如果做的是带大非对称载荷的飞机，必须把惯量积Ixz加回来，否则螺旋桨滑流带来的不对称力矩会算不准。

2.3 气动力和力矩怎么塞进方程：先做接口，再填数据

在写aero_force_moment之前，先定义清楚输入输出。我这里用的是最常见做法：输入空速、攻角、侧滑角、角速率和舵面偏角，输出体轴系下的力和力矩。这样做的原因是，后续无论你是用线性气动系数、风洞数据查表，还是CFD计算出来的数据集，都可以替换同一个函数而不动主程序。

function [Fa, Ma] = aero_force_moment(V, alpha, beta, p, q, r, param) % 气动力与气动力矩，输出在体轴系 % 动压 qbar = 0.5*rho*V^2 qbar = 0.5 * param.rho * V^2; S = param.S; b = param.b; c = param.c; % 风轴系气动力系数 CL = param.CL0 + param.CLalpha * alpha + param.CLq * q * c / (2*V); CD = param.CD0 + param.K * (CL - param.CLmin) .^ 2; % 极曲线近似 CY = param.CYbeta * beta + param.CYp * p * b / (2*V) + param.CYr * r * b / (2*V); % 风轴到体轴转换 Lwind = qbar * S * CL; Dwind = qbar * S * CD; Ywind = qbar * S * CY; Fx_a = -Dwind * cos(alpha) + Lwind * sin(alpha); Fy_a = Ywind; Fz_a = -Dwind * sin(alpha) - Lwind * cos(alpha); % 力矩系数 Cl = param.Clbeta * beta + param.Clp * p * b/(2*V) + param.Clda * param.da; Cm = param.Cm0 + param.Cmalpha * alpha + param.Cmq * q * c/(2*V) + param.Cmde * param.de; Cn = param.Cnbeta * beta + param.Cnr * r * b/(2*V) + param.Cnda * param.da; Ma = [Cl * qbar * S * b; Cm * qbar * S * c; Cn * qbar * S * b]; Fa = [Fx_a; Fy_a; Fz_a]; end

这段代码有两个参数容易被忽略。第一个是CD用极曲线近似：CD0 + K*(CL-CLmin)^2，真实飞机手册里一般给的是这个形式，而不是直接给一个固定CD。第二个是角速率产生的阻尼项：CLq、CYp、CYr、Cmq、Cnr这些必须无量纲化，用c/(2V)或b/(2V)乘进去，不处理单位的话，仿真的动态响应会明显偏“活跃”。

3. 气动数据从哪儿来：线性气动系数把第一版模型撑起来

3.1 线性气动模型够不够用：CFD和风洞数据到位之前的最佳选择

很多人在建模之初就卡在“没有气动数据”。实机厂家通常只给一个升力系数曲线，力矩系数给得很保守，CFD自己算又没把握。我的建议是，第一版不要追求高精度，先用线性气动系数把模型跑起来。固定翼在小攻角、未失速的范围内，升力线斜率、阻尼导数这些线性模型表现足够好，用于飞控设计验证完全够用。等你需要模拟失速、大侧滑这些非线性边界时再把系数表换成插值表，整体架构不用动。

线性气动模型的核心参数大概有十几项：升力线斜率CLalpha、零攻角升力CL0、俯仰力矩斜率Cmalpha、升降舵效率Cmde、滚转阻尼Clp、偏航阻尼Cnr、侧滑刚度Cnbeta、滚转刚度Clbeta等等。这些参数网上能找到一些通用小型无人机的参考值，但最好对照你手头飞机的翼型数据做一次量纲检查。我把它们维护在一个结构体里，作为唯一的“飞机事实来源”。

3.2 系数查表脚本：把线性参数和插值表统一成一个接口

如果你拿到的气动数据是表格形式，例如风洞报告里alpha从-10度到20度每隔2度一个CL值，那线性函数就不够用了。我在上级项目里的做法是写一个统一的查表函数，内部用MATLAB自带的interp1或者griddedInterpolant，支持一维表（攻角）和二维表（攻角+侧滑角），外部接口仍然返回力和力矩。

function [CL, CD, CY] = aero_coeff_table(alpha_deg, beta_deg, param) % 从气动数据库插值，alpha_deg和beta_deg单位是度 % 小于0.2度直接用线性模型，否则查表 if ~isfield(param, 'CL_alpha_table') % 没有表格就退回线性模型 CL = param.CL0 + param.CLalpha * deg2rad(alpha_deg); CD = param.CD0; CY = param.CYbeta * deg2rad(beta_deg); return; end % alpha表是列向量，CL表按alpha行排列 CL = interp1(param.CL_alpha_table(:,1), param.CL_alpha_table(:,2), ... clamp(alpha_deg, param.alpha_min, param.alpha_max), 'linear', 'extrap'); CD = interp1(param.CD_alpha_table(:,1), param.CD_alpha_table(:,2), ... clamp(alpha_deg, param.alpha_min, param.alpha_max), 'linear', 'extrap'); % 侧滑角用二阶多项式修正，避免二维表过大 CY = (param.CYbeta * deg2rad(beta_deg)) * ... (1 + param.CY_beta2 * abs(deg2rad(beta_deg))); end

注意我在这里用了一个小技巧：侧滑角方向的数据如果也是完整二维表，会导致查表数据量很大而且需要风洞给足够多的beta状态；很多项目其实只有beta一维导数，所以用CYbeta乘beta，加上一个二阶修正项模拟方向舵或机身侧向力非线性。逻辑说明：interp1的最后一个参数我写了'extrap'，这意味着即使仿真因为某些原因短暂计算出攻角25度，也不会返回NaN导致整个积分崩溃，而是线性外推一个值。但从工程角度这其实是“掩耳盗铃”，所以我在后面会加一个总攻角警告，一旦超过失速攻角就打印提示。

3.3 一组能直接用于第一版仿真的参数参考

下表是一架翼展约2米、起飞重量约8公斤的小型固定翼无人机的典型参数，适合用来跑通模型、验证配平，不能直接当真实飞机的气动数据用。量纲统一：角度相关导数用rad^{-1}，角速率阻尼也用rad^{-1}。

拿到这张表之后，第一步不是急着仿真，而是用直觉检查一遍：Cmalpha负值说明飞机低头方向有恢复力矩，这是静稳定飞机的基本特征；如果你在代码里随手填了一个正值，那配平出来的结果会是尾翼产生负升力来配平，不是不能用，但整个飞机的响应会非常怪异。

4. 让仿真闭环跑起来：配平、线性化与小扰动响应

4.1 配平初值怎么选：手算一个合理的“顺风起飞”状态

配平就是求一组状态和舵面偏角，使得dx/dt等于零。对常规平飞来说，意味着高度不变、速度不变、俯仰角不变。这个初值影响极大——给不好的话，数值求解器会一直找不到平衡点，或者找到一个完全没意义的解（比如负攻角飞行）。

我一般先手算粗值。平飞时升力等于重力：CL = (2mg)/(rhoSV^2)，然后根据升力线斜率反推alpha。假设巡航速度25m/s，空气密度1.225kg/m^3，按上表参数，CL约等于0.93，对应alpha约8度。这个攻角合理，因为小型低速无人机巡航攻角本来就在5到10度之间。然后推力估算成配平阻力，初始俯仰角取等于攻角。把这些塞进状态初值：

% 巡航配平手算初值 V0 = 25.0; % 巡航速度 m/s gamma0 = 0; % 飞行路径角 = 0 平飞 alpha0_deg = 8; % 由升力平衡手算得到，单位度 theta0 = gamma0 + deg2rad(alpha0_deg); % 俯仰角 % 状态初值：稳定平飞 x0 = [V0*cos(theta0); ... % u 0; ... % v V0*sin(theta0); ... % w 0; 0; 0; ... % p, q, r 0; 0; 100; ... % pN, pE, 高度100m 0; theta0; 0]; % phi, theta, psi % 舵面初值 de0 = 0; % 升降舵先给0，等配平求解器调 th0 = 0.5; % 油门，无量纲0~1，等求解器调

注意这里我没直接把升降舵设成“正确值”，而是让求解器去调。原因很简单：手算可以估算攻角，但你算不出舵偏角，因为Cm0和Cmde的影响必须一起去配平。如果你对MATLAB不熟，配平求解可以先用fsolve：变量选[theta, de, th]三个，方程选dudt=0, dwdt=0, dqdt=0，高度和速度由x0固定。三个方程三个未知数，非常好收敛。

4.2 用数值线性化把飞机变成A、B矩阵：状态空间才是控制器设计的地基

有了配平点，下一步要拿状态空间模型。很多教材直接摆A矩阵公式，但真到你的飞机参数，公式展开容易错。我的做法是有限差分求雅可比：在配平点给每个状态一个微小扰动（比如1e-5），跑一次动力学函数，记录dx的变化量，除以扰动量，得到A矩阵。B矩阵同理，扰动升降舵和油门。

这段代码就是把非线性模型“降级”成线性模型的关键，控制器设计的极点配置、LQR、PID参数初估值都从这里来。要注意，有限差分步长不能太大也不能太小：太大则把非线性失真带进来，太小则数值精度不够。1e-5到1e-6一般是安全区间，如果模型本身用了很多插值表，步长还要再大一点免得触发查表跳变。

4.3 小扰动仿真：给飞机一个“轻推”，看它回不回来

配平线性化做完后，第一个仿真实验不是给大指令，而是给一个小扰动。比如在配平基础上把俯仰角速率q抬高一点，或者给升降舵一个短暂的0.5度偏转再回中，看飞机的响应是否衰减。这是检验模型物理自洽性最有效的手段。

% 在配平点做小扰动仿真 dt = 0.01; t_end = 20; t = 0:dt:t_end; % 初始化状态为配平点 x = repmat(x_trim, 1, length(t)); x(11, 1) = x_trim(11) + deg2rad(2); % 初始俯仰角加2度扰动 % 用ode4固定步长积分，或直接调用ode45 options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6); [t_out, x_out] = ode45(@(tt, xx) fixedwing_dynamics(tt, xx, param_trim), ... [0 t_end], x(:,1), options); % 画俯仰角响应 plot(t_out, rad2deg(x_out(:,11))); xlabel('时间 (s)'); ylabel('俯仰角 (deg)');

正常情况下，静稳定飞机的俯仰角应该先有一个短周期振荡（大概1-2Hz），然后快速回到配平值；如果它缓慢发散或者高频振荡越来越猛，先别急着调控制器，回过去查气动阻尼导数的符号和量纲。我见过太多人在这步直接跳到PID调试，结果控制器参数怎么整定都不对，其实是飞机模型本身就是不自洽的。

5. 固定翼建模与仿真避坑手册：五个最容易翻车的现场

5.1 现象：俯仰角一推杆，欧拉角解算突然变成NaN

原因：我在2.2节代码里写dpsi的分母是cos(theta)。当俯仰角接近正负90度时，cos(theta)趋近0，偏航角速率算式出现无穷大，数值积分器直接炸掉，状态变成NaN。固定翼正常飞行很少到这，但做大机动仿真、或者初始随机状态设置得太野时一定会触发。解决：一是把欧拉角替换成四元数做姿态解算，通用性最好；二是在dpsi计算中加一个阈值保护，如max(cos(theta), 0.1)，但这样会引入明显误差。我的实际选择是：常规设计验证用阈值保护，因为够简单；但只要是做机动包线仿真，就老老实实上四元数，这个坑一次都不想踩第二次。

5.2 现象：配平求出来的升降舵偏角大得离谱，和手算对不上

原因：气动系数里角度单位混了。最常见的是alpha表用度，而CLalpha定义成1/度而不是1/rad，或者说明书里给的是CLalpha = 4.8 /rad，你在代码里直接把deg2rad(alpha)漏了，导致攻角看起来比实际大57倍。遇到配平结果反常，第一步就是检查单位。解决：把单位统一写在结构体字段里，比如param.alpha_unit = 'rad'，并在查表函数的入口强制转换一次。我会在测试脚本里故意输入alpha=0.1，手动算出预期的CL值做断言校验，这类问题一次就能拦下来。

5.3 现象：trim函数报错“no solution found”，或者收敛到一个负速度解

原因：配平初值选得离真实平衡点太远。特别是横向状态v、p、r不为零的时候，求解器会乱找。另一个原因是升降舵初值给了一个过大的偏角，导致气动力矩曲线进入非线性区，没有平衡点。解决：配平前先手工把横向状态清零，纵向用我4.1节的方法手算一轮。如果还不行，就扫一遍速度：从15m/s到35m/s每隔2m/s设一个初值跑配平，看哪个收敛结果最符合预期，再把那个解作为正式初值。这个“扫速度配平”的方法几乎不用动脑，但特别管用。

5.4 现象：同一个模型，用ode45和固定步长ode4跑出来的结论完全相反

原因：固定翼模型里存在高频气动阻尼模态，比如荷兰滚或者短周期，时间常数可能只有0.2秒。用默认的变步长求解器看似精确，但如果步长控制不严，高频模态可能被数值方法“吃掉”，或者相反被放大。我遇到过一次：ode45说飞机稳定，换成Simulink里固定步长1ms跑，结果轻微振荡不衰减。解决：在正式做控制器验证之前，固定一个步长跑一遍（Simulink里推荐0.01s起步，做电机动态仿真则要0.001s），然后逐步减半步长确认结果不变。如果结果变了，说明你在和一个“数值假象”在战斗。

5.5 现象：仿真越来越慢，原来卡在气动系数函数里那行循环

原因：气动数据从线性换成插值表之后，如果在查表函数里写了“遍历所有表格行找最近邻”的手写循环，那每次调用几千次，ode45又要调用上万次，累积下来速度会非常难看。我早期干过这种蠢事，用length+for循环一个个找插值区间，被同事看到直接笑了。解决：不要手写插值，直接用interp1或griddedInterpolant，后者在多次查同一个表时还能缓存分片，速度优势明显。另一个被忽视的地方是，把舵面偏角、推力这些输入从结构体里取字段时，尽量一次取出来放在局部变量，而不是在循环内反复写param.xxx，逻辑一样，速度差几倍。

6. 再往前走一步：用结构体把气动模型组织成可替换接口，顺便接上STK联合仿真

固定翼建模与仿真做到能配平、能线性化、能输出响应，其实才完成了最基础的一步。如果你想让它成为一套能复用两三年的工具，我建议尽快做一次“接口化”改造：把上面代码里所有直接依赖全局变量的地方清干净，统一用param结构体传递飞机参数。这样无论是从线性气动模型切换到插值表，还是把某型真机的风洞数据换成CFD数据，都只改param，不动方程和仿真框架。我自己的习惯是给param增加版本号字段param.version和一句来源备注，比如“来自XX项目气动报告v2.1”。另外一个值得花半天做的事，是写一个自动检查脚本：输入几个已知状态组合，输出对应气动力和力矩，把结果和手算值做对比，误差超过1%立刻报警。这部分自动化检查看起来不起眼，但它能让你在改参数后放心大胆地继续往下走。

再往下走，很多做任务规划的人会把MATLAB仿真和STK做联合仿真：MATLAB提供飞机六自由度状态，STK提供地形、视线、传感器覆盖范围分析。常见做法是用MATLAB的引擎接口在循环里调用STK场景，把无人机位置每0.1秒写入STK的MotionModel。这里有个坑：STK的接口调用本身有开销，如果每步积分都同步一次，仿真会慢得让你怀疑人生。我的做法是：先在MATLAB里跑完一整段轨迹，最后一次性把轨迹点批量写入STK，让STK只做覆盖分析和可视化。这个“先算完再联合”的模式，比逐帧同步实用得多。

我在这条路上最大的教训，是刚开始总想贪多：一上来就想把失速、螺旋桨滑流、弹性模型全塞进去。结果模型复杂度上去了，bug也跟着上去了，配平越来越难查。后来我改成“先线性再非线性、先刚体再弹性、先单机再编队”的顺序，每一步都留足验证时间，反而走得快。希望帮到你。

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