統計数理（2015）
第 63 巻 第 1 号 129–144
c 2015 統計数理研究所
特集「地震予測と統計モデル」［研究詳解］

地球潮汐と地震活動との相関を用いた
地震活動予測岩田 貴樹1,2（受付 2015 年 1 月 5 日；採択 2 月 9 日）

要

旨過去，多くの研究によって，地球潮汐に起因する周期的な応力変動と地震活動との相関が議論されている．本稿では，この種の相関を議論するためによく用いられている統計的手法と，それに関連した近年の知見について紹介する．特に，この種の相関の有無が，地殻の応力状態を知るための指標となり得ることを指摘し，地震予測との関連性について議論する．更に，相関を網羅的に調べることに向けた提言を行ない，この分野の今後について展望する．
キーワード：地震活動，地球潮汐，点過程解析，地震活動予測．

1.

はじめに

月や太陽が地球に及ぼす万有引力により地球の海面が周期的に上下することは身近な現象である．液体部分と同様に，地球の固体部分も周期的に変形し，応力／歪み変化が生じており，これらを「地球潮汐」と呼ぶ．
この地球潮汐に起因する周期的な応力／歪み変化に伴い，地震が起きる現象は，古くから研究者の興味を引き，多くの研究が知られている．相関を肯定的に示した研究例としては，ある
限られた地域の，主に小地震の活動を調べたもの（例えば，Klein, 1976; Rydelek et al., 1988; 岩田・中西, 1998; Wilcock, 2001, 2009; Tolstoy et al., 2002; Stroup et al., 2007; Chen et al., 2012）
や，全地球の中〜大地震を含んだ地震カタログを解析したもの（Tsuruoka et al., 1995; Tanaka etがある．また，大地震のような顕著なイ
al., 2002a; Cochran et al., 2004; Métivier et al., 2009）
ベントに伴う，いわゆる余震活動に着目した研究もある（例えば，Mohler, 1980; Souriau et al.,．さらに，通常の地震が起こ
1982; Iwata, 2002; Crockett et al., 2006; Zhang and Zhuang, 2011）
す波に比べ，卓越周波数が低い波を起こす「低周波微動」と呼ばれる特異なイベントと地球潮汐との相関に関する報告が，近年は多くある（例えば，Shelly et al., 2007; Rubinstein et al., 2008;．そNakata et al., 2008; Thomas et al., 2009, 2012; Gallego et al., 2013; Ide and Tanaka, 2014），Heaton（1982），Rydelek et al.（1992），Vidale（1998）のように，相の一方で，Knopoﬀ（1964）
関に対して否定的な結果を示した研究も存在する．
このように結果が分かれる理由として考えられるのは，地震活動と地球潮汐との相関は普遍的な現象ではなく，限られた期間・地域にのみ現れるという可能性である．詳しくは 2.1 節にや Tanaka et al.（2002b）は，1982 年のトンガ地震の発生て後述するが，例えば，鶴岡（1995）
12

常磐大学 コミュニティ振興学部：〒 310–8585 茨城県水戸市見和 1–430–1
統計数理研究所 客員：〒 190–8562 東京都立川市緑町 10–3

130

統計数理

第 63 巻

第 1 号 2015

直前にのみ，有意な相関が見られることを指摘している．同様のことは，Li and Xu（2012）やTanaka（2012）
が，それぞれ 2008 年四川大地震や，2013 年東北地方太平洋沖地震の直前に限って，その震源域において地震活動と地球潮汐との有意な相関が見られることを指摘している．
2.2 節で触れる Yin et al.（1995）
を始めとする彼らの一連の研究も，カリフォルニア・中国・日
本の関東地方などで起きた顕著なイベントの前に，その前兆として，潮汐による応力が増加する期間に地震が偏って起きることを示している．これらが示唆することは，相関の有無は，地殻の応力状態に依存しており，（大）地震が起こりやすいような臨界状態にある時にのみ相関が現れるという考えである．これに基づけば，あるデータセットに対しては肯定的な結果が現れ，別のデータセットに対しては否定的な結果が現れるということへの，一応の説明にはなる．
とはいえ，地震活動と地球潮汐の相関を調べた研究は，あくまでケーススタディのレベルに留まっており，網羅的に調べることで相関の普遍性を議論した研究はまだ未達成な段階である．
こういった現状を踏まえ，本稿においては，この種の研究で典型的に用いられてきた統計的手法について紹介し，これらを用いた代表的な研究についてレビューを行なう．そして，包括的な研究に向けて，考え得る今後の展開について提言する．本稿の主たる目的は，あくまで手法とそれを用いた研究成果の詳解であり，地震活動と地球潮汐との相関の有無そのものについて議論するものではないことに留意されたい．
なお，著者は数年前に，やはり周期的な応力／歪み変化と地震活動変化との相関に関する統計的解析に関する解説を，Web を介して統計地震学教育を促進するためのプロジェクトである
the Community Online Resource for Statistical Seismicity Analysis（CORSSA）に寄稿している（Iwata, 2012）．本稿の一部は，その内容と重なっていることを，あらかじめお断りしておく．
2.

代表的な手法とその応用例

2.1 Schuster’s test
地震活動の周期性を調べるにあたって，古くからあり，かつ広く知られている方法が Schuster’s
であり，この方法は，地震活動と潮汐との相関を調べることにも適用出来
test（Schuster, 1897）
る．この方法を用いるにあたっては，まず，解析したい地震発生時系列のデータセットにおけ（位相）θi へ変換する．変換の際は，地球潮汐による応力/歪
る i 番目の地震の発生時刻を，角度み変化が極小となる隣り合った時刻を −180◦ および 180◦ ，これらに挟まれた応力/歪み変化が極大となる時刻を 0◦ とし，その間の時刻は線形に位相に割り振る．こうして得られた位相 θi のセットに対し，次の統計量 R を計算する．
2  N
2
 N
 


(2.1)
R=
cos θ
+
sin θ
.
i

