JP5819222B2 - リアルタイム震度計測装置とその方法 - Google Patents

Description

本発明は、地震発生後のリアルタイム震度を測定可能な震度計測装置とその方法に係り、特に、リアルタイムでも測定精度を高めることが可能なリアルタイム震度計測装置とその方法に関する。

地震が頻繁に発生する我が国では、地震が発生すると気象庁等が地震情報を発信し、ラジオ番組の放送中やテレビ番組の放映中であっても地震速報として、音声や字幕で各地での震度、震源及び地震規模（マグニチュード）が流されている。
しかしながら、その地震情報は震度３以上の場合に最も早く発信される震度速報であっても通常２分程度要している。これは、地震震度が、地震波の加速度センサによって測定される地震加速度の時間軸上の波形データを一定時間蓄積した上で、ＦＦＴ（Fast Fourier Transform：高速フーリエ変換）演算を気象庁が定めるＦＦＴフィルターで処理して行った後、それを逆ＦＦＴ変換して得られた加速度波形データを震度計算に用いているためである。
すなわち、リアルタイムに震度が計測されるのではなく、一旦地震の加速度波形データを一定時間蓄積した上で解析されているため、時間遅れが生じるのである。このような震度計測に関しては、官報第１８３１号で告示された気象庁震度階級表（気象庁告示第４号）に記載されている。
そこで、従来このような時間遅れを解消してリアルタイムに震度を計測するための研究がおこなわれており、その成果が発明として出願され開示されている。

例えば、特許文献１には、「計測震度概算装置、それを用いた計測震度概算システム及び計測震度概算方法」という発明が開示されている。
特許文献１に開示された発明は、リアルタイム性と換算誤差の小さな計測震度概算装置等を提供するものであり、ＦＦＴ演算を行うことなく地震の加速度を計測して地動加速度の時系列を得る地動加速度時系列取得手段と、地動加速度の時系列をフィルタ処理する時間領域フィルタ手段と、時間領域フィルタによりフィルタ処理された地動加速度の時系列から計測震度の概算値を算出する算出手段を備えたことを特徴としている。
このような特許文献１に開示される発明では、地震加速度波形データを蓄積することなく、時間領域フィルタによってフィルタ処理するために、リアルタイム性に優れると共に、フィルタ係数の演算を効率的に行うなどしてある程度の精度を保持することができるので、計測震度を概算することができるのである。

また、特許文献２には、「震度測定装置」という発明が開示されている。
この発明では、加速度波形データを入力する加速度入力部と、加速度波形データにフィルタ処理を施し、補正加速度波形データを出力するフィルタと、補正加速度波形データから震度を算出する震度算出部とを備えて、フィルタは加算、減算、２のべき乗数の乗算及び２のべき乗の除算を組み合わせた演算を実行することで、加速度波形データを補正加速度波形データに変換することを特徴とするものである。
このような特許文献２に開示される発明では、回路構成を複雑化させることなく、低消費電力で気象庁定義のフィルタを設計可能とすることができるのである。

特開２００８−１０７３３４号公報 特開２０１１−４７６５７号公報

しかしながら、特許文献１に開示された発明においては、時間領域フィルタを用いることでリアルタイム性を改善することが可能となったものの、その計測震度は発明の名称のとおり、概算に留まっており、その精度に改善の余地があった。

また、特許文献２に開示された発明においても、消費電力の低減は可能かもしれないものの、特許文献２の図３に示されるように、気象庁定義のフィルタ特性との近似にはズレが生じており、その精度には改善の余地があるという課題があった。

本発明はかかる従来の事情に対処してなされたものであり、地震発生直後から、リアルタイム性と高精度の両方を備えた震度計測を可能として、震度速報の発信までの時間短縮を行いつつ、震源に近い被災地であっても地震対応や避難が迅速かつ適切に実行可能とするリアルタイム震度計測装置とその方法を提供することを目的とする。

上記目的を達成するため、本願発明のリアルタイム震度計測装置は、地震による加速度を検出する加速度検出手段と、この加速度検出手段から出力されるアナログ信号を第１のデジタル信号に変換するＡ／Ｄ変換器と、このＡ／Ｄ変換器から出力される前記第１のデジタル信号をフィルター処理して第２のデジタル信号を出力するＩＩＲ（Infinite Impulse Response:無限インパルス応答）フィルターと、このＩＩＲフィルターから出力される第２のデジタル信号に基づいてリアルタイム震度を計算する震度演算部と、このリアルタイム震度を出力する出力部と、を有し、前記ＩＩＲフィルターは、ゲインピークを含む周波数帯域において、（ｎ＋１）次ハイパスフィルター（ＨＰＦ）（ｎは正整数、以下同様）とｎ次ハイパスフィルター（ＨＰＦ）を組み合わせて構成される第１のＩＩＲフィルターと、を備えたことを特徴とするものである。
このように構成されたリアルタイム震度計測装置においては、ＩＩＲフィルターがリアルタイムにフィルター処理するという作用を有することはもちろんのこと、ＩＩＲフィルターが（ｎ＋１）次ＨＰＦとｎ次ＨＰＦを組み合わせて構成されるため、ｎ次ＨＰＦのみで構成されるよりも（ｎ＋１）次ほど自由度が増え、より高い精度のフィルタリング特性を発揮させる作用を有する。
請求項１に記載の発明は、このリアルタイム震度計測装置において、前記ゲインピークを含む周波数（以下、周波数をｆで表わす）帯域は０．３Ｈｚ＜ｆ≦１Ｈｚの範囲を含む帯域であって、前記ｎは１であり第１のＩＩＲフィルターは遮断周波数と減衰定数の２つの自由度を有する２次ＨＰＦと遮断周波数の１つの自由度を有する１次ＨＰＦを組み合わせて構成されることを特徴とするものである。
このように構成されたリアルタイム震度計測装置においては、ＩＩＲフィルター特性として、特に自由度を必要とするゲインピークを含む周波数帯域で、遮断周波数と減衰定数の２つの自由度を有する２次ＨＰＦと遮断周波数の１つの自由度を有する１次ＨＰＦを組み合わせて第１のＩＩＲフィルターを構成させることで、必要とされる高い自由度を提供するという作用を有する。
また、請求項１記載のリアルタイム震度計測装置は、前記ＩＩＲフィルターが、前記ゲインピークを含む周波数帯域以下の低周波数帯域においては、前記第１のＩＩＲフィルターを用い、前記ゲインピークを含む周波数帯域よりも高い１Ｈｚ＜ｆ≦３Ｈｚを含む中周波数帯域においては、遮断周波数の１つの自由度を有する２つの１次フィルターを組み合わせて周波数の平方根の逆数に比例する第２のＩＩＲフィルターを備え、前記中周波数帯域よりも高い３Ｈｚ＜ｆを含む高周波数帯域においては、遮断周波数と減衰定数の３つの自由度を有する２次ローパスフィルター（ＬＰＦ）を第３のＩＩＲフィルターとして備えていることを特徴とするものである。
このように構成されたリアルタイム震度計測装置においては、ＩＩＲフィルターによってすべての周波数帯域において、地震加速度波形の特徴に沿ってフィルター処理されるという作用を有する。

本発明の請求項１に記載のリアルタイム震度計測装置によれば、リアルタイムで高精度の震度計測が可能となるので、地震速報をはじめ緊急地震情報をより迅速に高い精度で発信することができる。また、従来困難であった震源に近い地域での迅速な対応、避難が可能となり、住民の安全確保や財産の保全がより確実に実行可能となる。

東日本大震災を引き起こした２０１１年東北地方太平洋沖地震の際には、短時間に精度の高い地震の規模（マグニチュード）が求まらず、そのために、緊急地震速報や津波警報が適切に発表できないという不具合があった。
これらの不具合は、本発明の請求項１記載のリアルタイム震度計測装置により、国土の震度分布を精度良くリアルタイムに計測することで、解決できる。
例えば、「地震のマグニチュード」は従来の震源情報と本発明によるリアルタイム震度分布を組み合わせることで、短時間に精度よく求めることが可能となる。「緊急地震速報」に適用すれば、震源に近い観測点のリアルタイム震度を用いて、従来は困難であった直下地震の震央周辺でも、緊急地震速報の利用が可能となり、震央から離れた地域でも情報の迅速化と高精度化が可能になる。「津波警報」については、マグニチュードに基づく津波警報の迅速化と高精度が実現する他、津波被害が想定される沿岸部のリアルタイム震度が求まるので、リアルタイム震度に基づく津波警報の要否・警報のレベルも判定できるので、津波警報の信頼性を高めることができる。
このように、本発明を応用すると、地震に関する諸情報の発表が迅速化でき、情報の精度が向上して、地震直後のすべての行動に対する時間的余裕が従来に比べると増加し、被災による被害を最小限に留める基本的な効果を発揮することができるのである。

