（日本建築学会構造系論文報告集、No.501、1997/11）

Several techniques for efficiently computing Green's functions and normal mode solutions for layered half-spaces are proposed. First, the reflection/transmision (R/T) matrix method is applied to computing normal mode solutions in order to obtain numerically stable results up to very high frequencies. Second, the phase velocities obtained from above and Filon's quadrature are used for the wavenumber integration of Green's functions; this method greatly reduces the number of total integration points. Finally, the cubic-spline interpolation is used to spacially interpolate Green's function, greatly reducing computational time to obtain many Green's functions. The validity of the methods are checked by comparing results using other methods.

Keywords: Green's function, normal mode solution, layered-half space, R/T matrix method, Filon's quadrature, Green function library, cubic-spline interpolation
グリーン関数、正規モード解、成層地盤、透過反射マトリックス法、フィロンの数値積分法、グリーン関数ライブラリー、スプライン補間法

* 工学院大学建築学科・工学博士　Kogakuin University, Dept. of Architecture, Dr. Eng.

1. はじめに

　成層地盤の正規モード解及び動的グリーン関数を効率的に計算することは、地震学・地震工学のみならず環境振動問題など建築・土木の分野でも重要である。特に近年、地震学・地震工学では震源過程のインバージョンや面的な強震動波形のシミュレーションに効率よく大量にグリーン関数を計算する必要に迫られている（例えば文献 1)）。成層地盤の動的グリーン関数は境界要素法などの数値解析法にも適用されており、その際には精度の高く数値的に安定した解が要求される 2)。
　著者はこれまで、成層地盤におけるグリーン関数を数値的に安定かつ効率的に計算する手法を提案してきた 3),4),5)。そこでは、成層を表記する伝達マトリックスにLuco and Apselの反射／透過マトリックス法（R/T Matrix法）6)を改良し、高振動数までの数値安定性を保証するとともに、振源の深さと観測点の深さが近いときに生ずる数値不安定の問題も、波数積分に解析的な漸近解を導入し解決した。
　本論文では著者の方法を更に発展させ、R/T Matrix法を用いて正規モード解を効率的に計算する方法、グリーン関数の波数積分を高速化する方法、さらに多くの数のグリーン関数をライブラリー値と補間法を組み合わせて計算する方法を提案する。なお、ここで示されるフォートランによるプログラムはインターネットを通じ研究目的に公開されている。

2. 定式

　本論文で使われる成層地盤のグリーン関数の基本的な定式は文献 3)を、振源と観測点の深さが近い場合の解析的な処理法に関しては文献 4)を、参照されたい。以下の２．１ではそれらの要点のみ記し、２．２で正規モード解への適用、２．３と２．４で波数積分法の効率化について、２．５で膨大な数のグリーン関数を計算する諸方法について述べる。

2.1 R/T Matrix法と振源と観測点の深さが近い場合の処理

　成層地盤を表記する伝達マトリックスにはThomson-Haskell法が多用されているが、この方法には高振動数で指数関数が発散し、数値不安定になるという欠点がある。このため、Haskell法を改良または全く別な定式による多くの手法が提案され、応用されている（例えば1),6),7),8),9)）。本論文では発散する指数関数を完全に除去した反射・透過マトリックス法（reflection /transmission matrix、R/T Matrix法）を用いている。その記述法に関してはKennett 7)とLuco and Apsel 6)の２つの方法があるが、ここではLuco and Apsel の方法を採用している。但し、原論文では振源を記述する項に発散する指数関数が残っているため、ここではHarkrider 10)やKennettら 7)が行ったように振源の深さに仮想の層境界を設け、そこにステップ応力条件を導入してその問題を解決している 3)。
　一方、グリーン関数を計算する際、振源と観測点の深さが同じか近い場合は、波数積分の被積分関数が波数とともに発散してしまう。このため何らかの対策が必要になるが、被積分関数の漸近解を導入し収束を早める方法 11),12)、振動する被積分関数の極小・極大値の平均値から近似解を求める方法 13)、積分路を複素積分路に拡大して収束する被積分関数に変換する方法 14)、などが提案されている。それぞれの手法には一長一短があるが（例えば文献 3),5)）、本論文では被積分関数に漸近解を導入する方法を採用している。その中でも多用されている方法は、波数が増大する（振動数が０に近づく）とグリーン関数が振源層と同じ物性の一様全無限体の静的解に漸近する性質を利用するものである。すなわち、その漸近解を被積分関数から引いて、被積分関数の収束を保証してから数値積分を行い、次にその積分値に漸近解の解析積分値を加えて最終的な解を求める方法である 11),12)。但し、この方法では振源が、自由表面や層境界に近い場合、境界上にある場合、または観測点が振源とは違う層にある場合などは、層境界から振幅の大きな反射波や透過波が生じるため、一様全無限体の静的解の収束性は悪くなる。そこで本研究では波数が大きい場合の反射・透過波を含む漸近解を解析的に求め、それを用いている 4)。

