Earthquake Location and Magnitude Estimation
with Graph Neural Networks

The basic computational unit of a standard GNN is the generic message passing layer. The message passing layer captures the powerful notion of graph convolution (e.g., [6, 1]), wherein the latent states of nodes on graphs are updated based on (learned) mappings from the latent states of neighboring nodes. The operator is generically defined as $\mathcal{G}(\boldsymbol{h}_{\mathcal{A}},\mathcal{E}_{\mathcal{A}^{\prime}\leftarrow\mathcal{A}}):\mathcal{A}\longrightarrow\mathcal{A}^{\prime}$, and is given as

In the forward call (1), on a per-node basis $(i)$, the latent features of each neighboring node $j\in\mathcal{N}(i)$ are collected, and transformed by the learnable Fully Connected Network (FCN) layer, $\phi^{msg}$, prior to being pooled over the entire neighborhood, $\mathcal{N}(i)$, where the pooling is over the node axis ($POOL:\mathbb{R}^{|\mathcal{N}(i)|\times K}\longrightarrow\mathbb{R}^{K}$), and is any global pooling operator such as sum-, mean, or max-pool. The pooled aggregate message is then concatenated with the self state of node $\boldsymbol{h}_{i}$, and transformed by an additional learnable FCN, $\phi^{agg}$, resulting in the updated state on node $i$. Inside the message passing layer $\phi^{msg}$, the optional term of $\boldsymbol{e}_{ij}$ represents edge data between two nodes $(i,j)$, such as Cartesian offsets of node coordinates (if nodes represent spatially located points), and $\boldsymbol{z}=\text{POOL}\{\phi^{glb}(\boldsymbol{h}_{j})\}$ is an optional ‘global’ summary feature, that can be extracted from the entire graph at each step, using an additional learnable FCN, $\phi^{glb}$. By combining these facets, the overall capability of the graph convolution cell is to learn local feature transformations that work well for a wide range of nodes, while also exploiting local information transfer through the adjacency matrix.

In our method there is one component (or graph) that specifies the station set, $\mathcal{S}$, and another that represents the model space, $\mathcal{H}$, where we consider either source region $\mathcal{H}=\mathcal{X}$ or magnitude axis $\mathcal{H}=\mathcal{M}$. For either model space, we simply need a (regular or irregular) grid that spans the region of interest. For $\mathcal{X}$, we take a collection of $\sim$100’s of sparsely distributed points covering the source region of interest, and for the magnitude axis we take a fine spacing of nodes in magnitude (e.g., [-3,7] M with 0.1 or 0.25 M increment).

Rather than leave the station set and model spaces as unordered sets, we give these sets structure by interpreting them as graphs. Specifically, $(\mathcal{S},\mathcal{E}_{\mathcal{S}})$ is the station graph, $(\mathcal{X},\mathcal{E}_{\mathcal{X}})$ is the source graph, and $(\mathcal{M},\mathcal{E}_{\mathcal{M}})$ is the magnitude graph, where we define edge sets, $\mathcal{E}_{\mathcal{S}},\mathcal{E}_{\mathcal{X}},\mathcal{E}_{\mathcal{M}}$, or equivalently adjacency matrices, that specify which elements in the respective sets are linked. While many edge construction schemes are possible [3], GNN’s often work well for a wide range of adjacency matrix choices [20] because the GNN can learn how to use whichever distribution of graphs it is trained over. For point cloud data, such as represented in $\mathcal{S}$ and $\mathcal{X}$, a natural choice of edge construction scheme is to use K-nearest-neighbor (K-NN) graphs or $\epsilon-$distance graphs [3], since these favor connecting nearby elements to one another. For $\mathcal{M}$, while it is simply a linear grid of nodes, $K-NN$ and $\epsilon-$distance graphs effectively represent 1D convolution kernels over the $\mathcal{M}$ axis, and are hence also a natural choice.

For our purposes we use $K-NN$ instead of $\epsilon-$distance graphs for all of $\mathcal{S}$, $\mathcal{X}$, $\mathcal{M}$, since these limit the incoming degree of all nodes to $K$, whereas $\epsilon-$distance graphs can have substantially different degrees across different nodes of the graph. While we did not optimize these choices fully, modest trial and error led to

We still need to define one essential graph type: the Cartesian product graph of $\mathcal{S}\times\mathcal{H}$, which has as nodes all pairs of station-model ($\boldsymbol{s}$, $\boldsymbol{h}$) space coordinates. As discussed earlier, and as highlighted in the next section, it is most naturally all pairs of $(\boldsymbol{s},\boldsymbol{h})\in\mathcal{S}\times\mathcal{H}$ stations, and model space coordinates, that we can use to parse the data. The graphs already defined can be used directly to construct the Cartesian product graph with two-edge types, given by

