6.1.2 計算手法

(1) ハイブリッド合成法

 ハイブリッド合成法は、短周期領域での地震動を統計的グリーン関数法で、長周期領域を3次元差分法で算定し、最終的に両者を重ね合わせて広帯域地震動を作成する方法である。この手法の特徴は、地盤や断層の不均質性の影響を受けにくく波動モデルで説明可能な長周期地震動を物理モデルに基づいて理論的に計算し、逆に地盤や断層の不均質性の影響を受け易く統計的にしか説明できない、もしくは理論的な計算効率が悪い短周期地震動を統計的に計算することである。両領域を接続する周期は、地下構造データや振動特性、差分法計算に用いる計算機の能力を考慮して決定する。短周期成分はハイパスフィルタ処理を、長周期域成分はローパスフィルタ処理を行った後、両者を重ね合わせて接続する(このフィルタをマッチングフィルタと呼ぶ)。両フィルタとも接続周期Tcで0.5となるコサイン型(-90°〜90°)で、Tc/1.4〜Tc/0.6の周期幅で1から0に減衰する。

 想定糸魚川−静岡構造線断層の近くに位置する松本周辺での深部地盤の一次周期は3秒程度である(((財)地震予知総合研究センター、2001))。松本盆地の一次モードによる地震動が3次元差分法で計算可能なように、併せて計算機の性能も考慮し、接続周期を1.5秒とした。図6.1-2に今回用いた接続周期1.5秒のマッチングフィルタを示す。

(2)統計的グリーン関数法

(a) 地震基盤における統計的グリーン関数の作成

 地震基盤における統計的グリーン関数は、佐藤(1994a、1994b)が仙台地域で観測された主に海溝型地震の記録から推定したパラメータを用いたスペクトルモデルと経時特性モデルを用いて、Boore(1983)と同様の手法で作成した。
 地震基盤における S 波主要動のスペクトルは、次式で表わされる地震動の加速度フーリエスペクトルのモデルによって表現した。

(6.1-3)

 ここに、pq は要素断層に関する添字で、Apq(f) は地震動の加速度フーリエスペクトル、F は地震波の放射特性、rpq および bpq は要素断層における地殻の密度および剪断波速度、M0pq は地震モーメント、fcpq は臨界振動数、fmaxpq は高周波遮断振動数、m は定数、rpq は震源距離、Q(f) は地殻の Q 値、rsb および bsb は地震基盤の密度および剪断波速度である。rpqbpq は、震源における密度3.03 g/cm3 とS波速度3.93 km/sである。最終項は,自由表面の影響および要素断層における地殻のインピーダンスと地震基盤のインピーダンスとの相違 (佐藤, 1978) を考慮した項である(壇・他、2000)。
 地殻の Q 値であるQ(f) は、

Q =110f 0.69  f ≧0.8Hz (6.1−4)
Q =110  f <0.8Hz (6.1−5)

で表わされる値とした。(6.1.2−4)式は、佐藤・他(1994b)の推定値であるが、最近の研究では、0.8Hz前後以下でQ値が頭打ちするかあるいは減少することが指摘されていることから、0.8Hz以下で一定値とした。

 要素断層の大きさは、4章で述べたように2km×2kmである。断層モデルの中心点付近で、地震モーメント M0pq を3.28×1023 dyn-cm、断層の長さ Lpq を 2 km 、幅 Wpq を 2 km 、実効応力 σpq を 10MPaとして、統計的グリーン関数を作成した。なお、fcpq は以下の式から算出できる。

(6.1-6)

 更に、強震動を評価する領域内にある 1地点を考えると、その地点は各要素断層から見て常に同じ方位に位置しているわけではないこと、および、短周期の地震記録には明確な放射特性が見られないことから、F は平均的な値を用いることとした。計算地点と断層面との幾何学的関係および断層の滑りのタイプを考慮して、Boore and Boatwright (1984)に示された F の値のうちSH 波 (S 波のうち重力の方向に直交する成分) とSV 波 (S 波のうち SH 波に直交する成分) の値の大きい方を採用した。具体的には、 北部1・2断層では、30°dip dip-slip のSV 波の値である 0.62を、中部1・2断層では、vertical strike-slipのSH波の値0.45を採用した。 fmaxpqm は、佐藤・他 (1994b)の推定値である fmaxpq = 13.5 Hz およびm = 4.2 である。
 一方、時刻歴波形の作成には経時特性あるいは位相特性が必要であるが、地震基盤におけるS波主要動の経時特性に関しては現在までに十分な研究成果が得られていない。しかし、本検討では、最終的に統計的グリーン関数を定義する位置を工学的基盤上とするため、佐藤・他(1994a)の仙台地域の工学的基盤における地震記録から求められた経験的経時特性モデルを準用した。

