長岡工業高等専門学校研究紀要

論第 54 巻（2018）

文

アンサンブルカルマンフィルタFEMに基づく浅水流推定解析に対する時変サンプル数の適用
倉橋 貴彦1・小林 雅輝1・衛藤 俊彦2
1 長岡技術科学大学（Department of Mechanical Engineering, Nagaoka University of Technology）
2環境都市工学科（Department of Civil Engineering, National Institute of Technology, Nagaoka College）

Application of Time Variant Sample Number for Estimation Analysis
in Shallow Water Flow Based on the Ensemble Kalman Filter FEM
Takahiko KURAHASHI1, Masaki KOBAYASHI1 and Toshihiko ETO2
Abstract
In this paper, we present the investigation of influence of time variant sample number for estimation
results of shallow water flow using the ensemble Kalman filter FEM. The shallow water equation isemployed as the governing equation, and the SUPG method is applied to discretize the governing equation.
In this study, estimation results are compared by changing the sample number in time.
Key Words : ensemble Kalman filter FEM, SUPG method, shallow water equation, sample number

１．はじめに
２．アンサンブルカルマンフィルタFEM
非線形モデルを対象とした状態推定法として
Evensenによりアンサンブルカルマンフィルタ 1), 2)が開発され，近年では天気予報に適用し，高精度な状

本論文では，支配方程式として，式(1)，(2)に示す浅水長波方程式を用いる．

態推定に関する研究 が行われている．筆者らもア
3)

ンサンブルカルマンフィルタと有限要素法を融合した方法に関する研究を行っており，浅水流のモデル

𝑢̇ 𝑖 + 𝑢𝑗 𝑢𝑖,𝑗 + 𝑔𝜂,𝑖 − 𝜈(𝑢𝑖,𝑗 + 𝑢𝑗,𝑖 ),𝑗 = 0(1)𝜂̇ + {(ℎ + 𝜂)𝑢𝑖 },𝑖 = 0(2)に対して状態推定解析のプログラムを構築し，サンプル数の設定が水位の値に対して与える影響について考察を行っている4), 5)．従来の検討では，時間で一定のサンプル数を使用しており，各時間ステップ

ここに，ui，h, η，g，νはx,y方向の流速成分，基準において適切にサンプル数を変更できれば，計算時水深，基準水深からの水位変動量，重力加速度，動間の削減にも繋がると考えられる．そこで，本論文

粘性係数を示す．式(1)，(2)についてSUPG法およびでは，サンプル数を時間で変えた場合について状態後退オイラー法を適用し離散化する 6)．結果として推定解析を行い，推定精度に与える影響について検アンサンブルカルマンフィルタFEMにおけるシステ討を行う．解析モデルとしては，矩形の水路を対象

ム方程式は式(3)のように書くことができる．とする．

- 1 -倉橋 貴彦・小林 雅輝・衛藤 俊彦

𝑛+1𝑛𝑛[𝐴1 (𝜙)]{𝜙(𝑘)(−)} = [𝐴2 (𝜙)]{𝜙(𝑘)(+)} + {𝑞(𝑘)}(k=1,….M)

𝑛+1𝑛𝑛[𝐴1 (𝜙)]{𝜙(𝑘)(−)} = [𝐴2 (𝜙)]{𝜙(𝑘)(+)} + {𝑞(𝑘)}(3)

ここに[𝐴1 (𝜙)]，[𝐴2 (𝜙)]は有限要素法による定式化(k=1,….M)
8.