i=1

i

i=1

ここで，N はデータセットに含まれている地震の個数である．幾何学的には，2 次元平面上において，原点を出発点として，方向（偏角）を θi とする長さ 1 のベクトル N 本を順に繋いでいった際の，終点の原点からの距離に R は相当する．また，統計科学においては，R2 /N は点過程
スペクトルの推定量として知られているものである（Bartlett, 1963）．
さて，もし，地震活動と地球潮汐の間に相関がない場合，即ち θi の分布が一様分布となる場
合は，2 次元平面において出鱈目な方向を持つ長さ 1 のベクトルを繋いでいった，いわゆる酔歩
問題と見做すことが出来，2R2 /N の分布は，近似的に自由度 2 のカイ 2 乗分布に従うことが知「地震活動と地球られている．よって，実際にデータから計算した R の値が R0 であった場合，潮汐の間に相関がない」という帰無仮説の下でこのような値が偶然得られる確率（有意確率）
は，(2.2)

P (R ≥ R0 ) = exp(−R2 /N )

地球潮汐と地震活動との相関を用いた地震活動予測131

図 1．2011 年東北地方太平洋沖地震の前 3000 日間のデータに対し，200 km × 200 km の空
間ウィンドウを，地震断層の辺に対して平行に 50 km ずつ動かしつつ，各々のウィンド（Tanaka,
ウ内の地震に対して Schuster’s test の P 値を計算した場合の，その空間分布
2012 より引用）．
Fig. 1. Spatial variation in P -values from Schuster’s test for earthquakes that occurred
during the 3000 days before the 2011 oﬀ the Paciﬁc coast of Tohoku Earthquake.
The spatial variation was computed with spatial windows of 200 km × 200 km
shifted by 50 km along both the strike and dip directions. This ﬁgure was takenfrom Takana（2012）, and is produced by permission of American GeophysicalUnion. Copyright 2012 American Geophysical Union.

で与えられる．以降，この値を P 値と呼ぶこととする．P 値が十分に小さければ，上の帰無仮説を棄却して，「地震活動と地球潮汐の間に相関がある」と結論づけることになる．
地球潮汐は，多数の異なる周期・振幅・位相を持つ成分（分潮）の足し合わせであるため，その時刻歴は非常に複雑である．それゆえ，地震発生時系列を位相に直すにあたって上記のようないくらか手間のかかる手続きを必要とする．しかし，例えば，周期 12.42 時間の太陰半日周のような顕著な振幅を持った特定の分潮にのみ着目し，単純に θi = 2πti /T（T は着潮（M2 分潮）
目した分潮の周期）として位相を計算することもしばしばある．
Schuster’s test を用いた研究例としては，Tsuruoka et al.（1995），岩田・中西（1998），Tanaka
et al.（2002a, 2002b）, Cochran et al.（2004），Stroup et al.（2007），Li and Xu（2012）などがあを紹介する．
るが，その中から Tanaka（2012）
この論文では，2011 年東北地方太平洋沖地震の発生前に，その震源域において，地震活動と地球潮汐との間に相関があるかどうかを調べている．まず，震源域全体をほぼカバーする 500 km（図 1 の黒四角）
に対し，「Global CMT カタログ」と呼ばれる，全地球の地震系
× 200 km の領域
列を含むカタログから，その地震検知能力を考慮して Mw ≥ 5.0 の地震を取り出した．解析期

132

統計数理

第 63 巻

第 1 号 2015

図 2．図 1 の “region A” に関する Schuster’s test の P 値の時間変化．3000 日間のタイム
ウィンドウを 500 日ごとに動かして計算したもの．本震発生時は縦の灰色点線で表されている（Tanaka, 2012 より引用）．
Fig. 2. Temporal variation in P -values from Schuster’s test for earthquakes that occurred in “region A” shown in Fig. 1. The temporal variation was computedwith time windows of 3000 days shifted by 500 days. This ﬁgure was taken from
Takana（2012）, and is produced by permission of American Geophysical Union.
Copyright 2012 American Geophysical Union.

間を本震発生前 3000 日間として Schuster’s test を行なったところ，P 値は 12%となり，有意な相関は見られない．しかし，同様のデータセットに対して，200 km ×200 km の空間領域（ウィンドウ）を，本震断層領域として設定された矩形の各辺の方向に対して平行に 50 km 間隔で動かしながら，各ウィンドウに対して P 値の計算を行なったところ，図 1 に示す通り，断層領域の北の部分で有意な相関が見い出された（図 1 では，各空間ウィンドウに対して得られた P 値を，
そのウィンドウの中心にあたる 50 km × 50 km の領域の色で示してある）．最も小さい P 値は0.34%と非常に強い相関が見られ，また，これを得た領域（図 1 中で，“Region A” として示された灰色の四角で囲まれた領域）は本震震源，即ちその破壊開始点や，本震発生の 2 日前に起きた最大前震の震源を含むものである．
さらに，この “Region A” については，P 値の時間変化についても調べた．やはり 3000 日間の時間幅（タイムウィンドウ）を設定し，これを 500 日の間隔で動かしつつ，各タイムウィンド
ウ内に含まれる地震に対して，P 値を計算した．図 2 に示す通り，30 年以上にわたる時間変化
を見ると，2000 年代の前半を含むような時期から P 値の減少が始まり，本震発生直前期を含む（図中，縦の灰色点線の右側），
タイムウィンドウで P 値が最小となっている．また，本震発生後
即ち，東北地方太平洋沖地震の余震活動に対しては，有意な相関は見られない．
以上のことから，本震発生に向けての応力集中，特に破壊開始点におけるそれに呼応して，地震活動と地球潮汐との有意な相関が見られるようになり，本震発生により蓄積された応力が解放されたため，相関が見られなくなったというシナリオが示唆される．また，他の大地震に対，Tanaka et al.（2002b）, Tanaka（2010）
は同様の現象を示しておする解析として，鶴岡（1995）
り，必ずしも東北地方太平洋沖地震に限られたものではないようである．
この小節の最後に，Schuster’s test を適切に用いるための留意点について簡単に述べる．まず，その適用にあたっては，地震のクラスター性を事前に取り除く必要がある．これは大地震
直後の激しい余震活動のような，非常に短い時間に集中して地震が起きるようなことがあれば，当然，その時間帯に対応する位相に地震が集中することになり，周期性とは無関係に P 値は小さくなるからである．これを避けるためには，地震のクラスター性を取り除く「デクラスタリング」と呼ばれる処理を，事前に施す必要がある．それにも関わらず，この点に全く頓着せずにSchuster’s test を用いた研究も散見され，こういった研究の結果には留意する必要がある．ま地球潮汐と地震活動との相関を用いた地震活動予測

