課題演習 DC（地震）「短周期微動源の推定」（４回）
課題 0．観測について
・ 期間，使用機器，場所など

課題 1．波形を解析してみよう
1. 波形をテキストデータに変換する．
2. 適当なチャネルおよび時間で切り出す．
・ 使うコマンド（例）： wed, dewin, cat, awk, sed
3. プロットしてみる．
・ 使うコマンド（例）： gnuplot など
・ 時系列（単成分，3 成分，複数の観測点）
・ パーティクルモーション（例：EW 成分を横軸，NS 成分を縦軸にしたプロット）
4. 微動の強度を計算してみよう
・ たとえば，30 分間のデータについて，1 分ごとの平均 2 乗振幅を計算してみよう．
・ 読み込んだ１分間 12000 点の波形データを 2 乗して足し合わせるプログラム．
http:// www-seis1.kugi.kyoto-u.ac.jp/dc/18/tremor/tremor.html
・ データを読み込むループと計算するループ
・ OPEN，READ，WRITE，FORMAT，DO，変数（配列）の定義，など
・ 読み込んだ波形データを 12000 点ごとに 2 乗平均するプログラムにする．
・ まずは火口西の自分の担当の地震計の波形から計算してみよう．
・ 微動に変化のありそうな部分を計算してみるとよい．
・ 計算する時間（解析期間，間隔）は変更してもよい．

1

課題 2．短周期微動源のセンブランス解析
1. 阿蘇火山で観測される短周期微動の波形を見つける
2. 波形の切り出し
3. センブランス解析
3.1

センブランス値

3.2

人工的な地震波形記録を作る

3.3

解析プログラムを作る

3.4

実データの解析1. 阿蘇火山で観測される短周期微動の波形を見つける
・ 複数の波形を並べてみて，そのなかからイベントを見つける．
・ なるべく振幅の大きなイベントで，全ての観測点（複数の観測点）で観測されているものが望ましい．

2. 波形の切り出し
・ 火山性微動が含まれている部分の波形を切り出してファイルに書き出す．
・使うコマンド（例）： wed, dewin, cat, paste, awk, sed

・ 切り出す時間幅は自由．ただし余り長くすると、後々の計算時間が長くなるので注意しよう．

3. センブランス解析
3.1 センブランス値
アレイで観測地震波について，時間軸をずらしながら足し合わせることを考える．今回のようなアレイで観測すれば，相関のよい波形が得られる．もし，正しい到来方向と速度の組を用いて時間軸をずらして足し合わせれば，それぞれの山谷どうしを足し合わせることになり，振幅の大きな波形となる．逆に，実際の到来方向と異なる方向から来たと仮定して時間軸をずらして足し合わせると，山谷が合わず振幅の小さな波形となる．
さまざまな方向（と速度の組）から波が到来したと仮定して，山谷の強め具合を評価することで，最適な到来方向と速度の組をさがす．このような解析に適した量として，ここではセンブランス値を用いる．
センブランス値は次の式で定義される．このセンブランス値は波形が同じだと１にな2る。波形が異なると、値が小さくなる。

ここで，u(xi, tk) は観測点（地震計） i における時刻 tk の振幅
・１つの観測点の波形は M 個（サンプリング周波数×時間）の時間成分（データ）・観測点は N 点

xi (d)：基準とする観測点 ( i0 ) と、i 番目の観測点までの距離．
・震源から観測点 i までの距離を Di とすると、xi = Di - Di0

s：スローネス（速度の逆数）
・右辺の分子はずらして足し合わせた波形の二乗振幅和，分母は波形の二乗振幅和の総和になっている．
なお，時間軸をずらすために，観測点間の距離 x が必要となる．地球楕円体に対し２点間の距離 D (m) を求めるために，ヒュベニの距離計算式を用いる（距離を求めるためには他にもたくさんの近似式がある）．

P：２点の平均緯度
ΔP：２点の緯度差
ΔR：２点の経度差

M：子午線曲率半径N：卯酉線曲率半径3.2

人工的な地震波形記録を作る

センブランス解析のプログラムをテストするための人工データをつくる．
サンプルプログラムと基準波形データ（kyo2:/home/seis/kuge/2018/）
・make_data.f90
・ricker.dat
このサンプルプログラムは，震源が火口湯溜まり付近にあったときに，指定した位置にあるアレイでの波形データを計算しファイルに出力する．
・コンパイルして実行すると ex1.dat,..., ex6.dat というファイルができる．
中身は，

