アンサンブルカルマンフィルタ（EnKF）の理論とPython実装：高次元データ同化のためのモンテカルロ共分散近似

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

粒子フィルタのPython実装では、非線形・非ガウスな状態空間モデルに対応できる一方で、「高次元問題では粒子数が指数的に増加する次元の呪い」に直面することを述べました。Unscented Kalman Filter（UKF）はシグマポイントという決定論的な代表点（状態次元 \(n\) に対して \(2n+1\) 個）で非線形性を扱うことでこの問題を緩和しますが、シグマポイントの数は次元 \(n\) に比例して増加し、共分散更新には \(O(n^3)\) のコレスキー分解が必要です。

アンサンブルカルマンフィルタ（Ensemble Kalman Filter, EnKF）は、この2つの手法の中間に位置する手法です。粒子フィルタと同じくモンテカルロ的な「アンサンブル」（サンプル集団）で確率分布を表現しますが、重み付けや重点サンプリングは行わず、アンサンブルの標本平均・標本共分散を使って線形カルマンフィルタと同じ更新式を適用します。アンサンブルサイズ \(N\) はユーザーが自由に選べるパラメータであり、UKFのシグマポイントのように次元 \(n\) に縛られません。この性質により、EnKFは気象予報・海洋データ同化のように状態次元が \(10^6\) 〜\(10^9\) に達するような超高次元問題でも実用可能な、標準的な**データ同化（data assimilation）**手法として使われています。

EnKFは1994年にGeir Evensenによって提案され、その後、粒子フィルタと同様にモンテカルロサンプリングを共分散近似に用いる点、UKFや非線形カルマンスムーザと同様にガウス近似のカルマン更新式を用いる点で、両者の折衷案として発展してきました。

状態空間モデル

一般の非線形状態空間モデルを考えます。

\(f\) は状態遷移関数、\(h\) は観測関数、\(w_{k-1}, v_k\) はそれぞれ平均ゼロ・共分散 \(Q, R\) のガウス過程ノイズ・観測ノイズです。粒子フィルタとは異なり、EnKFはノイズがガウス分布であることを前提とします（この点はUKFやカルマンフィルタと同様です）。\(f, h\) は非線形でも構いません。

アルゴリズム導出

アンサンブルによる予測ステップ（Forecast）

時刻 \(k-1\) の解析（更新後）アンサンブル \(\{x_{k-1}^{a,(i)}\}_{i=1}^{N}\) が与えられているとします。各メンバーを状態遷移モデルで独立に伝播させます。

アンサンブル平均と標本共分散（不偏推定量として \(N-1\) で除す）を計算します。

式(4)-(5)は、線形カルマンフィルタの予測共分散 \(P_{k|k-1} = A P_{k-1|k-1} A^T + Q\) （カルマンフィルタの記事の式4）を、アンサンブルのモンテカルロ標本統計量で近似的に置き換えたものです。\(f\) が非線形であっても、アンサンブルメンバーをそのまま \(f\) に通すだけでよく、EKFのようなヤコビ行列や、UKFのようなシグマポイントの生成は不要です。

摂動観測によるアンサンブル更新（Analysis）

EnKFの分析ステップの要点は、観測 \(y_k\) をそのまま全メンバーに適用するのではなく、**各メンバーに独立な観測ノイズを加えた「摂動観測（perturbed observations）」**を用いる点です。これはBurgers, van Leeuwen & Evensen (1998)が示した、更新後アンサンブルの共分散が真の解析共分散と一致するために必要な操作です。

各予測メンバーを観測関数で変換します。

状態と観測の交差共分散、観測の共分散（イノベーション共分散）をアンサンブル統計量として計算します。

アンサンブルカルマンゲイン：

各メンバーを個別に更新します。