133

た，デクラスタリングを行なっていたとしても，これまでのところ，確立されたデクラスタリ
ング手法はないことを鑑みると，例えば複数の手法を用いて，結果の頑健性を確かめるような慎重さが不可欠であろう．
加えて，もし真の位相分布に何らかの偏りがあった場合であっても，P 値の大小はデータ数にも大きく左右される．言い換えると，相関の時空間変化を調べるにあたって，Tanaka（2012）に代表される彼女（ら）の一連の研究では，同じ時間あるいは空間幅のウィンドウを動かして P
値の変化を見ているが，これが必ずしも適切とは言えない．例えば，同じ地震数を含むようにウィンドウの幅を可変にして解析を行なう工夫も必要かもしれない．また，ウィンドウ幅（あるいは各ウィンドウに含まれる地震数）の設定を変えると，どのように結果が変化するのかについても検証が必要と思われる．
2.2 Load/Unload Response Ratio（LURR）
前節で紹介した Schuster’s test の発想といくらか重なるところがある手法として，Load/Unload（Yin et al., 1995）
がある．LURR を用いた比較的近年までの研究は，
Response Ratio（LURR）
Yin et al.（2006）に一通りまとめられている．また，最近は，このアプローチの発展形として，
ETAS モデル（Ogata, 1988）と組み合わせることで，LURR を次節で述べる点過程解析の枠組み．
で扱う手法も提案されている（Zhang and Zhuang, 2011; 3 節でも触れる）
Yin et al.（1995）
で提示された手法について紹介する．1 節の冒頭で触れたように，地球潮汐による応力変化は周期的なものであるため，これにより解析対象とする地震の発生を促すようなと，逆に発生を抑制するような「負」の応力が生じる期「正」の応力が生じる期間（loading period）
が交互に訪れる．それゆえ，地震をその発生時刻に応じて，loading period
間（unloading period）に起きたものと unloading period に起きたものとの 2 群に分けることが出来る．そして，各々「正」の群に関する和を，「負」の群に関すの群に属する地震のエネルギーの m 乗の和を計算し，Y を求める．式で表すと以下の通りである．
るそれで割ったもの（即ち，両者の比）
N −
N+

 − m
+ m
(2.3)
Y =
(Ei )
(Ei ) .
i=1

i=1

ここで，Ei+ および Ei− はそれぞれ，「正」の群および「負」の群における i 番目の地震のエネル+−ギーであり，N および N はそれぞれの群に含まれる地震数である．また，m を 0 にすれば，2 群における地震数を比較していることになり，m = 1/2 とすれば，いわゆる Benioﬀ 歪みを考を始めとする，彼らの一連の研究では m = 1/2 が使われる
えることになる．Yin et al.（1995）
ことが多い．
また，地震が完全にランダムである（潮汐と相関がない）場合に，Y が従う分布を数値シミュに記されている．これを用いて，Schuster’s test
レーションで求める方法は，Yin et al.（2000）のように統計的検定を行なうことも可能である．
実際に LURR を用いた研究では，有意性の検定そのものより，ある一定の時間幅（タイムウィ
ンドウ）
を設け，そのタイムウィンドウを動かしつつ Y を計算することにより，その時間変化を調べ，大地震発生の前兆として使えるかどうかを議論することが多い．例えば，Yin et al.
（2000）は南カリフォルニアの地震データに対して LURR の時間変化を計算した．図 3 に示す通り，大地震発生（図中の矢印）の前に，Y の顕著な上昇が見られるとしている．また，その有意性は上記の検定で確認されている．
但し，この結果に関連しては，特に 1994 年ノースリッジ地震前の Y の変化を詳しく検証し直し，有意な前兆は見られないとする研究がある（Smith and Sammis, 2004）．さらに，Trotta
は，Yin との個人的なやり取りと，Smith and Sammis（2004）にある情報を元
and Tullis（2006）

134

統計数理

第 63 巻

第 1 号 2015

図 3．Load/Unload Response Ratio（LURR）の適用例として，南カリフォルニアの地震系列を解析したもの．式（2.3）で定義された Y の値の時間変化を見ると，いくつかの大地
震（図中の矢印）の前に有意な上昇が見られる（Yin et al., 2000 より引用）．
Fig. 3. Application of the Load/Unload Response Ratio（LURR）to several earthquake
sequences in southern California. The values of Y deﬁned in eq.（2.3）signiﬁcantly
increase prior to the occurrences of the large earthquakes indicated by arrows.
This ﬁgure was taken from Yin et al.（2000）.

に，それらの研究が示した Y の時間変化を再現しようと試みたが，再現できなかったと記している．そして，LURR に基づく解析は，論文などには詳らかに現れない細かいパラメータの設定に大きく左右されると結論している．
また，図 3 として引用した，Yin et al.（2000）の図にある元のキャプションには “all strong地球潮汐と地震活動との相関を用いた地震活動予測

135

