5.2　 詳細法（ハイブリッド合成法）について

　詳細法計算（ハイブリット合成法計算）とは、予め求めておいた２つの計算結果を合成する方法である。ここでは、短周期側を統計的グリーン関数法、長周期側として波数積分法または差分法によりもとめた時刻歴波形を川瀬・松島 (1998)が用いた手法と同様の手法で合成する。合成はある周期（接続周期）を中心とするマッチングフィルターを用い、時刻歴波形をフィルター処理して行う。

　5.2.1　三次元有限差分法

　１） 手法の概要　

　有限差分法は有限要素法とならんで解析領域全体を離散化して解く方法（領域法）の代表的なものである。具体的には、対象領域内に配置された離散化された節点(格子点)において、波動場を記述する波動方程式の変数にテイラー展開を適用し微係数を差分近似することで格子点での値に関する連立一次方程式（差分方程式）を作成し、これを逐次解いていく方法である。
　格子点毎に独立に物理定数が付与できるので不均質性が強い場の問題を解くのに適するほか、非線形の問題にも適用でき、応用範囲が広い。解き得る問題の規模は計算機の記憶容量と計算時間に依存するが、差分法の計算アルゴリズム自体がベクトル化に適しており、他の手法と比べて大規模な数値演算が容易に行える利点がある。また、他の領域型の離散化手法と同様にすべての格子点での応答が同時に得られるので、可視化の技法を使えば地動の面的分布は容易に表現できる。しかし、次に示すような領域法に属する離散化手法が避けて通れない難点がある。第一は、空間上の離散化された格子点すべてに地盤の弾性定数を与えなければならず扱うデータ量が膨大となることである。第二は空間上の離散化によって生じる数値誤差のため計算精度に留意する必要があり、適用できる周波数帯域は格子点間隔から決まる特定の周波数よりも低周波数側に限定されることである。第三は地盤の半無限性（無限境界）の扱いが完全ではないことである。これらに加え、次のような差分法特有の問題点があげられる。
　差分格子は一般的には等間隔の直交格子が基本となるため、格子点の配置に有限要素法のような融通性がなく複雑な境界条件の扱いも容易ではない。例えば、凹凸がある地表面の応力解放条件を満足させるのは容易ではない。時空間上の有限差分法では計算効率が陰解法と比べて格段に優れている陽解法を採用するのが普通である。陽解法とは各時間ステップ毎の各格子点における解をその点自身とその極く近傍の格子点の数ステップ前までの解を用いて、式を解かずに直接計算するだけで求める方法であり、陰解法に比べて計算に必要な記憶容量が非常に少なくて済み、計算時間も大幅に短縮されるという利点がある。しかし、陽解法を採用する場合、空間および時間の離散化において数値安定性を考慮する必要があり、空間的に変化する 値を導入する場合、近似的な扱いとならざるを得ない。上述したようないくつかの難点もあるが、有限要素法と比べて簡便で計算時間が速いので、最近では大規模な領域を対象とした地震波の波動伝播問題によく用いられている。

　２） 地盤および波動伝播特性のモデル化と精度

(1) 格子点の配置

　媒質の離散化における格子点（グリッド）の配置方法として、食い違い格子あるいはスタッガードグリッド（staggered grid）と呼ばれるものを用いる。これは、変位（あるいは速度）を評価する格子点と応力を評価する格子点を格子点間隔 の半分だけずらす方法である（例えば、Virieux and Madariaga, 1982）。変位勾配から歪み・応力が決まるという物理的な性質をよく表しており自然であり、数値安定性が向上する利点がある。ただし、変位（あるいは速度）の各成分の格子点が互いに異なる点にあるため、厳密には同一点での値を求めることができない不自由さがある。

(2) 離散化誤差と数値安定性

　微分演算の差分近似は、微分点の前後での値の変化が十分に小さいことを仮定してテイラー展開を用いて誘導されるので、離散化誤差を小さくするためには格子点間隔を対象波長に対して小さく選ぶ必要がある。格子点間隔を対象波長の何分の一以下にすべきかは差分方程式の近似度（あるいは精度ともいう）によって決まる。微分法の値は近似度４（４次精度）の場合、その周囲の４つの格子点の値から定まる。三次元波動場の計算に最近よく用いられる差分法は空間に関して近似度4のものが多い（例えば、Frankel and Vidale, 1992; Graves, 1996;Pitarka, 1999）。実際の計算上の精度を確保するための格子点間隔は、空間に関して近似度４の場合、対象波長の５分の１以下にする必要がある（Levander, 1988）。それよりも短い波長の波は減衰させられる（例えば、Buell, 1991; 竹中, 1993）と同時に位相速度の変化（数値分散）が顕著になる。波動方程式を陽解法型の差分方程式に近似して解く場合の差分スキームの安定性を確保するための必要条件は、単位の時間ステップ だけ経過したときの波面の進行距離が単位のグリッド内におさまるという条件である。この場合、グリッド間隔 の等間隔グリッドを用いた三次元差分方程式では、最速波の速度を とすると時間ステップ は以下の条件を満足しなければならない（Graves, 1996）。

