Microsoft Word - 特集準備資料１

圧密問題の逆解析－

1－

1 はじめに 地盤工学における逆解析の研究は，主として圧密沈下の問題を対象に進められてきた．この理由は定 かではないが，日本の都市・人口が軟弱な粘性土地盤で構成される沖積平野に集中しており，建設工事 にともなう地盤の（圧密）沈下量を精度よく予測したいというニーズがあったことと，地盤の沈下量が 比較的計測しやすく，逆解析の研究に必要な観測データが容易に得られたということが理由として考え られる．圧密沈下問題の逆解析に関する研究は 1970 年代後半にはじまり，コンピュータの発達により 1990 年代中ごろまで（ある程度）盛んに研究された．最近では，遺伝的アルゴリズム（Genetic Algorithm, GA）や粒子群最適化法（Particle Swarm Optimization (PSO) method）などの新しいアルゴリズムを圧密逆 問題に適用した研究などがみられる．しかしながら，圧密の逆問題に関する多くの研究例があるものの， これまで「あるアルゴリズムでパラメータが同定でできるかどうか」や「逆解析に基づいた将来予測が 観測データに合っているかどうか」といった議論が中心であり，圧密逆問題の性質について定量的に示 された，あるいは議論された研究例は少ない（例えば[1]，[2]）．そこで圧密 WG では，活動の第一ステ ップとして，最も基本的な一次元圧密問題を対象とし，その逆問題の性質について定量的に調べること とした．本稿ではその結果を報告する． 2 圧密現象の定式化 圧密とは，飽和した地盤材料，主として粘性土に荷重を載荷したときに，時間経過に伴いその変形， とくに鉛直変位（沈下）が進行する現象を指す．砂質土や礫質土と異なり，粘性土は透水性が低く間隙 中の水が瞬時に移動できないことから，時間遅れの変位が生じる．圧密現象は，Terzaghi [3]や三笠[4] による圧密方程式により記述することができる．ここでは（著者の土質力学の復習も兼ねて）一次元圧 密方程式の誘導方法を示す[5]． 2.1 飽和粘土の圧密挙動[5] Fig. 1 は，飽和粘土の一次元圧密の模式図を表している．容器の中に土の骨格を表すバネと間隙水を 表す水が入っている．上蓋に小さい孔が空いており，排水に時間がかかることを表している．図(a)は載 荷直後の状態を表しており，土骨格は縮まず排水は生じていない．よって釣合い条件から，間隙水圧 u は上載荷重p に等しい．載荷から時間が経過し，排水とバネの圧縮が進んだ状態が図(b)に示されている． 上載荷重の一部を土骨格（バネ）が受け持ち，残りを間隙水圧が受け持つ状態であり，次式が成り立つ．

u

p

'

(1) 図(c)は十分時間が経過した排水完了時点の状態を表しており， ' = p, u = 0 が成り立つ．この説明から分 かることは，飽和粘土の変形は，p によって規定されるのではなく， 'によって規定されるということ

である． P=pA _{P P} p p p

z

0

z

z

f 0 t p u , 0 ' t p u ' t 0 , ' p u drainage (a) (b) (c) Fig. 1 飽和粘土の一次元圧密[5] 2.2 応力－ひずみ関係[5] 先述した飽和粘土の一次元圧密問題において，載荷荷重（圧力）p を増加させると排水も増え，圧密 完了時の圧縮量や圧縮ひずみ （(z0‐zf) / z0）も大きくなる．載荷圧力p と圧縮ひずみ の関係の模式図を Fig. 2 に示す．

)

'

(

p

O

p

Fig. 2 一次元圧縮時のひずみと有効応力[5] ここで，ひずみ の代わりに間隙比e を用いて Fig. 2 の応力‐ひずみ関係を増分形式で表すと，

p

e

e

p

m

)

1

