Influence of boundary condition for result of flow field estimation in shallow water region based on Kalman filter FEM

010

Copyright @2016社団法人 日本機械学会

カルマンフィルタ

FEM

Influence of boundary condition for result of flow field estimation in shallow water region based on Kalman filter FEM

吉荒 太一・長岡技術科学大学大学院

Taichi Yoshiara, Graduate School of Nagaoka University of Technology 正 倉橋 貴彦・長岡技術科学大学大学院

Takahiko Kurahashi, Nagaoka University of Technology

正 小林 泰秀・長岡技術科学大学大学院

Yasuhide Kobayashi, Nagaoka University of Technology Key Words: Kalman filter, Finite element method, Boundary condition, Flow field estimation, Shallow water equation

Abstract: In this study, we present numerical examples for estimation of shallow water flow based on the Kalman filter FEM. The shallow water equation is employed as the governing equation, and the Galerkin method using the triangular element and the selective lumping method are applied to discretize the governing equation in space and time, respectively. We investigated difference of numerical results by changing position of the observation points. It was found that high accurate flow field estimation result is obtained in case that the observation points are set to the upstream side in the open channel in comparison with the case that the observation points are set to the downstream side in the open channel. Therefore, inflow boundary condition is introduced to estimation analysis of shallow water flow based on Kalman filter FEM, and examinations of inflow boundary condition on the estimation accuracy are carried out.

1. 序論

カルマンフィルタ⁽¹⁾⁽²⁾は誤差を含む観測値から対象となる システムの状態量を高精度に推定することができることか ら，様々な工学的分野で用いられている．空間的に広がりを 持ち，時間進展するモデルに対してカルマンフィルタを適用 する場合，空間における状態量分布の離散化式を導くために，

境界要素法や有限要素法が適用されてきた⁽³⁾⁽⁴⁾．本研究では，

カルマンフィルタ理論および有限要素法に基づき，ある水域 における流れ場の推定解析を行うことに焦点を当てる．現象 を表すモデルは平面2次元場の線形浅水長波方程式とし，空 間方向には有限要素法（ガラーキン法），時間方向に対して はセレクティブ・ランピング法により離散化を行った．これ までにもカルマンフィルタを用いて，水域内の流れ場を推定 する検討⁽⁵⁾が行われてきた．先の研究⁽⁶⁾では使用する観測変 数（流速および水位変動量），観測点の数および位置の条件 が推定精度に与える影響について開流路を対象に数値実験 が行われ，その関係性が示された．この検討から，カルマン フィルタFEMに基づく流れ場の推定を行う際，流入境界側 に観測点を配置しない場合には，流入境界から最初の観測点 までの区間における流れ場の推定は困難であることがわか っている．そこで本研究では，流入境界側に観測点を配置し ない場合であっても，適切な境界条件を与えることで推定が 可能となるのか検証し，境界条件が推定精度に与える影響に ついて開流路を対象に数値実験を行う．

2. 支配方程式の離散化

対象とする流れは式(1)および式(2)に示す平面2次元場の 線形浅水長波方程式で表す．

𝑢̇𝑖+ 𝑔𝜂,𝑖= 0 (1)

𝜂̇ + ℎ𝑢_𝑖,𝑖= 0 (2)

ここに，𝑢_𝑖は𝑥，𝑦方向に対する流速，𝜂は平均水位からの水

位変動量，𝑔は重力加速度，ℎは水深を表す．空間近似には三 角形1次要素を用いたガラーキン法，時間近似にはセレクテ ィブ・ランピング法を用いることにより，離散化された式(3) が得られる．すべての要素について足し合わせ，式(4)のよう に簡略的に記述する．ここに，𝛥𝑡は微小時間間隔，ℎ̅は要素 内の平均水深，[𝑆_𝑖]は圧力項に関するマトリックス，[𝑀]は質 量マトリックス，[𝑀̅]は質量マトリックスの各要素を対角項 に集中させた行列，[𝑀̃]は解の安定化のため，ランピングパ

ラメータ e (0≦e≦1)によって式(5)の様に変形した行列であ

る．

{ {𝑢_𝑥^𝑛+1} {𝑢_𝑦^𝑛+1} {𝜂^𝑛+1}

} = [

[[𝑀̅]] [0] −𝑔𝛥𝑡[𝑆_𝑥] [0] [𝑀̅] −𝑔∆𝑡[𝑆_𝑦]

−ℎ̅∆𝑡[𝑆_𝑥] −ℎ̅∆𝑡[𝑆_𝑦] [𝑀̅]

]

−1

{ {𝑢_𝑥^𝑛} {𝑢_𝑦^𝑛} {𝜂^𝑛}

} (3)

{𝜙̂^𝑛+1} = [𝐴]{𝜙̂^𝑛} (4)

[𝑀̃] = (1 − 𝑒)[𝑀̅] + 𝑒[𝑀] (5)

