特開 2023-82905 地震波解析装置、方法及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2023082905

(43)【公開日】2023-06-15

(54)【発明の名称】地震波解析装置、方法及びプログラム

(51)【国際特許分類】

   G01V 1/28 20060101AFI20230608BHJP

   G01V 1/00 20060101ALI20230608BHJP

【ＦＩ】

G01V1/28

G01V1/00 D

【審査請求】未請求

【請求項の数】6

【出願形態】ＯＬ

(21)【出願番号】P 2021196907

(22)【出願日】2021-12-03

【新規性喪失の例外の表示】特許法第３０条第２項適用申請有り ・一般社団法人 日本音響学会 ２０２１年 春季研究発表会 講演論文集（令和３年２月２４日発行） ・特定非営利活動法人 海洋音響学会 ２０２１年度 第３７回研究発表会 講演予稿集 ｈｔｔｐｓ：／／ｗｗｗ．ｍａｓｊ．ｊｐ／ｅｖｅｎｔｓ＿ｒｅｃｒｕｉｔ／ｒｐ－２０２１０５２７＿ｇｅｎｅｒａｌ／ｒｐ－２０２１０５２７／（令和３年５月２７日掲載） ・一般社団法人 日本音響学会 ２０２１年 秋季研究発表会 講演論文集（令和３年８月２４日発行） ・特定非営利活動法人 超音波エレクトロニクス協会 ２０２１年 第４２回Ｓｙｍｐｏｓｉｕｍ ｏｎ ＵｌｔｒａＳｏｎｉｃ Ｅｌｅｃｔｒｏｎｉｃｓ（ＵＳＥ２０２１）（超音波エレクトロニクスの基礎と応用に関するシンポジウム） 講演予稿集 ｈｔｔｐ：／／ｗｗｗ．ｓｅｎｋｙｏ．ｃｏ．ｊｐ／ｕｓｅ／ＵＳＥ２０２１／ｉｎｄｅｘ．ｈｔｍｌ（令和３年１０月２５日掲載）

(71)【出願人】

【識別番号】505363178

【氏名又は名称】菊池 年晃

(74)【代理人】

【識別番号】100080090

【弁理士】

【氏名又は名称】岩堀 邦男

(72)【発明者】

【氏名】菊池 年晃

【テーマコード（参考）】

2G105

【Ｆターム（参考）】

2G105AA03

2G105BB01

2G105DD01

2G105EE02

2G105GG04

2G105MM01

2G105MM03

(57)【要約】

【課題】地震予知に有益な情報を提供する。
【解決手段】地震波解析装置１０は、ある観測点で実測された地震データを逐次取得する地震データ取得手段１１と、地震データ取得手段１１で取得された地震データに対して時間周波数解析を行うことにより、地震波の最低周波数を算出する最低周波数算出手段１２と、最低周波数算出手段１２で算出された最低周波数の経時変化に基づき、最低周波数の上昇を検出する周波数上昇検出手段１３と、を備えたものである。
【選択図】図１

【特許請求の範囲】

【請求項1】

ある観測点で実測された地震データを逐次取得する地震データ取得手段と、
この地震データ取得手段で取得された前記地震データに対して時間周波数解析を行うことにより、地震波の最低周波数を算出する最低周波数算出手段と、
この最低周波数算出手段で算出された前記最低周波数の経時変化に基づき、前記最低周波数の上昇を検出する周波数上昇検出手段と、
を備えた地震波解析装置。

【請求項2】

複数の観測点で一定期間実測された本震及び本震前の地震データを取得する地震データ取得手段と、
前記複数の観測点で実測された本震の地震データに特異値分解法を適用して、震源地における地震波の放射軸方向を算出する放射軸方向算出手段と、
前記放射軸方向に位置する前記観測点における前記本震前の地震データに対して時間周波数解析を適用して、前記本震前の地震波の最低周波数を算出する最低周波数算出手段と、
前記最低周波数の経時変化と前記本震との関係を検出する前兆検出手段と、
を備えた地震波解析装置。

【請求項3】

ある観測点で実測された地震データを逐次取得し、
取得された前記地震データに対して時間周波数解析を行うことにより、地震波の最低周波数を算出し、
算出された前記最低周波数の経時変化に基づき前記最低周波数の上昇を検出する、
地震波解析方法。

【請求項4】

複数の観測点で一定期間実測された本震及び本震前の地震データを取得し、
前記複数の観測点で実測された本震の地震データに特異値分解法を適用して、震源地における地震波の放射軸方向を算出し、
前記放射軸方向に位置する前記観測点における前記本震前の地震データに対して時間周波数解析を適用して、前記本震前の地震波の最低周波数を算出し、
前記最低周波数の経時変化と前記本震との関係を検出する、
地震波解析方法。