(5-5)

(3) 吸収境界条件

　計算機の能力に限界があることから、有限差分法でモデル化できる媒質の範囲（計算領域）は限られる。計算領域の境界面に何らの境界条件も与えないとこれらの面は剛体壁（固定端）として振る舞うため、この面に達した波は完全反射して計算領域内に戻ってきてしまう。しかし、実際の地盤にこのような壁は存在しない。この計算上の反射波を防ぐには、十分大きな計算領域を設定し、反射波が着目地点に戻ってくる前に計算を打ち切る必要がある。しかし、この方法では計算する地震波の継続時間が長ければ長いほど計算領域を大きくとらなければならなくなり、記憶容量、計算時間とも膨大になってしまう。そこで、計算領域を必要以上に大きくとらない方法として、境界面に入射した波が吸収される吸収境界条件（absorbing boundary condition）が開発された。この吸収境界条件としてはClayton and Engquist （1977）のものが有名であるが、この吸収境界条件だけでは反射波の抑止が完全ではない。そこで最近では、境界面から20 あるいは30 点の格子点を帯状の波動吸収領域として用い、この中では時間ステップ毎に波の振幅に指数関数を乗じ振幅を徐々に減少させる方法（Cerjan et al., 1985）が用いられている。Graves (1996)、Pitarka (1999)、Aoi and Fujiwara (1999)の有限差分法では、Clayton and Engquist （1977）の吸収境界条件とCerjan et al. (1985)の帯状の吸収領域が併用されている。

(4) 地表面における応力解放条件

　直交格子を基本とする有限差分法で地表面における応力解放条件を満足させ、かつ数値安定性を確保することはそう簡単ではない。有限差分法で地表面を扱う方法の一つは自由地表面での応力零の条件を陽に定式化する方法である（Zero - stress formulation）。地表面が水平な場合には場の逆対称性を利用して応力零の条件を満足するために必要な地表面上および地表面より上の空中の格子点の変位や応力の値を容易に計算できる（例えば、Levander, 1988; Graves, 1996; Pitarka, 1999; Aoi and Fujiwara, 1999）。この場合には数値的な安定性が確保され精度もよい。しかし、この方法を凹凸がある地表面に適用することは難しい。

(5) 値

　媒質の一般的な非弾性減衰を考慮するには応力と歪みとのコンボリューションの操作が必要となるため単純な陽解法を維持できなくなり、計算に必要な記憶容量や計算時間が格段に増大する（例えば、Emmerich and Korn, 1987）。そのため、これまでの有限差分法の計算では 値の効果を無視する場合が多い（例えば、Yomogida and Etgen, 1993; Olsen and Archuleta, 1996）。最近、近似的扱いにより陽解法を維持したまま空間的に変化する 値の効果を有限差分法に導入するテクニックが提案され、その有効性の確認が行われている（Graves, 1996）。このテクニックでは、限られた帯域波の計算を前提に非弾性減衰による位相速度の分散性を無視し、P波とS波の 値を区別せずに周波数に比例する 値を仮定することで、各時間ステップの速度場および応力場に指数関数型の減衰項を乗じる操作を施す。強震記録から推定された 値に関する多くの既往の研究では 値が周波数にほぼ比例する結果となっている。これらの研究で用いた記録の分解能からは約1Hz よりも低周波数帯域で 値が周波数に比例するのかどうか定かではないが、約1Hz より高周波数帯域の観測結果に基づけば、Graves (1996) の周波数比例の 値の仮定は不自然ではない。ただし、この方法ではP波とS波の 値を区別できない。

　３） 震源のモデル化

　断層面上の滑り時間関数や破壊時刻といった運動学的パラメターを先験的に与える運動学的断層モデルの場合、有限差分法のような領域法に属する離散化手法で最も問題なのは、一つの格子点にダブルカップルを厳密に作用させるのが簡単ではないことにある。この問題点の解決法として、分布震源を用いる方法とソースボックス法の二通りがある（竹中, 1993）。前者は有限な広さを持つ複数の格子点にそれぞれ異なった向きのシングルフォースを作用させることで分布震源としてダブルカップルを近似的に表現する方法である（Aboudi, 1971）。この方法は簡便なことからよく用いられてきた（例えば、Frankel, 1993）。最近、Graves (1996)はこの方法と食い違い格子（スタッガードグリッド）との取り合わせがよいことに着目し、三次元場の食い違い格子におけるモーメントテンソルの表現式を導いている。Aoi and Fujiwara (1999)でもGraves (1996)と同様の方法でダブルカップルの点震源をモデル化している。

