日本地震工学会論文集 第 7 巻、第 2 号（特集号）、2007

2003 年十勝沖地震における地盤の非線形応答：KiK-net 鉛直アレー記録の活用

山本明夫 1)、笹谷 努 2)
1) 正会員 応用地質株式会社 東京本社 技術センター、上級専門職 博士（理学）e-mail: yamamoto-akio@oyonet.oyo.co.jp
2) 正会員 北海道大学大学院理学研究院、助教授 理博
e-mail: sasatani@ares.sci.hokudai.ac.jp

要約
本論文は、KiK-net 鉛直アレー観測記録から地盤の非線形応答を評価したものである。取り上げた KiK-net 観測点は北海道十勝支庁豊頃町にある TKCH07 で、通常の KiK-net 観測点とは異なり、表層 30m 程度まで S 波速度 200m/sec 以下の軟弱な地盤が分布している。
ここでは 2003 年十勝沖地震（MJ=8.0）で地表最大加速度(EW 成分)=403.9cm/sec2 が観測されたほか、弱震動も多数観測されている。本論文では、ほぼ直下に発生した 22 地震を用いて、地表／地中の S 波スペクトル比を算出し、それを用いて地中地震計以浅の地盤の S
波速度構造および Q 値を遺伝的アルゴリズム (GA) を用いて同定した。また、その同定モデルを用いて線形地震応答解析、等価線形地震応答解析および非線形地震応答解析を行い、その結果をもとに、観測記録と良い一致を示す改良型等価線形地震応答解析手法を提案する。

キーワード： KiK-net 鉛直アレー、地震時非線形性、地盤モデル同定、
遺伝的アルゴリズム、地震応答解析、2003 年十勝沖地震
１．はじめに
地表での地震動特性は、断層での震源特性、震源から地震基盤までの伝達特性および地震基盤から地表までのサイト特性に支配される。この中でサイト特性、特に工学的基盤以浅の表層地盤の非線形性については、
適切な鉛直アレー観測記録に基づいた研究があまり行われておらず、地表地震動に与える非線形性の影響を的確に予測することは現在でも難しい状況であると考えられる。過去の研究では、Satoh et al. (2001) 1) は、1995
年兵庫県南部地震の本震・余震の地表－地中鉛直アレー記録を用いて、地表と工学的基盤のスペクトル比から逆解析を行い地盤モデルを同定するとともに、強震記録から得られる剛性率比と減衰比のせん断ひずみ依存性が室内実験とよく一致していることを示している。Iai et al. (1995) 2) は、1993 年釧路沖地震の釧路港での地表－地中鉛直アレー記録を用い、近傍の地盤調査、凍結サンプリングおよび室内動的試験を行い、一次元
有効応力解析により密な砂地盤のサイクリック・モビリティによる加速度波形の変化を再現している。
一方、鉛直アレー観測網である KiK-net は、2000 年 8 月から正式公開となり、全国で約 670 点配置されて-144-

いる。この鉛直アレー観測網の記録を利用することは、非線形応答を評価するためには非常に有用と考えられるが、観測点がほとんど山地に位置するため現在までそれほど積極的には利用されてこなかった。本研究でとりあげた 2003 年十勝沖地震は、KiK-net 設置後最初のマグニチュード 8 級の地震であり、北海道内の
KiK-net 観測点 112 点のうち 106 点で記録が得られ、その中の最大地表加速度は、KSRH03（標茶北）の 806.2gal（NS 成分）である。本論文では、液状化しない地点での大ひずみ時の地盤の非線形応答の問題を対象とし、2003 年十勝沖地震などの北海道 KiK-net 観測点の記録を用いて、非線形応答を示す観測点を抽出するとともに、地震応答解析において結果に最も大きな影響があると考えられる地盤構造について、弱震動記録による
観測地表／地中 S 波スペクトル比を用いて逆解析により速度層構造を同定する。そして、その構造を用いて種々の手法で地震応答解析を行い、地盤の非線形応答が地震動に与える影響およびその予測方法について検討する。
２．2003 年十勝沖地震時に地盤の非線形応答を示す KiK-net 観測点の抽出
2003 年十勝沖地震は、9 月 26 日午前 4 時 50 分頃発生し、気象庁マグニチュード 8.0、震源の深さ約 42kmであった。北海道内の最大震度は 6 強であり、道東から道央南部にかけて大きな震度となっている。図 1 に
北海道南部の KiK-net 観測点と非線形応答を示す観測点の抽出に用いた地震の震央位置を示す。ここでは、2003 年十勝沖地震による KiK-net 鉛直アレー記録について、その非線形性を検討するために、S 波スペクトル比（地表／地中）を算出し、同様な処理により算出した弱震時のそれと比較する。弱震時の記録としては、
KiK-net で観測された記録数も考慮して、図 1 に示す 3 地震によるものを対象とした。
スペクトル比の算出には、理論スペクトル比の計算に SH 波による重複反射理論を用いることから、地表および地中の波形データをトランスバース成分の波形に変換した S 波を用いた。その後、弱震動 3 地震については、S 波主要動部を 10.24 秒間切り出し、その両端に 10%のコサインテーパーをかけ、それぞれの S 波スペクトルから地表／地中スペクトル比を計算した。
本震については、同様にS 波主要動部を20.48 秒切り出した。
非線形を示す観測点の抽出は、スペクトル比の高周波数側の低下、およびスペクトル比卓越周波数の低周
波数側への移動により判断した（Higashi & Sasatani (2000) 3) ）。図 2a～図 2c に TKCH07（豊頃）観測点の強震動（本震）と弱震動（図 1 の 3 地震）のスペクトル比の比較図を示した。異なる弱震動と比較しても同じように地盤の非線形性の特徴が現れているのが確認できる。このような強震動と弱震動のスペクトル比を北海道内の全 KiK-net 観測点で検討したが、山地部に設置されている観測点が多いため、スペクトル比の高周波数側の低下とスペクトル比卓越周波数の低周波数側への移動の両方が発生している観測点は、TKCH07（豊頃）、4)
TKCH08（大樹）および KSRH10（浜中）の 3 箇所であった（山本 (2006) ）。本研究では、これら 3 観測点
5)の中で、地質柱状図 において最も軟弱地盤が厚く堆積する TKCH07（豊頃）観測点について検討する。
KiK-net 豊頃観測点(TKCH07)は、北海道十勝支庁豊頃町にあり、通常の KiK-net 観測点とは異なり、十勝川の左岸堤防から約 450m 離れた旧湿地 6) に位置する。図 3 に TKCH07（豊頃）観測点付近の治水地形分類図を示す。この地点では、防災科学技術研究所の地盤情報 5) によると、表層から S 波速度 Vs=200m/sec 以下の
腐植土および砂質土が約 30m 分布している。
３．遺伝的アルゴリズム(GA)を利用した地盤モデルの同定
工学的基盤以浅の表層地盤では、層相および層構造が非常に複雑で、従来から行われてきた最小自乗法による逆解析では、適切な初期モデルの設定が困難な場合が多く、初期モデルによっては繰返し計算の途中で地盤の速度や層厚などの変数が負となる場合がある。このような問題を解決するために、最近、最適化問題の解法として注目されている遺伝的アルゴリズム(GA)を用いて、KiK-net 鉛直アレー観測記録の地表／地中ス
ペクトル比から、S 波速度構造および Q 値を同定することとした。これは、後で行う非線形応答計算において、正確な S 波速度構造が重要であることを反映したものである。
遺伝的アルゴリズム(GA)は、自然淘汰に基づく生物の進化過程を模擬した数理的なモデルであり、1980 年-145-

