水平成層構造における地中の近地理論地震記象

OKAYAMA University Earth Science Reports, Vol.27, No.1, 29–38 (2021)

水平成層構造における地中の近地理論地震記象

Synthetic nearfield seismograms at a subsurface position

in a horizontally layered half-space

竹中博士 (Hiroshi T

)*

渡邉禎貢

(Tomotsugu W

)*

Abstract

We extend the computational code of Takenaka and Sasatani (2000) for synthetic nearfield seismograms for horizontally layered elastic media, based on the reflection/transmission matrices and the discrete wavenumber summation method, to calculate seismic motion and its spatial derivatives at a subsurface position in the attenuative media. In this paper we describe the theory of this extension and show some numerical examples to verify the extended code.

Keywords: synthetic seismogram, reflectivity method, layered half-space

1. はじめに 水平成層構造における地震波動の計算には一般 に reflectivity 法と呼ばれる半解析的な手法がしば しば用いられる（纐纈・竹中, 1989）．reflectivity 法 は ， 地 震 波 動 場 の 支 配 方 程 式 を 時 間 に つ い て Fourier 変換，水平方向の空間座標について 2 次元 のFourier 変換（平面波展開）または Fourier-Bessel 変換（円筒波展開）を施して水平波数－周波数領域 の方程式にする．この方程式は，微分が鉛直座標に ついての微分のみ残っている連立の常微分方程式 である．この常微分方程式は行列を使った各種の 手法（総称して行列法と呼ばれる）で解かれる．行 列法には非常に多くの種類があり，中でも古典的 なpropagator matrix 法（Haskell 法はその一つ）が単 純さゆえに現在でもよく用いられるが，震源を含 む計算では数値的に不安定である．不安定を軽減 するため，例えばDunkin (1965)の delta matrix を P-SV 波の計算に導入することも行われている（例え ば，Kind, 1978; Kohketsu, 1985; 中村・竹中, 2005; Nakamura and Takenaka, 2006; 竹中・中村, 2010）が， propagator matrix に基づく限り，数値的な困難が付 きまとう．この困難を完全に克服した数値的に最

* 岡山大学大学院自然科学研究科，〒700-8530 岡山市北区津島中三丁目１－１

Graduate School of Natural Science and Technology, Okayama University, Okayama 700-8530, Japan.

も安定な手法の一つに，本論文で扱う（propagator matrix ではなく）反射行列・透過行列に基づく解法 がある．そこにも定式化が異なる様々なバリエー ションがあるが，地震学では Kennett and Kerry (1979)流（例えば，Yao and Harkrider, 1983; Sasatani, 1985; 武尾, 1985; Takenaka and Sasatani, 2000）と Luco and Apsel (1983)流（例えば，Chen, 1993; Hisada, 1994, 1995）の二つが主流である．本論文では， Kennett and Kerry (1979)流の解法のひとつである Takenaka and Sasatani (2000)の計算プログラムを以 下のように拡張したので，報告する．まず，媒質に 因果的な非弾性減衰を導入した．次に地表だけで なく地中にも観測点を置けるようにし，さらに地 動（変位・速度・加速度）だけでなくその勾配（空 間微分）も計算できるようにした．

Takenaka and Sasatani (2000)の手法は，Sasatani (1985)の拡張である．Sasatani (1985)は，Aki and Richards (1980)の motion-stress vector と layer matrix に基づいて，Kennett and Kerry (1979)流の反射・透 過行列を定式化した．ただし，方位角展開にAki and Richards (1980)のような複素 Fourier 級数ではなく， 実Fourier 級数を採用していた．Fourier-Bessel 逆変

30 竹中博士・渡邉禎貢

