四元数转化成欧拉角.docVIP

四元数转化成欧拉角 笔者奋斗 修改时间：2015.11.23 姿态解算（以匿名版程序为例） 首先，程序中一般用了两种求解姿态的方法，一种为欧拉角法，一种为四元数法。由四元数法将MPU6050的值融合成四元再由欧拉角法把四元转换成欧拉角送入串级PID控制器 （1）欧拉角法静止状态，或者总加速度只是稍微大于g时，由加计算出的值比较准确。使用欧拉角表示姿态，令Φ,θ和Φ代表ZYX 欧拉角，分别称为偏航角、俯仰角和横滚角。载体坐标系下的 加 速 度(axB,ayB,azB)和参考坐标系下的加速度(axN, ayN, azN)之间的关系可表示为(1)。其中 c 和 s 分别代表 cos 和 sin。axB,ayB,azB就是mpu读出来的三个值。 这个矩阵就是三个旋转矩阵相乘得到的，因为矩阵的乘法可以表示旋转。 (1) 飞行器处于静止状态，此时参考系下的加速度等于重力加速度，即（2） 把（2）代入（1）可以解得 （3） （4） 即为初始俯仰角和横滚角，通过加速度计得到载体坐标系下的加速度即可将其解出，偏航角可以通过电子罗盘求出。 （2）四元数法（通过处理单位采样时间内的角增量(mpu的陀螺仪得到的就是角增量)，为了避免噪声的微分放大，应该直接用角增量） void IMUupdate(float gx, float gy, float gz, float ax, float ay, float az) { float norm; // 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) 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*(q1q3 - q0q2); //四元素中xyz vy = 2*(q0q1 + q2q3); vz = q0q0 - q1q1 - q2q2 + q3q3 ; // error is sum of cross product between reference direction of fields and direction measured by sensors ex = (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 measurements gx = gx + Kp*ex + exInt; //将误差PI后补偿到陀螺仪，即补偿零点漂移 gy = gy + Kp*ey + 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 + (q