之前寫完了

David LEE:卡爾曼濾波:從入門到精通?

zhuanlan.zhihu.com
圖標

之後一直說要寫非線性情況下的推導。但後來想想,EKF不過就是做了一階泰勒展開而已,要寫成一篇文章實在是太短,感覺也沒有太多乾貨,就一直沒有動筆。最近在看 《概率機器人》,覺得關於貝葉斯濾波器確實也可以稍作總結,於是就寫了這篇卡爾曼濾波家族。

本文對於擴展卡爾曼濾波、無跡卡爾曼濾波僅僅做了一些簡要介紹,不再想上次的文章那樣做詳細地推導了。但只要看過之前寫的卡爾曼濾波,相信這篇文章對於你來說也是很好理解的。

本文配圖均來自《概率機器人》

擴展卡爾曼濾波

假設狀態轉移概率和測量概率分別由非線性函數g和h控制,而不再是一個線性變換:

egin{aligned}  mathbf{x}_t & = g(mathbf{u}_t, mathbf{x}{t-1}) + varepsilon_t \ mathbf{z}_t & = h(mathbf{x}_t) + delta_t  end{aligned}

這種情況下,由於線性變換的關係不在了,因此概率分布也不再是高斯分布。整個系統不再有閉式解,這是最讓人頭疼的。

而EKF 的主要思想就是線性化:通過一個在高斯函數的均值處與非線性函數g相切的線性函數來近似g。

線性化的主要優點就是效率,一旦對g和h進行了線性化,KEF和KF就是等效的。

EKF採用一階泰勒展開的方式來進行線性化,其根據g的值和斜率構造一個函數g的線性近似函數:

g(mathbf{u}_t, mathbf{x}_{t-1}) = frac{partial g(mathbf{u}_t, mathbf{x}_{t-1})}{partial mathbf{x}_{t-1}}

線性點的選擇依據的是在線性化點附近自變數最有可能的狀態。對於高斯函數,最可能的狀態自然就是後驗的均值 mu_{t-1}

egin{aligned}  g(mathbf{u}_t, mathbf{x}_{t-1}) & approx g(mathbf{u}_t, mu_{t-1}) + g(mathbf{u}_t, mu_{t-1}) (mathbf{x}_{t-1} - mu_{t-1}) \ & = g(mathbf{u}_t, mu_{t-1}) + mathbf{G}_t (mathbf{x}_{t-1} - mu_{t-1})  end{aligned}

同理,將測量函數h線性化,有

egin{aligned}  h(mathbf{x}_t) & approx h(ar mu_t) + h(ar mu_t) (mathbf{x}_t - armu_t) \ & = h(ar mu_t) + mathbf{H}_t (mathbf{x}_t - ar mu_t)  end{aligned}

最後,整個EKF演算法的流程如下:

  1. 運動更新:

egin{aligned}  ar mu_t & = g(mathbf{u}_t, mathbf{x}_{t-1}) \  ar Sigma_t & = mathbf{G}_t Sigma_{t-1} mathbf{G}_t^T + mathbf{R}_t  end{aligned}

2. 測量更新

egin{aligned}  mathbf{K}_t & = ar Sigma_t mathbf{H}_t^T (mathbf{H}_t ar Sigma_t mathbf{H}_t^T + mathbf{Q}_t)^{-1} \ mu_t & = ar mu_t + mathbf{K}_t (mathbf{z}_t - h(ar mu_t)) \ Sigma_t & = (mathbf{I} - mathbf{K}_t mathbf{H}_t) ar Sigma_t  end{aligned}

很顯然,EKF 能否成功應用取決於兩個因素:

  1. 被近似的函數的局部非線性化程度;
  2. 概率分布自身的不確定度(協方差)。

上兩圖就明確展示了非線性函數在近似點非線性程度越高、概率分布本身越不確定,所得到的近似結果就越差。此時,採用擴展卡爾曼濾波的效果往往很差,甚至會導致發散。

無跡卡爾曼濾波

不同於EKF使用線性化來近似非線性函數,UKF通過無損變換來近似一個高斯分布,它通過使用加權統計線性回歸過程來實現隨機線性化。

下面就介紹些UKF無損變換的思想

UKF明確地從高斯分布中提取 幾個sigma 點,並將這些 sigma 點經過函數g變換。這些 sigma 點分別位於均值處、對稱分布於主軸的協方差處。對於具有均值 mu 和協方差 Sigma 的n維高斯分布,對應的2n+1個 sigmachi^{[i]} 根據如下規則進行選擇:

egin{aligned}  chi^0 & = mu \ chi^i & = mu + (sqrt{(n + lambda)Sigma)_i} quad i = 1, ..., n \  chi^i & = mu - (sqrt{(n + lambda)Sigma)_i} quad i = n + 1, ..., 2n  end{aligned}

其中, lambda = alpha^2 (n + kappa) - nalpha, kappa 是確定 sigma 點分布在離均值多遠的範圍內的比例參數(就是你的採樣範圍)。

每個 sigma 有兩個與之相關的權值:

  • omega_m^i 在計算均值時使用;
  • omega_c^i 在計算高斯協方差的時候使用。

egin{aligned}  omega_m^0 & = frac{lambda}{n + lambda} \  omega_c^0 & = frac{lambda}{n + lambda} + (1 - alpha^2 + eta) \  omega_m^i & = omega_c^i = frac{1}{2(n + lambda)} quad i = 1, ..., 2n  end{aligned}

選擇參數 eta 對高斯分布的附加的高階分布信息進行編碼。如果概率分布是精確的高斯分布,則 eta=2

然後, sigma 點經過函數g變換,來探測非線性函數g如何改變了高斯分布的形狀:

mathbf{y}^i = g(chi^i)

所得的結果高斯分布的參數 {ar mu, ar Sigma} 由映射的 sigmamathbf{y}^i 獲得:

egin{aligned}  ar mu & = sum_{i=0}^{2n} omega_m^i mathbf{y}^i \  ar Sigma & = sum_{i=0}^{2n} omega_c^i (mathbf{y}^i - ar mu)(mathbf{y} - ar mu)^T  end{aligned}

將以上過程分別代入到運動更新和測量更新中就可以得到無跡卡爾曼濾波了。

運動更新

1. 計算 sigma 點:

chi_{t-1} = (mu_{t-1} quad mu_{t-1} + sqrt{(n + lambda)Sigma_{t-1}}  quad mu_{t-1} - sqrt{(n + lambda)Sigma_{t-1}})

2. 非線性函數映射

ar chi^*_t = g(mu_t, chi_{t-1})

3. 進行運動更新

egin{aligned}  ar mu_t & = sum_{i=0}^{2n} omega_m^i ar chi^{* i}_t \  ar Sigma_t & = sum_{i=1}^{2n} omega_c^i ( ar chi^{ i}_t - ar mu_t)( ar chi^{ i}_t - ar mu_t)^T + mathbf{R}_t  end{aligned}

測量更新

1. 依據運動更新的結果計算新的 sigma

ar chi_t = (ar mu_t quad ar mu_t + sqrt{(n + lambda) ar Sigma_t}  quad ar mu_t - sqrt{(n + lambda) ar Sigma_t})

2. 計算每個sigma 點所對應的觀測值的預測

mathbf{ar Z}_t^i = h(archi_t^i)

3.計算預測的觀測值的不確定性

egin{aligned}  mathbf{hat z}_t & = sum_{i=0}^{2n} omega_m^i mathbf{ar Z}^i_t \  mathbf{S}_t & = sum_{i=0}^{2n} omega_c^i (mathbf{ar Z}_t^i - mathbf{hat z}_t)(mathbf{ar Z}_t^i - mathbf{hat z}_t)^T + mathbf{Q}_t  end{aligned}

4. 計算 sigma 點集在狀態空間和測量空間之間的互協方差

ar Sigma_t^{x, z} = sum_{i=0}^{2n} omega^i_c (ar chi_t^i - ar mu_t)(mathbf{ar Z}_t^i - mathbf{hat z}_t)^T

5. 進行測量更新

egin{aligned}  mathbf{K}_t & = ar Sigma_t^{x,z} mathbf{S}^{-1} \  mu_t & = ar mu_t + mathbf{K}_t (mathbf{z}_t - mathbf{hat z}_t) \  Sigma_t & = ar Sigma_t - mathbf{K}_t mathbf{S}_t mathbf{K}_t^T  end{aligned}

以上就是擴展卡爾曼濾波和無跡卡爾曼濾波的全過程了。只要明白了貝葉斯濾波或者卡爾曼濾波器,這兩個擴展都是信手拈來。

SLAM中,KEF往往比UKF使用得更廣。我個人淺見,主要原因有一下幾點:

1. UKF是一種抽樣近似,無可避免地導致計算量較大; 2. 只有在高度不線性的情況下(或方差很大),UKF才有明顯的優勢; 3. UKF的優勢主要集中在開環的情況下,但SLAM中更常使用迴環檢測這樣的閉環來消除誤差; 4. 如果建圖與定位分開,UKF在建圖過程中的作用可能會大些(開環的情況下);在定位中,用其來提高運動模型作用不大,觀測模型中準確地和地圖進行匹配對定位精度影響更大。

信息濾波

信息濾波(IF)是卡爾曼濾波的對偶濾波演算法,二者的不同在於對高斯分布的表示方式。信息濾波以正則參數來表示高斯分布,由一個信息矩陣和信息向量組成。

egin{aligned}  Omega & = Sigma^{-1} \  xi & = Sigma^{-1} mu  end{aligned}

對於一個高斯分布

p(mathbf{x}) = det(2piSigma)^{-1/2} exp { -frac{1}{2}(mathbf{x} - mu)^T  Sigma^{-1} (mathbf{x} - mu) }

經過變換後

egin{aligned}  p(mathbf{x}) & = det(2piSigma)^{-1/2} exp { -frac{1}{2}mathbf{x}^T Sigma^{-1}mathbf{x} + mathbf{x}^T Sigma^{-1}mu - frac{1}{2}mu^T Sigma^{-1} mu } \  & = det(2piSigma)^{-1/2} exp { - frac{1}{2}mu^T Sigma^{-1} mu } exp{ -frac{1}{2}mathbf{x}^T Sigma^{-1}mathbf{x} + mathbf{x}^T Sigma^{-1}mu } \  & = etaexp{ -frac{1}{2}mathbf{x}^T Sigma^{-1}mathbf{x} + mathbf{x}^T Sigma^{-1}mu } \  & = etaexp{ -frac{1}{2}mathbf{x}^T Omegamathbf{x} + mathbf{x}^T xi }  end{aligned}

該式在一些情況下會比矩參數(卡爾曼濾波中均值和協方差分別是一階矩和二階矩)的表示更加簡潔。例如,信息濾波下高斯分布的負對數為

-log p(mathbf{x}) = const. + frac{1}{2}mathbf{x}^T Omegamathbf{x} - mathbf{x}^T xi

此外,IF表示全局不確定性更加簡單,設置 Omega =0 即可。

對應的信息濾波演算法和擴展信息濾波演算法如下。

信息濾波演算法

對應的運動方程和測量方程為

egin{aligned} mathbf{x}_t & = mathbf{A}_t mathbf{x}_{t-1} + mathbf{B}_t mathbf{u}_t + varepsilon \ mathbf{z}_t & = mathbf{C}_t mathbf{x}_t + delta_t end{aligned}

對應的信息濾波演算法為:

egin{aligned}  ar Omega_t & = (mathbf{A}t Omega^{-1}{t-1} mathbf{A}^T_t + mathbf{R}_t)^{-1} \  ar xi_t & = ar Omega_t (mathbf{A}t Omega^{-1}{t-1} xi_{t-1} + mathbf{B}_t mu_t) \  \  Omega_t & = mathbf{C}_t^T mathbf{Q}_t^{-1}mathbf{C}_t + ar Omega_t \  xi_t & = mathbf{C}_t^T mathbf{Q}_t^{-1} mathbf{z}_t + ar xi_t  end{aligned}

對比信息濾波演算法和卡爾曼濾波,二者在計算複雜度上有很大的不同。

信息濾波中運動更新涉及矩陣求逆,計算量大;但在卡爾曼濾波中運動更新只涉及狀態向量的一部分。而信息濾波中測量更新是增量的,而KF的測量更新涉及矩陣求逆而較為困難。這也展示了二者的特性是對偶的。

擴展信息濾波演算法

egin{aligned}  mu_{t-1} & = Omega_{t-1}^{-1} xi_{t-1} \ ar mu_t & = g(mathbf{u}t, mu{t-1}) \  ar Omega_t & = (mathbf{G}t Omega^{-1}{t-1} mathbf{G}^T_t + mathbf{R}_t)^{-1} \  ar xi & = ar Omega_t ar mu_t \  \  Omega_t & = ar Omega_t + mathbf{H}_t^T mathbf{Q}_t^{-1}mathbf{H}_t \  xi_t & = ar xi_t + mathbf{H}_t^T mathbf{Q}_t^{-1}[ mathbf{z}_t - h(armu_t) + mathbf{H}_t ar mu_t]  end{aligned}

推薦閱讀:

相关文章