カルマンフィルタ有限要素法を用いた河川の解析に関する研究

(1)

修士論文発表要旨

(2008

年度)

カルマンフィルタ有限要素法を用いた河川の解析に関する研究

Estimation of River Current Using Reduced Kalman Filter Finite Element Method

土木工学専攻１０号尾島康則

Yasunori OJIMA

1

はじめに

近年,水質汚染は非常に深刻な環境問題になっている.

水質汚染とは火山・海藻の異常発生や暴風雨及び地震と同様に水質と水系の生態系の状態における大きな変質を招くものなのだが,河川の水質がある基準を超えた時にこれを汚染という. そのため,人々が生活をするに当たり, 河川の状態を把握することは大変重要である. その状態を把握するために,様々な観測が行われるが,必要な観測データの入手が困難であることや, 観測する際に機械誤差, 人為誤差や自然誤差など様々な誤差が含まれるという問題点がある. そこで,本研究では有限要素法を用いたカルマンフィルタ理論を適用し,モデルの不確実性を確率過程に基づいて修正し,推定を行う. カルマンフィルタ理

論

^[1][2]

とは,確率過程に基づいた解析手法の一つであり

誤差を含む観測値から,その誤差を考慮しながら推定,同定,制御,予測,データ同化等様々な分野で用いられるフィルタリング理論の一つである.そして,このカルマンフィルタ理論と有限要素法を組み合わせることによって,時系列のみならず空間方向にも解析が可能となる.

本研究では,数値解析例として,日本でも有数の汚染河川である千葉県手賀沼を実在のモデルとして扱った. 特に河川の汚染を表す指標である溶存酸素量の推定に非常に影響をもつ拡散係数に着目し,拡散係数を同定し,同時に推定を行った. 同定をするためにフィルタリング理論である拡張カルマンフィルタと有限要素法を組み合わせた拡張カルマンフィルタ有限要素法を用いる. また, 多くの解析時間及び計算容量がかかってしまうという拡張カルマンフィルタ有限要素法の問題の削減を目的とした, 縮小カルマンフィルタ有限要素法を提案する. 新手法である縮小カルマンフィルタ有限要素法を河川に適用し,有効性を検証する. また有限要素法の基礎方程式として浅水長波方程式を用いる.

2

数値解析手法

2.1

カルマンフィルタ

カルマンフィルタとは, 確率過程に基づいた推定手法の一つである. 本手法を用いることによりモデルの不確

実性を確率過程に基づいて修正し, 状態量を逐次決定的に推定することができる. 同定問題は非線形問題となるため,カルマンフィルタを拡張する必要がある. 拡張カルマンフィルタの状態空間モデル,システム方程式

(1)

と観測方程式

(2)

をそれぞれ以下に示す.

x _k+1 = f _k (x _k ) + G _k w _t (1) y _k = h _k x _k + v _k (2)

ここで

x _k

は状態量,

f _k (x _k )

は状態遷移行列,

G _k

は駆動行列,

y _k

は観測値,

h _k

は観測行列である.

w _k

及び

v _k

はシステムノイズ及び観測ノイズと呼ばれ, 平均

0,

共分散行列に従う正規白色ノイズである.

x _k

の関数で表される状態遷移行列

f _k (x _k )

と観測行列

h _k (x _k )

を線形化するために以下に示す高次の項を無視したテイラー展開を適用する.

f _k (x _k ) f _k (x ^∗ _k/k ) + F _k (x _k − x ^∗ _k/k ) (3) h _k (x _k ) h _k (x ^∗ _k/k−1 ) + H _k (x _k − x ^∗ _k/k−1 ) (4)

この行列をそれぞれ式

(1), (2)

へ代入すると以下の線形化された状態空間モデルを得ることができる.

x _k+1 = F _k x _k + G _k w _k + f _k (x ^∗ _k/k ) − F _k x ^∗ _k/k (5) y _k = H _k x _k + v _k + h _k (x ^∗ _k/k−1 ) − H _k x ^∗ _k/k−1 (6)

上式の

F _k ,H _k

は感度行列と呼ばれ,以下のように表される.

