拡張カルマンフィルタFEMに基づく水位状況推定解析

JST新技術説明会 2023年 6月15日（木）倉橋 貴彦
長岡技術科学大学
技学研究院 機械系

1

発明のイメージ図

画像解析後、水位をデータ化
カルマンフィルタFEMによる浸水域の推定解析

カメラによる画像撮影

結果の可視化（※リアルタイムにおける浸水域の高精度把握）

浸水域：撮影点（観測点）（カルマンフィルタFEMの利点）
通常の数値シミュレーションでは，モデルの仮定に従った解は得られるが，観測値に含まれる要因を全て網羅はできない．
しかし，カルマンフィルタFEMでは，観測値に含まれる要因を，観測点のみならず流れ場に反映でき，シミュレーションの結果を観測値に合うように補正することができる．

東京湾モデルの流れ場推定解析の例の紹介

2

1. 拡張カルマンフィルタFEMについて観測値を用いた最適推定解析をする（最適未知パラメータ算定の）ためのアプローチ• 確定論的アプローチ (随伴変数法, 直接微分法 等)*** 問題 : シミュレーション値が観測値に一致するように「未知パラメータ」を算定する．
• 確率論的アプローチ (カルマンフィルタ, 粒子フィルタ 等)***問題：推定誤差を最小とするように「未知パラメータ」を算定する．
カルマンフィルタにおける“n+1”ステップの最適推定値
𝑛+1𝑛+1𝑢(+)= 𝑢(−)+ 𝐾 𝑛+1 𝑧 𝑛+1 − 𝑢𝑛+1−

推定値を算定するための修正項

𝐾 𝑛+1 : カルマンゲイン(※推定誤差を最小とする様に与えられる．)
𝑧 𝑛+1 : 観測値

𝑛+1= (1 − 𝐾 𝑛+1 )𝑢(−)+ 𝐾 𝑛+1 𝑧 𝑛+1u : 状態量(-) : シミュレーションの結果（同化前の状態量）(+) : 観測値zを用いて修正した値（同化後の状態量）

観測値 𝑧 𝑛+1
推定値

𝑛+1𝑢(+)
𝑛+1

シミュレーション値 𝑢(−)“n+1” ステップt : 時間Fig : 確率論的アプローチにおける推定値の計算のイメージ図

検討例の紹介
1. Kalman filterとFEMを融合した手法により，東京湾モデルに対する潮流推定解析を実施し推定精度に関する考察を実施．
2. 潮流発電に対する発電ポテンシャルを算定．

3

2. 有限要素法に基づく支配方程式の離散化
支配方程式(非線形浅水長波方程式)𝑢ሶ 𝑖 + 𝑢𝑗 𝑢𝑖,𝑗 + 𝑔𝜂,𝑖 − 𝜈 𝑢𝑖,𝑗 + 𝑢𝑗,𝑖𝜂ሶ +ℎ + 𝜂 𝑢𝑖 ,𝑖 = 0,𝑗+ 𝑓𝑢𝑖 = 0x2(=y)η(u1,u2)=(ux,uy)

ui : x,y方向流速
η : 水位g : 重力加速度h : 平均水深
ν : 渦動粘性係数
f : 摩擦係数

ガラーキン法
状態変数と重み関数に同様の補間関数を適用(N : 形状関数)

三角形一次要素33𝑢𝑖 𝑥, 𝑦 = ෍ 𝑁𝑘 𝑢𝑖𝑘
𝑘=1

Water surfaceΩeη

3x1(=x)1𝑢𝑖∗ 𝑥, 𝑦2

𝑘=1

h

3

Sea bottom

重み付き残差方程式∗= ෍ 𝑁𝑘 𝑢𝑖𝑘

Assumption:*The current velocity from the water surface to the sea bottom is uniform.
*Water pressure is represented by hydrostatic pressure.

流速ui・水位η
同様の補間関数を使用（同次補間）

𝜂 𝑥, 𝑦 = ෍ 𝑁𝑘 𝜂𝑘𝑘=1
3𝜂 𝑥, 𝑦 = ෍ 𝑁𝑘 𝜂𝑘∗∗𝑘=1