また，カルマンフィルタの理論では，駆動行列（単位行 列）[Γ]にシステムノイズによるベクトル{𝑞^𝑛}を乗じたベク トルを，式(5)に加えることにより真値を表すシステム方程 式となる．また，システムノイズは平均0，分散𝜎²により 表した正規乱数により表されるものとすると，真値の時間 進展の方程式は，式(6)のように与えられる．ここに，行列

[𝐴]は状態遷移行列と呼ばれる．また，観測値{𝑧^𝑛+1}につい

ては真値{𝜙^𝑛+1}に，観測ノイズ{𝑟^𝑛+1}を加えた値により表さ れるものとすると，式(7)に示す観測方程式が得られる．観 測ノイズもシステムノイズと同様に，平均0，分散𝜎²によ り表した正規乱数により表されるものとする．行列[𝐻]は観 測行列と呼ばれ，各観測値に対応する真値の成分に1，それ 以外には0が入る行列とする．

{𝜙^𝑛+1} = [𝐴]{𝜙̂^𝑛} + [Γ]{𝑞^𝑛} (6)

{𝑧^𝑛+1} = [𝐻]{𝜙^𝑛+1} + {𝑟^𝑛+1} (7)

3. 状態推定シミュレーションの計算アルゴリズム

本検討で使用した計算アルゴリズムを以下に示す．(－)と (＋)の添字はそれぞれデータ同化の前後を表す．

1.初期値の入力：

[𝐴], [𝑃₍₊₎⁰ ], {𝜙₍₊₎⁰ }, [Γ], [𝑄], [𝑅], 𝜀, {𝑧^𝑛+1}(𝑛 = 0~𝑖𝑚𝑎𝑥) 2.同化前の誤差共分散行列の計算：

[𝑃(−)] = [𝐴][𝑃(+)][𝐴]^𝑇+ [Γ][𝑄][Γ]^𝑇 3.カルマンゲイン行列の計算：

[𝐾₁] = [𝑃₍₋₎][𝐻]^𝑇([𝐻][𝑃₍₊₎][𝐻]^𝑇+ [𝑅])⁻¹ 4.同化後の誤差共分散行列の計算：

[𝑃(+)] = [𝑃(−)] − [𝐾1][𝐻][𝑃₍₋₎]

5.同化後の誤差共分散行列のフロベニウスノルムによる収束 判定：

もし、 √𝑡𝑟 (([𝑃₍₊₎^𝑘+1] − [𝑃₍₊₎^𝑘 ])([𝑃₍₊₎^𝑘+1] − [𝑃₍₊₎^𝑘 ])^𝑇) ≤ 𝜀 ならステップ6へ，そうでなければステップ2へ 6.支配方程式の有限要素方程式によるステップの更新：

{𝜙̂₍₋₎^𝑛+1} = [𝐴]{𝜙̂₍₊₎^𝑛 } 7.同化後の最適推定値の計算：

_{𝜙̂(+)^𝑛+1_{} = {𝜙̂(−)}^𝑛+1_{} + [𝐾1}]({𝑧^𝑛+1} − [𝐻]{𝜙̂₍₋₎^𝑛+1})

4. カルマンフィルタ FEM による流れ場推定に対する数値実 験

4-1 観測点配置に関する検討

まず，観測点の位置と数について以下のように検討を行う．

本検討では図1のような開流路において，左側の流入境界か ら周期2π秒，振幅1.0mの正弦関数にしたがう水位変動を与 え，システムノイズを考慮して得られた数値解析結果を真値，

さらに観測ノイズを付加したものを人工的な観測値として 取り扱う．ここでシステムノイズ，観測ノイズはそれぞれ平

均0.0，分散0.0001，平均0.0，分散0.1000とした正規乱数と

する．観測点は図2に示すように，Case Aでは流路中央線上 に一様に配置され，Case Bでは流出境界側にのみ配置される．

また，表1に計算条件を整理して示す．

Table1 Computational conditions.

Fig.1 Computational model.

Fig.2 Numerical Cases.

図3および図4に2.0s後の流路中央線上(y=0.4m)における 水位変動量の分布を示す．横軸は流入境界からの距離，縦軸 は水位変動量をそれぞれ表す．各Case とも観測点の配置さ れている区間においては，観測値に付加された観測ノイズが 除去され，推定値は真値に対して良い一致を示している．し かし，Case Bの観測点の配置されていない流入境界から5m までの区間においては水位変動量の推定値が一定となって しまい，水位変動量の推定が困難であることがわかる．

Fig.3 Water elevation on the center line in Case A at T=2.0s.

Fig.4 Water elevation on the center line in Case B at T=2.0s.

4-2 流入境界条件の検討

前節における推定結果より，観測点が流入境界側に配置さ れない場合，最初の観測点までの区間における推定は困難で 時間増分量 ∆𝑡, s. 0.001

ステップ数 2000

節点数 153

要素数 200

