3.4 粒子フィルタ

土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
山本真哉（清水建設）
3.4 粒子フィルタ
特長
粒子フィルタはアンサンブルカルマンフィルタと同様にアンサンブルベースの逐次ベイ
ズフィルタの一種である。しかし、アンサンブルカルマンフィルタでは状態変数や観測デ
ータがガウス分布に従わない場合や、観測モデルが線形関数でない場合に、妥当な計算結
果が得られないことがあるのに対し、粒子フィルタでは、そのような仮定を一切必要とせ
ず、より一般的な状態空間モデルを対象にできるのが最大の特長である。また、アルゴリ
ズムが単純なため、実装が容易であることも長所の一つとして挙げられる。
粒子フィルタは、非線形性や非ガウス性を有する状態空間モデルを対象としたフィルタ
手法[1]をモンテカルロ近似で解く方法として、モンテカルロ・フィルタ[2]、またはブート
ストラップフィルタ[3]という名前で提案されたのが始まりであるが、現在では一般的に粒
子フィルタと呼ばれ、シミュレーションモデルと観測データを逐次的に同化する手法とし
て認知されている。
アルゴリズム
いま、以下のような非線形性および非ガウス性を考慮した状態空間モデルを対象として
粒子フィルタによる状態推定を考える。
xt  F ( xt 1 , vt )
(1)
y t  H ( xt , wt )
(2)
ここで、システムモデル F および観測モデル H は非線形の関数である。また、モデル誤差
に相当するシステムノイズ v と観測ノイズ w の確率分布は必ずしもガウス分布に従うもの
ではない。
粒子フィルタにおける状態変数の事前確率密度分布（以降、事前分布と記す）を求める
手順はアンサンブルカルマンフィルタと共通であり、前時間ステップで事後分布として得
た個々のアンサンブルメンバー（慣習的に粒子と呼ばれる）について（1）式のモデル方程
式により時間更新を解くことで行う。いま、状態変数 x の事前分布を N 個の粒子により近
似すると、ディラックのデルタ関数を用いて次のように表わされる。
p ( xt ) 
1 N
 ( xt  xti )

N i 1
(3)
ただし、xi はモデル方程式より求めたある粒子の予測値であり、i は粒子のインデックスを
示している。
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
山本真哉（清水建設）
次に、新たに観測データ y が与えられた際の状態の事後分布はベイズの定理より次式のよ
うに与えられる。
p ( xt | y t ) 
p ( y t | xt ) p ( xt )
(4)
 p( yt | xt ) p( xt )dx
ここで、P(x|y)は条件付き確率を表し、y が与えられた条件の下で x が得られる確率を意味
する。一方、P(y|x)は尤度であり、観測 y が得られたときの x の尤もらしさを表す。さらに
上式に式(3)を代入すると事後分布は次のような形で表される。
N
p ( xt | yt )   i ( xt  xti )
(5)
i 1
ただし、
i 
p( yt | xti )
N
 p( y
i 1
t
(6)
i
t
|x )
である。ω は総和が 1 となるように正規化された各粒子の尤度であり、すなわち、個々の
粒子の予測結果と観測結果との近さを表している。これより、式(5)の事後分布は尤度によ
って重みづけられた事前分布の粒子で近似できることが分かる。例えば、観測ノイズが平
均 0、共分散行列 R のガウス分布に従うと仮定した場合、尤度は次式により算出できる[4]。
p ( y t | xi ) 




T
 1

exp  yt  H ( xti ) R 1 y  H ( xti ) 
 2