F _k = ∂f _k

∂x _k

x k =x ^∗ _k/k

, H _k = ∂h _k

∂x _k

x k =x ^∗ _k/k−1

ここで

∂x _k/k

は観測データ

y ₀ · · · y _t

が与えられたときの

x _k/k

の最小分散推定値である. 以下にカルマンフィルタのアルゴリズムを示す.

1. [ˆΓ _0|−1 ] = [Γ _x ₀ ] , { x ˆ _0|−1 } = { ¯ x ₀ }

2. [K _k ] = [ˆΓ _k|k−1 ][H _k ] ^T ([R _v _k ] + [H _k ][ˆΓ _k|k−1 ][H _k ] ^T ) ⁻¹ 3. [ ˆ P _k|k ] = ([I] − [K _k ][H _k ])[ˆΓ _k|k−1 ]

4. [ˆΓ _k+1|k ] = [F _k ][ ˆ P _k|k ][F _k ] ^T + [G _k ][Q _w _k ][G _k ] ^T 5. { x ˆ _k+1|k } = [F _k ] { x ˆ _k|k }

6. { x ˆ _k|k } = { x ˆ _k|k−1 } + [K _k ]( { y _k } − [H _k ] { x ˆ _k|k−1 } )

(2)

ここで

K _k

はカルマンゲイン, ˆ

P _k|k

は推定誤差共分散行列, ˆΓ

_k+1|k

は予測誤差共分散行列,

Q _wk , R _vk

は共分散行列である.

2.2

縮小カルマンフィルタ有限要素法

前述した拡張カルマンフィルタには物理モデルが含まれていない, そこで有限要素法で近似した有限要素方程式を物理モデルとして拡張カルマンフィルタに組み込むことで拡張カルマンフィルタ有限要素法を導く. カルマンフィルタに物理モデルを組み込む方法は大別して

2

種類ある. 一つは状態方程式に組み込む方法,もう一つは観測方程式に組み込む方法である. 前者の方は,時系列過程を追従することができるので, 前述したカルマンフィルタのアルゴリズムをそのまま用いることができる.

本手法を拡張カルマンフィルタ有限要素法と呼ぶ. 本手法を用いることにより, 状態量の時間的分布のみにならず空間的分布においても推定が可能となる. また,本研究では共分散行列に

4

次スプライン関数,

S = ¯ S{1 . 0 − 6 . 0

„ r i

d m

« ₂ + 8 . 0

„ r i

d m

« ₃

− 3 . 0

„ r i

d m

« ₄ } (7)

を提案する. ここで,Sは共分散の空間分布を,

r _i , d _m

は節点間距離,最大節点間距離を表している. また,観測誤差

v _k ,

システム誤差

w _k

は観測値と推定値の差,最適推定値と推定値の差によってそれぞれ求めることができる.

共分散の空間分布を節点周りに定義することによってカルマンフィルタ有限要素法が持つ解析時間及び解析容量が増加する問題点を解決することが出来る.そこで本研究では縮小カルマンフィルタ有限要素法を提案し, 解決する. 観測点付近の情報のみでカルマンフィルタ理論を適用するために式

(1)

を変化させ,

x _k+1 z _k+1

=

F _k L _k M _k N _k

x _k z _k

+

Gw _k

0 (8)

このように, 観測点付近とそれ以外に分離させる. これによりカルマンゲイン

K

を求めるために用いる行列が

F _k

のみとなり計算要領と時間の削減を行うことが可能となる.

dm ri

main point related point

図

1;

節点及び節点周りの領域

0 0.2 0.4 0.6 0.8 1

Quartic Spline Function

図

2; 4

次スプライン関数

2.3

状態方程式

流れの状態を表す方程式として以下に示す非線形浅水長波方程式と移流拡散方程式を用い, Galerkin法を適用することにより,有限要素近似を行う.

＜非線形浅水長波方程式＞

u ˙ _i + u _j u _i,j + g(η + H + h), _i −ν(u _i,j + u _j,i ) _,j = 0 (9) η ˙ + { (η + H)u _i } _,i = 0, (10)