Spectral Ratio (Suf./Base)10050

TKCH07 2001-0814-M6.2
Weak-motion
Strong-motion(Main Shock)

10
5
1
0.5
0.1
0.1

0.5

1

5

10

Frequency

10050

TKCH07 2003-0929-M6.5
Weak-motion
Strong-motion(Main Shock)Spectral Ratio (Suf./Base)

Spectral Ratio (Suf./Base)

図 1 北海道 KiK-net 観測点位置および検討対象地震震央 図 2a 強震動と弱震動のスペクトル比の比較(1)

10
5
1
0.5
0.1
0.1

0.5

1

5

TKCH07 2003-1011-M6.1
Weak-motion
Strong-motion(Main Shock)

10
5
1
0.5
0.1
0.1

10

0.5

1

5

10

Frequency

Frequency

図 2b 強震動と弱震動のスペクトル比の比較(2)10050

図 2c 強震動と弱震動のスペクトル比の比較(3)

TKCH07十勝川

図 3 TKCH07（豊頃）観測点付近の治水地形分類図-146-

代後半からシステム設計などと組み合わせ、最適化問題の解法としての可能性が検討されてきた。この方法
は、地震学での逆解析問題にも広く利用されており、山中・石田 (1995) 7) はこの方法を用いて微動などの位
相速度データから S 波速度構造を同定している。本論文では、野口・笹谷 (2004) 8) による遺伝的アルゴリズ
ムと一次元重複反射理論を組み合わせ、観測スペクトル比と理論スペクトル比の差を最小とする S 波速度構
造と Q 値を同定した。S 波構造モデルの同定では、Q 値を暫定的に固定して GA を用いて層厚と S 波速度を同定し、次に、得られた最適な速度構造モデルを用いて再び GA で Q 値パラメータを同定した。このような
2 段階評価方法は、速度構造と Q 値を一度に同定すると自由度が高すぎて結果が安定しないためである。また、Q 値の場合は、Q=Q0・f n（f：周波数(Hz)）として、Q0 と n のトレードオフを避けるため、n=0.5, 0.6, 0.7に固定して Q0 を同定した。
地盤モデルの同定に用いた地震は、S 波の鉛直入射の仮定を満たすために、震源の深さが 50km 以深で、震央距離が震源深さの 2 倍以内の図 4 に示す 22 地震とした。これらの地震の観測波形よりトランスバース成分を抽出し、S 波主要動部を約 10 秒間切り出した。その波形の両端に 10%のコサインテーパーをかけ、それぞれのフーリエスペクトルを算出し地表／地中のスペクトル比を求めた。図 5 に示すように算出したスペクト
ル比からその相乗平均を求め、それを地盤モデルの同定に用いるターゲットスペクトル比とした。なお、このスペクトル比は、地表地震計の設置方向が+40 度ずれていることが判明しているためそれを補正している（前田ほか (2005) 9) ）。
GA を用いた速度構造の同定は、次の拘束条件を用いた。まず理論スペクトル比(Rcal)と観測スペクトル比(Robs)の差を misfit として式(1)のように定義した。

misfit =1 N ⎡ log Robs ( f i ) − log Rcal ( f i ) ⎤∑⎥ + Tdiff
N i =1 ⎢⎣f i 0.5⎦(1)

この中で Tdiff は、仮定した速度モデルの地中から地表までの S 波の片道走時が観測された 0.41～0.42 秒以外の場合に加算されるペナルティである。計算は乱数の種を変えて 5 回行い、その中の最小の misfit のモデルを最適解とした。次に、この最適な構造を基にして Q 値を推定した。Q 値の同定解析における misfit には式(2)を用いた。第 2 項の W は重み関数で、ある層の Q0 がそのすぐ上の層の Q0 より小さい場合にのみ、その上の層との差を misfit に加算することを示している。

