Deep learning · Part A: from a neuron to a U-Net picker

A seismic network records tens of thousands of earthquakes a year at the Geysers, and every one of them needs its P and S arrival times picked on every station before it can be located. An analyst picks a few hundred records a day. How can a computer learn to pick them as an analyst does — and when it has learned, what has it learned?

Every method in this course so far has fitted a model that we wrote down: a straight line for the magnitude–frequency relation (Session 1), travel times through a layered Earth (Session 2), the radiation pattern of a double couple (Session 3). A neural network is different. It is a very flexible function with many adjustable numbers, called parameters, and nobody writes down what they should be. They are learned from examples: thousands of seismograms, each with the arrival times an analyst picked. Fitting a network is still what we have done all along — choose the parameters that minimise a loss — but with so many parameters that two new things are needed: a cheap way to compute the gradient, and a way to tell learning from memorising.

This notebook builds a phase picker from the ground up, one idea at a time, and trains it on real Geysers seismograms. The picker is a simplified copy of PhaseNet (Zhu & Beroza 2018), the network our group built and that much of the field now uses, with the same names and the same design.

The questions this notebook answers

1. What is a neural network, and what can it represent that a linear model cannot? (A.1)
2. Where does the gradient come from, and why is it cheap even with millions of parameters? (A.2)
3. How is a network trained? (A.3)
4. How does a network say "noise, P or S" at a sample, and how is that scored? (A.4)
5. What are the Geysers seismograms and their labels, and how must they be split to test a picker honestly? (A.5)
6. Why use convolutions for a seismogram, and how far does a network see? (A.6–A.7)
7. How does a network output a probability at every sample, and how is it trained without memorising? (A.8–A.9)
8. How good are its picks, what do incomplete labels teach it, and does a probability of 0.5 mean a 50 % chance? (A.10–A.12)

Part B runs the published PhaseNet on Geysers earthquakes and compares it with the network trained here.

A.1 · A network is a function with adjustable parts

Start with a curve every seismologist knows: the first-arrival time of P against distance, for a source at the surface of a crust with two layers — a slow layer, 3 km thick, over a faster half-space. Session 2, Location & relocation, computed such curves. Near the source the first arrival is the direct wave, which travels in the top layer, so its time grows as distance divided by the top-layer speed $v_1$. Beyond the crossover distance the first arrival is the head wave, which runs along the top of the faster layer at speed $v_2$; its time grows more slowly with distance. The curve is two straight lines that meet at a kink.

We make 200 noisy observations of this curve and ask the simplest model we know, a straight line fitted by least squares (Session 1), to describe it.

crossover distance 11.7 km; head-wave intercept 1.39 s
straight-line fit: rms misfit 0.257 s against 0.05 s of noise

The line misses by several times the noise: no single slope fits both branches. We could write the right model down — we just did — but suppose we did not know the physics. We want a model that can bend where the data bend.

A neuron, and a layer of them

The building block is the neuron. It takes an input $x$, multiplies it by a weight $w$, adds a bias $b$, and passes the result through a fixed non-linear function called the activation. The activation used almost everywhere is the rectified linear unit, ReLU:

$$\mathrm{relu}(z) = \max(0, z).$$

It is zero for negative $z$ and equal to $z$ for positive $z$, so a neuron's output, $\mathrm{relu}(wx + b)$, is flat up to $x = -b/w$ and a straight line after it: one kink, at a place set by the bias. A hidden layer is $H$ neurons side by side, all reading the same input. The network's output adds their outputs, each with its own weight $a_j$, plus a final bias $c$:

$$\hat t(x) = c + \sum_{j=1}^{H} a_j\, \mathrm{relu}(w_j x + b_j).$$

This is a neural network with one hidden layer, also called a multilayer perceptron (MLP). Its parameters are the $3H + 1$ numbers $w_j, b_j, a_j, c$. Because each neuron contributes one kink, the network is a piecewise-linear function with up to $H$ kinks. The word "hidden" only means that the layer's outputs are not observed; the network computes them on the way.

Our curve has one kink, so two neurons are enough. One neuron, $\mathrm{relu}(x)$, gives the slope $1/v_1$ from the origin; a second, starting at the crossover distance, subtracts the difference in slope. We can set the six parameters by hand from the physics:

two neurons set by hand: rms misfit 0.048 s

Two neurons fit the curve to the noise level. Their parameters are physics: $a_1 = 1/v_1$ is the slowness of the top layer, $a_1 + a_2 = 1/v_2$ the slowness of the half-space, and $-b_2/w_2$ the crossover distance. Add more neurons and the network can follow more kinks; with enough of them, a one-hidden-layer network can approximate any continuous curve on an interval as closely as we want. That flexibility is the whole point — and it is also the danger, because a network with many kinks can follow the noise as well as the signal.

What we cannot do for a seismogram is set the parameters by hand. PhaseNet has about 270,000 of them, and nobody knows what any one should be. They have to be learned: found by minimising a loss, exactly as least squares found the line. The loss here is still the mean squared misfit, but the model is no longer linear in its parameters (they sit inside the ReLU), so there is no normal equation to solve. We need the gradient of the loss with respect to every parameter, and a way to walk downhill. A.2 computes the gradient; A.3 walks.

Optional: why stack layers

A one-hidden-layer network can represent any continuous function, but it may need very many neurons to do so. Feeding one layer's outputs into another layer of neurons composes kinks: a second layer can fold the piecewise-linear output of the first, and the number of linear pieces can grow exponentially with the number of layers while the number of parameters grows only linearly. A deep network therefore represents complicated functions with far fewer parameters than a wide, shallow one. This is the usual argument for depth (Prince, Understanding Deep Learning, ch. 4); in practice depth also builds features in stages, which A.7 makes concrete for seismograms.

A.2 · The gradient by the chain rule: backpropagation

Take the smallest network that shows everything: one neuron and a squared misfit on one observation $(x, t)$. The computation runs forward in four steps, each a simple operation whose derivative we know:

1. $z = w x + b$, whose derivatives are $\partial z/\partial w = x$ and $\partial z/\partial b = 1$;
2. $u = \mathrm{relu}(z)$, with $\partial u/\partial z = 1$ if $z > 0$ and $0$ otherwise;
3. $\hat t = a u + c$, with $\partial \hat t/\partial a = u$, $\partial \hat t/\partial u = a$ and $\partial \hat t/\partial c = 1$;
4. $L = (\hat t - t)^2$, with $\partial L/\partial \hat t = 2(\hat t - t)$.

Walking the graph backwards

Such a chain of operations is a computational graph. We want $\partial L/\partial w$ and the other three derivatives. The chain rule says that the derivative of $L$ with respect to anything is the derivative of $L$ with respect to the next step, times the derivative of that step with respect to the thing. So we start at the end, where $\partial L/\partial L = 1$, and walk backwards through the graph, multiplying by one local derivative per step:

$$\frac{\partial L}{\partial \hat t} = 2(\hat t - t), \qquad \frac{\partial L}{\partial a} = \frac{\partial L}{\partial \hat t}\,u, \qquad \frac{\partial L}{\partial c} = \frac{\partial L}{\partial \hat t},$$

$$\frac{\partial L}{\partial u} = \frac{\partial L}{\partial \hat t}\,a, \qquad \frac{\partial L}{\partial z} = \frac{\partial L}{\partial u}\,\mathbb{1}[z > 0], \qquad \frac{\partial L}{\partial w} = \frac{\partial L}{\partial z}\,x, \qquad \frac{\partial L}{\partial b} = \frac{\partial L}{\partial z}.$$

Each line reuses the one before it: the derivative arriving from later in the graph (the upstream gradient) times the local derivative of one step. That is all backpropagation is — the chain rule, evaluated from the loss backwards, reusing every intermediate product. We do it with numbers, for an observation at 20 km.

forward:  z = 14.00, u = 14.00, t_hat = 4.700 s, t = 4.726 s, L = 0.0007
by hand:  dL/dw = -0.3088, dL/db = -0.0154, dL/da = -0.7205, dL/dc = -0.0515

Two checks. The first is a finite difference: nudge one parameter by a tiny amount and see how much the loss changes. It needs one extra forward pass per parameter, which is why nobody trains a network this way, but it is the honest test of a gradient. The second is automatic differentiation in PyTorch, the library we use from here on. PyTorch records every operation as it runs the forward pass — it builds the computational graph for us — and loss.backward() walks it backwards exactly as we just did by hand.

finite difference: -0.3088, -0.0154, -0.7205, -0.0515
PyTorch autograd:  -0.3088, -0.0154, -0.7205, -0.0515

All three agree. Notice what the backward pass cost: one multiplication per step, the same work as the forward pass. That holds for any network, however large. The gradient with respect to every parameter costs about as much as one more evaluation of the network. For PhaseNet's 270,000 parameters a finite difference would need 270,000 extra forward passes; backpropagation needs one backward pass.

This is why networks are trained by gradient descent rather than by the Gauss–Newton and Levenberg–Marquardt methods of Session 2. Those methods need the full Jacobian, the derivative of every residual with respect to every parameter: a matrix of (number of data) × (number of parameters) entries, here millions × hundreds of thousands. Backpropagation never forms it. It returns only the gradient of one number, the loss — which is all gradient descent needs.

Optional: reverse mode, forward mode and the adjoint

The chain rule can be evaluated in either direction. Going forwards from one parameter gives the derivative of every output with respect to that one parameter (forward mode), so a full gradient needs one pass per parameter. Going backwards from one output gives the derivative of that output with respect to every parameter (reverse mode), so a full gradient needs one pass per output. A loss is a single number, so reverse mode wins by a factor of the number of parameters. The price is memory: every intermediate value of the forward pass must be stored until the backward pass uses it. Geophysicists met reverse mode long before deep learning, as the adjoint-state method used in waveform and travel-time tomography: one forward simulation and one adjoint (backward) simulation give the gradient of a misfit with respect to every cell of a velocity model. Session 8, Tomography, returns to it.

A.3 · Training: walking downhill

With the gradient in hand, gradient descent repeats one update: move every parameter a small step against its gradient,

$$\theta \leftarrow \theta - \eta\, \frac{\partial L}{\partial \theta},$$

where $\theta$ stands for all the parameters together and $\eta$, the learning rate, sets the step. Each pass of update is one iteration. The loss $L$ is now the mean squared misfit over all 200 observations.

Unlike A.1 we give the network no help: eight neurons, with weights drawn at random. One preparation comes first, and it matters more than it looks. Distances run from 0 to 40 km and times from 0 to 8 s, while random starting weights are of order one. We standardise both: subtract the mean and divide by the standard deviation, so input and target are each of order one and centred on zero. Every network in this notebook is fed standardised numbers; for seismograms (A.5) that means dividing each window by its own standard deviation.

The learning rate is the first thing to get right. We try four values on the same starting weights.

learning rate 0.01: rms misfit after 2000 steps 0.150 s
learning rate 0.1 : rms misfit after 2000 steps 0.048 s
learning rate 1.0 : rms misfit after 2000 steps 7.131 s

learning rate 2.0 : rms misfit after 2000 steps diverged

Too small a step (0.01) and the network is still far from the data after 2,000 iterations. A step ten times larger (0.1) reaches the noise level. At 1.0 the first steps overshoot so far that every neuron's input becomes negative for all the data: the neurons die and pass no gradient, and what is left is a single constant, the output bias, which each step throws from one side of the data to the other, so the misfit never falls. At 2.0 the steps grow without bound and the loss diverges. There is no formula for the right learning rate — it is found by trying, which is why the first thing to do with any new network is to watch its loss curve.

Mini-batches: stochastic gradient descent

Computing the gradient over all the data at every step is affordable for 200 points but not for a million seismograms. Stochastic gradient descent (SGD) estimates the gradient from a small random mini-batch of examples at each step. The estimate is noisy, but it is unbiased and costs a fraction of the full one, so far more steps fit in the same time. One pass through all the training examples is an epoch.

Momentum and Adam