(

₀ v (2) となり，mvは体積圧縮係数と呼ばれる．図から理解できるように，mv は有効応力が増加すれば小さく なる．

2.3 一次元圧密理論[5] 粘土層下端を原点として上向きにz 軸をとり，点 z で断面積が A，高さ z の（円筒）要素を考える（Fig. 3 (a)）．粘性土の間隙水の排水速度を v とすると，要素下面ではv( zt, )であり，上面ではv ( v/ z) zで あるため， t中の排水（脱水）体積は， Vw ( v/ z)A t zと表せる．一方，要素の鉛直ひずみを ( zt, ) とすると， t中の要素の体積変化は， V ( / t)A t zと表せる．土粒子と水の圧縮性を無視すると， w V V であるから，

x

v

t

(3) が成り立つ． 一方，点z の静水圧，過剰間隙水圧をそれぞれ us，u とすれば，全水頭 h = z + us/ w + u/ wである．右 辺の1 項と 2 項の和は地下水面と粘土層底面間の鉛直距離であり，一定である．よって，Darcy の法則 から，

z

u

k

z

h

k

ki

v

w (4) が成り立つ．ここに，k は透水係数，i は動水勾配である． 0 z z Dz t A z z v v t vA 0 z z Dz A z z ' ' A ' ) (j i w z jA (a) (b) Fig. 3 一次元脱水と作用力[5] 次に，Fig. 3 (b)から，土粒子に働く有効力のつり合いを考える．図中の 'は有効応力であり，j は透水 力（浸透力）を表す．土粒子に働く有効物体力を無視すると， '/ z i wであり，式(4)から，

v

k

i

z

w w

'

(5) となる．一次元圧縮条件に対しては， mv 'であるので，mvがe に依存するとしても z に無関係な らば式(5)は，

z

m

k

v

w v (6) となり，式(6)を式(3)に代入すると，

z

m

k

z

t

_{v w} (7) が導かれる．これは三笠の圧密方程式[4]と呼ばれる． 一方，式(7)に

m

_v

'

を代入すると，

z

k

z

t

m

'

'

w v (8) が得られる．ここで，載荷圧p = p0が一定であれば，式(5)は間隙水圧に基づいた次式で表される． 2 2 v

_z

u

c

t

u

(9) ここにcvは圧密係数と呼ばれ， w v v

m

k

c

(10) で表される．式(9)は Terzaghi の圧密方程式[3]と呼ばれる．

3 1 次元圧密方程式を対象とした逆問題 ここでは，Fig. 4 に示すモデル地盤を対象に「地盤内の間隙水圧を測定し，初期の間隙水圧分布を求 める」という逆問題を考える．ここでは圧密方程式を対象とした逆問題の性質を分析することが目的で あるため，問題設定にはこだわらないことにする．なお，ここで対象とする拡散方程式（圧密方程式） の逆問題は小國[6]に取り上げられており，本稿における計算も小國[6]を参考にした． 5. 0m observation Dx = 0.125m

200

sec

000

,

40

~

0

200

)

kPa

(

50

'

)

cm/sec

(

10

1

8

.

0

6

.

1

T 5 0

N

T

t

p

k

C

e

c 10.0kPa Fig. 4 モデル地盤 3.1 一次元圧密方程式の離散化 Terzaghi の圧密方程式（支配方程式）と初期条件，境界条件は以下の通りである。 支配方程式： 2 2 v

_z

u

c

t

u

(

0

m

z

5

m

,

0

t

)

(11) 境界条件：

0

)

,

( t

z

u

(at z

0m,

5m)

(12) 初期条件：

)

(

)

0

,

(

z

f

x

u

(13)

ここでcvは先述した圧密係数である． 式(11) ~ (13)で表される初期値境界値問題を数値的に解くために，時間 t と空間 x を有限の等間隔の格 子で離散化する．離散化の間隔について，時間については t，空間については x とし，前進差分により 式(11)を離散化すると，

z

z

t

z

z

u

t

z

u

z

t

z

u

t

z

z

u

c

t

t

z

u

t

t

z

u

(

,

)

(

,

)

(

,

)

(

,

)

(

,

)

(

,

)

v 2 v

)

,

(

)

,

(

2

)

,

