Lorenz96モデルのデータ同化:Local Ensemble Transform Kalman Filter

はじめに

ずいぶん間があいてしまいましたが前回のEnKFの基礎の記事の続きです。今回はETKFとその局所化を導入した形であるLETKFの解説をします。EAKFの解説は飛ばすことにしました。申し訳ありません。
LETKFは並列計算による効率化が図れることとサンプル数が少なくても動作するという性質から大規模計算に強く、現実の天気予報でも用いられています。続きものの記事なので、データ同化の基礎やExtended Kalman Filterの記事に書いてあることは全て前提とします。わからない用語はそちらを参照してください。

EnKF:Ensemble Transform Kalman Filter

直接$P^{a},P^{f}$で計算するとサイズ$N\times N$の行列計算をする必要があります。
サンプルから作った行列$\delta X^{f} = (x^{f,1}-\bar{x}^{f},\ldots,x^{f,m}- \bar{x}^{f})$で計算すれば
サイズ$N \times m$の行列計算でよく、一般にサンプル数はモデルの次元より小さいのでこちらの方が計算が効率的です。
そこで$P^{a} = (I-KH)P^{f}$の代わりに$\delta X^{a} = \delta X^{f} T$という形にして、
平均からのずれ$\delta X^{f}$に変換$T$をかけて計算すればいいという考え方に基づいたのがETKFです．

導出

けっこう面倒なので読み飛ばしても大丈夫だと思います。
$\delta X^{a} = (x^{a,1}-\bar{x}^{a},\ldots,x^{a,m}-\bar{x}^{a}) $に対して
サイズ$m\times m$の変換行列$T$を用いて$\delta X^{a} =\delta X^{f} T$という形を仮定し、
$P^a = \frac{1}{m-1} \delta X^{a} (\delta X^{a})^{\top} = \frac{1}{m-1} \delta X^{f} TT^{\top}(\delta X^{f})^{\top}$が$(I-KH)P^{f}$に等しくなるように$T$を定めます。
$K = P^{f}H^{\top}(HP^{f}H^{\top}+R)^{-1}$をサンプル行列$\delta X^{f}$を用いて表して変形することで

\begin{align}
(I-KH)P^{f} &= \left(I -\frac{1}{m-1} \delta X^{f} (\delta X^{f})^{\top}H^{\top} \left(H \frac{1}{m-1} \delta X^{f} (\delta X^{f})^{\top} H^{\top} + R \right)^{-1} H \right) \frac{1}{m-1} \delta X^{f} (\delta X^{f})^{\top}\\
&= \frac{1}{m-1} \left(I -\delta X^{f} (\delta X^{f})^{\top}H^{\top} \left(H \delta X^{f} (\delta X^{f})^{\top} H^{\top} + (m-1)R \right)^{-1} H \right)\delta X^{f} (\delta X^{f})^{\top}\\
&= \frac{1}{m-1} \left(I -\delta X^{f} (\delta Y^{f})^{\top} \left(\delta Y^{f} (\delta Y^{f})^{\top} + (m-1)R \right)^{-1} H \right)\delta X^{f} (\delta X^{f})^{\top}\\
&= \frac{1}{m-1} \delta X^{f} \left(I -(\delta Y^{f})^{\top} \left(\delta Y^{f} (\delta Y^{f})^{\top} + (m-1)R \right)^{-1} H \delta X^{f} \right) (\delta X^{f})^{\top}
\end{align}

となります。したがって

\begin{align}
TT^{\top} &= I -(\delta Y^{f})^{\top} \left(\delta Y^{f} (\delta Y^{f})^{\top} + (m-1)R \right)^{-1} \delta Y^{f} \\
&= I - \left((m-1)I +  (\delta Y^{f})^{\top}R^{-1} \delta Y^{f} \right)^{-1} (\delta Y^{f})^{\top} R^{-1}\delta Y^{f}\\
&= \left((m-1)I +  (\delta Y^{f})^{\top}R^{-1} \delta Y^{f} \right)^{-1} \left( (m-1)I +  (\delta Y^{f})^{\top}R^{-1} \delta Y^{f} - (\delta Y^{f})^{\top}R^{-1} \delta Y^{f}\right)\\
&= (m-1) \left((m-1)I +  (\delta Y^{f})^{\top}R^{-1} \delta Y^{f} \right)^{-1}
\end{align}

