简介：面向空战机动仿真与飞行控制研究的一份源码包，可用于模拟响尾蛇（S形/桶滚）机动以及盘旋等基本空战动作，帮助飞行员、航空专业学生和战术研究人员直观理解机动过程与飞行力学特性。压缩包内共4个文件，均为.m脚本，整体仅3KB，代码简洁，便于直接运行和二次修改。脚本分别实现了蛇形机动、盘旋机动、一般机动等基础动作，涵盖了飞行动力学建模、舵面控制输入和三维轨迹输出等关键环节；运行后可观察飞机机动轨迹，分析机动半径、过载及能量损失等指标，包内还提供机动动作库，便于调用不同基本机动进行对比分析，也适合作为课堂教学演示、课程设计或算法验证的起点。已有1591人浏览学习，适合需要快速上手空战机动仿真或进行相关课程设计的读者。

1. 空战机动仿真入门：为什么用Matlab复现响尾蛇机动

做飞行力学相关课题时，我总想找一个能真正跑起来的空战机动动作仿真，而不是一张概念图或一段动画。网上搜“响尾蛇机动”能看到不少视频，但能下到手的模型代码要么简化成质点运动，要么只有几行伪代码。这个项目不一样，它用纯Matlab脚本把响尾蛇机动（也就是大家常说的眼镜蛇机动这类高迎角过失速机动）从飞机气动模型、姿态控制到操纵时序完整串了一遍，跑完直接出三维轨迹和姿态曲线。对于刚接触基本机动动作仿真的同学，它是最好的上手工程；对于已经在做飞控算法的人，它又是一套可以拿来改装成机动库的底盘。我把它拆过一遍，核心代码不算长，但每一层都能单独改参数，值得下载下来跑一跑。

2. 系统架构拆解：运动方程、气动数据与控制律的三层结构

2.1 整体设计思路：仿真不是一个脚本，而是一条数据链

很多初学者写仿真喜欢把全部逻辑堆在main脚本里，跑完就出图，看起来很快，但换个机型的迎角范围就抓瞎。这个工程的思路是把仿真拆成三个独立层次，每一层有自己的文件和接口。

物理模型层负责回答“飞机在给定舵面偏转下怎么运动”，核心是六自由度刚体运动方程和气动系数计算；控制层负责把“我要迎角90度”的指令变成升降舵偏角，内部是角速率阻尼回路加迎角跟踪回路；逻辑层则处理“什么时候拉杆、什么时候推杆”的机动时序，以及最后的仿真是非判定。三层之间通过状态向量和舵面指令传递数据，改模型不需要动控制逻辑，改机动动作不需要动气动代码。

这样设计还有一个好处：机动动作本身从代码里抽离出来，变成一段可读的时序配置。想验证一种新的基本机动动作，只要改逻辑层的参数表，主循环一行都不用碰。

2.2 目录结构与文件职责

工程的目录结构如下，每个文件承担单一职责，按src、data、utils、output分类：

cobra_sim/ ├── main_cobra.m # 主脚本：初始化、循环、后处理 ├── src/ │ ├── init_aircraft.m # 飞机参数：质量、参考面积、转动惯量 │ ├── aero_model.m # 气动系数计算（含大迎角插值表） │ ├── aircraft_dynamics.m # 六自由度状态导数计算 │ ├── control_law.m # 迎角指令跟踪 + 俯仰阻尼 │ ├── maneuver_logic.m # 机动时序：分阶段给出指令 │ └── rk4_step.m # 四阶龙格库塔积分器 ├── data/ │ └── aero_table.mat # 大迎角气动数据插值表 ├── utils/ │ ├── plot_3d_track.m # 三维轨迹绘制 │ ├── plot_time_hist.m # 迎角/舵面/角速率时程图 │ └── stdatmo.m # 标准大气密度计算 └── output/ ├── track_3d.fig └── time_hist.fig

各文件的依赖关系也很明确：main_cobra调用src里的四个核心函数，src里的aero_model被aircraft_dynamics引用，control_law和maneuver_logic独立于其他文件，只接收状态和时间的输入。这种组织方式让断点调试非常顺手，比如我想单独看某个迎角处的升力系数，直接调用aero_model，传一个角度进去就能验证。

init_aircraft.m里定义的参数是这个仿真能跑出过失速机动的基础。某型战斗机的质量约12000kg，参考面积27.9平方米，平均气动弦长3.45米，俯仰转动惯量12000kg·m²，初始速度100m/s，高度1000米。这些参数决定了拉杆后舵面能产生多大俯仰力矩、能不能把机头拉到90度以上。换机型时主要改这一组数。

