Ｃ４−４

第 10 回日本地震工学シンポジウム(1998)

不連続格子を用いた４次精度差分法による波形合成
3-D Fourth-Order Finite-Difference Seismograms Using Discontinuous Grids青井 真 1)、藤原 広行 2)
Shin AOI1) and Hiroyuki FUJIWARA2)
1) 科学技術庁防災科学技術研究所、研究員、博士（理学）National Research Institute for Earth Science and Disaster Prevention, Researcher, Ph. D.
2) 科学技術庁防災科学技術研究所、主任研究官、博士（理学）National Research Institute for Earth Science and Disaster Prevention, Senior Researcher, Ph. D.

要約
大規模な平野において有限差分法による波形合成を行うさいに、地震波速度構造に合わせた不均質な格子を用いるのが非常に有効である。青井・藤原(1997) では食い違い格子による空間・時間２次精度の有限差分法に関して不連続な格子を用いた定式化を行い、十分な計算精度があることを確認した。今回我々は、空間４次・時間２次の有限差分法で同様の不連続格子を用いた場合について定式化を行った。４次精度の有限差分近似を用いることにより、２次精度の有限差分近似を用いるのに比べ格子間隔を 1.5-2 倍程度大きくすることが出来るため、3-8 分の１程度にメモリー量を節約することが出来る。本手法による合成波形は均質格子を用いた有限差分法によるものと１７％以下の誤差で一致しており、十分な精度があるものと考えられる。

SUMMARY
We formulated and checked the accuracy of the finite-difference method (FDM) of second-order andfourth-order approximation for the time and spatial derivatives, respectively, with discontinuous grids where thenumber of grids is changed by eliminating or inserting them.

Discontinuous grids were proved to have sufficient

calculation accuracy concerning the second-order finite-difference operators for space and time (Aoi and Fujiwara,1997).

The use of these discontinuous grids in the FDM with the fourth-order operators for space enabled theapplication of this method to problems of larger scale.

Levander (1988)、Graves (1996)）が一般的に用いられて１．はじめにいる。
地震波のシミュレーションを行うために、運動方程式を有限電子計算機の性能の向上により、今やかなり現実的なモデ差分近似することにより数値的に解く手法は昔から用いられてルを用いて関東平野やロサンジェルスなどの大規模な平野におり（例えば、Boore (1972), Kelly et al. (1976)）、現在ではおける３次元波動場のシミュレーションを行うことが可能になり食い違い格子による定式化（例えば、Virieux (1984, 1986)、つつある（例えば、Sato(1997)）。しかしながら、ごく表層部の- 875 -Ｃ４−４

第 10 回日本地震工学シンポジウム(1998)

低速度の層に関しては波形に大きな影響をおよぼすにも関わxらず、導入することが出来ないのが現状である。また、構造インバージョン（例えば Aoi et al. (1995, 1997)）や震源インバ

yージョンへの応用を考えると、計算機に対する負担の軽い波z形合成法の開発が重要である。
均質な大きさの格子を用いる限り、格子サイズは計算すべき最低波長により決定されるため、低速度の層がごく表層にのvxvyvz

tii
txytxztyz

み存在する場合でも計算領域全体を小さな格子に分割せざるを得ず、大規模な計算を行うのには自ずと限界がある。地震Fig. 1. A unit cell for staggered grids.

波形のシミュレーションを行う場合、地表に近いほど極端に地震波速度が遅い構造を取り扱うことが多いため、地震波速度た方法に十分な精度があることが検証された。

構造に合わせた不均質な格子を用いる必要がある。このよう現在、有限差分法による波形合成のさいには、格子間隔をな不均質な格子を用いることにより、計算機に対する負荷を軽1.5-2 倍程度大きくすることが出来る、空間４次精度の有限差減する試みは２次元問題に関してはすでにいくつかの研究が分近似を用いて計算するのが一般的である（例えば、なされている（例えば、Moczo(1989)、Pitarka(1994)）。

Levander (1988)）。このことにより、メモリー量を３ないし８分不均質格子を用いることによる効果がより顕著である３次元の１に、計算時間も大幅に軽減することが出来る。そこで今回問題に関して、青井・藤原(1997)、Aoi and Fujiwara (1998)

我々は、不連続格子を用いた空間４次・時間２次精度の有限では不連続格子による時間・空間２次精度の有限差分法の定差分法の定式化を行う。

式化が行われ、他の手法との比較により、不連続格子を用いTop of the Rigion IIA
BCDE：Vx：Vy：Vz：τii：τxy：τxz：τyz