を満たせば十分です。途中で$A^{\top}(AA^{\top}+B)^{-1} = (I+A^{\top}B^{-1}A)^{-1}A^{\top}B^{-1}$を用いました。
このような$T$にはユニタリ行列分の自由度が残りますが、

T = \sqrt{(m-1) \left((m-1)I +  (\delta Y^{f})^{\top}R^{-1} \delta Y^{f} \right)^{-1}}

とするのが計算の安定性、アンサンブルが持つべき性質の観点などから推奨されています。
解析値のアンサンブルの更新は、$(M)_{k}$で行列$M$の$k$列目を表すことにして、

\begin{align}
x^{a,k} &= \bar{x}^{a} + (\delta X^{a})_{k}\\
&= \bar{x}^{f} + K(y^{o} - H \bar{x}^{f}) + (\delta X^{f} T)_{k} \\
&= \bar{x}^{f} + \delta X^{f} \left((m-1)I +  (\delta Y^{f})^{\top}R^{-1} \delta Y^{f} \right)^{-1}
(\delta Y^{f})^{\top} R^{-1} (y^{o} - H\bar{x}^{f})+ (\delta X^{f} T)_{k}\\
&= \bar{x}^{f} + \delta X^{f} \left( \tilde{P}^{a} C (y^{o}-\bar{y}^{f}) + (T)_{k} \right)
\end{align}

によって求まります。ここで最後の行への変形で

\begin{align}
\tilde{P}^{a} &= \left((m-1)I +  (\delta Y^{f})^{\top}R^{-1} \delta Y^{f} \right)^{-1}\\
C &= (\delta Y^{f})^{\top} R^{-1}\\
\bar{y}^{f} &= H \bar{x}^{f}
\end{align}

としました。

アルゴリズム

工夫の余地は沢山あるのですが、例えば"Brian R.Hunt, Eric J.Kostelich,Istvan Szunyogh,Efficient data assmilation for spatiotemporal chaos: A local ensemble transform Kalman filter,Physica D, 2007"に従ってアルゴリズムを示します。
Forecast Step

\begin{align}
x^{f,k} &= M(x^{a,k})\\
\delta X^f &= \left(x^{f,1}-\bar{x}^f, \ldots , x^{f,m} - \bar{x}^f \right) \\
\bar{y}^{f} &= H \bar{x}^{f} \\
\delta Y^f &= H \delta X^{f}
\end{align}

Analysis Step

\begin{align}
C &= (\delta Y^{f})^{\top} R^{-1} \\
\tilde{P}^{a} &= \left((m-1)I + C \delta Y^{f} \right)^{-1}\\
T &= \sqrt{(m-1)\tilde{P}^{a}} \\
\bar{w}^{a} &= \tilde{P}^{a} C (y^{o} - \bar{y}^{f})\\
w^{a,k} &= \bar{w}^{a} + (T)_{k} \\
x^{a,k} &= \bar{x}^{f} + \delta X^{f}w^{a,k}

\end{align}

実装と解説

例のごとく全体を示してから解説することにします。

using LinearAlgebra
using Statistics
using Random
using DelimitedFiles


##L96modelの右辺
function L96(u; F = 8.0, N = 40)
    f = fill(0.0, N)
    for k = 3:N-1
        f[k] = (u[k+1] - u[k-2]) * u[k-1] - u[k] + F
    end
    f[1] = (u[2] - u[N-1]) * u[N] - u[1] + F
    f[2] = (u[3] - u[N]) * u[1] - u[2] + F
    f[N] = (u[1] - u[N-2]) * u[N-1] - u[N] + F

    return f
end

