开源制作小四旋翼全过程

肖久

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 - q1q1 - 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 = (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滤波，其实就是互补滤波

7．

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;

龙格库塔法。。。就是方程的数值解法。。近似解。。一阶解法

0736fac7228d4220f912874ee8cee5e5_21.png (0 Bytes, 下载次数: 0)

下载附件

2010-12-14 22:54 上传


这个跟四元数的微分方程对应有兴趣的看看书。。。。

8．

norm = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);

q0 = q0 / norm;

q1 = q1 / norm;

q2 = q2 / norm;

q3 = q3 / norm;

对四元数进行规范化，即化为模长为1 ，因为只有规范化的四元数才能表示刚体旋转。

9．

Q_ANGLE.Y = asin(-2 * q1 * q3 + 2 * q0* q2)* 57.3; // pitch

Q_ANGLE.X = atan2(2 * q2 * q3 + 2 * q0 * q1, -2 * q1 * q1 - 2 * q2* q2 + 1)* 57.3; // roll

仍旧一一对应关系发现2（q1q3 -q0q2）刚好跟欧拉角对应，由此利用自带库函数即可求得俯仰角，横滚角类似，偏航角由于没有罗盘进行校正求没有意义，控制中采用采用PD控制。

补充，由于陀螺仪会有零点漂移开始一定要进行补偿。这段是在mpu6050.c中程序，对直流偏执进行了补偿。

MPU6050_ACC_LAST.X=((((int16_t)mpu6050_buffer[0]) << 8) | mpu6050_buffer[1]) - ACC_OFFSET.X;

MPU6050_ACC_LAST.Y=((((int16_t)mpu6050_buffer[2]) << 8) | mpu6050_buffer[3]) - ACC_OFFSET.Y;

MPU6050_ACC_LAST.Z=((((int16_t)mpu6050_buffer[4]) << 8) | mpu6050_buffer[5]) - ACC_OFFSET.Z;

MPU6050_GYRO_LAST.X=((((int16_t)mpu6050_buffer[8]) << 8) | mpu6050_buffer[9]) - GYRO_OFFSET.X;

MPU6050_GYRO_LAST.Y=((((int16_t)mpu6050_buffer[10]) << 8) | mpu6050_buffer[11]) - GYRO_OFFSET.Y;

MPU6050_GYRO_LAST.Z=((((int16_t)mpu6050_buffer[12]) << 8) | mpu6050_buffer[13]) - GYRO_OFFSET.Z;

这里还要说一点，这里加速计的数据用的是滑动平均值滤波法，一定要有这个。。不然由于机械振动造成的影响非常大。。。

void repare_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_BUF[filter_cnt] = MPU6050_ACC_LAST.X;

ACC_Y_BUF[filter_cnt] = MPU6050_ACC_LAST.Y;

ACC_Z_BUF[filter_cnt] = MPU6050_ACC_LAST.Z;

for(i=0;i<FILTER_NUM;i++)

{

temp1 += ACC_X_BUF;

temp2 += ACC_Y_BUF;

temp3 += ACC_Z_BUF;

}

选择716电机，720的转速跟716的差不多，注意电机孔直径一定要大一点。。不然塞不进去。。然后补充一下MPU6050的摆放位置没有关系。。同一坐标系下测的的加速度角速度都是没有关系的。。。

关于电源还有QFN的芯片。。。注意这个其实很好焊，用烙铁沾点锡加点焊锡膏在对准的引脚上拖焊就行。。可以拿个灯照着看反光比较容易对准。

关于软件最关建的说说。。只有姿态解算部分。。PID部分我的算法还得改。。。。这个网上有开源的就是串机PID。。。额。。本人菜鸟。。还没看懂。。大二还没学自控。。回去会看的。。。注意这里MPU最好用硬件IIc，因为小四轴的姿态更新频率是1000HZ比较快，这里的IIC只是一个器件。。目前还没出什么问题。
这里有一部分是我之前写的总结

（1）欧拉角法静止状态，或者总加速度只是稍微大于g时，由加计算出的值比较准确。

使用欧拉角表示姿态，令Φ,θ和Φ代表ZYX 欧拉角，分别称为偏航角、俯仰角和横滚角 。 载体坐标系下的 加 速 度(axB,ayB,azB)和参考坐标系下的加速度(axN, ayN, azN)之间的关系可表示为(1)。其中 c 和 s 分别代表 cos 和 sin。axB,ayB,azB就是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 + (q0*gz + q1*gy - q2*gx)*halfT;