A
BCDERegion IIRegion I

Bottom of the Rigion I
Fig. 2. (center)3-D discontinuous grid system. (left) Two transections on the top and at the bottom ofthe overlapping region of Regions I and II, where the elimination or the insertion of grids are necessary.
(right) Two profiles of the discontinuous grids.

- 876 -Ｃ４−４

第 10 回日本地震工学シンポジウム(1998)

２．手法
異なる格子点間隔を持つ二つの食い違い格子（格子点配列はそれぞれ Fig. 1 のとおり）から成る不連続格子（Fig. 2 中央）を用いる。領域 I は細かい格子点間隔（ x , y , z 方向の格子点間隔はそれぞれ ∆x , ∆y , ∆z ）を持つ領域、領域 IIはその３倍の粗い格子点間隔（ x , y , z 方向の格子点間隔はそれぞれ 3∆x , 3∆y , 3∆z ）を持つ領域であり、２つの領域は3∆z / 2 分の領域だけオーバーラップしている。領域 I と領域II の内部における偏微分の項に対しては、４次精度の有限差分近似（例えば、Levander (1988)）、fi =1{ c 0 ( f i +1 / 2 − f i −1 / 2 ) − c1 ( f i +3 / 2 − f i −3 / 2 )}hTable 1: Physical parameters of the structure
VpVs
Density Thickness[km/s] [km/s] [g/cm3][km]
Sediment
2.4
0.8
1.8
1.0Rock4.3
2.5
2.5
∞
堆積層は直径 10km、深さ 1.0km の円柱（中心はｘ、ｙ方向に−50m,0m シフトしている）
Table 2: Source parameters
震源位置 x=-1.2[km] y=-1.2[km] Z=9.4[km]時間関数 Ricker wavelet ( fc = 0.8[Hz] )力の向き single force ( x-direction )（ただし、 c0 = 9 / 8 、 c1 = 1 / 24 ）、を用いている（ h は格子点間隔）。ただし、領域 I と領域 II のＣ面内および領域 II のＤ・Ｅ面内では（Fig. 2 右）、ｚに関する接続領域付近では、ｚに関する偏微分の項に対して４次精度偏微分の項に対しては２次精度の有限差分近似（例えば、の有限差分近似を用いることが出来ない。そこで、領域 I のＢ・Virieux (1984)）、FDM(discontinuous grids)

FDM(homogeneous grids)▼▼▼▼▼▼ Vp=2.4km/s▼ Vs=0.8km/s▼
3▼ ρ =1.8g/cm
▼▼
▼
1km
▼
▼
Region I1.5km▼▼▼ Region II▼▼▼▼
Vp=4.3km/s▼
Vs=2.5km/s▼▼
3ρ =2.5g/cm
▼▼
▼
▼
▼
▼
▼
▼
▼▼
▼
▼▼▼
▼▼▼
▼
▼
▼
▼
▼
▼
▼
▼
▼

X-comp.
0816

24 time[s]single force(-1.2,-1.2,9.4) km

0816

24 time[s]Fig. 3: Velocity waveform of x-component at the observation points on the z-axis ( 0 km ≤ z ≤ 4.9km ,with the interval of 100 m. This interval is 300 m in Region II) in 3-D basin structure, calculated(right) by the FDM using discontinuous grids and (left) by the FDM using uniform grids. On the center,the schematic of the structure employed is shown. The arrow indicates the boundary between the tworegions.

- 877 -Ｃ４−４fi =第 10 回日本地震工学シンポジウム(1998)
３．計算例

1( f i +1 / 2 − f i −1 / 2 )hを用いる。また、４次精度の有限差分近似を用いる場合に一本手法の精度を検証するために、簡単な３次元盆地モデル

般的に行われるように、領域 I の上面、領域 II の底面およびを用いて、不連続格子による計算結果を、通常の均質格子を各領域の側面においても適宜２次精度の有限差分近似を用用いた有限差分法による計算結果と比較する。

いる。

ここで用いるモデルは、直径 20km・深さ 1.0km の円柱状のさらに、領域 I のＤ面内および領域 II の A 面内では（Fig. 2堆積層と岩からなる盆地モデルであり（Fig. 3 中央）、それぞ右）、２次精度の有限差分近似を用いても速度や応力を求めれの媒質は Table 1 に示す物理常数を持つ均質完全弾性体ることは出来ない。これらの値は、他方の領域の値から内挿おである。震源はｘ方向へのシングルフォース（ x = −1.2 km、よび間引くことにより求める。Ｄ面内における領域 I の値は、領y = −1.2 km、 z = 9.4 km）、震 源時間関 数は中心周波数

