前説

著作権表明

• ソースコードはGPLとして公開します。簡単に言えば著作者名と改変者名の連続性を保持すれば無償で誰でも利用（改変やコピー）でき、商用に使っても手数料を取って配布してもかまいません。

言葉の規定

計測震度	気象庁が公開しているアルゴリズム（計算法）によって計算した値
ガル	加速度　1cm/s2　が　1ガル
重力加速度	地域によって少し異なりますが、国際標準では980.665cm/s2として扱います
	京都大学地質学鉱物学教室重力室（国際其準点）におけるgの値＝9．7970727m/s2 東京大学理学部化学館地下原点室におけるgの値＝9．7978872m/s2（理科年表2000）

原典

計測震度は、震度計内部で以下のようなディジタル処理によって計算されます。　気象庁の解説WEB　URL ディジタル加速度記録３成分（水平動２成分、上下動１成分）のそれぞれに、フーリエ変換・フィルター処理・逆フーリエ変換の手順で、以下に示す特性のフィルターを掛ける。 　計測震度の計算で使われているフィルター処理は、周波数 0.5 - 10Hzの範囲で地震動の加速度と速度の中間の波形を求めていることに相当します。つまり、両対数のグラフ上で見ると、このフィルターの特性曲線の傾きが上の周波数範囲で -1/2 となっています。また、1.0Hz で倍率が１となるよう定数が選ばれています。また、フィルターの式は、以下の３つの部分からなります。 SQRT(1/F) 1/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) SQRT(1-EXP(-(F/0.5)**3)) 　ここで、Fは周波数（Hz）、XはF/10です。１は上で述べた加速度と速度の中間の特徴を表すフィルター、２はハイカットフィルター、３はローカットフィルターです。 　従来から用いられてきた、最大加速度を震度に換算するいわゆる河角の式との違いは、加速度記録に低周波数側を強調する上記のようなフィルターを施したうえ、最大値そのものではなく0.3秒以上継続する値を使う点です。 得られたフィルター処理済みの記録３成分から、ベクトル波形を合成する。 ベクトル波形の絶対値がある値 a 以上となる時間の合計を計算したとき、これがちょうど 0.3秒となるような a を求める。 この a から　I = 2 log a + 0.94　により計測震度 Ｉ を計算する。

「 震度を知る」　−知識とその活用−　気象庁監修　Ａ５版・２３８ページ　定価3,000円

実装

条件

言語：　C++のクラスライブラリとして提供します。他の言語への変換はご自分でしてください
コンパイラ：　ボーランドのC++ビルダ　Ver5
ライブラリ：　数理設計研究所の行列・ベクトル計算パッケージ（vecmath.h)をテンプレートとし、FFTパッケージも使います

構成

• 計測震度の計算　sindo.h へのリンク sindo.cpp へのリンク
  • フーリエ変換　fft.h へのリンク FFT.CPP へのリンク
  • ベクトル支援　vecmath.h へのリンク

解説

sindo.h　sindo.cpp

class Sindo_Base 基本クラス　気象庁の計測震度をそのまま実現している
class Sindo : public Sindo_Base　実用するクラス　周囲雑音を適当に排除するアルゴリズムを追加したもの

　どちらを使っても計測震度としては問題ない。計測震度の実証試験は公表されている新潟中越地震と阪神大震災の地震波形で同じ値になることを確認した。

関数の使い方：

データ投入→100SPS観測してV3D（３元ベクトル）にgalの単位にしてひとつづつ、またはn個のパッケージにしてset_galを呼び出す。
帰り値が１なら必ず計測震度を引き取ること。。計測震度は呼び出すたびにゼロクリアされる。
int set_gal(V3D gal);
int set_gal(V3D *gal, int n);

上記の帰り引数が１ならば、計測震度を得ることができる、計算そのものは終わっている
double get_sindo(void);

class Sindo : public Sindo_Base　には周囲環境の微細振動を差し引くことができるので、記憶をリセットする関数もある
void reset_noise(void);

sindo.cpp 内部計算

既定値

• １００SPS　毎秒100サンプル
• 内部バッファのサイズ＝FFT計算のサイズ＝２０４８であり、１００SPSならば20.48秒
• 計測震度の算出間隔＝５０サンプル単位＝0.5秒毎に計算する

set_gal　内部アルゴリズム　Sindo_Base判

1. 持ち込まれた観測値が算出間隔50個より小さければ、そのサイズ、大きければ50個、過去のバッファを動かして空きを作り、そこに新しい値をコピーする。残り数がなくなるまで以下の処理をする。X,Y,Zそれぞれのバッファを計算用にコピーしてCalcを呼び出す
2. Calc内の処理　　このCalc内で計測震度の計算をしている
   1. X,Y,ZそれぞれについてFFTして周波数軸に変換
   2. 気象庁ご推薦のフィルタ係数（あらかじめ初期化してある）を乗算
   3. X,Y,ZそれぞれについてIFFTして時間軸に戻す
   4. X,Y,Zの絶対値を作り、IFFTした際の規格係数を乗算
   5. 継続時間を考慮した振幅の決定　→　最大値を調べて、最大値をてがかりに分解能４０００のヒストグラムを作成
   6. ヒストグラムの上位から探索し、0.3秒分＝30個になるgal値を求め
   7. KeisokuSindo = 2.0 * log10(kg) + 0.94;　で計測震度値を更新する
   8. round(KeisokuSindo, 1); 　これは小数点１位で四捨五入するテンプレートである。
      小数点３位にしたいのなら　round(KeisokuSindo, 3); 　と変更せよ
3. Calcで計算した計測震度が、set_galに入ってきたときより大きければ１を持ち帰る・・・・つまり新しい計測震度！

Calc　には　Sindo_Base判とSindo判の2種がある。上記解説はSindo_Base判である。Sindo判は周囲雑音をエネルギーとして蓄積し、これを差し引く操作が入っている。実用的には気象庁発表の震度情報に良くあった値を示すが、計測震度の計算としては標準ではない。しかし、気象庁の観測拠点のように静かな場所に置く事ができない一般向けの計測震度計算機として良い手法であるかもしれないと思って実験中である。

..end