式(9)-(12)は、カルマンフィルタの更新式（式5-9）と本質的に同じ構造です。異なるのは、\(H P_{k|k-1} H^T\) や \(H P_{k|k-1}\) を解析的に計算する代わりに、アンサンブルの標本共分散 \(P_k^{zz}, P_k^{xz}\) で近似する点だけです。観測関数 \(h\) が非線形であっても、式(7)でアンサンブルメンバーを \(h\) に通すだけでよく、この点はUKFのシグマポイント変換と同じ発想ですが、シグマポイントが決定論的な \(2n+1\) 点であるのに対し、EnKFのアンサンブルは確率的な \(N\) 点であり、\(N\) は次元 \(n\) と独立に選べます。

Python実装

式(3)-(12)をそのままnumpyで実装します。

import numpy as np

def enkf_step(ensemble, f, h, Q, R, z, rng, inflation=1.0):
    """摂動観測（perturbed observations）によるEnKFの1ステップ（式3-12）

    ensemble : (N, n) 各行がアンサンブルメンバー
    f, h     : 状態遷移関数・観測関数（1メンバーを受け取り1メンバーを返す）
    Q, R     : プロセスノイズ・観測ノイズの共分散
    z        : 観測ベクトル (m,)
    inflation: 乗法的共分散膨張の係数（1.0で膨張なし）
    """
    N, n = ensemble.shape
    m = R.shape[0]

    # ---- 予測ステップ（式3-5） ----
    Xf = np.array([f(ensemble[i]) for i in range(N)])
    Xf += rng.multivariate_normal(np.zeros(n), Q, size=N)
    xf_mean = Xf.mean(axis=0)

    # 乗法的共分散膨張：アンサンブル平均のまわりで偏差を拡大する
    if inflation != 1.0:
        Xf = xf_mean + inflation * (Xf - xf_mean)

    # ---- 分析ステップ（式6-12） ----
    Zf = np.array([h(Xf[i]) for i in range(N)])
    zf_mean = Zf.mean(axis=0)

    Ex = Xf - xf_mean
    Ez = Zf - zf_mean
    Pxz = (Ex.T @ Ez) / (N - 1)          # 式9
    Pzz = (Ez.T @ Ez) / (N - 1) + R      # 式10

    K = Pxz @ np.linalg.inv(Pzz)         # 式11（アンサンブルカルマンゲイン）

    Y = z + rng.multivariate_normal(np.zeros(m), R, size=N)  # 式6：摂動観測
    Xa = Xf + (Y - Zf) @ K.T             # 式12

    return Xa


def run_enkf(x0_ensemble, f, h, Q, R, observations, rng, inflation=1.0):
    """観測系列に対してEnKFを逐次実行し、各時刻のアンサンブル平均を返す"""
    ensemble = x0_ensemble.copy()
    n = x0_ensemble.shape[1]
    means = np.zeros((len(observations), n))
    for k, z in enumerate(observations):
        ensemble = enkf_step(ensemble, f, h, Q, R, z, rng, inflation=inflation)
        means[k] = ensemble.mean(axis=0)
    return means

比較対象として、解析解を与える通常の（線形）カルマンフィルタも実装します（カルマンフィルタの記事と同一実装）。

class KalmanFilter:
    """線形カルマンフィルタ（解析解、EnKFの収束先の基準として使用）"""

    def __init__(self, A, H, Q, R, x0, P0):
        self.A, self.H, self.Q, self.R = A, H, Q, R
        self.x, self.P = x0.copy(), P0.copy()

    def predict(self):
        self.x = self.A @ self.x
        self.P = self.A @ self.P @ self.A.T + self.Q
        return self.x.copy()

    def update(self, z):
        y = z - self.H @ self.x
        S = self.H @ self.P @ self.H.T + self.R
        K = self.P @ self.H.T @ np.linalg.inv(S)
        self.x = self.x + K @ y
        self.P = (np.eye(len(self.x)) - K @ self.H) @ self.P
        return self.x.copy()

実験1：低次元線形ガウス問題でのKF収束検証

