# Assessing the predicting power of GPS data for aftershocks forecasting

###### Abstract

We present a machine learning approach for the aftershock forecasting of
Japanese earthquake catalogue from 2015 to 2019. Our method takes as sole
input the ground surface deformation as measured by Global Positioning System
(GPS) stations at the day of the mainshock, and processes it with a
Convolutional Neural Network (CNN), thus capturing the input’s spatial
correlations. Despite the moderate amount of data the performance of this new
approach is very promising. The accuracy of the prediction heavily relies on
the density of GPS stations: the predictive power is lost when the mainshocks
occur far from measurement stations, as in offshore regions.

Enter journal name here

Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France Université Paris-
Saclay, CNRS, INRIA, Laboratoire Interdisciplinaire des Sciences du Numérique,
TAU team, 91190 Gif-sur-Yvette, France The Institute of Statistical
Mathematics, Research Organization of Information and Systems, Tokyo, Japan

Vincenzo Maria<EMAIL_ADDRESS><EMAIL_ADDRESS>
Building a dataset based on GPS data for aftershocks forecasting.

Using machine learning algorithms for aftershock forecasting.

Evaluating the usefulness of GPS data for seismic prediction purposes.

## Plain Language Summary

Forecasting large aftershocks is a challenge of great importance for human
security. Today we dispose of statistical predictive models called Epidemic
Type Aftershock Sequence (ETAS) tuned on the earthquake catalogue of the past
seismicity. This catalogues contains basic information such as the location,
the time and the magnitude of an earthquake. However we dispose of much richer
data set about the crust dynamics, such as the daily displacement of the
ground surface, that is nowadays measured by numerous GPS stations, devices
that send their absolute position everyday to sattellites, thus telling us
about how the ground deforms. In this study, we propose to forecast the
Japanese aftershocks by means of a machine learning study of the GPS data
alone. Our results show that this method is very promising and relies on the
quality and the quantity of the available data.

##  1 Introduction

A large earthquake triggers numerous aftershocks with a frequency that
decreases with time according to the celebrated Omori law [Utsu . (1995)]. The
aftershocks magnitude is typically small, but some of them can be very large
and dangerous because they have occurred in areas that already weakened.
Forecasting their spatial location remains an outstanding open problem. Today,
for practical purposes this task is carried out using a statistical model
called Epidemic-Type Aftershock Sequence (ETAS) model [Ogata (1988), Ogata
(1998), Helmstetter Sornette (2003), Ogata Zhuang (2006)]. The ETAS is a self-
exciting point process in which aftershocks occurrence is influenced by past
events. Few free parameters define the model and they are calibrated by the
simple catalog data containing magnitude, time and location of the events
[Lippiello . (2014), Lombardi (2015), Molkenthin . (2022), Ross (2021)]. This
task is challenging due to the overlapping of seismic coda waves, which makes
difficult to detect small events in the initial stages of seismic sequences,
leading to a significant bias in the estimation of ETAS parameters [Kagan
(2004), Peng . (2007), Omi . (2013), Hainzl (20161), Hainzl (20162), De
Arcangelis . (2018), Lippiello, Petrillo, Godano . (2019)]. Recently, some
fitting procedures have been developed to address this data incompleteness,
but their effectiveness relies on the number of recorded aftershocks
[Lippiello, Petrillo, Godano . (2019), Omi . (2015), Omi . (2019)]. Another
limitation of these forecasting methods is the need to convert the estimated
aftershock rate into the observed ground shaking, via accurate empirical
attenuation functions [Boore Joyner (1982), Fukushima Tanaka (1990)].

A second strategy of aftershock forecasting relies on the study of more rich
and structured data such as the ground displacement map generated by the
mainshock. This map is estimated from the inversion of the waveform data by
means of finite-fault rupture models designed by geophysicists for each
seismic fault. The Coulomb fault stress change obtained from the map provides
a simple and geophysically inspired criterion for aftershock forecasting [King
. (1994), Toda . (1998), Parsons . (1999), Reasenberg Simpson (1992), Jacques
. (1996), Nostro . (1997)]. However its validity is still controversial
[Hardebeck . (1998), Mallman Zoback (2007), Felzer Brodsky (2005), Petrillo .
(2022)] and cannot compete with the performance of ETAS based predictions (in
particular ETAS systematically overperform on the standard CSEP test [Savran .
(2022), Savran . (2021), Khawaja . (2023)]). More recently machine learning
approaches have been explored with more promising results. In this context
neural networks, namely models defined by many free parameters, are trained
using this more rich data set. Specifically the authors in [DeVries . (2018)]
trained a Multilayer Perceptron (MLP) with the ground deformation maps to
forecast aftershock locations of different large mainshocks around the world.
As discussed above, these maps are not the outcome of a direct measurement but
rely on the inversion with complex rupture kernels. In this work we explore
the possibility to train a neural network using surface displacement maps
recorded by GPS stations. We focus on the data of Japan (1998-2018) where the
GPS stations are dense, with an average separation of $11$ km, and a
relatively good amount of independent sequences are available in the Japan
Meteorological Agency Earthquake Catalog (JMA)
111https://www.data.jma.go.jp/svd/eqev/data/bulletin.

The paper is structured as follows: Section II provides an introduction to the
data and the methods utilized. The results are presented in Section III, and
the final section is dedicated to conclusions.

##  2 Data and methods

Figure 1: The location of the GPS stations (black dots) in Japan. The training
set contains $48$ mainshocks occurred before 2015/01/01 (blue circles). The
test set contains $19$ mainshocks occurred after 2015/01/01 (red circles)
Insets: The $250$ km $\times\,250$ km maps centered around two mainshocks: the
2019 Tsuruoka earthquake (label 17)(upper left) and the 2017 offshore
Hachinohe earthquake (label 7) (lower right). The aftershocks with magnitude
larger than $4$ within $45$ days are depicted in yellow. The hatched regions
correspond to masked cells, far from GPS stations.

