《新手四旋翼算法总结》由会员分享，可在线阅读，更多相关《新手四旋翼算法总结（6页珍藏版）》请在金锄头文库上搜索。
1、 新手四旋翼算法总结一姿态结算（匿名版程序）首先，程序中一般用了两种求解姿态的方法，一种为欧拉角法，一种为四元数法（1）欧拉角法静止状态，或者总加速度只是稍微大于 g 时，由加计算出的值比较准确。使用欧拉角表示姿态，令 , 和 代表 ZYX 欧拉角，分别称为偏航角、俯仰角和横滚角 。 载体坐标系下的 加 速 度(a xB,ayB,azB)和参考坐标系下的加速度 (axN, ayN, azN)之间的关系可表示为(1)。其中 c 和 s 分别代表 cos 和 sin。axB,ayB,azB 就是 mpu 读出来的三个值。这个矩阵就是三个旋转矩阵相乘得到的，因为矩阵的乘法可以表示旋转。(1) axc
2、ssaxNyBcs cyz z 飞行器处于静止状态，此时参考系下的加速度等于重力加速度，即（2）0xNyzag把（2）代入（1）可以解的（3）2()xByzarctg（4）yBzarctg即为初始俯仰角和横滚角，通过加速度计得到载体坐标系下的加速度即可将其解出，偏航角可以通过电子罗盘求出。（2）四元数法（通过处理单位采样时间内的角增量(mpu 的陀螺仪得到的就是角增量)，为了避免噪声的微分放大，应该直接用角增量-抄的书）上匿名的程序void IMUupdate(float gx, float gy, float gz, float ax, float ay, float az)float no
3、rm;/ float hx, hy, hz, bx, bz;float vx, vy, vz;/ wx, wy, wz;float ex, ey, ez; / 先把这些用得到的值算好float q0q0 = q0*q0;float q0q1 = q0*q1;float q0q2 = q0*q2;/ float q0q3 = q0*q3;float q1q1 = q1*q1;/ float q1q2 = q1*q2;float q1q3 = q1*q3;float q2q2 = q2*q2;float q2q3 = q2*q3;float q3q3 = q3*q3;if(ax*ay*az=0)r
4、eturn;norm = sqrt(ax*ax + ay*ay + az*az); /acc 数据归一化ax = ax /norm;ay = ay / norm;az = az / norm;/ estimated direction of gravity and flux (v and w) vx = 2*(q1q3 - q0q2); /四元素中 xyz 的 vy = 2*(q0q1 + q2q3);vz = q0q0 - q1q1 - q2q2 + q3q3 ;/ error is sum of cross product between reference direction of fi
5、elds and direction measured by sensorsex = (ay*vz - az*vy) ; /向量外积在相减得到差分就是误差ey = (az*vx - ax*vz) ;ez = (ax*vy - ay*vx) ;exInt = exInt + ex * Ki; /对误差进行积分eyInt = eyInt + ey * Ki;ezInt = ezInt + ez * Ki;/ adjusted gyroscope measurementsgx = gx + Kp*ex + exInt; /将误差 PI后补偿到陀螺仪，即补偿零点漂移gy = gy + Kp*ey +
6、eyInt;gz = gz + Kp*ez + ezInt; /这里的 gz 由于没有观测者进行矫正会产生漂移，表现出来的就是积分自增或自减 / integrate quaternion rate and normalise /四元素的微分方程q0 = q0 + (-q1*gx - q2*gy - q3*gz)*halfT;q1 = q1 + (q0*gx + q2*gz - q3*gy)*halfT;q2 = q2 + (q0*gy - q1*gz + q3*gx)*halfT;q3 = q3 + (q0*gz + q1*gy - q2*gx)*halfT;/ normalise quate
7、rnionnorm = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);q0 = q0 / norm;q1 = q1 / norm;q2 = q2 / norm;q3 = q3 / norm;/Q_ANGLE.Yaw = atan2(2 * q1 * q2 + 2 * q0 * q3, -2 * q2*q2 - 2 * q3* q3 + 1)* 57.3; / yawQ_ANGLE.Y = asin(-2 * q1 * q3 + 2 * q0* q2)* 57.3; / pitchQ_ANGLE.X = atan2(2 * q2 * q3 + 2 * q0 * q1,
8、-2 * q1 * q1 - 2 * q2* q2 + 1)* 57.3; / roll逐条解释。姿态矩阵可以由以下两种方式表示第一个就是上图所说的欧拉角法（式（1） ） ，还有一个就是四元数法021232(q103)2(q1302) CbR(q0)q 注意！1 这里是 CbR,假设 b 为四旋翼固连坐标系， R 为参考坐标系，那么CbR 表示 b 系到 R 系的坐标变换矩阵，由于 (1)式表示的为 R 系到 b 系的坐标矩阵，要用上式表示，则要对四元数法矩阵求逆，又因为该矩阵为正交阵，逆等于转置，则描述 R 系到 b 系的四元数矩阵为 021232(q103)2(q1302)C(q0)q 此
9、时矩阵跟 1 式矩阵一一对应。1.float q0q0 = q0*q0;float q0q1 = q0*q1;float q0q2 = q0*q2;/ float q0q3 = q0*q3;float q1q1 = q1*q1;/ float q1q2 = q1*q2;float q1q3 = q1*q3;float q2q2 = q2*q2; float q2q3 = q2*q3;float q3q3 = q3*q3;这段程序就是为了把需要用到的姿态矩阵的元素求出来给出的。2.vx = 2*(q1q3 - q0q2); /vy = 2*(q0q1 + q2q3);vz = q0q0 - q1
10、q1 - q2q2 + q3q3 ;可以看到 vx，vy，vz 为 CRb 的最后一列的三项，四元数矩阵带入（1）式得 vx，vy，vz分别是 axB，ayB，azB 每一项 g 前的系数。且静止情况下 vx，vy，vz 组成向量模长基本可以认为为 1.3.norm = sqrt(ax*ax + ay*ay + az*az); /acc 数据归一化ax = ax /norm;ay = ay / norm;az = az / norm;以上已说，由四元数倒推回去的加速度，向量模长为 1，为了比较误差进行归 1 化，算的由加计得出的向量。4.ex = (ay*vz - az*vy) ; ey =
11、(az*vx - ax*vz) ;ez = (ax*vy - ay*vx) ;接着可以通过叉乘（向量外积）计算误差5.exInt = exInt + ex * Ki;eyInt = eyInt + ey * Ki;ezInt = ezInt + ez * Ki;对误差进行积分6.gx = gx + Kp*ex + exInt; gy = gy + Kp*ey + eyInt;gz = gz + Kp*ez + ezInt;进行 pi 滤波7q0 = q0 + (-q1*gx - q2*gy - q3*gz)*halfT;q1 = q1 + (q0*gx + q2*gz - q3*gy)*hal
12、fT;q2 = q2 + (q0*gy - q1*gz + q3*gx)*halfT;q3 = q3 + (q0*gz + q1*gy - q2*gx)*halfT;对四元数进行跟新，这里用的是方程的数值解法，求得的解释近似解，总之就是跟新了四元数8norm = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);q0 = q0 / norm;q1 = q1 / norm; q2 = q2 / norm;q3 = q3 / norm;对四元数进行规范化，即化为模长为 1 ，因为只有规范化的四元数才能表示刚体旋转。9Q_ANGLE.Y = asin(-2 * q1 * q3
13、+ 2 * q0* q2)* 57.3; / pitchQ_ANGLE.X = atan2(2 * q2 * q3 + 2 * q0 * q1, -2 * q1 * q1 - 2 * q2* q2 + 1)* 57.3; / roll仍旧一一对应关系发现 2（q1q3 -q0q2）刚好跟欧拉角法的 对应，由此利用自带库sin函数即可求得俯仰角，横滚角类似，偏航角由于没有罗盘进行校正求没有意义，控制中采用采用 PD 控制。补充，由于陀螺仪会有零点漂移开始一定要进行补偿。这段是在 mpu6050.c 中程序，对直流偏执进行了补偿。MPU6050_ACC_LAST.X=(int16_t)mpu605
14、0_buffer0) 8) | mpu6050_buffer1) - ACC_OFFSET.X;MPU6050_ACC_LAST.Y=(int16_t)mpu6050_buffer2) 8) | mpu6050_buffer3) - ACC_OFFSET.Y;MPU6050_ACC_LAST.Z=(int16_t)mpu6050_buffer4) 8) | mpu6050_buffer5) - ACC_OFFSET.Z;MPU6050_GYRO_LAST.X=(int16_t)mpu6050_buffer8) 8) | mpu6050_buffer9) - GYRO_OFFSET.X;MPU6
15、050_GYRO_LAST.Y=(int16_t)mpu6050_buffer10) 8) | mpu6050_buffer11) - GYRO_OFFSET.Y;MPU6050_GYRO_LAST.Z=(int16_t)mpu6050_buffer12) 8) | mpu6050_buffer13) - GYRO_OFFSET.Z;这里还要说一点，这里加速计的数据用的是滑动平均值滤波法void Prepare_Data(void)static uint8_t filter_cnt=0;int32_t temp1=0,temp2=0,temp3=0;uint8_t i;MPU6050_Read();MPU6050_Dataanl();ACC_X_BUFfilter_cnt = MPU6050_ACC_LAST.X;ACC_Y_BUFfilter_cnt = MPU6050_ACC_LAST.Y;ACC_Z_BUFfilter_cnt = MPU6050_ACC_LAST.Z;f