【請求項5】

ある観測点で実測された地震データを逐次取得する地震データ取得手段、
この地震データ取得手段で取得された前記地震データに対して時間周波数解析を行うことにより、地震波の最低周波数を算出する最低周波数算出手段、及び、
この最低周波数算出手段で算出された前記最低周波数の経時変化に基づき、前記最低周波数の上昇を検出する周波数上昇検出手段、
としてコンピュータを機能させるための地震波解析プログラム。

【請求項6】

複数の観測点で一定期間実測された本震及び本震前の地震データを取得する地震データ取得手段、
前記複数の観測点で実測された本震の地震データに特異値分解法を適用して、震源地における地震波の放射軸方向を算出する放射軸方向算出手段、
前記放射軸方向に位置する前記観測点における前記本震前の地震データに対して時間周波数解析を適用して、前記本震前の地震波の最低周波数を算出する最低周波数算出手段、及び、
前記最低周波数の経時変化と前記本震との関係を検出する前兆検出手段、
としてコンピュータを機能させるための地震波解析プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、微弱な地震による地震データに基づいて地震波を解析し、これにより地震前兆を観測する地震波解析装置、方法及びプログラムに関する。

【背景技術】

【0002】

地震に関する研究は様々な観点から進められている（特許文献１など）。例えば、地震の震源構造を解明するために、地震波を雑音性信号と捉えて、幾つかの信号処理法を駆使して、震源構造を解明することが考えられる。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特許第５６１６３９２号公報

【非特許文献】

【0004】

【非特許文献1】T.Kikuchi K.Mizutani USE2019 Vibration Structure and radiation waves of active fault

【非特許文献2】菊池、水谷 海洋音響学会 2020.5 活断層の周辺環境によるラテラル波の変動

【非特許文献3】馬杉正男 信号解析 森北出版 P81

【非特許文献4】菊池、水谷 音響学会 2018.5 特異値分解法による震源振動の放射特性

【非特許文献5】菊池、水谷 音響学会 2020.9 活断層振動の時間周波数解析とその方位性

【非特許文献6】Nai-chyuan Yen et al. JASA,87,2359(1990)

【非特許文献7】陳 延偉 : “独立成分分析法(ICA)のパターン認識・画像処理への応用とMATLABシミュレーション”、株式会社トリケップス

【発明の概要】

【発明が解決しようとする課題】

【0005】

しかしながら、震源内部で発生する地震固有の振動と伝搬環境に内在する雑音成分とは、周波数的に差異がないため容易には分離できない。そこで、本発明では、特異値分解法、時間周波数解析法及び独立成分分析法を用いて、震源の全体構造を把握した上で、その内部振動構造を解明することにした。その結果、防災上重要な結果が得られた。

【課題を解決するための手段】

【0006】

本発明に係る第一の地震波解析装置は、
ある観測点で実測された地震データを逐次取得する地震データ取得手段と、
この地震データ取得手段で取得された前記地震データに対して時間周波数解析を行うことにより、地震波の最低周波数を算出する最低周波数算出手段と、
この最低周波数算出手段で算出された前記最低周波数の経時変化に基づき、前記最低周波数の上昇を検出する周波数上昇検出手段と、
を備えたものである。

【0007】

本発明に係る第二の地震波解析装置は、
複数の観測点で一定期間実測された本震及び本震前の地震データを取得する地震データ取得手段と、
前記複数の観測点で実測された本震の地震データに特異値分解法を適用して、震源地における地震波の放射軸方向を算出する放射軸方向算出手段と、
前記放射軸方向に位置する前記観測点における前記本震前の地震データに対して時間周波数解析を適用して、前記本震前の地震波の最低周波数を算出する最低周波数算出手段と、
前記最低周波数の経時変化と前記本震との関係を検出する前兆検出手段と、
を備えたものである。

【0008】

本発明に係る地震波解析方法は、本発明に係る第一及び第二の地震波解析装置の動作に対応するものである。本発明に係る地震波解析プログラムは、本発明に係る第一及び第二の地震波解析装置の各手段をコンピュータに機能させるためのものである。

【発明の効果】

【0009】

本発明によれば、実測した地震データに基づいて活断層の存在、その方位性に加えて、活断層の発達を監視することができ、地震の前兆を正確に監視することができる。

【図面の簡単な説明】

【0010】

【図1】本発明の実施形態に係る地震波解析装置の構成を示すブロック図であり、［Ａ］は実施形態１、［Ｂ］は実施形態２である。

【図2】本発明の実施形態に係る地震波解析方法を説明するフローチャートであり、［Ａ］は実施形態１、［Ｂ］は実施形態２である。