#4-4Runge-Kutta
function Model(u; dt = 0.05)
    du = u
    s1 = L96(u)
    s2 = L96(u + s1 * dt / 2)
    s3 = L96(u + s2 * dt / 2)
    s4 = L96(u + s3 * dt)
    du += (s1 + 2 * s2 + 2 * s3 + s4) * (dt / 6)
    return du
end

function main()
    Time_Step = 14600
    member = 30
    rho = 1.05
    N = 40
    M = 40
    F = 8.0
    IN = Matrix(1.0I, N, N)

    u_true = readdlm("L96_truestate.txt")
    u_obs = readdlm("L96_observation.txt")


    H = Matrix(1.0I, M, N)
    R = Matrix(1.0I, M, M)

    ua = fill(0.0, (N,member))
    uf_mean = fill(0.0, N)
    ua_mean = fill(0.0, N)
    uf = fill(0.0, (N,member))
    dXf = fill(0.0, (N,member))
    dYf = fill(0.0, (M, member))
    dXa = fill(0.0, (N,member))

    open("L96_ETKF_output.txt","w") do output
        for m = 1:member
            ua[:, m] = rand(N) .+ F
            for i = 1:Time_Step
                ua[:,m] = Model(view(ua,:,m))
            end
        end

        for i in 1:Time_Step

            ##forecast step
            for m in 1:member
                uf[:,m] = Model(view(ua,:,m))
            end
            uf_mean = mean(uf, dims = 2)
            yf_mean = H*uf_mean
            dXf = uf .- uf_mean
            dYf = H * dXf

            ##analysis step
            C = transpose(dYf)*inv(R)
            Pa = inv((member-1)/rho * I + C*dYf)
            dWa = sqrt(Symmetric((member-1)*Pa))
            wa = Pa * C * (u_obs[i,2:N+1] - yf_mean)
            W = dWa .+ wa
            ua_mean = uf_mean + dXf * wa
            ua = dXf*W .+ uf_mean

            writedlm(output,[(i / 40) (norm(u_true[i, 2:N+1] - ua_mean)/sqrt(N)) (norm(u_obs[i, 2:N+1] -u_true[i, 2:N+1])/sqrt(N))])
        end
    end
end

main()

Forecast Step

\begin{align}
x^{f,k} &= M(x^{a,k})\\
\delta X^f &= \left(x^{f,1}-\bar{x}^f, \ldots , x^{f,m} - \bar{x}^f \right) \\
\bar{y}^{f} &= H \bar{x}^{f} \\
\delta Y^f &= H \delta X^{f}
\end{align}

            ##forecast step
            for m in 1:member
                uf[:,m] = Model(view(ua,:,m))
            end
            uf_mean = mean(uf, dims = 2)
            yf_mean = H*uf_mean
            dXf = uf .- uf_mean
            dYf = H * dXf

ほぼそのまま書けます。$H$が非線形の場合は$\bar{y}^{f} = H(\bar{x}^{f})$および$\delta Y^f = \left(H(x^{f,1})-\bar{y}^f, \ldots , H(x^{f,m}) - \bar{y}^f \right)$としてあとは同様にすることが提案されています。

Analysis Step

\begin{align}
C &= (\delta Y^{f})^{\top} R^{-1} \\
\tilde{P}^{a} &= \left((m-1)I + C \delta Y^{f} \right)^{-1}\\
T &= \sqrt{(m-1)\tilde{P}^{a}} \\
\bar{w}^{a} &= \tilde{P}^{a} C (y^{o} - \bar{y}^{f})\\
w^{a,k} &= \bar{w}^{a} + (T)_{k} \\
x^{a,k} &= \bar{x}^{f} + \delta X^{f}w^{a,k}

\end{align}

            ##analysis step
            C = transpose(dYf)*inv(R)
            Pa = inv((member-1)/rho * I + C*dYf)
            dWa = sqrt(Symmetric((member-1)*Pa))
            wa = Pa * C * (u_obs[i,2:N+1] - yf_mean)
            W = dWa .+ wa
            ua_mean = uf_mean + dXf * wa
            ua = dXf*W .+ uf_mean