（ａ）は、本発明の実施の形態に係るリアルタイム震度計測装置の構成図であり、（ｂ）はＩＩＲフィルターの構成図である。 気象庁が定めるＦＦＴフィルターの周波数特性図である。 気象庁が定めるＦＦＴフィルターの特性を周波数帯域毎に分解して説明する図であり、（ａ）は低周波数成分の遮断フィルターJ₁(f )の周波数特性図、（ｂ）は周波数の平方根に逆比例する特性を備えたフィルターJ₂(f )の周波数特性図、（ｃ）は高周波成分の遮断フィルターJ₃(f )の周波数特性図である。 本実施の形態に係るリアルタイム震度計測装置の第１のＩＩＲフィルターについて、入力時系列とフィルター出力時系列との関係を示すブロック図である。 本実施の形態に係るリアルタイム震度計測装置の第２のＩＩＲフィルターについて、入力時系列とフィルター出力時系列との関係を示すブロック図であり、（ａ）はy₁(t )、（ｂ）はy₂(t)（ｃ）はy₃(t )について示している。 本実施の形態に係るリアルタイム震度計測装置の第３のＩＩＲフィルターについて、入力時系列とフィルター出力時系列との関係を示すブロック図である。 本実施の形態に係るリアルタイム震度計測装置のＩＩＲフィルターG (z)の周波数特性図である。 本実施の形態に係るリアルタイム震度計測装置のＩＩＲフィルターG (z)の周波数特性と、気象庁の定義するフィルターJ の周波数特性とを比較して示す周波数特性図である。 本実施の形態に係るリアルタイム震度計測装置のＩＩＲフィルターG (z)と特許文献１に開示されている従来技術のフィルタの周波数特性とを比較して示す周波数特性図である。 本実施の形態に係るリアルタイム震度計測装置のＩＩＲフィルターを用いて計算される計測震度相当値と気象庁計測震度との比較を示す評価図である。 本実施の形態に係るリアルタイム震度計測装置のＩＩＲフィルターを用いた場合の震度誤差の分布図である。 本実施の形態に係るリアルタイム震度計測装置のＩＩＲフィルターを用いた場合の震度誤差と震央距離の関係を示す評価図である。 本実施の形態に係るリアルタイム震度計測装置のＩＩＲフィルターを用いた場合の震度誤差と計測震度との関係を示す評価図である。 本実施の形態に係るリアルタイム震度計測装置のＩＩＲフィルターを用いた場合の震度誤差と地震のマグニチュードとの関係を示す評価図である。 特許文献１に開示された技術によるＩＩＲフィルターを用いた場合の震度誤差の分布図である。 特許文献１に開示された技術によるＩＩＲフィルターを用いた場合の震度誤差と震央距離の関係を示す評価図である。 特許文献１に開示された技術によるＩＩＲフィルターを用いた場合の震度誤差と計測震度との関係を示す評価図である。 特許文献１に開示された技術によるＩＩＲフィルターを用いた場合の震度誤差と地震のマグニチュードとの関係を示す評価図である。

以下に、本発明の第１の実施の形態に係るリアルタイム震度計測装置について図１乃至図１４を参照しながら説明する。
図１（ａ）は、本発明の第１の実施の形態に係るリアルタイム震度計測装置の構成図であり、（ｂ）はＩＩＲフィルターの構成図である。
図１（ａ）において、リアルタイム震度計測装置１は、大きくは演算部２と計時部３及び出力部４から構成されている。
また、演算部２は、加速度計５、Ａ／Ｄ変換器６、ＩＩＲフィルター７、震度演算部８から構成されている。
加速度計５は、地震による揺れから加速度成分を検出するものであり、この加速度計５で検出されたアナログの加速度波形信号５Ｓ（アナログ信号）を、計時部３から発信される連続した時間信号３Ｓ１で一定時間ごとにサンプリングすることにより離散化された加速度信号６Ｓ（デジタル信号）に変換するＡ／Ｄ変換器６を備えている。
また、このＡ／Ｄ変換器６から出力されるデジタルの加速度信号を周波数毎に重み付けを行ってデジタルの重み付け加速度信号７Ｓを出力するＩＩＲフィルター７を備えており、本実施の形態では、図１（ｂ）に示されるとおり、このＩＩＲフィルター７はさらに、第１のＩＩＲフィルター７１、第２のＩＩＲフィルター７２及び第３のＩＩＲフィルター７３から構成されている。すなわち、ＩＩＲフィルター７は気象庁が定めるＦＦＴフィルターの近似フィルターであり、第１のＩＩＲフィルター７１、第２のＩＩＲフィルター７２、第３のＩＩＲフィルター７３はそれぞれＩＩＲフィルター７の一部を構成するものである。
第１のＩＩＲフィルター７１は、周波数が１Ｈｚよりも低い周波数帯域において処理するフィルターであり、第２のＩＩＲフィルター７２は１Ｈｚ以上で３Ｈｚ未満の周波数帯域において処理するフィルターであり、第３のＩＩＲフィルター７３は、３Ｈｚ以上の周波数帯域で処理するフィルターである。これらのフィルターの特性については詳細に後述する。

ＩＩＲフィルター７で周波数毎に重み付けをされた重み付け加速度信号７Ｓを用いて、震度演算部８では震度を算出し、震度信号８Ｓを出力する。
計時部３は、前述のとおりサンプリングするための時間信号３Ｓ１を発信するのみならず、演算部２を構成する加速度計５、Ａ／Ｄ変換器６、ＩＩＲフィルター７、震度演算部８において出力される信号に対して、互いに共通の時間軸として関連付けるための時刻信号３Ｓ２を出力する。
従って、加速度計５、Ａ／Ｄ変換器６、ＩＩＲフィルター７及び震度演算部８によってそれぞれ出力される加速度波形信号５Ｓ、加速度信号６Ｓ、重み付け加速度信号７Ｓ及び震度信号８Ｓにおいても時刻信号３Ｓ２に重畳して出力され、最終的に得られる震度信号８Ｓによって示されるリアルタイム震度は、時刻に沿って把握することが可能である。
出力部４は、少なくとも震度演算部８から出力される震度信号８Ｓを入力して外部へ出力するものであるが、加速度計５から出力される加速度波形信号５Ｓ、Ａ／Ｄ変換器６から出力される加速度信号６Ｓ及びＩＩＲフィルター７から出力される重み付け加速度信号７Ｓについても入力して外部へ出力すると、それぞれの構成要素におけるデータ処理の検証も可能であるため便利である。

特に、加速度計５から出力される加速度波形信号５Ｓは地震の加速度波形を直接的に把握することができ、例えば震度信号８Ｓと重畳させて出力部４から出力させることで、地震の加速度とリアルタイムの震度について比較を行うことが可能であり、地震加速度と震度のリアルタイムの相関を把握することが可能であり、防災計画を立てる際にも有益な情報となる。
出力部４としては、具体的にはＣＲＴ、液晶、プラズマあるいは有機ＥＬなどによるディスプレイ装置、あるいはプリンタ装置などの出力装置、さらには外部装置への伝送を行うためのトランスミッタなどの発信装置などが考えられる。もちろん、外部装置への伝送のための出力に対するインターフェースであってもよい。
なお、図１（ａ）、（ｂ）を用いて本発明のリアルタイム震度計測装置１を説明したが、このリアルタイム震度計測装置１をリアルタイム震度計測方法というような方法発明として捉えることも可能である。
すなわち、加速度計５は、地震による加速度を検出してアナログ信号である加速度波形信号５Ｓを出力する加速度検出工程と捉えることができる。また、Ａ／Ｄ変換器６は、加速度検出工程から出力される加速度波形信号５Ｓをデジタル信号である加速度信号６Ｓに変換するＡ／Ｄ変換工程と捉えることができる。さらに、ＩＩＲフィルター７は、Ａ／Ｄ変換工程から出力される加速度信号６Ｓをフィルター処理してデジタル信号である重み付け加速度信号７Ｓを出力するＩＩＲフィルタリング工程と捉えることができる。震度演算部８は、ＩＩＲフィルタリング工程から出力される重み付け加速度信号７Ｓに基づいてリアルタイム震度を計算する震度演算工程と捉えることができ、この震度演算工程で計算されたリアルタイム震度は震度信号８Ｓとして出力され、出力部４は、この震度信号８Ｓを受信して外部に出力する出力工程として捉えることができるのである。
また、リアルタイム震度計測装置１における計時部３から発信される時間信号３Ｓ１と時刻信号３Ｓ２に関する作用、効果についても、リアルタイム震度計測装置１を構成する加速度計５、Ａ／Ｄ変換器６、ＩＩＲフィルター７、震度演算部８、出力部４などの要素に対応する、それぞれ加速度検出工程、Ａ／Ｄ変換工程、ＩＩＲフィルタリング工程、震度演算工程、出力工程において同様に作用、効果が発揮される。