E(t) = (t- ta) / (tb -ta )2 (tattb) (6.1−7)
E(t) =1 (tbttc) (6.1−8)
E(t) = exp(-ln10 (t -tc )/ (td -tc )) (tcttd ) (6.1−10)

ここに、

log (tb - ta ) = 0.229 MJMApq -1.112 (6.1−11)
log (tb - ta ) = 0.433 MJMApq-1.936 (6.1−12)
log (tb - ta ) = 0.778 logrpq-0.340 (6.1−13)

である。気象庁マグニチュードは、佐藤(1989)に基づいて、地震モーメントM0pqから算出した。

(b)深部地盤上面での統計的グリーン関数の作成

 各計算ポイント直下の三次元地盤モデルから、各計算ポイント毎に1次元地盤モデルを作成し、(a)において作成された地震基盤における統計的グリーン関数を入射波とし、S波の1次元重複反射理論により、3次元地盤構造モデル上面での統計的グリーン関数を計算する。Qs値は、(6.1−5)式、(6.1−5)式と同じである。これは、佐藤ほか(1994b)において、仙台で観測された強震記録のスペクトルインバージョンに基づき求められた経験的地盤増幅率が Qs=110f 0.69を用いた1次元地盤モデルによりほぼ説明できることが示されているためである。

(c)深部地盤上面における統計的グリーン関数を用いた波形合成

 (b)で作成された深部地盤上面における統計的グリーン関数を用いて、壇・佐藤(1998)の断層の非一様滑り破壊を考慮した半経験的波形合成法により波形合成を行なう。この際,大地震の要素断層の滑り量と応力降下量が再現されるように、震源スペクトルを補正した。このようにして算定された要素断層の波形を大地震の断層面全体の破壊過程を考慮して合成した。

(3) 3次元差分法

A. 手法の概要

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

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

・ 格子点の配置

 媒質の離散化における格子点(グリッド)の配置方法として、食い違い格子あるいはスタッガードグリッド(staggered grid)と呼ばれるものを用いる。これは、変位(あるいは速度)を評価する格子点と応力を評価する格子点を格子点間隔hの半分だけずらす方法である(例えば、Virieux and Madariaga、 1982)。変位勾配から歪み・応力が決まるという物理的な性質をよく表しており自然であり、数値安定性が向上する利点がある。ただし、変位(あるいは速度)の各成分の格子点が互いに異なる点にあるため、厳密には同一点での値を求めることができない不自由さがある。図6.1-3 にGraves(1996)の論文に記載されている食い違い格子の模式図を示す。図中、速度の各成分は黒印の点で計算され、応力の各成分は、その黒印の点から半グリッド(h/2)だけずれた白印の点で計算される。また、媒質を規定するパラメータである密度およびラーメの定数(密度、P波速度、S波速度を与えればラーメの定数は一意に定まる)は白丸の点毎に独立に与えられる。

・ 離散化誤差と数値安定性

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

Δt< 0 495 h/ Vmax (6.1-14)

・ 吸収境界条件

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

・ 地表面における応力解放条件

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

・ Q値

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

C. 震源のモデル化

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

D. 不等間隔格子による有限差分法

 本検討では、不等間隔格子を用いたPitarka(1999)の有限差分法により計算を行う。不等間隔格子の模式図を図6.1-4に示す。有限差分法による計算を行う際には、計算の対象とする媒体を直方体の格子点に離散化する。この時、計算の安定条件を満たすために、格子間隔は(6.1.2-10)式と次の条件を満たす必要がある。

f < Vmin /(6h) (6.1-15)

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

E.  想定糸魚川−静岡構造線断層帯地震における三次元差分法の計算条件

 接続周期は1.5秒であるため、マッチングフィルタのカットオフ周期は約1秒となる。従って差分法で解析する最短周期を1秒とする。
 最小グリッド間隔hは最も遅い弾性波速度Vminが1km/sであり、(6.1-14)式をみたすようにh=166mとした。また時間刻みΔtは(6.1-15)式をみたすようにΔt=0.005秒とした。