四旋翼飞行动力学建模：AirSim 源码深度解析

一、为什么要读 AirSim 的动力学代码

Microsoft AirSim 是无人机仿真领域最知名的开源项目之一。它发布于 2017 年，核心动力学代码（AirLib 库）约 15 万行 C++，定位在 Unreal Engine/Unity 之上，为多旋翼提供高保真物理仿真。

虽然 AirSim 已于 2022 年进入维护模式（不再新增功能），但其多旋翼飞行动力学实现——特别是转子推力模型、气动阻力和传感器噪声建模——至今仍是 C++ 开源仿真中相当完整的参考，值得逐行阅读。

本文聚焦 AirSim 动力学内核（AirLib），不涉及 Unreal Engine 渲染层、Python API 和 ROS 桥接，只看 AirLib/include/ 下的物理建模代码。

分析范围：

模块	所在文件	代码行数
物理引擎主循环	FastPhysicsEngine.hpp	~465 行
多旋翼动力学体	MultiRotorPhysicsBody.hpp	~230 行
螺旋桨参数	RotorParams.hpp	~65 行
螺旋桨执行器	RotorActuator.hpp	~147 行
多旋翼参数配置	MultiRotorParams.hpp	~550 行
环境大气模型	Environment.hpp	~130 行
运动学状态	Kinematics.hpp	~110 行
物理体基类	PhysicsBody.hpp	~270 行
IMU 传感器	ImuSimple.hpp	~118 行

全部文件加在一起约 2000 行，是 AirSim 动力学内核的精华所在。

---

二、整体架构：四层设计

AirLib 的动力学代码呈四层架构：

┌──────────────────────────────────────────────────────┐
│  VehicleApiBase / MultirotorApiBase                  │
│  └─ MavLinkMultirotorApi / SimpleFlightApi           │  ← API 层
├──────────────────────────────────────────────────────┤
│  MultiRotorPhysicsBody                               │
│  └─ 继承自 PhysicsBody                               │  ← 动力学体
├──────────────────────────────────────────────────────┤
│  FastPhysicsEngine                                   │
│  └─ 继承自 PhysicsEngineBase                         │  ← 物理引擎
├──────────────────────────────────────────────────────┤
│  RotorActuator × N                                   │  ← 力顶点(转子)
│  PhysicsBodyVertex × 6 (DragFace)                    │  ← 力顶点(阻力面)
│  SensorCollection (IMU/Baro/GPS/Mag)                 │  ← 传感器
└──────────────────────────────────────────────────────┘

每一层只做一件事：

• 力顶点层（RotorActuator + 6 个 DragVertex）：接受控制信号，产生力/力矩
• 物理引擎层（FastPhysicsEngine）：遍历所有物理体，做时间积分和碰撞响应
• 动力学体层（MultiRotorPhysicsBody）：组装转子参数、阻力顶点、传感器
• API 层：飞控固件的控制逻辑（SimpleFlight 或 PX4/ArduCopter）

这个设计的巧妙之处在于：动力学体和物理引擎完全解耦。物理引擎只看到 PhysicsBody 的抽象接口（wrenchVertexCount()、dragVertexCount()），不关心具体是四旋翼还是汽车。如果要添加新机型，只需实现 PhysicsBody 子类。

关键还在于，所有代码都是 header-only。没有 .cpp 文件，模板和内联函数全在头文件里——好处是编译期优化充分，坏处是编译慢（每个包含头文件的编译单元都要重新实例化）。

---

三、物理引擎主循环（FastPhysicsEngine）

FastPhysicsEngine::update() 是每一帧的入口：

virtual void update() override {
    PhysicsEngineBase::update();
    for (PhysicsBody* body_ptr : *this) {
        updatePhysics(*body_ptr);
    }
}

它遍历所有注册的 PhysicsBody，为每个调用 updatePhysics。核心逻辑在 getNextKinematicsNoCollision() 中：

3.1 时间步长

TTimeDelta dt = clock()->updateSince(body.last_kinematics_time);

