JP2014228373A - 計測震度概算装置及びそれを用いた計測震度概算システム - Google Patents

Abstract

【課題】即時性を維持しつつ、さらに換算の誤差の小さい計測震度概算装置及びそれを用いた計測震度概算システムを提供することを目的とする。【解決手段】 地震の加速度を計測し、計測した加速度から計測震度を概算する計測震度概算装置において、地動加速度の時系列を得る地動加速度時系列取得手段７１と、地動加速度の時系列をフィルタ処理する時間領域フィルタ手段７２と、時間領域フィルタによりフィルタ処理された地動加速度の時系列から計測震度の概算値を算出する算出手段７３とを備え,時間領域フィルタ手段は、式（１１）乃至（１６）において、演算１乃至演算７を実行することを特徴とする。【選択図】図１

Description

本発明は、地震の加速度を計測し、その計測値から震度を概算する計測震度概算装置及びそれを用いた計測震度概算システムに関する。

従来から、地震が発生した際の警報発報や鉄道等の制御の指標として、気象庁が官報第１８３１号気象庁告示第４号気象業務法施行規則に関する告示で示した計測震度が用いられているのが現状である。計測震度を算出するためには、地震記録を一定時間蓄積してから、ＦＦＴを用いてフィルタ処理をすることが必要なため、刻々と測定される地震記録に対して値を算出することができず、即時性が要求される用途に適していなかった。

計測震度は、現在最も広く認知されている強震動指標であるが、演算に周波数領域のフィルタ処理を必要とすることから速報性が求められる用途には適していない。これを解決するため、本発明者は、特許文献１で、周波数領域でのフィルタ演算を時間領域の近似フィルタで代用することによって実時間で震度の近似値を得るリアルタイム演算法を提案した。

図１５は、特許文献１に記載された振幅周波数特性を示す図である。実線Ａは特許文献１に記載された振幅周波数特性であり、破線Ｂは気象庁告示に基づくフィルタ振幅周波数特性を、実線Ｃは特許文献１に記載されたフィルタの振幅周波数特性を気象庁告示に基づくフィルタの振幅周波数特性で除した相対的な振幅周波数特性をそれぞれ表している。

特許文献１に記載されたフィルタは、主要な帯域０．１Ｈｚ〜１０Ｈｚで気象庁告示に基づくフィルタとの差違が０．８２８倍から１．０８４倍に収まるように構成されている。これは相対的な振幅周波数特性である実線Ｃから読み取ることができる。このように、特許文献１に記載されたフィルタは、即時性があり、換算誤差の小さい計測震度の概算をすることができた。

特許第４２２９３３７号公報

本発明は、従来のフィルタと比較して、即時性を維持しつつ、さらに換算の誤差の小さい計測震度概算装置及びそれを用いた計測震度概算システムを提供することを目的とする。

前記課題を解決するために、本発明の計測震度概算装置は、
地震の加速度を計測し、計測した加速度から計測震度を概算する計測震度概算装置において、
地動加速度の時系列を得る地動加速度時系列取得手段と、
前記地動加速度の時系列をフィルタ処理する時間領域フィルタ手段と、
前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列から計測震度の概算値を算出する算出手段と、
を備え、
前記算出手段は、
前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列を振幅時系列に変換する振幅時系列変換手段と、
前記振幅時系列を震度換算振幅時系列に換算する震度換算振幅時系列換算手段と、
震度換算振幅時系列換算手段により求めた震度換算振幅時系列を離散化する離散化手段と、
前記離散化手段により求めた各離散化震度換算振幅値に継続時間カウントを実行する継続時間カウント手段と、
を有し、
前記時間領域フィルタ手段は、以下の式（１１）乃至（１６）において、演算１乃至演算７を実行する
ことを特徴とする計測震度概算装置。

・フィルタ1
y(k)=[-α₁y(k-1) -α₂y(k-2)+β₀x(k) +β₁x(k-1) +β₂x(k-2)]/α₀ (15)
ただし、
α₀=8/(ΔT)²+(4ω_a1＋2ω_a2 )/(ΔT)+ω_a1ω_a2 (11)
α₁=2ω_a1ω_a2-16/(ΔT)²
α₂=8/(ΔT)²-(4ω_a1＋2ω_a2 )/(ΔT)+ω_a1ω_a2
β₀=4/(ΔT)²+2ω_a2 /(ΔT)
β₁=-8/(ΔT)²
β₂=4/(ΔT)²-2ω_a2 /(ΔT)
ω_a1=2πf_a1, ω_a2=2πf_a2
x(k)はフィルタの入力時系列
y(k)はフィルタの出力時系列
kは時間ステップ数
ΔTはサンプリング間隔
πは円周率
である。

・フィルタ2
y(k)=[-α₁y(k-1) -α₂y(k-2)+β₀x(k) +β₁x(k-1) +β₂x(k-2)]/α₀ (15)
ただし、
α₀=16/(ΔT)²+17ω_a3/(ΔT)+ ω_a3 ² (12)
α₁=2ω_a3 ²-32/(ΔT)²
α₂=16/(ΔT)²-17ω_a3/(ΔT)+ ω_a3 ²
β₀=4/(ΔT)²+8.5ω_a3 /(ΔT) + ω_a3 ²
β₁=2ω_a3 ²-8/(ΔT)²
β₂=4/(ΔT)²-8.5ω_a3 /(ΔT) + ω_a3 ²
ω_a3=2πf_a3
である。