misfit =nl
1 N
log Robs ( f i ) − log Rcal ( f i ) + W ∑ Q0 k − Q0 k −1∑
N i =1K =2(2)GA を用いて TKCH07（豊頃）観測点の地盤モデルを同定した結果を図 6 に示した。この結果は、清水・高
井・笹谷 (2005) 10) による豊頃町中央若葉地区(WKB; TKCH07 を含む)の地震記録および微動探査により求められた速度構造モデル（深さ 400m）の浅部をさらに再解析したものである。図 6 には、ターゲットスペクト
ル比、GA により同定した速度層構造から求めた理論スペクトル比および PS 検層による速度層構造から求めた理論スペクトル比を示している。GA により同定した理論スペクトル比は、ターゲットスペクトル比と比べて、第 4 次ピークまでよく一致している。しかし、5Hz よりも高周波数側では、若干の過大評価となっている。
これは、式(1)の misfit の評価を低周波数側に重みを置いた結果である。一方、PS 検層の速度層構造から求めた理論スペクトル比は、第 4 次ピークが大きくずれていることがわかる。表 1 に同定解析から得られた速度
層構造および Q 値を示した。また、図 7 には同定した速度層構造を示した。
Q 値については、Q=Q0・f n の指数値(n)を 0.7, 0.6, 0.5 に固定して、GA による Q0 の同定を行った。表 1 により、第 1 層～第 4 層までは、n=0.7～0.5 と変化してもそれほど定数値(Q0)に変化がなく、第 5 層と第 6 層の変化が大きいことがわかる。また、n=0.7, 0.6, 0.5 の同定解析の misfit は、n=0.5 の場合が一番小さな値となった。

-147-

図 4 TKCH07（豊頃）観測点の地盤モデル同定解析に用いた 22 地震の震央位置とマグニチュード

Spectral Ratio (Suf./Base)

100

10

1

0.1
0.1

1

10

Frequency (Hz)図 5 TkCH07（豊頃）観測点における 22 地震による S 波スペクトル比とその相乗平均値(●)
Vs(m/sec)

100
0

GA-Model
10

200

400

600

0
PS-ModelPS-Model-20

depth(m)Spectral Ratio (Suf./Base)

Target1

GA-Model

-40

-60

-80

0
0.1

1.0

10.0

Frequency (Hz)-100

図 6 速度層構造の同定結果を基に算出した理論スペクトル比-148-

図 7 同定した速度層構造800

表 1 GA を用いた速度層構造および Q 値の同定結果
観測点 層No.

TKCH07

1
2
3
4
5
6

速度構造[PS検層結果]
Q値
S波速度
層厚
Q0 h0(%)
Q0 h0(%)
nn(m/sec)(m)
51[80]
5.9[2.0]
4.0 12.5
4.9 10.2
102[110]
7.1[2.0]
7.2
6.9
5.7
8.8
158[110]
7.8[10.0]
7.2
6.9
10.3
4.9
0.7
0.6
207[200]
12[24.0] 26.1
1.9
60.7
0.8
347[350]
16[10.0] 140.5
0.4
147.6
0.3
700[530] 51.2[52.0] 188.1
0.3
183.3
0.3
misfit=0.01313
misfit=0.01271
合計層厚
100.0

Q0

h0(%)

n

6.7
7.5
7.2
6.9
13.4
3.7
0.5
32.4
1.5
53.0
0.9
88.1
0.6
misfit=0.01220

これまでの GA による解析結果から、TKCH07（豊頃）観測点の最表層の S 波速度は、Vs=51m/sec であり、従来の PS 検層結果等で観測される S 波速度と比較して、かなり低い速度である。図 3 に示した治水地形分類
図から、本地点は「旧湿地」に分類され、軟弱な腐植土が分布していることが示唆される。これを確認するために、本地点でさらにサウンディング および微動アレー探査（高井ほか (2005) 11) ）を行っている。それ
らの結果から、最表層に非常に軟弱な腐植土が分布していることが確認され、微動アレー観測においては
Vs=60m/sec の S 波速度が深度 5m まで分布する構造が得られている。
４．弱震時の線形地震応答解析
GA を用いた同定解析より得られた S 波速度構造モデルおよび Q 値モデルを用いて弱震動時の線形地震応
答解析を行い、同定結果の有効性を検証する。
線形地震応答解析は、一次元重複反射理論に基づいた応答計算プログラム DYNEQ（吉田・末富 (2003) 12) ）を用いた。線形計算の場合、計算パラメータは、地盤層構造モデルとして S 波速度、層厚、密度および Q 値（または減衰比 h）だけとなる。この中で、S 波速度および層厚は GA による同定結果または PS 検層結果を用い、密度は層ごとに粘土～1.7tf/m3、砂～1.8tf/m3、礫～2.0tf/m3 の値を用いた。また、Q 値については、同定解析で求めた Q=Q0・f n の形の散乱減衰(n=0.5)を用いたものと、福島・翠川 (1994) 13) による Q0 の平均値を採用して、Vs=500m/sec 以上の層で Q0=36.2 (h=1.38%)、Vs=500m/sec 未満の層で Q0=10.8 (h=4.65%)のように周波数に依存しない履歴減衰（Q=Q0：一定）を用いた地震応答解析を行った。
図 8 に Q 値の設定の違いによる線形地震応答解析結果の比較を示した。周波数依存させた散乱減衰を用いた計算結果は、高周波数側（6Hz 以上）の減衰比(h)が小さくなりすぎ、最大加速度が観測値と比較して大きくなっていることがわかる。この他の 4 つの弱震動でも同様な傾向が確かめられたため、今回は周波数に依存しない履歴減衰を使用することとした。この履歴減衰については、原・清田 (1977) 14) による室内動土質試
験より、周波数依存性がないことが確かめられており、周波数依存性を持つ散乱減衰とどのように組み合わせるかは今後の課題と考えられる。検討対象とした弱震動は、TKCH07（豊頃）観測点で記録された地表最大
加速度 50cm/sec2 以上の 7 つの弱震動と、2003 年十勝沖地震の本震記録を用いた。また、線形地震応答解析に
は、KiK-net 地中地震計のトランスバース成分記録の S 波初動部分から約 40 秒間を切り出して使用した。
表 2 に示す弱震動および 2003 年十勝沖地震本震の記録の 8 ケースで線形地震応答解析を行い、GA 同定モ
デルと PS 検層モデルで計算された地表最大加速度について比較した。図 9 に観測地表最大加速度との関係を示したが、弱震動では GA 同定モデルの計算地表最大加速度値が観測値とよく合致している。PS 検層モデルの場合は、計算地表最大加速度値が観測値より小さくなっている。この傾向は、最大加速度値だけなく、フーリエスペクトルなどの周波数領域でも同様であった。一方、本震では両モデルとも計算値が 700cm/sec2 以上と算出されているのに対し、観測値はその約半分の 355cm/sec2 であり強い非線形性が現れている。

