JP4229337B2 - 計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法 - Google Patents
計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法 Download PDFInfo
- Publication number
- JP4229337B2 JP4229337B2 JP2007249485A JP2007249485A JP4229337B2 JP 4229337 B2 JP4229337 B2 JP 4229337B2 JP 2007249485 A JP2007249485 A JP 2007249485A JP 2007249485 A JP2007249485 A JP 2007249485A JP 4229337 B2 JP4229337 B2 JP 4229337B2
- Authority
- JP
- Japan
- Prior art keywords
- seismic intensity
- time series
- filter
- value
- acceleration
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Landscapes
- Geophysics And Detection Of Objects (AREA)
Description
本発明は、地震の加速度を計測し、その計測値から震度を概算する計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法に関する。
従来から、地震が発生した際の警報発報や鉄道等の制御の指標として、気象庁が官報第１８３１号気象庁告示第４号気象業務法施行規則に関する告示で示した計測震度が用いられているのが現状である。計測震度を算出するためには、地震記録を一定時間蓄積してから、ＦＦＴを用いてフィルタ処理をすることが必要なため、刻々と測定される地震記録に対して値を算出することができず、即時性が要求される用途に適していなかった。
また、最大加速度、最大速度、ＳＩ値などの地震動指標を測定し、その地震動指標に基づいて簡単な演算処理をすることで震度相当の値を推定する方法がある（特許文献１参照）。
ところが、この方法は、即時性があるが、換算の誤差が大きく、また換算式設計に際し、過去の大量の地震データを使用した解析を必要とする。
本発明は、前記課題を解決するために、即時性があり、換算の誤差の小さい計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法を提供することを目的とする。
前記課題を解決するために、本発明の計測震度概算装置は、地震の加速度を計測し、計測した加速度から計測震度を概算する計測震度概算装置において、地動加速度の時系列を得る地動加速度時系列取得手段と、前記地動加速度の時系列をフィルタ処理する時間領域フィルタ手段と、前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列から計測震度の概算値を算出する算出手段とを備えたことを特徴とする。
また、前記算出手段は、前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列を振幅時系列に変換する振幅時系列変換手段と、前記振幅時系列を震度換算振幅時系列に換算する震度換算振幅時系列換算手段と、震度換算振幅時系列換算手段により求めた震度換算振幅時系列を離散化する離散化手段と、前記離散化手段により求めた各離散化震度換算振幅値に継続時間カウントを実行する継続時間カウント手段と、を有することを特徴とする。
前記時間領域フィルタ手段は、下記の式（１１）、（１２）及び（１４―１）乃至（１６）において、演算１乃至演算6を実行することを特徴とする。
y(k)=-α1/α0×y(k-1)+β0/α0×x(k) +β1/α0×x(k-1) ・・・（１２）
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) ・・・（１５）
+β1/α0×x(k-1)+β2/α0×x(k-2) ・・・（１５）
fb、hbは周波波数特性を規定するパラメータ
y(k)=gc×x(k) ・・・（１６）
ただし、gcはゲイン
ただし、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はゲイン調整パラメータ
演算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) ・・・（１２）
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)
・・・（１５）
・・・（１５）
fb、hbは周波波数特性を規定するパラメータ
y(k)=gc×x(k) ・・・（１６）
ただし、gcはゲイン
ただし、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はゲイン調整パラメータ
演算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]に変換する（ＳＴ４）。
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]=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]のどの要素がいつの時間ステップの値であるかを記憶しておく。これで初期化が完了する。
ステップ２０１において、初期化が行われていないと判断された場合、ステップ２０２で、初期化を実行する（ＳＴ２０２）。初期化について説明する。まず、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=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次で減衰する周波数領域の位置を規定するパラメータである。
×(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つのフィルタを順次作用させることに等しい。
(ω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/πと置き換えてもよい。
ω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） ・・・（９）
を得た。図１０には、このパラメータを用いたときの総合的な近似フィルタの振幅特性と本来の振幅特性が示されている。
（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) ・・・（１０）
双一次変換、s-1=(ΔT/2)×(1+ z-1)/( 1- z-1)を適用すれば
F(z)=(β0+z-1×β1)/(α0+z-1×α1) ・・・（１０）
入力時系列x(k)にフィルタF(z)を作用させた出力時系列y(k)を得るには、
y(k)=-α1/α0×y(k-1)+β0/α0×x(k) +β1/α0×x(k-1) ・・・（１２）
なる演算を行えばよい。ここで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) ・・・（１３）
G(z)=(β0+z-1×β1+z-2×β2)/(α0+z-1×α1+z-2×α2) ・・・（１３）
また、ラプラス変換形式で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) ・・・（１３）
G(z)=(β0+z-1×β1+z-2×β2)/(α0+z-1×α1+z-2×α2) ・・・（１３）
入力時系列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) ・・・（１５）
なる演算を行えばよい。
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) ・・・（１６）
なる演算を行えばよい。
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)の演算
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を用いない本方法のアルゴリズムは計算負荷を軽減する上でも有用である。
１…計測震度概算装置、２…処理部、３…制御・演算部、４…計測震度算出部、５…ＳＩ値算出部、６…応答値算出部、７…計測震度概算部、７１…地動加速度時系列取得手段、７２…時間領域フィルタ手段、７３…算出手段、７４…振幅時系列変換手段、７５…震度換算振幅時系列換算手段、７６…離散化手段、７７…継続時間カウント手段、８…電源装置、９…計測装置、１０１…計測震度概算システム、１０２…地震計、１０３…通信ネットワーク、１０４…地震データ処理装置
Claims (5)
- 地震の加速度を計測し、計測した加速度から計測震度を概算する計測震度概算装置において、地動加速度の時系列を得る地動加速度時系列取得手段と、前記地動加速度の時系列をフィルタ処理する時間領域フィルタ手段と、前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列から計測震度の概算値を算出する算出手段とを備え、
前記算出手段は、前記時間領域フィルタによりフィルタ処理された前記地動加速度の時系列を振幅時系列に変換する振幅時系列変換手段と、前記振幅時系列を震度換算振幅時系列に換算する震度換算振幅時系列換算手段と、震度換算振幅時系列換算手段により求めた震度換算振幅時系列を離散化する離散化手段と、前記離散化手段により求めた各離散化震度換算振幅値に継続時間カウントを実行する継続時間カウント手段と、を有することを特徴とする計測震度概算装置。 - 前記時間領域フィルタ手段は、下記の式（１１）、（１２）及び（１４−１）乃至（１６）において、演算１乃至演算６を実行することを特徴とする請求項１に記載された計測震度概算装置。
y(k)=-α1/α0×y(k-1)+β0/α0×x(k)+β1/α0×x(k-1) ・・・（１２）
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) ・・・（１５）
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) ・・・（１２）
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) ・・・（１５）
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はゲイン調整パラメータ - 地震動を計測する複数の地震計と、通信ネットワークと、前記複数の地震計で観測された地震動を、前記通信ネットワークを通じ、離れた場所で前記計測震度概算装置による概算の震度値算出を実行することを特徴とする請求項１乃至請求項３のいずれか１項に記載された計測震度概算装置を用いた計測震度概算システム。
- 地震の加速度を計測し、計測した加速度から計測震度を概算する計測震度概算方法において、地動加速度時系列を取得するステップと、前記地動加速度時系列を時間領域でフィルタ処理し、フィルタ処理地動加速度時系列を取得するステップと、前記フィルタ処理地動加速度の時系列から計測震度の概算値を算出するステップとからなり、
前記フィルタ処理地動加速度の時系列から計測震度の概算値を算出するステップは、前記フィルタ処理地動加速度時系列から振幅時系列を取得するステップと、前記振幅時系列を震度換算振幅時系列に換算するステップと、前記震度換算振幅時系列を離散化し、各離散化震度換算振幅値を取得するステップと、前記各離散化震度換算振幅値について継続時間をカウントするステップと、前記各離散化震度換算振幅値と継続時間のカウント数から計測震度概算値を算出するステップとからなることを特徴とする計測震度概算方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2007249485A JP4229337B2 (ja) | 2006-09-28 | 2007-09-26 | 計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法 |
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2006265192 | 2006-09-28 | ||
| JP2007249485A JP4229337B2 (ja) | 2006-09-28 | 2007-09-26 | 計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2008107334A JP2008107334A (ja) | 2008-05-08 |
| JP4229337B2 true JP4229337B2 (ja) | 2009-02-25 |
Family
ID=39440774
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2007249485A Active JP4229337B2 (ja) | 2006-09-28 | 2007-09-26 | 計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP4229337B2 (ja) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2014178226A (ja) * | 2013-03-15 | 2014-09-25 | National Research Institute For Earth Science & Disaster Provention | 計測震度概算システム及び計測震度概算方法 |
| JP2014228373A (ja) * | 2013-05-22 | 2014-12-08 | 独立行政法人防災科学技術研究所 | 計測震度概算装置及びそれを用いた計測震度概算システム |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4551955B2 (ja) * | 2008-11-25 | 2010-09-29 | シャープ株式会社 | 温水器 |
| JP5313706B2 (ja) * | 2009-01-27 | 2013-10-09 | 大成建設株式会社 | 地震動波形推定方法 |
| JP5507903B2 (ja) * | 2009-06-19 | 2014-05-28 | 白山工業株式会社 | 震度推定方法及び装置 |
| JP5375435B2 (ja) * | 2009-08-25 | 2013-12-25 | 株式会社ホームサイスモメータ | 震度測定装置 |
| JP5819222B2 (ja) * | 2012-03-01 | 2015-11-18 | 株式会社シグネット | リアルタイム震度計測装置とその方法 |
| JP6457276B2 (ja) * | 2015-01-21 | 2019-01-23 | 国立研究開発法人防災科学技術研究所 | 地震動補正装置、それを用いた地震動補正システム、及び地震動補正方法 |
| JP6872794B2 (ja) * | 2017-09-26 | 2021-05-19 | 国立研究開発法人防災科学技術研究所 | 地震動計測装置、それを用いた地震動計測システム、及び地震計の傾斜補正方法 |
| CN117607967B (zh) * | 2024-01-19 | 2024-03-26 | 中国建筑西南设计研究院有限公司 | 一种基于遗传算法的地震动基线校正方法及电子设备 |
-
2007
- 2007-09-26 JP JP2007249485A patent/JP4229337B2/ja active Active
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2014178226A (ja) * | 2013-03-15 | 2014-09-25 | National Research Institute For Earth Science & Disaster Provention | 計測震度概算システム及び計測震度概算方法 |
| JP2014228373A (ja) * | 2013-05-22 | 2014-12-08 | 独立行政法人防災科学技術研究所 | 計測震度概算装置及びそれを用いた計測震度概算システム |
Also Published As
| Publication number | Publication date |
|---|---|
| JP2008107334A (ja) | 2008-05-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP4229337B2 (ja) | 計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法 | |
| Noël et al. | Subspace-based identification of a nonlinear spacecraft in the time and frequency domains | |
| Wan et al. | Structural response reconstruction based on the modal superposition method in the presence of closely spaced modes | |
| US6968296B2 (en) | Cable detector with decimating filter and filtering method | |
| US9161126B2 (en) | Systems and methods for protecting a speaker | |
| JP3771195B2 (ja) | 重量測定用ノイズ除去装置および重量測定用ノイズ除去方法 | |
| Calabrese et al. | Adaptive constrained unscented Kalman filtering for real‐time nonlinear structural system identification | |
| US20150330950A1 (en) | Structural fatigue crack monitoring system and method | |
| CN103547899B (zh) | 振动监视系统 | |
| US20140254804A1 (en) | Systems and methods for protecting a speaker | |
| CN110794170B (zh) | 一种加速度计两自由度动态模型参数辨识的方法 | |
| JP5375435B2 (ja) | 震度測定装置 | |
| Carbajo et al. | ASDAH: An automated structural change detection algorithm based on the Hilbert–Huang transform | |
| Spanos et al. | Numerical treatment of seismic accelerograms and of inelastic seismic structural responses using harmonic wavelets | |
| CN111024214B (zh) | 一种实时获取声共振混合机运行过程中固有频率的方法 | |
| Cooper et al. | On-line physical parameter estimation with adaptive forgetting factors | |
| Spanos et al. | Deterministic and stochastic analyses of a nonlinear system with a Biot visco‐elastic element | |
| JP2867477B2 (ja) | オンライン設備の寿命予測方法 | |
| JP5946067B2 (ja) | 計測震度概算装置及びそれを用いた計測震度概算システム | |
| JP2023139362A (ja) | 時間領域フィルタ装置、計測震度概算装置、計測震度概算システム及び時間領域フィルタ方法 | |
| Li et al. | Parametric time‐domain identification of multiple‐input systems using decoupled output signals | |
| Tuhta | OMA of model chimney using Bench-Scale earthquake simulator | |
| Prawin et al. | Extraction of opening and closing states of cracked structure using adaptive volterra filter model | |
| Ghrib et al. | An adaptive filtering-based solution for the Bayesian modal identification formulation | |
| Bicchi et al. | Optimal design of dynamic multi-axis force/torque sensor |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20080618 |
|
| A871 | Explanation of circumstances concerning accelerated examination |
Free format text: JAPANESE INTERMEDIATE CODE: A871 Effective date: 20080618 |
|
| A975 | Report on accelerated examination |
Free format text: JAPANESE INTERMEDIATE CODE: A971005 Effective date: 20080716 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20080910 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20081024 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20081126 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20081127 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 4229337 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111212 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111212 Year of fee payment: 3 |
|
| S111 | Request for change of ownership or part of ownership |
Free format text: JAPANESE INTERMEDIATE CODE: R313117 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111212 Year of fee payment: 3 |
|
| R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121212 Year of fee payment: 4 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20131212 Year of fee payment: 5 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| S533 | Written request for registration of change of name |
Free format text: JAPANESE INTERMEDIATE CODE: R313533 |
|
| R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |