5iMX宗旨:分享遥控模型兴趣爱好

5iMX.com 我爱模型 玩家论坛 ——专业遥控模型和无人机玩家论坛(玩模型就上我爱模型,创始于2003年)
查看: 7579|回复: 27
打印 上一主题 下一主题

随手科普无人机核心技术(含代码)

  [复制链接]
跳转到指定楼层
楼主
发表于 2018-1-25 15:52 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
点击查看详情
本帖最后由 空中飞鹅 于 2018-2-2 19:04 编辑

2018年1月30日修改。
把最开始放出来的代码和伪码去掉了。
直接给出一个GPS飞控的工作逻辑框架。
为方便网站审核,尽量使用中文而避免使用非中文。因为非中文一多,帖子发了一天也未必审核通过。全是中文的话,几秒钟就审核完了。
——————————————————————————————————————————————————————————
   一、通讯逻辑

为简单起见,暂时不考虑地面站(PC端或移动端),而只采用遥控器控制。地面站功能跟遥控器工功能有很多重叠的地方,而不重叠的地方对飞控主体程序影响不大。为了简单起见,很多地方的细节都暂时忽略,而主要针对躯干部分进行描述。


遥控器有8个通道,其中摇杆占四个通道,多段开关占两个通道,另外有两个按键通道。


摇杆通道定义为美国手模式,即左手油门。注:按照“左西右东”的原则,日本为东,美国为西,很容易记得美国手和日本手的区别。

通道1-通道4分别为:横滚、俯仰、油门和航向。左手摇杆左右方向为航向通道,打到最左端对应1000us,落在中间为1500us,打到最右端对应2000us。左手摇杆上下方向为油门通道,打到最下段为1000us,中间为1500us,最上端为2000us。右手摇杆左右方向为横滚通道,打到最左端为1000us,中间为1500us,最右端为2000us。右手摇杆上下方向为俯仰通道,打到最上端为1000us,中间为1500us,最下段为2000us。注:俯仰通道和油门通道不一样,位时映射关系正好相反。

有头状态下,并且在摇杆控制有效的情况下:航向通道打到左边,机头逆时针转动,航向通道打到右边,机头顺时针转动;油门通道打到上边,飞机往上走,油门通道打到下边,飞机往下走;俯仰通道打到上边,飞机往前走,俯仰通道打到下边,飞机往后走;横滚通道打到左边,飞机往左走,横滚通道打到右边,飞机往右走。


遥控器左上角三段开关定义为通道5,右上角三段开关定义为通道6

假定飞控支持手动、增稳、定高、自动起飞/自动降落,无头/有头、悬停、返航、绕点和跟随。通道5和通道6搭配起来实现:手动、增稳、定高、悬停、绕点、跟随和返航。通道7实现自动起飞/自动降落,通道8实现无头/有头。

通道5在上段,通道6在上段,为手动模式;通道5在上段,通道6在中段,为增稳模式,通道5在上段,通道6在下段,为定高模式。通道5在中段,通道6在上段,为悬停模式,通道5在中段,通道6在中段,为绕点模式,通道5在下段,通道6在下段,为跟随模式。通道5在下段,为返航模式。


两个按键开关定义为通道7-8

飞机尚未起飞的时候,长按遥控器中间偏左的通道7按键,飞机为自动起飞状态,飞机已经在空中并非是自动起飞状态,这时候按下通道7按键,飞机为自动降落状态。飞机上电后,飞机默认为有头状态,上电后长按遥控器中间偏右的通道8按键,飞机为无头状态,再长按返回有头状态。


二、控制逻辑


一共有12个控制环。其中位置3个,速度3个,角度3个,角速度3个。

GPS定位之后,位置环和速度环参与运作。GPS丢失或未定位,位置环和速度环失效。

假定飞机的姿态和位置都已经算好了,控制部分代码可以分为四个大函数:水平控制函数、高度控制函数、航向控制函数和电机控制函数。


1、水平控制函数又分速度计算、角度计算、舵机输出计算三大块。

速度计算部分,分为两段:有GSP和无GPS

1)GPS时,速度环失效,速度环PID参数清零。

2)GPS时,又分速度打杆、速度刹车和位置锁定三种情况。

设横滚俯仰死区都是50。为了严格区分三种情况,需要设定一个位置锁定标志。

情况一(速度打杆):横滚杆量(或俯仰杆量)大于1450,或小于1550。位置锁定标志为0

此时杆量直接对应机体坐标系的目标速度。当杆量小于1450时,目标速度为10米每秒乘以(杆量-1450)再除以(1450-1000)。当杆量大于1550时,目标速度为10米每秒乘以(杆量-1550)再除以(2000-1550)。

机体坐标系目标速度经矢量分解,变到NED坐标系,得到东向目标速度和北向目标速度。

情况二(速度刹车):横滚杆量(或俯仰杆量)在14501550之间。位置锁定标志为零。

此时直接令目标速度为0。飞机进入刹车状态。