EnKFは線形ガウスモデルでは、\(N \to \infty\) の極限で通常のカルマンフィルタの解析解に収束するはずです。これを4次元（位置・速度 × 2軸）の等速直線運動モデルで検証します。状態遷移行列 \(A\) 、観測行列 \(H\) （位置のみ観測）はカルマンフィルタの記事と同じ構造の2次元版です。

観測系列を1つ固定して生成し、通常のKF（決定論的な解析解）を1回計算した上で、EnKFをアンサンブルサイズ \(N \in \{10, 20, 50, 200, 1000\}\) ごとに40通りの乱数シードで実行し、EnKF推定とKF解析解とのRMSEを平均しました。

Ns = [10, 20, 50, 200, 1000]
n_trials = 40
results = {}
for N in Ns:
    rmses = []
    for trial in range(n_trials):
        rng = np.random.default_rng(1000 * N + trial)
        ens0 = x0 + rng.multivariate_normal(np.zeros(n), P0, size=N)
        enkf_means = run_enkf(ens0, f, h, Q, R, observations, rng)
        rmse = np.sqrt(np.mean((enkf_means[:, :2] - kf_means[:, :2]) ** 2))
        rmses.append(rmse)
    results[N] = (np.mean(rmses), np.std(rmses))
    print(f"N={N}: EnKF-vs-KF position RMSE = {np.mean(rmses):.5f}")

# RMSE ~ C / sqrt(N) をフィット
C = np.mean([results[N][0] * np.sqrt(N) for N in Ns])

実行結果は以下の通りです（位置成分のRMSE、40試行平均）。

フィットした係数は \(C = 0.9065\) です。\(N=10\) から \(N=1000\) （100倍）にかけてRMSEは \(0.358 \to 0.025\) （約 14.1分の1）に減少しており、これは \(\sqrt{100} = 10\) に近いオーダーで、サンプリング誤差が \(O(1/\sqrt{N})\) で減衰するというモンテカルロ推定の理論的性質と整合しています（\(N\) を100倍しても誤差は10分の1にしかならない、収束の遅さもこの実験から確認できます）。なお通常のKF（解析解）自体の真の状態に対するRMSEは0.4990で、EnKFの \(N=1000\) でのKFとの残差（0.0254）はこれよりずっと小さく、\(N\) が大きいアンサンブルは実質的にKFと見分けがつかないレベルまで収束していることが分かります。

実験2：Lorenz-96モデルによる高次元データ同化

EnKFの実用上の強みは高次元問題での挙動にあります。気象データ同化のテストベッドとして標準的なLorenz-96モデル（40変数、カオス的挙動）で検証します。

\(F=8\) （カオス的レジーム）、\(n=40\) とし、4次のルンゲ・クッタ法（\(\Delta t = 0.05\) ）で離散化します。観測は40次元のうち偶数インデックスの20成分のみ（部分観測、\(R = I_{20}\) ）とし、真値は \(Q = 0.01 I_{40}\) のプロセスノイズを加えて生成しました。

def lorenz96_step(x, dt=0.05, F=8.0):
    """4次のルンゲ・クッタ法によるLorenz-96モデルの1ステップ積分（式13）"""
    def rhs(x):
        return (np.roll(x, -1) - np.roll(x, 2)) * np.roll(x, 1) - x + F
    k1 = rhs(x)
    k2 = rhs(x + 0.5 * dt * k1)
    k3 = rhs(x + 0.5 * dt * k2)
    k4 = rhs(x + dt * k3)
    return x + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)

EnKF（アンサンブルサイズ \(N=40\) 、乗法的共分散膨張1.10）、UKF（\(n=40\) なので \(2n+1=81\) 個のシグマポイント）、そして粒子フィルタ（粒子フィルタの記事と同じブートストラップ型、系統的リサンプリング）を同一の観測系列（200ステップ）に適用しました。

実測結果（全200ステップに対する全状態成分のRMSE）：

