真値

土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
珠玖隆行（岡山大）
拡張カルマンフィルター
1
はじめに
最も基本的な非線形フィルタリングである拡張カルマンフィルター（Extended Kalman filter,
EKF）について説明する．はじめに，カルマンフィルター（Kalman filter, KF）の基礎となる線形
最小分散推定（linear minimum mean square estimate, LMMSE）の理論について概説し，カルマンフ
ィルターを誘導した後，拡張カルマンフィルターを誘導する．なお，ここでの説明や誘導方法は
文献[1]に負うところが大きい．
2
線形最小分散推定（linear minimum mean square estimate, LMMSE）
カルマンフィルターの誘導にはいくつかの方法があるが（例えば，[2] [3] [4]），ここでは著者が
分かりやすいと感じた線形最小分散推定に基づいた誘導方法を示す．よって，はじめにカルマン
フィルターの基礎となる線形最小分散推定について簡単に説明する．
はじめに，xt を時刻 t のある地点におけるある物理量の予測結果とし，yt を時刻 t のある地点に
おけるある物理量の観測データと定義する．ここで，最適な推定値 xtest を 2 つのデータ xt と yt か
ら推定する問題を考える．
xt と yt の真値を xttrue とし，xt と yt の誤差をそれぞれ pt，vt とすると，xt，yt は以下のように表す
ことができる．
xt
xttrue
yt
true
t
pt
x
(1a)
vt
(1b)
最適な xtest を推定するにあたり， 誤差 pt と vt が，
E pt
E vt
0
(2)
E pt pt
Pt , E vt vt
E pt vt
E vt pt
Rt
0
(3)
(4)
を満たし，かつ各誤差が正規分布（ガウス分布）に従うと仮定する．式(3)中の Pt，Rt は誤差分散
を表す．
以上の準備から，最適解は，推定値 xtest を 2 つのデータ xt と yt の線形結合で表し，この推定値
の誤差分散を最小にする重み x，
y を求めることで得られる．
xtest
x
x t
y
yt
(5)
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
珠玖隆行（岡山大）
式(5)の期待値をとると，
E xtest
E
x
E
x t
true
x t
yt
y
true
y t
x
x
(6)
が得られ，
1
x
(7)
y
のとき，xtest は不偏推定値となる．式(7)を式(5)に代入すると，
xtest
が得られ，xtest の誤差分散
2
2
E ( xtest
x
) yt
(8)
xttrue ) 2
E{
x t
x
E{
x
(1
P
x
( xttrue
E pt2
2
x t
2
(1
は次式で表される．
2
x
推定値の誤差分散
x
x t
pt ) (1
2
x
(1
が最小となる重み
xttrue }2
) yt
x
(1
x
x
) 2 Rt
)( xttrue
) E p t vt
(1
y
xttrue }2
vt )
(1
x
) 2 Pt
2
y
) 2 E vt2
Rt
(9)
y は，式(9)を微分することにより，以下のように求められ
る．
2
d(
d
)
2 Pt
2
Pt
以上より，最適な推定値 xtest とその誤差分散
Rt
Pt
2
y
Rt
0
(10)
Pt
y
xtest
P
y t
y
Rt
2
2
は次式で表される．
Pt
xt
Pt
(1
y
Rt
yt
) 2 Pt
2
y
Rt
Pt
Pt
xt
Pt
Rt
( yt
xt )
(12)
Rt
2
Rt
Pt
(11)
Rt
2
Pt
Pt
Rt
Rt
Rt Pt
Pt Rt
(13)
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
珠玖隆行（岡山大）
これまで示したように，推定値を予測結果 xt と観測データ yt の線形結合で表し，その誤差分散を
最小にする推定法は線形最小分散推定と呼ばれる．
Fig.1 に，推定値 xtest と予測結果 xt，観測データ yt それぞれの関係を示す．予測結果 xt と観測デ
ータ yt は真値からそれぞれ距離 pt，vt に位置する．xt と yt は無相関であるため，真値 xttrue から xt，
yt へ向かうベクトルは直交関係にある．
xtest
xt
(|
x
t
est
t
x
t
yt
(|
| pt )
y
t
| vt )
xttrue
Fig.1 線形最小分散推定の幾何学的解釈[1]
最適な推定値は，予測結果 xt，観測データ yt を通る直線上の，真値からの距離が最小となる点と
なる．また，Fig.1 から式(13)を幾何学的に導くことも可能である．図より最適推定値の誤差は元
の誤差の推定値の誤差よりも減少していることがわかる．
3
カルマンフィルター（Kalman filter , KF）
前節の線形最小分散推定では，ある一時点のある一地点における予測結果と観測データから，
最適な推定値を求める問題を考えた．ここでは時系列データに対する最適な状態推定法であるカ
ルマンフィルターを，線形最小分散推定に基づき導出する．
対象とするシステム
（問題）において，
シミュレーションによって得られた予測結果を x1, x2, x3,...,
xT とする．xt は時刻 t における対象とするシステムの解析結果（状態変数）を表す多次元のベク
トルを意味する．時刻 t の関数である状態変数 x の時間発展を与える線形作用素を Ft-1(・)とおく
と，式(1a)は以下のように書き換えられる．
xt
Ft 1 ( xt 1 ) wt
Ft 1 ( xttrue1
1
pt 1 ) wt
1
(14)
ここで wt-1 はシミュレーションによる時間発展に伴う誤差（システムノイズ，システム誤差）ベ
クトルである．式(14)のように状態の動的な性質を表すモデルは，一般的にシステムモデルと呼ば
れる．
一方，対象とするシステムで観測されたデータを y1, y2, y3,..., yT とし，yt は時刻 t において対象と
するシステムの種々の地点で観測される観測データを表す多次元のベクトルとする．観測データ
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
珠玖隆行（岡山大）
に関しては以下の表現を用いる．
H t ( xttrue ) v t
yt
(15)
ここで，Ht は観測行列と呼ばれ，全物理量（状態変数）から一部のデータを取り出すためのマト
リックスである．このように，システム方程式と観測方程式を連立させたシステムの表現（モデ
ル）は状態空間表現（状態空間モデル，state space model）と呼ばれ，カルマンフィルターをはじ
めとする各種フィルタリング手法を説明する際に用いられる．
以上の準備より，解析誤差共分散行列 Pt を最小化する問題を考える．推定値 xtest は，線形最小
分散推定における説明と同じように，式(11)おける Pt/(Pt+Rt)を重みマトリックス Wt とおくと，
x test
xt
Wt ( yt
Ft 1 ( xttrue1
xttrue
H t xt )
pt 1 ) wt
Ft 1 pt
1
wt
Wt {( H t xttrue
1
1
Wt ( v t
H t Ft 1 pt
vt ) H t ( Ft 1 ( xttrue
1
pt 1 ) wt 1 )}
H t wt 1 )
(16)
と表されるため，式(9)は以下のように書き換えられる．
Pt
E ( xtest
xttrue ) 2
E ( Ft 1 pt
1
wt
1
Wt vt Wt H t Ft 1 pt
Ft 1 Pt 1Ft T1 2 Ft 1 Pt Ft T1 H tTWt T
Qt
1
1
Wt H t w t 1 ) 2
2Qt 1 H tTWt T
Wt RtWt T Wt H t Ft 1 Pt 1 Ft T1 H tTWt T Wt H t Qt 1 H tTWt T
(17)
最適な重みマトリックスを求めるために，行列の微分公式を用いて上式を Wt で微分し，Wt に関
してまとめると次式が得られる．
tr ( Pt )
Wt
2 Ft 1 Pt 1Ft T1 H tT
Wt
2Qt 1 H tT
2Wt Rt
( Ft 1 Pt 1Ft T1 Qt 1 ) H tT {Rt
2Wt H t Ft 1 Pt 1Ft T1 H tT
2Wt H t Qt 1 H tT (18)
H t ( Ft 1 Pt 1Ft T1 Qt 1 ) H tT }
1
(19)
式(19)は式(12)をマトリックス表現したものであり，最適な重み Wt がシステムノイズと観測ノイ
ズの和に対するシステムノイズの比を表している．カルマンフィルターでは行列 Wt をとくにカル
マンゲイン Kt と呼び，多くの教科書では式(19)の代わりに次式が用いられている．
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
珠玖隆行（岡山大）
Pt H tT ( Rt
Kt
Pt
H t Pt H tT )
Ft 1 Pt 1Ft T1
Qt
1
(20)
(21)
1
ここで，Pt は時刻 t における解析誤差共分散行列である．式(21)から明らかなように，Pt は時刻 t-1
における諸量から計算することができ，観測データが得られていない状態での解析誤差共分散行
列の時間発展を表している．このような解析誤差分散の時間発展を表す式はリヤプノフ方程式
（Lyapunov equation）と呼ばれている．
これまでの説明を整理するために，時刻 t = 0～T でシステムの状態推定を行う場合のカルマン
フィルターのフローを Fig.2 に示す．
Start
Initial parameters
xt , Pt (t
0)
t
t 1
Time updating
xt
Pt
Ft 1 ( xt 1 ) wt
T
t 1
Ft 1 Pt 1F
Qt
1
1
Observation data
yt
Observation updating
K t Pt H tT ( Rt H t Pt H tT )
xtest xt K t ( yt H t xt )
No
1
t T?
Yes
End
Fig.2 カルマンフィルターのフロー
4
拡張カルマンフィルター（extended Kalman filter , EKF）
カルマンフィルターでは，線形なシステム方程式（シミュレーションモデル）に対して最適な
推定値を求めたが，実際に取り扱う問題の多くは非線形である．ここでは非線形システムに適用
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
珠玖隆行（岡山大）
可能な拡張カルマンフィルターについて簡単に説明する．
以下に示す非線形状態空間モデルに対して，カルマンフィルターを適用することを考える．
xt
f t ( xt 1 )
wt
yt
H t ( xttrue ) v t
(22a)
1
(22b)
ここに，ft (・)は非線形作用素を表す．カルマンフィルターを非線形モデルに適用するために，非
線形作用素 ft を状態ベクトル xt-1 の推定値 x test1 の周りで線形化（微分）し，式(14)と同様な線形作
用素 F̂t として扱うことを考える．
Fˆt
ft ( x)
x
(23)
x
xtest1
F̂t はいわゆるヤコビ行列である．以上の準備から，カルマンフィルター説明する際に用いた式
(13)と(20)の代わりに
xt
f t ( x test1 )
Pt
Fˆt 1 Pt 1Fˆt T1
wt
(14)’
1
Qt
1
(21)’
を用いることで，カルマンフィルターのアルゴリズムが非線形システムの最適推定問題に適用で
きる．これを拡張カルマンフィルターという．なお，非線形観測モデル ht を扱う場合でも，シス
テムモデルと同様に線形化してカルマンフィルターの更新式を使うことができる．
最後に，拡張カルマンフィルターの問題点について触れる．拡張カルマンフィルターでは各時
点で微分による線形化を行うため，解析的に微分できない場合には数値微分を行う必要がある．
また，高次項を無視しているため，非線形性の強いシステムに対しては計算が不安定になること
が知られている．このような問題に対しては，アンセンテッドカルマンフィルター（unscented
Kalman filter，UKF）やアンサンブルカルマンフィルター（ensemble Kalman filter, EnKF）が有効と
なる．
参考文献
[1] 淡路敏之，蒲池政文，池田元美，石川洋一：データ同化 観測・実験とモデルを融合するイノ
ベーション，京都大学出版会，pp.15-19，2009．
[2] 樋口知之，上野玄太，中野慎也，中村和幸，吉田
亮：データ同化入門－次世代のシミュレー
ション技術－，朝倉書店，pp.47-77，2011．
[3] 片山 徹：新版 応用カルマンフィルタ，朝倉書店，pp.83-107，2000．
[4] 谷萩隆嗣：カルマンフィルタと適応信号処理，コロナ社，pp.1-47，2005．