当飞机实际速度足够小(在正负0.1米之间),把位置锁定标志设为1,把当前位置设为目标位置。

情况三(位置锁定):横滚杆量(或俯仰杆量)在14501550之间。位置锁定标志为1

A)当飞机处于悬停模式,由目标位置和当前位置做PID计算,得到NED坐标系的目标速度。

B)当飞机处于返航模式。切入返航瞬间,把目标位置设为当前位置。先进入返航三状态中的原地爬高状态。

先判断飞机当前位置与回家点之间的距离。如果距离超过20米,把目标高度设为20米,等飞机高度到达20米后才进入原高返航状态。如果距小于20米,直接进入原高返航状态。

在原高爬高状态,由目标位置和当前位置做PID计算,得到北向目标速度和东向目标速度。由目标高度得到垂直目标速度。

在原高返航状态,目标位置改为回家点。同样由位置环得到NED目标速度。

当飞机与回家点之间的距离小于1米时,进入自动降落状态。

在自动降落状态,目标位置还是回家点,目标高度为零。由位置环得到NED目标速度。由高度环得到垂直目标速度。

C) 当飞机处于绕点模式时,刚切入时目标位置就是当前位置。当右边摇杆打杆时,临时执行打杆控制操作,松杆刹车后飞机停住。这个位置到之前的目标位置之间的距离,如果小于3米,就认为绕点无效,飞机原地悬停。如果距离大于3米,则飞机之后的目标位置,就是以刚切入时的目标位置为圆心,以距离为半径的圆。利用圆的方程,以固定的角速度改变从圆心出发到飞机的方位角,就能让飞机沿着圆形轨迹飞行。飞行过程中,目标航向角就是从飞机出发到圆心的方位角。依然由目标位置和当前位置得到NED目标速度。

D) 当飞机处于跟随模式时,刚切入时计算一个飞机与人之间的距离(人的位置由遥控器上的GPS给出),然后,当人的位置发生改变时,把人的位置改变量加到飞机的位置坐标上去,目标位置就出来了。目标位置有了,当前位置有了,NED目标速度也就有了。航向角控制上,让目标航向角等于从飞机到人的方位角。


角度计算部分。


还是分为有GPS和无GPS两种情况。


GPS的时候,目标速度与当前速度做PID计算,可以得到北向的目标倾角和东向的目标倾角。经过矢量分解,得到机体坐标系的目标倾角,即目标横滚角和目标俯仰角。


GPS情况下。

如果是手动模式,目标角速度由杆量给出。映射关系与杆量到速度的非常类似。

如果是手动模式,目标角度由杆量给出。映射关系与杆量到速度的非常类似。


舵机输出计算部分。

经过角度角速度串级PID(或单角速度环),可以得到横滚俯仰两个通道的舵机输出。


2、高度控制

  摇杆设死区。死区操作如前。

  如果飞机没有起飞,自动起飞/自动降落按键有值,就给一个目标高度,经过高度环和垂直速度环的串级PID计算,得到油门通道的舵机输出。

  如果飞机已经在空中,自动起飞/自动降落按键有值,同样操作。

  如果飞机已经在空中,并且工作在手动模式、增稳模式,对高度环、垂直速度环PID参数清零,其他情况,杆量对应于高度增量。经高度环和垂直速度环,可以得到舵机输出。


3、航向控制


摇杆设死区。死区操作如前。

杆量对应角度增量。经角度环和角速度环,可得航向舵机输出。


4、电机控制函数

把电机状态分为关闭和打开两种状态。

当飞机没有起飞,并且自动起飞按钮未按下,飞机掰杆解锁(外八或内八)。

飞机解锁后,由四个通道的舵机输出(通道控制量)进行动力分配,得到四个电机对应的控制量。

这里给出动力分配矩阵的数学推导过程。


设机身右上角的电机为电机1(逆时针转),右下角的电机为电机2(顺时针转),左下角的电机为电机3(逆时针转),左上角的电机为电机4(顺时针转)。

当俯仰角增加时,电机1和电机4转速加大,电机2和电机3转速减小。

当横滚角增加时,电机1和电机2转速减小,电机3和电机4转速增大。

当高度增加时,四个电机同时增加。

当航向角增加时,电机1和电机3转速增大,电机2和电机4转速减小。

得:俯仰控制量=电机1控制量-电机2控制量-电机3控制量+电机4控制量

    横滚控制量=-电机1控制量-电机2控制量+电机3控制量+电机4控制量

    油门控制量=电机1控制量+电机2控制量+电机3控制量+电机4控制量

    航向控制量=电机1控制量-电机2控制量+电机3控制量-电机4控制量

写成矩阵形式,求逆矩阵,可得:

电机1控制量=油门控制量-横滚控制量+俯仰控制量+航向控制量

电机2控制量=油门控制量-横滚控制量-俯仰控制量-航向控制量

电机3控制量=油门控制量+横滚控制量-俯仰控制量+航向控制量

电机4控制量=油门控制量+横滚控制量+俯仰控制量-航向控制量


