1、新手四旋翼算法总结一.姿态结算(匿名版程序)首先，程序中一般用了两种求解姿态的方法，一种为欧拉角法，一种为四元数法(1)欧拉角法静止状态，或者总加速度只是稍微大于g时，由加计算出的值比较准确。使用欧拉角表示姿态，令，。和代表ZYX欧拉角，分别称为偏航角、俯仰角和横 滚角。载体坐标系下的力口速 度SB,ayB,azB)和参考坐标系下的加速度(axN, ayN, a,N)之 间的关系可表示为(1)。其中C和S分别代表cos和sino axB,ayB,azB就是mpu读出来的 三个值。这个矩阵就是三个旋转矩阵相乘得到的，因为矩阵的乘法可以表示旋转。axBcccs-saxNayB=-cs + ssc
2、cc + sssscayNazBss + csc -sc + cssccazN(1)飞行器处于静止状态，此时参考系下的加速度等于重力加速度，即把(2)代入(1)可以解的 = arctg (c-=) (3)ayB + azB - arctg(4)即为初始俯仰角和横滚角，通过加速度计得到载体坐标系下的加速度即可将其解出，偏 航角可以通过电子罗盘求出。(2)四元数法(通过处理单位采样时间内的角增量(mpu的陀螺仪得到的就是角增量)，为 了避免噪声的微分放大，应该直接用角增量抄的书)上匿名的程序void IMUupdate(float gx, float gy, float gz, float ax,
3、 float ay, float az)(float norm;/ float hx, hy, hz, bx, bz;float vx, vy, vz;/ wx, wy, wz;float ex, ey, ez;/先把这些用得到的值算好float qq = q*q;float qql = q*ql;float qq2 = q*q2;/ float qq3 = q*q3;float qlql = ql*ql;/ float qlq2 = ql*q2;float qlq3 = ql*q3;float q2q2 = q2*q2;float q2q3 = q2*q3;float q3q3 = q3*q
4、3;if(ax*ay*az=O)return;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*(qlq3 - qq2);四元素中 xyz 的vy = 2*(q0ql + q2q3);vz = qq - qlql - q2q2 + q3q3 ;/ error is sum of cross product between reference direct
5、ion of fields and direction measured by sensorsex = (ay*vz - az*vy);向量外积在相 减得到差分就是误差ey = (az*vx - ax*vz); ez =(ax*vy - ay*vx);对误差进行积分/将误差PI/这里的gz由于exlnt = exlnt + ex * Ki; eylnt = eylnt + ey * Ki; ezlnt = ezlnt + ez * Ki;/ adjusted gyroscope measurements gx = gx + Kp* ex + exlnt;后补偿到陀螺仪，即补偿零点漂移gy =
6、gy + Kp*ey + eylnt; gz = gz + Kp*ez + ezlnt;没有观测者进行矫正会产生漂移，表现出来的就是积分自增或自减/四元素的微分方程/ integrate quaternion rate and normalise q = q + (-ql*gx - q2*gy - q3*gz)*half; ql = ql + (q*gx + q2*gz - q3*gy)*half; q2 = q2 + (q*gy - ql*gz + q3*gx)*half; q3 = q3 + (q*gz + ql*gy - q2*gx)*half;/ normalise quaternio
7、nnorm = sqrt(q*q + ql*ql + q2*q2 + q3*q3);q = q / norm;ql = ql / norm;q2 = q2 / norm;q3 = q3 / norm;/Q_ANGLE.Yaw = atan2(2 * ql * q2 + 2 * q * q3, -2 * q2*q2 - 2 * q3* q3 + 1)* 57.3; / yawQ_ANGLE.Y = asin(-2 * ql * q3 + 2 * q* q2)* 57.3; / pitchQ_ANGLE.X = atan2(2 * q2 * q3 + 2 * q * ql, -2 * ql * q
8、l - 2 * q2* q2 + 1)* 57.3;/ roll)逐条解释。姿态矩阵可以由以下两种方式表示第一个就是上图所说的欧拉角法(式(1),还有一个就是四元数法%0八2 + 少2-夕2八2-夕3八2 2(qlql-q0q3)2(qlq3+q0q2)CbR=2(qlq2 + q0q3)1 2(ql 八 2 + q3 八 2)2(q2q3-q0ql)2(qlq3-q0q2)2(q2q3 + q0ql) 02-12-22 + 32注意！ 1这里是CbR,假设b为四旋翼固连坐标系，R为参考坐标系，那么CbR表示b系到R系的坐标变换矩阵，由于式表示的为R系到b系的坐标矩阵，要用上式表示， 则要对四
9、元数法矩阵求逆，又因为该矩阵为正交阵，逆等于转置，则描述R系到b系的四 元数矩阵为%0八2 + 少2-夕2八2-叱22(qlql-q0q3)2(qlq3-q0q2)CRb=2(qlq2 + q0q3)1-2(ql 八 2 + q3 八 2)2(q2q3 + q0ql)2(qlq3 + q0q2)2(q2q3-q0ql)如八 2-必八 2-q2八 2 + 夕3 八 2此时矩阵跟1式矩阵一一对应。1.float qq = q*q;float qql= q*ql;float qq2 = q*q2;/ float qq3 = q*q3;float qlql = ql*ql;/ float qlq2 =
10、 ql*q2;float qlq3 = ql*q3;float q2q2 = q2*q2;float q2q3 = q2*q3;float q3q3 = q3*q3;这段程序就是为了把需要用到的姿态矩阵的元素求出来给出的。2.Vx = 2*(qlq3 - qq2);/Vy = 2*(q0ql + q2q3);vz = qq - qlql - q2q2 + q3q3 ;可以看到vx, vy, VZ为CRb的最后一列的三项，四元数矩阵带入(1)式得vx, vy, VZ分 别是axB, ayB, azB每一项g前的系数。且静止情况下vx, vy, VZ组成向量模长基本可以 认为为1.3.norm =
11、 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 = (az*vx - ax*vz);ez = (ax*vy - ay*vx);接着可以通过叉乘(向量外积)计算误差5.exlnt = exlnt + ex * Ki;eylnt = eylnt + ey * Ki;ezlnt = ezlnt + ez * Ki;对误差进行积分6.g
12、x = gx + Kp* ex + exlnt;gy = gy + Kp*ey + eylnt;gz = gz + Kp*ez + ezlnt;进行Pi滤波7.q = q + (-ql*gx - q2*gy - q3*gz)*half;ql = ql + (q*gx + q2*gz - q3*gy)*half;q2 = q2 + (q*gy - ql*gz + q3*gx)*half;q3 = q3 + (q*gz + ql*gy - q2*gx)*half;对四元数进行跟新，这里用的是方程的数值解法，求得的解释近似解，总之就是跟新了四 元数8.norm = sqrt(q*q + ql*ql
13、 q2*q2 + q3*q3);q = q / norm;ql = ql / norm;q2 = q2 / norm;q3 = q3 / norm;对四元数进行规范化，即化为模长为1,因为只有规范化的四元数才能表示刚体旋转。9.Q_ANGLE.Y = asin(-2 * ql * q3 + 2 * q* q2)* 57.3; / pitchQ_ANGLE.X = atan2(2 * q2 * q3 + 2 * q * ql, -2 * ql * ql - 2 * q2* q2 + 1)* 57.3;/ roll仍旧一一对应关系发现2 (qlq3 -qq2)刚好跟欧拉角法的Sine对应，由此利
14、用自带库函 数即可求得俯仰角，横滚角类似，偏航角由于没有罗盘进行校正求没有意义，控制中采用采 用PD控制。补充，由于陀螺仪会有零点漂移开始一定要进行补偿。这段是在mpu6050.c中程序，对直 流偏执进行了补偿。MPU6050_ACC_LAST.X=(intl6_t)mpu6050_buffer0) 8) mpu6050.bufferl)- ACC_OFFSET.X;MPU6050_ACC_LAST.Y=(intl6_t)mpu6050_buffer2) 8) mpu6050.buffer3)- ACC_OFFSET.Y;MPU6050_ACC_LAST.Z=(intl6_t)mpu6050_
15、buffer4) 8) mpu6050.buffer5)- ACC_OFFSET.Z;MPU6050_GYRO_LAST.X=(intl6_t)mpu6050_buffer8) 8) mpu6050.buffer9)- GYRO_OFFSET.X;MPU6050_GYRO_LAST.Y=(intl6_t)mpu6050_buffer10) 8) mpu6O5O.bufferll) -GYRO_OFFSET.Y;MPU6050_GYRO_LAST.Z=(intl6_t)mpu6050_buffer12) 8) mpu6O5O.buffer13)- GYRO_OFFSET.Z;这里还要说一点，这里
16、加速计的数据用的是滑动平均值滤波法void Prepare_Data(void)(static uint8_t filter_cnt=O;int32_t temp 1=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;for(i=0;iFIETER_NUM;i+)tempi +=
17、ACC_X_BUFi;temp2 += ACC_Y_BUFi;temp3 += ACC_Z_BUFi;ACC_AVG.X = tempi / FILTER_NUM;ACC_AVG.Y = temp2 / FILTER_NUM;ACC_AVG.Z = temp3 / FILTER_NUM;filter_cnt+;if(filter_cnt=FILTER_NUM) filter_cnt=O;GYROJ.X += MPU6050_GYRO_LAST.X*Gyro_G*0.0001;GYROJ.Y += MPU6050_GYRO_LAST. Y*Gyro_G*0.0001;GYROJ.Z += MPU6050_GYRO_LAST.Z*Gyro_G*0.0001;)资料在附带文件中。
免费区 