換の波数積分は Bouchon (1981)の離散化波数法 (discrete wavenumber summation method)によって実 行 し て い る ． 震 源 の タ イ プ は 断 層 震 源 と 爆 発 （explosion）震源のみ扱えた．Takenaka and Sasatani (2000)では，方位角展開に Aki and Richards (1980)と 同じ複素フーリエ級数を採用し，震源のタイプを 単一力（シングル・フォース）とモーメント・テン ソルに拡張した．次の節では，今回のさらなる拡張 について説明する． 2. 計算手法 座標系として，Fig. 1 のようなデカルト座標系[𝑋 , 𝑋 , 𝑋 ]の原点を震央とし，𝑋 , 𝑋 , 𝑋 軸をそれぞれ 北向き，東向き，震源を含む鉛直下向きに設定する と，円筒座標r, θ, z はそれぞれ震央距離，方位角， 深さになる．その場合のデカルト座標𝑋 ，𝑋 , 𝑋 を 以降，x, y, z と記す．座標 r, θ に対応する基底ベク トルの向きをそれぞれ radial 方向，transverse 方向 と呼ぶ．このとき，変位場𝑼のデカルト座標系の水 平成分𝑈 （北向き成分）, 𝑈 （東向き成分）と円筒 座標系のRadial 成分𝑈 ，Transverse 成分𝑈 の間には 𝑈 cos 𝜃 𝑈 sin 𝜃 𝑈 , 𝑈 sin 𝜃 𝑈 cos 𝜃 𝑈 , 1 の関係がある． 時間依存を𝑒 とすると円筒座標系の変位成分 は 𝑈 𝑟, 𝜃, 𝑧, 𝜔 1 2π 𝑒 𝑟 𝑧 𝜕𝐽 𝑘𝑟 𝜕 𝑘𝑟 𝑙 𝑧 𝑖𝑚 𝑘𝑟𝐽 𝑘𝑟 𝑘𝑑𝑘, 𝑈 𝑟, 𝜃, 𝑧, 𝜔 1 2π 𝑒 𝑟 𝑧 𝑖𝑚 𝑘𝑟𝐽 𝑘𝑟 𝑙 𝑧 𝜕𝐽 𝑘𝑟 𝜕 𝑘𝑟 𝑘𝑑𝑘, 𝑈 𝑟, 𝜃, 𝑧, 𝜔 1 2π 𝑒 𝑟 𝑧 𝐽 𝑘𝑟 𝑘𝑑𝑘 , 2 とFourier-Bessel 展開表示できる．ここで，𝑖は虚数 単位，𝜔は角周波数，𝑡は時間，𝑘は水平波数，𝑚は 整数，𝐽 ∙ は𝑚次の Bessel 関数である．𝑟 𝑧 , 𝑟 𝑧 , 𝑙 𝑧 は Aki and Richards (1980)で定義され ているmotion-stress vector の変位成分である．その 成分の前二つがP-SV 波，最後が SH 波に当たる． (2)式中の𝑟 𝑧 , 𝑟 𝑧 , 𝑙 𝑧 を motion-stress vector のトラクション成分𝑟 𝑧 , 𝑟 𝑧 , 𝑙 𝑧 でこの順 に置き換えると，水平面に対するトラクションに 対 応 す る 応 力 成 分𝜏 𝑟, 𝜃, 𝑧, 𝜔 , 𝜏 𝑟, 𝜃, 𝑧, 𝜔 , 𝜏 𝑟, 𝜃, 𝑧, 𝜔 の Fourier-Bessel 展開式となる．媒質が 水平成層構造の場合，変位場の方位依存性は震源 に置く力源にのみ依存するので，方位角依存度を 表す整数𝑚も力源の種類によって決まる．単一力 （シングル・フォース）の場合，𝑚 0, 1に限ら Fig. 1. 本論文で採用する座標系．原点を震央とす るデカルト座標系の𝑋 ，𝑋 , 𝑋 軸をそれぞれ北向 き，東向き，鉛直下向きにとると，円筒座標r, θ, z はそれぞれ震央距離，方位角，深さになる．本文で は，その場合のデカルト座標𝑋 ，𝑋 , 𝑋 を x, y, z と 表す．

水平成層構造における地中の近地理論地震記象 31

れ，モーメント・テンソルの場合，𝑚 0, 1, 2 に 限られるので，𝑚についての総和は実際には 2か ら 2までである．表記を簡単にするため，これ以 降特に必要がある場合を除いて，添え字の𝑚 を省 いて記述することにする．Takenaka and Sasatani (2000)ではこの Fourier-Bessel 展開表示式の右辺の 1/ 2π の項が1/ 4π𝜌𝜔 という因子になっていた ことに注意してほしい．ここでの𝜌は震源における 媒質の密度である．この因子は，震源に関わる項で 全成分に共通して現れるので，Takenaka and Sasatani (2000)ではそれを括り出して陽に表示していた．