＜移流拡散方程式＞

c ˙ + u _i c, _i −κc, _ii = 0 (11)

ここで,u

_i

は

x,y

方向の流速,gは重力加速度,ηは水位変動量,H は基準水深,hは河床勾配,νは渦動粘性係数,cは溶存酸素量,κは拡散係数を表している. そして空間方向の離散化に有限要素法を用い,時間方向の離散化に陽的

Euler

法を用いることで,以下の有限要素方程式を得る.

M ¯ _αiβj u ⁿ⁺¹ _βj = ˜ M _αiβj u ⁿ _βj − ∆t(K _αβγj u ⁿ _βj u ⁿ _γi (12) +D _αiβj u ⁿ _βj − H _αiλ η _λ ⁿ + Ω ⁿ _αi )

M ¯ _λµ η ⁿ⁺¹ _µ = ˜ M _λµ η _µ ⁿ − ∆t(A _λµβj u ⁿ _βj η ⁿ _µ + B _λβj u ⁿ _βj ) (13)

M ¯ _αβ c ⁿ⁺¹ _β = ˜ M _αβ c ⁿ _β − ∆t(K _αβγj c ⁿ _β u ⁿ _γj (14) +D _αβ c _β + Ψ _α )

ここで,質量集中行列

M ¯ _αβ

と

selective lumping parame- ter

を用いた混合質量行列

M ˜ _αβ

は以下のように表される.

M ¯ _αβ = ∆ 3

⎡

⎢ ⎢

⎣

1 0 0 0 1 0 0 0 1

⎤

⎥ ⎥

⎦ (15)

M ˜ _αβ = e M ¯ + (1 − e)M (16)

ここで,

e

は

lumping parameter

で本研究では

0.9

を用いている.

(3)

2.4

同時推定同定手法

溶存酸素量に非常に関係のある拡散係数のパラメータ同定を行う事で最適な拡散係数を導出し, 得られた拡散係数の値で推定を同時に行う新しい手法を提案した. 同時推定同定手法として, 定数である拡散係数を以下のように展開を行う.

κ ^k+1 _m = κ ^k _m + κ ^k _m,β (c ^k _β − c ^k+1 _β ) + · · · (17)

次に,展開して得られた式

(17)

を式

(14)

に代入し展開することで以下の式を得ることができる.

» M ¯ αβ 0 0 E mn

– j c β

κ m

ﬀ _k+1

= »

H αβ G αn

A βn E mn

– j c β

κ n

ﬀ _k

− B βm

ここで,H

_αβ

は状態遷移行列,G

_αn

は拡散行列,A

_βn

は感度行列,E

_mn

は単位行列,

B _βm

は既知項である. 本研究では以上の式を用いて,推定と同定を同時に行った.

2.5

新アルゴリズム

以上に基づいて前述したアルゴリズムの改良を以下に示す.

⎧ ⎪

⎪ ⎪

⎪ ⎨

⎪ ⎪

⎩ u v η c

⎫ ⎪

⎪ ⎪

⎪ ⎬

⎪ ⎪

⎭

k+1

= [D _k ]

⎧ ⎪

⎪ ⎪

⎪ ⎨

⎪ ⎪

⎩ u v η c

⎫ ⎪

⎪ ⎪

⎪ ⎬

⎪ ⎪

⎭

k

+ { A _k } , (18)

ここで,D

_k

を式

(8)

のように

4

分割することで状態遷移行列

F _k

を得ることができる. また

A _k

は既知項として扱う. 以下に新アルゴリズムを示す.

＜縮小カルマンフィルタ有限要素法＞

1 . [ˆ Γ _0|−1 ] = [Γ _x ₀ ] , { ˆ x 0|−1 } = { x ¯ 0 }

2 . [ K k ] = [ˆ Γ _k|k−1 ][ H k ] ^T ([ R v _k ] + [ H k ][ˆ Γ _k|k−1 ][ H k ] ^T ) ⁻¹ 3 . [ ˆ P k|k ] = ([ I ] − [ K k ][ H k ])[ˆ Γ _k|k−1 ]

