2.3 系统辨识

0. 前言

​ 在之前的课程中，我们已经对四旋翼的各种物理参数和数学模型有了较为详细的了解。这些参数在后续控制算法设计中至关重要，其准确性直接决定了算法的稳定性。然而，对于手中的四旋翼，我们该如何较为精准地确定这些参数呢？有些参数如质量很容易精确测量，但像某些转动惯量、电机延迟、反扭系数等参数则不容易直接测得。这时，系统辨识通常是一种方便且常用的系统参数拟合方法。所谓系统辨识，就是通过对系统施加特定输入，并记录其输出或响应，从而在输入输出过程中体现系统的物理参数。我们可以采用如最小二乘等方法，从这些输入输出数据中反推并拟合出这些参数。在大多数情况下，这种方法能够达到较高的精度，完全可以直接用于后续控制算法的模型设计与仿真。

读完本文后，你将能够：

• 理解系统辨识在整个模型依赖控制体系中的位置：它是连接"理论方程"和"真实飞机"的桥梁
• 设计 3~4 次简单飞行实验，分别激励推力、转动惯量和扭矩系数对应的动力学
• 使用 sysid.tools 在线平台完成从 ULog 飞行日志到完整模型参数的自动辨识
• 理解每个辨识步骤背后的数学原理（最小二乘、EMA 延迟模型、联合优化）
• 获得一组可直接用于后续 LQR/MPC/NMPC 控制器设计的四旋翼物理参数

---

目录

1. 为什么要做系统辨识
2. 你需要准备什么
3. 飞行实验设计
4. 使用 sysid.tools 进行辨识
5. 基本原理

---

1. 为什么要做系统辨识

1.1 PID 的局限与模型依赖控制的需求

​ PX4 默认使用的 PID 级联控制器有一个很大的优点：你不需要知道飞行器的质量、转动惯量、推力系数等物理参数。调参就是反复试，找到一组让飞机不抖且能跟得上指令的 P/I/D 增益。这项工作虽然门槛不高，但是流程较为繁琐。然而PID因为不清楚模型的具体参数，只能在控制系统出现偏差之后来”进行补救“，并且其补救控制量对控制模型来说可能也不是最优的，导致控制系统的跟随性能可能较差，并且十分依赖参数调节效果。相比之下，基于模型的控制算法则可以根据目标控制量来计算准确的前馈输入，并能在出现外部扰动的时候给出合理的反馈控制量。对用户来说，调节参数的过程更加简单和形象。之前PID的参数调节决定了系统是否能够保持稳定，而模型控制算法的参数调节则是直接在保证稳定的基础上，进一步实现能量和性能的权衡调节，大幅度降低了参数调节时间和对经验的要求。

• LQR（线性二次型调节器）：需要状态空间模型，即 $\dot{x}=Ax+Bu$ 中的 $A$ 和 $B$ 矩阵
• MPC（模型预测控制）：需要在每个控制周期预测未来若干步的系统响应
• NMPC（非线性 MPC）：需要完整的非线性动力学模型

​ 以上所有方法都依赖一个前提：你有一个准确描述飞行器行为的数学模型。系统辨识的任务，就是通过飞行数据反推出这个模型的各个参数。

1.2 本教程能帮你得到什么

通过约 3~4 次简单飞行（每次约 40 秒），不借助推力台等额外设备或者测量方式，你可以辨识出以下参数：

参数	符号	作用
推力二次曲线参数	$a,b,c$	归一化油门（0~1）到实际推力（N）的二次映射
电机时间常数	$T_m$	电机响应的延迟特性
转动惯量	$I_{xx},I_{yy},I_{zz}$	绕三轴的转动惯量
扭矩系数	$K_{\tau }$	推力到反扭矩的比例系数

这些参数同一些能够轻松测量的数据例如质量和尺寸等，组成了四旋翼完整动力学模型中的所有物理参数。

---

2. 你需要准备什么

2.1 硬件

​ 你仅需准备一个电子秤用于测量飞行器总起飞质量（精度 0.1g 即可）。HumpRotor 的电机坐标、推力方向、扭矩方向等几何参数已由提供的模型文件给出，你无需手动测量。

2.2 软件

系统辨识使用在线工具完成，无需安装任何软件：

sysid.tools：https://sysid.tools/

