粒子フィルタと固定ラグスムーザ

概要

このコラムでは、非線形・非ガウスな動的システムにおける状態推定手法として広く用いられる粒子フィルタ(Particle Filter, PF)と、 推定値の時間的な一貫性を高めるための固定ラグ平滑化(Fixed-lag Smoothing)を組み合わせた実装を紹介します。

ベイズの定理

PFは、ベイズフィルタのモンテカルロ近似です。ベイズフィルタは以下の2ステップを通じて、逐次的に状態を更新します:

1. 予測ステップ:

\[ p(x_t \mid z_{1:t-1}) = \int p(x_t \mid x_{t-1}, u_t) \cdot p(x_{t-1} \mid z_{1:t-1})\,dx_{t-1} \] これは過去の状態 \(x_{t-1}\) の分布に対して、システム入力 \(u_t\) に基づく予測分布を畳み込む操作です。 ここで,

\( z_{1:t-1} = \{ z_1,z_2,...,z_{t-1}\}\)
は、時刻1から\(t-1\)までに得られた観測の履歴を意味します。 この式は、「過去の観測履歴をもとに、時刻 \(t\) における状態 \(x_t\) がどれくらい確からしいか」を計算するために、 前時刻 \(x_{t-1}\) の確率と、そこから現在への遷移確率 \(p(x_t \mid x_{t-1}, u_t)\) を組み合わせて、すべての \(x_{t-1}\) にわたって統合(積分)していることを意味します。

2. 更新ステップ:

\[ p(x_t \mid z_{1:t}) \propto p(z_t \mid x_t) \cdot p(x_t \mid z_{1:t-1}) \] これは観測値 \(z_t\) に基づき、予測分布を重み付けによって修正する操作です。

粒子フィルタ

粒子フィルタは、状態分布を明示的な関数ではなく、粒子(サンプル)と重みの集合 \(\{x_t^{(i)}, w_t^{(i)}\}\) により表現します。 これにより非線形・非ガウスの状態空間でも近似的な推論が可能となります。実装上は以下の手順で構成されます:

  1. 各粒子 \(x_{t-1}^{(i)}\) に対して、状態遷移モデルに基づいて新しい粒子 \(x_t^{(i)}\) をサンプリング(propagate()
  2. 観測値 \(z_t\) に対する尤度 \(p(z_t \mid x_t^{(i)})\) を計算し、重み \(w_t^{(i)}\) を更新(computeLikelihood()
  3. 重みの正規化後、サンプリング(resample())により低確率粒子を除去

これにより、ベイズ更新 \[ p(x_t^{(i)} \mid z_{1:t}) \propto p(z_t \mid x_t^{(i)}) \cdot p(x_t^{(i)} \mid z_{1:t-1}) \] を粒子集合として近似することが可能になります。

固定ラグスムーザ

粒子フィルタは時間方向に向かう逐次処理なので,過去の状態に対する修正能力は持ちません。 これに対し、固定ラグスムーザは、\(t-L\) から \(t\) までの履歴を保持し、 未来情報を用いて過去の粒子分布を再評価することができます。

具体的には以下のような後退処理を行います:

  1. 時刻 \(k+1\) のスムーズ粒子の平均 \(\mu_{k+1}\) と分散(あるいは共分散行列)\(\Sigma_{k+1}\) を評価
  2. 時刻 \(k\) の粒子 \(x_k^{(i)}\) を1ステップ予測し、\(x_{k+1}\) における一致度(マハラノビス距離など)を評価
  3. その一致度から粒子の重みを更新し、リサンプリング

この手法により、粒子は未来の情報に整合的なものが選ばれ、時系列的な一貫性を向上することができます。

特徴的なのは、スムーザで行われる予測および修正の処理が、フィルタで用いる操作と本質的に共通している点です。 状態遷移や尤度評価の操作が共通であることは、実装と理論の両面で大きな利点をもたらします。

このように、共通のロジックに基づいて構築されたスムーザは、コードとしても理論としても自然で扱いやすく、 将来的な応用にも強い設計となっていると考えます。

デモ内容

シミュレーションでは、\(z\) 方向に単振動する対象を想定し、 真の加速度からIMU観測を生成し、 GNSSからは一定間隔で位置を観測します。 これに対し粒子フィルタにより逐次推定を行い、さらに固定ラグスムーザで後方の状態推定を改善します。

スムーザでは、状態空間 \([z, \dot{z}]\) 全体を使ったもの(デフォルト)と、 位置 \(z\) のみを用いた簡易的なものが実装されている. 状態空間全体を使った手法では,尤度評価において位置成分と速度成分の重み係数 \(w_z, w_\dot{z} \) を設定可能としており(ハードコーディング)、 粒子のばらつきに応じた柔軟なスムージングを可能にしています。

このコラムを通じて,ベイズフィルタに基づく推定の全体像と、粒子ベースによる実装の具現化を確認してもらえることを期待します.

デモ

プログラムコード

粒子フィルタおよびスムーザ部分の実装コード(Javascript)を以下に示します. htmlを含むコード全体は,ここからダウンロードできます.

Fuijta@PARI(2025)