以下の小節では，Takenaka and Sasatani (2000)の計 算プログラムを拡張するために今回新たに導入し た非弾性減衰，地中の地震波形及び変位勾配の計 算スキームについて説明する． 2.1 非弾性減衰の導入 reflectivity 法に周波数一定の Q 値を用いて因果 律を満たす非弾性減衰を導入する方法は，例えば 武尾(1985)によって提案されている．武尾(1985)は 媒質におけるP 波と S 波の速度を角周波数がゼロ 付近（𝜔 1/𝜏 ; 𝜏 2×104 _{(s)，逆数はほぼ地球の} 自由振動の最小の角振動数に対応）を除いて以下 の分散性を持つ複素数𝑣 𝜔 𝑣 𝜔 𝑣 1 1 𝜋𝑄ln 𝜔 𝜔 𝑖 2𝑄 for 𝜔 ≫ 1 𝜏 , 3

で表し（Aki and Richards, 1980），周波数がゼロ付近 で Liu et al. (1976)の提案する連続的な緩和スペク トルを持つレオロジーモデルに基づいた 𝑣 𝜔 0 𝑣 / 1 1 𝜋𝑄ln 𝜔 𝜏 , 4 に置き換えることで因果律を満たす減衰を導入し た．ここで，𝜔 は𝜔 2π 𝑓 であり，𝑓 は速度構造 モデルの参照周波数 (reference frequency)，𝑣 は入 力として与える実数の速度で周波数𝑓 における P 波またはS 波の位相速度である．なお，武尾(1985) では𝑓 =1 (Hz)に固定しているが，本論文の計算プロ グラムでは入力パラメータとしている． 2.2 水平成層構造中の変位場 震源の深さを𝑧 ，観測点の深さを𝑧 とする．観測 点におけるmotion-stress vector の変位成分ベクトル 𝒓 𝑧 は，観測点が地表にある場合（𝑧 0），反射・ 透過行列を用いて 𝒓 𝑧 0 𝑾 𝑰 𝑹 𝑹 𝑻 𝑰 𝑹 𝑹 𝚺 𝑹 𝚺 , 𝑾 𝑴 𝑴 𝑹, 5 によって計算できる(Kennett, 1983; Sasatani, 1985; Takenaka and Sasatani, 2000)．ここで，各種の反射行 列𝑹と透過行列𝑻の表記の仕方は Kennett (1983)に従 っており（Fig. 2(a)参照），その計算は Sasatani (1985) の式に基づいて行う．𝑴 , 𝑴 は，下降波，上昇 波 を そ れ ぞ れ の 変 位 成 分 に 変 換 す る 行 列 で ， Sasatani (1985)の(11), (13)式（Aki and Richards (2002) では(7.49), (7.55)式）で定義された E 行列の上半分 に位置する二つの小行列である．また，震源に置か れた力源（単一力またはモーメント・テンソル）か ら放射された地震波の震源の深さ𝑧 における下降 波 成 分 𝚺 と上昇波成分 𝚺 の式は Takenaka and Sasatani (2000)にその導出過程とともに記されてい る（ただし，(2)式のところで述べた因子の扱いには 注意）．なお，上昇波成分 𝚺 の符号の定義は， Kennett and Kerry (1979)ではなく Kennett (1983)に基 づいているため，上式では 𝚺 の前の符号が Kennett and Kerry (1979)や Yao and Harkrider (1983), 武尾 (1985)と反対になっているので，注意していただき たい．

観測点が地中にある場合は，その深さが震源よ りも浅い場合（𝑧 𝑧 ）と深い場合（𝑧 𝑧 ）で

