發布日期:2022-10-09 點擊率:182
展開全文
卡爾曼濾波器是傳感器融合工程師用于自動駕駛汽車的工具。想象一下,你有一個雷達傳感器,告訴你另一輛車距離15米,一個激光傳感器說車輛距離20米。你如何協調這些傳感器測量?這就是卡爾曼濾波器的功能。卡爾曼濾波在自動駕駛汽車上的應用十分廣泛,本文講述卡爾曼濾波算法,希望對你有所幫助。
卡爾曼濾波算法在控制領域有極廣泛的應用,在發動機燃油噴射控制中,可以應用擴展的卡爾曼濾波理論研究瞬態工況下發動機循環進氣量的最優估計算法,在雷達中,人們感興趣的是跟蹤目標,但目標的位置、速度、加速度的測量值往往在任何時候都有噪聲。卡爾曼濾波利用目標的動態信息,設法去掉噪聲的影響,得到一個關于目標位置的好的估計。
為了以后更好的工程實踐應用卡爾曼濾波算法,今天小編帶領著大家了解卡爾曼濾波算法的理論,及其在自動駕駛多傳感器融合算法中的應用。ー 1 ー什么是卡爾曼濾波?
你可以在任何含有不確定信息的動態系統中使用卡爾曼濾波,對系統下一步的走向做出有根據的預測,即使伴隨著各種干擾,卡爾曼濾波總是能指出真實發生的情況。
在連續變化的系統中使用卡爾曼濾波是非常理想的,它具有占用內存小的優點(除了前一個狀態量外,不需要保留其它歷史數據),并且速度很快,很適合應用于實時問題和嵌入式系統。
在Google上找到的大多數關于實現卡爾曼濾波的數學公式看起來有點晦澀難懂,這個狀況有點糟糕。實際上,如果以正確的方式看待它,卡爾曼濾波是非常簡單和容易理解的,下面我將用漂亮的圖片和色彩清晰的闡述它,你只需要懂一些基本的概率和矩陣的知識就可以了。
ー 2 ー我們能用卡爾曼濾波做什么?
用玩具舉例:你開發了一個可以在樹林里到處跑的小機器人,這個機器人需要知道它所在的確切位置才能導航。我們可以說機器人有一個狀態: ,表示位置和速度:
注意這個狀態只是關于這個系統基本屬性的一堆數字,它可以是任何其它的東西。在這個例子中是位置和速度,它也可以是一個容器中液體的總量,汽車發動機的溫度,用戶手指在觸摸板上的位置坐標,或者任何你需要跟蹤的信號。這個機器人帶有GPS,精度大約為10米,還算不錯,但是,它需要將自己的位置精確到10米以內。樹林里有很多溝壑和懸崖,如果機器人走錯了一步,就有可能掉下懸崖,所以只有GPS是不夠的。或許我們知道一些機器人如何運動的信息:例如,機器人知道發送給電機的指令,知道自己是否在朝一個方向移動并且沒有人干預,在下一個狀態,機器人很可能朝著相同的方向移動。當然,機器人對自己的運動是一無所知的:它可能受到風吹的影響,輪子方向偏了一點,或者遇到不平的地面而翻倒。所以,輪子轉過的長度并不能精確表示機器人實際行走的距離,預測也不是很完美。GPS 傳感器告訴了我們一些狀態信息,我們的預測告訴了我們機器人會怎樣運動,但都只是間接的,并且伴隨著一些不確定和不準確性。但是,如果使用所有對我們可用的信息,我們能得到一個比任何依據自身估計更好的結果嗎?回答當然是YES,這就是卡爾曼濾波的用處。ー 3 ー卡爾曼濾波是如何看到你的問題的?
下面我們繼續以只有位置和速度這兩個狀態的簡單例子做解釋。我們并不知道實際的位置和速度,它們之間有很多種可能正確的組合,但其中一些的可能性要大于其它部分: 卡爾曼濾波假設兩個變量(位置和速度,在這個例子中)都是隨機的,并且服從高斯分布。每個變量都有一個均值μ,表示隨機分布的中心(最可能的狀態),以及方差,表示不確定性。在上圖中,位置和速度是不相關的,這意味著由其中一個變量的狀態無法推測出另一個變量可能的值。下面的例子更有趣:位置和速度是相關的,觀測特定位置的可能性取決于當前的速度: 這種情況是有可能發生的,例如,我們基于舊的位置來估計新位置。如果速度過高,我們可能已經移動很遠了。如果緩慢移動,則距離不會很遠。跟蹤這種關系是非常重要的,因為它帶給我們更多的信息:其中一個測量值告訴了我們其它變量可能的值,這就是卡爾曼濾波的目的,盡可能地在包含不確定性的測量數據中提取更多信息! :這種相關性用協方差矩陣來表示,簡而言之,矩陣中的每個元素, 表示第 i 個和第 j 個狀態變量之間的相關度。(你可能已經猜到協方差矩陣是一個對稱矩陣,這意味著可以任意交換 i 和 j)。協方差矩陣通常用來表示,其中的元素則表示為ー 4 ー使用矩陣來描述問題
我們基于高斯分布來建立狀態變量,所以在時刻 k 需要兩個信息:最佳估計
(即均值,其它地方常用 μ 表示),以及協方差矩陣。
(當然,在這里我們只用到了位置和速度,實際上這個狀態可以包含多個變量,代表任何你想表示的信息)。接下來,我們需要根據當前狀態(k-1 時刻)來預測下一狀態(k 時刻)。記住,我們并不知道對下一狀態的所有預測中哪個是“真實”的,但我們的預測函數并不在乎。它對所有的可能性進行預測,并給出新的高斯分布。我們可以用矩陣來表示這個預測過程:
它將我們原始估計中的每個點都移動到了一個新的預測位置,如果原始估計是正確的話,這個新的預測位置就是系統下一步會移動到的位置。那我們又如何用矩陣來預測下一個時刻的位置和速度呢?下面用一個基本的運動學公式來表示:現在,我們有了一個預測矩陣來表示下一時刻的狀態,但是,我們仍然不知道怎么更新協方差矩陣。此時,我們需要引入另一個公式,如果我們將分布中的每個點都乘以矩陣 A,它的協方差矩陣系會怎樣變化呢?很簡單,下面給出公式:結合方程(4)和(3)得到: ー 5 ー外部控制量我們并沒有捕捉到一切信息,可能存在外部因素會對系統進行控制,帶來一些與系統自身狀態沒有相關性的改變。 以火車的運動狀態模型為例,火車司機可能會操縱油門,讓火車加速。相同地,在我們機器人這個例子中,導航軟件可能會發出一個指令讓輪子轉向或者停止。如果知這些額外的信息,我們可以用一個向量來表示,將它加到我們的預測方程中做修正。假設由于油門的設置或控制命令,我們知道了期望的加速度,根據基本的運動學方程可以得到:
以矩陣的形式表示就是: 稱為控制矩陣,稱為控制向量(對于沒有外部控制的簡單系統來說,這部分可以忽略)。讓我們再思考一下,如果我們的預測并不是100%準確的,該怎么辦呢?ー 6 ー外部干擾如果這些狀態量是基于系統自身的屬性或者已知的外部控制作用來變化的,則不會出現什么問題。但是,如果存在未知的干擾呢?例如,假設我們跟蹤一個四旋翼飛行器,它可能會受到風的干擾,如果我們跟蹤一個輪式機器人,輪子可能會打滑,或者路面上的小坡會讓它減速。這樣的話我們就不能繼續對這些狀態進行跟蹤,如果沒有把這些外部干擾考慮在內,我們的預測就會出現偏差。 在每次預測之后,我們可以添加一些新的不確定性來建立這種與“外界”(即我們沒有跟蹤的干擾)之間的不確定性模型: 原始估計中的每個狀態變量更新到新的狀態后,仍然服從高斯分布。我們可以說的每個狀態變量移動到了一個新的服從高斯分布的區域,協方差為
。換句話說就是,我們將這些沒有被跟蹤的干擾當作是協方差的的噪聲來處理。這產生了具有不同協方差(但是具有相同的均值)的新的高斯分布。我們通過簡單地添加得到的擴展的協方差,下面給出預測步驟的完整表達式:由上式可知,新的最優估計是根據上一個最優估計預測的道德,并加上已知外部控制量的修正。而新的不確定性由上一不確定性預測得到,并加上外部環境的干擾。好了,我們對系統可能的動向有了一個模糊的估計,用和來表示。如果結合傳感器的數據會怎樣呢?
ー 7 ー用測量值來修正估計值我們可能會有多個傳感器來測量系統當前的狀態,哪個傳感器具體測量的是哪個狀態變量并不重要,也許一個是測量位置,一個是測量速度,每個傳感器間接地告訴了我們一些狀態信息。注意,傳感器讀取的數據的單位和尺度有可能與我們要跟蹤的狀態額單位和尺度不一樣,我們用矩陣來表示傳感器的數據。我們可以計算出傳感器讀數的分布,用之前的表示方法如下式所示: 卡爾曼濾波的一大優點就是能處理傳感器噪聲,換句話說,我們的傳感器或多或少都有點不可靠,并且原始估計中的每個狀態可以和一定范圍內的傳感器讀數對應起來。
從測量到的傳感器數據中,我們大致能猜到系統當前處于什么狀態。但是由于存在不確定性,某些狀態可能比我們得到的讀數更接近真實狀態。
我們將這種不確定性(例如:傳感器噪聲)用協方差表示,該分布的均值就是我們讀取到的傳感器數據,稱之為:傳感器噪 現在我們有了兩個高斯分布,一個是在預測值附近,一個是在傳感器讀數附近。我們必須在預測值(粉紅色)和傳感器測量值(綠色)之間找到最優解。 那么,我們最有可能的狀態是什么呢?對于任何可能的度數有兩種情況:(1)傳感器的測量值;(2)由前一狀態得到的預測值。如果我們想知道這兩種情況都可能發生的概率,將這兩個高斯分布相乘就可以了。剩下的就是重疊部分了,這個重疊部分的均值就是兩個估計最可能的值,也就是給定的所有信息中的最優估計。瞧!這個重疊的區域看起來像另一個高斯分布。如你所見,把兩個具有不同均值和方差的高斯分布相乘,你會得到一個新的具有獨立均值和方差的高斯分布!下面用公式講解。ー 8 ー融合高斯分布先以一維高斯分布來分析比較簡單點,具有方差 和 μ 的高斯曲線可以用下式表示:
如果把兩個服從高斯分布的函數相乘會得到什么呢?
將式(9)代入到式(10)中(注意重新歸一化,使總概率為1)可以得到: 將式(11)中的兩個式子相同的部分用 k 表示:
下面進一步將式(12)和(13)寫成矩陣的形式,如果 Σ 表示高斯分布的協方差, 表示每個維度的均值,則:
矩陣稱為卡爾曼增益,下面將會用到。放松!我們快要完成了!
將所有公式結合起來我們有兩個高斯分布,預測部分,和測量部分,將它們放到式(15)中算出它們之間的重疊部分:
由式(14)可得卡爾曼增益為: 將式(16)和式(17)的兩邊同時左乘矩陣的逆(注意里面包含了 )將其約掉,再將式(16)的第二個等式兩邊同時右乘矩陣 的逆得到以下等式: 上式給出了完整的更新步驟方程。就是新的最優估計,我們可以將它和放到下一個預測和更新方程中不斷迭代。 以上所有公式中,你只需要用到式(7)、(18)、(19)。(如果忘了的話,你可以根據式(4)和(15)重新推導一下)
我們可以用這些公式對任何線性系統建立精確的模型,對于非線性系統來說,我們使用擴展卡爾曼濾波,區別在于EKF多了一個把預測和測量部分進行線性化的過程。ー 9 ー自動駕駛中的傳感器融合算法
追蹤靜止和移動的目標是自動駕駛技術領域最為需要的核心技術之一。來源于多種傳感器的信號,包括攝像頭,雷達,以及激光雷達(基于脈沖激光的測距設備)等傳感器組合的組合體來估計位置,速度,軌跡以及目標的種類,例如其他車輛和行人。
你可能會問——為什么我們需要這么多的傳感器?
這是因為每種傳感器提供了追蹤物體所需要的不同精度和類型的信息,尤其是在不同天氣條件下。比如,以激光雷達為基礎的傳感器能很好地解決位置的問題,但是在糟糕的天氣條件下其精度和性能都會有很大程度的下降。另一方面,雷達的空間解決方案也相對便宜很多,與此同時,在惡劣的天氣條件下也能提供更準確的數據。
同樣的,不像激光雷達傳感器,雷達可以提供目標的速度和方位。雷達數據也是計算密集型的,因為一束激光發射非常多包含每個獨立的激光點的范圍的數據,它使得你必須理解你的算法。
組合來自不同傳感器信息的技術稱之為傳感器融合技術。之所以較早的討論這個,是因為應用在傳感器融合之上的算法必須處理短暫的,充滿雜訊的輸入,生成可靠的運動狀態估計的概率。
下面我們將展示在位置的追蹤和估計中最通用的算法,卡爾曼濾波器的變種——‘擴展卡爾曼濾波器’。
1.使用激光雷達數據的基礎卡爾曼濾波器:
卡爾曼濾波器的歷史已經超過半個世紀,但是對于輸入數據的噪聲信息和狀態估計的平滑來說仍然是最有效的傳感器融合算法之一。它假定位置參數是符合高斯分布的,即完全可以被均值和協方差參數化:X~N(μ, σ2)
當傳感器的信息流開始的時候,卡爾曼濾波器使用一系列的狀態信息來預測和衡量更新步驟去更新被追蹤目標的信心值(概率)。預測和更細心的步驟如下圖所示:
狀態預測:
我們會使用一個簡化的線性狀態空間模型(鏈接:去闡述濾波器的工作方式. 一個系統在t時刻的線性狀態可以從t-1時刻根據以下等式被估計:
衡量更新:
卡爾曼濾波器的下一部分則是去使用實測參數z去更新預測狀態'x',通過縮放因子(通常稱之為卡爾曼增益)成比例的計算估計值和測量值之間的誤差。
你可以在以下鏈接中找到衡量更新等式的求導:
舉個例子
理論部分結束!讓我們嘗試用一些代碼去表示一些卡爾曼濾波器的基礎過程。
對此,我們模擬一個目標,它的狀態被四維向量x=[px py vx vy]所描述。
在這個例子中測量傳感器是返回位置數據但是沒有速度信息的激光傳感器。為了觀測到速度信息我們需要使用雷達傳感器數據。在接下來的章節中,當我們討論擴展卡爾曼濾波的時候我們將會涉及到這些。
讓我們以一些假設開始:
偽代碼:
基礎版本的卡爾曼濾波器代碼步驟列在了下面。你可以找到一個基礎的例子:
最后的迭代步驟通過測量和應用預測以及濾波器的更新步驟如下所示:
plot_position_variance(x,P,edgecolor='r') #plot initial position and covariance in red
for z in measurements:
x,P =predict(x, P)
x,P = update(x, P,z)
plot_position_variance(x,P,edgecolor='b') #plot updates in blue
print(x)
print(P)
卡爾曼濾波器迭代:濾波器在迭代之后向真實值收斂
上方的圖闡述了濾波器在每次迭代中狀態向量的px,py維度和位置的協方差發生了哪些變化。紅圈表示初始過程不確定性。隨著預測和測量更新,我們開始有了較小的誤差(半徑),狀態估計開始越來越精確。
如同你看到的,最終估計的狀態向量x[11.99,2.05]非常接近最終的觀測值,誤差最小值也縮小到了0.05。
2.擴展卡爾曼濾波器——使用雷達數據
雷達數據帶來了一些更困難的挑戰。雷達返回數據值的是基于極坐標系,其由三個部分構成:
- ρ /Range(從原點到此的距離)
- ? / bearing ( ρ 和 x的夾角),
- ρ˙:接近率/距離變化率
由于沒有H矩陣將狀態向量映射到雷達的測量空間,我們需要一個函數h(x)來將狀態空間銀蛇到測量空間以測量更新步驟。這個方法是映射極坐標到笛卡爾坐標的方法的衍生方法,定義如下:
這個映射闡述了一個非線性函數,它將使得卡爾曼濾波器的過程和測量符合高斯分布的假設無效。擴展卡爾曼濾波器使用局部線性模型來逼近非線性模型,然后使用卡爾曼濾波應用到逼近值上。局部線性逼近是通過計算當前狀態估計的一階泰勒展開得出的。一階的逼近也叫雅克比矩陣。我們不會在這里過多介紹雅克比矩陣的推到過程。相關內容在網絡中可以找到優秀的介紹,但是如果你想要直接使用這些東西,你可以在以下的github引用中找到代碼實例:Github link
現實參考
你可以在github存儲庫中找到關于Kalman過濾器的c++編寫代碼:
代碼示例:
來源:牛喀網
本篇文章借鑒了加拿大勞里埃大學的Simon D. Levy 教授博客對擴展卡爾曼濾波的理解思路。謹慎地對其中一些冗余的基礎內容進行刪減,并按照自己的理解進行增添改寫。主要是為了自己以后查閱,原博客是很有意思的交互式頁面,建議有英文閱讀能力的同仁閱讀原博主的文章:Kalman Tutorial - Simon D. Levy
1.一個基本模型1.1 機器人小M
為了逐步解釋擴展卡爾曼濾波的原理,我們先建立一個最基本的機器人小M,它只能在直線上做 方向的一維運動,如圖所示:
1.2 位移狀態預測值
我們用如下規律預測機器人的運動:其當前時刻的位移 是上一時刻位移 的 倍。即位移滿足如下遞推公式:
這個公式很好理解,假如小M的初始位移 ,那么第 時刻的位移就依次為 ,是一條關于時刻的指數曲線。
但實際上,小M的運動總是存在各種誤差,這種誤差來源于機器人運動預測模型本身的不精確,并不是觀測引起的(即使用絕對精確的傳感器進行觀測,這種誤差依然存在)。我們用估計值 來替代上式中的精確值 。另外,運動過程必然存在噪聲,用 表示:
一般來說,我們將運動噪聲視為一個均值為 ,方差為 的高斯分布。因此我們在接下來的討論中忽略 ,將上式寫為:
(1)式說明,如果僅僅依靠運動模型進行預測,那么每一步都會得到含有誤差的估計值。我們用 表示第 時刻的預測誤差。則顯然預測誤差和運動模型一樣,也呈不精確的遞推關系:
可以注意到,預測誤差的遞推式中,系數變成了 ,這是因為我們采用了平方差來表示誤差,其中的邏輯后面會講到。 是過程噪聲的平方差,前面已經說到了。
從上面兩式可以看出,每一步遞推得到的狀態值和預測誤差都是不準確的估計值。如果每一步都只包含以上計算,必然會導致誤差累積。為了修正這種誤差,我們就要引入傳感器測量值。
1.3 傳感器測量值
假設我們安裝了一個傳感器,用以測量小M的精確位移 。我們用 表示 時刻傳感器所測量到的位移。理想情況下,傳感器測到的值就是真實的位移 。但實際上,由于傳感器精度有限,我們得到的測量值總是含有部分噪聲。我們用 表示傳感器噪聲。則傳感器測量值可以表示如下:
改寫上式,我們可以用傳感器測量值和噪聲得到位移狀態:
按道理來說,有這個公式,我們就可以計算狀態準確值 了。但實際上,傳感器的噪聲 是未知的(雖然從廠家樣本中包含噪聲的分布狀況,但無法獲知某一時刻精確的噪聲值)。因此我們只能另辟蹊徑,采用其它方法以避開精確的傳感器噪聲。
1.4 狀態準確值 的另一種計算方法
我們注意到 1.2 中已經得到了一個位移狀態的估計值 。而 1.3 中得到了傳感器的測量值 。我們將二者結合起來:
這個公式的含義是:利用運動模型所得到的估計值 和當前時刻的傳感器測量值 來計算當前時刻的準確狀態值 。其中,我們引入了一個增益系數 來表示本次估計的采信偏向。
增益系數
為了理解公式(3)中增益系數的含義,我們可以將增益系數 在 到 之間移動。可以發現,當增益系數等于 時,當前時刻的狀態準確值 等于估計值 。反之,若增益系數等于 ,狀態準確值 等于傳感器測量值 。因此,增益系數的大小實際上表示了對運動模型估計值和傳感器測量值的采信程度(更相信哪一方)。
如何確定增益系數呢?可以做一個感性的理解:既然增益系數代表了對傳感器測量值的采信程度。那么顯然,如果傳感器噪聲很大,傳感器測量值的可信程度就很小。雖然我們無法獲知傳感器噪聲的精確值,但廠家所提供的產品樣板中通常會包含傳感器噪聲的分布狀況,例如給出了噪聲的方差 。則增益系數可以通過如下公式確定:
分析上式,如果傳感器噪聲方差 越大,則增益系數 越小,代表傳感器測量值越不可信,邏輯上沒有問題。 的含義在1.2中已經說過了,是本次預測誤差的估計值。
我們通過公式(3)已經得到了狀態的準確值 ,遞推公式(1)已經可以成功運行了。但遞推公式(2)中還缺少了預測誤差的準確值 有待計算。
1.5 預測誤差準確值 的計算方法
我們已經知道了預測誤差的估計值 ,利用增益系數 可以計算預測誤差準確值:
分析上式,當增益系數 時,預測誤差的準確值就等于估計值 ,這個邏輯是合理的,因為當增益系數為0時,我們的狀態準確值也是直接采信了狀態估計值,說明我們認為傳感器誤差太大,直接采信了運動模型的預測結果;當增益系數 時,我們的預測誤差準確值將直接歸零,這個邏輯也是合理的,從公式(4)可以看出,只有當噪聲方差為0(傳感器絕對準確)時,增益系數才可能為1,而且此時我們將完全采用傳感器的測量值作為準確值(公式3)。既然傳感器絕對準確,那么本步的預測當然不存在任何誤差。
當然,以上兩類均是不可能出現的極端情況,實際情況總是介于兩者之間(沒有絕對準確的運動模型,也沒有絕對準確的傳感器),否則傳感器融合也就失去了意義。
1.6 總結
我們將上面標出序號的式子重新按照計算順序排列出來:
我們把前兩步叫做預測,后三步叫做更新。所謂預測,就是利用上一步的狀態準確值和預測誤差準確值,結合運動模型得到下一步的狀態估計值和預測誤差估計值;所謂更新,就是融入傳感器信息,計算增益系數,將估計值再次轉換為準確值。
至此,我們已經利用一個基本模型小M理解了卡爾曼濾波的雛形。雖然真正的卡爾曼濾波形式與此有一些區別,但基本的邏輯是一致的。接下來我們逐步完善這一公式。
2 推廣到線性卡爾曼濾波
上面的基本模型和線性卡爾曼濾波的標準模型還有一定的差異,接下來我們分步填補這些差異,得到真正的線性卡爾曼濾波模型。
2.1 加入控制量
首先我們看狀態估計值 的遞推方程:
在這個遞推方程中,下一時刻的位移狀態完全取決于上一時刻的位移狀態(和系數 的乘積)。但實際上,一個機器人的狀態變化往往不會這么簡單,我們通常會對機器人施加一個人為的影響,以改變其運動過程(例如對機器人施加一個外力推動其前進)。這種人為的影響無法直接用狀態量來表達,因此我們需要額外引入控制量 及其增益 ,改寫狀態方程如下:
上式表示,某一時刻的狀態估計值由上一時刻的狀態準確值與控制量共同決定。
2.2 引入觀測系數
在基本模型中,我們假設傳感器的測量值 在絕對精確的時候,就等于狀態準確值 。但實際上,我們的傳感器往往并不會直接測量機器人的狀態量,而需要經過一個變換,我們用乘上一個觀測系數 來表示:
其中 仍然為噪聲精確值。為了理解變換系數 的含義,我們還是拿小M舉例。假如我們的傳感器是一臺放在終點的位移傳感器,其正向與小M的位移正向剛好相反,那么其測得的測量值就是真實位移的相反數,此時 ,如圖所示:
又或者,我們用一個碼盤來測量履帶所轉的圈數,那碼盤的測量值就是真實位移除以履帶的周長,即 ,其中 為履帶周長,如圖所示:
引入觀測系數 以后,整個更新過程的形式修改如下:
我們可能會疑惑其中的 是按照什么邏輯插入的,但只要深入分析就能理解其合理性。還是按照極端情況分析,假設傳感器噪聲 ,那么此時增益系數 ,代入公式第二行中,得到 ,說明由于傳感器噪聲為0,是絕對精確的,因此狀態準確值 直接采信了傳感器測量值 ,且其間存在一個觀測系數 的轉換關系;若將增益系數代入第三行,得到 ,說明由于傳感器絕對精確,且直接采信了傳感器測量值,因此該步驟的預測誤差歸零。以上邏輯全部與基本模型中所推導的一致,說明公式(8)的修改是合理的。
2.3 擴展到多維
現在,我們的公式已經修改為了如下形式:
在真實的機器人中,我們所需要的狀態量往往不止一個。例如對小M,我們可能還需要描述其速度 。于是,我們先將狀態從一個維度擴展到多個維度,用狀態向量 來表示:
為了將公式(9)完整地多維化,我們需要設計各個系數和變量的多維形式。
第一行:狀態向量預測
由于引入了多維狀態向量,因此 和 不再代表簡單的位移遞推關系,而要代表各個狀態分量之間的基于運動模型的計算關系。對小M來說,就是要代表本時刻的位移和速度,是如何通過上一次時刻的位移和速度(以及控制量)計算得到的。顯然,這種多維的計算關系也不能用一個變量來表示,而要將 轉換為矩陣。
乍一看很復雜,好在對小M來說,這種計算關系是比較容易找到的,就是勻加速直線運動的位移和速度計算公式:
利用線性代數知識,我們很容易將上式寫為矩陣形式:
對照上式與公式(9)的第一行,我們很容易看出其間的一一對應關系。因此第一行的多維形式如下:
第二行:誤差預測
1.2 節中提到了,我們用平方差來表示預測誤差。在多維情況下,我們用協方差矩陣來替代方差變量,其中主對角線上的值依次代表每個狀態變量本身的誤差(用方差表示),副對角線代表一個狀態變量的誤差隨著另一個狀態變量的誤差的變化,一般為0。例如我們小M的位移誤差為0.1,速度誤差也為0.1,則協方差矩陣為:
關于更詳細的協方差矩陣定義請自行查找資料。
因此第二行可改寫為:
其中 為狀態變量的協方差矩陣, 為噪聲的協方差矩陣。
為 的矩陣乘法形式。
第三到五行:更新
按照上述邏輯,我們可以把更新過程全部多維化。只需要將各變量改為向量形式,系數改為矩陣形式,常數1改為單位矩陣 ,數的倒數改為矩陣的逆即可。
將公式(9)完全轉化為多維形式,我們就得到了完整的線性卡爾曼濾波公式,如下:
真誠說明
接下來我們要做的是將卡爾曼濾波公式進行非線性化,得到擴展卡爾曼濾波公式,并通過一個C++程序來實現。本來想一篇文章搞定,但知乎編輯器寫太多以后變得十分卡頓,體驗很惡劣,因此我先把線性部分發出來,剩下的盡快補上。
我很努力地按照原博主的文章和自己的理解盡可能地寫清楚,但原博主的文章本身存在一些缺陷,而我自己對卡爾曼濾波的理解也相當有限,所以文章中必然存在一些不足之處。如果發現了,歡迎大家友善地指出~
關注:決策智能與機器學習,深耕AI脫水干貨
作者:vincentqin ?來源:計算機視覺SLAM
轉載請聯系作者
譯者注:這恐怕是全網有關卡爾曼濾波最簡單易懂的解釋,如果你認真的讀完本文,你將對卡爾曼濾波有一個更加清晰的認識,并且可以手推卡爾曼濾波。原文作者使用了漂亮的圖片和顏色來闡明它的原理(讀起來并不會因公式多而感到枯燥),所以請勇敢地讀下去!
背景
關于濾波
首先援引來自知乎大神的解釋。
“一位專業課的教授給我們上課的時候,曾談到:filtering is weighting(濾波即加權)。濾波的作用就是給不同的信號分量不同的權重。最簡單的loss pass filter, 就是直接把低頻的信號給1權重,而給高頻部分0權重。對于更復雜的濾波,比如維納濾波, 則要根據信號的統計知識來設計權重。
從統計信號處理的角度,降噪可以看成濾波的一種。降噪的目的在于突出信號本身而抑制噪聲影響。從這個角度,降噪就是給信號一個高的權重而給噪聲一個低的權重。維納濾波就是一個典型的降噪濾波器。”
關于卡爾曼濾波
Kalman Filter 算法,是一種遞推預測濾波算法,算法中涉及到濾波,也涉及到對下一時刻數據的預測。Kalman Filter 由一系列遞歸數學公式描述。它提供了一種高效可計算的方法來估計過程的狀態,并使估計均方誤差最小。卡爾曼濾波器應用廣泛且功能強大:它可以估計信號的過去和當前狀態,甚至能估計將來的狀態,即使并不知道模型的確切性質。
Kalman Filter 也可以被認為是一種數據融合算法(Data fusion algorithm),已有50多年的歷史,是當今使用最重要和最常見的數據融合算法之一。Kalman Filter 的巨大成功歸功于其小的計算需求,優雅的遞歸屬性以及作為具有高斯誤差統計的一維線性系統的最優估計器的狀態。
Kalman Filter 只能減小均值為0的測量噪聲帶來的影響。只要噪聲期望為0,那么不管方差多大,只要迭代次數足夠多,那效果都很好。反之,噪聲期望不為0,那么估計值就還是與實際值有偏差[3]。
問題描述
上面的描述可能把大家繞暈了,下面來點輕松的。
我們會有一個疑問:卡爾曼濾波到底是如何解決實際問題的呢?
我們以機器人為例介紹卡爾曼濾波的原理,我們的任務是要預測機器人的狀態
,包括位置
與速度
,即可表示為:
某個時刻,我們不知道真實的位置與速度到底是多少,二者存在一個所有可能性的組合,大致如下圖所示。
卡爾曼濾波假設狀態所有的變量都是隨機的且都服從高斯分布,每個變量都有其對應的均值
以及方差
(它代表了不確定性)。
在上圖中,位置和速度是不相關(uncorrelated)的,這意味著某個變量的狀態不會告訴你其他變量的狀態是怎樣的。即,我們雖然知道現在的速度,但無法從現在的速度推測出現在的位置。但實際上并非如此,我們知道速度和位置是有關系的(correlated),這樣一來二者之間的組合關系變成了如下圖所示的情況。
這種情況是很容易發生的,例如,如果速度很快,我們可能會走得更遠,所以我們的位置會更大。如果我們走得很慢,我們就不會走得太遠。
這種狀態變量之間的關系很重要,因為它可以為我們提供更多信息:One measurement tells us something about what the others could be。這就是卡爾曼濾波器的目標,我們希望從不確定的測量中盡可能多地獲取信息!
這種狀態量的相關性可以由協方差矩陣表示。簡而言之,矩陣
的每個元素是第i個狀態變量和第j個狀態變量之間的相關度。(顯然地可以知道協方差矩陣是對稱的,這意味著交換i和j都沒關系)。協方差矩陣通常標記為“?
”,因此我們將它們的元素稱為“
”。
狀態預測
問題的矩陣形式表示
我們把狀態建模成高斯分布(Gaussian blob,由于二維高斯分布長得像一個個小泡泡,之所以長這個樣子,可參考鏈接[2])。我們需要求解/估計在時間
時刻的兩個信息:1. 最優估計
以及它的協方差矩陣
,我們可以寫成下面矩陣形式:
(當然,這里我們僅使用位置和速度,但是請記住狀態可以包含任意數量的變量,并且可以表示所需的任何變量)
接下來,我們需要某種方式來查看當前狀態(
時刻)并預測在時刻
處的狀態。請記住,我們不知道哪個狀態是“真實”狀態,但是這里提到的預測(prediction)并不在乎這些。
我們可以用一個矩陣
來表示這個預測過程:
這個矩陣
將原始估計中的每個點移動到新的預測位置。
那么問題來了,應該如何使用上述矩陣來預測下一時刻的位置和速度呢?為了闡述這個過程,我們使用了一個非常基礎的運動學公式(初中物理中就學過)進行描述:
寫成矩陣形式:
現在我們有了一個預測矩陣或者叫做狀態轉移矩陣,該矩陣可以幫助我們計算下一個時刻的狀態。但我們仍然不知道如何更新狀態的協方差矩陣,其實過程也是很簡單,如果我們將分布中的每個點乘以矩陣
,那么其協方差矩陣
會發生什么?
將公式(3)代入公式(4)我們可以得到:
External influence
不過我們并沒有考慮到所有的影響因素。可能有一些與狀態本身無關的變化——如外界因素可能正在影響系統。
例如,我們用狀態對列車的運動進行建模,如果列車長加大油門,火車就加速。同樣,在我們的機器人示例中,導航系統軟件可能會發出使車輪轉動或停止的命令。如果我們很明確地知道這些因素,我們可以將其放在一起構成一個向量
,我們可以對這個量進行某些“處理”,然后將其添加到我們的預測中對狀態進行更正。
假設我們知道由于油門設置或控制命令而產生的預期加速度
。根據基本運動學原理,我們可以得到下式:
將其寫成矩陣形式:
其中
被稱為控制矩陣,
被稱為控制向量。(注意:對于沒有外部影響的簡單系統,可以忽略這個控制項)。
如果我們的預測并不是100%準確模型,這會發生什么呢?
External uncertainty
如果狀態僅僅依賴其自身的屬性進行演進,那一切都會很好。如果狀態受到外部因素進行演進,我們只要知道這些外部因素是什么,那么一切仍然很好。
但在實際使用中,我們有時不知道的這些外部因素到底是如何被建模的。例如,我們要跟蹤四軸飛行器,它可能會隨風搖晃;如果我們跟蹤的是輪式機器人,則車輪可能會打滑,或者因地面顛簸導致其減速。我們無法跟蹤這些外部因素,如果發生任何這些情況,我們的預測可能會出錯,因為我們并沒有考慮這些因素。
通過在每個預測步驟之后添加一些新的不確定性,我們可以對與“世界”相關的不確定性進行建模(如我們無法跟蹤的事物):
這樣一來,由于新增的不確定性原始估計中的每個狀態都可能遷移到多個狀態。因為我們非常喜歡用高斯分布進行建模,此處也不例外。我們可以說
的每個點都移動到具有協方差
的高斯分布內的某個位置,如下圖所示:
這將產生一個新的高斯分布,其協方差不同(但均值相同):
所以呢,我們在狀態量的協方差中增加額外的協方差
,所以預測階段完整的狀態轉移方程為:
換句話說:新的最佳估計是根據先前的最佳估計做出的預測,再加上對已知外部影響的校正。
新的不確定度是根據先前的不確定度做出的預測,再加上來自環境額外的不確定度。
上述過程描繪了狀態預測過程,那么當我們從傳感器中獲取一些測量數據時會發生什么呢?
狀態更新
利用測量進一步修正狀態
假設我們有幾個傳感器,這些傳感器可以向我們提供有關系統狀態的信息。就目前而言,測量什么量都無關緊要,也許一個讀取位置,另一個讀取速度。每個傳感器都告訴我們有關狀態的一些間接信息(換句話說,傳感器在狀態下運作并產生一組測量讀數)。
請注意,測量的單位可能與狀態量的單位不同。我們使用矩陣
對傳感器的測量進行建模。
所以我們期望傳感器的度數可以被建模成如下形式:
卡爾曼濾波器的偉大之處就在于它能夠處理傳感器噪聲。換句話說,傳感器本身的測量是不準確的,且原始估計中的每個狀態都可能導致一定范圍的傳感器讀數,而卡爾曼濾波能夠在這些不確定性存在的情況下找到最優的狀態。
根據傳感器的讀數,我們會猜測系統正處于某個特定狀態。但是由于不確定性的存在,某些狀態比其他狀態更可能產生我們看到的讀數:
我們將這種不確定性(如傳感器噪聲)的協方差表示為
,讀數的分布均值等于我們觀察到傳感器的讀數,我們將其表示為
這樣一來,我們有了兩個高斯分布:一個圍繞通過狀態轉移預測的平均值,另一個圍繞實際傳感器讀數。
因此,我們需要將基于預測狀態(粉紅色)的推測讀數與基于實際觀察到的傳感器讀數(綠色)進行融合。
那么融合后最有可能的新狀態是什么?對于任何可能的讀數
,我們都有兩個相關的概率:(1)我們的傳感器讀數
是
的測量值的概率,以及(2)先前估計值的概率認為
是我們應該看到的讀數。
如果我們有兩個概率,并且想知道兩個概率都為真的機會,則將它們相乘。因此,我們對兩個高斯分布進行了相乘處理:
兩個概率分布相乘得到的就是上圖中的重疊部分。而且重疊部分的概率分布會比我們之前的任何一個估計值/讀數都精確得多,這個分布的均值就是兩種估計最有可能配置(得到的狀態)。
事實證明,兩個獨立的高斯分布相乘之后會得到一個新的具有其均值和協方差矩陣的高斯分布!下面開始推公式。
合并兩個高斯分布
首先考慮一維高斯情況:一個均值為
,方差為
的高斯分布的形式為:
我們想知道將兩個高斯曲線相乘會發生什么。下圖中的藍色曲線表示兩個高斯總體的(未歸一化)交集:
將公式(9)代入公式(10),我們可以得到新的高斯分布的均值和方差如下所示:
我們將其中的一小部分重寫為
:
這樣一來,公式的形式就簡單多了!我們順勢將公式(12)和(13)的矩陣形式寫在下面:
其中
表示新高斯分布的協方差矩陣,
是每個維度的均值,
就是大名鼎鼎的“卡爾曼增益”(Kalman gain)。
公式匯總
我們有兩個高斯分布,一個是我們預測的觀測,另外一個是實際的觀測(傳感器讀數)
,我們將這兩個高斯分布帶入公式(15)中就可以得到二者的重疊區域:
我們有兩個高斯分布,一個是我們預測的觀測,另外一個是實際的觀測(傳感器讀數)
,我們將這兩個高斯分布帶入公式(15)中就可以得到二者的重疊區域:
從公式(14)我們可以知道,卡爾曼增益是:
然后我們將公式(16)與公式(17)中的
去除,同時將
后面的
去除,我們可以得到最終的化簡形式的更新方程:
圖說
大功告成,
就是更新后的最優狀態!接下來我們可以繼續進行預測,然后更新,重復上述過程!下圖給出卡爾曼濾波信息流。
總結
在上述所有數學公式中,你需要實現的只是公式(7)(18)和(19)。(或者,如果你忘記了這些,可以從等式(4)和(15)重新推導所有內容。)
這將使你能夠準確地對任何線性系統建模。對于非線性系統,我們使用擴展卡爾曼濾波器,該濾波器通過簡單地線性化預測和測量值的均值進行工作。
參考資料
[1]: How a Kalman filter works, in pictures, 圖解卡爾曼濾波是如何工作的:
?
[2]: A geometric interpretation of the covariance matrix, 協方差矩陣的幾何解釋:?
[3]: Kalman Filter 卡爾曼濾波:?
歷史精華好文
專輯1:AI工程落地
專輯2:AI核心算法
專輯3:AI優質資源
交流合作
請加微信號:yan_kylin_phenix,注明姓名+單位+從業方向+地點,非誠勿擾。
我目前認為的卡爾曼濾波器是這樣的:
1.設置狀態轉移矩陣
2.[預測]通過 t-1 時的狀態預測 t 時的狀態
3.[更新]用傳感器的測量值與 2 中的預測值估計出最優值,并更新卡爾曼增益等常數
4.重復 2~3 步驟
(如上正是opencv庫中自帶的卡爾曼濾波器例子的步驟,也是初學者能最直觀理解到的,它只涉及到一個測量值)
那么請問,在某些涉及多傳感器融合的問題中,例如 陀螺儀(累積誤差)和電子羅盤(干擾)的數據融合,此時有多個測量值存在,卡爾曼濾波器又是怎樣工作的呢?
貼一個數據融合的結果,用陀螺儀和GPS數據進行數據融合,濾波估計車輛前輪轉角即行業內所說的慣導角位移。LZ精準農業領域,器件型號不表,用霍爾傳感器(精度0.1度)標定,精度在0.17度,用于農機自動駕駛,真他娘的好用。。。。。。
補一個IMU姿態測量系統數據融合的結果。陀螺和加表數據融合,濾波估計系統俯仰角和橫滾角。器件型號不表。
這個是實驗室桌面實驗,手動翻轉IMU
下面是室外實驗,用GPS雙天線的俯仰角進行精度標定,GPS主副天線基線2米,俯仰角精度0.1度。IMU剛性鏈接在GPS工裝上,兩手托舉行走。
更一個濾波融合的例子。
霍爾傳感器安裝在車輛輪子處測量車輛前輪的角度,但是霍爾傳感器測量的是0到360度的絕對角度,安裝在車輪轉向節處時,在車輛前輪打正的情況下角度不是0而是一個隨機的固定的偏置,采用濾波估計的方法估計其安裝偏置,精度可以達到0.1度:
再更一個IMU或者GPS在載體上安裝角估計的應用。IMU或者雙天線GPS安裝在載體上,由于安裝誤差,IMU的前向或GPS主副天線的航向與載體的前向有一個固定的安裝誤差角,采用濾波估計的方法進行估計,精度可以達到0.2度:
在更一個組合導航的例子。
單天線GNSS+IMU做組合導航,提供位置速度姿態,IMU型號是ADI某較低成本器件,位置速度精度和RTK狀態下GNSS結果保持一致,航向精度0.12度(RMS)。
產品是農機自動駕駛儀,用來輔助農機的自動駕駛進行作業,目前主要應用在插秧機、拖拉機、割草機等四輪小車。
組合航向精度和雙天線GNSS航向精度基本一致,且靜止情況下航向不漂。
農機橫向控制精度<2.5cm RMS。
下一篇: PLC、DCS、FCS三大控
上一篇: 電氣控制線路圖控制原