2.2 反射・透過マトリックス法による正規モード解の計算

　Kennett法による反射・透過マトリックス法による表面波の分散曲線やモード形、正規モード解などの計算法は文献 15) などに示されているが、ここではChen 16)に習い、R/T法によるHisada 3)の記述による定式を行う。
　初めに表面波の特性方程式を導く。Hisada 3)の(A-1)式を用いると、自由表面で応力０の条件からLove波・Rayleigh波ともに次式を得る。

...............(1)

記号の詳細はHisada 3)を参照されたい。上式の C_d¹ は、第一層内での下方に向かう波動の振幅を表すが、これが存在するためには上式の中カッコ内の固有値が０でなければならない。無次元化を行うと、Love波に対しては

...............(2)

Rayleigh波に対しては

...............(3)

を得る。ここで μ₁ は第一層の剪断弾性係数であり、また、 と はHisada 3)に示されているように、

...............(4)

である。 と は、第一層のＰ波・Ｓ波速度に該当する分岐点で０となるが、(2)、(3)式の分子も同じ分岐点で０となり、これを相殺するために上式には導入している。さらに分子が波数とともに振幅が増大して発散するのを防ぐ目的もある。
　(2)、(3)式が表面波の存在条件を決める特性方程式であり、周期ごとに各モードの固有値（位相速度）が決定される。一方、深さ方向の振幅分布である固有ベクトル、及び群速度は以下のように計算される。まず固有変位ベクトルはHisada 3)の(A-1)式左辺のＤより、固有応力ベクトルは同じくＳより、容易に計算される。一方、群速度は求められた固有ベクトルを深さ方向にエネルギー積分を求め、文献 17)の(7.70)式で計算する。この方法は、各層ごとに固有ベクトルが解析的に積分でき、効率的に群速度が計算される。また微分ではなく積分を用いるため数値的に安定している利点もある。

2.3 グリーン関数の効率的な波数積分法（その１：極の扱い）

　グリーン関数を波数積分して求める際、二つの点に注意しなければならない。まずは被積分関数が表面波の位相速度に相当する波数（極）で発散することである。通常、層構造の速度に虚数を導入して減衰を表すため発散は無くなる。しかし減衰が小さな場合は鋭いピーク値を示すため特別な配慮が必要になる。
　これに対処するために良く用いられるのはPhinneyによる方法 18)である。これは振動数に仮想の虚数 (ω_I) を導入し、被積分関数をなだらかな関数にして波数積分し、時間領域に戻す際にexp(ω_It) を掛けて仮想減衰による減少分をもとに戻す方法である。この方法は非常に簡便であるが、後の計算例に示すように時間領域で exp(ω_It) を掛けた際、数値誤差の蓄積で波形の尾部の精度が悪化する場合があるので注意を要する。
　本研究では、2.2 で計算した表面波の位相速度を利用する。すなわちグリーン関数に先だって位相速度を計算し、各モードの極の位置を調べておき、積分する際に極回りで密な積分点を配置する方法である。この方法では全継続時間で精度良い解が得られる。しかし、構造の減衰が非常に小さい場合には、極回りの積分評価にはCauchyの主値などによる厳密な評価をしないと、精度が悪くなる可能性があることに注意を要する。

2.4 グリーン関数の効率的な波数積分法（その２：波数が大きい場合）

　グリーン関数を波数積分して求める際の二つめの問題は、振源と観測点の水平距離と波数の積が大きい場合、被積分関数が激しく振動し、数値積分が困難になることである。このような関数を効率的に積分する方法が各種提案されている 13),14)。本研究では、被積分関数は調和関数に収束することを利用し、Filonの積分法を採用する。これは、ある小さな積分区間ごとに調和関数の包絡線を高次関数で近似し、その高次関数と調和関数との積を解析的に積分する方法である。必要な積分点は、調和振動する包絡形状のみであるため、通常の積分法に比べて積分点数を著しく減らせる。包絡線を１次近似した方法は文献 9)などによって用いられているが、ここではより精度の高い二次近似式を用いる。
　グリーン関数に必要な波数積分は次式でまとめられる。

...............(5)

ここで、Ｊはベッセル関数、ｋは波数、ｒは振源から観測点までの水平距離である。上式右辺で、kr（波数 S水平距離）が３程度より小さい場合、すなわち第１項の数値積分にはSimpson法を用いる。
　一方、第二項にはFilonの積分を用いる。 k₀×r の値が３程度以上とすると、ベッセル関数は以下の式で非常に精度良く近似される（例えば文献 19)）。

...............(6)

　(6)式を(5)式に代入し、sin項とcos項で整理すると(5)式の第二項は以下の二つの式にまとめられる。