背景技術の欄でも述べたとおり、地震震度は、加速度センサによって測定される地震加速度の時間軸上の波形データを一定時間蓄積した上で、ＦＦＴ（高速フーリエ変換）演算を気象庁が定めるＦＦＴフィルターで処理して行った後、それを逆ＦＦＴ変換して得られた加速度波形データを用いて計算される。
この気象庁が定めるＦＦＴフィルターの周波数特性図を図２に示し、図３に、このＦＦＴフィルターの特性を周波数帯域毎に分解して説明する図を示す。図３（ａ）は低周波数成分の遮断フィルターJ₁(f )の周波数特性図、（ｂ）は周波数の平方根に逆比例する特性を備えたフィルターJ₂(f )の周波数特性図、（ｃ）は高周波成分の遮断フィルターJ₃(f )の周波数特性図である。
図２において、横軸は地震波形の周波数（Ｈｚ）であり、縦軸は相対利得（倍）を示している。このように定められるＦＦＴフィルターの周波数特性は、地震震度という自然現象の程度を物理的に示すためのものであり、それ自体には物理的な意味があるものである。
この気象庁が定めるＦＦＴフィルターに代えて、近似した特性を発揮するＩＩＲフィルター７を備えることでリアルタイム震度計測装置１が実現されるのである。

ここで、このＦＦＴフィルターの定義について説明する。
加速度記録から気象庁計測震度を求めるときに用いられるフィルターJ (f )は、周波数帯域で以下のように定義される。
J(f)=J₁(f)J₂(f)J₃(f) （１）
ここで、J₁(f)、J₂(f)、J₃(f)は以下のとおりである。
J₁(f)=[1-exp{-(f/0.5)³}]^1/2 （２）
J₂(f)=(1/f)^1/2 （３）
J₃(f)=1/{1+Sc_k(f/10)^2k}^1/2,k=1,...,6 （４）
ｆは周波数、c_kは定数で、c₁=0.694, c₂=0.241, c₃=0.0557, c₄=0.009664,
c₅=0.00134, c₆=0.000155

さらに、これらJ₁(f )、J₂(f )、J₃(f)及びJ (f)の特性について詳細に説明する。
（１）J₁(f)の特性（図３（ａ）参照）
フィルターJ₁(f)は、遮断周波数０．５Ｈｚの１．５次のＨＰＦである。およそ０．３Ｈｚより低周波数側では、振幅は周波数の１．５乗に比例し、０．３〜１Ｈｚ の周波数帯域では、振幅の増加率は徐々に滑らかに低くなり、およそ１Ｈｚより高周波数側では平坦となる。
（２）J₂(f)の特性（図３（ｂ）参照）
フィルターJ₂(f)は、振幅が周波数の平方根の逆数に比例する。
（３）J₃(f)の特性（図３（ｃ）参照）
フィルターJ₃(f)は、遮断周波数１０Ｈｚの６次のＬＰＦである。およそ３Ｈｚより低周波数側では、振幅は平坦であり、３〜３０Ｈｚの周波数帯域では、振幅は徐々に滑らかに減少し、およそ３０Ｈｚより高周波数側では周波数の６乗に逆比例する。

（４）J (f)の特性（図２参照）
フィルターJ (f)の特性は、以下のとおりである。
（ａ）ｆ≦０．３Ｈｚの低い低周波数帯域
ｆ≦０．３Ｈｚの低い周波数帯域では、フィルターJ₃(f)は平坦な特性をもち、１．５次ＨＰＦ（高域通過フィルター）であるフィルターJ₁(f)による特性と周波数の平方根に逆比例するフィルターJ₂(f)による特性の重畳となることから、利得が周波数に比例する特性をもつ。すなわち、1次のＨＰＦと同様の特性をもつ。
（ｂ）０．３＜ｆ≦１Ｈｚの帯域
０．３＜ｆ≦１Ｈｚの帯域では、フィルターJ₃(f)は平坦な特性、フィルターJ₂(f)は周波数の平方根に逆比例する特性と、どちらも単純な特性をもつ。この帯域において、周波数に比例して増加する振幅特性から、徐々に滑らかに増加率が低下し、およそ０．６Ｈｚで極大値をとり、０．６Ｈｚより高周波数側では徐々に滑らかに振幅減少率が高くなり、およそ１Ｈｚで周波数の平方根に逆比例する特性となる特徴は、フィルターJ₁(f)によるものである。すなわち、この帯域にゲインピークが存在している。
（ｃ）１＜ｆ≦３Ｈｚの中周波数帯域
１＜ｆ≦３Ｈｚの中周波数帯域では、フィルターJ₁(f)、及びフィルターJ₃(f)はともに平坦な特性をもつことから、フィルターJ₂(f)により周波数に逆比例する特性をもつ。
（ｄ）３＜ｆ≦３０Ｈｚの帯域
３＜ｆ≦３０Ｈｚの帯域では、フィルターJ₁(f)は平坦な特性、フィルターJ₂(f)は周波数の平方根に逆比例する特性と、どちらも単純な特性をもつ。この帯域において、徐々に滑らかに振幅の減少率が高くなり、およそ３０Ｈｚで周波数の６．５乗に逆比例する特性となる特徴は、フィルターJ₃(f)によるものである。
（ｅ）３０Ｈｚ＜ｆの高周波数帯域
３０Ｈｚより高い周波数帯域では、フィルターJ₁(f)は平坦な特性をもち、周波数の平方根に逆比例するフィルターJ₂(f)による特性と６次ＬＰＦ（低域通過フィルター）であるフィルターJ₃(f)の特性の重畳となることから、６．５次のＬＰＦと同様の特性をもつ。なお、本願においては、周波数帯域の上下限を表現する際に、境界を明確化すべく、便宜上Ａ＜ｆ≦Ｂという表現を用いたが、Ａ≦ｆ＜Ｂであってもよく、上下限のいずれの数値を範囲に含めるようにするかは特に限定するものではない。

以上のような特性を有する気象庁のＦＦＴフィルターに近似させるべく、本実施の形態に係るリアルタイム震度計測装置１のＩＩＲフィルター７では以下に説明するように３つの周波数帯域に分けて周波数特性を持たせた。この３つの周波数帯域に分けた周波数特性を持つそれぞれのフィルターが、図１（ｂ）に示す第１のＩＩＲフィルター７１、第２のＩＩＲフィルター７２、第３のＩＩＲフィルター７３に対応している。
（ａ）ｆ≦０．３Ｈｚの低周波数帯域及び０．３＜ｆ≦１Ｈｚの周波数帯域
０．３Ｈｚ以下の低周波数帯域においては、フィルターJ (f)は1次で減衰する特性をもつ。この帯域においては１次で減衰するＨＰＦで近似することができる。１次で減衰するＨＰＦは、遮断周波数より十分低周波数側では振幅は１次で減衰し、遮断周波数周辺より十分高周波数側では平坦であり、遮断周波数周辺の帯域では、遷移的特性をもつ。このことを利用して、ｆ≦０．３Ｈｚの低い周波数帯域及び０．３＜ｆ≦１Ｈｚの帯域における特性を、１次で減衰するＨＰＦで近似する。
１次で減衰するＨＰＦとして、第１のＩＩＲフィルター７１を本実施の形態では、（ｎ＋１）次ＨＰＦとｎ次ＨＰＦとの比（ｎは正整数）で実現させた。１次ＨＰＦそのものは自由度が１であり、自由に設定できる条件は遮断周波数のみである。従って、遮断周波数周辺の帯域での特性は、１次ＨＰＦ固有の特性に固定される。震度計算に用いられる前述のフィルターJ (f)の遮断周波数周辺における周波数特性は、１次ＨＰＦとは大きく異なり、遮断周波数周辺におけるフィルターJ (f)の振幅は、１次ＨＰＦの振幅より有意に大きい。周波数特性の再現を考慮すると、１次ＨＰＦ単独での近似は不適切であると考えられる。