域 II における値から線形補間（Appendix を参照）を用いて求0.8Hz の Ricker wavelet である （Table 2）。中心周波数めている。また、Ａ面内において領域 II の値を求めるべき格

0.8Hz の Ricker wavelet はおおむね 1.6Hz の周波数帯ま子点は全て領域 I の格子点であるため、領域 I の変数を間引でパワーを持っており、計算領域における最短波長は 600mくことによりその値をそのまま用いることが出来る。注目するべ程度であることから、１波長あたり 6 格子点である。Fig. 3 にきことは、すべての変数の内挿および間引きはそれぞれ１つ本手法（右）および均質な格子を用いた空間４次精度の有限の水平面内でのみ行われており、その内挿以外は通常の有差分法（左）による計算結果を示した。観測点は原点から深さ限差分近似により変数の更新が行われていることである。

方向に 100m 間隔で 40 点（ただし、不連続格子の領域 II に計算領域の境界からの人工的な反射波を最小限に押さえるあっては 300m）である。

ため、Cerjan et al. (1985)による吸収境界条件と、Clayton

Fig. 3 の結果から、両手法でほぼ合成波形が一致しているand Engquist (1977)による無反射境界条件を併用する。

ことが分かり、また、領域 I と領域 II の接続部分を通過することで波形が乱れていないことから、変数の内挿の精度は十分であることが分かる。領域 I と領域 II を接続している部分における領域 II の格子点数は、1.6Hz のＳ波の場合１波長あたり 6格子点程度、Ricker wavelet の中心周波数 0.8Hz の場合１波長あたり 12 格子点程度である。Aoi and Fujiwara (1998)
において、１波長あたりの格子点数と線形内挿の精度に関して考察がなされている。
さらに詳しく誤差を評価するために、

誤差 =両手法の結果の差のＲＭＳ均質格子による結果のＲＭＳを用いて、３つの深さにおける両手法の波形を比較する（Fig.
4）。全ての観測点において、誤差が 17%以下におさまっており、不連続格子による合成波形は非常に精度よく計算出来ることが分かる。
ここで示した例の場合、不連続格子を用いて計算すること
により、メモリー量および計算時間を 4.5 分の１程度に軽減することに成功した。さらに大規模な計算においては、十分の一程度まで計算機に対する負荷を軽減することが出来る。このことはメモリー量を増やすことなく格子サイズを半分以下に出来ることを意味する。

４．まとめFig. ４: The results of computation on z-axis atz = 0.0km, 0.8km and 1.6km by (a) the FDM
using discontinuous grids and (b) the FDM
using uniform grids. The waveforms in thelower part of the figure show the difference
between them.

- 878 -

現実的なモデルを用いた有限差分法による３次元波形シミュレーションを行う上で最も大きな問題は、計算時間およびメモリー要求量であるが、地震波速度構造に合わせた不均質な格子と空間４次精度の有限差分法を組み合わせることにより、より大規模な地震波シミュレーションを行うことの出来る手法を

Ｃ４−４第 10 回日本地震工学シンポジウム(1998)
705-708.

提案した。
細かい格子とその３倍の格子点間隔を持つ粗い格子の二つから成る不連続格子を用い、内挿や間引きにより２つの領域を接続する手法は空間４次精度の有限差分法においても有効であることを示した。通常の均質格子を用いた場合との残差でClayton, R. and B. Engquist, (1977).

Absorbingboundary conditions for acoustic and elastic wave
equations. Bull. Seism. Soc. Am. 67, 1529-1540.
Graves, R. W. (1996).

誤差を評価すると、不連続格子を用いることによる誤差は１propagation７％以下となり、十分な計算精度があることが分かった。４次精

staggered-grid finite difference, Bull. Seism. Soc. Am.

度の有限差分法を用いることにより、２次精度の場合に比べ格子間隔を 1.5-2 倍程度大きくすることが出来るため、メモリー量を３ないし８分の１程度に、計算時間も大幅に軽減することが出来る。また、不連続格子を用いることによりメモリー量および計算時間を数分の１から十分の１程度に軽減することに成功たが、このことはメモリー量を増やすことなく格子サイズを半分

inSimulating seismic wave

3Dmediausing