这是一个基于浏览器的交互式系统辨识平台（由原论文作者团队开发并开源）。你只需上传 PX4 的 ULog 飞行日志，选择数据区间，即可自动完成辨识。其源代码托管在 GitHub：

https://github.com/arplaboratory/data-driven-system-identification

2.3 飞控日志配置

系统辨识需要飞控记录以下数据项（在 PX4 的 SDLOG_PROFILE 参数中启用高速率 + 系统辨识主题即可）：

ULog 字段	含义
actuator_motors_mux_control[0~3]	四个电机的归一化输出（0~1）
vehicle_angular_velocity_xyz[0~2]	三轴角速度
vehicle_angular_velocity_xyz_derivative[0~2]	三轴角加速度
vehicle_acceleration_xyz[0~2]	三轴加速度（机体坐标系）

---

3. 飞行实验设计

辨识需要 3~4 次飞行，分别激励不同的动力学特性。每次飞行约 40 秒。

3.1 飞行一：上下飞行（激励线加速度）

目的：辨识推力曲线参数 $a,b,c$ 和电机延迟时间常数 $T_m$。

操作：让飞行器在垂直方向上来回运动，油门从接近零推到最大，再拉回零。关键是 覆盖尽可能宽的油门范围，因为推力曲线需要在整个工作区间上拟合。

激励的物理量：z 轴加速度（机体坐标系）。

3.2 飞行二：Roll 大角度摇摆（激励横滚转动）

目的：辨识 $I_{xx}$（绕 x 轴的转动惯量）。

操作：在遥控器上快速、大幅度地推拉 Roll 摇杆，让飞行器绕 x 轴做明显的角加速运动。不需要飞得很远，重点是角度变化的剧烈程度（要有较多的角加速度过程）。

激励的物理量：x 轴角加速度。

3.3 飞行三：Pitch 大角度摇摆（激励俯仰转动）

目的：辨识 $I_{yy}$（绕 y 轴的转动惯量）。

操作：与 Roll 飞行类似，但改为快速、大幅度地推拉 Pitch 摇杆，让飞行器绕 y 轴做明显的角加速运动。

激励的物理量：y 轴角加速度。

说明：如果 Roll 和 Pitch 的激励在一次飞行中均足够充分，可以将飞行二与飞行三合并为一次飞行。分开飞行的好处在于每轴的数据更"干净"，辨识结果更稳定。

3.4 飞行四：Yaw 旋转（激励偏航转动）

目的：辨识扭矩系数 $K_{\tau }$。

操作：让飞行器原地来回旋转，快速改变偏航方向。Yaw 的动力学通常较慢，因此需要足够的激励强度才能获得有效数据。

激励的物理量：z 轴角加速度。

3.5 操作要点

• 所有飞行在手动/自稳模式下进行，不需要 GPS 或自主飞行
• 飞行中尽量让单个轴的运动幅度大一些，避免所有轴同时小幅晃动
• 确保日志中没有大的数据缺口
• 只有飞机完全处于空中的数据才是有效的，要严格排除飞机未起飞时或者碰撞地面时的数据

---

4. 使用 sysid.tools 进行辨识

本教程使用 FLU 坐标系（与 PX4 飞控内部一致）：x 轴指向前方（Front），y 轴指向左方（Left），z 轴指向上方（Up）。

4.1 模型文件准备

在开始辨识前，你需要准备一个描述飞行器基本几何的 JSON 模型文件。HumpRotor160 的参数如下：

电机位置（以下为 HumpRotor160 的数值，仅供参考，请以你自己飞行器的实际尺寸为准）：

电机	x (m)	y (m)	位置	扭矩方向
M1	0.0525	-0.0525	右前	+
M2	-0.0525	-0.0525	左前	−
M3	-0.0525	0.0525	左后	+
M4	0.0525	0.0525	右后	−

完整模型 JSON（以下为 HumpRotor160 GPS版本 的示例，请根据你的飞机实际情况来修改数值）：

