データ同化

データ同化とは

---

データ同化とは、確率・統計の考え方を用いて、数値シミュレーションと実験データ・観測データを統合することにより、 シミュレーションでより良い予測を行うための手法です。

データ同化は、元々は気象学や海洋学の分野で用いられていましたが、近年では数値流体解析でも活用が進んでいます。

例えば、Katoら¹は角柱周りの流れに関する風洞実験に対してデータ同化を適用しており、 角柱の表面圧力の実験データをシミュレーションに同化させることで、 図1のようにシミュレーションで得られる流速分布が実験値に近づいたと報告しています。

図1 Katoらによる角柱を用いた風洞実験にデータ同化を適用した研究¹

データ同化により、このようにシミュレーション精度を向上させることが可能である他、未知パラメータの値を推定することもできます。

本記事では、代表的なデータ同化手法の一つであるアンサンブルカルマンフィルタ(EnKF)について紹介します。

EnKFは、アルゴリズムの性質上計算コストが大きくなりやすいですが、 シミュレーションのソルバー部分のコードをほとんど書き換える必要がないため、 実装が容易であるというメリットがあります。

EnKFの理論を説明するにあたっては、文献²を参考にしています。

ベイズ統計

---

本記事では、ベイズ統計の枠組みからEnKFの理論を解説します。

本節では、そのために必要な知識について簡単に説明します。 $P(\cdot)$は確率、$p(\cdot)$は確率密度を表すとします。

確率変数・確率分布

---

確率的に値が変動する変数を確率変数といいます。

例えば、サイコロの目は、1, 2, 3, 4, 5, 6のいずれかの値をランダムにとるので確率変数です。

確率変数には、サイコロの目のように離散的な値をとる離散型確率変数と、連続的な値をとる連続型確率変数があります。

平らな面の上に立てた針が倒れたときに鉛直線となす角は、$0 \sim 360$の実数値をランダムにとるので連続型確率変数です。

確率変数の取りえる値と、その値をとる確率の対応関係を確率分布といいます。 例えば、サイコロの目の確率分布は表1のようになります。

• 表1　サイコロの目の確率分布

  目	1	2	3	4	5	6
  確率	1/6	1/6	1/6	1/6	1/6	1/6

連続型確率変数では、取りえる値の個数は無限なので、ある1つの値を取る確率は0としないと、（確率の総和）= 1が成り立たなくなります。

値ではなく、値の範囲に対してなら確率を考えることができます。

次式のように、区間$[a,b]$で積分すると$a<X<b$となる確率$P(a<X<b)$が得られる関数$p(x)$を確率密度関数と呼びます。

$$ P(a\leq X \leq b) = \int_a^b p(x) dx \tag{1} $$

確率の性質から、$p(x)$は次式を満たします。

$$ p(x) \geq 0 \tag{2} $$

$$ \int_{-\infty}^\infty p(x) dx = 1 \tag{3} $$

確率密度は、その値の相対的な生じやすさを表します。

連続型確率変数がある1つの値$a$しかとらない場合は、確率密度関数はデルタ関数$\delta(x-a)$で与えられます。 デルタ関数の定義は以下の通りです。