AirSim 采用变步长模式——dt 由实际系统时钟决定（通常是 Unreal Engine 的渲染帧间隔），而不是固定值。这种模式的好处是易于与游戏引擎集成，坏处是帧率抖动时物理也不稳定。

3.2 力/力矩合成

// 1. 从力顶点收集 wrench（转子推力/力矩）
const Wrench body_wrench = getBodyWrench(body, current.pose.orientation);

// 2. 计算平均速度用于阻力
avg_linear = current.twist.linear + current.accelerations.linear * (0.5f * dt_real);
avg_angular = current.twist.angular + current.accelerations.angular * (0.5f * dt_real);

// 3. 计算阻力
const Wrench drag_wrench = getDragWrench(body, current.pose.orientation, avg_linear, avg_angular, wind);

// 4. 合成总 wrench
next_wrench = body_wrench + drag_wrench;

getBodyWrench() 遍历所有 wrenchVertex（即 rotor），用 τ = r × F 计算每个转子对质心产生的力矩：

for (uint i = 0; i < body.wrenchVertexCount(); ++i) {
    const auto& vertex_wrench = vertex.getWrench();
    wrench += vertex_wrench;
    wrench.torque += vertex.getPosition().cross(vertex_wrench.force);  // τ = r × F
}

注意：力矩在体坐标系合成，但力会被转换到世界坐标系。

3.3 刚体运动方程（Verlet 积分）

AirSim 使用半隐式 Verlet 积分，不是常见的 RK4：

加速度计算：

// 线加速度：F = ma 加上重力
next.accelerations.linear = (next_wrench.force / body.getMass()) + gravity;

// 角加速度：Euler 旋转方程
// L = I·ω, dL/dt = τ - ω × L
const Vector3r angular_momentum = body.getInertia() * avg_angular;
const Vector3r angular_momentum_rate = next_wrench.torque - avg_angular.cross(angular_momentum);
next.accelerations.angular = body.getInertiaInv() * angular_momentum_rate;

其中 angular_momentum_rate = τ - ω × (I·ω) 正是 Euler 旋转方程（体坐标系形式）：

$$I \dot{\boldsymbol{\omega}} + \boldsymbol{\omega} \times (I \boldsymbol{\omega}) = \boldsymbol{\tau}$$

速度更新：

// Verlet 风格——用当前和下一个加速度的平均值
next.twist.linear = current.twist.linear + (current.accelerations.linear + next.accelerations.linear) * (0.5f * dt_real);
next.twist.angular = current.twist.angular + (current.accelerations.angular + next.accelerations.angular) * (0.5f * dt_real);

位姿更新：

// 位置
next.pose.position = current_pose.position + avg_linear * dt_real;

// 姿态四元数
real_T angle_per_unit = avg_angular.norm();
if (angle_per_unit > 0) {
    AngleAxisr angle_dt_aa(angle_per_unit * dt_real, avg_angular / angle_per_unit);
    Quaternionr angle_dt_q(angle_dt_aa);
    next.pose.orientation = current_pose.orientation * angle_dt_q;
    next.pose.orientation.normalize();
}

对四元数更新，核心逻辑是将角速度向量转化为四元数增量，右乘到当前姿态上。证明如下：