{
  "gravity": 9.81,
  "mass": 0.2436,
  "rotor_positions": [
    [ 0.05657, -0.05657, 0],
    [-0.05657, -0.05657, 0],
    [-0.05657,  0.05657, 0],
    [ 0.05657,  0.05657, 0]
  ],
  "rotor_thrust_directions": [
    [0, 0, 1],
    [0, 0, 1],
    [0, 0, 1],
    [0, 0, 1]
  ],
  "rotor_torque_directions": [
    [0, 0,  1],
    [0, 0, -1],
    [0, 0,  1],
    [0, 0, -1]
  ]
}

​ 请将 mass、rotor_positions 替换为你自己的实测值。gravity、rotor_thrust_directions、rotor_torque_directions 通常无需修改。

4.2 上传日志与选取数据区间

1. 打开 https://sysid.tools/
2. 上传所有飞行的 ULog 文件（3~4 个）
3. 为每个辨识目标选取有效数据区间：
   • Thrust（推力辨识）：选取上下飞行的区间
   • Roll/Pitch Inertia（转动惯量）：选取大角度摇摆的区间
   • Yaw（偏航）：选取旋转飞行的区间

选取区间时注意：只保留有明显激励的段落，去掉起飞前和着陆后的无效数据。

4.3 辨识结果解读

电机推力曲线

​ 电机模型使用二次多项式，将归一化油门 $\omega \in [0,1]$ 映射为实际推力（N）。sysid.tools 输出的格式为 $a+bx+c{x}^2$，输出值按 [a, b, c] 顺序排列。在控制器代码中对应参数名 motor_a, motor_b, motor_c。

$$f(\omega )=a+b\cdot \omega +c\cdot {\omega }^2$$

​ 其中 $\omega$ 是飞控输出的归一化电机控制量（actuator_motors_mux_control，范围 0~1），$f$ 是该电机产生的推力（单位：牛顿）。三个系数的物理含义：

符号	含义
$a$（零次项）	截距修正，理想情况下应为 0（零油门 → 零推力），拟合值通常为很小的负数
$b$（一次项）	推力的线性分量，低油门区域主要由该项贡献
$c$（二次项）	推力的二次分量，高油门区域主导推力增长（空气动力学中推力近似与转速平方成正比）

推力 $\leftrightarrow$ 油门双向转换公式

​ 在实际控制器中，因为很多控制算法直接在推力上进行控制，所以需要频繁地进行推力和归一化油门之间的相互转换。

油门 → 推力（直接代入二次式）：

$$f(\omega )=a+b\cdot \omega +c\cdot {\omega }^2$$

推力 → 油门（反解二次方程，取正根）：

$$c\cdot {\omega }^2+b\cdot \omega +(a-f)=0$$$$\omega (f)=\frac{-b+\sqrt{b^2-4ac+4c\cdot f}}{2c}$$

电机时间常数

​ 一阶电机延迟模型：${\omega }_m(t+\mathrm{\Delta }t)=\alpha \cdot {\omega }_m(t)+(1-\alpha )\cdot {\omega }_{sp}(t)$，其中 $\alpha =e^{-\mathrm{\Delta }t/T_m}$。这个延迟如果不考虑，推力预测会有明显偏差。

转动惯量

​ 惯性比值（经验公式，仅供参考）：$C_{xy\to z}\approx 1.832$ 是基于多机型统计数据的回归值（见原论文图 1 及表 II）。该值对常见 X 型四旋翼具有较好的普适性，但对于异形机架或特殊质量分布的飞行器可能有较大偏差。由此推算的 $I_{zz}$ 和 $K_{\tau }$ 应视为近似估计。

扭矩系数

​ 扭矩系数的估计通常噪声较大（偏航动力学慢，信噪比低），但其精度对整体控制效果影响有限。

---

5. 基本原理

本节在保持可读性的前提下，给出每个辨识步骤的核心公式与推导逻辑。想深入了解完整推导的读者请参阅原论文 Data-Driven System Identification of Quadrotors Subject to Motor Delays 第 III 节。

5.1 总体辨识流程

整个辨识过程是顺序依赖的：前一步的结果是后一步的输入。流程图如下：

flowchart TD
    A[输入: ULog 飞行日志 + 基本模型 JSON] --> B[步骤一: 辨识 T_m 与推力曲线]
    B --> C[步骤二: 辨识 Ixx, Iyy]
    C --> D[步骤三: 估计 Izz]
    D --> E[步骤四: 辨识 Kτ]
    E --> F[输出: 完整动力学模型]

各步骤的输入/输出关系：