...............(7)

　本研究では被積分関数の収束性はどのような条件下でも2.1 によって保証されているため、積分範囲の上限を k_2n で打ち切っている。
　 k₀ から k_2n までの積分範囲を間隔 Δk で２ｎ個に分割し、小さい方から波数を k₀， k₁， k₂，．．．、 k_2n とする。この時、２次のFilon積分法を用いると(7)式は以下のように表現できる 20)。

...............(8)

ここで、

...............(9)

2.5 膨大な数のグリーン関数を計算するために有効な諸手法

　上に示した諸手法によってグリーン関数は、どのような条件下でも高振動数まで高速に計算される。しかし多重震源モデルの震源インバージョンを行う場合（例えば 21)）、境界要素法の基本解に用いる場合（例えば 2)）、有限要素法や差分法への３次元入射場を計算する場合、などには数千から数百万のグリーン関数の計算が必要となることがある。この場合、以下に示される２つの手法が有効である。
　まず初めは、Hisada 3) (A-1)式で各層の下降・上昇波の振幅を表す係数 C_d^j と C_u^j を、予めライブラリーとしてプログラムの中で計算しておく方法である。実際の波数積分は、このライブラリー値を用いて計算する。この方法は、 C_d^j と C_u^j の計算に最も演算時間をとられること、この係数はある周波数では振源深さhと波数kのみの関数であり、観測点とは独立していること、などから、点振源に対して観測点数がたくさんある場合に有効である。それに対してダブルカップル震源の場合は、Hisada 3) の(22)式に示されるように相反定理を用いていないため、hが観測点深さとなり、逆に震源の数がたくさんある場合に有効となる。
　次の方法はグリーン関数ライブラリーと補間法を組み合わせて用いる方法である（例えば 21)）。すなわち、予めある任意面の荒いグリッド上に、厳密なグリーン関数のライブラリーを作成しておき、任意点のグリーン関数値はライブラリー値を補間して計算する方法である。本研究では補間法としてcubic-spline法 22)を用いている。

3. ソースプログラムの公開

　本論文で示された正規モード解、グリーン関数、補間法などのソースプログラム（フォートラン）及びマニュアル等はインターネットを通じ、研究用に公開されている。anonymous ftpサイトは 133.80.156.38で、ディレクトリーは pub/hisada である。最新情報はそこのREADMEファイルを参照されたい。

4. 計算結果及び考察

　以下の解析では、すべて表１に示す地殻モデルを用いて本諸手法の妥当性の確認を行った。

表１：本研究で用いる地殻構造モデル

密度(g/cm3)	P波速度(m/sec)	Qp	S波速度	Qs	厚さ(m)
2.50	3900.0	100.0	2200.0	50.0	2500.0
2.60	5100.0	200.0	2700.0	100.0	1000.0
2.70	6000.0	300.0	3500.0	150.0	12500.0
2.90	6800.0	500.0	3800.0	250.0	16000.0
3.10	7600.0	500.0	4250.0	250.0	-

4.1　分散曲線と固有モードの計算

　図１に2.2 に示した方法による表面波の位相速度と群速度の分散曲線を示す。図２は、１ HZ の時の特性関数で、○が位相速度に相当する。ここでは特性関数の絶対値（実線）を用いて、その振幅が０になるところを位相速度としている。図に見られるように振動数が大きくなると固有モードの数が増大し、特に各層の剪断波速度より少し大きな所に近傍に集中する。従って位相速度を求めるアルゴリズムには、この部分にはより細かな分割数を用いている。図３には１HZの時の深さ５０ｋｍまでの変位と応力のモード形を示す。低次モードでは表層付近に、高次になるほど深い層まで振幅が分布している。当然ながら、応力は地表面では０となる。

4.2 グリーン関数の計算

　表１の地盤モデルと表２に示す震源モデルと観測点位置を用いて、グリーン関数の計算を行う。表のＸ方向が北、Ｙが東に相当し、振動数は１ HZ まで計算している。まず手法のチェックとして、図４に観測点が地表面にある場合の本手法とSaikiaのプログラムを用いた結果の比較を示す。Saikia9)は Propagator Matrix として Haskell 法を改良した Compound Matrix 法を用い、波数積分法はFrequency - Wavenumber 積分法と１次 Filon 法を組み合わせて用いている。図４に示されるように両者の結果は非常に良く一致している。但し、子細に見ると後続波形部分に少し位相のずれが見られ、またSaikiaの結果には波形の尾部のノイズが見える。これは 2.3 に述べたようにPhinney法 18)に起因すると考えられる。因みに演算時間は、日立のワークステーション3050RXを用いてSaikiaでは波形を得るのに約20秒、本手法では約100秒であり、Saikiaの方が５倍程度早い。すなわち計算時間に関してはPhinney法は積分点を大きく減らせるため有効である。
　後続波の精度、及び観測点深さが震源深さにごく近い場合の精度をチェックするため、Harkrider10)のアルゴリズムで正規モード解を計算して比較する。但し、HarkriderはHaskell 法を用いており、高振動数で解が数値的に発散してしまうため、振動数は0.488 HZ まで落としている。Ｑ値は表１の全層で100として、正規モード解には文献 17)の (7.93)式 を用いる。図５に本手法による結果と正規モード解による結果を示す。このうち(b)と(c)は、観測点深さが震源深さの1000mに比べ、10ｍ浅い場合と深い場合である。どの結果も、実体波に相当する初動部分を除き、両手法による結果は完全に一致している。従って、本手法による結果の妥当性は実証されたと考えられる。

