1、卡爾曼濾波算法及 C 語言實現摘要：本文著重討論了卡爾曼濾波器的原理，典型算法以及應用領域。清晰地闡述了kalman filter 在信息估計方面的最優性能。著重介紹簡單 kalman filter algorithm 的編程，使用 kalman filter 的經典 5filter 的最優估計能力。現最優化遞歸公式來編程。通過 c 語言編寫程序實現 kalman：kalman filter；最優估計；C 語言1 引言Kalman Filter 是一個高效的遞歸濾波器，它可以實現從一系列的噪聲測量中，估計動態于 Rudolf Emil Kalman 在 1960 年的博士和的A New問題的新

2、方法）。并系統的狀態。Approach to Linear Eiltering and Prediction Problems（線性濾波與上應用，此后 kalman filter 取得且最先在阿登月計劃軌跡發展和完善。它的廣泛應用已經超過 30 年，包括人導航，。傳感器數據融合甚至在軍事方面的雷達系統以及追蹤等等，近年來更被廣泛應用于計算機圖像處理，例如頭臉識別，圖像分割，圖像邊緣檢測等等。2 kalman filter 最優化遞歸估計Kalman filter 是一個“optimal recursive data processing algorithm（最優化遞歸數據處理方法）”。對于解決

3、很大部分的問題，他是最優，效率最高甚至是最有用的方法。而 kalman filter 最為 的內容是體現它最優化估計和遞歸特點的 5 條公式。舉一個例子來詳細說明 5 條公式的物理意義。假設我們要研究的對象是某一個房間的溫度信號。對于室溫來說，一分鐘內或一小段時間內的值是基本上不變的或者變化范圍很小。也就是說t1 時刻的溫度T1 和t2 時刻的溫度T2基本不變，即T2 = T1 。在這個過程中，因為畢竟溫度還是有所改變的，設有幾度的偏差。我們把這幾度的偏差看成是高斯白噪聲 w(t) ，也就是說 Ew(t) = 0 ，Dw(t) = s 2 。除此之外我們在用一個溫度計來實時測量房間的溫度值 Z

4、 ，但由于量具本身的誤差，所測得的溫度值也是確的，也會和實際值偏差幾度，把這幾度的偏差看成是測量噪聲v(t) 。即滿足 Ev(t) = 0 ， Dv(t) = s 。21此時我們對于這個房間的溫度就得到了兩個數據。一個是你根據經驗得到的經驗值1T2 = T1 ，一個是從溫度計上得到的測量值 Z ，以及各自引入的高斯白噪聲。下面就具體講解 kalman filter 來估計房間溫度的原理與步驟。要估計 K 時刻的實際溫度值，首先要根據 K-1 時刻的溫度值K 時刻的溫度，按照之前我們討論的 T = T，若 k-1 時刻的溫度值是 T= 230C ，那么 此時的 k -121T = T= 230C

5、 ，假如該值的噪聲是 w(k) =50C ，5°是這樣得到的，若果 k-1 時刻估算kk -1出的最優溫度值的噪聲是 w(k)' =30C ，的噪聲是 v(k)' =40C ，所以總體的噪聲為w(k) =w(k)'2 + v(k)'2 =50C 。此時再從溫度計上得到 K 時刻的溫度值為T = 25oC ，設kz該測量值的噪聲是40C 。現在發現問題了，在 k 時刻我們就有了兩個溫度值T = 230C 和T = 25oC ，要信那個kkz呢，簡單的求平均已經不能滿足精度的要求了。我們可以用他們的協方差 covariance 來。協方差本身就能體現兩個信

6、號的相關性，通過它就能到底真值更逼近于值還是測量值。引入 kalman gain（ kg ），有公式計算kg ，222'2222kg = w(k) /(w(k) + v(k) ) = 5 /(5 + 4 )所以kg =0.78。我們可以估算出 K 時刻的實際溫度值是，（1）T = T + k ´(T -T ) = 23 + 0.78´ (25 - 23) = 24.56oC（2）kgkzk可以看出這個值接近于溫度計測量到的值，所以估算出的最優溫度值偏向溫度計的值。這時我們已經得到了 K 時刻的最優溫度值，接下來估計 K+1 時刻的最優溫度值。既然參數kgkalman