步骤	输入数据	待求参数	所依赖的前序结果
步骤一	上下飞行日志	$T_m$, $a,b,c$	无（但 $T_m$ 和 $a,b,c$ 需联合优化）
步骤二	Roll/Pitch 飞行日志	$I_{xx},I_{yy}$	$T_m$, $a,b,c$
步骤三	（经验公式计算得出）	$I_{zz}$	$I_{xx},I_{yy}$
步骤四	Yaw 飞行日志	$K_{\tau }$	$T_m$, $a,b,c$, $I_{zz}$

下面逐一展开每步的数学细节。

5.2 步骤一：推力曲线与电机时间常数

5.2.1 不考虑延迟时的线性方程组

四旋翼的线运动方程（机体坐标系下，FLU 约定）：

$$m\cdot {\mathbf{a}}_{\text{acc}}=\sum _{i=1}^4{\mathbf{r}}_{f_i}\cdot f_i$$

其中 ${\mathbf{r}}_{f_i}$ 是第 $i$ 个电机的推力方向（垂直机架配置下均为 $(0,0,1)$），$f_i$ 是该电机产生的推力，${\mathbf{a}}_{\text{acc}}$ 是加速度计读数。质量 $m$ 已由电子秤实测得到。

假设推力 $f_i$ 是电机转速 ${\omega }_{m_i}$ 的二次函数：

$$f_i=a_i+b_i{\omega }_{m_i}+c_i{\omega }_{m_i}^2$$

代入线运动方程并写成矩阵形式。对于第 $k$ 个采样时刻，z 轴方向（推力主方向）：

$$m\cdot a_{\text{acc},k}^{(z)}=\sum _{i=1}^4(a_i+b_i{\omega }_{m_i,k}+c_i{\omega }_{m_i,k}^2)$$

将所有 $N$ 个采样时刻堆叠，得到超定线性方程组：

$$\mathbf{A}\mathbf{x}=\mathbf{b}$$

其中各矩阵/向量的构造为：

$$\mathbf{b}={\left[\begin{array}{c}m\cdot a_{\text{acc},1}^{(z)}\\ m\cdot a_{\text{acc},2}^{(z)}\\ ⋮\\ m\cdot a_{\text{acc},N}^{(z)}\end{array}\right]}_{N\times 1},\mathbf{x}={\left[\begin{array}{c}{a}_1\\ b_1\\ c_1\\ ⋮\\ a_4\\ b_4\\ c_4\end{array}\right]}_{12\times 1}$$$$\mathbf{A}={\left[\begin{array}{ccccccc}1& {\omega }_{m_1,1}& {\omega }_{m_1,1}^2& \cdots & 1& {\omega }_{m_4,1}& {\omega }_{m_4,1}^2\\ 1& {\omega }_{m_1,2}& {\omega }_{m_1,2}^2& \cdots & 1& {\omega }_{m_4,2}& {\omega }_{m_4,2}^2\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 1& {\omega }_{m_1,N}& {\omega }_{m_1,N}^2& \cdots & 1& {\omega }_{m_4,N}& {\omega }_{m_4,N}^2\end{array}\right]}_{N\times 12}$$

通过最小二乘直接求解：

$${\mathbf{x}}^{\ast }=\arg \underset{\mathbf{x}}{min}\parallel \mathbf{A}\mathbf{x}-\mathbf{b}{\parallel }_2^2=({\mathbf{A}}^T\mathbf{A}{)}^{-1}{\mathbf{A}}^T\mathbf{b}$$

实际中四个电机的推力曲线通常取均值 $\overline{a}=\frac{1}{4}\sum _{i=1}^4{a}_i$（$b,c$ 同理），因为同型号电机和电调的推力特性几乎一致，取均值可以提高鲁棒性。

5.2.2 电机延迟的 EMA 模型

上述推导要求已知各时刻的电机转速 ${\omega }_{m_i}$，对于很多无法获取电机转速数据的飞控来说，其日志中只有电机设定值 ${\omega }_{s{p}_i}$（归一化到 0~1 的油门指令）。而电机又无法瞬时响应设定值变化，二者之间存在一阶动态关系：

$${\dot{\omega }}_{m_i}(t)=\frac{1}{T_m}({\omega }_{s{p}_i}(t)-{\omega }_{m_i}(t))$$

