Identification of Dynamic Motion of the Ground using the Kalman Filter Finite Element Method

(1)

修士論文発表要旨

(2009

年度)

カルマンフィルタ有限要素法を用いた地盤動的運動の同定に関する研究

Identification of Dynamic Motion of the Ground using the Kalman Filter Finite Element Method

土木工学専攻３９号山本智史

Satoshi YAMAMOTO

1

はじめに

ダムを建設するにあたり,材料費用削減の為に近隣の地山を切り土し建設材料の一部として利用することがある. 地山の作業現場においての安全性を確保する為に,事前に軟弱地盤の特定や地山全体の状況を把握することは大変重要あり,その把握方法の一つとして,発破によって発生した振動を加速度計,速度計を用いて観測する方法がある. しかしながら,この方法では地山全体の状態を把握する為に,多数の観測機械を設置する必要があり,膨大な費用,時間が必要となる問題や,最も知りたい点での観測が困難であることや,観測する際に機械,人為的誤差等の誤差が含まれてしまうという問題がある. そこで,本研究では有限要素法を用いたカルマンフィルタ理論を適用し,モデルの不確実性を確率過程に基づいて修正し,推定,及び同定を行う. カルマンフィルタ理論とは,確率過程に基づいた解析手法の一つであり誤差を含む観測値から,その誤差を考慮しながら推定,同定,予測,データ同化等様々な分野で用いられるフィルタリング理論の一つである. 更に,カルマンフィルタ理論と有限要素法を組み合わせることで,時系列のみならず空間的に解析することが可能になる.

本研究では,数値解析例として,宮城県箕輪山を実在のモデルとして扱い,発破時の振動による速度を観測データとして加速度の推定を行う. ここで,実際の発破外力の大きさを決定することは非常に困難である. それゆえに,加速度の推定と同時に発破外力の同定を行う. また,多くの解析時間及び,計算容量がかかってしまうというカルマンフィルタ有限要素法の問題の解決を目的とした, 縮小カルマンフィルタ有限要素法を提案する. この縮小カルマンフィルタを実際の箕輪山に適用し,有効性を検証する.

2

数値解析手法

2.1

カルマンフィルタ

カルマンフィルタとは, 確率過程に基づいた推定手法の一つである. 本手法を用いることによりモデルの不確実性を確率過程に基づいて修正し, 状態量を逐次決定的に推定することができる. カルマンフィルタの状態空間モデルを以下に示す.

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

ここで,

x k

はシステムの状態量,

f k

は状態遷移行列,

G k

は駆動行列,

y k

は観測データ,

h k

は観測行列である.

v k

及び

w k

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

0 ,

共分散行列

Q w k , R v k

に従う正規白色ノイズであ

る. 以下にカルマンフィルタのアルゴリズムを示す.

1. [ ˆ P 0|−1 ] = [P x 0 ] , {ˆ x 0|−1 } = {¯ x 0 } 2. {ˆ x k+1|k } = [F k ]{ˆ x k|k }

3. [K k ] = [ ˆ P k|k−1 ][H k ] ^T ([R _vk ] + [H k ][ ˆ P k|k−1 ][H k ] ^T ) ⁻¹ 4. [ ˆ P k|k ] = ([I] − [K k ][H k ])[ ˆ P k|k−1 ]

5. {ˆ x k|k } = {ˆ x k|k−1 } + [K k ]({y k } − [H k ]{ˆ x k|k−1 }) 6. [ ˆ P k+1|k ] = [F k ][ ˆ P k|k ][F k ] ^T + [G k ][Q _wk ][G k ] ^T

ここで, ˆ

x k|k

は観測データ

y 0 ,

…,

y k

が与えられたときの

x k

の最小分散推定値である.

K k

はカルマンゲインである.

2.2

状態方程式

トンネルの周辺の地山の状態を表す状態方程式として弾性体の方程式を用いる. 応力の釣り合い方程式,変位とひずみ関係式, 応力とひずみの関係を表す土の構成則に対して四面体要素を用い,

Galerkin

法を適用することにより,有限要素近似を行う. ここで,空間方向に離散化した有限要素方程式は以下の式で与えられる.

M iαkβ u ¨ kβ + C iαkβ u ˙ kβ + K iαkβ u kβ = ˆ γ αi (3)

また, 時間方向の離散化には

N ewmarkβ

法を用いた.

N ewmarkβ

法は

