Stevens《飞行器控制与仿真》读书笔记（一）：飞行动力学建模——从6-DOF到状态空间方程

本文是 Stevens、Lewis 与 Johnson《Aircraft Control and Simulation》（第3版）读书笔记系列的第一篇，覆盖教材第1-2章的核心内容：六自由度运动方程的完整建立与状态空间表达。阅读本文前需掌握微分方程基础（尤其是一阶ODE方程组）。

---

第一篇：飞行动力学建模——从6-DOF到状态空间方程

一、概念引入——飞控设计的三个核心问题

1.1 问题

设计一个无人机的飞行控制系统，需要回答三个层层递进的问题：

1. 系统怎么动？ — 建立被控对象的数学模型（6-DOF 方程）
2. 系统本身稳不稳？ — 分析系统固有动态（模态分析、特征值）
3. 怎么让系统按期望动？ — 设计控制器（PID、状态反馈、制导律）

Stevens 教材的第1-2章回答第一个问题。它的建模方法与已有的《无人机飞行物理学》文章有重叠，但 Stevens 的视角更面向控制——它不满足于”推导出6-DOF方程”，而是进一步将方程整理成控制系统设计需要的标准形式：状态空间方程。

1.2 历史的五分钟

时间	人物	贡献
1687	牛顿	三大运动定律——飞行器动力学的物理基础
1750	欧拉	刚体转动方程
1903	莱特兄弟	首次动力飞行——催生飞行动力学这门学科
1911	布赖恩（Bryan）	《稳定性与操纵性》——首次用线性化方法分析飞机运动
1930s	冯·卡门	提出”小扰动假设”——将非线性的6-DOF方程在平衡点线性化
1960s	卡尔曼	状态空间方法——控制系统设计的现代框架

---

二、坐标系——飞行动力学的语言

2.1 概念引入

飞行动力学中坐标系是最关键的”基础设施”。同一个物理量在不同坐标系下的表达完全不同——搞错坐标系是飞控开发中最常见的错误之一。

Stevens 教材第1章详细定义了三个坐标系：

坐标系	缩写	原点	x轴	z轴	用途
地轴系	NED	地面某点	北（North）	地心方向（Down）	位置、导航
机体轴系	Body	重心	机头方向	“向下”腹部方向	角速度、力/力矩
稳定轴系	Stability	重心	投影到水平面的机体x轴	垂直于x轴向下	小扰动分析

2.2 NED 地轴系 vs 导航常用的 ENU

注意：航空标准使用 NED（北-东-地），这与机器人/无人机行业中常用的 ENU（东-北-天）不同。

特性	NED (Stevens)	ENU (常见ROS/四旋翼)
x轴	北	东
y轴	东	北
z轴	向下（正方向指向地心）	向上（正方向指向天空）
偏航角	$\psi$ 从北测量	$\psi$ 从东测量

如果你从 Stevens 教材转向 PX4/ROS 开发，需要注意这个差异带来的航向角符号变化。

2.3 欧拉角

机体系相对于地轴系的姿态用三个欧拉角描述：

角度	符号	范围	正方向	物理含义
偏航角	$\psi$	$[0^\circ, 360^\circ)$	顺时针（从上往下看）	机头指向
俯仰角	$\theta$	$(-90^\circ, 90^\circ)$	抬头	机头上下
滚转角	$\phi$	$(-180^\circ, 180^\circ)$	右侧抬起	机翼倾斜

从地轴系到机体轴的旋转顺序（Stevens 使用的 3-2-1 旋转，即偏航→俯仰→滚转）：

$$\mathbf{R}_{BE} = \mathbf{R}_x(\phi)\mathbf{R}_y(\theta)\mathbf{R}_z(\psi)$$

展开：