により得られる係数行列であり，{𝜙}は各節点の状状態量に対するアンサンブル行列の計算．
̃ 𝑛+1
̃ 𝑛+1
̃ 𝑛+1
̃ 𝑛+1 ] = [{ϕ[Φ(1)(+) } {ϕ(2)(+) } ⋯ {ϕ(𝑀)(+) }](+)

態変数ui，ηからなるベクトル，{𝑞}はシステムノイズを示す．また，nは時間ステップ，kはサンプル数

9.の番号，Mはサンプル数を示す．(+)および(－)の記号は同化前，同化後を示しており，式(3)の右辺の同化前の状態変数ベクトル{φn(k)(+)}は，カルマンゲ
イン行列[K]，観測値のベクトル{z}，観測行列[H]，観測ノイズ {r}を使用することにより与えられる

時間ステップnをn+1に更新し，状態量に対する𝑛+1
̃ (+)]→ア ン サ ン ブ ル 行 列 を 同 化 前 も の （ [Φ𝑛̃ (−) ] ）とし，ステップ2へ戻る．
[Φ３．数値実験（式(4)）．
図-1に解析モデルを示す．流れ場の推定解析にお𝑛𝑛{𝜙(𝑘)(+)} = {𝜙(𝑘)(−)}

ける観測点と推定点の位置を表-1に，またその他計算条件を表-2に示す．観測点は観測値を与える点で𝑛𝑛+[𝐾 𝑛 ] ({𝑧 𝑛 } − ([𝐻]{𝜙(𝑘)(−)} + {𝑟(𝑘)}))(k=1,….M)あり，推定点は，観測値をもとに推定した結果を表(4)示する点を示す．本検討では，図-2に示す水位変動に対する境界条件を図-1の左端の境界に与え，観測ここで，{𝜙̅ 𝑛 } および{𝑧̅𝑛 }を各サンプル数に対するベクトル{𝜙(𝑘) } および{𝑧(𝑘) }の平均値とし，各サンプル𝑛𝑛と 平 均 値 の 差 を { 𝜙̃(𝑘) } = { 𝜙(𝑘) } − {𝜙̅ 𝑛 } お よ び

値および真値をあらかじめ算出しておき，この解析

𝑛𝑛{ 𝑧̃(𝑘)} = { 𝑧(𝑘)} − {𝑧̅ 𝑛 }により表すと計算の流れは以

サンプル数の時間変化を考え，流入境界条件を時間で 微 分 し た cos 波 を 用 い て ， サ ン プ ル 数 は Case1

下の様に書くことができる

5), 7)．

により得られた観測値をもとに流れ場の推定解析を行う．流入境界条件（sin波）の変化量に関連した

M(t)=900|cos(2πt)|+100，Case2 M(t)=100|cos(2πt)|+100
1.

2.
3.