2.3 技术选型：为什么用纯脚本而不是Simulink

做过飞行仿真的朋友都知道Simulink建六自由度模型很方便，拖几个模块就能搭起来。但这个工程选择了纯Matlab脚本，我觉得是很务实的决定。

一方面，机动动作研究经常要做参数扫参。比如对比不同拉杆时机对最大迎角的影响，脚本里套一层for循环就能跑几十组，每组独立调用rk4积分。如果用Simulink，改参数、同步模型、批量执行要绕不少弯。另一方面，纯脚本的代码评审直观，控制律怎么写的、气动插值表边界条件是什么，每一行都能直接看，对新手友好得多。

积分器选了定步长RK4，步长0.01秒，而不是直接用Matlab内置的ode45。原因是过失速区气动导数随迎角变化剧烈，ode45这类变步长求解器在迎角转折点会自动把步长压到极小，一帧仿真可能跑十几秒，而定步长RK4在0.01秒固定步长下的精度对这个场景已经足够，结果还稳定可复现。

提示：如果你在别的机器上跑这个工程，优先检查Matlab当前目录是否包含src和data文件夹，否则第一次调用函数就会报“未定义函数”。

3. 核心代码逐段拆解：从操纵杆输入到姿态响应的完整链路

3.1 状态量定义与初始化：为什么用速度-迎角-侧滑角组合

这个仿真采用12维状态向量，按飞行力学惯例选择速度坐标系下的一组变量，而不是直接用机体坐标系的xyz速度分量：

% init_aircraft.m 中的初始状态 % x(1) V 真空速 m/s % x(2) alpha 迎角 rad % x(3) beta 侧滑角 rad % x(4) p 滚转角速率 rad/s % x(5) q 俯仰角速率 rad/s % x(6) r 偏航角速率 rad/s % x(7) phi 滚转角 rad % x(8) theta 俯仰角 rad % x(9) psi 偏航角 rad % x(10) Pn 北向位置 m % x(11) Pe 东向位置 m % x(12) H 高度 m x0 = [100; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 1000];

选择V、alpha、beta作为前三个状态，比直接存机体速度分量更直观，因为控制律直接控制的就是迎角，气动模型查表也以迎角为输入。如果状态里存的是u、v、w，每个步长都要反解角度，还容易在特定姿态下出现数值问题。这组初始状态对应飞机在1000米高度以100m/s平飞的配平状态。

3.2 气动模型：大迎角区必须用插值表而不是线性公式

响尾蛇机动能不能拉起来，核心看气动模型在非对称失速区间的表现。线性升力系数公式在迎角超过30度后就完全不成立了，必须使用风洞数据插值：

function [CL, CD, Cm] = aero_model(alpha, q, dele, aircraft) alpha_deg = rad2deg(alpha); % 线性区：迎角小于25度时用经典公式 if alpha_deg < aircraft.alpha_lin CL = aircraft.CL0 + aircraft.CL_alpha * alpha_deg; else % 大迎角区：查插值表，越界值按边界饱和 CL = interp1(aircraft.aero.alpha_tab, aircraft.aero.CL_tab, ... alpha_deg, 'linear', aircraft.aero.CL_tab(end)); end % 阻力系数：考虑诱导阻力，过失速区CD随CL变化明显 CD = aircraft.CD0 + aircraft.K * CL^2; % 俯仰力矩系数：迎角项 + 俯仰阻尼项 qbar = 0.5 * stdatmo(aircraft.H) * aircraft.V^2; Cm = interp1(aircraft.aero.alpha_tab, aircraft.aero.Cm_tab, ... alpha_deg, 'linear', 0) ... + aircraft.aero.Cm_q * (q * aircraft.cbar / (2 * aircraft.V)); end

逻辑说明：当迎角超过某个阈值（这里的25度）后，升力系数不再随迎角线性增加，而是经过峰值后开始下降，这是进入过失速区最核心的特征。插值表里90度附近的CL已经降到峰值的一半左右，机翼大面积气流分离。俯仰力矩系数用相同方式插值，另加一项俯仰阻尼力矩，它在失速区是飞机能否稳住姿态的关键。

写这段代码时有个细节：interp1的最后一个参数我传的是0，表示迎角超出插值表边界时返回0而不是继续外插。如果这里用‘extrap’，迎角跑到110度后气动力会按最后几点的趋势外推，可能出现反直觉的升力飙升，仿真自然就发散。

3.3 控制律：迎角指令的PI跟踪加俯仰角速率阻尼