一方、（ｎ＋１）次ＨＰＦとｎ次ＨＰＦとの比は自由度が２ｎ＋１であり、ｎ（整数）を適切に選択することで、フィルターJ (f)の遮断周波数周辺における周波数特性を近似させることが可能である。発明者は、フィルターJ (f)の遮断周波数周辺における特性を分析したうえで、ここでは一例としてｎ＝１を採用し、自由度３の２次と１次の比で十分な近似が得られることを確認した。すなわち第１のＩＩＲフィルター７１としてG_a(f)を遮断周波数f_a2、減衰定数h_aの２次ＨＰＦと遮断周波数f_a1の１次ＨＰＦを組み合わせて以下のようにする。
G_a(f)=H₂(f;f_a2,h_a,g_a)/H₁(f;f_a1,g_a),
（５）
ここで、f は周波数、f_a1、f_a2は遮断周波数、h_aは減衰定数、g_a1、g_a2は利得で、g_a1=g_a2=１とする。式（５）の中のH₂(f;f_a2,h_a,g_a)は、遮断周波数f₀、減衰定数h、利得gの２次ＨＰＦ、H₁(f;f_a1,g_a)は、遮断周波数f₀、利得gの１次ＨＰＦで、以下のように表される。
H₂(f;f₀,h,g)=g/{(f₀/jf)²+2h(f₀/jf)+1} （６）
H₁(f;f₀,g)=g/{(f₀/jf)+1} （７）
ここで、j=（-1）＾(1/2)は虚数単位、すなわち、（−１）の平方根である。
このようにｎ＝１を採用して自由度が３となったことで、２次ＨＰＦの遮断周波数f_a2、減衰定数h_aと１次ＨＰＦの遮断周波数f_a1という３つの定数をとることが可能であり、これによってＦＦＴフィルターにおける０．３＜ｆ≦１Ｈｚの帯域で示される特性、具体的には特に、およそ０．６Ｈｚで極大値、すなわちゲインピークをとるという特性が表現できるのである。
従って、リアルタイム震度計測装置１をリアルタイム震度計測方法として捉えた場合には、ＩＩＲフィルタリング工程で、ゲインピークを含む周波数帯域において、（ｎ＋１）次ハイパスフィルター（ＨＰＦ）（ｎは正整数、以下同様）とｎ次ハイパスフィルター（ＨＰＦ）を組み合わせてフィルタリングすることを意味している。また、その際には、ｎ＝１を採用して自由度を３としてもよい。
さらに、ＩＩＲフィルター７が第１のＩＩＲフィルター７１、第２のＩＩＲフィルター７２及び第３のＩＩＲフィルター７３を組み合わせたフィルターとして表現できることから、ＩＩＲフィルタリング工程においては、これらの３つのフィルターを組み合わせて用い、このうち、第１のＩＩＲフィルター７１において、ゲインピークを含む周波数帯域において、（ｎ＋１）次ハイパスフィルター（ＨＰＦ）（ｎは正整数、以下同様）とｎ次ハイパスフィルター（ＨＰＦ）を組み合わせてフィルタリングし、その他の第２のＩＩＲフィルター７２と第３のＩＩＲフィルター７３については、これから述べる特性を備えたフィルターとしてフィルタリングするようにするとよい。

（ｂ）１＜ｆ≦３Ｈｚの中周波数帯域
周波数１＜ｆ≦３Ｈｚ程度の範囲で周波数の平方根の逆数に比例するフィルターJ (f )の特性を実現するために、第２のＩＩＲフィルター７２としてG_b(f )を以下のように1次HPFを組み合わせる。
G_b(f)=G_b1(f)G_b2(f)G_b3(f) （８）
G_b1(f)=G_b12(f)/G_b11(f) （９）
G_b2(f)=G_b22(f)/G_b21(f) （１０）
G_b3(f)=G_b32(f)/G_b31(f) （１１）
G_b11(f)=H₁(f;f_b,g_b1) （１２）
G_b12(f)=H₁(f;f_b/η,g_b1/η) （１３）
G_b21(f)=H₁(f;γf_b,g_b2) （１４）
G_b22(f)=H₁(f;f_b/γη,g_b2/γη) （１５）
G_b31(f)=H₁(f;f_b/γ,g_b3) （１６）
G_b32(f)=H₁(f;γf_b/η,γg_b3/η) （１７）
ここで、f_bは遮断周波数、γ、ηは減衰率と帯域幅を決定する定数、またg_b1、g_b2、g_b3は利得で、g_b1=g_b2=g_b3=1とする。このフィルターの振幅は、遮断周波数より十分低周波数側、及び十分高周波数側では平坦である。

（ｃ）３＜ｆ≦３０Ｈｚの帯域及び３０Ｈｚ＜ｆの高周波数帯域
フィルターJ (f )は高周波数帯域において６．５次で減衰する特性をもつ。この帯域においては、単純には６．５次のＬＰＦが考えられるが、後述するように、周波数fと変数zとの変数変換の近似の影響により、高周波数帯域においては、デジタルＬＰＦフィルターの減衰率はＬＰＦの次数以上となる。このことを利用して、高周波数帯域においては２次ＬＰＦを採用することとする。また、実際の地震動は伝播に伴う減衰のために高周波数帯域では振幅は小さいことから、２次ＬＰＦによる近似で実用上問題がないことが確認される。
２次ＬＰＦは遮断周波数より十分低周波数側では振幅は平坦であり、遮断周波数周辺においては徐々に滑らかに減少する。従って、フィルターJ(f)の３＜ｆ≦３０Ｈｚにおける特性も、本実施の形態に係るリアルタイム震度計測装置１の第３のＩＩＲフィルター７３としてのG_c(f)では２次ＬＰＦで近似することとする。
G_c(f)=L₂(f;f_c,h_c,g_c) （１８）
ここで、f_cは遮断周波数、h_cは減衰定数、g_cは利得で、g_c=1とする。
L₂(f;f_c,h_c,g_c)は２次ＬＰＦで、以下のように表される。
L₂(f;f_c,h_c,g_c)=g_c(f_c/jf)²/{(f_c/jf)²+2h_c(f_c/jf)+1} （１９）

次に、上記の第１のＩＩＲフィルター７１（G_a(f)）、第２のＩＩＲフィルター７２（G_b(f)）、第３のＩＩＲフィルター７３（G_c(f)）に現れる１次ＨＰＦ、２次ＨＰＦ、及び２次ＬＰＦは、以下のようにラプラス変換の表現で書き換えられる。
H₁(s;f₀,g)=g/{(2πf₀/s)+1} （２０）
H₂(s;f₀,h,g)=g/{(2πf₀/s)²+4πhf₀/s+1} （２１）
L₂(s;f₀,h,g)=g(2πf₀/s)²/{(2πf₀/s)²+4πhf₀/s+1} （２２）
ここで、sはラプラス変数で、
s=2πjf （２３）
である。
次に、上記の式（２０）、（２１）、（２２）をｚ変換による表現に書き換える。ラプラス変数ｓと変数ｚとの関係は以下で表される。Δtはサンプリング時間間隔である。
z=exp(sΔt) （２４）
あるいは、
s=(1/Δt)lnz （２５）
上式の ｌｎｚを級数展開することにより以下の関係が得られる。
s^-1=(Δt/2)(1+z^-1)/(1-z^-1) （２６）
s^-2=(Δt²/12)(1+10z^-1+z^-2)/(1-z^-1)² （２７）
１次ＨＰＦ（２０）は、式（２６）を用いて、以下のように書き換えられる。
H₁(z;f₀,g)=g(1-z^-1)/{(1+πf₀Δt)+(-1+πf₀Δt)z^-1} （２８）
２次ＨＰＦ（２１）、２次ＬＰＦ（２２）は、以下の２通りの方法でｚ変換に書き換えられる。
(A) 式（２６）で示される関係
(B) 式（２６）及び式（２７）で示される関係
２次ＨＰＦの低周波数帯域での特性については（Ａ）あるいは（Ｂ）いずれを用いても周波数特性の相違はほとんどないが、２次ＬＰＦの高周波数帯域での特性については、（Ａ）の方法によるほうが遮断周波数以上の帯域での減衰率が（Ｂ）の方法によるよりも大きい。このことを利用して、１０Ｈｚ以上の帯域での気象庁計測震度フィルターの特性をよりよく近似できるように、ここでは（Ａ）の方法を採用することとする。２次ＨＰＦ（２１）は以下のように書き換えられる。
H₂(z;f₀,h,g)=g(1-2z^-1+z^-1)/H(z;f₀,h) （２９）
H(z;f₀,h)={1+2πhf₀Δt+(πf₀Δt)²}+{-2+2(πf₀Δt)²}z^-1
+{1-2πhf₀Δt+(πf₀Δt)²}z^-2 （３０）
また、２次ＬＰＦ（２２）は以下のように書き換えられる。
L₂(z;f₀,h,g)=g(πf₀Δt)²(1+2z^-1+z^-2)/H(z;f₀,h) （３１）

以下、第１のＩＩＲフィルター７１（G_a(f)）、第２のＩＩＲフィルター７２（G_b(f)）、第３のＩＩＲフィルター７３（G_c(f)）のｚ変換による表現について説明する。

（ａ）第１のＩＩＲフィルター７１（G_a(f)）のｚ変換による表現
式（５）に式（２０）、（２１）、（２６）を代入し整理すると次式を得る。
G_a(z)={β_a0+β_a1z^-1+β_a2z^-2}/{α_a0+α_a1z^-1+α_a2z^-2} （３２）
ここで、β_a0=1+ζ_a1，β_a1=-2，β_a2=1-ζ_a1，α_a0=1+2h_aζ_a2+ζ_a2 ²，
α_a1=2ζ_a2 ²-2，α_a2=1-2h_aζ_a2+ζ_a2 ²，ζ_a1=πf_a1Δt，ζ_a2=πf_a2Δt

