无人机物理学新突破

背景

本文《Drone Physics》源自一篇发表于 Hacker News 的技术文章，原作者基于其在自适应控制领域的研究成果以及于 AIAA DASC 2023 会议上发表的关于多旋翼无人机 Python 仿真框架 multirotor 的工作，系统性地介绍了描述无人机飞行物理所需的基础数学与力学模型。文章面向具备初等线性代数、初等微积分和初等经典力学知识的读者，旨在为无人机建模、仿真与控制提供清晰的坐标系、状态变量、力与力矩、运动方程以及螺旋桨控制逻辑的理论框架。文中符号约定在很大程度上借鉴了 Charles Tytler 的出色工作。

核心内容

飞行器描述

多旋翼 UAV 具有六个自由度：三个平动轴（$x, y, z$）和三个转动轴（$\phi, \theta, \psi$）。为了使用坐标，需要约定一个坐标系。通常采用北-东-地（North-East-Down, NED）系统：正 $x$ 轴方向为“前向”/北，正 $y$ 轴为“右向”/东，正 $z$ 轴为“向下”。这是一个右手坐标系：绕轴的正旋转方向为拇指指向轴正向时其余手指弯曲的方向。例如，正 $z$ 轴旋转从 $+x$ 转向 $+y$。

使用两个参考坐标系表示无人机的状态：

• 惯性参考系 $n$：静态参考系，轴与世界坐标系对齐，表示为列向量 $\hat{n} = [\hat{x}^n, \hat{y}^n, \hat{z}^n]^T$。
• 机体固连非惯性参考系 $b$：轴与运动中的刚体（无人机）重心对齐，随机体移动和旋转，表示为列向量 $\hat{b} = [\hat{x}^b, \hat{y}^b, \hat{z}^b]^T$。机体坐标系的原点固定在无人机质心上。

平动表示

机体坐标系 $b$ 固定在质心处，因此无人机相对于机体坐标系的位置位移始终为 0。无人机的位置由机体坐标系相对于惯性坐标系的位置矢量 $\hat{r}^n$ 表示。无人机的速度是机体坐标系相对于惯性坐标系的速度，这里选择在机体坐标系中表示该速度。

转动表示

无人机在惯性参考系中的姿态遵循 Tait-Bryan 角约定，即姿态可由三个顺序旋转描述：偏航（$\psi$）、俯仰（$\theta$）、滚转（$\phi$），且顺序不可颠倒。从惯性系 $\hat{n}$ 开始，先绕惯性 $z$ 轴旋转 $\psi$（偏航），得到 $\hat{n}\psi$；再绕新 $y$ 轴旋转 $\theta$（俯仰），得到 $\hat{n}{\psi,\theta}$；最后绕新 $x$ 轴旋转 $\phi$（滚转），得到 $\hat{n}_{\psi,\theta,\phi}$，即机体参考系。

无人机的角速度 $\hat{w}^T = [\omega_x, \omega_y, \omega_z]$ 是机体坐标系绕自身轴的瞬时旋转速率。各姿态角定义在各自旋转阶段的参考系中（$\hat{n}, \hat{n}\psi, \hat{n}{\psi,\theta} = \hat{b}$）。滚转旋转最后进行，结果即为机体坐标系，因此 $\omega_x$ 等滚转速率 $\dot{\phi}$。俯仰在初始偏航之后、滚转之前进行，因此 $\omega_y$ 是经过滚转 $\phi$ 旋转到机体坐标系下的俯仰速率。偏航首先在惯性系中进行，因此 $\omega_z$ 是经过俯仰 $\theta$ 和滚转 $\phi$ 旋转到机体坐标系下的偏航速率。可用矩阵方程表示：

$$ \begin{bmatrix} \omega_x \ \omega_y \ \omega_z \end{bmatrix} = R(\phi)\cdot R(\theta) \begin{bmatrix} 0 \ 0 \ \dot{\psi} \end{bmatrix} + R(\phi) \begin{bmatrix} 0 \ \dot{\theta} \ 0 \end{bmatrix} + \begin{bmatrix} \dot{\phi} \ 0 \ 0 \end{bmatrix} $$