The earthquakes used in this study lie in the longitude interval
$[123^{\circ},148^{\circ}]$ and latitude interval $[22^{\circ},46^{\circ}]$
and are extracted from the JMA catalogue. We consider as a mainshock an
earthquake with magnitude larger than $6$ and associate to each of them the
list of the aftershocks that occurred within the next $45$ days, in the
surrounding region (square box of $250$ km of width). Only big aftershocks
with magnitude larger than $4$ are retained. Note that, in the case of two
earthquakes with magnitude larger than $6$ lying within the same spatio-
temporal window, only the first is labeled as a mainshock while the second is
labeled as an aftershock of the former. The choice of $45$ days comes as a
compromise between missing aftershocks (if the window was chosen shorter) or
discarding sequences (if the window was chosen longer). In Fig.1 we indicate
the position of all mainshocks: in blue are earthquakes that occurred before
1st January 2015 and used for training the neural network, in red are the
earthquakes that occurred after 1st January 2015 and build the test set to
assess the quality of our results.

GPS data are extracted from the Nevada Geodetic Laboratory repository (see
Blewitt2018). We use 24-hour (daily) final solutions in the frame ’IGS14’,
format ’tenv3’. Each GPS station, provides its daily position
$\vec{r}(t)=(e(t),n(t))$ (here $e(t)$ stays for the longitude (eastward) and
$n(t)$ for latitude (northward) at the day $t$). Altitude measurements are
also present, but are more noisy and thus discarded from our study. Hence, the
daily displacement of the station in the east-north plane is simply obtained
as $\vec{v}(t)=\vec{r}(t)-\vec{r}(t-1)$.

The input of the neural network are surface displacement maps of $250$ km
$\times\,250$ km centered at the mainshock and discretized in cells of $5$ km
$\times\,5$ km. This choice of $5$ km is identical to that of [DeVries .
(2018)] and allows to avoid interpolation problems that arise when choosing
overly fine-grained discretization. We take a box of $250$ km to avoid missing
out too many aftershocks, while avoiding an absurdly large window where most
cells would be inactive. The maps are obtained from the interpolation of the
GPS measurements available only at the station positions. How to perform this
interpolation is a classical issue in geophysics where the study of long time
series of GPS data is important to estimate the strain rates induced by
relative fault motion. Classical interpolation schemes are based on least
squares methods (weighted linear regression or biharmonic spline, see [Shen .
(2015)]). Here we use instead a different approach, recommended in [Sandwell
Wessel (2016)] and based on the solution of the 2D Navier-Cauchy equations at
the equilibrium, describing the displacement of a (2D) solid at rest. In Japan
the network of GPS is dense, but many earthquakes occur in the small islands
or close to the sea where GPS measurements are only partially available.
Because an interpolation cannot extrapolate displacements far away from
stations, we do not include remote regions in our study. Concretely, cells for
which the closest GPS station is further away than $d^{*}=86\rm{km}$ are
masked, i.e. they are excluded both from learning as well as forecasting. This
choice is a compromise between quality of data (better quality with smaller
$d^{*}$) and quantity of data (larger amount with larger $d^{*}$). In Fig.2
top show some surface displacements maps at the day of the mainshock, the
input of the neural network. For a given cell $(x,y)$, the displacement
$\vec{v}(x,y)=(v_{e}(x,y),v_{n}(x,y))$ is represented by an arrow (for the
direction) and a color (for the modulus). Each cell has an associated binary
label $A(x,y)$ which is $1$ if at least one aftershock with magnitude $m>4$
occurs in the $45$ days following the mainshock.

We disregard sequences in remote regions, far from GPS stations222We keep the
surface displacement maps with at least 3 stations employed for interpolation
and with at least one inside the $250\rm{km}\times 250\rm{km}$ window., and
select $67$ seismic sequences divided in $N_{\rm tr}=48$ sequences for the
training set and $19$ for the test set. For each sequence, the goal is to
predict the binary labels $A(x,y)$ given the surface displacement maps
$\vec{v}(x,y)$ as input. Because the displacement induced by the earthquake is
sometimes recorded on the next day, for each mainshock we actually provide two
surface displacement maps: the one at the day of the mainshock as well as the
one of the following day. Thus, aftershocks occurring the immediate day after
the mainshock are discarded from the prediction as well. This prediction
problem is clearly cast as a standard machine learning task. We deal with an
input image of $50\times 50$ cells of 4 channels each (the two components of
the GPS displacement and the two days), and attempt to perform image
segmentation, or in other words, binary classification of each cell. We use a
Convolutional Neural Network (CNN), which is a neural network particularly
suited for Computer Vision tasks, as it is translationally equivariant by
construction [Lecun . (1998), Krizhevsky . (2012)]: this means that a shift in
the input image induce an identical shift in the output. The loss function is
the Binary Cross Entropy loss, a classical choice for binary classification.
The input vectors of the displacement map are encoded using polar coordinates
rather than cartesian ones, and the first layers of the CNN deal with the
radial and angular part of the signal separately. The last layers combine
these inputs to produce a scalar output, i.e. a probability map
$\pi(x,y)\in[0,1]$ of aftershock occurrence. To face the issue of having such
a small training set, we use two techniques. First we perform data
augmentation, a classical technique: to each displacement map we apply all the
plane transformations that leave the map a square, namely all possible
distinct combinations of $4$ rotations of
$0^{\circ},90^{\circ},180^{\circ},270^{\circ},$ and axial flip (8 possible
variations in total). Second and most importantly, instead of using a single
neural network trained with all the maps we construct an ensemble of neural
networks. This ensemble procedure is used to have a more robust prediction, as
the training set is small [Dietterich (2000)]. In practice, all the networks
have identical architecture, but they are trained with $N_{\rm tr}-1$
independent mainshocks (hence two networks of the ensemble differ by a single
mainshock). As a result, for each mainshock $i$ on the test set, we obtain an
ensemble of $N_{\rm tr}$ predictions, indexed by $k=1...N_{\rm tr}$ :

