Pythonで地震計の地震加速度のデータから地震震度を計算する - ServersMan@VPS（CentOS）でお気楽サーバー運営 (^^♪　　(忘れっぽいので個人メモ用)

興味本意で^^;地震情報関係のサイトを今作っているのですが、地震計のデータから地震の計測震度を計算する必要が出てきて、Pythonで地震計の地震加速度のデータから計測震度を計算する関数を作ったので、ちょっとメモ。

防災科学技術研究所強震観測網(K-NET,KiK-net) さんで公開されてる地震計の地震加速度のデータファイルから、地震の計測震度をPythonで計算してみる。

K-NETさんで公開されてる地震計のデータファイルは↓こんな感じのフォーマットなので、

----+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
Origin Time       2006/01/18 23:25:00<LF>         (1) 地震発生時刻
Lat.              37.798<LF>                      (2) 震央北緯
Lon.              142.200<LF>                     (3) 震央東経
Depth. (km)       36<LF>                          (4) 震源深さ
Mag.              5.7<LF>                         (5) マグニチュード
Station Code      MYG002<LF>                      (6) 観測点コード
Station Lat.      38.7262<LF>                     (7) 観測点北緯
Station Lon.      141.5109<LF>                    (8) 観測点東経
Station Height(m) 79<LF>                          (9) 観測点標高
Record Time       2006/01/10 23:25:46<LF>         (10) 記録開始時刻
Sampling Freq(Hz) 100Hz<LF>                       (11) サンプリング周波数
Duration Time(s)  108<LF>                         (12) 計測時間
Dir.              N-S<LF>                         (13) チャンネル
Scale Factor      3920(gal)/6182761<LF>           (14) スケールファクタ
Max Acc. (gal)    48.060<LF>                      (15) 最大加速度
Last Correction   2006/01/18 23:25:31<LF>         (16) 最終校正時刻
Memo.             <LF>                            (17) 備考
    8009     9006     7998     8001     8009     8007     8002     8005<LF>
    8006     8003     8003     8005     8005     8003     8001     8004<LF>
 :
 :
 :
    8006     8006     8000     8000     8005     8009     8004     7999<LF>
    8003     8007     8004     8000     8000     8005     8009     8005<LF>
[EOF]
----+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8

東西、南北、上下動のアスキーフォーマットの地震計のテキストデータ（つまり文字列）を読み込んでそれぞれ、dataEW、dataNS、dataUDという変数にセットしてあげて、calc(dataEW, dataNS, dataUD)と呼び出すと、地震計が記録したデータから計測震度が計算出来る。

ちなみに、計測震度の計算方法は気象庁 | 計測震度の算出方法にあるとおり、以下。（詳しくはネットで見つけたこちらのPDFが参考になります。『地震動強さ指標と新しい気象庁震度との対応関係』）

ディジタル加速度記録３成分（水平動２成分、上下動１成分）のそれぞれのフーリエ変換を求める。

地震波の周期による影響を補正するフィルターを掛ける。

逆フーリエ変換を行い、時刻歴の波形にもどす。

得られたフィルター処理済みの３成分の波形をベクトル的に合成をする。

ベクトル波形の絶対値がある値 a 以上となる時間の合計を計算したとき、これがちょうど 0.3秒となるような a を求める。

5.で求めたaを、　I = 2 log a + 0.94　により計測震度Ｉを計算する。計算されたIの小数第３位を四捨五入し、小数第２位を切り捨てたものを計測震度とする。

#!/usr/bin/env python
# -*- coding: utf-8 -*- 

import numpy as np 
from scipy import sqrt, power, exp, log10
from scipy.fftpack import fft, ifft, fftfreq
import tarfile


def getEarthquake(data):

    tokens = data.split() 
    
    # Scale Factor の値を取り出す 
    (Scale, Factor) = tokens[tokens.index('Factor')+1].split('(gal)/') 
    
    # サンプリング周波数
    fs = tokens[tokens.index('Freq(Hz)')+1].rstrip('Hz')

    # 波形データを取り出す 
    items=tokens[tokens.index('Memo.')+1:] 

    w0 = np.array(items, dtype=np.float64) 
    w0 = float(Scale)  / float(Factor) * (w0 - w0[0])
    maxgal = max(abs(w0)) 

    return maxgal, fs, w0
    


def calc(dataEW, dataNS, dataUD):

    # 波形データを取り出す 
    (maxgalE, fsE, sigE) = getEarthquake(dataEW)
    (maxgalN, fsN, sigN) = getEarthquake(dataNS)
    (maxgalU, fsU, sigU) = getEarthquake(dataUD)

    # サンプリング周波数
    fs = int(fsE)

    # サンプル数
    L = len(sigE)

    # 地震波形をフーリエ変換してフーリエスペクトルを求める
    spectrumE = fft(sigE)
    spectrumN = fft(sigN)
    spectrumU = fft(sigU)

    # 周波数軸の値を計算 
    freqList = abs( fftfreq(L, d=1.0/fs) )

    # フィルターの作成 
    # 周期効果フィルター用の関数窓を作成
    winX = sqrt( 1 / freqList[1:] )
    winX=np.concatenate(([0],winX))
    # ハイカットフィルター用の関数窓を作成
    y = freqList / 10.
    winY = 1.0/sqrt(1.0 + 0.694*power(y,2) + 0.241*power(y,4) + 0.0557*power(y,6) + 0.009664*power(y,8) + 0.00134*power(y,10) + 0.000155*power(y,12))
    # ローカットフィルター用の関数窓を作成
    winZ = sqrt(1.0 - exp( power( -1.0 * (freqList / 0.5), 3) ))
    # ３つの関数窓を合成
    win = winX * winY * winZ

    # スペクトルにフィルターをかける
    spectrum_winE = win * spectrumE 
    spectrum_winN = win * spectrumN
    spectrum_winU = win * spectrumU

    # フィルターをかけたスペクトルを逆フーリエ変換して時系列の波形に戻す。
    resyn_sigE = ifft(spectrum_winE)
    resyn_sigN = ifft(spectrum_winN)
    resyn_sigU = ifft(spectrum_winU)

    # 3成分のベクトル合成波形を求める。
    S2 = sqrt(power(resyn_sigE,2) + power(resyn_sigN,2) + power(resyn_sigU,2))

    # 0.3秒間の地震波形のサンプル数は「int(0.3*fs)」個なので、合成波形データを降順でソートして、地震加速度の大きい方から数えて「int(0.3*fs)」個目の加速度から地震計の計測震度を求める
    maxS2 = sorted(S2, reverse=True)
    I = 2*log10(maxS2[int(0.3*fs)]) + 0.94
    I = int( round(I.real,2) * 10 ) / 10.0
    
    return I