轉(zhuǎn)載請(qǐng)注明出處：http://blog.csdn.net/hongbin_xu 或 http://hongbin96.com/
文章鏈接：http://blog.csdn.net/hongbin_xu/article/details/56846490 或 http://hongbin96.com/111

上次我們討論了姿態(tài)解算基礎(chǔ)理論以及幾個(gè)比較重要的公式的一些推導(dǎo)，如果有不熟悉的請(qǐng)點(diǎn)擊這里打開鏈接。這次來(lái)介紹一些實(shí)際的姿態(tài)解算算法吧！
一般在程序中，姿態(tài)解算的方式有兩種：一種是歐拉角法，一種是四元數(shù)法。這里不介紹歐拉角法，只介紹四元數(shù)法，如有興趣可以去查找相關(guān)資料。

互補(bǔ)濾波算法：
顧名思義，是多組數(shù)據(jù)結(jié)合互補(bǔ)，并進(jìn)行濾波處理穩(wěn)定輸出，得到姿態(tài)的算法。而我們使用的傳感器就是加速度計(jì)和陀螺儀。加速度計(jì)用于測(cè)量加速度，陀螺儀用于測(cè)量角速度。 加速度計(jì)的靜態(tài)穩(wěn)定性更好，而在運(yùn)動(dòng)時(shí)其數(shù)據(jù)相對(duì)不可靠；陀螺儀的動(dòng)態(tài)穩(wěn)定性更好，但是靜止時(shí)數(shù)據(jù)相對(duì)不可靠。所以，我們可以通過加速度計(jì)的輸出來(lái)修正陀螺儀的漂移誤差，換句話說，通過加速度計(jì)來(lái)修正陀螺儀。

這個(gè)是我在網(wǎng)上找到的說明互補(bǔ)濾波法的框圖：（原圖下載：點(diǎn)擊這里打開)

首先，我們?nèi)《▽?dǎo)航坐標(biāo)系n中標(biāo)準(zhǔn)重力加速度g，定義為，那么將導(dǎo)航坐標(biāo)系n下的 轉(zhuǎn)換為載體坐標(biāo)系b下的：。
這里用到了，前面在推導(dǎo)基礎(chǔ)公式時(shí)導(dǎo)出了使用四元數(shù)表示的旋轉(zhuǎn)矩陣。
公式如下：

但是我們要用的是，所以還要對(duì) 做一個(gè)矩陣逆變換。由于它是正交矩陣，對(duì)于正交矩陣有這個(gè)性質(zhì)：正交矩陣的逆矩陣等于其的轉(zhuǎn)置。所以我們很容易得到：

將上式代入 ，而且我們已知，所以得到：

接著再定義載體坐標(biāo)系b中加速度計(jì)輸出為a，由于前面計(jì)算導(dǎo)航坐標(biāo)系時(shí)我們采用的重力加速度是標(biāo)準(zhǔn)重力加速度，所以還需要對(duì)其進(jìn)行歸一化，才能繼續(xù)運(yùn)算。
設(shè)加速度計(jì)三個(gè)軸的值分別是ax，ay，az。
首先求模：。
歸一化：

根據(jù)框圖的說明再整理一下前面得到的結(jié)果：
標(biāo)準(zhǔn)重力加速度從n系轉(zhuǎn)到b系中的矩陣表示：

b系下加速度計(jì)測(cè)量得到的加速度的矩陣表示：
（備注：這是歸一化之后的值）

對(duì) 和 做向量叉乘，即可得到給陀螺儀的校正補(bǔ)償值e。

表示成矩陣形式更為直觀：

然后再使用PI控制器進(jìn)行濾波，準(zhǔn)確地說事消除漂移誤差，只要存在誤差控制器便會(huì)持續(xù)作用，直至誤差為0。控制的效果取決于P和I參數(shù)，分別對(duì)應(yīng)比例控制和積分控制的參數(shù)。

這里給出PI控制的公式：

是我們要負(fù)反饋給陀螺儀進(jìn)行校正補(bǔ)償?shù)闹?#xff0c; 是比例控制項(xiàng)， 是積分控制項(xiàng)，在程序中采用離散累加計(jì)算。關(guān)于PID控制理論的東西，這里不做贅述。

如框圖中所寫，接下來(lái)將前面得到的補(bǔ)償值加在陀螺儀輸出的數(shù)據(jù)上進(jìn)行校正。