こちらもほぼそのままです。Pa = inv((member-1)/rho * I + C*dYf)におけるrhoはInflation parameterで、$\delta X^{f}$を$\sqrt{\rho}\delta X^{f}$に置き換えることと等価です。

結果

サンプル数を30とし、rho=1.05とした結果を下に示します。

localizationなしでもうまくいっていることがわかります。サンプル数は少なければ少ないほど計算が楽なのでもっと減らしたいのですが、このままサンプル数を減らしていくと10個あたりで動作しなくなります。そこで局所化をうまく導入したいのですが、サンプルの次元で計算を済ませるには少し工夫が必要で、それが次節のLETKFです。

EnKF:Local ETKF

サンプル数をさらに削減したいので、先ほどのETKFに対して局所化を導入します。共分散行列の対角成分を小さくする（すなわち、遠くにある予報データの誤差の相関を強制的に小さくする）ように成分ごとの積をとる、という方法を取ろうとすると$P^{f}$や$P^{a}$の各成分ごとの計算が必要になり、$\delta X^{a}$などで閉じた計算ができません。
そこで、ある分割領域の予測データに対してその近傍にある観測のみを取り込む、という形で局所化を行います。このような局所化のアルゴリズム(LEKF)自体は"Ott et al,Exploiting Local Low Dimensionality of the Atmospheric Dynamics for Efficient Ensemble Kalman Filtering,2002"の時点ですでに存在していますが、明示的にETKFと組み合わせてLETKFという名称を提唱したのはETKFのところでも引用したHunt et al,2007だと思います。

###導出
Lorenz96モデルでは変数が40個あるので、これをいくつかに分割します。分割の仕方は自由ですが、例えば$(x_{1},\ldots,x_{40})$を$(x_{1},\ldots,x_{5}),(x_{6},\ldots,x_{10}),\ldots,(x_{36},\ldots,x_{40})$のように5成分ずつに分割したとしましょう。各分割に対して、どの観測を同化するのかを決めます。たとえば分割$(x_{1},\ldots,x_{5})$に対しては観測$(y_{1},\ldots,y_{40})$のうち$(y_{1},\ldots,y_{5})$から遠くにある成分はあまり関係がないでしょう。そこで$(y_{1},\ldots,y_{5})$を同化することにします。もちろん近い成分なら関係があると思われるので$(y_{40},y_{1},y_{2},\ldots,y_{5},y_{6})$を同化しても良いかもしれません。

分割$(x_{1},\ldots,x_{5})$に対応するサンプル$x_{l}^{f,k} =(x_{1}^{f,k},\ldots,x_{5}^{f,k})$を取ってきて、$\delta X_{l}^{f}$を作ります。同様に同化することに決めた観測$(y_{1},\ldots,y_{5})$に対応する観測$y_{l}^{o} =(y_{1}^{o},\ldots,y_{5}^{o})$と共分散$R_{l}$を取ってきます。あとは通常のETKFと同様の計算を行って、解析アンサンブル$x_{l}^{a}$を求めます。

この計算できちんと同化ができるということを一応示しておきます。面倒なので読み飛ばしても問題ないと思います。
$P_{l} \colon \mathbb{R}^{N} \to \mathbb{R}^{d_{l}}$を全変数に対して分割の各成分を対応させる行列とし、
$Q_{l} \colon \mathbb{R}^{d_{o}} \to \mathbb{R}^{d_{l}^{o}}$を分割に対して同化する観測を対応させる行列とします。
局所化は$P^{f}$を$\sum_{l}P_{l} P^{f} P_{l}^{\top}$とし、$R$を$\sum_{l} Q_{l} R Q_{l}^{\top}$とすることによって行います。
例のような5成分ずつの分割では$5\times 5$の行列を斜めに並べたブロック対角化で局所化していることに相当します。