earthquakes that occurred in southern California from 1980 to 1994”（の前の Y の時間変化を調べた）とある．しかし，ここで取り上げられた地震の近傍では，他にもいくつか M 6 程度の地震があり，“strong earthquakes” の選択基準がよく分からない．何れにしても，LURR が地震予測に有効であるか，またそもそも地震活動と地球潮汐との相関を測るのに適切な手段であるかど
うか，詳しい検討が必要であろう．
2.3 点過程解析に基づくアプローチ
地震の発生事象のように，時間軸上（実際には，空間または時空間にも拡張出来るが，ここでは単純化のため，時間のみを扱うことにする）にその発生時刻を「点」として抽象的に表し得る確率的事象は点過程と呼ばれる（例えば Daley and Vere-Jones, 2003．地震学への応用に関しては，
Ogata, 1999, 2013 の解説論文も参照）．こういったデータの解析においては，その点パターン，
具体的には単位時間当たりの点の個数（即ち，密度）を表す「強度関数」（intensity function）を推定することが，まず重要となる．なお，地震データの解析においては，強度関数は「地震活動度」（seismicity rate）と実質上同義と考えて差し支えない．
今，何らかの仮定に基づき，時刻 t における強度関数 λ(t|θ) が与えられており，この強度関数のモデルパラメータ θ の推定を行ないたいとする．解析する地震データ（発生時系列）におい（N はデータ内の地震て，ti が i 番目の地震の発生時刻とすると，λ(t|θ) の {ti ; i = 1, 2, . . . , N }数）
に対する当てはまり具合いを表す対数尤度 ln L(θ) は，(2.4)

ln L(θ) =N

i=1

ln λ(ti |θ) −T

λ(t|θ)dt
S

で与えられる．ここで，S および T はそれぞれ，解析期間の始まりと終わりの時刻である．
さて，この枠組みで，地震活動に何らかの周期的な変化があるかどうかを調べたいとする．こういった場合の強度関数として，Ogata（1983）は，(2.5)

λ(t) =（トレンド）+（クラスター）+（周期性）というものを用いることを提案している．地震活動の場合，「トレンド」
は長期的な地震活動の変化，クラスターは余震活動に代表される地震活動の群れに相当する．
周期性が有意であるかどうかについては，（2.5）式の強度関数から周期性成分を取り除いたモ（Akaike, 1974）を計算し，比較デルとそうでないモデル各々について，赤池情報量規準（AIC）
することで判断できる（いわゆる「モデル選択」）．よく知られているように，AIC はモデルパラ
メータ数の違いを考慮した上で，モデルとデータとの合い具合いを比較するための統計的指標
である．
Iwata and Katao（2006）（以下，IK2006）
は，この考え方に基づいて，丹波山地の微小地震活動について調べた．この地域は，1995 年兵庫県南部地震の震源域に隣接しており，兵庫県南部地震の発生から数年間にわたって，地震活動と月齢との間に相関が見られることが，定性的に片尾（2002）によって指摘されている．
IK2006 で用いた強度関数においては，トレンド成分を多項式で表すこととした．クラスター
では，当時の計算機における計算負荷を考慮して，Laguerre 型
成分については，Ogata（1983）
多項式を用いることが提案されているが，その後の計算機能力の向上を鑑みて，より複雑ながらも地震学では現在のところ標準的な地震活動モデルとされている ETAS モデル（Ogata, 1988）を導入した．そして，周期性成分については，三角関数を用いることとした．
よって，具体的な強度関数の式は，以下の通りになる．

統計数理

136

第 63 巻

第 1 号 2015

表 1．式
（2.6）に基づく点過程解析の結果．調べた 4 つのモデルに対する AIC の値と，選択されたトレンド成分に使われた多項式の次数 J ．4 つのモデルの中で，AIC 最小となった
ものを太字で示す．
Table 1. Akaike’s Information Criterion（AIC）values and the values of J, the order ofthe polynomials in the trend component obtained for the intensity function asshown in eq.（2.6）. The minimum AIC values for the examined four models arehighlighted in bold characters.

(2.6)

λ(t|θ) = μ +

J

k=1

ak t k + K exp(α(Mi − Mz ))i;ti <t(t − ti + c)p+A1 sin θ(t) + B1 cos θ(t) + A2 sin(2θ(t)) + B2 cos(2θ(t)).

ここで，θ(t) は時刻 t を月齢に関する位相に変換する関数である．即ち，ある連続する新月の時刻を 0◦ または 360◦ ，その間の満月の時刻を 180◦ とし，その間の時刻は線形に角度へと割り振ったものである．これにより，θ(t) および 2θ(t) が，それぞれ月の満ち欠けの 1 周期（例えば，（例えば，ある新月から次ある新月から次の新月まで）に相当する 1 朔望月周期およびその半分の満月まで）となる半朔望月周期に対応することになり，以下の 4 つのモデルを調べることが出
1 朔望月周期のみあり，（iii）半朔望月周期のみあり，来る:（i）月齢に関連した周期性なし，（ii）
1 朔望月・半朔望月周期ともあり．そして，それぞれの場合に対応するモデルパラメータへ（iv）
A1 = A2 = B1 = B2 = 0，（ii）A2 = B2 = 0，（iii）A1 = B1 = 0，（iv）制の制約はそれぞれ，（i）
約なし，となる．また，トレンド成分に相当する多項式については，その次数 J を定める必要があるが，モデル選択の考え方と同様に，AIC に基づいて適切な次数を選ぶことにした．
IK2006 では，京都大学防災研究所阿武山系観測網により決定された，丹波山地の地震を解析が扱ったものと同様とし，兵庫県南部地震発生から約 2した．具体的な地域設定は片尾（2002）
年間（片尾（2002）とほぼ同じ 1995 年 1 月 17 日正午から 1996 年 12 月 11 日正午まで）の M ≥ 1.2の地震データを取り出し，式（2.6）の強度関数に基づくモデル比較を行った結果を表 1 に示す．
AIC が小さい方が，よりよいモデルであることから，（iv）のケース，即ち「1 朔望月・半朔望月
周期ともあり」とするモデルが最もよいことが分かる．また，他の 3 つのモデルとの AIC の差
は 2.45 またはそれ以上であることから，両周期が存在することは統計的に十分有意であること．さらに，比較のも言える（より詳しくは，IW2006 に示したカイ 2 乗分布に基づく検定を参照）（1993 年 1 月 17 日 0 時から 1995 年 1 月 16 日 24 時）におけため，兵庫県南部地震発生前 2 年間
る同じ地域の地震データを解析したところ，「月齢に関連した周期なし」に相当する（i）のケースが選択された（表 1）．即ち，月齢に関連した地震活動の周期性は，兵庫県南部地震発生後のみに見られることになる．