（参考：真の軌道の気候学的標準偏差は3.6792。RMSEがこれを上回る手法は「観測を使わず気候値を答えるより悪い」ことを意味します。）

EnKF（\(N=40\) ）は開始直後の過渡応答を除けば良好に追跡でき、終盤100ステップのRMSEは0.6392まで下がり、UKF（RMSE 0.8782）と同等以上の精度を達成しました。しかも必要なサンプル数はUKFの \(81\) 個のシグマポイントに対し \(40\) 個（約半分）であり、実行時間もUKFの \(1.30\) 秒に対しEnKFは \(0.58\) 秒（約2.3倍高速）でした。一方、同じ乗法的膨張を適用してもアンサンブルサイズを \(N=20\) まで減らすとRMSEは \(4.33\) まで悪化し、気候学的標準偏差（3.68）に近づきます。これは、共分散膨張だけでは補いきれない「アンサンブルの標本ランク不足によるスプリアス相関」の影響であり、実運用のEnKFでは本記事では扱わない**共分散局所化（covariance localization）**を併用してこの問題に対処します。

対照的に、粒子フィルタは \(N=20\) でも \(N=1000\) でもほぼ気候学的標準偏差と同水準のRMSE（5.29、5.24）にとどまり、実質的に追跡に失敗しています。有効サンプルサイズ（ESS）を計測すると、\(N=20\) では平均ESS \(3.48\) （ほぼ全ての重みが数個の粒子に集中）、\(N=1000\) に増やしても平均ESS \(20.12\) にとどまりました。40次元の観測尤度に対して重みの分散が急激に大きくなり、粒子のほとんどが無視される「重みの縮退」が起きていることが分かります。これはまさに粒子フィルタの記事で述べた「次元の呪い」の実例であり、EnKFがアンサンブル平均・共分散という要約統計量だけで更新を行う（個々のメンバーの尤度を評価して重み付けしない）ことが、高次元での頑健性の源泉であることが実測からも裏付けられます。

実験3：共分散膨張とフィルタ発散

アンサンブルサイズが小さいと、標本共分散がサンプリング誤差によって過小評価され、フィルタが自身の予測を過信して観測を軽視するようになります。これが進行すると、共分散が実際の誤差より小さいまま推定が真の状態から乖離し続ける**フィルタ発散（filter divergence）**が起こります。この現象を、次元を落とした10変数のLorenz-96モデル（\(n=10\) 、観測5成分）とアンサンブルサイズ \(N=8\) （次元の8割程度、標本ランク不足になりやすい設定）で再現しました。

for k, z in enumerate(observations):
    ensemble = enkf_step(ensemble, f, h, Q, R, z, rng, inflation=inflation)
    # inflation=1.0（膨張なし） と inflation=1.05（乗法的膨張）を比較

膨張なし（inflation=1.0）と、乗法的共分散膨張1.05（式(3)直後にアンサンブル偏差を1.05倍に拡大）を適用した場合の、20ステップおきのRMSEの推移は以下の通りです。

膨張なしの場合、\(k=60\) 付近まではRMSE 0.5前後まで良好に追跡できていたにもかかわらず、\(k=120\) 以降にRMSEが単調に増加し、\(k=180\) では4.875まで悪化しています（この系の気候学的標準偏差は3.7134であり、これを上回るということは「観測なしで気候値を答える」より悪い状態、すなわちフィルタ発散です）。標本ランク不足によるアンサンブル共分散の過小評価が、時間とともに自己強化的に進行し、フィルタが観測を信用しなくなっていく典型的なパターンです。一方、乗法的共分散膨張1.05を適用した場合はRMSEが終始2.0以下にとどまり、発散は見られませんでした。