32 竹中博士・渡邉禎貢 計算式が異なる．観測点の深さが震源よりも浅い 場合（Fig. 2(b)）， 𝒓 𝑧 𝑧 𝑾 𝑰 𝑹 𝑹 𝑻 𝑰 𝑹 𝑹 𝚺 𝑹 𝚺 , 𝑾 𝑴 𝑴 𝑹 , 6 観測点の深さが震源よりも深い場合（Fig. 2(c)）， 𝒓 𝑧 𝑧 𝑾 𝑰 𝑹 𝑹 𝑻 𝑰 𝑹 𝑹 𝚺 𝑹 𝚺 , 𝑾 𝑴 𝑴 𝑹 , 7 によって観測点におけるmotion-stress vector の変位 成分を求めることができる．次の小節で扱う変位 の鉛直勾配の計算には，観測点における motion-stress vector のトラクション成分も必要になる．ト ラ ク シ ョ ン 成 分 は(5) ～ (7) 式 に お け る𝑾 , 𝑾 , 𝑾 の式の右辺にある𝑴 , 𝑴 , 𝑴 , 𝑴 をそ れ ぞ れ ト ラ ク シ ョ ン に 該 当 す る𝑵 , 𝑵 , 𝑵 , 𝑵 に置き換えるだけで求めることができる． 2.2 変位勾配 歪場などの変形を求める場合には，変位勾配（変 位の空間微分）を計算する必要がある．今回拡張し た計算コードでは，変位の空間微分も計算して出 力できるようにした．本節では，その計算の仕方を 記述する. 2.2.1 水平勾配 変位場の円筒座標系成分の radial 方向の偏微分 係数は，(2)式を震央距離𝑟で偏微分した 𝜕𝑈 𝜕𝑟 1 2π 𝑒 𝑟 𝑧 𝑘 𝜕 𝐽 𝑘𝑟 𝜕 𝑘𝑟 𝑙 𝑧 𝑖𝑚 𝑟 𝜕𝐽 𝑘𝑟 𝜕 𝑘𝑟 1 𝑘𝑟𝐽 𝑘𝑟 𝑘𝑑𝑘, Fig. 2. 震源からの上昇波・下降波の放射と各種の 反射・透過行列の作用を表す模式図．(a) 地表観測 点（𝑧 0）の場合．(b) 震源よりも浅い地中観測 点（𝑧 𝑧 ）の場合．(c) 震源よりも深い地中観測 点（𝑧 𝑧 ）の場合．

水平成層構造における地中の近地理論地震記象 33 𝜕𝑈 𝜕𝑟 1 2π 𝑒 𝑟 𝑧 𝑖𝑚 𝑟 𝜕𝐽 𝑘𝑟 𝜕 𝑘𝑟 1 𝑘𝑟𝐽 𝑘𝑟 𝑙 𝑧 𝑘𝜕 𝐽 𝑘𝑟 𝜕 𝑘𝑟 𝑘𝑑𝑘, 𝜕𝑈 𝜕𝑟 1 2π 𝑒 𝑟 𝑧 𝑘 𝜕𝐽 𝑘𝑟 𝜕 𝑘𝑟 𝑘𝑑𝑘, 8 から求められる．Bessel 関数の微分は，例えば 𝑑𝐽 𝑥 𝑑𝑥 𝐽 𝑥 𝑚 𝑥 𝐽 𝑥 𝑚 𝑥 𝐽 𝑥 𝐽 𝑥 1 2 𝐽 𝑥 𝐽 𝑥 , 9 などの恒等式を用いることにより，もとのBessel 関 数から容易に計算できる． 一方，変位場の円筒座標系成分の方位角方向の 偏微分係数（𝜃による偏微分）は，(2)式の複素指数 関数𝑒 をその𝜃微分𝑖𝑚𝑒 に置き換えるだけで 求まる． デカルト座標系における水平動の水平勾配の成 分は，偏微分の連鎖律による以下の式 𝜕𝑈 𝜕𝑥 𝜕𝑈 𝜕𝑟 𝜕𝑟 𝜕𝑥 𝜕𝑈 𝜕𝜃 𝜕𝜃 𝜕𝑥 , 𝜕𝑈 𝜕𝑥 𝜕𝑈 𝜕𝑟 𝜕𝑟 𝜕𝑥 𝜕𝑈 𝜕𝜃 𝜕𝜃 𝜕𝑥 , 𝜕𝑈 𝜕𝑦 𝜕𝑈 𝜕𝑟 𝜕𝑟 𝜕𝑦 𝜕𝑈 𝜕𝜃 𝜕𝜃 𝜕𝑦 , 𝜕𝑈 𝜕𝑦 𝜕𝑈 𝜕𝑟 𝜕𝑟 𝜕𝑦 𝜕𝑈 𝜕𝜃 𝜕𝜃 𝜕𝑦 , 𝜕𝑟 𝜕𝑥 cos 𝜃, 𝜕𝑟 𝜕𝑦 sin 𝜃, 𝜕𝜃 𝜕𝑥 sin 𝜃 𝑟 , 𝜕𝜃 𝜕𝑦 cos 𝜃 𝑟 , 10 と，(1)式から得られる 𝜕𝑈 𝜕𝑟 cos 𝜃 𝜕𝑈 𝜕𝑟 sin 𝜃 𝜕𝑈 𝜕𝑟 , 𝜕𝑈 𝜕𝑟 sin 𝜃 𝜕𝑈 𝜕𝑟 cos 𝜃 𝜕𝑈 𝜕𝑟 , 𝜕𝑈 𝜕𝜃 sin 𝜃 𝑈 cos 𝜃 𝑈 cos 𝜃𝜕𝑈 𝜕𝜃 sin 𝜃 𝜕𝑈 𝜕𝜃 , 𝜕𝑈 𝜕𝜃 cos 𝜃 𝑈 sin 𝜃 𝑈 sin 𝜃𝜕𝑈 𝜕𝜃 cos 𝜃 𝜕𝑈 𝜕𝜃 , 11 及び上述の円筒座標系成分の radial 方向及び方位 角方向の偏微分係数の式から求めることができる． 2.2.2 鉛直勾配 (1)式から 𝜕𝑈 𝜕𝑧 cos 𝜃 𝜕𝑈 𝜕𝑧 sin 𝜃 𝜕𝑈 𝜕𝑧 , 𝜕𝑈 𝜕𝑧 sin 𝜃 𝜕𝑈 𝜕𝑧 cos 𝜃 𝜕𝑈 𝜕𝑧 , 12 の関係がある．変位場の円筒座標系成分の𝑧方向の 偏微分係数は，(2)式の右辺にある motion-stress vector の変位成分𝑟 𝑧 , 𝑟 𝑧 , 𝑙 𝑧 をその𝑧微分 で置き換えることにより求まる．𝑟 𝑧 , 𝑟 𝑧 , 𝑙 𝑧 の𝑧微分は，motion-stress vector のトラクショ ン成分𝑟 𝑧 , 𝑟 𝑧 , 𝑙 𝑧 を含む

