アンサンブルカルマンフィルタによる 浅水長波流れの推定解析

日本機械学会［No.187-1］北陸信越支部 第 55 期総会・講演会 講演論文集 [2018.3.3 福井県福井市]

［No.187-1］日本機械学会北陸信越支部 第 55 期総会・講演会 講演論文集 [2018.3.3 福井県福井市]

C012

アンサンブルカルマンフィルタによる 浅水長波流れの推定解析

齋藤 浄^*1，倉橋 貴彦^*2

Shallow water flow estimation analysis using the ensemble Kalman filter FEM

Kiyora SAITO

and Takahiko KURAHASHI

*1*2 Nagaoka University of Technology. Dept. of Mechanical Engineering Kamitomioka, Nagaoka-shi, Niigata, 1603-1 Japan

In this study, we present the estimation of shallow water flow field based on the ensemble Kalman filter FEM using the SUPG method. As the governing equation, the shallow water equation is introduced, and the Stream Upwind Petrov Gelerkin method (SUPG method), one of discretization method in the FEM, is applied to discretize the governing equation. The open channel model is emploed as the computational model.

Key Words : Shallow water flow, Ensemble Kalman filter FEM, SUPG method, Sample number.

1． 緒 言

アンサンブルカルマンフィルタは，Evensen⁽¹⁾⁽²⁾ により開発された初期条件の集合を用いたカルマンフィルタの 計算理論である．この理論では，リカッチ方程式を直接計算する必要がなく，カルマンゲインを近似的に解くこ とや，誤差共分散行列の計算において，システム方程式内の状態遷移行列を含まないため，非線形モデルの適用 も可能となる．また，一般にシステム方程式は非線形の方程式となることが多いこともあり，従来のカルマンフ ィルタに比べて，汎用的な計算理論であると考えられる．本研究では，開水路流れを対象とし，流れ場の推定シ ミュレーションを行う．数値実験では水路中央線上における真値と推定値の一致率に関する考察を行う．

2．アンサンブルカルマンフィルタFEMによる定式化および計算の流れ

2・1浅水長波方程式の離散化とアンサンブルカルマンフィルタ理論による計算アルゴリズム 本研究で対象とする非線形浅水長波方程式を式（1），式（2）に示す．

+ , + , − , + , , = 0 (1)

+ ℎ + , = 0 (2)

ここで は ， 方向流速m s ， は水位変動量m， は重力加速度m s ， ℎは基準水深m， はカルマ ン定数， は粘性係数（式（3））である．なお， ^∗はマニングの粗度係数である．

= 6 ∙ ^∗ ! " "

ℎ + ^#$ ℎ + (3)

*1 学生員，長岡技術科学大学 工学部

*2 正員，長岡技術科学大学

次に，重み付き残差法とθ法を用い空間方向に離散化する（式（4），（5））．

% ^∗+ & " ,"∗ + , + , − , + , , 'Ω

)* = 0 (4)

% ^∗+ & " ,"∗ + ℎ + , '+

,_* = 0 (5)

また，安定化パラメータ&を式（6）に示す．また式（6）における ^.は式（7）に示す．

τ = 012

Δ45 + 62 ^. ℎ. 7 8

9#

(6)

^. = : ; _.,

<

=#

(7)

式（6）においてΔ4は時間増分，ℎ.は要素長さである．ここで，要素領域を三角形とし，ある要素おける状 態変数 ， および重み関数 ^∗， ^∗に対して補間関数を適用すると式（8），（9）に示す有限要素方程式が得ら れる．

> ? @ +B CD E? @ + F D ? @ + G CH EA + CH EI J

+ & K" >L" ? @ + K"+ K CH"E? @ + KM H" ? @ = ?N @ (8)

> ? @ + ℎ +KKKKKKKKKK O ? @ +B O ? @ +P B O ?ℎ@ P

+& K"C>_L"E? @ + K"KKKKKKKKKK Hℎ + " ? @ + K"K H" ? @ = 0 (9)

ここでK は ， 方向要素内平均流速m s ，ℎKは要素内平均基準水深mである．また，時間方向に対してはθ 法を用いて離散化すると式（10）,（11）となる．