（ｂ）第２のＩＩＲフィルター７２（G_b(f)）のｚ変換による表現
式（６）に式（２０）、（２６）を代入し整理すると次式を得る。
G_b1(z)=G_b12(z)/G_b11(z)={ β_b10+β_b11z^-1}/{α_b10+α_b11z^-1} （３３）
G_b2(z)=G_b22(z)/G_b21(z)={ β_b20+β_b21z^-1}/{α_b20+α_b21z^-1} （３４）
G_b3(z)=G_b32(z)/G_b31(z)={ β_b30+β_b31z^-1}/{α_b30+α_b31z^-1} （３５）
ここで、β_b10=ζ_b+1，β_b11=ζ_b-1，α_b10=ζ_b+η，α_b11=ζ_b-η，
β_b20=γζ_b+1，β_b21=γζ_b-1，α_b20=ζ_b+γη，α_b21=ζ_b-γη，
β_b30=ζ_b/γ+1，β_b31=ζ_b/γ-1，α_b30=ζ_b+η/γ，α_b31=ζ_b-η/γ，ζ_b=πf_bΔt

（ｃ）第３のＩＩＲフィルター７３（G_c(f)）のｚ変換による表現
式（７）に式（２２）、（２６）を代入し整理すると次式を得る。
G_c(z)={β_c0+β_c1z^-1+β_c2z^-2}/{α_c0+α_c1z^-1+α_c2z^-2} （３６）
ここで、β_c0=ζ_c ²，β_c1=2ζ_c ²，β_c2=ζ_c ²，α_c0=1+2h_cζ_c+ζ_c ²，
α_c1=2ζ_c ²-2，α_c2=1-2h_cζ_c+ζ_c ²，ζ_c=πf_cΔt

次に、第１のＩＩＲフィルター７１（G_a(f)）、第２のＩＩＲフィルター７２（G_b(f)）、第３のＩＩＲフィルター７３（G_c(f)）のｚ変換による表現を利用して、それぞれの入力時系列x(t)とＩＩＲフィルター７の出力y(t)の関係について説明する。
入力時系列をx(t)、フィルター出力をy(t)とすると、
y(z)=G(z)x(z) （３７）
また、
y(z)z^-1 −＞ y(t-Δt) （３８）
式（３２）、（３３）〜（３５）、（３６）と、式（３７）に基づいて、また式（３８）を用いて、ＩＩＲフィルター７の時間領域における表現を得る。

（ａ）第１のＩＩＲフィルター７１（G_a(f)）
式（３２）と、式（３７）、及び式（３８）から、
y(t)=b_a0x(t)+b_a1x(t-Δt)+b_a2x(t-2Δt)-a_a1y(t-Δt)-a_a2y(t-2Δt) （３９）
を得る。
ここで、b_a0=β_a0/α_a0，b_a1=β_a1/α_a0，b_a2=β_a2/α_a0，a_a1=α_a1/α_a0，
a_a2=α_a2/α_a0
本実施の形態に係るリアルタイム震度計測装置１の第１のＩＩＲフィルター７１の入力時系列とフィルター出力時系列との関係は式（３９）のように表現でき、これらの関係は図４に示されるとおりである。

（ｂ）第２のＩＩＲフィルター７２（G_b(f)）
式（３３）〜（３５）と、式（３７）、及び式（３８）から、
y₁(t)=b_b10x(t)+b_b11x(t-Δt)-a_b11y₁(t-Δt) （４０）
y₂(t)=b_b20y₁(t)+b_b21y₁(t-Δt)-a_b21y₂(t-Δt) （４１）
y₃(t)=b_b30y₂(t)+b_b31y₂(t-Δt)-a_b31y₃(t-Δt) （４２）
を得る。
ここで、b_b10=β_b10/α_b10，b_b11=β_b11/α_b10，a_b11=α_b11/α_b10，
b_b20=β_b20/α_b20，b_b21=β_b21/α_b20，a_b21=α_b21/α_b20，
b_b30=β_b30/α_b30，b_b31=β_b31/α_b30，a_b31=α_b31/α_b30
本実施の形態に係るリアルタイム震度計測装置１の第２のＩＩＲフィルター７２の入力時系列とフィルター出力時系列との関係は、式（４０）から（４２）のように表現でき、これらの関係は図５（ａ）〜（ｃ）に示すとおりである。すなわち、（ａ）はy₁(t)、（ｂ）はy₂(t)（ｃ）はy₃(t)についてそれぞれ示すブロック図である。

（ｃ）第３のＩＩＲフィルター７３（G_c(f)）
式（３６）と、式（３７）、及び式（３８）から、
y(t) = b_c0x(t)+b_c1x(t-Δt)+b_c2x(t-2Δt)-a_c1y(t-Δt)-a_c2y(t-2Δt) （４３）
を得る。
ここで、b_c0=β_c0/α_c0，b_c1=β_c1/α_c0，b_c2=β_c2/α_c0， a_c1=α_c1/α_c0，
a_c2=α_c2/α_c0
本実施の形態に係るリアルタイム震度計測装置１の第３のＩＩＲフィルター７３の入力時系列とフィルター出力時系列は式（４３）のように表現でき、この関係は図６に示すとおりである。

これまで説明したことをまとめると、計測震度を求めるときに用いられるフィルターを近似するＩＩＲフィルター７（G(f)）は、ｚ返還で表現すれば、Ga(z)、Gb1(z)、Gb2(z)、Gb3(z)、Gc(z)及び利得g_Tの積として式（４４）のように表現することができる。ここでは、演算の順序も含めて式（４４）のように表現している。
G(z)=Gc(z)Ga(z)Gb1(z)Gb3(z)g_TGb2(z) （４４）
フィルターG (z)の周波数特性は、f_a1、f_a2、h_a、f_b、γ、η、g_T、f_c、h_cの９つの係数で決定される。これらの係数については、式（６）、（７）、（１２）〜（１７）、（１９）に記載されるとおりである。
フィルターG (z)が、周波数０．０１〜２０Ｈｚの帯域で、フィルターJ を最もよく近似するように、これら９つの係数を最小二乗法により決定した。この決定に関しては後述する。

図７に、本実施の形態に係るＩＩＲフィルター７、すなわち、第１のＩＩＲフィルター７１（G_a(z)）、第２のＩＩＲフィルター７２（G_b(z)）、第３のＩＩＲフィルター７３（G_c(z)）から構成されるG (z )の周波数特性を示す。横軸は周波数（Ｈｚ）を示し、横軸は相対利得（倍）を示す。図７から理解されるとおり、図２に示される気象庁が定義するＦＦＴフィルターの特性と精度よく近似されていることがわかる。

また、図８に本実施の形態に係るリアルタイム震度計測装置１のＩＩＲフィルター７としてのG(z)と気象庁の定義するフィルターJの周波数特性との比較を示す。横軸と縦軸は図７と同一である。周波数０．０１〜２０Ｈｚの帯域で、フィルターG(z)の特性はフィルターJのものとほぼ一致している。
さらに、図９に本実施の形態に係るリアルタイム震度計測装置１のＩＩＲフィルター７としてのG(z)と特許文献１で開示されている従来技術のフィルタの周波数特性との比較を示す。図９において、実線は本実施の形態に係るリアルタイム震度計測装置１のＩＩＲフィルター７の周波数特性を示しており、点線は特許文献１に係る従来技術のフィルタの周波数特性を示している。これらを比較すると、０．６Ｈｚ近傍での差が最も大きい。
これは、従来技術がゲインピークを含む周波数帯域において、単に１次で減衰するＨＰＦを用いて、遮断周波数を０．６Ｈｚ近傍に持たせることで、この遮断周波数より十分低周波数側では振幅は１次で減衰し、遮断周波数周辺より十分高周波数側では平坦であり、遮断周波数周辺の帯域では、遷移的特性を持つようにしたものの、この０．６Ｈｚ近傍で気象庁の定義するフィルターJの周波数特性を十分に表現できず、結局本願実施の形態に係るリアルタイム震度計測装置１のＩＩＲフィルター７との周波数特性とも差異が大きく出てしまったものである。
一方、本実施の形態に係るＩＩＲフィルター７においては、ゲインピークを含む周波数帯域、すなわち０．３Ｈｚ＜ｆ≦１Ｈｚの範囲を含む帯域において、第１のＩＩＲフィルター７１を２次ハイパスフィルター（ＨＰＦ）と１次ハイパスフィルター（ＨＰＦ）を組み合わせて構成している。この第１のＩＩＲフィルター７１では、このような構成によって自由度が３となったことで、２次ＨＰＦの遮断周波数f_a2、減衰定数h_aと１次ＨＰＦの遮断周波数f_a1という３つの定数をとることが可能であり、これによってＦＦＴフィルターにおける０．３＜ｆ≦１Ｈｚの帯域で示される特性、具体的には特に、およそ０．６Ｈｚで極大値をとる特性が気象庁の定義するフィルターJとほぼ同一の特性として表現できている。従って、ゲインピークを含む周波数帯域において、図８ではほぼ同一の周波数特性として表れているのである。