表２：本研究で用いる点震源モデルのパラメータと観測点位置

震源位置(x,y,z:m) 地震moment(dyne*cm) strike角(度) dip角(度) (0,0,1000) 5.0e+23 220 50
rake角(度) rise time(sec) 20 0.4	観測点位置(x,y:m) 観測点深さ(m) (17430.0,-199239.0) 0,990 or 1010

4.3 グリーン関数ライブラリーとSpline補間によるグリーン関数の空間的補間法

　2.5 で述べたように、cubic-spline法を用いてグリーン関数を空間的に補間することを試みる。ここではライブラリーのベースとなるグリーン関数のグリッド間隔が、どの程度の波長間隔であれば、精度良く補間できるかを調べる。
　図６に使用するモデルを示す。傾斜角が90度の2.2 km×2.2 kmの平面に矩形のグリッドを置き、厳密に計算したグリーン関数とspline補間による関数との比較を図７に示す。観測点は平面の対角線上に採る。図７の破線は、最上層内で最小周期１秒の１／８波長に相当するグリッド間隔（275ｍ）で計算した厳密なグリーン関数、実線はspline補間による関数である。図７(a)の実線は、2.2ｋｍ間隔グリッド（すなわち最上層内で最短周期１秒の１波長に相当し、図で△の付いている場所）でグリーン関数を計算し、その間の275ｍ間隔の波形をspline補間した波形である。破線の厳密解と比べると長周期の波は良く再現されているが、短周期波は全く再現されず、最大振幅も過小に評価されている。次に図７(b)の実線は、1.1ｋｍ間隔グリッド（すなわち最小周期１秒の1/2波長に相当し、図で△の付いている場所）でグリーン関数を計算し、その間を補間した波形である。図７(a)に比べ、補間波形は改善されており、実用上はこの程度の荒いグリッド間隔でも支障はなさそうである。但し、子細にみると短周期波形の振幅を過小評価している。最後に図７(c)の実線は、0.55ｋｍ間隔グリッド（すなわち最小周期１秒の1/4波長に相当し、図で△の付いている場所）でグリーン関数を計算し、その間を補間した波形である。破線の厳密な波形と完全に一致していることが分かる。すなわち、cubic-spline補間を用いてグリーン関数を補間して計算する場合、ベースとなるグリーン関数のグリッド間隔は、最小周期の1/4波長で充分である。

5. 結論

　本論文では、成層地盤のグリーン関数を計算するHisadaの方法をさらに発展させ、正規モード解計算への応用、正規モード解及び高次 Filon 積分を利用した波数積分の高速化、大量のグリーン関数を計算するためにcubic-spline法を用いた補間法、などを提案し、その諸手法のチェックを行った。その結果、１）本手法がR/Tマトリックス法を用いているため、分散曲線やモード形なども高振動数まで安定して計算すること、２）観測点深さが震源深さに近い場合にも、精度良いグリーン関数が計算されること、３）最小周期の1/4波長のグリッド間隔でベースとなるグリーン関数を計算しておけば、spline補間によって精度良いグリーン関数が計算されること、などが明らかになった。また本研究で使用したソースコードは、研究用にすべてインターネットにて公開した。

謝辞

　R/Tマトリックス法のグリーン関数と正規モード解への応用は、著者が南カルフォルニア大学に滞在中、X. Chen 氏（現、北京大学）から助言を頂き、プログラム化を行いました。C. Saikia 氏（Woodward Clyde 社）にはグリーン関数計算のプログラム使用を許可して頂き、使用に際しては佐藤俊明氏（清水建設）と供に助言を頂きました。また吉田一博氏（清水建設）及び査読者には貴重なご意見を頂きました。本研究の一部は、米国NSF（ contract Num. ACS-9318163）と Southern California Earthquake Center（NSF Num. EAR - 8920136 and USGS Num. 14-08-0001-A0899）、平成８年度文部省科研費奨励研究、及び工学院大学総合研究所プロジェクト研究によって補助を頂きました。