u ⁽ⁿ⁺¹⁾ _iβ ,¨ u ⁽ⁿ⁺¹⁾ _iβ

を以下のように仮定する.

u ⁿ⁺¹ _βk = u ⁿ _βk + Δ t

2 ( ˙ u ⁿ⁺¹ _βk + ˙ u ⁿ _βk ) , (4) u ¨ ⁿ⁺¹ _βk = −¨ u ⁿ _βk + 2

Δ t ( ˙ u ⁿ⁺¹ _βk − u ˙ ⁿ _βk ) , (5)

上式で求めた

u ⁽ⁿ⁺¹⁾ _iβ , ˙ u ⁽ⁿ⁺¹⁾ _iβ

を以下の式に代入することで加速度を求めることができる. 以下に加速度

u ¨ βk ,

速度

u ˙ βk ,

変位

u βk

に関する有限要素方程式を示す.

u ˙ ⁿ⁺¹ _βk = F αiβk u ˙ ⁿ _βk − L ⁻¹ _αiβk ( M αiβk u ˙ ⁿ + K αiβk u ⁿ + γ αi ) , (6)

ここで,

γ αi

は外力項であり. 発破によって発生するボアホール圧力は一気に上昇し,また時間を掛けて伝播していくため,外力をそのように与える必要があり,その圧力の履歴を追えるように与える必要がある. 以下のボアホール圧を表す式で定義される.

ˆ γ αi = M m=1

A ^(m) _αi ( e ^−ηt − e ^−ξt ) (7)

上式で求めた

A ^(m) _αi

は,外力の最大振幅を表し,本研究では

x, y, z

方向の外力の最大振幅である

A ^(m) _αi

をそれぞれ同時に同定する.

(2)

2.3

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

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

2

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

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

4

次スプライン関数,

S = ¯ S{1.0 − 6.0

„ r i

d m

« ₂ + 8.0

„ r i

d m

« ₃

− 3.0

„ r i

d m

« ₄ } (8)

を提案する. ここで,

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 (9)

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

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.4

同時推定同定手法

推定を行うにあたり,振動の起因となる発破外力をどのように定めるがということは, 非常に重要な問題である.

それゆえ,振動の起因となる発破外力の最大振幅

A

の同定を行う事で最適な発破外力を導出し, 得られた発破外力の値で推定を行う新しい手法を提案した. 同時推定同定手法として,以下のように式を展開する.

u ˙ ⁿ⁺¹ _αi γ _βj ⁿ⁺¹

=

F αiβj −L ⁻¹ _αiβj k

u ˙ ⁿ _αi γ _βj ⁿ

+ B αi (10)

k = γ βj n+1

γ _βj ⁿ = e ^ξ(t ⁿ ^+Δt) − e ^η(t ⁿ ^+Δt)

e ^ξ(t ⁿ ⁾ − e ^η(t ⁿ ⁾ (11)

ここで,

F αβ

は,カルマンフィルタ有限要素法を用いて推定のみ行う場合に適用される状態遷移行列であり,

B αi

は既知項である. そして,発破外力同定を行う上で,重要な

k

は,外力が発生している微小間内におけるある時点の外力項とその次の時点における外力項の比を表している. 本研究では,

F αβ ,−L ⁻¹ _αiβj , k

を,新しい状態遷移行列

F n

として以下の縮小カルマンフィルタアルゴリズムへ適用することによって,推定と同定を同時に行った.

2.5

縮小カルマンフィルタアルゴリズム

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

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

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 ) ⁻¹ 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 } = {B 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)

における発破外力によって引き起こされた状態量の推定及び発破外力の同定を同時に行う. そして,この手法において発破外力の最大振幅を表す

A

の初期値を決定する必要があるため,2つのケースに分けて解析を行い本手法の有効性を検討する. 図

4

に有限要素分割図を示す. 総接点数は

5603,

総要素数は

26102

である. 本研究では

2005

年

9

月

20

日において観測された速度

u ˙ i

を観測値として用いた. 表

1

は,各ケースにおける,

x, y, z

方向それぞれの

A

の初期値を表す.

(3)

図

3;

箕輪山山頂における二つ石現場

図

4;

有限要素分割図

図

5;

観測点配置図

本研究では,Nos.1,2,4の速度データを観測値を使い,No.2 における加速度を推定する. 外力条件として図

8

における白い楕円エリアが外力を与えたエリアであり,

x, y, z

方向にそれぞれ外力を加える. 以下に各方向へ加えた外力のイメージ図を示す.

