カルマンフィルタの導出

はじめに

wiki: 「カルマンフィルタ (Kalman Filter、KFと略す) は、誤差のある観測値を用いて、ある動的システムの状態を推定あるいは制御するための、無限インパルス応答フィルタの一種である」。観測値に観測雑音、状態予測値にシステム雑音があって場合により時間とともに誤差がドンドン蓄積してそのまま使えないので、観測と予測のガウス性を利用した線形センサデータフュージョンがカルマンフィルタの原点である。

状態空間モデル

時系列解析の中で、予測値と観測値の間何らかの因果関係を見つけて、何らかの方法でそれらのデータを絡んでモデル化して状態を推定していく。ここで汎用的な状態方程式観測方程式は以下の式にする。
$$\begin{eqnarray*}&& x_t = F_{t-1}(x_{t-1} ) + q_{t-1}\\&& y_t = H_t(x_t) + r_t\end{eqnarray*}$$
ただし、システム雑音、観測雑音\((q_t, r_t)\)は期待値が\(0\)、分散が\((Q_t, R_t)\)の独立正規分布(ガウス分布)にする。
$$\begin{eqnarray*}&& p(q_t) ~ N(0,Qt)\\&& p(r_t) ~ N(0,Rt)\\&& E[q_t r_t]=0 \end{eqnarray*}$$

基本原理

カルマンフィルタの原理は以下のイメージに示すように、最小解析誤差共分散推定(最小二乗)から予測値、観測値をそれぞれの割合\(I-K_t H_t , K_t\)で線形合成した値を状態推定値にする。

カルマンフィルタの原理
カルマンフィルタの原理

式の導出

時系列データに適用して、以下予測Predict→更新Correctまたは、時間更新Time Update→測定更新 Measurement Updateの繰り返しによって、状態の数学期待最小解析誤差共分散を求める手順になる。

Step0=パラメータ初期化(t=0)

状態推定(数学期待)、解析誤差共分散の初期化
$$\begin{eqnarray*}&& \hat{x}_0 = E[x_0]\\&& \hat{P_0} = E[(x_0-\hat{x_0})(x_0-\hat{x_0})^T]\end{eqnarray*}$$

Step1=状態値、共分散行列予測(t>0)

$$\begin{eqnarray*}&& \bar{x_t} = F_{t-1}(\hat{x}_{t-1})\\&& \bar{P_t} = F_{t-1}\hat{P}_{t-1}F_{t-1}^T + Q_{t-1}\end{eqnarray*}$$
ただし、\(\bar{x_t}\)は予測値、\(\bar{P_t}\)は事前共分散行列、\(Q_{t-1}=E[q_{t-1}q^T_{t-1}]\)

Step2=カルマンゲイン、状態値、共分散行列更新

$$\begin{eqnarray*}&& K_t = \bar{P_t}H_t^T(H_t\bar{P_t}H_t^T + R_t)^{-1}\\
&& \hat{x_t} = \bar{x_t} + K_t(y_t – H_t\bar{x_t})\\&&\hat{P_t} = (I-K_tH_t)\bar{P_t} \end{eqnarray*}$$
ただし、\(K_t\)はカルマンゲイン、\(\hat{x_t}\)は状態推定値、\(\hat{P_t}\)は事後共分散行列、\(R_{t}=E[r_{t}r^T_{t}]\)

KFフローチャート

Kalman_Filter_Flowchart_2
Kalman_Filter_Flowchart_2

拡張カルマンフィルタ

拡張カルマンフィルタ(Extended Kalman Filter、EKFと略す)は、非線形フィルタリングである。前述した状態方程式、観測方程式より、以下の状態空間モデルの\(f(⋅)\)または\(h(⋅)\)が非線形関数であり、拡張カルマンフィルタが適用される。テイラー展開より、2次微分以降の項目を省略して、非線形である\(f(⋅), h(⋅)\)の1次微分を線形化とし、前述したカルマンフィルタのアルゴリズムが適用可能となる。しかし、\(f(⋅), h(⋅)\)の微分では\(f(⋅), h(⋅)\)の一部しか表現できず、この線形化処理(1次微分)が誤差を大きく招く可能性がある。
$$\begin{eqnarray*}&& F_{t-1} =\frac{\partial f_{t-1}(x)}{\partial x}|_{x=\hat{x}_{t-1}}\\&& H_{t} =\frac{\partial h_{t}(x)}{\partial x}|_{x=\hat{x}_{t}}\end{eqnarray*}$$