4 . [ˆ Γ _k+1|k ] = [ F k ][ ˆ P k|k ][ F k ] ^T + [ G k ][ Q w k ][ G k ] ^T 5 . { x ˆ k+1|k } = [ F k ] { x ˆ k|k } + {M k }

6 . {M k } = {A k−1 } + {C k }

7 . {C k } = {L k ( M k−1 x k−1 + N k−1 z k−1 + V k−1 ) } + {U k } 8 . { x ˆ k|k } = { x ˆ k|k−1 } + [ K k ]( {y k } − [ H k ] { x ˆ k|k−1 } )

3

数値解析例

3.1

解析モデル

数値解析例として, カルマンフィルタ有限要素法を用いて千葉県手賀沼

(図 3)

における流況を推定する. モデルである千葉県手賀沼は全長

4km,

面積

6.5km ² ,

最大水

深

3.8m,

平均水深

0.9m

である. 解析領域である有限要

素分割図は有限要素分割図は総節点数

1316,

総要素数は

2378

であり,図

4

に示す. 本研究では

2003

年

8

月

16

日

における流速

u, v

と水位変動量

η

及び溶存酸素量

c

を観測値として用いた.

図

3;

手賀沼全景

提供;国土交通省

No.2 No.1

No.3 No.4

図

4;

有限要素分割図

流入条件に観測地点

1

番で得られた観測データを与えました. 図

5〜8

に観測地点

1

番における流速

u, v

と水位変動量

η

c

の観測値をそれぞれ示す.

-0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3

4 8 12 16 20 24

X-Velocity

Time Observed Data

図

5; x

方向の流速

-0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2

4 8 12 16 20 24

Y-Velocity

Time Observed Data

図

6; y

方向の流速

1.24 1.26 1.28 1.3 1.32 1.34 1.36 1.38

4 8 12 16 20 24

Water Elevation

Time Observed Data

図

7;

水位変動量

4 4.5 5 5.5 6 6.5 7

4 8 12 16 20 24

Dissolved Oxygen

Time Observed Data

図

8;

溶存酸素量提供;利根川下流河川事務所また,3日前からの観測データを解析の初期条件とする. 図

9

及び

10

に解析で用いた流速及び溶存酸素量の初期条件を示す.

(4)

X (m)

Y(m)

500 1000 1500 2000

200 400 600 800

0.1 m/s

図

9;

流速の初期条件

X (m)

Y(m)

500 1000 1500 2000

200 400 600 800

C (ppm) 8 7.375 6.75 6.125 5.5 4.875 4.25 3.625 3 2.375 1.75 1.125 0.5

図

10;

溶存酸素量の初期条件

3.2

解析

1(推定)

解析

1

では,本来一定でない定数である,拡散係数

κ

を

0.2

としてカルマンフィルタ有限要素法を用た河川の流況の推定を行った. 観測点

1

番,2番及び

4

番の観測値を用いて,観測値

3

番の流況を自作のプログラムにより推定した. その結果として,図

11〜14

に推定点

(観測地点 3

番) における流速

u, v

と水位変動量

η

c

の推定値の時刻暦と観測値の時刻暦との比較をそれぞれ示す.

また表

1

より, 観測値とほぼ同等の推定値が得られていることがわかる.

-0.03 -0.02 -0.01 0 0.01 0.02 0.03

4 8 12 16 20 24

(M/SEC)

8/16 estimated velocity observed velocity

図

11; x

方向の流速

-0.03 -0.02 -0.01 0 0.01 0.02 0.03

4 8 12 16 20 24

(M/SEC)

8/16 estimated velocity observed velocity

図

12; y

方向の流速

0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84

4 8 12 16 20 24

(M)

8/16 estimated water elevation observed water elevation

図

13;

水位変動量

3 4 5 6 7 8 9 10

4 8 12 16 20 24

Contaminant Concentration (ppm)

Time (hour) identified value (k=0.2) observed data at No.3

図

14;

溶存酸素量

相関係数

x

方向の流速

8.3 y

方向の流速

8.5

水位変動量

9.2

溶存酸素量