【図3】２０１２年１月２８日に富士山付近で発生したＭ５．４の地震（以下「対象地震」という。）に関し、震源地から各観測点を見た方位角と各観測点における主要成分周波数との関係を示すグラフである。

【図4】図３に示すデータに各観測点の地域性を加味して補正したグラフである。

【図5】図４に示す方位角２８８°の駒ケ根における受信波のＷＤＦ分布を示すグラフである。

【図6】図４に示す方位角２８７°の芦安における受信波のＷＤＦ分布を示すグラフである。

【図7】図４に示す方位角９°の皆野における受信波のＷＤＦ分布を示すグラフである。

【図8】図４に示す方位角６７°の府中における受信波のＷＤＦ分布を示すグラフである。

【図9】独立成分分析を示す概念図である。

【図10】図４及び図１１に示す駒ケ根における受信波のＩＣＡ処理結果を示すグラフであり、上図が第一分離成分、中図が第二分離成分、下図が第三分離成分である。

【図11】活断層の軸線上に沿った観測点を示す概念図である。

【図12】図４及び図１１に示す上松における受信波のＩＣＡ処理結果を示すグラフであり、上図が第一分離成分、中図が第二分離成分、下図が第三分離成分である。

【図13】図４及び図１１に示す下呂北における受信波のＩＣＡ処理結果を示すグラフであり、上図が第一分離成分、中図が第二分離成分、下図が第三分離成分である。

【図14】活断層のダクトモデルを示す概念図である。

【図15】図１４に示すダクトから放射されたパルスを示すグラフである。

【図16】対象地震に関し、前兆地震の発生を示すグラフであり、縦軸は上図が緯度、中図が緯度、下図がマグニチュードであり、各図の横軸は観測開始からの経過日数である。

【図17】対象地震に関し、前兆地震から本震に至るまでのマグニチュード及びＷＤＦの最低周波数の変化を示すグラフであり、横軸は観測開始からの経過日数である。

【図18】活断層周辺の音速変化による放射パルスの変化を示すグラフである。

【図19】活断層周辺の音速変化による放射パルスのスペクトル変化を示すグラフである。

【発明を実施するための形態】

【0011】

＜本発明の基礎となる知見＞
本発明者は、地震予知に関する研究を進めた結果、以下のような知見を得た。この知見を図面に沿って説明する。

【0012】

１．地震波の放射構造
1.1地震波の起源
地震波の起源は活断層内部に生じた亀裂による短いパルスであり、そのパルスは活断層特有の狭い層を伝搬した後に放射される。その際、狭い層の両側面内に生じた二つのラテラル波が干渉しながら伝搬する。その伝搬の際に、ラテラル波の干渉によって衝撃波と周期的パルスが生成される。また、活断層の長さは周期的波の波長に比して非常に大きいため、活断層から放射される波の指向性は狭く複雑になる。詳細は非特許文献１、２を参照。

【0013】

1.2特異値分解法
特異値分解とは、任意のｍ×ｎの行列Ｘが与えられたときに、Ｕ^TＵ＝Ｖ^TＶ＝Ｅ（Ｅ：単位行列)を満たす正規直交行列Ｕ（ｍ×ｒ）、Ｖ（ｎ×ｒ）を用いて分解するための手段であり、次の関係が成立する。

【数1】

…(1)
ここで、

【数2】

は対角成分以外が０となる対角行列である。対角成分

【数3】

は特異値と呼ばれる。(1)式に関して、

【数4】

…(2)
とおくと、

【数5】

…(3)
と表現される。これは行列のスペクトル分解と呼ばれる。
観測で得られた地震波Ｘを成分の大きさ順に級数展開し、その第一項を主要成分とする。詳細は非特許文献３、４を参照。

【0014】

1.3地震波の方位角分布
対象地震に対して検証し、その結果を示す。解析には上記の特異値分解法を用いる。求めた放射波の方位角分布を図４に示す。

【0015】

図４から明らかなように、西側と東側で異なる分布が見られる。全体的にみると双極性放射特性である。図４において、西側と東側の境界線（２０°から２１０°）は、気象庁が求めた震源球の節面と一致する。この節面は活断層の方位を表すと言われているが、通常の音源の指向性からかけ離れている。したがって、実際の活断層の方位は、もう一つの節面、すなわち３００°方向と考えられる。詳細は非特許文献５参照。

【0016】

更に詳しく説明する。地殻の亀裂によって生じた活断層の振動によって、地震波が放射される。しかしながら、活断層の振動構造はまだ十分に解析されていない。そこで、活断層の振動モデルの理論解析と実際の地震とに基づいて、放射波を解析する。ここでは地震データの解析法について述べる。前述したように、対象地震の震源を囲む観測点で得られたＰ波に、特異値分解法を適用して周波数成分を求めた。

【0017】