34 竹中博士・渡邉禎貢 𝑑𝑙 𝑧 𝑑𝑧 1 𝜇 𝑧 𝑙 𝑧 , 𝑑𝑟 𝑧 𝑑𝑧 𝑘𝑟 𝑧 1 𝜇 𝑧 𝑟 𝑧 , 𝑑𝑟 𝑧 𝑑𝑧 1 𝜆 𝑧 2𝜇 𝑧 𝑘𝜆 𝑧 𝑟 𝑧 𝑟 𝑧 , 13 の関係から計算することができる．ここで, 𝜆 𝑧 , 𝜇 𝑧 は深さ𝑧における媒質の Lamé 定数である．こ の関係は，motion-stress vector についての同次常微 分方程式（Aki and Richards (2002)の(7.24), (7.28)式; Sasatani (1985)の(6)～(8)式）の一部である． 3. 計算例 本節では今回拡張した計算プログラムを検証す るために行った計算の例（地中観測点の場合）を示 す．理論記象を計算するために用いた地下構造モ デル（Table 1）は，以下のすべての計算例で共通で ある．このモデルは，Sasatani (1985)及び Takenaka and Sasatani (2000)で計算例に使用された速度構造 モデルに非弾性減衰のパラメータ QP,QSを追加し たものである．以下，地中観測点が震源よりも浅い 場合と深い場合に分けて計算例を示す． 3.1 観測点の深さが震源よりも浅い場合 震源を深さ𝑧 10 (km)，走行 70 (°)，傾斜 85 (°)， すべり角 83 (°)，地震モーメント9.46 10 (Nm) の断層震源とし，震源時間関数に Takenaka and Sasatani (2000)と同じ幅 1.0 秒のベル型パルスを用 いる．観測点を震央距離𝑟 3 (km)，方位角𝜃 0 (°)， 深さ𝑧 1 (km)の地点に置いた場合の計算結果を Fig. 3 に示す．検証のため，propagator matrix を用い て地中の波形を計算する竹中・中村 (2010)の計算 プログラムで計算した波形も一緒に示している． 両者ともに波形が完全に一致していることから， 本論文の計算プログラムが正しい波形を計算でき ていることが分かる．ここでは示さないが，より複 雑な構造モデルや他の断層メカニズムや単一力な どの場合についても両プログラムで計算結果を比 較し，すべてのケースにおいて本論文の計算プロ グラムが正しい波形を与えることを確認している． α は P 波速度，β は S 波速度，ρ は密度．reference frequency は 1 Hz. α (km/s) β(km/s) ρ (g/cm3_{) thickness (km) QP QS} 3.5 2.0 2.4 2.0 150 100 6.0 3.5 2.7 ∞ 300 200 Table 1. 地下構造モデル. Fig. 3. 震源よりも浅い地中観測点（𝑧 𝑧 ）の場 合の理論記象（地動速度）．赤色の線は本研究の計 算プログラムによる波形，緑色の破線は竹中・中村 (2010)の計算プログラムによる波形. 両線は完全 に重なっている.