・フィルタ3
y(k)=[-α₁y(k-1) -α₂y(k-2)+β₀x(k) +β₁x(k-1) +β₂x(k-2)]/α₀ (15)
ただし、
α₀=12/(ΔT)²+12h_b2ω_b /(ΔT)+ω_b ² (13)
α₁=10ω_b ²-24/(ΔT)²
α₂=12/(ΔT)²-12h_b2ω_b /(ΔT)+ω_b ²
β₀=12/(ΔT)²+12h_b1ω_b /(ΔT)+ω_b ²
β₁=10ω_b ²-24/(ΔT)²
β₂=12/(ΔT)²-12h_b1ω_b /(ΔT)+ω_b ²
ω_b=2πf_b
である。

・フィルタ4
y(k)=[-α₁y(k-1) -α₂y(k-2)+β₀x(k) +β₁x(k-1) +β₂x(k-2)]/α₀ (15)
ただし、
α₀=12/(ΔT)²+12h_cω_c /(ΔT)+ω_c ² (14)
α₁=10ω_c ²-24/(ΔT)²
α₂=12/(ΔT)²-12h_cω_c /(ΔT)+ω_c ²
β₀=ω_c ²
β₁=10ω_c ²
β₂=ω_c ²
ω_c=2πf_c
である。

・フィルタ5
y(k)=g_d x(k) (16)
ただし、g_d はゲインである。

演算1： f_a1＝f₀, f_a2＝f₁として，式(11)，(15)の演算
演算2： f_a3＝f₁として，式(12)，(15)の演算
演算3： h_b1 = h_2a, h_b2 = h_2b, f_b＝f₂ として，式(13)，(15)の演算
演算4： h_c = h₃, f_c＝f₃ として，式(14)，(15)の演算
演算5： h_c = h₄, f_c＝f₄ として，式(14)，(15)の演算
演算6： h_c = h₅, f_c＝f₅ として，式(14)，(15)の演算
演算7： g_d＝ g として，式(16)の演算
ただし、
f₀はフィルタのカットオフ周波数を規定するパラメータ
f₁は0.5次で減衰する周波数領域の位置を規定するパラメータ
f₂、h_2a、h_2bは補正フィルタの特性を規定するパラメータ
f₃、h₃はフィルタのカットオフ周波数およびダンピングを規定するパラメータ
f₄、h₄はフィルタのカットオフ周波数およびダンピングを規定するパラメータ
f₅、h₅はフィルタのカットオフ周波数およびダンピングを規定するパラメータ
gはゲイン調整パラメータ
である。

また、本発明の計測震度概算装置を用いた計測震度概算システムは、地震動を計測する複数の地震計と、通信ネットワークと、前記複数の地震計で観測された地震動を、前記通信ネットワークを通じ、離れた場所で前記計測震度概算装置による概算の震度値算出を実行することを特徴とする。

このような計測震度概算装置及び計測震度概算システムにより、従来のフィルタと比較して、即時性を維持しつつ、さらに換算の誤差の小さい計測震度の概算をすることができる。

本実施形態の計測震度概算装置のブロック図である。 本実施形態の計測震度概算部のブロック図である。 本実施形態の計測震度概算値を求めるフローチャート図である。 本実施形態の離散化を示す図である。 本実施形態の継続時間カウントの考え方を示す図である。 本実施形態の継続時間カウントのフローチャート図である。 第２の継続時間カウントの考え方を示す図である。 第２の継続時間カウントのフローチャート図である。 本実施形態の計測震度概算値を示す図である。 本実施形態のフィルタ特性を示す図である。 本実施形態のフィルタを用いて計算されたリアルタイム震度（Ir）と気象庁震度（Ijma）の比較を示す図である。 本実施形態の（ΔI＝Ijma−Ir）のヒストグラムを示す図である。 本実施形態の（ΔI＝Ijma−Ir）のヒストグラムを示す図である。 本実施形態の計測震度概算システムを示す図である。 特許文献１に記載されたフィルタ特性を示す図である。

本発明の実施の形態を図により説明する。図１は、本発明の実施形態の計測震度概算装置の主要構成を示す図である。図中、１は計測震度概算装置、２は処理部、３は制御・演算部、４は計測震度算出部、５はＳＩ値算出部、６は応答値算出部、７は計測震度概算部、８は電源装置、９は計測装置である。

計測震度概算装置１は、計測装置９で計測した加速度等を制御・演算部３に入力することで、震度を概算する。処理部２は、電源装置８により駆動される制御・演算部３を有する。本実施形態の制御・演算部３は、よりよい精度を求めるため、現行の計測震度計に計測震度概算部７を追加した形態としたが、即時性をより重視する場合等では、計測震度算出部４、ＳＩ値算出部５及び応答値算出部６を備えなくてもよい。また、計測装置９を処理部２内に備えてもよい。

図２は、本発明の実施形態の計測震度概算部７の主要構成を示す図である。図中、７１は地動加速度時系列取得手段、７２は時間領域フィルタ手段、７３は算出手段、７４は振幅時系列変換手段、７５は震度換算振幅時系列換算手段、７６は離散化手段、７７は継続時間カウント手段である。

地動加速度時系列取得手段７１は直交３成分の地動加速度の時系列を取得し、時間領域フィルタ手段７２に送る。時間領域フィルタ手段７２は地動加速度時系列取得手段７１が取得した地動加速度の時系列をフィルタ処理する。算出手段７３は時間領域フィルタ手段７２によりフィルタ処理された地動加速度の時系列から計測震度の概算値を算出する。

