1. 模型推导的整体思路为什么四旋翼的数学模型让人又爱又恨四旋翼无人机这玩意儿飞行控制圈里人人都说“先建模再调参”可真到了自己动手推公式的时候不少人直接被那堆矩阵、角度、力矩给劝退了。我一开始做飞控时也踩过这个坑照着论文抄了一版动力学方程结果扔到真机上姿态稳不住稍微打杆就乱飘。后来才明白不是方程抄错了是我压根没搞懂每个项对应飞机上的哪个物理部分、哪些项能忽略、哪些项必须精确标定。聊四旋翼数学模型之前先说一个最核心的认知四旋翼本质上是一个欠驱动、强耦合、非线性的系统。它有六个自由度三个位置方向、三个姿态角却只有四个电机输入这意味着你想让飞机往前飞并不能直接给一个“前进力”只能通过改变姿态、让升力倾斜用水平分量去驱动位移。这个“通过姿态间接控制位置”的逻辑就是整篇数学模型里贯穿始终的主线。所以推导模型时我们通常分两层来拆位置动力学描述飞机在空间里的平移运动和姿态动力学描述飞机绕自身三个轴的旋转运动。两层之间通过姿态矩阵耦合在一起。写出来就是标准的牛顿-欧拉方程体系。这套模型能解决什么问题往小了说它能告诉你“给定四个螺旋桨转速飞机会产生多大加速度”这是仿真模拟、飞控算法开发的基础往大了说后面的状态估计、姿态解算、PID控制器设计、甚至故障检测全都长在这棵树上。适合谁来参考自己做四旋翼毕设的学生、刚入门飞控开发的工程师、以及想搞无人机仿真验证的爱好者这篇文章应该能把你们从“能飞”带到“知道为什么能飞”的那一步。2. 坐标系定义建模前的必修课90%的新手都挂在这一步2.1 机体坐标系与惯性坐标系到底怎么选做四旋翼建模第一步不是列方程而是把“在哪里描述运动”这件事说清楚。这里涉及两个坐标系地面惯性坐标系和机体坐标系。地面系一般取起飞点为原点X轴指向正北或任意固定方向Z轴垂直向下或向上满足右手定则即可。机体系则是固连在飞机上的坐标架原点在飞机重心X轴指向机头方向Y轴指向机右翼Z轴垂直机体向下。这里特别提醒一下不同教材Z轴取向可能相反一种是“北西天”坐标一种是“东北地”坐标推导公式前必须统一否则后面算力矩时符号全乱。为什么要两个坐标系配合用因为传感器测量的量分散在两个坐标系里加速度计和陀螺仪测得的是机体系下的物理量而GPS、光流传感器测得的是地面系下的位置。控制要的是地面系下的位移误差而执行机构电机是在机体系下出力。没有坐标变换这两个世界就是割裂的。2.2 旋转矩阵连接两个坐标系的桥从机体系到地面系的旋转一般用三个欧拉角来描述滚转角绕X轴记为、俯仰角绕Y轴记为、偏航角绕Z轴记为。三个单轴旋转矩阵相乘就得到了完整的旋转矩阵常用的顺序是Z-Y-X先偏航再俯仰最后滚转写成教科书里常见的形式$$ R R_z(\psi) \cdot R_y(\theta) \cdot R_x(\phi) $$展开之后是长这样的东西$$ R \begin{bmatrix} c\psi c\theta c\psi s\theta s\phi - s\psi c\phi c\psi s\theta c\phi s\psi s\phi \ s\psi c\theta s\psi s\theta s\phi c\psi c\phi s\psi s\theta c\phi - c\psi s\phi \ -s\theta c\theta s\phi c\theta c\phi \end{bmatrix} $$这个矩阵的意义很直观你把它乘上机体系下的任意向量得到的就是地面系下的表达。比如机体系下的升力方向是垂直机体平面的Z轴负方向乘以旋转矩阵后就能算出这个升力在地面系里沿三个轴的分量——这部分直接进入位移方程。我用一个生活类比帮助理解机身朝哪个方向就好比你斜端着一碗水地面系的Z轴是竖直向下而机体系的Z轴是垂直于碗底向外。水受重力是沿地面系Z轴的但你手托碗底感受到的力却是沿机体系Z轴的这中间的“角度差”就是旋转矩阵在数学上干的事情。2.3 欧拉角速率和机体角速度之间的转换关系建模里另一个容易让人头大的点是欧拉角随时间的变化率和机体系下的角速度并不相等。机体系下的角速度p, q, r是陀螺仪直接测量的量而我们要的姿态角导数是两者之间必须有转换矩阵$$ \begin{bmatrix} \dot{\phi} \ \dot{\theta} \ \dot{\psi} \end{bmatrix} \begin{bmatrix} 1 \sin\phi\tan\theta \cos\phi\tan\theta \ 0 \cos\phi -\sin\phi \ 0 \sin\phi/\cos\theta \cos\phi/\cos\theta \end{bmatrix} \begin{bmatrix} p \ q \ r \end{bmatrix} $$注意这个矩阵在俯仰角接近正负90度时出现奇异值——这就是欧拉角表示法的“万向节锁死”问题。做工程时如果飞机做特技飞行可能出现倒飞或大角度机动用欧拉角建模就会失效这时需要转成四元数表示。四元数没有奇异性运算也更快只是直观性稍差。实操建议初学者先用欧拉角理解物理意义写代码时尽早切四元数。我在PX4开源固件里见过姿态估计的代码内部清一色四元数运算只在最外层输出欧拉角给上层控制或地面站显示用得着。还有一个细节姿态解算的0.00001度的误差经过长时间积分就变成好几米的定位漂移这也是为什么建模时不能简单用“角度”本身而要用“角度的变化率”配合传感器融合算法来估计姿态。3. 力和力矩分量的建模每一项都对应飞机的真实物理部件3.1 拉力模型的主导项但计算远比想象复杂四旋翼垂直运动的本质是螺旋桨旋转产生升力。每个螺旋桨产生的拉力F与转速的平方成正比经验公式写作$$ F_i C_T \cdot \rho \cdot A \cdot R^2 \cdot \omega_i^2 $$如果飞机尺寸和飞行高度变化不大空气密度可以算作常数把所有常系数合并就得到简化的轴流模型$$ F_i k_T \cdot \omega_i^2 $$其中就是拉力系数。看起来很简单但这个系数怎么来是建模时第一个大坑。理论上可以用螺旋桨叶素理论推导但实际桨叶形状、桨距、翼型分布都是非线性的纯理论计算误差在10%以上。我实测过两个不同厂家的1045桨同样是3000转每分钟拉力差了近15%这直接影响后续仿真与实飞的一致性。更稳妥的方法是台架标定把电调和电机固定在力传感器上给定PWM值测转速和拉力用最小二乘法拟合k_T。具体方法后面会用一整节来讲这里先记住一个结论想靠模型精确拉力系数一定不能用仿真库里的默认值。3.2 反扭力矩和陀螺效应两个总被忽略却真实存在的角运动影响电机旋转除了产生升力还给机身一个反方向的扭矩。这个反扭力矩同样与转速平方成正比$$ M_i C_M \cdot \omega_i^2 k_M \cdot \omega_i^2 $$注意四旋翼的电机是成对反向旋转的两个顺时针、两个逆时针。当四个电机转速相等时反扭力矩整体互相抵消飞机纹丝不动一旦某对角电机转速差变大偏航力矩就出现了。这就是四旋翼偏航控制的基本原理。再看陀螺效应。高速旋转的螺旋桨相当于一个陀螺当机身旋转时桨叶会产生抵抗这个旋转的陀螺力矩。完整的高精度模型里每一项都是$$ M_{gyro,i} I_{prop} \cdot (\Omega_{body} \times \omega_i) $$也就是螺旋桨的角动量与机体角速度的叉积。很多人做悬停仿真的忽略这一项但在做快速横滚或剧烈的偏航机动时陀螺力矩会让飞机出现明显的“进动”现象姿态响应变得非常奇怪。3.3 空气阻力与风扰悬停仿真可以忽略高速飞行就必须建模空气在机身平面上产生的阻力用气动阻力模型来描述$$ F_{drag} \frac{1}{2} \rho C_D A v^2 $$这个力的方向与速度方向相反。低速小于3m/s悬停飞行时阻力项数值很小可以当作扰动处理但如果做高速巡航或抗风试验这个力就是决定稳态速度上限的关键项忽略它会导致仿真里的最大速度明显大于实际。还有一个容易漏掉的项是平动阻力对姿态回路的影响。机身高速飞行时气流会在机身上产生一个气动力矩这个力矩在某些攻角下会不稳定导致飞机抬头或低头。想要做精确的高速飞行模型就需要查表或风洞数据这部分已经属于气动辨识的专业领域了。4. 牛顿-欧拉方程推导从物理守恒定律到最终的六自由度模型4.1 平动动力学方程让位置变化与力直接挂钩现在可以写核心的动力学方程了。根据牛顿第二定律在地面惯性坐标系下四旋翼的平动满足$$ m \begin{bmatrix} \ddot{x} \ \ddot{y} \ \ddot{z} \end{bmatrix} R \cdot \begin{bmatrix} 0 \ 0 \ -T \end{bmatrix} \begin{bmatrix} 0 \ 0 \ mg \end{bmatrix} F_{drag} $$其中是四个螺旋桨升力之和那个带着旋转矩阵的项就是“机体系升力投影到地面系”的结果。重力项是沿地面系Z轴按我的坐标约定Z轴向下为正重力项为正。加上气动阻力项之后就构成了完整的位置动力学。这个方程解释了四旋翼所有的平移行为水平方向没有直接的输入力只能靠倾斜机身改变R矩阵的俯仰和滚转角让升力分解出水平分量来驱动运动。这就是为什么四旋翼的“前进”被戏称为“低头往前冲”。4.2 转动动力学方程姿态变化的根本原因姿态动力学用欧拉方程描述在机体系下$$ I \begin{bmatrix} \dot{p} \ \dot{q} \ \dot{r} \end{bmatrix} - \begin{bmatrix} p \ q \ r \end{bmatrix} \times I \begin{bmatrix} p \ q \ r \end{bmatrix} \begin{bmatrix} \tau_\phi \ \tau_\theta \ \tau_\psi \end{bmatrix} M_{gyro} $$这里I是转动惯量矩阵。对于四旋翼这类近似对称的结构可以近似为对角阵也就是I_xx、I_yy、I_zz三个独立参数。其中的叉积项代表陀螺力矩的耦合——当飞机横滚角速度大时即使没有输入俯仰力矩俯仰方向也会受到耦合影响。三个轴的力矩输入分别是横滚力矩来自左右电机升力差俯仰力矩来自前后电机升力差偏航力矩来自反扭力矩差。写成控制关系$$ \tau_\phi l \cdot k_T \cdot (\omega_2^2 - \omega_4^2) $$$$ \tau_\theta l \cdot k_T \cdot (\omega_1^2 - \omega_3^2) $$$$ \tau_\psi k_M \cdot (\omega_1^2 \omega_3^2 - \omega_2^2 - \omega_4^2) $$其中l是电机到飞机重心质心的水平距离。这里的电机编号要和你实际装机保持一致否则就是“控制分配矩阵填错飞控输出全反”的悲剧。4.3 从单机功率到模型整体六自由度的完整状态表达把平动方程和转动方程组合起来加上姿态角速率的转换方程就得到了完整的6自由度模型。用状态空间的形式写出来状态量是位置x、y、z速度、姿态角、和角速度p、q、r一共12个状态。完整模型写出来很长但表达的信息其实就一句话给定四个电机的转速就能推算出飞机的位置和姿态随时间的变化。这是仿真器最核心的动力学引擎。反过来想让飞机按指定轨迹飞行就需要设计控制律根据期望的位置和姿态反解所需的转速——这是后面控制器设计的任务。5. 实操建模实物参数测量与辨识技巧5.1 质量、重心与转动惯量的测量别再用CAD里的数模型里的十几个参数最影响真实度的不是公式形式而是数值准不准。很多初学者直接用SolidWorks里的估算值结果仿真和实机飞行状态差了一大截。我个人的经验是质量和重心一定要实测转动惯量至少用悬线法测一次。质量用电子秤直接称重心位置用两个秤支撑法做一个“找支点”的操作把飞机平放在一根细杆上调节细杆位置直到飞机保持平衡细杆位置所在平面即为重心平面。三个维度都用这个方法做一遍重心的空间位置就出来了。转动惯量的测法比较有趣最常用的是三线摆法用三根等长的线吊起飞机让飞机绕垂直轴做小角度扭转摆动测量摆动周期T根据公式$$ I_{zz} \frac{m g r^2}{4\pi^2 L} \cdot T^2 $$就能算出绕垂直轴的转动惯量。绕另外两个轴的测量需要重新调整挂架方向。如果条件不允许也可以用CAD算出的值做初值再飞一轮测试根据实际响应调整。5.2 拉力系数与扭矩系数的台架标定一个下午能做完标定k_T和k_M需要一个简易的台架力传感器、转速计或者电调遥测数据、直流电源、飞控板。把单个电机-电调-桨叶组合固定在力传感器上依次给不同油门值记录转速和拉力。我标定一个四轴一般测10个点从10%油门到90%油门每个点稳定采集5秒取平均值。再用线性回归去拟合和的关系。扭矩系数k_M的标定稍微麻烦一点需要测量电机座上的反扭力。我是用一个L型支架把电机轴向力传感器压在支架上记录侧向力乘上力臂得到扭矩然后同样拟合与转速平方的关系。这里分享一个注意电池电压的变化会显著影响拉力系数的标定结果。同一个电机3S电池满电12.6V和接近没电11.1V时相同转速下输出拉力并不会有太大差别因为转速是闭环控制但电调的输出特性会变。所以标定时尽量保持输入电压稳定或者记录电压值在后续控制里做电压补偿。5.3 混控矩阵把期望力与力矩变成四个电机的转速命令有了参数以后可以写出完整的控制分配矩阵。常规X型四旋翼期望的升力T、横滚力矩、俯仰力矩、偏航力矩与四个电机转速平方的关系是$$ \begin{bmatrix} T \ \tau_\phi \ \tau_\theta \ \tau_\psi \end{bmatrix} \begin{bmatrix} k_T k_T k_T k_T \ -k_T l k_T l k_T l -k_T l \ k_T l k_T l -k_T l -k_T l \ -k_M k_M -k_M k_M \end{bmatrix} \begin{bmatrix} \omega_1^2 \ \omega_2^2 \ \omega_3^2 \ \omega_4^2 \end{bmatrix} $$注意矩阵的具体数值和正负号完全取决于你的电机序号排列和旋向定义。装机之后第一件事就是拿着这个矩阵去对照你的飞机实际构型把每一行每一列都核对一遍。我见过不少飞友炸机的根源就是飞控里选的机型跟实际X型还是十字型不一致导致混控矩阵直接反了。6. 从数学到控制模型的具体工程化应用6.1 在仿真环境里验证你的模型给控制器当“训练场”模型建好之后第一件能做的事就是仿真。以我常用的MATLAB/Simulink为例可以把第4节的12个状态方程封装进一个S-Function或Stateflow模块里外面接上控制器作为闭环。这样可以在不上真机的情况下先验证控制器的稳定性、调节PID参数、做故障模拟某个电机突然失效。做仿真时重要的一点是加入噪声和误差传感器模型要叠加高斯白噪声执行器要有响应延迟和饱和限制否则仿真里效果完美的控制器搬上真机照样乱来。很多开源仿真器如AirSim、Gazebo内部也是用这类动力学模型只是把空气动力学项做得更细。6.2 模型在姿态控制器设计中的关键作用PID参数可以不用再瞎试有了模型姿态控制器的设计从“盲调”变成“有据可查”。以角度环为例外环是姿态角控制内环是角速度控制。根据模型中的转动方程角速度回路的被控对象可以近似成一阶惯性环节增益就是力矩系数除以转动惯量。这告诉你p、q、r三个角速度环比例增益的初值怎么定$$ K_p \approx \frac{I_{xx}}{\tau_\phi} $$用这个初值起调比纯靠脸盲调省了一个晚上的时间。我自己调试时通常会把模型仿真出来的增益直接作为真机PID的初值然后根据实际飞行做小范围微调一般两三轮就能稳住。6.3 模型验证与误差源分析为什么模型算出来和真机飞起来不一致建好模型不等于万事大吉。模型和真机之间的误差来源主要有几类。首先是参数误差质量估计不准、重心不在几何中心、转动惯量测得不准都会导致模型输出与实机响应偏差。这时可以做一个“模型校准飞行”记录实机舵面激励下的姿态响应再与仿真对比用对比结果反过来修正参数。其次是未建模动态比如机架的弹性形变、螺旋桨在高速下的气动失速、电机的响应延迟这些都很难用简单公式描述。处理思路是先把简单模型建好跑通再逐个叠加更精细的修正项千万不要一上来就搞CFD级别的建模工程上性价比极低。最后是延迟从控制器输出到电机转动再到产生拉力的过程存在几十到上百毫秒的延迟这个延迟会极大影响控制稳定性。模型里至少要加入一阶惯性环节模拟电机延迟否则仿真的相位裕度算出来和实际差距很大。7. 常见问题排查建模时容易踩的坑和解决办法7.1 坐标方向搞反升力符号不对飞机直接扎地板这个错误极其常见。很多人在旋转矩阵里用了“Z轴向上”的约定把升力方向写成正Z结果跟重力方向一致仿真里飞机直接加速砸向地面。检查方法非常简单把飞机水平放置欧拉角全为0给一个正升力看z轴加速度是不是向上负方向或正方向取决于你的坐标约定。如果一个简单的悬停命令都让飞机往地上加速八成是坐标方向或升力符号的问题。7.2 混控矩阵符号不对打杆方向全反翻滚就是翻机混控矩阵的符号错了表现是给一个正的横滚指令飞机往反方向滚。这类问题一定要在系留测试或仿真里先验证——把飞机绑在测试台架上给一个小的偏航指令观察四个电机的转速变化是否符合预期。我习惯在混控代码里打印每个电机的转速命令做一次“地面仿真输出检查”这比上天试错安全得多。7.3 转动惯量差太多姿态振荡的隐形原因如果你的姿态响应在调参时表现出奇怪的振荡——某些频率下特别容易共振无论怎么降PID增益都不行——大概率是转动惯量估计偏差太大。此时可以用频率响应法做辨识给一个扫频激励测量姿态响应通过波特图估算实际转动惯量。这一步做完再调PID你会感觉像是换了架飞机。8. 一些实操心得和个人的小建议写到这里模型推导的全部主线已经走完。最后分享几点我在实际项目中摸爬滚打出来的体会。关于建模精度别追求“完全精确”追求“够用且一致”。控制上需要的不是真实系统的完美复制而是一个能反映主要动态特性的简化模型。你舍掉陀螺力矩悬停性能照样很好你忽略空气阻力低速飞行问题也不大。但如果你用的模型跟实机的趋势都相反那再精妙的控制律也白搭。关于工具链我强烈建议用一套PythonNumPy/SciPy先快速搭建模型原型验证模型行为和控制器设计然后再移植到MATLAB/Simulink做正式仿真或直接生成C代码。Python处理矩阵运算和可视化都比MATLAB轻量不少而且后期接真机数据做参数辨识也方便。等你需要更复杂的仿真环境时再考虑Gazebo或AirSim不迟。关于模型文档化所有推导的变量定义、坐标约定、参数标定记录一定要形成文档并写明日期、环境温度、电压等级因为建模参数跟使用条件强相关。我吃过一个亏某次做仿真用了一套夏天标定的参数冬天温度低了电池内阻变大实际转速-拉力关系变了飞行性能和仿真对不上排查了一天才发现是参数标定的环境变了。四旋翼的数学模型推导说白了就是三步定义坐标系、列出力和力矩、代入牛顿欧拉方程。难的是在每个环节都不出错并把模型参数标得接近真实。只要把这篇文章里面提到的坑都绕过去后面的控制器设计、状态估计、路径规划再复杂也有一个靠谱的“物理底座”托着。