　４） 不等間隔格子による有限差分法

　本検討では、不連続隔格子を用いたAoi and Fujiwara (1999)の有限差分法により計算を行う。有限差分法による計算を行う際には、計算の対象とする媒体を直方体の格子点に離散化する。この時、計算の安定条件を満たすために、格子間隔は（5-6）式と次の条件を満たす必要がある。

(5-6)

ここで、 は計算したい最大周波数(Hz)、 (km/s) はモデル内の媒質の伝播速度のうち最も遅い速度、 (km) は最小格子間隔である。
　また、用いる格子点数はすなわち必要とする計算機のメモリ量となるため、同じメモリ量を使用できる場合はより大きな領域を計算できる。
　以上の理由により三次元盆地構造を考慮した地震動の計算に不連続間隔格子による有限差分法を用いた。

　５） 三次元有限差分法による計算条件

　以下の条件で計算を実施することとした。

　領域内の最小弾性波速度700m/sと最小グリッド間隔0.1kmを用いると(5-6)式により0.86秒までの計算を精度よく行うことが可能である。このような理由から、後述するハイブリッド法によるマッチングフィルターの中心周期1秒としている。

　5.2.2　統計的グリーン関数法

　壇・佐藤 (1998) による統計的グリーン関数法は、断層面を小断層に分割し、小断層ごとにBoore (1983) の統計的震源モデルを分布させ、Irikura (1986) の方法で モデルに従うように重ね合わせる方法に準拠するが、断層の非一様すべりを破壊モデルを規定する量のうち、低振動数の地震波の放出量に対応する非一様すべり量 と高震動数の地震波の放出量に対応する非一様すべり量 とを考慮した合成方法である。小地震波である半経験的グリーン関数が使用できない場合にも使用できるが、位相をランダムと仮定しており、特に短周期で有効な方法である。

　１）　統計的グリーン関数（加速度フーリエスペクトル）の作成

　断層の食い違い理論(Aki and Richards, 1980)によれば、無限媒質のときの遠方場におけるS波の変位波形 　は、

(5-7)

と表される。ここに、 はS波の放射特性(radiation pattern)、 は媒質の密度、 は媒質のS波速度、 は震源距離、 を断層面上の点、 は断層面上の点 におけるすべり速度時間関数、 は断層面である。本文では、(5-7)式の(5-8)式に相当する部分を震源時間関数部といい、そのフーリエ変換を震源スペクトルと呼ぶ。

(5-8)

　今、すべり速度時間関数 が、断層面上で によらず同形とすると、点 が破壊する時刻を として、(5-8)式は下のように書くことができる。

(5-9)

　以上により地震基盤におけるS波の主要動の統計的グリーン関数をBoore (1983)の統計モデルに準拠して作成した。この統計モデルは、

(5-10)

で表される地震動の加速度フーリエスペクトルのモデルである。
　ここに、 は要素断層に関する添え字で、 は地震動の加速度フーリエスペクトル、 は地震動の放射特性、 および は要素断層における地殻の密度およびせん断波速度、 は地震モーメント、 は臨界振動数、 は高周波遮断振動数、 は定数、 は震源距離、 は地殻の 値、 および は地震基盤の密度およびせん断波速度である。最終項は、自由表面の影響および要素断層における地殻のインピーダンスと地震基盤のインピーダンスとの比較との相違を考慮したものである。
　また、震源の大きさに関する量である および については、以下の式で求める。

(5-11)

ここに は地殻の剛性率、 は実効応力である。

　２） 　時刻歴（統計的経時特性）の作成

　時刻歴の作成には経時特性もしくは位相特性が必要であるが、地震基盤におけるS波の主要動の経時特性に関しては、現在までに研究があまりされていない。このプログラムでは、最終的に統計的グリーン関数を定義する位置を工学的基盤上と考え、気象庁マグニチュード と震源距離 とで規定される佐藤ほか (1994a)による仙台地域の工学的基盤におけるS波の主要動の地震記録から求められた統計的経時特性を準用している。

　経時特性の式は以下の式で表される。

(5-12)

ここで、

である。気象庁マグニチュードは、佐藤 (1989)による式で地震モーメント より(5-13)で計算する。

(5-13)

　３） 要素地震波形（種地震）の作成

　　１） および ２） の結果より要素地震波形（種地震）を作成する。種地震の作成は １） で作成した加速度フーリエスペクトルにランダムな位相を与え(理論的なスペクトルに[-π,π]の一様な乱数で位相を与える)、ランダムな位相を与えた後の処理は以下の通り。

