六轴融合算法

先说什么叫六轴融合?

在3Dof姿态追踪功能中,最主要的传感器就是陀螺仪(Gyroscope),它可以提供3个轴的角加速度,对时间进行积分,就可以得出物体旋转的方向角度。但是因为硬件精度等各方面原因,会产生误差,随着时间的累积,计算得到的角度误差就会越来越大,即产生漂移。

为了防止漂移,这就引入了另一个传感器,加速度计(Accelerometer)。在一般的3Dof运动中,由运动产生的加速度较少,对物体影响最大的是重力,而重力始终是垂直向下。引入Acc的本质就是利用重力来做参考,对Gyro计算的数据进行校正。

所谓六轴融合,就是Gyro三个轴+Acc三个轴,计算得到3Dof姿态。

值得注意的是:

1. 重力只能校准Roll和Pitch方向,Yaw方向依然会有误差,还要再引入地磁计(Magnetic),那就是九轴融合算法了。在Android上,SensorManager提供了注册监听的类型,其中TYPE_GAME_ROTATION_VECTOR就是用的六轴融合,TYPE_ROTATION_VECTOR用的就是九轴融合

2. 因为硬件本身的精度,还有设备组装等方面的原因,传感器必然存在原始的误差,需要做零漂校准。也就是尽量保证设备在水平静置状态下,Gyro数据都是0,Acc只有垂直水平面的轴数据为重力值,其余为0。

(很惭愧,这两样我都没有做过)

下面是摘自网上的一个简单算法,还有好多不懂的地方,先记下来,慢慢再研究。

public class IMUTest {
    private final float Kp = 4.0f;   //Kp越大代表重力校正的权重越大,过大可能产生抖动。
    private final float Ki = 0.001f; //Ki应该是零漂校准产生的,没仔细研究过还不确定,先用个较小的值,没啥影响
    private float q0, q1, q2, q3; //依次对应四元数的w,x,y,z
    private float exInt, eyInt, ezInt;
    private long lastTime;

    public IMUTest() {
        q0 = 1.f;
        q1 = 0.f;
        q2 = 0.f;
        q3 = 0.f;
        exInt = 0.f;
        eyInt = 0.f;
        ezInt = 0.f;
    }

    public void update(float GyroX, float GyroY, float GyroZ, float AccX, float AccY, float AccZ, long timeMS) {
        long diff = timeMS - lastTime;
        lastTime = timeMS;
        float halfT = diff * 0.001f * 0.5f;//单位毫秒转为秒

        //加速度计处于自由落体状态时,不进行姿态解算,因为三个轴都为0,会产生分母无穷大的情况
        if ( Math.abs(AccX * AccY * AccZ) < 0.000001f) return;

        float norm; //矢量的模或四元数的模

        //归一化加速度计,这样变更量程也不需要修改Kp参数
        norm = (float)Math.sqrt(AccX * AccX + AccY * AccY + AccZ * AccZ);
        float ax = AccX / norm;
        float ay = AccY / norm;
        float az = AccZ / norm;
        //用当前姿态计算出重力在三个轴上的分量,
        //参考坐标n系转化到载体坐标b系,用四元数表示的方向余弦矩阵第三列即是,为啥?
        float vx = 2.f * (q1 * q3 - q0 * q2);
        float vy = 2.f * (q0 * q1 + q2 * q3);
        float vz = q0 * q0 - q1 * q1 - q2 * q2 + q3 * q3;
        //计算传感器测量重力与姿态计算的重力之间的误差,叉乘可以表示这一误差,又是为啥?
        float ex = ay * vz - az * vy;
        float ey = az * vx - ax * vz;
        float ez = ax * vy - ay * vx;
        //对误差进行积分
        exInt = exInt + ex * Ki;
        eyInt = eyInt + ey * Ki;
        ezInt = ezInt + ez * Ki;
        //将误差PI后补偿到陀螺仪,即补偿零点漂移
        //这里的水平旋转方向由于没有观测者进行矫正会产生漂移,表现出来的就是积分自增或自减,所以应该直接用GyroX?
        float gx = GyroX + Kp * ex + exInt;
        float gy = GyroY + Kp * ey + eyInt;
        float gz = GyroZ + Kp * ez + ezInt;

        //下面进行姿态的更新,也就是四元数微分方程的求解
        float q0temp = q0;
        float q1temp = q1;
        float q2temp = q2;
        float q3temp = q3;
        //采用一阶毕卡解法,相关知识可参见《惯性器件与惯性导航系统》P212
        q0 = q0temp + (-q1temp * gx - q2temp * gy - q3temp * gz) * halfT;
        q1 = q1temp + (q0temp * gx + q2temp * gz - q3temp * gy) * halfT;
        q2 = q2temp + (q0temp * gy - q1temp * gz + q3temp * gx) * halfT;
        q3 = q3temp + (q0temp * gz + q1temp * gy - q2temp * gx) * halfT;
        //对四元数做归一化,单位化四元数在空间旋转时不会拉伸,仅有旋转角度,这类似线性代数里的正交变换
        norm = (float)Math.sqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3);
        q0 /= norm;
        q1 /= norm;
        q2 /= norm;
        q3 /= norm;
    }
}

