Shindo_tanaka.exe の説明書

http://www.seisvol.kishou.go.jp/eq/kyoshin/kaisetsu/calc_sindo.htm を参考にした。

●2011.04 田中が作成。（exe およびこの説明書。）

●3成分（水平2成分、上下1成分）の時刻暦加速度データから、
  気象庁の計測震度および震度階を求めるプログラム。

●加速度のデータの単位は（cm/s/s）であること。

●加速度のデータ数の上限値は70000。これ以上の場合、強制終了する。

●加速度のデータのフォーマットは、下記の通り。

Time(s)	NS	EW	UD    
0	0.004	0.003	-0.004
0.01	0.008	0	-0.009
0.02	0.011	-0.004	-0.007

※1行目のヘッダーを必ず付けること。
  時間および加速度データはフリーフォーマットであるが、
  基本的にはCSVまたはテキストデータとする。

●計算結果の整合性は確認済み。

●計算はすべて倍精度で行っている。

************************************************************************
計算方法

計測震度は、平成8 年2 月15 日気象庁告示第4 号に定められている。
具体的な計算手順を以下に示す。

(1)ディジタル加速度記録3 成分（水平動2 成分、上下動1 成分）
   のそれぞれの フーリエ変換を求める。
(2)地震波の周期による影響を補正する3 種類のフィルターを掛ける。
   3 種類のフィルターは以下のようである。
   ただし、f は地震波の周波数(Hz)、y はf に10分の1
   を乗じた値である( y = f/10 )。
   �@周期の効果を表すフィルター
     (1/F)^0.5
   �Aハイカットフィルター
     1/ (1+ 0.694y^2 + 0.241y^4 + 0.0557y^6 + 0.009664y^8 + 0.00134y^10 + 0.000155y^12)^0.5
   �Bローカットフィルター
     (1- exp(-( f / 0.5)^3)^0.5
(3)逆フーリエ変換を行い、時刻歴の波形にもどす。
(4)得られたフィルター処理済みの3 成分の波形をベクトル的に合成する。
      p=( px^2 + py^2 + pz^2 )^0.5

(5)ベクトル波形の絶対値がある値 a 以上となる時間の合計を計算したとき、
   これがちょうど 0.3 秒となるような a を求める。
   ※具体的には、デジタル記録のサンプリング時間間隔をΔt としたとき、
   ベクトル波形p を絶対値の大きい順に並べて、0.3/Δt 番目の値をa とする。
(6)(5)で求めたa を、 I = 2 log a + 0.94 により計測震度 Ｉ を計算する。
(7)計算されをたI の小数第３位を四捨五入し、
   小数第２位を切り捨てたもの計測震度とする。


   |--------------------------|
   |震度階級 | 計測震度       |
   |--------------------------|
   |０       | 0.5未満        |
   |１       | 0.5以上 1.5未満|
   |２       | 1.5以上 2.5未満|
   |３       | 2.5以上 3.5未満|
   |４       | 3.5以上 4.5未満|
   |５弱     | 4.5以上 5.0未満|
   |５強     | 5.0以上 5.5未満|
   |６弱     | 5.5以上 6.0未満|
   |６強     | 6.0以上 6.5未満|
   |７       | 6.5以上        |
   |--------------------------|

                                                                   以上