状態空間モデル

$$\begin{eqnarray*}&& x_t = f_{t-1}(x_{t-1}) + q_{t-1}\\&& y_t = h_t(x_t) + r_t\end{eqnarray*}$$

Step0=パラメータ初期化(t=0)

$$\begin{eqnarray*}&& \hat{x}_0 = E[x_0]\\&& \hat{P_0} = E[(x_0-\hat{x_0})(x_0-\hat{x_0})^T]\end{eqnarray*}$$

Step1=状態値、共分散行列予測(t>0)

$$\begin{eqnarray*}&& \bar{x_t} = f_{t-1}(\hat{x}_{t-1})\\&& \bar{P_t} = F_{t-1}\hat{P}_{t-1}F_{t-1}^T + Q_{t-1}\end{eqnarray*}$$
ただし、\(\bar{x_t}\)は予測値、\(\bar{P_t}\)は事前共分散行列、\(Q_{t-1}=E[q_{t-1}q^T_{t-1}]\)

Step2=カルマンゲイン、状態値、共分散行列更新

$$\begin{eqnarray*}&& K_t = \bar{P_t}H_t^T(H_t\bar{P_t}H_t^T + R_t)^{-1}\\
&& \hat{x_t} = \bar{x_t} + K_t(y_t – h_t\bar{x_t})\\&&\hat{P_t} = (I-K_tH_t)\bar{P_t} \end{eqnarray*}$$
ただし、\(K_t\)はカルマンゲイン、\(\hat{x_t}\)は状態推定値、\(\hat{P_t}\)は事後共分散行列、\(R_t=E[r_t r^T_t]\)

EKFフローチャート

Extended_Kalman_Filter_Flowchart_2
Extended_Kalman_Filter_Flowchart_2

拡張カルマンフィルタを6軸IMUへの適用

6軸IMU~拡張カルマンフィルタ

カルマンフィルタの再考

カルマンフィルタでは、状態推定値を予測結果\(x_t\)(実装例ではジャイロセンサーデータ)と観測データ\(y_t\)(実装例では加速度センサデータ)の線形結合で合成し,その誤差分散を最小にする推定法だと分かる。これはシステム誤差、観測誤差の数学期待が0の正規分布との前提条件から由来した推定法である。しかし、カルマンフィルタをかけることで、状態推定値は予測結果と観測データの間にあるのは、真値からかなり乖離してしまう場合にあるのか。これはシステム誤差と観測誤差が無相関かつ直交という前提から、勿論真値が推定値と観測値の間にある結論を結ぶ考えである。また、ジャイロセンサーデータと、加速度センサデータとも観測値\(y_t\)にする方法がある。それぞれのパーフォーマンスの検証は、比較的精確な実験環境(比較用の高精度ジャイロセンサ、加速度センサ、エンコーダ、モータ)がないと、実は容易ではない。というよりも、シミュレーションをかけてカルマンフィルタのアルゴリズムを検証するのが、確実に可能である。

参考文献

1-wikipedia: カルマンフィルタオイラー角
2-Greg Welch氏、Gary Bishop氏: An Introduction to the Kalman Filter
3-田島洋氏: マルチボディダイナミクスの基礎―3次元運動方程式の立て方

関連記事

6軸IMU MPU-6050をロボットに組み込もう
YDLIDAR G4=16m 薄型 ROS対応SLAM LIDAR
研究開発用 台車型ロボット キット
オイラー角~ジンバルロック~クォータニオン
SLAM~拡張カルマンフィルタ
SLAM~Unscentedカルマンフィルタ

ロボット・ドローン部品お探しなら
ROBOT翔・電子部品ストア