どすえのブログ

ソフトウェア開発ブログ

非線形カルマンフィルタ (2) - アンサンブルカルマンフィルタ -

目次

本シリーズについて

非線形現象に対するカルマンフィルタを本で勉強したのでメモも兼ねてまとめます.

前回の記事 dosuex.hatenablog.com では状態遷移、観測がともに線形である(行列演算で記述できる)ケースに適用できるカルマンフィルタを紹介しました。

今回は非線形システムを扱うことができるアンサンブルカルマンフィルタ(Ensemble Kalman Filter; EnKF)を紹介します. アンサンブルカルマンフィルタは, モンテカルロ近似を用いて, 非線形のシステムモデルをそのまま使って分散共分散行列の推定を行います。

コード( Jupyter Notebook )はこちらです.

github.com

今回はシステムモデルに関しては非線形モデルを, 観測モデルに関しては線形モデルを仮定する.

 
\begin{align}
x_{t} &= f_t(x_{t-1}, v_t) \tag{1} \\
y_t &= H_t x_t + w_t \tag{2} \\
\end{align}

x_t状態ベクトル,  y_t は観測データを表す. 観測ノイズは w_t \sim N(0,R_t)としてガウス分布からサンプルされるノイズを仮定する. 一方, システムノイズ v_t \sim q_t(v_t)は任意の分布で良い.

EnKFではアンサンブル近似を用いて分布の近似を行う. アンサンブル近似とは, 以下の図のようにある確率変数 xの確率分布を N個のサンプルの集合を用いて近似することである.

f:id:dosuex:20210517222636p:plain

このサンプルの集合 \{x^{(i)}\}_{i=1}^N をアンサンブルを呼ぶ. またそれぞれの x^{(i)} のことを粒子と呼ぶ.

EnKFでは, (1)式, (2)式で与えられる状態空間モデルのもとで, カルマンフィルタと同様に1期先予測とフィルタの手順を繰り返す.

カルマンフィルタでは平均値と共分散行列を逐次的に計算したが, EnKFでは予測分布とフィルタ分布のアンサンブルを逐次的に計算することになる.

これを改めて式で表す.

時刻t-1でのフィルタ分布p(x_{t-1}|y_{1:t-1})を近似するN個の粒子から構成されるアンサンブル \{x^{(i)}_{t-1|t-1}\}_{i=1}^N が準備されているものとする. これは

 
\begin{align}
p(x_{t-1}|y_{1:t-1}) \approx \frac{1}{N}\sum_{i=1}^N \delta (x_{t-1} - x^{(i)}_{t-1|t-1})
\end{align}

と書ける. このフィルタアンサンブルを使って時刻tでの観測y_tを取り込んで

 
\begin{align}
p(x_{t}|y_{1:t}) \approx \frac{1}{N}\sum_{i=1}^N \delta (x_{t} - x^{(i)}_{t|t})
\end{align}

となる \{x^{(i)}_{t|t}\}_{i=1}^N を求める一連の手続きが, EnKFである.

以下では, 添字の  t|t-1 は時刻  t-1 の情報を用いて時刻tの状態を推定する操作(予測)を表し,  t|tは現在の時刻の情報を用いて状態を最適化する操作(フィルタリング)を表す.

またプライム記号\primeは行列の転置を表す.

アンサンブルカルマンフィルタ

  1. 初期状態のアンサンブル \{x_{0|0}^{(i)}\}_{i=1}^N を生成する. t=1とする.

  2. 以下を行う.

    1. 1期先予測

      • システムノイズのアンサンブル \{ v_t^{(i)} \}_{i=1}^N を生成する.

      • iに対して x_{t|t-1}^{(i)} = f_t (x_{t-1|t-1}^{(i)}, v_t^{(i)}) を計算して, 1期先予測アンサンブルを得る.

    2. フィルタ

      • システムノイズのアンサンブル \{ w_t^{(i)} \}_{i=1}^N

      • 分散共分散行列 \hat{V}_{t|t-1} ,  \hat{R}_t ならびに \hat{K}_t を計算する

      • iに対して,  x_{t|t}^{(i)} = x_{t|t-1}^{(i)} + \hat{K}_t (y_t + \check{w}_t^{(i)} - H_t x_{t|t-1}^{(i)}) を計算し, フィルタアンサンブル \{x_{t|t}^{(i)}\}_{i=1}^N を得る.

3: t=Tならば停止. それ以外は t\leftarrow t+1として2に戻る.

以上の手続きで使用した変数の定義を示す.

  • 状態変数ベクトルの偏差
 
\begin{align}
\check{x}_{t|t-1}^{(i)} = x_{t|t-1}^{(i)} - \frac{1}{N} \sum_{j=1}^N x_{t|t-1}^{(j)}
\end{align}
  • 状態変数ベクトルの分散共分散行列
 
\begin{align}
\hat{V}_{t|t-1} = \frac{1}{N-1} \sum_{j=1}^N x_{t|t-1}^{(j)} x_{t|t-1}^{(j)^{\prime}}
\end{align}
  • 観測ノイズの偏差
 
\begin{align}
\check{w}_{t}^{(i)} = w_{t}^{(i)} - \frac{1}{N} \sum_{j=1}^N w_{t}^{(j)}
\end{align}
  • 観測ノイズの分散共分散行列
 
\begin{align}
\hat{R}_{t} = \frac{1}{N-1} \sum_{j=1}^N w_{t}^{(j)} w_{t}^{(j)^{\prime}}
\end{align}
  • カルマンゲイン
 
\begin{align}
\hat{K}_{t} = \hat{V}_{t|t-1} H_t^{\prime} (H_t \hat{V}_{t|t-1} H_t^{\prime} + \hat{R}_t)^{-1}
\end{align}

適用

非線形の減衰を受けた振動子の運動を記述するVan der Pol 方程式

 
\begin{align}
\frac{d^{2} x}{d t^{2}}-\mu\left(1-x^{2}\right) \frac{d x}{d t}+x=0
\end{align}

にEnKFを適用する. 減衰の強さを制御するパラメータである \muを未知数とし, このパラメータの真値を同定することで, 真の力学モデルを同定する問題を解く.

下図に適用した結果を示す. 赤線が真値であり, 緑色の三角が観測ノイズの乗った観測である. EnKFによる推定結果である青線が、次第に真値に近づいていく様子が見られる. f:id:dosuex:20210518015005p:plain

こちらは減衰パラメータの時間ステップごとの推定結果を示している. 青線がそれぞれのアンサンブル粒子(状態変数ベクトル)の保持する減衰パラメータに相当する. 全ての粒子が次第に正しい減衰パラメータの値に収束していく様子が見られる. f:id:dosuex:20210518224912p:plain

参考書籍