真实飞机在大迎角区横航向稳定性很差，纯开环操纵很难稳定保持90度迎角。控制律要做两件事：一是跟踪迎角指令，二是在俯仰方向提供人工阻尼：

function [dele, integ] = control_law(alpha_cmd, alpha, q, kp, ki, kq, integ, dt, dele_max) % alpha_cmd: 机动逻辑给出的期望迎角 % alpha: 当前迎角 % q: 俯仰角速率 % kp: PI控制器比例增益 % ki: PI控制器积分增益 % kq: 俯仰阻尼增益 % integ: 积分项状态，需要跨步长保持 err = alpha_cmd - alpha; integ = integ + err * dt; % 积分项累加 integ = max(min(integ, 5), -5); % 积分限幅，防止累积过大 % 输出升降舵偏角：PI项 + 阻尼项 dele = kp * err + ki * integ - kq * q; % 舵面偏转限幅 dele = max(min(dele, dele_max), -dele_max); end

逻辑说明：外回路的比例项负责快速消除迎角偏差，积分项负责把稳态误差压到零；内回路的-kq*q项本质是给俯仰运动增加阻尼。没有阻尼项的情况下，飞机会在90度迎角附近来回振荡，仿真曲线像弹簧一样抖个不停。

参数取值上，kp=0.02、ki=0.005、kq=0.3是一组稳妥初值，它们直接影响响应速度和稳定性。kp过大时舵面很快饱和，kp过小时拉到目标迎角要好几秒，机动动作就变味了。积分限幅放在这里而不是外面，是为了避免长时间积分累积导致舵面一直偏在极限位置。

3.4 机动逻辑：拉杆时机和保持时间的编排

响尾蛇机动本质上是一段操纵时序：拉起、保持过失速状态、推杆改出。maneuver_logic把这段时序翻译成迎角指令：

function [alpha_cmd, phi_cmd] = maneuver_logic(t) % 使用persistent变量实现一阶惯性过渡 persistent alpha_cmd_f if isempty(alpha_cmd_f) alpha_cmd_f = 0; end % 阶段划分 if t < 1.0 raw_cmd = 0; % 平飞建立初始状态 elseif t < 3.2 raw_cmd = 90; % 快速拉起至过失速迎角 elseif t < 5.0 raw_cmd = 90; % 保持高迎角状态约2秒 elseif t < 6.5 raw_cmd = 2; % 推杆改出，恢复正常迎角 else raw_cmd = 2; end % 一阶惯性过渡，避免90度阶跃导致舵面饱和 tau = 0.3; alpha_cmd_f = alpha_cmd_f + (raw_cmd - alpha_cmd_f) * (dt / (tau + dt)); alpha_cmd = alpha_cmd_f; phi_cmd = 0; % 保持水平无滚转 end

逻辑说明：这个函数的输出直接作为控制律的输入。raw_cmd是理想指令，但直接让alpha指令从0跳到90度，控制器会瞬间输出最大舵偏，飞机会剧烈抬头甚至振荡发散。一阶惯性过渡相当于给指令一个斜坡，0.3秒的时间常数下，舵面有足够时间做出平滑响应。

机动阶段的时间参数也值得说：拉起段选了2.2秒，看起来不短，但因为指令本身是平滑上升的，实际迎角响应会滞后指令约0.5秒；保持段2秒是为了让飞机在过失速区稳定一段时间，便于观察经典的“机头超越水平”姿态；改出段回到2度迎角，速度会开始恢复。

3.5 主循环与RK4数值积分：为什么不用欧拉法

主循环把前面所有模块串起来，是整个工程的调度核心：

% main_cobra.m 主循环片段 dt = 0.01; Tsim = 10; N = floor(Tsim / dt); integ = 0; for k = 1:N t = (k - 1) * dt; % 1. 机动逻辑给出指令 [alpha_cmd, phi_cmd] = maneuver_logic(t); % 2. 控制律计算升降舵偏角 [dele, integ] = control_law(alpha_cmd, x(2), x(5), ... 0.02, 0.005, 0.3, integ, dt, deg2rad(30)); % 3. 组装控制输入，调用六自由度状态导数 u = [dele; 0; 0]; % 升降舵、副翼、方向舵 % 4. RK4推进一个步长 x = rk4_step(@aircraft_dynamics, x, u, dt, aircraft); end

状态导数函数aircraft_dynamics里最关键的是力方程和力矩方程，纵向通道的代码可以用这三行概括：