地球潮汐と地震活動との相関を用いた地震活動予測

137

表 2．式
（2.7）に基づく点過程解析の結果．調べた 4 つのモデルに対する AIC の値と，選択されたトレンド成分および周期性成分に使われた多項式の次数 J, L1 , L2 ．4 つのモデルの中で，AIC 最小となったものを太字で示す．
Table 2. Akaike’s Information Criterion（AIC）values and the values of J, L1 , and L2 ,the orders of the polynomials in the trend and periodic components obtained forthe intensity function as shown in eq.（2.7）. The minimum AIC values for theexamined four models are highlighted in bold characters.

次に，周期性がどのように時間変化しているかを調べることを試みた．式（2.6）では，周期性成分に対応する三角関数の係数を定数としている．しかし，この係数が時間変化出来るようにモデルを拡張することで，周期性の強度の時間変化に関する解析を行ない得る．IK2006 では単純ではあるが，各三角関数の係数が，トレンド成分と同様に多項式で表されるものとした．即ち，強度関数は以下のように拡張される．
J K exp(α(Mi − Mz ))
ak t k +λ(t) = μ +(2.7)(t − ti + c)pi;t <t
k=1
i

+

L1

k=1

+

L2

k=1

A1k tk−1 · sin θ(t) +L1


B1k tk−1 · cos θ(t)k=1

A2k tk−1 · sin(2θ(t)) +L2


B2k tk−1 · cos(2θ(t)).

k=1

1 朔望月および半朔望月周期に対応する項に現れる多項式の次数 L1 および L2 は，トレンド成（2.6）を強度関数に分に対応する多項式の次数 J と同様に，AIC に基づいて定める．但し，式用いた解析の場合と同様に，1 朔望月および半朔望月に対応する周期性の時間変化がある・な
しの場合に応じて，4 つのモデルを調べることとなり，各モデルは次に示すように，L1 および L2 に対する制約条件として表すことが出来る:（i）L1 = L2 = 0，（ii）L1 ≥ 1, L2 = 0，（iii）
L1 = 0, L2 ≥ 1，（iv）L1 ≥ 1, L2 ≥ 2．
周期性の時間変化を詳らかにするために，解析期間を 2 年間から 4 年間（1995 年 1 月 17 日正午から 1999 年 1 月 17 日正午まで）と延長した上で，この解析を行なった．まず，AIC による 4
つのモデルの比較結果を表 2 に示す．モデル（iv），即ち，1 朔望月および半朔望月周期とも，その強度に時間変化があるとするモデルが最適のものとして選ばれた．さらに，他の 3 つのモデルとの AIC の差は 4 程度またはそれ以上であり，その有意性も確かと言える．
図4（a）に，モデル（iv）
より得られた周期性成分の時間変化を示す．そして，周期性成分の強度の時間変化を分かりやすくするために，式（2.7）にある三角関数の振幅の時間変化を次の式で

統計数理

138

第 63 巻

第 1 号 2015

図 4．（a）兵庫県南部地震発生後 4 年間
（1995 年 1 月 17 日 – 1999 年 1 月 17 日）の丹波山地の地震活動を，式（2.7）の強度関数を用いた点過程解析で調べた結果，得られた周期性成（a）に示した周期性成分の時間変化を，式（2.8）の gi (t) により三角分の時間変化．
（b）関数の振幅，即ち周期性成分の強度に変換したもの．実線および点線がそれぞれ 1 朔望
（i = 2）に対応する．また，横軸は解析期間の始点からの月（i = 1）および半朔望月周期
経過日数である（Iwata and Kanao, 2006 より引用）．
Fig. 4.（a）Temporal variation in the intensity function of the trigonometric terms ineq.（2.7）that were obtained by analyses of microearthquakes that occurred from
17 January 1995 to 17 January 1999.（b）Temporal variation in the amplitudesof the trigonometric terms, which are given as gi (t) in eq.（2.8）. The solid anddotted lines correspond to the periodicities of synodic（i = 1）and half-synodic（i = 2）months, respectively. This ﬁgure is taken from Iwata and Kanao（2006）,
and is produced by permission of American Geophysical Union. Copyright 2006
American Geophysical Union.

計算した．
(2.8)


2  L
2
 Li
i
 


gi (t) =Aik tk−1+
Bik tk−1 (i = 1, 2).
k=1

k=1

i = 1 および 2 は，それぞれ 1 朔望月，半朔望月周期の三角関数に対応する．計算された g1 (t)および g2 (t) を合わせて図 4（b）に示す．
図に示した通り，1 朔望月・半朔望月周期とも，その強度は本震発生直後に最も強く，その後，時間を追うにつれ減衰している．また，兵庫県南部地震発生前 4 年間（1991 年 1 月 17 日 0 時かのデータについても，同じ解析を行なったが，これについては，ら 1996 年 1 月 16 日 24 時まで）
，即ち 1 朔望月・半朔望月周期ともなしとするものが最適であった（表 2）．以上をまモデル（i）
とめると，兵庫県南部地震前には，地震活動と地球潮汐との相関はなかったが，その発生直後に地震活動と地球潮汐との相関が最も強く，時間の経過と共に相関が失われていったと言える．
IK2006 で解析対象とした丹波山地は，兵庫県南部地震の断層運動により正の応力変動（クーロン応力変化）を受けたことが指摘されている（Toda et al., 1998）．よって，この正の応力変動を原因として，丹波山地の応力レベルが臨界状態となることで，一時的に強い地震活動と地球潮地球潮汐と地震活動との相関を用いた地震活動予測