Two refinements are used almost everywhere. Momentum keeps a running average of recent gradients and steps along it, so the noise of the mini-batches averages out and the steps grow along directions where the gradient keeps the same sign. Adam (Kingma & Ba 2015) adds a second running average, of the squared gradient, and divides each parameter's step by its square root, so every parameter moves at a similar rate whatever the scale of its gradient. Adam's default learning rate, $10^{-3}$, works for a wide range of networks, which is why it is the usual first choice; PhaseNet is trained with it. We compare the three on mini-batches of 16 observations. One run can be lucky, so each curve is the median over ten random starts.

SGD               : median rms misfit after 1500 steps 0.189 s; never below 0.06 s

SGD with momentum : median rms misfit after 1500 steps 0.053 s; below 0.06 s after 1205 steps

Adam              : median rms misfit after 1500 steps 0.054 s; below 0.06 s after 786 steps

With the same learning rate, momentum and Adam reach the noise level in a fraction of the steps plain SGD needs. The curves are jagged because every step sees a different mini-batch.

Where it starts matters

The loss of a network has many valleys, as the likelihood of a Gaussian mixture had many local optima in Session 4, Clustering. Training from different random starting weights can end in different places. With ReLU neurons there is a specific way to fail: a neuron whose input $w_j x + b_j$ is negative for every observation outputs zero everywhere, so its gradient is zero and it never recovers — a dead neuron. With only two neurons, losing one leaves a straight line. We train the two-neuron network of A.1, and an eight-neuron one, from ten random starts each.

2 neurons, rms misfit from 10 random starts (s): 0.26, 0.24, 0.26, 0.05, 2.14, 0.05, 0.26, 0.05, 0.23, 0.26  -> 3 of 10 at the noise level

8 neurons, rms misfit from 10 random starts (s): 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05  -> 10 of 10 at the noise level

Two neurons are exactly enough to represent the curve, yet most random starts fail: they end on the straight-line misfit of A.1 (about 0.26 s), a poorer valley, or at a constant because both neurons died. Eight neurons succeed far more often, although six of them are not needed to represent the curve. Two things make failures rare in practice. Networks are given more neurons than strictly needed, so a few dead ones or a poor start do not matter, and their starting weights are drawn at a carefully chosen scale: He initialisation draws each weight with variance $2/n_\text{in}$, where $n_\text{in}$ is the number of inputs to the neuron, so that signals neither shrink nor grow as they pass through many layers. PyTorch's layers start from a scale of this kind by default.

A.4 · From numbers to classes: softmax and cross-entropy

A phase picker does not predict a number. At every sample of a seismogram it must say which of three classes the sample belongs to: noise, P or S. Session 3, Focal mechanisms, met the two-class version: a polarity is up or down, and the logistic loss scores a predicted probability against the observed sign. Here is the version with any number of classes.

The network outputs one number per class at every sample, called a logit: $s_\text{noise}$, $s_P$, $s_S$. Any real values are allowed. The softmax turns them into probabilities that are positive and add up to one:

$$p_k = \frac{e^{s_k}}{e^{s_\text{noise}} + e^{s_P} + e^{s_S}}.$$

The largest logit gets the largest probability, and adding the same constant to all three logits changes nothing. With two classes the softmax is the logistic function of Session 3.

The loss is the cross-entropy. If the true class at a sample is $k$, the loss there is $-\ln p_k$: zero when the network gives the true class probability one, and growing without bound as that probability goes to zero. It is the negative log-likelihood of the labels, so minimising it is maximum likelihood, as least squares was for Gaussian noise in Session 1.

logits [2.0, 0.5, -1.0] -> probabilities noise 0.79, P 0.18, S 0.04;  cross-entropy if the truth is P: 1.74
logits [0.0, 0.0, 0.0] -> probabilities noise 0.33, P 0.33, S 0.33;  cross-entropy if the truth is P: 1.10
logits [-1.0, 4.0, 0.0] -> probabilities noise 0.01, P 0.98, S 0.02;  cross-entropy if the truth is P: 0.02
ln 3 = 1.10

The middle case is worth remembering. A network that knows nothing gives every class the same probability, one third, and its cross-entropy is $\ln 3 \approx 1.10$. So the loss of a freshly initialised three-class network should be close to 1.10; if it is not, something is wrong with the data or the code before any training has happened. This is the first check in Karpathy's widely used recipe for training networks (Karpathy 2019), and we use it in A.8.

Soft labels

An analyst's pick is not exact: two careful analysts can differ by a few hundredths of a second. PhaseNet therefore does not label a single sample as P. It spreads the label over a short Gaussian centred on the pick: the target probability of P is $\exp[-(i - i_P)^2 / 2\sigma^2]$ at sample $i$, with $\sigma = 10$ samples (0.1 s), and the probability of noise is whatever is left. The target is then a probability distribution $q$ rather than a single class, and the cross-entropy becomes $-\sum_k q_k \ln p_k$. PhaseNet subtracts the part that does not depend on the network, $-\sum_k q_k \ln q_k$, so that a perfect prediction scores exactly zero. The result is the Kullback–Leibler divergence between the target and the prediction, and it is the loss we train with:

perfect logits: 0.0000;  uninformed logits: 0.686

A.5 · The Geysers seismograms, their labels, and an honest split

A picker learns from examples, so before any network we look hard at the examples. They come from the Northern California Earthquake Data Center (NCEDC): for every earthquake in the catalogue, the seismogram recorded at each station, with the P and S arrival times picked for that catalogue. We keep the earthquakes of 2016 inside the producing field and the stations inside the 40 km study circle — the same region as the project notebook — and only velocity sensors:

• BG.DP, the Geysers' own borehole geophones, operated for the field and densest near the steam wells;
• NC.EH, short-period seismometers of the regional network, almost all recording only the vertical component;
• NC.HH and BK.HH, broadband seismometers.

Every example is a window of 30 s, sampled 100 times a second (3,000 samples), on three components: east, north and vertical (E, N, Z). A station that records only Z has zeros on E and N. Amplitudes range over four orders of magnitude between a small distant earthquake and a large near one, so each window is divided by its own standard deviation, the standardisation of A.3.

The file also holds the arrival times — the labels — as sample indices inside each window, and a table that says, for every window, which earthquake and station it is, the distance between them, the magnitude, and which part of the experiment the window belongs to.

