(57)【特許請求の範囲】
【請求項１】
地震波のＰ波を抽出し、前記Ｐ波に基づいてＰ波震度Ｖｐを計算すると共に、前記地震波を複数のバンドパスフィルタでフィルタ処理してそれぞれの最大加速度Ａｉを検出し、前記Ｐ波震度Ｖｐの過去分データ及び推定されたＳ波震度Ｖｓの過去分データに基づいて、多変数の最小二乗法を用いて分散が最小になる定数αｉ及びβｉを求め、観測点に依存する定数Ｃｊ、前記定数αｉ及びβｉ、前記Ｐ波震度Ｖｐ、前記最大加速度Ａｉから、ｆｉ（）を関数として、Ｖｓ＝Ｃｊ＋Σαｉ・ｆｉ（Ａｉ）＋βｉ・Ｖｐに基づいて前記Ｓ波震度Ｖｓを推定することを特徴とする震度推定方法。
【請求項２】
地震波のＰ波を抽出し、前記Ｐ波に基づいてＰ波震度Ｖｐを計算すると共に、前記地震波を低周波バンドパスフィルタ及び高周波バンドパスフィルタでフィルタ処理してそれぞれの最大加速度Ａ１及びＡ２を検出し、前記Ｐ波震度Ｖｐの過去分データ及び推定されたＳ波震度Ｖｓの過去分データに基づいて、多変数の最小二乗法を用いて分散が最小になる定数α１、α２及びβを求め、観測点に依存する定数Ｃｊ、前記定数α１、α２及びβ、前記Ｐ波震度Ｖｐ、前記最大加速度Ａ１及びＡ２から、ｆ１（）及びｆ２（）を関数として、Ｖｓ＝Ｃｊ＋α１・ｆ１（Ａ１）＋α２・ｆ２（Ａ２）＋β・Ｖｐに基づいて前記Ｓ波震度Ｖｓを推定することを特徴とする震度推定方法。
【請求項３】
地震波のＰ波を抽出するＰ波抽出部と、前記Ｐ波に基づいてＰ波震度Ｖｐを計算するＰ波震度計算部と、前記地震波を複数のバンドパスフィルタでフィルタ処理するバンドパスフィルタと、前記バンドパスフィルタからの出力の最大加速度Ａｉを検出する最大加速度検出部と、前記Ｐ波震度Ｖｐの過去分データ及び推定されたＳ波震度Ｖｓの過去分データに基づいて、多変数の最小二乗法を用いて分散が最小になる定数αｉ及びβｉを求める定数演算部と、観測点に依存する定数Ｃｊ、前記定数αｉ及びβｉ、前記Ｐ波震度Ｖｐ、前記最大加速度Ａｉから、ｆｉ（）を関数として、Ｖｓ＝Ｃｊ＋Σαｉ・ｆｉ（Ａｉ）＋βｉ・Ｖｐに基づいて前記Ｓ波震度Ｖｓを推定するＳ波震度推定部とを具備したことを特徴とする震度推定装置。
【請求項４】
地震波のＰ波を抽出するＰ波抽出部と、前記Ｐ波に基づいてＰ波震度Ｖｐを計算するＰ波震度計算部と、前記地震波を低周波バンドパスフィルタ処理する低周波バンドパスフィルタと、前記地震波を高周波バンドパスフィルタ処理する高周波バンドパスフィルタと、前記低周波バンドパスフィルタの出力信号の最大加速度Ａ１を検出する第１最大加速度検出部と、前記高周波バンドパスフィルタの出力信号の最大加速度Ａ２を検出する第２最大加速度検出部と、前記最大加速度Ａ１及びＡ２の比率Ａ１／Ａ２を求める比率計算部と、前記Ｐ波震度Ｖｐの過去分データ及び推定されたＳ波震度Ｖｓの過去分データに基づいて、多変数の最小二乗法を用いて分散が最小になる定数α及びβを求める定数演算部と、観測点に依存する定数Ｃｊ、前記定数α及びβ、前記Ｐ波震度Ｖｐ、前記比率Ａ１／Ａ２から、ｆ（）を関数として、Ｖｓ＝Ｃｊ＋α・ｆ（Ａ１／Ａ２）＋β・Ｖｐに基づいてＳ波震度Ｖｓを推定するＳ波震度推定部とを具備したことを特徴とする震度推定装置。
【請求項５】
地震波のＰ波を抽出するＰ波抽出部と、前記Ｐ波に基づいてＰ波震度Ｖｐを計算するＰ波震度計算部と、前記地震波を低周波バンドパスフィルタ処理する低周波バンドパスフィルタと、前記地震波を高周波バンドパスフィルタ処理する高周波バンドパスフィルタと、前記低周波バンドパスフィルタの出力信号の最大加速度Ａ１を検出する第１最大加速度検出部と、前記高周波バンドパスフィルタの出力信号の最大加速度Ａ２を検出する第２最大加速度検出部と、前記Ｐ波震度Ｖｐの過去分データ及び推定されたＳ波震度Ｖｓの過去分データに基づいて、多変数の最小二乗法を用いて分散が最小になる定数α１、α２及びβを求める定数演算部と、観測点に依存する定数Ｃｊ、前記定数α１、α２及びβ、前記Ｐ波震度Ｖｐ、前記最大加速度Ａ１及びＡ２から、ｆ１（）及びｆ２（）を関数として、Ｖｓ＝Ｃｊ＋α１・ｆ１（Ａ１）＋α２・ｆ２（Ａ２）＋β・Ｖｐに基づいて前記Ｓ波震度Ｖｓを推定するＳ波震度推定部とを具備したことを特徴とする震度推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は単独観測点による地震の震度推定方法及び装置に関し、特に地震発生から最初に到達するＰ波（初期微動）成分を用いて、Ｓ波（主要動）が到達する前に観測点における震度を精度良く推定するための方法及び装置に関する。
【背景技術】
【０００２】
従来、単独観測点方式で、初動部分であるＰ波の情報からＳ波（主要動）の震度を推定する方法では、Ｐ波の加速度の値が閾値を超えた場合に、その後に大きな主要動が来ると想定してトリガ信号を出すレベルトリガ法が主に採用されている。このレベルトリガ法は主要動が到達する前にトリガ信号を出すため、大きく揺れる前に設備の終了動作を始めるなど、地震被害を軽減させることができる。
【０００３】
しかしながら、発震機構の影響により、Ｐ波の振幅とＳ波の振幅の相関のばらつきが大きいため、Ｐ波のレベルでトリガ信号を出したにも拘わらず大きなＳ波主要動が来なかったり、反対にＰ波でトリガ信号が出なかったが大きなＳ波主要動が来たりする誤報が非常に多かった。また、単独観測点で地震波を観測するときに、その初動であるＰ波の情報を基に到達地震波の震度を推定する別の方法が気象庁の緊急地震速報で採用されており、以下の方法が提案されている。
【０００４】
特開２００２−２７７５５７号公報（特許文献１）では、先ず単独観測点において観測される地震波のＰ波の立ち上がりからの数秒間の情報を基に、震央距離とマグニチュードを求め、それらの情報から観測地点での震度を推定している。即ち、特許文献１では、地震計から得られる地震波初動部分の波形形状の特徴に注目し、その波形形状をパラメータが数個の簡易な関数でフィッティングすることで定量化し、得られたパラメータから震央距離及びマグニチュードを推定している。
【０００５】
しかしながら、特許文献１に開示された方法では計算された震度の誤差が大きいため、その情報を基に設備機器などを制御することは事実上難しい。従って、気象庁の地震観測ネットワークにおいては基本的に、複数の観測点で地震波を観測してから判定する手法が採用されている。
【０００６】
一方、自前で地震計を設置して地震観測を行う場合には、広域に多点の地震観測を展開することが難しいため、単独観測点での情報に基づいて、震度をできるだけ正確に推定することが求められる。現況は、上述のように精度の良くない震度推定を行うか、特開２００８−１０７３３４号公報（特許文献２）に開示されているように、計測震度相当値をリアルタイムで計算する方法を採らざるを得ない。即ち、特許文献２の装置は、地動加速度の時系列を得る地動加速度時系列取得手段と、地動加速度の時系列をフィルタ処理する時間領域フィルタ手段と、時間領域フィルタ手段によりフィルタ処理された地動加速度の時系列から計測震度の概算値を算出する算出手段とを備えている。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２００２−２７７５５７号公報
【特許文献２】特開２００８−１０７３３４号公報
【非特許文献】
【０００８】
【非特許文献１】日本物理探査学会誌「緊急地震速報における震度マグニチュードの有効性」、第６０巻第５号、ｐ４０７〜４１７
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、特許文献２の装置によるリアルタイム計測は、実際に観測点で揺れた事実を基に迅速に判定するものであり、Ｐ波の初期微動を基に震度を推定するものではないため、Ｓ波主要動が到達する前の制御は事実上不可能である。
【００１０】
実際に断層面がずれて地震波が発生した場合、観測点と断層面の位置関係により、到達する地震波のＰ波の振幅とＳ波の振幅の比が異なってしまう。これを発震機構といい、理論的にはＰ波部分の振幅が大きい方向でＳ波部分の振幅が小さく、Ｐ波部分の振幅が小さい方向でＳ波部分の振幅が大きい。そのため、Ｐ波の振幅情報のみを使用してＳ波の最大振幅や震度を推定すると必然的に誤差を含んでしまう。発震機構は、断層により地震が発生したときに断層形状により理論的に求めることができる。そのため、断層を囲むように多点の観測点があれば理論的に求めることが可能であるが、単独観測点における情報では求めることは不可能である。
【００１１】
発震機構以外にも、地震波の減衰による影響もある。地震波は地中を伝わるときに高周波の減衰が特に大きいため、この減衰特性を補正する必要がある。そのため、地震の震度計測においては、全体の簡易性から単独観測点方式が望ましく、単独の地震観測点で観測若しくは計測される震度を、地震波のＰ波部分の情報を基に、発震機構及び減衰による補正を行って、地震波の主要動であるＳ波の最大振幅や震度をより正確に計測することが要請されている。
【００１２】
本発明は上述のような事情よりなされたものであり、本発明の目的は、地震波を発生させた発震機構自体が分からなくても、地震波が地中を伝播する際に発生する散乱の影響を考慮することによって、Ｐ波の発震機構の影響を想定して補正すると共に、地震波減衰の影響も補正し、Ｐ波の震度を高速度に計算することにより、Ｓ波が実際に到達する前に震度を精度良く推定する方法及び装置を提供することにある。
【課題を解決するための手段】
【００１３】
本発明は震度推定方法に関し、本発明の上記目的は、地震波のＰ波を抽出し、前記Ｐ波に基づいてＰ波震度Ｖｐを計算すると共に、前記地震波を複数のバンドパスフィルタでフィルタ処理してそれぞれの最大加速度Ａｉを検出し、前記Ｐ波震度Ｖｐの過去分データ及び推定されたＳ波震度Ｖｓの過去分データに基づいて、多変数の最小二乗法を用いて分散が最小になる定数αｉ及びβｉを求め、観測点に依存する定数Ｃｊ、前記定数αｉ及びβｉ、前記Ｐ波震度Ｖｐ、前記最大加速度Ａｉから、ｆｉ（）を関数として、Ｖｓ＝Ｃｊ＋Σαｉ・ｆｉ（Ａｉ）＋βｉ・Ｖｐに基づいて前記Ｓ波震度Ｖｓを推定することにより達成される。
また、本発明は震度推定装置に関し、本発明の上記目的は、地震波のＰ波を抽出するＰ波抽出部と、前記Ｐ波に基づいてＰ波震度Ｖｐを計算するＰ波震度計算部と、前記地震波を複数のバンドパスフィルタでフィルタ処理するバンドパスフィルタと、前記バンドパスフィルタからの出力の最大加速度Ａｉを検出する最大加速度検出部と、前記Ｐ波震度Ｖｐの過去分データ及び推定されたＳ波震度Ｖｓの過去分データに基づいて、多変数の最小二乗法を用いて分散が最小になる定数αｉ及びβｉを求める定数演算部と、観測点に依存する定数Ｃｊ、前記定数αｉ及びβｉ、前記Ｐ波震度Ｖｐ、前記最大加速度Ａｉから、ｆｉ（）を関数として、Ｖｓ＝Ｃｊ＋Σαｉ・ｆｉ（Ａｉ）＋βｉ・Ｖｐに基づいて前記Ｓ波震度Ｖｓを推定するＳ波震度推定部とを設けることにより達成される。
【発明の効果】
【００１４】
従来は、単独観測点での地震波を観測してＳ波の震度を推定する方法としては、Ｐ波の情報を基にした精度の悪いＳ波震度を推定するか、実際に現地で大きい揺れを観測してから、その時の計測震度相当値をいち早く出力するしかなかったが、本発明によれば、単独観測点におけるＰ波（初期微動）の情報を基に、Ｓ波（主要動）が実際に到達する前に、従来より精度良く震度を推定することができる。
【図面の簡単な説明】
【００１５】
【図１】発震機構によるＰ波及びＳ波の振幅例を示す分布図である。
【図２】地震波の一例を示す波形図である。
【図３】本発明に係る震度推定装置の一例を示すブロック構成図である。
【図４】Ｐ波震度計算部の構成例を示すブロック図である。
【図５】地震信号の水平動及び上下動の加速度信号の一例を示す波形図である。
【図６】定数演算部の構成例を示すブロック図である。
【図７】本発明の全体の動作例を示すフローチャートである。
【図８】Ｐ波震度計算の処理例を示すフローチャートである。
【図９】定数演算の処理例を示すフローチャートである。
【発明を実施するための形態】
【００１６】
地震の規模が同じでも発震機構（ラディエーションパターン）の影響で、Ｐ波（初期微動）及びＳ波（主要動）の各振幅は方位により変わる。理論的な振幅は、Ｐ波の振幅が大きい方向でＳ波の振幅が小さく、逆にＰ波の振幅が小さい方向でＳ波の振幅が大きくなる。発震機構による振幅の方位変化を補正項で表すと、Ｐ波に対応する補正項が大きいところでＳ波に対応する補正項は小さく、逆にＰ波に対応する補正項が小さいところでＳ波に対応する補正項が大きくなる。即ち、地震波は地中を伝播するとき、高周波成分は散乱の影響を受け、振幅の方位依存性が観測されないことが知られており、低周波成分の振幅は地中伝播時の散乱の影響を受け難いため、発震機構による振幅の方位依存性が見られる。方位依存性が見られるデータと、方位依存性が見られないデータとの振幅の違いを測定することにより、Ｐ波の発震機構による補正項を求めることができる。
【００１７】
図１は発震機構によるＰ波及びＳ波の各振幅の分布例を示しており、図１（Ａ）がＰ波の分布例を、図１（Ｂ）がＳ波の分布例をそれぞれ示している。クローバ形状の羽が大きい方向に振幅が大きいことを示している。図１が示すように、Ｐ波とＳ波とでは振幅が大きくなる方向が異なり、Ｐ波の振幅が大きい方向でＳ波の振幅が小さく、逆にＳ波の振幅が大きい方向でＰ波の振幅が小さい。また、地震波は周波数が高いほど地中伝播で散乱の影響を受け、Ｓ波よりＰ波の方が周波数が高いために、より散乱の影響を受ける。Ｐ波のうち、数１０Ｈｚ程度の成分（高周波）は散乱の影響で、発震機構の分布が見られずほぼ一様になるのに対し、数Ｈｚ程度の成分（低周波）は散乱の影響を余り受けずに発震機構の影響を受ける。従って、Ｐ波成分の数１０Ｈｚ成分（高周波成分）と数Ｈｚ成分（低周波成分）の比率若しくは相関関係をとることにより、その方位からの発震機構による影響度を推定することができる。推定された影響度の値を使って、Ｓ波の振幅を補正してＳ波震度を推定することが可能となる。
【００１８】
このようにＳ波振幅の発震機構による補正項は、Ｐ波の補正項が大きい方向で小さくなっていることから、Ｐ波の発震機構による補正項が求まれば、Ｓ波の補正項を求めることが可能である。
【００１９】
本発明は高周波成分の振幅と低周波成分の振幅の関係を基礎データとして、Ｐ波とＳ波の各補正項を求め、Ｐ波の振幅に基づいて正確なＳ波の振幅を推定する。その方法として、例えばＰ波において、伝播時の散乱の影響を余り受けない低周波成分と散乱の影響を受ける高周波成分の比をとることが考えられ、これにより観測されるＰ波の発震機構の影響を求めることができる。
【００２０】
ここで、Ｐ波に基づいて、到来するＳ波の震度を推定できるという本発明の原理を説明する。
【００２１】
独立行政法人防災科学技術研究所が運用する全国強震観測網（K-net）の観測点の震度別積算回数を調べてみると、表１のようになる。
【００２２】
【表１】
0005507903.tif 000002
表１から分かるように、震度が“１”異なると発生頻度は大体１０倍異なる。この情報を基に、緊急地震速報を制御（特定の震度が到来したときに、設備機器を止めるなどの操作をすること）に用いる場合の適切な制御の割合を想定した場合、震度がある閾値以上になると９５％の確率で制御できることと、震度の推定誤差が正規分布に従うこととを仮定すると、表２のようになる。「適切な制御」とは、制御を予定される震度に対して実際に制御を行うことであり、制御しなくても良い状態で制御を行ってしまうことを「不適切な制御」としている。
【００２３】
気象庁のサイトでは、予測される震度の±１階級以内に９４％が収まるとあり、仮にある施設で、「震度６弱でシステムの電源の終了操作など、何か制御を始める」ということを緊急地震速報で行わせる場合、緊急地震速報で震度6が来るという情報では、震度5強以上6強弱の範囲の地震が来ることが想定される。しかし、実際には、震度６弱の地震よりも震度５強の地震の方が約１０倍多いため、適切に震度６弱で制御をしたいのに、震度５強で動いてしまうケースが１０倍多くなってしまう。
【００２４】
【表２】
0005507903.tif 000003
表２において、例えば震度推定のＲＭＳ（root-mean-square）が±０．５である場合、９５％の確率で制御を可能にさせるならば、１０回のうち９回は小さい地震を大きい地震と判断してしまうということである。つまり、大きい地震を小さいと判断してしまう確率は５％である。緊急地震速報においては、震度推定のＲＭＳ±１の範囲で９４％が収まるようになっているが、この誤差が±１より更に小さくなれば、適切な制御が増えるというのが表２の意味である。
【００２５】
また、様々な指標を基に制御を行うとすると、表３のようになる。表３中の割合は、ＲＭＳが“０．４７の時は、８．７回に1回適切な制御になるということを示し、同様に気象庁マグニチュードのＲＭＳが”０．５６“では１３．２回に1回の割合で、P波センサのＲＭＳが”０．５４“では９．１回に1回の割合でそれぞれ適切な制御をすることを示している。
【００２６】
【表３】
0005507903.tif 000004
表３から分かるように気象庁震度マグニチュードの割合が一番良いが、それでもＲＭＳが“０．４７”である。本発明では、図１に示すようなＰ波とＳ波の発震機構の違いに着目し、図２に示すような地震波に対して、単独観測点で計測若しくは観測されるＰ波の震度から、地震波全体の震度（Ｓ波主要動）を推定する精度を上げる。発震機構の影響は地震波の低周波成分に顕著に見られ、高周波成分は地震波が伝播する際に散乱してしまい、顕著に見られない。従って、観測される高周波成分と低周波成分の比率若しくは相関を用いて補正することにより、Ｓ波震度推定の精度が向上させている。
【００２７】
以下に本発明の実施形態の一例を、図面を参照して説明する。
【００２８】
図３は本発明の実施形態の構成例を示しており、センサ及びＡ／Ｄ変換器等で成る震度計１０が単独観測点に設置され、図２に示すような地震波が震度計１０で検知されてディジタル化され、ディジタル化された地震信号Ｖｅが震度計１０から通信回線等で伝送されてＰ波抽出部２０に入力される。Ｐ波抽出部２０はノイズを除去してＰ波成分のみを抽出するが、これは、地震波のうち高周波成分は地中伝播時に直ぐに減衰すること、ノイズ源（数１００Ｈｚ）は近距離の高周波ノイズ成分（数１０Ｈｚ）と比較して周波数オーダーが１桁相違することから、数１００Ｈzサンプリングで地震波データを取得することで、単独観測点においても１秒程度で地震波（Ｐ波）であるかノイズであるかを容易に識別して抽出することができる。なお、地震信号Ｖｅの伝送をアナログで行い、Ｐ波抽出部２０に設けたＡ／Ｄ変換器でディジタル信号に変換してから、Ｐ波の抽出を行っても良い。
【００２９】
Ｐ波抽出部２０で抽出されたＰ波信号ＶｅｐはＰ波震度計算部３０に入力されると共に、低周波用バンドパスフィルタ（ＢＰＦ）２１及び高周波用バンドパスフィルタ（ＢＰＦ）２２に入力され、低周波用ＢＰＦ２１からのＰ波信号Ｖｅｐ１は最大加速度検出部２３に入力され、高周波用ＢＰＦ２２からのＰ波信号Ｖｅｐ２は最大加速度検出部２４に入力される。低周波用ＢＰＦ２１は通過帯周波数が数Ｈｚ（１〜５Ｈｚ）幅のバンドパスフィルタであり、高周波用ＢＰＦ２２は通過帯周波数が数１０Ｈｚ（１０〜５０Ｈｚ）幅のバンドパスフィルタである。
【００３０】
なお、バンドパスフィルタ２１及び２２は、通過域で最大平坦な振幅特性を有しているバターワースバンドパスフィルタが最適である。また、高周波、低周波のほかに中周波（５〜１０Ｈｚ）のバンドパスフィルタを設け、３つのパラメータとすることも可能である。
【００３１】
Ｐ波震度計算部３０は、Ｐ波信号Ｖｅｐに基づいてＰ波震度Ｖｐを計算する。気象庁のサイトで示されているように一旦周波数領域に変換してから計算しても良いが、フィルタで時間領域に戻す処理が必要なため時間がかかってしまう。このため、Ｓ波が到来する前に、高速度にＰ波震度を計算する必要がある本発明には適用できず、本発明では時間領域で逐次計測震度相当値を計算できる時間領域フィルタを使用することで、リアルタイムにＰ波震度Ｖｐを計算するようにしている。Ｐ波震度Ｖｐは、例えば非特許文献１に記載されているような手法で計算して求める。非特許文献１では、即時処理用に、速度から加速度への変換、センサの特性補正及び計測震度計算用のフィルタ処理を、時間領域で同時に行うことができる５０次の漸化フィルタを使用している。
【００３２】
地震情報などにより発表される震度階級は、観測点における揺れの強さの程度を数値化した計測震度から換算されており、表４のように階級化されている。
【００３３】
【表４】
0005507903.tif 000005
図４はＰ波震度計算部３０の構成例を示しており、Ｐ波信号Ｖｅｐは加速度計３１１（Ｘ方向（東西））、加速度計３１２（Ｙ方向（南北））及び加速度計３１３（Ｚ方向（上下））に入力され、図５（Ａ）〜（Ｃ）に示すような水平動＃１、水平動＃２及び上下動の３方向の加速度信号が得られ、各加速度信号はベクトル合成部３２で合成され、その後、地震波の周期による影響を補正するために、時間領域の例えば５０次の漸化フィルタ３３に入力され、フィルタ３３でフィルタ処理された合成加速度信号ＡＣ０は変数算出部３４に入力される。変数算出部３４は合成加速度信号ＡＣ０の波形の絶対値がある値a以上となる時間の合計をτ（ａ）とし、このτ（ａ）が下記数１を満たす値ａの最大値をａ０とする。
（数１）
τ（ａ）＞０．３秒
そして、変数算出部３４で算出された変数a０は震度計算部３５に入力され、下記数２に従ってＰ波震度Ｖｐが０．１単位で計算される。本例では、a０＝１２７．８５galである。
（数２）
Ｖｐ＝２log  a０  +0.94  
震度計算部３５で数２に従って計算されたＰ波震度Ｖｐは、Ｓ波震度推定部４０及び定数α、βを演算するための定数演算部５０に入力される。
【００３４】
一方、低周波用ＢＰＦ２１を通過したＰ波信号Ｖｅｐ１は最大加速度検出部２３に入力され、高周波用ＢＰＦ２２を通過したＰ波信号Ｖｅｐ２は最大加速度検出部２４に入力され、それぞれ所定時間内での最大加速度Ａ１及びＡ２が検出される。最大加速度Ａ１及びＡ２は比率計算部２５に入力され、比率Ａ１／Ａ２が計算される。計算された比率Ａ１／Ａ２はＳ波震度推定部４０に入力され、定数設定部４１に設定されている観測点に依存した定数Ｃｊ、定数演算部５０で演算された定数α、βもＳ波震度推定部４０に入力される。
【００３５】
Ｓ波震度推定部４０はＰ波震度Ｖｐ、比率Ａ１／Ａ２、定数Ｃｊ及びα、βを用いて、下記数３に従ってＳ波震度Ｖｓを推定する。
（数３）
Ｖｓ＝Ｃｊ＋α・ｆ（Ａ１／Ａ２）＋β・Ｖｐ
ただし、ｆ（）は関数であり、例えばlog関数である。
ここで、定数α、βを計算する定数演算部５０の構成を図６に示して説明する。
【００３６】
定数演算部５０は定数α、βを計算して地震発生の都度修正を行うものである。即ち、新しく地震が観測されたときに震度Ｖｐ，Ｖｓが観測されるので、その震度データは過去の地震データから作成できるＶｐ−Ｖｓ相関図上に1点追加し、その分布を基に最小二乗法で定数α、βを計算して修正を行う。
【００３７】
Ｐ波震度計算部３０で計算されたＰ波震度Ｖｐはメモリ５１に入力されると共に演算部５６に入力され、Ｓ波震度推定部４０で推定されたＳ波震度Ｖｓはメモリ５２に入力されると共に演算部５６に入力される。メモリ５１及び５２にはそれぞれ過去分のＰ波震度Ｖｐ及びＳ波震度Ｖｓが累積して記憶されており、メモリ５１及び５２におけるＰ波震度Ｖｐ及びＳ波震度Ｖｓの過去分データが演算部５５に入力され、演算部５５及び５６でそれぞれ前記数３を演算することにより震度Ｖｐ，Ｖｓの相関を求め、その相関の中から分散（標準偏差）が最も小さい回帰直線で定数α、βを求めて出力し、新しく求められた定数α、βをＳ波震度推定部４０に入力する。演算部５５，５６及び修正部５７では、多変数の最小二乗法を用いて分散が最小になるような定数α、βを計算若しくは修正する。具体的な計算の方法は、以下の通りである。
【００３８】
上記数３より、入力変数はＣｊ、α、βであり、観測から得られる値は震度Ｖｓ，関数ｆ（Ａ１／Ａ２），震度Ｖｐであるので、データは次のように表せる。
（数４）
Ｖｓ１＝Ｃｊ＋α・ｆ（Ａ１１／Ａ２１）＋β・Ｖｐ１
Ｖｓ２＝Ｃｊ＋α・ｆ（Ａ１２／Ａ２２）＋β・Ｖｐ２
Ｖｓ３＝Ｃｊ＋α・ｆ（Ａ１３／Ａ２３）＋β・Ｖｐ３
・・・
これを行列で表記すると、Ｖｓｎのベクトルを[Ｖｓ]とし、パラメータの転置ベクトルをＰ（Ｃｊ，α，β）としたとき、
（数５）
[Ｖｓ]＝Ｍ・Ｐ（Ｃｊ，α，β）
なる行列が成立する。ここで、最良の解を求めるためにＭの転置行列Ｍｔを考え、両辺の左から乗算することにより
（数６）
Ｍｔ・[Ｖｓ]＝（Ｍｔ・Ｍ）・Ｐ（Ｃｊ，α，β）
が得られ、この数６より下記数７となって、転置ベクトルＰ（Ｃｊ，α，β）が求まる。
（数７）
Ｐ（Ｃｊ，α，β）＝（Ｍｔ・Ｍ）−  １・Ｍｔ・[Ｖｓ]
なお、新たな地震データが観測された場合には、Ｖｓｘ＝Ｃｊ＋α・ｆ（Ａ１ｘ／Ａ２ｘ）＋β・Ｖｐｘの行が前記数４に１つ追加されることになり、上述と同様な計算で転置ベクトルＰ（Ｃｊ，α，β）が求まる。
【００３９】
このような構成において、その動作を図７のフローチャートを参照して説明する。
【００４０】
地震発生により図２に示すような地震波が発生し、この地震波が震度計１０で検知されるが（ステップＳ１）、震度計１０は周波数に基づいてノイズではない地震波のみを検知する。震度計１０からの地震信号Ｖｅは伝送されてＰ波抽出部２０に入力され、所定時間だけＰ波部分を抽出する（ステップＳ２）。地震波で最初に到達するのはＰ波であるので、Ｐ波の立ち上がりから所定時間（例えば３秒）だけＰ波を抽出し、抽出されたＰ波信号ＶｅｐをＰ波震度計算部３０及び低周波用ＢＰＦ２１、高周波用ＢＰＦ２２に入力するが、Ｐ波震度計算部３０は前述のような方法でＰ波震度Ｖｐの計算を行い（ステップＳ１０）、低周波用ＢＰＦ２１及び高周波用ＢＰＦ２２はフィルタ処理を行い（ステップＳ２０）、最大加速度検出部２３及び２４はそれぞれ最大加速度Ａ１及びＡ２を検出し（ステップＳ２１）、比率計算部２５は比率Ａ１／Ａ２を計算する（ステップＳ２２）。また、定数演算部５０は、定数α及びβを前述の方法で演算し（ステップＳ３０）、Ｓ波震度推定部４０はＰ波震度Ｖｐ、比率Ａ１／Ａ２、定数α、β及び設定されている定数Ｃｊに基づいて前記数３からＳ波震度Ｖｓを推定する（ステップＳ４０）。
【００４１】
なお、Ｐ波震度の計算（ステップＳ１０）、フィルタ処理（ステップＳ２０）〜比率計算（ステップＳ２２）及び定数演算（ステップＳ３０）の順番は任意に可変である。
【００４２】
ところで、上記ステップＳ１０の詳細は図８に示すようになっており、Ｐ波信号Ｖｅｐが加速度計３１１〜３１３に入力され（ステップＳ１１）、各加速度計３１１〜３１３でＰ波信号Ｖｅｐの３方向の加速度（水平動＃１、＃２及び上下動）が検出され（ステップＳ１２）、ベクトル合成部３２でベクトル合成され（ステップＳ１３）、合成された合成加速度が時間領域のフィルタ３３に入力されてフィルタ処理される（ステップＳ１４）。フィルタ処理された合成加速度ＡＣ０は変数算出部３４に入力され、変数算出部３４で前記所定式に従って変数ａ０が算出され（ステップＳ１５）、震度計算部３５は前記数２に基づいてＰ波震度Ｖｐを計算する（ステップＳ１６）。
【００４３】
また、上記ステップＳ３０の定数α、βの演算処理は、過去の地震データなどから観測点に依存しない定数α、βを求めるものであり、詳細は図９に示すようになっている。
【００４４】
先ずＰ波震度Ｖｐを入力して（ステップＳ３１）メモリ５１に記憶し（ステップＳ３２）、同様にＳ波震度Ｖｓを入力して（ステップＳ３３）メモリ５２に記憶し（ステップＳ３４）、演算部５５は数２を実行する（ステップＳ３５）。また、最新の観測データであるＰ波震度Ｖｐ及びＳ波震度Ｖｓは演算部５６に入力され（ステップＳ３６）、演算部５６は数２を実行し（ステップＳ３７）、修正部５７は演算部５５及び５６の演算結果から前述の計算方法を用いて定数α、βを求めて出力する（ステップＳ３８）。
【００４５】
上述のようにＳ波震度推定部３０で推定したＳ波震度Ｖｓを使用し、目的の震度に合わせた接点信号（ＯＮ／ＯＦＦ）や地震情報を出力させ、設備の制御などに使用する。現場地震計単独でも可能であるし、緊急地震速報と並列させることにより、緊急地震速報からの情報を補完すると共に、緊急地震速報では間に合わない近距離で発生した地震の震度予測にも有効である。
【００４６】
なお、上述では最大加速度Ａ１、Ａ２を検出してその比率Ａ１／Ａ２を関数ｆ（）に入力し、数３に従ってＳ波震度Ｖｓを推定しているが、最大加速度Ａ１、Ａ２の各関数ｆ１（）、ｆ２（）から下記数８に従って推定することも可能である。
（数８）
Ｖｓ＝Ｃｊ＋α１・ｆ１（Ａ１）＋α２・ｆ２（Ａ２）＋β・Ｖｐ
また、上述では２若しくは３個のバンドパスフィルタを用いて加速度を求め、Ｓ波震度Ｖｐの推定を行っているが、複数のバンドパスフィルタを用いて加速度Ａｉを求め、各加速度Ａｉに対応した定数αｉ、βｉ、関数ｆｉを用いて下記数９のような一般式でＳ波震度Ｖｓを推定することができる。
（数９）
Ｖｓ＝Ｃｊ＋αｉ・ｆｉ（Ａｉ）＋βｉ・Ｖｐ
【実施例】
【００４７】
Ｐ波に基づく震度ＶｐからＳ波の震度Ｖｓを推定するための補正として、観測結果を基に以下の数１０を求めた。
【００４８】
（数１０）
Ｖｓ＝0.824−0.307×log(A3hz/A40hz)+0.909×Ｖｐ
ただし、A3hzはＰ波加速度の３Hz成分の最初の２秒間の絶対値の和（平均的な振幅
）であり、A40hzはＰ波加速度の４０Hz成分の最初の２秒間の絶対値の和（平均的
な振幅）である。
観測点毎の補正として、定数項のみを変えて下記数１１を求めた。
（数１１）
Ｖｓ＝Ｃｊ−0.307×log(A3hｚ/A40hz)+0.909×Ｖｐ
この補正をした時と、補正をしなかった時のＲＭＳ残差を比較すると、次の表５のようになる。即ち、表５は、Ｐ波震度ＶｓからＳ波震度Ｖｓを推定したときのＲＭＳ残差の比較を示している。
【００４９】
【表５】
0005507903.tif 000006
観測点補正、発震機構の補正の両方を行った場合が、ＲＭＳ値が一番低い。観測点補正をするだけよりも、ＲＭＳが０．０５程度良くなっている。また、ボアホール（地中の地震計）でのＰ波震度→ボアホールのＳ波震度推定の方が、地表のＰ波震度→地表のＳ波震度推定よりもＲＭＳ値が低い。因みに地下のＰ波震度→地表のＳ波震度のＲＭＳ値はこれらよりも相対的に高い。即ち、地表の震度を推定するためには、地表のデータを使用するのが一番良い。これは、地下から地表までの地盤による増幅の非線形性によるもので、その補正が一筋縄ではいかないことを示している。
【符号の説明】
【００５０】
１０                    震度計
２０                    Ｐ波抽出部
２１                    低周波用バンドパスフィルタ（ＢＰＦ）
２２                    高周波用バンドパスフィルタ（ＢＰＦ）
２３、２４              最大加速度検出部
２５                    比率演算部
３０                    Ｐ波震度計算部
３１１〜３１３          加速度計
３２                    ベクトル合成部
３３                    フィルタ
３４                    変数算出部
３５                    震度計算部
４０                    Ｓ波震度推定部
５０                    定数演算部
５１，５２              メモリ
５５，５６              演算部
５７                    修正部
【図１】
0005507903.tif 000007
【図２】
0005507903.tif 000008
【図３】
0005507903.tif 000009
【図４】
0005507903.tif 000010
【図５】
0005507903.tif 000011
【図６】
0005507903.tif 000012
【図７】
0005507903.tif 000013
【図８】
0005507903.tif 000014
【図９】
0005507903.tif 000015