化简得：

$$ \begin{bmatrix} \omega_x \ \omega_y \ \omega_z \end{bmatrix} = \begin{bmatrix} \dot{\phi} + \dot{\psi} \sin{\theta}\ - \dot{\psi} \sin{\phi} \cos{\theta} + \dot{\theta}\ \dot{\psi} \cos{\phi} \cos{\theta} \end{bmatrix} $$

惯性系与机体系的转换

机体系与惯性系可能因（1）位移和（2）旋转而不同。机体系原点固定于 UAV 原点，因此机体系相对于惯性系的位移即为 UAV 的惯性位置：$\hat{r}^n = [x, y, z]^T$。

定义在机体系原点处的矢量，若平移到惯性系原点，则会表现为旋转。给定惯性系中的矢量 $\hat{\mathcal{V}}^n = [x^n, y^n, z^n]^T$ 和同一矢量在机体系原点处的表示 $\hat{\mathcal{V}}^b = [x^b, y^b, z^b]^T$，从机体系到惯性系的旋转矩阵 $R_b^n$ 定义如下（每个基本旋转矩阵为）：

$$ R(\phi) = \begin{bmatrix} 1 & 0 & 0 \ 0 & c\phi & -s\phi \ 0 & s\phi & c\phi \end{bmatrix}, \quad R(\theta) = \begin{bmatrix} c\theta & 0 & s\theta \ 0 & 1 & 0 \ -s\theta & 0 & c\theta \end{bmatrix}, \quad R(\psi) = \begin{bmatrix} c\psi & -s\psi & 0 \ s\psi & c\psi & 0 \ 0 & 0 & 1 \end{bmatrix} $$

其中 $c \mid s$ 表示余弦和正弦。从惯性系到机体系的变换为：

$$ \hat{\mathcal{V}}^b = R(\phi)\cdot R(\theta) \cdot R(\psi) \cdot \hat{\mathcal{V}}^n = R_n^b \hat{\mathcal{V}}^n $$

即：

$$ \begin{bmatrix} x^b \ y^b \ z^b \end{bmatrix} = \begin{bmatrix} c\psi c\theta & - s\psi c\theta & s\theta\ s\phi s\theta c\psi + s\psi c\phi & - s\phi s\psi s\theta + c\phi c\psi & - s\phi c\theta\ s\phi s\psi - s\theta c\phi c\psi & s\phi c\psi + s\psi s\theta c\phi & c\phi c\theta \end{bmatrix} \begin{bmatrix} x^n \ y^n \ z^n \end{bmatrix} $$

此变换可用于从惯性视角观察作用在机体系上的力。

旋转参考系中的导数协调

如果矢量本身在参考系旋转的同时发生变化，则其时间变化率并不简单。变化率由两部分组成：（1）矢量在机体系中的变化；（2）坐标系坐标的变化。

即在以瞬时角速度 $\hat{\omega} = [\omega_x, \omega_y, \omega_z]$ 绕自身轴旋转的机体系中，在共位的惯性系中测量到的机体系中矢量 $\hat{\mathcal{V}}^b = \hat{b} \cdot \hat{\mathcal{V}}$ 的时间导数由 Transport 定理给出：

$$ \frac{d \hat{\mathcal{V}}^b}{d t} = \hat{b} \cdot \frac{d \hat{\mathcal{V}}}{d t} + \frac{d \hat{b}}{d t} \cdot \hat{\mathcal{V}} = \hat{b} \cdot \frac{d \hat{\mathcal{V}}}{d t} + \hat{\omega} \times \hat{\mathcal{V}}^b = \hat{b} \cdot \frac{d \hat{\mathcal{V}}}{d t} + \begin{bmatrix} 0 & -\omega_z & \omega_y \ \omega_z & 0 & -\omega_x \ -\omega_y & \omega_x & 0 \end{bmatrix} \hat{\mathcal{V}} $$

Transport 定理与叉积的关系在此处为后续推导运动方程奠定了基础。

关键要点