図

6 x,y

方向の外力

図

7 z

方向の外力

図

8〜10

に観測点

No.3

での

x,y,z

方向の速度観測デー

タを示す. その他の点についての観測データは本稿では割愛する.

図

8; x

方向の速度

図

9; y

方向の速度

図

10; z

方向の速度

表

1;

各ケースにおける

x, y, z

方向の初期値

Case X

方向初期値

Y

方向初期値

Z

方向初期値解析

1 2 . 925

×

10 ⁹ 2 . 925

×

10 ⁹ 3 . 696

×

10 ⁸

解析

2 2 . 925

×

10 ⁸ 2 . 925

×

10 ⁸ 3 . 696

×

10 ⁷

(4)

3.2

解析

1

結果

図

11

に解析

1

における同定された発破外力の時刻暦, 図

12〜図 14

に観測点

No.2

における同時推定同定手法より推定された各方向の加速度を表す. ここで,解析

1

における各方向初期値は推定のみの検証を行いある程度定められた値である. 推定結果に着目すると,各方向共に観測値と推定値が高い相関を持っていることがわかる. しかしながら,このケースでは初期値

A

を,推定のみの検証で精度の良い結果が得られた

A

として設定しているため良い結果が得られたと考えられる. それゆえ,初期値を変えたケースの解析結果を解析

2

において検証する.

図

11;

同定発破外力図

12; x

方向加速度

図

13; y

方向加速度

図

14; z

方向加速度

3.2

解析

2

結果

図

15

に解析

1

における同定された発破外力の時刻暦, 図

16〜図 18

に同時推定同定手法より推定された各方向の加速度を表す. ここで,解析

2

における

A

の初期値は解析

1

と比べて

1

オーダー下げた値である. 解析

1

と

2

の同定結果を比べると,各方向において解析

1

の同定値よりも小さな値を示している. 次に推定結果に着目すると,Z方向におけるピークでの現象を追えていないが,各方向ともに挙動を追えていることがわかる.

図

15;

同定発破外力図

16; x

方向の加速度

図

17; y

方向加速度図

18; z

方向加速度

5

まとめ

本研究では箕輪山における, 地山の発破外力の同定及び,加速度の推定手法として縮小カルマンフィルタ有限要素法を提案し, その妥当性を実際の観測データとの比較によって確認した. 本研究では

3

つの限られた観測点から推定点における加速度を推定できたことから, 限られた観測値からの推定に成功したと言えよう. 本研究の最大の目的である,誤差を含む観測データを用いて推定と同定を同時に行い,同定対象

A

の初期値が推定精度にどの程度影響するのかを

2

つのケースに分けて検証し,同時同定推定に成功した. また,今までカルマンフィルタ有限要素法を扱う上で,最大の問題となっていた解析時間と解析容量の大幅な削減を縮小カルマンフィルタ有限要素法を用いることで可能とした. その結果,以前から課題であった

3

次元解析が可能となった.

しかしながら,ケース

2

の結果に着目すると,

x, y

方向においては精度の良い推定に成功したが, z方向の結果に着目するとピーク値が約

40

％ずれてしまっている. その為,本手法が各方向初期値

A

に依存してしまうという問題点が大きな検討課題であるといえる..

参考文献

1 R.E.Kalman, A New Approach to Linear Filtering and Prediction Problems, Trans.ASME,J. Basic Eng., vol182D, no.1,34-45,(1960)

2 A.Murakami and T.Hasegawa, Back Analysis by Kalman Filter Finite Elements and a Determina- tion of Optimal Observed Points Location. Jour- nals of the Japan Society of Civil Engineers Vol.388,227- 235,(1987)

3

船越雄二,’拡張カルマンフィルタ有限要素法を用いた物質濃度の推定’,中央大学修士論文,(2001)

4

樋川敦,’逆解析によるトンネル切羽前方のヤング係

数の同定’,中央大学修士論文,(2004)

5

加藤祐介,’拡張カルマンフィルタ有限要素法を用いた弾性係数の同定’,中央大学修士論文,(2006)

6

尾島康則,’カルマンフィルタ有限要素法を用いた河

川の解析’,中央大学修士論文,(2008)

Identification of Dynamic Motion of the Ground using the Kalman Filter Finite Element Method

