简介:本资源是一套面向嵌入式姿态解算开发者的完整工程实现,聚焦于基于四元数的扩展卡尔曼滤波(EKF)与双级卡尔曼滤波在STM32F103ZET6平台上的落地应用,专为解决MPU9250九轴传感器在动态环境中姿态漂移、温漂补偿不足等核心难题而设计,适用于无人机、机器人、可穿戴设备等对实时性与精度要求较高的姿态估计场景。压缩包共1031个文件,含549个C源码(含ARM CMSIS-DSP数学库调用如arm_dct4_init_f32.c)、244个头文件(定义状态模型、协方差矩阵结构及传感器校准参数)、64个依赖描述文件(.d)及编译输出文件(.o/.axf/.hex),整体大小36.6MB,体现完整的Keil MDK工程架构与底层驱动-算法-校准全链路实现。已有173人学习下载,读者可直接获取已适配FPU加速的四元数EKF主循环、双级滤波器解耦架构代码、MPU9250椭球拟合磁力计校准模块、多速率数据同步逻辑及I2C底层驱动,大幅降低从理论到嵌入式部署的门槛。 拿STM32F103ZET6加一块MPU9250做姿态解算,听起来像是入门的标配,但真正要把姿态做到稳、准、快,里面的门道比想象中多得多。这篇文章我把实际跑通的方案完整梳理了一遍:基于四元数的扩展卡尔曼滤波(EKF)做主姿态融合,再在外面套一级卡尔曼滤波做航向修正与平滑,也就是标题里说的双级卡尔曼结构。我会先把方案选型的逻辑讲透,再把滤波器推导、参数整定、STM32工程实现、实测数据全部展开,最后整理调试过程中踩过的一堆坑。无论你是要做自平衡车、四轴、机械臂云台,还是搞体感交互,这套框架都能直接拿来用。
1. 项目概述与核心思路
1.1 这个系统到底在解决什么问题
姿态解算的本质,是用传感器的数据去估计一个物体在三维空间中的朝向。MPU9250里面有一颗三轴陀螺仪、一颗三轴加速度计和一颗三轴磁力计,正好组成了九轴方案。陀螺仪积分角度很平滑,但零偏会随时间累积漂移;加速度计能感知重力方向,但受运动加速度干扰严重,而且无法给出航向;磁力计能提供绝对航向,但在室内、电机旁边容易被磁场干扰。
这三类传感器单独用都不靠谱,正确的做法是融合。问题在于怎么融合:算法复杂度要照顾STM32F103ZET6这颗72MHz主频、没有硬件浮点单元的Cortex-M3,精度又要尽量向大疆、飞控板看齐。这篇项目的核心思路就是用四元数表达姿态,用扩展卡尔曼滤波把陀螺仪动态数据和加速度计、磁力计的绝对参考结合起来,最后再用一级卡尔曼滤波专门处理航向修正,形成一个双级结构。
1.2 为什么选“四元数+EKF+双级卡尔曼”这个组合
传感器姿态融合的常用方案就几个:互补滤波、Mahony梯度下降、Madgwick算法、卡尔曼滤波家族。互补滤波实现简单,代码几十行,适合刚入门跑通Demo,但动态性能和抗干扰能力有限,在剧烈运动或者低帧率情况下姿态容易拉偏。Madgwick算法在校准好的情况下效果不错,计算量也小,但它是确定性算法,没有噪声模型,无法根据传感器质量动态调节信任度。
扩展卡尔曼滤波的想法很直接:既然陀螺仪积分在短时间内是可靠的,那就以陀螺仪数据驱动状态预测;而加速度计和磁力计的观测虽然各自有缺陷,但长期来看是绑定了真实姿态的,于是把二者作为量测持续修正预测结果。EKF的代价是需要推导非线性方程的雅可比矩阵,计算量比互补滤波大不少,但换来的是精度和噪声自适应能力。在STM32F103上,经过优化的EKF跑200Hz完全可行,这也是我实测过的结论。
至于为什么还要再加一级卡尔曼滤波,这是我在这类项目里踩过最深的一个坑。单级EKF如果同时把九轴数据都塞进去,动态响应和航向修正确实都能兼顾,但磁力计的干扰信号会通过量测更新直接污染整个姿态四元数。比如你在电机旁边转一下,roll和pitch会跟着抖一下,这就是典型的磁力计干扰串进了整个姿态。把航向修正单独拆出来放到第二级,EKF第一级只负责陀螺仪+加速度计的高频姿态融合,第二级再专门处理磁力计航向修正,两级各管各的,问题就彻底分开了。
1.3 方案选型对比表
| 方案 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|
| 互补滤波 | 代码量小、计算量低、容易调通 | 动态性能一般、抗干扰能力弱 | 入门学习、玩具级控制 |
| Madgwick | 计算量小、静态精度尚可 | 没有噪声建模、参数固定 | 资源极受限场景 |
| 单级EKF(九轴) | 精度高、能建模噪声 | 磁力计干扰污染所有姿态 | 磁场环境干净、无电机干扰 |
| 双级EKF+KF | 第一级专注横滚俯仰,第二级专注航向,抗干扰性强 | 结构稍复杂、参数调教工作量增加 | 有电机、磁场不稳定的移动平台 |
从表格可以看出来,双级结构不是性能指标上全面碾压单级EKF,而是在实际工程环境里多了一堵保护墙,让传感器故障或者干扰不至于一波带走整体姿态。这也是我最终选择它的根本原因。
2. 四元数基础与姿态表示
2.1 从欧拉角到四元数:为什么要换坐标系表示
很多初学者习惯了用pitch、roll、yaw描述姿态,那为什么不直接拿欧拉角做卡尔曼滤波?原因有两个。第一个是万向锁问题:当俯仰角接近±90°时,滚转和偏航会退化成同一个旋转轴,系统失去一个自由度,姿态解算会变得极不稳定。第二个问题更实际:欧拉角不适合做连续旋转积分,因为角速度不是欧拉角的简单导数,需要在更新过程中反复做三角函数运算,这么一套在72MHz的M3上跑实在不太划算。
四元数用一个标量加上三个虚部描述旋转,数学上相当于四维空间里的单位向量。它的最大好处是:没有奇异点,在任意姿态下都能稳定更新;乘法运算可以简洁地表示旋转复合,只需要一堆乘法和加法,没有三角函数;归一化操作可以持续修正数值误差,积分过程整体非常稳定。工程实现上,四元数也是现代飞控的主流选择。
2.2 坐标系约定与初始四元数计算
做姿态解算的第一步是统一坐标系。我用的当地坐标系是北东地(NED),X轴指向北,Y轴指向东,Z轴指向地心。机体坐标系是固定的:X轴朝机头,Y轴朝右翼,Z轴朝下。MPU9250输出的加速度和角速度都是机体坐标系下的量,而姿态四元数表达的是机体坐标到当地坐标系之间的旋转关系。
初始四元数通常需要根据传感器静止时刻的姿态来生成。如果开机时设备水平放置,初始欧拉角就是(0, 0, 0),对应的四元数是(1, 0, 0, 0)。如果设备开机时不是水平的,可以先用加速度计估算初始横滚和俯仰角:
roll = atan2(accel_y, accel_z) pitch = atan2(-accel_x, sqrt(accel_y^2 + accel_z^2))
再结合磁力计得出初始航向角yaw,然后用欧拉角转四元数的公式生成初始值:
q0 = cos(roll/2)*cos(pitch/2)*cos(yaw/2) + sin(roll/2)*sin(pitch/2)*sin(yaw/2) q1 = sin(roll/2)*cos(pitch/2)*cos(yaw/2) - cos(roll/2)*sin(pitch/2)*sin(yaw/2) q2 = cos(roll/2)*sin(pitch/2)*cos(yaw/2) + sin(roll/2)*cos(pitch/2)*sin(yaw/2) q3 = cos(roll/2)*cos(pitch/2)*sin(yaw/2) - sin(roll/2)*sin(pitch/2)*cos(yaw/2)这里要特别注意欧拉角的旋转顺序。不同行业习惯不一样,有的用ZYX、有的用ZXY,四元数转欧拉角的公式是按顺序约定的。我全程用的是Z-Y-X内旋顺序,也就是先偏航、再俯仰、最后横滚,这是无人机领域最主流的约定,转出来之后可以直接对接后续的控制环路。
2.3 四元数运动学方程与离散化实现
四元数的连续时间运动学方程是:
q_dot = 0.5 * q ⊗ ω
其中ω是机体坐标系下的角速度四元数 (0, ωx, ωy, ωz),⊗表示四元数乘法。工程上我们用的是离散时间更新,最常见的一阶近似是把q_dot乘以dt后加到原四元数上:
float w_norm = sqrt(wx*wx + wy*wy + wz*wz); float half_dt = 0.5f * dt; float q_new[4]; q_new[0] = q[0] - half_dt * (q[1]*wx + q[2]*wy + q[3]*wz); q_new[1] = q[1] + half_dt * (q[0]*wx + q[2]*wz - q[3]*wy); q_new[2] = q[2] + half_dt * (q[0]*wy - q[1]*wz + q[3]*wx); q_new[3] = q[3] + half_dt * (q[0]*wz + q[1]*wy - q[2]*wx);这个一阶近似在采样频率足够高(一般要求角增量小于30°)时精度完全够用。如果追求更高精度,可以用角增量积分或四阶龙格库塔,但在F103上没必要,一阶近似配合200Hz采样已经能把误差压得很低。
每次更新完后必须做一次归一化:把四个分量除以它们的模长。这是四元数算法最容易忽略的细节,单位四元数如果不归一化,累积误差会让姿态渐渐畸变,最终导致卡尔曼滤波发散。
3. 扩展卡尔曼滤波:状态方程与量测方程
3.1 为什么标准卡尔曼滤波不够用,必须扩展
标准卡尔曼滤波的假设是系统状态方程和观测方程都是线性的,也就是形如 x' = Ax + Bu + w 这种形式。但四元数运动学方程里包含了四元数和角速度的乘法,而加速度计观测姿态又是通过旋转矩阵计算出重力分量,方程里全是q1q3、q0q2这种二次项,完全是非线性关系。强行用线性卡尔曼会导致模型失配,滤波器会不稳定甚至发散。
扩展卡尔曼滤波的思路是:把非线性方程在当前估计值附近做一阶泰勒展开,求出雅可比矩阵,再用这个线性化后的模型走标准卡尔曼滤波流程。只要误差不大,一阶近似是足够准确的。姿态解算恰恰满足这个前提,因为EKF的预测循环以200Hz跑,相邻两帧之间姿态变化很小。
3.2 状态向量与误差模型设计
我采用的状态向量是7维:
x = [q0, q1, q2, q3, bgx, bgy, bgz]^T
前四个分量是姿态四元数,后三个是陀螺仪的三轴零偏。把零偏放进状态里去估计,是陀螺仪漂移校正的关键。陀螺仪测量值可以写成:
ω_measured = ω_true + b + v
其中b就是缓慢变化的零偏,v是测量白噪声。EKF在每次量测更新时,会根据加速度计观测到的姿态误差来调整零偏估计,这样陀螺仪的长期漂移就得到了补偿。
状态方程的离散形式是:
q_{k+1|k} = q_k + 0.5 * dt * q_k ⊗ (ω_measured - b_k) b_{k+1|k} = b_k
也就是说,四元数随着角速度积分而演化,零偏假设近似不变,由EKF通过观测持续修正。这个模型写作代码就是:
// 预测步骤 float wx = gx - bgx; float wy = gy - bgy; float wz = gz - bgz; float half_dt = 0.5f * dt; x[0] += half_dt * (-x[1]*wx - x[2]*wy - x[3]*wz); x[1] += half_dt * ( x[0]*wx + x[2]*wz - x[3]*wy); x[2] += half_dt * ( x[0]*wy - x[1]*wz + x[3]*wx); x[3] += half_dt * ( x[0]*wz + x[1]*wy - x[2]*wx); // bgx, bgy, bgz 保持不变3.3 量测方程与雅可比矩阵
加速度计在机体坐标系里测得的是重力向量在当地坐标系的投影经旋转矩阵转回机体系的值。当地重力向量(NED下)是 (0, 0, 1),所以观测方程是:
a_hat = R(q)^T * [0, 0, 1]
展开旋转矩阵后,三个分量为:
a_hat_x = 2*(q1q3 - q0q2) a_hat_y = 2*(q0q1 + q2q3) a_hat_z = q0^2 - q1^2 - q2^2 + q3^2
EKF的预测输出是q,我们拿这个q计算出来的a_hat应该与加速度计实测值a_meas接近。两者的偏差乘以卡尔曼增益就产生状态修正量。
为了执行卡尔曼更新,需要对观测向量关于状态向量求雅可比矩阵H。这里很容易算错,建议用符号推导工具验证。H矩阵是3行7列,前四列是对四个四元数分量的偏导,后三列是零偏相关项(加速度观测不直接依赖零偏,所以为0)。前四列的偏导形式是:
h11 = -2*q2 h12 = 2*q3 h13 = -2*q0 h14 = 2*q1 h21 = 2*q1 h22 = 2*q0 h23 = 2*q3 h24 = 2*q2 h31 = 2*q0 h32 = -2*q1 h33 = -2*q2 h34 = 2*q3这个雅可比矩阵是在每一次量测更新时用当前q实时计算的,不是固定的。正因为H矩阵随q变化,整个滤波过程才是“扩展”的线性化滤波。
3.4 噪声矩阵与协方差的整定
EKF有两个核心噪声矩阵:过程噪声矩阵Q和量测噪声矩阵R。Q代表状态方程的信任程度,R代表传感器的信任程度。Q越大,滤波器越相信量测、响应越快但噪声越大;Q越小,滤波器越平滑但滞后明显。R同理,方向相反。
经过实际调试,我用的一组比较稳定的初值是:
- Q对角线:四元数相关部分取 1e-4 到 1e-3,零偏相关部分取 1e-6 到 1e-5
- R对角线:加速度计三个轴取 1e-2 到 1e-1(不同量程下要缩放)
- 初始协方差矩阵P0:取单位阵放大到 1e-2 量级,代表初始姿态不确定度
为什么加速度计R取这么大?因为运动过程中加速度计不仅有重力分量,还叠加了运动加速度,运动加速度就是观测噪声的一部分。R取大一点,相当于告诉滤波器别太相信加速度计,用陀螺仪积分先扛住,这样在快速翻转时姿态不会跟着加速度剧烈抖动。后面我会讲如何基于这个思想做自适应R。
4. 双级卡尔曼滤波的设计
4.1 第一级和第二级各自管什么
双级卡尔曼的核心原则是责任分离。第一级EKF以200Hz频率运行,输入是陀螺仪和加速度计,输出是高频姿态四元数。它负责解决姿态的动态跟踪和横滚、俯仰精度。第二级卡尔曼以20到50Hz频率运行,输入是磁力计推算的航向角以及第一级输出的姿态,输出是经过滤波、平滑、修正后的最终航向。
为什么航向要单独拉出来处理?因为在EKF量测方程里直接加入磁力计,航向误差会通过增益传导到四元数的所有分量上。现实中电机、电源线都会产生磁场畸变,哪怕磁力计异常只持续几百毫秒,姿态也会被拉偏一度甚至更多。把航向修正放到第二级后,第一级的横滚俯仰完全不受磁力计干扰,磁力计异常最多影响航向,而航向异常又会被第二级卡尔曼的低通特性过滤掉一部分,整个系统稳定性天然高出一截。
4.2 第二级卡尔曼滤波器的具体形式
第二级滤波器可以设计成一个轻量的线性卡尔曼滤波器。状态量取两维:
x2 = [ψ, ω_ψ]
分别是航向角和航向角速度。状态方程是一个带角速率的随机游走模型:
ψ_{k+1} = ψ_k + ω_ψ * dt + w1 ω_ψ_{k+1} = ω_ψ_k + w2
观测是磁力计推算出的航向角 ψ_mag。量测方程是线性的:
z_k = [1, 0] * [ψ, ω_ψ]^T + v
这个滤波器是完全线性的,所以直接用标准卡尔曼滤波更新即可。它的作用有三个:一是平滑磁力计噪声带来的航向抖动;二是当磁力计出现野值时,通过卡方检验可以把异常观测挡在门外;三是可以把第一级输出的、带陀螺仪积分的航向与磁力计航向做最优加权,让航向既不漂移也不跟手迟钝。
4.3 两级之间的衔接与最终姿态输出
两级滤波最终如何拼装?第一级EKF输出的是四元数q_ekf,其中roll和pitch已经接近真实值,yaw会缓慢漂移。第二级输出的是修正后的航向角ψ_final。最终对外输出的姿态,是把q_ekf的四元数转成欧拉角,提取出roll和pitch,再用第二级的ψ_final替换掉q_ekf里的yaw,然后重新合成一个四元数。
这个“替换航向再合成四元数”的操作要从欧拉角到四元数的公式反推。具体流程是:
// 1. 从q_ekf提取roll和pitch roll = atan2(2*(q0*q1 + q2*q3), 1 - 2*(q1*q1 + q2*q2)) pitch = asin(2*(q0*q2 - q3*q1)) // 2. 用第二级滤波后的ψ_final替代原始yaw // 3. 用新的roll, pitch, ψ_final重新合成四元数这样做的好处是最终四元数既包含了EKF对横滚俯仰的精确估计,又包含了第二级对航向的长期校正,同时保证四元数依然是单位四元数,方便对后续控制环路直接使用。
4.4 两级滤波器参数协作的实操经验
双级结构调参时,第一级和第二级是相对独立的。通常先把第一级EKF调到横滚俯仰在动态旋转时响应灵敏、静态时平滑稳定,然后再去调第二级航向滤波。第二级的Q和R相对好理解:Q大一点航向跟随快但噪声大,R大一点航向平滑但滞后明显。一个典型的起始值是:
- 第二级Q对角线取 [0.001, 0.001]
- 第二级R取 5.0(对应磁力计航向在强磁干扰下的噪声方差)
在实际调试中,你用MCU的串口把第一级的yaw原始值和第二级的yaw滤波值都打印出来对比,就能直观看到航向曲线从“在真实值附近抖动”变成“平滑收敛”。这个调参过程是整个双级方案最花时间的地方,但没有捷径,只能一遍一遍调Q、R的比值。
5. STM32F103ZET6 + MPU9250 的工程实现
5.1 硬件连接与通信方式选择
STM32F103ZET6是LQFP144封装的Cortex-M3,外设资源非常充足,USB、CAN、多路串口、SPI、I2C全都有,跑姿态解算完全浪费不完。MPU9250支持I2C和SPI两种接口。我的建议是优先用SPI连接IMU的加速度计和陀螺仪,因为SPI的速率更高(我实际配置在1MHz)且没有I2C的ACK时序抖动,读取一致性更好。
这里有一个新手很容易踩的坑:MPU9250内部的磁力计AK8963是通过内部I2C通路访问的,当SPI模式使用时,磁力计数据需要通过MPU9250的I2C master模式或者旁路模式来读取。也就是说,你在SPI总线上读出的是MPU9250内部的ACCEL和GYRO数据,磁力计数据是通过MPU9250帮你“转发”出来的。实际代码里要把MPU9250设置为bypass模式,让外部主控可以直接用I2C读取AK8963的寄存器,或者用SPI往MPU9250发命令,由它通过内部I2C读取磁力计再通过SPI返回。这个细节如果不注意,很可能会卡在磁力计数据一直为零的问题上。
5.2 驱动初始化与传感器校准流程
传感器驱动初始化的顺序不能乱。我总结了一套稳定的流程:
- 复位MPU9250,等待100ms(写PWR_MGMT_1的reset位)
- 配置电源管理寄存器,选择内部时钟源
- 配置陀螺仪量程为±1000dps,加速度计量程为±2g
- 配置数字低通滤波器(DLPF),带宽设为差不多20Hz
- 配置SPI接口为单传感器模式(I2C禁止,但保留磁力计的bypass通路)
- 磁力计独立初始化,连续测量模式,16位输出
- 传感器静止校准
校准的核心操作是静态求平均。把设备静置在水平桌面,连续采集陀螺仪和加速度计各500个样本,然后取平均,得到陀螺仪零偏初始值和加速度计零偏(顺带判断传感器安装是否水平)。陀螺仪零偏会随温度漂移,所以EKF里面仍然需要把零偏作为状态估计,校准只是给EKF一个不起眼但靠谱的初值。
5.3 主循环任务的调度与时序设计
姿态解算是一个固定频率的实时任务,在F103上建议用定时器中断驱动,而不是在while(1)里读取+计算。我用TIM4产生1kHz的中断,中断里做计时,每5个中断触发一次200Hz的传感器读取和EKF更新,每20个中断触发一次50Hz的磁力计读取和第二级航向滤波。
这种分层调度的好处是:EKF的预测更新有固定的dt,不会因为主循环里其他逻辑阻塞而时间漂移。卡尔曼滤波最怕的就是时间步长不一致,一旦dt抖动,Q矩阵的等效放大倍数就变了,滤波器会表现出时好时坏的行为。
5.4 关键代码框架
下面是主循环和两级滤波的骨架,可以对照自己的工程改造:
#define IMU_READ_HZ 200 #define MAG_READ_HZ 50 float x[7]; // EKF状态:q0-q3, bgx, bgy, bgz float P[7*7]; // EKF协方差矩阵 float x2[2]; // 第二级状态:yaw, yaw_rate float P2[2*2]; // 第二级协方差 float q_final[4]; // 最终输出四元数 void TIM4_IRQHandler(void) { static uint16_t imu_cnt = 0; static uint16_t mag_cnt = 0; if (imu_cnt++ >= 5) { imu_cnt = 0; imu_read_accel_gyro(&ax, &ay, &az, &gx, &gy, &gz); ekf_predict(gyro, dt); ekf_update_accel(accel); normalize_q(x); if (mag_cnt >= 4) { mag_cnt = 0; float yaw_mag = get_yaw_from_mag(q, mag); kf_yaw_update(yaw_mag); } fusing_output_quaternion(q_final, x, x2); } }如果想减少代码规模,要特别注意避免在主循环里用动态内存和malloc,嵌入式卡尔曼滤波器用固定数组足够,动态分配会导致内存碎片,也可能拖慢执行速度。
5.5 计算量评估与浮点优化
STM32F103ZET6没有FPU,所有浮点运算都是以软件浮点形式执行。7状态EKF的一轮预测加量测更新,包括四元数更新、雅可比矩阵计算、矩阵乘法和求逆,我统计过大约需要1.2ms左右。在200Hz下每5ms执行一次,占用大约24%的CPU时间,后续再跑PID控制、通信、日志打印都绰绰有余。
如果没有加优化直接编译,可能单次EKF要跑到2.5ms,这会挤占其他任务的时间。几个关键优化手段是:
- 编译器开-O2优化,这是最划算的优化
- 所有变量声明为float,避免double,double在M3上是模拟运算,速度慢一倍以上
- 矩阵运算用固定展开,不要用for循环处理小矩阵,比如3x3矩阵乘法直接写出9个表达式
- 雅可比矩阵中很多项重复,可以存成临时变量复用
实测下来,优化后单次EKF能压到1ms以内。这个优化前后的对比非常明显,建议读者在自己板子上跑一遍基准测试,你会明显感受到为什么嵌入式里不能像PC上那样随意写代码。
6. 实测结果与分析
6.1 静态漂移测试
把整个系统静止放在桌面上,记录5分钟的姿态输出。纯陀螺仪积分模式下,横滚和俯仰的漂移很容易就超过2°,航向漂移更是分钟级别可达十几度。而双级EKF输出的姿态在5分钟内横滚和俯仰稳定在±0.3°以内,航向漂移在±0.8°以内。这组数据说明陀螺零偏估计和磁力计航向修正都正常工作。
6.2 动态测试:快速旋转与剧烈加速
用手拿着板子快速转90°、机身翻滚、快速推拉,观察姿态是否跟手和是否反弹。纯互补滤波在快速旋转时容易出现“转不回来”的滞后现象,因为加速度计的低通特性拖慢了响应。单级EKF动态响应好很多,但在电机启动时磁场突变会导致所有轴都被拉偏一点。双级结构下,横滚俯仰不受磁场突变影响,航向虽受影响但被二级滤波抑制,表现最稳。
| 测试场景 | 纯陀螺积分 | 互补滤波 | 单级EKF | 双级EKF+KF |
|---|---|---|---|---|
| 静态漂移(5分钟) | 大 | 中等 | 较小 | 最小 |
| 快速旋转响应 | 跟手但漂移大 | 明显滞后 | 跟手 | 跟手 |
| 磁场突变 | 不影响 | 影响小 | 横滚俯仰被拉偏 | 只有航向微抖 |
| 剧烈运动加速度干扰 | 不受影响 | 明显过冲 | 轻微过冲 | 轻微过冲 |
6.3 对第一级第二级的观察
调试过程中最有意思的是把第一级EKF输出的yaw和第二级输出的yaw同时打出来对比。第一级的yaw在来回摆动时很跟手,但是会缓慢飘,磁力计偶尔出现的相位跳变也会让它抖一下。第二级的yaw是一条很平滑的曲线,长期来看不飘,短期来看响应滞后大概在100ms以内,这个滞后对大多数机器人控制完全可接受。这就是双级结构的效果:速度由第一级负责,长期精度由第二级负责,两级合起来才是一个完整的姿态。
7. 常见问题与避坑指南
7.1 陀螺仪零偏随时间漂移
现象:静态时pitch和roll慢慢往一个方向偏,速度不快但一直在动。原因就是温度升高导致陀螺零偏缓慢变化。解决办法是EKF状态里加零偏估计项,这是我在状态向量里包含bgx、bgy、bgz的原因。同时可以每几十秒在系统静止时自动触发一次零偏更新,进一步压住长时间漂移。
7.2 磁力计硬磁干扰与软磁干扰
磁力计放在PCB上,附近的铁磁元器件和电源线都会产生干扰。硬磁干扰表现为磁力计数据在空间中形成一个偏心的球,软磁干扰让球变成椭球。简单的做法是拿着板子绕“8”字旋转,采集一批数据,拟合出中心偏移(硬磁补偿)和缩放矩阵(软磁补偿)。如果不想做得太重,至少要做硬磁补偿,把三轴偏移值减去。没有补偿的航向数据在旋转时误差可能达到十几度。
7.3 四元数忘记归一化导致发散
这是新手最容易犯的错误。四元数如果长期不归一化,模长会逐渐偏离1,旋转矩阵会变形,最终姿态彻底乱掉。归一化的操作必须在每次四元数更新之后、任何输出之前执行。看起来是小事,但影响是致命的。我习惯在EKF更新函数最后强制调用归一化,同时断言模长误差不超过1e-6级别。
7.4 卡尔曼滤波的P矩阵发散
如果协方差矩阵P的元素变得特别大或者出现负值,滤波器就会出现数值发散。常见原因是Q矩阵设得过大,导致估计误差协方差持续膨胀,或者浮点运算中出现了NaN。另一个常见原因是量测更新时R矩阵过小,导致增益计算中除数为接近零的数。排查方法是在串口里周期打印P矩阵对角线元素检查合理性,一旦发现异常立即寻找是哪一步引入的问题。
7.5 MPU9250 SPI模式下磁力计读不出数据
这个问题在第五部分提过,这里再强调一次:MPU9250的磁力计AK8963挂在内部的I2C总线上,SPI模式下外部主控不能直接读磁力计寄存器,必须通过MPU9250的I2C Master或bypass模式访问。很多移植失败的工程都卡在这里。最简单的做法是把MPU9250的I2C_MST_EN关闭、I2C_BYPASS_EN打开,让外部主控通过I2C接口直连AK8963,这样虽然要开两路总线,但稳定性最好。
7.6 采样频率不稳定导致姿态抖动
卡尔曼滤波依赖固定的dt。如果采样时间不稳定,预测步噪声方差会被放大或缩小,姿态会表现出随机抖动。解决方案是用定时器中断驱动读取,而不是在串口回调里读取。我在F103上用TIM4产生固定1kHz中断计数,确保每次EKF更新的dt都严格等于5ms。
我把这些常见问题整理成一个速查表,方便遇到问题时快速对照:
| 现象 | 可能原因 | 检查与解决办法 |
|---|---|---|
| 姿态静止时缓慢漂移 | 陀螺零偏变化 | 检查EKF状态是否有零偏估计,重校准 |
| 姿态快速旋转时过冲 | 加速度计R过小 | 增大R,或做自适应R |
| 姿态缓慢滞后 | 加速度计R过大或Q过小 | 减小R,增大Q |
| 航向抖动大 | 磁力计噪声大、未校准 | 做硬磁软磁校准,加大二级滤波R |
| 卡尔曼发散、数值爆炸 | P矩阵异常、NaN | 检查归一化、Q/R比值、打印P对角线 |
| 航向在某个角度跳变 | 欧拉角转换顺序不一致 | 统一旋转顺序为ZYX,并核对四元数公式 |
8. 实操心得与扩展方向
这套双级卡尔曼滤波方案我先后在两个平台上验证过,一次是四轮小车,一次是二轴云台。两次经历最大的体会是:滤波器设计的核心不是算法越复杂越好,而是每一级都能用最合适的手段解决自己那部分问题。第一级的EKF解决了动态跟踪和横滚俯仰修正,第二级的KF解决航向长期精度和异常扰动隔离,各司其职,调参时思路也清晰很多。
在实际调试中,我给自己的建议是先把EKF单独调到满意,再加第二级。很多朋友一上来就想双级一把调通,结果两级参数互相干扰,根本分不清是哪里出了问题。串口调试输出的帮助很大,我建议把第一级的原始四元数、欧拉角、第二级的航向、磁力计原始航向都打印出来,用可视化工具画成曲线,这样调试效率比对着串口数字要高一个量级。
再分享一个小技巧:EKF的自适应R矩阵并不难加。在每个量测更新周期,计算加速度计模长与重力模长的偏差,如果偏差超过一个阈值,就把R放大几倍。这样在剧烈运动时滤波器自动降低对加速度计的信任,回到静止时再恢复正常R,姿态抖动会明显减少。这个改动大概只增加不到20行代码,但对动态场景的收益很可观。
后续想扩展的话,几个方向都是现成的。这套双级结构的第二级滤波器可以直接替换成融合GPS航向或视觉里程计的最优估计器,第一级也可以无缝移植到无刷云台的电机控制上。四轴飞行器的话,直接把最终四元数输出给姿态控制PID环即可。总之这个骨架是可以持续往上搭东西的,核心的滤波逻辑和工程经验不会过时。
本文还有配套的精品资源,点击获取