测量部分后面展开。此处从略。

至此,一个GPS飞控的简单框架完全确立。






评分

参与人数 1威望 +2 收起 理由
怪怪 + 2 很给力!

查看全部评分

欢迎继续阅读楼主其他信息

主题

  • 没有相关信息
  • 没有相关信息
  • 没有相关信息
沙发
发表于 2018-1-25 17:36 | 只看该作者
牛逼,开源顶级代码

3
发表于 2018-1-25 19:13 | 只看该作者
本帖最后由 xxdcq 于 2018-1-25 23:40 编辑

姿态计算里的各加速度需要先归一化处理,而且缺少积分运算,只有比例运算是不稳定的我给你贴一个完全正确的

#define pi 3.14159265f    //圆周率                          
#define Kp 0.8f           //比例系数                        
#define Ki 0.001f         //积分系数                        
#define halfT 0.004f      //8毫秒运算周期,所以半周期就是4毫秒=0.004秒           
float idata q0=1,q1=0,q2=0,q3=0;  //四元数初始化赋值   
float idata exInt=0,eyInt=0,ezInt=0;
float idata q0q0 = 1;
float idata q0q1 = 0;
float idata q0q2 = 0;
float idata q0q3 = 0;
float idata q1q1 = 0;
float idata q1q2 = 0;
float idata q1q3 = 0;
float idata q2q2 = 0;
float idata q2q3 = 0;
float idata q3q3 = 0;
float idata DCMgb02=0;
float idata DCMgb12=0;
float idata DCMgb22=1;

//函数名:invSqrt(void)
//描述:求平方根的倒数
//该函数是经典的Carmack求平方根算法,效率极高,使用魔数0x5f375a86
static float invSqrt(float number)
{
   volatile long i;
   volatile float x, y;
   volatile const float f = 1.5F;
   x = number * 0.5F;
   y = number;
   i = * (( long * ) &y);
   i = 0x5f375a86 - ( i >> 1 );
   y = * (( float * ) &i);
   y = y * ( f - ( x * y * y ) );
   return y;
}

void IMUupdate(float gx, float gy, float gz, float ax, float ay, float az)
{
  float idata recipNorm;
  float idata vx, vy, vz;
  float idata ex, ey, ez;
  recipNorm = invSqrt(ax*ax + ay*ay + az*az);  
  ax *= recipNorm;     //归一化处理
  ay *= recipNorm;     //归一化处理
  az *= recipNorm;     //归一化处理         
  vx = DCMgb02;                                                                                                
  vy = DCMgb12;
  vz = DCMgb22;
  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;  //误差积分
  gx = gx + Kp*ex + exInt;                                                                                                   
  gy = gy + Kp*ey + eyInt;
  gz = gz + Kp*ez + ezInt;
//得到新的四元数                                                                                                                                      
  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;
//归一化处理
  recipNorm = invSqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);
  q0 *= recipNorm;
  q1 *= recipNorm;
  q2 *= recipNorm;
  q3 *= recipNorm;
  q0q0 = q0*q0;
  q0q1 = q0*q1;
  q0q2 = q0*q2;
  q0q3 = q0*q3;
  q1q1 = q1*q1;
  q1q2 = q1*q2;
  q1q3 = q1*q3;
  q2q2 = q2*q2;
  q2q3 = q2*q3;
  q3q3 = q3*q3;
//得到地理坐标系转机体坐标系的旋转矩阵
//机体坐标系转地理坐标系只需要将行列交换就行了
  DCMgb02 = 2*(q1q3 - q0q2);                        // 13
  DCMgb12 = 2*(q2q3 + q0q1);                        // 23
  DCMgb22 = q0q0 - q1q1 - q2q2 + q3q3;                // 33

  Angle=asin(2*(q0q2-q1q3 ))* 57.2957795f; // 俯仰角
  Angley=asin(2*(q0q1+q2q3 ))* 57.2957795f; // 横滚角
}




4
发表于 2018-1-25 21:11 | 只看该作者
好厉害虽然大部分看不懂
来自安卓客户端来自安卓客户端
5
发表于 2018-1-25 22:13 | 只看该作者
厉害,如果结合图形案例讲解就更具体了
6
 楼主| 发表于 2018-1-26 02:03 | 只看该作者

两个坐标系与三个姿态角