算出手段７３は、振幅時系列変換手段７４、震度換算振幅時系列換算手段７５、離散化手段７６及び継続時間カウント手段７７を有する。振幅時系列変換手段７４は時間領域フィルタ手段７２によりフィルタ処理された地動加速度の時系列を振幅時系列に変換し、震度換算振幅時系列換算手段７５は振幅時系列を震度換算振幅時系列に換算する。離散化手段７６は震度換算振幅時系列換算手段により求めた震度換算振幅時系列を離散化し、継続時間カウント手段７７は離散化手段により求めた各離散化震度換算振幅値に継続時間カウントを実行する。

次に、図３に示すフローチャートを用いて、このような計測震度概算装置１の計測震度概算方法について説明する。まず、ステップ１で、地震動を観測する（ＳＴ１）。

次に、ステップ２で、時間間隔dt（例えば0.01秒間隔）でサンプリングされた、直交３成分の地動加速度時系列xa[k]、ya[k]、za[k]を得る。ここで、kは時間ステップ数とする（ＳＴ２）。本方法の入力として必要なのは、地動加速度の時系列であるが、地震計の特性補正や速度の微分演算等を施して地動加速度に変換することができれば、必ずしも加速度計により得られたものでなくてもよい。また、地震動の観測は東西方向と南北方向等の水平方向２成分と鉛直方向成分についてそれぞれ行い３成分の地動加速度の時系列を得る。ただし、本方法はいずれかの２成分もしくは１成分の地動加速度の時系列に関しても適用可能である。地震動の観測が水平方向２成分と鉛直方向成分でなされていない場合は、必要な変換を行い３成分の地動加速度の時系列を得る。なお、この変換は、次に示すステップ３のフィルタ処理の後に実行してもよい。

次に、ステップ３で、直交３成分の地動加速度時系列（xa[k]、ya[k]、za[k]）に後述する演算Ａで定めるフィルタ処理を行い、３成分のフィルタ処理された地動加速度の時系列（xf[k]、yf[k]、zf[k]）を得る（ＳＴ３）。なお、フィルタの特性に関しては、後述する。また、地震観測の成分数が３成分に満たない場合は、存在しない（xf[k]、yf[k] 、zf[k]）を0として、ステップ３以降の処理を実行することができる。例えば、x,y,z,の順番は無関係なので、地震観測の成分数が2成分の場合は、xf[k]、yf[k] 、zf[k] のうちいずれかの1つを0、地震観測の成分数が１成分の場合は、xf[k]、yf[k] 、zf[k] のうちいずれかの2つを0とする。

次に、ステップ４で、フィルタ処理された３成分地動加速度の時系列を
af[k]=(xf[k]×xf[k]+yf[k]×yf[k]+zf[k]×zf[k])^1/2 ・・・（Ａ）
なる演算で振幅時系列af[k]に変換する（ＳＴ４）。

気象庁告示第4号の本来の定義に従えば、この振幅時系列から、一定の時間区間内において積算継続時間が0.3秒以上となる最大の振幅レベルa₀を求め、換算式I=2×log(a₀)+0.94で計測震度相当に換算することになるが、計算の効率化のため、計測震度相当に変換してから積算継続時間に関する処理を行う。

まず、ステップ５で、ステップ４において得られた振幅時系列を
If[k]=2×log(af[k])+0.94 ・・・（Ｂ）
なる換算式で、震度換算した震度換算振幅時系列If[k]に変換する（ＳＴ５）。

次に、ステップ６で、図４に示すように、If[k]を必要な精度の離散化幅dIで離散化する（ＳＴ６）。本実施形態では、計測震度が0.1刻みで算出されるため、最終的な結果の十分な精度を保持するためdIは0.1の10分の1である0.01で離散化した。

具体的には、If[k]から下限値Ii（本実施例では-8）を減じた値を離散化幅dIで除し整数化する。このとき整数化した値が負になった場合は0に置き換える。こうして得られる離散化震度換算振幅値の時系列をId[k]とする。次にId[k]の値を上限値N未満に制限する。実施例ではN=1800とし、Id[k]の値がN以上となる場合はN-1に置き換える。これによりId[k]の値は0からN-1の値をとる。実施例においては離散化後の震度換算振幅値が0から1799になるが、これは離散化前の値で-8から9.99に相当する値の制限を行ったことになる。

値の制限範囲は、求めようとする計測震度の範囲を含む必要がある。範囲の上限については、震度7の下限である計測震度6.5よりも十分に大きい値にする必要がある。実施例では、9.99とした。これは計測された加速度が10000galを上回る場合に相当し、地震波センサ自体の計測上限を超える値であるため、十分に大きい。下限については、震度0の上限である計測震度0.4よりも小さい値にする必要がある。実施例では-8とした。これは地震の無い時の微弱な地面の揺れを地震動とみなして計算された計測震度よりも小さい値として設定した値であり、十分に小さい値である。

次に、ステップ７で、各離散化震度換算振幅値について継続時間のカウントを行う（ＳＴ７）。

ここで、図５を用いて第１の継続時間のカウント方法についての基本的な考え方を説明する。図５は、時間ステップ数kと離散化震度換算振幅値との関係を示すもので、図５（ａ）は積算継続時間のカウント範囲を時間ステップ数の２〜４とした場合であり、図５（ｂ）は積算継続時間のカウント範囲を時間ステップ数の３〜５とした場合である。この時、時間ステップ数２に対応する離散化震度換算振幅値はカウント範囲を外れ、時間ステップ数５に対応する離散化震度換算振幅値はカウント範囲に加わることとなる。従って、時間ステップ数３〜４の離散化震度換算振幅値に関しては、時間ステップ数が１増えても変わらないので、そのまま残すようにする。