7、 filter 是一個最優化的遞歸處理方法，那么遞歸就體現在該算法的一個上，由公式（1） kg 的算法可知每次計算時的kg 是不一樣的。這樣我們要估計 K+1 時刻的最優溫度值，就得先算出 K 時刻的kg ，然后才能利用公式（2）估計 K+1 時刻的最優溫度值。由此可以看出我們只需知道初始時刻的值和它所對應的協方差以及測量值，就可以進行kalman 估計了。3 Kalman Filter Algorithm首先以一個離散程來描述。過程為例討論 kalman filter algorithm。該系統可用一個線性微分方X (k) = A× X (k -1) + B ×U (k)

8、 +W (k) （3）2Z(k) = H × X (k) +V (k) （4）（3）式和（4）式中， X (k) 是 K 時刻的系統狀態，U (k) 是 K 時刻對系統的量，A 和 B 是系統參數，對于多模型系統，它們為矩陣。 Z (k) 是 K 時刻的測量值，H 是測量系統的參數，對于多測量系統，H 為矩陣。W (k) 和V (k) 分別表示系統和測量過程中的噪聲，使用kalman filter 估計時，我們認為噪聲滿足高斯白噪聲模型，設W (k) 和V (k) 的 covariance分別為 Q 和R。討論kalman filter algorithm 的 5 個經典公式。第一步

9、，現在的狀態：X (k | k -1) = A× X (k -1| k -1) + B ×U (k)（5）式（5）中 X (k | k -1) 是利用上一狀態的結果， X (k -1| k -1) 是上一時刻的最優值，U (k) 為現在狀態的量，如果沒有，可以為 0。經過公式（5）后系統結果已經更新了，對應于 X (k | k -1) 的 covariance 還沒有更新，用 P 表示 covariance，P(k | k -1) = A× P(k -1| k -1) AT + Q（6）式（6）中 P(k | k -1) 是 X (k | k -1) 對應的 co