න 𝑢𝑖∗ (𝑢ሶ 𝑖 + 𝑢𝑗 𝑢𝑖,𝑗 + 𝑔𝜂,𝑖 − 𝜈 𝑢𝑖,𝑗 + 𝑢𝑗,𝑖Ω𝑒න 𝜂 ∗ (𝜂ሶ +Ω𝑒,𝑗+ 𝑓𝑢𝑖 ൰ = 0

有限要素方程式ℎ + 𝜂 𝑢𝑖 ,𝑖 )𝑑Ω = 0

𝑀𝑒 𝑈ሶ 𝑒𝑖 + 𝑢ത 𝑒𝑗 𝑆𝑒𝑗 𝑈𝑒𝑖 + 𝑔 𝑆𝑒𝑖 𝐸𝑒 − 𝜈 𝐻𝑒𝑗𝑗 𝑈𝑒𝑖 + 𝐻𝑒𝑖𝑗 𝑈𝑒𝑗+ 𝑓 𝑀𝑒 𝑈𝑒𝑖 = {0}𝑀𝑒 𝐸ሶ 𝑒𝑖 + 𝑢ത 𝑒𝑖 𝑆𝑒𝑖 𝐸𝑒𝑖 + 𝑢ത 𝑒𝑖 𝑆𝑒𝑖 ℎ𝑒𝑖 + (ℎത 𝑒 + 𝜂ҧ𝑒 ) 𝑆𝑒𝑖 𝑈𝑒𝑖 = {0}4支配方程式に対する有限要素方程式
𝑀𝑒 𝑈ሶ 𝑒𝑖 + 𝑢ത 𝑒𝑗 𝑆𝑒𝑗 𝑈𝑒𝑖 + 𝑔 𝑆𝑒𝑖 𝐸𝑒 − 𝜈 𝐻𝑒𝑗𝑗 𝑈𝑒𝑖 + 𝐻𝑒𝑖𝑗 𝑈𝑒𝑗+ 𝑓 𝑀𝑒 𝑈𝑒𝑖 = {0}𝑀𝑒 𝐸ሶ 𝑒𝑖 + 𝑢ത 𝑒𝑖 𝑆𝑒𝑖 𝐸𝑒𝑖 + 𝑢ത 𝑒𝑖 𝑆𝑒𝑖 ℎ𝑒𝑖 + (ℎത 𝑒 + 𝜂ҧ𝑒 ) 𝑆𝑒𝑖 𝑈𝑒𝑖 = {0}

時間方向の離散化：セレクティブ・ランピング法の適用

時間方向に離散化された有限要素方程式
1𝑛
෩ 𝑈𝑒𝑖 𝑛 ) + 𝑢ത 𝑒𝑗( 𝑀𝑒 𝑈𝑒𝑖 𝑛+1 − 𝑀𝑆𝑒𝑗 𝑈𝑒𝑖 𝑛 + 𝑔 𝑆𝑒𝑖 𝐸𝑒 𝑛 − 𝜈Δ𝑡𝐻𝑒𝑗𝑗 𝑈𝑒𝑖 𝑛 + 𝐻𝑒𝑖𝑗 𝑈𝑒𝑗𝑛+ 𝑓 𝑀𝑒 𝑈𝑒𝑖 𝑛 = {0}1𝑛𝑛෩ 𝐸𝑒𝑖 𝑛 ) + 𝑢ത 𝑒𝑖( 𝑀𝑒 𝐸𝑒𝑖 𝑛+1 − 𝑀𝑆𝑒𝑖 𝐸𝑒𝑖 𝑛 + 𝑢ത 𝑒𝑖
𝑆𝑒𝑖 ℎ𝑒𝑖 𝑛 + ℎത 𝑒𝑛 + 𝜂ҧ𝑒𝑛 𝑆𝑒𝑖 𝑈𝑒𝑖 𝑛 = {0}Δ𝑡整合質量行列, 集中質量行列

2 1 1 
1 0 0
Ae Ae M  = 1 2 1, M = 0 1 0
123
1 1 20 0 1~M = (1 − e ) M + eM    e: ランピングパラメータ （0 ≦ e≦1)

