修士論文発表要旨
(2010
年度)カルマンフィルタ有限要素法を用いた河川の流れ及び COD の解析に関する研究
Parameter Identification of River Current Using Reduced Kalman Filter Finite Element Method
土木工学専攻
14
号 金井 靖Osamu KANAI
1
はじめに近年,水質汚染は非常に深刻な環境問題になっている.
水質汚染とは火山・海藻の異常発生や暴風雨及び地震と 同様に水質と水系の生態系の状態における大きな変質を 招くものなのだが,河川の水質がある基準を超えた時にこ れを汚染という. そのため,人々が生活をするに当たり, 河川の状態を把握することは大変重要である. その状態 を把握するために,有限要素法は解析手段の一つとして非 常に有効である，しかしながら計算に用いる際にはある 程度、実際の観測値が必要となってくる，様々な観測が 行われるが, 必要な観測データの入手が困難であること や,観測する際に機械誤差,人為誤差や自然誤差など様々 な誤差が含まれるという問題点がある. そこで,本研究で は有限要素法を用いたカルマンフィルタ理論を適用し,モ デルの不確実性を確率過程に基づいて修正し,推定を行 う. カルマンフィルタ理論とは,確率過程に基づいた解析 手法の一つであり誤差を含む観測値から,その誤差を考慮 しながら推定,同定,制御,予測,データ同化等様々な分野 で用いられるフィルタリング理論の一つである.そして, このカルマンフィルタ理論と有限要素法を組み合わせる ことによって,時系列のみならず空間方向にも解析が可能 となる.
本研究では,数値解析例として,日本でも有数の汚染河 川である千葉県手賀沼を実在のモデルとして扱った. そ してこれまでの研究では,河川の汚染を表す指標である
COD
の推定結果について満足な結果を得られていなかっ たため，推定に非常に影響をもつ拡散係数に着目し,拡 散係数を同定し,同時に推定を行った. 同定をするために フィルタリング理論である拡張カルマンフィルタと有限 要素法を組み合わせた拡張カルマンフィルタ有限要素法 を用いる. また,多くの解析時間及び計算容量がかかって しまうという拡張カルマンフィルタ有限要素法の問題の 削減を目的とした, 縮小カルマンフィルタ有限要素法を 提案する. 新手法である縮小カルマンフィルタ有限要素 法を河川に適用し,有効性を検証する. また有限要素法の 基礎方程式として河川の流れを再現するために浅水長波 方程式，CODの流れを再現するために移流拡散方程式を 用いる.今回は同定手法の有効性を示すために異なった境界条 件での計算と、領域分割した場合と一つの領域として計 算した時の計算時間や計算の正確性についてのいくつか 場合に分け計算を行った.
2
数値解析手法2.1
状態方程式今回は手賀沼川という水深の非常に浅く，COD濃度の
高い河川をモデルとして計算するため，河川の水流を表 現するために浅水長波方程式を用い
COD
の移流拡散を 表現するために移流拡散方程式を状態方程式として用い る.＜浅水長波方程式＞
˙
u i + u j u i,j + g(η + H + h), i −ν(u i,j + u j,i ) ,j = 0 (1)
˙
η + {(η + H)u i } ,i = 0, (2)
＜移流拡散方程式＞
˙
c + u i c, i −κc, ii = 0 (3)
ここで,ui
はx,y
方向の流速,gは重力加速度,ηは水位変 動量,H は基準水深,hは河床勾配,νは渦動粘性係数,cはCOD,κ
は拡散係数を表している.2.2
カルマンフィルタカルマンフィルタとは, 確率過程に基づいた推定手法 の一つである. 本手法を用いることによりモデルの不確 実性を確率過程に基づいて修正し, 状態量を逐次決定的 に推定することができる. 同定問題は非線形問題となる ため,カルマンフィルタを拡張する必要がある. 拡張カル マンフィルタの状態空間モデル,システム方程式
(4)
と観 測方程式(5)
をそれぞれ以下に示す.{x n+1 } = [F n ]{x n } + [G]{w n } (4) {y} n = [H]{x} n + [G]{v n } (5)
ここでx n
は状態量,F n
はx n
の関数で表される状態遷 移行列,G n
は駆動行列,y n
は観測値,H
は観測行列であ る.w n
及びv n
はシステムノイズ及び観測ノイズと呼ば れ,平均0,
共分散行列に従う正規白色ノイズである.以下にカルマンフィルタのアルゴリズムを示す.
1. [ˆ Γ
0|−1] = [Γ
x0] , {ˆ x
0|−1} = {¯ x
0} 2. [K
n] = [ˆ Γ
n][H ]
T([R] + [H][ˆ Γ
n][H ]
T)
−13. [ ˆ P
n] = ([I] − [K − n][H ])[ˆ Γ
n]
4. [ˆ Γ
n] = [F
n][ ˆ P
n][F
n]
T+ [G
n][Q][G
n]
T5. {ˆ x
n|n−1} = [F
n]{ˆ x
x−1/x} + {B}
6. {ˆ x
n|n} = {ˆ x
n|n−1} + [K
n]({y
n} − [H
n]{ˆ x
n|n−1})
ここで
K n
はカルマンゲイン, ˆP n
は推定誤差共分散行 列, ˆΓ n
は予測誤差共分散行列,Q, R
は共分散行列である.本研究ではカルマンフィルタのアルゴリズムに有限要素 法で近似した有限要素方程式を物理モデルとして組み込 む. そして,浅水長波方程式と移流拡散方程式を
Galerkin
法と陽的Euler
法によって離散化した有限要素近似を用いることにより河川の流速,水位変動量及び
COD
を推定 する. ここで,空間方向及び時間方向に離散化した有限要 素方程式は以下の式で与えられる.＜浅水長波方程式＞
1
∆t M ¯ αiβj u n+1 = 1
∆t M ˜ αiβj u n βj − (K αβγj U βj u γi (6) +K αβγj u βj U γi + D αiβj u βj + H αiλ η λ + ˆ Ω αi )
1
∆t M ¯ λµ η n+1 µ = 1
∆t M ˜ λµ η µ n − (A λµβj u βj E µ (7) +A λµβj U βj η µ + B λβj u βj = ˆ Σ λ )
＜移流拡散方程式＞
1
∆t M αβ c n+1 β = 1
∆t M αβ c n β − K αβγj c γ (U βj − u βj ) (8) +D αβ c β − Ψ ˆ α
ここで
M αiβj = M αβ δ ij K αβγj = Z
V
(Φ α Φ β Φ γ,j )dV
M αβ = Z
V
(Φ α Φ β )dV H αiλ = g Z
V
(Φ α,i Φ λ )dV
D αiβj = ν Z
V
(Φ α,j Φ β,i )dV + ν Z
V
(Φ α,k Φ β,k )δ ij dV
A λµβj = Z
V
(Φ λ Φ µ Φ β,j )dV + Z
V
(Φ λ Φ β Φ µ,j )dV
B λβj = A λµβj (H µ + h µ )dV
D αβ = κ Z
V
(Φ α,k Φ β,k )dV Σ ˆ λ = ˆ Σ u λ
Ψ ˆ α = Z
S
(Φ α ˆ b)dS + ˆ Ψ u α Ω αi = Z
S
(Φ α s ˆ i )dS + ˆ Ω u αi
まず有限要素法で近似した有限要素方程式を物理モデ ルとしてカルマンフィルタに組み込むことでカルマンフィ ルタ有限要素法を導く. カルマンフィルタに物理モデル を組み込む方法は大別して2
種類ある. 一つは状態方程 式に組み込む方法,もう一つは観測方程式に組み込む方法 である. 前者の方は,時系列過程を追従することができる ので, カルマンフィルタのアルゴリズムをそのまま用い ることができる.本手法をカルマンフィルタ有限要素法と呼ぶ. 本手法を 用いることにより, 状態量の時間的分布のみにならず空 間的分布においても推定が可能となる. また,本研究では 共分散行列に
4
次スプライン関数,S = ¯ S {1 . 0 − 6 . 0
„ r
id
m«
2+ 8 . 0
„ r
id
m«
3− 3 . 0
„ r
id
m«
4} (9)
dm ri
main point related point
図
1;
節点及び節点周りの領域0 0.2 0.4 0.6 0.8 1
0 0.2 0.4 0.6 0.8 1
Quartic Spline Function
Quartic Spline Function
図
2; 4
次スプライン関数を提案する. ここで,Sは共分散の空間分布を,
r i , d m
は節点間距離,最大節点間距離を表している. また,観測 誤差
v k ,
システム誤差w k
は観測値と推定値の差,最適推 定値と推定値の差によってそれぞれ求めることができる.共分散の空間分布を節点周りに定義することによってカ ルマンフィルタ有限要素法が持つ解析時間及び解析容量 が増加する問題点を解決することが出来る.そこで本研究 では観測点付近の情報のみでカルマンフィルタ理論を適 用する,観測点周りの領域とそれ以外の領域の領域に分け 観測点に近いところ程重みを大きくして計算領域を縮小 する．(4)式は次のように変形される．
x
n+1z
n+1ff
=
» F
nL
nM
nN
n– x
nz
nff +
U
nV
nff +
Gw
n0
ff (10)
{y} n = [H ]{x} n + [G]{v n } (11)
ここでx k
は状態量,F n {x n }
は状態遷移行列,G
は駆 動行列,y n
は観測値,H
は観測行列である.U n ,V n
及びv n
はシステムノイズ及び観測ノイズと呼ばれ,平均0,
共 分散行列Q, R
に従う正規白色ノイズである.このように,計算領域を 観測点付近とそれ以外に分離 させるとシステム方程式を四分割する事ができ，これに よりカルマンゲイン
K
を求めるために用いる行列がF n
のみとなり計算要領と時間の削減を行うことが可能とな る.以下に計算時間の短縮を行う時のアルゴリズムを示す．
1. [ˆ Γ
0/−1] = [Γ
x0] , {ˆ x
0/−1} = {¯ x
0} 2. [K
k] = [ˆ Γ
n][H ]
T([R] + [H][Γ
n][H ]
T)
−13. [P
n= ([I] − [K
n][H])[Γ
n]
4. [Γ
n+1] = [F
n][P
n][F
n]
T+ [G
n][Q][G
n]
T5. {C
n} = [L
n]{[M
n−1]{x}
n−1+ [N
n−1]{z}
n−1+ {V }
n−1} + {U
n} 6. {ˆ x
n/n−1} = [F
n]{ˆ x
n−1/n−1} + {C}
7. {ˆ x
n/n} = {ˆ x
n/n−1} + [K
k]({y
n} − [H ]{ˆ x
n/n−1}
ここで
x ˆ k/k
は観測データy 0 · · · y t
が与えられたときのx k/k
の最小分散推定値である.2.3
同時推定同定手法本研究では
COD
に非常に関係のある拡散係数のパラ メータ同定を行う事で最適な拡散係数を導出し, 得られ た拡散係数の値で推定を同時に行う新しい手法を提案し た. 同時同定手法として,拡散係数を深さの異なる場所に おいて異なり，また時間により変化すると仮定して定数 である拡散係数をカルマンフィルタを用いて逐次的に求 めていく,まず拡散係数を以下のように展開する.κ n+1 m = κ n m + κ n m,β (c k β − c n+1 β ) + · · · (12)
この時，mは分割された領域の番号であり，以下がその イメージ図である．κ
3κ
4κ
1κ
2κ
3κ
4κ
1κ
2図
3;
領域分割1
∆t M ¯ αβ c n+1 β = 1
∆t M ˜ αβ c n β − K αβγj c γ u βj (13) +κ λ D αβ λ c β + ˆ Ψ α
次に, (12)式
(13)
式より以下の式を得ることができる.» M ¯
αβ0 0 E
mn– » c
βκ
m–
n+1=
» H
αβG
αnA
βnE
mn– » c
βκ
n–
n(14)
−
" Φ ˆ
α− A
λµcn−1 β
#
nここで,H
αβ
は状態遷移行列,Gαn
は拡散行列,Aβn
は感 度行列,Emn
は単位行列, ˆΦ α
は既知項である. 本研究で は以上の式を用いて,推定と同定を同時に行った.3
数値解析例3.1
解析モデル数値解析例として,カルマンフィルタ有限要素法を用い て千葉県手賀沼における流況を推定する. 観測値は
2003
年8
月1
日から5
日における流速u, v
と水位変動量η
及 びCOD
が4
点の観測地点において利根川下流河川事務 局によって得られている. 推定には観測地点No.1, No.2, No4
の3
点を用い, No3での観測値と推定値の比較を行 い,本手法の有効性を示した.解析領域と観測地点を図4
に示す. 解析領域は観測点No1
から手賀沼右端までとし，有限要素分割図は総節点数
4044，総要素数 7552
の分割 として図5
に示す．また全ての解析において，時間増分 量を∆t
を0.01，重力加速度 g
を9.8(m/s 2 )
としている.図
4;
解析領域図
5;
有限要素分割図流入条件に観測地点
1
番で得られた観測データを与え ました. 図6〜9
に観測地点1
番における流速u, v
と水 位変動量η
及びCOD
の観測値をそれぞれ示す. ,観測値 はいずれも利根川下流河川事務所より提供していただい たものである．-5 0 5 10 15 20 25
0 20 40 60 80 100 120
x-velocity
TIME(hour) observe
図
6; x
方向の流速-10 -5 0 5 10 15
0 20 40 60 80 100 120
y-velocity
TIME(hour) observe
図
7; y
方向の流速-0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25
0 20 40 60 80 100 120
Water elevation(m)
TIME(hour) observe
図
8;
水位変動量4 5 6 7 8 9 10 11
0 20 40 60 80 100 120
COD
TIME(hour) observe
図
9; COD
また,前日の観測値を用いて定常になるまで順解析を 行い，その値を初期条件とする．図
10
に解析で用いた流 速，水位変動量及びCOD
の初期条件を示す.0.1m/s
H
-0.01 -0.0095 -0.009 -0.0085 -0.008 -0.0075 -0.007 -0.0065 -0.006 -0.0055 -0.005 -0.0045 -0.004 -0.0035 -0.003 -0.0025 -0.002
C
8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3 2.5 2 1.5 8.5 9
図
10;
初期条件3.2
解析1(推定)
解析
1
では,本来一定でない定数である,拡散係数κ
を0.2
としてカルマンフィルタ有限要素法を用た河川の流況 の推定を行った. 観測点1
番,2番及び4
番の観測値を用 いて,観測値3
番の流況を自作のプログラムにより推定し た. その結果として,図11〜14
に推定点(観測地点 3
番) における流速u, v
と水位変動量η
及びCOD
の推定値の 時刻暦と観測値の時刻暦との比較をそれぞれ示す. また 表1
より, 観測値とほぼ同等の推定値が得られているこ とがわかる.-8 -6 -4 -2 0 2 4 6 8
0 20 40 60 80 100 120
x-velocity
TIME(hour) observe estimate
図
11; x
方向の流速-12 -10 -8 -6 -4 -2 0 2 4
0 20 40 60 80 100 120
y-velocity
TIME(hour) observe estimate
図
12; y
方向の流速-0.1 -0.05 0 0.05 0.1 0.15
0 20 40 60 80 100 120
Water elevation
TIME(hour) observe estimate
図
13;
水位変動量6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 11.5 12
0 20 40 60 80 100 120
COD
TIME(hour) observe estimate
図
14; COD
相関係数
x
方向の流速0.71 y
方向の流速0.69
水位変動量0.71
COD 0.63
表
1;
相関係数3.3
解析2(推定及び同定)
解析
1
では,本来一定でない値であるはずの拡散係数κ
を一定として解析を行ってしまっている. この問題を解 決するために解析2
として, カルマンフィルタ有限要素 法の同定手法を用いることで拡散係数の同定を行い, 同 定を行うと同時に推定を行うことでより実際の現象に近 いCOD
の推定を行った. 拡散係数の同定結果を図15
に, 同定された拡散係数を用いたCOD
の推定結果を図16
に それぞれ示す. また表2
より,解析1
より良い精度の解析 を行えた.0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.21
0 20 40 60 80 100 120
k
TIME(hour) Diffusion coefficient
図
15;
拡散係数6.5 7 7.5 8 8.5 9 9.5 10 10.5 11 11.5 12
0 20 40 60 80 100 120
COD
TIME(hour) observe estimate
図
16; COD
相関係数
COD 0.73
表
2;
相関係数4
まとめ本研究では河川における,流況の推定及び同定手法とし て拡張カルマンフィルタ有限要素法を提案し,その妥当性 の検証を観測データとの比較によって確認した. 誤差を 含む観測データを用いて推定と同定を同時に行い, 拡散 係数を変数として扱うことでより実際の現象に近い解析 を行った. また,今までカルマンフィルタ有限要素法を扱 う上で一番の問題となっていた解析時間と解析容量の大 幅な削減を縮小カルマンフィルタ有限要素法を用いるこ とで可能とした. 今後の課題として,より大規模な実在モ デルへの適用と, 未来を推定する予測問題へ発展させる ことが挙げられる.
参考文献