10、variance，P(k -1| k -1) 是 X (k -1| k -1)對應的 covariance， AT 是 A 的轉置矩陣。Q 是系統的噪聲，（5）和（6）式便是 kalman filter中的前兩個公式。對系統的。有了系統的，接下來就要參考測量值進行估計了。X (k | k) = X (k | k -1) + kg (k) ×(Z(k) - H × X (k | k -1)（7）由上面分析可知為了實現遞歸，每次的kg 都是實時更新的。kg (k) = P(k | K -1) × H /(H × P(k | k -1) × H +

11、R)TT（8）P(k | k) = (1- kg (k) × H ) × P(k | k -1)（9）這樣每次 P(k | k) 和 kg (k) 都需要前一時刻的值來更新，遞歸的估計下去。（5）（9）式便是 kalman filter algorithm 的五條公式。4 利用 C 語言編程實現 Kalman Filter Algorithm3要求是給定一個固定量，然后由測量值來使用 kalman filter 估計系統真實值。為了編程簡單，我將（5）式中的 A=1，U (k) =0，（5）式改寫為下面的形式，X (k | k -1) = X (k -1| k -1)（10

12、）式（6）改寫為，P(k | k -1) = P(k -1| k -1) + Q（11）再令 H=1，式（7），（8），（9）可改寫為，X (k | k) = X (k | k -1) + kg (k) ×(Z(k) - X (k | k -1)（12）kg (k) = P(k | K -1) /(P(k | k -1) + R)（13）P(k | k) = (1- kg (k) × P(k | k -1)算法）。（14）使用C 語言編程實現（x_mid=x_last;/x_last=x(k-1|k-1),x_mid=x(k|k-1)p_mid=p_last+Q;/p_mi

13、d=p(k|k-1),p_last=p(k-1|k-1),Q=噪聲kg=p_mid/(p_mid+R); /kg 為 kalman filter，R 為噪聲z_measure=z_real+frand()*0.03;/測量值x_now=x_mid+kg*(z_measure-x_mid);/估計出的最優值p_now=(1-kg)*p_mid;/最優值對應的 covariancep_last = p_now;/更新 covariance 值x_last = x_now;/更新系統狀態值5 算法測試為了測試kalman filter algorithm，我設計了一個簡單實驗，來驗證 kalman

14、filter 的優良性。程序中給定一個初值，然后給定一組測量值，驗證kalman filter 估值的準確性。根據kalman filter algorithm，我們需要給定系統初始值 x_last，系統噪聲 Q 和測量噪聲R，以及初始值所對應的協方差 P_last。為了驗證優劣性，還需要給定真實值 z_real 來計算kalman filter 誤差error_kalman 以及測量誤差error_measure 以及它們在有限次的計算中的累積誤差，累積kalman 誤差 sumerror_kalman 和累積測量誤差sumerror_measure。實驗中給定 x_last=0，p_las

15、t=0,Q=0.018,R=0.0542.實驗中可以通過適當改變 Q 和 R 來獲得更好的估計結果。也可以改變 p_last 和 x_last 的值，由于kalman filter 是對協方差的遞歸算法來估計信號數據的，所以 p_last 對算法結果的影響很大，圖 3 就說明了這一情況，由于在初始時就有協方差，所以在運行過程中算法累積誤差相比初始時沒有誤差的就比較大。給定值為 z_real=0.56 時運行結果如圖 1 所示：4圖 1 真值為 0.56 的運行結果給定值 z_real=3.32 時的運行結果如圖 2圖 2 真值為 3.32 的運行結果圖 3 為 Q，R 不變，p_last=0.

16、02，x_last=0，z_real=0.56 時的測試結果。通過和前兩次結果比較發現p_last 對估計結果影響較大，分析可知這種現象是符合kalman filter 的，通過改變 Q 和R 的值也能算法的性能，但是實際操作中我們并不能這兩個量。5圖 3 改變 p_last 的測試結果6 結論本文通過對 kalman filter algorithm 的深入探討，對kalman filter 有了更深刻的認識，理解了的 5 條公式的物理意義，以及kalman filter 的思想，并通過理解算法編程實踐，驗證了kalman filter 在數據處理方面的優良性能。并且通過實驗結果分析了kal

17、man filter algorithm 的本質對每次估計產生的協方差遞歸結合當前測量值來估計系統當前的最佳狀態。如要算法的性能就必須要盡可能的減小系統噪聲和測量噪聲，優化程序，減小估計的協方差。參考文獻1譚浩強.C 程序設計(第三版)M.北京:,2005,91130.2崔平遠,黃曉瑞.基于卡爾曼濾波的多傳感器信息融合算法及其應用J.電機與學報,2001,9(5): 204-207.3黨宏社,韓崇昭,段戰勝.基于多卡爾曼濾波器的自適應傳感器融合J.系統工程與電子技術,2004,5(26):311-313.4文貢堅,. 一種穩健的直線提取算法J.軟件學報,2001,11(11):1660-166

18、5.6附錄：源程序#include "stdio.h" #include "stdlib.h" #include "math.h"double frand()return 2*(rand()/(double)RAND_MAX) - 0.5);/隨機噪聲void main()float float float float float float float float float float float float floatint i;x_last=0; p_last=0.02; Q=0.018; R=0.542;kg; x_mid;

19、 x_now; p_mid; p_now;z_real=0.56;/0.56 z_measure; sumerror_kalman=0;sumerror_measure=0;x_last=z_real+frand()*0.03; x_mid=x_last; for(i=0;i<20;i+) x_mid=x_last;/x_last=x(k-1|k-1),x_mid=x(k|k-1) p_mid=p_last+Q; /p_mid=p(k|k-1),p_last=p(k-1|k-1),Q=噪聲kg=p_mid/(p_mid+R); /kg 為kalman filter，R 為噪聲z_measure=z_real+frand()*0.03;/測量值x_now=x_mid+kg*(z_measure-x_mid);/估計出的最優值p_now=(1-kg)*p_mid;/最優值對應的 covariancepri