简介：这份基于MATLAB实现的飞机机动动作轨迹仿真工程，覆盖盘旋、蛇形、眼镜蛇等典型机动动作的建模与可视化，适合飞行器仿真入门者、MATLAB学习者以及需要快速生成飞行轨迹的科研用户。压缩包共5个文件，包括4个m脚本和1份md使用说明文档，m文件里含主函数与多个机动模块，按功能拆分便于调用，文档则说明运行流程，整体包仅16KB，结构精简单方便携。目前已有140人浏览学习，可直接用于课程设计、毕业设计或飞行控制预研。借助说明文档，用户能快速理解各机动模型的组织方式，按需调用盘旋、蛇形、眼镜蛇等模块，替换数据或参数即可运行并输出轨迹效果图；对新手而言，无需深究底层推导也能复现仿真结果，适合与MATLAB基础学习结合，边看边练快速积累飞行仿真经验。
1. 基于 MATLAB 实现的飞机机动动作轨迹仿真：盘旋、蛇形与眼镜蛇的建模差异
拿到标题里「基于 MATLAB 实现的飞机机动动作轨迹仿真 + 使用说明文档」这类工程包，最常见的困局是代码能出图，结果却经不起检查：盘旋绕一整圈后起点和终点对不齐，蛇形的横向振幅和理论值差一倍，眼镜蛇拉到接近 90° 时高度曲线直接发散。问题几乎都不在 MATLAB 本身，而在三种机动在运动学上根本不是同一类问题。盘旋是常值坡度下的定常转弯，误差在长时间积分里慢慢累积；蛇形是周期控制输入驱动的跟踪问题，振幅由坡度幅值和频率共同决定；眼镜蛇则要处理大迎角下升力、阻力与速度的强耦合，点质量模型里必须加进迎角到气动系数的映射，才能做出速度骤减、高度近似不变的效果。下面沿「状态定义 → 微分方程 → ODE45 积分 → 可视化与参数扫描 → 解析校验」这条链，把三种机动放进同一个 MATLAB 框架实现，并给出参数表、数值坑和扩展方向。
2. 盘旋、蛇形与眼镜蛇机动的运动学模型与参数映射
2.1 坐标系与状态量定义：NEU 地面系下的六维状态
先定坐标系，否则后面所有公式都会串位。仿真里取北东上（NEU）地面系：x 指北，y 指东，h 是海拔高度（向上为正）。飞机位置由 (x, y, h) 描述，V 为真空速，γ 是航迹倾角（爬升为正），χ 是航向角（自北顺时针量取），六维状态向量写成
s = [x, y, h, V, γ, χ]ᵀ
对应运动学微分方程为
ẋ = V cosγ cosχ
ẏ = V cosγ sinχ
ḣ = V sinγ
V̇ = (T − D)/m − g sinγ
γ̇ = (L cosφ − m g cosγ) / (m V)
χ̇ = L sinφ / (m V cosγ)
其中 φ 是坡度，T 是推力，D 是阻力，L 是升力。这套点质量三自由度方程的控制量是坡度和升力（或迎角），比 6DOF 少了三个姿态角速率方程，但已经能正确描述航迹几何和时间历程。盘旋和蛇形只在水平面内讨论，γ 恒为 0，L 与 m g 平衡，χ̇ 简化成 g tanφ/V；眼镜蛇必须保留 V̇ 和 γ̇，并把 L、D 写成迎角的函数。状态量全部用 SI 单位、角度用弧度，是后面所有代码能对齐的前提。
2.2 盘旋机动：定常转弯的半径、坡度与角速度约束
水平盘旋即协调转弯：L cosφ = m g，γ = 0，由 χ̇ = g tanφ / V 得到
转弯半径 R = V / χ̇ = V² / (g tanφ)
转弯角速度 ω = g tanφ / V
半径与速度平方成正比、与坡度正切成反比，这意味着调速度对半径的影响远大于调坡度。常见做法有两种：直接给常值航向角速率，或者给常值坡度加速度。我一般用后者，因为改坡度即可直观看到半径变化，且与真实飞行的协调转弯约束一致。以下是一组可复现的对照数据。
| 真空速 V (m/s) | 坡度 φ | 理论半径 R (m) | 角速度 ω (°/s) | 转一整圈时间 (s) |
|---|---|---|---|---|
| 100 | 20° | 2800 | 2.05 | 176 |
| 100 | 30° | 1766 | 3.25 | 111 |
| 150 | 30° | 3973 | 2.16 | 166 |
| 200 | 45° | 4077 | 2.81 | 128 |
注意转一整圈的时间是 2πR/V，仿真时长必须按这个公式对齐，否则画出来的圆永远不会闭合。若想更真实，可以在盘旋入口加一个坡度建立过程 φ(t) = φ_cmd(1 − e^(−t/τ))，τ 取 2 s 左右，入口瞬态就不那么突兀。
2.3 蛇形机动：正弦坡度指令驱动下的周期地面轨迹
蛇形机动（serpentine maneuver）的俯视航迹是一条周期曲线，建模方式是让坡度按正弦变化
φ(t) = φ_max sin(ωt)
代入合作转弯关系 χ̇ = (g/V) tan(φ_max sinωt)。在小坡度、小角度假设下 tanφ ≈ φ、sinχ ≈ χ，可以得到解析解
y(t) = y0 − (g φ_max / ω²) sin(ωt)
横向振幅与波长分别为
A = g φ_max / ω²
λ = 2π V / ω
这里有个反直觉的结论：蛇形振幅不显含速度，速度只决定波长；想拉开横向幅度，降 ω 是最有效的，因为振幅按 1/ω² 增长。若改用横向加速度指令 a_y = a_{y,max} sinωt，则 A = a_{y,max}/ω²，与坡度写法等价，因为 a_y = g tanφ。下表给出三组典型参数。
| 速度 V (m/s) | 最大坡度 φ_max | 角频率 ω (rad/s) | 理论振幅 A (m) | 波长 λ (m) |
|---|---|---|---|---|
| 100 | 25° | 0.1 | 428 | 6283 |
| 100 | 25° | 0.2 | 107 | 3142 |
| 100 | 25° | 0.3 | 48 | 2094 |
| 150 | 25° | 0.2 | 107 | 4712 |
大坡度时解析式会偏离数值解，因为 tanφ 比线性更快增长，所以仿真代码里必须用完整的 g tan(φ)/V，而不是直接把 φ 当弧度塞进线性公式。
2.4 眼镜蛇机动：大迎角下拉力、阻力与速度的强耦合
眼镜蛇机动的本质是：高度基本不变，飞机快速拉起，迎角越过失速区到 90°–110°，随后推杆恢复，速度从进入值掉到一半以下。轨迹仿真关心的耦合在两点：迎角 α = θ − γ 同时进入升力和阻力公式，阻力剧增使 V 下降，V 下降又改变动压、拉低气动力。右端项必须同时保留 V̇ 和 γ̇。
点质量模型用平板块气动模型就能抓住主要特征
C_L = sin(2α)
C_D = C_D0 + 2 sin²α
α = 45° 时升力最大，α = 90° 时升力为零、阻力最大，正好对应眼镜蛇末尾「速度骤减、轨迹平缓」的现象。姿态侧用如下平滑俯仰剖面，该函数的 MATLAB 形式在 3.4 节直接复用。
function theta = cobra_pitch_profile(t, p) % p.t_pull 拉杆时间, p.t_hold 保持时间, p.theta_max 最大俯仰角 if t <= p.t_pull theta = p.theta_max * 0.5 * (1 - cos(pi * t / p.t_pull)); elseif t <= p.t_pull + p.t_hold theta = p.theta_max; else tt = min(t - p.t_pull - p.t_hold, p.t_pull); theta = p.theta_max * 0.5 * (1 + cos(pi * tt / p.t_pull)); end end这个剖面由拉杆、保持、推杆三段组成，余弦过渡保证俯仰角连续，平均俯仰速率约等于 θ_max / t_pull。θ_max = 110°、t_pull = 2 s 时约 55°/s，与真实眼镜蛇的拉起速率同量级。需要说明适用边界：点质量模型能看到高度、速度、航迹倾角和过载历程，但看不到机头转动的惯性，真实的俯仰振荡和迎角过冲要刚体 6DOF 加俯仰力矩系数才能复现。
三种机动的参数映射汇总如下，这一表也建议直接写进工程包的使用说明文档。
| 机动 | 控制输入 | 核心关系 | 轨迹特征 |
|---|---|---|---|
| 盘旋 | 常值坡度 φ | R = V²/(g tanφ) | 水平圆，半径由 V、φ 决定 |
| 蛇形 | 正弦坡度 φ_max sinωt | A = g φ_max / ω² | 周期横向偏移，波长 2πV/ω |
| 眼镜蛇 | 俯仰剖面 θ(t) | α = θ − γ，C_L、C_D(α) | 高度近似不变，速度骤降 |
3. 用 MATLAB 搭建三种机动轨迹仿真器：ODE45 与统一入口函数
3.1 顶层统一入口：mode 参数分工的 simulate_maneuver
把所有工况放进一个函数，用 mode 字符串分发，使用说明文档里只需要约定一个 API 和一份参数结构体。状态向量、单位、输出字段全在注释里写死，避免换人维护时把速度和角度单位搞混。
function res = simulate_maneuver(mode, params, s0) % mode - 'circle' / 'serpentine' / 'cobra' % params - 结构体, 字段按工况在调用方设置 % s0 - 6 维列向量 [x0; y0; h0; V0; gamma0; chi0] % 长度单位 m, 速度单位 m/s, 角度单位 rad tspan = 0:params.dt:params.t_end; opt = odeset('RelTol', 1e-9, 'AbsTol', 1e-9, 'MaxStep', 1.0); switch mode case 'circle' [t, s] = ode45(@(t,s) rhs_circle(t, s, params), tspan, s0, opt); case 'serpentine' [t, s] = ode45(@(t,s) rhs_serpentine(t, s, params), tspan, s0, opt); case 'cobra' [t, s] = ode45(@(t,s) rhs_cobra(t, s, params), tspan, s0, opt); otherwise error('未知机动模式 %s, 请检查 mode 字符串', mode); end res.t = t; % 时间列向量 res.s = s; % 状态矩阵, 各列对应 x y h V gamma chi end逻辑说明：tspan 用 0:dt:t_end 生成，这里 dt 只决定输出点密度，不影响 ode45 内部自适应步长；RelTol 和 AbsTol 都压到 1e-9，是为了让长时间盘旋的闭合误差不累积；MaxStep 1.0 防止眼镜蛇这种强非线性段积分器大步长绕过剧烈变化。params 里每个字段在调用前赋值，mode 写错时 error 直接给出提示。
3.2 盘旋工况：常值坡度右端项与半径对照
右端项是六个分量的导数，唯一有物理量的是航向角速率 χ̇ = g tanφ/V。
function ds = rhs_circle(t, s, p) % 状态解包 x = s(1); y = s(2); h = s(3); V = s(4); gam = s(5); chi = s(6); g = 9.80665; chi_dot = g * tan(p.phi) / V; % 协调转弯, 坡度恒定, 定速定高 ds = zeros(6,1); ds(1) = V * cos(gam) * cos(chi); ds(2) = V * cos(gam) * sin(chi); ds(3) = V * sin(gam); % gam 初值为 0 时高度不变 ds(4) = 0; % 定速假设: 推力等于阻力 ds(5) = 0; % 航迹倾角不变 ds(6) = chi_dot; end参数说明：p.phi 必须是弧度，这是最常见的报错来源；V 取常值是工程简化，真实盘旋需要推力配平。调用驱动脚本时，把 t_end 对齐到整圈时间 2πV/(g tanφ)，再算半径和闭合误差。
s0 = [0; 0; 1000; 100; 0; 0]; % 起点原点, 高度 1 km, 速度 100 m/s p.dt = 0.1; p.V = 100; p.phi = deg2rad(25); p.t_end = 2*pi*p.V / (9.81*tan(p.phi)); % 刚好一整圈 res = simulate_maneuver('circle', p, s0); R_sim = mean(sqrt(res.s(:,1).^2 + res.s(:,2).^2)); R_ana = p.V^2 / (9.81*tan(p.phi)); fprintf('仿真半径 %.2f m, 理论半径 %.2f m\n', R_sim, R_ana);这里 R_sim 用数值解每个时刻到原点的距离取平均，R_ana 用解析式；两者偏差在 0.5% 以内说明右端项正确。闭合误差直接算终点与起点的位置差，留到第 5 章统一校验。
3.3 蛇形工况：正弦坡度指令与相位对齐
蛇形右端项和盘旋几乎一样，只是坡度从常值变成正弦函数，因此要额外传入角频率 p.w 和最大坡度 p.phi_max。
function ds = rhs_serpentine(t, s, p) x = s(1); y = s(2); h = s(3); V = s(4); gam = s(5); chi = s(6); g = 9.80665; phi_cmd = p.phi_max * sin(p.w * t); % 坡度指令为纯正弦 chi_dot = g * tan(phi_cmd) / V; % 合作转弯关系 ds = zeros(6,1); ds(1) = V * cos(gam) * cos(chi); ds(2) = V * cos(gam) * sin(chi); ds(3) = V * sin(gam); % 定高飞行 ds(4) = 0; ds(5) = 0; ds(6) = chi_dot; end参数说明：p.w 的单位是 rad/s，300 rad/s 之类纯属没转换单位；从 t = 0 起 sin 初值为 0、斜率为正，坡度指令连续，不需要额外加 ramp。初始状态把 χ0 设 0，横向偏移就关于 y0 对称。若要验证理论振幅，直接对 y 列取 (max − min)/2，与 g·φ_max/w² 对比，第 5 章给出完整校验代码。
3.4 眼镜蛇工况：平板块气动模型与推重比假设
眼镜蛇右端项最复杂，需要气动系数、动压和重力投影。推力按慢车处理取 0，密度按海平面常数 1.225 kg/m³。
function ds = rhs_cobra(t, s, p) x = s(1); y = s(2); h = s(3); V = s(4); gam = s(5); chi = s(6); g = 9.80665; rho = 1.225; theta = cobra_pitch_profile(t, p); % 俯仰角指令, rad alpha = max(-pi/2, min(pi/2, theta-gam));% 迎角限幅到 [-90, 90] CL = sin(2*alpha); % 平板块升力模型 CD = p.CD0 + 2*sin(alpha)^2; % 大迎角阻力模型 qbar = 0.5 * rho * V^2; % 动压 L = qbar * p.S * CL; D = qbar * p.S * CD; ds = zeros(6,1); ds(1) = V * cos(gam) * cos(chi); ds(2) = V * cos(gam) * sin(chi); ds(3) = V * sin(gam); ds(4) = (-D)/p.m - g*sin(gam); % 推力取 0 ds(5) = (L - p.m*g*cos(gam))/(p.m*V); % V 极小时需保护, 见说明 ds(6) = 0; % 保持进入航向 end参数说明：迎角限幅到 ±90° 是工程简化，真实眼镜蛇会略过垂直再回来，点质量模型在越过垂直的区间没有俯仰惯性支撑，限幅反而让速度、高度的时间历程更稳。ds(5) 里 V 在分母，仿真中速度会从 110 m/s 掉到 50 m/s 左右，不会趋零；若改了参数导致 V 接近 0，应在调用处加if V < 1e-3的保护分支。驱动参数见下表。
| 参数 | 取值 | 物理含义 |
|---|---|---|
| p.m | 15000 kg | 飞机质量 |
| p.S | 50 m² | 参考面积 |
| p.CD0 | 0.02 | 零升阻力系数 |
| p.theta_max | 110° | 最大俯仰角 |
| p.t_pull / p.t_hold | 2 s / 1 s | 拉杆与保持时间 |
| p.t_end | 8 s | 仿真时长 |
这套参数跑出来，速度曲线会在 2–3 s 内显著下降，高度变化不超过 50 m，符合眼镜蛇「先掉速、不掉高」的定性特征。注意真实机型需要查气动数据库，表里的数值只用于演示运动学结构。
4. 跑通三种机动算例：MATLAB 可视化与参数影响
4.1 三维轨迹与四联图绘制脚本
统一入口的好处是绘图脚本也只需要写一份。用 tiledlayout 把三维轨迹、水平投影、高度和速度时间历程放在同一张图里，水平投影必须加 axis equal，否则圆会被拉伸成椭圆，这是盘旋仿真里最常见的主观误判。
figure('Color','w'); tiledlayout(2,2,'TileSpacing','compact'); nexttile; plot3(res.s(:,1), res.s(:,2), res.s(:,3), 'b', 'LineWidth', 1.2); grid on; axis equal; xlabel('北 x (m)'); ylabel('东 y (m)'); zlabel('高度 h (m)'); title('三维轨迹'); view(-30, 25); nexttile; plot(res.s(:,1), res.s(:,2), 'b'); grid on; axis equal; xlabel('北 x (m)'); ylabel('东 y (m)'); title('水平投影'); nexttile; plot(res.t, res.s(:,3)); grid on; xlabel('时间 t (s)'); ylabel('h (m)'); title('高度时间历程'); nexttile; plot(res.t, res.s(:,4)); grid on; xlabel('时间 t (s)'); ylabel('V (m/s)'); title('速度时间历程');逻辑说明：plot3 第一二三个参数分别是 x、y、h 列；view(-30, 25) 给出斜视角，方便看盘旋的圆形和眼镜蛇的高度保持。高度和速度子图用于检查蛇形是否定高、眼镜蛇速度是否如预期骤降。
4.2 参数扫描：坡度、角频率与拉起时间的敏感度
三种机动各有最敏感的参数：盘旋是速度（半径按 V² 变），蛇形是角频率（振幅按 1/ω² 收缩），眼镜蛇是最大俯仰角（决定峰值阻力和速度损失）。写一个循环扫描坡度，观察半径变化：
phi_list = deg2rad([15 20 25 30 35]); figure('Color','w'); hold on; for k = 1:numel(phi_list) p.phi = phi_list(k); p.t_end = 2*pi*p.V / (9.81*tan(p.phi)); % 每圈时长随坡度变化 res = simulate_maneuver('circle', p, s0); plot(res.s(:,1), res.s(:,2), 'LineWidth', 1.2); end hold off; grid on; axis equal; legend(arrayfun(@(d) sprintf('坡度 %d°', rad2deg(d)), phi_list, 'UniformOutput', false)); xlabel('北 x (m)'); ylabel('东 y (m)'); title('坡度对盘旋半径的影响');这里的重点是 p.t_end 必须跟随 φ 重新计算，否则小坡度圆还没画完、大坡度圆已经超过一圈，对比图会误导。参数影响规律整理如下。
| 机动 | 参数变化 | 轨迹响应 |
|---|---|---|
| 盘旋 | φ 增大 | 半径按 1/tanφ 减小，周期缩短 |
| 盘旋 | V 增大 | 半径按 V² 增大，比调坡度灵敏得多 |
| 蛇形 | ω 增大 | 振幅按 1/ω² 收缩，波长按 1/ω 缩短 |
| 蛇形 | φ_max 增大 | 振幅线性增大，大坡度时解析式偏差变大 |
| 眼镜蛇 | θ_max 增大 | 峰值阻力增大，速度损失更多 |
| 眼镜蛇 | t_pull 增大 | 平均俯仰速率降低，过载减小，速度下降变缓 |
4.3 三个高发数值坑与对应处理
第一个是度与弧度混用。参数表里写 25，代码里直接 tan(25)，得到的半径和理论值差一个数量级，报错形态是「仿真半径 780 m，理论 2800 m」。统一用 deg2rad 包装输入，输出角度如需显示再 rad2deg。
第二个是 ODE45 默认容差导致长时间积分漂移。盘旋转 176 s，默认 RelTol 1e-3 下位置误差可能累积到数十米，圆末端出现可见的开口。把 RelTol、AbsTol 压到 1e-9，闭合误差通常在 1e-6 m 量级。代价是积分步数增加，对单圈仿真完全可接受。
第三个是欧拉角奇异性。点质量模型的状态量是 γ 和 χ，不是机体姿态欧拉角，所以正常仿真不会触发万向锁；但只要你想输出俯仰角 θ 并把它画成时间历程，就会遇到 θ 接近 90° 的问题。处理办法有两个：模型内用四元数传播姿态，或者绘图时用 atan2 类函数避免主值跳变。第 5 章给出四元数扩展的最小代码。
5. 轨迹仿真结果的解析校验与两步扩展技巧
5.1 用解析解闭环校验盘旋、蛇形与眼镜蛇
仿真的意义在于可复现，所以每一步都要有独立于数值积分的判据。盘旋看整圈航向角变化和位置漂移，蛇形看振幅，眼镜蛇看能量守恒。以下代码直接贴在主脚本尾部即可。
% 盘旋: 航向角应恰好转 2*pi, 位置漂移应在 1e-6 量级 omega = 9.81*tan(p.phi)/p.V; drift = norm(res.s(end,1:2)' - res.s(1,1:2)'); fprintf('盘旋 航向角变化 %.6f rad, 位置漂移 %.3e m, 半径 %.2f m\n', ... omega*res.t(end), drift, R_ana); % 蛇形: 仿真振幅对比解析振幅 A = g*phi_max/w^2 p.w = 0.2; p.phi_max = deg2rad(25); p.t_end = 60; p.dt = 0.1; res = simulate_maneuver('serpentine', p, s0); A_sim = (max(res.s(:,2)) - min(res.s(:,2))) / 2; A_ana = 9.81 * p.phi_max / p.w^2; fprintf('蛇形 仿真振幅 %.2f m, 解析振幅 %.2f m, 偏差 %.2f%%\n', ... A_sim, A_ana, abs(A_sim-A_ana)/A_ana*100); % 眼镜蛇: 动能损失应远大于势能变化, 且总机械能损失等于阻力做功 res = simulate_maneuver('cobra', p, s0); dEk = 0.5*p.m*(res.s(1,4)^2 - res.s(end,4)^2); dEp = p.m*9.81*(res.s(end,3) - res.s(1,3)); fprintf('眼镜蛇 动能损失 %.1f kJ, 势能变化 %.1f kJ\n', dEk/1e3, dEp/1e3);逻辑说明：蛇形振幅的解析式是线性近似，仿真用完整 tan 函数，所以偏差随 φ_max 增大而增大，但 φ_max 在 30° 以内时偏差通常不超过 3%，超过这个范围就该以数值解为准。眼镜蛇的能量校验收尾时再补一步数值积分计算阻力做功 ∫D·V dt，与 dEk + dEp 对比，闭合误差在 1% 内说明气动模型和积分器没有系统性错误。
5.2 扩展技巧：四元数姿态传播与 6DOF 升级路径
如果下一步要把点质量轨迹升级成带姿态的 6DOF 仿真，最直接的动作是把俯仰角剖面换成四元数传播，绕开欧拉角在 ±90° 俯仰时的奇异问题。最小实现如下。
% 四元数 q = [qw qx qy qz], omega_body 是机体轴角速度 (rad/s) function qd = quat_rate(q, omega_body) w = q(1); v = q(2:4); qd = 0.5 * [ -v' * omega_body; w * omega_body + cross(v, omega_body) ]; end参数说明：这个函数对应四元数微分方程 q̇ = 0.5 q ⊗ ω，omega_body 是机体坐标系下的角速度向量。升级到 6DOF 时，给俯仰角速度设限、加进俯仰力矩方程和 Cm(α) 曲线，即可看到真实眼镜蛇的迎角过冲和机头抖动。四元数转欧拉角再交给绘制接口处理，姿态显示仍用 rad2deg 输出，但内部积分永远不碰奇异点。
另一个工程化技巧是把三种机动的参数结构体统一写成函数返回，并生成一份 CSV 参数表。这样使用说明文档里的参数表、MATLAB 里的 params 结构体、仿真输出三处永远一致，改参数不会出现「文档改了、代码没改」的经典事故。校验脚本和参数表放在同一个目录，别人拿到工程包后跑一遍主脚本，三种机动就能全部复现并自动打印误差，比对着文档逐条核对可靠得多。
本文还有配套的精品资源，点击获取