このような考え方を用いて継続時間のカウントをサブルーチンとして実行する。なお、ここでいう継続時間とは、気象庁告示第４号の本来の定義から明らかであるが、離散化震度換算振幅値の時系列がある一定の振幅に等しい状態が継続した時間ではなく、ある一定の振幅以上となった状態が継続した時間をいう。

図６は第１の継続時間のカウント方法を示すサブルーチンのフローチャートである。あらかじめ、N個の離散化震度換算振幅値の値ごとに継続時間カウント用の記憶領域を確保しておく。この配列は0に初期化しておき、この記憶領域をIc[i]とする。配列の添え字iは0からN-1とする。

まず、ステップ１０１で、上記のId[k]、Ic[i]を用いて、継続時間のカウントを、Id[k]の値以下の添え字iをもつカウント用記憶領域Ic[i]の値をすべて＋1することで実行する（ＳＴ１０１）。

次に、ステップ１０２で、積算継続時間をカウントする時間区間範囲を外れたものについてカウントした値を元にもどす必要がある。これは、Id[k-m]の値以下の添え字iをもつIc[i]の値をすべて-1することで効率よく実行できる（ＳＴ１０２）。ここでmは積算継続時間をカウントする時間区間範囲に対応するサンプル数である。実施例の場合は100Hzサンプリングで積算継続時間をカウントする時間区間範囲を10秒とした。この場合mは1000となる。これを実行するためには、m個の記憶領域を確保し、常にId[k]からId[k-m]までのm個の値を記憶しておく必要がある。このカウント方法により、時間ステップ毎にm個の継続時間に関する処理を行わずにすむ。

次に、ステップ１０３で、積算継続時間0.3秒に対応する総カウント数M（実施例のように100Hzサンプリングの場合は30）以上の値となるIc[i]のうち最大の添え字のものをIc[p]とすれば、pが求まり（ＳＴ１０３）、サブルーチンを終了する。このように、サブルーチンを実行することでpが求まり、図３のフローチャートに戻る。

また、図７を用いて第２の継続時間のカウント方法についての基本的な考え方を説明する。図７（ａ）は、サンプリング周波数を１０Ｈｚ、積算継続時間をカウントする時間区間を１秒として、積算継続時間をカウントする時間区間内の震度換算振幅値を値の大きい順にならびかえたものである。なお、図中で上段の数字は、震度換算振幅値の値を、下段の数字はその値が得られた、時間ステップ数を表している。この場合、積算継続時間が0.3×10Hz=3となる震度換算振幅値は、上から３番目に大きい値６となる。

次に、図７（ｂ）に示すように、時間ステップ数が１つ進んだ場合、図７（ａ）における時間ステップ数０の離散化震度換算振幅値が積算継続時間をカウントする時間区間から外れる。そして、図７（ｃ）に示すように、新たに時間ステップ数１０の離散化震度換算振幅値が積算継続時間をカウントする時間区間に入る。この場合、積算継続時間が0.3×10Hz=3となる震度換算振幅値は、上から３番目に大きい値５となる。

図８は第２の継続時間のカウント方法を示すサブルーチンのフローチャートである。

まず、ステップ２０１で、配列Is[j]の初期化が行われているか判断する（ＳＴ２０１）。
ステップ２０１において、初期化が行われていないと判断された場合、ステップ２０２で、初期化を実行する（ＳＴ２０２）。初期化について説明する。まず、m個の記憶領域を持つ配列Is[j]を用意する。配列の添字jは1からmまでとする。ここで、mは図６のステップ１０２で用いた第１の方法ものと同様に、積算継続時間をカウントする時間区間範囲に対応するサンプル数である。本第２の方法の場合は100Hzサンプリングで、時間区間範囲10秒なので、1000である。配列Is[j]の初期化として、最新の時間ステップ（時間ステップ数をkとする）の1つ前の時間ステップからm個過去にさかのぼり、離散化震度換算振幅値の値Id[k-m]〜Id[k-1]を配列Is[j]にコピーする。次にIs[j]を値の大きい順に並び替える。このときm個の配列It[n]を用意し（配列の添字nは1からm）、並び替えられたIs[j]のどの要素がいつの時間ステップの値であるかを記憶しておく。これで初期化が完了する。

次に、時間ステップ数kにおいて、積算継続時間をカウントする時間区間範囲から外れる値Id[k-m]と、新たに加わる値Id[k]に関する操作のみを行って、Is[j]を値の大きい順に並んだ配列に保つようにする。

具体的には、ステップ２０１において、初期化が行われていると判断された場合、又は、ステップ２０２で初期化が行われた場合、ステップ２０３で、It[q]=k-mとなるqを見つけ、時間区間範囲から外れる値であるId[k-m]に対応しているIs[q]を配列Is[j]から削除する。削除した部分は詰めて配列Is[j]をm-1個の配列とする。配列It[n]も同様にIt[q]を削除し、削除した部分は詰める（ＳＴ２０３）。

次に、ステップ２０４で、Id[k]の値と、配列Is[j]の値を比較し、Is[j]が大きい順にならぶような場所を見つけて挿入し、Is[j]を再びm個の配列とする。配列It[n]も同じ場所に要素を挿入し、その値はkとし追加された値に対応する時間ステップ数を記憶させる（ＳＴ２０４）。

以上の操作が終了すると、Is[j]は大きい順にならんだ配列となっており、It[n]も更新されたIs[j]について、どの要素がいつの時間ステップの値であるかを記憶している。

