１． はじめに
hypoTD_ver2 は、過去に発生した地震を用いた相対震源決定法であり、ver.1 に比べ、ユーザーがより簡単に操作できるよう変更されている。
メインプログラム（ateq_td）の実行で、テンプレート地震の登録、指定された回数のイタレーションによる相対震源決定が行われるようにした。
震源の緯度、経度、深さの推定誤差は、horiuchi et al.(2025) に記載されている方法で求めるようにした。
走時計算は、１）２層傾斜構造、２）多層構造、３）気象庁走時表 の３種類から選べるようにした。
本ソフトウェア利用で必要な入力ファイルは、
１）到着時刻データと初期震源、
２）速度構造、
３）観測点座標、
４）使用するファイル名を指定するファイル etc/wave_dir、
５）地震とノイズとを識別する判定パラメータのテーブル（../etc/hypo_remove.dat）である。
出力データは、hypoTDで計算された震源リスト（arriv_dd/hypo**）である。
２．Directory構造
(1) JMA_measure、win_measure、hypoDD_measure
気象庁検測値データ、Winフォーマット検測値ファイル、hypoDDによるフォーマットを、hypoTD用データフォーマットに変換するソフト
(2) auto_td：hypoTDによる相対震源決定（ateq_td）のソースコード
(3) arriv：hypoTDのフォーマットに変換した検測値データ
(4) etc：速度構造、観測点の座標、入力ファイル名（wave_dir.dat）、地震とノイズを識別する判定パラメータのテーブル（../etc/hypo_remove.dat）
(5) arriv_dd：hypoTDによる震源決定結果の震源リスト
３．気象庁検測値データ、Winフォーマット検測値ファイル、hypoDDフォーマットを hypoTD 用データフォーマットに変換するソフト
以下を実行することにより、検測値データを hypoTD のフォーマットで arriv ディレクトリーに格納される。
以下の例での 1984 は、1984年のデータを表す。
arriv/arv_1984.dat（バイナリー形式の、年毎の到着時刻データ）
arriv/hypo_1984.dat（フォーマット付き、ダイレクトアクセス形式の年毎の震源リスト）
arriv/event.dat（フォーマット付き、ダイレクトアクセス形式の全ての地震の震源リスト）
注）event.dat は、ateq_td を実行すると、計算された震源リストに置き換えられる。
（１）気象庁一元化震源到着時刻データ
1）コンパイル
cd JMA_measure
make
２）入力データ
検測値データのファイル名を print.dat に入れる。複数個ある場合は複数個。
print.dat の例：
measure_20240102_1.txt
３）実行
./read_jma_arriv
（２）Winフォーマット
1）コンパイル
cd win_measure
make
２）データ
print.dat に読み込むピックファイルのディレクトリー名を入れる。複数個ある場合は複数個。
print.dat の例：
picks/1611
３）実行
./read_win
（３）hypoDDフォーマット
1）コンパイル
cd hypoDD_measure
make
２）データ
print.dat に読み込むファイル名を入れる。複数個ある場合は複数個。
print.dat の例：
test1.pha
３）実行
./read_lv
４．入力データ
到着時刻以外の入力データは、
１）観測点座標
２）速度構造
３）データ選択ファイル
４）地震とノイズを識別するテーブル
（１）観測点座標
例１：
etc/jma_stn.dat
観測点 番号 緯度 経度 海抜(km)
EHTO 1 36.5354 140.5136 0.041
ETKM 2 36.4386 140.5687 0.009
EMTH 3 36.3396 140.5182 0.025
EIBR 4 36.2575 140.4601 0.027
ENGK 5 36.1456 140.4295 0.020
例２：
etc/stn_dd.dat
NCAAR 39.275936 -121.026962
NCAAS 38.430138 -121.109589
NCABJ 39.165771 -121.192993
NCABR 39.138126 -121.488167
・観測点番号は、入れなくてもよい。
・海抜が書かれていない場合は、0.0km が指定される。
・観測点ファイル名を間違えた場合は、観測点コードを書いて計算が終了する。
（２）速度構造モデル
etc/vel_model.dat に速度モデルを入れる。
例：
# はコメント行
# Component Name
UD NS EW AZ AX AYMZ MX MY SZ SX SY TZ TX TY !削除しないこと
# Select type of velmodel. 1; two layer, 2; JMA table, 3; N Layer N<50
1
# 速度構造例（コメント付き）
# 5.40 -59.0 7.75 -2.3 3.15 -59.0 4.35 -2.3 31.0 !Tohoku Univ.
# 5.50 -53.49 7.80 -2.13 3.25 -44.68 4.41 -1.57 32.0 !Ukawa
6.00 -30.0 8.0 -2.3 3.53 -30.0 4.50 -2.3 40.0
# Number of layers
# depth Vp Vs(km)
6
-2.00 4.800 2.844
0.50 4.950 2.931
1.00 6.000 3.521
24.00 6.648 3.858
40.01 8.10 4.50
700.00 12.0 6.66
＃を除く行
１行目 未使用。 消さないこと。
２行目 1 速度モデルの選択 1,2,3を入れる。
１；２層傾斜モデル
２；気象庁走時表
３；N層モデル（層内の速度は一定）
３行目 ２行目を、１にした場合の速度モデル
6.00 -30.0 8.0 -2.3 3.53 -30.0 4.50 -2.3 40.0 とした場合の速度は以下で与えられる。
２層傾斜モデルの速度計算例：
Vp = 6.0 * {(R - dep) / R} -30.0
Vs = 3.53 * {(R - dep) / R} -30.0
２層目の深さ：40km
２層目の速度：
Vp = 8.0 * {(Rm - dep) / Rm} -2.3
Vs = 4.5 * {(Rm - dep) / Rm} -2.3
R = 6371.05
Rm = 6371.05 - 40.0
と表せる。
４行目 ２行目を３とした場合の速度モデルの層数
５行目以降
N層モデルの深さ、P波速度、S波速度
この例だと、深さ-2km(標高２ｋｍ) から、0.5kmまでの、Ｐ波速度は、4.8km/sec, S波は2.844 km/sである。
N層モデルでは、平行層モデルが使われており、地球の半径は、考慮されていない。
（３）データ選択ファイル
etc/wave_dir.dat
例：
# Input files of ateq_td
3 !1 Number of iteration
# station code file
../etc/stn_dd.dat !2 Station for hypoDD sample data
# File name of Event list for re-location
../arriv/event.dat !3 Relocation list including template EQ
# lat1,lat2,lon1,lon2,dep1,dep2
-90.0 90.0 -180.0 180.0 -2.0 700.0 !4 Relocation area
# Magnitude range for the selection of Template EQ
1.0 6.0 !5 Magnitude range for template EQ selection
# Minimum numbers of P and S arrival times for the selection of Template EQ
0 0 !6 Minimum P and S readings for template EQ selection
# Shallowest depth for templete EQ selection
0.1 !7 Shallowest depth(km) for template EQ selection
説明：
以下の行番号では、コメント行は除く。
１行目 イタレーションの回数
２行目 観測点データのファイル名
３行目 全地震の震源リストのファイル名
震源リストには、地震番号が書かれており、検測値データ（バイナリーデータ）と、震源リストが紐付けされている。地震番号は、年毎に指定されている。
４行目 震源決定を行う地震の緯度、経度、深さ範囲。
５行目 テンプレート地震として、登録する地震のマグニチュードの範囲。
６行目 テンプレート地震として、登録する地震のP波、S波読み取り数の下限。
２０ １０ とすると、 P波、S波読み取り数が２０，１０以上で、且つ、指定されたマグニチュードの範囲の地震を登録する。
７行目 テンプレート地震として登録する地震の上限の深さ（ｋｍ）。
0.1 とすると、震源の深さが地表に決定された地震は、テンプレート地震として、登録されない。
（４）地震とノイズを識別する判定パラメータのテーブル
../etc/hypo_remove.dat
例：
0.1 !1 err_min Minimum sd_cross
2.0 !2 rat_stn_sd2 remove error picks sd2=2.0*sd
0.4 !3 sd_remove remove hypo when large sd (sec)
10.0 !4 dl_remove remove hypo when large location error (km)
0.2 !5 rat_err_eq remove hypo if 60% of large sd eq.
70.0 !6 del_lim_sd_change km
5.0 !7 magnitude to remove S arrival
説明：
１行目 ；相対震源決定の走時残差の最小値
２行目 ；読み取りデータを削除する場合の基準
１行目が 0.1 ２行目が2.0の場合が、
50個の地震の走時残差のRMSをsd, 走時残差の平均値からのずれを av とすると
Sd が0.1秒より小さければ、sd=0.1 とし、av>2sd であれば、その観測点データを削除
３行目 50個の地震の走時残差の中で、RMSがこの値以上であれば、そのテンプレート地震との震源決定を行わない
４行目 50個の相対震源決定を行った場合の、震源位置のばらつきが、この値以上であれば、震源決定を行わない。
５行目 50個の相対震源決定で、震源決定できた地震の割合が、この値以下であれば、震源決定を行わない。
６行目 平均的観測点間距離の閾値。 例えば、70km で、平均的観測点間距離が、Lとすると、３行目、５行目のデータは、
Lが70kmを超えた場合に、３行目、５行目のデータを、L/70で、除して、判定に使う。
７行目 マグニチュードがこの値以上であれば、S波到着時刻を使用しない。
震源決定できない地震の割合が数％以下であれば、上記パラメータは変更不要。
５. hypoTDによる震源決定
（1）コンパイル
auto_td ディレクトリーでmake とすると、ateq_tdがコンパイルされる。
（2）実行
cd auto_td
./ateq_td
６．hypoTDの出力
arriv_dd に、震源リストが作成される。
（1）ファイル名：arriv_dd/hypotd_asc_ita.dat
ita：iteration の回数
ファイル形式：アスキーフォーマット
例：hypotd_asc_001.dat
write(38,1302) iyear,im,id,ih,imin,t,
1 phi,dphi,ram,dram,dep,ddep,fmag,nobsp,nobss,sd_all_cross
1302 format(i4,4i2,f7.2,f8.4,f7.4,f10.4,f7.4,2f7.2,f5.2,2i4,f6.2)
Date origin Lat err Long err Depth err Mag Np Ns RMS
1984 1 31931 5.80 36.8646 0.0104 -121.4131 0.0155 3.91 1.52 2.52 93 0 0.16
1984 1 32058 34.57 36.8755 0.0007 -121.6182 0.0015 3.37 0.29 2.39 62 0 0.14
1984 1 51756 17.10 36.7045 0.0008 -121.3740 0.0008 1.25 0.23 1.44 19 0 0.13
震源決定できなかったイベントは、初期震源が入っており、sd_all_cross = 9.9 になっている。
（2）ファイル名：arriv_dd/hypotd_dir_ita.dat
ita：iteration の回数
ファイル形式：フォーマット付きダイレクトアクセス、レコード長 120 byte
例：hypotd_dir_001.dat
322
84 1 31931 0.00 5.80 0.00 36.8646 0.0104 -121.4131 0.0155 3.91 1.52 2.52 93 0 93 3 0.16
84 1 32058 0.00 34.57 0.00 36.8755 0.0007 -121.6182 0.0015 3.37 0.29 2.39 62 0 62 4 0.14
84 1 51756 0.00 17.10 0.00 36.7045 0.0008 -121.3740 0.0008 1.25 0.23 1.44 19 0 20 5 0.13
１行目：地震数
２行目以降：
（年-1900）, 月, 日, 時, 分, 秒, 発震時, 誤差, 緯度, 誤差, 経度, 誤差, 深さ, 誤差, マグニチュード, P波読み取り数, S波読み取り数, 観測点数, イベント番号, 相対震源決定の走時残差の標準誤差
character*1 cret, c_temp, c_hypo_flag
open(7,file=file_hypo,access='direct',form='formatted',
1 recl=120,status='unknown')
write(7,1303,rec=irece) cret,iy,im,id,ih,imin,sec,t,dt,
1 phi,dphi,ram,dram,dep,ddep,fmag,nobsp,nobss,nobsall,
1 c_temp,c_hypo_flag,irec_arv,sd_all_cross
1303 format(a1,1x,5i2,f6.2,f8.2,f6.2,f8.4,f7.4,f9.4,f7.4,2f7.2,f5.2,3i4,1x,2a1,i9,f6.2)
cret 改行
iy,im,id,ih,imin,sec 時刻
t,dt origin time とその誤差
phi,dphi,ram,dram,dep,ddep 緯度・経度・深さとその誤差
fmag マグニチュード
nobsp,nobss,nobsall P波・S波読み取り数・観測点数
c_temp,c_hypo_flag 気象庁フラグ
irec_arv イベント番号
sd_all_cross 走時残差の標準誤差
震源決定できなかったイベントは、気象庁震源が入っており、sd_all_cross = 9.9 になっている。
７．気象庁走時表の変更方法
１）走時表を変更する場合は、3.入力データ ２) etc/vel_model.dat に速度構造モデルを入れる。
用いる走時表のタイプを ２ にする。
２）../etc/jma_trav_del_dep.dat を変更する。
３）深さ、距離方向のグリッド数が、２３６，１０６を超える場合は、
common /c_trav_jma/ パラメータの引数を増やす。
以下は sub_trv.f の最後の、走時表読み取りプログラムである。
subroutine read_jma_trv
common /c_trav_jma/idelg(236),idepg(106),travp(106,236),
1 travs(106,236)
open(7,file='../etc/jma_trav_del_dep.dat')
ndel=236
ndep=106
do kdep=1,ndep
do kdel=1,ndel
read(7,100) travp(kdep,kdel),travs(kdep,kdel),
1 idepg(kdep),idelg(kdel)
100 format(1x,f9.3,2x,f9.3,i4,i7)
enddo
enddo
close(7)
return
end
８．リアルタイム相対震源決定の方法
１）etc/wave_dir.dat の変更
２）ateq_td.f の中の以下の system call をコメントにする。
text='cp ../arriv_dd/hypotd_dir.dat '//event_list
ccc call system(text) !cp calculated hypo for the next iteration
３）ateq_td.f の 179行目 subroutine read_new_event_test(iflag) を変更し、
新しい地震が発生したら、その地震を震源決定するようにする。
９．ダウンロード
以下のリンクをクリックすると、hypoTD Ver.2 をダウンロードできます。
hypoTD Ver.2 ソースファイルをダウンロード(1.08MB)
ご質問のある方は、「お問い合せ」ページ記載のメールアドレス宛にご連絡ください。
|