【2012年10月29日記】
2012年10月28日にカナダ西海岸沖でMw7.7の地震が起きました。
以下にはフィリピンの東方沖の地震について記してありますが、
演習には、カナダ西岸の地震を調べてくれてもよいです。
カナダ西岸の地震の概要と本堂との位置関係は
ココ
8月の阿蘇火山地震観測実習で、本堂のトンネルに 広帯域地震計を設置しました。実習の後も、その地震計はずっと 地震動をとり続けています。
8月31日現地時間20時47分、 フィリピンの東沖合でマグニチュードMw7.6の地震が起きました。 この地震の発生により、太平洋津波警報センター（ホノルル）は 一時、フィリピン、インドネシアなどに津波警報を発令しました。 日本でも、気象庁が、東北から九州にかけての太平洋沿岸と沖縄などに 津波注意報をだしましたが、1日午前0時過ぎに解除しました。
この地震は、 8月31日12時47分33.4秒（世界標準時）（日本時間では31日21時47分33.4秒）に 発生しました。 規模（マグニチュード）は、 奥尻島に大きな津波の被害を及ぼした 1993年北海道南西沖地震（Mw7.7）より少し小さいです。
本堂の広帯域地震計は、この地震からのＰ波を記録しました。 ４回の演習で、そのＰ波の波形（"かたち"）をモデリングする（計算機で つくる）ことで、 地震の断層やずれの向き、断層のずれの継続時間や時間関数、震源の深さ などを推測、考察しましょう。
フィリピンは地震が多く起こっている国で、 とても複雑なテクトニクスをしています。 南部では、東側のフィリピン海溝からフィリピン海プレートが西向きに沈み込む一方、 北部では、西側のマニラ海溝からユーラシアプレートが東向きに沈み込んでいます。 フィリピン諸島の中央には、左横ずれ断層のフィリピン断層が 南北にはしっています。 今回の地震は、フィリピン海プレートが沈み込む、 フィリピン海溝のプレート境界付近で発生したと考えられます。
フィリピンのテクトニクス
（http://volcano.oregonstate.edu/vwdocs/volc_images/southeast_asia/philippines/tectonics.html
から引用）。
本堂の震央距離は、 22.3° です。 震央距離が 30-100°の範囲を、地震学では「遠地（teleseismic distance）」と 呼んでいます。 この範囲のＰ波の波形は、 地球内部を伝播するときに上部マントルやマントル遷移層の 複雑な部分の構造の影響を受けにくいために、 地震の震源域で起こったことを推測するときによく使用します。 今回は、この遠地よりも短い距離で記録がとられているので、 P波の波形が遷移層や他の経路から来る波の影響を受けている 可能性があります。 他の観測点の波形とも比較しながら、少し注意深く使うことにします。
地球の内部の速度構造を仮定すると、Ｐ波などの地震波が 伝播するにかかる時間を推測することができます。 "ak135"という世界標準１次元速度モデルを用いると、 Ｐ波は本堂に295.1秒かかって到達し、 Ｓ波は本堂に539.0秒かかって到達すると予想されます。 また、あわせて、 本堂に伝播するＰ波が、震源から鉛直下向きから上方に 39°の角度で射出したＰ波であることも推測されます。 （この角度を射出角take-off angleと呼んでいます）
また、遷移層での急激な速度変化によって 少し異なった経路を通った複数のP波の到来が予想され、 その走時(震源から観測点までの所要時間） は、ak135モデルから、295.1, 296.7, 297.9秒と予想されています。 更に、 PP波の走時は318.1秒、PPP波の走時は326.0, 387.8, 387.9秒 と予想されています。
なお、地震の震源からみた本堂の方向（方位 azimuth）は、 10°、本堂からみた地震の震源の方向（逆方位 back-azimuth）は、-168°です。 いずれも、北をゼロとして時計回りを正に測られています。
★が地震、太い矢印が波の伝播方向、ihが射出角、φsが方位。
地動を見るためには、取得されたデータから、地震計の応答関数を 取り除く必要があります。 今回は、この作業を、SACというソフトウェアを使って行います。
阿蘇のデータがあるサーバにログインして、
作業を行ってください。
作業の詳細はココ
サーバへのログインやwin形式のデータを扱うコマンド等は、
第１回目の演習のメモにあります。
ココ
方位が0〜20°にある、他の観測点で記録されたＰ波の変位波形
ココ
と比較して、波形のどの特徴が安定してみえているかを
みます。その部分の特徴をモデリングします。
（参考に、変位に変換する前のもとのデータの波形は
ココ）
遠地P波の変位波形（変位の時間関数）は、波線理論をもとに、
で、近似的に計算できます。
ここで、
S(t)は地震の断層運動によって決まる時間関数、
E(t)は地下構造によって生じる種々の波の到着を与える時間関数
（ここでは、直達P波とともに、震源そばの地表での
反射波や変換波を与える時間関数）、
P(t)は震源から観測点まで伝わる時の非弾性の効果を与える時間関数。
演習では、これらを計算するプログラムを作成します。
内容を説明するpdfノートはココ
【1-1】 断層運動から出るＰ波とＳ波の放射特性を計算する
断層の走向、傾斜角、すべり角（すべり方向）、Ｐ波の射出角と方位の５つの 変数を与えたときに、Ｐ波放射特性（Rp）を計算するプログラムをつくる。
プログラムの動作確認として、
(1) 横ずれ断層（走向90°、傾斜角90°、すべり角0°）で射出角が45°のとき、 地震波の伝播する方位を0から360°まで増やしていった場合の Ｐ波放射特性（Rp）をgnuplotで図にしてみる。
(2) 逆断層（走向0°、傾斜角45°、すべり角90°）で 地震波の伝播する方位が90°のとき、 射出角を0から180°まで増やしていった場合のＰ波放射特性（Rp）を gnuplotで図にしてみる。
Ｐ波が出来たら、ＳＶ波放射特性(Rsv)も作成して、同じように確認して。
【1-2】地表での反射波（pP）と変換波（sP）を考慮して、E(t)を作成する
E(t)をgnuplotで図にして、出てくる結果を確認。
【2-1】震源の時間関数S(t)を計算するプログラムをつくり、動作を確認
S(t)をgnuplotで図にして、出てくる結果を確認。
【2-2】出来たS(t)を計算するプログラムを使って、E(t)*S(t)を計算させる
E(t)*S(t)をgnuplotで図にして、出てくる結果を確認。
時間の余裕にあわせて、以下２つから、どちらかを選択して行う。
◇時間に余裕のある人（３回目の演習や４回目開始時にこの項目に達した人）
【ステップ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) その結果をもとにしながら、フィリピン東方沖の地震において、 断層の走向、傾斜角、すべり角、ずれの継続時間、 深さがどのぐらいの値であるかを、 計算波形と観測波形の比較から推測する。 どのパラメタが決まりやすいのか、決まりにくいのかも考える。
2012年11月20日（火）締め切り
Wordファイルかpdfファイルにまとめたものを久家宛にメールで送ってください。
メールの件名（サブジェクト）を「ＤＣレポート」にしてください。
レポート内には、名前を忘れずにいれてください。