2,540 windows of 30 s from 930 earthquakes at 55 stations

                         earthquakes  earthquake windows  noise-only windows
training                         508               1,119                 163
validation                       210                 454                  70
test                             207                 452                  68
five case earthquakes              5                 214                   0

Each row is one part of the data, and A.5's last subsection explains why there are three plus the cases. An earthquake window holds the arrivals of an earthquake at one station. A noise-only window was cut from the same recording before the earthquake arrived, and checked against the catalogue to contain no other earthquake; a picker must learn to stay silent on it. The five case earthquakes are kept whole — every station of each — to look at closely, and are used for nothing else. We begin with the largest of the five, the M 2.9 of 31 March 2016 — the reference event of the project notebook — drawn as a record section: the vertical component of every station, placed at its distance from the epicentre.

50 stations within 40 km; P labelled on 50, S labelled on 18; farthest S label at 13.3 km

The P arrivals (red) line up along a curve that steepens with distance — the moveout of A.1 — and the S arrivals (blue) follow later on a steeper one. Every station has a P label. The S labels stop well inside the circle, although the S wave is plainly visible on the farther traces. The analyst did not pick them: S is harder to pick than P, because it arrives inside the P coda, and an S pick far away adds little to a location that the near stations already constrain.

Look also at the far traces after the P: their oscillations have flat tops of constant height. Those are NC.EH short-period stations, whose signals reach the data centre over analog telemetry with a limited range; when the ground moves more than that range allows, the record clips — it saturates at its maximum and stays there. The P onset survives clipping, but the shape of what follows, including the S, is lost. We count windows whose largest component has a run of flat samples at its maximum:

NC.EH windows that clip, by magnitude: M<0.5 0%, M0.5-1.5 0%, M1.5-2.5 7%, M>=2.5 26%
stations of the reference event that clip: 10, all NC.EH

A quarter of the short-period records of earthquakes above M 2.5 clip, and none below M 1.5; the BG geophones and the broadband stations record the same earthquakes without clipping. We keep clipped windows — a picker in operation meets them — but they are one more reason the S labels thin out where they do. This is the most important fact about these labels, and it holds for the whole year:

share of records with an S label, all reviewed 2016 field events:
   0-2  km: 49.8% of 10,348 records
   2-5  km: 30.9% of 38,023 records
   5-10 km: 11.2% of 43,719 records
  10-15 km:  1.6% of 21,308 records
  15-20 km:  1.0% of  5,585 records
  20-30 km:  0.5% of  6,408 records
  30-40 km:  0.1% of  4,876 records
by sensor: DP 20.9% (vertical only: 0%), EH 0.6% (vertical only: 94%), HH 26.1% (vertical only: 0%)

Half the records within 2 km carry an S label, and almost none beyond 15 km do. On the NC.EH stations, nearly all of which record only the vertical component, almost no S is ever picked: S shakes the ground mostly horizontally. A missing S label does not mean there is no S. If a network is told that everything without a label is noise, it will be taught to call distant S arrivals noise. A.11 measures what that costs and how to avoid it.

Which events are reviewed

The catalogue's picks are of two kinds: manual, reviewed by an analyst, and automatic, made by the network's own picking software and never checked. Which events get reviewed depends on their size:

  M -2.0 to 0.5 : 1,385 of 3,441 events reviewed
  M  0.5 to 1.0 :   818 of 5,735 events reviewed
  M  1.0 to 1.5 : 1,358 of 2,487 events reviewed
  M  1.5 to 2.0 :   606 of   613 events reviewed
  M  2.0 to 6.0 :   207 of   214 events reviewed

Every event above M 1.5 is reviewed; most between M 0.5 and 1.0 are not. We use reviewed events only, so that the labels are an analyst's, but keep reviewed small events deliberately, so the picker sees the whole range of sizes. (Part B looks at an automatic-only event.)

Splitting so that the test is honest

A network must be scored on examples it never trained on: a training set to fit the parameters, a validation set to make decisions during training (when to stop, which settings to keep), and a test set, touched once at the end, to report how good the picker is. Session 1 split data at random. Seismograms cannot be split at random, for three physical reasons.

1. The records of one earthquake are not independent. They share the source, the magnitude and the origin time. So we split whole earthquakes, never single records.
2. Earthquakes repeat. At the Geysers, half of all earthquakes have another earthquake within about 110 m of them (below). Two earthquakes that close give nearly identical seismograms at the same station, because the waves travel the same path. With a random split of earthquakes, the test set would be full of near-copies of training examples, and would measure how well the network remembers paths, not how well it picks. So we cut the field into cells of 1 km by 1 km and assign whole cells: every earthquake in a cell goes to the same set.
3. Performance depends on distance and magnitude. Small and distant earthquakes are harder. A test set drawn in natural proportions would be dominated by small, near records and say little about the rest. So the validation and test sets are balanced: the same number of records in every combination of four magnitude ranges and three distance ranges, with the sensor types mixed in each. Every score can then be reported per combination, and an average over combinations weights each equally.

distance from each 2016 earthquake to its nearest neighbour: median 110 m, 10 % of earthquakes within 44 m

val: records per magnitude and distance
dbin      0-5 km  5-15 km  15-40 km
mbin                               
M<0.5         40       40        17
M0.5-1.5      40       40        40
M1.5-2.5      40       40        40
M>=2.5        39       40        38

test: records per magnitude and distance
dbin      0-5 km  5-15 km  15-40 km
mbin                               
M<0.5         40       40        12
M0.5-1.5      40       40        40
M1.5-2.5      40       40        40
M>=2.5        40       40        40

The three sets interleave across the field, cell by cell, so each covers the whole field while no test earthquake has a near-copy in training. The validation and test tables hold forty records in every cell but one: earthquakes below M 0.5 are seldom recorded beyond 15 km, so that combination has fewer — which is itself a fact about detection, not a flaw in the selection.

Finally we turn the labels into training targets, exactly as PhaseNet does (A.4): a Gaussian of width $\sigma = 10$ samples centred on each P and each S, and noise as whatever probability is left.

A.6 · Why convolution