その主要成分の周波数を震源からの方位角に対して表した結果を、図３に示す。方位によって放射形態が大きく異なっていることが分かる。そして、震源から観測点までの距離は、数キロから数百キロまで大きく分散している。したがって、図３のデータには、地殻の地域的変化や拡散の影響が包含されていると考えられる。そこで、それらの影響を低減するためにタイムリバーサル法を適用する（特許文献１参照）。つまり、観測点での受信波にタイムリバーサル処理を施して、震源付近に置換した波を求める。得られた波に特異値分解法を適用した結果を、図４に示す。図４では、図３と比較すると対称性が顕著である。そして、方位角２１０°を境として放射パターンが大きく異なっている。すなわち、方位角２１０°以上では方位角に対する放射周波数の変動が急峻（方位角が少し変わっただけで周波数が大きく変動する）であるが、方位角２１０°以下では放射周波数の変動は緩かである。全体的にみると、変動が急峻な部分と緩やかな部分との双極性放射である。

【0018】

その他の地震に関しても同様な処理を行った結果、双極性放射であることが示された。更に、前兆地震の放射特性も双極性放射であることも明らかにした。それは、活断層の方位が本震と前兆地震とで同一であることを示している。その方位同一性を利用して、それらの時系列変化を調べることにより予知情報を得ることができる。

【0019】

２．方位角成分の分析
2.1時間周波数解析
複雑な信号は、その時間密度や周波数密度のフーリエ変換だけでは、その性質を完全には表すことができない。それを同時に解析する方法として時間周波数解析が用いられる。時間周波数解析の中で最も基本的な方式はウイグナー分布関数（Wigner distribution function：WDF）で、そのデジタル表示は次式で表される。詳細は非特許文献６を参照。

【数6】

…(4)
ここで、

【数7】

また

【数8】

は対象とする関数である。詳細は非特許文献５、６を参照。

【0020】

2.2方位角と時間周波数解析
図４における方位角分布で、３００°方向に節面、すなわち活断層の方位があることが指摘された。それらを確認するために、ここでは地震波に対してウイグナー分布関数を適用する。適用した結果を図５、図６、図７及び図８に示す。これらの図において、横軸は周波数、縦軸は時間である。図５から明らかなように、四つの周波数成分があり、それらが周波数を変えずに時間的に変化している。すなわち、それらは互いに独立である。これにより、図５の観測点は、活断層の真正面に位置することが分かる。一方、図７では、８秒付近に全ての周波数成分が重なっている。これは、多数の成分が重合して観測される、活断層の真横方向に位置することを表している。この観測点である皆野は図４の節面方向と一致している。

【0021】

更に詳しく説明する。１節で活断層放射の双極性放射が明らかになった。この双極性は５節による理論解析の結果、活断層の振動方位に関することが示された。そこで、この項では活断層の振動方位を求める。その解析法として時間周波数解析法を用いる。観測点における受信波にＷＤＦを適用して、時間周波数解析を行う。図４の双極性の節に近い皆野（方位角９°、距離６５．６ｋｍ）に対する結果を、図７に示す。時間８秒付近に、周波数０．５から１０Ｈｚまでのほぼ全周波数成分が分布している。ここでは、活断層の長さに比して観測点までの距離が非常に大きいため、活断層上の全ての点における振動が同時刻に到達する。これは、パルスの進行方向に直角の位置であることを表している。

【0022】

パルスの進行方向は活断層の亀裂の方向と一致することが知られている。その延長方向に位置する観測点は、駒ケ根（方位角２８８．１°、距離９５．９ｋｍ）である。その受信波を解析した結果を、図５に示す。主な周波数成分は、１．５、３．５、５．９、７．８及び１０．５Ｈｚである。そして、それら振動成分はそれらの周波数を変えずに時間的に継続している。各成分とも他成分との変調はなく一定の周波数を維持しており、いわゆる音源の本来の振動を表す縦型パターンである。駒ケ根を中心とした方位角２１０°から１０°の範囲では、周波数は異なるがほぼ縦型パターンである。活断層の延長線上で駒ケ根の正反対方向に位置する横須賀（方位角１１５．０°、距離７２．３ｋｍ）でも、縦型パターンが見られたが周波数成分は０．５Ｈｚのみである。また、図５ではパルスの先頭部分のレベルが高いが、横須賀ではパルスの中間から後半でレベルが高くなる。すなわち、駒ケ根方向に進行したパルスが反射して横須賀方向へ進行するが、その間に高周波成分は減衰したと考えられる。一方、方位角が１０°から２１０°では散乱パターンになる。例えば八王子（方位角５２．２°、距離３２．９ｋｍ）の散乱パターンは、活断層内の進行波と反射波が混合した状態を表している。すなわち、これらの散乱パターンは双極性特性の緩やかな分布が形成される主因である。

