[S14P-02] Enhancing Statistical and Machine Learning Methods for Early Estimation of GR Law and MO Formula Parameters in Aftershock Activity
1. 現状の余震予測
2016年に発生した熊本地震では，一連の地震活動の中で震度6弱以上の地震が7回発生し，甚大な被害を与えた．そのため，大地震発生の際には迅速に余震活動についての予測を行い，一般に公開してゆく必要がある．余震の確率評価の手法は地震調査研究推進本部地震調査委員によって提案されている．具体的には，地震の規模別頻度分布を表すGutenberg-Richter則（以降“GR則”と称する）と，余震活動の単調減少を表す改良大森公式（以降“MO式”と称する）の二つの経験則を用いて確率論的に余震活動を評価している．
2. 余震観測の問題点
余震の観測にあたって，マグニチュードの大きい地震と比べて小さい地震の方が地震計での観測が難しいという観測技術上の問題点が存在し，この問題は本震発生直後に，より顕著にみられる．小さい規模の地震の記録が脱落した状態ではGR則のbやMO則のKなどの予測式のパラメータが過小評価されてしまう．これにより，本震発生直後に正確な余震予測情報を公開することは困難である．
3. 使用データと研究手法
データセットには，気象庁震源データを使用する．1997年10月から2023年3月を対象期間とし，本震はマグニチュード6.0以上の地震，余震は本震発生後90日間のマグニチュード2.5以上の地震として機械的に抽出し，解析に用いる本震余震系列データを作成した．予測式の各パラメータの算出方法については，最尤法に基づく式によってGR則のbを推定した．MO式のパラメータは地震調査研究推進本部地震調査委員会の手法を元に算出した．上記の手法により，抽出された35系列の余震系列においてGR則のbと，MO式のKが本震発生直後に実際に過小評価されていること，またその過小評価の度合いと本震のマグニチュード（以降“M0”と称する）との間に相関があることを確認した．本研究では，本震直後のパラメータの過小評価の度合いを表す値b’, K’を以下の様に定義した．ここで下付き文字3hは本震発生時を起点として，3時間分のデータを使用して算出されたパラメータであること，1dは1日分のデータを使用して算出されたパラメータであることを表している．
b' = b1d/b3h
K' = K1d/K3h
b’, K’が大きいほど，本震後3時間時点でのパラメータが本震後 1日後のパラメータから過小評価されるという関係を表す．本研究では，このb’, K’を説明変数として，単回帰モデル，lassoモデル，サポートベクター回帰モデルの3つのモデルを構築した．この時，単回帰モデルにおいては前述の関係よりM0を説明変数とした．lassoモデルおよび，サポートベクター回帰モデルにおいてはM0に加えて，余震系列のマグニチュードの時系列データをPythonライブラリのtsfreshを使用して特徴量化したものを使用した．時系列特徴量はlassoモデルのアルゴリズムにより説明変数選択を行った．本研究では，これらのモデルを使用して，本震発生後3時間時点で1日後までの余震予測式のパラメータを推定し，予測精度の向上を目指す．
4. 結果と考察
本研究で作成したモデルを使用して推定したb', k'を元に算出されたb1d, K1dと実測値との誤差を比較したものが添付の図である．現行の手法における，b1dの推定値を過大評価，K1dの推定値を過小評価してしまうといった系統的な誤差が本研究の手法では改善されていることが分かる．算出したb1d, K1dを使用して実際に過去の本震余震系列において本震後3時間後時点で1日後までの余震数を予測し実測値との誤差を比較したものも，本研究の手法は現行の手法と比較して平均的に精度が良いことが分かる．また，200回のモンテカルロリサンプリングによるb’の推定を行うlassoモデルの学習によって選択された説明変数のうち最も選択された変数は「余震マグニチュードの上位70%の範囲にある値の尺度」となっている．これが大きいことは余震のマグニチュードが大きい傾向を表しており，活動全体が活発であることを間接的に表し，小さな規模の余震の脱落をM0よりも正確に表していられるのではないかと推測ができる．M0は全体で6番目に重要な説明変数であった．選択された説明変数を解釈するうえでは，tsfreshで作成された説明変数は機械的に作成された時系列特徴量であり，必ずしも小さな規模の余震の脱落を直接的に表現する値ではないということに留意しなければならない．また，lassoモデルにより選択された説明変数はサポートベクター回帰モデルにおいても精度の良い予測結果をもたらすことも分かった．
2016年に発生した熊本地震では，一連の地震活動の中で震度6弱以上の地震が7回発生し，甚大な被害を与えた．そのため，大地震発生の際には迅速に余震活動についての予測を行い，一般に公開してゆく必要がある．余震の確率評価の手法は地震調査研究推進本部地震調査委員によって提案されている．具体的には，地震の規模別頻度分布を表すGutenberg-Richter則（以降“GR則”と称する）と，余震活動の単調減少を表す改良大森公式（以降“MO式”と称する）の二つの経験則を用いて確率論的に余震活動を評価している．
2. 余震観測の問題点
余震の観測にあたって，マグニチュードの大きい地震と比べて小さい地震の方が地震計での観測が難しいという観測技術上の問題点が存在し，この問題は本震発生直後に，より顕著にみられる．小さい規模の地震の記録が脱落した状態ではGR則のbやMO則のKなどの予測式のパラメータが過小評価されてしまう．これにより，本震発生直後に正確な余震予測情報を公開することは困難である．
3. 使用データと研究手法
データセットには，気象庁震源データを使用する．1997年10月から2023年3月を対象期間とし，本震はマグニチュード6.0以上の地震，余震は本震発生後90日間のマグニチュード2.5以上の地震として機械的に抽出し，解析に用いる本震余震系列データを作成した．予測式の各パラメータの算出方法については，最尤法に基づく式によってGR則のbを推定した．MO式のパラメータは地震調査研究推進本部地震調査委員会の手法を元に算出した．上記の手法により，抽出された35系列の余震系列においてGR則のbと，MO式のKが本震発生直後に実際に過小評価されていること，またその過小評価の度合いと本震のマグニチュード（以降“M0”と称する）との間に相関があることを確認した．本研究では，本震直後のパラメータの過小評価の度合いを表す値b’, K’を以下の様に定義した．ここで下付き文字3hは本震発生時を起点として，3時間分のデータを使用して算出されたパラメータであること，1dは1日分のデータを使用して算出されたパラメータであることを表している．
b' = b1d/b3h
K' = K1d/K3h
b’, K’が大きいほど，本震後3時間時点でのパラメータが本震後 1日後のパラメータから過小評価されるという関係を表す．本研究では，このb’, K’を説明変数として，単回帰モデル，lassoモデル，サポートベクター回帰モデルの3つのモデルを構築した．この時，単回帰モデルにおいては前述の関係よりM0を説明変数とした．lassoモデルおよび，サポートベクター回帰モデルにおいてはM0に加えて，余震系列のマグニチュードの時系列データをPythonライブラリのtsfreshを使用して特徴量化したものを使用した．時系列特徴量はlassoモデルのアルゴリズムにより説明変数選択を行った．本研究では，これらのモデルを使用して，本震発生後3時間時点で1日後までの余震予測式のパラメータを推定し，予測精度の向上を目指す．
4. 結果と考察
本研究で作成したモデルを使用して推定したb', k'を元に算出されたb1d, K1dと実測値との誤差を比較したものが添付の図である．現行の手法における，b1dの推定値を過大評価，K1dの推定値を過小評価してしまうといった系統的な誤差が本研究の手法では改善されていることが分かる．算出したb1d, K1dを使用して実際に過去の本震余震系列において本震後3時間後時点で1日後までの余震数を予測し実測値との誤差を比較したものも，本研究の手法は現行の手法と比較して平均的に精度が良いことが分かる．また，200回のモンテカルロリサンプリングによるb’の推定を行うlassoモデルの学習によって選択された説明変数のうち最も選択された変数は「余震マグニチュードの上位70%の範囲にある値の尺度」となっている．これが大きいことは余震のマグニチュードが大きい傾向を表しており，活動全体が活発であることを間接的に表し，小さな規模の余震の脱落をM0よりも正確に表していられるのではないかと推測ができる．M0は全体で6番目に重要な説明変数であった．選択された説明変数を解釈するうえでは，tsfreshで作成された説明変数は機械的に作成された時系列特徴量であり，必ずしも小さな規模の余震の脱落を直接的に表現する値ではないということに留意しなければならない．また，lassoモデルにより選択された説明変数はサポートベクター回帰モデルにおいても精度の良い予測結果をもたらすことも分かった．