跟著框圖往下走：

將前面的陀螺儀數(shù)據(jù)通過四元數(shù)微分方程轉(zhuǎn)換為四元數(shù)輸出。
因?yàn)橛袔讉€(gè)地方我也沒搞懂，所以就簡(jiǎn)單介紹一下四元數(shù)微分方程，詳細(xì)步驟請(qǐng)查閱秦永元的慣性導(dǎo)航一書（第三篇 9.2.3節(jié)）：
由于載體的運(yùn)動(dòng)，四元數(shù)Q是變量，其參數(shù)可以表示成關(guān)于時(shí)間的函數(shù)。
使用四元數(shù)的三角形式：
為剛體瞬時(shí)繞軸轉(zhuǎn)過的角度， 為歸一化后的位置矢量。
設(shè)角速度為：
對(duì)四元數(shù)的三角表示形式求導(dǎo)：

因?yàn)?， 且
所以：

將上式表示成矩陣形式：

或者
上面的兩個(gè)公式被稱為四元數(shù)微分方程。利用陀螺儀的數(shù)據(jù)進(jìn)行離散累積便可得到四元數(shù)的值，最后再轉(zhuǎn)換成歐拉角形式輸出即可。

再附上代碼：

//////////////////////////////////////////////////////////////////////////////// #define Kp 10.0f // proportional gain governs rate of convergence to accelerometer/magnetometer #define Ki 0.008f // integral gain governs rate of convergence of gyroscope biases #define halfT 0.001f // half the sample period采樣周期的一半float q0 = 1, q1 = 0, q2 = 0, q3 = 0; // quaternion elements representing the estimated orientation float exInt = 0, eyInt = 0, ezInt = 0; // scaled integral error 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數(shù)據(jù)歸一化ax = ax /norm;ay = ay / norm;az = az / norm;// estimated direction of gravity and flux (v and w) 估計(jì)重力方向和流量/變遷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 sensorsex = (ay*vz - az*vy) ; //向量外積在相減得到差分就是誤差ey = (az*vx - ax*vz) ;ez = (ax*vy - ay*vx) ;exInt = exInt + ex * Ki; //對(duì)誤差進(jìn)行積分eyInt = eyInt + ey * Ki;ezInt = ezInt + ez * Ki;// adjusted gyroscope measurementsgx = gx + Kp*ex + exInt; //將誤差PI后補(bǔ)償?shù)酵勇輧x，即補(bǔ)償零點(diǎn)漂移gy = gy + Kp*ey + eyInt;gz = gz + Kp*ez + ezInt; //這里的gz由于沒有觀測(cè)者進(jìn)行矯正會(huì)產(chǎn)生漂移，表現(xiàn)出來(lái)的就是積分自增或自減// 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 quaternionnorm = 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, -2 * q1 * q1 - 2 * q2* q2 + 1)* 57.3; // roll }

這段代碼網(wǎng)上很多地方都可以看見，匿名四軸（老版本的）的程序也是用的這段。
上面也添加了一些注釋，把互補(bǔ)濾波算法的理論部分過了一遍后再來(lái)對(duì)著代碼讀一遍，應(yīng)該不會(huì)覺得有多難了吧。程序中就只是把前面推導(dǎo)的公式一個(gè)個(gè)用上去了而已，如果不懂原理直接使用也不會(huì)有太大問題。
程序中的積分運(yùn)算是采用離散累積的方法運(yùn)算的。PID控制器只是起到了一個(gè)穩(wěn)定數(shù)據(jù)消除漂移誤差的作用。

希望這些筆記能給更多的人提供幫助，歡迎交流指正。

我的前一篇博文講的是姿態(tài)解算的基礎(chǔ)知識(shí)和一些相關(guān)公式的推導(dǎo)，有興趣的可以去查看：http://blog.csdn.net/hongbin_xu/article/details/55667899

最后再分享一些我在網(wǎng)上找到的資料，個(gè)人覺得對(duì)學(xué)習(xí)很有幫助。
鏈接：
http://blog.csdn.net/wkdwl/article/details/52119163
http://blog.csdn.net/nemol1990/article/details/16924745
http://blog.csdn.net/super_mice/article/details/45619945
http://www.openedv.com/thread-42657-1-1.html

總結(jié)

以上是生活随笔為你收集整理的四旋翼姿态解算——互补滤波算法及理论推导的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。