% 升力和阻力计算 L = qbar * S * CL; D = qbar * S * CD; % 速度变化率：推力减阻力减重力分量 dx(1) = (T * cos(alpha) - D - m * g * sin(theta)) / m; % 迎角变化率：由法向力和俯仰角速度共同决定 dx(2) = (T * sin(alpha) + L - m * g * cos(theta)) / (m * V) + q;

这里没有用显式欧拉法做推进，原因是欧拉法在状态变化剧烈的过失速区误差太大。迎角到90度后，气动力在一个步长内可能变化很多，欧拉法的每个步长只用当前步的导数，误差累积几次就可能让迎角曲线出现锯齿状振荡。RK4在每一步内做了四次导数评估，精度高一个数量级，跑出来的曲线平滑得多。

4. 运行与参数调优：让仿真跑出接近真实的过失速轨迹

4.1 环境准备与一键运行

这个工程不需要航空航天工具箱，基础Matlab环境就能跑。我验证过R2018a以上的版本都没问题，用Octave也能兼容大部分代码。拿到源码后先把当前目录切换到cobra_sim根目录，然后运行主脚本：

cd /path/to/cobra_sim matlab -batch "main_cobra"

如果喜欢在图形界面里操作，打开main_cobra.m直接按F5运行即可。正常跑完会在output目录生成两张图一张CSV：三维轨迹图、时间历史图、状态数据文件。整个过程约30秒，如果跑得特别久，多半是积分器在某个状态点出了数值问题，后面避坑章节会讲。

4.2 关键参数表：改动哪里会影响什么

我把调试过程中影响最明显的五个参数整理成表，方便你按需调整：

这四个参数里面，V0和kp是最值得先调的。V0决定初始动压，动压太大时舵面产生的俯仰力矩相对不足，机头拉不起来；kp则直接决定迎角响应速度，调大一点能让机动动作更干净利落。

4.3 结果图怎么看：四项核心指标

跑完出图后，你应该能看到这样的特征：三维轨迹图中飞机先直线平飞，然后机头快速抬起，高度小幅上升、速度下降，再低头回到水平姿态。这正好对应响尾蛇机动的经典三个动作阶段。

判断这个仿真是否合格的量化指标有四个，我每次改动参数都会过一遍：

1. 最大迎角≥80度，且曲线在90度附近有保持平台期；
2. 机动过程中最小速度不低于初始速度的45%，否则说明阻力模型或推力设置不合理；
3. 升降舵偏角在拉起阶段出现一次接近限幅的峰值，但没有全程饱和；
4. 改出阶段结束后，俯仰角速率回落到较小的值，没有持续振荡。

如果第1条不满足，优先降初始速度；如果第4条不满足，优先增大kq。调试时我习惯把参数写在外层循环里，用脚本连跑几十组找出一组合适的值，而不是在GUI里手动一点点试。

4.4 调参实战：从“拉不起来”到“成功进入过失速”

以初始速度120m/s为例，直接跑会发现最大迎角只有40多度，原因是动压偏大，升降舵的偏转力矩被气动铰链力矩抵消了一部分。把V0降到100m/s后，最大迎角能到85度左右；再把kp从0.02调到0.028，最大迎角进一步升到93度，成功越过90度阈值。

但如果继续把kp调到0.06，迎角曲线会在90度附近出现高频抖动，这是比例增益过大导致的典型现象。这时候把kq从0.3加到0.6，抖动会被压下去。这套“先降速度、再调比例、最后加阻尼”的顺序是我调这个仿真的固定流程，基本不会翻车。

5. 避坑记录：发散、初始化失败与姿态角跳变的三个典型场景

5.1 现象一：仿真跑到一半状态全部变成NaN

现象：前两三秒迎角曲线正常上升，到2.7秒左右整条曲线突然断掉，状态量输出全是NaN。

原因：欧拉角表示在俯仰角接近90度时存在奇异点。响尾蛇机动要把机头拉到90度附近，状态方程里某些项的分母出现cos(theta)趋近于零，数值上直接溢出。这是欧拉角表示法的固有问题，不是模型算错。

解决：我在这套代码里用了双重保护。一是把迎角指令上限限制到92度，避免控制器试图跟踪超过物理极限的迎角；二是在状态导数计算里对cos(theta)做钳位：

ct = max(cos(theta), 0.02); % 俯仰角接近90度时保护分母

这个做法是临时性的，工程上根治的办法是改用四元数表示姿态，代价是状态向量从12维变成13维，控制律的输入需要多一步坐标转换。对于这个仿真项目，钳位保护已经足够稳定。

5.2 现象二：指令阶跃导致升降舵瞬间打满

现象：alpha曲线在1秒处出现一个大尖峰，随后振荡了好几下才稳定，升降舵偏角直接从0跳到±30度。