気象庁の定義するフィルターJをよい精度で近似する本実施の形態に係るリアルタイム震度計測装置１のＩＩＲフィルター７G (z)を用いて計算される計測震度相当値と気象庁による計測震度との比較を図１０に示す。図１０における横軸は計測震度、縦軸はＩＩＲフィルター７を採用した本実施の形態に係るリアルタイム震度計測装置１において計算された計測震度相当値である。すなわち、リアルタイムに得られた震度である。また、計測震度相当値の気象庁計測震度からのずれの分布を図１１に示す。図１１における横軸は震度誤差、すなわちＩＩＲフィルター７を用いて得られた計測震度相当値から気象庁計測震度を差し引いた値である。縦軸は相対頻度を示している。
図１０、図１１では、近年日本国内および近海で発生した２５２の地震を国内の地震観測点で観測した４３，３９０の地震記録を用いた。本実施の形態に係るＩＩＲフィルター７を用いて計算される計測震度相当値は、計測震度の広い範囲で気象庁計測震度と十分な精度で一致していることが確認できる。ずれの平均値はほぼ０であり、ずれの幅（分布の標準偏差）は非常に小さい。

図１２、図１３、図１４に、震度誤差と、それぞれ震央距離、計測震度（気象庁計測震度）、地震のマグニチュード（気象庁によるマグニチュード）との関係を示す。図１２における横軸は震央距離、図１３における横軸は計測震度、図１４における横軸は気象庁によるマグニチュード値であり、図１２から図１４の縦軸はいずれも図１１で横軸に示されていた震度誤差である。これらの図からわかるように、気象庁の定義するフィルターJを高精度で近似するＩＩＲフィルター７G (z)を用いて計算される計測震度相当値は、震央距離、計測震度、地震のマグニチュードの広い範囲で、計測震度を十分な精度で近似できる。

図９に示すように、本発明の実施形態と特許文献１に開示された技術を用いたＩＩＲフィルターの特性の差はわずかに見えるが、図１１乃至図１４に対応する震度誤差（計測震度相当値と気象庁計測震度の差）を特許文献１に開示された技術を用いて求めてグラフ化すると、図１５乃至図１８となる。すなわち、図１５乃至図１８における震度誤差は特許文献１に開示された技術を用いて得られた計測震度相当値と気象庁計測震度の差を示している。
図１５では誤差分布のマイナス側に裾を引いている。これは、特許文献１に開示されたＩＩＲフィルタが０．６Ｈｚの周辺で気象庁の定義するＦＦＴフィルターの特性を正しく近似できず、利得を下回ったことが主たる原因と考えられる。
図１６の横軸は震央距離であり、１００ｋｍ以上の遠地地震で卓越周波数が０．６ＨｚになるとＩＩＲフィルターの近似誤差がマイナス側に偏った震度誤差として顕在化したことが示されている。
図１７の横軸は計測震度であり、図１５でマイナス側に偏った震度誤差分布は広く分散している。
図１８の横軸はマグニチュードであり、図１５でマイナス側に偏った震度誤差分布はＭ５以上にまとまっている。
前述のとおり、本発明の実施形態と特許文献１に開示された技術を用いたＩＩＲフィルター特性の差は図９に示すようにわずかであるように見えるが、図１１乃至図１４と図１５乃至図１８を対比すれば、本発明の実施の形態に係るリアルタイム震度計測装置１の震度測定精度が大幅に向上していることが確認できる。

以上説明したとおり、計測震度を求めるときに用いられる気象庁のＦＦＴフィルターを近似するＩＩＲフィルター７は、G_a(z)、G_b1(z)、G_b2(z)、G_b3(z)、G_c(z)、及び利得g_Tの積で与えられる。以下に、f_a1、f_a2、h_a、f_b、γ、η、g_T、f_c、h_cの９つの係数を定めるべく、最小二乗法を用いて行った解析について説明する。
G(z)=g_TG_a(z)G_b1(z)G_b2(z)G_b3(z)G_c(z) （４５）
ここで、
G_a(z)={β_a0+β_a1z^-1+β_a2z^-2}/{α_a0+α_a1z^-1+α_a2z^-2}， （４６）
β_a0=1+ζ_a1，β_a1=-2，β_a2=1-ζ_a1，α_a0=1+2h_aζ_a2+ζ_a2 ²，α_a1=2ζ_a2 ²-2，
α_a2=1-2h_aζ_a2+ζ_a2 ²，ζ_a1=πf_a1Δt，ζ_a2=πf_a2Δt

G_b1(z)=G_b12(z)/G_b11(z)={ β_b10+β_b11z^-1}/{α_b10+α_b11z^-1}， （４７）
G_b2(z)=G_b22(z)/G_b21(z)={ β_b20+β_b21z^-1}/{α_b20+α_b21z^-1}， （４８）
G_b3(z)=G_b32(z)/G_b31(z)={ β_b30+β_b31z^-1}/{α_b30+α_b31z^-1}， （４９）
β_b10=ζ_b+1，β_b11=ζ_b-1，α_b10=ζ_b+η，α_b11=ζ_b-η，β_b20=γζ_b+1，
β_b21=γζ_b-1，α_b20=ζ_b+γη，α_b21=ζ_b-γη，β_b30=ζ_b/γ+1，β_b31=ζ_b/γ-1，
α_b30=ζ_b+η/γ，α_b31=ζ_b-η/γ，ζ_b=πf_bΔt

G_c(z)={β_c0+β_c1z^-1+β_c2z^-2}/{α_c0+α_c1z^-1+α_c2z^-2}， （５０）
β_c0=ζ_c ²，β_c1=2ζ_c ²，β_c2=ζ_c ²，α_c0=1+2h_cζ_c+ζ_c ²，α_c1=2ζ_c ²-2，
α_c2=1-2h_cζ_c+ζ_c ²，ζ_c=πf_cΔt
である。式（４５）は、ＩＩＲフィルターのラプラス変換に対して、次式（５１）を用いてラプラス変数を近似して得たものである。
s^-1=(Δt/2)(1+z^-1)/(1-z^-1) （５１）
ラプラス変数とzとの関係は、本来次式（５２）で与えられる。
z=exp(sΔt) （５２）
変数sと周波数fとの関係は、次式（５３）で与えられる。
ｓ=2πjf （５３）

ＩＩＲフィルター７の周波数特性は、式（４６）、（４７）、（４８）、（４９）、及び（５０）に式（５２）および（５３）を代入して得られる。ＩＩＲフィルター７の振幅を、
A(f)=|G(f )|，A_a(f)=|G_a(f)|，A_b(f)=|G_b(f)|，A_c(f)=|G_c(f)|
として、以下のように表わされる。
A(f)=g_TA_a(f)A_b(f)A_c(f)， （５４）
A_a(f)=C_aP^1/2(ζ_a1)Q^-1/2(ζ_a2,h_a2)， （５５）
A_b(f)=C_b(η)P^1/2(ζ_b)P^-1/2(ζ_b/η)P^1/2(ζ_bγ)P^-1/2(ζ_b/ηγ)P^1/2(ζ_b/γ)
×P^-1/2(ζ_bγ/η)，（５６）
A_c(f)=C_pζ_c ²Q^-1/2(ζ_c,h_c) ， （５７）
ここで、
P(ζ)=C_pζ²+C_m， （５８）
Q(ζ,h)=C_p ²ζ⁴+2C_pC_mH(h)ζ²+C_m ²， （５９）
H(h)=2h²-1，f=2πfΔt，C_a=(1-cosφ)^1/2，C_b(η)=η^-3，
C_p=1+cosφ，C_m=1-cosφ， （６０）
である。