> + &K"C>_L"E ? ^<Q#@ + Δ4R K CO E + &K"K CH"E ? ^<Q#@

+ Δ4GR CH E ^<Q# Δ4 R O + &K" H" ? ^<Q#@ − > + &K"C>_L"E ? ^<@ + Δ4 1 − R K CO E + &K"K CH"E ? ^<@ + Δ4G 1 − R CH E ^<

+ Δ4 1 − R O + &K" H" ? ^<@ = ?0@

(10)

> + &K"C>_L"E + K O ? ^<Q#@ + Δ4 ℎ +KKKKKKKKKKR O + &K" H" ? ^<Q#@ + Δ4 O ?ℎ^<@

− > + &B BTH T ? ^<@ + Δ4ℎK 1 − R O + &K" H" ? ^<@ = ?0@ (11)

2・2アンサンブルカルマンフィルタ FEM による計算の流れ

アンサンブルカルマンフィルタFEMの計算手順を以下に示す． ここで， と>は時間ステップとサンプル番 号を示す． さらに（+）および（−）符号は，同化後の値および同化前の値を示す．推定値の初期値行列は，

状態変数の行ベクトルにより現れている．状態変数の行ベクトルは式（12）により計算する．

UVK^<₉ W = X"YUVZ^<_Q W[ + U\^< W (12) ここで推定値より推定値平均からの誤差を式（13），（14）により計算する．

IV^<₉J = 1

> :UVK^<9 W

]

=#

(13)

UV^^<₉ W = UVK^<₉ W − IV^<₉J T = 1, ⋯ , > (14)

結果として，誤差アンサンブル行列は，式（15）に示すように，状態推定誤差の行ベクトルによって表される．

C̀₉E = bUV^^{< #}₉ WUV^^<₉ W ⋯ UV^^{< ]}₉ Wc (15) また，観測値は式（16）のベクトルで示される．

Ud̀^<₉ W = H^< UVK^<₉ W + Uf^< W (16) ここで観測値より，観測値平均からの誤差を式（17），（18）と計算する．また式（17）より予測誤差アンサン ブル行列は式（19）のように与えられる．

Id^<₉J = 1

> :Ud̀^<9 W

]

=#

(17)

Ud̀^<₉ W = Ud̀ W − Id<9J (18)

