気象庁 計測震度の算出方法例との比較
2000年１０月６日に発生した鳥取県西部地震の米子市（計測震度＝５．１）を使用する。
震度算出
>> SeismicIntensity('AA06EA01.csv',',',7,1,2,3,100)
ans =
5.1000震度結果は一致する。
オリジナルの加速度波形
入力データとしては同一と確認できる。
オリジナル加速度波形の周波数スペクトル
これも同一。
震度計算のためのフィルター特性
バンドパスフィルタをするためのハイカットとローカット。
黒線が合成して実際に使用するフィルター特性。
周期効果について。
調べてみたが詳細は不明。
推測としては、地震の特性上の補正係数で、周波数の高い方が早々に減衰するため、震度の情報から除去していると思われる。
フィルター補正後の加速度波形
フィルター後の３成分合成加速度
ここまで再現できているので、今回のアルゴリズムは一定の正しさがあると言える。
その他のサンプルによる震度算出
(2011年)東北地方太平洋沖地震のデータセットから、
以下を使用
L311E4E1.csv (東京都 東京千代田区大手町 ５強 5.1 )
L3114B91.csv (宮城県 大崎市古川三日町 ６強 6.2)
結果
>> SeismicIntensity('L311E4E1.csv',',',7,1,2,3,100)
ans =
5.1000
>> SeismicIntensity('L3114B91.csv',',',7,1,2,3,100)
ans =
6.2000
>>気象庁公表値と一致した。
まとめ
- 複雑な演算がある場合、再現がメンドクサイ
- 計算式があっていても、演算順序やFPUの違い、FPUモードの違いで演算誤差がでる。
- 上記に伴い、正しいことの証明がブラックボックスアプローチになる。
- 今回は、気象庁が途中のデータ特性を公表してくれたおかげでステップごとに検証できた。
- よって、中間データを如何に取得するかが難易度を下げるコツとなる。
- FFTが入ると、MATLAB coderがイマイチな感じになる。
- 可変要素数不可
- 要素数が2のべき乗であること必須(性能面を考えると必須の方が良いが)
コメント