8.5

表

1;

相関係数

3.3

解析

2(推定及び同定)

解析

1

では,本来一定でない値であるはずの拡散係数

κ

を一定として解析を行ってしまっている. この問題を解決するために解析

2

として, カルマンフィルタ有限要素法の同定手法を用いることで拡散係数の同定を行い, 同定を行うと同時に推定を行うことでより実際の現象に近い溶存酸素量の推定を行った. 拡散係数の同定結果を図

15

に,同定された拡散係数を用いた溶存酸素量の推定結果を図

16

にそれぞれ示す. また表

2

より,解析

1

より良い精度の解析を行えた.

0 0.05 0.1 0.15 0.2 0.25 0.3

4 8 12 16 20 24

Diffusion Coefficient

Time (hour) identified diffusion coefficient

図

15;

拡散係数

5 6 7 8 9

4 8 12 16 20 24

Contaminant Concentration (ppm)

Time (hour) estimated contaminant observed contaminant

図

16;

溶存酸素量

相関係数溶存酸素量

9.5

表

2;

相関係数

4

まとめ

本研究では河川における,流況の推定及び同定手法として拡張カルマンフィルタ有限要素法を提案し,その妥当性の検証を観測データとの比較によって確認した. 誤差を含む観測データを用いて推定と同定を同時に行い, 拡散係数を変数として扱うことでより実際の現象に近い解析を行った. また,今までカルマンフィルタ有限要素法を扱う上で一番の問題となっていた解析時間と解析容量の大幅な削減を縮小カルマンフィルタ有限要素法を用いることで可能とした. 今後の課題として,より大規模な実在モデルへの適用と, 未来を推定する予測問題へ発展させることが挙げられる.

カルマンフィルタ有限要素法を用いた河川の解析に関する研究