139

汐との相関が生じた．そして，その後，丹波山地およびその周辺域における地震活動による応力再配分，あるいは地殻の粘弾性的な振る舞いによる応力緩和などにより，応力レベルが徐々に低下し，これに伴って，相関が弱まっていたと考え得る．
3.

今後への展開と地震予測

1 節でも述べた通り，地球潮汐と地震活動との相関は存在するのか，また存在するにしても，あくまで限られた期間・地域にのみ現れるのか，さらにそうであればどういった期間・地域に現れるのか，こういった疑問は未解決であり，これに向けての網羅的な解析が望まれる．
，鶴岡（1995），Tanaka et al.（2002b），Tanaka（2010, 2012）な
特に，Yin et al.（1995, 2000）
どが示す大地震発生の前兆現象としての地球潮汐と地震活動との相関検出は，地震予測の観点から鑑みても魅力的である．重要なことは，上記の一連の研究のような「後予知」ではなく，大地震発生前に検出できる体勢を整えることである．例えば，Schuster’s test に依るのであればその P 値，LURR であれば Y の計算を常時行ない，モニタリングすることが必要である．
但し，例えば Schuster’s test であれば，Tanaka（2012）が指摘している通り，対象とする地震によって，その前に現れる P 値低下の期間はかなり異なっている．言い換えると，前兆現象としての P 値低下を検出できるかどうかは，タイムウィンドウの幅を適切に設定できるかどうかに大きく依存するが，この設定に関する具体的な指針がないのが現状である．同じ問題は，空間ウィンドウの大きさの設定に関しても生じるし，これは LURR においても共通の課題である．
勿論，様々な大きさの時間・空間ウィンドウを用いたモニタリングを同時に行なうことは，こういった問題点に対する一策ではあるし，こういったモニタリングを介して，適切なウィンドウ設定に関する情報が得られる可能性もある．
他方，2.3 節に示した点過程解析によるモデリングは，時間・空間ウィンドウの設定が必要な
いため，このような問題は回避することが出来る．また，2.3 節に示した例では，相関の空間（Ogata, 1988）が時空間版（Ogata, 1998; Ogata
変動までは考えていないが，元々の ETAS モデル
へと拡張されているように，式（2.5）のような，地震活動の周期性を調べるた
and Zhuang, 2006）
めの強度関数も，時空間へと拡張することは，理屈の上ではさほど難しいことではない．
但し，具体的にどういう関数を用いて強度関数を書き表すかについて，考える必要があるし，それ以前に，地震活動を構成する各成分をどう組み合わせるかについても考察の余地がある．
2.3 節では，式（2.5）に基づいて，式（2.7）のような強度関数を設定して，相関の時間変化についのように各成分の和を取る形ではなく，て調べた．しかし，式
（2.5）
(3.1)

λ(t) = （トレンド）+（クラスター） ×（周期性）のように，周期性の効果を，それ以外の成分に積の形で与えることも考えられる．
実際，岩田・片尾（2006）
では，式（2.5）に替え，式（3.1）に従う形で強度関数を構成し直して，
IK2006 で扱ったものと同じデータを解析した．クラスター成分は ETAS モデルで表し，トレン
ド成分及び周期性成分の三角関数の係数部分を多項式で表した点は，IK2006 と同様である．再解析の結果，IK2006 で行なった解析に比べ，AIC の値は 5 またはそれ以上よくなり，モデルの改善が示された．
LURR に関連しては，2.2 節で触れた Zhang and Zhuang（2011）
が，LURR の統計量 Y（式（2.3））に ETAS モデルを掛けたものを強度関数とする手法を提案しており，四川地震の余震活と同様の発想に基づくものである．
動を解析している．このモデルは式（3.1）
このような試みは行なわれつつあるにしても，点過程解析の土俵で，地震活動と潮汐との相関を扱うための種々のモデル提案と，実際のデータ解析に基づくその検証は単純な作業ではな

統計数理

140

第 63 巻

第 1 号 2015

い．特に時空間へと問題を展開した場合には，かなり複雑なものとなり，一層の慎重な考察が求められるであろう．また，Schuster’s test や LURR はその統計量の計算が比較的簡単であるのに対し，点過程解析における対数尤度の計算とその最大化（モデルパラメータ推定）
は，往々にして数値計算上の困難を伴うことになる点に留意する必要がある．但し，こういった困難は，今後の計算機能力の向上とともに，ある程度は克服出来ることが期待される．
最終的には，CSEP Japan（Nanjo et al., 2011）の枠組みに乗るような，地震活動と地球潮汐との相関の時空間変動を情報として取り入れた予測モデルの構築が望まれる．しかしながら，上述の通り，そこに至るまでには取り組むべき多くの課題があることから，自分自身も含め，この分野に取り組む研究者の奮起が必要である．
謝辞

本稿を書く機会を設けて下さいました本特集号オーガナイザーの尾形良彦氏，有益なコメン
トを下さった匿名の査読者に感謝致します．また，本稿執筆の一部は，科学研究費補助金・基（課題番号 26330056）
を用いて行ないました．
盤研究（C）

参考文

献