全200ステップでの平均を見ても、**膨張なしの全体RMSEは2.313、膨張1.05では1.658（約28%改善）、終盤50ステップに限れば3.686 → 1.783（約52%改善）**と、共分散膨張がフィルタ発散を明確に抑制していることが数値で確認できます。ただし膨張率を上げすぎても改善するとは限らず（本実験では1.02〜1.05付近が最良で、1.10ではむしろ発散が再発しました）、膨張率はチューニングが必要なハイパーパラメータです。

フィルタ手法の使い分け

これまで本ブログで扱ってきたカルマン系フィルタとEnKFを、線形性・次元・計算コストの観点で比較します。

EnKFは「共分散の伝播をアンサンブルのモンテカルロ近似に置き換える」という点で粒子フィルタと発想を共有しますが、更新式自体は線形カルマンフィルタと同じくガウス近似に基づきます。この折衷によって、多峰性への対応力は失う代わりに、UKFのようにシグマポイント数が次元に比例して増える制約や、粒子フィルタの次元の呪いを回避し、状態次元が数百万に達する数値天気予報のような超高次元問題でも扱えるようになります。逆に、低次元かつ強い非線形性・非ガウス性を持つ問題では、本記事の実験1・2が示す通りEnKFにもサンプリング誤差が存在するため、UKFや粒子フィルタの方が有利な場面もあります。

おすすめ書籍

カルマンフィルタの導出から実装例までを日本語で丁寧に解説した定番書です。EnKFの更新式（式9-12）が線形カルマンフィルタの更新式（式5-9）とどう対応しているかを理解する土台になります。

※ 上記は Amazon アソシエイトのリンクです。

関連記事

• カルマンフィルタの理論とPython実装 - EnKFの更新式が近似する、線形ガウス系での解析解を解説しています。
• Unscented Kalman Filter（UKF）の理論とPython実装 - シグマポイントという決定論的な代表点で非線形性を扱う手法。EnKFの確率的なアンサンブルとの対比になります。
• 粒子フィルタのPython実装：リサンプリング手法の比較 - 同じくモンテカルロサンプリングを使いますが、重み付けとリサンプリングにより任意の分布に対応する一方、高次元では次元の呪いに直面します。本記事の実験2でその挙動を直接比較しています。
• 非線形カルマンスムーザ（Extended RTS Smoother・Unscented RTS Smoother） - EKF/UKFをスムージングに拡張した手法。EnKFにも同様のアンサンブルスムーザ（EnKS）が存在します。
• 拡張カルマンフィルタ（EKF）の理論とPython実装 - ヤコビ行列による線形化アプローチ。EnKFはヤコビ行列を必要としません。
• カルマンスムーザの比較 - カルマンフィルタ系のスムージング手法群を横断比較しています。
• MCMC入門：メトロポリス・ヘイスティングス法とギブスサンプリング - EnKFや粒子フィルタと同じくモンテカルロサンプリングを推定に用いる手法です。
• 機械学習による時系列予測・分類・異常検知ハブ - 線形カルマン・EKF・UKF・粒子フィルタを含む系列モデル群のハブ記事。EnKFは超高次元データ同化という追加の軸を提供します。

参考文献

• Evensen, G. (1994). “Sequential Data Assimilation with a Nonlinear Quasi-Geostrophic Model Using Monte Carlo Methods to Forecast Error Statistics.” Journal of Geophysical Research, 99(C5), 10143-10162.
• Burgers, G., van Leeuwen, P. J., & Evensen, G. (1998). “Analysis Scheme in the Ensemble Kalman Filter.” Monthly Weather Review, 126(6), 1719-1724.
• Evensen, G. (2009). Data Assimilation: The Ensemble Kalman Filter (2nd ed.). Springer.
• Lorenz, E. N. (1996). “Predictability: A Problem Partly Solved.” Proceedings of the Seminar on Predictability, ECMWF.
• Houtekamer, P. L., & Mitchell, H. L. (1998). “Data Assimilation Using an Ensemble Kalman Filter Technique.” Monthly Weather Review, 126(3), 796-811.