$$ \delta (x-1) = \{ \begin{array}{cc} 0 & (x \not= a)\\ \infty & (x = a) \end{array} \tag{4} $$

$$ \int_{-\infty}^\infty \delta(x-a) dx = 1 \tag{5} $$

同時確率・条件付き確率

---

2つの確率変数$X,Y$について、$X=i$と$Y=j$が同時に成り立つ確率を同時確率といい、$P(X=i,Y=j)$と書きます。

また、$Y=j$という条件のもとでの$X=i$となる確率を条件付き確率といい、$P(X=i│Y=j)$と書きます。

定義より、次式が成り立ちます。

$$ P(X=i | Y=j) = \frac{P(X=i, Y=j)}{P(Y=j)} \tag{6} $$

例えば、生徒数$100$人の学校で、各学年の男女別の生徒数が表2のように与えられるとします。

• 表2　生徒数100人の学校の学年、男女別の人数

  	1年	2年	3年	合計
  男子	18	16	28	62
  女子	12	14	12	38
  合計	30	30	40	100

このとき、$P(男子, 1年) = 18/100$, $P(1年) = 30/100$なので$P(男子|1年) = 18/30 = 3/5$となります。

また、$P(男子, 1年) + P(男子, 2年) + P(男子, 3年) = P(男子)$となっていることが分かります。

一般に、2つの確率変数$X, Y$について次式が成り立ちます。

$$ \sum_j P(X=i, Y=j) = P(X=i) \tag{7} $$

$X,Y$が連続型確率変数のときは、$Y=y$という条件のもとでの$X$の確率密度関数を条件付き確率密度関数といい、$p(x|y)$と書きます。

ベイズの定理

---

式(6)から、同時確率$P(X=i,Y=j)$は以下のように2通りに表せます。

$$ P(X=i,Y=j)=P(X=i│Y=j)P(Y=j)=P(Y=j│X=i)P(X=i) \tag{8} $$

式(7), (8)から次に示すベイズの定理が導かれます。

$$ \begin{array}{ll} P(X=i│Y=j) &= \frac{P(Y=j│X=i)P(X=i)}{P(Y=j)} \\ &= \frac{P(Y=j│X=i)P(X=i)}{\sum_k P(Y=j│X=k)P(X=k)} \end{array} \tag{9} $$

ベイズの定理では、$Y=j$という確定情報（観測データ）が新たに得られたときに、 それを反映した$X=i$の確率$P(X=i│Y=j)$を、既知の確率$P(X=i)$から計算します。

$P(X=i)$は観測データが得られる前の確率なので、事前確率と呼びます。

$P(Y=j│X=i)$は$X=i$という条件のもとで$Y=j$となる確率ですが、 これは$Y=j$という観測データに対する$X=i$の尤もらしさを表すため、尤度と呼ばれます。

$P(X=i│Y=j)$は、観測データ$Y=j$が得られた後の$X=i$の確率なので、事後確率と呼びます。

ベイズの定理について次の例題を使って説明します。

• 例題: 赤玉と青玉が計10個入った箱から、玉を1個取り出して色を確認し元に戻す試行を10回行うと、 表3に示す結果が得られた。試行の結果から、箱の中にある赤玉の個数をベイズの定理により推定する。

• 表3　試行の結果

  試行の回数	取り出した玉の色
  1	赤
  2	赤
  3	赤
  4	青
  5	青
  6	赤
  7	赤
  8	青
  9	赤
  10	赤

この例題にベイズの定理を適用すると、$X$は箱の中にある赤玉の個数で取りえる値は0, 1, 2, …, 10、 $Y$は取り出した玉の色で取りえる値は赤, 青となります。

観測データが何もない初期では、箱の中にある赤玉の個数は$0 \sim 10$の値を等しく取りえると仮定し、 $X$の事前確率は次式で与えます。

$$ P(X=i) = \frac{1}{11} \hspace{0.5cm}(i = 0, 1, 2,\cdots, 10) \tag{10} $$

1回目の試行についてベイズの定理を実際に計算してみます。

$$P(X=i) = 1/11 \hspace{0.5cm}(仮定より)$$

$$ \begin{array}{ll} P(Y=赤|X=i) &= (箱の中にある赤玉がi個であるときに赤玉を取り出す確率) \\ &= i/10 \end{array} $$

したがって、

$$ P(X=i│Y=赤) = \frac{P(Y=赤│X=i)P(X=i)}{\sum_{k=0}^{10} P(Y=赤│X=k)P(X=k)} $$

$$ = \frac{\frac{i}{10} \cdot \frac{1}{11}}{\sum_{k=0}^{10} \frac{k}{10} \cdot \frac{1}{11}}=\frac{i}{55} \tag{11} $$

これより、ベイズの定理によって赤玉の個数の確率分布は図2のように更新されます。

図2 試行1の結果による赤玉の個数の確率分布の更新

$n(2 \leqq n \leqq 10)$回目の試行についてベイズの定理を適用するときは、$n-1$回目の試行による事後確率を事前確率として用います。

取り出した玉の色の情報をこのように逐次的に反映させていくと、赤玉の個数の確率分布は図3のように更新されていきます。

10回目の試行後の確率分布を見ると、赤玉の個数は7である確率が最も高いことが分かります。 10回の試行によって赤玉が7回、青玉が3回出ていることを考えると、これは妥当な結果と考えられます。

このように、ベイズの定理を用いると、既知の確率分布に観測データの情報を取り入れて更新することができます。

図3 試行1～10の結果による赤玉の個数の確率分布の更新

なお、データ同化では連続型確率変数を扱うため、実際に用いるのは次式に示す連続型のベイズの定理となります。

$$ p(x|y) = \frac{p(y│x)p(x)}{\int_x p(y│x)p(x)dx} \tag{12} $$

ここで、$p(y│x)$は尤度であり、$p(x)$, $p(x|y)$はそれぞれ事前分布、事後分布といいます。

状態ベクトル・状態空間モデル

---

状態ベクトル

---

データ同化では、図4のようにシミュレーションで時間変化を計算するすべての変数をひとつのベクトルに入れて扱います。

このベクトルを状態ベクトルと呼び、$\mathbf{x}_t$と書きます。

状態ベクトルに推定したいパラメータも含めることで、そのパラメータを推定することが可能です。

図4 状態ベクトルの定義

システムモデル

---

状態ベクトルを用いると、数値シミュレーションとは次式のように状態ベクトルを更新する操作であると考えることができます。

$$ \mathbf{x}_{t+1} = f_t(\mathbf{x}_t) \tag{13} $$

ここで、$f_t$は状態ベクトルの時間発展を表す関数で、 大雑把にいうと数値シミュレーションで解く離散化された物理モデル（微分方程式）です。

さて、シミュレーションの解と実際の自然現象の間には、以下の要因で誤差が存在すると考えられます。

• シミュレーションで解く微分方程式は、そもそも自然現象を完璧に表現できるものではない
• 微分方程式を離散化する際に誤差が生じる

データ同化では、次式のように確率変数$\mathbf{v}_t$を加えることで、 シミュレーションの解と実際の現象とのギャップを考慮します。

$$ \mathbf{x}_{t+1} = f_t(\mathbf{x}_t)+\mathbf{v}_t \tag{14} $$

上式をシステムモデル、$\mathbf{v}_t$をシステムノイズといいます。

確率変数$\mathbf{v}_t$を毎ステップ加えるため、状態ベクトル$\mathbf{x}_t$も確率変数であり、 幅を持った確率分布に従うことになります。

これは、通常のシミュレーションでは初期条件を決めれば任意の時刻でシミュレーション結果が一つに定まるのと大きく異なる点です。

通常のシミュレーションとデータ同化で行うシミュレーションを比較すると、図5のようになります。

図5 通常のシミュレーションとデータ同化で行うシミュレーションの比較

数値シミュレーションに確率の考え方を取り入れている例として、 図6に示す台風の進路予測があります。

台風の進路予測は数値シミュレーションによって行われておりますが、 台風の予報円は、台風の中心が70%の確率でその中に入ることを表しています。

図6 台風の進路予測³

観測モデル

---

観測データも一つのベクトル$\mathbf{y}_t$に入れて扱い、観測データ$\mathbf{y}_t$と状態ベクトル$\mathbf{x}_t$の関係を次式で表現します。

$$ \mathbf{y}_t = h_t(\mathbf{x}_t) + \mathbf{w}_t \tag{15} $$

上式を観測モデル、$\mathbf{w}_t$を観測ノイズといいます。

観測ノイズは、観測データとシミュレーション結果の差を表す項で、データの測定に伴う測定誤差に加えて、 物理モデルの不完全性に起因する誤差も含んでいます。 システムモデルと観測モデルを合わせて、状態空間モデルと呼びます。

$h_t$は状態ベクトルを観測データと比較可能な形に変換する演算子で、観測演算子といいます。

観測演算子が線形で行列により表される場合は、その行列を観測行列と呼び、$\mathbf{H}_t$と表記します。

例えば、観測点がシミュレーションの格子点と重なっており、 単に観測点での変数の値を状態ベクトルから抜き出すだけでいい場合は、 $h_t$は線形となります。

逐次データ同化

---

EnKFは、時系列の観測データを1個ずつ同化する逐次データ同化に分類されます。

逐次データ同化では、予測分布$p(\mathbf x_t│\mathbf y_{1:t-1})$とフィルタ分布$p(\mathbf x_t│\mathbf y_{1:t})$を交互に計算します。

ここで、$\mathbf y_{1:t}=${$\mathbf y_1,\mathbf y_2,\cdots,\mathbf y_t$}であり、$\mathbf y_{1:0}$は空集合とします。

予測分布は時刻$t-1$までの観測データ$\mathbf y_{1:t-1}$を反映した$\mathbf x_t$の確率分布、 フィルタ分布は時刻$t$までの観測データ$\mathbf y_{1:t}$を反映した$\mathbf x_t$の確率分布を表します。

フィルタ分布から予測分布を求めることを一期先予測、予測分布からフィルタ分布を求めることをフィルタリングといいます。

図で表すと以下のようになります。

図7 逐次データ同化の基本的な流れ

カルマンフィルタ

---

システムモデルと観測モデルが次式のようにともに線形である場合を考えます。

$$ \mathbf{x}_{t+1} = \mathbf{F}_t\mathbf{x}_t + \mathbf{v}_t \tag{16} $$

$$ \mathbf{y}_{t} = \mathbf{H}_t\mathbf{x}_t + \mathbf{w}_t \tag{17} $$

ここで、$p(\mathbf{v}_t) = N(\mathbf{0}, \mathbf{Q}_t)$, $p(\mathbf{w}_t)=N(\mathbf{0}, \mathbf{R}_t )$とします。

$N(\mathbf \mu, \mathbf \Sigma)$は平均ベクトル$\mathbf \mu$、分散共分散行列$\mathbf \Sigma$の正規分布であり、$\mathbf{0}$は零ベクトルを表します。

分散共分散行列は、対角成分がベクトルの各成分の分散、非対角成分がベクトルの異なる成分同士の共分散である行列で、多次元正規分布の形状を決定します。

このとき、フィルタ分布$p(\mathbf x_{t-1}│\mathbf y_{1:t-1})$が正規分布で $p(\mathbf x_{t-1}│\mathbf y_{1:t-1})=N(\mathbf \mu_{t-1|t-1}, \mathbf V_{t-1|t-1})$なら、 予測分布$p(\mathbf x_t│\mathbf y_{1:t-1})$も正規分布となり、一期先予測の計算式は以下で与えられます。

$$ p(\mathbf x_t│\mathbf y_{1:t-1}) = N(\mathbf \mu_{t|t-1}, \mathbf V_{t|t-1}) \tag{18} $$

$$ \mathbf \mu_{t|t-1} = \mathbf F_t \mathbf \mu_{t-1|t-1} \tag{19} $$

$$ \mathbf V_{t|t-1} = \mathbf F_t \mathbf V_{t-1|t-1} \mathbf F_t^T+\mathbf Q_t \tag{20} $$

また、予測分布$p(\mathbf x_t│\mathbf y_{1:t-1})$が正規分布で $p(\mathbf x_t│\mathbf y_{1:t-1})=N(\mathbf \mu_{t|t-1}, \mathbf V_{t|t-1})$なら、 フィルタ分布$p(\mathbf x_t│\mathbf y_{1:t})$も正規分布となり、フィルタリングの計算式は以下で与えられます。

$$ p(\mathbf x_t│\mathbf y_{1:t}) = N(\mathbf \mu_{t|t}, \mathbf V_{t|t}) \tag{21} $$

$$ \mathbf \mu_{t|t} = \mathbf \mu_{t|t-1} + \mathbf K_t(\mathbf y_t - \mathbf H_t \mathbf \mu_{t|t-1}) \tag{22} $$

$$ \mathbf K_t = \mathbf V_{t|t-1} \mathbf H_t^T(\mathbf H_t \mathbf V_{t|t-1} \mathbf H_t^T + \mathbf R_t)^{-1} \tag{23} $$

$$ \mathbf V_{t|t} = \mathbf V_{t|t-1} - \mathbf K_t\mathbf H_t \mathbf V_{t|t-1} \tag{24} $$

したがって、初期分布$p(\mathbf x_0│\mathbf y_{1:0})=p(\mathbf x_0)$を正規分布とすれば、 上記の2つの式を交互に用いて任意の時刻の予測分布とフィルタ分布を求めることができます。

このアルゴリズムを、カルマンフィルタ(KF)といいます。

フィルタリングの式は、ベイズの定理で(事前分布) = (予測分布), (事後分布) = (フィルタ分布)とすれば得られます。

したがって、フィルタリングとはベイズの定理により観測データを反映してシミュレーション結果の確率分布を更新する操作となります。

このとき、尤度は観測モデルより$p(\mathbf y_t│\mathbf x_t ) = N(\mathbf H_t \mathbf x_t, \mathbf R_t)$とします。

フィルタリングの式の意味を理解するため、状態ベクトルと観測データがともに1次元で$H_t=1$とすると、以下のようになります。

$$ \mu_{t|t}=\mu_{t|t-1}+K_t(y_t-\mu_{t|t-1})=\frac{\sigma_{t|t-1}^2 y_t+\sigma_{obs}^2 \mu_{t|t-1}}{\sigma_{t|t-1}^2+\sigma_{obs}^2} \tag{25} $$

$$ K_t=\frac{\sigma_{t|t-1}^2}{\sigma_{t|t-1}^2+\sigma_{obs}^2} \tag{26} $$

$$ V_{t|t}=\frac{\sigma_{obs}^2}{\sigma_{t|t-1}^2+\sigma_{obs}^2} V_{t|t-1} \tag{27} $$

ここで、$\sigma_{t|t-1}^2$は予測分布$p(x_t│y_{1:t-1} )$の分散、$\sigma_{obs}^2$は観測ノイズの分散を表します。 平均ベクトルの式は、元の平均ベクトル$\mu_{t|t-1}$と観測データ$y_t$の加重平均の形になっており、$y_t$の重みは$\sigma_{t|t-1}^2$, $\mu_{t|t-1}$の重みは$\sigma_{obs}^2$です。

したがって、$\sigma_{t|t-1}^2$が$\sigma_{obs}^2$より大きいほど$\mu_{t|t}$は$y_t$に近づき、 逆に$\sigma_{t|t-1}^2$が$\sigma_{obs}^2$より小さいほど$\mu_{t|t}$は$\mu_{t|t-1}$から変化しないことが分かります。 また、$\sigma_{obs}^2/(\sigma_{t|t-1}^2+\sigma_{obs}^2)\leq 1$であるため、フィルタ分布の分散$V_{t|t}$は予測分布の分散$V_{t|t-1}$より小さくなります。

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

---

カルマンフィルタ(KF)は、システムモデル・観測モデルが線形で、 なおかつシステムノイズ・観測ノイズが正規分布に従う場合のみ適用可能です。

しかしながら、数値シミュレーションのほとんどではシステムモデルが非線形となるため、カルマンフィルタは使えません。

アンサンブルカルマンフィルタ（EnKF）は、以下のようにシステムモデルが非線形な場合にもKFを適用できるようにしたアルゴリズムです。

$$ \mathbf x_{t+1} = f_t(\mathbf x_t) + \mathbf v_t \tag{28} $$

$$ \mathbf y_t = \mathbf H_t \mathbf x_t + \mathbf w_t \tag{29} $$

ここで、KFと同様に、$p(\mathbf v_t) = N(\mathbf 0, \mathbf Q_t)$, $p(\mathbf w_t) = N(\mathbf 0, \mathbf R_t)$とします。

EnKFでは、図8のようにシミュレーション結果の確率分布を多数のサンプルによって近似します。

そのために、条件の異なる多数のシミュレーションを並列に実行します。

状態ベクトルの確率分布を近似するためのサンプルの集合をアンサンブル、その中の個々のサンプルをアンサンブルメンバーといいます。

$N$個のシミュレーションを並列実行するとき、$i (i = 1, 2, 3, \cdots, N)$番目のシミュレーションの状態ベクトルを$\mathbf x_t^{(i)}$と表記します。

図8 多数のシミュレーション結果による確率分布の近似

EnKFでは、以下の方法で一期先予測とフィルタリングを行います。

• (一期先予測)

フィルタ分布のアンサンブル$\{\mathbf x_{t-1|t-1}^{(i)}\}_{i=1}^N$から、 予測分布のアンサンブル$\{x_{t|t-1}^{(i)}\}_{i=1}^N$を計算します。 次式のように、システムモデルを用いて$N$個の状態ベクトルを更新します。

$$ \mathbf x_{t|t-1}^{(i)} = f_t(\mathbf x_{t-1|t-1}^{(i)})+\mathbf v_t^{(i)} \tag{30} $$

ここで、$\mathbf v_t^{(i)}$はシステムノイズのサンプルです。 $N$個のシミュレーションの各々で、ソルバーの計算を1ステップ分行った後、 システムノイズのサンプル$\mathbf v_t^{(i)}$を疑似乱数により生成し、状態ベクトルに加えます。

• (フィルタリング) 予測分布のアンサンブル$\{x_{t|t-1}^{(i)}\}_{i=1}^N$から、 フィルタ分布のアンサンブル$\{x_{t|t}^{(i)}\}_{i=1}^N$を以下の手順で計算します。

[1] 予測分布のアンサンブルから、サンプル分散共分散行列$\hat {\mathbf V}_{t|t-1}$を計算します。

$$\hat{\mathbf V}_{t|t-1}=\frac{1}{N-1}\sum_{j=1}^N \check{\mathbf x}^{(j)}_{t|t-1} \check{\mathbf x}^{(j)T}_{t|t-1} \tag{31}$$

$$\check{\mathbf x}_{t|t-1}^{(i)}=\mathbf x_{t|t-1}^{(i)}-\frac{1}{N}\sum_{j=1}^N\mathbf x_{t|t-1}^{(j)} \tag{32}$$

[2] 疑似乱数により観測ノイズのサンプル$\mathbf w_t^{(i)}$を生成し、 それらの偏差$\check{\mathbf w}_t^{(i)}$と観測ノイズのサンプル分散共分散行列$\hat{\mathbf R}_t$を計算します。

$$ \check{\mathbf w}_t^{(i)} = \mathbf w_t^{(i)}-\frac{1}{N}\sum_{j=1}^N \mathbf w_t^{(j)} \tag{33} $$

$$ \hat{\mathbf R}_t = \frac{1}{N-1}\sum_{j=1}^N \check{\mathbf w}_t^{(j)}\check{\mathbf w}_t^{(j)T} \tag{34} $$

[3] カルマンゲイン$\hat{\mathbf K}_t$を計算します。

$$ \hat{\mathbf K}_t = \hat{\mathbf V}_{t|t-1} \mathbf H_t^T(\mathbf H_t^T \hat{\mathbf V}_{t|t-1}\mathbf H_t^T + \hat{\mathbf R}_t)^{-1} \tag{35} $$

観測データ$\mathbf y_t$の次元がアンサンブルメンバー数より大きいとき、 $\mathbf H_t^T \hat{\mathbf V}_{t|t-1}\mathbf H_t^T + \hat{\mathbf R}_t$の逆行列は存在しないため、Moore-Penroseの疑似逆行列を用います。

[4] 次式のように、KFのフィルタリングにおける平均ベクトルの更新式(式(22))を各アンサンブルメンバーに適用することで、 フィルタ分布のアンサンブルを求めます。

$$ \mathbf x_{t|t}^{(i)}=\mathbf x_{t|t-1}^{(i)} + \hat{\mathbf K}_t(\mathbf y_t + \check{\mathbf w}_t^{(i)}-\mathbf H_t \mathbf x_{t|t-1}^{(i)}) \tag{36} $$

システムモデルと観測モデルがともに線形で、 システムノイズと観測ノイズがともに正規分布に従うとき、 上式で求まるフィルタ分布のアンサンブル$\{\mathbf x_{t|t}^{(i)}\}_{i=1}^N$のサンプル平均とサンプル分散共分散行列は、 KFのフィルタリング後の平均ベクトルと分散共分散行列に近似的に一致します。

上式でKFのフィルタリングの式には無かった項$\check{\mathbf w}_t^{(i)}$を$\mathbf y_t$に加えていますが、 これは擾乱付き観測と呼ばれます。

KFのフィルタリングの式をそのまま用いると、フィルタ分布のサンプル分散共分散行列がKFの分散共分散行列より小さくなってしまうため、擾乱付き観測を用いる必要があります。

上記の説明は、一期先予測とフィルタリングを交互に行う前提で書いていますが、 観測データが無い時刻ではフィルタリングをスキップすることも可能です。

EnKFのアルゴリズムをまとめると、図9のようになります。

観測データが存在する時刻になるまで一期先予測を行い、時間ステップを進めます。

観測データが存在する時刻に到達したら、フィルタリングによってシミュレーション結果に観測データを同化し、シミュレーション結果を補正します。

すべての実験・観測データの同化が完了するまで、これらの処理を繰り返し行います。

図9 EnKFのアルゴリズム

• cfd
• データ同化
• データサイエンス
• アンサンブルカルマンフィルタ