$\displaystyle\\{\pi_{i}^{(1)}(x,y),\dots,\pi_{i}^{(k)}(x,y),\dots,\pi_{i}^{(N_{\rm
tr})}(x,y)\\}$ (1)

To determine the final prediction $\pi_{i}(x,y)$ we use three protocols:

  * •

The most cautious protocol is MAX, where, for each cell $(x,y)$, we select the
maximal probability of the ensemble,
$\pi_{i}(x,y)=\max_{k}\left[\pi_{i}^{(k)}(x,y)\right]$.

  * •

The less cautious protocol is MIN, where we select the minimal probability of
the ensemble, $\pi_{i}(x,y)=\min_{k}\left[\pi_{i}^{(k)}(x,y)\right]$.

  * •

Finally, at the intermediate level, we define MEDIAN where $\pi_{i}(x,y)$ is
the median value among the $N_{\rm tr}$ sorted values of the ensemble.

We discuss in the next section how these different protocols perform on the
test set. In the Appendices we provide the details of the interpolation method
based on the 2D Navier-Cauchy equations and on the network architecture and
training procedure. Moreover we provide the catalogue of the earthquakes
employed in this study and all the Python codes for downloading and
manipulating GPS data as well as training and testing the neural network.

##  3 Results

Figure 2: Top: From left to right, surface displacement map of the mainshocks
$1,15,17$ (see Tab.(1)). This is the input of the neural network. Down: the
final prediction $\pi_{i}(x,y)$. This is the output of the neural network,
obtained with the MAX protocol.

The goal of this study is to use the GPS surface displacement map of the
mainshock $i$ to predict if, within $45$ days, at least an aftershock has
occurred inside each cell $(x,y)$ of area $25\rm{km}^{2}$. To convert the
continuous probability $\pi_{i}(x,y)$ in a deterministic event with
probability $0$ or $1$ we introduce a global threshold $\pi_{\rm th}$. Hence,
we call positive a cell with $\pi_{i}(x,y)>\pi_{\rm th}$. We divided the
positive cells in two categories:

  * (i)

the true positive are the positive cells where the prediction is correct. It
is convenient to normalize their number by the total number of cells where at
least an aftershock has occurred and obtain the true positive rate (TPR).

  * (ii)

Similarly, the false positive are the positive cells where the prediction is
wrong. Their number should be normalized with the total number of cells
without aftershocks to obtain the false positive rate (FPR).

To assess the quality of our classifier we use the receiver operating
characteristic (ROC) curve [Ian T. Jolliffe (2012), Zechar (2010)], obtained
by varying the parameter $\pi_{\rm th}$. The score associated to a ROC curve
is usually given by the area under the curve (AUC). For example, the ROC curve
of a random classifier is the line $\rm{FPR}=\rm{TPR}$ with an AUC of $0.5$;
an AUC above $0.75$ is considered a very good classifier. In Fig. 3 left the
ROC curves obtained from the test set are shown. All protocols outperform the
random classifier. The best protocol is the MAX, with an AUC of $0.7$. A
fingerprint of the robustness of our results is that the MEDIAN protocol has a
very similar ROC curve, proving the majority of the networks in the ensemble
tends to agree. Another fingerprint of such robustness is that the output
probabilities $\pi_{i}(x,y)$ after each mainshock are broadly distributed with
cells with a very low probabilities $\sim 0.2$ and cells with high
probabilities $\sim 0.8$. This means that the GPS displacement maps have a
great forecasting potential because in the case of non-informative input data
the output probabilities will have values narrow distributed around $\sim
0.5$.

To get further understanding in the inset of Fig. 3 left the AUC obtained for
each mainshock of the test set are shown. Most of them have a very good score,
but a fraction of them perform worst than a random classifier. A
representative mainshock belonging to this class is the 2017 Hachinohe (label
$7$) shown in Fig. 3 right. In this case the mainshock is offshore quite close
to the masked cells and thus in a region which is quite far from the GPS
stations. By inspecting Tab.1 we observe that a bad performance occurs either
in maps where the number of GPS stations used for the interpolation is small
(labels $2$, $9$, $10$, $19$) either in situations where is located close to
the masked cells (labels $7$ and $18$).

Figure 3: Left: the ROC curves obtained with the three protocols MAX, MEDIAN,
MIN, averaged over the $19$ mainshocks of the test set. Inset: the AUC scores
computed for each sample of the test set (MAX protocol). Right top: Surface
displacement map of the mainshock $7$ (see Tab.(1)). Right down: final
prediction $\pi(x,y)$ (MAX protocol). In this case the network fails to
produce a satisfying prediction. Even though 13 GPS stations were used to
construct the surface displacement map, most of the events happen in the sea,
where the signal quality is poor.

##  4 Conclusion

In conclusion our results show that this machine learning method based solely
on GPS data is very promising. The principal limitation of the performance
concerns the density of the GPS station which is high in Japan but much
smaller in other seismic regions such as California, Italy, Greece, etc.
Our study is based only on the displacement map of the day of the mainshock.
However, recent studies [Perfettini Avouac (2004), Perfettini Avouac (2007),
Perfettini . (2005)] have shown that large earthquakes are followed by an
increased afterslip in the regions where the aftershock rate is very high as
also confirmed by numerical models [Petrillo . (2020), Lippiello, Petrillo,
Landes Rosso (2019)]. It will be interesting to explore if the GPS data
available after the mainshock can be used to improve the aftershock
forecasting.
Furthermore, for concrete improvement of the current aftershock prediction
methods (ETAS models), one should consider a combination of past catalog data
(ETAS approach) and recent surface deformation (GPS data processed with a
CNN).

###### Acknowledgements.