-149-

(gal)
120.0(gal*s)1e2

観測地表波形
75.3

1e1
1e0
0.0
1e-1
1e-2
-120.0Acceleration (cm/s2)

(gal)
120.0

1e-3
0.1
(gal*s)1e2

計算地表波形（Q 値散乱減衰）

105.7

1.0

10.0
(Hz)

1e1
1e0
0.0
1e-1
1e-2
1e-3
0.1
(gal*s)1e2-120.0(gal)
120.0

計算地表波形（Q 値履歴減衰）

87.3

1.0

10.0
(Hz)

1e1
1e0

0.0
1e-1
1e-2
1e-3
0.1

-120.0

1.0

10.0
(Hz)

0.0

2.0

4.0

6.0

8.0

10.0

12.0

14.0

16.0

Time (sec)

図 8 Q 値の設定の違いによる線形地震応答解析結果の比較 （2003-0316 M4.6 の弱震動）【GA-Model】【PS-Model】
1000
WeakMotion2

Observed PGA (cm/sec )WeakMotion2

Observed PGA (cm/sec )1000

StrongMotion(MainShock)

100

10

StrongMotion(MainShock)

100

10
10

100

1000

10

Linear Response Analysis PGA (cm/sec2)

100

1000

Linear Response Analysis PGA (cm/sec2)

図 9 GA 同定モデルおよび PS 検層モデルでの線形地震応答解析結果
表2

TKCH07（豊頃）観測点で線形地震応答解析に用いた地震

No.

Earthquake occurrence
date

1
2
3
4

2003-0316
2003-0926（Max. aftershock）
2003-1011
2004-0720

epi- PGA
focal
central obs.
Earthquake occurrence
Mag. depth
No.
dist. (cm/
date
(km)
(km) sec2)
4.6 100
55 75.3 5 2004-1129
7.1
21
123 67.6 6 2004-1206
6.1
28
130 64.4 7 2005-0118
5.0
98
47 51.4 8 2003-0926（Main shock）-150-

epi- PGA
focal
central obs.
Mag. depth
dist. (cm/(km)(km) sec2)
7.1
48
144 125.9
6.9
46
149 105.5
6.4
50
121 53.6
8.0
42
123 354.9

５．強震時の等価線形地震応答解析
強震動時の非線形性を考慮した地震応答解析は、重複反射理論を基に等価線形手法を取り入れた SHAKE(1972) 15) に始まったと考えてよい。本章では、従来の SHAKE による手法、その後提案された等価線形手法の改良手法およびここで提案する手法を用いて、TKCH07（豊頃）で同定した地盤モデルに対して計算を行い、それらの結果を評価した。
5.1 手法SHAKE による手法は、繰返し計算時の最大せん断ひずみと有効ひずみの関係によって特徴づけられる。式(3)に示すように係数αは、周波数によって変化せずα=0.65 の一定値である。また、吉田 (1994) 16) により、
SHAKE の適用性はとりまとめられている。

γ eff = αγ max(3)

ここで、γeff は有効せん断ひずみ、γmax は最大せん断ひずみである。計算の際に用いる動的変形特性曲線は、原位置で採取された不攪乱試料から室内試験を行って求める方法が最も精度が高いと考えられるが、本調査地ではこのような調査・試験は行われていないため、広く用いられている旧土木研究所の標準曲線（岩崎ら(1980) 17) 、横田ら (1982) 18) および岩崎ら(1977) 19) ）を用いた。図 10a と図 10b に旧土木研究所の沖積粘土(Ac)、洪積粘土(Dc)および沖積砂質土(As)の動的変形特性曲線を示した。ただし、Ac および As については平均有効主応力σ’m=98.1kPa の場合を表示した。計算においては Ac および As 層の動的変形特性曲線は、有効上載圧により変化させて用いている。
SHAKE の係数αが一定値であることにより、応答波形の高周波成分の増幅度が過小評価されることは、従来から指摘されてきており、これを改良するために杉戸ら(1994) 20) 、Yoshida et al. (2002) 21) および Kausel &Assimaki (2002) 22) により係数αを周波数依存型にすることが試みられている。まず、杉戸ら(1977)の係数αの定義を式(4)に示す。

γ eff ( ω ) = αF( ω )
γ max
Fmax

(4)

ここで、F(ω)は対象とする層内でのせん断ひずみ波形のフーリエスペクトル、Fmax はその最大値であり、原論文ではα=0.65 としている。この方法は、周波数の高いところで減衰を大きめに評価しすぎることを改良するために考案されたものであるが、せん断ひずみ波形のフーリエスペクトルを直接γeff とγmax を関係づける係数に用いており、せん断ひずみ波形のフーリエスペクトルの凹凸に影響を受けることが予想される。
Yoshida et al. (2002)は、強震記録からせん断ひずみの時刻歴を計算し、ゼロクロス法を用いてせん断ひずみと周波数の関係を求めた。これから、有効せん断ひずみと最大せん断ひずみの関係を式(5)のように設定した。

γ eff = γ max

f < fp⎛ ⎛ log f − log f ⎞ m ⎞p⎟ ⎟ × γ max
γ eff = ⎜ 1 − ⎜⎜⎜log f e − log f p ⎟⎠ ⎟⎝ ⎝⎠f p ≤ f ≤ fe