其中 $T_m$ 是一阶时间常数（单位：秒）。该微分方程在离散时间下的解为指数移动平均（EMA）递推式：

$${\omega }_{m_i}[k+1]=\alpha \cdot {\omega }_{m_i}[k]+(1-\alpha )\cdot {\omega }_{s{p}_i}[k]$$

其中 $\alpha =e^{-\mathrm{\Delta }t/T_m}$，$\mathrm{\Delta }t$ 为采样间隔。初始条件 ${\omega }_{m_i}[0]={\omega }_{s{p}_i}[0]$（假设初始时刻电机已达到设定值）。

5.2.3 $T_m$ 与推力曲线的联合优化

$T_m$ 本身也是未知的。辨识策略是在合理范围内（小型四旋翼通常 20ms~100ms，一般不会超过 100ms，具体范围应根据所用电机-电调组合的特性合理设定）做一维扫描 + 嵌套最小二乘：

对于每个候选 $T_m$ 值：

1. 用 EMA 递推式从 ${\omega }_{s{p}_i}$ 估计 ${\omega }_{m_i}$
2. 用 ${\omega }_{m_i}$ 构造 $\mathbf{A}$ 矩阵
3. 求解 ${\mathbf{x}}^{\ast }=({\mathbf{A}}^T\mathbf{A}{)}^{-1}{\mathbf{A}}^T\mathbf{b}$
4. 计算拟合均方根误差：$\text{RMSE}(T_m)=\sqrt{\frac{1}{N}\parallel \mathbf{A}{\mathbf{x}}^{\ast }-\mathbf{b}{\parallel }_2^2}$

最优 $T_m$ 为使 RMSE 最小的值：

$$T_m^{\ast }=\arg \underset{T_m\in [T_{min},T_{max}]}{min}\text{RMSE}(T_m)$$

这等价于在均匀先验下的最大后验估计（MAP）：

$$T_m^{\ast }=\arg \underset{T_m}{max}p(T_m\mid \mathcal{D})=\arg \underset{T_m}{min}\sum _{k=1}^N\parallel {\mathbf{b}}_k-{\mathbf{A}}_k\mathbf{x}{\parallel }_2^2$$

RMSE-$T_m$ 曲线通常在最优值处呈现清晰的最小值（参考原论文图 2）。找到 $T_m^{\ast }$ 后，对应的 ${\mathbf{x}}^{\ast }$ 即为最终推力曲线参数。

5.3 步骤二：转动惯量 $I_{xx},I_{yy}$ 的辨识

有了 $T_m$ 和推力曲线，就可以计算每个采样时刻每个电机的实际推力 $f_i(t)$。现在转向角运动方程。

5.3.1 x 轴转动惯量

绕 x 轴的角动力学（考虑纯垂直推力配置）：

$${\dot{\omega }}_{b_x}{I}_{xx}={\tau }_x$$

其中 ${\dot{\omega }}_{b_x}$ 是 x 轴角加速度（由 vehicle_angular_velocity_xyz_derivative[0] 直接提供，或对角速度做数值差分得到），${\tau }_x$ 是绕 x 轴的总力矩。

在小角速度假设下忽略陀螺进动项，${\tau }_x$ 仅由四个电机的推力对质心的力矩贡献：

$${\tau }_x=\sum _{i=1}^4({\mathbf{r}}_{p_i}\times {\mathbf{r}}_{f_i}{)}_x\cdot f_i$$

其中 ${\mathbf{r}}_{p_i}$ 是电机位置向量（由模型 JSON 提供），${\mathbf{r}}_{f_i}$ 是推力方向 $(0,0,1)$，叉积的 x 分量即为力矩臂。

将所有采样时刻堆叠为最小二乘形式：

$$\underset{{\mathbf{A}}_x}{\underbrace{\left[\begin{array}{c}{\dot{\omega }}_{b_x}[1]\\ {\dot{\omega }}_{b_x}[2]\\ ⋮\\ {\dot{\omega }}_{b_x}[N]\end{array}\right]}}\underset{{\mathbf{x}}_x}{\underbrace{\left[\begin{array}{c}{I}_{xx}\end{array}\right]}}=\underset{{\mathbf{b}}_x}{\underbrace{\left[\begin{array}{c}{\tau }_x[1]\\ {\tau }_x[2]\\ ⋮\\ {\tau }_x[N]\end{array}\right]}}$$