$$\mathbf{R}_{BE} = \begin{bmatrix} \cos\theta\cos\psi & \cos\theta\sin\psi & -\sin\theta \\ \sin\phi\sin\theta\cos\psi - \cos\phi\sin\psi & \sin\phi\sin\theta\sin\psi + \cos\phi\cos\psi & \sin\phi\cos\theta \\ \cos\phi\sin\theta\cos\psi + \sin\phi\sin\psi & \cos\phi\sin\theta\sin\psi - \sin\phi\cos\psi & \cos\phi\cos\theta \end{bmatrix}$$

---

三、刚体六自由度运动方程

3.1 平动方程（在机体系下表达）

在机体系下，牛顿第二定律的形式稍微复杂一些——因为机体系本身在旋转：

$$m(\dot{\mathbf{V}}_B + \boldsymbol{\omega}_B \times \mathbf{V}_B) = \mathbf{F}_B + m\mathbf{g}_B$$

其中：

• $\mathbf{V}_B = [u, v, w]^T$ — 空速在机体系的分量
• $\boldsymbol{\omega}_B = [p, q, r]^T$ — 角速度
• $\mathbf{F}_B$ — 气动力+推力在机体系的分量
• $\mathbf{g}_B$ — 重力在机体系的分量

展开为标量方程：

$$\dot{u} = rv - qw + \frac{F_x}{m} + g_x$$ $$\dot{v} = pw - ru + \frac{F_y}{m} + g_y$$ $$\dot{w} = qu - pv + \frac{F_z}{m} + g_z$$

3.2 转动方程（在机体系下表达）

欧拉转动方程：

$$\mathbf{I}\dot{\boldsymbol{\omega}}_B + \boldsymbol{\omega}_B \times (\mathbf{I}\boldsymbol{\omega}_B) = \mathbf{M}_B$$

其中 $\mathbf{I}$ 是惯性张量：

$$\mathbf{I} = \begin{bmatrix} I_x & -I_{xy} & -I_{xz} \\ -I_{xy} & I_y & -I_{yz} \\ -I_{xz} & -I_{yz} & I_z \end{bmatrix}$$

展开（假设飞行器对称 $I_{xy}=I_{yz}=0$，这是大多数固定翼和四旋翼适用的近似）：

$$\dot{p} = \frac{I_{xz}(I_x - I_y + I_z)pq + (I_z - I_y)I_zqr + I_zL + I_{xz}N}{I_xI_z - I_{xz}^2}$$ $$\dot{q} = \frac{(I_z - I_x)pr - I_{xz}(p^2 - r^2)}{I_y} + \frac{M}{I_y}$$ $$\dot{r} = \frac{(I_x - I_y + I_z)I_{xz}qr + (I_x - I_y)I_xpq + I_xN + I_{xz}L}{I_xI_z - I_{xz}^2}$$

其中 $L, M, N$ 是绕 x, y, z 轴的气动力矩。

3.3 姿态运动学

欧拉角的微分方程：

$$\dot{\phi} = p + q\sin\phi\tan\theta + r\cos\phi\tan\theta$$ $$\dot{\theta} = q\cos\phi - r\sin\phi$$ $$\dot{\psi} = q\sin\phi\sec\theta + r\cos\phi\sec\theta$$

注意：$\theta = \pm 90^\circ$ 时 $\sec\theta$ 发散——这就是万向锁问题。实际飞控在姿态解算中使用四元数来避免这个奇异性。

3.4 位置运动学

位置在 NED 坐标系中的变化率：

$$\dot{P}_N = u\cos\theta\cos\psi + v(\sin\phi\sin\theta\cos\psi - \cos\phi\sin\psi) + w(\cos\phi\sin\theta\cos\psi + \sin\phi\sin\psi)$$ $$\dot{P}_E = u\cos\theta\sin\psi + v(\sin\phi\sin\theta\sin\psi + \cos\phi\cos\psi) + w(\cos\phi\sin\theta\sin\psi - \sin\phi\cos\psi)$$ $$\dot{P}_D = -u\sin\theta + v\sin\phi\cos\theta + w\cos\phi\cos\theta$$

---

四、Stevens 的12状态标准形式

4.1 状态向量

综合以上所有方程，Stevens 将飞行器的状态定义为12维向量：