本帖最后由 空中飞鹅 于 2018-1-26 02:07 编辑

       绝大部分的飞控系统(含开源和商用)都会使用北东地坐标系作为静态参考。动态的自然就是机体坐标系了。把飞机放平,机头对准正北。此时机头方向就是X轴,机体右侧就是Y轴,垂直机身向下就是Z轴。通常认为北东地坐标系是不会动的,而机体坐标系时刻在动,这时候,两个坐标系之间就会产生一些夹角。最常用的夹角就是三个姿态角:横滚角、俯仰角和航向角。

        严格意义上讲,姿态角只有横滚和俯仰,航向角和姿态角并列看待。但对于一般人来说,没必要分得太细,直接认为三个都是姿态角就好了。

        俯仰角就是机头与水平面的夹角,机头高于水平面为正,低于水平面为负。横滚角非严格意义下看,就是机体右侧与水平面的夹角,机体右侧低于水平面为正,高于水平面为负。航向角就是机头与正北的夹角,顺时针转动为正,逆时针转动为负。

       可以笼统地认为,俯仰角就是由机身绕Y轴转动产生的,横滚角就是由机身绕X轴转动产生的,航向角就是由机身绕Z轴转动产生的。而X轴陀螺测量的就是机身绕X轴转动的角速度,Y轴陀螺测量的就是机身绕Y轴转动的角速度,Z轴陀螺测量的就是机身绕Z轴转动的角速度。

        我们可以用右手法则判断每个轴角速度的正方向。把右手拇指指向机头方向,另外四指弯曲的方向就是X角速度的正方向(顺时针转动)。把右手拇指指向机体右侧,另外四指弯曲的方向就是Y轴角速度的正方向(顺时针转动)。把右手拇指指向立轴方向,另外四指弯曲的方向就是Z轴角速度的正方向(顺时针转动)。

        机身放平时,横滚角为零。只绕X轴顺时针转动机身,X轴角速度基本都是正数,横滚角在变大,当机体右侧指向地面时,横滚角为90度,当机体翻过肚子,肚子平行于地面,机体右侧再次与水平面重合时,横滚角为180(正向最大值)。只绕X轴逆时针转动机身,X轴角速度基本都是负数,横滚角在变小,当机体右侧指向天空时,横滚角为-90度,当机体翻过肚子,肚子平行于地面,机体右侧再次与水平面重合时,横滚角为-180度(负向最大值)。

        机身放平时,俯仰角为零。只绕Y轴顺时针转动机身,Y轴角速度基本都是正数,俯仰角在变大,当机头指向天空时,俯仰角为90度(正向最大值),然后继续顺时针转动机身,俯仰角开始减小,70,60,40,20,当机头再次与水平面重合时,俯仰角为零。只绕Y轴逆时针转动机身,Y轴角速度基本都是负数,俯仰角在减小,当机头指向地面时,俯仰角为-90度(负向最小值),然后继续逆时针转动机身,俯仰角开始增大,-70,-60,-40,-20,当机头再次与水平面重合时,俯仰角为零。

        机身放平,机头对准正北时,航向角为零。只绕Z轴顺时针转动机身,Z轴角速度基本都是正数,航向角在变大,当机头指向正东时,航向角为90度,当机头指向正南时,航向角为180度(正向最大值)。只绕Z轴逆时针转动机身,Z轴角速度基本都是负值,航向角在变小,当机头指向正西时,航向角为-90度,当机头指向正南时,航向角为-180度(负向最小值)。






7
 楼主| 发表于 2018-1-26 09:04 | 只看该作者
xxdcq 发表于 2018-1-25 19:13
姿态计算里的各加速度需要先归一化处理,而且缺少积分运算,只有比例运算是不稳定的我给你贴一个完全正确的 ...

1、多谢指点。

2,加速度计归一化可以在代入参数之前就做了。放到函数里头做也可以。这个代码是发帖的时候随手写的,限于大家都知道的原因,不可能一下子就把真正在用的代码放出来(尽管几乎就是一样的,但还是换个表达方式稳妥一些)。

3、积分项加不加没有绝对性的影响。
反正我是上飞机验证过了。多了积分项,KI从-100调到正100,一点变化都看不出来。
当然,可能在别的飞机上情况不一样,但目前试过几个机型,还真的没有发现任何差异。


8
 楼主| 发表于 2018-1-26 10:27 | 只看该作者

坐标系对准

    加速度计测量的是重力加速度的分量,而非运动加速度的分量。
    假定陀螺仪的坐标系和机体坐标系完全重合,飞机水平放置的时候,Z轴加速度计的读数为+1g,把飞机翻过来肚子朝上,Z轴加速度计的读数为-1g。飞机水平放置(无论肚子朝上朝下),X轴加速度计和Y轴加速度计的读数都接近于零。

类似的,假定陀螺仪的坐标系和机体坐标系完全重合,把机头垂直向下,X轴加速度计的读数为1g,把机头垂直朝上,X轴加速度计读数为-1g。假定陀螺仪的坐标系和机体坐标系完全重合,把机体右侧垂直向下,Y轴加速度计的读数为1g,把机体右侧垂直朝上,Y轴加速度计的读数为-1g

这样,可以提供一种对准坐标系的方法。

实际使用中,陀螺仪的坐标系和机体坐标系不一定完全一致。当然,两个坐标系的三个轴是肯定可以相互覆盖的(不然正交误差太大,陀螺仪没法用)。按标准机体坐标系的坐标轴顺序(X,Y,Z)来看,陀螺仪(加速度计部分)的坐标轴顺序可能是:(x,y,z),(-x,y,-z),(-x,-y,z),等等,等等。

