(57)【要約】【課題】本震直後の限られた地震観測データであってもリアルタイムに余震確率を推定することができる、余震予測装置、余震予測方法、および、プログラムを提供することを目的とする。【解決手段】本発明は、少なくとも本震から任意の経過時間までの地震観測データを記憶し、時間ｔにおける５０％の検出率のマグニチュードを示すパラメータμ（ｔ）を少なくとも含む検出率関数とグーテンベルグ・リヒター則の積による関数を、各学習期間における地震観測データに逐次フィッティングさせ、時間ｔに応じて変動するパラメータμ（ｔ）を含むパラメータの値を推定し、推定されたパラメータの値を導入した地震発生率モデルを用いて、余震発生率を予測する際に、μ（ｔ）の時間変動を、状態空間法に基づいて非パラメトリックに推定することを特徴とする。【選択図】図３
2014219315.tif 000002
【特許請求の範囲】
【請求項１】
記憶部と制御部を少なくとも備えた余震予測装置であって、
上記記憶部は、
少なくとも本震から任意の経過時間までの地震観測データを記憶し、
上記制御部は、
時間ｔにおける５０％の検出率のマグニチュードを示すパラメータμ（ｔ）を少なくとも含む検出率関数とグーテンベルグ・リヒター則の積に基づく関数を、各学習期間における上記地震観測データに逐次フィッティングさせ、時間ｔに応じて変動するパラメータμ（ｔ）を含むパラメータの値を推定するパラメータ推定手段と、
上記パラメータ推定手段により推定されたパラメータの値を導入した地震発生率モデルを用いて、余震発生率を予測する余震発生率予測手段と、
を備えた余震予測装置において、
上記パラメータ推定手段は、
上記μ（ｔ）の時間変動を、状態空間法に基づいて非パラメトリックに推定すること
を特徴とする余震予測装置。
【請求項２】
請求項１に記載の余震予測装置において、
上記μ（ｔ）の時間変動は、
上記状態空間法に基づいて離散化した以下の式により表されること
【数１】
2014219315.tif 000028
を特徴とする余震予測装置。
【請求項３】
請求項２に記載の余震予測装置において、
上記μ（ｔ）に対する事前分布として以下の式を導入し、上記μ（ｔ）の時間変動に対してペナルティーを課すこと
【数２】
2014219315.tif 000029
（ここで、Ｖは変動の度合いを示す。）
を特徴とする余震予測装置。
【請求項４】
記憶部と制御部を少なくとも備えたコンピュータにおいて実行される余震予測方法であって、
上記記憶部は、
少なくとも本震から任意の経過時間までの地震観測データを記憶し、
上記制御部において実行される、
時間ｔにおける５０％の検出率のマグニチュードを示すパラメータμ（ｔ）を少なくとも含む検出率関数とグーテンベルグ・リヒター則の積による関数を、各学習期間における上記地震観測データに逐次フィッティングさせ、時間ｔに応じて変動するパラメータμ（ｔ）を含むパラメータの値を推定するパラメータ推定ステップと、
上記パラメータ推定ステップにより推定されたパラメータの値を導入した地震発生率モデルを用いて、余震発生率を予測する余震発生率予測ステップと、
を含む余震予測装置において、
上記パラメータ推定ステップは、
上記μ（ｔ）の時間変動を、状態空間法に基づいて非パラメトリックに推定すること
を特徴とする余震予測方法。
【請求項５】
少なくとも本震から任意の経過時間までの地震観測データを記憶する記憶部と制御部を少なくとも備えたコンピュータに実行させるためのプログラムであって、
上記制御部において、
時間ｔにおける５０％の検出率のマグニチュードを示すパラメータμ（ｔ）を少なくとも含む検出率関数とグーテンベルグ・リヒター則の積による関数を、各学習期間における上記地震観測データに逐次フィッティングさせ、時間ｔに応じて変動するパラメータμ（ｔ）を含むパラメータの値を推定するパラメータ推定ステップと、
上記パラメータ推定ステップにより推定されたパラメータの値を導入した地震発生率モデルを用いて、余震発生率を予測する余震発生率予測ステップと、
を実行させるためのプログラムであって、
上記パラメータ推定ステップは、
上記μ（ｔ）の時間変動を、状態空間法に基づいて非パラメトリックに推定すること
を特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、余震予測装置、余震予測方法、および、プログラムに関する。
【背景技術】
【０００２】
従来、地震の発生頻度等を表現する数理モデルが提案されている。
【０００３】
例えば、マグニチュードの大きさと、各マグニチュードの地震（余震）の発生頻度の対数との間には、比例関係があることを示す、グーデンベルク・リヒター則が発見されている。
【０００４】
また、１８９４年、大森房吉は、本震後の余震の回数が時間に関して双曲関数に従って減衰することを示す式（大森則）を発表しており、１９５７年には、宇津徳治が、これを改良し、一般にべき乗則に従って減衰することを示す大森・宇津則（改良大森公式）を発表している。
【０００５】
１９８９年、リーゼンベルクとジョーンズは、グーテンベルク・リヒター則と大森・宇津則を組み合わせて、将来の任意の時間、各マグニチュードにおける、余震の確率を予測する方法を発表した。
【０００６】
これら地震発生率モデルの各パラメータは、地震発生率モデルを地震観測データにフィッティングさせることにより求めることができるが、特に本震から１日以内では、複数の地震波が重なり合う影響で個々の余震を良く区別することができず、余震の検出率が低下するため、得られるデータは非常に不完全である。そのためフィッティングさせる際に検出率を加味する必要がある。
【０００７】
非特許文献１では、本震からの経過時間に従ってパラメトリックに変動する検出率の関数を導入し、地震観測データの不完全性を考慮に入れて余震確率を推定することが開示されている。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】Ｙｏｓｉｈｉｋｏ  Ｏｇａｔａ，  Ｋｏｉｃｈｉ  Ｋａｔｓｕｒａ  “Ｉｍｍｅｄｉａｔｅ  ａｎｄ  ｕｐｄａｔｅｄ  ｆｏｒｅｃａｓｔｉｎｇ  ｏｆ  ａｆｔｅｒｓｈｏｃｋ  ｈａｚａｒｄ”  Ｇｅｏｐｈｙｓｉｃａｌ  Ｒｅｓｅａｒｃｈ  Ｌｅｔｔｅｒｓ，  Ｖｏｌｕｍｅ  ３３，  Ｉｓｓｕｅ  １０，  Ｍａｙ  ２００６
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、非特許文献１等に記載の従来の方法では、本震直後の地震観測データを用いると、検出率関数をフィッティングさせることができずパラメータが発散してしまうことがあるという問題点を有していた。
【００１０】
本発明は、上記問題点に鑑みてなされたもので、本震直後の限られた地震観測データであってもリアルタイムに余震確率を推定することができる、余震予測装置、余震予測方法、および、プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１１】
このような目的を達成するため、本発明の余震予測装置は、記憶部と制御部を少なくとも備えた余震予測装置であって、上記記憶部は、少なくとも本震から任意の経過時間までの地震観測データを記憶し、上記制御部は、時間ｔにおける５０％の検出率のマグニチュードを示すパラメータμ（ｔ）を少なくとも含む検出率関数とグーテンベルグ・リヒター則の積に基づく関数を、各学習期間における上記地震観測データに逐次フィッティングさせ、時間ｔに応じて変動するパラメータμ（ｔ）を含むパラメータの値を推定するパラメータ推定手段と、上記パラメータ推定手段により推定されたパラメータの値を導入した地震発生率モデルを用いて、余震発生率を予測する余震発生率予測手段と、を備えた余震予測装置において、上記パラメータ推定手段は、上記μ（ｔ）の時間変動を、状態空間法に基づいて非パラメトリックに推定することを特徴とする。
【００１２】
また、本発明の余震予測装置は、上記記載の余震予測装置において、上記μ（ｔ）の時間変動は、離散化した以下の式により表されることを特徴とする。
【数１】
2014219315.tif 000003
【００１３】
また、本発明の余震予測装置は、上記記載の余震予測装置において、上記μ（ｔ）に対する事前分布として以下の式を導入し、上記μ（ｔ）の時間変動に対してペナルティーを課すことを特徴とする。
【数２】
2014219315.tif 000004
（ここで、Ｖは変動の度合いを示す。）
【００１４】
また、本発明の余震予測方法は、記憶部と制御部を少なくとも備えたコンピュータにおいて実行される余震予測方法であって、上記記憶部は、少なくとも本震から任意の経過時間までの地震観測データを記憶し、上記制御部において実行される、時間ｔにおける５０％の検出率のマグニチュードを示すパラメータμ（ｔ）を少なくとも含む検出率関数とグーテンベルグ・リヒター則の積による関数を、各学習期間における上記地震観測データに逐次フィッティングさせ、時間ｔに応じて変動するパラメータμ（ｔ）を含むパラメータの値を推定するパラメータ推定ステップと、上記パラメータ推定ステップにより推定されたパラメータの値を導入した地震発生率モデルを用いて、余震発生率を予測する余震発生率予測ステップと、を含む余震予測装置において、上記パラメータ推定ステップは、上記μ（ｔ）の時間変動を、状態空間法に基づいて非パラメトリックに推定することを特徴とする。
【００１５】
また、本発明のプログラムは、少なくとも本震から任意の経過時間までの地震観測データを記憶する記憶部と制御部を少なくとも備えたコンピュータに実行させるためのプログラムであって、上記制御部において、時間ｔにおける５０％の検出率のマグニチュードを示すパラメータμ（ｔ）を少なくとも含む検出率関数とグーテンベルグ・リヒター則の積による関数を、各学習期間における上記地震観測データに逐次フィッティングさせ、これによって時間ｔに応じて変動するパラメータμ（ｔ）を含むパラメータの値を推定するパラメータ推定ステップと、上記パラメータ推定ステップにより推定されたパラメータの値を導入した地震発生率モデルを用いて、余震発生率を予測する余震発生率予測ステップと、を実行させるためのプログラムであって、上記パラメータ推定ステップは、上記μ（ｔ）の時間変動を、状態空間法に基づいて非パラメトリックに推定することを特徴とする。
【００１６】
また、本発明は記録媒体に関するものであり、上記記載のプログラムを記録したことを特徴とする。
【発明の効果】
【００１７】
この発明によれば、本震直後の限られた地震観測データであってもリアルタイムに余震確率を推定することができるという効果を奏する。より具体的には、先行技術で提案された検出率のパラメトリック関数は、地震の検出率がなだらかに上昇するパラメトリックな関数を想定していたため、大きな本震が発生してから１日以内の地震観測データにおいて、検出率に揺らぎがある場合、柔軟にフィッティングさせることができずパラメータが発散してしまうことがあるという問題があったが、本発明によれば、検出率の時間変動を、非パラメトリックに推定するので、本震直後の不均質な地震観測データであっても、最適にフィッティングさせることが可能になる。
【００１８】
また、本発明によれば、時間を変数とする検出率の関数を離散化させ、それを状態空間法で推定することにより、特定の関数形を仮定することなくデータに対して適応的な推定が可能になるという効果を奏する。
【図面の簡単な説明】
【００１９】
【図１】図１は、グーテンベルグ・リヒター則を示すグラフである。
【図２】図２は、正規分布の累積分布関数による検出率関数を示すグラフである。
【図３】図３は、本実施の形態が適用される余震予測装置１００の一例を示すブロック図である。
【図４】図４は、本実施の形態における余震予測装置１００の処理の一例を示すフローチャートである。
【図５】図５は、東北地方太平洋沖地震前後の震源地分布を示す図である。
【図６】図６は、東北地方太平洋沖地震前後のマグニチュードの時系列データを示す図である。
【図７】図７は、検出率関数とグーテンベルグ・リヒター則の積に基づく式（６）を、地震観測データにフィッティングさせたグラフである。
【図８】図８は、各学習期間（本震から３時間、６時間、１２時間、２４時間）により得られたμ（ｔ）の値を示すグラフである。
【図９】図９は、各期間（３−６ｈ，６−１２ｈ，１２−２４ｈ）において、観測されたマグニチュードの正規化ヒストグラムと、各地震が起きた時刻における推定された瞬間マグニチュード分布の関係を示す図である。
【図１０】図１０は、本震後１日までの学習期間で推定されたパラメータを用いて、１日後から３０日までの期間、マグニチュード５以上の余震の頻度を予測した結果を示す図である。
【図１１】図１１は、本震後１日までの学習期間で推定されたパラメータを用いて、１日後から３０日までの期間、マグニチュード６以上の余震の頻度をそれぞれ予測した結果を示す図である。
【図１２】図１２は、短い学習期間において推定されたパラメータを用いて、その後の期間、マグニチュード５以上の余震の頻度を予測した結果を示す図である。
【図１３】図１３は、短い学習期間において推定されたパラメータを用いて、その後の期間、マグニチュード６以上の余震の頻度を予測した結果を示す図である。
【発明を実施するための形態】
【００２０】
以下に、本発明にかかる余震予測装置、余震予測方法、および、プログラムの実施の形態を図面に基づいて詳細に説明する。なお、この実施の形態によりこの発明が限定されるものではない。
【００２１】
以下、本発明の実施の形態の概要について説明し、その後、本実施の形態の構成および処理等の詳細について説明する。
【００２２】
［本発明の実施の形態の概要］
本発明の実施の形態の概要として、以下に、本実施の形態が開発された背景、先行技術例、本発明の実施の形態の概要について順に説明する。
【００２３】
［背景］
大きな地震の後にはおびただしい数の余震が起こる。特に本震後一日間には、追加の被害を引き起こしうる大きな余震が起こる確率が高く、すべての大きな余震のうち約半数がこの期間に起こる。そのため本震後できるだけ早く、余震の予測を行うことが望ましい。
【００２４】
その一方で、本震の直後に正確な予測を立てることは非常に難しく、実際、気象庁は本震後一日以上たってから予測を行っている。これは、本震とその後の多くの余震の地震波が重なることにより、それぞれの地震の発生場所、時間、マグニチュードを決めることが難しく、実際に起こった地震のうちごく一部の地震しか観測できないからである。
【００２５】
この困難を克服するために、観測からもれた余震も含め実際に起こった余震の活動を、余震の検出率を考慮して推定・予測することが提案されている（非特許文献１）。
【００２６】
この先行技術例（非特許文献１）では、まず余震の検出率を推定し、それと既存の予測モデルが組み合わされることにより、本震の直後から精度の高い予測が可能であることが２００３年の宮城県沖地震の余震活動を例にして示されている。以下に、非特許文献１の先行技術例について説明し、その後、本実施の形態による改良方法について説明する。
【００２７】
［先行技術例］
先行技術例（非特許文献１）では、実際に起こった余震のモデルとしてＲｅａｓｅｎｂｅｒｇ−Ｊｏｎｅｓのモデル（Ｓｃｉｅｎｃｅ  １９８９，Ｖｏｌ．２４３，ｎｏ．４８９５，ｐｐ．１１７３−１１７６）を用いている。このモデルでは、本震からの経過時間ｔにおけるマグニチュードＭの余震の発生率が、以下の式で与えられる。
【数３】
2014219315.tif 000005
ここで、Ｋ，ｐ，ｃは余震の減衰を表す大森・宇津則（改良大森公式）のパラメータであり、ｂはマグニチュードの頻度を表すグーテンベルグ・リヒター則のパラメータである。ここで、図１は、グーテンベルグ・リヒター則を示すグラフである。横軸はマグニチュード（Ｍ）を示し、縦軸は発生頻度ｎ（Ｍ）の対数を示す。
【００２８】
図１に示すように、グーテンベルグ・リヒター則は、マグニチュード（Ｍ）の大きさと、各マグニチュードの発生頻度（ｎ（Ｍ））の対数との間には、比例関係があることを示し、ｌｏｇｎ（Ｍ）＝ａ−ｂＭで表される（以下の式（５）参照）。すなわち、図１のグラフの傾きは−ｂであり、縦軸の切片はａである。
【００２９】
式（１）のモデルからマグニチュードＭｆ以上の余震の発生率は、以下の式で与えられる。
【数４】
2014219315.tif 000006
【００３０】
実際に起こった余震のモデル（１）のパラメータを、観測された余震のデータから推定するために地震の検出率を導入する。余震の検出率は、マグニチュードが小さいほど低く、マグニチュードが大きいほど高いというようにマグニチュードに依存し、その関係は、以下の正規分布の累積分布関数でモデル化される。
【数５】
2014219315.tif 000007
ここでパラメータμは検出率が５０％であるマグニチュードを表し、σは地震が部分的に観測されるマグニチュードの幅を表している。ここで、図２は、正規分布の累積分布関数による検出率関数を示すグラフである。横軸はマグニチュードを示し、縦軸は検出率を示す。
【００３１】
図２に示すように、μは、検出率０．５（５０％）のときのマグニチュード（この例ではμ＝４程度）である。
【００３２】
なお、余震の検出率には、このマグニチュードへの依存性に加えて、時間への依存性も存在する。すなわち、検出率は本震後に非常に低く、時間が経つにつれて通常値に戻っていくという時間への依存性がある。それを表現するために、パラメータμが時間の関数μ（ｔ）であるとする。
【００３３】
尾形・桂による先行技術例（非特許文献１）では、以下のように、μ０，μ１，α，γというパラメータをもつパラメトリックな時間依存性の検出率関数μ（ｔ）のモデルが提案されている。
【数６】
2014219315.tif 000008
【００３４】
検出率関数のパラメータは、観測された余震のマグニチュードのデータからフィッティングにより推定可能である。モデル（１）では、実際に起こっている地震のマグニチュードの分布は、以下のグーテンベルグ・リヒター則に従っていると仮定されている。
【数７】
2014219315.tif 000009
【００３５】
検出される地震のマグニチュードの確率分布は、実際に起こった地震のマグニチュードの分布（５）と検出率（３）の積に比例し、以下の式で与えられる。
【数８】
2014219315.tif 000010
【数９】
2014219315.tif 000011
【００３６】
以上の手続きにより、検出率関数にて時間変動する検出率が推定できる。
【００３７】
次に推定された検出率を用いて、予測モデル（１）の残りのパラメータＫ，ｐ，ｃの推定を行う。検出された地震の発生率は、実際に起こっている地震の確率強度（１）と検出率関数（３）の積である以下の式で表される。
【数１０】
2014219315.tif 000012
【００３８】
ここから、大森・宇津則のＫ，ｐ，ｃに対する対数尤度が以下の式で表され、これを最大化することによって、先行技術例ではパラメータ値を得ている。
【数１１】
2014219315.tif 000013
【００３９】
［本実施の形態による方法］
上述した先行技術例（非特許文献１）では、時間変動する余震の検出率を推定するために、式（４）にて示したように、パラメトリックな統計モデルが用いられている。しかしながら、このモデルのパラメータを短い余震時系列から推定しようとするとき、特に本震直後の限られた時間の地震観測データしか得られていない場合に、パラメータが発散してしまい最適なパラメータ値が得られないことがあった。
【００４０】
本願発明者らは、鋭意検討の結果、安定した検出率の推定を行うために、パラメトリックな統計モデルが用いることなく、状態空間法に基づく検出率の推定手法を導入することが適切であると想起するに至った。
【００４１】
本実施の形態では、安定して検出率関数を推定するために状態空間法を用いて、μ（ｔ）を推定する。状態空間法を用いることにより、特定の関数形を仮定することなく地震観測データに対して適応的な推定が可能になる。
【００４２】
本実施の形態では、解析を単純にするために時間連続関数であるμ（ｔ）を次式のように離散化する。ここで、ｔｉはｉ番目に検出された余震の発生時刻である。
【数１２】
2014219315.tif 000014
【００４３】
これにより、μ（ｔ）の推定問題は、有限個のパラメータμ＝［μ１，μ２，．．．，μＮ］Ｔを推定する問題になる。
【００４４】
状態空間法では、μ（ｔ）の変動に対してペナルティーを課して推定を行う。このために、μ（ｔ）に対する事前分布を次式のように導入する。
【数１３】
2014219315.tif 000015
【００４５】
ここで、初期分布はＰ（μ１，μ２）＝ｃｏｎｓｔ（一定）である。これによりμ（ｔ）の変動が制限され、予測に最適なスムーズな関数の推定が可能になる。変動の度合いはパラメータＶによって決まる。そして与えられたμからマグニチュード系列Ｍ＝［Ｍ１，Ｍ２，．．．，ＭＮ］Ｔが発生する確率は、式（６）から、次式で与えられる。
【数１４】
2014219315.tif 000016
【００４６】
このとき推定の目標はパラメータＶ，β，σの最適な推定のもとで、マグニチュード系列Ｍから、パラメータμの事後確率Ｐβ，σ，Ｖ（μ｜Ｍ）を最大化するμ＝μ＊を見つけることである。事後確率は、ベイズの公式から次式で与えられる。
【数１５】
2014219315.tif 000017
【００４７】
事後分布を最大にするには、μ＊をＮｅｗｔｏｎ法等により見つける。ただし、これらの推定はパラメータＶ，β，σの推定に大きく依存する。これらのパラメータの事後確率Ｐ（Ｖ，β，σ│Ｍ）もベイズの公式から次式で与えられる。
【数１６】
2014219315.tif 000018
【００４８】
事前分布Ｐ（Ｖ，β，σ）は、過去のデータから決めておいてもよい。一例として、Ｐ（Ｖ，β，σ）は、βの正規分布であるとする。Ｐ（Ｖ，β，σ│Ｍ）を最大にするＶ，β，σは、ＥＭ法（期待値最大化法）等を用いてフィッティングさせることにより求められる。
【００４９】
以上が本実施の形態の概要である。なお、以上は、本発明が具現化される一形態を、本実施の形態の概要として記載したものであり、本発明は、上述または後述する実施の形態に限定されるものではない。以下の実施の形態においては、上述した本実施の形態をハードウェアとソフトウェアの協働により実現させるため、本実施の形態の構成および処理等の詳細について説明する。
【００５０】
［余震予測装置の構成］
次に、本実施の形態の余震予測装置の構成について図３を参照して説明する。図３は、本実施の形態が適用される余震予測装置１００の一例を示すブロック図であり、該構成のうち本実施の形態に関係する部分を概念的に示している。
【００５１】
図３に示すように、本実施の形態における余震予測装置１００は、概略的に、制御部１０２と記憶部１０６を少なくとも備え、本実施の形態において、更に、入出力制御インターフェース部１０８と通信制御インターフェース部１０４を備える。ここで、制御部１０２は、余震予測装置１００の全体を統括的に制御するＣＰＵ等である。また、通信制御インターフェース部１０４は、通信回線等に接続されるルータ等の通信装置（図示せず）に接続されるインターフェースであり、入出力制御インターフェース部１０８は、入力装置１１２や出力装置１１４に接続されるインターフェースである。また、記憶部１０６は、各種のデータベースやテーブルなどを格納する装置である。これら余震予測装置１００の各部は任意の通信路を介して通信可能に接続されている。更に、この余震予測装置１００は、ルータ等の通信装置および専用線等の有線または無線の通信回線を介して、ネットワーク３００に通信可能に接続されている。
【００５２】
記憶部１０６に格納される各種のデータベースやテーブル（リアルタイム観測ファイル１０６ａ、パラメータファイル１０６ｂ、予測結果ファイル１０６ｃ等）は、固定ディスク装置等のストレージ手段である。例えば、記憶部１０６は、各種処理に用いる各種のプログラム、テーブル、ファイル、データベース、および、ウェブページ等を格納する。
【００５３】
これら記憶部１０６の各構成要素のうち、リアルタイム観測ファイル１０６ａは、地震観測データを記憶する地震観測データ記憶手段である。リアルタイム観測ファイル１０６ａは、少なくとも本震から任意の経過時間までのデータを記憶し、好適には、本震からの経過時間に応じてリアルタイムに地震観測データが更新される。例えば、余震予測装置１００の制御部１０２が、ネットワーク３００を介して、気象庁等の外部システム２００から、随時更新される地震観測データ（気象庁カタログデータ、気象庁一元化震源データ等）をダウンロードすることにより、リアルタイム観測ファイル１０６ａに記憶された地震観測データを更新してもよい。また、余震予測装置１００の制御部１０２が、地震計等の入力装置１１２から、あるいは、地震観測網を構成する地震計等の外部システム２００から、リアルタイムに地震波形データを取得して、発震時刻と震源地の決定とマグニチュードの算出を行うことにより地震観測データを取得して、リアルタイム観測ファイル１０６ａの地震観測データを更新してもよい。
【数１７】
2014219315.tif 000019
【００５４】
また、パラメータファイル１０６ｂは、各種パラメータ（係数値等）を記憶するパラメータ記憶手段である。例えば、パラメータファイル１０６ｂは、グーテンベルグ・リヒター則のパラメータ（例えば、式（５）のβ（＝ｂｌｎ１０）等）、離散化されたμ（ｔ）のパラメータ（例えば、式（９）のインデックスｉと経過時間ｔの関連付けデータ、μｉ（ｉ＝１，２，．．．，Ｎ））、余震の減衰に関するパラメータ（例えば、式（１）中の大森・宇津則に関するパラメータＫ，ｐ，ｃ）、マグニチュードに関する検出率関数（正規分布の累積分布関数等）のパラメータ（例えば、式（３）のσ）、ガウス分布に関するパラメータ（例えば、式（１０）のＰ（μ１，μ２），Ｖ）等の各種パラメータの値（初期値、更新値等）を記憶してもよい。
【００５５】
また、予測結果ファイル１０６ｃは、余震発生率予測部１０２ｂによる余震発生率の予測結果を記憶する予測結果記憶手段である。例えば、予測結果ファイル１０６ｃは、パラメータの値が導入された地震発生率モデル（大森・宇津則）、Ｒｅａｓｅｎｂｅｒｇ−Ｊｏｎｅｓのモデル等）に基づく計算結果（例えば、３日以内にマグニチュード５以上の余震が発生する確率○○％など）を記憶してもよい。
【００５６】
入出力制御インターフェース部１０８は、入力装置１１２や出力装置１１４の制御を行う。ここで、出力装置１１４としては、モニタ（家庭用テレビを含む）の他、スピーカを用いることができる（なお、以下においては出力装置１１４をモニタとして記載する場合がある）。また、入力装置１１２としては、キーボード、マウス、およびマイク等のほか、地震計や３軸加速度センサーなど地震の動きを計測する手段であってもよい。なお、複数台の入力装置１１２が余震予測装置１００に接続されてもよい。
【００５７】
また、図３において、制御部１０２は、ＯＳ（Ｏｐｅｒａｔｉｎｇ  Ｓｙｓｔｅｍ）等の制御プログラムや、各種の処理手順等を規定したプログラム、および、所要データを格納するための内部メモリを有する。そして、制御部１０２は、これらのプログラム等により、種々の処理を実行するための情報処理を行う。制御部１０２は、機能概念的に、パラメータ推定部１０２ａ、および、余震発生率予測部１０２ｂを備える。概略的に、パラメータ推定部１０２ａは、数理モデルを実測値（現在までの地震観測データ）にあてはめてパラメータの値を推定するのに対し、余震発生率予測部１０２ｂは、推定されたパラメータの値を導入した数理モデルを用いて、将来の余震発生率を予測する機能を有する。
【００５８】
すなわち、パラメータ推定部１０２ａは、パラメータμ（ｔ）を少なくとも含む検出率関数とグーテンベルグ・リヒター則の積に基づく関数を、各学習期間における地震観測データに逐次フィッティングさせ、時間ｔに応じて変動するパラメータμ（ｔ）を含むパラメータの値を推定するパラメータ推定手段である。なお、パラメータ推定部１０２ａは、推定したパラメータの値や、新たに取得された地震観測データに基づいて更新されたパラメータの値を、パラメータファイル１０６ｂに格納する。
【００５９】
ここで、本実施の形態において、パラメータ推定部１０２ａは、μ（ｔ）の時間変動を、状態空間法に基づいて非パラメトリックに推定する。より具体的には、パラメータ推定部１０２ａは、時間を変数とするμ（ｔ）の関数を、状態空間法に基づいて以下の式により離散化させてもよい。例えば、パラメータ推定部１０２ａは、地震観測データから得られる有限個（Ｎ）のマグニチュード系列Ｍ＝［Ｍ１，Ｍ２，．．．，ＭＮ］Ｔに基づいて、対応する有限個（Ｎ）のパラメータμ＝［μ１，μ２，．．．，μＮ］Ｔを推定する。
【数１８】
2014219315.tif 000020
ここで、ｉは検出された地震のインデックスであり、ｔは本震からの経過時間である。
【００６０】
なお、パラメータ推定部１０２ａは、ここで、μ（ｔ）に対する事前分布として以下の式を導入してもよい。これにより、μの時間変動に対して適切なペナルティーを課して推定することができる。
【数１９】
2014219315.tif 000021
（ここで、Ｖは変動の度合いを示す。）
【００６１】
ここで、パラメータ推定部１０２ａは、各種のパラメータを推定するために、２段階でフィッティングを行ってもよい。すなわち、第１段階では、パラメータ推定部１０２ａは、マグニチュードに依存する関数（パラメータμ（ｔ）を少なくとも含む検出率関数とグーテンベルグ・リヒター則の積に基づくマグニチュード分布を表す式（６）等）を、各学習期間における地震観測データにフィッティングさせることにより、学習期間ごとに、マグニチュードに依存する関数のパラメータを推定してもよい。なお、この第１段階では、μ（ｔ）は定数として扱われる。
【００６２】
そして、第２段階では、パラメータ推定部１０２ａは、時間に依存する関数（時間に依存する検出率に応じた、余震の発生確率の減衰に関する関数等）を、本震発生から現在得られている経過時間までの地震観測データにフィッティングさせることにより、時間に依存する関数に関するパラメータを推定してもよい。例えば、パラメータ推定部１０２ａは、過去の地震観測データから定めた事前分布Ｐ（Ｖ，β，σ）から、事後分布Ｐ（Ｖ，β，σ│Ｍ）が最大になるＶ，β，σを、ＥＭ法（期待値最大化法）等を用いてフィッティングさせることにより求めてもよい。
【００６３】
また、余震発生率予測部１０２ｂは、パラメータ推定部１０２ａにより推定されたパラメータの値を導入した地震発生率モデルを用いて、余震発生率を予測する余震発生率予測手段である。例えば、余震発生率予測部１０２ｂは、推定されたパラメータの値を導入した地震発生率モデル（大森・宇津則、Ｒｅａｓｅｎｂｅｒｇ−Ｊｏｎｅｓのモデル等）をグラフ化して将来の余震発生率をグラフに表してもよく、３日以内にマグニチュード５以上の余震発生率を積分等により算出してもよい。なお、余震発生率予測部１０２ｂは、予測結果データを、予測結果ファイル１０６ｃに格納してもよく、モニタ等の出力装置１１４に出力（表示等）させてもよく、外部システム２００（利用者端末等）に送信してもよい。
【００６４】
以上が、本実施の形態における余震予測装置１００の構成の一例である。なお、余震予測装置１００は、ネットワーク３００を介して外部システム２００に接続されてもよい。この場合、通信制御インターフェース部１０４は、余震予測装置１００とネットワーク３００（またはルータ等の通信装置）との間における通信制御を行う。すなわち、通信制御インターフェース部１０４は、他の端末と通信回線を介してデータを通信する機能を有する。また、ネットワーク３００は、余震予測装置１００と外部システム２００とを相互に接続する機能を有し、例えば、インターネット等である。
【００６５】
また、外部システム２００は、ネットワーク３００を介して、余震予測装置１００と相互に接続され、各種パラメータに関する外部データベースや、接続された情報処理装置に余震予測方法を実行させるためのプログラム等を提供する機能を有する。
【００６６】
ここで、外部システム２００は、ＷＥＢサーバやＡＳＰサーバ等として構成していてもよい。また、外部システム２００のハードウェア構成は、一般に市販されるワークステーション、パーソナルコンピュータ等の情報処理装置およびその付属装置により構成していてもよい。また、外部システム２００の各機能は、外部システム２００のハードウェア構成中のＣＰＵ、ディスク装置、メモリ装置、入力装置、出力装置、通信制御装置等およびそれらを制御するプログラム等により実現される。
【００６７】
以上で、本実施の形態の構成の説明を終える。
【００６８】
［余震予測装置１００の処理］
次に、このように構成された本実施の形態における余震予測装置１００の処理の一例について、以下に図４を参照して詳細に説明する。図４は、本実施の形態における余震予測装置１００の処理の一例を示すフローチャートである。
【００６９】
図４に示すように、まず、パラメータ推定部１０２ａは、最新の地震観測データを取得する（ステップＳＡ−１）。例えば、パラメータ推定部１０２ａは、ネットワーク３００を介して、気象庁等の外部システム２００から、随時更新される地震観測データ（気象庁カタログデータ、気象庁一元化震源データ等）をダウンロードすることにより、リアルタイム観測ファイル１０６ａに最新の地震観測データを格納してもよい。また、パラメータ推定部１０２ａは、地震観測網を構成する、地震計等の入力装置１１２や外部システム２００から、リアルタイムに地震波形データを取得して、発震時刻と震源地の決定とマグニチュードの算出を行うことにより地震を特定して地震観測データを取得してもよい。
【数２０】
2014219315.tif 000022
【００７０】
そして、図３のパラメータ推定部１０２ａは、数理モデルを地震観測データにフィッティングさせる処理を行う（ステップＳＡ−２）。例えば、パラメータ推定部１０２ａは、パラメータμ（ｔ）を少なくとも含む検出率関数とグーテンベルグ・リヒター則の積に基づく以下の式（６）を、各学習期間における地震観測データであるマグニチュード系列Ｍ＝［Ｍ１，Ｍ２，．．．，ＭＮ］Ｔに、逐次フィッティングさせることにより、各パラメータμ＝［μ１，μ２，．．．，μＮ］Ｔを推定する。
【数２１】
2014219315.tif 000023
【００７１】
なお、このとき、パラメータ推定部１０２ａは、μ（ｔ）の変動を、ガウス分布に従うように制限する。すなわち、パラメータ推定部１０２ａは、μ（ｔ）の変動に対してペナルティーを課しながら推定を行う。μ（ｔ）に対する事前分布は、次式で表現される。
【数２２】
2014219315.tif 000024
初期分布Ｐ（μ１，μ２）は一定であり、Ｖは変動の度合いを示す。
【００７２】
与えられたμからマグニチュード系列Ｍ＝［Ｍ１，Ｍ２，．．．，ＭＮ］Ｔが発生する確率は、次式で表現される。
【数２３】
2014219315.tif 000025
【００７３】
パラメータ推定部１０２ａは、パラメータＶ，β，σの推定を最適化し、マグニチュード系列Ｍから、パラメータμの事後確率Ｐβ，σ，Ｖ（μ｜Ｍ）を最大化するμ＝μ＊を求めるため、ベイズ法を用いてもよい。事後確率は、ベイズの公式から次式で表現される。
【数２４】
2014219315.tif 000026
【００７４】
パラメータ推定部１０２ａは、事後分布を最大にするために、μ＊をＮｅｗｔｏｎ法等により求めてもよい。一方で、パラメータＶ，β，σの推定を最適化するため、パラメータ推定部１０２ａは、ベイズ法による次式の事後確率Ｐ（Ｖ，β，σ│Ｍ）を最大にするＶ，β，σを、ＥＭ法（期待値最大化法）等を用いてフィッティングさせることにより求める。ここで、Ｐ（Ｖ，β，σ）は、βの正規分布であるとする。
【数２５】
2014219315.tif 000027
【００７５】
再び図４に戻り、図３のパラメータ推定部１０２ａは、ステップＳＡ−２において新たに得られたパラメータの値を、パラメータファイル１０６ｂに格納することにより更新する（ステップＳＡ−３）。
【００７６】
そして、余震発生率予測部１０２ｂは、パラメータファイル１０６ｂから更新されたパラメータ値を地震発生率モデルに代入し、この地震発生率モデルに基づいて余震発生率を予測する（ステップＳＡ−４）。例えば、余震発生率予測部１０２ｂは、推定されたパラメータの値を導入した地震発生率モデル（大森・宇津則、Ｒｅａｓｅｎｂｅｒｇ−Ｊｏｎｅｓのモデル等）をグラフ化して、本震から現在を経て将来までの余震発生率をグラフに表してもよい。また、余震発生率予測部１０２ｂは、推定されたパラメータの値を導入した地震発生率モデルから、３日以内にマグニチュード５以上の余震が発生する確率を、積分等により算出してもよい。なお、余震発生率予測部１０２ｂは、予測結果データを、予測結果ファイル１０６ｃに格納してもよく、モニタ等の出力装置１１４に出力（表示等）させてもよく、外部システム２００（利用者端末等）に送信してもよい。
【００７７】
そして、余震予測装置１００の制御部１０２は、処理を終えるか否かを判定する（ステップＳＡ−５）。例えば、余震予測装置１００の制御部１０２は、本震から一定期間経過した場合に処理を終えると判定してもよく、利用者により入力部１１２を介して処理の終了が指示された場合に、処理を終了してもよい。
【００７８】
処理を終了しないと判定された場合（ステップＳＡ−５，Ｎｏ）、余震予測装置１００の制御部１０２は、上述したステップＳＡ−１〜ＳＡ−５の処理を繰返す。一方、処理を終了すると判定された場合（ステップＳＡ−５，Ｙｅｓ）、余震予測装置１００の制御部１０２は、処理を終了する。
【００７９】
以上が、本実施形態における余震予測装置１００の処理の一例である。
【００８０】
［実施例］
つづいて、上述した実施の形態を、２０１１年の東北地方太平洋沖地震（Ｍ９．０）の余震データに適用した実施例について、以下に図５〜図１３を参照して説明する。
【００８１】
ここで、図５は、東北地方太平洋沖地震前後の震源地分布を示す図である。なお、本震の震源地は星印で表されている。また、図６は、東北地方太平洋沖地震前後のマグニチュードの時系列データを示す図である。横軸は、本震の発生時点を０とする時刻であり、縦軸はマグニチュードを示す。これら地震観測データは、ＰＤＥ／ＮＥＩＣのカタログデータから取得した。
【００８２】
図７は東北地方太平洋沖地震の前の地震データ化に対して、検出率関数とグーテンベルグ・リヒター則の積に基づく式（６）を、μ（ｔ）が一定として、フィッティングさせたもののグラフである。横軸はマグニチュードであり、縦軸は発生頻度である。実線は、検出率関数とグーテンベルグ・リヒター則の積に基づく式を表し、黒丸は、マグニチュードの頻度を表す。モデルと観測が良く一致していることから検出率関数（３）の妥当性がわかる。
【００８３】
図８は、各学習期間（本震から３時間、６時間、１２時間、２４時間）の地震観測データに対して、状態空間法で推定されたμ（ｔ）の値を示すグラフである。黒点は観測された各地震を表す。図８の点線は学習期間３時間、一点鎖線は６時間、二点鎖線は１２時間、実線は２４時間で得られたμ（ｔ）の値の推移を示す。灰色の丸は、検出された各地震を表している。図８に示すように、本震直後の検出率の変動に揺らぎがある場合でも、適切に検出率を推定できることが分かった。
【００８４】
ここで、図９は、各期間（３−６ｈ，６−１２ｈ，１２−２４ｈ）において、観測されたマグニチュードの正規化ヒストグラムと、各地震が起きた時刻における推定された瞬間マグニチュード分布の関係を示す図である。横軸はマグニチュードであり、縦軸は確率密度を表す。黒丸は、観測された各地震のマグニチュードの正規化ヒストグラムを示し、灰色の曲線は、各地震が起きた瞬間のマグニチュード分布を示し、黒い曲線は、その平均曲線である。このように、各地震が起きた瞬間の推定値も実測値とよく一致していることが確かめられた。
【００８５】
次に推定された検出率を用いて、予測モデル（式（１））のパラメータの推定を行い、予測実験を行った。図１０および図１１は、本震後１日までの学習期間で推定されたパラメータを用いて、１日後から３０日までの期間、マグニチュード５以上、および、６以上の余震の頻度を、それぞれ予測した結果を示す図である。実線は、本実施の形態により推定された検出率関数とグーテンベルグ・リヒター則の積に基づく関数をもとに、大森・宇津則をフィッティングさせ、そこから予測された余震の頻度を示し、破線は、検出率を考慮せず大森・宇津則を直接フィッティングさせた曲線を示す。
【００８６】
また、図１２および図１３は、更に短い学習期間において推定されたパラメータを用いて、その後の期間、マグニチュード５以上、および、６以上の余震の頻度を、それぞれ予測した結果を示す図である。一点線は学習期間３時間、破線は６時間、二点鎖線は１２時間、実線は２４時間の学習期間で得られた予測結果を示す。その結果、本実施の形態による予測と、実際の観測がよく一致していることが確かめられた。
【００８７】
［他の実施の形態］
さて、これまで本発明の実施の形態について説明したが、本発明は、上述した実施の形態以外にも、特許請求の範囲に記載した技術的思想の範囲内において種々の異なる実施の形態にて実施されてよいものである。
【００８８】
例えば、余震予測装置１００がスタンドアローンの形態で処理を行う場合を一例に説明したが、余震予測装置１００は、クライアント端末からの要求に応じて処理を行い、その処理結果を当該クライアント端末に返却するようにしてもよい。
【００８９】
また、実施の形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部または一部を手動的に行うこともでき、あるいは、手動的に行われるものとして説明した処理の全部または一部を公知の方法で自動的に行うこともできる。
【００９０】
このほか、上記文献中や図面中で示した処理手順、制御手順、具体的名称、各処理の登録データや検索条件等のパラメータを含む情報、画面例、データベース構成については、特記する場合を除いて任意に変更することができる。
【００９１】
また、余震予測装置１００に関して、図示の各構成要素は機能概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。
【００９２】
例えば、余震予測装置１００の各装置が備える処理機能、特に制御部１０２にて行われる各処理機能については、その全部または任意の一部を、ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）および当該ＣＰＵにて解釈実行されるプログラムにて実現してもよく、また、ワイヤードロジックによるハードウェアとして実現してもよい。なお、プログラムは、後述する記録媒体に記録されており、必要に応じて余震予測装置１００に機械的に読み取られる。すなわち、ＲＯＭまたはＨＤなどの記憶部１０６などには、ＯＳ（Ｏｐｅｒａｔｉｎｇ  Ｓｙｓｔｅｍ）として協働してＣＰＵに命令を与え、各種処理を行うためのコンピュータプログラムが記録されている。このコンピュータプログラムは、ＲＡＭにロードされることによって実行され、ＣＰＵと協働して制御部を構成する。
【００９３】
また、このコンピュータプログラムは、余震予測装置１００に対して任意のネットワーク３００を介して接続されたアプリケーションプログラムサーバに記憶されていてもよく、必要に応じてその全部または一部をダウンロードすることも可能である。
【００９４】
また、本発明に係るプログラムを、コンピュータ読み取り可能な記録媒体に格納してもよく、また、プログラム製品として構成することもできる。ここで、この「記録媒体」とは、メモリーカード、ＵＳＢメモリ、ＳＤカード、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＥＰＲＯＭ、ＥＥＰＲＯＭ、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤ、および、Ｂｌｕ−ｒａｙ（登録商標）  Ｄｉｓｃ等の任意の「可搬用の物理媒体」を含むものとする。
【００９５】
また、「プログラム」とは、任意の言語や記述方法にて記述されたデータ処理方法であり、ソースコードやバイナリコード等の形式を問わない。なお、「プログラム」は必ずしも単一的に構成されるものに限られず、複数のモジュールやライブラリとして分散構成されるものや、ＯＳ（Ｏｐｅｒａｔｉｎｇ  Ｓｙｓｔｅｍ）に代表される別個のプログラムと協働してその機能を達成するものをも含む。なお、実施の形態に示した各装置において記録媒体を読み取るための具体的な構成、読み取り手順、あるいは、読み取り後のインストール手順等については、周知の構成や手順を用いることができる。
【００９６】
記憶部１０６に格納される各種のデータベース等（リアルタイム観測ファイル１０６ａ、パラメータファイル１０６ｂ、予測結果ファイル１０６ｃ等）は、ＲＡＭ、ＲＯＭ等のメモリ装置、ハードディスク等の固定ディスク装置、フレキシブルディスク、および、光ディスク等のストレージ手段であり、各種処理やウェブサイト提供に用いる各種のプログラム、テーブル、データベース、および、ウェブページ用ファイル等を格納する。
【００９７】
また、余震予測装置１００は、既知のパーソナルコンピュータ、ワークステーション等の情報処理装置として構成してもよく、また、該情報処理装置に任意の周辺装置を接続して構成してもよい。また、余震予測装置１００は、該情報処理装置に本発明の方法を実現させるソフトウェア（プログラム、データ等を含む）を実装することにより実現してもよい。
【００９８】
更に、装置の分散・統合の具体的形態は図示するものに限られず、その全部または一部を、各種の付加等に応じて、または、機能負荷に応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。すなわち、上述した実施形態を任意に組み合わせて実施してもよく、実施形態を選択的に実施してもよい。
【産業上の利用可能性】
【００９９】
以上詳述に説明したように、本発明によれば、本震直後の限られた地震観測データであってもリアルタイムに余震確率を推定することができる、余震予測装置、余震予測方法、および、プログラムを提供することができ、気象庁等による余震予測の公表や、民間企業による災害情報の予測を有料で配信するサービスなど、防災や地震予知などの様々な分野において極めて有用である。
【符号の説明】
【０１００】
１００  余震予測装置
１０２  制御部
１０２ａ  パラメータ推定部
１０２ｂ  余震発生率予測部
１０４  通信制御インターフェース部
１０６  記憶部
１０６ａ  リアルタイム観測ファイル
１０６ｂ  パラメータファイル
１０６ｃ  予測結果ファイル
１０８  入出力制御インターフェース部
１１２  入力装置
１１４  出力装置
２００  外部システム
３００  ネットワーク
【図１】
2014219315.tif 000030
【図２】
2014219315.tif 000031
【図３】
2014219315.tif 000032
【図４】
2014219315.tif 000033
【図７】
2014219315.tif 000034
【図１０】
2014219315.tif 000035
【図１１】
2014219315.tif 000036
【図１２】
2014219315.tif 000037
【図１３】
2014219315.tif 000038
【図５】
2014219315.tif 000039
【図６】
2014219315.tif 000040
【図８】
2014219315.tif 000041
【図９】
2014219315.tif 000042