γ eff = 0fe < f-151-(5)ここで、f は周波数、fp はゼロクロス法で求めた最大せん断ひずみ波形の最大振幅に対応する周波数、fe はそれ以上の周波数で弾性挙動を仮定する境界の周波数で fe=6Hz を提唱している。
以上、 式(3)の係数αに周波数依存性を取り入れる手法について記した。杉戸らの方法は、式(4)からも推察できるように、係数α・F(ω)/Fmax が 1 よりも小さくなり、強震時については観測値より大きい最大加速度を与えると予測されること、また、Kausel & Assimaki (2002) および Yoshida et al. (2002)で示されたせん断ひずみの周波数依存性が、周波数に対してほぼ単調減少の傾向にあること、さらに本研究では液状化しない地点での大ひずみ時の非線形応答の問題を対象としていることから、Yoshida et al. (2002)による式(5)のように不連続な関数で係数αを近似せず、連続した関数で係数αの周波数依存性を近似した。ここでは、Yoshida, Miura &Kobayashi (2000) 23) が提案した式(6)に示す対数型の関係式を用い、せん断ひずみ波形を繰返し計算ごとに算出
し、せん断ひずみ波形のフーリエスペクトル振幅と周波数の関係から、式(6)の係数 A を最小自乗法で求めた。
また係数 m については、m=1.0～3.0 で観測地表波形に最も合致するように設定した。

γ eff = γ max

f < fplog γ eff = A(log f − log f p )m + log γ max

f ≥ fp(6)

なお、式(3)と式(6)より式(3)の係数αは、式(7)で与えられる。

α = 10A(log f −log f p )m

0.3

0.8

Damping ratio h

Shear modulus ratio G/G0

1.0

0.6
0.4
0.2

(7)

Ac: σ'm=98.1kPa
Dc

0.2

Ac: σ'm=98.1kPa
Dc
As: σ'm=98.1kPa

0.1

As: σ'm=98.1kPa

0.0
1.0E-06

0.0

1.0E-05

1.0E-04
Shear strain

1.0E-03

1.0E-06

1.0E-02

図 10a 旧土木研究所の動的変形特性曲線(G/G0)

1.0E-05

1.0E-04
Shear strain

1.0E-03

1.0E-02

図 10b 旧土木研究所の動的変形特性曲線(h)5.2 結果
2003 年十勝沖地震の本震の地中記録を用いた等価線形地震応答解析による計算地表波形と観測地表波形との比較を図 11 に示す。同図で、従来の SHAKE による結果は、高周波成分の増幅が小さく評価されてしまい、
観測波形の最大加速度に比べて約 70%の小さな値となっている。等価線形地震応答解析による最表層の最大
せん断ひずみは 1%程度であり、等価線形解析の有効範囲限界に近い。フーリエスペクトルも高周波数の落ち込みが著しい。一方、杉戸ら (1994) による手法では、最大加速度が過大評価されてしまい、観測値の 2 倍程度となっている。フーリエスペクトルの高周波成分も観測値より大きくなっている。なお、杉戸ら(1994)の手法による地震応答解析は、吉田・末富 (2003)による応答計算プログラム DYNEQ の中のモジュールを用いた。
その際、式(4)のαは原論文で推奨されているα=0.65 に、強震動に最も影響すると考えられる地表から深度13m までの粘土層の減衰の下限値(hmin)については hmin=0.02 としている。
Yoshida et al. (2002) による手法では、
fe=6Hz 以上の高周波数で弾性挙動としているため、非常に高周波成分の多い波形となってしまっている。ま-152-

た、フーリエスペクトルにおいて 6Hz で線形計算のフーリエスペクトルに移行している。提案する手法では、
最大加速度においてほぼ観測値と対応しており、フーリエスペクトルも従来の SHAKE の手法による高周波数の落ち込みもない。また、せん断ひずみと周波数の関係を連続関数で回帰しているため、Yoshida et al. (2002) の

(gal)
700.0

(gal*s)1e4

Obs. ground motion at surface (TR comp.)

1e3
1e2
0.0
1e1
1e0

-354.91e-1
0.1(gal*s)1e4-700.0Acceleration (cm/s2)

(gal)
700.0

Calc. ground motion at surface (SHAKE)

1.0

10.0
(Hz)

1.0

10.0
(Hz)

1.0

10.0
(Hz)

1.0

10.0
(Hz)

1.0

10.0
(Hz)

1e3
1e2
0.0
1e1
-237.81e0-700.0(gal)
700.0

Calc. ground motion at surface (Sugito et al. (1994))

1e-1
0.1
(gal*s)
1e4
1e3
1e2

0.0
1e1
1e0
-700.01e-1
0.1-682.9

0.0

20.0

40.0

Time (sec)Calc. ground motion at surface (Yoshida et al. (2002))

(gal)
700.0

472.7

(gal*s)
1e4
1e3
1e2

0.0

Acceleration (cm/s2)

1e1
1e0

-700.0(gal)
700.0

Calc. ground motion at surface (Proposed method m=2.7)
342.4

1e-1
0.1
(gal*s)
1e4
1e3
1e2

0.0
1e1
1e0
-700.0
0.0

1e-1
0.1
20.0

40.0

Time (sec)

図 11 等価線形地震応答解析結果 波形およびフーリエスペクトル-153-

手法のようにフーリエスペクトルの線形計算への移行が現れない。このような検討結果から、本研究で提案した手法は、現状提案されている手法に比べて有用性があるのではないかと考えられる。
６．強震時の非線形地震応答解析
非線形地震応答解析については、応力－ひずみ関係を双曲線モデルまた R-O モデルなどで定式化して用いるのが一般的である（石原 (1976) 24) ）。しかしながら、今回は等価線形地震応答解析結果と含めて評価を行うため、等価線形地震応答解析の動的変形曲線をそのまま使用できる吉田 (2003) 25) による手法を用いた。
時刻歴解析の非線形地震応答解析では、計算を安定させる目的を含めて、Rayleigh 減衰の設定が必要である。
これについては等価線形地震応答解析の Vs=500m/sec 以上の線形層で用いた h=1.38%の値が、0.5Hz と 6Hz の間で最大 1.38%になるように Rayleigh 減衰の係数を設定した。
図 12 に非線形地震応答解析結果を示した。中段の図が等価線形地震応答解析と共通した動的変形特性曲線を用いた結果であるが、計算最大加速度は、観測最大加速度に比べて約 70%の小さな値となっている。ただし、フーリエスペクトルについては、SHAKE による等価線形計算ほどの高周波数の落ち込みは現れていない。
本解析では、地表応答波形が動的変形特性曲線の採用値に依存することが考えられたため、TKCH07（豊頃）
観測点から南東へ約 1.4km 離れた十勝川堤防左岸堤内側地盤上で実施された第 1 層(Vs=51m/sec)の乱さない試料のサンプリングおよび動的変形特性試験結果（藤井ほか (2005) 26) ）を用いてさらに解析を行った。図 12の下段にその結果を示したが、最大加速度の出現時刻およびフーリエスペクトルの高周波数部について、原位置での動的変形特性曲線を用いた方が、観測値によく適合している。この結果は、原位置での動的変形特性試験結果を解析に取り込むことが重要であることを示している。
(gal)
400.0