The MLP of A.1 reads its whole input at once, with one weight for every input value. For a window of 3 × 3,000 samples that is 9,000 weights per neuron, and — the real problem — every weight is tied to one position in the window. A neuron that learns what a P arrival looks like at sample 1,000 knows nothing about a P arrival at sample 1,500. An arrival can be anywhere in a window, so the network would have to learn the same pattern again at every position.

We can see this happen. The task: decide whether a 4 s piece of the vertical component (400 samples) contains a P arrival. We train on pieces where the P sits in the first half and test on pieces where it sits in the second half, with equal numbers of pure-noise pieces in both.

training pieces 4175, of which 2237 contain a P; test pieces 827

A convolution reads the input differently. It slides one short kernel — a handful of weights, here 7 — along the trace and at every position takes the weighted sum of the samples under it:

$$y_i = \sum_{j=0}^{K-1} w_j\, x_{i+j} + b.$$

The same $K$ weights are used at every position: weight sharing. Two things follow. The layer has $K + 1$ parameters however long the trace is, instead of one per sample. And if the input shifts, the output shifts by the same amount and is otherwise unchanged: the layer is translation-equivariant. A kernel that responds to the onset of a P wave responds to it wherever it occurs. A convolutional layer applies several kernels at once, each producing one output trace called a channel, and reads several input channels (the three components) by giving each kernel one set of weights per input channel. Layers stack like the MLP's: convolution, then ReLU, then the next convolution.

To turn a trace of feature channels into one yes/no answer, the convolutional network takes the maximum of each channel over the whole piece — "did this feature occur anywhere?" — and feeds those maxima to a final linear layer. We train both networks the same way: the MLP, and a small convolutional network (CNN).

MLP: 25,729 parameters; accuracy on training pieces 98%, new pieces with P in the first half 81%, P in the second half 65%
CNN:  3,761 parameters; accuracy on training pieces 95%, new pieces with P in the first half 93%, P in the second half 93%

Both networks learn the training pieces, and both do well on new pieces with the P where they saw it in training. When the P moves to the second half, the MLP drops most of the way to guessing (50 %): its weights for the second half never saw an arrival. The CNN, with several times fewer parameters, does as well on the shifted pieces as on the others. That is the case for convolution in one figure: the right structure — here, that an arrival looks the same wherever it occurs — lets a network learn from far fewer examples, and generalise to cases it never saw. Building such an assumption into a network's design is called an inductive bias.

STA/LTA is a convolution

Seismologists used convolutions to detect arrivals long before neural networks. The classic trigger, STA/LTA (Allen 1978), divides the average energy in a short window (the short-term average, STA, here 0.5 s) by the average in a long window before it (the long-term average, LTA, here 5 s), and declares an arrival where the ratio jumps. Each average is a convolution of the squared trace with a box-shaped kernel of fixed weights. A CNN does the same kind of computation, but learns its kernels from the labels instead of having them chosen by hand — and stacks many of them.

STA/LTA peaks at 23.39 s; P label at 23.11 s

The ratio jumps at the P arrival, and again, less clearly, at the S. STA/LTA detects that something arrived; it does not say whether it is P or S, and its peak is late by a fraction of the short window. Telling the phases apart needs context — what came before — which is the next question.

A.7 · How far the network sees: the receptive field

An output sample of a convolutional network depends only on a limited stretch of the input around it: its receptive field. One kernel of 7 samples sees 7 samples. A second layer of kernel 7 on top sees 13: each of its 7 inputs saw 7. In general $n$ stacked layers of kernel $K$ see $n(K-1) + 1$ samples — the field grows by only $K - 1$ per layer.

That is too slow. To call an arrival S rather than P, the network must see the P that came before it, so its receptive field must span the S–P time. At the Geysers that is mostly under 1.5 s (150 samples), but it reaches several seconds at 40 km, and to call a window noise-free the network should see much more. The fix is downsampling: a convolution with stride $s$ moves its kernel $s$ samples at a time, so its output is $s$ times shorter, and every layer after it sees $s$ times further in input samples. PhaseNet downsamples by 4 at each of four stages.

The receptive field can be measured rather than derived. Ask how one output sample depends on every input sample: the gradient of that output with respect to the input, by backpropagation (A.2). Wherever the gradient is non-zero, the input can affect the output.

four layers of kernel 7: receptive field 25 samples (formula 4 x 6 + 1 = 25)
the same four layers, the last three with stride 4: 133 samples

A.8 · A probability at every sample: the U-Net

A picker must answer at every sample: noise, P or S. That is segmentation — labelling each element of the input — and the network built for it is the U-Net (Ronneberger et al. 2015), designed for images and adapted to seismograms by PhaseNet. It has two halves.

• The encoder (the left arm of the U) is a stack of convolutions with ReLU, with a strided convolution after each stage that shortens the trace four times. Each stage doubles the number of channels: as the trace gets shorter, each sample describes a longer stretch of the input in more ways. After four stages a trace of 3,000 samples is 12 samples long and every one of them sees a large part of the window — the context needed to tell P from S.
• The decoder (the right arm) undoes the shortening with transposed convolutions, which spread each sample back over four, until the trace is 3,000 samples long again.

Downsampling throws away timing: after four stages one sample stands for 256 input samples, 2.56 s, far too coarse to place a pick to a hundredth of a second. The U-Net's answer is the skip connection: at each stage, the decoder receives the encoder's output from the same stage, before it was shortened, and stacks it with its own channels. The coarse path says what is there; the skip connection says exactly where. A last convolution of kernel 1 turns the channels into the three logits of A.4, at every sample.

Two further ingredients sit after every convolution. ReLU is A.1's activation. Batch normalisation standardises each channel over the mini-batch — subtract its mean, divide by its standard deviation, then rescale with two learned numbers — which, like the standardisation of A.3, keeps the numbers inside a deep network of order one and makes training far less sensitive to the learning rate.

The class below is PhaseNet's UNet2018 written for one-dimensional traces, with its names, its kernel of 7 samples, its stride of 4, its 5 stages and its 8 starting channels. The published code pads its strided and transposed convolutions exactly as the original TensorFlow model did, which is needed to reuse the published weights; we train our own weights, so we use plain symmetric padding.

