地震計でとられたデータで使えるものは、 地震波が到達した時刻だけではありません。 地震波の波形（”かたち”）にはたくさんの情報が含まれています。 理論をもとに、観測された地震波の波形を計算機でモデリングする（復元する） ことによって、 地震波が伝わってきた地下構造や、地震波を放射した地震そのものの断層運動など の特性を知ることができます。

今回の４回の演習では、実際に、簡単な理論から Ｐ波の波形（"かたち"）をモデリングする（計算機でつくる）ことで、 地震の断層・ずれの向き、断層のずれの継続時間、震源の深さ などの一部を推測してみましょう。

9月の阿蘇火山地震観測実習で、本堂のトンネルに広帯域地震計を1台設置しました。 今回の演習では、この地震計が取得したデータを用いて、 2019年11月14日UTC（世界標準時）にインドネシアで起きたマグニチュードMw7.1の地震を調べます。

●インドネシアの地震と阿蘇・本堂との位置関係

地震の規模（モーメントマグニチュード）は Mw7.1と推定されています。

本堂までの震央距離は、　31.393° です。 距離は、地球中心を見込む角度で表しています。 震央距離が　30-100°の範囲を、地震学では「遠地（teleseismic distance）」と 呼んでいます。 この範囲のＰ波の波形は、 地球内部を伝播するときに、地震波速度が急変する 上部マントルやマントル遷移層の影響を受けにくいために、 地震の震源域で起こったことを推測するときによく使用します。

地球の内部の速度構造を仮定すると、Ｐ波などの地震波が 伝播するにかかる時間を推測することができます。 "iasp91"と呼ばれる世界標準１次元地震波速度モデルを用いると、 Ｐ波は本堂に377.79秒かかって到達し、 Ｓ波は本堂に683.06秒かかって到達すると予想されます。 あわせて、この本堂に到達するＰ波のray parameterが、8.79 [s/deg]と見積もられ、 この値からP波の射出角（take-off angle）が推測できます。

なお、地震の震央からみた本堂の方向（方位 azimuth）は、 7.5°、本堂からみた地震の震央の方向（逆方位 back-azimuth）は、-171.0°です。 いずれも、北をゼロとして時計回りを正に測られています。

★が地震、太い矢印が波の伝播方向、ihが射出角、φsが方位。

●演習０●　どのようなところで起こった地震かをみる

- どのような場所で起こった地震かをみてください。
- USGSが決めた地震の発生時刻をみてください。
- この地震後に発生した余震の分布を調べてください。
- W-phaseで推定したmoment tensor解の節面の方向をみてください。
※ W-phaseは、P波とS波の間に到達する超長周期の振動です。気象庁のページに簡単な説明があります（こちら）

●演習１●　観測波形をみる
本堂で記録されたＰ波の変位波形をみる

地動を見るためには、取得されたデータから、地震計の応答関数を 取り除く必要があります。 今回は、この作業を、SACというソフトウェアを使って行います。

阿蘇のデータがあるサーバにログインして、 作業を行ってください。
作業の詳細はココ

●演習２●　計算（理論）波形をつくる

遠地P波の変位波形（変位の時間関数）は、波線理論をもとに、

で、近似的に計算できます。

ここで、tは時間で、
S(t)は地震の断層運動によって決まる時間関数、
E(t)は地下構造によって生じる種々の波の到着を与える時間関数 （ここでは、直達P波とともに、震源そばの地表での 反射波や変換波を与える時間関数）、
P(t)は震源から観測点まで伝わる時の非弾性の効果を与える時間関数。

演習では、これらを計算するプログラムを作成します。

内容を説明するpdfノートはココ

【ステップ1】構造による時間関数E(t)を作る

【1-1】 断層運動から出るＰ波とＳ波の放射特性を計算する

断層の走向、傾斜角、すべり角（すべり方向）、Ｐ波の射出角と方位の５つの 変数を与えたときに、Ｐ波放射特性（Rp）を計算するプログラムをつくります。

プログラムの動作確認として、