Obs. ground motion at surface (TR comp.)(gal*s)
1e4
1e3
1e2

0.0
1e1
1e0
-400.0Acceleration (cm/s2)

(gal)
400.0

Calc. ground motion at surface(Non-linear)-354.91e-1
0.1(gal*s)
1e4

237.6

1.0

10.0
(Hz)

1.0

10.0
(Hz)

1.0

10.0
(Hz)

1e3
1e2

0.0
1e1
1e0

-400.0(gal)
400.0

Calc. ground motion at surface(Non-linear using in-situ
G/G0 & h data)

1e-1
0.1
(gal*s)
1e4

316.9

1e3
1e2
0.0
1e1
1e0
-400.0
0.0

1e-1
0.1
20.0

40.0

Time (sec)

図 12 非線形地震応答解析結果 波形およびフーリエスペクトル

-154-

1000

Linear calc.

10

1

0

0

y = 5067.999 x-5.528
R2 = 0.962

100

strain * sec （*10-4）

strain * sec （*10-4）100

Strain * sec (*10-4)

1000

1000

y = 917.761 x-3.185
R2 = 0.875

SHAKE

10

1

0

10

100

Proposedmethod m=2.7

10

1

0

0

0
1

y = 96.241 x-2.243
R2 = 0.815100

1Frequency (Hz)

10

1

100

10

100

Frequency (Hz)Frequency (Hz)20

4Sugito et al.
(1994)

3

Yoshida et al(2002)

2

Proposedmethod1

SHAKE

Damping ratio h (%)

52

Shear modulus G (MN/m )図 13 第 1 層（深度 2.95m, Vs=51m/sec)での応答計算結果のせん断ひずみ波形のフーリエスペクトル振幅と周波数の関係
Sugito et al.
(1994)

15

Yoshida et al.
(2002)

10

Proposedmethod

5

SHAKE

0

0
0.1

1.0

10.0

0.1

100.0

1.0

10.0

100.0

Frequency (Hz)Frequency (Hz)図 14 第 1 層での剛性率(G)および減衰比(h)の周波数依存性
7．考察
7.1 せん断ひずみ波形のスペクトル
提案した改良型等価線形地震応答解析結果に対して、せん断ひずみ波形の周波数特性について考察する。
まず、地震応答解析を行った全層の中でせん断ひずみの一番大きな最表層（深度 2.95m, Vs=51m/sec）のせん
断ひずみ波形を取り出し、せん断ひずみ波形のフーリエスペクトルを求めた。その周波数に対する傾きを線形計算、等価線形計算(SHAKE)および改良型等価線形計算結果で比較し図 13 に示した。せん断ひずみ波形の周波数特性は、どの解析手法においても周波数に対して 2Hz 付近からほぼ単調減少を示している。また、
SHAKE による傾きが一番大きく、高周波数での増幅度の過小評価と対応している。
7.2 剛性率および減衰比の周波数依存性
式(3)のαに周波数依存性を導入した杉戸ら (1994) の手法、Yoshida et al. (2002) の手法、本研究による提案
手法および SHAKE の計算結果から、剛性率 (G) と減衰比 (h) の周波数依存性を算出し図 14 に示した。杉戸らの手法では、剛性率と減衰比の周波数依存性に凹凸が現れており、結果として高周波数側で小さな減衰比になっていることが、過大な最大加速度が計算される主因だと考えられる。Yoshida et al.の手法では、弾性挙動を示す 6Hz より高周波数側の剛性率および減衰比がフラットになっており、これが応答計算結果の加速度
波形のフーリエスペクトルが線形計算のフーリエスペクトルに移行する原因となっている。提案した手法では、剛性率および減衰比は周波数に対して単調増加および単調減少しており、Yoshida et al. (2002)がゼロクロス法で求めたせん断ひずみ波形の振幅と周波数の関係と整合している。
7.3 係数 m と計算地表最大加速度の関係
式(6)および式(7)における係数 m と設定した m を用いた改良型等価線形計算の地表最大加速度については、
今回の解析の場合は図 15 に示す関係がある。m が増えるにつれて、高周波数が減少する計算となるため地表最大加速度は小さくなる。現時点では、m を事前に決める手法が見つけられていないが、TKCH07 のような

-155-

軟弱地盤の場合、m=1.0～3.0 が適切であると考える。
7.4 提案手法の工学的な利用例
ここまでは工学的基盤から地表までの計算の例を示したが、逆に地表から工学的基盤までの引き戻し計算
は、入力地震動を求めるために必要である。この場合、従来の SHAKE による計算では、高周波側の増幅率を過小評価するため、引き戻し計算では、逆に高周波数側が極端に大きくなる計算結果となり、非現実的な入力地震動となる事例が多くみられた。ここで提案したせん断ひずみの周波数依存性を考慮した等価線形地震
応答解析手法で引き戻し計算を行った例を図16 に示した。
これは、地表でスペクトルフィッテイングした1968
年十勝沖地震の八戸 NS 波を八戸観測点の工学的基盤（野田・上部 (1974) 27) ）に引き戻した計算である。従来の SHAKE による計算では、高周波数側が極端に大きくなる計算結果（図 16 中段）となる。一方、提案した手法（図 16 下段）では入射最大加速度は地表最大加速度の約 1/2.5 で、現実的な数値で引き戻し計算が可能となっている。
Fitting wave to the standard Sa
spectrum at surface