这时候,可以按照这种方法确定加速度计的正确顺序。假定姿态计算函数为ATTI(GXGYGZAXAYAZ),GXGYGZ为机体坐标系下的三轴角速度(单位为弧度),AXAYAZ为机体坐标系下的三轴重力加速度分量(归一化到-1~1之间)。这时定义陀螺仪(典型的就是MPU6050)的角速度读数(单位为弧度)为gx,gy,gz,加速度计读数(归一化到-1~1之间)为ax,ay,az

首先确定加速度计的坐标系。

把飞机水平放置,看哪个轴读数绝对值最大(另外两个轴读数接近于零),这个轴的读数就是机体坐标系的Z轴的读数 。如:ax=-1,ay=0.05,az=-0.1。这时候可以确认,AZ=fabs(ax)。再看ax为负,就可以确认AZ=-ax

把飞机机头垂直朝下,看哪个轴读数绝对值最大(另外两个轴读数接近于零),这个轴的读数就是机体坐标系的X轴的读数。如(已经和上面的坐标系不一样),ax=0,ay=1az=0.01。这时候,可以确认,AX=fabs(ay)。再看ay为正,就可以确认AX=ay

把机体右侧垂直向下,看哪个轴读数绝对值最大(另外两个轴读数接近于零),这个轴的读数就是机体坐标系的Y轴的读数。如(已经和上面的两个例子不一样),ax=0,ay=-1,az=0。这时候可以确认,AY=fabs(ay)。再看ay为负,就可以确认AY=-ay

再来确定陀螺仪的坐标系。

这时候要用到前面讲的机体绕坐标轴转动顺序(右手法则)。

把飞机水平放置,把飞机绕X轴顺时针转动(另外两个轴尽量保持原样),看哪个轴读数绝对值最大(另外两个轴读数接近于零),这个轴的读数就是机体坐标系的X轴的读数。如,gx=10,gy=0.5,gz=-0.1。这时候,可以确认GX=fabs(gx)。再看gx为正,就可以确认GX=gx

把飞机水平放置,把飞机绕Y轴顺时针转动(另外两个轴尽量保持原样),看哪个轴读数绝对值最大(另外两个轴读数接近于零),这个轴的读数就是机体坐标系的Y轴的读数。如(和前一个例子的坐标系不同),gx=-30,gy=1,gz=-2。这时候可以确认GY=fabs(gx)。再看gx为负,就可以确认GY=-gx

把飞机水平放置,把飞机绕Z轴顺时针转动(另外两个轴尽量保持原样),看哪个轴读数绝对值最大(另外两个轴读数接近于零),这个轴的读数就是机体坐标系的Z轴的读数。如(和前面两个例子的坐标系不同),gx=0,gy=0.5,gz=15。这时候可以确认GZ=fabs(gz)。再看gz为正,就可以确认GZ=gz



9
 楼主| 发表于 2018-1-26 10:32 | 只看该作者

磁强计的坐标系对准

本帖最后由 空中飞鹅 于 2018-1-26 22:47 编辑

       前面讲了陀螺仪和加速度计的坐标对准,这里来讲一下磁强计的坐标对准。

        磁强计,测量的是地磁的分量。

        和陀螺仪类似,磁强计自身的坐标系不一定与机体坐标系完全一致(但必须能彼此覆盖减小正交误差),有可能是磁强计自身x轴,就是机体坐标系的Y轴,等等。

        假设地磁在机体坐标系上的分量为MX,MY,MZ,磁强计自身坐标系的读数为mx,my,mz。      

        把飞机放平,机头对准北方,这时候,机体坐标系X轴的读数为正的最大值,机体坐标系Y轴的读数比其他两个读数明显小得多(接近于零),剩下一个就是机体坐标系的Z轴。

        例如:mx=50,my=400,mz=100,此时,可以很快判断出来:MX=my,MY=mx,MZ=MZ。

10
 楼主| 发表于 2018-1-27 00:07 | 只看该作者

四元数与姿态计算

本帖最后由 空中飞鹅 于 2018-1-27 00:18 编辑

    一、四元数是个什么东西?

    四元数纯粹就是数学上的东西,为了简化姿态相关计算而引入的一套规格,没有太多实质的物理意义。在民用无人机领域,基本不用去深究理论的东西,只要知道有些什么东西,然后怎么用就好了。研究多了,隔几天就忘了,以后还得重新研究,浪费时间浪费精力。

    二、四元数与姿态矩阵

    北东地坐标系(或叫导航坐标系)与机体坐标系之间的坐标变换,就是通过姿态矩阵来实现的。这个姿态矩阵除了可以用三角函数来表示,还可以用四元数来表示。

Tn2b=[T11,T12,T13;T21,T22,T23,T31,T32,T33]