(2009

カルマンフィルタ有限要素法を用いた地盤動的運動の同定に関する研究

Identification of Dynamic Motion of the Ground using the Kalman Filter Finite Element Method

Satoshi YAMAMOTO

1

2

2.1

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

x k

f k

G k

y k

h k

v k

w k

0 ,

Q w k , R v k

1. [ ˆ P 0|−1 ] = [P x 0 ] , {ˆ x 0|−1 } = {¯ x 0 } 2. {ˆ x k+1|k } = [F k ]{ˆ x k|k }

3. [K k ] = [ ˆ P k|k−1 ][H k ] T ([R vk ] + [H k ][ ˆ P k|k−1 ][H k ] T ) −1 4. [ ˆ P k|k ] = ([I] − [K k ][H k ])[ ˆ P k|k−1 ]

5. {ˆ x k|k } = {ˆ x k|k−1 } + [K k ]({y k } − [H k ]{ˆ x k|k−1 }) 6. [ ˆ P k+1|k ] = [F k ][ ˆ P k|k ][F k ] T + [G k ][Q wk ][G k ] T

x k|k

y 0 ,

y k

x k

K k

2.2

Galerkin

M iαkβ u ¨ kβ + C iαkβ u ˙ kβ + K iαkβ u kβ = ˆ γ αi (3)

N ewmarkβ

N ewmarkβ

u (n+1) iβ ,¨ u (n+1) iβ

u n+1 βk = u n βk + Δ t

2 ( ˙ u n+1 βk + ˙ u n βk ) , (4) u ¨ n+1 βk = −¨ u n βk + 2

Δ t ( ˙ u n+1 βk − u ˙ n βk ) , (5)

u (n+1) iβ , ˙ u (n+1) iβ

u ¨ βk ,

u ˙ βk ,

u βk

u ˙ n+1 βk = F αiβk u ˙ n βk − L −1 αiβk ( M αiβk u ˙ n + K αiβk u n + γ αi ) , (6)

γ αi

ˆ γ αi = M m=1

A (m) αi ( e −ηt − e −ξt ) (7)

A (m) αi

x, y, z

A (m) αi

2.3

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 } (8)

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

(9)

K

F k

1;

2; 4

3. [K k ] = [ ˆ P k|k−1 ][H k ] ^T ([R _vk ] + [H k ][ ˆ P k|k−1 ][H k ] ^T ) ⁻¹ 4. [ ˆ P k|k ] = ([I] − [K k ][H k ])[ ˆ P k|k−1 ]

5. {ˆ x k|k } = {ˆ x k|k−1 } + [K k ]({y k } − [H k ]{ˆ x k|k−1 }) 6. [ ˆ P k+1|k ] = [F k ][ ˆ P k|k ][F k ] ^T + [G k ][Q _wk ][G k ] ^T

u ⁽ⁿ⁺¹⁾ _iβ ,¨ u ⁽ⁿ⁺¹⁾ _iβ

u ⁿ⁺¹ _βk = u ⁿ _βk + Δ t

2 ( ˙ u ⁿ⁺¹ _βk + ˙ u ⁿ _βk ) , (4) u ¨ ⁿ⁺¹ _βk = −¨ u ⁿ _βk + 2

Δ t ( ˙ u ⁿ⁺¹ _βk − u ˙ ⁿ _βk ) , (5)

u ⁽ⁿ⁺¹⁾ _iβ , ˙ u ⁽ⁿ⁺¹⁾ _iβ

u ˙ ⁿ⁺¹ _βk = F αiβk u ˙ ⁿ _βk − L ⁻¹ _αiβk ( M αiβk u ˙ ⁿ + K αiβk u ⁿ + γ αi ) , (6)

A ^(m) _αi ( e ^−ηt − e ^−ξt ) (7)

A ^(m) _αi

A ^(m) _αi

« ₂ + 8.0

« ₃

« ₄ } (8)

u ˙ ⁿ⁺¹ _αi γ _βj ⁿ⁺¹

F αiβj −L ⁻¹ _αiβj k

u ˙ ⁿ _αi γ _βj ⁿ

γ _βj ⁿ = e ^ξ(t ⁿ ^+Δt) − e ^η(t ⁿ ^+Δt)

e ^ξ(t ⁿ ⁾ − e ^η(t ⁿ ⁾ (11)

F αβ ,−L ⁻¹ _αiβj , k

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 }

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 })

10 ⁹ 2 . 925

10 ⁹ 3 . 696

10 ⁸

10 ⁸ 2 . 925

10 ⁸ 3 . 696

10 ⁷