課題演習DC「地球の鼓動を探る」:モデリング編
◆◇◆遠地P波波形のモデリング◆◇◆
2009年11月16日〜12月21日(4回)久家慶子担当
(Last Modified: December 3, 2009)

9月おわりの阿蘇火山地震観測実習で、本堂のトンネルに 広帯域地震計を設置しました。実習の後も、その地震計はずっと 地震動をとり続けています。データは、刻々と、 京都のPCに送られています。

今年は、設置してから、世界でたくさんの大地震が起きました。

9月30日早朝(日本時)には、マグニチュード8.1の大地震が 南太平洋サモア諸島で発生しました。 この地震は海の下の浅い深さで起こったために、 大きな津波を発生し、サモア諸島とその周辺の島々に大きな被害をもたらしました。

同日夜(日本時)には、インドネシアのスマトラ島で、 マグニチュード7.5の地震が起きました。 この地震は、深さ約80kmで発生しました。 強震動による建物の倒壊による被害が深刻でした。

そして、10月8日(日本時)には、南太平洋バヌアツ諸島で、 マグニチュード7.7 7.8 7.4の3つの大地震が、 最初の地震から約15分後、約1時間後と、立て続けに起こりました。 日本でも、一時、津波注意報がだされました。

サモア諸島とスマトラ島の2つの大地震については、 本堂トンネルの広帯域地震計の上下動成分の記録は、 データを送るためのケーブル結線が悪く、 残念ながら、正しくとることができませんでした。 このため、P波の波形解析に使えません。 水平動成分に記録されたSH波の波形も、 地震の震源の解析に使うことができるのですが、 残念ながら、これらの2つの大地震については、 日本列島の方向と距離は、SH波の波形解析をするには、 あまり適当ではありませんでした。

一方、スマトラ島の地震のあと、大倉先生がケーブルの結線を 直してくれたおかげで、バヌアツ諸島で起こった大地震では、 上下動成分のデータもきれいにとることができました。 そこで、今回は、この上下動成分のデータを用いて、 バヌアツ諸島で起こった最初の大地震のP波の波形解析を 行います。

この地震は、 10月7日22時3分15.8秒(世界標準時)(日本時間では8日7時3分15.8秒)に 発生しました。 規模は、奥尻島に大きな津波被害をもたらした 1993年北海道南西沖地震(Mw7.7)と同じ程度です。

4回の演習で、そのP波の波形("かたち")をモデリングする(計算機上で 復元する)ことで、 地震の断層やずれの向き、断層のずれの継続時間、震源の深さ を推測・考察します。

●バヌアツ地震の概要と阿蘇・本堂との位置関係

アメリカ地質調査所・地震情報センター(USGS NEIC)によって決められている 地震の震源の位置は、
南緯 13.05° 東経 166.34°
です。

震源があるところでは、インド・オーストラリアプレートが東に向かって 太平洋プレートの下に沈み込んでいます。 バヌアツ地震は、太平洋プレートとインド・オーストラリアプレートの プレート境界の地震と思われます。

本堂の震央距離は、 57° です。 震央距離が 30-100°の範囲を、地震学では「遠地(teleseismic distance)」と 呼んでいます。 この範囲のP波の波形は、 伝播するときに上部マントルや遷移層の複雑な構造の影響を受けにくく、 震源で起こったことを推測するときによく使用します。

地球の内部の速度構造を仮定すると、P波などの地震波が 伝播するにかかる時間を推測することができます。 "iasp91"という世界標準1次元速度モデルを用いると、 P波は本堂に579.2秒かかって到達し、 S波は本堂に1049.4秒かかって到達すると予想されます。 また、あわせて、 本堂に伝播するP波が、震源から鉛直下向きから上方に 31°の角度で射出したP波であることも推測されます。 (この角度を射出角take-off angleと呼んでいます)

地震の震源からみた本堂の方向(方位 azimuth)は、 -36°、本堂からみた地震の震源の方向(逆方位 back-azimuth)は、138°です。 いずれも、北をゼロとして時計回りを正に測られています。

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

●観測波形を準備する:本堂で記録されたP波の変位波形をつくる

地震計で記録されたデータは、地動そのものではありません。 地震計固有の応答関数が含まれたデータになっています。
これらの説明はココ

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

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

サーバへのログインやwin形式のデータを扱うコマンド等は、 第1回目の演習のメモにあります。
ココ

それから、今回は解析には使いませんが、せっかくなので、 サモア諸島とスマトラ島の地震のS波(SH波、SV波)も 見ておきましょう。
作業の詳細はココ

●計算(理論)波形をつくる

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

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

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

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

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

【演習1】構造による時間関数E(t)を作る

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

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

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

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

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

P波が出来たら、SV波放射特性(Rsv)も作成して、同じように確認して。

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

E(t)をgnuplotで図にして、出てくる結果を確認。

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

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

S(t)をgnuplotで図にして、出てくる結果を確認。

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

E(t)*S(t)をgnuplotで図にして、出てくる結果を確認。

【演習3】非弾性減衰の効果P(t)をいれる

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

◇時間に余裕のある人(3回目の演習や4回目開始時にこの項目に達した人)

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

フーリエ変換のプログラムは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
です。

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

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

上のFFTのプログラムでは、周波数領域でi番目の複素データ(前半部)は、
(i-1)/データ長
の周波数に対応します。
(ここで、データ長は、時間ステップ(△t)かける2のn乗です。)

◇時間に余裕のない人

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

●モデリングの課題

(1) 断層の走向、傾斜角、すべり角、ずれの継続時間、深さを変えると、 本堂において予測される計算波形はどのように変化するかを調べる。

(2) その結果をもとにしながら、断層の走向、傾斜角、すべり角、ずれの継続時間、 深さをどのように変えると、計算波形と本堂での観測波形が似るようになるかを調べる。

【12月3日追加】 あわせて、本堂以外の2つの観測点PFO KAPIのP波の波形でも、 パラメタから推測される波形がおおむね説明できることを確認してください。
2つの観測点のデータの情報は ココ

●レポート

モデリングの課題(1)〜(2)についてまとめ、考察を加える。
(全部できないときには、出来たところまででも提出)

2010年1月14日(木)締め切り

Wordファイルかpdfファイルにまとめたものを久家宛にメールで送るか (メールの件名(サブジェクト)を「DCレポート」にすること)、 あるいは理1号館255号室前の封筒へ提出。
レポートには、名前を忘れずにいれること。