1. 与えたフーリエ振幅と乱数位相で、フーリエ逆変換を行って、時刻歴を作成する。
2. この時刻歴に（時間領域で）統計的経時特性をかける。フーリエ変換をして、位相情報を残す。（振幅の情報は捨てる。）
3. 再度、与えた加速度フーリエ振幅と残しておいた位相で、フーリエ逆変換を行って、時刻歴を作成する。
4. 位相情報を変化させないようにエンベロープ処理を行う。

図5.3にランダムな位相を与えた種地震（要素地震）波形例を示す。

本検討で用いた種地震作成パラメータを表5.1に示す。

　４） 波形合成法（統計的グリーン関数法）アルゴリズム

●壇・佐藤(1998)の方法

　壇・佐藤 (1998)は、断層の非一様すべり破壊モデルを規定する量のうち、低振動数の地震波の放出量に対応する非一様すべり と高振動数の地震波の放出量に対応する非一様すべり速度 とを考慮した合成方法を提案した。壇・佐藤 (1998)の合成方法は、大地震の 番目の要素断層の震源スペクトル および小地震の震源スペクトル を下式のように変更したものである。

(5-14)

ここに、大地震の要素断層の大きさと小地震の大きさは等しいとしている 。また、 は虚数単位、 は小地震の臨界円振動数で、 は媒質のせん断剛性率、 は媒質のS波速度、 は実効応力、 は要素断層を面積が等価な円としたときの半径である。上の大地震の 番目の要素断層の震源スペクトル と小地震の震源スペクトル との比率をとると、スケールファクター は、

(5-15)

となる。壇・佐藤 (1998)の合成方法によって得られる合成波形のフーリエ変換の低振動数領域および高振動数領域における振幅は、距離補正項を無視すると、

(5-16)

となる。ただし、高振動数領域における振幅はランダム和となることを考慮した。ここで非一様な実効応力 が全ての要素断層において に等しいとして、小地震の実効応力 に対する比率を と等しいとして、小地震のすべり量 に対する比率を とおくと、

(5-17)

となる。さらに小地震と大地震との間で巨視的断層パラメータの相似則 が成り立っていると仮定すると、

(5-18)

となる。

　壇・佐藤 (1998)は、この合成方法を、Wald and Somerville (1995)の測地データと周期４秒以上地震記録から同定した。計算の対象としている周期は0.067秒〜4秒であり、やや短周期帯域の地震動を主体にしている。
　壇・佐藤 (1998)は、これらの合成波形より算定した計測震度は、報告されている気象庁の震度と良く対応した値となって強震動の予測問題という観点から、断層の非一様すべり破壊モデルによる結果を、従来の巨視的断層モデルおよび(5-19)で定義される断層面全体の短周期レベルが非一様すべり破壊モデルと等しくなるような等価一様すべり破壊モデルの２つを考え、３つのモデルによる合成結果を相互比較した。

(5-19)

　その結果、計算対象となっている６地点のうち５地点で、非一様すべり破壊モデルによる合成結果と等価一様すべり破壊モデルによる合成結果はほぼ同じとなり、従来の巨視的断層モデルによる合成結果はやや小さくなった。以上より距離補正項を入れたスケールファクター の式を以下に示す。

(5-20)

　5.2.3　ハイブリッド合成法

　ハイブリッド合成法は、短周期領域と長周期領域においてそれぞれ求めておいた2 つの計算結果を合成して広帯域地震動を評価する方法である（例えば、川瀬・松島, 1998 ;入倉・釜江, 1999）。今回は、短周期側では統計的グリーン関数法により求めた = 550m/s層上面での時刻歴波形を用いる。長周期側では三次元有限差分法により求めた = 550m/s層上面での時刻歴波形を用いる。合成は、ある周期（接続周期）を中心とするコサイン型フィルタの組み合わせによるマッチングフィルタを用いて時刻歴波形をフィルタ処理し、時刻歴で重ね合わせて行う。三次元有限差分法により周期0.87秒までの計算を精度よく行うことが可能であることから、ハイブリッド合成法によるマッチングフィルタの中心周期を1.0秒とした。図5.4に、用いたマッチングフィルタを示す。統計的グリーン関数法により計算された波形のロールオフ周期、カットオフ周期はぞれぞれ、0.71秒、1.6 秒であり、三次元有限差分法により計算された波形のロールオフ周期、カットオフ周期はぞれぞれ、1.6 秒、0.71 秒である。 = 550m/s層上面でのハイブリッド合成波を入射波とした。計算波形を表示する10地点では、 = 550m/s層上面でのハイブリッド合成波を「詳細法工学的基盤」での合成波と呼ぶ。６章で図2.1に示す10地点での時刻歴速度波形及び減衰定数5％の疑似速度応答スペクトルを表示する。