(cache)加速度から計測震度を計算してみる

ふと気象庁が発表している震度を計算したくなったので、最近よく使うtypescriptで実装しようと思います。

前提

この記事ではある程度プログラムができそこそこ地震観測に詳しい人向けです。

計算方法

気象庁の震度は、1分長の3成分加速度情報をフーリエ演算して各成分ごとにフィルターを掛け合わせてからの逆フーリエ演算、その後各成分を合成し3成分合成加速度上位0.3秒目の値から震度算出と、文字に起こすとややこしいものになっています。

実装

今回、FFT(高速フーリエ演算)のライブラリはfft-jsを使うことにします。
データは気象庁強震観測データを利用します。

2011/03/11 14:46:18に発生した地震で、大崎市古川三日町の観測所のデータを使います。
110311_4B9.json に3成分加速度情報を2次元配列にしてJSON形式で保存します。
※上下、東西、北南の順番はどの順番でもいいです。
※データは3分長ですが計算結果に影響はほぼないので大丈夫。

import { fft, ifft } from 'fft-js';
import gals from './110311_4B9.json';

加速度をFFTに通すために、2のべき乗個にしなくてはいけないので、2のべき乗個になるように後ろに0を追加していき、その後FFT関数を実行します。

const fftGals: [number, number][][] = gals.map(gal => fft(arrayPadding(gal)));

function arrayPadding(arr: number[]) {
    let dips = 2;

    while (dips <= arr.length) {
        dips *= 2;
    }

  return arr.concat(new Array(dips - arr.length).fill(0));
}

次に地震波の周期による影響を補正するフィルターを各成分ごとにかけます。

// サンプリング間隔(seconds)
const dt = 0.01;

const NF = fftGals[0].length;
const NF2 = NF / 2;

for (let i = 1; i <= NF2; i++) {
    // 周波数
    const frq = i / (NF * dt);

    const x = frq / 10;
    const fc = Math.sqrt(1 / frq);
    const fh =
        1 /
        Math.sqrt(
            1 +
            0.694 * x ** 2 +
            0.241 * x ** 4 +
            0.0557 * x ** 6 +
            0.009664 * x ** 8 +
            0.00134 * x ** 10 +
            0.000155 * x ** 12
        );
    const fl = (1 - Math.exp((-(frq / 0.5)) ** 3)) ** 0.5;
    const fa = fc * fh * fl;

    // 計算回数を抑える
    for (let lk = 0; lk < fftGals.length; lk++) {
        const fs: [number, number] = [fftGals[lk][i][0] * fa, fftGals[lk][i][1] * fa];

        fftGals[lk][i] = fs;
        fftGals[lk][NF - i] = [fs[0], -fs[1]];
    }
}

フィルターを通したらIFFT（逆フーリエ変換）をして3成分合成加速度を得ます。

const filterGals = fftGals.map(fftGal => ifft(fftGal));

const lenGals: number[][] = [];

filterGals.forEach((gal, i) => {
    gal.forEach(([gas], t) => {
        (lenGals[t] || (lenGals[t] = [])).push(gas);
    });
});

// 各成分を合成
const syntheticGal: number[] = lenGals.map(([s1, s2, s3]) => Math.sqrt(s1 ** 2 + s2 ** 2 + s3 ** 2));

合成加速度を得たら震度を算出していきます。

// 合成加速度を降順に直す。
syntheticGal.sort((a, b) => b - a);

// 合成加速度上位0.3秒目のindex値を得る。
const top03 = Math.floor(0.3 / dt) - 1;

// 震度計算
const int = 2 * (Math.log(syntheticGal[top03]) / Math.LN10) + 0.94;

// 気象庁の告示に従い計測震度に直す。
const kInt = (Math.floor(Math.round(int * 100) / 10) / 10).toFixed(1);

console.log('計測震度: ' + kInt);
// < 計測震度: 6.2

完成です。

参考
平成8年気象庁告示第4号