参考自https://blog.csdn.net/asr9k/article/details/53258308?utm_medium=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.control&dist_request_id=ccb9df40-57cb-4714-873f-727e69a96dee&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.control

  • 0
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
根据提供的引用内容,EKF C语言六轴姿态融合算法是一种使用扩展卡尔曼滤波器(EKF)进行六轴姿态融合算法。该算法使用加速度计和磁力计来估计姿态,并使用EKF来融合两个感器的数据。在加速度计为0的条件下,该算法可以准确地估计姿态。在后续的融合算法中,加速度计和磁力计只起到修正的作用,影响不大。 以下是一个简单的EKF C语言六轴姿态融合算法的示例: ```c #include <stdio.h> #include <math.h> #define PI 3.14159265358979323846 // 定义加速度计和磁力计的数据结构 typedef struct { double x; double y; double z; } Vector3; // 定义姿态的数据结构 typedef struct { double roll; double pitch; double yaw; } Attitude; // 定义EKF的数据结构 typedef struct { double q0; double q1; double q2; double q3; double p[3]; double k[4]; double h[4]; double f[16]; double q[16]; double r[4]; double x[4]; double z[4]; double y[4]; double s[4]; double t[16]; } EKF; // 初始化EKF void EKF_Init(EKF *ekf) { ekf->q0 = 1.0; ekf->q1 = 0.0; ekf->q2 = 0.0; ekf->q3 = 0.0; ekf->p[0] = 0.0; ekf->p[1] = 0.0; ekf->p[2] = 0.0; ekf->k[0] = 0.0; ekf->k[1] = 0.0; ekf->k[2] = 0.0; ekf->k[3] = 0.0; ekf->h[0] = 0.0; ekf->h[1] = 0.0; ekf->h[2] = 0.0; ekf->h[3] = 0.0; ekf->f[0] = 1.0; ekf->f[1] = 0.0; ekf->f[2] = 0.0; ekf->f[3] = 0.0; ekf->f[4] = 0.0; ekf->f[5] = 1.0; ekf->f[6] = 0.0; ekf->f[7] = 0.0; ekf->f[8] = 0.0; ekf->f[9] = 0.0; ekf->f[10] = 1.0; ekf->f[11] = 0.0; ekf->f[12] = 0.0; ekf->f[13] = 0.0; ekf->f[14] = 0.0; ekf->f[15] = 1.0; ekf->q[0] = 0.01; ekf->q[1] = 0.0; ekf->q[2] = 0.0; ekf->q[3] = 0.0; ekf->q[4] = 0.0; ekf->q[5] = 0.01; ekf->q[6] = 0.0; ekf->q[7] = 0.0; ekf->q[8] = 0.0; ekf->q[9] = 0.0; ekf->q[10] = 0.01; ekf->q[11] = 0.0; ekf->q[12] = 0.0; ekf->q[13] = 0.0; ekf->q[14] = 0.0; ekf->q[15] = 0.01; ekf->r[0] = 0.1; ekf->r[1] = 0.0; ekf->r[2] = 0.0; ekf->r[3] = 0.1; ekf->x[0] = 1.0; ekf->x[1] = 0.0; ekf->x[2] = 0.0; ekf->x[3] = 0.0; } // 更新EKF void EKF_Update(EKF *ekf, Vector3 acc, Vector3 mag, double dt) { double ax = acc.x; double ay = acc.y; double az = acc.z; double mx = mag.x; double my = mag.y; double mz = mag.z; double q0 = ekf->q0; double q1 = ekf->q1; double q2 = ekf->q2; double q3 = ekf->q3; double p0 = ekf->p[0]; double p1 = ekf->p[1]; double p2 = ekf->p[2]; double k0 = ekf->k[0]; double k1 = ekf->k[1]; double k2 = ekf->k[2]; double k3 = ekf->k[3]; double h0 = ekf->h[0]; double h1 = ekf->h[1]; double h2 = ekf->h[2]; double h3 = ekf->h[3]; double f0 = ekf->f[0]; double f1 = ekf->f[1]; double f2 = ekf->f[2]; double f3 = ekf->f[3]; double f4 = ekf->f[4]; double f5 = ekf->f[5]; double f6 = ekf->f[6]; double f7 = ekf->f[7]; double f8 = ekf->f[8]; double f9 = ekf->f[9]; double f10 = ekf->f[10]; double f11 = ekf->f[11]; double f12 = ekf->f[12]; double f13 = ekf->f[13]; double f14 = ekf->f[14]; double f15 = ekf->f[15]; double q0q0 = q0 * q0; double q0q1 = q0 * q1; double q0q2 = q0 * q2; double q0q3 = q0 * q3; double q1q1 = q1 * q1; double q1q2 = q1 * q2; double q1q3 = q1 * q3; double q2q2 = q2 * q2; double q2q3 = q2 * q3; double q3q3 = q3 * q3; double norm; double hx; double hy; double hz; double bx; double bz; double vx; double vy; double vz; double wx; double wy; double wz; double ex; double ey; double ez; double qa; double qb; double qc; double qd; double s0; double s1; double s2; double s3; double s4; double s5; double s6; double s7; double s8; double s9; double s10; double s11; double s12; double s13; double s14; double s15; double t0; double t1; double t2; double t3; double t4; double t5; double t6; double t7; double t8; double t9; double t10; double t11; double t12; double t13; double t14; double t15; // 计算加速度计和磁力计的模 norm = sqrt(ax * ax + ay * ay + az * az); ax /= norm; ay /= norm; az /= norm; norm = sqrt(mx * mx + my * my + mz * mz); mx /= norm; my /= norm; mz /= norm; // 计算磁力计的方向 hx = mx * q0q0 - 2.0 * q0 * my * q3 + 2.0 * q0 * mz * q2 + mx * q1q1 + 2.0 * q1 * my * q2 + 2.0 * q1 * mz * q3 - mx * q2q2 - my * q1q1 + 2.0 * q2 * mz * q3 - mz * q0q0 - mz * q1q1; hy = 2.0 * q0 * mx * q3 + my * q0q0 - 2.0 * q0 * mz * q1 + 2.0 * q1 * mx * q2 - my * q2q2 + my * q3q3 + 2.0 * q2 * mz * q3 - mz * q1q1 - mz * q2q2; hz = 2.0 * q0 * my * q1 + 2.0 * q0 * mz * q0 + mz * q1q1 - my * q2q2 + 2.0 * q1 * mx * q3 - mz * q3q3 + 2.0 * q2 * my * q3 - mx * q1q1 - mx * q2q2; bx = sqrt((hx * hx) + (hy * hy)); bz = hz; // 计算预测状态 vx = 2.0 * (q1q3 - q0q2); vy = 2.0 * (q0q1 + q2q3); vz = q0q0 - q1q1 - q2q2 + q3q3; wx = 2.0 * bx * (0.5 - q2q2 - q3q3) + 2.0 * bz * (q1q3 - q0q2); wy = 2.0 * bx * (q1q2 - q0q3) + 2.0 * bz * (q0q1 + q2q3); wz = 2.0 * bx * (q0q2 + q1q3) + 2.0 * bz * (0.5 - q1q1 - q2q2); ex = (ay * vz - az * vy) + (my * wz - mz * wy); ey = (az * vx - ax * vz) + (mz * wx - mx * wz); ez = (ax * vy - ay * vx) + (mx * wy - my * wx); ekf->x[0] += dt * (q1 * ex + q2 * ey + q3 * ez); ekf->x[1] += dt * (q0 * ex + q2 * ez - q3 * ey); ekf->x[2] += dt * (q0 * ey - q1 * ez + q3 * ex); ekf->x[3] += dt * (q0 * ez + q1 * ey - q2 * ex); // 计算雅可比矩阵 s0 = q0q0 + q1q1 - q2q2 - q3q3; s1 = 2.0 * (q1 * q2 + q0 * q3); s2 = 2.0 * (q1 * q3 - q0 * q2); s3 = 2.0 * (q1 * q2 - q0 * q3); s4 = q0q0 - q1q1 + q2q2 - q3q3; s5 = 2.0 * (q2 * q3 + q0 * q1); s6 = 2.0 * (q1 * q3 + q0 * q2); s7 = 2.0 * (q2 * q3 - q0 * q1); s8 = q0q0 - q1q1 - q2q2 + q3q3; t0 = 2.0 * (-q2q2 - q3q3); t1 = 2.0 * (q1q2 - q0q3); t2 = 2.0 * (q1q3 + q0q2); t3 = 2.0 * (q1q2 + q0q3); t4 = 2.0 * (-q1q1 - q3q3); t5 = 2.0 * (q

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值