Ch<9E = bUd̀^{< #}₉ WUd̀^<₉ W ⋯ Ud̀^{< ]}₉ Wc (19) 式（15）に示す推定値ベクトルの予測誤差アンサンブル行列と式（19）に示す観測ベクトルの予測誤差アンサ ンブル行列より，共分散行列はそれぞれ式（20），（21）のように計算される．

CiB<9E = 1

> − 1 C̀^<9EC̀<9E^k (20) ClK^<₉E = 1

> − 1 Ch^<9EC ^<₉E^k (21)

式（20），（21）を用いてカルマンゲイン行列を式（22）のように計算する．

n^< = bi<9

c bl<9

c^9# (22)

カルマンゲイン行列を用いて最適推定値を算出する（式（23），（24））．

UV^<₉ W = UV^<₉ W + n^< Y?d^<@ − Ud^<₉ W[ T = 1, ⋯ , > (23)

UV^<Q#₉ W = o^< UV^<₉ W + Ip^< J T = 1, ⋯ , > (24)

以降，この計算フローを繰り返すことにより推定値を算出する．

3． 矩形水路モデルに対する数値解析（マニング粗度係数の違いが解析結果へ与える影響）

粘性が解析結果に与える影響について検討するために図１の水路において図2に示す波形を左端の境界条 件から与え，式（10），（11）を用い順解析を行う．また，解析条件を表1に示す．

以上の条件においてマニングの粗度係数を変え計算を行う．各ケースにおけるマニングの粗度係数の値を 表2，点Dにおける水位変動量の推定結果を図3に示す．結果として，粘性項におけるマニングの粗度係数 ^∗ を変えた解析結果ではマニングの粗度係数の影響により水位変動量の最大値が減少し，水位変動量のピーク時 間が後方にシフトしていることを確認できる．

Fig. 1 Computation model

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 0.2 0.4 0.6 0.8 1

Water elevation, m

Time, s Sin wave

Fig. 2 Boundary condition for water elevation

Table 1 Computational conditions

Time increments 0.001

Number of time step 1000 Standard water depth ℎ， m 10 Gravitational acceleration ， m2/s 9.81

Parameter θ 1.0

Initial condition for and

Average 0

Variance 0.001

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 0.2 0.4 0.6 0.8 1

Water elevation, m

Time, s case1

case2 case3 case4 case5

Fig. 3 Time history of water elevation at point D

Table 2 Set of Manning roughness coefficient

case1 ∗= 0.0 s/m^#$

case2 ∗= 0.01 s/m^#$

case3 ∗= 0.025 s/m^#$

case4 ∗= 0.05 s/m^#$

case5 ∗= 0.1 s/m^#$

4． アンサンブルカルマンフィルタを用いた解析結果の検討

4・1 計算条件と観測値の作成

アンサンブルカルマンフィルタFEMによる流れ場の推定計算を行うにあたり3章で求めた順解析の結果を用 い，人工的に観測値を作成する．表1に示す平均・分散の値のもと，各節点に異なる正規乱数を計算し，流速・

水位変動量にそれぞれ付加する．観測点は図1に示すA，B，Cの三点であり，観測点Cでの観測値を横軸時 間，縦軸水位変動量とした場合のグラフを図4に示す．

Fig.4 Observed value at point C 4・2 波形推定

4・1節で作成した観測値を使用し，水路全体において波形が伝播する様子を推定する．サンプル数と推定 精度の関係を確認するためにサンプル数を表3に示すよう変え，確認点D，E，Dにおいて水位変動量と時間 変化の関係を図5～7に示す．

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 0.2 0.4 0.6 0.8 1

Water elevation, [m]

Time, [s]

observed value

Table 3 Calculation conditions case1 sample number > = 100 case2 sample number > = 175 case3 sample number > = 250

case4 sample number > = 1000 ^-2

-1.5 -1 -0.5 0 0.5 1 1.5 2

0 0.2 0.4 0.6 0.8 1

Water elevation, m

Time, s case1

case2 case3 case4 True value

Fig. 5 Estimated value at point D for each ensemble number

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 0.2 0.4 0.6 0.8 1

Water elevation, [m]

Time, [s]

case1 case2 case3 case4 True value

Fig. 6 Estimated value at point E for each ensemble number

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 0.2 0.4 0.6 0.8 1

Water elevation, [m]

Time, [s]

case1 case2 case3 case4 True value

Fig. 7 Estimated value at point F for each ensemble number

4・3 推定結果

図5～図7より各確認点において，サンプル数を1000とした場合の推定値と真値は良く一致していること がわかる．サンプル数が減少するに従い，細かなノイズと周期的な振動が発生し，サンプル数100の場合の推 定結果は，最終的に発散してしまっている．ここで，各確認点において推定値の真値の一致率の違いを検証す るために，各時刻ステップにおける推定値と真値の差に関する二乗値を計算し，0～1.4秒間において和を計算 したものを表4に整理する．結果として，流入点に近いほど，推定値は真値に近くなることを確認した．

Table 4 Square sum of residual between true value and estimated value in case of sample number M=1000

Point D 1.7164

Point E 1.7278

Point F 2.6578

5． 結 語

非線形浅水長波方程式を適用し水路内の流れ場推定を行い，マニングの粗度係数の値が水位変動量の推定結果 へ与える影響を確認した．また，本検討モデルではサンプル数は1000程度用意することにより，流れ場の推定を 適切に行うことが出来ることを確認した．

謝 辞

本論文を執筆するにあたり科学研究費助成金（基盤（C））15K5786 の援助を受けた．ここに謝意を表す．ま た，本論文中に示した数値解析を行うにあたり，九州大学情報基盤研究開発センターの高性能演算サーバーシス テムを利用させていただいた．九州大学情報基盤研究開発センターのスタッフに方々に対して謝意を表す．

文 献

(1) G. Evensen，“Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics”, Journal of Geophysical Research,Vol. 99, (1986), pp. 10143-10162.

(2) G. Evensen，“The Ensemble Kalman filter: theoretical formulation and practical implementation” Ocean Dynamics, Vol. 53, (2003), pp. 343-367.