【0023】

以上のとおり、振動形態すなわち時間周波数分布から、活断層の主軸の放射の方位角が求められた。この方位角に位置する観測点（ここでは駒ケ根）における時間周波数解析を前兆地震に適用すると、地震予知に関するデータが得られる。地震波のビームは非常に狭いので、主軸のビームも狭く観測方位は限られる。

【0024】

３．活断層内部振動の抽出
２節で活断層の方位角が明らかになった。しかしながら、図５が示すように、各周波数分布はある幅を持っている、すなわち、雑音成分も含まれている。そこで、活断層内部の固有の振動を取得するために、次に述べる独立成分分析法を用いる。

【0025】

3.1独立成分分析法
独立成分分析法（Independence Component Analysis：ICA）は、新しい多次元信号解析法として、１９８０年代に提案され、１９９０年代に入ってから、その理論と応用が盛んに研究されるようになっている。ＩＣＡは、互いに独立な信号が重なり合った混合信号を幾つかの異なる条件で観測し、それを基に原信号（独立成分）を分離する問題として、定式化される。詳細は非特許文献７を参照。

【0026】

ＩＣＡの基本概念図を図９に示す。観測信号は、幾つかの統計的独立な原信号の線形和からなる。原信号と混合プロセスが未知で、観測信号のみで原信号を分離・推定する。
図９に示すように、ｎ個の未知原信号を次のベクトルで与えられたとする。

【数9】

ｔ＝１，２，…，Ｎ …(5)
各信号は統計的に互いに独立で、平均が０であるとする。これらの信号が線形的に混合され、ｍ個のセンサーで観測した信号を、

【数10】

ｔ＝１，２，…，Ｎ …(6)
とする。ｓ（ｔ）とｘ（ｔ）との間に次のような線形関係が成り立つ。
ｓ（ｔ）＝Ａｘ（ｔ） …(7)
ここで、Ａは混合行列で、ｍ×ｎのfull column rank行列（ｍ≧ｎ）であるとする。係数ａ_ijは、原信号ｓ_iが観測信号ｘ_jへの寄与又は度合いを表す。混合行列Ａが既知の場合、原信号ｓ（ｔ）の最小二乗解は、Ａのpseudo inverse Ａ^＋を用いて、Ａ^＋ｘ（ｔ）と簡単に求めることができるが、一般の場合Ａも未知なので観測信号ｘのみから原信号ｓ又は混合行列Ａを推定しなければならない。すなわち、観測信号ｘのみから、あるｍ×ｎの分離行列Ｗを用いて、

【数11】

…(8)
なる分離行列ｙを計算する。ＩＣＡでは、ｙの各成分が互いに統計的に独立となるＷを求めることが目的である。理想的には、

【数12】

…(9)
となれば、ｙとｓは一致する。しかし、ｙの各成分が互いに順序及び大きさは独立性に影響しないことから、

【数13】

…(10)
となる分離行列を求めることになる。ここで、Ｄは対角行列であり、Ｐは置換行列（各行各列に１となる要素が１つだけある行列）である。結局、(10)式を最適化問題として解くことになり、いろいろな方法が提案されている。

【0027】

3.2軸方向におけるＩＣＡ解
1.3項と同様に、対象地震に対して検証する。通常、地震計は東西動、南北動及び上下動の三成分を同時に計測・記録している。観測点の駒ヶ根において受波した信号の三成分を原信号として、ＩＣＡ処理した結果を、図１０に示す。上図の第一分離成分は、通常観測される地震波形に類似している。中図の第二分離成分は、先頭部分に周期的成分が形成されている。下図の第三分離成分は、しばしば観測される通常の衝撃波に類似している。

【0028】

２節で示したように、地震波の放射波形は、方位角に対する変化が鋭敏（狭ビーム放射）である。この性質は、方位角すなわち水平方向の変化のみならず、垂直方向の変化も鋭敏と考えられる。大多数の観測点は地上にあるので、垂直方向の変化は距離方向の変化として現れる。駒ヶ根に近い方位角、すなわち、活断層の軸上線で、距離を変えた点（図１１）におけるＩＣＡ結果を求めた。

【0029】

駒ヶ根（図１０）、上松（図１２）及び下呂北（図１３）の各図の第一分離波の波形を見ると、明らかに上松の波形が他の波形とかけ離れている。上松の波形は、正弦的周期性を持ったパルスが継続している。この波形は、活断層をモデル化した薄い層内を伝搬したラテラル波の放射に一致している。したがって、図１２の第一分離成分は活断層内部の固有振動を表していることが明らかである。

【0030】