\begin{align}
\mathrm{tr}P^{a} &= \mathrm{tr} (I-KH)P^{f}(I-KH)^{\top} + \mathrm{tr}KRK^{\top} \\
&= \sum_{l} \left(\mathrm{tr} (I-KH)P_{l} P^{f} P_{l}^{\top} (I-KH)^{\top} + \mathrm{tr}KQ_{l}R Q_{l}^{\top}K^{\top}\right)\\
&=\sum_{l} \left(\mathrm{tr} (P_{l}-KHP_{l})P^{f} (P_{l}-KHP_{l})^{\top}
 + \mathrm{tr}KQ_{l}R Q_{l}^{\top} K^{\top}\right)
\end{align}

の最小化をすればよく、

\begin{equation}
K = \sum_{l} (P_{l} P^{f}P_{l}^{\top} H^{\top})(HP_{l}P^{f}P_{l}^{\top}H^{\top} + Q_{l} R Q_{l}^{\top})^{-1}
\end{equation}

とすればよいことになります。これは各分割上で通常のETKFのアルゴリズムを実行することと等価です。

アルゴリズム

Forecast Stepは特に変更はありません。Analysis Stepにおいて局所化を行います。

Analysis Step
各$l$に対して$\delta Y_{l}^{f},y_{l}^{o},R_{l}$および$\bar{y}_{l}^{f}$を全体から各分割の成分を抜いてきて作ります。あとは通常のETKFと同様にして

\begin{align}
C_{l} &= (\delta Y_{l}^{f})^{\top} R_{l}^{-1} \\
\tilde{P_{l}}^{a} &= \left((m-1)I + C_{l} \delta Y_{l}^{f} \right)^{-1}\\
T_{l} &= \sqrt{(m-1)\tilde{P}_{l}^{a}} \\
\bar{w}_{l}^{a} &= \tilde{P}_{l}^{a} C_{l} (y_{l}^{o} - \bar{y}_{l}^{f})\\
w_{l}^{a,k} &= \bar{w}_{l}^{a} + (T_{l})_{k} \\
x_{l}^{a,k} &= \bar{x}_{l}^{f} + \delta X_{l}^{f}w_{l}^{a,k}

\end{align}

として解析値を求めます。全ての$l$に対して更新が終われば解析値$x^{a}$が求まっています。どの順番で計算をしても構わないことが重要で、実際の大規模計算ではこの部分を並列計算で行うことで計算の効率化を図っているとのことです。

実装と解説

全体を示してからAnalysis Stepのところだけ解説をします。

using LinearAlgebra
using Statistics
using Random
using DelimitedFiles


##L96modelの右辺
function L96(u; F = 8.0, N = 40)
    f = fill(0.0, N)
    for k = 3:N-1
        f[k] = (u[k+1] - u[k-2]) * u[k-1] - u[k] + F
    end
    f[1] = (u[2] - u[N-1]) * u[N] - u[1] + F
    f[2] = (u[3] - u[N]) * u[1] - u[2] + F
    f[N] = (u[1] - u[N-2]) * u[N-1] - u[N] + F

    return f
end

#4-4Runge-Kutta
function Model(u; dt = 0.05)
    du = u
    s1 = L96(u)
    s2 = L96(u + s1 * dt / 2)
    s3 = L96(u + s2 * dt / 2)
    s4 = L96(u + s3 * dt)
    du += (s1 + 2 * s2 + 2 * s3 + s4) * (dt / 6)
    return du
end