求めるべきＩＩＲフィルター７の係数はf_a1、f_a2、h_a、f_b、γ、η、g_T、f_c、h_cである。ただし、ここではγ=4、η=2として固定し、f_a1、f_a2、f_b、f_cに関しては、
ζ_a1=πf_a1Δt、ζ_a2=πf_a2Δt、ζ_b=πf_bΔt、ζ_c=πf_cΔtによって変換されたζ_a1、ζ_a2、ζ_b、ζ_cを求める係数とする。式（５４）〜（５９）で与えられるＩＩＲフィルター７の特性はこれらの変数に関して非線形である。フィルター係数を求めるにあたって、先ず係数のそれぞれについて初期値を設定し、初期値からの微小な量の変動を考え、微小量の３次以上の項を無視して、フィルター特性を与える式を初期値の周りで微小量のべきで展開することにする。
x_i=x_i0+δx_i，
A(f)=A₀(f)+Σ_i∂A(f)/∂x_iδx_i+1/2Σ_iΣ_j∂²A(f)/∂x_i∂x_jδx_iδx_j （６１）
ここで、x_iは、ζ_a1、ζ_a2、h_a、ζ_b、g_T、ζ_c、h_cのうちいずれかであり、x_i0はx_iの初期値、dx_iは初期値からの微小な変動量、また、A₀(f)は係数が初期値のときのA(f )の値である。気象庁の定義するＦＦＴフィルター特性を広い周波数帯域で最もよく近似するＩＩＲフィルター７の最適解は、以下の量のいずれかを最小にするフィルター係数の組み合わせで求められる。
ρ_N=Σ_fw(f){|J(f)|-A(f)}²， （６２）
ρ_G=Σ_fw(f){log|J(f)|-logA(f)}²， （６３）
ここで、Σ_fは近似の程度を評価する周波数帯域での総和、w(f)は総和をとるときの重みである。また、ρ_N及びρ_Gは重み付き残差平方和である。最適解は次式を満たす。
∂ρ_N/∂x_m=Σ_f(-2)w(f){|J(f)|-A(f)}∂A/∂x_m=0， （６４）
∂ρ_G/∂x_m=Σ_f(-2)w(f){log|J(f)|-logA(f)}A^-1(f)∂A/∂x_m=0， （６５）
式（６４）では、フィルターのスペクトル振幅の差、式（６５）では比を、評価する周波数帯域内で平均として小さくする係数が解として得られることになる。式（６４）あるいは（６５）に式（６１）を代入し、∂x_iの２次以上の項を無視すると以下のように整理される。
Σ_iΣ_fw(f){∂A/∂x_mΣ_i∂A(f)/∂x_i-R_N(f)∂²A(f)/∂x_m∂x_i}δx_i
=Σ_fw(f)R_N(f)∂A/∂x_m，
R_N(f)=|J(f)|-A₀(f)， （６６）
Σ_iΣ_fw(f)[A^-2(f){1+R_G(f)}∂A/∂x_mΣ_i∂A(f)/∂x_i-R_G(f)A^-1(f)∂²A(f)/∂x_m∂x_i]δx_i
=Σ_fw(f)R_G(f)A^-1(f)∂A/∂x_m，
R_G(f )=log|J(f)|-logA₀(f) ， （６７）
式（６６）、及び式（６７）は、以下のように行列を用いて表すことができる。
Σ_iK_miδx_i＝D_m， （６８）
ただし、式（６６）を用いたときは、
K_mi=Σ_fw(f){∂A/∂x_m∂A(f)/∂x_i-R_N(f)∂²A(f)/∂x_m∂x_i}，
D_m=Σ_fw(f)R_N(f)∂A/∂x_m， （６９）
であり、式（６７）を用いたときは、
K_mi=Σ_f w(f)[A^-2(f){1+R_G(f)}∂A/∂x_m∂A(f)/∂x_i-R_G(f)A^-1(f)∂²A(f)/∂x_m∂x_i]，
D_m=Σ_fw(f)R_G(f)A^-1(f)∂A/∂x_m， （７０）
である。実際に解を求めるときは、安定性を考慮して式（６８）の代わりに以下の式を用いる。
Σ_I(K_mi+kδ_mi)δx_i＝D_m， （７１）
ここで、kは安定して解を得るための定数である。ひとつの初期値x_i0に対して式（７１）を用いて得られたδx_iから、次の初期値x_i0+δx_iが設定される。この操作を繰り返して、すべてのδx_iが十分に小さくなったときのx_iをフィルター係数の解とする。