86, 1091-1106.
Kelly, K. R., R. W. Ward, S. Treitel, and R. M. Alford(1976).

Synthetic seismograms a finite-difference

approach, Geophysics 41, 2-27.
Levandar, A. R. (1988). Fourth-order finite-difference
P-SV seismograms, Geophysics 53, 1425-1436.
Moczo, P. (1989).

以下に出来ることを意味する。

elastic

Finite-difference technique for

SH-waves in 2-D media using irregular grids
application

totheseismicresponse

problem,Geophys. J. Int. 99, 321-329.

謝辞：
Geophysical Institute Slovak Academy of Sciences のPitarka, A. (1994), SH wave propagation in anelastic

Peter Moczo 博士には有益な議論をしていただきました。計

and heterogeneous media using 2D Finite difference

算に科学技術庁防災科学技術研究所のスーパーコンピュー

method with irregular grid, 日本地震学会講演予稿タを使用させていただきました。記して感謝いたします。

集,1994 年秋季大会
Sato, T., R. W. Graves and P. G. Somerville (1997).
Preliminary analysis of long-period basin responsein the Tokyo metropolitan area from 1990 Odwara

参考文献

earthquake (M5.1) using 3D finite difference method,
青井・藤原 (1997)、不連続格子による３次元差分法を用いた波形合成（その１）、日本地震学会 1997 年度秋季大会予地球惑星科学関連学会 1997 年合同大会予稿集、103、
B41-10Virieux,

稿集、C41
Aoi, S., T. Iwata, K. Irikura, and F. J. Sanchez-Sesma(1995). Waveform inversion for determining theboundary shape of the basin structure, Bull. Seism.

Soc. Am., 85, 1445-1455
shape

two-dimensional

waveformbasin

three-component

array(1984).

heterogeneous

inversion

structure
datawithSH-wave

propagation

media:

inVelocity-stressfinite-difference method, Geophysics 49, 1933-1957.
Virieux, J. (1986).
heterogeneous

Aoi, S., T. Iwata, H. Fujiwara and K. Irikura (1997).
BoundaryJ.

P-SV wave propagation inmedia:

Velocity-stressfinite-difference method, Geophysics 51, 889-901.

for
using

Appendix

obliquely

azimuthal plane incident wave, Bull. Seism. Soc.

Am., 87, 222-233

領域 I と領域 II の接続部で必要となる内挿法を具体的に述Aoi, S. and H. Fujiwara (1998). 3-D finite-difference
method using discontinuous grids , Bull. Seism. Soc.

Am., 投稿中
Boore, D. M. (1972).
seismicwave

Finite-difference methods forpropagationinべる。Fig. A-1 に示すように、

a 0 ( x) = 1 − x
a 1 ( x) = x(0 ≤ x ≤ 1)

heterogeneous

で与えられる内挿関数を用いて線形の内挿を行う。これらのmaterials, in Methods in Computational Physics,

内挿関数を用いて、1/3 間隔で内挿するための重みを Table

Vol.11, B. A. Bolt (Editor), Academic Press, New

A-1 に示す。

York.

領域 I の最下面（本文 Fig. 2 の D 面内）において、各変数

Cerjan, C., D. Kosloff, R. Kosloff, and M. Reshef (1985).

は x-y 平面上で内挿を行う必要がある（本文 Fig. 2 左下）。

A nonreflecting boundary condition for discrete

内挿が必要となる各変数の格子点は Fig. A-2 のような配置とacoustic and elastic wave equations, Geophysics 50,なっている。図中の (I,J)、 (i,j) は内挿を行うための局所的- 879 -Ｃ４−４

第 10 回日本地震工学シンポジウム(1998)

Table A-1: Weights for linear interpolation

な番号付けであり、内挿によって得られる値は
110

1

2

30

1/3

2/3

1ai0 = 1 − xiai1 = xi

12/3

1/3

0

0

1/3

2/3

1i1

uij = ∑∑ 1α iI, ,jJU I ,J (i, j = 0,1,2,3; I , J = 0,1)xi

I =0 J =0α iI, ,jJ = aiI ⋅ a Jjとなる。ただし、 U I , J 、 ui, j はそれぞれ領域 I、領域 II の速度または応力を表す。

aI(x)0j I 0J0 01231101/3

2/312xFig. A-1: Interpolation functions for linearinterpolation.

- 880 -3 1Fig. A-2: Grid locations on the plane forinterpolation, where (I, J) and (i, j) are local
numberings for the interpolation.