function main()
    Time_Step = 14600
    member = 5
    local_dim = 5
    local_obs_dim = 40
    rho = 1.3
    N = 40
    M = 40
    F = 8.0
    IN = Matrix(1.0I, N, N)

    u_true = readdlm("L96_truestate.txt")
    u_obs = readdlm("L96_observation.txt")


    H = Matrix(1.0I, M, N)
    R = Matrix(1.0I, M, M)

    ua = fill(0.0, (N,member))
    uf_mean = fill(0.0, N)
    ua_mean = fill(0.0, N)
    uf = fill(0.0, (N,member))
    dXf = fill(0.0, (N,member))
    dYf = fill(0.0, (M, member))
    dXa = fill(0.0, (N,member))

    open("L96_LETKF_output.txt","w") do output
        for m = 1:member
            ua[:, m] = rand(N) .+ F
            for i = 1:Time_Step
                ua[:,m] = Model(view(ua,:,m))
            end
        end

        for i in 1:Time_Step

            ##forecast step
            for m in 1:member
                uf[:,m] = Model(view(ua,:,m))
            end
            uf_mean = mean(uf, dims = 2)
            yf_mean = H*uf_mean
            dXf = uf .- uf_mean
            dYf = H * dXf

            ##analysis step

            ##devide into local patch
            patch_num::Int64 = N / local_dim

            for local_ind in 1:patch_num
                uf_local_begin::Int64 = (local_ind -1)*local_dim + 1
                uf_local_end::Int64 = local_ind*local_dim
                obs_local_begin = uf_local_begin
                obs_local_end = uf_local_end
                local_obs_dim = obs_local_end-obs_local_begin+1

                uf_mean_local = uf_mean[uf_local_begin:uf_local_end]
                dXf_local = dXf[uf_local_begin:uf_local_end,:]

                dYf_local = dYf[obs_local_begin:obs_local_end,:]
                yf_mean_local = yf_mean[obs_local_begin:obs_local_end]
                u_obs_local = u_obs[i,obs_local_begin+1:obs_local_end+1]
                R_local = Matrix(1.0I,local_obs_dim,local_obs_dim)


                C = transpose(dYf_local)*inv(R_local)
                Pa_local = inv((member-1)/rho * I + C*dYf_local)
                dWa_local = sqrt(Symmetric((member-1)*Pa_local))
                wa_mean = Pa_local*C*(u_obs_local - yf_mean_local)

                wa = dWa_local .+ wa_mean
                ua_mean[uf_local_begin:uf_local_end] = uf_mean_local + dXf_local*wa_mean
                ua[uf_local_begin:uf_local_end,:] = dXf_local*wa .+ uf_mean_local
            end

            writedlm(output,[(i / 40) (norm(u_true[i, 2:N+1] - ua_mean)/sqrt(N)) (norm(u_obs[i, 2:N+1] -u_true[i, 2:N+1])/sqrt(N))])
        end
    end
end

main()

Analysis Step

Forecast Stepについては特に変更がないので、Analysis Stepについてのみ解説をします。考えるのが面倒になって最初に思いついた書き方でそのまま実装したので、観測数を減らすなどするとこのままでは動作しなくなります。各自で工夫してより良い書き方をして教えてくださると助かります。抜き出したものが次になります。

            ##analysis step

            ##devide into local patch
            patch_num::Int64 = N / local_dim

            for local_ind in 1:patch_num
                uf_local_begin::Int64 = (local_ind -1)*local_dim + 1
                uf_local_end::Int64 = local_ind*local_dim
                obs_local_begin = uf_local_begin
                obs_local_end = uf_local_end
                local_obs_dim = obs_local_end-obs_local_begin+1

                uf_mean_local = uf_mean[uf_local_begin:uf_local_end]
                dXf_local = dXf[uf_local_begin:uf_local_end,:]

                dYf_local = dYf[obs_local_begin:obs_local_end,:]
                yf_mean_local = yf_mean[obs_local_begin:obs_local_end]
                u_obs_local = u_obs[i,obs_local_begin+1:obs_local_end+1]
                R_local = Matrix(1.0I,local_obs_dim,local_obs_dim)


                C = transpose(dYf_local)*inv(R_local)
                Pa_local = inv((member-1)/rho * I + C*dYf_local)
                dWa_local = sqrt(Symmetric((member-1)*Pa_local))
                wa_mean = Pa_local*C*(u_obs_local - yf_mean_local)

                wa = dWa_local .+ wa_mean
                ua_mean[uf_local_begin:uf_local_end] = uf_mean_local + dXf_local*wa_mean
                ua[uf_local_begin:uf_local_end,:] = dXf_local*wa .+ uf_mean_local
            end