461.1
*

(gal)
500.0

PGA (cm/sec2)

1000
0.0

500
-500.0
0.0

1.0

1.5
2.0
2.5
Coefficient m

3.0

図 15 係数 m と地表最大加速度の関係

8．結論

3.5
Acceleration (cm/sec2)

100
0.5

20.0

40.0

463.6
*

(gal)
500.0

60.0

80.0

Time (sec)

100.0

120.0

140.0

160.0
(sec)

Calc. incident motion at base by SHAKE

0.0-500.0
0.0

20.0

40.0

60.0

80.0

100.0

120.0

140.0

160.0

(sec)
2003 年十勝沖地震の KiK-net 鉛直アレー観
Time (sec)
測記録を用いて地盤の非線形応答の評価を行(gal)
Calc. incident motion at base500.0
by proposed methodった。まず、北海道における KiK-net 全観測
187.9*
記録の地表／地中 S 波スペクトル比を算出した。強震時の非線形性の特徴として、スペク
0.0
トル比の高周波数側の低下およびスペクトル
比卓越周波数の低周波数側への移動の両方が-500.0
現れた観測点として TKCH07（豊頃）
、
0.0
20.0
40.0
60.0
80.0
100.0
120.0
140.0
160.0
TKCH08（大樹）および KSRH10（浜中）を(sec)Time (sec)
抽出し、その中で軟弱層の一番厚い TKCH07
を解析対象とした。
次に TKCH07 観測点のほぼ直下に発生し
図 16 提案手法の工学的な利用例
た 22 地震を用いて、弱震動に対する地表／地中S 波スペクトル比を算出し、その平均値をターゲットとして、遺伝的アルゴリズム (GA) を用いて地盤モデルの同定解析を行った。同定解析で得られた S 波速度構造および Q 値構造から求めた理論スペクトル比は、観測スペクトル比の平均値と良い一致をみた。同定解析した S 波速度構造を用いて弱震動による線形地震応答解析を実施した。この結果から同定した S 波速度構造の妥当性が確認できた。
さらに、2003 年十勝沖地震本震の記録を用いて、従来の SHAKE による等価線形地震応答解析、せん断ひ-156-

ずみの周波数依存性を考慮した等価線形地震応答解析および非線形全応力地震応答解析を行い、観測地表波
形と計算地表波形の比較を行った。SHAKE による等価線形解析は、従来から指摘されていた地表加速度波形の高周波成分を過小評価することを確認し、これを改良するためにせん断ひずみの周波数依存性を考慮した改良型等価線形地震応答解析手法を提案した。
非線形全応力地震応答解析については、TKCH07 観測点近傍の原位置での動的変形特性試験結果を用いて同解析を行った結果、標準的な曲線を用いるより良い結果が得られた。これは、非線形全応力地震応答解析
については、原位置での動的変形特性試験結果を解析に取り込むことが重要であることを示唆している。
今後、工学的基盤以浅の表層地盤においても、今回のような KiK-net 鉛直アレー観測記録を基礎として、同
定解析を用いた地盤構造のチューニングが地震応答解析の予測精度向上に欠かせないと考えられる。また、解析事例を増やして、提案した改良型等価線形地震応答解析のパラメータ m を的確に決定できるよう改善していきたい。特に土質との関係を明らかにすることが重要と考えられる。地盤の非線形応答が地震動に与える影響を検討するためには、地震発生頻度の高い地域で、軟弱地盤に多点の鉛直アレー強震計を設置し、長期間にわたって良質の記録を得ることが必要だと考える。TKCH07（豊頃）観測点はこのような意味で、大変貴重な観測点であり、今後とも継続的な観測が望まれる。
謝 辞本研究では、独立行政法人防災科学研究所の KiK-net による地震記録および地盤構造資料を使用させていただいた。また、非線形地震応答解析で使用した原位置の動的変形特性試験結果は、北海道開発局帯広開発建設部よりデータを提供いただいたものである。東北学院大学の吉田 望教授には、DYNEQ および DYNES3Dの使用法についてご教授いただいた。また、本研究の地震応答解析には、吉田 望・末富岩雄 (2003) による
DYNEQ および吉田 望 (2003) による DYNES3D を使用した。編集委員長の中村 晋博士、特集号担当編集者および 3 人の査読者によるコメントにより原稿が大幅に改善された。記して、ここに感謝の意を表す。なお、本論文は、著者の一人の学位論文の一部であり、第 3 章のサウンディングの詳細については、山本 (2006)
を参照されたい。

参考文献
1) Satoh T., M., Fushimi and Y. Tatsumi: Inversion of strain-dependent nonlinear characteristics of soil using weak andstrong motions observed by borehole sites in Japan, Bull. Seism. Soc. Am., Vol.91, 2001, pp.365-380.
2) Iai S., Morita T., kameoka T., Matsunaga Y. and K. Abiko: Response of a dense sand deposit during 1993 Kushiro-oki
Earthquake, Soils and Foundations, Vol.35-1, 1995, pp.115-131.
3) Higashi, S. and T. Sasatani: Nonlinear site response in Kushiro during 1994 Hokkaido Toho-oki Earthquake, Bull.
Seism. Soc. Am., Vol.90, 2000, pp.1082-1095.
4) 山本明夫：地盤の非線形応答が地震動に与える影響及びその予測に関する研究、北海道大学博士論文、2006
年、pp.17-40, 138 pp.
5) 防災科学技術研究所：http://www.bosai.go.jp/kik/6) 国土地理院：治水地形分類図－茂岩, 技術資料 D1-No.441、1978 年.
7) 山中浩明、石田 寛：遺伝的アルゴリズムによる位相速度の逆解析、日本建築学会構造系論文集、第 46 号、1995 年、pp.9-17.
8) 野口科子、笹谷 努：KiK-net 観測点 TKCH08 における S 波速度構造の評価、北海道大学地球物理学研究
報告、No.67、2004 年、pp.181-195.
9) 前田宜浩、笹谷 努、高井伸雄、清水 学：北海道内の KiK-net 観測点における地表地震計の設置方向の推定、北海道大学地球物理学研究報告、第 68 号、2005 年、pp.141-152.
10) 清水 学、高井伸雄、笹谷 努：十勝川下流域における地下構造と強震動特性（２）浅部 S 波速度構造の推定と増幅特性の検討、日本建築学会大会学術講演梗概集（近畿）、2005 年、pp.267-268.
11) 高井伸雄、清水 学、前田宜浩、新屋雅之、三輪田吾郎、齋藤誠治、山本明夫、笹谷 努：十勝支庁豊頃町における地下構造探査と地震動特性の検討、日本地震工学会・大会 2005 梗概集、2005 年、pp.452-453.