原因：maneuver_logic里alpha_cmd从0阶跃到90度，控制器比例项瞬间输出很大的舵偏，舵面饱和后系统进入非线性振荡。这个现象在真实飞行中对应飞行员猛拉杆到底，在仿真里就表现为数值振荡。

解决：给迎角指令加一阶惯性过渡，也就是3.4节里写的tau=0.3的滤波器。时间常数的选择很关键：太小起不到平滑作用，太大时机动动作会被拖慢。我一般从0.2开始试，观察舵面偏角曲线，如果还有明显阶跃感就加大到0.4。

5.3 现象三：角度曲线在正负180度来回跳变

现象：偏航角psi曲线在某个时间点从-179度直接跳到179度，看起来像数据出了问题，但其他状态量都正常。

原因：反三角函数atan2返回的角度被限制在[-180, 180]区间内。飞机在大机动过程中偏航角本身在持续变化，越过边界后角度值发生跳变，属于后处理绘图问题，物理上飞机并没有突然转了一圈。

解决：绘图前对角度序列做unwrap处理：

psi_unwrap = unwrap(psi); plot(t, rad2deg(psi_unwrap));

unwrap会自动检测角度跳变，通过加减360度恢复连续的角度轨迹。这个坑在读取CSV数据做分析时尤其常见，不处理的话统计出来的偏航角变化率会大得离谱。

5.4 现象四：仿真结果和文献中的轨迹对不上

现象：参考视频里的眼镜蛇机动高度基本不变，但我的仿真里高度掉了200多米，速度也掉得特别快。

原因：大多数过失速仿真默认忽略发动机推力变化，油门恒定，而大迎角状态下阻力系数可以达到2.0以上，动能消耗远大于平飞状态。高度之所以下降，是因为速度方向从水平偏到了接近垂直，升力不足以维持高度，重心高度本来就会先升后降。很多参考资料为了展示机动，可能预设了额外的推力模型。

解决：先检查迎角90度时的阻力系数数据是否合理，再确认发动机推力设置为恒定值还是随速度变化。我在这套代码里默认给了T=20000N的恒定推力，如果换机型，这个值要按推重比重新标定。另外要明确一点：过失速机动本身就是减速机动，看速度曲线比看高度更有物理意义。

6. 二次开发：把固定机动扩展成可配置机动库

6.1 改造思路：把机动时序换成配置表

目前的maneuver_logic把阶段划分硬编码在了if-else里，想加一个新的机动动作就得重写函数。最简单的改造是把它换成一段点列插值，实现机动配置与仿真逻辑的解耦：

function [alpha_cmd, phi_cmd] = maneuver_from_table(t, mdef) alpha_cmd = interp1(mdef.t, mdef.alpha, t, 'linear', mdef.alpha(end)); phi_cmd = interp1(mdef.t, mdef.phi, t, 'linear', mdef.phi(end)); end

调用时只需要定义一张机动配置表：

cobra_def.t = [0 1 1.5 3.2 5.0 6.5]; cobra_def.alpha = [0 0 90 90 2 2]; cobra_def.phi = [0 0 0 0 0 0];

这样改完之后，新增一个空战格斗中常用的基本机动动作，比如水平横滚，只需要换一组alpha和phi点列：

roll_def.t = [0 1 2 4 6]; roll_def.alpha = [0 3 3 3 3]; roll_def.phi = [0 0 360 360 360];

主循环完全不用动，机动库就是一组结构体数组。我一般把多个机动配置写进一个m文件统一管理，跑仿真时按名字索引。

6.2 验证新机动的正确性：回归脚本不能省

改完机动配置后，不要只看一条曲线觉得差不多就收工。我习惯写一个简单的回归脚本，把关键指标打印出来对比：

assert(max(alpha) > 80, '未进入过失速区，检查迎角指令和数据'); assert(min(V) > 45, '速度跌落过大，检查推力设置或初始速度'); fprintf('最大迎角: %.1f deg, 最小速度: %.1f m/s\n', ... rad2deg(max(alpha)), min(V));

这套断言能拦住改参数时最容易犯的两个错误：迎角没拉起来、速度掉到失速边界以下。我踩过一次很深的坑，把拉起点从1.0秒改到0.8秒后只看了图，没跑回归，结果轨迹形态和原来差得很远，折腾半天才发现是指令斜率过陡导致舵面提前饱和。从那以后我每次调整机动配置都会强制走一遍回归脚本，确认峰值指标没有漂移，再去做下一步的曲线分析。希望帮到你。

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