Akaike, H.（1974）. A new look at the statistical model identiﬁcation, IEEE Transactions on AutomaticControl, AC-19, 716–723.
Bartlett, M. S.（1963）. The spectral analysis of point processes, Journal of the Royal Statistical Society,
Series B, 25, 264–296.
Chen, H. J., Chen C. C., Tseng, C. Y. and Wang, J. H.（2012）. Eﬀect of tidal triggering on seismicity
in Taiwan revealed by the empirical mode decomposition method, Natural Hazards and Earth
System Sciences, 12, 2193–2202.
Cochran, E. S., Vidale, J. E. and Tanaka, S.（2004）. Earth tides can trigger shallow thrust fault earthquakes, Science, 306, 1164–1166.
Crockett, R. G. M., Gillmore, G. K., Phillips, P. S., Denman, A. R. and Groves-Kirkby, C. J.（2006）.
Tidal synchronicity of the 26 December 2004 Sumatran earthquake and its aftershocks, Geophysical Research Letters, 33, L19302, 10.1029/2006GL027074.
Daley, D. J. and Vere-Jones, D.（2003）. An Introduction to the Theory of Point Processes, Vol. I,
Springer, New York.
Gallego, A., Russo, R. M., Comte, D., Mocanu, V., Murdie, R. E. and Vandecar, J. C.（2013）. Tidalmodulation of continuous nonvolcanic seismic tremor in the Chile triple junction region, Geochemistry, Geophysics, Geosystems, 14, 851–863.
Heaton, T. H.（1982）. Tidal triggering of earthquakes, Bulletin of the Seismological Society of America,
72, 2180–2200.
Ide, S. and Tanaka, Y.（2014）. Controls on plate motion by oscillating tidal stress: Evidence from deep
tremors in western Japan, Geophysical Research Letters, 41, 3842–3850.
Iwata, T.（2002）. Tidal stress/strain and acoustic emission activity at the Underground Research Laboratory, Canada, Geophysical Research Letters, 29, 1126, 10.1029/2001GL014277.
Iwata, T.（2012）. Earthquake triggering caused by the external oscillation of stress/strain changes, Community Online Resource for Statistical Seismicity Analysis, 1–27, doi:10.5078/corssa-65838518.
Iwata, T. and Katao, H.（2006）. Correlation between the phase of the moon and the occurrences of
microearthquakes in the Tamba region through point-process modeling, Geophysical Research
Letters, 33, L07302, 10.1029/2005GL025510.

地球潮汐と地震活動との相関を用いた地震活動予測141

岩田貴樹, 片尾 浩（2006）. 月齢と丹波山地の微小地震発生の相関—改良した強度関数による解析—, 日本
地震学会 2006 年度秋期大会講演予稿集, p.212.
岩田貴樹，中西一郎（1998）. 長野県松代における地球潮汐と地震発生の関係, 地震, 51, 51–59.
片尾 浩（2002）. 月齢と丹波山地の微小地震発生の相関について, 地学雑誌, 111, 248–255.
Klein, F. W.（1976）. Earthquake swarms and the semidiurnal solid earth tide, Geophysical Journal of
the Royal Astronomical Society, 45, 245–295.
Knopoﬀ, L.（1964）. Earth tides as a triggering mechanism for earthquakes, Bulletin of the SeismologicalSociety of America, 54, 1865–1870.
Li, Q. and Xu, G. M.（2012）. Tidal triggering of earthquakes in Longmen Shan region: The relation to
the fold belt and basin structures, Earth, Planets and Space, 64, 771–776.
Métivier, L., de Viron, O., Conrad, C. P., Renault, S., Diament, M. and Patau, G.（2009）. Evidenceof earthquake triggering by the solid earth tides, Earth and Planetary Science Letters, 278,
370–375.
Mohler, A. S.（1980）. Earthquake/earth tide correlation and other features of the Susanville, California,
earthquake sequence of June-July 1976, Bulletin of the Seismological Society of America, 70,
1583–1593.
Nakata, R., Suda, N. and Tsuruoka, H.（2008）. Non-volcanic tremor resulting from the combined eﬀectof Earth tides and slow slip events, Nature Geoscience, 1, 676–678.
Nanjo, K. Z., Tsuruoka, H., Hirata, N. and Jordan, T. H.（2011）. Overview of the ﬁrst earthquake
forecast testing experiment in Japan, Earth Planets Space, 63, 159–169.
Ogata, Y.（1983）. Likelihood analysis of point processes and its applications to seismological data: A
prospect of earthquake prediction research, Bulletin of the International Statistical Institute,
50（Book 2）
, 943–961.
Ogata, Y.（1988）. Statistical models for earthquake occurrence and residual analysis for point processes,
Journal of the American Statistical Association, 83, 9–27.
Ogata, Y.（1998）. Space-time point-process models for earthquake occurrences, Annals of the Instituteof Statistical Mathematics, 50, 379–402.
Ogata, Y.（1999）. Seismicity analysis through point-process modeling: A review, Pure and AppliedGeophysics, 155, 471–507.
Ogata, Y.（2013）. A prospect of earthquake prediction research, Statistical Science, 28, 521–541.
Ogata, Y. and Zhuang, J.（2006）. Space-time ETAS models and an improved extension, Tectonophysics,
413, 13–23.
Rubinstein, J. L., La Rocca, M., Vidale, J. E., Creager, K. C. and Wech, A. G.（2008）. Tidal modulationof nonvolcanic tremor, Science, 319, 186–189.
Rydelek, P. A., Davis, P. M. and Koyanagi, R. Y.（1988）. Tidal triggering of earthquake swarms at
Kilauea Volcano, Hawaii, Journal of Geophysical Research, 93, 4401–4411.
Rydelek, P. A., Sacks, I. S. and Scarpa, R.（1992）. On tidal triggering of earthquakes at Campi Flegrei,
Italy, Geophysical Journal International, 109, 125–135.
Schuster, Y.（1897）. On lunar and solar periodicities of earthquakes, Proceedings of the Royal Societyof London, 61, 455–465.
Shelly, D. R., Beroza, G. C. and Ide, S.（2007）. Complex evolution of transient slip derived from precise tremor locations in western Shikoku, Japan, Geochemistry, Geophysics, Geosystems, 8,
Q10014, doi:10.1029/2007GC001640.
Smith, S. W. and Sammis, C. G.（2004）. Revisiting the tidal activation of seismicity with a damage
mechanics and friction point of view, Pure and Applied Geophysics, 161, 2393–2404.
Souriau, M., Souriau, A. and Gagnepain, J.（1982）. Modeling and detecting interactions between earth
tides and earthquakes with application to an aftershock sequence in the Pyrenees, Bulletin of

