加速度计和陀螺仪数据融合
生活随笔
收集整理的這篇文章主要介紹了
加速度计和陀螺仪数据融合
小編覺得挺不錯的,現在分享給大家,幫大家做個參考.
本帖翻譯自
IMU（加速度計和陀螺儀設備）在嵌入式應用中使用的指南。
這篇文章主要介紹加速度計和陀螺儀的數學模型和基本算法，以及如何融合這兩者，側重算法、思想的討論
介紹
本指南旨在向興趣者介紹慣性MEMS（微機電系統）傳感器，特別是加速度計和陀螺儀以及其他整合IMU（ 慣性測量單元 ）設備。
? 2012-8-22 16:38 上傳 下載附件?(16.48 KB)
IMU單元例子 ：上圖中MCU頂端的 ACC Gyro 6DOF ，名為 USBThumb ，支持USB/串口通信
在這篇文章中我將概括這么幾個基本并且重要的話題：
- 加速度計（accelerometer）檢測什么
- 陀螺儀（gyroscope，也稱作 gyro）檢測什么
- 如何將傳感器ADC讀取的數據轉換為物理單位（加速度傳感器的單位是g，陀螺儀的單位是 度/秒）
- 如何結合加速度傳感器和陀螺儀的數據以得到設備和地平面之間的傾角的準確信息
在整篇文章中我盡量將數學運算降低到最少。如果你知道什么是正弦、余弦、正切函數，那無論你的項目使用哪種平臺你應該都會明白和運用這篇文章中的思想，這些平臺如Arduino、Propeller、Basic Stamp、Ateml芯片、PIC芯片等等。總有些人認為使用IMU單元需要復雜的數學運算（復雜的FIR或IIR濾波，如卡爾曼濾波，Parks-McClellan濾波等）。你如果研究這些會得到很棒且很復雜的結果。我解釋事情的方式，只需要基本的數學。我非常堅信簡單的原則。我認為一個簡單的系統更容易操作和監控，另外許多嵌入式設備并不具備能力和資源去實現需要進行矩陣運算的復雜算法。
我會用我設計的一個新IMU模塊—— Acc_Gyro Accelerometer + Gyro IMU 作為例子。在下面的例子中我們會使用這個設備的參數。用這個模塊作為介紹非常合適，因為它由3個設備組成：
- LIS331AL ( datasheet ) – 3軸 2G 模擬加速度計
- LPR550AL ( datasheet ) – 雙軸（俯仰、翻滾） 500°/s 加速度傳感器
- LY550ALH ( datasheet ) –單軸（偏航）陀螺儀 最后這個設備在這篇介紹中不使用，不過他在DCM Matrix implementation 中有重要作用
它們一起組成了一個6自由度的慣性測量單元。這是個花哨的名字！然而，在花哨的名字后面是個非常有用的設備組合，接下來我們會詳細介紹之。
第一部分 加速度計
要了解這個模塊我們先從加速度計開始。當我們在想象一個加速度計的時候我們可以把它想作一個圓球在一個方盒子中。你可能會把它想作一個餅干或者甜圈，但我就把它當做一個球好了：
? 2012-8-22 16:38 上傳 下載附件?(15.85 KB)
我們假定這個盒子不在重力場中或者其他任何會影響球的位置的場中，球處于盒子的正中央。你可以想象盒子在外太空中，遠離任何天體，如果很難想象，那就當做盒子在航天飛機中，一切東西都處于無重力狀態。在上面的圖中你可以看到我們給每個軸分配了一對墻（我們移除了Y+以此來觀察里面的情況）。設想每面墻都能感測壓力。如果我們突然把盒子向左移動（加速度為1g=9.8m/s^2），那么球會撞上X-墻。然后我們檢測球撞擊墻面產生的壓力，X軸輸出值為-1g。
? 2012-8-22 16:38 上傳 下載附件?(19.23 KB)
請注意加速度計檢測到得力的方向與它本身加速度的方向是相反的。這種力量通常被稱為 慣性力或假想力 ?。在這個模型中你你應該學到加速度計是通過間接測量力對一個墻面的作用來測量加速度的，在實際應用中，可能通過彈簧等裝置來測量力。這個力可以是加速度引起的，但在下面的例子中，我們會發現它不一定是加速度引起的。
如果我們把模型放在地球上，球會落在Z-墻面上并對其施加一個1g的力，見下圖：
? 2012-8-22 16:38 上傳 下載附件?(19.54 KB)
在這種情況下盒子沒有移動但我們任然讀取到Z軸有-1g的值。球在墻壁上施加的壓力是由引力造成的。在理論上，它可以是不同類型的力量 - 例如，你可以想象我們的球是鐵質的，將一個磁鐵放在盒子旁邊那球就會撞上另一面墻。引用這個例子只是為了說明加速度計的本質是檢測力而非加速度。只是加速度所引起的慣性力正好能被加速度計的檢測裝置所捕獲。
雖然這個模型并非一個MEMS傳感器的真實構造，但它用來解決與加速度計相關的問題相當有效。實際上有些類似傳感器中有金屬小球，它們稱作傾角開關，但是它們的功能更弱，只能檢測設備是否在一定程度內傾斜，卻不能得到傾斜的程度。
到目前為止，我們已經分析了單軸的加速度計輸出，這是使用單軸加速度計所能得到的。三軸加速度計的真正價值在于它們能夠檢測全部三個軸的慣性力。讓我們回到盒子模型，并將盒子向右旋轉45度。現在球會與兩個面接觸：Z-和X-，見下圖：
? 2012-8-22 16:38 上傳 下載附件?(34.03 KB)
0.71g這個值是不是任意的，它們實際上是1/2的平方根的近似值。我們介紹加速度計的下一個模型時這一點會更清楚。
在上一個模型中我們引入了重力并旋轉了盒子。在最后的兩個例子中我們分析了盒子在兩種情況下的輸出值，力矢量保持不變。雖然這有助于理解加速度計是怎么和外部力相互作用的，但如果我們將坐標系換為加速度的三個軸并想象矢量力在周圍旋轉，這會更方便計算。
? 2012-8-22 16:38 上傳 下載附件?(7.19 KB)
請看看在上面的模型，我保留了軸的顏色，以便你的思維能更好的從上一個模型轉到新的模型中。想象新模型中每個軸都分別垂直于原模型中各自的墻面。矢量R是加速度計所檢測的矢量（它可能是重力或上面例子中慣性力的合成）。RX，RY，RZ是矢量R在X，Y，Z上的投影。請注意下列關系：
，R ^ 2 = RX ^ 2 + RY ^ 2 + RZ ^ 2（ 公式1 ）
此公式等價于 三維空間勾股定理 。
還記得我剛才說的1/2的平方根0.71不是個隨機值吧。如果你把它們代回上式，回顧一下重力加速度是1g，那我們就能驗證：
1 ^ 2 =（SQRT（1/2））^ 2 + 0 ^ 2 +（SQRT（1/2））^ 2
在 公式1 中簡單的取代： R=1, Rx = -SQRT(1/2), Ry = 0 , Rz = -SQRT(1/2)
經過一大段的理論序言后，我們和實際的加速度計很靠近了。RX，RY，RZ值是實際中加速度計輸出的線性相關值，你可以用它們進行各種計算。
在我們運用它之前我們先討論一點獲取加速度計數據的方法。大多數加速度計可歸為兩類：數字和模擬。數字加速度計可通過I2C，SPI或USART方式獲取信息，而模擬加速度計的輸出是一個在預定范圍內的電壓值，你需要用ADC（模擬量轉數字量）模塊將其轉換為數字值。我將不會詳細介紹ADC是怎么工作的，部分原因是這是個很廣的話題，另一個原因是不同平臺的ADC都會有差別。有些MCU具有內置ADC模塊，而有些則需要外部電路進行ADC轉換。不管使用什么類型的ADC模塊，你都會得到一個在一定范圍內的數值。例如一個10位ADC模塊的輸出值范圍在0 .. 1023間，請注意，1023 = 2 ^ 10 -1。一個12位ADC模塊的輸出值范圍在0 .. 4095內，注意，4095 = 2 ^ 12-1。
我們繼續，先考慮下一個簡單的例子，假設我們從10位ADC模塊得到了以下的三個軸的數據：
AdcRx = 586
AdcRy = 630
AdcRz = 561
每個ADC模塊都有一個參考電壓，假設在我們的例子中，它是3.3V。要將一個10位的ADC值轉成電壓值，我們使用下列公式：
VoltsRx = AdcRx * VREF / 1023
小注：8位ADC的最大值是255 = 2 ^ 8 -1，12位ADC最大值是4095 = 2 ^ 12 -1。
將3個軸的值代入上式，得到：
VoltsRx = 586 * 3.3 / 1023 =～1.89V（結果取兩位小數）
VoltsRy = 630 * 3.3 / 1023 =～2.03V
VoltsRz = 561 * 3.3 / 1023 =～1.81V
每個加速度計都有一個零加速度的電壓值，你可以在它的說明書中找到，這個電壓值對應于加速度為0g。通過計算相對0g電壓的偏移量我們可以得到一個有符號的電壓值。比方說，0g電壓值 VzeroG= 1.65V，通過下面的方式可以得到相對0g電壓的偏移量:
DeltaVoltsRx = 1.89V - 1.65V = 0.24V
DeltaVoltsRy = 2.03V - 1.65V = 0.38V
DeltaVoltsRz = 1.81V - 1.65V = 0.16V
現在我們得到了加速度計的電壓值，但它的單位還不是g（9.8m/s^2），最后的轉換，我們還需要引入加速度計的靈敏度（Sensitivity），單位通常是 mV/g。比方說，加速度計的靈敏度 Sensitivity= 478.5mV / g = 0.4785V /g。靈敏度值可以在加速度計說明書中找到。要獲得最后的單位為g的加速度，我們使用下列公式計算：
RX = DeltaVoltsRx /Sensitivity
RX = 0.24V / 0.4785V / G =～0.5g
RY = 0.38V / 0.4785V / G =～0.79g
RZ = 0.16V / 0.4785V / G =～0.33g
當然，我們可以把所有的步驟全部放在一個式子里，但我想通過介紹每一個步驟以便讓你了解怎么讀取一個ADC值并將其轉換為單位為g的矢量力的分量。
Rx = (AdcRx * Vref / 1023 – VzeroG) / Sensitivity ?(公式2)
Ry = (AdcRy * Vref / 1023 – VzeroG) / Sensitivity
Rz = (AdcRz * Vref / 1023 – VzeroG) / Sensitivity
現在我們得到了慣性力矢量的三個分量，如果設備除了重力外不受任何外力影響，那我們就可以認為這個方向就是重力矢量的方向。如果你想計算設備相對于地面的傾角，可以計算這個矢量和Z軸之間的夾角。如果你對每個軸的傾角都感興趣，你可以把這個結果分為兩個分量：X軸、Y軸傾角，這可以通過計算重力矢量和X、Y軸的夾角得到。計算這些角度比你想象的簡單，現在我們已經算出了Rx，Ry，Rz的值，讓我們回到我們的上一個加速度模型，再加一些標注上去：
? 2012-8-22 16:38 上傳 下載附件?(9.9 KB)
我們感興趣的角度是向量R和X，Y，Z軸之間的夾角，那就令這些角度為Axr，Ayr，Azr。觀察由R和Rx組成的直角三角形：
cos(Axr) = Rx / R , 類似的：
cos(Ayr) = Ry / R
cos(Azr) = Rz / R
從 公式1 我們可以推導出 R = SQRT( Rx^2 + Ry^2 + Rz^2)
通過arccos()函數（cos()的反函數）我們可以計算出所需的角度：
Axr = arccos(Rx/R)
Ayr = arccos(Ry/R)
Azr = arccos(Rz/R)
我們花了大段的篇幅來解釋加速度計模型，最后所要的只是以上這幾個公式。根據你的應用場合，你可能會用到我們推導出來的幾個過渡公式。我們接下來要介紹陀螺儀模塊，并向大家介紹怎么融合加速度計和陀螺儀的數據以得到更精確的傾角值。
但在此之前，我們再介紹幾個很常用的公式：
cosX = cos(Axr) = Rx / R
cosY = cos(Ayr) = Ry / R
cosZ = cos(Azr) = Rz / R
這三個公式通常稱作 方向余弦 ?，它主要表達了單位向量（長度為1的向量）和R向量具有相同的方向。你可以很容易地驗證：
SQRT（cosX ^ 2 + COSY ^ 2 + cosZ ^ 2）= 1
這是個很好的性質，因為它避免了我們一直檢測R向量的模（長度）。通常如果我們只是對慣性力的方向感興趣，那標準化模長以簡化其他計算是個明智的選擇。
第二部分陀螺儀
對于陀螺儀我們將不會像加速度計一樣介紹它的等價盒子模型，而是直接跳到加速度計的第二個模型，通過這個模型我們會向大家介紹陀螺儀是怎么工作的。
? 2012-8-22 16:38 上傳 下載附件?(12.55 KB)
陀螺儀的每個通道檢測一個軸的旋轉。例如，一個2軸陀螺儀檢測繞X和Y軸的旋轉。為了用數字來表達這些旋轉，我們先引進一些符號。首先我們定義：
Rxz – 慣性力矢量R在XZ平面上的投影
Ryz – 慣性力矢量R在YZ平面的上投影
在由Rxz和Rz組成的直角三角形中，運用勾股定理可得：
Rxz^2 = Rx^2 + Rz^2 ，同樣：
Ryz^2 = Ry^2 + Rz^2
同時注意：
R^2 = Rxz^2 + Ry^2 ，這個公式可以 公式1 和上面的公式推導出來，也可由R和Ryz所組成的直角三角形推導出來
R ^ 2 = Ryz ^ 2 + RX ^ 2
在這篇文章中我們不會用到這些公式，但知道模型中的那些數值間的關系有助于理解。
相反，我們按如下方法定義Z軸和Rxz、Ryz向量所成的夾角：
AXZ - Rxz（矢量R在XZ平面的投影）和Z軸所成的夾角
AYZ - Ryz（矢量R在YZ平面的投影）和Z軸所成夾角
現在我們離陀螺儀要測量的東西又近了一步。陀螺儀測量上面定義的角度的變化率。換句話說，它會輸出一個與上面這些角度變化率線性相關的值。為了解釋這一點，我們先假設在t0時刻，我們已測得繞Y軸旋轉的角度（也就是Axz），定義為Axz0，之后在t1時刻我們再次測量這個角度，得到Axz1。角度變化率按下面方法計算：
RateAxz = (Axz1 – Axz0) / (t1 – t0).
如果用度來表示角度，秒來表示時間，那這個值的單位就是 度/秒。這就是陀螺儀檢測的東西。
在實際運用中，陀螺儀一般都不會直接給你一個單位為度/秒的值（除非它是個特殊的數字陀螺儀）。就像加速度計一樣，你會得到一個ADC值并且要用類似 公式2 的式子將其轉換成單位為 度/秒的值。讓我們來介紹陀螺儀輸出值轉換中的ADC部分（假設使用10位ADC模塊，如果是8位ADC，用1023代替255，如果是12為ADC用4095代替1023）。
RateAxz = (AdcGyroXZ * Vref / 1023 – VzeroRate) / Sensitivity? 公式3
RateAyz = (AdcGyroYZ * Vref / 1023 – VzeroRate) / Sensitivity
AdcGyroXZ，AdcGyroYZ - 這兩個值由ADC讀取，它們分別代表矢量R的投影在XZ和YZ平面內里的轉角，也可等價的說，旋轉可分解為單獨繞Y和X軸的運動。
Vref – ADC的參考電壓，上例中我們使用3.3V
VzeroRate – 是零變化率電壓，換句話說它是陀螺儀不受任何轉動影響時的輸出值，對 Acc Gyro 板來說，可以認為是1.23V（此值通常可以在說明書中找到——但千萬別相信這個值，因為大多數的陀螺儀在焊接后會有一定的偏差，所以可以使用電壓計測量每個通道的輸出值，通常這個值在焊接后就不會改變，如果有跳動，在設備使用前寫一個校準程序對其進行測量，用戶應當在設備啟動的時候保持設備靜止以進行校準）。
Sensitivity –陀螺儀的靈敏度，單位mV/(deg/s)，通常寫作mV/deg/s,它的意思就是如果旋轉速度增加1°/s,陀螺儀的輸出就會增加多少mV。 Acc_Gyro 板的靈敏度值是2mV/deg/s或0.002V/deg/s
讓我們舉個例子，假設我們的ADC模塊返回以下值：
AdcGyroXZ = 571
AdcGyroXZ = 323
用上面的公式，在代入 Acc Gyro 板的參數，可得：
RateAxz = (571 * 3.3V / 1023 – 1.23V) / ( 0.002V/deg/s) =~ 306 deg/s
RateAyz = (323 * 3.3V / 1023 – 1.23V) / ( 0.002V/deg/s) =~ -94 deg/s
換句話說設備繞Y軸（也可以說在XZ平面內）以306°/s速度和繞X軸（或者說YZ平面內）以-94°/s的速度旋轉。請注意，負號表示該設備朝著反方向旋轉。按照慣例，一個方向的旋轉是正值。一份好的陀螺儀說明書會告訴你哪個方向是正的，否則你就要自己測試出哪個旋轉方向會使得輸出腳電壓增加。最好使用示波器進行測試，因為一旦你停止了旋轉，電壓就會掉回零速率水平。如果你使用的是萬用表，你得保持一定的旋轉速度幾秒鐘并同時比較電壓值和零速率電壓值。如果值大于零速率電壓值那說明這個旋轉方向是正向。
第三部分 將它們綜合起來。融合加速度計和陀螺儀的數據。
如果你在閱讀這篇文章你可能已經有了或準備購買一個IMU設備，或者你準備用獨立的加速度計和陀螺儀搭建一個。
注：具體的代碼實現和算法測試，請閱讀這篇文章：
http://starlino.com/imu_kalman_arduino.html
在使用整合了加速度計和陀螺儀的IMU設備時，首先要做的就是統一它們的坐標系。最簡單的辦法就是將加速度計作為參考坐標系。大多數的加速度計技術說明書都會指出對應于物理芯片或設備的XZY軸方向。例如，下面就是Acc Gyro板的說明書中給出的XYZ軸方向：
? 2012-8-22 16:39 上傳 下載附件?(35.18 KB)
接下來的步驟是：
- 確定陀螺儀的輸出對應到上述討論的RateAxz，RateAyz值。
- 根據陀螺儀和加速度計的位置決定是否要反轉輸出值
不要設想陀螺儀陀的輸出有XY，它會適應加速度計坐標系里的任何軸，盡管這個輸出是IMU模塊的一部分。最好的辦法就是測試。
接下來的示例用來確定哪個陀螺儀的輸出對應RateAxz。
- 首先將設備保持水平。加速度計的XY軸輸出會是零加速度電壓（ Acc Gyro 板的值是1.65V）
- 接下來將設備繞Y軸旋轉，換句話說就是將設備在XZ平面內旋轉，所以X、Z的加速度輸出值會變化而Y軸保持不變。
-當以勻速旋轉設備的時候，注意陀螺儀的哪個通道輸出值變化了，其他輸出應該保持不變。
- 在陀螺儀繞Y軸旋轉（在XZ平面內旋轉）的時候輸出值變化的就是AdcGyroXZ，用于計算RateAxz
-最后一步，確認旋轉的方向是否和我們的模型對應，因為陀螺儀和加速度的位置關系，有時候你可能要把RateAxz值反向
-重復上面的測試，將設備繞Y軸旋轉，這次查看加速度計的X軸輸出（也就是AdcRx）。如果AdcRx增大（從水平位置開始旋轉的第一個90°），那AdcGyroXZ應當減小。這是因為我們觀察的是重力矢量，當設備朝一個方向旋轉時矢量會朝相反的方向旋轉（相對坐標系運動）。所以，如果你不想反轉RateAxz，你可以在 公式3 中引入正負號來解決這個問題：
RateAxz = InvertAxz * (AdcGyroXZ * Vref / 1023 – VzeroRate) / Sensitivity ,其中InvertAxz= 1 或-1
同樣的方法可以用來測試RateAyz，將設備繞X軸旋轉，你就能測出陀螺儀的哪個輸出對應于RateAyz，以及它是否需要反轉。一旦你確定了InvertAyz，你就能可以用下面的公式來計算RateAyz：
RateAyz = InvertAyz * (AdcGyroYZ * Vref / 1023 – VzeroRate) / Sensitivity
如果對Acc Gyro板進行這些測試，你會得到下面的這些結果：
- RateAxz的輸出管腳是GX4，InvertAxz = 1
- RateAyz輸出管腳是GY4，InvertAyz = 1
從現在開始我們認為你已經設置好了IMU模塊并能計算出正確的Axr，Ayr，Azr值（在第一部分加速度計中定義）以及RateAyz，RateAyz（在第二部分陀螺儀中）。下一步，我們分析這些值之間的關系并得到更準確的設備和地平面之間的傾角。
你可能會問自己一個問題，如果加速度計已經告訴我們Axr，Ayr，Azr的傾角，為什么還要費事去得到陀螺儀的數據？答案很簡單：加速度計的數據不是100%準確的。有幾個原因，還記加速度計測量的是慣性力，這個力可以由重力引起（理想情況只受重力影響），當也可能由設備的加速度（運動）引起。因此，就算加速度計處于一個相對比較平穩的狀態，它對一般的震動和機械噪聲很敏感。這就是為什么大部分的IMU系統都需要陀螺儀來使加速度計的輸出更平滑。但是怎么辦到這點呢？陀螺儀不受噪聲影響嗎？
陀螺儀也會有噪聲，但由于它檢測的是旋轉，因此對線性機械運動沒那么敏感，不過陀螺儀有另外一種問題，比如漂移（當選擇停止的時候電壓不會回到零速率電壓）。然而，通過計算加速度計和陀螺儀的平均值我們能得到一個相對更準確的當前設備的傾角值，這比單獨使用加速度計更好。
接下來的步驟我會介紹一種算法，算法受卡爾曼濾波中的一些思想啟發，但是它更簡單并且更容易在嵌入式設備中實現。在此之前，讓我們先看看我們需要算法計算什么值。所要算的就是重力矢量R=[Rx,Ry,Rz]，它可由其他值推導出來，如Axr，Ayr，Azr或者cosX，cosY，cosZ，由這些值我們能得到設備相對地平面的傾角值，這些關系我們在第一部分已經討論過。有人可能會說-根據第一部分的 公式2 我們不是已經得到Rx，Ry，Rz的值了嗎？是的，但是記住，這些值只是由加速度計數據推導出來的，如果你直接將它們用于你的程序你會得到難以忍受的噪聲。為了避免進一步的混亂，我們重新定義加速度計的測量值：
Racc – 是由加速度計測量到得慣性力矢量，它可分解為下面的分量（在XYZ軸上的投影）：
RxAcc = (AdcRx * Vref / 1023 – VzeroG) / Sensitivity
RyAcc = (AdcRy * Vref / 1023 – VzeroG) / Sensitivity
RzAcc = (AdcRz * Vref / 1023 – VzeroG) / Sensitivity
現在我們得到了一組只來自于加速度計ADC的值。我們把這組數據叫做“vector”，并使用下面的符號：
Racc = [RxAcc,RyAcc,RzAcc]
因為這些Racc的分量可由加速度計數據得到，我們可以把它當做算法的輸入。
請注意Racc測量的是重力，如果你得到的矢量長度約等于1g那么你就是正確的：
|Racc| = SQRT(RxAcc^2 +RyAcc^2 + RzAcc^2),
但是請確定把矢量轉換成下面的矢量非常重要：
Racc(normalized) = [RxAcc/|Racc| , RyAcc/|Racc| , RzAcc/|Racc|].
這可以確保標準化Racc始終是1。
接來下我們引進一個新的向量：
Rest = [RxEst,RyEst,RzEst]
這就是算法的輸出值，它經過陀螺儀數據的修正和基于上一次估算的值。
這是算法所做的事：
-加速度計告訴我們：“你現在的位置是Racc”
我們回答：“謝謝，但讓我確認一下”
-然后根據陀螺儀的數據和上一次的Rest值修正這個值并輸出新的估算值Rest。
-我們認為Rest是當前設備姿態的“最佳值”。
讓我們看看它是怎么實現的。
數列的開始，我們先認為加速度值正確并賦值：
Rest(0) = Racc(0)
Rest和Racc是向量，所以上面的式子可以用3個簡單的式子代替，注意別重復了：
RxEst（0）= RxAcc（0）
RyEst（0）= RyAcc（0）
RzEst（0）= RzAcc（0）
接下來我們在每個等時間間隔T秒做一次測量，得到新的測量值，并定義為Racc（1），Racc（2），Racc（3）等等。同時，在每個時間間隔我們也計算出新的估算值Rest（1）,Rest（2），Rest（3），等等。
假設我們在第n步。我們有兩列已知的值可以用：
Rest(n-1) – 前一個估算值，Rest(0) = Racc(0)
Racc(n) – 當前加速度計測量值
在計算Rest（n）前，我們先引進一個新的值，它可由陀螺儀和前一個估算值得到。
叫做Rgyro，同樣它是個矢量并由3個分量組成：
Rgyro = [RxGyro,RyGyro,RzGyro]
我們分別計算這個矢量的分量，從RxGyro開始。
? 2012-8-22 16:39 上傳 下載附件?(12.55 KB)
首先觀察陀螺儀模型中下面的關系，根據由Rz和Rxz組成的直角三角形我們能推出：
tan(Axz) = Rx/Rz => Axz = atan2(Rx,Rz)
你可能從未用過atan2這個函數，它和atan類似，但atan返回值范圍是（-PI/2,PI/2），atan2返回值范圍是（-PI，PI），并且他有兩個參數。它能將Rx，Rz值轉換成360°（-PI，PI）內的角度。更多信息請閱讀? atan2 .
所以，知道了RxEst（n-1）和RzEst（n-1）我們發現：
Axz(n-1) = atan2( RxEst(n-1) , RzEst(n-1) ).
記住，陀螺儀測量的是Axz角度變化率，因此，我們可以按如下方法估算新的角度Axz（n）：
Axz(n) = Axz(n-1) + RateAxz(n) * T
請記住，RateAxz可由陀螺儀ADC讀取得到。通過使用平均轉速可由得到一個更準確的公式：
RateAxzAvg =（RateAxz（N）+ RateAxz（N-1））/ 2
Axz(n) = Axz(n-1) + RateAxzAvg * T
同理可得：
Ayz(n) = Ayz(n-1) + RateAyz(n) * T
好了，現在我們有了Axz（n），Ayz（n）。現在我們如何推導出RxGyro/RyGyro?根據 公式1 我們可以把Rgyro長度寫成下式：
| Rgyro | = SQRT（RxGyro ^ 2 + RyGyro ^ 2 + RzGyro ^ 2）
同時，因為我們已經將Racc標準化，我們可以認為它的長度是1并且旋轉后保持不變，所以寫成下面的方式相對比較安全：
| Rgyro | = 1
我們暫時采用更短的符號進行下面的計算：
x =RxGyro , y=RyGyro, z=RzGyro
根據上面的關系可得：
x = x / 1 = x / SQRT(x^2+y^2+z^2)
分子分母同除以SQRT（X ^ 2 + Z ^ 2）
x = ( x / SQRT(x^2 + z^2) ) / SQRT( (x^2 + y^2 + z^2) / (x^2 + z^2) )
注意x / SQRT(x^2 + z^2) = sin(Axz), 所以:
x = sin(Axz) / SQRT (1 + y^2 / (x^2 + z^2) )
將SQRT內部分式的分子分母同乘以z^2
x = sin(Axz) / SQRT (1 + y^2??* z ^2 / (z^2 * (x^2 + z^2)) )
注意 z / SQRT(x^2 + z^2) = cos(Axz)， y / z = tan(Ayz), 所以最后可得:
x = sin(Axz) / SQRT (1 + cos(Axz)^2 * tan(Ayz)^2 )
替換成原來的符號可得：
RxGyro = sin(Axz(n)) / SQRT (1 + cos(Axz(n))^2 * tan(Ayz(n))^2 )
同理可得：
RyGyro = sin(Ayz(n)) / SQRT (1 + cos(Ayz(n))^2 * tan(Axz(n))^2 )
提示：這個公式還可以更進一步簡化。分式兩邊同除以sin（axz（你））可得：
RxGyro =??1??/ SQRT (1/ sin(Axz(n))^2??+ cos(Axz(n))^2 / sin(Axz(n))^2??* tan(Ayz(n))^2 )
RxGyro =??1??/ SQRT (1/ sin(Axz(n))^2??+ cot(Axz(n))^2??* sin(Ayz(n))^2??/ cos(Ayz(n))^2 )??
現在加減? ?cos(Axz(n))^2/sin(Axz(n))^2? ?= cot(Axz(n))^2?
RxGyro =??1??/ SQRT (1/ sin(Axz(n))^2??-??cos(Axz(n))^2/sin(Axz(n))^2? ?+ cot(Axz(n))^2??* sin(Ayz(n))^2??/ cos(Ayz(n))^2??+ cot(Axz(n))^2 )
綜合條件1、2和3、4可得：
RxGyro =??1??/ SQRT (1??+? ?cot(Axz(n))^2 * sec(Ayz(n))^2 ),? ???其中??cot(x) = 1 / tan(x)??， sec(x) = 1 / cos(x)
這個公式只用了2個三角函數并且計算量更低。如果你有Mathematica程序，通過使用 FullSimplify [Sin[A]^2/ ( 1 + Cos[A]^2 * Tan[B]^2)]你可以驗證這個公式。
現在我們發現：
RzGyro??=??Sign(RzGyro)*SQRT(1 – RxGyro^2 – RyGyro^2).
其中，當 RzGyro>=0時，Sign(RzGyro) = 1 , 當 RzGyro<0時，Sign(RzGyro) = -1 。
一個簡單的估算方法：
Sign(RzGyro) = Sign(RzEst(n-1))
在實際應用中，當心RzEst（n-1）趨近于0。這時候你可以跳過整個陀螺儀階段并賦值：Rgyro=Rest（n-1）。Rz可以用作計算Axz和Ayz傾角的參考，當它趨近于0時，它可能會溢出并引發不好的后果。這時你會得到很大的浮點數據，并且tan()/atan()函數得到的結果會缺乏精度。
現在我們回顧一下已經得到的結果，我們在算法中的第n步，并計算出了下面的值：
Racc – 加速度計讀取的當前值
Rgyro –根據Rest(-1)和當前陀螺儀讀取值所得
我們根據哪個值來更新Rest（n）呢？你可能已經猜到，兩者都采用。我們會用一個加權平均值，得：
Rest(n) = (Racc * w1 + Rgyro * w2 ) / (w1 + w2)
分子分母同除以w1，公式可簡化成：
Rest(n) = (Racc * w1/w1 + Rgyro * w2/w1 ) / (w1/w1 + w2/w1)
令w2=w1=wGyro，可得：
Rest(n) = (Racc + Rgyro * wGyro ) / (1 + wGyro)
在上面的公式中，wGyro表示我們對加速度計和陀螺儀的相信程度。這個值可以通過測試確定，根據經驗值5-20之間會得到一個很好的結果。
此算法和卡爾曼濾波最主要的差別是它的權重是相對固定的，而卡爾曼濾波中的權重會隨著加速度計讀取的噪聲而改變。卡爾曼濾波注重給你一個“最好”的理論結果，而此算法給你的是實際項目中“夠用”的結果。你可以實現一個算法，它能根據測量的噪聲而改變wGyro值，但對大部分應用來說固定的權重也能工作的很好。
現在得到最新的估算值還差一步：
RxEst(n) = (RxAcc + RxGyro * wGyro ) / (1 + wGyro)
RyEst(n) = (RyAcc + RyGyro * wGyro ) / (1 + wGyro)
RzEst(n) = (RzAcc + RzGyro * wGyro ) / (1 + wGyro)
現在，再次標準化矢量：
R = SQRT(RxEst(n) ^2 + RyEst(n)^2 +??RzEst(n)^2 )
RxEst(n) = RxEst(n)/R
RyEst(n) = RyEst(n)/R
RzEst(n) = RzEst(n)/R
這篇文章主要介紹加速度計和陀螺儀的數學模型和基本算法，以及如何融合這兩者，側重算法、思想的討論
介紹
本指南旨在向興趣者介紹慣性MEMS（微機電系統）傳感器，特別是加速度計和陀螺儀以及其他整合IMU（ 慣性測量單元 ）設備。
? 2012-8-22 16:38 上傳 下載附件?(16.48 KB)
IMU單元例子 ：上圖中MCU頂端的 ACC Gyro 6DOF ，名為 USBThumb ，支持USB/串口通信
在這篇文章中我將概括這么幾個基本并且重要的話題：
- 加速度計（accelerometer）檢測什么
- 陀螺儀（gyroscope，也稱作 gyro）檢測什么
- 如何將傳感器ADC讀取的數據轉換為物理單位（加速度傳感器的單位是g，陀螺儀的單位是 度/秒）
- 如何結合加速度傳感器和陀螺儀的數據以得到設備和地平面之間的傾角的準確信息
在整篇文章中我盡量將數學運算降低到最少。如果你知道什么是正弦、余弦、正切函數，那無論你的項目使用哪種平臺你應該都會明白和運用這篇文章中的思想，這些平臺如Arduino、Propeller、Basic Stamp、Ateml芯片、PIC芯片等等。總有些人認為使用IMU單元需要復雜的數學運算（復雜的FIR或IIR濾波，如卡爾曼濾波，Parks-McClellan濾波等）。你如果研究這些會得到很棒且很復雜的結果。我解釋事情的方式，只需要基本的數學。我非常堅信簡單的原則。我認為一個簡單的系統更容易操作和監控，另外許多嵌入式設備并不具備能力和資源去實現需要進行矩陣運算的復雜算法。
我會用我設計的一個新IMU模塊—— Acc_Gyro Accelerometer + Gyro IMU 作為例子。在下面的例子中我們會使用這個設備的參數。用這個模塊作為介紹非常合適，因為它由3個設備組成：
- LIS331AL ( datasheet ) – 3軸 2G 模擬加速度計
- LPR550AL ( datasheet ) – 雙軸（俯仰、翻滾） 500°/s 加速度傳感器
- LY550ALH ( datasheet ) –單軸（偏航）陀螺儀 最后這個設備在這篇介紹中不使用，不過他在DCM Matrix implementation 中有重要作用
它們一起組成了一個6自由度的慣性測量單元。這是個花哨的名字！然而，在花哨的名字后面是個非常有用的設備組合，接下來我們會詳細介紹之。
第一部分 加速度計
要了解這個模塊我們先從加速度計開始。當我們在想象一個加速度計的時候我們可以把它想作一個圓球在一個方盒子中。你可能會把它想作一個餅干或者甜圈，但我就把它當做一個球好了：
? 2012-8-22 16:38 上傳 下載附件?(15.85 KB)
我們假定這個盒子不在重力場中或者其他任何會影響球的位置的場中，球處于盒子的正中央。你可以想象盒子在外太空中，遠離任何天體，如果很難想象，那就當做盒子在航天飛機中，一切東西都處于無重力狀態。在上面的圖中你可以看到我們給每個軸分配了一對墻（我們移除了Y+以此來觀察里面的情況）。設想每面墻都能感測壓力。如果我們突然把盒子向左移動（加速度為1g=9.8m/s^2），那么球會撞上X-墻。然后我們檢測球撞擊墻面產生的壓力，X軸輸出值為-1g。
? 2012-8-22 16:38 上傳 下載附件?(19.23 KB)
請注意加速度計檢測到得力的方向與它本身加速度的方向是相反的。這種力量通常被稱為 慣性力或假想力 ?。在這個模型中你你應該學到加速度計是通過間接測量力對一個墻面的作用來測量加速度的，在實際應用中，可能通過彈簧等裝置來測量力。這個力可以是加速度引起的，但在下面的例子中，我們會發現它不一定是加速度引起的。
如果我們把模型放在地球上，球會落在Z-墻面上并對其施加一個1g的力，見下圖：
? 2012-8-22 16:38 上傳 下載附件?(19.54 KB)
在這種情況下盒子沒有移動但我們任然讀取到Z軸有-1g的值。球在墻壁上施加的壓力是由引力造成的。在理論上，它可以是不同類型的力量 - 例如，你可以想象我們的球是鐵質的，將一個磁鐵放在盒子旁邊那球就會撞上另一面墻。引用這個例子只是為了說明加速度計的本質是檢測力而非加速度。只是加速度所引起的慣性力正好能被加速度計的檢測裝置所捕獲。
雖然這個模型并非一個MEMS傳感器的真實構造，但它用來解決與加速度計相關的問題相當有效。實際上有些類似傳感器中有金屬小球，它們稱作傾角開關，但是它們的功能更弱，只能檢測設備是否在一定程度內傾斜，卻不能得到傾斜的程度。
到目前為止，我們已經分析了單軸的加速度計輸出，這是使用單軸加速度計所能得到的。三軸加速度計的真正價值在于它們能夠檢測全部三個軸的慣性力。讓我們回到盒子模型，并將盒子向右旋轉45度。現在球會與兩個面接觸：Z-和X-，見下圖：
? 2012-8-22 16:38 上傳 下載附件?(34.03 KB)
0.71g這個值是不是任意的，它們實際上是1/2的平方根的近似值。我們介紹加速度計的下一個模型時這一點會更清楚。
在上一個模型中我們引入了重力并旋轉了盒子。在最后的兩個例子中我們分析了盒子在兩種情況下的輸出值，力矢量保持不變。雖然這有助于理解加速度計是怎么和外部力相互作用的，但如果我們將坐標系換為加速度的三個軸并想象矢量力在周圍旋轉，這會更方便計算。
? 2012-8-22 16:38 上傳 下載附件?(7.19 KB)
請看看在上面的模型，我保留了軸的顏色，以便你的思維能更好的從上一個模型轉到新的模型中。想象新模型中每個軸都分別垂直于原模型中各自的墻面。矢量R是加速度計所檢測的矢量（它可能是重力或上面例子中慣性力的合成）。RX，RY，RZ是矢量R在X，Y，Z上的投影。請注意下列關系：
，R ^ 2 = RX ^ 2 + RY ^ 2 + RZ ^ 2（ 公式1 ）
此公式等價于 三維空間勾股定理 。
還記得我剛才說的1/2的平方根0.71不是個隨機值吧。如果你把它們代回上式，回顧一下重力加速度是1g，那我們就能驗證：
1 ^ 2 =（SQRT（1/2））^ 2 + 0 ^ 2 +（SQRT（1/2））^ 2
在 公式1 中簡單的取代： R=1, Rx = -SQRT(1/2), Ry = 0 , Rz = -SQRT(1/2)
經過一大段的理論序言后，我們和實際的加速度計很靠近了。RX，RY，RZ值是實際中加速度計輸出的線性相關值，你可以用它們進行各種計算。
在我們運用它之前我們先討論一點獲取加速度計數據的方法。大多數加速度計可歸為兩類：數字和模擬。數字加速度計可通過I2C，SPI或USART方式獲取信息，而模擬加速度計的輸出是一個在預定范圍內的電壓值，你需要用ADC（模擬量轉數字量）模塊將其轉換為數字值。我將不會詳細介紹ADC是怎么工作的，部分原因是這是個很廣的話題，另一個原因是不同平臺的ADC都會有差別。有些MCU具有內置ADC模塊，而有些則需要外部電路進行ADC轉換。不管使用什么類型的ADC模塊，你都會得到一個在一定范圍內的數值。例如一個10位ADC模塊的輸出值范圍在0 .. 1023間，請注意，1023 = 2 ^ 10 -1。一個12位ADC模塊的輸出值范圍在0 .. 4095內，注意，4095 = 2 ^ 12-1。
我們繼續，先考慮下一個簡單的例子，假設我們從10位ADC模塊得到了以下的三個軸的數據：
AdcRx = 586
AdcRy = 630
AdcRz = 561
每個ADC模塊都有一個參考電壓，假設在我們的例子中，它是3.3V。要將一個10位的ADC值轉成電壓值，我們使用下列公式：
VoltsRx = AdcRx * VREF / 1023
小注：8位ADC的最大值是255 = 2 ^ 8 -1，12位ADC最大值是4095 = 2 ^ 12 -1。
將3個軸的值代入上式，得到：
VoltsRx = 586 * 3.3 / 1023 =～1.89V（結果取兩位小數）
VoltsRy = 630 * 3.3 / 1023 =～2.03V
VoltsRz = 561 * 3.3 / 1023 =～1.81V
每個加速度計都有一個零加速度的電壓值，你可以在它的說明書中找到，這個電壓值對應于加速度為0g。通過計算相對0g電壓的偏移量我們可以得到一個有符號的電壓值。比方說，0g電壓值 VzeroG= 1.65V，通過下面的方式可以得到相對0g電壓的偏移量:
DeltaVoltsRx = 1.89V - 1.65V = 0.24V
DeltaVoltsRy = 2.03V - 1.65V = 0.38V
DeltaVoltsRz = 1.81V - 1.65V = 0.16V
現在我們得到了加速度計的電壓值，但它的單位還不是g（9.8m/s^2），最后的轉換，我們還需要引入加速度計的靈敏度（Sensitivity），單位通常是 mV/g。比方說，加速度計的靈敏度 Sensitivity= 478.5mV / g = 0.4785V /g。靈敏度值可以在加速度計說明書中找到。要獲得最后的單位為g的加速度，我們使用下列公式計算：
RX = DeltaVoltsRx /Sensitivity
RX = 0.24V / 0.4785V / G =～0.5g
RY = 0.38V / 0.4785V / G =～0.79g
RZ = 0.16V / 0.4785V / G =～0.33g
當然，我們可以把所有的步驟全部放在一個式子里，但我想通過介紹每一個步驟以便讓你了解怎么讀取一個ADC值并將其轉換為單位為g的矢量力的分量。
Rx = (AdcRx * Vref / 1023 – VzeroG) / Sensitivity ?(公式2)
Ry = (AdcRy * Vref / 1023 – VzeroG) / Sensitivity
Rz = (AdcRz * Vref / 1023 – VzeroG) / Sensitivity
現在我們得到了慣性力矢量的三個分量，如果設備除了重力外不受任何外力影響，那我們就可以認為這個方向就是重力矢量的方向。如果你想計算設備相對于地面的傾角，可以計算這個矢量和Z軸之間的夾角。如果你對每個軸的傾角都感興趣，你可以把這個結果分為兩個分量：X軸、Y軸傾角，這可以通過計算重力矢量和X、Y軸的夾角得到。計算這些角度比你想象的簡單，現在我們已經算出了Rx，Ry，Rz的值，讓我們回到我們的上一個加速度模型，再加一些標注上去：
? 2012-8-22 16:38 上傳 下載附件?(9.9 KB)
我們感興趣的角度是向量R和X，Y，Z軸之間的夾角，那就令這些角度為Axr，Ayr，Azr。觀察由R和Rx組成的直角三角形：
cos(Axr) = Rx / R , 類似的：
cos(Ayr) = Ry / R
cos(Azr) = Rz / R
從 公式1 我們可以推導出 R = SQRT( Rx^2 + Ry^2 + Rz^2)
通過arccos()函數（cos()的反函數）我們可以計算出所需的角度：
Axr = arccos(Rx/R)
Ayr = arccos(Ry/R)
Azr = arccos(Rz/R)
我們花了大段的篇幅來解釋加速度計模型，最后所要的只是以上這幾個公式。根據你的應用場合，你可能會用到我們推導出來的幾個過渡公式。我們接下來要介紹陀螺儀模塊，并向大家介紹怎么融合加速度計和陀螺儀的數據以得到更精確的傾角值。
但在此之前，我們再介紹幾個很常用的公式：
cosX = cos(Axr) = Rx / R
cosY = cos(Ayr) = Ry / R
cosZ = cos(Azr) = Rz / R
這三個公式通常稱作 方向余弦 ?，它主要表達了單位向量（長度為1的向量）和R向量具有相同的方向。你可以很容易地驗證：
SQRT（cosX ^ 2 + COSY ^ 2 + cosZ ^ 2）= 1
這是個很好的性質，因為它避免了我們一直檢測R向量的模（長度）。通常如果我們只是對慣性力的方向感興趣，那標準化模長以簡化其他計算是個明智的選擇。
第二部分陀螺儀
對于陀螺儀我們將不會像加速度計一樣介紹它的等價盒子模型，而是直接跳到加速度計的第二個模型，通過這個模型我們會向大家介紹陀螺儀是怎么工作的。
? 2012-8-22 16:38 上傳 下載附件?(12.55 KB)
陀螺儀的每個通道檢測一個軸的旋轉。例如，一個2軸陀螺儀檢測繞X和Y軸的旋轉。為了用數字來表達這些旋轉，我們先引進一些符號。首先我們定義：
Rxz – 慣性力矢量R在XZ平面上的投影
Ryz – 慣性力矢量R在YZ平面的上投影
在由Rxz和Rz組成的直角三角形中，運用勾股定理可得：
Rxz^2 = Rx^2 + Rz^2 ，同樣：
Ryz^2 = Ry^2 + Rz^2
同時注意：
R^2 = Rxz^2 + Ry^2 ，這個公式可以 公式1 和上面的公式推導出來，也可由R和Ryz所組成的直角三角形推導出來
R ^ 2 = Ryz ^ 2 + RX ^ 2
在這篇文章中我們不會用到這些公式，但知道模型中的那些數值間的關系有助于理解。
相反，我們按如下方法定義Z軸和Rxz、Ryz向量所成的夾角：
AXZ - Rxz（矢量R在XZ平面的投影）和Z軸所成的夾角
AYZ - Ryz（矢量R在YZ平面的投影）和Z軸所成夾角
現在我們離陀螺儀要測量的東西又近了一步。陀螺儀測量上面定義的角度的變化率。換句話說，它會輸出一個與上面這些角度變化率線性相關的值。為了解釋這一點，我們先假設在t0時刻，我們已測得繞Y軸旋轉的角度（也就是Axz），定義為Axz0，之后在t1時刻我們再次測量這個角度，得到Axz1。角度變化率按下面方法計算：
RateAxz = (Axz1 – Axz0) / (t1 – t0).
如果用度來表示角度，秒來表示時間，那這個值的單位就是 度/秒。這就是陀螺儀檢測的東西。
在實際運用中，陀螺儀一般都不會直接給你一個單位為度/秒的值（除非它是個特殊的數字陀螺儀）。就像加速度計一樣，你會得到一個ADC值并且要用類似 公式2 的式子將其轉換成單位為 度/秒的值。讓我們來介紹陀螺儀輸出值轉換中的ADC部分（假設使用10位ADC模塊，如果是8位ADC，用1023代替255，如果是12為ADC用4095代替1023）。
RateAxz = (AdcGyroXZ * Vref / 1023 – VzeroRate) / Sensitivity? 公式3
RateAyz = (AdcGyroYZ * Vref / 1023 – VzeroRate) / Sensitivity
AdcGyroXZ，AdcGyroYZ - 這兩個值由ADC讀取，它們分別代表矢量R的投影在XZ和YZ平面內里的轉角，也可等價的說，旋轉可分解為單獨繞Y和X軸的運動。
Vref – ADC的參考電壓，上例中我們使用3.3V
VzeroRate – 是零變化率電壓，換句話說它是陀螺儀不受任何轉動影響時的輸出值，對 Acc Gyro 板來說，可以認為是1.23V（此值通常可以在說明書中找到——但千萬別相信這個值，因為大多數的陀螺儀在焊接后會有一定的偏差，所以可以使用電壓計測量每個通道的輸出值，通常這個值在焊接后就不會改變，如果有跳動，在設備使用前寫一個校準程序對其進行測量，用戶應當在設備啟動的時候保持設備靜止以進行校準）。
Sensitivity –陀螺儀的靈敏度，單位mV/(deg/s)，通常寫作mV/deg/s,它的意思就是如果旋轉速度增加1°/s,陀螺儀的輸出就會增加多少mV。 Acc_Gyro 板的靈敏度值是2mV/deg/s或0.002V/deg/s
讓我們舉個例子，假設我們的ADC模塊返回以下值：
AdcGyroXZ = 571
AdcGyroXZ = 323
用上面的公式，在代入 Acc Gyro 板的參數，可得：
RateAxz = (571 * 3.3V / 1023 – 1.23V) / ( 0.002V/deg/s) =~ 306 deg/s
RateAyz = (323 * 3.3V / 1023 – 1.23V) / ( 0.002V/deg/s) =~ -94 deg/s
換句話說設備繞Y軸（也可以說在XZ平面內）以306°/s速度和繞X軸（或者說YZ平面內）以-94°/s的速度旋轉。請注意，負號表示該設備朝著反方向旋轉。按照慣例，一個方向的旋轉是正值。一份好的陀螺儀說明書會告訴你哪個方向是正的，否則你就要自己測試出哪個旋轉方向會使得輸出腳電壓增加。最好使用示波器進行測試，因為一旦你停止了旋轉，電壓就會掉回零速率水平。如果你使用的是萬用表，你得保持一定的旋轉速度幾秒鐘并同時比較電壓值和零速率電壓值。如果值大于零速率電壓值那說明這個旋轉方向是正向。
第三部分 將它們綜合起來。融合加速度計和陀螺儀的數據。
如果你在閱讀這篇文章你可能已經有了或準備購買一個IMU設備，或者你準備用獨立的加速度計和陀螺儀搭建一個。
注：具體的代碼實現和算法測試，請閱讀這篇文章：
http://starlino.com/imu_kalman_arduino.html
在使用整合了加速度計和陀螺儀的IMU設備時，首先要做的就是統一它們的坐標系。最簡單的辦法就是將加速度計作為參考坐標系。大多數的加速度計技術說明書都會指出對應于物理芯片或設備的XZY軸方向。例如，下面就是Acc Gyro板的說明書中給出的XYZ軸方向：
? 2012-8-22 16:39 上傳 下載附件?(35.18 KB)
接下來的步驟是：
- 確定陀螺儀的輸出對應到上述討論的RateAxz，RateAyz值。
- 根據陀螺儀和加速度計的位置決定是否要反轉輸出值
不要設想陀螺儀陀的輸出有XY，它會適應加速度計坐標系里的任何軸，盡管這個輸出是IMU模塊的一部分。最好的辦法就是測試。
接下來的示例用來確定哪個陀螺儀的輸出對應RateAxz。
- 首先將設備保持水平。加速度計的XY軸輸出會是零加速度電壓（ Acc Gyro 板的值是1.65V）
- 接下來將設備繞Y軸旋轉，換句話說就是將設備在XZ平面內旋轉，所以X、Z的加速度輸出值會變化而Y軸保持不變。
-當以勻速旋轉設備的時候，注意陀螺儀的哪個通道輸出值變化了，其他輸出應該保持不變。
- 在陀螺儀繞Y軸旋轉（在XZ平面內旋轉）的時候輸出值變化的就是AdcGyroXZ，用于計算RateAxz
-最后一步，確認旋轉的方向是否和我們的模型對應，因為陀螺儀和加速度的位置關系，有時候你可能要把RateAxz值反向
-重復上面的測試，將設備繞Y軸旋轉，這次查看加速度計的X軸輸出（也就是AdcRx）。如果AdcRx增大（從水平位置開始旋轉的第一個90°），那AdcGyroXZ應當減小。這是因為我們觀察的是重力矢量，當設備朝一個方向旋轉時矢量會朝相反的方向旋轉（相對坐標系運動）。所以，如果你不想反轉RateAxz，你可以在 公式3 中引入正負號來解決這個問題：
RateAxz = InvertAxz * (AdcGyroXZ * Vref / 1023 – VzeroRate) / Sensitivity ,其中InvertAxz= 1 或-1
同樣的方法可以用來測試RateAyz，將設備繞X軸旋轉，你就能測出陀螺儀的哪個輸出對應于RateAyz，以及它是否需要反轉。一旦你確定了InvertAyz，你就能可以用下面的公式來計算RateAyz：
RateAyz = InvertAyz * (AdcGyroYZ * Vref / 1023 – VzeroRate) / Sensitivity
如果對Acc Gyro板進行這些測試，你會得到下面的這些結果：
- RateAxz的輸出管腳是GX4，InvertAxz = 1
- RateAyz輸出管腳是GY4，InvertAyz = 1
從現在開始我們認為你已經設置好了IMU模塊并能計算出正確的Axr，Ayr，Azr值（在第一部分加速度計中定義）以及RateAyz，RateAyz（在第二部分陀螺儀中）。下一步，我們分析這些值之間的關系并得到更準確的設備和地平面之間的傾角。
你可能會問自己一個問題，如果加速度計已經告訴我們Axr，Ayr，Azr的傾角，為什么還要費事去得到陀螺儀的數據？答案很簡單：加速度計的數據不是100%準確的。有幾個原因，還記加速度計測量的是慣性力，這個力可以由重力引起（理想情況只受重力影響），當也可能由設備的加速度（運動）引起。因此，就算加速度計處于一個相對比較平穩的狀態，它對一般的震動和機械噪聲很敏感。這就是為什么大部分的IMU系統都需要陀螺儀來使加速度計的輸出更平滑。但是怎么辦到這點呢？陀螺儀不受噪聲影響嗎？
陀螺儀也會有噪聲，但由于它檢測的是旋轉，因此對線性機械運動沒那么敏感，不過陀螺儀有另外一種問題，比如漂移（當選擇停止的時候電壓不會回到零速率電壓）。然而，通過計算加速度計和陀螺儀的平均值我們能得到一個相對更準確的當前設備的傾角值，這比單獨使用加速度計更好。
接下來的步驟我會介紹一種算法，算法受卡爾曼濾波中的一些思想啟發，但是它更簡單并且更容易在嵌入式設備中實現。在此之前，讓我們先看看我們需要算法計算什么值。所要算的就是重力矢量R=[Rx,Ry,Rz]，它可由其他值推導出來，如Axr，Ayr，Azr或者cosX，cosY，cosZ，由這些值我們能得到設備相對地平面的傾角值，這些關系我們在第一部分已經討論過。有人可能會說-根據第一部分的 公式2 我們不是已經得到Rx，Ry，Rz的值了嗎？是的，但是記住，這些值只是由加速度計數據推導出來的，如果你直接將它們用于你的程序你會得到難以忍受的噪聲。為了避免進一步的混亂，我們重新定義加速度計的測量值：
Racc – 是由加速度計測量到得慣性力矢量，它可分解為下面的分量（在XYZ軸上的投影）：
RxAcc = (AdcRx * Vref / 1023 – VzeroG) / Sensitivity
RyAcc = (AdcRy * Vref / 1023 – VzeroG) / Sensitivity
RzAcc = (AdcRz * Vref / 1023 – VzeroG) / Sensitivity
現在我們得到了一組只來自于加速度計ADC的值。我們把這組數據叫做“vector”，并使用下面的符號：
Racc = [RxAcc,RyAcc,RzAcc]
因為這些Racc的分量可由加速度計數據得到，我們可以把它當做算法的輸入。
請注意Racc測量的是重力，如果你得到的矢量長度約等于1g那么你就是正確的：
|Racc| = SQRT(RxAcc^2 +RyAcc^2 + RzAcc^2),
但是請確定把矢量轉換成下面的矢量非常重要：
Racc(normalized) = [RxAcc/|Racc| , RyAcc/|Racc| , RzAcc/|Racc|].
這可以確保標準化Racc始終是1。
接來下我們引進一個新的向量：
Rest = [RxEst,RyEst,RzEst]
這就是算法的輸出值，它經過陀螺儀數據的修正和基于上一次估算的值。
這是算法所做的事：
-加速度計告訴我們：“你現在的位置是Racc”
我們回答：“謝謝，但讓我確認一下”
-然后根據陀螺儀的數據和上一次的Rest值修正這個值并輸出新的估算值Rest。
-我們認為Rest是當前設備姿態的“最佳值”。
讓我們看看它是怎么實現的。
數列的開始，我們先認為加速度值正確并賦值：
Rest(0) = Racc(0)
Rest和Racc是向量，所以上面的式子可以用3個簡單的式子代替，注意別重復了：
RxEst（0）= RxAcc（0）
RyEst（0）= RyAcc（0）
RzEst（0）= RzAcc（0）
接下來我們在每個等時間間隔T秒做一次測量，得到新的測量值，并定義為Racc（1），Racc（2），Racc（3）等等。同時，在每個時間間隔我們也計算出新的估算值Rest（1）,Rest（2），Rest（3），等等。
假設我們在第n步。我們有兩列已知的值可以用：
Rest(n-1) – 前一個估算值，Rest(0) = Racc(0)
Racc(n) – 當前加速度計測量值
在計算Rest（n）前，我們先引進一個新的值，它可由陀螺儀和前一個估算值得到。
叫做Rgyro，同樣它是個矢量并由3個分量組成：
Rgyro = [RxGyro,RyGyro,RzGyro]
我們分別計算這個矢量的分量，從RxGyro開始。
? 2012-8-22 16:39 上傳 下載附件?(12.55 KB)
首先觀察陀螺儀模型中下面的關系，根據由Rz和Rxz組成的直角三角形我們能推出：
tan(Axz) = Rx/Rz => Axz = atan2(Rx,Rz)
你可能從未用過atan2這個函數，它和atan類似，但atan返回值范圍是（-PI/2,PI/2），atan2返回值范圍是（-PI，PI），并且他有兩個參數。它能將Rx，Rz值轉換成360°（-PI，PI）內的角度。更多信息請閱讀? atan2 .
所以，知道了RxEst（n-1）和RzEst（n-1）我們發現：
Axz(n-1) = atan2( RxEst(n-1) , RzEst(n-1) ).
記住，陀螺儀測量的是Axz角度變化率，因此，我們可以按如下方法估算新的角度Axz（n）：
Axz(n) = Axz(n-1) + RateAxz(n) * T
請記住，RateAxz可由陀螺儀ADC讀取得到。通過使用平均轉速可由得到一個更準確的公式：
RateAxzAvg =（RateAxz（N）+ RateAxz（N-1））/ 2
Axz(n) = Axz(n-1) + RateAxzAvg * T
同理可得：
Ayz(n) = Ayz(n-1) + RateAyz(n) * T
好了，現在我們有了Axz（n），Ayz（n）。現在我們如何推導出RxGyro/RyGyro?根據 公式1 我們可以把Rgyro長度寫成下式：
| Rgyro | = SQRT（RxGyro ^ 2 + RyGyro ^ 2 + RzGyro ^ 2）
同時，因為我們已經將Racc標準化，我們可以認為它的長度是1并且旋轉后保持不變，所以寫成下面的方式相對比較安全：
| Rgyro | = 1
我們暫時采用更短的符號進行下面的計算：
x =RxGyro , y=RyGyro, z=RzGyro
根據上面的關系可得：
x = x / 1 = x / SQRT(x^2+y^2+z^2)
分子分母同除以SQRT（X ^ 2 + Z ^ 2）
x = ( x / SQRT(x^2 + z^2) ) / SQRT( (x^2 + y^2 + z^2) / (x^2 + z^2) )
注意x / SQRT(x^2 + z^2) = sin(Axz), 所以:
x = sin(Axz) / SQRT (1 + y^2 / (x^2 + z^2) )
將SQRT內部分式的分子分母同乘以z^2
x = sin(Axz) / SQRT (1 + y^2??* z ^2 / (z^2 * (x^2 + z^2)) )
注意 z / SQRT(x^2 + z^2) = cos(Axz)， y / z = tan(Ayz), 所以最后可得:
x = sin(Axz) / SQRT (1 + cos(Axz)^2 * tan(Ayz)^2 )
替換成原來的符號可得：
RxGyro = sin(Axz(n)) / SQRT (1 + cos(Axz(n))^2 * tan(Ayz(n))^2 )
同理可得：
RyGyro = sin(Ayz(n)) / SQRT (1 + cos(Ayz(n))^2 * tan(Axz(n))^2 )
提示：這個公式還可以更進一步簡化。分式兩邊同除以sin（axz（你））可得：
RxGyro =??1??/ SQRT (1/ sin(Axz(n))^2??+ cos(Axz(n))^2 / sin(Axz(n))^2??* tan(Ayz(n))^2 )
RxGyro =??1??/ SQRT (1/ sin(Axz(n))^2??+ cot(Axz(n))^2??* sin(Ayz(n))^2??/ cos(Ayz(n))^2 )??
現在加減? ?cos(Axz(n))^2/sin(Axz(n))^2? ?= cot(Axz(n))^2?
RxGyro =??1??/ SQRT (1/ sin(Axz(n))^2??-??cos(Axz(n))^2/sin(Axz(n))^2? ?+ cot(Axz(n))^2??* sin(Ayz(n))^2??/ cos(Ayz(n))^2??+ cot(Axz(n))^2 )
綜合條件1、2和3、4可得：
RxGyro =??1??/ SQRT (1??+? ?cot(Axz(n))^2 * sec(Ayz(n))^2 ),? ???其中??cot(x) = 1 / tan(x)??， sec(x) = 1 / cos(x)
這個公式只用了2個三角函數并且計算量更低。如果你有Mathematica程序，通過使用 FullSimplify [Sin[A]^2/ ( 1 + Cos[A]^2 * Tan[B]^2)]你可以驗證這個公式。
現在我們發現：
RzGyro??=??Sign(RzGyro)*SQRT(1 – RxGyro^2 – RyGyro^2).
其中，當 RzGyro>=0時，Sign(RzGyro) = 1 , 當 RzGyro<0時，Sign(RzGyro) = -1 。
一個簡單的估算方法：
Sign(RzGyro) = Sign(RzEst(n-1))
在實際應用中，當心RzEst（n-1）趨近于0。這時候你可以跳過整個陀螺儀階段并賦值：Rgyro=Rest（n-1）。Rz可以用作計算Axz和Ayz傾角的參考，當它趨近于0時，它可能會溢出并引發不好的后果。這時你會得到很大的浮點數據，并且tan()/atan()函數得到的結果會缺乏精度。
現在我們回顧一下已經得到的結果，我們在算法中的第n步，并計算出了下面的值：
Racc – 加速度計讀取的當前值
Rgyro –根據Rest(-1)和當前陀螺儀讀取值所得
我們根據哪個值來更新Rest（n）呢？你可能已經猜到，兩者都采用。我們會用一個加權平均值，得：
Rest(n) = (Racc * w1 + Rgyro * w2 ) / (w1 + w2)
分子分母同除以w1，公式可簡化成：
Rest(n) = (Racc * w1/w1 + Rgyro * w2/w1 ) / (w1/w1 + w2/w1)
令w2=w1=wGyro，可得：
Rest(n) = (Racc + Rgyro * wGyro ) / (1 + wGyro)
在上面的公式中，wGyro表示我們對加速度計和陀螺儀的相信程度。這個值可以通過測試確定，根據經驗值5-20之間會得到一個很好的結果。
此算法和卡爾曼濾波最主要的差別是它的權重是相對固定的，而卡爾曼濾波中的權重會隨著加速度計讀取的噪聲而改變。卡爾曼濾波注重給你一個“最好”的理論結果，而此算法給你的是實際項目中“夠用”的結果。你可以實現一個算法，它能根據測量的噪聲而改變wGyro值，但對大部分應用來說固定的權重也能工作的很好。
現在得到最新的估算值還差一步：
RxEst(n) = (RxAcc + RxGyro * wGyro ) / (1 + wGyro)
RyEst(n) = (RyAcc + RyGyro * wGyro ) / (1 + wGyro)
RzEst(n) = (RzAcc + RzGyro * wGyro ) / (1 + wGyro)
現在，再次標準化矢量：
R = SQRT(RxEst(n) ^2 + RyEst(n)^2 +??RzEst(n)^2 )
RxEst(n) = RxEst(n)/R
RyEst(n) = RyEst(n)/R
RzEst(n) = RzEst(n)/R
現在，可以再次進行下一輪循環了。
注：關于此算法的具體實現和測試，請閱讀這篇文章：
http://starlino.com/imu_kalman_arduino.html
加速度計和陀螺儀IMU融合的其他資源：
http://www.mikroquad.com/pub/Res ... ryFilter/filter.pdf
http://stackoverflow.com/questio ... -accelerometer-data
http://www.dimensionengineering.com/accelerometers.htm
總結
以上是生活随笔為你收集整理的加速度计和陀螺仪数据融合的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Docker容器解决没有Vim命令
- 下一篇： 【Go】panic: reflect: