はじめに
データ同化は、数値モデルによる予報と観測データを組み合わせ、システムの状態を推定する手法です。数値天気予報をはじめ、さまざまな分野で利用されています。
今回は、カオス的な振る舞いを示す Lorenz 63モデルに 粒子フィルタ(particle filter; PF)を適用し、状態推定の過程をインタラクティブに可視化するシミュレーターを作成しました。
ブラウザ上で動作し、粒子の時間発展、観測による重みの更新、リサンプリングの様子を確認できます。
シミュレーター
Lorenz 63 — Particle Filter | climatecontour
JavaScriptで実装しており、Pythonなどの実行環境を用意する必要はありません。
1. Lorenz 63モデル
Lorenz 63モデルは、次の3変数の常微分方程式で表されます。
\begin{aligned}
\frac{dx}{dt} &= \sigma(y-x),\\
\frac{dy}{dt} &= x(\rho-z)-y,\\
\frac{dz}{dt} &= xy-\beta z.
\end{aligned}
今回は、標準的なパラメータ
\sigma=10,\quad \rho=28,\quad \beta=\frac{8}{3}
を用いています。
この条件では、軌道が二つのローブの間を不規則に遷移する、カオス的な振る舞いを示します。
数値積分には4次のRunge–Kutta法（RK4法）を使用し、時間刻みは0.01としています。
2. 粒子フィルタによる状態推定
粒子フィルタ（particle filter; PF）は、逐次モンテカルロ法（Sequential Monte Carlo; SMC）に基づき、状態の確率分布を多数の重み付き粒子で近似する状態推定手法です。
今回のシミュレーターでは、ブートストラップ粒子フィルタ（bootstrap particle filter）を実装しました。これは逐次重点サンプリング・リサンプリング（sequential importance resampling; SIR）に基づく手法です。予報に状態遷移モデルを使い、観測尤度で粒子の重みを更新します。リサンプリングは毎回ではなく、ESSに応じて実行します。
予報ステップ（粒子の時間発展）
各粒子をLorenz 63モデルで時間積分します。
各 予報粒子 には、独立な 加法性白色ガウス雑音（additive Gaussian process noise） を加えます。一方、シミュレーション上の真値（以下、真値）は雑音を加えず、決定論的なLorenz 63モデルで生成します。このため、真値を生成するモデルとフィルタが仮定する確率的な状態遷移モデルは、厳密には一致しません。
\mathbf{x}_{n+1}^{f,(i)}
=
M_{\Delta t}(\mathbf{x}_n^{(i)})
+
q\sqrt{\Delta t}\,\boldsymbol{\xi}_n^{(i)}
ここで、$M_{\Delta t}$ はLorenz 63モデルの時間積分写像、$\boldsymbol{\xi}_n^{(i)}\sim\mathcal{N}(\mathbf{0},I_3)$ は粒子ごと・時間ステップごとに独立に生成する標準正規乱数ベクトル、$q$ は雑音振幅です。
観測の生成
観測値は、真値に 観測演算子（observation operator） を適用し、ガウス分布に従う観測誤差を加えて生成します。観測誤差は観測時刻・観測成分ごとに独立とし、その標準偏差を同じ値に設定しています。
\mathbf{y}_k=H(\mathbf{x}_k^{\mathrm{true}})+\boldsymbol{\epsilon}_k
標準設定では $H(x,y,z)=(x,y)^\mathsf{T}$ とし、$x,y$ の2成分を観測します。観測演算子は切り替え可能で、3変数から1成分を返す線形観測演算子と非線形観測演算子も選択できます。非線形観測では、観測空間での平均 $E[H(X)]$ と、状態空間での平均を観測空間へ写像した $H(E[X])$ は一般に一致しません。
重みの更新とリサンプリング
観測が得られると、各予報粒子について 観測尤度（observation likelihood） を計算し、重要度重み（importance weights） を更新して正規化します。次式の左辺は正規化後の重みです（比例定数は重みの総和が1になるよう定めます）。
w_k^{(i)}
\propto
w_{k-1}^{(i)}
p(\mathbf{y}_k\mid\mathbf{x}_k^{f,(i)})
更新後の正規化重みから、**有効サンプルサイズ（effective sample size; ESS）**の推定量
N_{\mathrm{eff}}
=
\frac{1}{\sum_{i=1}^{N}(w_k^{(i)})^2}
を計算します。ESSが閾値（この実装では $0.55N$）を下回った場合、**系統的リサンプリング（systematic resampling）**を実行します。ESSはリサンプリング前の正規化重みから求め、解析推定値もリサンプリング前の重み付き粒子群から計算しています。
3. シミュレーターで可視化できること
3次元状態空間（位相空間）
3次元状態空間には、真値、予報・解析粒子群、粒子の**重み付き平均（解析推定値）**を表示しています。解析推定値は、観測を反映して更新された事後分布の平均の近似です。
時間発展に伴って粒子群が広がり、観測による重み更新とリサンプリングを通じて粒子分布が変化する様子を確認できます。
特に、Lorenzアトラクターの二つのローブの間を真値が遷移する際、粒子群がどのように状態を推定するかを観察できます。
状態変数と推定の不確かさ
$x,y,z$ の時間発展を表示し、真値と推定値を比較できます。
各状態成分の粒子分布から計算した 重み付き5–95%分位点区間（各周辺分布の中央90%区間）も表示します。これは3次元状態ベクトル全体の同時90%信用領域ではありません。また、区間はフィルタが仮定したモデルに基づく推定分布を表します。
解析誤差とESS
真値が既知の数値実験なので、解析推定値と真値との差を直接評価できます。画面の解析RMSEは、各観測時刻で3状態成分 $x,y,z$ について計算した二乗平均平方根誤差（root mean square error; RMSE）です。時間平均や複数の独立試行にわたるRMSEではありません。
\mathrm{RMSE}_k^{a}=\sqrt{\frac{1}{3}\sum_{j=1}^{3}\left(\hat{x}_{k,j}^{a}-x_{k,j}^{\mathrm{true}}\right)^2}
また、ESSの時間変化を確認することで、粒子の重みの集中とリサンプリングの関係を観察できます。
4. パラメータを変更して実験
以下は本記事作成時の標準設定です。
| 項目 | 標準設定 |
|---|---|
| 粒子数 | 1200 |
| 真値の初期状態 | (−8, 8, 27) |
| 観測成分 | x, y |
| 観測誤差の標準偏差 | 5.0 |
| 観測時間間隔 | 0.50 |
| 予報粒子の過程雑音振幅 | 1.8 |
設定画面では、以下のパラメータを変更できます。
- 粒子数
- 観測演算子
- 観測誤差の標準偏差
- 予報粒子に加える過程雑音の振幅
- 観測時間間隔
- 再生速度
例えば、観測誤差を大きくしたり、観測間隔を長くしたりすると、状態推定が難しくなる場合があります。実際の変化は観測演算子や初期粒子分布などにも依存します。
一方で、粒子数や過程雑音の設定も、粒子分布の広がりや推定精度に影響します。
これらのパラメータを変更することで、粒子フィルタの挙動を比較できます。
5. 実装について
シミュレーターはHTMLとJavaScriptを使用して実装しました。
Lorenz 63モデルの時間積分、重要度重みの更新、ESSの計算、系統的リサンプリングなどの計算をブラウザ内で実行しています。
可視化にはCanvasを使用し、3次元状態空間の視点をマウス操作で変更できるようにしています。
おわりに
Lorenz 63モデルに粒子フィルタを適用し、観測によって状態推定が更新される様子を可視化するシミュレーターを作成しました。
数式だけでは把握しにくい、粒子群の時間発展やリサンプリングの挙動を観察できるようにしています。
データ同化や粒子フィルタの学習などに利用していただければ幸いです。
シミュレーターはこちら：
Lorenz 63 — Particle Filter | climatecontour
また、同じカオスモデルを用いたブラウザゲームを作成しました。
市民の予報官｜4D-Varデータ同化で挑む気象予報ゲーム
参考文献
- 矢野浩一（2014）「粒子フィルタの基礎と応用：フィルタ・平滑化・パラメータ推定」『日本統計学会誌』44(1), 189–216. 本文（J-STAGE）
- van Leeuwen, P. J., et al. (2019). “Particle filters for high-dimensional geoscience applications: A review.” Quarterly Journal of the Royal Meteorological Society, 145, 2335–2365. https://doi.org/10.1002/qj.3551
- Douc, R., Cappé, O., and Moulines, E. (2005). “Comparison of Resampling Schemes for Particle Filtering.” Proceedings of the 4th International Symposium on Image and Signal Processing and Analysis (ISPA 2005), 64–69. https://doi.org/10.1109/ISPA.2005.195385
- Schenk, N., Potthast, R., and Rojahn, A. (2022). “On Two Localized Particle Filter Methods for Lorenz 1963 and 1996 Models.” Frontiers in Applied Mathematics and Statistics. https://doi.org/10.3389/fams.2022.920186