なお、ダクトモデルと伝搬を以下に要約する。図１４に示すダクトモデルで、幅Ｗを３５０ｍ、底部速度Ｃ_Bを３０００ｍ／ｓ、ダクトの長さを１１ｋｍとする。このとき、ダクトの基点から単周期パルスを放射し、ダクトの軸上の距離４０ｋｍで受けた波を図１５に示す。図１５において、Ａはダクト本来の放射波（先頭衝撃波）、Ｂは地表面からの反射波（地表反射波）、Ｃはダクト端面で反射した波（活断層内の反射波）である。Ａは、単パルスに続くＶ型ダクトの特有の指数関数減衰型である。Ｂは、ダクト内を通らないことからダクトの分散性を受けないので、発振パルスと同形である。Ｃは、ダクトの基点と先端とを往復したダクト内の反射波である。

【0031】

これらのラテラルパルスと図１２の第一分離波とを比較すると、先端の周期的パルス及び後方の周期的パルスの形成が両者で一致していることが明らかである。以上述べたように、活断層の内部の震動構造と放射構造を一体化して明らかにすることができた。

【0032】

活断層の主軸ビームの観測方位角が狭いので、観測点は限られる。しかし、５節に示す理論的解析によると、活断層内部の振動は活断層の幅によって決まる特有の振動性パルスを形成する。それを観測によって検出するために、独立成分分析法を適用する。

【0033】

４．活断層震動の経時変化
３節までで、活断層の内部の振動構造と放射構造を一体化して説明した。ここではそれらの応用について述べる。

【0034】

一般に、大きな地震の前後には小さな地震が多発することが多くみられる。今回解析した対象地震を本震と考えて、それ以前に発生した地震のマグニチュードと位置との関係を図１６に示す。約３年間にわたる発生位置は、ほぼ本震を含む限られた層内で発生したことが分かる。これらの地震に対して特異値分解法を適用して放射構造を求めたところ、本震（図４）に近いことが明らかになった。それに伴って、放射軸方向も本震の場合とほぼ同一方向であることが明らかになった。

【0035】

そこで、前兆地震の軸方向においてウイグナー分布解析を行い、その最低周波数を求めた。その結果を図１７に示す。マグニチュード（破線）は、ばらつきはあるが平均的には変化が見られない。一方、最低周波数（実線）の下部は、４５０日の１．１Ｈｚから１０５０日の２．５Ｈｚまで上昇している。これは、前兆現象の一つと考えられる。なお、地殻歪が増大すると、同時に周辺圧力が上昇する。圧力が上昇すると、ラテラル波の速度が上昇して、周波数が上昇することは、別途示す。

【0036】

更に詳しく説明する。３節に示したように、時間周波数解析法と独立成分分析法を併用することにより、活断層の変動を監視して、それらの急激な変動から地震の予知情報を得ることができる。図１７において、横軸は２０１２年１月２８日以前に発生した前兆地震の日々であり、破線はマグニチュード、実線はＷＤＦの最低周波数である。マグニチュードはばらつきがあるが平均的には変化があるとは認められない。一方、最低周波数は本震の５００日前付近から上昇が見られる。この現象をシミュレーションによって示したのが図１９である。

【0037】

５．活断層モデルの放射
４節までは、観測結果から得られた地震振動現象の特性を解析している。ここでは、観測結果を裏付けるために理論シミュレーションによって振動現象を解析する。

【0038】

一般に、活断層周辺の地殻の歪が増大すると、その媒質の音速が増大する。その音速が増大するとラテラル波の速度が上昇する。それによって活断層から放射されるパルスの周波数が増大すると考えられる。図１９は図１５のＣパルスについてのシミュレーション結果である。図１９において、最上の曲線は、周辺の音速が６１００ｍ／ｓのときの放射パルスのスペクトルである。その下の曲線は周辺の音速を１００ｍ／ｓずつ増加して、最下の曲線は周辺の音速が７０００ｍ／ｓの時の放射パルスのスペクトルである。周辺の音速の増大に伴って、放射パルスの周波数が上昇している。

【0039】

更に詳しく説明する。活断層の歪が増大すると圧力も増大する。その圧力の増大に伴って活断層のラテラル波の速度が増大する。固体の場合、圧力が上がると硬さが増し、剛性率が上がるために音速が増すからである。その結果、活断層から放射されるパルスが速くなりスペクトルが変化する。したがって、それらの変化により活断層の変化を認識できる。これらのメカニズムをシミュレーションによって確認する。

【0040】

活断層の放射機構の概念図を図1４に示す。左に速度の深度分布を示す。ダクト内部の速度構造は幾つか考えられるがここではβの様な分布とする。各部の値は前述したように以下のとおりである。
音源の深度 ２０ｋｍ
活断層の長さ ２０ｋｍ
活断層の幅 ４４０ｍ
起振源 インパルス
活断層の中心の速度 Ｃo：３０００ｍ／ｓ