(2 ) R
1
l
(7)
以上に述べた手順を観測結果が得られる度に実行することで、逐次的に状態を推定する
のが粒子フィルタの方法の流れとなる。
ところで，粒子フィルタの実行を繰り返すと、特定の一つの粒子の重みが 1 に近い値と
なり、それ以外の粒子の重みがほぼ 0 を示すようになる。この結果、重みが 0 に近い粒子
の計算結果がフィルタリング時に全く意味をなさなくなるとともに、事後分布を単一の粒
子のみで近似することになるため、状態の推定精度が大きく低下してしまう。このような
現象はフィルタの退化(degeneracy)と呼ばれる。粒子フィルタを適用する際には、退化の抑
制を目的として粒子のリサンプリングを行うのが一般的であり、このアルゴリズムを特に
SIR（Sequential Importance Resampling）法と呼ぶ。具体的には、式(5)において、粒子
の重みに従って粒子の数だけ復元抽出を行う。これにより、重みが小さい、言い換えれば、
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
山本真哉（清水建設）
観測結果から遠い粒子を破棄し、代わりに重みの大きい粒子を複製することで退化を抑え
ることができる（図 1）。
粒子フィルタのアルゴリズム（SIR）をフローとして図 2 に示す。このフローにおいて、
全粒子に関する繰り返し処理内部で、適用対象とする数値解析のソルバーをサブルーチン
として実行するのが一般的である。また、この繰り返し処理は独立して実行できるため、
容易に並列化が可能であり、その効果も高い。粒子フィルタの実装上の注意点やリサンプ
リングのアルゴリズム等については文献[4]が参考となる。
t-1
t
t+1
時刻
事後分布
事前分布
p(yt|xt)
尤度
x
p(xt|yt-1)
x
x
x
p(xt-1|yt-1)
p(xt|yt)
p(xt+1|yt)
事後分布
事前分布
図 1 SIR 法における事前分布と事後分布の関係
事前分布 p(xt|yt-1)は、直前の事後分布 p(xt-1|yt-1)を構成する粒子が非線形モデル方程式によ
り時間発展（予測）したもので表現される。また、SIR では観測結果から得られる尤度分
布 p(yt|xt)に基づいて p(xt|yt-1)の粒子をリサンプリングしたものが同時刻の事後分布 p(xt|yt)
となる。このとき、復元抽出であるため、尤度が大きい粒子ほど複数回選ばれる可能性が
高くなる。この例で t+1 において同じ粒子が異なる予測結果を示しているのは、システム
ノイズを加えているためである。
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
山本真哉（清水建設）
図 2 粒子フィルタの処理の流れ（N は総粒子数）
粒子フィルタの短所
粒子フィルタはアンサンブルカルマンフィルタと比較して取り扱うことができる状態空
間モデルの自由度は高いものの、状態空間の次元数の増加と共に確率分布を適切に近似す
るために必要な粒子数が飛躍的に増加するため[5]、高次元の問題に対しては粒子フィルタ
の適用が困難となる。また、一般的にアンサンブルカルマンフィルタよりも多くの粒子数
が必要となることから、計算コストの高いモデル方程式には適さない。
さらに、粒子フィルタには退化の問題が存在するため、リサンプリングを実施するのが
一般的であるが、リサンプリングによりごく限られた粒子だけが複製されて残り続けるこ
とになり、結果的にアンサンブルの多様性が失われる。この場合、システムノイズを適切
に与えないと結局は粒子フィルタの退化に繋がる。
したがって、非線形フィルタによる状態推定を成功させるためには対象とする問題の性
質や規模に応じて適切な手法を選択することが重要となる。
土木学会応用力学委員会逆問題小委員会ホームページ逆問題副読本
山本真哉（清水建設）
【参考文献】
[1] Kitagawa, G.: Non-Gaussian State—Space Modeling of Nonstationary Time Series,
Journal of the American statistical association, 82(400), pp.1032-1041. 1987.
[2] 北川源四郎：モンテカルロ・フィルタおよび平滑化について，統計数理，第 44 巻第１
号，pp.31-48, 1996.
[3] Gordon, Neil J., David J. Salmond, and Adrian FM Smith.: Novel approach to
nonlinear/non-Gaussian Bayesian state estimation, IEE Proceedings F (Radar and
Signal Processing), Vol. 140, No. 2, IET Digital Library, 1993.
[4] 樋口知之，上野玄太，中野慎也，中村和幸，吉田亮：データ同化入門－次世代のシミ
ュレーション技術－，朝倉書店，pp.47-77，2011.
[5] Snyder C, Bengtsson T, Bickel P, Anderson J.: Obstacles to High-Dimensional
Particle Filtering, Mon Weather Rev, 136(12), pp.4629–4640, 2008.

Download Report