-157-12) 吉田 望、末富岩雄：DYNEQ, A computer program for dynamic response analysis of level ground by equivalentlinear method, version3.23, 2003, pp.1-107.
13) 福島美光、翠川三郎：周波数依存性を考慮した表層地盤の平均的な Q-1 値とそれに基づく地盤増幅率の評
価、日本建築学会論文集、第 460 号、1994 年、pp.37-46.
14) 原 昭夫, 清田芳治：地盤震動解析のための土の動的性質の研究－せん断弾性定数, 減衰定数の周波数依
存性－、第 14 回土質工学研究発表会、1977 年、pp.533-536.
15) Schnabel, P. B., Lysmer, J. and H. B. Seed: SHAKE A Computer program for earthquake response analysis ofhorizontally layered sites, Report No. EERC72-12, University of California Berkeley, 1972.
16) 吉田 望：実用プログラム SHAKE の適用性、軟弱地盤の地震動増幅シンポジウム発表論文集、土質工学会、1994 年、pp.14-31.
17) 岩崎敏男、常田賢一、吉田清一：沖積粘性土の動的変形特性・強度特性について、第 15 回土質工学研究発表会、1980 年、pp.625-628.
18) 横田耕一郎、龍岡文夫：不攪乱洪積粘土のせん断変形係数について、土木学会第 32 回年次学術講演概要
集、第 3 部、1982 年、pp.257-258.
19) 岩崎敏男、龍岡文夫、高木義和：砂のせん断変形係数と減衰の歪依存性について、第 12 回土質工学研究発表会、1977 年、pp.417-420.
20) 杉戸真太、会田尚義、増田民夫：周波数特性を考慮した等価ひずみによる地盤の地震応答解析法に関する一考察、土木学会論文集、No. 493/Ⅲ-27、1994 年、pp.49-58.
21) Yoshida, N., Kobayashi, S., Suetomi, I. and Miura, K.: Equivalent linear method considering frequency dependent
characteristics of stiffness and damping, Soil Dynamics and Earthquake Engineering, vol. 22, 2002, pp.205-222.
22) Kausel, E. and D., Assimaki: Seismic simulation of inelastic soils via frequency-dependent moduli and damping, J.
Eng. Mechanics, vol. 128, January, 2002, pp.34-47.
23) Yoshida, N., Miura, K. and S., Kobayashi: Equivalent linear analysis considering large strains and frequency
dependent characteristics, 12WCEE, 2000, No.1832.
24) 石原研而：土質動力学の基礎、鹿島出版会、1976 年、pp.7-34.
25) 吉田 望: DYNES3D, A computer program for DYNamic response analysis of level ground by Effective Stressnonlinear method, version 2.41, 2003.
26) 藤井紀之、渡辺雅俊、澤田俊一、中田智広：2003 年十勝沖地震による十勝川河川堤防被害の解析的検証、
第 40 回地盤工学研究発表会、平成 17 年度発表講演集、2005 年、pp.2009-2010.
27) 野田節男、上部達生：東京湾沿岸地域の大地震時における地盤加速度、港湾技研資料、No.177、1974 年、
p.12.
（受理：2006 年 3 月 28 日）（掲載決定：2007 年 1 月 29 日）

Nonlinear Soil Response during the 2003 Tokachi-oki Earthquake: Application of KiK-net Vertical Array Records
YAMAMOTO Akio1) and SASATANI Tsutomu2)
1) Technical Manager, Earthquake Geotechnical Engineering Group, Tokyo Technical Center, OYO corp., Dr. Sci.
2) Associate Professor, Faculty of Science, Hokkaido University, Dr. Sci.
ABSTRACTThe nonlinear soil response is examined by using KiK-net vertical array records from a large (the 2003 Tokachi-oki
earthquake; M8) and small (M~6) events. First we search stations which show nonlinear soil response by comparing the
S-wave spectral ratio of the Tokachi-oki earthquake with those of three small events; the S-wave spectral ratio is takenfrom the surface spectrum divided by the borehole one. We find three stations from 106 KiK-net stations in Hokkaido that-158-

show the nonlinear soil response. These stations’ spectral ratios show amplitude decay at high frequencies and peak shift tolower frequency for the spectral ratio from the Tokachi-oki earthquake. The TKCH07 station at Toyokoro in
Tokachi-district, Hokkaido, is chosen to evaluate the nonlinear soil response. This site has very soft ground different fromanother KiK-net stations; S-wave velocity (Vs) is less than 200m/s down to a depth of 30m. This station observed peakground acceleration (PGA) of 403.9cm/sec2 during the Tokachi-oki earthquake. Next we construct the soil model (Vs andQs (quality factor for S-wave) structures) at TKCH07 based on the S-wave spectral ratios from 22 small earthquakes byusing the genetic algorithms (GA) inversion method. Finally we evaluate the nonlinear soil response at TKCH07 bycomparing the surface accelerogram observed during the Tokachi-oki earthquake with those calculated using the borehole
accelerogram and the constructed soil model; the linear response analysis, equivalent linear response analyses andnonlinear response analysis are used in the calculations. The frequency dependent equivalent linear response analysis is
newly proposed; this method well explains the observed surface accelerogram.
Key Words: KiK-net Vertical Array, Non-linearity of the Ground, Inversion Analysis, Genetic Algorithm,
Seismic Response Analysis, The 2003 Tokachi-oki Earthquake-159-