142

統計数理

第 63 巻

第 1 号 2015

the Seismological Society of America, 72, 165–180.
Stroup, D. F., Bohnenstiehl, D. R., Tolstoy, M., Waldhauser, F. and Weekly, R. T.（2007）. Pulse of
the seaﬂoor: Tidal triggering of microearthquakes at 9◦ 50 N East Paciﬁc Rise, GeophysicalResearch Letters, 34, L15301, 10.1029/2007GL03008.
Tanaka, S.（2010）. Tidal triggering of earthquakes precursory to the recent Sumatra megathrust earthquakes of 26 December 2004（Mw 9.0）
, 28 March 2005（Mw 8.6）and 12 September 2007（Mw
8.5）, Geophysical Research Letters, 37, L02301, 10.1029/2009GL041581.
Tanaka, S.（2012）. Tidal triggering of earthquakes prior to the 2011 Tohoku-Oki earthquake（Mw 9.1）,
Geophysical Research Letters, 39, L00G26, 10.1029/2009GL051179.
Tanaka, S., Ohtake, M. and Sato, H.（2002a）. Evidence for tidal triggering of earthquakes as revealed from statistical analysis of global data, Journal of Geophysical Research, 107, 2211,
10.1029/2001JB001577.
Tanaka, S., Ohtake, M. and Sato, H.（2002b）. Spatio-temporal variation of the tidal triggering effect on earthquake occurrence associated with the 1982 South Tonga earthquake of Mw 7.5,
Geophysical Research Letters, 29, 1756, 10.1029/2002GL015386.
Thomas, A. M., Nadeau, R. M. and Bürgmann, R.（2009）. Tremor-tide correlations and near-lithostatic
pore pressure on the deep San Andreas fault, Nature, 462, 1048–1051.
Thomas, A. M., Bürgmann, R., Shelly, D. R., Beeler, D. R. and Rudolph, M. L.（2012）. Tidal triggering of low frequency earthquakes near Parkﬁeld, California: Implications for fault mechanics within the brittle-ductile transition, Journal of Geophysical Research, 117, B05301,
doi:10.1029/2011JB009036.
Toda, S., Stein, R. S., Reasenberg, P. A., Dieterich, J. M. and Yoshida, A.（1998）. Stress transferredby the 1995 Mw = 6.9 Kobe, Japan, shock: Eﬀect on aftershocks and future earthquake probabilities, Journal of Geophysical Research, 103, 24543–24565.
Tolstoy, M., Vernon, F. L., Orcutt, J. A. and Wyatt, E. K.（2002）. Breathing of the seaﬂoor: Tidal
correlations of seismicity at Axial Volcano, Geology, 30, 503–506.
Trotta, J. E. and Tullis, T. E.（2006）. An independent assessment of the load/unload response ratio（LURR）proposed method of earthquake prediction, Pure and Applied Geophysics, 163,
2375–2387.

鶴岡 弘（1995）. 地震発生における地球潮汐の影響とその解釈, 博士論文, 東北大学大学院理学研究科.
Tsuruoka, H., Ohtake, M. and Sato, H.（1995）. Statistical test of the tidal triggering of earthquakes:
Contribution of the ocean tide loading eﬀect, Geophysical Journal International, 122, 183–194.
Vidale, J. E., Agnew, D. C., Johnston, M. J. S. and Oppenheimer, D. H.（1998）. Absence of earthquake correlation with Earth tides: An indication of high preseismic fault stress rate, Journalof Geophysical Research, 103, 24567–24572.
Wilcock, W. S. D.（2001）. Tidal triggering of microearthquakes on the Juan de Fuca Ridge, GeophysicalResearch Letters, 28, 3999–4002.
Wilcock, W. S. D.（2009）. Tidal triggering of earthquakes in the Northeast Paciﬁc Ocean, GeophysicalJournal International, 179, 1055–1070.
Yin, X. C., Chen, X. Z., Song, Z. P. and Yin, C.（1995）. A new approach to earthquake prediction: The
Load/Unload Response Ratio（LURR）Theory, Pure and Applied Geophysics, 145, 701–715.
Yin, X. C., Wang, Y. C., Peng, K. Y., Bai, Y. L., Wang, H. T. and Yin, X. F.（2000）. Development of
a new approach to earthquake prediction: Load/unload response ratio（LURR）theory, Pureand Applied Geophysics, 157, 2365–2383.
Yin, X. C., Zhang, L. P., Zhang, H. H., Yin, C., Wang, Y., Zhang, Y., Peng, K., Wang, H., Song, Z.,
Yu, H. and Zhuang, J.（2006）. LURR’s twenty years and its perspective, Pure and AppliedGeophysics, 163, 2317–2341.

地球潮汐と地震活動との相関を用いた地震活動予測

143

Zhang, L. and Zhuang, J.（2011）. An improved version of the Load/Unload Response Ratio method
for forecasting strong earthquakes, Tectonophysics, 509, 191–197.

144

Proceedings of the Institute of Statistical Mathematics Vol. 63, No. 1, 129–144 (2015)

Earthquake Forecasting Based on the Correlation between Earth Tides andEarthquake Occurrences
Takaki Iwata1,2
1 College of Community Development, Tokiwa University
2 The Institute of Statistical Mathematics

Many studies have investigated and discussed the correlation between the occurrence of earthquakes and periodic stress/strain changes due to earth tides. This article
introduces typical statistical methods that can be used to examine such correlation andreviews recent achievements made in past studies. In particular, it is proposed that thesignificance/non-significance of correlation can be used to measure the stress state of theearth’s crust and that this information may be valuable for earthquake forecasting. Additionally, this article suggests possible development of new statistical approaches in the
future that would allow for more comprehensive evaluation of such correlation.

Key words: Seismicity, earth tides, point-process analysis, earthquake forecasting.