=[cosPitch*cosYAW,cosPitch*sinYAW,-sinPitch;sinRoll*sinPitch*cosYAW-cosRoll*sinYAW,sinPitch*sinRoll*sinYAW+cosRoll*cosYAW,cosPitch*sinRoll;sinPitch*cosRoll*cosYAW+sinRoll*sinYAW,sinPitch*cosRoll*sinYAW-sinRoll*cosYAW,cosRoll*cosPitch]

=[q0*q0+q1*q1-q2*q2-q3*q3,2*(q1*q2+q0*q3),2*(q1*q3-q0*q2);
  2*(q1*q2-q0*q3),q0*q0-q1*q1+q2*q2-q3*q3,2*(q2*q3+q0*q1);
2*(q1*q3+q0*q2),2*(q2*q3-q0*q1),q0*q0-q1*q1-q2*q2+q3*q3]

    利用矩阵元素对应相等的规则,可以推知,

YAW=atan(2*(q1*q2+q0q3),q0*q0+q1*q1-q2*q2-q3*q3)

Pitch=asin(2*(q1*q3-q0*q2))

Roll=atan(2*(q2*q3+q0*q1),q0*q0-q1*q1-q2*q2+q3*q3))

    三、四元数微分方程

     代数式:Qdot = 1/2 *Q * w,其中w在实际应用中就是机体坐标系的角速度矢量。Q是四元数矢量(q0,q1,12,q3),Qdot是Q的导数。

     写成矩阵形式:

   [q0_dot,q1_dot,q2_dot,q3_dot]

=  1/2*[q0,-q1,-q2,-q3;q1,q0,-q3,q2;a2,a3,q0,-q1;q3,-q2,q1,q0]*[0,wx,wy,wz]

=  1/2*[q,-wx,-wy,-wz;wx,0,wz,-wy;wy,-wz,0,wx;wz,wy,-wx,0]*[q0,q1,q2,q3]

    使用一阶龙格库塔数值算法,可得:

  q0 = q0 + (-q1*wx - q2*wy - q3*wz)*0.5*dt;
   q1 = q1 + (q0*wx + q2*wz - q3*wy)*0.5*dt;
   q2 = q2 + (q0*wy - q1*wz + q3*wx)*0.5*dt;
   q3 = q3 + (q0*wz + q1*wy - q2*wx)*0.5*dt;

    于是,只要知道一个初始四元数,然后规定一个递推的时间间隔dt,再结合实时测量得到的机体坐标系角速度,就能得到实时的四元数。

    飞机起飞前有一段时间几乎就是完全静止的。这时可以使用加速度计和磁强计来计算得到初始姿态角。然后利用下面的公式来得到初始四元数。


  Q0=cos(YAW/2)*cos(Pitch/2)*cos(Roll/2)+sin(YAW/2)*sin(Pitch/2)*sin(Roll/2)
Q1=cos(YAW/2)*cos(Pitch/2)*sin(Roll/2)-sin(YAW/2)*sin(Pitch/2)*cos(Roll/2)
Q2=cos(YAW/2)*sin(Piitch/2)*cos(Roll/2)+sin(YAW/2)*cos(Roll/2)*sin(Roll/2)
Q3=cos(YAW/2)*cos(Pitch/2)*cos(Roll/2)-cos(YAW/2)*sin(Pitch/2)*sin(Roll/2)

    之后再代回上面姿态矩阵元素表示姿态角的公式,就能得到三个实时的姿态角。

11
 楼主| 发表于 2018-1-27 01:17 | 只看该作者

常见姿态算法

本帖最后由 空中飞鹅 于 2018-1-27 08:57 编辑

          常见姿态算法主要有:传统型互补滤波、改进型互补滤波、PI互补、梯度下降、姿态角做量测值的卡尔曼滤波,大型卡尔曼滤波(7阶以上)。

1、传统互补滤波

       使用陀螺积分,可以得到三个姿态角,使用加速度计和磁强计,也可以得到三个姿态角,这时候把对应角加权相加就好。
      这种方法在早期用得比较多,但现在几乎没人用了。小角度的情况下,精度还可以,但角度一大,误差明显加大。