$$\mathbf{x} = [u, v, w, p, q, r, \phi, \theta, \psi, P_N, P_E, P_D]^T$$

写成 $\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, \mathbf{u})$：

$$\begin{cases} \dot{u} = rv - qw - g\sin\theta + \frac{X + T}{m} \\ \dot{v} = pw - ru + g\cos\theta\sin\phi + \frac{Y}{m} \\ \dot{w} = qu - pv + g\cos\theta\cos\phi + \frac{Z}{m} \\ \dot{p} = (\Gamma_1 pq - \Gamma_2 qr) + \frac{1}{2}\rho V_a^2 S b(C_l + \Delta C_l) \\ \dot{q} = (\Gamma_5 pr - \Gamma_6(p^2 - r^2)) + \frac{1}{2}\rho V_a^2 S \bar{c} C_m / I_y \\ \dot{r} = (\Gamma_7 pq - \Gamma_1 qr) + \frac{1}{2}\rho V_a^2 S b(C_n + \Delta C_n) \\ \dot{\phi} = p + q\sin\phi\tan\theta + r\cos\phi\tan\theta \\ \dot{\theta} = q\cos\phi - r\sin\phi \\ \dot{\psi} = q\sin\phi\sec\theta + r\cos\phi\sec\theta \\ \dot{P}_N = u\cos\theta\cos\psi + v(\cdots) + w(\cdots) \\ \dot{P}_E = u\cos\theta\sin\psi + v(\cdots) + w(\cdots) \\ \dot{P}_D = -u\sin\theta + v\sin\phi\cos\theta + w\cos\phi\cos\theta \end{cases}$$

其中 $X, Y, Z$ 是气动力在机体系的分量，$C_l, C_m, C_n$ 是无量纲气动力矩系数。

4.2 输入向量

控制输入通常为：

• 固定翼：$\mathbf{u} = [\delta_e, \delta_a, \delta_r, \delta_T]^T$（升降舵、副翼、方向舵、油门）
• 四旋翼：$\mathbf{u} = [\omega_1^2, \omega_2^2, \omega_3^2, \omega_4^2]^T$（四个电机转速平方，正比于推力）

4.3 气动力与气动力矩

气动力和力矩通过无量纲系数表示：

$$\begin{aligned} X &= \bar{q} S C_X \\ Y &= \bar{q} S C_Y \\ Z &= \bar{q} S C_Z \\ L &= \bar{q} S b C_l \\ M &= \bar{q} S \bar{c} C_m \\ N &= \bar{q} S b C_n \end{aligned}$$

其中 $\bar{q} = \frac{1}{2}\rho V_a^2$ 是动压，$S$ 是机翼面积，$b$ 是展长，$\bar{c}$ 是平均气动弦长。

这些系数是状态变量和控制输入的复杂函数——Stevens 第2章用气动导数（如 $C_{L_\alpha}$、$C_{m_\alpha}$、$C_{m_q}$）来线性化它们。

---

五、气动导数与小扰动线性化

5.1 小扰动假设

Stevens 的核心思想：在平衡飞行状态（trim condition）附近，将非线性的6-DOF方程线性化。

假设平衡状态为 $\mathbf{x}_0, \mathbf{u}_0$，小扰动量为 $\Delta\mathbf{x} = \mathbf{x} - \mathbf{x}_0$。

将 $f$ 在平衡点附近做泰勒展开，取一阶项：

$$\Delta\dot{\mathbf{x}} = \frac{\partial\mathbf{f}}{\partial\mathbf{x}}\bigg|_{\mathbf{x}_0,\mathbf{u}_0} \Delta\mathbf{x} + \frac{\partial\mathbf{f}}{\partial\mathbf{u}}\bigg|_{\mathbf{x}_0,\mathbf{u}_0} \Delta\mathbf{u}$$

得到线性时不变（LTI）状态空间方程：

$$\boxed{\Delta\dot{\mathbf{x}} = \mathbf{A}\Delta\mathbf{x} + \mathbf{B}\Delta\mathbf{u}}$$