【0041】

活断層から水平に４０ｋｍの点で受波したパルスを図１５に示す。パルスＣは活断層内を１．５往復した波であるから活断層の影響を多く受けている。このパルスＣが周辺の音速Ｃoによって如何に変化するかを検証する。図１８は、周囲音速変化によるパルスの変化を示すシミュレーション結果である。図１８の最上のパルスはＣoが６０００ｍ／ｓのときのパルスで、それから下のパルスは速度を順次１００ｍ／ｓ増加させた場合のパルスを示している。そして、図１９は図１８に示したパルスのパワースペクトルである。図１８及び図１９から、周辺の音速の増大に伴って、放射パルスの周波数が上昇していることがわかる。以上、外周の圧力変化よって活断層の放射パルスが変化することを確認した。

【0042】

６．まとめ
活断層の内部振動を検出する手法を確立した。この手法は、活断層による地震が危惧されている地域や原子力発電所の様な重要な施設において、活断層の活動状況を監視するために活用することができる。本手法では、受信波の主要成分のみを扱い、雑音成分を除去している。その結果、明確な双極性特性を得ることができた。この明確な双極性特性から活断層の主軸がほぼ算出できるが、更に時間周波数解析で主軸放射角を確定することができた。更に独立成分分析法を適用することにより、活断層の内部振動を検出できた。その結果、シミュレーションによる予測解析と実測データとの比較が可能になり、地震の現況監視から予知につなげることが可能になる。

【0043】

本手法の処理手順の一例を述べれば、次のとおりである。
１）震源の全範囲にわたる観測点で地震波（Ｐ波）を観測し、その観測波の最大特徴を特異値分解法で求め、それらの方位分布を求める。この方法は、震源球と違ってマグニチュードが１．５程度の微弱な地震にも適用できるので、多くの前兆地震のデータを取り込むことができる。
２）次に振動軸方向を求める。ここでは、地震波に時間周波数分析法を適用し、１）でのおおよその放射構造を更に詳細に分析する。この必要性は地震波のビームが極端に狭いために起きる。ＷＤＦの放射パターンが縦型配列になる方位角を詳細に求める。方位角は、水平面のみならず、垂直面内でも同様である。垂直面内の角度変化は地上の観測点に対しては距離変化として現れる。このように、ＷＤＦパターンを方位と距離で探索し軸上ビームが地上に達する点（パラメトリックスポット）を求める。
３）独立成分分析。求められたパラメトリックスポットにおいて地震波に独立成分分析を適用し、活断層の内部振動を求める。
４）時系列変化。前兆地震に対して、２）と３）で得られたデータを時系列的に調べることにより、この地域の歪が安定しているか、又は歪が増大しているかを調べ、過去の同様な規模の地震と対比して危険度を取得する。

【0044】

本発明は上述の知見に基づきなされたものである。以下、本発明の実施形態を図面に沿って説明する。

【0045】

＜実施形態１＞
図１［Ａ］に示す本実施形態１の地震波解析装置１０は、ある観測点で実測された地震データを逐次取得する地震データ取得手段１１と、地震データ取得手段１１で取得された地震データに対して時間周波数解析を行うことにより、地震波の最低周波数を算出する最低周波数算出手段１２と、最低周波数算出手段１２で算出された最低周波数の経時変化に基づき、最低周波数の上昇を検出する周波数上昇検出手段１３と、を備えたものである。

【0046】

この地震データは、例えば防災科学技術研究所が管理しているデータを利用できる。地震データ取得手段１１は、地震が起きる度に、防災科学技術研究所のＨｉ－ｎｅｔから地震データをダウンロードして取得すればよい。最低周波数算出手段１２は、観測点における地震データに例えばウイグナー分布解析を適用して、その最低周波数を求める。周波数上昇検出手段１３は、最低周波数算出手段１２で逐次得られる最低周波数の経時変化に基づき、その最低周波数が上昇しているか否かを判断する。最低周波数が上昇すれば、図１７に示す実測例のように、大きな地震が起こる可能性が高まる。

【0047】

その理由は次のとおりである。断層周辺の地殻の歪が増大すると、断層周辺に大きな圧力が加わって、媒質の音速が増加する。圧力が上がると硬さが増し、剛性率が上がるために音速が増すからである。そして、音速が増加すると最低周波数が増加する。音速が増加すると最低周波数が増加することは、図１９に示すシミュレーション結果から明らかである。

【0048】

地震波解析装置１０によれば、逐次得られる地震データに対して時間周波数解析を適用して地震波の最低周波数を算出し、その最低周波数の上昇を検出することにより、地震予知に有益な情報を提供できる。

【0049】