(1) 横ずれ断層（走向0°、傾斜角90°、すべり角0°）で射出角が30°のとき、 地震波の伝播する方位を0から360°まで増やしていった場合の Ｐ波放射特性（Rp）をgnuplotで図にしてみます。

(2) 逆断層（走向0°、傾斜角45°、すべり角90°）で 地震波の伝播する方位が90°のとき、 射出角を0から180°まで増やしていった場合のＰ波放射特性（Rp）を gnuplotで図にしてみます。

Ｐ波が出来たら、ＳＶ波放射特性(Rsv)も作成して、同じように確認してください。

【1-2】地表での反射波（pP）と変換波（sP）を考慮して、E(t)を作成する

[1-1]のRpとRsvを使って、 深さが与えられたときのE(t)を作成するプログラムを作ります。 できたE(t)をgnuplotで図にして、出てくる結果を確認してください。

【ステップ2】震源の時間関数S(t)を作り、本堂でのE(t)*S(t)を計算する

【2-1】震源の時間関数S(t)を計算するプログラムをつくり、動作を確認

継続時間が与えられたときのS(t)を作成するプログラムを作ります。 S(t)をgnuplotで図にして、出てくる結果を確認します。

【2-2】出来たS(t)を計算するプログラムを使って、E(t)*S(t)を計算させる

E(t)*S(t)を導くプログラムを作ります。 結果をgnuplotで図にして、出てくる結果を確認します。

【ステップ3】非弾性減衰の効果P(t)をいれる

時間の余裕にあわせて、以下２つから、どちらかを選択して行う。

◇時間に余裕のある人

【ステップ2】の出来上がりの波形にフーリエ変換を行って、 周波数領域の中で、非弾性減衰の効果を掛け、 結果を逆フーリエ変換で時間領域に戻すように、 【ステップ2】で作ったプログラムを改良してください。 t*は１秒を使ってください。

フーリエ変換のプログラムはfcoolr.f

メインプログラムから

integer:: n
complex:: f(2048)
real:: rind

call fcoolr( n, f, rind )

のようにして呼び出すことができます。

ここで、配列fに含まれるデータは(2のn乗)個。 フーリエ変換ではrind=-1、逆フーリエ変換ではrind=1となります。 フーリエ変換のとき、変換された結果は、(2のn乗)倍になっていますので、 注意してください。

コンパイルの仕方は
gfortran メインプログラムのファイル名 fcoolr.f
です。

逆フーリエ変換した後、実数部にだけデータがある（虚数部がゼロ）ためには、 周波数領域の複素データが、折り返し点をはさんで、 後半部と前半部で共役になっていることが必要です。 非弾性減衰を前半部の複素データに作用させ、 後半部をその共役な値で埋めることにより、これを実現します。

例えば、以下のような周波数領域で１６の複素データをもつ場合、 ９つ目のデータが折り返し点になります。 この折り返し点をはさんで前後のデータを共役にします。

上のFFTのプログラムでは、周波数領域でi番目の複素データ（前半部）は、
(i-1)/データ長
の周波数に対応します。
（ここで、データ長は、時間ステップ（△t）かける２のｎ乗です。)
折り返し点が、ナイキスト周波数(振動数）になります。

◇時間に余裕のない人

非弾性減衰の応答関数として、デルタ関数に対する t*=1秒の場合にえられる時間関数（時間ステップ△t=0.2秒）を ココ におきました。このデータをプログラムで読み込み、 時間領域で畳み込み積分して、非弾性減衰の効果をいれてください。

●モデリングの課題●

(2) （1）で仮定した値から、 走向、傾斜角、すべり角、ずれの継続時間、深さを変えると、 本堂の理論変位波形がどのように変化するかを調べてください。

(3)（2）の結果と、本堂の観測変位波形との比較から、 走向、傾斜角、すべり角、ずれの継続時間、深さがどのぐらいの値になるか、 推測してください。
どのパラメタが決まりやすいのか、決まりにくいのかも、考察してください。

●レポート●

2020年2月5日（水）締め切り

Wordファイルかpdfファイルにまとめたものを久家宛にメールで送ってください。 メールの件名（サブジェクト）を「ＤＣレポート」にしてください。
レポート内には、名前を忘れずにいれてください。