そして、ステップ２０５で、値の大きい方からM番目の値である積算継続時間0.3秒に対応する離散化震度換算振幅値Is[M]を取得し（ＳＴ２０５）、サブルーチンを終了する。Is[M]は図６のステップ１０３で用いた第１の方法のpに対応する。また、Mは第１の方法と同様に、第２の方法の場合30である。このように、サブルーチンを実行することでpが求まり、図３のフローチャートに戻る。

なお、第２の継続時間のカウント方法の場合、カウントに関する処理は離散化震度換算振幅値の大小比較による。よって、整数値化された離散化震度換算振幅値ではなく、震度換算振幅値そのものを用いて処理することもできる。

また、継続時間のカウント方法は、同様の効果が得られるのであれば、第１、第２の方法以外であってもよい。さらに、データの保持は、配列でなくてもよい。

サブルーチンを終了した後、図３のフローチャートに戻り、最後に、ステップ８で、計測震度の概算値を算出する（ＳＴ８）。pの値に、離散化幅dIを乗じ下限値Iiを加えた値がその時点での計測震度の概算値Ia[k]となる。

以上の処理を、時間ステップ毎に行う。したがって、３成分の地動加速度（xa[k]、ya[k]、za[k]）を得る毎にその時点での震度の概算値Ia[k]を得ることができる。図９に示すように、本方法で求められる震度の概算値Ia[k]は、地震動が継続している間に時間変化する時系列として得られる。この時間変化する値の最大値Ixを、地震動が継続している区間に対して一つの値が決まる本来の定義の計測震度の概算値と定める。地震動が継続している区間を定義しがたい場合は、例えば地震検知後1分間の最大値のように定めればよい。このようにIa[k]は時々刻々と時間変化することから、最大値Ixに達する前に警戒レベルを超えた場合に即時警報を発することも可能になる。

なお、前述した各ステップの順序は、これに限定されるものではなく、計測震度概算値が変わらない限り、ステップの順序を変更してもよい。

次に、気象庁の告示した計測震度本来の定義のフィルタを近似した、時間領域フィルタの設計法について述べる。

近似フィルタの設計は、特許文献１で提案されたフィルタに次の３点の改良を加えることで行う。
（ａ）過大となっている１Ｈｚから１０Ｈｚおよび直流から０．１Ｈｚでの特性を改善するため、ゲイン調整のパラメータを変更する。
（ｂ）過小となっている０．５Ｈｚ付近の帯域での特性を改善するため、補正フィルタを付加する。
（ｃ）１０Ｈｚ以上の高周波域での特性を改善するために、２次ローパスフィルタを２段追加するとともに、既存の２次ローパスフィルタの特性も変更する。

なお、（ｂ）で用いた特性補正フィルタは、固有周期を同じくする２つの２次減衰振動系の比のラプラス変換表現に対し、z変換を適用したものである。また、特許文献１の１次４段のフィルタは、等価な２次２段のフィルタに変換した。上記をまとめると、総合的なフィルタは下記のラプラス変換で表される８つのフィルタと最終的な出力をｇ倍するゲイン調整を直列に組み合わせたものとなる。

1/(ω₀ s^-1+1) (1)
ただし、ω₀=2πf₀
(ω₁s^-1+1)/(ω₁s^-1+2) (2)
(4ω₁s^-1+1)/(ω₁s^-1+8) (3)
(0.25ω₁s^-1+1)/(ω₁s^-1+0.5) (4)
ただし、ω₁=2πf₁
(1+2h_2aω₂ s^-1+ω₂ ² s^-2)/(1+2h_2bω₂ s^-1+ω₂ ² s^-2) (5)
ただし、ω₂=2πf₂
ω₃ ² s^-2/(1+2h₃ω₃ s^-1+ω₃ ² s^-2) (6)
ただし、ω₃=2πf₃
ω₄ ² s^-2/(1+2h₄ω₄ s^-1+ω₄ ² s^-2) (7)
ただし、ω₄=2πf₄
ω₅ ² s^-2/(1+2h₅ω₅ s^-1+ω₅ ² s^-2) (8)
ただし、ω₅=2πf₅である。

このとき、フィルタの特性を規定するパラメータは、（f₀，f₁，f₂，f₃，f₄，f₅，h_2a，h_2b，h₃，h₄，h₅，g）の12個となる。これらの値を本来の振幅特性との一致の度合いや、最終的な計算精度を参考にして決定したところ、近似の度合いの良いパラメータの組み合わせとして、（f₀＝0.45[Hz]，f₁＝7.0[Hz]，f₂＝0.5[Hz]，f₃＝12.0[Hz]，f₄＝20.0[Hz]，f₅＝30.0[Hz] ，h_2a＝1.0，h_2b＝0.75，h₃＝0.9，h₄＝0.6，h₅＝0.6，g＝1.262）を得た。このパラメータを用いたときの総合的な近似フィルタの振幅特性が図１０に示されている。

フィルタのラプラス変換表示が与えられれば、z変換を用いてデジタルフィルタを構成することができる。ΔＴをサンプリング間隔、zを遅延演算子として、
s^-1=(ΔT/2)(1+ z^-1)/(1- z^-1)
及び
s^-2=(ΔT/12^1/2)²(1+ 10z^-1+ z^-2)/(1- z^-1) ²
というz変換を用いると、ラプラス変換表示で以下のように表される。
E(s)＝(aω_a s^-1+1)/(ω_a s^-1+b)
ただし、ω_a=2πf_aである。