(

z

t

z

z

u

t

z

u

t

z

z

u

c

(14) となる．上式より，各計算点（格子点）の間隙水圧ui,jの時間発展式が次のように得られる． _ij i j i j i j

t

u

_ij

z

u

u

u

c

u

1, ₂, 1, _, v 1 ,

2

j i j i j i

u

z

t

c

u

z

t

c

u

z

t

c

, 1 2 v , 2 v , 1 2 v

₍

₁

₂

₎

j i j i j i

u

u

u

_1,

(

1

2

)

_{, 1,} (15) ここに， = cv t / ( z2)であり，0 < < 0.5 を満たす必要がある．式(15)はマトリックス表現を用いて以 下のように表すことができる。 j j

AU

U

₁ (16) ここに，U j + 1は時刻j + 1 における節点 xiでの間隙水圧を1 列にならべたベクトルで， 1 , 1 , 2 1 , 1 1 p j N j j j

u

u

u

U

(17) と表される．ここにNpは全節点数を表す．なお，Ujに関しても同様の形で表現される． また，マトリックスA は次式で表される正方対称行列である．

2

1

0

0

0

0

2

1

0

0

2

1

0

0

2

1

A

(18)

式(16)に時刻 t = 0 での U，すなわち U0を代入し，各時間ステップで行列A をアップデートしていくこ とにより，n 番目の時間ステップにおける任意の節点での水圧が次式で計算できる． 0 1

A

U

U

n j (19) Fig. 5 に Fig. 4 の条件で計算した間隙水圧の消散過程を示す．図から明らかなように，載荷直後（計算 初期）では地盤の間隙水圧は上載圧に等しいが，時間の経過に伴い，水圧が消散していく様子が再現で きていることがわかる． 0.0 2.0 4.0 6.0 8.0 10.0 12.0 -5.0 -4.0 -3.0 -2.0 -1.0 0.0 : 0 : 2000sec : 10000sec : 30000sec : 60000sec : 100000sec

D

ep

th

(

m

)

Excess pore pressure (kPa)

Fig. 5 間隙水圧の消散過程 3.2 一次元圧密の逆問題 設定した逆問題は，「地盤内の間隙水圧を測定し，初期の間隙水圧分布を求める」であったが，具体 的には，式(19)における係数行列 Anの性質を調べる問題に相当する． 以上の準備に基づき，ここで取り扱う逆問題を次式で表す．

w

Mv

(20) ここに，

p 2 1 N

u

u

u

M

(21),

v

U

₀ (22), t p t , obs , 2 obs 1 , 2 obs 1 , 1 N N N

u

u

u

u

w

(23) n n n

HA

A

u

1

0

0

0

0

1

0

0

0

1

(24) で表され， obs , j i u は節点i における時刻 j の間隙水圧の観測データ，Ntは解析時間ステップ，H は観測行列 を表す．例えばFig. 4 で示したモデル地盤のように，境界から 1 節点内側の節点 2 点の観測データを用 いた場合，式(23)，(24)は次式のように書き換えられる． t p p p , 1 obs 2 , 1 obs 2 , 1 obs 1 , 1 obs 1 , 1 N N N N

u

u

u

u

u

w

(25),

u

_n

A

n

1

0

0

0

0

0

0

1

(26) 以上より，ここで対象とする逆問題は，連立一次方程式からv を求める問題，すなわち行列 M の性質を 調べる問題ということがわかる． 3.3 特異値分解 ここでは行列M の性質を調べる道具として特異値分解を用いる．そこでまず，特異値分解について概 説する．n×m 行列 M が与えられたとする．MTM という半正定値行列をつくり，左右から直交行列 V を掛け対角化したものを次式で表す． 2 2 3 2 2 2 1

0

0

0

0

0

0

0

0

0

0

0

0

0

0

)

