Parameter Identiﬁcation of Elastic Modulus of Rock Based on Blasting Waves in Tunnel Excavation

(1)

修士論文発表要旨

(2011

年度

)

トンネル掘削における発破振動を用いた弾性係数の同定に関する研究

Parameter Identiﬁcation of Elastic Modulus of Rock Based on Blasting Waves in Tunnel Excavation

土木工学専攻

34

号本山雄斗

Yuto MOTOYAMA

1 _はじめに

トンネルやダム，道路などの土木工事を行う上で地盤の性質や状態を予め調べることは重要である．地盤に外力が作用した場合，それがどのような振舞いをみせるか事前に把握しておくことは，工事を進める上での安全管理や費用対策，環境対策に直接的に関わってくる．今日の工事設計にあたり，ボーリング調査や物理探査，岩石試験等の地質調査方法が地山性状を調査する一般的なものとして採用されている．しかし，これら従来の方法は調査に費やす時間やコストを考えると，改善すべき点が多々挙げられる．そこで本研究では，トンネル掘削現場を対象として，数値解析を用いた新しい地質調査方法を提案する．数値解析を用いた地質調査は，安価で安全に行え，かつ汎用性が高いという点で有益である．トンネル掘削に限定して言えば，軟弱地盤層や断層破砕帯への進入は人的被害や工期の遅延を招く恐れがあるため回避しなければならない．数値解析により切羽前方の地山性状を把握できる本手法は，費用や安全面を考慮した革新的な調査方法と言えるだろう．

解析対象には，広島県に位置する尾道・松江線高野工区大万木トンネル工事現場

(以下，大万木トンネル)

を採用した．大万木トンネルは総延長

4, 878m

で，NATM工法を用いた山岳トンネルである．また，事前の調査で地山が

3

層に分類されていることがわかっている．従って，

パラメータ同定により地山を構成する各層の弾性係数を同定することを本研究の目的とする．

2 パラメータ同定

2.1 一次随伴法

添字記法と総和規約を用いて，パラメータ同定手法を示す．本稿で述べているパラメータとは弾性係数のことである. 弾性係数を同定することは，以下の拡張評価関数

J

^∗を最小にするような最小化問題と等価である．

J

^∗

= J + Λ + Ξ, (1)

式

(1)

における

J

は評価関数と呼ばれ，以下の形で定義される．

J = 1 2

∫ t

_f

t

0

( ˙ u αi − η αi )Q αiβj ( ˙ u βj − η βj )dt, (2)

ここに，˙

u αi

と

η αi

はそれぞれ，観測点

α

における

i

方向の計算速度と観測速度である．Q

αiβj

は観測点数の自由度を持った重み係数で，t

0

を

t f

はそれぞれ，初期時間と終端時間を意味する．支配方程式には式

(3)

に示す，変位で表したナビエの式を適用する．本研究では地盤を線形弾性体と仮定し，四面体要素を用いて動的解析を行う．

D _ijkl u _k,lj + ρb _i − ρ¨ u _i = 0, (3)

ここに，u

i

，b

i

，ρはそれぞれ，変位，物体力そして密度を意味し，従属変数の上部ドット記号は時間微分を表す．

D ijkl

は弾性係数行列と呼ばれ，地盤を構成する層の数だけ与えられる

(図 1

参照)．それぞれの層の弾性係数とポアソン比に基づいて，D

ijkl

は以下のように導くことができる．

D _ijkl =

∑ M

m=1

D ^(m) _ijkl , (4)

ここに，mは地山を構成する層の数を意味し，Mはその最大数である．式

(5)

は支配方程式を有限要素法により離散化し，Lagrange乗数法を適用したものである．

Λ =

∫ t

f

t

0

λ αi ( ˆ Ω αi − M αiβj u ¨ βj

− C _αiβj u ˙ _βj − K _αiβj u _βj )dt, (5)

式

(5)

における

λ αi

は

Lagrange

乗数を表し，係数行列

M _αiβj

，C

αiβj

，K

αiβj

はそれぞれ，質量，レイリー減衰，

剛性を意味する．

Ω ˆ αi

は外力項である. 式

(1)