水平成層構造における地中の近地理論地震記象 35

Fig. 4. 震源よりも深い地中観測点（𝑧 𝑧 ）の場合の理論記象（地動速度）の相反定理による 確認．赤色の線は(15)式の左辺のグリーン・テンソル，緑色の破線は(15)式の右辺のグリーン・ テンソル．両線は完全に重なっている. 波形の縦軸に示した振幅目盛の値の単位はモーメント の値が1.0 10 (Nm)のとき, cm/s になる.

36 竹中博士・渡邉禎貢 3.2 観測点の深さが震源よりも深い場合 竹中・中村 (2010)の計算プログラムでは，一旦地 表のmotion-stress vector を求めてそれを地中の観測 点の深さまでpropagator matrix によって移すことに よって地中の波動場を計算するスキームを採用し ている．このスキームでは，震源よりも深い観測点 については数値的な不安定性のために精度の良い 解を得ることが困難である．そこで，ここでは，竹 中・中村 (2010)の計算プログラムを用いる代わり に，以下のように相反定理を適用して観測点と震 源の位置を交換して計算した波形と比較すること で検証する． 単 一 力 に 対 す る グ リ ー ン 関 数 を𝐺 𝒙 ; 𝒙 ， 𝑛, 𝑝 𝑥, 𝑦, 𝑧 と表記する．ここで，グリーン関数の 引数（括弧の中）の最初のベクトルが観測点の位置， 次のベクトルが震源の位置，そして第１添字が変 位の成分，第２添字が力の向きを表わしている．こ のとき，モーメント・テンソルの成分𝑀 ， 𝑞 𝑥, 𝑦, 𝑧 に対するグリーン・テンソルの成分は， 𝐺 , 𝒙 ; 𝒙 と表せる．第３添字の前のコンマは震 源位置𝒙 の第３添字の座標で偏微分することを意 味している．グリーン関数についての相反定理か ら，観測点位置と震源位置を交換したグリーン関 数との関係式 𝐺 𝒙 ; 𝒙 𝐺 𝒙 ; 𝒙 , 14 が成り立つ．第１添字（変位の成分）と第２添字（力 の成分）も両辺で反対になっていることに注意さ れたい．この両辺を𝒙 の𝑞座標で偏微分すると 𝐺 , 𝒙 ; 𝒙 𝐺 , 𝒙 ; 𝒙 , 15 が得られる．𝑧 𝑧 の場合には，左辺の観測点の深 さは震源よりも深いが，右辺では逆になる．また， 左辺のグリーン・テンソルは力源がモーメント・テ ンソルの場合の変位であるが，右辺のグリーン・テ ンソルは力源が単一力の場合の変位勾配である． したがって，本論文の計算プログラムで両辺のグ リーン・テンソルをそれぞれ計算し，両者の波形を 比較して(15)式が成立することを確認することに より，計算プログラムが観測点の深さが震源より も深い場合の変位波形を正しく計算できることの みならず，地中観測点における変位勾配波形も正 しく計算できることも確かめることができると考 えられる． Fig. 4 は, (15)式左辺の震源の深さを𝑧 3 (km), 観測点の位置𝒙 を震央距離𝑟 3 (km)，方位角𝜃 30 (°)，深さ𝑧 10 (km) にして計算した(15)式の 両辺のグリーン・テンソル波形の比較である．この 計算ではモーメント速度関数（震源時間関数）を前 計算例の断層震源と同じベル型パルスを用いてお り，波形は地動速度に当たる．この図から(15)式の 両辺の波形が完全に一致していることが確認でき る．なお，ここでは示さないが，𝑧 と𝑧 の値を逆に したケースや地表震源（𝑧 0）の場合などでも計 算し，同様に(15)式や(14)式の両辺の波形が完全に 一致することを確認している． 4. 結論