V.M.S. thanks the valuable feedback from the conference ”1ères Journées du GDR
IDE ”Interaction Désordre Élasticité” - 2022” where this work in preparation
was presented. G.P. acknowledges support from MEXT Project for Seismology
TowArd Research innovation with Data of Earthquake (STAR-E Project), Grant
Number: JPJ010217. The GPS data used in the study is automatically downloaded
from the repositories of the Nevada Geodetic Laboratory [Blewitt . (2018)]
(http://geodesy.unr.edu/). We thank each independent work that processed the
data from the individual GPS stations employed. The earthquake catalog used in
this study is produced by the Japan Meteorological Agency, in cooperation with
the Ministry of Education, Culture, Sports, Science and Technology. The
catalog is based on seismic data provided by the National Research Institute
for Earth Science and Disaster Resilience, the Japan Meteorological Agency,
Hokkaido University, Hirosaki University, Tohoku University, the University of
Tokyo, Nagoya University, Kyoto University, Kochi University, Kyushu
University, Kagoshima University, the National Institute of Advanced
Industrial Science and Technology, the Geographical Survey Institute, Tokyo
Metropolis, Shizuoka Prefecture, Hot Springs Research Institute of Kanagawa
Prefecture, Yokohama City, and Japan Agency for Marine-Earth Science and
Technology.

##  Appendix A Computation of the surface displacement map via the 2D Navier-
Cauchy interpolation

GPS measurements are available only at the station positions
$\\{x_{s},y_{s}\\}_{s=1}^{N}$ ($x_{s}$ is the longitude, $y_{s}$ is the
latitude). Hence, at the day of the mainshock we dispose for each station of
its displacement vector
$\vec{v}(x_{s},y_{s})=(v_{e}(x_{s},y_{s}),v_{n}(x_{s},y_{s}))$. To construct
the map we need to compute the displacement vector of each cell $(x,y)$.
Following [Sandwell Wessel (2016)] we model the surface of the earth as $2D$
thin elastic sheet. The displacement satisfies the 2D Navier-Cauchy equations:

$\displaystyle\frac{2}{1-\nu}\partial_{x}^{2}v_{e}+\frac{2\nu}{1-\nu}\partial_{xy}^{2}v_{n}+\partial_{y}^{2}v_{e}+\partial_{xy}^{2}v_{n}$
$\displaystyle=-f_{e}\delta(x-x_{0})\delta(y-y_{0})$
$\displaystyle\partial_{xy}^{2}v_{e}+\partial_{x}^{2}v_{n}+\frac{2\nu}{1-\nu}\partial_{xy}^{2}v_{e}+\frac{2}{1-\nu}\partial_{y}^{2}v_{n}$
$\displaystyle=-f_{n}\delta(x-x_{0})\delta(y-y_{0})$ (2)

Here $\nu$ is the Poisson ratio which we set $1/2$ assuming the medium
isotropic, $f_{e}$ and $f_{n}$ are the point forces acting on the sheet at a
location $(x_{0},y_{0})$ (the forces are rescaled in the unit of the bulk
modulus which we assume constant in space). The solution of these linear
equations can be written as (see [Sandwell Wessel (2016)]):

$\displaystyle v_{e}(x,y)$ $\displaystyle=$ $\displaystyle
q(x-x_{0},y-y_{0})f_{e}+w(x-x_{0},y-y_{0})f_{n}$ $\displaystyle v_{n}(x,y)$
$\displaystyle=$ $\displaystyle
w(x-x_{0},y-y_{0})f_{e}+p(x-x_{0},y-y_{0})f_{n}$ (3)

where $q,p$ and $w$ are the Green’s functions:

$\displaystyle q(x,y)$ $\displaystyle=$
$\displaystyle(3-\nu)\ln(r)+(1+\nu)\frac{y^{2}}{r^{2}}$ (4) $\displaystyle
p(x,y)$ $\displaystyle=$
$\displaystyle(3-\nu)\ln(r)+(1+\nu)\frac{x^{2}}{r^{2}}$ (5) $\displaystyle
w(x,y)$ $\displaystyle=$ $\displaystyle-(1+\nu)\frac{xy}{r^{2}}$ (6)

with $r=\sqrt{x^{2}+y^{2}}$. To avoid the logarithmic divergence the Green’s
functions are regularized when $r=0$ by replacing $\ln(r)\to\ln(r+a)$ with the
constant $a=10$ km.

The interpolation strategy proposed in [Sandwell Wessel (2016)] is composed by
two steps:

  * •

First we assume that the point forces acting on the sheet are located at the
station positions. To determine their value, we invert the following $2N\times
2N$ linear system

$\begin{bmatrix}v_{e}(x_{s},y_{s})\\\
v_{n}(x_{s},y_{s})\end{bmatrix}=\underbrace{\begin{bmatrix}&q(x_{s}-x_{s^{\prime}},y_{s}-y_{s^{\prime}})&w(x_{s}-x_{s^{\prime}},y_{s}-y_{s^{\prime}})\\\
&w(x_{s}-x_{s^{\prime}},y_{s}-y_{s^{\prime}})&p(x_{s}-x_{s^{\prime}},y_{s}-y_{s^{\prime}})\end{bmatrix}}_{G}\begin{bmatrix}f_{e}(x_{s^{\prime}},y_{s^{\prime}})\\\
f_{n}(x_{s^{\prime}},y_{s^{\prime}})\end{bmatrix}$ (7)

To invert the linear system matrix we use SVD decomposition $G=U\Sigma V^{*}$
and keep the $80\%$ of the eigenvalues of the matrix $\Sigma$. This allows to
a smooth the final map and avoid singularities in the inversion.

  * •

Second, we use equation (A) with the forces computed in the first step to find
the displacement vector $(v_{e}(x,y),v_{n}(x,y))$ at the cell positions
$(x,y)$.

##  Appendix B Network architecture and training procedure

The building block of Convolutional Neural Networks (CNNs) is indeed the
convolution operation. Given an image $I(x,y)$, the convolution with a filter
of size $2k+1\times 2k+1$, denoted by $W_{ij}$, produces a new ’image’
$(W\star I)$:

$(W\star I)(x,y)=\sum_{i,j=-k}^{k}W_{i,j}I(x+i,y+j)$ (8)

In practice, we use mostly $k=2$, i.e. a filter of $5\times 5$. The parameters
$W_{i,j}$ are learned during training. As the convolution operation is linear,
one must add a non-linearity after the convolution operation in order to make
the network expressive. A popular choice is a ReLU (rectified linear unit)
function, which is nothing but $\textrm{ReLU}(z)=\max(z,0)$. In practice the
ReLU is applied, pixel by pixel, on the output of the convolution. A
convolution followed by a non linearity, $\textrm{ReLU}(W\star I)$, defines a
Convolution Block. A CNN consists in a series of Convolution Blocks stacked
one after the other.

As mentioned in the main text, we use the polar representation of surface
displacement maps as input, obtained as following:

$\displaystyle v_{r}(x,y)$ $\displaystyle=$
$\displaystyle\sqrt{v_{e}(x,y)^{2}+v_{n}(x,y)^{2}}$ $\displaystyle s(x,y)$
$\displaystyle=$ $\displaystyle v_{n}(x,y)/v_{r}(x,y)$ $\displaystyle c(x,y)$
$\displaystyle=$ $\displaystyle v_{e}(x,y)/v_{r}(x,y)$

$v_{r}(x,y)$ is thus the magnitude of the displacement and $s(x,y)$ and
$c(x,y)$ are the sine and cosines of the displacement vector on the $2D$
plane. To process the inputs, we use two CNNs made of 3 convolutional blocks
each. The radial input is fed to one, the angular to the other. The outputs
from the two networks are put together and fed to a final CNN with 3 blocks.
The last block in this case does not use a ReLU function but a sigmoid
$\sigma(z)=1/(1+e^{-z})$ that has domain on $[0,1]$ (i.e. represents a
probability). Given all these details, we can describe using three equations
our network:

$\displaystyle o_{r}(x,y)$ $\displaystyle=$
$\displaystyle\textrm{ReLU}(W_{3,r}\star\textrm{ReLU}(W_{2,r}\star\textrm{ReLU}(W_{1,r}\star
v_{r}(x,y))))$ (9) $\displaystyle o_{a}(x,y)$ $\displaystyle=$
$\displaystyle\textrm{ReLU}(W_{3,a}\star\textrm{ReLU}(W_{2,a}\star\textrm{ReLU}(W_{1,a}\star\\{s(x,y),c(x,y)\\})))$
(10) $\displaystyle\pi(x,y)$ $\displaystyle=$
$\displaystyle\sigma(W_{3}\star\textrm{ReLU}(W_{2}\star\textrm{ReLU}(W_{1}\star\\{o_{r}(x,y),o_{a}(x,y)\\})))$
(11)

The training of a neural network is nothing but the minimization of a Loss
function by means of gradient descent, or some refined version of this classic
optimization method, where the variables to be optimized are the weights of
the network (the coefficients of the filters $W$). The appropriate loss
function for a binary classification is the binary cross entropy loss:

$\displaystyle\mathcal{L}=-\sum_{x,y,i}A_{i}(x,y)\ln(\pi_{i}(x,y))+(1-A_{i}(x,y))\ln(1-\pi_{i}(x,y))$
(12)

where $i$ is the sample index (a sequence is a sample), $A_{i}(x,y)$ is the
binary label to predict, 1 if cell $(x,y)$ has an aftershock and 0 elsewhere,
and $\pi_{i}(x,y)$ is the predicted output probabilities (i.e. the whole
neural network output for the sample $i$). We mask the cells that are too far
from any station (further away than $d^{*}=86km$) by taking them out from the
sum. The whole architecture is summarized in Fig. 4.

Figure 4: Scheme for the network architecture. The two convolutional blocks on
the left are concatenated and passed to a last convolutional block that
produces the output probabilities.

We now provide details on how to search the best CNN architecture, i.e. the
procedure that led us to chose our machine learning model. A priori, numerous
architectures (different combinations of layers) can be considered, and the
goal is to pick whichever performs best. However, because of the scarcity of
data points, it’s not trivial to measure the notion of “works best”. A naive
approach would simply split the training set in a train and validation part,
using the classification score (as measured e.g. by the AUC) and loss function
values measured on the validation data as estimates of the performance of a
given architecture. However, this is very much not robust: given the small
number of train samples, and thus also of validation samples, this would
result in an essentially random choice, that would crucially depend on which
samples where picked for the validation set.
Thus, for each architecture or set of hyper-parameters we want to try, we
perform leave-one-out cross-validation: we build $N_{tr}$ train/validation
splits, for each one we train the model using $N_{tr}-1$ points and measure
the validation score on the left-out point. The average performance over all
splits is our measure of validation score, and the spread of scores between
splits is a measure of the robustness of the model: a large spread means the
model is not very robust.

Figure 5: Validation scores of all splits, for each event in the
train/validation set. Each event corresponds to an abscissa and is assigned a
random color for clarity. We note that for a given event (fixed abscissa), the
spread of scores of different models is small, with respect to the spread
between events.

In Fig. 5 we show that most of the randomness between splits comes from that
data points being different, not from the models. In our search of optimal
model, we tried to get both large validation score and small spread between
scores.

To reproduce the results in this work we refer to the code and data stored on
GitHub https://github.com/vicioms/gps_aftershocks_ml. Moreover in a README
file we provide the technical details of the training procedure.

##  Appendix C Catalogue of the test set

In the GitHub provided https://github.com/vicioms/gps_aftershocks_ml the full
list of mainshocks and aftershocks employed is available. For simplicity we
list in Table 1 the $19$ mainshocks of the test test. The blue color
corresponds to mainshocks for which the forecasting works well, the color red
is used for the others (see also the inset of Fig.2 left) .

Table 1: List of the $19$ earthquakes employed for the test. Index | Code Name | Data | Mw | Epicenter | #GPS Stations
---|---|---|---|---|---
1 | Kumamoto | 2016/04/14 | 6.5 | 32.74N 130.81E | 66
2 | Yilan | 2016/05/12 | 6.5 | 24.66N 121.92E | 1
3 | Yonakuni | 2016/06/24 | 6.2 | 23.50N 123.33E | 1
4 | Kurayoshi | 2016/10/21 | 6.6 | 35.38N 133.86E | 41
5 | Namie | 2016/11/22 | 7.4 | 37.35N 141.60E | 23
6 | Daigo | 2016/12/28 | 6.3 | 36.72N 140.57E | 46
7 | Hachinohe | 2017/09/27 | 6.1 | 40.27N 142.45E | 13
8 | Misawa | 2018/01/24 | 6.3 | 41.01N 142.45E | 6
9 | Hualien | 2018/02/04 | 6.5 | 24.12N 121.70E | 1
10 | Hualien | 2018/02/08 | 6.2 | 24.05N 121.76E | 2
11 | Matsue | 2018/04/09 | 6.1 | 35.18N 132.59E | 31
12 | Osaka | 2018/06/18 | 6.1 | 34.84N 135.62E | 75
13 | Chiba | 2018/07/07 | 6.0 | 35.16N 140.59E | 29
14 | Chitose | 2018/09/06 | 6.7 | 42.69N 142.01E | 31
15 | Nishinoomote | 2019/01/08 | 6.0 | 30.57N 131.16E | 15
16 | Miyazaki | 2019/05/10 | 6.3 | 31.80N 131.97E | 20
17 | Tsuruoka | 2019/06/18 | 6.7 | 38.61N 139.48E | 28
18 | Namie | 2019/08/04 | 6.4 | 37.71N 141.63E | 25
19 | Yilan | 2019/08/08 | 6.4 | 24.37N 121.87E | 3

## References

  * Blewitt . (2018) Blewitt2018Blewitt, G., Hammond, W. Kreemer, C.  201809\.  Harnessing the GPS Data Explosion for Interdisciplinary Science Harnessing the GPS data explosion for interdisciplinary science.  Eos99.  https://doi.org/10.1029/2018eo104623 10.1029/2018eo104623
  * Boore Joyner (1982) boore1982empiricalBoore, DM. Joyner, WB.  1982\.  The empirical prediction of ground motion The empirical prediction of ground motion.  Bulletin of the Seismological Society of America726BS43–S60.
  * De Arcangelis . (2018) de2018overlapDe Arcangelis, L., Godano, C. Lippiello, E.  2018\.  The overlap of aftershock coda waves and short-term postseismic forecasting The overlap of aftershock coda waves and short-term postseismic forecasting.  Journal of Geophysical Research: Solid Earth12375661–5674.
  * DeVries . (2018) DeVries2018DeVries, PMR., Viégas, F., Wattenberg, M. Meade, BJ.  201808\.  Deep learning of aftershock patterns following large earthquakes Deep learning of aftershock patterns following large earthquakes.  Nature5607720632–634.  https://doi.org/10.1038/s41586-018-0438-y 10.1038/s41586-018-0438-y
  * Dietterich (2000) dietterich2000ensembleDietterich, TG.  2000\.  Ensemble methods in machine learning Ensemble methods in machine learning.  Multiple Classifier Systems: First International Workshop, MCS 2000 Cagliari, Italy, June 21–23, 2000 Proceedings 1 Multiple classifier systems: First international workshop, mcs 2000 cagliari, italy, june 21–23, 2000 proceedings 1 ( 1–15).
  * Felzer Brodsky (2005) felzer2005testingFelzer, KR. Brodsky, EE.  2005\.  Testing the stress shadow hypothesis Testing the stress shadow hypothesis.  Journal of Geophysical Research: Solid Earth110B5.
  * Fukushima Tanaka (1990) fukushima1990newFukushima, Y. Tanaka, T.  1990\.  A new attenuation relation for peak horizontal acceleration of strong earthquake ground motion in Japan A new attenuation relation for peak horizontal acceleration of strong earthquake ground motion in japan.  Bulletin of the seismological Society of America804757–783.
  * Hainzl (20161) hainzl2016apparentHainzl, S.  20161\.  Apparent triggering function of aftershocks resulting from rate-dependent incompleteness of earthquake catalogs Apparent triggering function of aftershocks resulting from rate-dependent incompleteness of earthquake catalogs.  Journal of Geophysical Research: Solid Earth12196499–6509.
  * Hainzl (20162) hainzl2016rateHainzl, S.  20162\.  Rate-dependent incompleteness of earthquake catalogs Rate-dependent incompleteness of earthquake catalogs.  Seismological Research Letters872A337–344.
  * Hardebeck . (1998) hardebeck1998staticHardebeck, JL., Nazareth, JJ. Hauksson, E.  1998\.  The static stress change triggering model: Constraints from two southern California aftershock sequences The static stress change triggering model: Constraints from two southern california aftershock sequences.  Journal of Geophysical Research: Solid Earth103B1024427–24437.
  * Helmstetter Sornette (2003) HS03Helmstetter, A. Sornette, D.  2003\.  Importance of direct and indirect triggered seismicity in the ETAS model of seismicity Importance of direct and indirect triggered seismicity in the etas model of seismicity.  Geophysical Research Letters30.
  * Ian T. Jolliffe (2012) forecastVerIan T. Jolliffe, DBS.  2012\.  Forecast Verification: A Practitioner’s Guide in Atmospheric Science Forecast verification: A practitioner’s guide in atmospheric science.  Wiley.
  * Jacques . (1996) jacques1996seismicJacques, E., King, G., Tapponnier, P., Ruegg, J. Manighetti, I.  1996\.  Seismic activity triggered by stress changes after the 1978 events in the Asal Rift, Djibouti Seismic activity triggered by stress changes after the 1978 events in the asal rift, djibouti.  Geophysical research letters23182481–2484.
  * Kagan (2004) kagan2004shortKagan, YY.  2004\.  Short-term properties of earthquake catalogs and models of earthquake source Short-term properties of earthquake catalogs and models of earthquake source.  Bulletin of the Seismological Society of America9441207–1228.
  * Khawaja . (2023) kh2023Khawaja, AM., Hainzl, S., Schorlemmer, D., Iturrieta, P., Bayona, JA., Savran, WH.Marzocchi, W.  202301\.  Statistical power of spatial earthquake forecast tests Statistical power of spatial earthquake forecast tests.  Geophysical Journal International23332053-2066.  https://doi.org/10.1093/gji/ggad030 10.1093/gji/ggad030
  * King . (1994) king1994staticKing, GC., Stein, RS. Lin, J.  1994\.  Static stress changes and the triggering of earthquakes Static stress changes and the triggering of earthquakes.  Bulletin of the Seismological Society of America843935–953.
  * Krizhevsky . (2012) NIPS2012_c399862dKrizhevsky, A., Sutskever, I. Hinton, GE.  2012\.  ImageNet Classification with Deep Convolutional Neural Networks Imagenet classification with deep convolutional neural networks.  F. Pereira, C. Burges, L. Bottou K. Weinberger (), Advances in Neural Information Processing Systems Advances in neural information processing systems ( 25).  Curran Associates, Inc.  https://proceedings.neurips.cc/paper_files/paper/2012/file/c399862d3b9d6b76c8436e924a68c45b-Paper.pdf
  * Lecun . (1998) 726791Lecun, Y., Bottou, L., Bengio, Y. Haffner, P.  1998\.  Gradient-based learning applied to document recognition Gradient-based learning applied to document recognition.  Proceedings of the IEEE86112278-2324.  10.1109/5.726791
  * Lippiello . (2014) Lippiello2014Lippiello, E., Giacco, F., d. Arcangelis, L., Marzocchi, W. Godano, C.  201403\.  Parameter Estimation in the ETAS Model: Approximations and Novel Methods Parameter estimation in the ETAS model: Approximations and novel methods.  Bulletin of the Seismological Society of America1042985–994.  https://doi.org/10.1785/0120130148 10.1785/0120130148
  * Lippiello, Petrillo, Godano . (2019) lippiello2019forecastingLippiello, E., Petrillo, G., Godano, C., Tramelli, A., Papadimitriou, E. Karakostas, V.  2019\.  Forecasting of the first hour aftershocks by means of the perceived magnitude Forecasting of the first hour aftershocks by means of the perceived magnitude.  Nature communications1012953.
  * Lippiello, Petrillo, Landes Rosso (2019) lippiello2019faultLippiello, E., Petrillo, G., Landes, F. Rosso, A.  2019\.  Fault Heterogeneity and the Connection between Aftershocks and AfterslipShort Note Fault heterogeneity and the connection between aftershocks and afterslipshort note.  Bulletin of the Seismological Society of America10931156–1163.
  * Lombardi (2015) lombardi2015Lombardi, A.  201502\.  Estimation of the parameters of ETAS models by Simulated Annealing Estimation of the parameters of etas models by simulated annealing.  Scientific reports58417.  10.1038/srep08417
  * Mallman Zoback (2007) mallman2007assessingMallman, EP. Zoback, MD.  2007\.  Assessing elastic Coulomb stress transfer models using seismicity rates in southern California and southwestern Japan Assessing elastic coulomb stress transfer models using seismicity rates in southern california and southwestern japan.  Journal of Geophysical Research: Solid Earth112B3.
  * Molkenthin . (2022) molkenthin2022Molkenthin, C., Donner, C., Reich, S., Zöller, G., Hainzl, S., Holschneider, M. Opper, M.  202204\.  GP-ETAS: semiparametric Bayesian inference for the spatio-temporal epidemic type aftershock sequence model Gp-etas: semiparametric bayesian inference for the spatio-temporal epidemic type aftershock sequence model.  Statistics and Computing32.  10.1007/s11222-022-10085-3
  * Nostro . (1997) nostro1997staticNostro, C., Cocco, M. Belardinelli, ME.  1997\.  Static stress changes in extensional regimes: an application to southern Apennines (Italy) Static stress changes in extensional regimes: an application to southern apennines (italy).  Bulletin of the Seismological Society of America871234–248.
  * Ogata (1988) ogata1988Ogata, Y.  1988\.  Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes Statistical models for earthquake occurrences and residual analysis for point processes.  Journal of the American Statistical Association834019–27.
  * Ogata (1998) ogata1998Ogata, Y.  1998\.  Space-time point-process models for earthquake occurrences Space-time point-process models for earthquake occurrences.  Annals of the Institute of Statistical Mathematics502379–402.
  * Ogata Zhuang (2006) OZ98Ogata, Y. Zhuang, J.  2006\.  Space–time ETAS models and an improved extension Space–time etas models and an improved extension.  Tectonophysics413113-23.
  * Omi . (2013) omi2013forecastingOmi, T., Ogata, Y., Hirata, Y. Aihara, K.  2013\.  Forecasting large aftershocks within one day after the main shock Forecasting large aftershocks within one day after the main shock.  Scientific reports311–7.
  * Omi . (2015) omi2015intermediateOmi, T., Ogata, Y., Hirata, Y. Aihara, K.  2015\.  Intermediate-term forecasting of aftershocks from an early aftershock sequence: Bayesian and ensemble forecasting approaches Intermediate-term forecasting of aftershocks from an early aftershock sequence: Bayesian and ensemble forecasting approaches.  Journal of Geophysical Research: Solid Earth12042561–2578.
  * Omi . (2019) omi2019implementationOmi, T., Ogata, Y., Shiomi, K., Enescu, B., Sawazaki, K. Aihara, K.  2019\.  Implementation of a real-time system for automatic aftershock forecasting in Japan Implementation of a real-time system for automatic aftershock forecasting in japan.  Seismological Research Letters901242–250.
  * Parsons . (1999) parsons1999stressParsons, T., Stein, RS., Simpson, RW. Reasenberg, PA.  1999\.  Stress sensitivity of fault seismicity: A comparison between limited-offset oblique and major strike-slip faults Stress sensitivity of fault seismicity: A comparison between limited-offset oblique and major strike-slip faults.  Journal of Geophysical Research: Solid Earth104B920183–20202.
  * Peng . (2007) peng2007seismicityPeng, Z., Vidale, JE., Ishii, M. Helmstetter, A.  2007\.  Seismicity rate immediately before and after main shock rupture from high-frequency waveforms in Japan Seismicity rate immediately before and after main shock rupture from high-frequency waveforms in japan.  Journal of Geophysical Research: Solid Earth112B3.
  * Perfettini Avouac (2004) perfettini2004postseismicPerfettini, H. Avouac, JP.  2004\.  Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake, Taiwan Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the chi-chi earthquake, taiwan.  Journal of Geophysical Research: Solid Earth109B2.
  * Perfettini Avouac (2007) perf2007Perfettini, H. Avouac, JP.  2007\.  Modeling afterslip and aftershocks following the 1992 Landers earthquake Modeling afterslip and aftershocks following the 1992 landers earthquake.  Journal of Geophysical Research: Solid Earth112B7.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2006JB004399 https://doi.org/10.1029/2006JB004399
  * Perfettini . (2005) perf2005Perfettini, H., Avouac, JP. Ruegg, JC.  2005\.  Geodetic displacements and aftershocks following the 2001 Mw = 8.4 Peru earthquake: Implications for the mechanics of the earthquake cycle along subduction zones Geodetic displacements and aftershocks following the 2001 mw = 8.4 peru earthquake: Implications for the mechanics of the earthquake cycle along subduction zones.  Journal of Geophysical Research: Solid Earth110B9.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2004JB003522 https://doi.org/10.1029/2004JB003522
  * Petrillo . (2020) petrillo2020influencePetrillo, G., Lippiello, E., Landes, FP. Rosso, A.  2020\.  The influence of the brittle-ductile transition zone on aftershock and foreshock occurrence The influence of the brittle-ductile transition zone on aftershock and foreshock occurrence.  Nature Communications1113010.
  * Petrillo . (2022) petrillo2022testingPetrillo, G., Rosso, A. Lippiello, E.  2022\.  Testing of the seismic gap hypothesis in a model with realistic earthquake statistics Testing of the seismic gap hypothesis in a model with realistic earthquake statistics.  Journal of Geophysical Research: Solid Earth1276e2021JB023542.
  * Reasenberg Simpson (1992) reasenberg1992responseReasenberg, PA. Simpson, RW.  1992\.  Response of regional seismicity to the static stress change produced by the Loma Prieta earthquake Response of regional seismicity to the static stress change produced by the loma prieta earthquake.  Science25550521687–1690.
  * Ross (2021) ross2021Ross, GJ.  202105\.  Bayesian Estimation of the ETAS Model for Earthquake Occurrences Bayesian Estimation of the ETAS Model for Earthquake Occurrences.  Bulletin of the Seismological Society of America11131473-1480.  https://doi.org/10.1785/0120200198 10.1785/0120200198
  * Sandwell Wessel (2016) sandwell2016interpolationSandwell, DT. Wessel, P.  2016\.  Interpolation of 2-D vector data using constraints from elasticity Interpolation of 2-d vector data using constraints from elasticity.  Geophysical Research Letters432010–703.
  * Savran . (2022) savran2022Savran, WH., Bayona, JA., Iturrieta, P., Asim, KM., Bao, H., Bayliss, K.Werner, MJ.  202207\.  pyCSEP: A Python Toolkit for Earthquake Forecast Developers pyCSEP: A Python Toolkit for Earthquake Forecast Developers.  Seismological Research Letters9352858-2870.  https://doi.org/10.1785/0220220033 10.1785/0220220033
  * Savran . (2021) csepSavran, WH., Werner, MJ., Schorlemmer, D. Maechling, PJ.  202104\.  pyCSEP - Tools for Earthquake Forecast Developers. pyCSEP \- Tools for Earthquake Forecast Developers.
  * Shen . (2015) shen2015optimalShen, ZK., Wang, M., Zeng, Y. Wang, F.  2015\.  Optimal interpolation of spatially discretized geodetic data Optimal interpolation of spatially discretized geodetic data.  Bulletin of the Seismological Society of America10542117–2127.
  * Toda . (1998) toda1998stressToda, S., Stein, RS., Reasenberg, PA., Dieterich, JH. Yoshida, A.  1998\.  Stress transferred by the 1995 Mw= 6.9 Kobe, Japan, shock: Effect on aftershocks and future earthquake probabilities Stress transferred by the 1995 mw= 6.9 kobe, japan, shock: Effect on aftershocks and future earthquake probabilities.  Journal of Geophysical Research: Solid Earth103B1024543–24565.
  * Utsu . (1995) utsu1995centenaryUtsu, T., Ogata, Y. .  1995\.  The centenary of the Omori formula for a decay law of aftershock activity The centenary of the omori formula for a decay law of aftershock activity.  Journal of Physics of the Earth4311–33.
  * Zechar (2010) zechar2010evaluatingZechar, JD.  2010\.  Evaluating earthquake predictions and earthquake forecasts: A guide for students and new researchers Evaluating earthquake predictions and earthquake forecasts: A guide for students and new researchers.