• 当前姿态为 $q_0$，世界坐标系向量 $v$ 体轴表示为 $q_0' v q_0$
• 体轴角速度增量 $q_1$ 作用后为 $q_1(q_0' v q_0) q_1'$
• 再转回世界系：$q_0(q_1(q_0' v q_0) q_1')q_0' = (q_0 q_1) v (q_0 q_1)'$
• 所以新姿态为 $q_0 q_1$

这是经典的四元数姿态递推方法。

速度剪裁： AirSim 对速度做了上限保护——当线速度超过光速时会被截断。这个保护主要防止控制器 bug 导致数值爆炸。

3.4 碰撞响应

AirSim 的碰撞响应参考了 Chris Hecker 的经典论文（GDMPhys3.pdf）和 Gaffer On Games 的文章。使用冲量法：

// 冲量公式: j = -(1 + R) V·N / (1/m + [I⁻¹(r×N)×r]·N)
const real_T impulse_mag = -contact_vel_body.dot(normal_body) * (1 + restitution) / impulse_mag_denom;

// 更新速度
next.twist.linear  = vcur_avg + collision_info.normal * (impulse_mag / body.getMass());
next.twist.angular = angular_avg + r.cross(normal_body) * impulse_mag;

同时还实现了库仑摩擦力模型，在切向方向施加摩擦冲量。地面碰撞时会把弹性系数降为 0（防止弹跳）、摩擦系数升为 1，并锁定姿态（pitch/roll = 0）。

---

四、螺旋桨模型（RotorActuator + RotorParams）

4.1 推力/扭矩公式

AirSim 的螺旋桨模型基于 UIUC 风洞数据库（propDB），使用经典的动量理论公式：

// RotorParams.hpp
// F = C_T · ρ · n² · D⁴
// τ = C_P · ρ · n² · D⁵ / (2π)
max_thrust = C_T * air_density * revolutions_per_second² * propeller_diameter⁴;
max_torque = C_P * air_density * revolutions_per_second² * propeller_diameter⁵ / (2 * M_PI);

默认参数对应 GWS 9×5 三叶桨在 6396 RPM 下的风洞数据：

参数	数值	来源
$C_T$	0.109919	UIUC 风洞实测
$C_P$	0.040164	UIUC 风洞实测
转速	6396.667 RPM	UIUC 测试工况
桨直径	0.2286 m（9 inch）	DJI Phantom 2 规格
最大推力	4.179 N	计算值
最大扭矩	0.0556 N·m	计算值

空气密度修正通过 RotorActuator::updateEnvironmentalFactors() 完成：

air_density_ratio_ = environment_->getState().air_density / air_density_sea_level_;
wrench.force = normal * output_.thrust * air_density_ratio_;
wrench.torque = normal * output_.torque_scaler * air_density_ratio_;

高海拔推力损失会自动模拟。

4.2 控制信号处理链

控制信号从飞控到推力的完整路径：

控制信号 (0~1) → 一阶低通滤波 → speed = √(u · max_speed²) → thrust = u · max_thrust

其中 FirstOrderFilter 的传递函数为：

control_signal_filter_.initialize(params_.control_signal_filter_tc, 0, 0);

时间常数 control_signal_filter_tc = 0.005s（5ms），对应一阶低通滤波的截止频率约 32Hz。这个滤波器模拟了 ESC（电子调速器）和电机的动态响应。

关键问题：推力被建模为控制信号的线性函数（thrust = u · max_thrust），转速是开平方关系（speed = √(u · max_speed²)）。这意味着假设了 PWM 与 RPM 为线性关系——实际中，无刷电机的 PWM-RPM 曲线存在非线性和死区。

4.3 转子布局

MultiRotorParams.hpp 内置了多种机架配置模板：

• QuadX（4 旋翼 X 型）：标准 F450 框架，轴距 227.5mm
• HexX（6 旋翼 X 型）：轴距 227.5mm
• OctoX（8 旋翼 X 型）：轴距 227.5mm
• Flamewheel（F450 兼容）、Blacksheep（TBS Discovery）

QuadX 的布局采用经典 ArduPilot 编号方式：

        x-axis
   (2)  |   (0)
        |
---------------- y-axis
        |
   (1)  |   (3)

RotorPose 包括位置、法线方向（NED 坐标系，始终 (0, 0, -1)）和旋转方向（CW/CCW）。

4.4 转子推力的局限性

1. C_T 为常数，不随前进比 J 变化：实际中当飞行器快速前飞时，来流速度增大会减小螺旋桨的有效攻角，使 C_T 显著下降。AirSim 的常数模型会高估高速飞行时的推力。

2. 忽略了桨叶挥舞动力学：螺旋桨高速旋转时的桨叶挥舞（flapping）会产生陀螺效应和周期变距效应，这些在高速大机动时不可忽略。

3. 无滑流模型：螺旋桨滑流会冲刷机翼/机身，改变当地气流，这在 AirSim 中完全没有——因为四旋翼模型也没有机翼。

---

五、阻力模型（DragFaces）

AirSim 的阻力模型是一个六面体盒模型，非对称四旋翼的机体被简化为长方体：

// MultiRotorPhysicsBody::createDragVertices()
void createDragVertices() {
    // 机身六面投影面积
    real_T top_bottom_area = params.body_box.x() * params.body_box.y();
    real_T left_right_area = params.body_box.x() * params.body_box.z();
    real_T front_back_area = params.body_box.y() * params.body_box.z();

    // 加上桨盘对侧向的贡献
    real_T propeller_xsection = M_PI * diameter * height;

    // 阻力系数因子
    Vector3r drag_factor_unit = Vector3r(
        front_back_area + N_rotors * propeller_xsection,   // X 方向
        left_right_area + N_rotors * propeller_xsection,   // Y 方向
        top_bottom_area + N_rotors * propeller_area         // Z 方向
    ) * params.linear_drag_coefficient / 2;

    // 六个 DragVertex，每个在对应面的中心位置
    drag_faces_.emplace_back(Vector3r(0, 0, -z/2), Vector3r(0, 0, -1), factor_z);
    drag_faces_.emplace_back(Vector3r(0, 0,  z/2), Vector3r(0, 0,  1), factor_z);
    // ... X 和 Y 方向类似
}

阻力力方向始终与面法线反向，大小与相对速度平方成正比：

// FastPhysicsEngine::getDragWrench()
if (vel_comp > kDragMinVelocity) {
    Vector3r drag_force = vertex.getNormal() * (-vertex.getDragFactor() * air_density * vel_comp²);
    Vector3r drag_torque = vertex.getPosition().cross(drag_force);
    wrench.force += drag_force;    // 体坐标系
    wrench.torque += drag_torque;
}

最后转换到世界坐标系。注意：力矩通过 r × F 计算，r 是 DragVertex 的位置矢量。也就是说，即使阻力力在机头方向，如果作用点不在质心，也会产生俯仰力矩——这是合理的。

默认的阻力系数 linear_drag_coefficient = 1.3/4 = 0.325。文档注释提到参考了 klsin.bpmsg.com 的数据，除以 4 是为了”更激进的流线型设计，追求更高的极速”。Flamewheel 机架则会乘以 4 来让高速更真实。

核心缺陷：这个阻力模型完全没有升力面。四旋翼在高速前飞时，机身会因入射角产生升力（或负升力），影响俯仰平衡和总拉力需求。AirSim 的阻力面是纯耗散的（只产生阻力不产生升力），因此无法模拟：

• 前飞时的机身升力
• 侧飞时的侧向气动
• 角速率产生的气动阻尼力矩（除 r × F 带来的间接力矩外）

对低速机动（< 10 m/s）影响不大，但截击机的 80-100 m/s 飞行速度下，机身气动升力约占总重力的 5-15%，忽略它将导致显著误差。

如果要改进，需要在 createDragVertices() 中给每个面同时设置升力和阻力系数（类似 JSBSim 的 FGAerodynamics 表格），而非只有 drag_factor。

---

六、环境与大气模型（Environment）

Environment 实现了 MSL 标准大气模型：

void updateState(State& state) {
    // NED 位置 → 经纬高
    geodetic_converter_.ned2Geodetic(state.position, state.geo_point);

    // US 标准大气
    real_T geo_pot = EarthUtils::getGeopotential(altitude_km);
    state.temperature = EarthUtils::getStandardTemperature(geo_pot);
    state.air_pressure = EarthUtils::getStandardPressure(geo_pot, temperature);
    state.air_density = EarthUtils::getAirDensity(air_pressure, temperature);

    // 重力随高度变化
    state.gravity = Vector3r(0, 0, EarthUtils::getGravity(altitude));
}

重力：考虑了高度对引力加速度的影响（虽然变化量很小）：

// EarthUtils
real_T getGravity(double altitude_m) {
    return G_0 * pow(EARTH_RADIUS / (EARTH_RADIUS + altitude_m), 2);
}

标准重力 $g_0 = 9.80665$ m/s²，海拔 10 km 时 $g ≈ 9.774$，仅差 0.3%。

温度：使用 US Standard Atmosphere 的分层模型（对流层温度直减率 -6.5 K/km）。

空气密度：影响转子推力和阻力，根据温度和气压计算。

位置：使用 GeodeticConverter（经纬高 ↔ NED 局部坐标），支持通过 setHomeGeoPoint() 设置原点。

这些环境参数在物理引擎中影响两处：

1. 转子推力 → air_density_ratio_ 修正
2. 空气阻力 → air_density 参与阻力计算

值得注意的是，AirSim 的环境模型是静态的（不随时间变化），没有风切变、湍流、雷暴等气象效应，只有一个可设置的 wind_ 常量向量。

---

七、传感器模型

AirSim 的传感器模型通过 SensorBase 统一接口实现。每个传感器持有指向 Kinematics::State 和 Environment 的指针（Ground Truth），然后在真值基础上添加噪声。

7.1 IMU（ImuSimple）

IMU 噪声模型参考了 Oliver J. Woodman 的惯性导航入门论文：

void addNoise(Vector3r& linear_acceleration, Vector3r& angular_velocity) {
    real_T sqrt_dt = sqrt(max(dt, min_sample_time));

    // 陀螺仪：角度随机游走 (ARW) + 偏置漂移
    real_T gyro_sigma_arw = params_.gyro.arw / sqrt_dt;
    angular_velocity += gauss_dist.next() * gyro_sigma_arw + state_.gyroscope_bias;
    real_T gyro_sigma_bias = gyro_bias_stability_norm * sqrt_dt;
    state_.gyroscope_bias += gauss_dist.next() * gyro_sigma_bias;

    // 加速度计：速度随机游走 (VRW) + 偏置漂移
    real_T accel_sigma_vrw = params_.accel.vrw / sqrt_dt;
    linear_acceleration += gauss_dist.next() * accel_sigma_vrw + state_.accelerometer_bias;
    real_T accel_sigma_bias = accel_bias_stability_norm * sqrt_dt;
    state_.accelerometer_bias += gauss_dist.next() * accel_sigma_bias;
}

这是三项噪声的叠加：

噪声项	物理含义	数学形式	参数配置
角度随机游走 (ARW)	高速率白噪声	$\sigma/\sqrt{dt} \cdot \mathcal{N}(0,1)$	gyro.arw
偏置稳定性	偏置的低频漂移	$(\sigma/\sqrt{\tau})\sqrt{dt} \cdot \mathcal{N}(0,1)$	gyro.bias_stability, gyro.tau
开机偏置	每次上电的初始偏置	常数 $+\mathcal{N}(0, \sigma_{turn\_on})$	gyro.turn_on_bias

陀螺仪真实输出 = 真值角速度 + ARW 噪声 + 当前偏置（随机的偏置持续漂移）。

加速度计则输出体轴加速度（扣除重力）——这是真正的 IMU 输出，因为 IMU 无法区分加速度和重力：

output.linear_acceleration = ground_truth.kinematics->accelerations.linear 
                             - ground_truth.environment->getState().gravity;
// 再变换到体坐标系
output.linear_acceleration = VectorMath::transformToBodyFrame(output.linear_acceleration, orientation, true);

7.2 其他传感器

AirSim 还实现了：

• BarometerSimple：基于标准大气模型的气压高度
• GpsSimple：位置 + 速度 + 协方差，带接收机噪声
• MagnetometerSimple：地磁场模型，带硬铁/软铁校正参数

每个传感器的接口相同——initialize() 时传入 Ground Truth 指针，update() 时取真值加噪声。

传感器模型的共同缺陷：所有传感器都是零延迟的。真值计算发生在当前帧，噪声立即叠加，没有模拟：

• 传感器自身采样延迟（IMU 通常有 1-5 ms 的组延迟）
• 传输延迟（SPI/I2C 总线延迟）
• 不同传感器的时间戳对齐（通常用硬件同步）

对频率 < 50 Hz 的控制系统影响不大，但高频姿态控制（> 200 Hz）就需要考虑这些延迟。

---

八、飞控固件集成方式

AirSim 支持三种飞控模式：

模式	API 类	说明
SimpleFlight	SimpleFlightApi	内置简单飞控，纯高度+速度PID
PX4 MAVLink	MavLinkMultirotorApi	通过UDP串口外接 PX4 HITL/SITL
ArduCopter	MavLinkMultirotorApi	同上，连 ArduPilot SITL

最关键的是 Px4MultiRotorParams，它根据 connection_info_.model 选择机架参数，并创建 MavLinkMultirotorApi：

virtual std::unique_ptr<MultirotorApiBase> createMultirotorApi() override {
    auto api = make_unique<MavLinkMultirotorApi>();
    api->initialize(connection_info_, &getSensors(), true);
    return api;
}

MavLinkMultirotorApi 通过 MAVLink 接收 PX4 的 HIL_ACTUATOR_CONTROLS 消息，提取 0~1 的混控器输出，传递给 RotorActuator.setControlSignal()。反向路径，传感器数据从 getSensors() 取到后，封装成 HIL_SENSOR、HIL_GPS 等消息发回 PX4。

这建立了 AirSim 动力学 → PX4 飞控（实际固件）→ 气动响应 → AirSim 更新传感器 → 再喂给 PX4 的闭环。AirSim 只提供”真实”物理，飞控算法完全由 PX4 执行。

---

九、架构设计中的取舍

AirSim 的动力学设计展现了一系列明确的精度 vs 性能 vs 可维护性权衡：

选择正确的抽象层次

RotorActuator 继承 PhysicsBodyVertex 是核心抽象——每个转子是一个力/力矩产生点，对物理引擎透明。类似地，六个 drag face 也是 PhysicsBodyVertex 的子类。物理引擎不区分转子力和阻力，统一通过 getWrench() 收集。

这个模式定义了通用的力顶点接口（位置、法线、wrench），使得添加新的力产生器（如帆、舵面、推进器）不需要改动物理引擎。

简化的代价

简化	对低速四旋翼的影响	对高速飞行的影响
常量 C_T	< 5% 误差	10-30% 推力高估
长方体阻力模型	< 5% 差异	升力面缺失，姿态误差
无桨叶挥舞	悬停稳定，外环 OK	高速转弯时姿态震荡被低估
零延迟传感器	> 50 Hz 控制略微乐观	不明显
Verlet 积分	100 Hz 更新率下误差可忽略	大角速度时误差累积

为什么这个设计对仿真足够

AirSim 的主要设计目标是验证外环自主算法（路径规划、避障、航点跟踪），而不是飞控级动力学。在这些场景中：

• 飞行速度 < 15 m/s
• 加速度 < 2g
• 姿态角 < 30°

在这些条件下，上述简化带来的误差小于传感器噪声和风扰，完全不影响算法验证结论。

但如果目标是截击机级仿真（速度 > 80 m/s，机动到 90° 俯仰，法向过载 > 3g），每个简化都需要逐一升级。

---

十、总结：能学到什么

值得借鉴的设计

1. 力顶点模式：PhysicsBodyVertex 将力和力矩的生成统一抽象，新增执行器/舵面无需改动物理引擎
2. 双顶点体系：wrenchVertex（主动力）和 dragVertex（耗散力）分离设计
3. Header-only 架构：对中大型模板库来说是一种切实可行的发布方式（Eigen 也这么做）
4. 环境对象分离：大气模型独立于动力学体，可在仿真中动态注入风向/温度/气压
5. 碰撞响应实装：完整的冲量模型 + Coulomb 摩擦 + ground lock，有生产级参考价值

需要留意的局限

1. C_T 常量模型：不适用于高速前飞场景
2. 无升力面：机身、机翼、平尾的升力完全缺失
3. 积分精度：Verlet 在 100 Hz 以下时精度低于 RK4
4. 零延迟传感器：忽略采样/传输延迟
5. 静态大气：无湍流、风切变等扰动场

AirSim 的动力学代码编写质量相当高——注释详细（说明公式来源和论文引用）、边界条件处理完整（NaN 检测、速度截断、四元数归一化）。对需要在 C++ 中实现多旋翼仿真的开发者来说，AirLib/include/ 下的约 2000 行代码是优秀的参考起点。