１次フィルタのz変換E(z)は、
s^-1=(ΔT/2)(1+ z^-1)/(1- z^-1)
を適用すれば、
E(z)=(β₀+z^-1β₁)/(α₀+z^-1α₁) (9)
α₀=ω_a +2b/(ΔT)
α₁=ω_a -2b/(ΔT)
β₀=ω_a a+2/(ΔT)
β₁=ω_a a-2/(ΔT)
ω_a=2πf_a
である。

ここでは、１次のデジタルフィルタを２つずつ組み合わせて、２次のフィルタとして再構成する。一般に、２つの一次デジタルフィルタ
E₁(z)=(β₀₁+z^-1β₁₁)/(α₀₁+z^-1α₁₁)
E₂(z)=(β₀₂+z^-1β₁₂)/(α₀₂+z^-1α₁₂)
の積 E₁(z)・E₂(z)は、２次のデジタルフィルタとなる。その係数は、
F(z) = E₁(z)・E₂(z)= (β₀+z^-1β₁+z^-2β₂)/(α₀+z^-1α₁+z^-2α₂) (10)
α₀=α₀₁α₀₂
α₁=α₀₁α₁₂+α₀₂α₁₁
α₂=α₁₁α₁₂
β₀=β₀₁β₀₂
β₁=β₀₁β₁₂+β₀₂β₁₁
β₂=β₁₁β₁₂
と表される。

以上を用いれば（１），（２）式のフィルタを組み合わせたデジタルフィルタとして、（９）式と（１０）式およびa=0.0,b=1.0, f_a=f_a1（（１）式に対応するフィルタに対して）、a=1.0,b=2.0, f_a=f_a2（（２）式に対応するフィルタに対して）より、
α₀=8/(ΔT)²+(4ω_a1＋2ω_a2 )/(ΔT)+ ω_a1ω_a2 (11)
α₁=2ω_a1ω_a2-16/(ΔT)²
α₂=8/(ΔT)²-(4ω_a1＋2ω_a2 )/(ΔT)+ ω_a1ω_a2
β₀=4/(ΔT)²+2ω_a2 /(ΔT)
β₁=-8/(ΔT)²
β₂=4/(ΔT)²-2ω_a2 /(ΔT)
ω_a1=2πf_a1, ω_a2=2πf_a2
が得られる。

また、（３），（４）式のフィルタを組み合わせたデジタルフィルタとして、（９）式と（１０）式およびa=4.0,b=8.0, f_a=f_a3（（３）式に対応するフィルタに対して）、a=0.25,b=0.5, f_a=f_a3（（４）式に対応するフィルタに対して）より、
α₀=16/(ΔT)²+17ω_a3/(ΔT)+ ω_a3 ² (12)
α₁=2ω_a3 ²-32/(ΔT)²
α₂=16/(ΔT)²-17ω_a3/(ΔT)+ ω_a3 ²
β₀=4/(ΔT)²+8.5ω_a3 /(ΔT) + ω_a3 ²
β₁=2ω_a3 ²-8/(ΔT)²
β₂=4/(ΔT)²-8.5ω_a3 /(ΔT) + ω_a3 ²
ω_a3=2πf_a3
が得られる。

（５）式の補正フィルタに関しては、ラプラス変換表示で以下のように表される。
G(s)=（1+2h_b1ω_b s^-1+ω_b ²s^-2）/(1+2h_b2ω_b s^-1+ω_b ²s^-2)
ただし、ω_b=2πf_bである。

２次フィルタのz変換G(z)として、
s^-1=(ΔT/2)(1+ z^-1)/(1- z^-1)
及び
s^-2=(ΔT/12^1/2)²(1+ 10z^-1+ z^-2)/(1- z^-1) ²
を適用すれば、
G(z)=(β₀+z^-1β₁+z^-2β₂)/(α₀+z^-1α₁+z^-2α₂) (13)
α₀=12/(ΔT)²+12h_b2ω_b /(ΔT)+ω_b ²
α₁=10ω_b ²-24/(ΔT)²
α₂=12/(ΔT)²-12h_b2ω_b /(ΔT)+ω_b ²
β₀=12/(ΔT)²+12h_b1ω_b /(ΔT)+ω_b ²
β₁=10ω_b ²-24/(ΔT)²
β₂=12/(ΔT)²-12h_b1ω_b /(ΔT)+ω_b ²
ω_b=2πf_b
がデジタルフィルタとして得られる。

また、（６）、（７）、（８）式の２次ローパスフィルタに関しては、ラプラス変換表示で以下のように表される。
H(s)=ω_c ²s^-2/(1+2h_cω_c s^-1+ω_c ²s^-2)
ただし，ω_c=2πf_cである。

２次フィルタのz変換H(z)として、
s^-1=(ΔT/2)(1+ z^-1)/(1- z^-1)
及び
s^-2=(ΔT/12^1/2)²(1+ 10z^-1+ z^-2)/(1- z^-1) ²
を適用すれば、
H(z)=(β₀+z^-1β₁+z^-2β₂)/(α₀+z^-1α₁+z^-2α₂) (14)
α₀=12/(ΔT)²+12h_cω_c /(ΔT)+ω_c ²
α₁=10ω_c ²-24/(ΔT)²
α₂=12/(ΔT)²-12h_cω_c /(ΔT)+ω_c ²
β₀=ω_c ²
β₁=10ω_c ²
β₂=ω_c ²
ω_c=2πf_c
がデジタルフィルタとして得られる。