(

M

M

V

V

T T (27)

上式のV の縦ベクトル vjは，MTM の固有値 j2に対応する固有ベクトルである． ここで，VTMT = (MV)Tであり，W = MV とおく，すなわち， j j

Mv

w

(28) とおけば，i ≠ j のとき wiTwj = 0，i = j のとき wiTwj = i2となる．よって，直交行列W により， 2 T

_MV

_Γ

W

(29) と表すことができる．ここで， j j j

w

/

u

(30) とし，ujからなる行列U = (u1,…, um)を用いると，

0

0

0

0

0

0

0

0

0

0

0

0

0

0

3 2 1

MV

U

T (31) となる．このような分解を特異値分解M = U VT_という． 実際に特異値分解を行う場合，理論的には（原理的には）m×m の正方対称行列 MMT_{（もしくはn×n} の正方対称行列 MT_{M）の固有値解析により特異値が求められるが，数値計算上は好ましくないとされ} ている．その理由として，MT_{M を作る段階で多くの情報が失われてしまうこと[7]や，M が条件の良い} 行列であってもMT_{M も条件が良いとは限らない[8]などが挙げられている．よってここでも正方対称行} 列の固有値分解を行うのではなく，M を直接特異値分解し行列の性質を調べることとした．特異値分解 は種々の方法が提案されているが，ここではハウスホルダー法による2 重対角化と Givens 回転による追 い込みを組み合わせた方法を採用する．具体的な計算方法については文献[7]や[9]を参照されたい． 3.4 結果 はじめに，Fig. 4 の問題における行列 M を特異値分解し，得られた特異値の絶対値の分布を Fig. 6 に 示す．特異値はn の増加に伴い指数関数的に減少し，n = 20 あたりでほぼ 0 となる．これは，Fig. 4 の 問題に対して初期の間隙水圧分布を完全に求めるためには39 個の項が必要となるが，20 項程度しか使 用できないことを意味している．Fig. 7 ~ 10 に特異値分解で得られた全ての固有ベクトルを示す．n = 1 ~

7 程度の低次の固有ベクトルは比較的単調な（なめらかな）カーブを示しているが，次数の増加に伴っ て複雑な形状になっていくことが分かる．逆問題の解はこれら固有ベクトルの線形和で表されることか ら，低次の固有ベクトルのみ有効な場合には解のおおまかな形しか推定できず，逆に高次の固有ベクト ルも有効な場合は解の詳細な形まで推定できることを意味している．

0

10

20

30

40

0.0

0.2

0.4

0.6

0.8

1.0

S

in

gu

la

r

va

lu

e

|

n

|

n

Fig. 6 一次元圧密方程式の特異値（絶対値） 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 1 n = 2 n = 3 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 4 n = 5 n = 6 Fig. 7 固有ベクトル-1

0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 7 n = 8 n = 9 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 10 n = 11 n = 12 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 13 n = 14 n = 15 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 16 n = 17 n = 18 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 19 n = 20 n = 21 Fig. 8 固有ベクトル-2

0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 22 n = 23 n = 24 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 25 n = 26 n = 27 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 28 n = 29 n = 30 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 31 n = 32 n = 33 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 34 n = 35 n = 36 Fig. 9 固有ベクトル-3

0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n 0 10 20 30 40 -1.0 -0.5 0.0 0.5 1.0 n n = 37 n = 38 n = 39 Fig. 10 固有ベクトル-4 3.5 観測時間・観測間隔・観測数の影響 3.4 では Fig. 4 の条件で得られた M 行列の特異値分解により一次元圧密方程式の逆問題としての性質 を分析した．M 行列は観測の方法に影響されるため，ここでは観測時間や観測間隔，観測点数の違いに より，逆問題（特異値の分布）がどのように変わるかを調べる． Fig. 11 は観測間隔の違いによる特異値の絶対値の分布の違いを表している． t は観測データが得られ る時間間隔を表しており，100sec は 200sec の 2 倍の観測データが得られていることを意味している．観 測データが2 倍であるにも関わらず，逆問題の性質は改善しておらず，n の増加に伴い特異値は急激に 減少している．よって，計測間隔を短くすることは，一次元圧密の逆問題の改善にはつながらないこと がこの結果からわかる．

0

10

20

30

40

0.0

0.2

0.4

0.6

0.8

1.0

: t = 200 sec

: t = 100 sec

S

in

gu

la

r

va

lu

e

|

n

|

n

Fig. 11 観測間隔の影響 Fig. 12 は観測点数の違いによる特異値の絶対値の分布の違いを表している．図中の Nobsは観測データ 点数を表しており，観測点数以外の条件はFig. 4 と同じである．Nobs = 1 のケースは，地表面から 1 節点 内側の1 点のみの観測データを用いた場合の結果であり，Nobs = 2 のケースは Fig. 4 と同じ条件，Nobs = 3 のケースはFig. 4 に示された観測点に加え，地盤中心における観測点 1 点を加えたものである．図から

明らかなように，観測点の増加に伴い，特異値の減少傾向が改善されていることがわかる。しかしなが ら，指数関数的な減少傾向は改善されなかった．

0

10

20

30

40

0.0

0.2

0.4

0.6

0.8

1.0

: N

_obs

= 2

:

N

_obs

= 1

:

N

_obs

= 3

S

in

gu

la

r

va

lu

e

|

n

|

n

Fig. 12 観測データ数の影響 Fig. 13 は観測位置の違いによる特異値の絶対値の分布の違いを表している．図中の Nobs=1 のケースは

Fig. 12 における Nobs = 1 と同じ結果であり，Nobs = 1 at center は地盤中心の観測データ 1 点を用いた結果

を示している．地盤中心のデータを用いた場合，境界付近のデータを用いた場合よりも特異値の減少傾 向がわずかに改善されているが，両結果ともにn = 10 あたりで特異値がほぼゼロとなっている．

0

10

20

30

40

0.0

0.2

0.4

0.6

0.8

1.0

: N

_obs

= 1

:

N

_obs

= 1 at center

S

in

gu

la

r

va

lu

e

|

n

|

n

Fig. 13 観測データ位置の影響

Fig. 14 は観測期間の影響を表しており，観測期間以外の条件は Fig. 4 と同じである．T = 60,000 sec の

向は改善されていないことがわかる．よって，一次元圧密問題においては長期間の観測は逆問題の改善 に有効ではないことがわかる．

0

10

20

30

40

0.0

0.2

0.4

0.6

0.8

1.0

: T = 40,000 sec

: T = 60,000 sec

S

in

gu

la

r

va

lu

e

|

n

|

n

Fig. 14 観測期間の影響 Fig. 15 は観測開始時期の影響を表している．これまでと同様に，観測開始時期以外の条件は Fig. 4 と 同じである．両ケースともに時間40000sec の間に 200sec 間隔で計測された計 200 個の観測データを用 いてM マトリックスが構成されていることになるが，後半の観測データのみを用いたケースでは，有効 な解がほとんどないことがわかる．これは一次元圧密の逆問題を解くにあたり，圧密初期の観測データ が重要であり，ある一定時間を過ぎて観測されるデータは，逆問題においては有効ではないことを意味 している．

0

10

20

30

40

0.0

0.2

0.4

0.6

0.8

1.0

:

T = 0 ~ 40,000

:

T = 40,000 ~ 80,000 sec

Si

ng

ul

ar

v

al

ue

|

n

|

n

Fig. 15 観測開始時期の影響

これまで示したように，一次元圧密の逆問題は，観測計画を変更してもその性質を改善することは難 しく，n の増加に伴って特異値が指数関数的に減少するという特徴を持つことが明らかとなった．今回 は単純な一次元圧密方程式を対象としたが，このような性質が弾塑性構成モデルを用いた土‐水連成有 限要素解析を対象とした逆問題にも当てはまるのかどうか，今後検証が必要である．

参考文献

[1] Honjo, Y., Liu, W. T. and Guha, S.: Inverse analysis of an embankment on soft clay by extended Bayesian method, Int. J. Numer. Metho. Geomech., Vol. 18, pp. 709-734, 1994.

[2] Honjo, Y., Liu, W. T. and Sakajo, S.: Application of Akaike information criterion statistics to geotechnical inverse analysis: The extended Bayesian method, Structural Safety, Vol. 14, pp. 5-29, 1994.

[3] Terzaghi, K.: Theoretical Soil Mechanics, John Wiley & Sons, Inc., 1943. [4] 三笠正人：軟弱粘土の圧密，鹿島出版会，1963． [5] 山口柏樹：土質力学，技報堂出版，1984. [6] 小國健二：応用例で学ぶ逆問題と計測，オーム社，2011. [7] 村田健朗：線形代数と線形計算法序説，サイエンス社，1986． [8] 齋藤正徳：特異値分解と最小二乗法，物理探査，Vol. 36，No. 3，pp. 36 - 48，1983. [9] 中川 徹，小柳義男：最小二乗法による実験データ解析，東京大学出版会，1982. 2014. 7. 28