前半部分で全体のデータから局所化された部分へと射影をします。
今回は40変数をlocal_dim成分ずつに等分割しました。

                uf_local_begin::Int64 = (local_ind -1)*local_dim + 1
                uf_local_end::Int64 = local_ind*local_dim
                obs_local_begin = uf_local_begin
                obs_local_end = uf_local_end
                local_obs_dim = obs_local_end-obs_local_begin+1

                uf_mean_local = uf_mean[uf_local_begin:uf_local_end]
                dXf_local = dXf[uf_local_begin:uf_local_end,:]

                dYf_local = dYf[obs_local_begin:obs_local_end,:]
                yf_mean_local = yf_mean[obs_local_begin:obs_local_end]
                u_obs_local = u_obs[i,obs_local_begin+1:obs_local_end+1]
                R_local = Matrix(1.0I,local_obs_dim,local_obs_dim)

uf_local_beginとuf_local_endが予報値の分割がどの成分から始まって終わるかを表します。obs_local_beginとobs_local_endが観測の分割がどの成分から始まって終わるかを表します。今回は簡単に予報と観測で同じ成分を取ることにしました。ここはもっと凝ってもいいと思います。あとは対応する成分を全体のuf_meanやdXfから抜いてきて局所的な観測や予報値を作っています。

後半部分で局所的に解析値を求めます。

\begin{align}
C_{l} &= (\delta Y_{l}^{f})^{\top} R_{l}^{-1} \\
\tilde{P_{l}}^{a} &= \left((m-1)I + C_{l} \delta Y_{l}^{f} \right)^{-1}\\
T_{l} &= \sqrt{(m-1)\tilde{P}_{l}^{a}} \\
\bar{w}_{l}^{a} &= \tilde{P}_{l}^{a} C_{l} (y_{l}^{o} - \bar{y}_{l}^{f})\\
w_{l}^{a,k} &= \bar{w}_{l}^{a} + (T_{l})_{k} \\
x_{l}^{a,k} &= \bar{x}_{l}^{f} + \delta X_{l}^{f}w_{l}^{a,k}

\end{align}

                C = transpose(dYf_local)*inv(R_local)
                Pa_local = inv((member-1)/rho * I + C*dYf_local)
                dWa_local = sqrt(Symmetric((member-1)*Pa_local))
                wa_mean = Pa_local*C*(u_obs_local - yf_mean_local)

                wa = dWa_local .+ wa_mean
                ua_mean[uf_local_begin:uf_local_end] = uf_mean_local + dXf_local*wa_mean
                ua[uf_local_begin:uf_local_end,:] = dXf_local*wa .+ uf_mean_local
 

ほぼそのままです。ua_localを用いる代わりに直接 ua[uf_local_begin:uf_local_end,:]として全体の解析値の一部分を更新するようにしています。一つ一つの行列のサイズは高々局所的な次元×サンプル数程度なので計算自体は簡単です。大規模計算ではこの部分を並列化して計算することで効率化を図っているとのことです。

結果

各分割の要素数を5とし、サンプル数を5として、rho=1.3とした結果を示します。

サンプル数が5でも（精度は劣りますが）動作することがわかります。サンプル数8ではRMSEは0.2前後まで下がります。サンプル数が少なくて済む点、並列計算によって計算効率の向上を測れる点でLETKFは優れた手法で、現在でも用いられています。一方でどのように局所的に分割するか、各分割に対してどの観測を取り込んで同化するのかといったところは一般的な指針のようなものはなく¹、今のところモデルごとに手でチューニングするしか無いようです。

まとめ

今回はETKFとその局所化を取り入れた形であるLETKFの解説をしました。LETKFではサンプル数5でRMSEは0.4前後でしたが、分割の仕方や取り込む観測の工夫次第ではもっとRMSEを落とせるかもしれません。この記事を読んで挑戦してくださるとうれしいです。次回は3次元変分法と4次元変分法の解説になると思いますが、その間に接線形コード、アジョイントコードやそれを用いた最適化について一つ記事を挟むことになると思います。

1. もちろん大きな需要はあるようで、前述のHunt et al,2007の参考文献やBred Vector,Singular Vector,Assimilation in Unstable Subspaceといったキーワードで調べてみると沢山文献が出てきます。よい参考文献があれば教えていただけると助かります。 ↩