3

１列目：時刻
２列目：WES?や SUN?（指定した位置による）での波形データ
gnuplot 等で波形を確認しておこう．
3.3

解析プログラムを作る

・3.１「センブランス値」を参考にして解析プログラムを作成する．
・センブランス値が最大になる場所が，求める震源となる．
・震源の位置の分かっている人工的な記録（3.2 で作成したもの）を使って，実際どれくらいよく震源が決まっているかを確認しよう．
・ある点，あるスローネスを仮定（決め打ち）したときのセンブランス値を計算するプ
ログラムをつくる．
・波形を作るとき使ったプログラムを参考にする．
・震源 (elon, elat) = (131.085245, 32.884672)、スローネス s = 1.0 の時のセンブランス値を求めるプログラムを作り，「3.2 人工的な地震波形記録を作る」で作った波形を解析してみる．
・センブランス値が１になっているかどうか確認する（プログラムが正しければなるはず）．
・(elon, elat) や s をずらしたとき，センブランス値がどう変化するかも確認してみよう．
(波形が説明できない)
0(波形が説明できる)＜－－ [センブランス値] －－＞1・つぎに，(elon, elat)を固定しスローネスを未知数として，センブランス値が最大になるようなスローネスを求めるように変更する．たとえば s を 0.01 刻みに変化させて，その都度センブランス値を計算し，センブランス値が最大となるような s とそのときのセンブランス値を求める．
・人工的な記録（3.2 で作成したもの）を解析して，やはり s = 1.0 の時に最適解が出ることを確認する。
・これができたら任意の (elon, elat) についてセンブランス値とスローネス値の組を計算するようにプログラムを修正する．出力は，「経度」 「緯度」 「センブランス値」 「スローネス」の順にしておくと，後でのお絵かきに便利．
3.3

実データの解析・実際に「砂千里」「火口西」で記録された火山性微動について，上で用いたプログラ
ムをつかって実際の微動源を決定する．
・それぞれのアレイ毎に解析を行う

4

・２つのアレイで，同じ時間帯に得られた微動の波形を使う．
・解析する時間が長いと、計算時間が掛かるので注意しよう．
・センブランスの最大値を与える場所が，震源と推定される．
・その時のスローネスも求めよう．逆数をとることで，地震波の速度がわかる．

5

■ 計算機の使用方法などに関するメモ ■※セキュリティ上、ホームページ版では IP アドレスは伏せてあります
データ等の場所
サーバー kyo2 (XXX.XX.XX.X)に置いてある．
ディレクトリ構造/home/seis/

データ等の共通のファイルの置き場/home/seis/arraydata

アレイ観測のデータ（win フォーマット．観測点，日付ごと）/home/seis/cmgdata

広帯域地震計のデータ（〃）

波形画像
アレイ観測のデータを 30 分ごとの画像にして，以下にある．観測点ごとになっている．
http:// www-seis1.kugi.kyoto-u.ac.jp/dc/18/tremor/tremor.html

波形の振幅を物理値になおす方法
チャネル表に記載されている値を使う（単位は[9 列目] ）= データ x [13 列目] ÷ [8 列目] ÷ 10[12 列目]/20
パラメータファイル
kyo2:/home/seis/kuge/2018/dc18.tbl
このファイルには観測点のチャンネル番号[１列目]や位置（緯度経度）[14,15 列目]も入っている。
※WIN システムについては，http://eoc.eri.u-tokyo.ac.jp/WIN/ を参照のこと．

6

kyo2 (XXX.XX.XX.X) からファイルを転送する方法
Linux から
scp もしくは sftp コマンドを利用% scp kuge@XXX.XX.XX.X:/home/seis/kuge/2018/test.f.

% sftp kuge@XXX.XX.XX.X
フォルダの場合は –r オプション
Windows からWinSCP（プログラム->通信->WinSCP）で kyo2(XXX.XX.XX.X)にログイン

awk (gawk)
たとえば，asc.dat というファイルの 5 行目から 100 行目までの 1,2 列目だけを切り出して，asc2.dat というファイルに書き出したいとき，% awk '{if (NR >= 5 && NR <= 100) print $1, $2}' asc.dat > asc2.dat
1 列目が時刻歴で，3 秒から 11.5 秒までを切り出して，2 列目だけを asc3.dat というフ
ァイルに書き出したいとき，% awk '{if ($1 >= 3.0 && $1 <= 11.5) print $2}' asc.dat > asc3.dat

7