UNet2018: 268,443 parameters; output shape for one window (1, 3, 3000)
receptive field: 3,331 samples = 33.3 s

One output sample can be affected by about half a minute of input, which spans any S–P time in the field many times over. Every 3,000-sample window is therefore seen almost whole by every output sample.

Training

Training is A.3's loop with three additions. The data come in mini-batches of 32 windows. After every epoch the network is scored on the validation set, which it never trains on. And we keep the parameters from the epoch with the lowest validation loss — early stopping: if the network begins to memorise its training windows, the validation loss stops falling and starts to rise, and the weights we keep are those from before that happened.

First, Karpathy's check from A.4: before any training, the loss should be what a network that knows nothing would score.

validation loss of the untrained network 0.918; of a network that gives every class 1/3: 1.087

Both are close to $\ln 3 \approx 1.10$, as A.4 predicted for a network that knows nothing, so the data, the labels and the loss are wired up correctly. (The uniform network scores slightly less than $\ln 3$ because near each pick the target itself is uncertain, and the KL divergence does not charge for that; the untrained network differs a little more because its random weights are not exactly uninformed.)

Now the training run: 30 epochs, a few minutes on a GPU and longer on a processor alone (see the setup cell).

epoch  5: training loss 0.2158, validation loss 0.1929  (4 s)

epoch 10: training loss 0.0664, validation loss 0.0640  (7 s)

epoch 15: training loss 0.0316, validation loss 0.0337  (11 s)

epoch 20: training loss 0.0194, validation loss 0.0235  (15 s)

epoch 25: training loss 0.0137, validation loss 0.0197  (19 s)

epoch 30: training loss 0.0101, validation loss 0.0177  (22 s)
kept the weights of epoch 30

The loss curves and one validation window, before we score anything:

Training and validation losses fall together, and the predicted probabilities (thin lines) sit on the labels (wide pale lines). A.10 turns these probabilities into picks and scores them.

A.9 · Learned or memorised?

The network has about 270,000 parameters and only about 1,300 training windows. With so many parameters it could memorise every training window — learn the noise of each one by heart — and still pick nothing new. Why did it not? To see what memorising looks like, we train the same network on just 40 windows, in mini-batches of 8, for 300 epochs.

40 training windows: final training loss 0.0032, validation loss 0.0404 (lowest 0.0391)
all 1282 training windows: training loss 0.0101, validation loss 0.0177

With 40 windows the training loss keeps falling, far below anything the full training run reached, while the validation loss stops improving and stays more than twice as high as with all the windows (dashed). The gap between the two curves is the network fitting the particular noise of its 40 windows, which tells it nothing about new ones. That is overfitting, the bias–variance trade-off of Session 1 in a network: the training loss measures memory, only the validation loss measures learning.

The best remedy for overfitting is more data, and the second best is data that are varied. A network trained on a few thousand windows, from earthquakes all over the field, of every size and at every distance, cannot memorise them all, so it has to learn what is common to them: what an arrival looks like. The common remedies, besides early stopping, are:

• data augmentation: make new training examples by changes that do not change the answer. Our windows place the P at a random position, so the network cannot learn a position, and PhaseNet also adds noise and stacks two earthquakes in one window;
• weight decay: add $\lambda \sum \theta^2$ to the loss, the ridge penalty of Session 1, which keeps the weights small (Adam's variant is called AdamW);
• dropout: during training, set a random fraction of the neurons' outputs to zero at every step, so no neuron can rely on any single other one. UNet2018 can apply it after every layer.

A.10 · From probabilities to picks, and how good they are

A pick is a time, not a curve. PhaseNet turns each probability trace into picks with one rule: a pick is a local maximum of the probability that is higher than a threshold, and is the highest value within half a second either side. Its time is the sample of the maximum, and its probability is the pick's score. PhaseNet's default threshold is 0.3.

A pick is correct if it lies within 0.1 s (10 samples) of a label, the tolerance used to score PhaseNet. Counting over windows:

• a true positive is a label with a correct pick;
• a false negative is a label with no pick near it — a missed arrival;
• a false positive is a pick near no label — an arrival that is not there.

Two numbers summarise them. Precision is the fraction of picks that are correct, TP / (TP + FP); recall is the fraction of labelled arrivals that were found, TP / (TP + FN). Raising the threshold makes fewer, surer picks: precision rises and recall falls. The F1 score, $2PR/(P+R)$, is high only when both are.

One subtlety comes from A.5. An S pick on a window with no S label cannot be scored: the S may be there, unlabelled. So S is scored only on windows that carry an S label, and on noise windows, where any pick is false.

threshold 0.1: P precision 0.75 recall 0.81 F1 0.78 | S precision 0.79 recall 0.82 F1 0.81

threshold 0.3: P precision 0.84 recall 0.80 F1 0.82 | S precision 0.82 recall 0.80 F1 0.81

threshold 0.5: P precision 0.87 recall 0.79 F1 0.83 | S precision 0.84 recall 0.79 F1 0.81

threshold 0.7: P precision 0.92 recall 0.76 F1 0.83 | S precision 0.86 recall 0.72 F1 0.78

The F1 score changes little between 0.3 and 0.5, so we keep PhaseNet's 0.3. Now, once, the test set: earthquakes from cells the network never saw, balanced over magnitude and distance.

P, all test windows: precision 0.87, recall 0.82, F1 0.85 over 482 labelled arrivals; false picks on the 68 noise windows: 14
S, all test windows: precision 0.85, recall 0.84, F1 0.84 over 268 labelled arrivals; false picks on the 68 noise windows: 8

The picker is best for larger earthquakes recorded nearby, and degrades towards small earthquakes and large distances, where the arrival is weak against the noise. For S, the far column is empty: almost no S is labelled beyond 15 km, so there is nothing to score — not evidence that the picker fails there. Because the test set is balanced, each cell is measured with comparable precision, and the overall numbers weight all cells similarly. A test set drawn in natural proportions, mostly small and near, would report a single number that hides the far column.

How precise are the correct picks?

The analyst's pick is not the truth either. Each P label carries the analyst's quality remark: impulsive (i), a sharp onset, or emergent (e), a gradual one whose start is uncertain.

P residuals (correct picks): median 0 ms, standard deviation 38 ms

Picks agree with impulsive labels more closely than with emergent ones. Part of that spread is the label's own uncertainty: when the analyst could not see the onset clearly, no picker can be scored better than the analyst's guess. A network cannot be more precise than its labels.

Against STA/LTA

The baseline every seismic network ran for decades is the STA/LTA trigger of A.6. As a P picker we take the first time the ratio exceeds a threshold, chosen on the validation set for the best F1, and score it exactly like the network.

STA/LTA   P F1 by distance: 0-5 km 0.69, 5-15 km 0.67, 15-40 km 0.28;  overall 0.56
UNet2018  P F1 by distance: 0-5 km 0.91, 5-15 km 0.92, 15-40 km 0.71;  overall 0.85
(STA/LTA trigger level chosen on validation: 4)

STA/LTA finds the strong, near arrivals, but it times them by when the energy ratio crosses a level, which is after the onset, and it fires on noise bursts. The network, trained on analysts' picks, times arrivals the way the analysts did and keeps its precision at distance.

A.11 · The labels decide what is learnable

A.5 showed that most S arrivals are not labelled, and that the gaps are not random: they grow with distance and with vertical-only sensors. In A.8 we trained with the simplest rule — every sample without a label is noise — so on a window with no S label, the network was taught that the S arrival in it is noise. What did that teach?

The fix is to teach only what we know. On a window without an S label, over the stretch where the S could arrive — from just after the P until $0.15\,\mathrm{s/km} \times$ distance $+ 1$ s later, which covers the S–P time at any distance in the field — we do not know whether a sample is noise or S. We do know it is not P: the P was picked, and it came earlier. So there we replace the loss by the likelihood of what we know, a partial label:

$$L_i = -\ln\left(p_\text{noise} + p_S\right),$$

which is zero whenever the network puts all its probability on noise or S, and does not care how it divides it between them. Everywhere else, including the coda after that stretch, the full label stands. (Simply leaving the stretch out of the loss — masking it, as PhaseNet-Plus's phase_mask can — also stops the network learning that S is noise, but then nothing tells it that an arrival in that stretch is not a P, and it starts calling distant S arrivals P.) We retrain with the partial label, everything else the same.

every gap is noise : P F1 0.85, S F1 on labelled S 0.84

partial label      : P F1 0.85, S F1 on labelled S 0.88

On the arrivals we can score, the two networks are close for P, and the partial label scores somewhat better on S. The real difference is where we cannot score: the S arrivals nobody labelled. There is still a physical test. For a wave path of length $R$, the S–P time is $R\,(1/v_S - 1/v_P)$: it grows with distance, and it can be no shorter than for the fastest rock the waves could have crossed. In crust no faster than $v_P = 6$ km/s, with the usual ratio $v_P/v_S = 1.73$, that floor is about 0.12 s per kilometre of epicentral distance (the path is never shorter than the epicentral distance). A pick far below the floor cannot be an S. We take every test window without an S label and plot each network's S picks against the floor, beside the analysts' S picks on labelled windows.

S-P floor: 0.122 s/km, so at least 3.6 s at 30 km; analysts' S picks below 0.7 x floor: 0%
every gap is noise: S picks on 54 of 191 unlabelled test windows; within 15 km 41 (0% too early to be S), beyond 15 km 13 (31% too early)
partial label: S picks on 103 of 191 unlabelled test windows; within 15 km 48 (0% too early to be S), beyond 15 km 55 (64% too early)

The dashed line is the floor and the shaded wedge lies well below it; the analysts' S picks sit above it, as they must. Within 15 km both networks pick S on unlabelled windows, the partial-label network more of them, and almost none fall in the wedge: they behave like real S arrivals the analysts did not pick. The network trained with every gap as noise picks far fewer: it learned that much of what follows a P is noise — exactly what its labels said. A network learns its labels, including their omissions.

Beyond 15 km the physics gives a sharper verdict. The partial-label network now makes many "S" picks there, but many sit about a second after the P, inside the wedge: at 30 km a real S comes at least several seconds after the P. These picks are not S. They are most likely a later P arrival — at these distances the first arrival is a head wave, and the direct P that A.1 described follows it. The partial label stopped the network learning that distant S is noise, but it could not teach it what a distant S looks like, because no training window had a labelled S beyond about 15 km. A loss can stop a network learning the wrong thing from missing labels; only labels can teach it the right thing. No score on labelled arrivals could have shown this; the S–P line did.

A.12 · Does a probability of 0.5 mean a 50 % chance?

Each pick carries a probability. Is it a probability in the everyday sense — of all the picks scored 0.5, are half correct? A predictor whose probabilities mean this is calibrated. Session 2 asked the same of an error ellipse: does the 95 % ellipse contain the truth 95 % of the time?

The check is a reliability diagram. Take every P pick on the test set, at a low threshold of 0.05 so that weak picks are included, group the picks by their probability, and plot the fraction correct in each group against the group's mean probability. A calibrated picker lies on the diagonal. Each group is a sample of picks, so its fraction correct has an uncertainty: we draw the binomial standard error, $\sqrt{f(1-f)/n}$ for a fraction $f$ of $n$ picks. The expected calibration error (ECE) is the average distance from the diagonal, each group weighted by its number of picks.

674 P picks; 46% have probability above 0.9, and 98% of those are correct; picks between 0.3 and 0.7: 45, 33% correct

The grey numbers are the picks in each group. Almost half the picks are confident, above 0.9, and nearly all of those are correct: at the top the curve meets the diagonal. Below about 0.6 it lies under the diagonal, beyond its error bars: a pick scored 0.5 is right only about a quarter of the time. The weak picks are overconfident. The ECE is nevertheless small, because it weights each group by its number of picks and few picks are weak — a single number that hides where the calibration fails. So the answer to the question is no: here a 0.5 means roughly a one-in-four chance, and only the strong picks mean what they say.

Even where it holds, calibration is a measurement, not a guarantee. A network's softmax outputs are not probabilities by construction, because training rewards being right, not being honest about doubt, and large modern networks are often overconfident (Guo et al. 2017). Calibration also changes when the data change: a picker trained on one region and run on another may keep its ranking of picks while its probabilities drift. So a threshold such as 0.5 is a number on a scale the network invented, and before a probability is used as one — to weight a pick in a location, to decide what an analyst should review — its reliability diagram must be measured on data the network did not train on. Part B asks this of the published PhaseNet on Geysers data.

Optional: recalibrating

When a reliability curve does lie off the diagonal, it can be corrected after training without changing which picks are found. Temperature scaling divides every logit by one number $T$ before the softmax, fitted on the validation set: $T > 1$ softens the probabilities and $T < 1$ sharpens them (Guo et al. 2017). Fitting a monotone function of the pick score to whether picks were correct on the validation set — a logistic curve (Platt scaling) or a step function (isotonic regression) — does the same at the level of picks. Because each is monotone, a stronger pick stays stronger: the ranking, and so the precision–recall trade-off, is unchanged, and only what a number means moves.

A.13 · The programs the field runs

• STA/LTA (Allen 1978) and its descendants, the trigger in almost every network's real-time system for four decades. Hand-designed convolutions (A.6).
• PhaseNet (Zhu & Beroza 2018), the U-Net of A.8, trained on northern California picks. Our UNet2018 has its layers, but not the original's padding; the published weights need that padding to reproduce PhaseNet's picks, so Part B adds it before running them.
• GPD (Ross et al. 2018), a convolutional network that classifies short windows as P, S or noise — the detector of A.6 grown up — slid along the trace.
• EQTransformer (Mousavi et al. 2020), which adds attention layers and detects the earthquake and picks its phases together.
• SeisBench (Woollam et al. 2022), the common library that packages these networks, their published weights and benchmark datasets in one interface, so that pickers can be compared on the same data. It is the obvious tool for a student who wants to use a picker today.
• PhaseNet-Plus, our group's current code, which adds first-motion polarity and event detection heads and the label mask of A.11; the functions in this notebook are distilled from it.

Takeaways

Machine learning

• A neural network is a flexible function of many parameters — with ReLU activations, a piecewise-linear one — fitted by minimising a loss, as least squares is.
• Backpropagation is the chain rule applied from the loss backwards, and it gives the gradient with respect to every parameter for about the cost of one forward pass.
• Gradient descent needs a learning rate found by watching the loss; mini-batches, momentum and Adam make it fast, and more neurons than strictly needed make it reliable.
• Softmax and cross-entropy turn a network's outputs into class probabilities and score them by likelihood; soft Gaussian labels make the loss PhaseNet's KL divergence.
• Convolution shares one small kernel across all positions, so a convolutional network recognises an arrival wherever it occurs and generalises where a fully connected network does not.
• A U-Net shortens the trace to see context and restores timing through skip connections, so it can give a probability at every sample.
• Overfitting shows as a validation loss that rises while the training loss falls; more varied data, augmentation and early stopping prevent it.
• A network learns its labels including their omissions: a loss that scores only what is known (here, "not P") stops it learning that unlabelled arrivals are noise, but only labels can teach it what those arrivals look like.
• A softmax output is not a calibrated probability by construction; whether it is must be measured with a reliability diagram on held-out data.

Seismology

• A phase picker turns a three-component seismogram into P and S arrival times, the input to association, location and everything after them in this course.
• Records of one earthquake share its source, and nearby earthquakes share paths, so a picker is tested on earthquakes from places it never trained on, split by space, not at random.
• Picking skill depends on magnitude and distance, so a picker is scored on a test set balanced over both, and reported per combination.
• Catalogue S picks are missing mostly at distance and on vertical-only sensors, so an unlabelled S is not an absent S.
• A picker's timing can be no better than its labels: picks agree more closely with impulsive than with emergent analyst picks.
• S–P time grows in proportion to distance, which gives a physical check on picks that no label can score, and it shows that distant "S" picks learned without distant S labels are not S.

References

• Allen, R. V. (1978). Automatic earthquake recognition and timing from single traces. Bull. Seismol. Soc. Am., 68, 1521–1532. https://doi.org/10.1785/BSSA0680051521
• Guo, C., Pleiss, G., Sun, Y., & Weinberger, K. Q. (2017). On calibration of modern neural networks. ICML. https://arxiv.org/abs/1706.04599
• He, K., Zhang, X., Ren, S., & Sun, J. (2015). Delving deep into rectifiers. ICCV. https://arxiv.org/abs/1502.01852
• Ioffe, S., & Szegedy, C. (2015). Batch normalization. ICML. https://arxiv.org/abs/1502.03167
• Karpathy, A. (2019). A recipe for training neural networks. http://karpathy.github.io/2019/04/25/recipe/
• Kingma, D. P., & Ba, J. (2015). Adam: A method for stochastic optimization. ICLR. https://arxiv.org/abs/1412.6980
• Mousavi, S. M., Ellsworth, W. L., Zhu, W., Chuang, L. Y., & Beroza, G. C. (2020). Earthquake transformer. Nat. Commun., 11. https://doi.org/10.1038/s41467-020-17591-w
• Prince, S. J. D. (2023). Understanding Deep Learning. MIT Press. https://udlbook.github.io/udlbook/
• Ronneberger, O., Fischer, P., & Brox, T. (2015). U-Net: Convolutional networks for biomedical image segmentation. MICCAI, 234–241. https://doi.org/10.1007/978-3-319-24574-4_28
• Ross, Z. E., Meier, M.-A., Hauksson, E., & Heaton, T. H. (2018). Generalized seismic phase detection with deep learning. Bull. Seismol. Soc. Am., 108, 2894–2901. https://doi.org/10.1785/0120180080
• Woollam, J., Münchmeyer, J., Tilmann, F., Rietbrock, A., et al. (2022). SeisBench — a toolbox for machine learning in seismology. Seismol. Res. Lett., 93, 1695–1709. https://doi.org/10.1785/0220210324
• Zhu, W., & Beroza, G. C. (2018). PhaseNet: a deep-neural-network-based seismic arrival-time picking method. Geophys. J. Int. https://doi.org/10.1093/gji/ggy423