(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2025-05-15
(45)【発行日】2025-05-23
(54)【発明の名称】時間領域フィルタ装置、計測震度概算装置、計測震度概算システム及び時間領域フィルタ方法
(51)【国際特許分類】
G01V 1/28 20060101AFI20250516BHJP
【FI】
G01V1/28
(21)【出願番号】P 2022044852
(22)【出願日】2022-03-22
【審査請求日】2024-07-18
(73)【特許権者】
【識別番号】501138231
【氏名又は名称】国立研究開発法人防災科学技術研究所
(74)【代理人】
【識別番号】100139103
【氏名又は名称】小山 卓志
(74)【代理人】
【識別番号】100139114
【氏名又は名称】田中 貞嗣
(74)【代理人】
【識別番号】100214260
【氏名又は名称】相羽 昌孝
(72)【発明者】
【氏名】▲功▼刀 卓
(72)【発明者】
【氏名】青井 真
(72)【発明者】
【氏名】中村 洋光
【審査官】鴨志田 健太
(56)【参考文献】
【文献】特開2014-228373(JP,A)
【文献】特許第4229337(JP,B1)
【文献】特開2011-047657(JP,A)
【文献】特開2014-181915(JP,A)
(58)【調査した分野】(Int.Cl.,DB名)
G01V 1/28
(57)【特許請求の範囲】
【請求項1】
地動加速度の時系列をフィルタ処理する時間領域フィルタ装置であって、
以下の式(11)において、フィルタ処理する第1フィルタと、
以下の式(12)において、フィルタ処理する第2フィルタと、
以下の式(13)において、フィルタ処理する第3フィルタと、
以下の式(14)において、フィルタ処理する第4フィルタと、
以下の式(15)において、入力時系列x(k)にこれらのフィルタを作用させて出力時系列y(k)を取得する出力時系列取得部と、
以下の式(16)において、ゲインを調整するゲイン調整部と、
演算11乃至演算17を実行する演算部と、
を備える
ことを特徴とする時間領域フィルタ装置。
・第1フィルタ
α
0=8/(ΔT)
2+(4ω
a1+2ω
a2 )/(ΔT)+ω
a1ω
a2 (11)
α
1=2ω
a1ω
a2-16/(ΔT)
2
α
2=8/(ΔT)
2-(4ω
a1+2ω
a2 )/(ΔT)+ω
a1ω
a2
β
0=4/(ΔT)
2+2ω
a2 /(ΔT)
β
1=-8/(ΔT)
2
β
2=4/(ΔT)
2-2ω
a2 /(ΔT)
ω
a1=2πf
a1, ω
a2=2πf
a2
・第2フィルタ
α
0=16/(ΔT)
2+17ω
a3/(ΔT)+ ω
a3
2
(12)
α
1=2ω
a3
2-32/(ΔT)
2
α
2=16/(ΔT)
2-17ω
a3/(ΔT)+ ω
a3
2
β
0=4/(ΔT)
2+8.5ω
a3 /(ΔT) + ω
a3
2
β
1=2ω
a3
2-8/(ΔT)
2
β
2=4/(ΔT)
2-8.5ω
a3 /(ΔT) + ω
a3
2
ω
a3=2πf
a3
・第3フィルタ
α
0=1+ h
b2ω
bΔT +γ
b2ω
b
2 (ΔT)
2 (13)
α
1=-2+ω
b
2(ΔT)
2(1-2γ
b2)
α
2=1- h
b2ω
bΔT +γ
b2ω
b
2 (ΔT)
2
β
0=1+ h
b1ω
bΔT +γ
b1ω
b
2 (ΔT)
2
β
1=-2+ω
b
2(ΔT)
2(1-2γ
b1)
β
2=1- h
b1ω
bΔT +γ
b1ω
b
2 (ΔT)
2
ω
b=2πf
b
・第4フィルタ
α
0=1+ h
cω
cΔT +γ
c2ω
c
2 (ΔT)
2 (14)
α
1=-2+ω
c
2(ΔT)
2(1-2γ
c2)
α
2=1- h
cω
cΔT +γ
c2ω
c
2 (ΔT)
2
β
0=γ
c1ω
c
2(ΔT)
2
β
1=ω
c
2(ΔT)
2(1-2γ
c1)
β
2=γ
c1ω
c
2(ΔT)
2
ω
c=2πf
c
・出力時系列取得部
y(k)=[-α
1y(k-1) -α
2y(k-2)+β
0x(k) +β
1x(k-1) +β
2x(k-2)]/α
0 (15)
ただし、
x(k)はフィルタの入力時系列、
y(k)はフィルタの出力時系列、
kは時間ステップ数、
ΔTはサンプリング間隔、
πは円周率、
とする。
・ゲイン調整部
y(k)=g
d x(k) (16)
ただし、g
d はゲインである。
・演算部
演算11:f
a1 = f
0, f
a2 = f
1として、式(11),(15)の演算
演算12:f
a3 = f
1として、式(12),(15)の演算
演算13:h
b1= h
2a, h
b2= h
2b, f
b= f
2, γ
b1= γ
2a, γ
b2= γ
2b として、式(13),(15)の演算
演算14:h
c= h
3, f
c= f
3, γ
c1= γ
3a, γ
c2= γ
3b として、式(14),(15)の演算
演算15:h
c = h
4, f
c = f
4, γ
c1 = γ
4a, γ
c2 = γ
4b として、式(14),(15)の演算
演算16:h
c = h
5, f
c = f
5, γ
c1 = γ
5a, γ
c2 = γ
5b として、式(14),(15)の演算
演算17:g
d = g として、式(16)の演算
ただし、
f
0はフィルタのカットオフ周波数を規定するパラメータ、
f
1は0.5次で減衰する周波数領域の位置を規定するパラメータ、
f
2、h
2a、h
2b、は補正フィルタの特性を規定するパラメータ、
f
3、h
3はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
f
4、h
4はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
f
5、h
5はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
gはゲイン調整パラメータ、
γ
2a、γ
2b、γ
3a、γ
3b, γ
4a、γ
4b、γ
5a、γ
5b はフィルタの二次積分の近似方法を
規定するパラメータ、
であり、
下記の式(21)乃至(24)を満たす。
(1/4-1/(2πf
2ΔT)
2)≦γ
2b (21)
(1/4-1/(2πf
3ΔT)
2)≦γ
3b (22)
(1/4-1/(2πf
4ΔT)
2)≦γ
4b (23)
(1/4-1/(2πf
5ΔT)
2)≦γ
5b (24)
【請求項2】
前記フィルタの二次積分の近似方法を規定するパラメータにおいて、下記の式(25)乃至(28)のうち少なくとも一つを満たす
ことを特徴とする請求項1に記載された時間領域フィルタ装置。
γ
2a = γ
2b (25)
γ
3a = γ
3b (26)
γ
4a = γ
4b (27)
γ
5a = γ
5b (28)
【請求項3】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(1)を満たす
ことを特徴とする請求項1または請求項2に記載された時間領域フィルタ装置。
条件(1)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz]の場合、
f
s≧ 80[Hz]のとき
γ
2b = 1/12, γ
3b = 1/12, γ
4b = 1/12, γ
5b = 1/12
80[Hz] > f
s≧ 70[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/12, γ
5b = 1/8
70[Hz] > f
s≧ 60[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/12, γ
5b = 1/6
60[Hz] > f
s≧ 50[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/10, γ
5b = 1/4
50[Hz] > f
s≧ 40[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/6, γ
5b = 1/4
40[Hz] > f
s≧ 30[Hz]のとき
γ
2b= 1/12, γ
3b = 1/10, γ
4b = 1/4, γ
5b = 1/4
30[Hz] > f
s≧ 5[Hz]のとき
γ
2b= 1/12, γ
3b = 1/4, γ
4b = 1/4, γ
5b = 1/4
5[Hz] > f
s≧ 1[Hz]のとき
γ
2b= 1/6, γ
3b = 1/4, γ
4b = 1/4, γ
5b = 1/4
【請求項4】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(2)を満たす
ことを特徴とする請求項1または請求項2に記載された時間領域フィルタ装置。
条件(2)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz] の場合、
f
s ≧ 80[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/12, γ
5b = 1/12
80[Hz] > f
s≧ 70[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/12, γ
5b = 1/8
70[Hz] > f
s≧ 60[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/12, γ
5b = 1/6
60[Hz] > f
s≧ 50[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/8, γ
5b = 1/4
50[Hz] > f
s≧ 40[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/6, γ
5b = 1/4
40[Hz] > f
s≧ 30[Hz]のとき
γ
2b= 1/12, γ
3b = 1/8, γ
4b = 1/4, γ
5b = 1/4
30[Hz] > f
s≧ 5[Hz]のとき
γ
2b= 1/12, γ
3b = 1/4, γ
4b = 1/4, γ
5b = 1/4
5[Hz] > f
s≧ 1[Hz]のとき
γ
2b= 1/6, γ
3b = 1/4, γ
4b = 1/4, γ
5b = 1/4
【請求項5】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(3)を満たす
ことを特徴とする請求項1または請求項2に記載された時間領域フィルタ装置。
条件(3)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz] の場合、
f
s ≧ 80[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/12, γ
5b = 1/12
80[Hz] > f
s≧ 60[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/12, γ
5b = 1/6
60[Hz] > f
s≧ 40[Hz]のとき
γ
2b= 1/12, γ
3b = 1/12, γ
4b = 1/6, γ
5b = 1/4
40[Hz] > f
s≧ 30[Hz]のとき
γ
2b= 1/12, γ
3b = 1/6, γ
4b = 1/4, γ
5b = 1/4
30[Hz] > f
s≧ 5[Hz]のとき
γ
2b= 1/12, γ
3b = 1/4, γ
4b = 1/4, γ
5b = 1/4
5[Hz] > f
s≧ 1[Hz]のとき
γ
2b= 1/6, γ
3b = 1/4, γ
4b = 1/4, γ
5b = 1/4
【請求項6】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(4)を満たす
ことを特徴とする請求項1または請求項2に記載された時間領域フィルタ装置。
条件(4)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz] の場合、
f
s ≧ 80[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/12
80[Hz] > f
s ≧ 60[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/4
60[Hz] > f
s ≧ 40[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/4, γ
5b = 1/4
40[Hz] > f
s ≧ 5[Hz]のとき
γ
2b = 1/12, γ
3b= 1/4, γ
4b = 1/4, γ
5b = 1/4
5[Hz] > f
s ≧ 1[Hz]のとき
γ
2b = 1/4, γ
3b = 1/4, γ
4b= 1/4, γ
5b = 1/4
【請求項7】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(5)を満たす
ことを特徴とする請求項1または請求項2に記載された時間領域フィルタ装置。
条件(5)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz] の場合、
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/12
【請求項8】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(6)を満たす
ことを特徴とする請求項1または請求項2に記載された時間領域フィルタ装置。
条件(6)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz] の場合、
γ
2b = 1/4,
γ
3b = 1/4,
γ
4b= 1/4,
γ
5b = 1/4
【請求項9】
前記サンプリング間隔ΔTをステップ毎に異なる値とする
ことを特徴とする請求項1乃至請求項8のいずれか1項に記載された時間領域フィルタ装置。
【請求項10】
地動加速度の時系列を取得する地動加速度時系列取得部と、
前記地動加速度時系列取得部が取得した前記地動加速度の時系列をフィルタ処理する請求項1乃至9のいずれか1項に記載された時間領域フィルタ装置と、
前記時間領域フィルタ装置によりフィルタ処理された前記地動加速度の時系列から計測震度の概算値を算出する算出部と、
を備え、
前記算出部は、
前記時間領域フィルタ装置によりフィルタ処理された前記地動加速度の時系列を振幅時系列に変換する振幅時系列変換部と、
前記振幅時系列変換部が変換した振幅時系列を震度換算振幅時系列に換算する震度換算振幅時系列換算部と、
前記震度換算振幅時系列換算部により求めた震度換算振幅時系列を離散化する離散化部と、
前記離散化部により求めた各離散化震度換算振幅値に継続時間カウントを実行する継続時間カウント部と、
を有する
ことを特徴とする計測震度概算装置。
【請求項11】
請求項10に記載された計測震度概算装置と、
地震動を計測する複数の地震計と、
通信ネットワークと、
前記複数の地震計で観測された地震動を、前記通信ネットワークを通じ、離れた場所で前記計測震度概算装置による概算の震度値算出を実行する地震データ処理装置と、
を備える
ことを特徴とする計測震度概算システム。
【請求項12】
地動加速度の時系列をフィルタ処理する時間領域フィルタ方法であって、
以下の式(11)乃至(16)において、演算11乃至演算17の少なくとも一つを実行する
ことを特徴とする時間領域フィルタ方法。
・第1フィルタ
α
0=8/(ΔT)
2+(4ω
a1+2ω
a2 )/(ΔT)+ω
a1ω
a2 (11)
α
1=2ω
a1ω
a2-16/(ΔT)
2
α
2=8/(ΔT)
2-(4ω
a1+2ω
a2 )/(ΔT)+ω
a1ω
a2
β
0=4/(ΔT)
2+2ω
a2 /(ΔT)
β
1=-8/(ΔT)
2
β
2=4/(ΔT)
2-2ω
a2 /(ΔT)
ω
a1=2πf
a1, ω
a2=2πf
a2
・第2フィルタ
α
0=16/(ΔT)
2+17ω
a3/(ΔT)+ ω
a3
2
(12)
α
1=2ω
a3
2-32/(ΔT)
2
α
2=16/(ΔT)
2-17ω
a3/(ΔT)+ ω
a3
2
β
0=4/(ΔT)
2+8.5ω
a3 /(ΔT) + ω
a3
2
β
1=2ω
a3
2-8/(ΔT)
2
β
2=4/(ΔT)
2-8.5ω
a3 /(ΔT) + ω
a3
2
ω
a3=2πf
a3
・第3フィルタ
α
0=1+ h
b2ω
bΔT +γ
b2ω
b
2 (ΔT)
2 (13)
α
1=-2+ω
b
2(ΔT)
2(1-2γ
b2)
α
2=1- h
b2ω
bΔT +γ
b2ω
b
2 (ΔT)
2
β
0=1+ h
b1ω
bΔT +γ
b1ω
b
2 (ΔT)
2
β
1=-2+ω
b
2(ΔT)
2(1-2γ
b1)
β
2=1- h
b1ω
bΔT +γ
b1ω
b
2 (ΔT)
2
ω
b=2πf
b
・第4フィルタ
α
0=1+ h
cω
cΔT +γ
c2ω
c
2 (ΔT)
2 (14)
α
1=-2+ω
c
2(ΔT)
2(1-2γ
c2)
α
2=1- h
cω
cΔT +γ
c2ω
c
2 (ΔT)
2
β
0=γ
c1ω
c
2(ΔT)
2
β
1=ω
c
2(ΔT)
2(1-2γ
c1)
β
2=γ
c1ω
c
2(ΔT)
2
ω
c=2πf
c
・入力時系列x(k)にこれらのフィルタを作用させた出力時系列y(k)
y(k)=[-α
1y(k-1) -α
2y(k-2)+β
0x(k) +β
1x(k-1) +β
2x(k-2)]/α
0 (15)
ただし、
x(k)はフィルタの入力時系列、
y(k)はフィルタの出力時系列、
kは時間ステップ数、
ΔTはサンプリング間隔、
πは円周率、
とする。
・ゲイン調整
y(k)=g
d x(k) (16)
ただし、g
d はゲインである。
演算11:f
a1 = f
0, f
a2 = f
1として、式(11),(15)の演算
演算12:f
a3 = f
1として、式(12),(15)の演算
演算13:h
b1= h
2a, h
b2= h
2b, f
b= f
2, γ
b1= γ
2a, γ
b2= γ
2b として、式(13),(15)の演算
演算14:h
c= h
3, f
c= f
3, γ
c1= γ
3a, γ
c2= γ
3b として、式(14),(15)の演算
演算15:h
c = h
4, f
c = f
4, γ
c1 = γ
4a, γ
c2 = γ
4b として、式(14),(15)の演算
演算16:h
c = h
5, f
c = f
5, γ
c1 = γ
5a, γ
c2 = γ
5b として、式(14),(15)の演算
演算17:g
d = g として、式(16)の演算
ただし、
f
0はフィルタのカットオフ周波数を規定するパラメータ、
f
1は0.5次で減衰する周波数領域の位置を規定するパラメータ、
f
2、h
2a、h
2b、は補正フィルタの特性を規定するパラメータ、
f
3、h
3はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
f
4、h
4はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
f
5、h
5はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
gはゲイン調整パラメータ、
γ
2a、γ
2b、γ
3a、γ
3b, γ
4a、γ
4b、γ
5a、γ
5b はフィルタの二次積分の近似方法を規定するパラメータ、
であり、
γ
2b、γ
3b、γ
4b、γ
5b は、以下の式(21)乃至(24)を満たす。
(1/4-1/(2πf
2ΔT)
2)≦γ
2b (21)
(1/4-1/(2πf
3ΔT)
2)≦γ
3b (22)
(1/4-1/(2πf
4ΔT)
2)≦γ
4b (23)
(1/4-1/(2πf
5ΔT)
2)≦γ
5b (24)
【請求項13】
前記フィルタの二次積分の近似方法を規定するパラメータにおいて、以下の式(25)乃至(28)のうち少なくとも一つを満たす
ことを特徴とする請求項12に記載された時間領域フィルタ方法。
γ
2a = γ
2b (25)
γ
3a = γ
3b (26)
γ
4a = γ
4b (27)
γ
5a = γ
5b (28)
【請求項14】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(1)を満たす
ことを特徴とする請求項12又は請求項13に記載された時間領域フィルタ方法。
条件(1)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz]の場合、
f
s≧ 80[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/12
80[Hz] > f
s ≧ 70[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/8
70[Hz] > f
s ≧ 60[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/6
60[Hz] > f
s ≧ 50[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/10, γ
5b = 1/4
50[Hz] > f
s ≧ 40[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/6, γ
5b = 1/4
40[Hz] > f
s ≧ 30[Hz]のとき
γ
2b = 1/12, γ
3b= 1/10, γ
4b = 1/4, γ
5b = 1/4
30[Hz] > f
s ≧ 5[Hz]のとき
γ
2b = 1/12, γ
3b= 1/4, γ
4b = 1/4, γ
5b = 1/4
5[Hz] > f
s ≧ 1[Hz]のとき
γ
2b = 1/6, γ
3b = 1/4, γ
4b= 1/4, γ
5b = 1/4
【請求項15】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(2)を満たす
ことを特徴とする請求項12又は請求項13に記載された時間領域フィルタ方法。
条件(2)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz] の場合、
f
s ≧ 80[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/12
80[Hz] > f
s ≧ 70[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/8
70[Hz] > f
s ≧ 60[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/6
60[Hz] > f
s ≧ 50[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/8, γ
5b = 1/4
50[Hz] > f
s ≧ 40[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/6, γ
5b = 1/4
40[Hz] > f
s ≧ 30[Hz]のとき
γ
2b = 1/12, γ
3b= 1/8, γ
4b = 1/4, γ
5b = 1/4
30[Hz] > f
s ≧ 5[Hz]のとき
γ
2b = 1/12, γ
3b= 1/4, γ
4b = 1/4, γ
5b = 1/4
5[Hz] > f
s ≧ 1[Hz]のとき
γ
2b = 1/6, γ
3b = 1/4, γ
4b= 1/4, γ
5b = 1/4
【請求項16】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(3)を満たす
ことを特徴とする請求項12又は請求項13に記載された時間領域フィルタ方法。
条件(3)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz] の場合、
f
s ≧ 80[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/12
80[Hz] > f
s ≧ 60[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/6
60[Hz] > f
s ≧ 40[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/6, γ
5b = 1/4
40[Hz] > f
s ≧ 30[Hz]のとき
γ
2b = 1/12, γ
3b= 1/6, γ
4b = 1/4, γ
5b = 1/4
30[Hz] > f
s ≧ 5[Hz]のとき
γ
2b = 1/12, γ
3b= 1/4, γ
4b = 1/4, γ
5b = 1/4
5[Hz] > f
s ≧ 1[Hz]のとき
γ
2b = 1/6, γ
3b = 1/4, γ
4b= 1/4, γ
5b = 1/4
【請求項17】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(4)を満たす
ことを特徴とする請求項12又は請求項13に記載された時間領域フィルタ方法。
条件(4)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz] の場合、
f
s ≧ 80[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/12
80[Hz] > f
s ≧ 60[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/4
60[Hz] > f
s ≧ 40[Hz]のとき
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/4, γ
5b = 1/4
40[Hz] > f
s ≧ 5[Hz]のとき
γ
2b = 1/12, γ
3b= 1/4, γ
4b = 1/4, γ
5b = 1/4
5[Hz] > f
s ≧ 1[Hz]のとき
γ
2b = 1/4, γ
3b = 1/4, γ
4b= 1/4, γ
5b = 1/4
【請求項18】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(5)を満たす
ことを特徴とする請求項12又は請求項13に記載された時間領域フィルタ方法。
条件(5)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz] の場合、
γ
2b = 1/12, γ
3b= 1/12, γ
4b = 1/12, γ
5b = 1/12
【請求項19】
前記フィルタの二次積分の近似方法を規定するパラメータ、前記フィルタの特性を規定するパラメータ、前記フィルタのカットオフ周波数を規定するパラメータにおいて、以下の条件(6)を満たす
ことを特徴とする請求項12又は請求項13に記載された時間領域フィルタ方法。
条件(6)
サンプリング周波数をf
s = 1/ΔTとして、
f
2=0.5[Hz],f
3=12.0[Hz],f
4=20.0[Hz],f
5=30.0[Hz] の場合、
γ
2b = 1/4,
γ
3b = 1/4,
γ
4b = 1/4,
γ
5b= 1/4
【請求項20】
前記サンプリング間隔ΔTをステップ毎に異なる値とする
ことを特徴とする請求項12乃至請求項19のいずれか1項に記載された時間領域フィルタ方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、地震の加速度の時系列をフィルタ処理する時間領域フィルタ装置及び時間領域フィルタ方法、時間領域フィルタ装置を用いてフィルタ処理した地震の加速度の時系列から震度を概算する計測震度概算装置及び計測震度概算システムに関する。
【背景技術】
【0002】
従来から、地震が発生した際の警報発報や鉄道等の制御の指標として、気象庁が官報第1831号気象庁告示第4号気象業務法施行規則に関する告示で示した計測震度が用いられているのが現状である。計測震度を算出するためには、地震記録を一定時間蓄積してから、FFTを用いてフィルタ処理をすることが必要なため、刻々と測定される地震記録に対して値を算出することができず、即時性が要求される用途に適していなかった。
【0003】
計測震度は、現在最も広く認知されている強震動指標であるが、演算に周波数領域のフィルタ処理を必要とすることから速報性が求められる用途には適していない。これを解決するため、本発明者は、特許文献1で、周波数領域でのフィルタ演算を時間領域の近似フィルタで代用することによって実時間で震度の近似値を得るリアルタイム演算法を提案した。また、本出願人は、特許文献1の近似フィルタを改良することにより、即時性を維持しつつ、さらに換算の誤差の小さいリアルタイム演算法を特許文献2で提案した。
【0004】
特許文献2に記載されたフィルタは、主要な帯域0.1Hz~10Hzで気象庁告示に基づくフィルタとの差違が0.974倍から1.029倍に収まるように構成されている。これは相対的な振幅周波数特性である実線Cから読み取ることができる。このように、特許文献2に記載されたフィルタは、即時性があり、換算誤差の小さい計測震度の概算をすることができた。
【先行技術文献】
【特許文献】
【0005】
【文献】特許第4229337号公報
【文献】特許第5946067号公報
【発明の概要】
【発明が解決しようとする課題】
【0006】
しかしながら、従来の計測震度の概算は、低いサンプリング周波数(約77Hz以下)において、安定したフィルタを構成できず、100Hzのサンプリングに補間した上で、近似フィルタを用いていた。
【0007】
本発明は、低いサンプリング周波数においても、補間を行うことなく、そのまま処理できる時間領域フィルタ装置及び時間領域フィルタ方法、並びに、時間領域フィルタ装置を用いて、即時性を維持しつつ、換算の誤差の小さい計測震度概算装置及び計測震度概算システムを提供することを目的とする。
【課題を解決するための手段】
【0008】
前記課題を解決するために、本発明の時間領域フィルタ装置は、
地動加速度の時系列をフィルタ処理する時間領域フィルタ装置であって、
以下の式(11)において、フィルタ処理する第1フィルタと、
以下の式(12)において、フィルタ処理する第2フィルタと、
以下の式(13)において、フィルタ処理する第3フィルタと、
以下の式(14)において、フィルタ処理する第4フィルタと、
以下の式(15)において、入力時系列x(k)にこれらのフィルタを作用させて出力時系列y(k)を取得する出力時系列取得部と、
以下の式(16)において、ゲインを調整するゲイン調整部と、
演算11乃至演算17を実行する演算部と、
を備える。
・第1フィルタ
α0=8/(ΔT)2+(4ωa1+2ωa2 )/(ΔT)+ωa1ωa2 (11)
α1=2ωa1ωa2-16/(ΔT)2
α2=8/(ΔT)2-(4ωa1+2ωa2 )/(ΔT)+ωa1ωa2
β0=4/(ΔT)2+2ωa2 /(ΔT)
β1=-8/(ΔT)2
β2=4/(ΔT)2-2ωa2 /(ΔT)
ωa1=2πfa1, ωa2=2πfa2
・第2フィルタ
α0=16/(ΔT)2+17ωa3/(ΔT)+ ωa3
2
(12)
α1=2ωa3
2-32/(ΔT)2
α2=16/(ΔT)2-17ωa3/(ΔT)+ ωa3
2
β0=4/(ΔT)2+8.5ωa3 /(ΔT) + ωa3
2
β1=2ωa3
2-8/(ΔT)2
β2=4/(ΔT)2-8.5ωa3 /(ΔT) + ωa3
2
ωa3=2πfa3
・第3フィルタ
α0=1+ hb2ωbΔT +γb2ωb
2 (ΔT)2 (13)
α1=-2+ωb
2(ΔT)2(1-2γb2)
α2=1- hb2ωbΔT +γb2ωb
2 (ΔT)2
β0=1+ hb1ωbΔT +γb1ωb
2 (ΔT)2
β1=-2+ωb
2(ΔT)2(1-2γb1)
β2=1- hb1ωbΔT +γb1ωb
2 (ΔT)2
ωb=2πfb
・第4フィルタ
α0=1+ hcωcΔT +γc2ωc
2 (ΔT)2 (14)
α1=-2+ωc
2(ΔT)2(1-2γc2)
α2=1- hcωcΔT +γc2ωc
2 (ΔT)2
β0=γc1ωc
2(ΔT)2
β1=ωc
2(ΔT)2(1-2γc1)
β2=γc1ωc
2(ΔT)2
ωc=2πfc
・出力時系列取得部
y(k)=[-α1y(k-1) -α2y(k-2)+β0x(k) +β1x(k-1) +β2x(k-2)]/α0 (15)
ただし、
x(k)はフィルタの入力時系列、
y(k)はフィルタの出力時系列、
kは時間ステップ数、
ΔTはサンプリング間隔、
πは円周率、
とする。
・ゲイン調整部
y(k)=gd x(k) (16)
ただし、gd はゲインである。
演算11:fa1 = f0, fa2 = f1として、式(11),(15)の演算
演算12:fa3 = f1として、式(12),(15)の演算
演算13:hb1= h2a, hb2= h2b, fb= f2, γb1= γ2a, γb2= γ2b として、式(13),(15)の演算
演算14:hc= h3, fc= f3, γc1= γ3a, γc2= γ3b として、式(14),(15)の演算
演算15:hc = h4, fc = f4, γc1 = γ4a, γc2 = γ4b として、式(14),(15)の演算
演算16:hc = h5, fc = f5, γc1 = γ5a, γc2 = γ5b として、式(14),(15)の演算
演算17:gd = g として、式(16)の演算
ただし、
f0はフィルタのカットオフ周波数を規定するパラメータ、
f1は0.5次で減衰する周波数領域の位置を規定するパラメータ、
f2、h2a、h2b、は補正フィルタの特性を規定するパラメータ、
f3、h3はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
f4、h4はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
f5、h5はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
gはゲイン調整パラメータ、
γ2a、γ2b、γ3a、γ3b, γ4a、γ4b、γ5a、γ5b はフィルタの二次積分の近似方法を
規定するパラメータ、
であり、
下記の式(21)乃至(24)を満たす。
(1/4-1/(2πf2ΔT)2)≦γ2b (21)
(1/4-1/(2πf3ΔT)2)≦γ3b (22)
(1/4-1/(2πf4ΔT)2)≦γ4b (23)
(1/4-1/(2πf5ΔT)2)≦γ5b (24)
【発明の効果】
【0009】
このような時間領域フィルタ装置及び時間領域フィルタ方法によれば、低いサンプリング周波数においても、補間を行うことなく、そのまま処理することできる。また、計測震度概算装置及び計測震度概算システムによれば、時間領域フィルタ装置を用いることで、即時性を維持しつつ、換算の誤差を小さくすることができる。
【図面の簡単な説明】
【0010】
【
図1】本実施形態の計測震度概算装置のブロック図である。
【
図2】本実施形態の計測震度概算部のブロック図である。
【
図3】本実施形態の計測震度概算値を求めるフローチャート図である。
【
図5】本実施形態の継続時間カウントの考え方を示す図である。
【
図6】本実施形態の継続時間カウントのフローチャート図である。
【
図7】第2の継続時間カウントの考え方を示す図である。
【
図8】第2の継続時間カウントのフローチャート図である。
【
図9】本実施形態の時間領域フィルタ装置のブロック図である。
【
図10】地震波形と従来の補間を行う方法の波形と本実施形態の計測震度概算装置を用いた波形を示す図である。
【
図11】本実施形態の計測震度概算システムを示す図である。
【発明を実施するための形態】
【0011】
本発明の実施の形態を図により説明する。
図1は、本発明の実施形態の計測震度概算装置の主要構成を示す図である。図中、1は計測震度概算装置、2は処理部、3は制御・演算部、4は計測震度算出部、5はSI値算出部、6は応答値算出部、7は計測震度概算部、8は電源装置、9は計測装置である。
【0012】
計測震度概算装置1は、計測装置9で計測した加速度等を制御・演算部3に入力することで、震度を概算する。処理部2は、電源装置8により駆動される制御・演算部3を有する。本実施形態の制御・演算部3は、よりよい精度を求めるため、現行の計測震度計に計測震度概算部7を追加した形態としたが、即時性をより重視する場合等では、計測震度算出部4、SI値算出部5及び応答値算出部6を備えなくてもよい。また、計測装置9を処理部2内に備えてもよい。
【0013】
図2は、本発明の実施形態の計測震度概算部7の主要構成を示す図である。図中、71は地動加速度時系列取得部、72は時間領域フィルタ装置、73は算出部、74は振幅時系列変換部、75は震度換算振幅時系列換算部、76は離散化部、77は継続時間カウント部である。
【0014】
地動加速度時系列取得部71は直交3成分の地動加速度の時系列を取得し、時間領域フィルタ装置72に送る。時間領域フィルタ装置72は地動加速度時系列取得部71が取得した地動加速度の時系列をフィルタ処理する。算出部73は時間領域フィルタ装置72によりフィルタ処理された地動加速度の時系列から計測震度の概算値を算出する。
【0015】
算出部73は、振幅時系列変換部74、震度換算振幅時系列換算部75、離散化部76及び継続時間カウント部77を有する。振幅時系列変換部74は時間領域フィルタ装置72によりフィルタ処理された地動加速度の時系列を振幅時系列に変換し、震度換算振幅時系列換算部75は振幅時系列を震度換算振幅時系列に換算する。離散化部76は震度換算振幅時系列換算部により求めた震度換算振幅時系列を離散化し、継続時間カウント部77は離散化部により求めた各離散化震度換算振幅値に継続時間カウントを実行する。
【0016】
図3は、本実施形態の計測震度概算値を求めるフローチャート図である。
図4は、本実施形態の離散化を示す図である。
【0017】
次に、
図3に示すフローチャートを用いて、このような計測震度概算装置1の計測震度概算方法について説明する。まず、ステップ1で、地震動を観測する(ST1)。
【0018】
次に、ステップ2で、時間間隔dt(例えば0.01秒間隔)でサンプリングされた、直交3成分の地動加速度時系列xa[k]、ya[k]、za[k]を得る。ここで、kは時間ステップ数とする(ST2)。本方法の入力として必要なのは、地動加速度の時系列であるが、地震計の特性補正や速度の微分演算等を施して地動加速度に変換することができれば、必ずしも加速度計により得られたものでなくてもよい。また、地震動の観測は東西方向と南北方向等の水平方向2成分と鉛直方向成分についてそれぞれ行い3成分の地動加速度の時系列を得る。ただし、本方法はいずれかの2成分もしくは1成分の地動加速度の時系列に関しても適用可能である。地震動の観測が水平方向2成分と鉛直方向成分でなされていない場合は、必要な変換を行い3成分の地動加速度の時系列を得る。なお、この変換は、次に示すステップ3のフィルタ処理の後に実行してもよい。
【0019】
次に、ステップ3で、直交3成分の地動加速度時系列(xa[k]、ya[k]、za[k])に後述する演算Aで定めるフィルタ処理を行い、3成分のフィルタ処理された地動加速度の時系列(xf[k]、yf[k]、zf[k])を得る(ST3)。なお、フィルタの特性に関しては、後述する。また、地震観測の成分数が3成分に満たない場合は、存在しない(xf[k]、yf[k] 、zf[k])を0として、ステップ3以降の処理を実行することができる。例えば、x,y,z,の順番は無関係なので、地震観測の成分数が2成分の場合は、xf[k]、yf[k] 、zf[k] のうちいずれかの1つを0、地震観測の成分数が1成分の場合は、xf[k]、yf[k] 、zf[k] のうちいずれかの2つを0とする。
【0020】
次に、ステップ4で、フィルタ処理された3成分地動加速度の時系列を
af[k]=(xf[k]×xf[k]+yf[k]×yf[k]+zf[k]×zf[k])1/2 ・・・(A)
なる演算で振幅時系列af[k]に変換する(ST4)。
【0021】
気象庁告示第4号の本来の定義に従えば、この振幅時系列から、一定の時間区間内において積算継続時間が0.3秒以上となる最大の振幅レベルa0を求め、換算式I=2×log(a0)+0.94で計測震度相当に換算することになるが、計算の効率化のため、計測震度相当に変換してから積算継続時間に関する処理を行う。
【0022】
まず、ステップ5で、ステップ4において得られた振幅時系列を
If[k]=2×log(af[k])+0.94 ・・・(B)
なる換算式で、震度換算した震度換算振幅時系列If[k]に変換する(ST5)。
【0023】
次に、ステップ6で、
図4に示すように、If[k]を必要な精度の離散化幅dIで離散化する(ST6)。本実施形態では、計測震度が0.1刻みで算出されるため、最終的な結果の十分な精度を保持するためdIは0.1の10分の1である0.01で離散化した。
【0024】
具体的には、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に相当する値の制限を行ったことになる。
【0025】
値の制限範囲は、求めようとする計測震度の範囲を含む必要がある。範囲の上限については、震度7の下限である計測震度6.5よりも十分に大きい値にする必要がある。実施例では、9.99とした。これは計測された加速度が10000galを上回る場合に相当し、地震波センサ自体の計測上限を超える値であるため、十分に大きい。下限については、震度0の上限である計測震度0.4よりも小さい値にする必要がある。実施例では-8とした。これは地震の無い時の微弱な地面の揺れを地震動とみなして計算された計測震度よりも小さい値として設定した値であり、十分に小さい値である。
【0026】
次に、ステップ7で、各離散化震度換算振幅値について継続時間のカウントを行う(ST7)。
【0027】
ここで、
図5を用いて第1の継続時間のカウント方法についての基本的な考え方を説明する。
図5は、時間ステップ数kと離散化震度換算振幅値との関係を示すもので、
図5(a)は積算継続時間のカウント範囲を時間ステップ数の2~4とした場合であり、
図5(b)は積算継続時間のカウント範囲を時間ステップ数の3~5とした場合である。この時、時間ステップ数2に対応する離散化震度換算振幅値はカウント範囲を外れ、時間ステップ数5に対応する離散化震度換算振幅値はカウント範囲に加わることとなる。従って、時間ステップ数3~4の離散化震度換算振幅値に関しては、時間ステップ数が1増えても変わらないので、そのまま残すようにする。
【0028】
このような考え方を用いて継続時間のカウントをサブルーチンとして実行する。なお、ここでいう継続時間とは、気象庁告示第4号の本来の定義から明らかであるが、離散化震度換算振幅値の時系列がある一定の振幅に等しい状態が継続した時間ではなく、ある一定の振幅以上となった状態が継続した時間をいう。
【0029】
図6は第1の継続時間のカウント方法を示すサブルーチンのフローチャートである。あらかじめ、N個の離散化震度換算振幅値の値ごとに継続時間カウント用の記憶領域を確保しておく。この配列は0に初期化しておき、この記憶領域をIc[i]とする。配列の添え字iは0からN-1とする。
【0030】
まず、ステップ101で、上記のId[k]、Ic[i]を用いて、継続時間のカウントを、Id[k]の値以下の添え字iをもつカウント用記憶領域Ic[i]の値をすべて+1することで実行する(ST101)。
【0031】
次に、ステップ102で、積算継続時間をカウントする時間区間範囲を外れたものについてカウントした値を元にもどす必要がある。これは、Id[k-m]の値以下の添え字iをもつIc[i]の値をすべて-1することで効率よく実行できる(ST102)。ここでmは積算継続時間をカウントする時間区間範囲に対応するサンプル数である。実施例の場合は100Hzサンプリングで積算継続時間をカウントする時間区間範囲を10秒とした。この場合mは1000となる。これを実行するためには、m個の記憶領域を確保し、常にId[k]からId[k-m]までのm個の値を記憶しておく必要がある。このカウント方法により、時間ステップ毎にm個の継続時間に関する処理を行わずにすむ。
【0032】
次に、ステップ103で、積算継続時間0.3秒に対応する総カウント数M(実施例のように100Hzサンプリングの場合は30)以上の値となるIc[i]のうち最大の添え字のものをIc[p]とすれば、pが求まり(ST103)、サブルーチンを終了する。このように、サブルーチンを実行することでpが求まり、
図3のフローチャートに戻る。
【0033】
また、
図7を用いて第2の継続時間のカウント方法についての基本的な考え方を説明する。
図7(a)は、サンプリング周波数を10Hz、積算継続時間をカウントする時間区間を1秒として、積算継続時間をカウントする時間区間内の震度換算振幅値を値の大きい順にならびかえたものである。なお、図中で上段の数字は、震度換算振幅値の値を、下段の数字はその値が得られた、時間ステップ数を表している。この場合、積算継続時間が0.3×10Hz=3となる震度換算振幅値は、上から3番目に大きい値6となる。
【0034】
次に、
図7(b)に示すように、時間ステップ数が1つ進んだ場合、
図7(a)における時間ステップ数0の離散化震度換算振幅値が積算継続時間をカウントする時間区間から外れる。そして、
図7(c)に示すように、新たに時間ステップ数10の離散化震度換算振幅値が積算継続時間をカウントする時間区間に入る。この場合、積算継続時間が0.3×10Hz=3となる震度換算振幅値は、上から3番目に大きい値5となる。
【0035】
図8は第2の継続時間のカウント方法を示すサブルーチンのフローチャートである。
【0036】
まず、ステップ201で、配列Is[j]の初期化が行われているか判断する(ST201)。
【0037】
ステップ201において、初期化が行われていないと判断された場合、ステップ202で、初期化を実行する(ST202)。初期化について説明する。まず、m個の記憶領域を持つ配列Is[j]を用意する。配列の添字jは1からmまでとする。ここで、mは
図6のステップ102で用いた第1の方法ものと同様に、積算継続時間をカウントする時間区間範囲に対応するサンプル数である。本第2の方法の場合は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]のどの要素がいつの時間ステップの値であるかを記憶しておく。これで初期化が完了する。
【0038】
次に、時間ステップ数kにおいて、積算継続時間をカウントする時間区間範囲から外れる値Id[k-m]と、新たに加わる値Id[k]に関する操作のみを行って、Is[j]を値の大きい順に並んだ配列に保つようにする。
【0039】
具体的には、ステップ201において、初期化が行われていると判断された場合、又は、ステップ202で初期化が行われた場合、ステップ203で、It[q]=k-mとなるqを見つけ、時間区間範囲から外れる値であるId[k-m]に対応しているIs[q]を配列Is[j]から削除する。削除した部分は詰めて配列Is[j]をm-1個の配列とする。配列It[n]も同様にIt[q]を削除し、削除した部分は詰める(ST203)。
【0040】
次に、ステップ204で、Id[k]の値と、配列Is[j]の値を比較し、Is[j]が大きい順にならぶような場所を見つけて挿入し、Is[j]を再びm個の配列とする。配列It[n]も同じ場所に要素を挿入し、その値はkとし追加された値に対応する時間ステップ数を記憶させる(ST204)。
【0041】
以上の操作が終了すると、Is[j]は大きい順にならんだ配列となっており、It[n]も更新されたIs[j]について、どの要素がいつの時間ステップの値であるかを記憶している。
【0042】
そして、ステップ205で、値の大きい方からM番目の値である積算継続時間0.3秒に対応する離散化震度換算振幅値Is[M]を取得し(ST205)、サブルーチンを終了する。Is[M]は
図6のステップ103で用いた第1の方法のpに対応する。また、Mは第1の方法と同様に、第2の方法の場合30である。このように、サブルーチンを実行することでpが求まり、
図3のフローチャートに戻る。
【0043】
なお、第2の継続時間のカウント方法の場合、カウントに関する処理は離散化震度換算振幅値の大小比較による。よって、整数値化された離散化震度換算振幅値ではなく、震度換算振幅値そのものを用いて処理することもできる。
【0044】
また、継続時間のカウント方法は、同様の効果が得られるのであれば、第1、第2の方法以外であってもよい。さらに、データの保持は、配列でなくてもよい。
【0045】
サブルーチンを終了した後、
図3のフローチャートに戻り、最後に、ステップ8で、計測震度の概算値を算出する(ST8)。pの値に、離散化幅dIを乗じ下限値Iiを加えた値がその時点での計測震度の概算値Ia[k]となる。
【0046】
以上の処理を、時間ステップ毎に行う。したがって、3成分の地動加速度(xa[k]、ya[k]、za[k])を得る毎にその時点での震度の概算値Ia[k]を得ることができる。
図9に示すように、本方法で求められる震度の概算値Ia[k]は、地震動が継続している間に時間変化する時系列として得られる。この時間変化する値の最大値Ixを、地震動が継続している区間に対して一つの値が決まる本来の定義の計測震度の概算値と定める。地震動が継続している区間を定義しがたい場合は、例えば地震検知後1分間の最大値のように定めればよい。このようにIa[k]は時々刻々と時間変化することから、最大値Ixに達する前に警戒レベルを超えた場合に即時警報を発することも可能になる。
【0047】
なお、前述した各ステップの順序は、これに限定されるものではなく、計測震度概算値が変わらない限り、ステップの順序を変更してもよい。
【0048】
次に、気象庁の告示した計測震度本来の定義のフィルタを近似した、時間領域フィルタの設計法について述べる。
【0049】
特許文献2で提案されたフィルタは下記のラプラス変換で表される8つのフィルタと最終的な出力をg倍するゲイン調整を直列に組み合わせたものであった。
【0050】
1/(ω0 s-1+1) (1’)
ただし、ω0=2πf0
(ω1s-1+1)/(ω1s-1+2) (2’)
(4ω1s-1+1)/(ω1s-1+8) (3’)
(0.25ω1s-1+1)/(ω1s-1+0.5) (4’)
ただし、ω1=2πf1
(1+2h2aω2 s-1+ω2
2s-2)/(1+2h2bω2 s-1+ω2
2 s-2) (5’)
ただし、ω2=2πf2
ω3
2 s-2/(1+2h3ω3 s-1+ω3
2s-2) (6’)
ただし、ω3=2πf3
ω4
2 s-2/(1+2h4ω4 s-1+ω4
2s-2) (7’)
ただし、ω4=2πf4
ω5
2 s-2/(1+2h5ω5 s-1+ω5
2s-2) (8’)
ただし、ω5=2πf5である。
【0051】
このとき、フィルタの特性を規定するパラメータは、(f0,f1,f2,f3,f4,f5,h2a,h2b,h3,h4,h5,g)の12個となる。これらの値を本来の振幅特性との一致の度合いや、最終的な計算精度を参考にして決定したところ、近似の度合いの良いパラメータの組み合わせとして、(f0=0.45[Hz],f1=7.0[Hz],f2=0.5[Hz],f3=12.0[Hz],f4=20.0[Hz],f5=30.0[Hz] ,h2a=1.0,h2b=0.75,h3=0.9,h4=0.6,h5=0.6,g=1.262)が得られている。
【0052】
フィルタのラプラス変換表示が与えられれば、z変換を用いてデジタルフィルタを構成することができる。特許文献2では、ΔTをサンプリング間隔、zを遅延演算子として、
s-1=(ΔT/2)(1+ z-1)/(1- z-1) (30)
及び
s-2=(ΔT/121/2)2(1+ 10z-1+ z-2)/(1- z-1) 2 (31)
というz変換を用いて、以下の式(11’)乃至式(14’)で表されるデジタルフィルタを得た。
【0053】
・フィルタ1
α0=8/(ΔT)2+(4ωa1+2ωa2 )/(ΔT)+ ωa1ωa2 (11’)
α1=2ωa1ωa2-16/(ΔT)2
α2=8/(ΔT)2-(4ωa1+2ωa2 )/(ΔT)+ ωa1ωa2
β0=4/(ΔT)2+2ωa2 /(ΔT)
β1=-8/(ΔT)2
β2=4/(ΔT)2-2ωa2 /(ΔT)
ωa1=2πfa1, ωa2=2πfa2
【0054】
・フィルタ2
α0=16/(ΔT)2+17ωa3/(ΔT)+ ωa3
2
(12’)
α1=2ωa3
2-32/(ΔT)2
α2=16/(ΔT)2-17ωa3/(ΔT)+ ωa3
2
β0=4/(ΔT)2+8.5ωa3 /(ΔT) + ωa3
2
β1=2ωa3
2-8/(ΔT)2
β2=4/(ΔT)2-8.5ωa3 /(ΔT) + ωa3
2
ωa3=2πfa3
【0055】
・フィルタ3
α0=12/(ΔT)2+12hb2ωb /(ΔT)+ωb
2 (13’)
α1=10ωb
2-24/(ΔT)2
α2=12/(ΔT)2-12hb2ωb /(ΔT)+ωb
2
β0=12/(ΔT)2+12hb1ωb /(ΔT)+ωb
2
β1=10ωb
2-24/(ΔT)2
β2=12/(ΔT)2-12hb1ωb /(ΔT)+ωb
2
ωb=2πfb
【0056】
・フィルタ4
α0=12/(ΔT)2+12hcωc /(ΔT)+ωc
2 (14’)
α1=10ωc
2-24/(ΔT)2
α2=12/(ΔT)2-12hcωc /(ΔT)+ωc
2
β0=ωc
2
β1=10ωc
2
β2=ωc
2
ωc=2πfc
【0057】
入力時系列x(k)にこれらのフィルタを作用させた出力時系列y(k)を得るには、以下の式(15’)の演算を行えばよい。
y(k)=[-α1y(k-1) -α2y(k-2)+β0x(k) +β1x(k-1) +β2x(k-2)]/α0 (15’)
ただし、
x(k)はフィルタの入力時系列、
y(k)はフィルタの出力時系列、
kは時間ステップ数、
ΔTはサンプリング間隔、
πは円周率、
とする。
【0058】
また、ゲイン調整はゲインを gdとすれば、入力時系列をx(k)、出力時系列をy(k)として、以下の式(16’)を演算すればよい。
y(k)=gdx(k) (16’)
【0059】
以上をまとめると、特許文献2に記載された近似フィルタ演算としては、加速度記録に対して、下記の7つのデジタルフィルタ演算を直列に行えばよい。
演算1’: fa1=f0, fa2=f1として、式(11’),(15’)の演算
演算2’: fa3=f1として、式(12’),(15’)の演算
演算3’: hb1 = h2a, hb2 = h2b, fb=f2 として、式(13’),(15’)の演算
演算4’: hc = h3, fc=f3 として、式(14’),(15’)の演算
演算5’: hc = h4, fc=f4 として、式(14’),(15’)の演算
演算6’: hc = h5, fc=f5 として、式(14’),(15’)の演算
演算7’: gd= g として、式(16’)の演算
【0060】
特許文献2では、実用的な地震防災に役立てるために、地動加速度時系列のサンプリング間隔ΔTは、実際に運用されている震度計の標準である0.01秒を想定してフィルタを設計し、精度検証を行った。一方、研究等を目的として標準の0.01秒より長いΔTを利用する場合もある。この際に、ΔTと,f2,f3,f4,f5 の組み合わせによっては、演算3’乃至演算6’が不安定となり出力時系列が発散する場合がある。
【0061】
出力時系列の発散は、以下の式(31)に示す変換(Boxer-Thaler integrator)を用いたために生じている。
s-2=(ΔT/121/2)2(1+ 10z-1+ z-2)/(1- z-1) 2 (31)
式(31)に示す変換は、計算精度は高いが、ΔTの値によっては得られるデジタルフィルタが不安定性を生じる可能性がある。
【0062】
一方、以下の式(32)に示す変換(Tustin integrator)を用いると、ΔTの値によらず、得られるデジタルフィルタは安定する。
s-2= s-1 s-1= (ΔT/2)(1+ z-1)/(1- z-1) (ΔT/2)(1+ z-1)/(1- z-1)
=(ΔT/41/2)2(1+ 2z-1+ z-2)/(1- z-1) 2 (32)
しかし、式(32)に示す変換は、計算精度が低い。
【0063】
さらに、以下の式(33)に示す変換(Madwed integrator)も可能である。
s-2=(ΔT/61/2)2(1+ 4z-1+ z-2)/(1- z-1) 2 (33)
式(33)に示す変換は、安定性と計算精度に関して、前記二つの変換の中間の性質を持つ。
【0064】
これらの式(31)乃至式(33)の3つの変換は、
s-2=γ(ΔT)2(1+(1/γ-2)z-1+ z-2)/(1- z-1) 2 (34)
という形にまとめて表現できる。したがって、本実施形態では、式(34)の変換を用いて、デジタルフィルタを構成し最適なγを利用できるようにする。なお、本実施形態ではγの値は、1/12,1/6,1/4に限定しない。ここで、γはフィルタの二次積分の近似方法を
規定するパラメータとなる。
【0065】
式(34)の変換を用いれば、特許文献2に示した式(13’)と式(14’)で表されるデジタルフィルタの替わりに、本実施形態では下記の式(13)と式(14)で表されるデジタルフィルタを得る。なお、変換に際しては、置き換えるs-2毎に異なるγの値(式(13)と式(14)では、γb1、γb2、γc1、γc2)をとることが可能なものとする。
【0066】
・第3フィルタ
α0=1+ hb2ωbΔT +γb2ωb
2 (ΔT)2 (13)
α1=-2+ωb
2(ΔT)2(1-2γb2)
α2=1- hb2ωbΔT +γb2ωb
2 (ΔT)2
β0=1+ hb1ωbΔT +γb1ωb
2 (ΔT)2
β1=-2+ωb
2(ΔT)2(1-2γb1)
β2=1- hb1ωbΔT +γb1ωb
2 (ΔT)2
ωb=2πfb
【0067】
・第4フィルタ
α0=1+ hcωcΔT +γc2ωc
2 (ΔT)2 (14)
α1=-2+ωc
2(ΔT)2(1-2γc2)
α2=1- hcωcΔT +γc2ωc
2 (ΔT)2
β0=γc1ωc
2(ΔT)2
β1=ωc
2(ΔT)2(1-2γc1)
β2=γc1ωc
2(ΔT)2
ωc=2πfc
【0068】
図9は、本実施形態の時間領域フィルタ装置のブロック図である。
【0069】
以上をまとめると、本実施形態の地動加速度の時系列をフィルタ処理する時間領域フィルタ装置は、以下の式(11)において、フィルタ処理する第1フィルタ721と、以下の式(12)において、フィルタ処理する第2フィルタ722と、以下の式(13)において、フィルタ処理する第3フィルタ723と、以下の式(14)において、フィルタ処理する第4フィルタ724と、以下の式(15)において、入力時系列x(k)にこれらのフィルタを作用させて出力時系列y(k)を取得する出力時系列取得部725と、以下の式(16)において、ゲインを調整するゲイン調整部726と、演算11乃至演算17を実行する演算部727と、を備える。そして、 以下の式(11)乃至(16)において、演算11乃至演算17の少なくとも一つを実行する。
・第1フィルタ
α0=8/(ΔT)2+(4ωa1+2ωa2 )/(ΔT)+ωa1ωa2 (11)
α1=2ωa1ωa2-16/(ΔT)2
α2=8/(ΔT)2-(4ωa1+2ωa2 )/(ΔT)+ωa1ωa2
β0=4/(ΔT)2+2ωa2 /(ΔT)
β1=-8/(ΔT)2
β2=4/(ΔT)2-2ωa2 /(ΔT)
ωa1=2πfa1, ωa2=2πfa2
・第2フィルタ
α0=16/(ΔT)2+17ωa3/(ΔT)+ ωa3
2
(12)
α1=2ωa3
2-32/(ΔT)2
α2=16/(ΔT)2-17ωa3/(ΔT)+ ωa3
2
β0=4/(ΔT)2+8.5ωa3 /(ΔT) + ωa3
2
β1=2ωa3
2-8/(ΔT)2
β2=4/(ΔT)2-8.5ωa3 /(ΔT) + ωa3
2
ωa3=2πfa3
・第3フィルタ
α0=1+ hb2ωbΔT +γb2ωb
2 (ΔT)2 (13)
α1=-2+ωb
2(ΔT)2(1-2γb2)
α2=1- hb2ωbΔT +γb2ωb
2 (ΔT)2
β0=1+ hb1ωbΔT +γb1ωb
2 (ΔT)2
β1=-2+ωb
2(ΔT)2(1-2γb1)
β2=1- hb1ωbΔT +γb1ωb
2 (ΔT)2
ωb=2πfb
・第4フィルタ
α0=1+ hcωcΔT +γc2ωc
2(ΔT)2 (14)
α1=-2+ωc
2 (ΔT)2(1-2γc2)
α2=1- hcωcΔT +γc2ωc
2(ΔT)2
β0=γc1ωc
2(ΔT)2
β1=ωc
2 (ΔT)2(1-2γc1)
β2=γc1ωc
2(ΔT)2
ωc=2πfc
・出力時系列取得部
y(k)=[-α1y(k-1) -α2y(k-2)+β0x(k) +β1x(k-1) +β2x(k-2)]/α0 (15)
ただし、
x(k)はフィルタの入力時系列、
y(k)はフィルタの出力時系列、
kは時間ステップ数、
ΔTはサンプリング間隔、
πは円周率、
とする。
・ゲイン調整部
y(k)=gd x(k) (16)
ただし、gd はゲインである。
・演算部
演算11:fa1 = f0, fa2 = f1として、式(11),(15)の演算
演算12:fa3 = f1として、式(12),(15)の演算
演算13:hb1= h2a, hb2= h2b, fb= f2, γb1= γ2a, γb2= γ2b として、式(13),(15)の演算
演算14:hc= h3, fc= f3, γc1= γ3a, γc2= γ3b として、式(14),(15)の演算
演算15:hc = h4, fc = f4, γc1 = γ4a, γc2 = γ4b として、式(14),(15)の演算
演算16:hc = h5, fc = f5, γc1 = γ5a, γc2 = γ5b として、式(14),(15)の演算
演算17:gd = g として、式(16)の演算
ただし、
f0はフィルタのカットオフ周波数を規定するパラメータ、
f1は0.5次で減衰する周波数領域の位置を規定するパラメータ、
f2、h2a、h2b、は補正フィルタの特性を規定するパラメータ、
f3、h3はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
f4、h4はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
f5、h5はフィルタのカットオフ周波数およびダンピングを規定するパラメータ、
gはゲイン調整パラメータ、
γ2a、γ2b、γ3a、γ3b, γ4a、γ4b、γ5a、γ5b はフィルタの二次積分の近似方法を規定するパラメータ、
である。
【0070】
次に、安定なデジタルフィルタを得るために必要なγの条件を導出する。デジタルフィルタが安定であるためにはすべての極が単位円の内側にある必要があり、これから、以下の式(21)乃至式(24)が安定なデジタルフィルタであるための条件となる。
(1/4-1/(2πf2ΔT)2)≦γ2b (21)
(1/4-1/(2πf3ΔT)2)≦γ3b (22)
(1/4-1/(2πf4ΔT)2)≦γ4b (23)
(1/4-1/(2πf5ΔT)2)≦γ5b (24)
【0071】
したがって、fb、fcが与えられた場合に、式(21)乃至式(24)を満たすような、γ2b、γ3b、γ4b、γ5b をΔTに応じて選べばよい。このとき、γ2a、γ3a、γ4a、γ5aの選び方には安定性に関する条件に影響されない。
【0072】
また、時間領域フィルタ装置は、フィルタの二次積分の近似方法を規定するパラメータにおいて、下記の式(25)乃至(28)のうち少なくとも一つを満たしてもよい。
γ2a = γ2b (25)
γ3a = γ3b (26)
γ4a = γ4b (27)
γ5a = γ5b (28)
【0073】
実際のγの選び方に関しては、f2=0.5[Hz],f3=12.0[Hz],f4=20.0[Hz],f5=30.0[Hz]の場合は、サンプリング周波数 fs = 1/ΔTとして、として、下記のようにすればよい。
【0074】
(1)5種類のγを用いる場合、
fs ≧ 80[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/12, γ5b = 1/12
80[Hz] > fs ≧ 70[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/12, γ5b = 1/8
70[Hz] > fs ≧ 60[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/12, γ5b = 1/6
60[Hz] > fs ≧ 50[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/10, γ5b = 1/4
50[Hz] > fs ≧ 40[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/6, γ5b = 1/4
40[Hz] > fs ≧ 30[Hz]のとき
γ2b = 1/12, γ3b= 1/10, γ4b = 1/4, γ5b = 1/4
30[Hz] > fs ≧ 5[Hz]のとき
γ2b = 1/12, γ3b= 1/4, γ4b = 1/4, γ5b = 1/4
5[Hz] > fs ≧ 1[Hz]のとき
γ2b = 1/6, γ3b = 1/4, γ4b= 1/4, γ5b = 1/4
【0075】
(2)4種類のγを用いる場合、
fs ≧ 80[Hz]のとき
γ2b = 1/12, γ3b = 1/12, γ4b = 1/12, γ5b= 1/12
80[Hz] > fs ≧ 70[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/12, γ5b = 1/8
70[Hz] > fs ≧ 60[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/12, γ5b = 1/6
60[Hz] > fs ≧ 50[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/8, γ5b = 1/4
50[Hz] > fs ≧ 40[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/6, γ5b = 1/4
40[Hz] > fs ≧ 30[Hz]のとき
γ2b = 1/12, γ3b= 1/8, γ4b = 1/4, γ5b = 1/4
30[Hz] > fs ≧ 5[Hz]のとき
γ2b = 1/12, γ3b= 1/4, γ4b = 1/4, γ5b = 1/4
5[Hz] > fs ≧ 1[Hz]のとき
γ2b = 1/6, γ3b = 1/4, γ4b= 1/4, γ5b = 1/4
【0076】
(3)3種類のγを用いる場合、
fs ≧ 80[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/12, γ5b = 1/12
80[Hz] > fs ≧ 60[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/12, γ5b = 1/6
60[Hz] > fs ≧ 40[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/6, γ5b = 1/4
40[Hz] > fs ≧ 30[Hz]のとき
γ2b = 1/12, γ3b= 1/6, γ4b = 1/4, γ5b = 1/4
30[Hz] > fs ≧ 5[Hz]のとき
γ2b = 1/12, γ3b= 1/4, γ4b = 1/4, γ5b = 1/4
5[Hz] > fs ≧ 1[Hz]のとき
γ2b = 1/6, γ3b = 1/4, γ4b= 1/4, γ5b = 1/4
【0077】
(4)2種類のγを用いる場合、
fs ≧ 80[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/12, γ5b = 1/12
80[Hz] > fs ≧ 60[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/12, γ5b = 1/4
60[Hz] > fs ≧ 40[Hz]のとき
γ2b = 1/12, γ3b= 1/12, γ4b = 1/4, γ5b = 1/4
40[Hz] > fs ≧ 5[Hz]のとき
γ2b = 1/12, γ3b= 1/4, γ4b = 1/4, γ5b = 1/4
5[Hz] > fs ≧ 1[Hz]のとき
γ2b = 1/4, γ3b = 1/4, γ4b= 1/4, γ5b = 1/4
【0078】
(5)1種類のγを用いる場合で演算精度を高める場合、
γ2b = 1/12, γ3b= 1/12, γ4b = 1/12, γ5b = 1/12
このγは、f2=0.5[Hz],f3=12.0[Hz],f4=20.0[Hz],f5=30.0[Hz]、以外の場合にも適用可能である。
【0079】
(6)1種類のγを用いる場合で安定性の保証を行う場合、
γ2b = 1/4, γ3b = 1/4, γ4b= 1/4, γ5b = 1/4
このγは、f2=0.5[Hz],f3=12.0[Hz],f4=20.0[Hz],f5=30.0[Hz]、以外の場合にも適用可能である。
【0080】
式(13)及び式(14)を用いれば常に安定なデジタルフィルタをサンプリング間隔ΔTに応じて得ることができる。したがって、サンプリング間隔ΔTがステップ毎に異なる値とする地動加速度時系列に対しても時間領域フィルタ演算を行うことが可能である。
【0081】
地動加速度時系列のサンプリング周波数 fs が5[Hz]等の低い値の場合は、f5=30[Hz]などの、高いカットオフ周波数を持つフィルタ演算は無駄な処理となる。このような場合は、該当するフィルタ処理を迂回すればよい。このような場合でも本実施形態のフィルタ演算は動作可能である。
【0082】
図10は、地震波形と地震波形を補間し従来の方法を用いた波形と本実施形態の計測震度概算装置を用いた波形を示す図である。
図10(a)は50Hzの地震波形、
図10(b)は50Hzの地震波形を補間して100Hzとして従来の方法で計測震度を概算した波形、
図10(c)は本実施形態の計測震度概算装置を用いて50Hzのままで計算した波形である。
【0083】
図10(b)と
図10(c)を比較すると、ほぼ同じような結果を示している。しかしながら、
図10(b)に示した従来の装置は、補間を行うため、
図10(c)に示した本実施形態の計測震度概算装置と比較して、即時性が劣ることになる。
【0084】
図11は、本実施形態の計測震度概算装置をシステムとして構築した場合を示す。図中、101は計測震度概算システム、102は地震計、103は通信ネットワーク、104は地震データ処理装置である。このシステムでは、複数の地震計102で観測された地震動を、通信ネットワーク103を通じ他の機器やデータセンターに転送し、離れた場所で計測震度概算部7を有する地震データ処理装置104による概算の震度値算出を行うものである。ここで、演算処理を行う地震データ処理装置は複数存在してもよい。このようなシステムにおいてもいち早く震度の概算値を得られるという本方法の利点が活かされる。データセンターでは多くの観測点のデータを同時に処理しなければならない場合が多いため、FFTを用いない本方法のアルゴリズムは計算負荷を軽減する上でも有用である。
【符号の説明】
【0085】
1…計測震度概算装置、2…処理部、3…制御・演算部、4…計測震度算出部、5…SI値算出部、6…応答値算出部、7…計測震度概算部、71…地動加速度時系列取得部、72…時間領域フィルタ装置、721…第1フィルタ、722…第2フィルタ、723…第3フィルタ、724…第4フィルタ、725…第1フィルタ出力時系列取得部、726…ゲイン調整部、727…演算部、73…算出部、74…振幅時系列変換部、75…震度換算振幅時系列換算部、76…離散化部、77…継続時間カウント部、8…電源装置、9…計測装置、101…計測震度概算システム、102…地震計、103…通信ネットワーク、104…地震データ処理装置