(57)【要約】【課題】即時性があり、換算の誤差の小さい計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法を提供することを目的とする。【解決手段】地震の加速度を計測し、計測した加速度から計測震度を概算する計測震度概算装置において、地動加速度の時系列を得る地動加速度時系列取得手段７１と、前記地動加速度の時系列をフィルタ処理する時間領域フィルタ手段７２と、前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列から計測震度の概算値を算出する算出手段７３とを備えたことを特徴とする。【選択図】図１
2008107334.tif 000002
【特許請求の範囲】
【請求項１】
地震の加速度を計測し、計測した加速度から計測震度を概算する計測震度概算装置において、地動加速度の時系列を得る地動加速度時系列取得手段と、前記地動加速度の時系列をフィルタ処理する時間領域フィルタ手段と、前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列から計測震度の概算値を算出する算出手段とを備えたことを特徴とする計測震度概算装置。
【請求項２】
前記算出手段は、前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列を振幅時系列に変換する振幅時系列変換手段と、前記振幅時系列を震度換算振幅時系列に換算する震度換算振幅時系列換算手段と、震度換算振幅時系列換算手段により求めた震度換算振幅時系列を離散化する離散化手段と、前記離散化手段により求めた各離散化震度換算振幅値に継続時間カウントを実行する継続時間カウント手段と、を有することを特徴とする請求項１に記載された計測震度概算装置。
【請求項３】
前記時間領域フィルタ手段は、下記の式（１１）、（１２）及び（１４−１）乃至（１６）において、演算１乃至演算６を実行することを特徴とする請求項１又は請求項２に記載された計測震度概算装置。
y(k)=-α1/α0×y(k-1)+β0/α0×x(k)+β1/α0×x(k-1)         ・・・（１２）
【数１】
2008107334.tif 000010
x(k)はフィルタの入力時系列
y(k)はフィルタの出力時系列
kは時間ステップ数
ωa=2×π×fa
ΔTはサンプリング間隔
a、b、faは周波数特性を規定するパラメータ
πは円周率
y(k)=-α1/α0×y(k-1)-α2/α0×y(k-2)+β0/α0×x(k)
+β1/α0×x(k-1)+β2/α0×x(k-2)   ・・・（１５）
【数２】
2008107334.tif 000011
ωb=2×π×fb
fb、hbは周波波数特性を規定するパラメータ
y(k)=gc×x(k)                                              ・・・（１６）
ただし、gcはゲイン
演算1：a=0.0、    b=1.0、    fa＝f0  として、式(11),(12)の演算
演算2：a=1.0、    b=2.0、    fa＝f1  として、式(11),(12)の演算
演算3：a=4.0、    b=8.0、    fa＝f1  として、式(11),(12)の演算
演算4：a=0.25、   b=0.5、    fa＝f1  として、式(11),(12)の演算
演算5：hb = h、    fb＝f2             として、式(14-1),(15)の演算
演算6：gc＝g                         として、式(16)の演算
ただし、f0はフィルタのカットオフ周波数を規定するパラメータ、
f1は0.5次で減衰する周波数領域の位置を規定するパラメータ、
f2はフィルタのカットオフ周波数を規定するパラメータ、
hはダンピングを規定するパラメータ、
gはゲイン調整パラメータ
【請求項４】
前記時間領域フィルタ手段は、下記の式（１１）、（１２）及び（１４−２）乃至（１６）において、演算１乃至演算６を実行することを特徴とする請求項１又は請求項２に記載された計測震度概算装置。
y(k)=-α1/α0×y(k-1)+β0/α0×x(k) +β1/α0×x(k-1)          ・・・（１２）
【数３】
2008107334.tif 000012
x(k)はフィルタの入力時系列
y(k)はフィルタの出力時系列
kは時間ステップ数
ωa=2×π×fa
ΔTはサンプリング間隔
a、b、faは周波数特性を規定するパラメータ
πは円周率
y(k)=-α1/α0×y(k-1)-α2/α0×y(k-2)+β0/α0×x(k)
+β1/α0×x(k-1)+β2/α0×x(k-2)    ・・・（１５）
【数４】
2008107334.tif 000013
ωb=2×π×fb
fb、hbは周波波数特性を規定するパラメータ
y(k)=gc×x(k)                                               ・・・（１６）
ただし、gcはゲイン
演算1：a=0.0、    b=1.0、    fa＝f0  として、式(11),(12)の演算
演算2：a=1.0、    b=2.0、    fa＝f1  として、式(11),(12)の演算
演算3：a=4.0、    b=8.0、    fa＝f1  として、式(11),(12)の演算
演算4：a=0.25、   b=0.5、    fa＝f1  として、式(11),(12)の演算
演算5：hb = h、    fb＝f2             として、式(14-2),(15)の演算
演算6：gc＝g                        として、式(16)の演算
ただし、f0はフィルタのカットオフ周波数を規定するパラメータ、
f1は0.5次で減衰する周波数領域の位置を規定するパラメータ、
f2はフィルタのカットオフ周波数を規定するパラメータ、
hはダンピングを規定するパラメータ、
gはゲイン調整パラメータ
【請求項５】
地震動を計測する複数の地震計と、通信ネットワークと、前記複数の地震計で観測された地震動を、前記通信ネットワークを通じ、離れた場所で前記計測震度概算装置による概算の震度値算出を実行することを特徴とする請求項１乃至請求項４のいずれか１項に記載された計測震度概算装置を用いた計測震度概算システム。
【請求項６】
地震の加速度を計測し、計測した加速度から計測震度を概算する計測震度概算方法において、地動加速度時系列を取得するステップと、前記地動加速度時系列を時間領域でフィルタ処理し、フィルタ処理地動加速度時系列を取得するステップと、前記フィルタ処理地動加速度の時系列から計測震度の概算値を算出するステップとからなることを特徴とする計測震度概算方法。
【請求項７】
前記フィルタ処理地動加速度の時系列から計測震度の概算値を算出するステップは、前記フィルタ処理地動加速度時系列から振幅時系列を取得するステップと、前記振幅時系列を震度換算振幅時系列に換算するステップと、前記震度換算振幅時系列を離散化し、各離散化震度換算振幅値を取得するステップと、前記各離散化震度換算振幅値について継続時間をカウントするステップと、前記各離散化震度換算振幅値と継続時間のカウント数から計測震度概算値を算出するステップとからなることを特徴とする請求項６に記載された計測震度概算方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、地震の加速度を計測し、その計測値から震度を概算する計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法に関する。
【背景技術】
【０００２】
従来から、地震が発生した際の警報発報や鉄道等の制御の指標として、気象庁が官報第１８３１号気象庁告示第４号気象業務法施行規則に関する告示で示した計測震度が用いられているのが現状である。計測震度を算出するためには、地震記録を一定時間蓄積してから、ＦＦＴを用いてフィルタ処理をすることが必要なため、刻々と測定される地震記録に対して値を算出することができず、即時性が要求される用途に適していなかった。
【０００３】
また、最大加速度、最大速度、ＳＩ値などの地震動指標を測定し、その地震動指標に基づいて簡単な演算処理をすることで震度相当の値を推定する方法がある（特許文献１参照）。
【０００４】
【特許文献１】特開２００３−２５５０５２号公報
【発明の開示】
【発明が解決しようとする課題】
【０００５】
ところが、この方法は、即時性があるが、換算の誤差が大きく、また換算式設計に際し、過去の大量の地震データを使用した解析を必要とする。
【０００６】
本発明は、前記課題を解決するために、即時性があり、換算の誤差の小さい計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法を提供することを目的とする。
【課題を解決するための手段】
【０００７】
前記課題を解決するために、本発明の計測震度概算装置は、地震の加速度を計測し、計測した加速度から計測震度を概算する計測震度概算装置において、地動加速度の時系列を得る地動加速度時系列取得手段と、前記地動加速度の時系列をフィルタ処理する時間領域フィルタ手段と、前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列から計測震度の概算値を算出する算出手段とを備えたことを特徴とする。
【０００８】
また、前記算出手段は、前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列を振幅時系列に変換する振幅時系列変換手段と、前記振幅時系列を震度換算振幅時系列に換算する震度換算振幅時系列換算手段と、震度換算振幅時系列換算手段により求めた震度換算振幅時系列を離散化する離散化手段と、前記離散化手段により求めた各離散化震度換算振幅値に継続時間カウントを実行する継続時間カウント手段と、を有することを特徴とする。
【０００９】
前記時間領域フィルタ手段は、下記の式（１１）、（１２）及び（１４―１）乃至（１６）において、演算１乃至演算6を実行することを特徴とする。
【００１０】
y(k)=-α1/α0×y(k-1)+β0/α0×x(k) +β1/α0×x(k-1)          ・・・（１２）
【００１１】
【数５】
2008107334.tif 000003
x(k)はフィルタの入力時系列
y(k)はフィルタの出力時系列
kは時間ステップ数
ωa=2×π×fa
ΔTはサンプリング間隔
a、b、faは周波数特性を規定するパラメータ
πは円周率
【００１２】
y(k)=-α1/α0×y(k-1)-α2/α0×y(k-2)+β0/α0×x(k)
+β1/α0×x(k-1)+β2/α0×x(k-2)     ・・・（１５）
【００１３】
【数６】
2008107334.tif 000004
ωb=2×π×fb
fb、hbは周波波数特性を規定するパラメータ
【００１４】
y(k)=gc×x(k)                                                ・・・（１６）
ただし、gcはゲイン
【００１５】
演算1：a=0.0、    b=1.0、    fa＝f0  として、式(11),(12)の演算
演算2：a=1.0、    b=2.0、    fa＝f1  として、式(11),(12)の演算
演算3：a=4.0、    b=8.0、    fa＝f1  として、式(11),(12)の演算
演算4：a=0.25、   b=0.5、    fa＝f1  として、式(11),(12)の演算
演算5：hb = h、    fb＝f2             として、式(14-1),(15)の演算
演算6：gc＝g                        として、式(16)の演算
ただし、f0はフィルタのカットオフ周波数を規定するパラメータ、
f1は0.5次で減衰する周波数領域の位置を規定するパラメータ、
f2はフィルタのカットオフ周波数を規定するパラメータ、
hはダンピングを規定するパラメータ、
gはゲイン調整パラメータ
【００１６】
また、前記時間領域フィルタ手段は、下記の式（１１）、（１２）及び（１４―２）乃至（１６）において、演算１乃至演算6を実行することを特徴とする。
【００１７】
y(k)=-α1/α0×y(k-1)+β0/α0×x(k) +β1/α0×x(k-1)          ・・・（１２）
【００１８】
【数７】
2008107334.tif 000005
x(k)はフィルタの入力時系列
y(k)はフィルタの出力時系列
kは時間ステップ数
ωa=2×π×fa
ΔTはサンプリング間隔
a、b、faは周波数特性を規定するパラメータ
πは円周率
【００１９】
y(k)=-α1/α0×y(k-1) -α2/α0×y(k-2)+β0/α0×x(k) +β1/α0×x(k-1) +β2/α0×x(k-2)
・・・（１５）
【００２０】
【数８】
2008107334.tif 000006
ωb=2×π×fb
fb、hbは周波波数特性を規定するパラメータ
【００２１】
y(k)=gc×x(k)                                                ・・・（１６）
ただし、gcはゲイン
【００２２】
演算1：a=0.0、    b=1.0、    fa＝f0  として、式(11),(12)の演算
演算2：a=1.0、    b=2.0、    fa＝f1  として、式(11),(12)の演算
演算3：a=4.0、    b=8.0、    fa＝f1  として、式(11),(12)の演算
演算4：a=0.25、   b=0.5、    fa＝f1  として、式(11),(12)の演算
演算5：hb = h、    fb＝f2             として、式(14-2),(15)の演算
演算6：gc＝g                        として、式(16)の演算
ただし、f0はフィルタのカットオフ周波数を規定するパラメータ、
f1は0.5次で減衰する周波数領域の位置を規定するパラメータ、
f2はフィルタのカットオフ周波数を規定するパラメータ、
hはダンピングを規定するパラメータ、
gはゲイン調整パラメータ
【００２３】
また、本発明の計測震度概算装置を用いた計測震度概算システムは、地震動を計測する複数の地震計と、通信ネットワークと、前記複数の地震計で観測された地震動を、前記通信ネットワークを通じ、離れた場所で前記計測震度概算装置による概算の震度値算出を実行することを特徴とする。
【００２４】
さらに、本発明の計測震度概算方法は、地震の加速度を計測し、計測した加速度から計測震度を概算する計測震度概算方法において、地動加速度時系列を取得するステップと、前記地動加速度時系列を時間領域でフィルタ処理し、フィルタ処理地動加速度時系列を取得するステップと、前記フィルタ処理地動加速度の時系列から計測震度の概算値を算出するステップとからなることを特徴とする。
【００２５】
また、前記フィルタ処理地動加速度の時系列から計測震度の概算値を算出するステップは、前記フィルタ処理地動加速度時系列から振幅時系列を取得するステップと、前記振幅時系列を震度換算振幅時系列に換算するステップと、前記震度換算振幅時系列を離散化し、各離散化震度換算振幅値を取得するステップと、前記各離散化震度換算振幅値について継続時間をカウントするステップと、前記各離散化震度換算振幅値と継続時間のカウント数から計測震度概算値を算出するステップとからなることを特徴とする。
【発明の効果】
【００２６】
このような計測震度概算装置及び計測震度概算方法により、即時性があり、換算の誤差の小さい計測震度の概算をすることができる。
【発明を実施するための最良の形態】
【００２７】
本発明の実施の形態を図により説明する。図１は、本発明の実施形態の計測震度概算装置の主要構成を示す図である。図中、１は計測震度概算装置、２は処理部、３は制御・演算部、４は計測震度算出部、５はＳＩ値算出部、６は応答値算出部、７は計測震度概算部、８は電源装置、９は計測装置である。
【００２８】
計測震度概算装置１は、計測装置９で計測した加速度等を制御・演算部３に入力することで、震度を概算する。処理部２は、電源装置８により駆動される制御・演算部３を有する。本実施形態の制御・演算部３は、よりよい精度を求めるため、現行の計測震度計に計測震度概算部７を追加した形態としたが、即時性をより重視する場合等では、計測震度算出部４、ＳＩ値算出部５及び応答値算出部６を備えなくてもよい。また、計測装置９を処理部２内に備えてもよい。
【００２９】
図２は、本発明の実施形態の計測震度概算部７の主要構成を示す図である。図中、７１は地動加速度時系列取得手段、７２は時間領域フィルタ手段、７３は算出手段、７４は振幅時系列変換手段、７５は震度換算振幅時系列換算手段、７６は離散化手段、７７は継続時間カウント手段である。
【００３０】
地動加速度時系列取得手段７１は直交３成分の地動加速度の時系列を取得し、時間領域フィルタ手段７２に送る。時間領域フィルタ手段７２は地動加速度時系列取得手段７１が取得した地動加速度の時系列をフィルタ処理する。算出手段７３は時間領域フィルタ手段７２によりフィルタ処理された地動加速度の時系列から計測震度の概算値を算出する。
【００３１】
算出手段７３は、振幅時系列変換手段７４、震度換算振幅時系列換算手段７５、離散化手段７６及び継続時間カウント手段７７を有する。振幅時系列変換手段７４は時間領域フィルタ手段７２によりフィルタ処理された地動加速度の時系列を振幅時系列に変換し、震度換算振幅時系列換算手段７５は振幅時系列を震度換算振幅時系列に換算する。離散化手段７６は震度換算振幅時系列換算手段により求めた震度換算振幅時系列を離散化し、継続時間カウント手段７７は離散化手段により求めた各離散化震度換算振幅値に継続時間カウントを実行する。
【００３２】
次に、図３に示すフローチャートを用いて、このような計測震度概算装置１の計測震度概算方法について説明する。まず、ステップ１で、地震動を観測する（ＳＴ１）。
【００３３】
次に、ステップ２で、時間間隔dt（例えば0.01秒間隔）でサンプリングされた、直交３成分の地動加速度時系列xa[k]、ya[k]、za[k]を得る。ここで、kは時間ステップ数とする（ＳＴ２）。本方法の入力として必要なのは、地動加速度の時系列であるが、地震計の特性補正や速度の微分演算等を施して地動加速度に変換することができれば、必ずしも加速度計により得られたものでなくてもよい。また、地震動の観測は東西方向と南北方向等の水平方向２成分と鉛直方向成分についてそれぞれ行い３成分の地動加速度の時系列を得る。ただし、本方法はいずれかの２成分もしくは１成分の地動加速度の時系列に関しても適用可能である。地震動の観測が水平方向２成分と鉛直方向成分でなされていない場合は、必要な変換を行い３成分の地動加速度の時系列を得る。なお、この変換は、次に示すステップ３のフィルタ処理の後に実行してもよい。
【００３４】
次に、ステップ３で、直交３成分の地動加速度時系列（xa[k]、ya[k]、za[k]）に後述する演算Ａで定めるフィルタ処理を行い、３成分のフィルタ処理された地動加速度の時系列（xf[k]、yf[k]、zf[k]）を得る（ＳＴ３）。なお、フィルタの特性に関しては、後述する。また、地震観測の成分数が3成分に満たない場合は、存在しない（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とする。
【００３５】
次に、ステップ４で、フィルタ処理された3成分地動加速度の時系列を
af[k]=(xf[k]×xf[k]+yf[k]×yf[k]+zf[k]×zf[k])1/2             ・・・（１）
なる演算で振幅時系列af[k]に変換する（ＳＴ４）。
【００３６】
気象庁告示第4号の本来の定義に従えば、この振幅時系列から、一定の時間区間内において積算継続時間が0.3秒以上となる最大の振幅レベルa0を求め、換算式I=2×log(a0)+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が求まり、図３のフローチャートに戻る。
【００４８】
また、図７を用いて第２の継続時間のカウント方法についての基本的な考え方を説明する。図7は(a)は、サンプリング周波数を10Ｈｚ、積算継続時間をカウントする時間区間を1秒として、積算継続時間をカウントする時間区間内の震度換算振幅値を値の大きい順にならびかえたものである。なお、図中で上段の数字は、震度換算振幅値の値を、下段の数字はその値が得られた、時間ステップ数を表している。この場合、積算継続時間が0.3×10Hz=3となる震度換算振幅値は、上から３番目に大きい値６となる。
【００４９】
次に、図７（ｂ）に示すように、時間ステップ数が１つ進んだ場合、図7（ａ）における時間ステップ数0の離散化震度換算振幅値が積算継続時間をカウントする時間区間から外れる。そして、図7（ｃ）に示すように、新たに時間ステップ数10の離散化震度換算振幅値が積算継続時間をカウントする時間区間に入る。この場合、積算継続時間が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が求まり、図３のフローチャートに戻る。
【００５７】
なお、第2の継続時間のカウント方法の場合、カウントに関する処理は離散化震度換算振幅値の大小比較による。よって、整数値化された離散化震度換算振幅値ではなく、震度換算振幅値そのものを用いて処理することもできる。
【００５８】
また、継続時間のカウント方法は、同様の効果が得られるのであれば、第１、第２の方法以外であってもよい。さらに、データの保持は、配列でなくてもよい。
【００５９】
サブルーチンを終了した後、図３のフローチャートに戻り、最後に、ステップ８で、計測震度の概算値を算出する（ＳＴ８）。pの値に、離散化幅dIを乗じ下限値Iiを加えた値がその時点での計測震度の概算値Ia[k]となる。
【００６０】
以上の処理を、時間ステップ毎に行う。したがって、3成分の地動加速度（xa[k]、ya[k]、za[k]）を得る毎にその時点での震度の概算値Ia[k]を得ることができる。図９に示すように、本方法で求められる震度の概算値Ia[k]は、地震動が継続している間に時間変化する時系列として得られる。この時間変化する値の最大値Ixを、地震動が継続している区間に対して一つの値が決まる本来の定義の計測震度の概算値と定める。地震動が継続している区間を定義しがたい場合は、例えば地震検知後1分間の最大値のように定めればよい。このようにIa[k]は時々刻々と時間変化することから、最大値Ixに達する前に警戒レベルを超えた場合に即時警報を発することも可能になる。
【００６１】
なお、前述した各ステップの順序は、これに限定されるものではなく、計測震度概算値が変わらない限り、ステップの順序を変更してもよい。
【００６２】
次に、気象庁の告示した計測震度本来の定義のフィルタを近似した、時間領域フィルタの設計法について述べる。
【００６３】
計測震度本来のフィルタ特性は、官報第１８３１号気象庁告示第４号気象業務法施行規則に関する告示に示されている。図１０は、このフィルタの地動加速度（計測の単位はｃｍ／ｓ2)を入力とする振幅特性を図示したものである。この振幅特性の特徴は、領域Ｆ１では周波数が下がるに従って1次で減衰し、領域Ｆ２では周波数が上がるに従って０．５次で減衰し、領域Ｆ３では周波数が上がるに従って次第に減衰が強くなることである。近似フィルタ設計においては、この特徴を演算量が少ない再帰型フィルタを複数段用いて実現する。
【００６４】
１次の減衰を持つ再帰型フィルタとしては、例えば式（３）のラプラス変換で表されるフィルタが知られている。
【００６５】
1/(ω0×s-1+1)                                               ・・・(３)
ただし、ω0=2×π×f0
ここで、f0はフィルタのカットオフ周波数を規定するパラメータである。
【００６６】
一定の区間で0.5次の減衰を持ちその他の区間で平坦な特性をもつ再帰型フィルタとしては、例えば式（４）のラプラス変換で表されるフィルタが知られている。
【００６７】
(ω1×s-1+1)/(ω1×s-1+2)×(4×ω1×s-1+1)/(ω1×s-1+8)
×(0.25×ω1×s-1+1)/(ω1×s-1+0.5)     ・・・（４）
ただし、ω1=2×π×f1
ここで、f1は0.5次で減衰する周波数領域の位置を規定するパラメータである。
【００６８】
なお、式（４）で表されるフィルタは、
(ω1×s-1+1)/(ω1×s-1+2)                                     ・・・（５）
(4×ω1×s-1+1)/(ω1×s-1+8)                                  ・・・（６）
(0.25×ω1×s-1+1)/(ω1×s-1+0.5)                             ・・・（７）
と表される3つのフィルタを順次作用させることに等しい。
【００６９】
図１０において、領域Ｆ３におけるフィルタの特性は不要な高周波成分を除去するために必要なものである。通常の強震観測では３０Ｈｚ乃至４０Ｈｚより高い周波数は地震動の数値化の段階で除去されているため、この部分の近似精度は最終的な演算結果への影響が少ない。このためこの部分は２次で減衰するものとして近似する。この特性を実現する再帰型フィルタとしては、例えば式（８）のラプラス変換で表されるフィルタが知られている。
【００７０】
(ω2×ω2×s-2)/[1+2×h×ω2×s-1+(ω2×ω2×s-2)]              ・・・（８）
ω2=2×π×f2
ここで、f2はフィルタのカットオフ周波数、hはダンピングを規定するパラメータである。なお、f2を必要に応じて、f2＝tan(πf2ΔT)/ΔT/πと置き換えてもよい。
【００７１】
これら、式（３）、（５）、（６）、（７）及び（８）のフィルタとゲイン調整として最終的な出力をｇ倍する処理を組み合わせれば近似フィルタを構成できる。このとき、フィルタの特性をコントロールするパラメータは（ｆ0、ｆ1、ｆ2、ｈ、ｇ）の５個となる。これらの値は、本来の振幅特性との近似の度合いをグラフ上で比較することや、最小自乗法等を用いることで決定することができる。実施例ではグラフ上での比較から、近似の度合いの良い組み合わせの例として、
（f0＝0.45[Hz]、f1＝7.0[Hz]、f2＝11.0[Hz]、h＝0.9、g＝1.409） ・・・（９）
を得た。図１０には、このパラメータを用いたときの総合的な近似フィルタの振幅特性と本来の振幅特性が示されている。
【００７２】
再帰型フィルタのラプラス変換表示が与えられれば、ｚ変換の理論を用いて計算機での処理に必要となるデジタルフィルタを構成することができる。
【００７３】
ラプラス変換形式でF(s)=(a×ωa×s-1+1)/(ωa×s-1+b)、ωa=2×π×fa、と表される一次フィルタのz変換F(z)は、ΔTをサンプリング間隔、zを遅延演算子として、
双一次変換、s-1=(ΔT/2)×(1+ z-1)/( 1- z-1)を適用すれば
F(z)=(β0+z-1×β1)/(α0+z-1×α1)                        ・・・（１０）
【００７４】
【数９】
2008107334.tif 000007
である。
【００７５】
入力時系列x(k)にフィルタF(z)を作用させた出力時系列y(k)を得るには、
y(k)=-α1/α0×y(k-1)+β0/α0×x(k) +β1/α0×x(k-1)     ・・・（１２）
なる演算を行えばよい。ここでkは時間ステップ数とする。
【００７６】
また、ラプラス変換形式でG(s)= (ωb×ωb×s-2)/[1+2×hb×ωb×s-1+(ωb×ωb×s-2)]、ωb=2×π×fb、と表される二次フィルタのz変換G(z)は、ΔTをサンプリング間隔、zを遅延演算子として、双一次変換、s-1=(ΔT/2)×(1+ z-1)/( 1- z-1)、およびs-2=[(ΔT/2)×(1+ z-1)/( 1- z-1)]×[(ΔT/2)×(1+ z-1)/( 1- z-1)]を適用すれば
G(z)=(β0+z-1×β1+z-2×β2)/(α0+z-1×α1+z-2×α2)       ・・・（１３）
【００７７】
【数１０】
2008107334.tif 000008
である。
【００７８】
また、ラプラス変換形式でG(s)= (ωb×ωb×s-2)/[1+2×hb×ωb×s-1+(ωb×ωb×s-2)]、ωb=2×π×fb、と表される二次フィルタのz変換G(z)は、ΔTをサンプリング間隔、zを遅延演算子として、変換、s-1=(ΔT/2)×(1+ z-1)/( 1- z-1)、およびs-2=ΔT×ΔT/12×(1+10×z-1+ z-2)×(1-z-1)/( 1- z-1)を適用すれば
G(z)=(β0+z-1×β1+z-2×β2)/(α0+z-1×α1+z-2×α2)      ・・・（１３）
【００７９】
【数１１】
2008107334.tif 000009
とあらわすこともできる。
【００８０】
入力時系列x(k)にフィルタG(z)を作用させた出力時系列y(k)を得るには、
y(k)=-α1/α0×y(k-1) -α2/α0×y(k-2)+β0/α0×x(k)
+β1/α0×x(k-1) +β2/α0×x(k-2) ・・・（１５）
なる演算を行えばよい。
【００８１】
また、ゲイン調整はゲインをgcとすれば、入力時系列をx(k)、出力時系列をy(k)として、
y(k)=gc×x(k)                                           ・・・（１６）
なる演算を行えばよい。
【００８２】
したがって、時間領域の近似フィルタ処理としては下記の6つの操作を直列に行えばよい（演算Ａとする）。
【００８３】
1     a=0.0、    b=1.0、   fa＝f0  として、式(11)、(12)の演算
2     a=1.0、    b=2.0、   fa＝f1  として、式(11)、(12)の演算
3     a=4.0、    b=8.0、   fa＝f1  として、式(11)、(12)の演算
4     a=0.25、   b=0.5、   fa＝f1  として、式(11)、(12)の演算
5     hb = h、   fb＝f2           として、式(14-1)又は式(14-2)、(15)の演算
6     gc＝g                        として、式(16)の演算
【００８４】
1から6の演算は、地動加速度の時系列が得られるごとに行うことができる。従って、最終的なフィルタ処理波形も地動加速度の時系列が得られるごと(たとえば0.01秒ごと)に得ることができる。ここでフィルタ処理の入力は地動加速度である必要があるが、地震計の特性補正や微分演算等を施して地動加速度に変換することができれば、必ずしも加速度計により得られたものでなくてもよい。
【００８５】
１から６の演算はフィルタ係数の合成等によってより効率的に行うことも可能である。
【００８６】
また、演算Ａには他のフィルタ演算を追加してもよい。
【００８７】
次に、図１１に示すグラフから本実施形態の精度を検証する。本方法の精度を検証するために、独立行政法人防災科学技術研究所が運用する全国強震観測網（K-NET）で公開されている96370件の強震データを対象に、本来の方法に基づいて計算した計測震度Ijmaと、本方法で計算された計測震度の概算値（振動継続期間中最大値）Ixを比較した。図１１に示すように、標準偏差は0.061であり、SI値、最大加速度、最大速度等の指標から換算する方法よりも精度が高い。
【００８８】
このように、ＦＦＴを用いずに、時間領域のフィルタを用いることで、即時性があり、ある程度の精度を保持する計測震度を概算することができる。
【００８９】
図１２は本実施形態の計測震度概算装置をシステムとして構築した場合を示す。図中、１０１は計測震度概算システム、１０２は地震計、１０３は通信ネットワーク、１０４は地震データ処理装置である。このシステムでは、複数の地震計１０２で観測された地震動を、通信ネットワーク１０３を通じ他の機器やデータセンターに転送し、離れた場所で計測震度概算部７を有する地震データ処理装置１０４による概算の震度値算出を行うものである。ここで、演算処理を行う地震データ処理装置は複数存在してもよい。このようなシステムにおいてもいち早く震度の概算値を得られるという本方法の利点が活かされる。データセンターでは多くの観測点のデータを同時に処理しなければならない場合が多いため、FFTを用いない本方法のアルゴリズムは計算負荷を軽減する上でも有用である。
【図面の簡単な説明】
【００９０】
【図１】本実施形態の計測震度概算装置のブロック図である。
【図２】本実施形態の計測震度概算部のブロック図である。
【図３】本実施形態の計測震度概算値を求めるフローチャート図である。
【図４】本実施形態の離散化を示す図である。
【図５】本実施形態の継続時間カウントの考え方を示す図である。
【図６】本実施形態の継続時間カウントのフローチャート図である。
【図７】第２の継続時間カウントの考え方を示す図である。
【図８】第２の継続時間カウントのフローチャート図である。
【図９】本実施形態の計測震度概算値を示す図である。
【図１０】本実施形態のフィルタ特性を示す図である。
【図１１】本来の計測震度値と本実施形態の計測震度概算値との比較を示す図である。
【図１２】本実施形態の計測震度概算システムを示す図である。
【符号の説明】
【００９１】
１…計測震度概算装置、２…処理部、３…制御・演算部、４…計測震度算出部、５…ＳＩ値算出部、６…応答値算出部、７…計測震度概算部、７１…地動加速度時系列取得手段、７２…時間領域フィルタ手段、７３…算出手段、７４…振幅時系列変換手段、７５…震度換算振幅時系列換算手段、７６…離散化手段、７７…継続時間カウント手段、８…電源装置、９…計測装置、１０１…計測震度概算システム、１０２…地震計、１０３…通信ネットワーク、１０４…地震データ処理装置
【図１】
2008107334.tif 000014
【図２】
2008107334.tif 000015
【図３】
2008107334.tif 000016
【図４】
2008107334.tif 000017
【図５】
2008107334.tif 000018
【図６】
2008107334.tif 000019
【図７】
2008107334.tif 000020
【図８】
2008107334.tif 000021
【図９】
2008107334.tif 000022
【図１０】
2008107334.tif 000023
【図１１】
2008107334.tif 000024
【図１２】
2008107334.tif 000025