一、卡爾曼濾波器要解決的問題
首先說一下卡爾曼濾波器要解決的是哪一類問題,這類系統應該如何建模。這里說的是線性卡爾曼濾波器,顧名思意,那就是線性動態的離散系統。這類系統可以用如下兩個方程來表示:
\[\begin{array}{l}
x(n + 1) = {\bf{F}}(n + 1,n)x(n) + {v_1}(n) \\
y(n) = {\bf{C}}(n)x(n) + {v_2}(n) \\
\end{array}\]
其中:
x(n)表示系統的狀態
F(n+1,n)為狀態轉移矩陣,表示狀態隨時間的變化規律。通俗的講,從當前狀態到下一個狀態之間有什么關系。
C(n)表示觀測值與狀態的關系
y(n)表示狀態的觀測值
v1表示系統過程的噪聲
v2表示觀測過程中產生的噪聲
上面的兩個方程中,第一個方程是過程方程,它表示系統狀態x(n)隨時間的更新過程。第二個方程為測量方程,表示狀態x(n)與測量結果y(n)的關系。這里我們要先對這兩個方程中的概念做下解釋。
首先解釋下狀態這個概念。狀態是對系統特征進行的一個抽象,由預測系統未來特性時所需要的、與系統過去行為有關的最少數據組成。
這個概念不好理解吧!那么舉個例子。相信不少朋友在網上看到過有人拿房間溫度測量來講述卡爾曼濾波原理。這里房間里真實的溫度就是狀態,它可以是一個參數,也可以是多個參數。那么,用溫度計測出來的值,就是這里的觀測值y(n)。再說一個例子,假如我們要對一個運動的物體進行跟蹤,那么,物體的位移和速度完全可以表示這個運動物體所組成的系統的主要特征。這時的狀態就可以用一個具有位移和速度兩個特征的向量來表示。解釋到這里,相信很多朋友已經正確理解了狀態這個概念,它表示的是系統客觀存在的真實特征。
再說一下系統狀態與其觀測值之間為什么有C(n)的存在,這里它表示的是觀測值與狀態的關系。再拿室內測度測量來舉例子,室內客觀真實溫度(未知量)做為這個系統中的狀態,用溫度計來測量這個狀態。測出來的溫度就是我們的觀測溫度y(n)。這里很可能系統狀態跟其觀測值不是簡單的加個測量噪聲的關系。那么這種關系就可以用C(n)來建模。
上面的兩個方程,就是線性卡爾曼濾波器對要解決的系統的建模。這里再次對這兩個方程表示的模型做出更加通俗的解釋:
過程方程:它講述的是系統從一個狀態到另一個狀態應該隨時間如何變化
測量方程:它講述的是當前狀態與當前測量值之間的關系
也就是說,我們如果對一個系統感興趣,首先要找出這個系統的一個或者幾個主要特征(狀態),然后對這幾個特征進行觀測,並得到一組觀測值,通過不斷的觀測來認識這個系統,利用僅有的觀測數據找到最優的方式來求解過程方程和測量方程。這就是卡爾曼濾波器要解決的問題。為了便於理解,后面我們只討論只有一個特征的狀態所組成的觀測系統。
二、投影定理與新息過程
現在我們已經清楚了卡爾曼濾波器要解決什么樣的問題,對於這樣的問題應該如何建立模型。下面我們來看下具體用什么樣的方法來最優的求解過程方程和測量方程。卡爾曼濾波器的推導過程有兩種思路:
(1)正交投影定理(有的地方也稱為射影定理,這里后面簡稱為投影定理)
(2)新息過程
卡爾曼老先生最初介紹卡爾曼濾波器時是使用正交投影定理來推導的。后來Kailath又提出來基於新息過程的推導方法。這里打算兩個都介紹一下。上面說到,當對系統觀測時,會得到一組觀測值y(n),那好,我們首先考慮一個估計的問題,更具體的說是如何利用僅有的觀測數據來預測狀態x(n)的問題。由回歸分析的內容我們可以知道,由觀測值y求隨機變量x的線性最小二乘估計(線性最小方差估計)可以表示為:
\[\hat x = E[x] + Cov(x,y)D{(y)^{ - 1}}\{ y - E[y]\} \]
這時,隨機變量x與其估計的差值跟觀測值y正交(不相關),我們就稱x的估計值為x在y的上投影,記為:
\[\hat x = proj[x|y] = proj[x|y(1),y(2),...,y(k)]\]
其幾何意義表示如下:
好,再來思考一個問題,隨機變量x的真實值是我們最終想求的,但並不能真正得到。如果我們轉變一下思路,用前n-1個觀測值來預測當前n時刻的觀測值,會推出什么樣的結果。這里我們需要再定義一個新的概念--新息,它表示當前觀測值與其估計值的誤差。可以表示為下面的方式:
\[\alpha (n) = y(n) - proj[y(n)|y(1),y(2),...,y(n - 1)] = y(n) - \hat y(n|{Y_{n - 1}})\]
新息幾何意義表示如下:
再來看新息產生的過程,我們得到一個觀測值,就能同時得到一個對其對應的誤差。反之也是成立的,也就是說觀測值與新息,存在着一一對應的關系,那如果我們想要得到最終的系統狀態,也新息來代替觀測值來估計也是可以的,即:
\[\{ y(1),y(2),...,y(n)\} \mathbin{\lower.3ex\hbox{$\buildrel\textstyle\rightarrow\over
{\smash{\leftarrow}\vphantom{_{\vbox to.5ex{\vss}}}}$}} \{ \alpha (1),\alpha (2),...,\alpha (n)\} \]
因此
\[proj[x|y(1),y(2),...,y(n)] = proj[x|\alpha (1),\alpha (2),...,\alpha (n)]\]
寫成遞推的形式為:
\[\begin{array}{l}
proj[x|y(1),y(2),...,y(n)] = proj[x|\alpha (1),\alpha (2),...,\alpha (n)] \\
= proj[x|\alpha ] = E[x] + \sum\limits_{i = 1}^{n - 1} {E\{ x} {\alpha ^T}(i)\} E{\{ \alpha (i){\alpha ^T}(i)\} ^{ - 1}}\alpha (i) + E\{ x{\alpha ^T}(n)\} E{\{ \alpha (n){\alpha ^T}(n)\} ^{ - 1}}\alpha (n) \\
= proj[x|\alpha (1),\alpha (2),...,\alpha (n - 1)] + E\{ x{\alpha ^T}(n)\} E{\{ \alpha (n){\alpha ^T}(n)\} ^{ - 1}}\alpha (n) \\
= proj[x|y(1),y(2),...,y(n - 1)] + E\{ x{\alpha ^T}(n)\} E{\{ \alpha (n){\alpha ^T}(n)\} ^{ - 1}}\alpha (n) \\
\end{array}\]
看來,利用投影的思想,我們是可以遞歸的表達出來狀態的預測過程的。也可以說,在幾何上,kalman濾波器可以看做狀態在由觀測生成的線性空間上的投影。但如果只分析到這里,估計很多朋友可能還不能真正的理解卡爾曼濾波器的思想。那么我們再以新息過程來分析下如何來求解過系統狀態。
認真考慮下就會發現,整個過程有兩個誤差。一個是觀測過程中的預測誤差:新息;一個是狀態更新過程中的預測誤差(一直隱含在投影過程中沒有被注意到),我們來看下它們之間的關系。
由過去所有的觀測值y(1),...,y(n-1)和系統的測量方程,可以得到當前觀測值y(n)的最小均方估計為
\[\hat y(n|{Y_{n - 1}}) = {\bf{C}}(n)\hat x(n|{{\bf{Y}}_{n - 1}})\]
於是,新息過程還可以表示為:
\[\begin{array}{*{20}{c}}
{\alpha (n) = y(n) - \hat y(n|{{\bf{Y}}_{n - 1}})} \\
{ = y(n) - {\bf{C}}(n)\hat x(n|{{\bf{Y}}_{n - 1}})} \\
{ = {\bf{C}}(n)\varepsilon (n,n - 1) + {v_2}(n)} \\
\end{array}\]
這里ε(n,n-1)是狀態x(n)和其一步預測值之差
\[\varepsilon (n,n - 1) = x(n) - \hat x(n|{{\bf{Y}}_{n - 1}})\]
如果我們定義
\[\begin{array}{l}
{\bf{R}}(n) = E[\alpha (n){\alpha ^H}(n)] \\
{\bf{K}}(n,n - 1) = E[\varepsilon (n,n - 1){\varepsilon ^H}(n,n - 1)] \\
{{\bf{Q}}_2}(n) = {v_2}(n)v_2^H(n) \\
\end{array}\]
這里R(n)為新息誤差相關矩陣,代表觀測信息的質量優劣、K(n,n-1)為狀態誤差相關矩陣,代表狀態時間更新的質量優劣、Q2(n)表示測量噪聲向量的相關矩陣。再由觀測方向可得:
\[{\bf{R}}(n) = {\bf{C}}(n){\bf{K}}(n,n - 1){{\bf{C}}^H}(n) + {{\bf{Q}}_2}(n)\]
這就是狀態預測誤差與新息之間的二階統計關系。因此,只要我們盡可能的優化時間更新的質量,同時提高觀測信息的質量,那么,我們一定能估計出想要的系統狀態。好,問題反過來,再回到文初提出來的問題,我們有了一組觀測值,又得到了與比對應的新息序列。如何好好利用它們估計現我們想要的系統狀態。下面就要好好說下:
根據觀測值與新息之間的一一對應關系,如果我們以新息過程來估計狀的最小均方估計,則這個估計可以表示為新息過程序列的線性組合。即:
\[\hat x(i|{{\bf{Y}}_n}) = \sum\limits_{k = 1}^n {{{\bf{B}}_i}(k)\alpha (k)} \]
這里B為待定系數矩陣,那么,根據線性濾波的正交原理,預測狀態誤差向量與新息過程正交,即:
\[E[\varepsilon (i,n){\alpha ^H}(m)] = E\{ [x(i) - x(i|{{\bf{Y}}_n})]{\alpha ^H}(m)\} = 0\]
再將新息過程導出來狀態x(i)的最小均方估計代入,可得
\[E[x(i){\alpha ^H}(m)] = {{\bf{B}}_i}(m)E[\alpha (m){\alpha ^H}(m)] = {{\bf{B}}_i}(m){\bf{R}}(m)\]
這里的R(m)為新息過程的相關矩陣。上式兩邊同時乘以R(m)的逆矩陣,可得待定系數矩陣B的表達式為
\[{{\bf{B}}_i}(m) = E[x(i){\alpha ^H}(m)]{{\bf{R}}^{ - 1}}(m)\]
再將B代入新息過程導出來狀態x(i)的最小均方估計公式中,可以得到狀態x(i)的最小均方誤差估計為:
\[\begin{array}{l}
\hat x(i|{{\bf{Y}}_n}) = \sum\limits_{k = 1}^n {E[x(i){\alpha ^H}(k)]{{\bf{R}}^{ - 1}}(k)\alpha (k)} \\
= \sum\limits_{k = 1}^{n - 1} {E[x(i){\alpha ^H}(k)]{{\bf{R}}^{ - 1}}(k)\alpha (k)} + E[x(i){\alpha ^H}(k)]{{\bf{R}}^{ - 1}}(k)\alpha (k) \\
\end{array}\]
對於i=n+1,有
\[\hat x(n + 1|{{\bf{Y}}_n}) = \sum\limits_{k = 1}^{n - 1} {E[x(n + 1){\alpha ^H}(k)]{{\bf{R}}^{ - 1}}(k)\alpha (k)} + E[x(n + 1){\alpha ^H}(k)]{{\bf{R}}^{ - 1}}(k)\alpha (k)\]
我們看到,這個結果與投影定理推導的結果的含義是相同的,只是使用的工具和准則不同。
最后,我們再回想下觀測數據與新息、狀態之間的關系,以及根據這些關系導出的結果。不難得出結論,卡爾曼濾波器的出發點是盡可能的好好利用觀測值,用盡可能少的計算量來估計系統狀態。呵呵,是不是感覺數字的游戲還是很好玩的,並不一定很枯燥。好的,今天就分析到這里吧,后續將會使用這次分析的結果,逐漸對卡爾曼濾波器及其使用做一個說明。這里的分析內容來源於我在學習卡爾曼濾波器過程中的一些思考(參閱書籍:自適應濾波器原理、卡爾曼濾波原理與應用),水平有限,如有錯誤之處,還請各位大牛不吝賜教。