入力時系列x(k)にこれらのフィルタを作用させた出力時系列y(k)を得るには、
y(k)=[-α₁y(k-1) -α₂y(k-2)+β₀x(k) +β₁x(k-1) +β₂x(k-2)]/α₀ (15)
なる演算を行えばよい。なお、ここでkは時系列のステップ数を表す。また、ゲイン調整はゲインを g_dとすれば、入力時系列をx(k)、出力時系列をy(k)として、
y(k)=g_dx(k) (16)
なる演算を行えばよい。

以上をまとめると、改良された近似フィルタ演算としては、加速度記録に対して、下記の７つのデジタルフィルタ演算を直列に行えばよい（演算Ａとする）。
１．f_a1＝f₀, f_a2＝f₁として、式（１１），（１５）の演算
２．f_a3＝f₁として、式（１２），（１５）の演算
３．h_b1 = h_2a, h_b2 = h_2b, f_b＝f₂ として、式（１３），（１５）の演算
４．h_c = h₃, f_c＝f₃ として、式（１４），（１５）の演算
５．h_c = h₄, f_c＝f₄ として、式（１４），（１５）の演算
６．h_c = h₅, f_c＝f₅ として、式（１４），（１５）の演算
７．g_d＝ g として、式（１６）の演算

１から７の演算は、地動加速度の時系列が得られるごとに行うことができる。従って、最終的なフィルタ処理波形も地動加速度の時系列が得られるごと(たとえば0.01秒ごと)に得ることができる。ここでフィルタ処理の入力は地動加速度である必要があるが、地震計の特性補正や微分演算等を施して地動加速度に変換することができれば、必ずしも加速度計により得られたものでなくてもよい。

１から７の演算はフィルタ係数の合成等によってより効率的に行うことも可能である。

また、演算Ａには他のフィルタ演算を追加してもよい。

図１０は、フィルタの振幅周波数特性を示す図である。図１０において、実線Ａは、本実施形態のフィルタの振幅周波数特性、破線Ｂは気象庁告知に基づくフィルタの振幅周波数特性、実線Ｃは本実施形態のフィルタの振幅周波数特性を気象庁告知に基づくフィルタの振幅周波数特性で除した相対的な振幅周波数特性、をそれぞれ表している。

本実施形態のフィルタは、0.1〜50Hzの帯域で気象庁告知に基づくフィルタとの差違が0.974倍から1.029倍に収まっている。

図１１は、本実施形態のフィルタを用いて計算されたリアルタイム震度（Ir）と気象庁震度（Ijma）の比較を示す図である。また、図１２及び図１３は、（ΔI＝Ijma−Ir）のヒストグラムを示す図である。

図１１に示すように、本実施形態のフィルタを用いて計算されたリアルタイム震度は、図１５に示した特許文献１に記載されたフィルタを用いて計算されたリアルタイム震度と比較して、気象庁震度にかなり近い結果であることがわかる。

全データを対象としたΔIの平均値は−0.0055、標準偏差は0.0272である。震度４以上（Ijma＞＝3.495）のデータを対象にすると、ΔIの平均値は−0.0011、標準偏差は0.0270である。図１２及び図１３からわかるように、ΔIの分布の左右非対称性は緩和されている。

本実施形態のフィルタによる標準偏差は、特許文献１に記載されたフィルタによる標準偏差に対し、全データを対象とした場合と震度４以上のデータを対象とした場合に、それぞれ0.70倍、0.63倍と改善されており、近似フィルタ改良の効果が期待通りに現れている。

表１は、本実施形態のフィルタのΔIが一定の値以内に収まるデータ数とその率を震度域に応じて計算した表である。

本実施形態のフィルタを用いた計算では、ΔIの絶対値が0.1以内となる率がすべての震度域において99％を超えているのがわかる。さらに、本実施形態のフィルタを用いた計算では、震度5弱以上の領域でΔIの絶対値がすべて0.15以内に収まっている。通常、計測震度が0.1刻みの値で扱われることを考えると、本実施形態のフィルタを用いた計算の精度は、速報用途としては十分に高い。

図１４は、本実施形態の計測震度概算装置をシステムとして構築した場合を示す。図中、１０１は計測震度概算システム、１０２は地震計、１０３は通信ネットワーク、１０４は地震データ処理装置である。このシステムでは、複数の地震計１０２で観測された地震動を、通信ネットワーク１０３を通じ他の機器やデータセンターに転送し、離れた場所で計測震度概算部７を有する地震データ処理装置１０４による概算の震度値算出を行うものである。ここで、演算処理を行う地震データ処理装置は複数存在してもよい。このようなシステムにおいてもいち早く震度の概算値を得られるという本方法の利点が活かされる。データセンターでは多くの観測点のデータを同時に処理しなければならない場合が多いため、FFTを用いない本方法のアルゴリズムは計算負荷を軽減する上でも有用である。

１…計測震度概算装置、２…処理部、３…制御・演算部、４…計測震度算出部、５…ＳＩ値算出部、６…応答値算出部、７…計測震度概算部、７１…地動加速度時系列取得手段、７２…時間領域フィルタ手段、７３…算出手段、７４…振幅時系列変換手段、７５…震度換算振幅時系列換算手段、７６…離散化手段、７７…継続時間カウント手段、８…電源装置、９…計測装置、１０１…計測震度概算システム、１０２…地震計、１０３…通信ネットワーク、１０４…地震データ処理装置

Claims (2)