(2008

カルマンフィルタ有限要素法を用いた河川の解析に関する研究

Estimation of River Current Using Reduced Kalman Filter Finite Element Method

Yasunori OJIMA

1

[1][2]

2

2.1

(1)

(2)

x k+1 = f k (x k ) + G k w t (1) y k = h k x k + v k (2)

x k

f k (x k )

G k

y k

h k

w k

v k

0,

x k

f k (x k )

h k (x k )

f k (x k ) f k (x ∗ k/k ) + F k (x k − x ∗ k/k ) (3) h k (x k ) h k (x ∗ k/k−1 ) + H k (x k − x ∗ k/k−1 ) (4)

(1), (2)

x k+1 = F k x k + G k w k + f k (x ∗ k/k ) − F k x ∗ k/k (5) y k = H k x k + v k + h k (x ∗ k/k−1 ) − H k x ∗ k/k−1 (6)

F k ,H k

F k = ∂f k

∂x k

x k =x ∗ k/k

, H k = ∂h k

∂x k

x k =x ∗ k/k−1

∂x k/k

y 0 · · · y t

x k/k

1. [ˆΓ 0|−1 ] = [Γ x 0 ] , { x ˆ 0|−1 } = { ¯ x 0 }

2. [K k ] = [ˆΓ k|k−1 ][H k ] T ([R v k ] + [H k ][ˆΓ k|k−1 ][H k ] T ) −1 3. [ ˆ P k|k ] = ([I] − [K k ][H k ])[ˆΓ k|k−1 ]

4. [ˆΓ k+1|k ] = [F k ][ ˆ P k|k ][F k ] T + [G k ][Q w k ][G k ] T 5. { x ˆ k+1|k } = [F k ] { x ˆ k|k }

6. { x ˆ k|k } = { x ˆ k|k−1 } + [K k ]( { y k } − [H k ] { x ˆ k|k−1 } )

K k

P k|k

k+1|k

Q wk , R vk

2.2

2

4

S = ¯ S{1 . 0 − 6 . 0

„ r i

d m

« 2 + 8 . 0

„ r i

d m

« 3

− 3 . 0

„ r i

d m

« 4 } (7)

r i , d m

v k ,

w k

(1)

x k+1 z k+1

=

F k L k M k N k

x k z k

+

Gw k

0

(8)

K

F k

1;

2; 4

2.3

u ˙ i + u j u i,j + g(η + H + h), i −ν(u i,j + u j,i ) ,j = 0 (9) η ˙ + { (η + H)u i } ,i = 0, (10)

c ˙ + u i c, i −κc, ii = 0 (11)

i

x,y

Euler

^[1][2]

x _k+1 = f _k (x _k ) + G _k w _t (1) y _k = h _k x _k + v _k (2)

x _k

f _k (x _k )

G _k

y _k

h _k

w _k

v _k

x _k

f _k (x _k )

h _k (x _k )

f _k (x _k ) f _k (x ^∗ _k/k ) + F _k (x _k − x ^∗ _k/k ) (3) h _k (x _k ) h _k (x ^∗ _k/k−1 ) + H _k (x _k − x ^∗ _k/k−1 ) (4)

x _k+1 = F _k x _k + G _k w _k + f _k (x ^∗ _k/k ) − F _k x ^∗ _k/k (5) y _k = H _k x _k + v _k + h _k (x ^∗ _k/k−1 ) − H _k x ^∗ _k/k−1 (6)

F _k ,H _k

F _k = ∂f _k

∂x _k

x k =x ^∗ _k/k

, H _k = ∂h _k

∂x _k

x k =x ^∗ _k/k−1

∂x _k/k

y ₀ · · · y _t

x _k/k

1. [ˆΓ _0|−1 ] = [Γ _x ₀ ] , { x ˆ _0|−1 } = { ¯ x ₀ }

2. [K _k ] = [ˆΓ _k|k−1 ][H _k ] ^T ([R _v _k ] + [H _k ][ˆΓ _k|k−1 ][H _k ] ^T ) ⁻¹ 3. [ ˆ P _k|k ] = ([I] − [K _k ][H _k ])[ˆΓ _k|k−1 ]

4. [ˆΓ _k+1|k ] = [F _k ][ ˆ P _k|k ][F _k ] ^T + [G _k ][Q _w _k ][G _k ] ^T 5. { x ˆ _k+1|k } = [F _k ] { x ˆ _k|k }

6. { x ˆ _k|k } = { x ˆ _k|k−1 } + [K _k ]( { y _k } − [H _k ] { x ˆ _k|k−1 } )

K _k

P _k|k

_k+1|k

Q _wk , R _vk

« ₂ + 8 . 0

« ₃

« ₄ } (7)

r _i , d _m

v _k ,

w _k

x _k+1 z _k+1

F _k L _k M _k N _k

x _k z _k

Gw _k

F _k

u ˙ _i + u _j u _i,j + g(η + H + h), _i −ν(u _i,j + u _j,i ) _,j = 0 (9) η ˙ + { (η + H)u _i } _,i = 0, (10)

c ˙ + u _i c, _i −κc, _ii = 0 (11)

_i

M ¯ _αiβj u ⁿ⁺¹ _βj = ˜ M _αiβj u ⁿ _βj − ∆t(K _αβγj u ⁿ _βj u ⁿ _γi (12) +D _αiβj u ⁿ _βj − H _αiλ η _λ ⁿ + Ω ⁿ _αi )

M ¯ _λµ η ⁿ⁺¹ _µ = ˜ M _λµ η _µ ⁿ − ∆t(A _λµβj u ⁿ _βj η ⁿ _µ + B _λβj u ⁿ _βj ) (13)

M ¯ _αβ c ⁿ⁺¹ _β = ˜ M _αβ c ⁿ _β − ∆t(K _αβγj c ⁿ _β u ⁿ _γj (14) +D _αβ c _β + Ψ _α )

M ¯ _αβ

M ˜ _αβ

M ¯ _αβ = ∆ 3

M ˜ _αβ = e M ¯ + (1 − e)M (16)

κ ^k+1 _m = κ ^k _m + κ ^k _m,β (c ^k _β − c ^k+1 _β ) + · · · (17)

ﬀ _k+1

ﬀ _k

_αβ

_αn

_βn

_mn

B _βm

= [D _k ]

+ { A _k } , (18)

_k