図２［Ａ］に示す本実施形態１の地震波解析方法は、例えば地震波解析装置１０の動作に対応する。すなわち、本実施形態１の地震波解析方法は、ある観測点で実測された地震データを逐次取得し（ステップＳ１１）、取得された地震データに対して時間周波数解析を行うことにより、地震波の最低周波数を算出し（ステップＳ１２）、算出された最低周波数の経時変化に基づき最低周波数の上昇を検出する（ステップＳ１３）。

【0050】

なお、本実施形態１の地震波解析装置１０はハードウェアとして構築したが、これに限られるものではなく、コンピュータに地震波解析プログラムを実行させることにより、ソフトウェア上に、地震データ取得手段１１、最低周波数算出手段１２及び周波数上昇検出手段１３を構築するようにしてもよい。

【0051】

＜実施形態２＞
図１［Ｂ］に示す本実施形態２の地震波解析装置２０は、複数の観測点で一定期間実測された本震及び本震前の地震データを取得する地震データ取得手段２１と、複数の観測点で実測された本震の地震データに特異値分解法を適用して、震源地における地震波の放射軸方向を算出する放射軸方向算出手段２２と、放射軸方向に位置する観測点における本震前の地震データに対して時間周波数解析を適用して、本震前の地震波の最低周波数を算出する最低周波数算出手段２３と、最低周波数の経時変化と本震との関係を検出する前兆検出手段２４と、を備えている。

【0052】

この地震データは、実施形態１と同様に防災科学技術研究所が管理しているデータを利用できる。地震データ取得手段２１は、防災科学技術研究所のＨｉ－ｎｅｔから地震データをダウンロードして取得すればよい。放射軸方向算出手段２２は、複数の観測点で実測された本震の地震データに特異値分解法を適用して、図４に示すように主要成分周波数の方位角分布を求め、この方位角分布に基づき、震源地における地震波の放射軸方向を求める。このとき、精度向上のため、タイムリバーサル法を適用してもよい。最低周波数算出手段２３は、放射軸方向に位置する観測点における本震前の地震データに例えばウイグナー分布解析を適用して、それらの最低周波数を求める。前兆検出手段２４は、最低周波数算出手段２３で得られる最低周波数の経時変化に基づき、例えばその最低周波数が上昇しているか否かを判断する。最低周波数が上昇していれば、図１７に示す実測例のように、本震前の地震波が本震の前兆であると判断できる。その理由は、実施形態１で説明したとおりである。

【0053】

地震波解析装置２０によれば、本震前の地震データに対して時間周波数解析を適用して地震波の最低周波数を算出し、最低周波数の経時変化と本震との関係を検出することにより、地震予知に有益な情報を提供できる。

【0054】

図２［Ｂ］に示す本実施形態２の地震波解析方法は、例えば地震波解析装置２０の動作に対応する。すなわち、本実施形態２の地震波解析方法は、複数の観測点で一定期間実測された本震及び本震前の地震データを取得し（ステップＳ２１）、複数の観測点で実測された本震の地震データに特異値分解法を適用して、震源地における地震波の放射軸方向を算出し（ステップＳ２２）、放射軸方向に位置する観測点における本震前の地震データに対して時間周波数解析を適用して、本震前の地震波の最低周波数を算出し（ステップＳ２３）、最低周波数の経時変化と本震との関係を検出する（ステップＳ２４）。

【0055】

なお、本実施形態２の地震波解析装置２０はハードウェアとして構築したが、これに限られるものではなく、コンピュータに地震波解析プログラムを実行させることにより、ソフトウェア上に、地震データ取得手段２１、放射軸方向算出手段２２、最低周波数算出手段２３及び前兆検出手段２４を構築するようにしてもよい。

【0056】

以上のように、実施形態２によれば、実測した地震データに基づいて活断層の存在、その方位性に加えて、活断層における亀裂の発達を監視することができ、地震の前兆を正確に監視することができる。

【0057】

＜その他＞
以上、上記実施形態を参照して本発明を説明したが、本発明は上記実施形態に限定されるものではない。本発明の構成や詳細については当業者が理解し得るさまざまな変更を加えることができ、そのように変更された技術も本発明に含まれる。

【産業上の利用可能性】

【0058】

本発明は、実測した地震データに基づいて活断層の存在、その方位性に加えて、活断層における亀裂の発達を監視することにより、断層の実体に即した地震の前兆を的確に把握でき、地震予知に貢献できるものである。

【符号の説明】

【0059】

１０ 地震波解析装置
１１ 地震データ取得手段
１２ 最低周波数算出手段
１３ 周波数上昇検出手段
２０ 地震波解析装置
２１ 地震データ取得手段
２２ 放射軸方向算出手段
２３ 最低周波数算出手段
２４ 前兆検出手段

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】