// normalise quaternion

norm = 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; // yaw

Q_ANGLE.Y = asin(-2 * q1 * q3 + 2 * q0* q2)* 57.3; // pitch

Q_ANGLE.X = atan2(2 * q2 * q3 + 2 * q0 * q1, -2 * q1 * q1 - 2 * q2* q2 + 1)* 57.3; // roll

}

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 - q1q1 - 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 = (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滤波，其实就是互补滤波

7．

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;

龙格库塔法。。。就是方程的数值解法。。近似解。。一阶解法


这个跟四元数的微分方程对应有兴趣的看看书。。。。

8．

norm = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);

q0 = q0 / norm;

q1 = q1 / norm;

q2 = q2 / norm;

q3 = q3 / norm;

对四元数进行规范化，即化为模长为1 ，因为只有规范化的四元数才能表示刚体旋转。

9．

Q_ANGLE.Y = asin(-2 * q1 * q3 + 2 * q0* q2)* 57.3; // pitch

Q_ANGLE.X = atan2(2 * q2 * q3 + 2 * q0 * q1, -2 * q1 * q1 - 2 * q2* q2 + 1)* 57.3; // roll

仍旧一一对应关系发现2（q1q3 -q0q2）刚好跟欧拉角对应，由此利用自带库函数即可求得俯仰角，横滚角类似，偏航角由于没有罗盘进行校正求没有意义，控制中采用采用PD控制。

补充，由于陀螺仪会有零点漂移开始一定要进行补偿。这段是在mpu6050.c中程序，对直流偏执进行了补偿。

MPU6050_ACC_LAST.X=((((int16_t)mpu6050_buffer[0]) << 8) | mpu6050_buffer[1]) - ACC_OFFSET.X;

MPU6050_ACC_LAST.Y=((((int16_t)mpu6050_buffer[2]) << 8) | mpu6050_buffer[3]) - ACC_OFFSET.Y;

MPU6050_ACC_LAST.Z=((((int16_t)mpu6050_buffer[4]) << 8) | mpu6050_buffer[5]) - ACC_OFFSET.Z;

MPU6050_GYRO_LAST.X=((((int16_t)mpu6050_buffer[8]) << 8) | mpu6050_buffer[9]) - GYRO_OFFSET.X;

MPU6050_GYRO_LAST.Y=((((int16_t)mpu6050_buffer[10]) << 8) | mpu6050_buffer[11]) - GYRO_OFFSET.Y;

MPU6050_GYRO_LAST.Z=((((int16_t)mpu6050_buffer[12]) << 8) | mpu6050_buffer[13]) - GYRO_OFFSET.Z;

这里还要说一点，这里加速计的数据用的是滑动平均值滤波法，一定要有这个。。不然由于机械振动造成的影响非常大。。。

void repare_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_BUF[filter_cnt] = MPU6050_ACC_LAST.X;

ACC_Y_BUF[filter_cnt] = MPU6050_ACC_LAST.Y;

ACC_Z_BUF[filter_cnt] = MPU6050_ACC_LAST.Z;

for(i=0;i<FILTER_NUM;i++)

{

temp1 += ACC_X_BUF;

temp2 += ACC_Y_BUF;

temp3 += ACC_Z_BUF;

}


电路图文件：

四轴飞行第三版-改.PcbDoc (6.96 MB, 下载次数: 24)

2016-7-21 10:34 上传

点击文件名下载附件
描述：四轴第三版
教程
基于PI_PD控制器的四旋翼姿态控制_唐健杰.pdf (232.97 KB, 下载次数: 26)

2016-7-21 10:34 上传

点击文件名下载附件
描述：论文
源代码
四轴第三版-改.SchDoc (317.5 KB, 下载次数: 18)

2016-7-21 10:34 上传

点击文件名下载附件
描述：四轴飞行第三版
其他文件
小四旋翼.rar.part (1.95 MB, 下载次数: 25)

2016-7-21 10:34 上传

点击文件名下载附件
描述：参考资料等

我知道答案 目前已有8人回答

阅读全文

肖久

stary666 发表于 2016-7-21 11:24
小九，验证过吗？？

亲测啊

modou0

要开搞麽？

kingsings

墙裂收藏

zhugenping126

学习一下啦

华尔兹兹

6666666666666

华尔兹兹

虽然看不懂  不明觉厉 呱唧呱唧

chen00000009

谢谢分享，代码略乱啊

Y.fei

看着好高深