【気象庁公表】震度算出方法 【MATLABプログラミング】

はじめに
発端
結果
MATLABによる震度算出方法
ブロック線図
mスクリプト
気象庁データの取得
次のページへ

はじめに

国土交通省気象庁で震度算出方法が公開されている。
https://www.data.jma.go.jp/svd/eqev/data/kyoshin/kaisetsu/calc_sindo.htm

気象庁｜強震観測について｜計測震度の算出方法

計測震度は、震度計内部のディジタル処理によって計算されます

この計算をMATLABで再現するというのが今回の目的となる。

MATLAB、Scilab、Python関係の数値計算関連の話はこちら。
https://www.simulationroom999.com/blog/comparison-of-matlab-python-scilab/

発端

同僚からの相談。

現状は震度計算をMATLABでやっているが、他の言語(C/C++等)でやって欲しいと言われている。
MATLABのコード(mスクリプト)は貰えていない。
よって、「それっぽいmスクリプト書いてくれ！！」とは言われていないが、なんとなく察した。

mスクリプトがあれば、MATLAB coderで一発でC/C++のオートコード生成ができるが、無いとなると自作するしかない。

いろいろ調べると、算出方法の概要は分かったので、とりあえずmスクリプトを書いてMATLAB coderでオートコード生成を試みる。

結果

まず、オートコード生成はできると言えばできる。
しかし、FFT(高速フーリエ変換)関数に於いて、入力要素数固定じゃないとオートコード生成が利かないことが分かった。
さすがに、要素数固定では使い物にならない。
しかし、一応mスクリプトによるアルゴリズムの再現は出来たので、普通にC/C++でゴリゴリ書ける水準までは来たと言える。
それでOKとしよう。
場合によっては、オートコードをベースに手修正しても良い。

MATLAB coderによるFFTでは、静的配列変数として三角関数テーブルとビット反転テーブルを宣言する方式と取っている。
この部分をmalloc的な動的配列にして、各テーブルも演算前に要素数に合わせて作ってしまえばOK。

FFTでオートコード生成する場合は、要素数を２のべき乗にする必要あり。

以下、公式ページで紹介されているデータと比較しながらの再現方法解説。

MATLABによる震度算出方法

ブロック線図

時間領域、周波数領域、NS加速度、EW加速度、UP加速度、FFT、BPF、IFFT、ノルム化、合成加速度、降順ソート、0.3秒データ抜き出し

小数第3位を四捨五入し、小数第2位を切り捨てしたものを震度とする。

mスクリプト

まずは書き出したmスクリプト。

至る所にglobal変数を配置しているが、途中計算確認用に使用しているのみ。
global定義を全部削除しても動作する。

function y = SeismicIntensityBody(ns,ew,ud,fs)

    global fft_ns;
    global fft_ew;
    global fft_ud;

    % FFTで周波数領域に
    fft_ns = fft(ns);\n
    fft_ew = fft(ew);\n
    fft_ud = fft(ud);\n
    
    global num;
    num = size( ns, 2 );       % データ数(列数)
    n = 0:( num - 1 );         % プロット軸元ネタ    
    %f =  n / (num/fs);
    % 0割回避＆疑似0値のため微小値を加算
    f =  n / (num / fs) + 0.00000000000000000001; %関数窓用プロット軸
    
    % 乗除は行列乗除("*" , "/")ではなく要素単位の乗除(".*" , "./")
    
    % 周期効果関数窓
    global winX;
    winX = sqrt(1 ./ f);
    
    % ハイカット関数窓
    Y = f ./ 10;
    global winY;
    winY = 1 ./ sqrt(ones(1,num) + 0.694*Y.^2 + 0.241*Y.^4 + 0.0557*Y.^6 + 0.009664*Y.^8 + 0.00134*Y.^10 + 0.000155*Y.^12);
    
    % ローカット関数窓
    global winZ;
    winZ = sqrt(ones(1,num) - exp(-(f./0.5).^3));
    
    % 関数窓合成
    win1 = winX .* winY .* winZ;    % 前半分用
    win2 = fliplr(winX) .* fliplr(winY) .* fliplr(winZ); % 後半用
    
    nn = floor(num/2);
    global win;
    win = [win1(1:nn),win1(nn+1),win2(nn+2:num)];
    
    % バンドパスフィルタ
    global fft_ns_;
    global fft_ew_;
    global fft_ud_;
    fft_ns_ = win.*fft_ns;
    fft_ew_ = win.*fft_ew;
    fft_ud_ = win.*fft_ud;
    
    % 逆フーリエ変換
    global res_ns;
    global res_ew;
    global res_ud;    
    res_ns = ifft(fft_ns_);
    res_ew = ifft(fft_ew_);
    res_ud = ifft(fft_ud_);
    
    % 合成加速度
    global anorm;
    anorm = abs(sqrt(res_ns.^2 + res_ew.^2 + res_ud.^2));
    
    % 降順ソート
    sorted_anorm = sort(anorm,'descend');
    
    % 0.3秒後のデータサンプル確保(開始位置を0秒とする。)
    S = sorted_anorm(floor(0.3*fs)+1);
    
    % 気象庁指定の演算で震度算出
    I = 2 * log10(S) + 0.94;
    y = floor(10*(I+0.005)) / 10;   % 小数第3位を四捨五入し、小数第2位を切り捨て

※全コードはGitHubで公開

GitHub - KEIKEI999/SeismicIntensity: 気象庁で開示されている震度計算をMATLABスクリプトで書いた

気象庁で開示されている震度計算をMATLABスクリプトで書いた. Contribute to KEIKEI999/SeismicIntensity development by creating an account on GitHub.

気象庁データの取得

気象庁で公開されているデータはcsv形式で、先頭7行がヘッダで、8行目以降からが実データとなる。
よって、dlmread関数を使用して行列として抜き出すことができる。

data = dlmread([ファイル名],',',7,0);

1列目がNS(NorthSouth：南北)加速度、2列目がEW(EastWest：東西)加速度、3列目がUD(UpDown：上下)加速度となる。
よって、データをスライシングする。

ns = data(:,1)';     % NSデータ スライシング
ew = data(:,2)';     % EWデータ スライシング
ud = data(:,3)';     % UDデータ スライシング

これを上記関数に渡す。

y = SeismicIntensityBody(ns,ew,ud,100);

次のページへ

次のページから実際の利用方法と何となくの正しさの証明。

【気象庁公表】震度算出方法【MATLABプログラミング】