なお、P (ζ)、P(aζ)、Q (ζ,h)の一階及び二階の偏微分係数は以下のとおり表される。
∂P(ζ)/∂ζ=2C_pζ， （７２）
∂P(aζ)/∂ζ=2C_pa²ζ， （７３）
∂²P(ζ)/∂ζ_a1 ²=2C_p， （７４）
∂²P(aζ)/∂ζ_a1 ²=2C_pa²， （７５）
∂Q(ζ,h)/∂ζ=4{C_p ²ζ³+C_pC_mH(h)ζ}， （７６）
∂Q(ζ,h)/∂h=8C_pC_mζ²h， （７７）
∂²Q(ζ,h)/∂ζ²=4{3C_p ²ζ²+C_pC_mH(h)}， （７８）
∂²Q(ζ,h)/∂h²=8C_pC_mζ²， （７９）
∂²Q(ζ,h)/∂ζ∂h=16C_pC_mζh， （８０）
これらを用いて、A (f)のζ_a1、ζ_a2、h_a2、ζ_b、ζ_c、h_cについての一階及び二階の偏微分係数は以下のとおり表される。
∂A(f)/∂ζ_a1=(1/2)g_TA_b(f)A_c(f)C_aQ ^-1/2(ζ_a2,h_a2)P^-1/2(ζ_a1)∂P(ζ_a1)/∂ζ_a1，
∂A(f)/∂ζ_a2=(-1/2)g_TA_b(f)A_c(f)C_aP ^1/2(ζ_a1)Q^-3/2(ζ_a2,h_a2)∂Q(ζ_a2,h_a2)/∂ζ_a2，
∂A(f)/∂h_a2=(-1/2)g_TA_b(f)A_c(f)C_aP ^1/2(ζ_a1)Q ^-3/2(ζ_a2,h_a2)∂Q(ζ_a2,h_a2)/∂h_a2，
∂A(f)/∂ζ_b =g_TA_a(f)A_c(f)A_b(f)C_pC_mζ_bT(ζ_b,γ,η)，
∂A(f)/∂ζ_c=g_TA_a(f)A_b(f)C_p{2ζ_cQ^-1/2(ζ_c,h_c)-(1/2)ζ_c ²Q^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂ζ_c}，
∂A(f)/∂h_c=(-1/2)g_TA_a(f)A_b(f)C_pζ_c ²Q^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂h_c，
∂A(f)/∂g_T=A_a(f)A_b(f)A_c(f)，
∂²A(f)/∂ζ_a1 ² =(1/2)g_TA_b(f)A_c(f)C_aQ^-1/2(ζ_a2,h_a2)
P^-3/2(ζ_a1){(-1/2)(∂P(ζ_a1)/∂ζ_a1)²+P(ζ_a1)∂²P(ζ_a1)/∂ζ_a1 ² }，
∂²A(f)/∂ζ_a2 ²=(-1/2)g_TA_b(f)A_c(f)C_aP ^1/2(ζ_a1)
Q^-5/2(ζ_a2,h_a2){(-3/2)(∂Q(ζ_a2,h_a2)/∂ζ_a2)²+Q(ζ_a2,h_a2)∂²Q(ζ_a2,h_a2)/∂ζ_a2 ²}，
∂²A(f)/∂h_a2 ²=(-1/2)g_TA_b(f)A_c(f)C_aP^1/2(ζ_a1)
Q^-5/2(ζ_a2,h_a2){(-3/2)(∂Q(ζ_a2,h_a2)/∂h_a2)²+Q(ζ_a2,h_a2)∂²Q(ζ_a2,h_a2)/∂h_a2 ²}，
∂²A(f)/∂ζ_b ²=g_TA_a(f)A_c(f)A_b(f){C_p ²C_m ²ζ_b ²T ²(ζ_b,γ,η)
+C_pC_mT (ζ_b,γ,η)+C_pC_mζ_b∂T(ζ_b,γ,η)/∂ζ_b，
∂²A(f)/∂ζ_c ²=g_TA_a(f)A_b(f)C_p[2Q^-1/2(ζ_c,h_c)-2ζ_cQ^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂ζ_c
-(1/2) ζ_c ²{Q^-3/2(ζ_c,h_c)∂²Q(ζ_c,h_c)/∂ζ_c ²-(3/2)Q^-5/2(ζ_c,h_c)(∂Q(ζ_c,h_c)/∂ζ_c)²}]，
∂²A(f)/∂h_c ²=(-1/2)g_TA_a(f)A_b(f)C_pζ_c ²
{Q^-3/2(ζ_c,h_c)∂²Q(ζ_c,h_c)/∂h_c ²-(3/2)Q^-5/2(ζ_c,h_c)(∂Q(ζ_c,h_c)/∂h_c)²}，
∂²A(f)/∂ζ_a1∂ζ_a2=(-1/4)g_TA_b(f)A_c(f)C_a
{P^-1/2(ζ_a1)∂P(ζ_a1)/∂ζ_a1Q^-3/2(ζ_a2,h_a2)∂Q(ζ_a2,h_a2)/∂ζ_a2}，
∂²A(f)/∂ζ_a1∂h_a2=(-1/4)g_TA_b(f)A_c(f)C_a
{P^-1/2(ζ_a1)∂P(ζ_a1)/∂ζ_a1Q^-3/2(ζ_a2,h_a2)∂Q(ζ_a2,h_a2)/∂h_a2}，
∂²A(f)/∂ζ_a1∂ζ_b=(1/2)g_TA_c(f)C_aQ^-1/2(ζ_a2,h_a2)P^-1/2(ζ_a1)∂P(ζ_a1)/∂ζ_a1
×A_b(f)C_pC_mz_bT(ζ_b,γ,η)，
∂²A(f)/∂ζ_a1∂ζ_c=(1/2)g_TA_b(f)C_aQ^-1/2(ζ_a2,h_a2)P^-1/2(ζ_a1)∂P(ζ_a1)/∂ζ_a1
×C_p{2ζ_cQ^-1/2(ζ_c,h_c)-(1/2)ζ_c ²Q^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂ζ}，
∂²A(f)/∂ζ_a1∂h_c=(1/2)g_TA_b(f)C_aQ^-1/2(ζ_a2,h_a2)P^-1/2(ζ_a1)∂P(ζ_a1)/∂ζ_a1
×{(-1/2)C_pζ_c ²Q^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂h_c}，
∂²A(f)/∂ζ_a2∂h_a2=(-1/2)g_TA_b(f)A_c(f)C_aP^1/2(ζ_a1)
×{Q^-3/2(ζ_a2,h_a2)∂²Q(ζ_a2,h_a2)/∂ζ_a2∂h_a2
- (3/2)Q^-5/2(ζ_a2,h_a2)∂Q(ζ_a2,h_a2)/∂ζ_a2∂Q(ζ_a2,h_a2)/∂h_a2}，
∂²A(f)/∂ζ_a2∂ζ_b=(-1/2)g_TA_b(f)A_c(f)C_aP^1/2(ζ_a1)
×Q^-3/2(ζ_a2,h_a2)∂Q(ζ_a2,h_a2)/∂ζ_a2C_pC_mζ_bT(ζ_b,γ,η)，
∂²A(f)/∂ζ_a2∂ζ_c=(-1/2)g_TA_b(f)C_aP^1/2(ζ_a1)Q^-3/2(ζ_a2,h_a2)∂Q(ζ_a2,h_a2)/∂ζ_a2
×C_p{2ζ_cQ ^-1/2(ζ_c,h_c)-(1/2) ζ_c ²Q^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂ζ_c}，
∂²A(f)/∂ζ_a2∂h_c=(-1/2)g_TA_b(f)C_aP^1/2(ζ_a1)Q^-3/2(ζ_a2,h_a2)∂Q(ζ_a2,h_a2)/∂ζ_a2
×C_pζ_c ²Q^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂h_c，
∂²A(f)/∂h_a2∂ζ_b=(-1/2)g_TA_b(f)A_c(f)C_aP^1/2(ζ_a1)
×Q^-3/2(ζ_a2,h_a2)∂Q(ζ_a2,h_a2)/∂h_a2C_pC_mζ_bT(ζ_b,γ, η)，
∂²A(f)/∂h_a2∂ζ_c=(-1/2)g_TA_b(f)C_aP ^1/2(ζ_a1)Q^-3/2(ζ_a2,h_a2)∂Q(ζ_a2,h_a2)/∂h_a2
×C_p{2ζ_cQ^-1/2(ζ_c,h_c)-(1/2) ζ_c ²Q^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂ζ_c}，
∂²A(f)/∂h_a2∂h_c=(-1/2)g_TA_b(f)C_aP^1/2(ζ_a1)Q^-3/2(ζ_a2,h_a2)∂Q(ζ_a2,h_a2)/∂h_a2
×C_pζ_c ²Q^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂h_c，
∂²A(f)/∂ζ_b∂ζ_c=g_TA_a(f)A_b(f)C_pC_mζ_bT(ζ_b,γ, η)
×C_p{2ζ_cQ^-1/2(ζ_c,h_c)-(1/2)ζ_c ²Q^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂ζ_c}，
∂²A(f)/∂ζ_b∂h_c=g_TA_a(f)A_b(f)C_pC_mζ_bT(ζ_b,γ,η)
×(-1/2)C_pζ_c ²Q^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂h_c，
∂²A(f)/∂ζ_c∂h_c=g_TA_a(f)A_b(f)C_p[-ζ_cQ^-3/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂h_c
-(1/2)ζ_c ²Q^-3/2(ζ_c,h_c)∂²Q(ζ_c,h_c)/∂ζ_c∂h_c
+(3/4)ζ_c ²Q^-5/2(ζ_c,h_c)∂Q(ζ_c,h_c)/∂ζ_c∂Q(ζ_c,h_c)/∂h_c]
ここで、
T(ζ_b,γ,η)=(1-η^-2)P^-1(ζ_b)P^-1(ζ_b/η)+(γ²-γ^-2η^-2)P ^-1(γζ_b)P ^-1(ζ_b/γη)
+(γ^-2-γ²η^-2)P^-1(ζ_b/γ)P^-1(γζ_b/η)，
∂T(ζ_b,γ,η)/∂ζ_b=-(1-η^-2)P ^-2(ζ_b)P^-2(ζ_b/η)
{P(ζ_b/η)∂P(ζ_b)/∂ζ_b+P(ζ_b)∂P(ζ_b/η)/∂ζ_b}
-(γ²-γ^-2η^-2)P^-2(γζ_b)P^-2(ζ_b/γη){P(ζ_b/γη)∂P(γζ_b)/∂ζ_b+P(γζ_b)∂P(ζ_b/γη)/∂ζ_b}
-(γ^-2-γη^-2)P ^-2(ζ_b/γ)P ^-2(γζ_b/η){P(γζ_b/η)∂P(ζ_b/γ)/∂ζ_b+P(ζ_b/γ)∂P(γζ_b/η)/∂ζ_b}，
である。

以上、説明した本実施の形態に係るリアルタイム震度計測装置１は、０．３＜ｆ≦１Ｈｚの周波数帯域において、極大値を取る気象庁定義のＦＦＴフィルターの特性を、遮断周波数と減衰定数に自由度を有する２次ＨＰＦと遮断周波数に自由度を有する１次ＨＰＦを比として組み合わせた第１のＩＩＲフィルター７１を用いることで十分高精度な近似を得ることが可能となった。しかも、ＩＩＲフィルター７による出力を震度演算部８によって処理することで、リアルタイムの震度計算が可能であるので、地震直後から出力部４によって震度情報を提供することが可能であり、時間遅れの小さな防災活動が可能である。頻発する地震に対する安全対策を講じる上で顕著な効果を発揮し得るのである。

本発明の請求項１乃至請求項４に記載された発明は、地震防災措置を講ずる必要がある施設や地域において有効であると同時に、地震速報をはじめとして緊急地震速報を報知する機関や施設においても利用可能である。

１…リアルタイム震度計測装置 ２…演算部 ３…計時部 ３Ｓ１…時間信号 ３Ｓ２…時刻信号 ４…出力部 ５…加速度計 ５Ｓ…加速度波形信号 ６…Ａ／Ｄ変換器 ６Ｓ…加速度信号 ７…ＩＩＲフィルター ７Ｓ…重み付け加速度信号 ８…震度演算部 ８Ｓ…震度信号 ７１…第１のＩＩＲフィルター ７２…第２のＩＩＲフィルター ７３…第３のＩＩＲフィルター

Claims (1)

1. 地震による加速度を検出する加速度検出手段と、
   この加速度検出手段から出力されるアナログ信号を第１のデジタル信号に変換するＡ／Ｄ変換器と、
   このＡ／Ｄ変換器から出力される前記第１のデジタル信号をフィルター処理して第２のデジタル信号を出力するＩＩＲ（Infinite Impulse Response:無限インパルス応答）フィルターと、
   このＩＩＲフィルターから出力される第２のデジタル信号に基づいてリアルタイム震度を計算する震度演算部と、
   このリアルタイム震度を出力する出力部と、を有するリアルタイム震度計測装置において、
   前記ＩＩＲフィルターは、ゲインピークを含む周波数帯域において、遮断周波数と減衰定数の２つの自由度を有する２次ハイパスフィルター（ＨＰＦ）と遮断周波数の１つの自由度を有する１次ハイパスフィルター（ＨＰＦ）を組み合わせて構成される第１のＩＩＲフィルターを備え、
   前記ゲインピークを含む周波数（以下、周波数をｆで表わす）帯域は０．３Ｈｚ＜ｆ≦１Ｈｚの範囲を含む帯域であって、
   前記ＩＩＲフィルターは、前記ゲインピークを含む周波数帯域以下の低周波数帯域においては、前記第１のＩＩＲフィルターを用い、前記ゲインピークを含む周波数帯域よりも高い１Ｈｚ＜ｆ≦３Ｈｚを含む中周波数帯域においては、遮断周波数の１つの自由度を有する３つの１次フィルターを組み合わせて周波数の平方根の逆数に比例する第２のＩＩＲフィルターを備え、前記中周波数帯域よりも高い３Ｈｚ＜ｆを含む高周波数帯域においては、遮断周波数と減衰定数の２つの自由度を有する２次ローパスフィルター（ＬＰＦ）を第３のＩＩＲフィルターとして備えていることを特徴とするリアルタイム震度計測装置。

Images