有限要素方程式

𝜙෠ 𝑛+1 = 𝐴𝑛 𝜙෠ 𝑛 + 𝑁 𝑛5

3. Kalman filter 理論において使用する方程式
有限要素方程式（推定値の計算）
𝑛+1𝑛𝜙෠(−)= 𝐴𝑛 𝜙෠(+)+ 𝑁𝑛𝑛+1𝜙෠(+)= 𝜙෠ 𝑛+1+ 𝐾1𝑛+1 ( 𝑧 𝑛+1 − 𝐻 𝜙෠ 𝑛+1)−−

システム方程式（真値の計算）

システムノイズ

𝜙 𝑛+1 = 𝐴𝑛 𝜙 𝑛 + 𝑁 𝑛 + Γ 𝑞𝑛
駆動行列

状態遷移行列
観測方程式

推定値と真値の関係

𝜙෠ = 𝜙 + 𝑝
推定値

z = H  + r n +1(-)：同化前（観測値による修正前）(+)：同化後（観測値による修正後）

n +1

観測行列

真値

推定誤差

n +1

推定値は，真値と推定誤差の和により表される．

観測ノイズ

Kalman filter理論に基づく定式化

同化後の推定誤差を最小とするように「Kalman gain行列」を誘導

6

Kalman gain行列の誘導 （※展開の流れ程度）
同化後の推定値

推定値と真値の関係

𝑛+1 + 𝑝 𝑛+1
𝜙෠ 𝑛+1=𝜙(+)+移項𝑛+1𝑝(+)= 𝜙෠ 𝑛+1− 𝜙 𝑛+1+𝑛+1：同化後の推定誤差
𝑝(+)

同化後の推定値の計算式，観測方程式を代入

推定値

真値

推定誤差

推定値は，真値と推定誤差の和により表される．

= 𝜙෠ 𝑛+1+ 𝐾1𝑛+1 ( 𝑧 𝑛+1 − 𝐻 𝜙෠ 𝑛+1) − 𝜙 𝑛+1−−=

𝐼 − 𝐾1𝑛+1 𝐻

𝑝𝑛+1+ 𝐾1𝑛+1 𝑟 𝑛+1−

同化後の推定誤差の共分散行列
𝑃𝑛+1=< 𝑝𝑛+1
𝑝𝑛+1+++𝑇※同化後の推定誤差{p(+)n+1}をカルマンゲイン行列[K1n+1]の関数として表示>
同化後の推定誤差を最小とするように「Kalman gain行列」を誘導𝜕𝑡𝑟𝜕 𝐾1𝑛+1

𝑃𝑛+1+= −2 𝑃𝑛+1
𝐻 𝑇 + 2 𝐾1𝑛+1 𝐻 𝑃𝑛+1𝐻 𝑇 + 2 𝐾1𝑛+1 [𝑅𝑛+1 ]=[0]−−𝐾1𝑛+1𝑃𝑛+1−𝐻𝑇=

𝐻𝑃𝑛+1−

𝐻 𝑇+𝑅𝑛+1−17拡張Kalman filter FEMの計算の流れ（※非線形のシステム方程式を用いたKalman filter）[𝐹 𝑛 ] =𝜕 𝐴 𝜙෠ 𝑛 + 𝑁 𝑛อ𝜕 𝜙𝑛0
1 初期条件 {𝜙෠ (+) } およびその他計算条件の入力

𝜙𝑛 ={0}線形化した状態遷移行列

2 推定誤差共分散行列の計算

nn𝑃𝑛+1=𝐹𝑃+−

3 Kalman gain 行列の計算

𝐾1𝑛+14 予測誤差共分散行列の計算

𝑛+1𝑛+1𝑛+1𝑃𝑛+1=𝑃−𝐾[𝐻]𝑃1+−−

5 推定値の計算

𝑛n 𝜙𝑛෠𝜙෠ 𝑛+1=𝐴−+ + {𝑁 }6 最適推定値の計算

𝜙෠ 𝑛+1= 𝜙෠ 𝑛+1+ 𝐾1𝑛+1 ( 𝑧 𝑛+1 − 𝐻 𝜙෠ 𝑛+1)+−−=𝑃𝑛+1−

𝐹 n 𝑇 + 𝛤 𝑄 𝑛 [𝛤]𝐻𝑇𝐻𝑃𝑛+1−

𝐻 𝑇+𝑅𝑛+1−18

4. 数値解析例

東京湾モデル （CASE A）
東京
千葉

東京湾

横浜

https://ja.wikipedia.org/wiki/%E6%
9D%B1%E4%BA%AC%E6%B9%BE

・Position of observation points, and observation variable

横須賀

：観測点

Inflow boundary

Fig. Position of observation points and points for verification .

・東京・千葉・横浜

観測値： 検潮所における水位の履歴： 推定点
観測値と推定値を比べ、推定精度を検証する点

目的：
拡張Kalman filter FEMにより潮流推定解析を実施し，推定精度に関する考察を実施
・「線形浅水長波方程式，Kalman filter FEM」，「非線形浅水長波方程式，拡張Kalman filter FEM」の比較
・観測点位置が推定精度に与える影響
・潮流発電ポテンシャルの算定

9

計算条件

節点数; 1052
要素数; 1731

有限要素メッシュ図

平均水深の分布図
Fig. Finite element mesh and distribution of water depth.

10境界条件 および その他計算条件
陸境界（スリップ境界条件）

𝑢𝑠𝑙𝑖𝑝
cos 2 𝜃=𝑣𝑠𝑙𝑖𝑝
sin 𝜃 cos 𝜃cos 𝜃 sin 𝜃 𝑢𝑣
sin2 𝜃

𝑢𝑠𝑙𝑖𝑝
𝑣𝑠𝑙𝑖𝑝𝑢𝑣𝜃Flow velocity for X direction after slip treatment [m/s]Flow velocity for Y direction after slip treatment [m/s]Flow velocity for X direction before slip treatment [m/s]Flow velocity for Y direction before slip treatment [m/s]Tangent angle on land boundary [deg.]

Tab. 計算条件
流入境界条件主要４分潮4𝜂 𝑡 = ෍ 𝑎𝑖 sin
𝑖=12𝜋𝑡+ 𝜅𝑖
𝑇𝑖

Tab. Tidal harmonic constants for major four tidal components in Tateyama.

𝑖

各成分

振幅 𝑎𝑖 , m

位相差 𝜅𝑖 , deg.

周期 𝑇𝑖

1

M20.3605

146.4812h25m

2S20.1707

174.3612h00m

3K10.2318

177.1325h49m

4O10.1823

158.4023h50m※e=0.8と設定．

11データ同化の計算結果ランピングパラメータ e=0.8

観測値と推定値の比較
千葉

東京
観測値

東京

千葉

推定値

横浜

横浜

Inflow boundary

12横須賀における観測水位と推定水位の比較（CASE A）

観測値

ランピングパラメータ e=0.8

推定値（非線形浅水長波eq. + EKF FEM)

東京

千葉

横浜

横須賀
推定値
（線形浅水長波eq. + KF FEM)

Inflow boundary

推定値（システム方程式：線形浅水長波eq. )に比べて，推定値（システム方程式：非線形浅水長波eq. + EKF FEM)の方が観測値との位相差を小さくすることができている．（推定値の振幅には大きな差は見受けられない．）

13

Case A：流速ベクトルおよび水位の分布図 （EKF-FEMによる推定結果）
ランピングパラメータ e=0.8

T=4.0h

T=8.0h

T=12.0h

T=16.0h

T=20.0h

T=24.0h14

東京における観測水位と推定水位の比較（CASE B）（観測点：千葉，横浜，横須賀， 確認点：東京 とした場合）
推定値（非線形浅水長波eq. + EKF FEM)

観測値