重力加速度 𝑔, 𝑚/𝑠² 9.8 ランピングパラメータ 𝑒 0.8 同化後の誤差共分散の初期値 𝑃₍₊₎⁰ 1.0 最適推定値の初期値 𝜙̂₍₊₎⁰ 0

収束判定定数 ε 0.01

あることがわかった．そこで，流入境界条件として観測値を 生成した時と同じ周期2π秒，振幅1.0mの正弦関数にしたが う水位変動を与え，Case Bの観測点配置で再度推定を行う．

その他の計算条件は前節表2に示す同様の値を使用した．こ

の検討をCase Cとし，図5に2.0s後の流路中央線上(y=0.4m)

における水位変動量の分布を示す．図5より，適切な流入境 界条件を与えた場合，Case Bと比べて流入境界側の観測点が 配置されていない区間においても，水位変動の傾向を捉えた 推定が行われていることが見て取れる．この結果から，観測 点を対象領域に一定間隔で配置できない場合や流入境界側 に観測点を配置できない場合においてカルマンフィルタ FEM に基づく流れ場の推定を行う際には，適切な流入境界 条件を設定する必要があると思われる．

Fig.5 Water elevation on the center line in Case C at T=2.0s.

また，真値と観測値を生成した際とは異なる流入境界条件 を与えた場合にはどのように推定が行われるか検証するた

め，同じCase B の観測点配置において，流入境界条件とし

て周期2π秒，振幅1.0mの余弦関数にしたがう水位変動を与 える．その他計算条件は同じとし，この場合をCase Dとし たときの 2.0s 後の流路中央線上(y=0.4m)における水位変動 量の分布を図6に示す．真値と観測値を生成した際とは異な る流入境界条件を与えた場合，流入境界から最初の観測点ま での区間においては，正しく推定が行われない．また，最初 の観測点から流出境界までの区間においては，適切でない流 入境界条件の影響を受けず，良好な推定結果が得られた．

Fig.6 Water elevation on the center line in Case D at T=2.0s.

5. 結言

本研究では，カルマンフィルタ有限要素法に基づき，開流 路における流れ場の推定を行った．平面2次元場の浅水長波 方程式を支配方程式とし，空間方向にはガラーキン法，時間 方向にはセレクティブ・ランピング法を用いて離散化を行っ た．開流路の流れ場の推定精度について流入境界条件に関す る検証を行った．その結果，以下の知見が得られた．

(1)対象領域に一定間隔に観測点を配置した場合(Case A)，観 測値から観測ノイズが取り除かれ，真値に沿う良好な推定結 果が得られた．

(2)観測点を流入境界側に配置しない場合(Case B)，流入境界 から最初の観測点までの区間では水位変動量が一定となり，

流れ場の推定が困難となることがわかった．

(3)観測点を流入境界側に配置しない場合であっても，流入境 界条件が適切に与えられるとき(Case C)，流入境界から最初 の観測点までの区間の推定は概ね良好に行われることがわ かった．

謝辞

本論文を執筆するにあたり，科学研究費補助金（基盤(C)）

15K05786の援助を受けた．ここに謝意を表す．また，本論文

中に示した数値計算を行うにあたり，九州大学情報基盤研究 開 発 セ ン タ ー の 高 性 能 演 算 サ ー バ シ ス テ ム

PRIMERGYCX400を利用させて頂いた．九州大学情報基盤研

究開発センターのスタッフの方々に対して謝意を表す．

参考文献

(1) Kalman, R.E., A New Approach to Linear Filtering and Prediction Problems, Transactions for the ASME – Journal of Basic of Engineering, Vol. 82, No.1(1960), pp.35-45.

(2) Kalman, R. E. and Bucy, R. S., New Results in Linear Filtering and Prediction Theory, Transactions of the ASME - Journal of Basic of Engineering, Vol. 83, No. 1 (1961), pp.95–

108.

(3) Nakamura, M., Tanaka, M., Adachi, O. and Ishikawa, H., Optimum Design of Transient Heat Conduction Fields Using Boundary Element Inverse Analusis, JSME Int Journal. Ser. A, Mechanics and material engineering Vol.

38(1995), No.4, pp.480-486.

(4) 井上藍，宮崎保光，後藤信夫，脳磁図による電流源推定 に関する逆問題の有限要素法:カルマンフィルタ解析，

情 報 科 学 技 術 フ ォ ー ラ ム 一 般 講 演 論 文 集 2002(2), pp.327-328, 2002

(5) Heemink, A. W., Two-dimensional Shallow Water Flow Identification, Applied Mathematical Modelling, Vol.

12(1988), pp.109-118.

(6) 倉橋貴彦，吉荒太一，小林泰秀，カルマンフィルタ有限 要素法による浅水域における流れ場の推定精度の検証

（観測変数・観測点の数および位置が推定精度に与える 影響)，日本機械学会論文集, 第82巻，第835号 pp.1- 19, 2016.