の右辺第３項の

Ξ

は安定化項で，以下の形で定義される．

Ξ = 1 2

∫ t

f

t

₀

(X _λ − X ˜ _λ )W _λµ (X _µ − X ˜ _µ )dt, (6) X _λ

は本研究の同定対象である弾性係数を意味している．

また，W

αiβj

は計算の安定性を確保するための安定化重みである．式

(6)

は繰返し計算収束時にはゼロとなる．

(2)

図

1:

地盤の構成

2.2 二次随伴法

逆問題の定式化において，本研究では随伴変数法を用いる．また，Lagrange乗数法により拡張された拡張評価関数

J

^∗は以下に示す．

J

^∗

= J +

∫ t

_f

t

0

λ αi (M αiβj u ¨ βj + C αiβj u ˙ βj

+ K αiβj u βj − Ω ˆ αi )dt, (7)

各変数のテイラー展開

u _i = u ⁽⁰⁾ _i + u ⁽¹⁾ _i + u ⁽²⁾ _i + u ⁽³⁾ _i + ..., (8) λ _αi = λ ⁽⁰⁾ _αi + λ ⁽¹⁾ _αi + λ ⁽²⁾ _αi + λ ⁽³⁾ _αi + ..., (9)

各変数と同様に弾性係数

X _i

も以下のようにテイラー展開を行う．

X i = X _i ⁽⁰⁾ + X _i ⁽¹⁾ + X _i ⁽²⁾ + ..., (10)

展開された各変数を拡張評価関数

J

^∗に代入すると，J^∗ は以下のように書き表すことができる．

J

^∗

= J

^∗

⁽⁰⁾ + J

^∗

⁽¹⁾ + J

^∗

⁽²⁾ + J

^∗

⁽³⁾ + ..., (11)

ここで，J^∗は以下のように書き表すことができる．

J

^∗

= J

^∗

⁽⁰⁾ +

∫ t

f

t

0

(g i X i + 1

2 X i H ij X j + ..., )dt (12)

ここで，g

i

は弾性係数

X i

の勾配であり，

H ij

はヘシアンマトリックスである．また，拡張評価関数は以下のように各項を分類することができる．

J

^∗

= J ₁

^∗

⁽⁰⁾ + J ₁

^∗

⁽¹⁾ + J ₁

^∗

⁽²⁾

+ J ₂

^∗

⁽²⁾ + J ₁

^∗

⁽³⁾ + J ₂

^∗

⁽³⁾ + ... (13)

J ₁

^∗

⁽¹⁾ =

∫ t

_f

t

0