求解得到：

$$I_{xx}=\frac{\sum _{k=1}^N{\dot{\omega }}_{b_x}[k]\cdot {\tau }_x[k]}{\sum _{k=1}^N{\dot{\omega }}_{b_x}[k{]}^2}$$

5.3.2 y 轴转动惯量

完全对称地：

$$I_{yy}=\frac{\sum _{k=1}^N{\dot{\omega }}_{b_y}[k]\cdot {\tau }_y[k]}{\sum _{k=1}^N{\dot{\omega }}_{b_y}[k{]}^2}$$

其中 ${\tau }_y=\sum _{i=1}^4({\mathbf{r}}_{p_i}\times {\mathbf{r}}_{f_i}{)}_y\cdot f_i$。

5.4 步骤三：$I_{zz}$ 的估计

对于纯垂直推力（推力方向均为 $(0,0,1)$）的四旋翼，z 轴的力矩来源与 x、y 轴不同。x、y 轴力矩由推力对质心的几何叉积 ${\mathbf{r}}_{p_i}\times {\mathbf{r}}_{f_i}$ 产生，但由于推力方向平行于 z 轴，该叉积的 z 分量恒为零。z 轴力矩完全来自螺旋桨的反扭矩（rotor torque），由模型 JSON 中的 rotor_torque_directions 定义（${\mathbf{r}}_{{\tau }_i}$，z 分量为 $\pm 1$）：

$${\tau }_z=\sum _{i=1}^4{r}_{{\tau }_i,z}\cdot K_{\tau }\cdot f_i$$

因此 z 轴角动力学为：

$${\dot{\omega }}_{b_z}{I}_{zz}=K_{\tau }\sum _{i=1}^4{r}_{{\tau }_i,z}\cdot f_i$$

观察上式可以发现，只能辨识出比值 $I_{zz}/K_{\tau }$，而非各自的绝对值。为解决这个歧义性，利用了以下统计规律：

对大量四旋翼平台的实测数据（原论文表 II 及图 1）进行线性回归后发现，$(I_{xx}+I_{yy})/2$ 与 $I_{zz}$ 之间存在稳定的比例关系：

$$I_{zz}=\frac{I_{xx}+I_{yy}}{2}\cdot C_{xy\to z},C_{xy\to z}\approx 1.832$$

其物理直观：对于近似扁平对称的四旋翼，绕 z 轴的转动惯量约是绕水平轴平均值的 1.8 倍（质量分布在外围电机位置，离 z 轴较远）。

5.5 步骤四：扭矩系数 $K_{\tau }$ 的辨识

有了 $I_{zz}$，即可从偏航飞行数据中辨识 $K_{\tau }$。定义 z 轴力矩输入：

$${\tau }_z^{\text{(input)}}=\sum _{i=1}^4{r}_{{\tau }_i,z}\cdot f_i$$

z 轴角动力学：

$${\dot{\omega }}_{b_z}{I}_{zz}=K_{\tau }\cdot {\tau }_z^{\text{(input)}}$$

构造最小二乘问题：

$${\mathbf{A}}_z=\left[\begin{array}{c}{\tau }_z^{\text{(input)}}[1]\\ ⋮\\ {\tau }_z^{\text{(input)}}[N]\end{array}\right],{\mathbf{b}}_z=\left[\begin{array}{c}{\dot{\omega }}_{b_z}[1]\cdot I_{zz}\\ ⋮\\ {\dot{\omega }}_{b_z}[N]\cdot I_{zz}\end{array}\right]$$$$K_{\tau }=\frac{\sum _{k=1}^N{\tau }_z^{\text{(input)}}[k]\cdot ({\dot{\omega }}_{b_z}[k]\cdot I_{zz})}{\sum _{k=1}^N{({\tau }_z^{\text{(input)}}[k])}^2}$$

实践中 $K_{\tau }$ 的估计噪声通常较大（偏航动力学慢，信噪比低），但其精度对整体控制效果影响有限。若辨识结果偏差过大，可近似取同量级数值，不会显著影响最终控制效果。

本章参考资料与引申阅读已汇总至 参考资料。