The input features we pass to the GNN are defined for arbitrary input pick data sets, where for location, it includes pick times (with assigned phase types), and for magnitude it includes amplitude measurements (with assigned phase types). Each input can be any given station set, $\mathcal{S}$, and individual stations in this set can still have missing picks or false picks. For the source location problem our input tensor is defined analogously to the ‘pre-stack’ back-projection metric (e.g., [14]), which, for each spatial coordinate $\boldsymbol{x}\in\mathcal{X}$ is the measure of how close the nearest arrival, $\tau_{i}^{k}\in\mathcal{D}_{i}$, of phase type $k$ on station $\boldsymbol{s}_{i}\in\mathcal{S}$, is to the theoretical arrival time from that source coordinate. Hence, we compute

For the magnitude input, a similar input scheme is defined. We compute

The GNN architecture we propose is designed to be a simple, effective way to map the very high-dimensional input tensors $\boldsymbol{h}_{k}^{\mathcal{X}}(\boldsymbol{s}_{i},\boldsymbol{x})$, $\boldsymbol{h}_{k}^{\mathcal{M}}(\boldsymbol{s}_{i},\boldsymbol{m})$, which are feature vectors measured on all pairs of station-model space points, to explicit predictions of source likelihood $p(\boldsymbol{x})\subset\Omega_{\mathcal{X}}$, and magnitude likelihood $p(\boldsymbol{m})\subset\Omega_{\mathcal{M}}$, at any query points within the full continuous Euclidean model space, $\Omega_{\mathcal{H}}$. Each input can also be for a different number of stations, and a different station graph.

To achieve this overall mapping, we first (i) apply repeated graph convolutions on the Cartesian product $\mathcal{S}\times\mathcal{H}$, then (ii) embed into the model space by stacking over stations (in the latent space, $\mathcal{S}\times\mathcal{H}\rightarrow\mathcal{H}$). Then, (iii) we repeat graph convolutions on $\mathcal{H}$, and finally (iv) we make predictions of likelihood at arbitrary coordinates $p(\boldsymbol{h}_{q})\subset\Omega_{\mathcal{H}}$, by using local masked-attention (i.e., effectively graph-based interpolation [19]) of the queries $\boldsymbol{h}_{q}$ with respect to the fixed nodes in $\mathcal{H}$. This design is straightforward, and in a natural, physically motivated way, transforms the raw data measured over arbitrary station networks to predictions over arbitrary model space domains. The architecture is given in more detail in Algorithm (1).

To train the models, we generate a diverse set of synthetic data, for the assumed physical travel time and local magnitude scale models. In either case, for a batch of 150 samples, for each sample, we (i) sample an arbitrary subset of 50 - 450 stations of the NC network, (ii) pick a random source coordinate uniformly over the domain (and random magnitude uniformly over the axis), (iii) compute theoretical arrival times and amplitude measurements, (iv) remove arrivals of source-receiver distances greater than $d+\epsilon$ with $d\sim$ Uniform [10, 1000] km per event, and $\epsilon\sim\mathcal{N}$(0, 30 km) per event and per station, (v) add per-pick Laplacian noise of std. 1$-$5$\%$ the quantity for travel time, and 25$-$250$\%$ the quantity for amplitude (hence, uncertainties increase proportional to the observed quantity). Finally, we (vi) randomly delete 10$-$80$\%$ of arrivals, and randomly corrupt 10$-$30$\%$ of the remaining arrivals to have anomalous (uniformly random) arrival time and amplitude measurements. These simple steps can produce complex input datasets (including strong noise and anomalous picks), for variable seismic networks, with data that are still physically tied to the assumed travel time and magnitude scale models.

For the labels, we represent each target solution as a continuous Gaussian, wherein the target spatial prediction is a multivariate Gaussian of fixed width $\sigma_{\mathcal{X}}$ (= 25 km), localized at the correct source coordinate, and similarly, the target magnitude prediction is a 1D Gaussian, localized on the true magnitude, with width $\sigma_{\mathcal{M}}$ (= 0.5 M). If a particular input sample has $\leq 3$ non-corrupted arrivals following the synthetic scheme above, the labels are set to all zero for that sample. Models are implemented in PyTorch Geometric [5], and we train by computing L2-norm losses of the output with respect to the target over the batch, and optimizing with the Adam optimizer [10] for $\sim$10,000 update steps with a learning rate of $10^{-3}$. In $\sim$10’s of training runs, the GNN performed similarly and converged in training at similar rates, hence indicating that training was stable.