其中 $\mathbf{A} \in \mathbb{R}^{12\times12}$ 是系统矩阵，$\mathbf{B} \in \mathbb{R}^{12\times4}$ 是输入矩阵。

5.2 气动导数的参数化

气动力系数用气动导数表示为状态的线性函数。以升力系数为例：

$$C_L = C_{L_0} + C_{L_\alpha}\alpha + C_{L_{\dot{\alpha}}}\dot{\alpha} + C_{L_q} \frac{q\bar{c}}{2V_a} + C_{L_{\delta_e}}\delta_e$$

• $C_{L_\alpha}$ — 升力线斜率（最重要的气动导数）
• $C_{L_q}$ — 俯仰速率对升力的影响
• $C_{L_{\delta_e}}$ — 升降舵效率

所有气动导数都可以通过风洞实验或 CFD 获取，或者从飞机的几何参数估算。

5.3 纵向/横航向解耦

Stevens 的核心观察：在平衡直线飞行状态下，纵向运动（速度、俯仰）和横航向运动（滚转、偏航）近似解耦。

因此 12×12 的系统矩阵可以分解为两个独立的子系统：

纵向（5阶）：$[u, w, q, \theta, \dot{P}_D]^T$（或常用 $[V, \alpha, q, \theta, h]^T$）

横航向（4阶）：$[v, p, r, \phi]^T$（或带偏航角 $[v, p, r, \phi, \psi]^T$）

这个解耦为模态分析奠定了基础。

---

六、完整概念地图

12状态6-DOF系统
     │
     ├── 平动 (3): u, v, w ← 牛顿第二定律
     ├── 转动 (3): p, q, r ← 欧拉方程
     ├── 姿态 (3): φ, θ, ψ ← 姿态运动学
     └── 位置 (3): PN, PE, PD ← 位置运动学

非线性系统: ẋ = f(x, u)
     │
     ↓ 小扰动线性化
     │
线性状态空间: Δẋ = AΔx + BΔu
     │
     ├── 纵向 (5阶): [Δu Δw Δq Δθ Δh]
     └── 横航向 (4阶): [Δv Δp Δr Δφ Δψ]

气动导数: 连接几何参数→力和力矩
     │
     ↓ 模态分析
     │
模态分析 → 特征值 → 稳定性 + 响应特性

---

七、核心公式速查卡

公式	含义
$\mathbf{V}_B = [u,v,w]^T$	空速机体系分量
$\boldsymbol{\omega}_B = [p,q,r]^T$	角速度
$\mathbf{x} = [u,v,w,p,q,r,\phi,\theta,\psi,P_N,P_E,P_D]^T$	12状态向量
$\Delta\dot{\mathbf{x}} = \mathbf{A}\Delta\mathbf{x} + \mathbf{B}\Delta\mathbf{u}$	线性化状态空间
$\bar{q} = \frac12\rho V_a^2$	动压
$C_L = C_{L_0} + C_{L_\alpha}\alpha + C_{L_q}\frac{q\bar{c}}{2V_a} + C_{L_{\delta_e}}\delta_e$	升力系数的气动导数展开

---

参考文献

1. Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft Control and Simulation (3rd ed., Chapters 1-2). Wiley. ISBN: 978-1-118-87098-4.
2. Etkin, B., & Reid, L. D. (1995). Dynamics of Flight: Stability and Control (3rd ed.). Wiley.
3. Cook, M. V. (2012). Flight Dynamics Principles (3rd ed.). Butterworth-Heinemann.
4. Bryan, G. H. (1911). Stability in Aviation: An Introduction to Dynamical Stability as Applied to the Motions of Aeroplanes. Macmillan.

---

下一节：飞行器模态分析：纵向与横航向

通过计算系统矩阵A的特征值，分析飞行器的固有力学模态。纵向：短周期、长周期（浮沉）；横航向：荷兰滚、滚转收敛、螺旋模态——这些都是你在PID调参文章中见过的专业术语，本文系统地讲清楚它们的数学来源。