1. 地震の加速度を計測し、計測した加速度から計測震度を概算する計測震度概算装置において、
   地動加速度の時系列を得る地動加速度時系列取得手段と、
   前記地動加速度の時系列をフィルタ処理する時間領域フィルタ手段と、
   前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列から計測震度の概算値を算出する算出手段と、
   を備え、
   前記算出手段は、
   前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列を振幅時系列に変換する振幅時系列変換手段と、
   前記振幅時系列を震度換算振幅時系列に換算する震度換算振幅時系列換算手段と、
   震度換算振幅時系列換算手段により求めた震度換算振幅時系列を離散化する離散化手段と、
   前記離散化手段により求めた各離散化震度換算振幅値に継続時間カウントを実行する継続時間カウント手段と、
   を有し、
   前記時間領域フィルタ手段は、以下の式（１１）乃至（１６）において、演算１乃至演算７を実行する
   ことを特徴とする計測震度概算装置。
   
   ・フィルタ1
   y(k)=[-α₁y(k-1) -α₂y(k-2)+β₀x(k) +β₁x(k-1) +β₂x(k-2)]/α₀ (15)
   ただし、
   α₀=8/(ΔT)²+(4ω_a1＋2ω_a2 )/(ΔT)+ω_a1ω_a2 (11)
   α₁=2ω_a1ω_a2-16/(ΔT)²
   α₂=8/(ΔT)²-(4ω_a1＋2ω_a2 )/(ΔT)+ω_a1ω_a2
   β₀=4/(ΔT)²+2ω_a2 /(ΔT)
   β₁=-8/(ΔT)²
   β₂=4/(ΔT)²-2ω_a2 /(ΔT)
   ω_a1=2πf_a1, ω_a2=2πf_a2
   x(k)はフィルタの入力時系列
   y(k)はフィルタの出力時系列
   kは時間ステップ数
   ΔTはサンプリング間隔
   πは円周率
   である。
   
   ・フィルタ2
   y(k)=[-α₁y(k-1) -α₂y(k-2)+β₀x(k) +β₁x(k-1) +β₂x(k-2)]/α₀ (15)
   ただし、
   α₀=16/(ΔT)²+17ω_a3/(ΔT)+ ω_a3 ² (12)
   α₁=2ω_a3 ²-32/(ΔT)²
   α₂=16/(ΔT)²-17ω_a3/(ΔT)+ ω_a3 ²
   β₀=4/(ΔT)²+8.5ω_a3 /(ΔT) + ω_a3 ²
   β₁=2ω_a3 ²-8/(ΔT)²
   β₂=4/(ΔT)²-8.5ω_a3 /(ΔT) + ω_a3 ²
   ω_a3=2πf_a3
   である。
   
   ・フィルタ3
   y(k)=[-α₁y(k-1) -α₂y(k-2)+β₀x(k) +β₁x(k-1) +β₂x(k-2)]/α₀ (15)
   ただし、
   α₀=12/(ΔT)²+12h_b2ω_b /(ΔT)+ω_b ² (13)
   α₁=10ω_b ²-24/(ΔT)²
   α₂=12/(ΔT)²-12h_b2ω_b /(ΔT)+ω_b ²
   β₀=12/(ΔT)²+12h_b1ω_b /(ΔT)+ω_b ²
   β₁=10ω_b ²-24/(ΔT)²
   β₂=12/(ΔT)²-12h_b1ω_b /(ΔT)+ω_b ²
   ω_b=2πf_b
   である。
   
   ・フィルタ4
   y(k)=[-α₁y(k-1) -α₂y(k-2)+β₀x(k) +β₁x(k-1) +β₂x(k-2)]/α₀ (15)
   ただし、
   α₀=12/(ΔT)²+12h_cω_c /(ΔT)+ω_c ² (14)
   α₁=10ω_c ²-24/(ΔT)²
   α₂=12/(ΔT)²-12h_cω_c /(ΔT)+ω_c ²
   β₀=ω_c ²
   β₁=10ω_c ²
   β₂=ω_c ²
   ω_c=2πf_c
   である。
   
   ・フィルタ5
   y(k)=g_d x(k) (16)
   ただし、g_d はゲインである。
   
   演算1： f_a1＝f₀, f_a2＝f₁として，式(11)，(15)の演算
   演算2： f_a3＝f₁として，式(12)，(15)の演算
   演算3： h_b1 = h_2a, h_b2 = h_2b, f_b＝f₂ として，式(13)，(15)の演算
   演算4： h_c = h₃, f_c＝f₃ として，式(14)，(15)の演算
   演算5： h_c = h₄, f_c＝f₄ として，式(14)，(15)の演算
   演算6： h_c = h₅, f_c＝f₅ として，式(14)，(15)の演算
   演算7： g_d＝ g として，式(16)の演算
   ただし、
   f₀はフィルタのカットオフ周波数を規定するパラメータ
   f₁は0.5次で減衰する周波数領域の位置を規定するパラメータ
   f₂、h_2a、h_2bは補正フィルタの特性を規定するパラメータ
   f₃、h₃はフィルタのカットオフ周波数およびダンピングを規定するパラメータ
   f₄、h₄はフィルタのカットオフ周波数およびダンピングを規定するパラメータ
   f₅、h₅はフィルタのカットオフ周波数およびダンピングを規定するパラメータ
   gはゲイン調整パラメータ
   である。
2. 地震動を計測する複数の地震計と、通信ネットワークと、前記複数の地震計で観測された地震動を、前記通信ネットワークを通じ、離れた場所で前記計測震度概算装置による概算の震度値算出を実行することを特徴とする請求項１に記載された計測震度概算装置を用いた計測震度概算システム。

Images