Takenaka and Sasatani (2000) が報告している，単 一力またはモーメント・テンソル震源を含む弾性 体の水平成層構造の地表における理論地震記象を 計 算 す る た め の 反 射 ・ 透 過 行 列 に 基 づ い た reflectivity 法の計算プログラムに対して以下の拡 張を施した．まず，媒質に因果的な非弾性減衰を導 入し，次に地表だけでなく地中にも観測点を置け るようにし，さらに地動（変位・速度・加速度）だ けでなくその勾配（空間微分）も計算できるように した．そして，今回拡張した計算プログラムを検証 した際に行った計算の例を示した．この計算プロ グラムは, 地中の地動波形のみならず歪や応力を 計算する際にも有用と考えられる． 謝辞 小松正直博士（岡山大学）には論文原稿を読んで

水平成層構造における地中の近地理論地震記象 37

頂き，有益なコメントを頂きました．一部の作図は Generic Mapping Tools (Wessel and Smith, 1998)を使 用しました．藤原貴生氏（岡山大学）には本稿の体 裁を調整して頂きました．記して感謝致します．

引用文献

Aki, K. and P. G. Richards, 1980, Quantitative Seismology: Theory and Methods, Vol. 1, W. H. Freeman Company, San Francisco.

Aki, K. and P. G. Richards, 2002, Quantitative seismology (2nd edition), University Science Book, Sausalito.

Bouchon, M., 1981, A simple method to calculate Green's functions for elastic layered media, Bulletin of the Seismological Society of America, 71, 959–971. Chen, X., 1993, A systematic and efficient method of

computing normal modes for multilayered half-space, Geophysical Journal International, 115, 391–409. Dunkin, J. W., 1965, Computation of modal solutions in

layered, elastic media at high frequencies, Bulletin of the Seismological Society of America, 55, 335–358. Hisada, Y., 1994, An efficient method for computing

Green's functions for a layered half-space with sources and receivers at close depths, Bulletin of the Seismological Society of America, 84, 1456–1472. Hisada, Y., 1995, An efficient method for computing

Green's functions for a layered half-space with sources and receivers at close depths (Part 2), Bulletin of the Seismological Society of America, 85, 1080– 1093.

Kennett, B. N. L., 1983, Seismic Wave Propagation in Stratified Media, Cambridge University Press, Cambridge.

Kennett, B. L. N. and N. J. Kerry, 1979, Seismic waves in a stratified half space, Geophysical Journal of the Royal Astronomical Society, 57, 557–583.

Kind, R., 1978, The reflectivity method for a buried source, Journal of Geophysics, 44, 603–612.

Kohketsu, K., 1985, The extended reflectivity method for synthetic near-field seismogram, Journal of Physics of the Earth, 33, 121–131.

纐纈一起・竹中博士, 1989, 近地地震波の伝播に関 する理論, 地震（第 2 輯）, 42, 391–403.

Liu, H.–P., D. L. Anderson and H. Kanamori, 1976, Velocity dispersion due to anelasticity; implications for seismology and mantle convection, Geophysical Journal of the Royal Astronomical Society, 47, 41–58. Luco, J. E. and R. J. Apsel, 1983. On the Green’s

functions for a layered half-space, Part I, Bulletin of the Seismological Society of America, 73, 909–929. 中村武史・竹中博士, 2005, 断層破砕帯の異方性が

地震波形に及ぼす影響の検討, 地震（第 2 輯）,

57, 331–342.

Nakamura, T. and H. Takenaka, 2006, A numerical analysis of seismic waves associated with an anisotropic fault zone, Earth, Planets and Space, 58, 569–582.

Sasatani, T., 1985, Seismic wave propagation in a stratified half-space, Journal of the Faculty of Science Hokkaido University, Series VII (Geophysics), 7, 401– 420.

Takenaka, H. and T. Sasatani, 2000, Seismic waves in a horizontally layered half-space with a general point source, Memoirs of the Faculty of Sciences, Kyushu University, Series D, Earth and Planetary Sciences,

31, 19–28. 竹中博士・中村武史, 2010, 横等方的媒質における 理論地震記象への非弾性減衰の導入, 北海道大 学地球物理学研究報告, 73, 185–194. 武尾 実, 1985, 非弾性減衰を考慮した震源近傍で の地震波合成－堆積層での非弾性減衰の効果に ついて－, 気象研究所研究報告, 36, 245–257. Wessel, P. and W. H. F. Smith, 1998, New, improved

version of the Generic Mapping Tools released, EOS Transaction American Geophysical Union, 79, 579.

38 竹中博士・渡邉禎貢

Yao, Z. X. and D. G. Harkrider, 1983, A generalized reflection-transmission coefficient matrix and discrete wavenumber method for synthetic seismograms, Bulletin of the Seismological Society of