(g i X _i ⁽¹⁾ dt (14)

J ₁

^∗

⁽²⁾ =

∫ t

_f

t

0

(g i X _i ⁽²⁾ dt (15)

J ₂

^∗

⁽²⁾ = 1 2

∫ t

f

t

0

(X _i ⁽¹⁾ H ij X _j ⁽¹⁾ )dt (16)

J ₁

^∗

⁽³⁾ =

∫ t

f

t

₀

(g _i X _i ⁽³⁾ dt (17)

J ₂

^∗

⁽³⁾ = 1 2

∫ t

_f

t

₀

(X _i ⁽¹⁾ H _ij X _j ⁽²⁾

+ X _i ⁽²⁾ H _ij X _j ⁽¹⁾ )dt (18)

以上より，拡張評価関数が最少となる必要条件は

J ₁

^∗

⁽¹⁾ = 0 (19) J ₁

^∗

⁽²⁾ + J ₂

^∗

⁽³⁾ = 0 (20)

上記の

2

条件から，

J

^∗

− J

^∗

⁽⁰⁾ = J ₂

^∗

⁽²⁾ ≥ 0 (21)

となる．また，本研究で用いた最小化手法は準ニュートン法である．ヘシアンマトリックスを更新する手法としては

BFGS

法を用いた．

3 検証

同定手法の検証を行う為，地盤を単純化して考え

3

次元の矩形有限要素分割図を用いる．計算条件として，図

2

のようなものを適用させた．外力は両問題とも

X = 200[m]

の側面に等分布荷重を載荷し，任意の地点で速度データを採取し，それを観測データとして弾性係数の同定を行う．有限要素法の順解析に用いる外力は模擬発破外力として単位時間に外力を加え，その衝撃によって発生する振動を任意の位置で採取し，それを観測データとして用いる．実際の現象においては観測データのみによってパラメータ同定を行うが，検証では目的とするパラメータを設定して有限要素法で順解析を計算してその速度データを用いて任意に設定した初期値から同定を行う．有限要素分割は節点数

1,029，要素数 4,320

となっている．境界条件としては

X = 0.0[m]

の面の変位を固定し，片持ち梁を想定している．

表

1:

設定条件

弾性係数

(初期値)

密度ポアソン比

[kN/m

²

] [kg/m

³

]

層

1 1.0 × 10 ⁶ 2.3 × 10 ³ 0.3

図

2:

有限要素分割図

(3)

4 検証結果

図

5

と図

6

は，それぞれ一次随伴法と二次随伴法を用いた際の弾性係数の推移を表している．また，図

3

と図

4

はその評価関数の推移を示したものである．両手法とも評価関数が

0

に収束しパラメータを同定することができた．検証段階での観測速度は，アルゴリズムの有効性を確認するため，数値解析にて得られたものを観測速度として定義し同定を行った．そのため，評価関数が

0

に収束したということは計算速度と観測速度は完全なる一致をしたということである．つまり，本手法は

3

次元におけるパラメータ同定の方法として有効であることが確認できた．

No. of iteration

Performancefunction

20 40 60 80

0 0.2 0.4 0.6 0.8 1

図

3:

評価関数の推移

(一次

随伴)

Performancefunction

10 20 30 40 50 60 70

0 0.2 0.4 0.6 0.8 1

図

4:

(二次

随伴)

Elasticmodulus[kN/m**2]

20 40 60 80

2E+06 4E+06 6E+06 8E+06 1E+07

図

5:

弾性係数の推移

(一次

随伴)

10 20 30 40 50 60 70

2E+06 4E+06 6E+06 8E+06 1E+07

図

6:

(二次

随伴)

5 _{数値解析例}

解析モデルとして，広島県に位置する大万木トンネル付近の地形データから，有限要素分割図を作成した

(図 7

参照)．総節点数と要素数はそれぞれ

8039

および

42023

となっている．また，解析領域はそれぞれの方向に

X = 254 [m]，Y = 254 [m]，Z = 205 [m]

となっている．観測値としての地盤振動は，トンネル坑口付近に設けた観測点から得られたものを使用する

(図 8

参照)．トンネル切羽にダイナマイトを仮定した衝撃力を与え，それによって生じた弾性波を用いて各層の弾性係数を同時に同定する．

発破外力の設定方法であるが，切羽に与えた外力は以下の式

(22)

で表される．外力を与える場所としては，図

10

図

7:

有限要素分割図図

8:

観測点

に示すように芯抜きを想定して，切羽中心の

6

点に

iX， Y

そして

Z

方向それぞれ載荷した．

Ω ˆ _αi =

∫

Γ

₂

N _α A _i (e

⁻

^ξt − e

⁻

^ηt )dΓ, (22)

ここに，A

i

は

i

方向において与えられる発破外力の最大振幅を意味する．ξと

η

は爆破力の時間変化に依存するパラメータである．本研究では弾性係数の同定を行う前に最大振幅

A i

の同定を行っている．また，ξと

η

をそれ

ぞれ

1000，5000

とした．解析領域はその解析目的に応じ

図

9:

トンネル坑口

図

10:

切羽における外力

て，トンネルの水平方向および鉛直方向に十分な範囲を設定し，それぞれの境界には適切な境界条件を設定しなければならない．そのため本研究では，領域側面はそれぞれ法線方向のみ変位を固定，領域底面は水平および鉛直方向全ての変位を固定，自由表面境界には応力自由条件を適用した．以上の解析条件で，トンネル掘削による発破振動を用い各層の弾性係数の同定を行う．

6 数値解析結果

図

11

と図

13

は数値解析例の評価関数の繰返し計算履歴を示している．両手法において評価関数は初期値を用いて計算された値から，ある一定の値に収束していることが見てとれる．評価関数が完全に

0

に収束していない理由としては，計算値と計測値が完全に一致していないためであると考えられる．計測データは機械的誤差や人為的誤差を含んでいるために数値解析と完全に一致させ

(4)

ることは難しいと考えられる．つまり，評価関数はゼロにはならないが，最小となるものが最適なパラメータとして判断されるのである．このことより，評価関数をより減少させることができた二次随伴法の方が有効であると言える．更に，図

12

と図

14

はそれぞれ一次随伴法と二次随伴法を用いた際の弾性係数の推移を表している．

Performancefunction

20 40 60 80

0.4 0.5 0.6 0.7 0.8 0.9 1

図

11:

(一

次随伴)

20 40 60 80

8E+06 1E+07 1.2E+07

1.4E+07 Stratum 1

Stratum 2 Stratum 3

図

12:

(一

次随伴)

Performancefunction

10 20 30 40 50

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

図

13:

(二

次随伴)

10 20 30 40 50

8E+06 9E+06 1E+07 1.1E+07

stratum 1 stratum 2 stratum 3

図

14:

(二

次随伴)

表

2:

弾性係数の計算結果

弾性係数

[kN/m ² ]

初期値一次随伴二次随伴層

1 1.00 × 10 ⁷ 1.40 × 10 ⁷ 1.15 × 10 ⁷

層

2 1.00 × 10 ⁷ 7.78 × 10 ⁶ 8.26 × 10 ⁶

層

3 1.00 × 10 ⁷ 7.73 × 10 ⁶ 7.43 × 10 ⁶

Time [s]

Velocity[m/s]

0.1 0.2 0.3 0.4

-0.001 -0.0005 0 0.0005 0.001

Computed velocity in X-direction Observed velocity in X-direction

図

15:

速度の比較-X方向

(一次随伴)

Time [s]

Velocity[m/s]

0.1 0.2 0.3 0.4

-0.001 -0.0005 0 0.0005 0.001

Computed velocity in X-direction Observed velocity in X-direction

図

16:

速度の比較-X方向

(二次随伴)

Time [s]

Velocity[m/s]

0.1 0.2 0.3 0.4

-0.001 -0.0005 0 0.0005 0.001

Computed velocity in Y-direction Observed velocity in Y-direction

図

17:

速度の比較-Y方向

(一次随伴)

Time [s]

Velocity[m/s]

0.1 0.2 0.3 0.4

-0.001 -0.0005 0 0.0005 0.001

Computed velocity in Y-direction Observed velocity in Y-direction

図

18:

速度の比較-Y方向

(二次随伴)

Time [s]

Velocity[m/s]

0.1 0.2 0.3 0.4

-0.001 -0.0005 0 0.0005 0.001

Computed velocity in Z-direction Observed velocity in Z-direction

図

19:

速度の比較-Z方向

(一次随伴)

Time [s]

Velocity[m/s]

0.1 0.2 0.3 0.4

-0.001 -0.0005 0 0.0005 0.001

Computed velocity in Z-direction Observed velocity in Z-direction

図

20:

速度の比較-Z方向

（二次随伴）

7 結論

数値解析の結果から，地表面に設置し実際に得られた弾性波を用いて，弾性係数を三層同時に同定することができた．従来までは，一次随伴方程式法を用いて弾性係数の同定を行うことは現存していたが，二次随伴方程式法を用いた同定までには手が及んでいなかった．その現状から本研究に取り組み，パラメータの同定解析手法を確立することができた．

参考文献

[1] K.Arai, H. Ohta and T. Yasui : “Simple optimiza- tion techniques for evaluating deformation moduli from ﬁeld observations”, Soils Found., 23(1), pp107 - 113, 1983

[2] N.Koizumi and M.Kawahara : “Parameter Identiﬁ- cation Method for Determination of Elastic Modu- lus of Rock Based on Adjoint Equation and Blasting Wave Measurements”, Int. J. Numer. Analy. Meth.

Parameter Identiﬁcation of Elastic Modulus of Rock Based on Blasting Waves in Tunnel Excavation

(2011

)

トンネル掘削における発破振動を用いた弾性係数の同定に関する研究

Parameter Identiﬁcation of Elastic Modulus of Rock Based on Blasting Waves in Tunnel Excavation

34

Yuto MOTOYAMA

1 はじめに

(以下，大万木トンネル)

4, 878m

3

2 パラメータ同定

2.1 一次随伴法

J

J

= J + Λ + Ξ, (1)

(1)

J

J = 1 2

∫ t

t

( ˙ u αi − η αi )Q αiβj ( ˙ u βj − η βj )dt, (2)

u αi

η αi

α

i

αiβj

0

t f

(3)

D ijkl u k,lj + ρb i − ρ¨ u i = 0, (3)

i

i

D ijkl

(図 1

ijkl

D ijkl =

∑ M

m=1

D (m) ijkl , (4)

(5)

Λ =

∫ t

t

λ αi ( ˆ Ω αi − M αiβj u ¨ βj

− C αiβj u ˙ βj − K αiβj u βj )dt, (5)

(5)

λ αi

Lagrange

M αiβj

αiβj

αiβj

Ω ˆ αi

(1)

Ξ

Ξ = 1 2

∫ t

t

(X λ − X ˜ λ )W λµ (X µ − X ˜ µ )dt, (6) X λ

αiβj

(6)

1:

2.2 二次随伴法

J

J

= J +

∫ t

t

λ αi (M αiβj u ¨ βj + C αiβj u ˙ βj

+ K αiβj u βj − Ω ˆ αi )dt, (7)

u i = u (0) i + u (1) i + u (2) i + u (3) i + ..., (8) λ αi = λ (0) αi + λ (1) αi + λ (2) αi + λ (3) αi + ..., (9)

X i

X i = X i (0) + X i (1) + X i (2) + ..., (10)

J

J

= J

(0) + J

(1) + J

(2) + J

(3) + ..., (11)

1 _はじめに

D _ijkl u _k,lj + ρb _i − ρ¨ u _i = 0, (3)

D _ijkl =

D ^(m) _ijkl , (4)

− C _αiβj u ˙ _βj − K _αiβj u _βj )dt, (5)

M _αiβj

(X _λ − X ˜ _λ )W _λµ (X _µ − X ˜ _µ )dt, (6) X _λ

u _i = u ⁽⁰⁾ _i + u ⁽¹⁾ _i + u ⁽²⁾ _i + u ⁽³⁾ _i + ..., (8) λ _αi = λ ⁽⁰⁾ _αi + λ ⁽¹⁾ _αi + λ ⁽²⁾ _αi + λ ⁽³⁾ _αi + ..., (9)

X _i

X i = X _i ⁽⁰⁾ + X _i ⁽¹⁾ + X _i ⁽²⁾ + ..., (10)

⁽⁰⁾ + J

⁽¹⁾ + J

⁽²⁾ + J

⁽³⁾ + ..., (11)

⁽⁰⁾ +

= J ₁

⁽⁰⁾ + J ₁

⁽¹⁾ + J ₁

⁽²⁾

+ J ₂

⁽²⁾ + J ₁

⁽³⁾ + J ₂

⁽³⁾ + ... (13)

J ₁

⁽¹⁾ =

(g i X _i ⁽¹⁾ dt (14)

J ₁

⁽²⁾ =

(g i X _i ⁽²⁾ dt (15)

J ₂

⁽²⁾ = 1 2

(X _i ⁽¹⁾ H ij X _j ⁽¹⁾ )dt (16)

J ₁

⁽³⁾ =

(g _i X _i ⁽³⁾ dt (17)

J ₂

⁽³⁾ = 1 2

(X _i ⁽¹⁾ H _ij X _j ⁽²⁾

+ X _i ⁽²⁾ H _ij X _j ⁽¹⁾ )dt (18)

J ₁

⁽¹⁾ = 0 (19) J ₁

⁽²⁾ + J ₂

⁽³⁾ = 0 (20)

⁽⁰⁾ = J ₂

⁽²⁾ ≥ 0 (21)

1 1.0 × 10 ⁶ 2.3 × 10 ³ 0.3