アンサンブル行列および初期タイムステップ
n=0の入力．
0̃ 0(2)(−) } ⋯ {ϕ̃ 0(𝑀)(−) }]̃ 0(1)(−) } {ϕ̃ (−)[Φ] = [{ϕと設定し（ 図-3 ）， 表-1 に示す推定点D，E，F点における水位変動の値を真値， M(t)=1000の一定値の観測値に対するアンサンブル行列の計算．
𝑛𝑛𝑛𝑛} ⋯ {z̃(𝑀)(−)}]} {z̃(2)(−)[Z̃(−)] = [{z̃(1)(−)

る．

共分散行列の計算．

を図-4～図-5 に示す．また，各推定点における真値

推定点 D～F における Case1 および Case2 の結果

1𝑛𝑛 𝑇[Z̃(−)][Z̃(−)]
𝑀−1との残差二乗和を 表-3～表-5 に示す．結果として，1𝑛̃(−)̃ 𝑛 ][Z̃ 𝑛 ]𝑇[V]=[Φ𝑀 − 1 (−) (−)
カルマンゲイン行列の計算．の一致率が高いことが確認できる．また，Case1 の𝑛̃ (−)[U]=4.

サンプル数を 1000 と固定した場合は，真値の値と結果もサンプル数が 1000 の結果と比較すると真値との一致率は若干下がるが，水位の変化の傾向は捉えることができていると考えられる．Case2 におい

𝑛𝑛 𝑇
̃ (−)
̃(−)
[𝐾 𝑛 ] = [U][V]

5.ては，水位の推定結果がどの点においても振動して最適推定値の計算．
𝑛𝑛{𝜙(𝑘)(+)} = {𝜙(𝑘)(−)}

おり，サンプル数を非定常的に変える場合においても，十分なサンプル数でない場合，適切な推定結果

𝑛𝑛+[𝐾 𝑛 ] ({𝑧 𝑛 } − ([𝐻]{𝜙(𝑘)(−)} + {𝑟(𝑘)}))(k=1,….M)
6.

最適推定値の平均値の計算．
𝑀𝑛{𝜙̅(+)}=

1𝑛∑{𝜙(+)(𝑖)}
𝑀
𝑖=1

7.

場合，Case1，Case2における水位の推定値と比較すシステム方程式の計算．

が得られないことがわかる．また，計算時間の結果の比較を表-6 に示す．この結果より，サンプル数を減らすことで計算時間を短縮できるということがわ
かる．サンプル数が 1000 の場合と，Case1，Case2における T=1 秒後の水面形を図-10～図-12 に示す．
Case1 においても，大きな数値振動なく水面形が得られていることが確認できるが Case2 では数値振動が大きいことを確認できる．

- 2 -アンサンブルカルマンフィルタ FEM に基づく浅水流推定解析に対する時変サンプル数の適用

図-1 有限要素分割図

図-3 サンプル数の経時変化

図-2 真値を作成の際に使用した境界条件
表-1 観測点と推定点の位置
ABCDEF

(観測点)
(観測点)
(観測点)
(推定点)
(推定点)
(推定点)

(x,y) = (0.0,0.5)(x,y) = (1.0,0.5)(x,y) = (2.0,0.5)(x,y) = (3.0,0.5)(x,y) = (4.0,0.5)(x,y) = (5.0,0.5)

図-4 推定点Dにおける水位の経時変化（Case1）

表-2 計算条件
時間増分量 Δt, s
時間ステップ数
基準水深 h, m
重力加速度 g, m/s2
初期条件 平均
初期条件 分散システムノイズ 平均
システムノイズ 分散
観測ノイズ 平均
観測ノイズ 分散0.1
1000
10.0
9.81
0.00.001
0.00.001
0.00.1図-5 推定点 E における水位の経時変化（Case1）- 3 -倉橋 貴彦・小林 雅輝・衛藤 俊彦

図-9 推定点Fにおける水位の経時変化（Case2）
図-6 推定点 F における水位の経時変化（Case1）

表-3 真値との残差二乗和の比較 (D 点)
M=1000
Case1 : M(t)=900|cos(2πt)|+100
Case2 : M(t)=100|cos(2πt)|+100

0.5882
1.607534.0462

表-4 真値との残差二乗和の比較 (E 点)
M=1000
Case1 : M(t)=900|cos(2πt)|+100
Case2 : M(t)=100|cos(2πt)|+100

1.13392.297739.2606

表-5 真値との残差二乗和の比較 (F 点)
M=1000
Case1 : M(t)=900|cos(2πt)|+100
Case2 : M(t)=100|cos(2πt)|+100

図-7 推定点Dにおける水位の経時変化（Case2）

1.29012.238821.9137

表-6 計算時間の比較
M=1000
Case1 : M(t)=900|cos(2πt)|+100
Case2 : M(t)=100|cos(2πt)|+100

1961.310s
1335.100s
314.213sm

図-8 推定点 E における水位の経時変化（Case2）m
図-10 1秒後の波形(M=1000)- 4 -アンサンブルカルマンフィルタ FEM に基づく浅水流推定解析に対する時変サンプル数の適用

系列において変えた場合においても，流れ場の推定解析を行えることを確認できた．本検討では流入条件を参考しサンプル数を経時的に変え検討を行ったが，サンプル数の変え方は本論文の方法以外にも設m定できるため，今後更なる検討を要すると考えられる．

謝辞：本論文の解析結果は，九州大学情報基盤研究開発センターの高性能アプリケーションサーバを使m用している．ここに謝意を表す．

図-11 1秒後の水面形(Case1)

参考文献
1)

G. Evensen, “Sequential data assimilation with a nonlinear
quasi-geostrophic model using Monte Carlo methods tom

forecast error statistics”, J. Geophys. Research, Vol.99,
pp.10143-10162, 1994.
2)

G. Evensen, “The Ensemble Kalman filter: theoreticalformulation and practical implementation”, Ocean Dyn.,
Vol.53, pp.343-367, 2003.

3)

三好建正，“アンサンブル・カルマンフィルタ－データ同化とアンサンブル予報の接点－”，天気，52
巻 2 号，pp.93-104，2005．

4)

T. Kurahashi, K. Saito and M. Nogami, “Application of theensemble Kalman filter FEM for estimation of flow field inshallow water regions”, JSIAM Letters, Vol.10, pp.4-8,2018.

5)

T. Kurahashi and K. Saito, “Investigation of the numericalaccuracy of shallow water flow estimation based on theensemble

Kalmanfilterusing

theSUPG FEM”,

International Journal of Computational Fluid Dynamics,
Vol.32, pp.146-157, 2018.
6)

倉橋貴彦，野上雅人，今井伸哉，吉荒太一，衛藤俊彦，波動水槽内の浅水流再現シミュレーション－
SUPG法による解析結果と水位計測値の比較－，長岡工業高等専門学校紀要，Vol.52, pp.21-27, 2016.

m7)

片 山 徹, “非 線 形 カ ル マ ン フ ィ ル タ”, 朝 倉 書 店,
pp.127-129, 2011.

図-12 1秒後の水面形(Case2)(2018. 8. 1 受付)

４．おわりに
本論文では，水路モデルに対してアンサンブルカ
ルマンフィルタFEMを適用し，サンプル数を時間で変えた場合について状態推定解析を行い，推定精度に関する検討を行った．結果としてサンプル数を時- 5 -