伪码如下:

          roll_gyro +=gy*dt;
          pitch_gyro+=gx*dt;
          yaw_gyro+=gz*dt;

          if(roll_gyro>PI)roll_gyro-=PI;
          if(roll_gyro<-PI)roll_gyro+=PI;

          if(pitchl_gyro>PI)pitch_gyro-=PI;
          if(pitch_gyro<-PI)pitch_gyro+=PI;     

          if(yaw_gyro>PI)yaw_gyro-=PI;
          if(yaw_gyro<-PI)yaw_gyro+=PI;      

          roll_acc=atan2(ay,az);
          pitch_acc=atan2(-ax,sqrt(ay*ay+az*az));
          yaw=atan2(-(MY*cos(roll_acc)-MZ*sin(roll_acc),sin(pitch_acc)*(MY*sin(roll_acc)+MZ*cos(roll_acc)+MX*cos(pitch_acc));




         g=sqrt(ax*ax+ay*ay+az*az);
         k1=fabs(g-g0)/(g+g0);
         k2=1-k1;


         roll=roll_gyro*k2+roll_acc*k1;
         pitch=pitch_gyro*k2+pitch_acc*k1;
          yaw=yaw_gyro*k2+yaw_acc*k1;




2、改进型互补滤波


       这种算法可以认为是卡尔曼算法的一个变种,主要用在云台上。飞机上用得不多。


        roll+=gy*dt + k*(roll_acc-roll);
        pitch+=gx*dt + k*(pitch_acc-pitch);
        yaw+=gz*dt + k*(yaw_acc-yaw);




3、PI互补
   
      这种算法是目前全世界使用最多的。九成以上的小飞机都使用这种算法。优点是静态精度很高,缺点是容易受震动影响。
      这是一种四元数算法。首先设定初始姿态为(0,0,0),得到初始四元数(1,0,0,0)。
      接着用四元数表示的加速度矢量与加速度计测得的机体坐标系重力加速度分量叉乘,得到角度误差(时间间隔很小时也可以认为是角速度误差)。
      在10楼的《坐标变换与姿态角计算》中,有提到:


      [ax,ay,az]=Tn2b*[0,0,g]

     展开可得:ax=T11*g;ay=T21*g,az=T31*g;
    把ax,ay,az单位化,即:norn=sqrt(ax*ax+ay*ay+az*az),ax=ax/norm,ay=ay/norm,az=az/norm,
    重力加速度也被单位化,于是有,ax=T11,ay=T21,az=T31。此时的T11,T21和T31由四元数表示。
   
    加速度叉乘得到的误差,可以用来修正陀螺角速度的累积误差。修正过后的角速度,代入四元数微分方程中,就可以递推新的四元数,然后表示出姿态角。

4、梯度下降


     网上资料很多,代码也不列出来了。
      计算量比PI互补大,动态精度比PI互补高,静态精度比PI互补低,但受震动影响较小。
     这种算法在高端飞控上比较常见。它比卡尔曼算法运算量少,精度不比卡尔曼差,所以性价比很高。
     需要注意的是,参数中,gx,gy,gz是机体坐标系的角速度,单位为弧度,ax,ay,az和mx,my,mz可以使用任意量纲,因为进到函数之后会有归一化操作。
   
5、姿态角做量测值的卡尔曼滤波

      这种算法比较鸡肋。
      运算量比梯度下降大很多,精度又不比梯度下降高多少。
      除了写论文能玩一下,其他场合下纯属垃圾。

6、大型卡尔曼滤波(7阶以上)

     在低端民用领域,这种算法实用性也不大。
     至少要在F4及以上处理器,才能跑得动。
     后面有必要再去展开。
     此处从略。





12
 楼主| 发表于 2018-1-27 10:35 | 只看该作者

PI互补代码详解

        想了下,还是对整个PI互补滤波的代码做一次梳理。

        首先要知道,PI互补是一个老外提出来的。它的基础就是梯度下降。

        而梯度下降也是另一个老外提出来的。

        两个老外一开始就把代码给我们写好了,国人哪怕在上面玩出花来,也是用了别人的现成程序。

        我自己虽然用了别的程序(暂时不会贴出来),但也是在老外的代码的基础上简化而来,这里就不昧心装逼了。

         主要有两个函数,第一个只是用陀螺和加速度计融合,第二个用于陀螺、加速度计与磁场融合。

         第一个函数可以计算俯仰角和横滚角,第二个函数除了俯仰角和横滚角,还可以计算航向角。

         这里有个坑,很大的坑。

         没有实战经验的人一般都不会想到。

        所以绝大部分的网上资料没有提到。

        哪怕陀螺、加速度计和磁强计都全了,还是应该把两个函数单独使用。不要只用能算三个角的那个函数。因为,只要磁场一受干扰,你用那个三个角的函数算出来的姿态肯定是歪的,飞机直接就掉下来了,你其他方面再牛逼,也没得玩了。

        横滚角和俯仰角要单独计算,航向角也是要另外计算。

        你可以都是用PI互补来计算,也可以只使用PI互补来算横滚俯仰,而用梯度下降来算航向。用PI互补和梯度下降算航向的时候,确实会重复计算一次横滚俯仰的关联数据,因为没有横滚俯仰的数据,航向角是怎么也算不出来的,能算出来也是不准的。

         先说第一个函数。
         使用加速度计去修正陀螺仪。
         step1:计算初始四元数。飞机放平,初始姿态为(0,0,YAW),此时也可以认为YAW就是0。只算姿态角的时候,航向角其实是什么都不要紧。吃屎姿态角为(0,0,0),代入公式

Q0=cos(YAW/2)*cos(Pitch/2)*cos(Roll/2)+sin(YAW/2)*sin(Pitch/2)*sin(Roll/2)
Q1=cos(YAW/2)*cos(Pitch/2)*sin(Roll/2)-sin(YAW/2)*sin(Pitch/2)*cos(Roll/2)
Q2=cos(YAW/2)*sin(Piitch/2)*cos(Roll/2)+sin(YAW/2)*cos(Roll/2)*sin(Roll/2)
Q3=cos(YAW/2)*cos(Pitch/2)*cos(Roll/2)-cos(YAW/2)*sin(Pitch/2)*sin(Roll/2)


可以求得Q=(1,0,0,0)。


step2:求四元数表示的加速度矢量。


前面已经讲解过一次,这里复习一下。


[ax,ay,az]=Tn2b*[0,0,g]


令norm=sqrt(ax*ax+ay*ay+az*az),ax=ax/notm,ay=ay/norm,az=az/norm,可得:


[ax,ay,az]=Tn2b*[0,0,1],展开后可得:


ax=T13,ay=T23,az=T33。


用四元数表示,可得:


       T13= 2 * (q1 * q3 - q0 * q2);
       T23 = 2 * (q2 * q3 + q0 * q1);
        T33 = q3 * q3 - q2 * q2 - q1 * q1 + q0 * q0;



step3:加速度修正角速度

飞机动起来之后,T13,T23,T33和ax,ay,az不再相等,而是存在一个夹角。可以等效认为这个夹角就是由陀螺仪旋转误差产生的。利用矢量叉乘等于矢量之间的误差,得到:

       ex = (ay*T33 - az*T23);
        ey = (az*T13 - ax*T33);
        ez = (ax*T23 - ay*T13);

这时的误差还是角度量纲,令ex=ex/dt,ey=ey/dt,ez=ez/dt。
误差变成角速度上的误差,把这个误差加到陀螺仪角速度上:

可以修正角速度:

gx=gx+ex*kp;
gy=gy+ey*kp;
gz=gz+ez*kp;

有人说加上误差积分性能会更加稳定,那这里也加上
   exInt = exInt + ex * Ki;                                                                  
   eyInt = eyInt + ey * Ki;  
   ezInt = ezInt + ez * Ki;  

继续修正角速度:

gx=gx+exInt;
gy=gy+eyInt;
gz=gz+ezInt


step4:利用四元数求解姿态角

之后利用四元数微分方程和一阶龙格库塔数值解法,可以得到横滚角和俯仰角:

        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;

        norm = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);//必须要有四元数归一化的步骤,不然会引入很大误差

        q0 = q0 / norm;
        q1 = q1 / norm;
        q2 = q2 / norm;
        q3 = q3 / norm;

       T13 = 2 * (q1 * q3 - q0 * q2);
       T23 = 2 * (q2 * q3 + q0 * q1);
        T33 = q3 * q3 - q2 * q2 - q1 * q1 + q0 * q0;

        Pitch= _asin(-T13);
        Roll = _atan2(T23, T33);


step5:参数调整

需要说明,在叉乘得到误差的时候,误差量纲发生了变化,但除以的dt是个常数,后面误差值乘以的kp和ki,在程序正常跑动的时候也是常数,所以可以把这个dt合并到kp和ki中去。统一去调这个kp,ki就好。一开始只用比例项就好,积分项暂时屏蔽掉。

dt一般用400HZ或500HZ,kp可以在0.01到100间调整,建议从1开始单向调整,先调大,如变成2,变成3,等。大了不行再调小,如变0.1,变0.01,等。
很简单,改了一个kp,先看静态数值,如果一个参数在合适范围,飞机水平静止放置10分钟,两个姿态角基本还是在0度附近,波动不超过1度。静态没问题,再看动态,把飞机用力摇晃一会,再水平放好,静看一两分钟,如果横滚俯仰可以在很短的时间内(最多几秒)恢复到0度附近,说明这个参数可以用。

只有kp的时候,kp调大表示更加相信加速度计,kp调小表示更加相信陀螺仪。其实也是符合互补滤波加权融合的原理。

调到比较满意的时候,再上飞机飞一下看看(前面没问题就可以上飞机了只是效果不一定很好),进行更细致的调整。

比例项我自己是不加的。这里就不去描述调整的过程了。

临时有事,下一个帖子再去讲加上磁强计之后的PI互补算法。

13
 楼主| 发表于 2018-1-27 14:35 | 只看该作者

保留

本帖最后由 空中飞鹅 于 2018-1-29 08:48 编辑

        正在从头修改已经发出来的内容。还没修到这段。暂时保留。
14
发表于 2018-2-5 10:27 | 只看该作者
好贴,学习了
15
发表于 2018-4-5 00:02 | 只看该作者
真的是强 学习了学习了
16
发表于 2018-6-5 23:21 | 只看该作者
真是纯干货啊!
来自苹果客户端来自苹果客户端
17
发表于 2018-6-13 01:39 | 只看该作者
18
发表于 2018-7-26 10:59 | 只看该作者
楼主还来不来呀,等 你的讲解
19
发表于 2018-7-30 12:01 | 只看该作者
20
发表于 2018-7-30 20:31 | 只看该作者
仰望大神.........................................................
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

关闭

【站内推荐】上一条 /1 下一条

快速回复 返回顶部 返回列表