確認点
東京

2
Observed vlue

Water elevation, m

1.5

観測点 No.1
千葉

Estimated value(KF-FEM)

Estimated value(EKF-FEM)

1

観測点 No.2
横浜

0.5
0-0.5
-1-1.504

8121620Time, h

24

観測点 No.3横須賀

推定値（線形浅水長波eq. + KF FEM)

推定値（システム方程式：線形浅水長波eq. )に比べて，推定値（システム方程式：非線形浅水長波eq. + EKF FEM)の方が観測値との振幅の一致性が高いことを確認できる．（位相差には大きな差は見受けられない．）

15

Case B：流速ベクトルおよび水位の分布図 （EKF-FEMによる推定結果）
ランピングパラメータ e=0.8

T=4.0h

T=8.0h

T=12.0h

T=16.0h

T=20.0h

T=24.0h16

ランピングパラメータ e=0.8
Case A, B：流速ベクトルおよび水位の分布図の比較（EKF-FEMによる推定結果）（※左：Case A（観測点：東京，千葉，横浜）， 右：Case B （観測点：千葉，横浜，横須賀） ）T=4.0h

T=8.0h

T=12.0h

T=16.0h

T=20.0h

T=24.0h17従来技術 （従来技術における問題点）
参考：「線形」浅水長波eq.を用いた順解析結果（通常のFEM）の結果と観測値の比較
東京千葉

ランピングパラメータ e=0.8

観測値

通常のFEMによる結果（ 「線形」浅水長波eq. ）横浜横須賀

Yokosuka

観測値を用いたデータ同化解析（最適推定シミュレーション）を実施することにより，通常の流れ解析に比べより実測に近いシミュレーションを実施できる．

18想定される用途：潮流発電関連のニュース

環境省「令和４年度潮流発電による地域の脱炭素化モデル構築事業」の一環として長崎海洋産業クラスター形成推進協議会（長崎県長崎市）と共同で実施するもので、発電開始は2024年度を予定している。
https://project.nikkeibp.co.jp/ms/atcl/19/news/00001/02436/?ST=msb#:~:text=
https://energy-shift.com/news/bba31792-6f5c-4d99-8fbd-e96501c4d95f

19

潮流発電ポテンシャルの算定

計算式（※風力発電のエネルギー量の計算式）

表：各パラメータの設定

Sweep area A [m2]Parameter η T

1.0Time increment dt [s]1.0Time step imax

86400

3Density ρ [kg/m ]1.01.0※最終的にポテンシャル分布を最大値により正規化するため，各種パラメータは1.0と与えている．
（※エネルギー変換効率の係数 ηT<1)

観測点位置と推定点位置

20

東京湾における発電ポテンシャルの算定結果

東京湾における発電ポテンシャルの分布(カルマンフィルタ,24h)(正規化処理有り)

東京湾における発電ポテンシャルの分布(拡張カルマンフィルタ,24h)

21■ 実用化に向けた課題

・適切な観測点個数，適切な観測点配置に関する調査・解析時間の短縮（現状1日分を対象とした解析において，スーパーコンピュータを用いたとしても，1日以上の計算時間を要している．（メッシュの粗さや時間刻みにも依るが，行列と行列の積等，行列の演算も多いため，多くの計算時間を要する．））■ 企業様への期待
・「水域の流れ場予測を高精度に行いたい」と考えています企業様には，本技術の導入が有効であると考えられる．

22

本技術に関する知的財産権•発明の名称：水位状況推定プログラム、水位状況推定システムおよび水位状況推定方法•出願番号 ：特願2021- 52765
•出願人
•発明者

：国立大学法人長岡技術科学大学
独立行政法人国立高等専門学校機構：倉橋貴彦、剱地利昭，五十嵐 晃平 23

お問い合わせ先
長岡技術科学大学
産学連携・研究推進課
知的財産係

ご清聴頂きましてありがとうございました。

ＴＥＬ 0258-47-9279
ＦＡＸ 0258-47-9040
e-mail patent@jcom.nagaokaut.ac.jp
24