Template matching · Part A: a template and a neural network, side by side¶
Many earthquakes at the Geysers are too small to stand out of the noise, so they never reach the catalogue. Yet many of them happen almost exactly where an earlier, larger earthquake happened, and their seismograms then look almost the same. Can we find a small earthquake by searching the continuous record for a copy of a known one — and how does that compare with the neural-network picker of Session 5?
Session 5, Deep learning, trained a network to recognise P and S arrivals from thousands of labelled examples. Template matching takes the opposite route. It uses no training at all: the seismogram of one catalogued earthquake, the template, is slid along the continuous record, and wherever the record looks enough like the template we declare a new earthquake. The idea goes back to the radar engineers' matched filter, and in seismology it has found ten times more earthquakes than the routine catalogue in Southern California (Ross et al. 2019).
This notebook builds template matching from one idea — similarity as a cosine — and then sets it against the network on the same data, under conditions we control, so that we can say exactly what each method can and cannot find.
The questions this notebook answers¶
- How do we measure whether two seismograms are alike? (A.1)
- How is that measure computed along a long record, and why is it a convolution? (A.2)
- Which earthquakes look alike, and in which frequency band? (A.3)
- How does a correlation become a detection, and how many detections would noise alone give? (A.4)
- Why use many stations at once? (A.5)
- How precisely can correlation time an arrival? (A.6)
- Template matching against a neural network: which finds more, and which finds what? (A.7)
- Is there a method between the two? (A.8)
Part B applies both methods to one hour of continuous Geysers data from every station of the array.
import importlib.util
import sys
import time
import urllib.request
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.ndimage
import scipy.signal
import scipy.stats
import torch
import torch.nn as nn
import torch.nn.functional as F
BASE = "https://raw.githubusercontent.com/AI4EPS/EPS207_Observational_Seismology/main/docs/notebooks"
for module in ("geysers_data.py",):
if "google.colab" in sys.modules or importlib.util.find_spec(module[:-3]) is None:
for attempt in range(4):
try:
urllib.request.urlretrieve(f"{BASE}/{module}", module)
break
except Exception as exc:
if attempt == 3:
raise RuntimeError(f"could not fetch {module}: {exc}; re-run this cell") from None
time.sleep(2 * (attempt + 1))
import geysers_data as gd
rng = np.random.default_rng(207)
torch.manual_seed(207)
plt.rcParams.update({"figure.dpi": 120, "font.size": 9, "axes.grid": True,
"grid.alpha": 0.3, "axes.axisbelow": True})
DEVICE = ("cuda" if torch.cuda.is_available() else
"mps" if torch.backends.mps.is_available() else "cpu")
print(f"torch {torch.__version__}; PhaseNet (A.7) will run on: {DEVICE}")
torch 2.14.1; PhaseNet (A.7) will run on: mps
The data: 265 earthquakes from one cubic kilometre¶
We need many earthquakes that occurred close together, recorded on the same stations. The densest part of the 2016 Geysers catalogue is a ball of radius 1 km near the centre of the field: it holds 265 earthquakes of magnitude 1.0 and above whose picks an analyst reviewed. We have their seismograms on 10 stations between about 1 and 12 km away, on three components (east, north, vertical: E, N, Z), sampled 100 times a second. Each seismogram starts 1 s before the origin time and lasts 10 s, and the analyst's P and S picks (arrival times) come with it.
For one station, SQK, about 1.6 km from the centre of the ball, we also have the first minute of
every recording: 40 s of background noise before the earthquake, then the earthquake. Strung
together, those 265 minutes stand in for 4.4 hours of continuous data, in which we know where
265 catalogued earthquakes are. Finally, the file holds the SQK seismograms of 30 earthquakes from elsewhere
in the field — earthquakes unlike those of the ball, which we will need in A.7.
The amplitudes are ground velocity in micrometres per second, as the archive gives them. The table below also lists each station's median P amplitude and noise level: two stations, HER and RTB, are far weaker than stations at similar distances, which suggests a gain error in the archive rather than in the ground. Correlation ignores amplitude, so it does not affect what follows, but it is a reminder to look at amplitudes before trusting them.
pack = np.load(gd.asset("geysers_patch_2016.npz"))
X = pack["window"].astype(np.float64) * pack["window_scale"][..., None] # (event, station, E N Z, 1000)
REC = pack["record"].astype(np.float64) * pack["record_scale"][..., None] # (event, E N Z, 6000) at SQK
FIELD = pack["field"].astype(np.float64) * pack["field_scale"][..., None] # (30, E N Z, 900), from 1 s before P
P, S = pack["p_index"].astype(int), pack["s_index"].astype(int) # sample in the window; -1 = no pick
events = pd.read_csv(gd.asset("geysers_patch_2016_events.csv"))
stations = pd.read_csv(gd.asset("geysers_patch_2016_stations.csv"))
patch, field = events[events.group == "patch"].reset_index(drop=True), events[events.group == "field"]
OFFSET = pack["record_offset"].astype(int) # where each 10 s SQK window starts in REC
FS, PRE = 100, 1.0 # Hz; the windows start 1 s before the origin
NE, NS = X.shape[:2]
SQK = int(np.where(stations.station == "SQK")[0][0])
# event positions in km about the centre of the ball, for distances
xyz = np.c_[(patch.event_longitude - patch.event_longitude.mean()) * gd.KM_LON,
(patch.event_latitude - patch.event_latitude.mean()) * gd.KM_LAT, patch.event_depth_km]
sep = np.linalg.norm(xyz[:, None] - xyz[None], axis=-1) # inter-event distance (km)
stations["distance_km"] = np.hypot((stations.station_longitude - patch.event_longitude.mean()) * gd.KM_LON,
(stations.station_latitude - patch.event_latitude.mean()) * gd.KM_LAT)
print(f"{NE} earthquakes, M {patch.event_magnitude.min():.1f}-{patch.event_magnitude.max():.1f} "
f"(median {patch.event_magnitude.median():.1f}), depth {patch.event_depth_km.min():.1f}-"
f"{patch.event_depth_km.max():.1f} km, {pd.to_datetime(patch.event_time).min():%d %b} to "
f"{pd.to_datetime(patch.event_time).max():%d %b %Y}")
print(f"median distance to the nearest other earthquake: {np.median(np.sort(sep, 1)[:, 1]) * 1000:.0f} m; "
f"pairs closer than 100 m: {(sep[np.triu_indices(NE, 1)] < 0.1).sum()}")
Xz = scipy.signal.sosfiltfilt(scipy.signal.butter(4, [5, 40], "bandpass", fs=FS, output="sos"), X[:, :, 2], axis=-1)
stations["P peak (um/s)"] = [np.median([np.abs(Xz[i, s, P[i, s]:P[i, s] + 100]).max() for i in np.where(P[:, s] >= 10)[0]])
for s in range(NS)]
stations["noise (um/s)"] = [np.median([Xz[i, s, :50].std() for i in range(NE) if Xz[i, s, :50].any()]) for s in range(NS)]
print(stations[["network", "station", "instrument", "distance_km", "P peak (um/s)", "noise (um/s)"]].round(3).to_string(index=False))
print(f"P picked on {(P >= 0).mean():.0%} and S on {(S >= 0).mean():.0%} of the {NE * NS:,} seismograms")
265 earthquakes, M 1.0-3.2 (median 1.4), depth 1.6-3.3 km, 02 Jan to 26 Dec 2016
median distance to the nearest other earthquake: 115 m; pairs closer than 100 m: 71
network station instrument distance_km P peak (um/s) noise (um/s)
BG SB4 DP 0.940 129.270 0.364
BG SQK DP 1.573 123.994 0.055
NC GDXB HH 2.050 105.058 0.085
BG DRK DP 2.852 56.991 0.043
BG AL3 DP 3.792 36.652 0.013
BG HVC DP 5.080 18.754 0.013
BG MCL DP 4.880 17.287 0.013
BG FNF DP 6.402 15.356 0.015
BG HER DP 9.061 0.187 0.005
BG RTB DP 11.882 0.697 0.002
P picked on 94% and S on 35% of the 2,650 seismograms
A pick missing on a seismogram means only that the analyst did not mark that phase there; S is marked much less often than P, as Session 5 found. Where we need a pick that the analyst did not make, we will say how we placed it.
One more preparation runs through the whole notebook. Earthquakes of magnitude 1 radiate most of
their energy above a few hertz, while ocean waves and wind put energy below 1 Hz, so a seismogram
is first band-pass filtered: we keep 5 to 40 Hz, with a zero-phase Butterworth filter of
order 4, which does not shift arrivals in time (sosfiltfilt runs the filter forwards and then
backwards). Part B uses the same band.
def bandpass(x, lo=5.0, hi=40.0):
sos = scipy.signal.butter(4, [lo, hi], "bandpass", fs=FS, output="sos")
return scipy.signal.sosfiltfilt(sos, x, axis=-1)
Xf = bandpass(X)
t_win = np.arange(X.shape[-1]) / FS - PRE # time after the origin (s)
A.1 · Similarity is a cosine¶
Look first at three earthquakes of the ball on the vertical component at SQK, each aligned on its P pick: a pair less than 100 m apart in the catalogue, and a third earthquake about 0.8 km from the first.
def piece(i, s=SQK, comp=2, phase=P, length=100):
"""`length` samples of a filtered seismogram, from 0.1 s before the pick."""
a = phase[i, s] - 10
return Xf[i, s, comp, a:a + length]
def cosine(x, y):
x, y = x - x.mean(), y - y.mean()
return x @ y / (np.linalg.norm(x) * np.linalg.norm(y))
ok = np.where(P[:, SQK] >= 10)[0]
close = [(cosine(piece(i), piece(j)), i, j) for i in ok for j in ok if i < j and sep[i, j] < 0.1]
_, i0, near = max(close) # the most alike of the close pairs
cand = ok[(np.abs(sep[i0, ok] - 0.8) < 0.1)]
far_ = cand[np.argsort([cosine(piece(i0), piece(j)) for j in cand])[len(cand) // 2]] # a typical one
fig, ax = plt.subplots(figsize=(8, 3))
for k, (j, c) in enumerate(((i0, "k"), (near, "tab:red"), (far_, "tab:blue"))):
z = Xf[j, SQK, 2]
lab = "the first" if j == i0 else f"{sep[i0, j] * 1000:.0f} m from the first"
sh = (P[j, SQK] - P[i0, SQK]) / FS # align on the P pick
ax.plot(t_win - sh, z / np.abs(z).max() - 2.2 * k, c=c, lw=0.7, label=f"M {patch.event_magnitude[j]:.1f}, {lab}")
ax.set(xlim=(0.5, 4), xlabel="time after the first earthquake's origin (s), aligned on P", yticks=[],
title="Three earthquakes of the ball at SQK, vertical, 5-40 Hz")
ax.legend(loc="upper right", fontsize=7)
plt.show()
The close pair looks alike wiggle for wiggle; the third does not, although it is less than a kilometre away. We want a number for "look alike".
Write a stretch of seismogram of \(n\) samples as a vector \(\mathbf{x} = (x_1, \dots, x_n)\), and the same stretch of another seismogram as \(\mathbf{y}\). Two vectors point the same way when the angle \(\theta\) between them is small. The cosine of that angle is the inner product of the two divided by their lengths:
It is 1 when \(\mathbf{y}\) is \(\mathbf{x}\) multiplied by any positive number, \(-1\) when the number is negative, and near 0 when the two are unrelated. Dividing by the lengths is what makes it ignore size: a magnitude 1 and a magnitude 2 at the same place have the same shape but ground motions about ten times different, and the cosine still says they are alike.
If we first subtract from each vector its own mean, the cosine becomes the correlation coefficient of statistics,
where \(\bar x\) and \(\sigma_x\) are the mean and standard deviation of \(\mathbf{x}\). This is the
normalised cross-correlation, the measure of similarity everyone in the field uses. We check
both claims — it equals np.corrcoef, and it ignores amplitude — on the three earthquakes, using
1 s of the vertical component starting 0.1 s before the P pick (the functions piece and
cosine are defined in the cell above).
x = piece(i0)
for j in (near, far_):
y = piece(j)
print(f"{sep[i0, j] * 1000:4.0f} m apart: cosine {cosine(x, y):.2f}, np.corrcoef {np.corrcoef(x, y)[0, 1]:.2f}, "
f"cosine with y x 1000 {cosine(x, 1000 * y):.2f}")
47 m apart: cosine 0.74, np.corrcoef 0.74, cosine with y x 1000 0.74 876 m apart: cosine 0.05, np.corrcoef 0.05, cosine with y x 1000 0.05
A.2 · Sliding the template is a convolution¶
To search a continuous record for copies of a template, we compute \(\rho\) between the template and every stretch of the record of the same length: first the stretch starting at sample 1, then at sample 2, and so on. That gives a correlation trace \(c(\tau)\), one value per start sample \(\tau\). Write the template as \(\mathbf{t}\), of \(n\) samples, and the record as \(d\). Then
Two things in this formula deserve a closer look.
The numerator is a convolution. Session 5's A.6 wrote a convolutional layer as \(\sum_k w_k\, d_{\tau+k}\): a short kernel \(w\) slid along the input. The numerator above is exactly that, with the kernel fixed to the demeaned template. (The mean of the record stretch can be left out of the numerator, because \(\sum_k (t_k - \bar t) = 0\).) So template matching is a one-layer convolutional network whose single kernel nobody trained: it is a copy of one earthquake. We will hold on to that sentence in A.7.
The denominator changes along the record. \(\sigma_d(\tau)\) is the standard deviation of the \(n\) record samples under the template at start \(\tau\). Computing it afresh for every \(\tau\) costs \(n\) operations each; instead we use cumulative sums. Let \(S_1(m) = \sum_{k \le m} d_k\) and \(S_2(m) = \sum_{k \le m} d_k^2\). The sum over any window is a difference of two cumulative sums, \(\sum_{k=\tau+1}^{\tau+n} d_k = S_1(\tau+n) - S_1(\tau)\), so
two subtractions per \(\tau\) whatever \(n\) is. The numerator for every \(\tau\) at once is computed by the fast Fourier transform (FFT): a convolution of a record of \(N\) samples with a kernel of \(n\) samples costs about \(N n\) operations done directly, and about \(N \log N\) through the FFT.
Below is the whole method in one function, ncc — the name and the conventions of CCTorch, our
group's correlation code. Where the record is perfectly flat (a gap filled with zeros), the
standard deviation is zero and we return 0 rather than divide by it. We check it against
np.corrcoef at every lag, and against a PyTorch convolution layer whose kernel we set to the
template.
def ncc(template, data):
"""Normalised cross-correlation of a template with every same-length stretch of `data`.
Output length len(data) - len(template) + 1, values in [-1, 1]."""
n = len(template)
if template.std() == 0: # a dead channel matches nothing
return np.zeros(len(data) - n + 1)
t = (template - template.mean()) / (template.std() * n)
num = scipy.signal.fftconvolve(data, t[::-1], mode="valid") # the sliding inner product
s1, s2 = np.cumsum(np.r_[0.0, data]), np.cumsum(np.r_[0.0, data ** 2])
var = (s2[n:] - s2[:-n]) / n - ((s1[n:] - s1[:-n]) / n) ** 2 # running variance
sd = np.sqrt(np.clip(var, 0, None))
return np.where(sd > 0, num / np.where(sd > 0, sd, 1.0), 0.0)
tpl = piece(i0)
record = np.ascontiguousarray(Xf[near, SQK, 2])
c = ncc(tpl, record)
direct = np.array([np.corrcoef(tpl, record[k:k + len(tpl)])[0, 1] for k in range(len(c))])
conv = nn.Conv1d(1, 1, len(tpl), bias=False).double()
conv.weight.data[:] = torch.tensor((tpl - tpl.mean()) / (tpl.std() * len(tpl)))[None, None]
num_torch = conv(torch.tensor(record)[None, None]).detach().numpy()[0, 0]
num_fft = scipy.signal.fftconvolve(record, ((tpl - tpl.mean()) / (tpl.std() * len(tpl)))[::-1], mode="valid")
print(f"largest difference from np.corrcoef: {np.abs(c - direct).max():.1e}")
print(f"largest difference between the FFT numerator and torch's Conv1d: {np.abs(num_fft - num_torch).max():.1e}")
largest difference from np.corrcoef: 8.2e-13 largest difference between the FFT numerator and torch's Conv1d: 1.1e-14
The differences are rounding error. (PyTorch's Conv1d, like most deep-learning code, slides
the kernel without flipping it, which a mathematician would call a correlation; fftconvolve
flips, so we flip the template ourselves.)
How much does the FFT save? We time the numerator — the sliding inner product — of one template against 4.4 hours of record (the 265 minutes at SQK strung together), directly and by FFT.
def join(records, taper=0.5):
"""String one-minute records together, each demeaned and tapered at both ends, so the joints
make no steps."""
w = scipy.signal.windows.tukey(records.shape[-1], 2 * taper * FS / records.shape[-1])
return np.concatenate(list((records - records.mean(-1, keepdims=True)) * w), axis=-1)
cont = bandpass(join(REC)) # (E N Z, 265 x 6000 samples)
print(f"continuous stand-in: {cont.shape[-1] / FS / 3600:.1f} h at {FS} Hz on 3 components")
m = 200_000 # time the direct sum on part, then scale
for n in (100, 1000): # a 1 s and a 10 s template
t_ = cont[2, 50_000:50_000 + n]
t0 = time.perf_counter(); scipy.signal.fftconvolve(cont[2], t_[::-1], mode="valid"); t_fft = time.perf_counter() - t0
t0 = time.perf_counter(); np.correlate(cont[2, :m + n - 1], t_, "valid") # the same numerator, directly
t_dir = (time.perf_counter() - t0) * cont.shape[-1] / m
print(f"a {n / FS:.0f} s template over {cont.shape[-1] / FS / 3600:.1f} h: FFT {t_fft:.2f} s, direct sum {t_dir:.2f} s")
continuous stand-in: 4.4 h at 100 Hz on 3 components a 1 s template over 4.4 h: FFT 0.03 s, direct sum 0.02 s a 10 s template over 4.4 h: FFT 0.03 s, direct sum 0.15 s
For a short template the direct sum is as fast as the FFT, because \(n\) is small; for a long one the FFT wins, because its cost hardly depends on the template's length while the direct sum grows in proportion to it. Real scans correlate thousands of templates with years of data from hundreds of channels; QTM, the Southern California scan, did it on graphics processors (Ross et al. 2019).
A.3 · Which earthquakes look alike?¶
A.1 showed two earthquakes less than 100 m apart that look alike and one under a kilometre away that does not. Why? A seismogram is the source filtered by the path to the station. Two earthquakes at nearly the same place, with the same mechanism, share the path; their seismograms differ by the source time function, which for these small events is shorter than the periods we look at. As the two move apart, the paths differ, and the differences that matter are those comparable to a wavelength \(\lambda = v / f\), the wave speed \(v\) divided by the frequency \(f\). Geller & Mueller (1980) found four central California earthquakes whose seismograms were nearly identical below 5 Hz but not above, and argued that they were within about a quarter of a wavelength of each other. How long is a wavelength here? The P travel times of the ball give the speed: we fit a straight line to the analyst's P times (after the origin) against the distance from each earthquake to each station.
hypo = np.sqrt(((stations.station_longitude.values[None] - patch.event_longitude.values[:, None]) * gd.KM_LON) ** 2
+ ((stations.station_latitude.values[None] - patch.event_latitude.values[:, None]) * gd.KM_LAT) ** 2
+ (patch.event_depth_km.values[:, None] + stations.station_elevation_m.values[None] / 1000) ** 2)
has = P >= 0
tt = P / FS - PRE # P travel time (s)
slope, intercept = np.polyfit(hypo[has], tt[has], 1)
vp = 1 / slope
print(f"apparent P speed from {has.sum():,} picks: {vp:.1f} km/s")
for f in (2, 5, 15, 40):
print(f" {f:2d} Hz: wavelength {1000 * vp / f:5.0f} m, a quarter of it {250 * vp / f:4.0f} m")
apparent P speed from 2,491 picks: 5.3 km/s 2 Hz: wavelength 2644 m, a quarter of it 661 m 5 Hz: wavelength 1058 m, a quarter of it 264 m 15 Hz: wavelength 353 m, a quarter of it 88 m 40 Hz: wavelength 132 m, a quarter of it 33 m
A quarter wavelength is a few hundred metres at the lowest frequencies we use and a few tens of metres at the highest. So two earthquakes a few hundred metres apart should still look alike at 2 Hz but not at 40 Hz.
We test this on all pairs of the 265 earthquakes. For a pair \((i, j)\) we correlate 3 s of earthquake \(i\) — from 0.1 s before P, on all three components — with earthquake \(j\), allowing a shift of up to 0.3 s on either side of \(j\)'s P pick (the picks are not exact, and that is the point of A.6), and keep the largest value. Each of the three components gives a correlation; we average them, and then average over the 10 stations. We do this in three frequency bands.
To read the numbers we need a reference for "no similarity at all". A correlation that is the maximum over 61 shifts is not near zero even between unrelated signals, because the largest of many random values is large. So we repeat the computation with every template reversed in time: a reversed seismogram has the same frequencies and the same length as the original, but it cannot look like any earthquake. Its correlations show the level that similarity must beat.
def windows_at(Xb, s, ok, length=300, lag=30):
"""Templates (from 0.1 s before P, `length` samples) and search windows (+-`lag` samples) at station s."""
T = np.stack([Xb[i, s, :, P[i, s] - 10:P[i, s] - 10 + length] for i in ok])
W = np.stack([Xb[j, s, :, P[j, s] - 10 - lag:P[j, s] - 10 + length + lag] for j in ok])
return T, W
def unit(a):
a = a - a.mean(-1, keepdims=True)
nrm = np.linalg.norm(a, axis=-1, keepdims=True)
return np.where(nrm > 0, a / np.where(nrm > 0, nrm, 1), 0.0) # a dead channel stays zero
def pair_cc(Xb, length=300, lag=30, reverse=False, channels=(0, 1, 2), stations_=range(NS)):
"""Largest correlation over +-lag for every pair, averaged over components and stations."""
total, count = np.zeros((NE, NE)), np.zeros((NE, NE))
for s in stations_:
ok = np.where(P[:, s] >= 10 + lag)[0]
T, W = windows_at(Xb, s, ok, length, lag)
T, W = T[:, channels], W[:, channels]
if reverse:
T = T[..., ::-1]
Wv = unit(np.lib.stride_tricks.sliding_window_view(W, length, axis=-1)) # (j, comp, shift, length)
J, C, K, L = Wv.shape
# one matrix product: every template against every window at every shift, summed over components
cc = unit(T).reshape(len(ok), C * L) @ Wv.transpose(0, 2, 1, 3).reshape(J * K, C * L).T
cc = cc.reshape(len(ok), J, K).max(-1) / len(channels)
total[np.ix_(ok, ok)] += np.minimum(cc, cc.T) # symmetric: the worse of i->j and j->i
count[np.ix_(ok, ok)] += 1
out = total / np.where(count > 0, count, 1)
np.fill_diagonal(out, np.nan)
return out
bands = {"2-8 Hz": (2, 8), "5-40 Hz": (5, 40), "15-40 Hz": (15, 40)}
cc_band, null_band = {}, {}
for name, (lo, hi) in bands.items():
Xb = bandpass(X, lo, hi)
cc_band[name], null_band[name] = pair_cc(Xb), pair_cc(Xb, reverse=True)
iu = np.triu_indices(NE, 1)
edges = np.array([0, 0.05, 0.1, 0.2, 0.3, 0.5, 0.75, 1.0, 1.5, 2.0])
mid = 0.5 * (edges[1:] + edges[:-1])
fig, ax = plt.subplots(figsize=(6.5, 3.6))
for (name, cc), col in zip(cc_band.items(), ("tab:blue", "k", "tab:red")):
k = np.digitize(sep[iu], edges) - 1
med = [np.median(cc[iu][k == b]) if (k == b).sum() >= 10 else np.nan for b in range(len(mid))]
ax.plot(mid * 1000, med, "o-", c=col, label=f"{name}: earthquake pairs (median)")
ax.axhline(np.nanquantile(null_band[name][iu], 0.99), c=col, ls=":", lw=1)
ax.set(xlabel="distance between the two earthquakes (m)", ylabel="correlation, 10 stations",
title="Nearby earthquakes look more alike; dotted: 99th percentile of reversed templates")
ax.legend(fontsize=7)
plt.show()
for name in bands:
close, apart = cc_band[name][iu][sep[iu] < 0.1], cc_band[name][iu][sep[iu] > 1.0]
print(f"{name:>9}: median correlation {np.median(close):.2f} for pairs < 100 m, {np.median(apart):.2f} beyond 1 km; "
f"reversed templates: 99th percentile {np.nanquantile(null_band[name][iu], 0.99):.2f}")
2-8 Hz: median correlation 0.51 for pairs < 100 m, 0.38 beyond 1 km; reversed templates: 99th percentile 0.30 5-40 Hz: median correlation 0.30 for pairs < 100 m, 0.22 beyond 1 km; reversed templates: 99th percentile 0.17 15-40 Hz: median correlation 0.25 for pairs < 100 m, 0.21 beyond 1 km; reversed templates: 99th percentile 0.16
Three things are in the figure.
- Closer pairs look more alike in every band, as the wavelength argument says, and the correlation falls off gradually rather than at one distance. Part of that gradual fall may be the catalogue's own location errors, which blur the distances on the horizontal axis.
- The lowest band is the most similar at every distance, again as the argument says: its wavelengths are longest. But its reversed templates also correlate more. A narrow band has few independent wiggles in 3 s, so even unrelated signals line up by chance. What matters for a detector is the distance between the earthquake curve and the dotted line, not the height of the curve. We return to this in A.4.
- How fast the similarity falls depends on the band, as the wavelength argument predicts. At 2–8 Hz it falls steadily over the whole kilometre; at 15–40 Hz it is already low for the closest pairs and hardly falls further, because a quarter wavelength there (tens of metres) is smaller than almost every separation we have, and smaller than the catalogue's location errors.
- Every curve lies above its dotted line, even beyond 1 km. Earthquakes of one ball share much of their path to each station, so even distant pairs keep some similarity, and the average over 10 stations pushes the reference down (A.5).
Correlation therefore reaches only as far as the earthquakes that resemble the template. Session 4, Clustering, found that Geysers earthquakes are packed tightly in space; this is why template matching works well here, and why it can only ever find earthquakes near the ones we already know.
A.4 · From a correlation to a detection¶
Now we scan. A template is the seismogram of one catalogued earthquake at SQK, cut as Part B also cuts it: P on the vertical component and S on the two horizontals, each 1 s long from 0.1 s before the analyst's pick. The three correlation traces are shifted so that each is indexed by the time of P, and averaged. Where an earthquake like the template passes, all three peak together.
We scan the 4.4 hours at SQK with every one of the 138 earthquakes whose P and S were both picked there. Each scan returns one correlation trace, which is mostly noise with an occasional peak.
def template_at(i, s=SQK):
"""(piece, component, offset from P in samples) for P on Z and S on E and N at station s."""
p, q = P[i, s], S[i, s]
return [(Xf[i, s, 2, p - 10:p + 90], 2, 0), (Xf[i, s, 0, q - 10:q + 90], 0, q - p),
(Xf[i, s, 1, q - 10:q + 90], 1, q - p)]
def scan(template, data):
"""Average of the component correlation traces, each shifted to the time of P."""
n = data.shape[-1] - 100 - max(off for _, _, off in template)
return np.mean([ncc(t, data[comp])[off:off + n] for t, comp, off in template], axis=0)
def reverse(template):
return [(t[::-1].copy(), comp, off) for t, comp, off in template]
with_s = np.where((P[:, SQK] >= 10) & (S[:, SQK] > 0))[0]
p_in_cont = np.arange(NE) * REC.shape[-1] + OFFSET + P[:, SQK] - 10 # the scan is indexed by P - 0.1 s
first = with_s[0]
c = scan(template_at(first), cont)
hours = np.arange(len(c)) / FS / 3600
fig, ax = plt.subplots(2, 1, figsize=(8, 4.2))
ax[0].plot(hours, c, c="0.3", lw=0.3)
ax[0].plot(hours[p_in_cont[first]], c[p_in_cont[first]], "v", c="tab:red", ms=6, label="the template's own earthquake")
ax[0].set(xlabel="hours", ylabel="correlation", title="One template scanned along 4.4 hours at SQK")
ax[0].legend(fontsize=7, loc="lower right")
k0 = p_in_cont[first]
z = slice(k0 - 3000, k0 + 3000)
ax[1].plot((np.arange(len(c))[z] - k0) / FS, c[z], c="0.3", lw=0.6)
ax[1].set(xlabel="seconds from the template's own P", ylabel="correlation")
plt.tight_layout()
plt.show()
print(f"at the template's own P the correlation is {c[k0]:.2f}; its median over the record is "
f"{np.median(c):.3f} and its standard deviation {c.std():.3f}")
at the template's own P the correlation is 1.00; its median over the record is 0.000 and its standard deviation 0.098
The template finds itself with a correlation of 1, as it must, and a few other earthquakes stand out as tall spikes. Everywhere else the trace wanders about zero. (The two narrow gaps are minutes in which SQK recorded nothing; there the correlation is set to 0.) To call a peak a detection we need a threshold, and the threshold must be set by the noise: the same template scanned over a noisier station or a narrower band gives a wider trace.
The usual measure of the trace's spread is not the standard deviation but the median absolute deviation, \(\mathrm{MAD} = \mathrm{median}(|c - \mathrm{median}(c)|)\). The standard deviation squares every value, so the few large peaks we are looking for inflate it; the median ignores them. For a Gaussian, \(\mathrm{MAD} = 0.6745\,\sigma\), so \(\sigma \approx 1.4826\,\mathrm{MAD}\). We call a peak a detection when it exceeds 8 × MAD above the median and is the largest value within 2 s on either side, so that one earthquake gives one detection. Part B uses the same rule.
Is 8 × MAD a strict threshold? If the trace were Gaussian, it would be \(8 / 1.4826 = 5.4\sigma\), which Gaussian noise exceeds very rarely (the cell below prints how rarely). We check that against the null: the same scans with every template reversed in time, as in A.3. Reversed templates cannot match any earthquake, so every detection they make is a false alarm.
def peaks(trace, threshold, half=2.0):
"""Indices of local maxima above `threshold` that are the largest within `half` seconds."""
top = scipy.ndimage.maximum_filter1d(trace, 2 * int(half * FS) + 1)
return np.where((top == trace) & (trace > threshold))[0]
origin_in_cont = np.arange(NE) * REC.shape[-1] + OFFSET + int(PRE * FS)
rows, z_true, z_null = [], [], []
for i in with_s:
for kind, tp in (("template", template_at(i)), ("reversed", reverse(template_at(i)))):
c = scan(tp, cont)
med = np.median(c); mad = np.median(np.abs(c - med))
(z_true if kind == "template" else z_null).append((c - med) / mad)
for k in peaks(c, med + 8 * mad):
nearest_p = np.min(np.abs(p_in_cont[P[:, SQK] >= 10] - k))
rec, at = divmod(k + 10, REC.shape[-1]) # which minute, and where in it
before = at < OFFSET[rec] # before that minute's earthquake
where = ("its own earthquake" if kind == "template" and abs(k - p_in_cont[i]) <= 20 else
"another catalogued P" if nearest_p <= 20 else
"before the catalogued earthquake" if before else "after a catalogued origin")
rows.append((kind, i, where, c[k], k))
det = pd.DataFrame(rows, columns=["kind", "template", "where", "cc", "sample"])
table = det.groupby(["where", "kind"]).size().unstack(fill_value=0)
print(f"{len(with_s)} templates, each scanned over {cont.shape[-1] / FS / 3600:.1f} h; detections above 8 x MAD:")
print(table.loc[["its own earthquake", "another catalogued P", "after a catalogued origin", "before the catalogued earthquake"]
if "its own earthquake" in table.index else table].to_string())
138 templates, each scanned over 4.4 h; detections above 8 x MAD: kind reversed template where its own earthquake 0 138 another catalogued P 0 795 after a catalogued origin 19 113 before the catalogued earthquake 6 81
Read the table row by row.
- Its own earthquake is found by every template, as it must be.
- Another catalogued P: these are the neighbours, the copies template matching exists to find. On average each template finds a handful of the other 264 earthquakes in the ball — the few that resemble it closely, as A.3 predicts.
- After a catalogued origin: a peak in the 20 s after a catalogued earthquake's origin but not at its P — for example the template's P window lining up with the S, or with the coda, the long tail of scattered waves after S, which comes from the same source region as the template and so resembles it a little. These are mostly not new earthquakes. The reversed templates make far fewer of them, because a reversed template shares neither source nor path: the null measures false alarms caused by noise, not those caused by other earthquakes. Part B meets this on a large scale in the burst of 23:59.
- Before the catalogued earthquake: peaks in the 39 s of background before each minute's catalogued earthquake. The reversed templates make a few there too, so some of these are false alarms; the excess over the reversed column are candidates for earthquakes the catalogue does not contain. A.7 tests them a second way, with PhaseNet.
Now the rate of false alarms. We pool the scans and compare the tail of the reversed traces with the tail a Gaussian would have.
zn, zt = np.concatenate(z_null), np.concatenate(z_true)
levels = np.arange(2, 12.5, 0.5)
gauss = scipy.stats.norm.sf(levels / 1.4826)
fig, ax = plt.subplots(figsize=(6, 3.4))
ax.semilogy(levels, [(zn > l).mean() for l in levels], "o-", c="0.3", ms=3, label="reversed templates (the null)")
ax.semilogy(levels, [(zt > l).mean() for l in levels], "o-", c="tab:red", ms=3, label="templates")
ax.semilogy(levels, gauss, "--", c="tab:blue", label="Gaussian noise")
ax.axvline(8, c="k", lw=0.8)
ax.set(xlabel="threshold (multiples of MAD above the median)", ylabel="fraction of samples above",
ylim=(1e-9, 1), title="Beyond about 5 x MAD the null is heavier-tailed than a Gaussian")
ax.legend(fontsize=7)
plt.show()
per_scan_hours = cont.shape[-1] / FS / 3600
g8 = scipy.stats.norm.sf(8 / 1.4826)
print(f"Gaussian noise exceeds 8 x MAD once in {1 / g8:.1e} samples, once in {1 / g8 / FS / 86400:.1f} days at {FS} Hz")
false = det[(det.kind == "reversed")].shape[0] / len(with_s) / per_scan_hours
for l in (6, 8, 10, 12):
print(f"samples above {l:2d} x MAD: reversed {np.mean(zn > l):.1e}, Gaussian {scipy.stats.norm.sf(l / 1.4826):.1e}")
print(f"at 8 x MAD the null exceeds the Gaussian by a factor of {np.mean(zn > 8) / scipy.stats.norm.sf(8 / 1.4826):.1f}")
print(f"false alarms per template per hour (reversed templates, 2 s apart): {false:.3f}")
Gaussian noise exceeds 8 x MAD once in 2.9e+07 samples, once in 3.4 days at 100 Hz
samples above 6 x MAD: reversed 3.2e-05, Gaussian 2.6e-05
samples above 8 x MAD: reversed 1.5e-07, Gaussian 3.4e-08 samples above 10 x MAD: reversed 0.0e+00, Gaussian 7.7e-12
samples above 12 x MAD: reversed 0.0e+00, Gaussian 2.9e-16 at 8 x MAD the null exceeds the Gaussian by a factor of 4.4 false alarms per template per hour (reversed templates, 2 s apart): 0.041
Up to about 5 × MAD the null follows the Gaussian; beyond it, the null's tail is heavier — at 8 × MAD by the factor printed above. Real noise is not a steady Gaussian process: it comes in bursts, and correlation traces inherit them. The templates' own curve lies far above both: those extra samples are the earthquakes. A threshold in MAD units is therefore only a convention, and its false-alarm rate has to be measured with a null and quoted with the detections. At 8 × MAD the rate is small for one template at one station, but a real scan runs thousands of templates for years, and the false alarms add up in proportion.
A.5 · Many stations: why the network stack works¶
A.3 averaged over 10 stations and three components; A.4 over three components at one station. How much does averaging help? Take the correlation of a template with an unrelated stretch of noise at one channel: it is a random number with mean near zero and some spread \(\sigma_1\). On \(N\) channels whose noise is independent, the average of the \(N\) correlations has spread \(\sigma_1 / \sqrt{N}\) — the standard error of a mean, Session 1. A real copy of the template, by contrast, correlates on every channel at once, so its average stays where it was. The gap between earthquake and noise widens by \(\sqrt{N}\). Gibbons & Ringdal (2006) built a detector on exactly this: correlation traces stacked over an array.
We measure it on the 265 earthquakes of the ball. For each earthquake we find its most similar partner, the "best partner", and compare that correlation with the 99.9th percentile of the reversed templates, for one channel (SQK vertical, P window of 3 s), three channels (SQK), and 30 channels (all 10 stations).
setups = {"1 channel (SQK Z)": dict(channels=(2,), stations_=[SQK]),
"3 channels (SQK)": dict(stations_=[SQK]),
"30 channels (10 stations)": dict()}
fig, ax = plt.subplots(figsize=(6.5, 3.4))
for k, (name, kw) in enumerate(setups.items()):
cc, nl = pair_cc(Xf, **kw), pair_cc(Xf, reverse=True, **kw)
best = np.nanmax(cc, 1)
best = best[np.isfinite(best)]
thr = np.nanquantile(nl[iu], 0.999)
ax.boxplot(best, positions=[k], widths=0.5, showfliers=False)
ax.plot([k - 0.35, k + 0.35], [thr, thr], c="tab:red", lw=2)
print(f"{name:>26}: null spread {np.nanstd(nl[iu]):.3f}, null 99.9 % {thr:.2f}; "
f"best partner median {np.median(best):.2f}; earthquakes whose best partner beats the null: "
f"{(best > thr).mean():.0%}")
ax.set_xticks(range(len(setups)), list(setups), fontsize=8)
ax.set(ylabel="correlation", title="Each earthquake's best partner (boxes) against the null (red)")
plt.show()
1 channel (SQK Z): null spread 0.055, null 99.9 % 0.40; best partner median 0.57; earthquakes whose best partner beats the null: 97%
3 channels (SQK): null spread 0.034, null 99.9 % 0.25; best partner median 0.54; earthquakes whose best partner beats the null: 99%
30 channels (10 stations): null spread 0.016, null 99.9 % 0.19; best partner median 0.49; earthquakes whose best partner beats the null: 100%
The boxes come down a little from left to right: distant stations see the pair less alike than SQK does, so their correlations pull the average down. But the null comes down much faster. From one channel to three its spread falls by close to \(\sqrt{3}\), as independent noise predicts; from three to thirty it falls less than \(\sqrt{10}\) would say, because noise on nearby stations is not fully independent. What matters is the gap between box and red line, and it widens with every channel added: with 30 channels even the least similar best partner is above the null. Stacking needs one thing we have passed over: the correlation traces must be aligned before they are averaged, each shifted by the travel time from the source to its station, the template's moveout. Part B does this for every station of the array, using each template's own picks; A.7 does it for our 10.
A.6 · Timing by correlation, to a fraction of a sample¶
A correlation peak says not only whether two earthquakes are alike but when one is relative to the other: the lag of the peak is the difference of their arrival times. Session 2, Location & relocation, used such differences between pairs of events — differential times — for double difference relocation, but took them from the analysts' picks. How much better are correlation lags?
A lag is a whole number of samples, 10 ms at 100 Hz. We can do better by fitting a parabola through the peak value \(c_0\) and its two neighbours \(c_{-1}\) and \(c_{+1}\). Write the parabola as \(c(x) = a x^2 + b x + c_0\), with \(x\) in samples from the peak sample. Putting in \(x = -1\) and \(x = +1\):
Subtracting the first from the second gives \(b = (c_{+1} - c_{-1})/2\); adding them gives \(a = (c_{+1} + c_{-1})/2 - c_0\). The parabola peaks where its slope \(2 a x + b\) is zero:
a fraction of a sample between \(-\tfrac12\) and \(\tfrac12\). (Poupinet et al. (1984) reached a resolution an order of magnitude finer than the sampling interval with a related method, which works with the phase of the cross-spectrum in the frequency domain.)
First a check that the parabola recovers a shift we plant: we move seismograms by a known fraction of a sample (by multiplying their Fourier transform by \(e^{-2\pi i f \delta}\)) and measure the shift back.
We measure, at each station, the differential P time of every pair whose correlation there exceeds 0.7: the template is earthquake \(i\)'s P window (1 s, vertical), and the lag is searched within 0.3 s of earthquake \(j\)'s P pick.
def parabola(cm, c0, cp):
"""Offset (in samples, between -1/2 and 1/2) of the peak of the parabola through three values."""
den = cm - 2 * c0 + cp
return np.where(den < 0, 0.5 * (cm - cp) / np.where(den < 0, den, -1), 0.0)
ok = np.where(P[:, SQK] >= 40)[0]
true = rng.uniform(-0.5, 0.5, len(ok)) # planted shifts, in samples
err = []
for i, d in zip(ok, true):
x = Xf[i, SQK, 2]
f = np.fft.rfftfreq(len(x))
y = np.fft.irfft(np.fft.rfft(x) * np.exp(-2j * np.pi * f * d), len(x)) # x moved later by d samples
t_, w = unit(x[P[i, SQK] - 10:P[i, SQK] + 90]), y[P[i, SQK] - 40:P[i, SQK] + 120]
c = unit(np.lib.stride_tricks.sliding_window_view(w, 100)) @ t_
k = int(np.argmax(c))
err.append(k + parabola(c[k - 1], c[k], c[k + 1]) - 30 - d)
print(f"planted shifts recovered with error {np.std(err) * 10:.2f} ms (rms), {len(ok)} seismograms")
planted shifts recovered with error 0.25 ms (rms), 263 seismograms
dt_rows = []
for s in range(NS):
ok = np.where(P[:, s] >= 40)[0]
T = unit(np.stack([Xf[i, s, 2, P[i, s] - 10:P[i, s] + 90] for i in ok])) # (i, 100)
W = np.stack([Xf[j, s, 2, P[j, s] - 40:P[j, s] + 120] for j in ok]) # (j, 160)
Wv = unit(np.lib.stride_tricks.sliding_window_view(W, 100, axis=-1)) # (j, 61 shifts, 100)
cc = (T @ Wv.reshape(-1, 100).T).reshape(len(ok), len(ok), -1) # (i, j, shift)
k = np.clip(cc.argmax(-1), 1, cc.shape[-1] - 2)
take = lambda d: np.take_along_axis(cc, (k + d)[..., None], -1)[..., 0]
peak, lag = take(0), k + parabola(take(-1), take(0), take(1))
a_, b_ = np.triu_indices(len(ok), 1)
keep = peak[a_, b_] > 0.7
for x, y in zip(a_[keep], b_[keep]):
i, j = ok[x], ok[y]
# arrival of j minus arrival of i, each relative to its own origin (s): by correlation and by picks
dt_rows.append((s, i, j, peak[x, y], (P[j, s] - 40 + lag[x, y] - (P[i, s] - 10)) / FS, (P[j, s] - P[i, s]) / FS))
dts = pd.DataFrame(dt_rows, columns=["station", "i", "j", "cc", "dt_cc", "dt_pick"])
dts["diff_ms"] = 1000 * (dts.dt_cc - dts.dt_pick)
print(f"{len(dts):,} pairs-at-a-station with correlation above 0.7, from {dts[['i', 'j']].stack().nunique()} earthquakes")
print(f"correlation minus pick differential time: median {dts.diff_ms.median():.1f} ms, "
f"middle 68 % within +-{(dts.diff_ms.quantile(0.84) - dts.diff_ms.quantile(0.16)) / 2:.1f} ms")
2,505 pairs-at-a-station with correlation above 0.7, from 258 earthquakes correlation minus pick differential time: median -0.7 ms, middle 68 % within +-47.9 ms
The two measures agree on average but differ, pair by pair, by tens of milliseconds. Which is right? There is no true time to compare with, but there is a consistency test. For three earthquakes \(i\), \(j\), \(k\) at one station, the differential times must close: \(\Delta t_{ij} + \Delta t_{jk} = \Delta t_{ik}\). Differential times taken from picks close exactly, whatever the picks' errors, because each pick enters twice with opposite signs, so the test says nothing about picks. Correlation lags are measured pair by pair, independently, so their closure error \(\Delta t_{ij} + \Delta t_{jk} - \Delta t_{ik}\) measures their precision.
key = {(r.station, r.i, r.j): r.dt_cc for r in dts.itertuples()}
closure = []
for (s, i, j), d_ij in key.items():
for (s2, j2, k), d_jk in list(key.items()):
if s2 == s and j2 == j and (s, i, k) in key:
closure.append(1000 * (d_ij + d_jk - key[(s, i, k)]))
closure = np.array(closure)
spread_cc = (np.quantile(closure, 0.84) - np.quantile(closure, 0.16)) / 2
spread_pick = (dts.diff_ms.quantile(0.84) - dts.diff_ms.quantile(0.16)) / 2
fig, ax = plt.subplots(figsize=(6, 3.2))
bins = np.linspace(-60, 60, 121)
ax.hist(dts.diff_ms, bins, color="0.6", label="correlation minus picks")
ax.hist(closure, bins, histtype="step", color="tab:red", lw=1.5, label="closure error of correlation lags")
ax.set(xlabel="milliseconds", ylabel="count", title="Correlation lags are an order of magnitude more precise than picks")
ax.legend(fontsize=7)
plt.show()
print(f"{len(closure):,} triangles; closure error of correlation lags: middle 68 % within +-{spread_cc:.1f} ms "
f"(about +-{spread_cc / np.sqrt(3):.1f} ms per lag if the three errors are independent)")
print(f"correlation minus picks: +-{spread_pick:.1f} ms, so the picks' own errors are of that size")
5,115 triangles; closure error of correlation lags: middle 68 % within +-2.7 ms (about +-1.5 ms per lag if the three errors are independent) correlation minus picks: +-47.9 ms, so the picks' own errors are of that size
Three correlation lags close to within a couple of milliseconds — a fraction of the 10 ms sample — while they differ from the picks by many times that. Closure measures precision, not accuracy: a bias shared by every lag, or the same wrong cycle chosen consistently, would still close. The planted shifts above check accuracy for the parabola itself. Together they say that most of the scatter between correlation and picks is pick error. This is why relocation programs feed on correlation differential times: Schaff & Waldhauser (2005) measured them for about 225,000 northern California earthquakes, and HypoDD (Session 2) uses them in place of picks wherever events correlate. Families of correlated earthquakes can then be found from the matrix of pair correlations with the hierarchical linkage of Session 4; Part B does so for the hour.
A.7 · Template matching against a neural network¶
The answer first. Template matching finds a copy of its template at a signal several times weaker than PhaseNet needs, at any threshold; it does not find earthquakes unlike the template at all, which PhaseNet finds as easily as copies. The rest of this section shows how we know.
We now have two ways of finding earthquakes in a continuous record. Template matching (A.2–A.6) slides one known seismogram along the record. A picker such as PhaseNet (Session 5) is a convolutional network that learned from analyst picks across northern California what P and S arrivals look like in general. A.2 showed that template matching is itself a convolution with one fixed kernel, so the two differ in what the kernels are: one copy of one earthquake, against many small kernels learned in layers from labels. Which finds more earthquakes, and which earthquakes?
On real data we cannot answer this cleanly, because we do not know every earthquake that is there (Part B meets exactly that difficulty). So we plant earthquakes where we know the answer. We take real noise, the 37 s before each catalogued earthquake at SQK, and add to it at a random time a real seismogram of an earthquake, scaled down to a chosen signal-to-noise ratio (SNR): the peak amplitude of the vertical component in the 2 s after P, divided by the standard deviation of the noise's vertical component, both in the 5–40 Hz band. Then we ask both methods to find it.
The planted earthquake is one of four kinds, relative to a template earthquake \(i\):
- the template itself — earthquake \(i\) in new noise, the best case for template matching;
- its closest look-alike — the earthquake of the ball most similar to \(i\) at SQK;
- a typical neighbour — a random other earthquake of the ball, less than 2 km from \(i\);
- an earthquake elsewhere — one of the 30 earthquakes from other parts of the field.
Template matching scans with earthquake \(i\)'s template and the rules of A.4 (detection if the correlation within 0.2 s of the planted P exceeds 8 × MAD of that trace). PhaseNet runs on the raw three components with the published weights and the pipeline of Session 5 (detection if it picks P with probability at least 0.5 within 0.5 s of the planted P). Both see exactly the same 37 s of data. (PhaseNet looks for P within 0.5 s and template matching within 0.2 s; PhaseNet's looser window only helps it.) First the network, the same as in Session 5, Part B; you need not read the class again — it is PhaseNet's U-Net with the padding rules its published weights need.
def standardise(x, window=3000):
"""Running mean and standard deviation per component (PhaseNet's normalize_long)."""
nt, shift = x.shape[-1], window // 2
xp = np.pad(x.astype(np.float64), ((0, 0), (shift, shift)), mode="reflect")
c1, c2 = (np.pad(np.cumsum(a, -1), ((0, 0), (1, 0))) for a in (xp, xp * xp))
s = np.arange(0, xp.shape[-1] - window + 1, shift) # a window every `shift` samples
mean = (c1[:, s + window] - c1[:, s]) / window
std = np.sqrt(np.clip((c2[:, s + window] - c2[:, s]) / window - mean ** 2, 0, None))
k = (nt + shift - 1) // shift
knot = np.arange(k + 1) * shift; knot[-1] = nt; use = np.clip(np.arange(k + 1), 1, k - 1)
m = np.stack([np.interp(np.arange(nt), knot, a[use]) for a in mean])
sd = np.stack([np.interp(np.arange(nt), knot, a[use]) for a in std])
live = int((sd.sum(-1) > 0).sum()) # a dead component counts as missing:
gain = 3.0 / live if 0 < live < 3 else 1.0 # the live ones are scaled up to make up
return ((x - m) / np.where(sd == 0, 1, sd) * gain).astype(np.float32)
class PhaseNet2018(nn.Module):
"""PhaseNet's U-Net (Session 5's UNet2018 with TensorFlow's 'same' padding); logits (batch, 3, samples)."""
def __init__(self, depths=5, filters_root=8, kernel_size=7, pool_size=4):
super().__init__()
K, S = kernel_size, pool_size
self.K, self.S, self.depths = K, S, depths
self.input_conv = nn.Conv1d(3, filters_root, K, padding=K // 2)
self.input_bn = nn.BatchNorm1d(filters_root, eps=1e-3)
self.down_conv1, self.down_bn1 = nn.ModuleList(), nn.ModuleList()
self.down_conv3, self.down_bn3 = nn.ModuleList(), nn.ModuleList()
c = filters_root
for d in range(depths):
f = 2 ** d * filters_root # 8, 16, 32, 64, 128 channels
self.down_conv1.append(nn.Conv1d(c, f, K, padding=K // 2, bias=False)); self.down_bn1.append(nn.BatchNorm1d(f, eps=1e-3))
if d < depths - 1:
self.down_conv3.append(nn.Conv1d(f, f, K, stride=S, bias=False)); self.down_bn3.append(nn.BatchNorm1d(f, eps=1e-3))
c = f
self.up_conv0, self.up_bn0 = nn.ModuleList(), nn.ModuleList()
self.up_conv1, self.up_bn1 = nn.ModuleList(), nn.ModuleList()
for d in range(depths - 2, -1, -1):
f = 2 ** d * filters_root
self.up_conv0.append(nn.ConvTranspose1d(c, f, K, stride=S, bias=False)); self.up_bn0.append(nn.BatchNorm1d(f, eps=1e-3))
self.up_conv1.append(nn.Conv1d(2 * f, f, K, padding=K // 2, bias=False)); self.up_bn1.append(nn.BatchNorm1d(f, eps=1e-3))
c = f
self.output_conv = nn.Conv1d(filters_root, 3, 1)
def forward(self, x):
K, S = self.K, self.S
net = F.relu(self.input_bn(self.input_conv(x)))
skips = []
for d in range(self.depths):
net = F.relu(self.down_bn1[d](self.down_conv1[d](net)))
skips.append(net)
if d < self.depths - 1:
r = net.shape[-1] % S # rule 1: TensorFlow pads K - S
pad = K - S if r == 0 else max(K - r, 0) # samples (or K - r), the extra
net = F.pad(net, (pad // 2, pad - pad // 2)) # one on the right
net = F.relu(self.down_bn3[d](self.down_conv3[d](net)))
for i, d in enumerate(range(self.depths - 2, -1, -1)):
n = net.shape[-1]
net = self.up_conv0[i](net)[..., (K - S) // 2:(K - S) // 2 + n * S] # rule 2: keep n * S samples
net = F.relu(self.up_bn0[i](net))
skip = skips[d]; c0 = (net.shape[-1] - skip.shape[-1]) // 2
net = torch.cat([skip, net[..., c0:c0 + skip.shape[-1]]], 1) # rule 3: centre on the skip
net = F.relu(self.up_bn1[i](self.up_conv1[i](net)))
return self.output_conv(net)
weights = gd.CACHE / "phasenet2018.pth" # the published weights, 1.1 MB
for attempt in range(4):
try:
if not weights.exists():
gd.CACHE.mkdir(exist_ok=True)
weights.write_bytes(gd.fetch("https://github.com/AI4EPS/models/releases/download/PhaseNet-2018/phasenet2018.pth").content)
break
except Exception as exc:
if attempt == 3:
raise RuntimeError(f"could not fetch the PhaseNet weights: {exc}; re-run this cell") from None
time.sleep(2 * (attempt + 1))
state = torch.load(weights, map_location="cpu", weights_only=False)["model"]
state = {k.removeprefix("backbone."): (v.squeeze(2) if v.ndim == 4 else v) for k, v in state.items() if k.startswith("backbone.")}
phasenet = PhaseNet2018()
phasenet.load_state_dict(state)
phasenet = phasenet.eval().to(DEVICE)
@torch.no_grad()
def phasenet_prob(batch):
"""P and S probabilities (n, 3, samples) for raw three-component records (n, 3, samples)."""
z = np.stack([standardise(x) for x in batch])
pad = -z.shape[-1] % 16
z = torch.from_numpy(np.pad(z, ((0, 0), (0, 0), (0, pad)))).to(DEVICE)
out = [torch.softmax(phasenet(z[k:k + 256]), 1)[..., :batch.shape[-1]].cpu().numpy() for k in range(0, len(z), 256)]
return np.concatenate(out)
def picks(prob_phase, vmin=0.5, half=0.5):
"""Local maxima of one phase's probability, at least `vmin`, the largest within `half` s."""
top = scipy.ndimage.maximum_filter1d(prob_phase, 2 * int(half * FS) + 1)
return np.where((top == prob_phase) & (prob_phase >= vmin))[0]
print(f"PhaseNet: {sum(p.numel() for p in phasenet.parameters()):,} parameters; "
f"a template: {3 * 100} numbers, copied from one earthquake")
PhaseNet: 268,443 parameters; a template: 300 numbers, copied from one earthquake
Now the experiment: 40 template earthquakes, the four kinds of planted earthquake, and seven signal-to-noise ratios from 0.25 to 16.
NOISE = 3700 # 37 s, all before the origin
live = np.where((np.abs(REC[:, :, :NOISE]).max(-1) > 0).all(1))[0]
usable = np.where((P[:, SQK] >= 100) & (P[:, SQK] + 800 <= X.shape[-1]))[0]
tpl_ids = rng.choice(np.intersect1d(with_s, usable), 40, replace=False)
sim_sqk = pair_cc(Xf, stations_=[SQK]) # A.3's measure, at SQK only
z_noise = bandpass(REC[:, 2, :NOISE]).std(-1)
def seismogram(kind, i):
"""A raw three-component seismogram from 1 s before P to 8 s after, and its similarity to i."""
if kind == "the template itself":
j = i
elif kind == "its closest look-alike":
j = usable[np.nanargmax(sim_sqk[i, usable])]
elif kind == "a typical neighbour":
j = rng.choice(usable[(usable != i) & (sep[i, usable] < 2)])
else:
return FIELD[rng.integers(len(FIELD))], np.nan
return X[j, SQK][:, P[j, SQK] - 100:P[j, SQK] + 800], (1.0 if j == i else sim_sqk[i, j])
kinds = ["the template itself", "its closest look-alike", "a typical neighbour", "an earthquake elsewhere"]
snrs = [0.25, 0.5, 1, 2, 4, 8, 16]
trials, data = [], []
for i in tpl_ids:
for kind in kinds:
for snr in snrs:
k = rng.choice(live[live != i])
x = REC[k, :, :NOISE].copy()
sig, sim = seismogram(kind, i)
sig = sig - sig.mean(-1, keepdims=True)
a = snr * z_noise[k] / np.abs(bandpass(sig)[2, 100:300]).max()
at = int(rng.integers(300, NOISE - 900))
x[:, at:at + 900] += a * sig
trials.append(dict(template=i, kind=kind, snr=snr, similarity=sim, p=at + 100))
data.append(x)
trials, data = pd.DataFrame(trials), np.stack(data)
t0 = time.perf_counter()
score_cc = [] # the peak near the planted P, in MAD units
for r, x in zip(trials.itertuples(), data):
c = scan(template_at(r.template), bandpass(x))
med = np.median(c); mad = np.median(np.abs(c - med))
k = r.p - 10 # the scan is indexed by P - 0.1 s
score_cc.append((c[k - 20:k + 21].max() - med) / mad)
time_cc = time.perf_counter() - t0
t0 = time.perf_counter()
prob = phasenet_prob(data.astype(np.float32))
time_pn = time.perf_counter() - t0
trials["score_cc"] = score_cc
trials["score_pn"] = [pr[1, r.p - 50:r.p + 51].max() for r, pr in zip(trials.itertuples(), prob)]
trials["template matching"] = trials.score_cc > 8 # the rules of A.4 and of Session 5
trials["PhaseNet"] = trials.score_pn >= 0.5
print(f"{len(trials):,} planted earthquakes; median similarity to the template at SQK: " +
", ".join(f"{k}: {trials[trials.kind == k].similarity.median():.2f}" for k in kinds[1:3]))
recall = trials.groupby(["kind", "snr"])[["template matching", "PhaseNet"]].mean()
1,120 planted earthquakes; median similarity to the template at SQK: its closest look-alike: 0.51, a typical neighbour: 0.22
Before the totals, one case to see what the numbers mean: the template's own earthquake planted at SNR 2, where template matching finds it and PhaseNet does not. The dashed line marks the planted P; the red dot, the correlation peak.
case = trials[(trials.kind == kinds[0]) & (trials.snr == 2) & trials["template matching"] & ~trials["PhaseNet"]].index[0]
r, x = trials.loc[case], data[case]
c = scan(template_at(r.template), bandpass(x))
med = np.median(c); mad = np.median(np.abs(c - med))
tt = np.arange(x.shape[-1]) / FS - r.p / FS # seconds from the planted P
fig, axs = plt.subplots(3, 1, figsize=(8, 5.2), sharex=True)
axs[0].plot(tt, bandpass(x)[2], c="0.3", lw=0.4)
axs[0].set(ylabel="velocity (um/s)", title="An earthquake planted at SNR 2 (vertical, 5-40 Hz): can you see it?")
axs[1].plot(tt[10:10 + len(c)], c, c="k", lw=0.5)
kp = r.p - 10 - 20 + int(np.argmax(c[r.p - 10 - 20:r.p - 10 + 21]))
axs[1].plot(tt[10 + kp], c[kp], "o", c="tab:red", ms=5)
axs[1].axhline(med + 8 * mad, c="tab:red", lw=1, label="8 x MAD")
axs[1].set(ylabel="correlation", title="Template matching: one clear peak above the threshold")
axs[2].plot(tt, prob[case, 1], c="tab:green", lw=0.8)
axs[2].axhline(0.5, c="tab:red", lw=1, label="0.5")
axs[2].set(ylabel="P probability", ylim=(0, 1), xlabel="seconds from the planted P",
title="PhaseNet: the probability never reaches 0.5")
for ax in axs:
ax.axvline(0, c="tab:blue", lw=0.8, ls="--")
ax.legend(fontsize=7, loc="upper right") if ax is not axs[0] else None
plt.tight_layout()
plt.show()
print(f"correlation at the planted P {c[r.p - 10 - 20:r.p - 10 + 21].max():.2f} against a threshold of {med + 8 * mad:.2f}; "
f"PhaseNet's largest P probability within 0.5 s: {prob[case, 1, r.p - 50:r.p + 51].max():.2f}")
correlation at the planted P 0.68 against a threshold of 0.46; PhaseNet's largest P probability within 0.5 s: 0.26
Now all 1,120 planted earthquakes. Each panel is one kind of planted earthquake; the two curves are the two methods. The horizontal gap between the curves at 50 % is the improvement: how many times weaker an earthquake template matching can find.
def snr_half(curve):
"""SNR at which a recall curve first reaches 50 %, by interpolation in log SNR (nan if never)."""
v = np.asarray(curve)
k = np.argmax(v >= 0.5)
if v[k] < 0.5:
return np.nan
if k == 0:
return snrs[0]
f = (0.5 - v[k - 1]) / (v[k] - v[k - 1])
return float(np.exp(np.log(snrs[k - 1]) + f * (np.log(snrs[k]) - np.log(snrs[k - 1]))))
def no_minor(ax):
ax.xaxis.set_minor_locator(plt.NullLocator())
def gap_arrow(ax, lo, hi):
"""An arrow at 50 % from the more sensitive curve to the other, labelled with the ratio."""
if np.isfinite(lo) and np.isfinite(hi) and hi / lo >= 1.5:
ax.annotate("", xy=(hi, 50), xytext=(lo, 50), arrowprops=dict(arrowstyle="<->", color="tab:red", lw=1.2))
ax.text(np.sqrt(lo * hi), 53, f"x{hi / lo:.0f}", color="tab:red", ha="center", fontsize=9)
fig, axs = plt.subplots(1, 4, figsize=(11, 3.2), sharey=True)
for ax, kind in zip(axs, kinds):
for method, col in (("template matching", "k"), ("PhaseNet", "tab:green")):
ax.semilogx(snrs, 100 * recall.loc[kind, method], "o-", c=col, label=method)
gap_arrow(ax, snr_half(recall.loc[kind, "template matching"]), snr_half(recall.loc[kind, "PhaseNet"]))
ax.set(title=kind, xlabel="signal-to-noise ratio", xticks=snrs, xticklabels=[f"{s:g}" for s in snrs])
no_minor(ax)
axs[0].set_ylabel("planted earthquakes found (%)")
fig.legend(*axs[0].get_legend_handles_labels(), loc="lower center", ncol=2, fontsize=8, bbox_to_anchor=(0.5, -0.08))
fig.suptitle("One station: template matching finds its own earthquake at a fraction of the signal, and misses the others", y=1.03)
plt.tight_layout()
plt.show()
half = pd.DataFrame({m: {k: snr_half(recall.loc[k, m]) for k in kinds} for m in ("template matching", "PhaseNet")})
print("signal-to-noise ratio at which half the planted earthquakes are found:")
print(half.round(1).to_string())
ratio = half.loc["the template itself", "template matching"]
print(f"PhaseNet needs {half.PhaseNet.median() / ratio:.0f} times the signal that template matching needs for the template itself")
signal-to-noise ratio at which half the planted earthquakes are found:
template matching PhaseNet
the template itself 1.2 6.3
its closest look-alike 4.9 7.0
a typical neighbour NaN 6.1
an earthquake elsewhere NaN 7.7
PhaseNet needs 5 times the signal that template matching needs for the template itself
The two panels are the comparison this session is about.
- For the template itself, template matching wins by a wide margin. It finds the earthquake at a signal-to-noise ratio several times smaller than PhaseNet needs. This is the matched filter theorem in action: for a known signal in white noise, correlating with the signal itself is the most sensitive linear detector there is.
- For its closest look-alike it still finds more at low SNR, but it levels off: however clear the signal, the correlation can never exceed the two waveforms' similarity (shown above), and only a pair similar enough clears 8 × MAD.
- For a typical neighbour or an earthquake elsewhere it finds almost nothing, whatever the SNR. The template is blind to earthquakes that do not resemble it.
- PhaseNet's curves lie almost on top of each other. What it needs is a clear enough arrival; it does not care which earthquake made it, because it learned what P waves look like from earthquakes everywhere in northern California.
So the two methods do not compete so much as divide the work. A general detector finds the earthquakes above its threshold; template matching then finds, below that threshold, copies of the earthquakes already found. Ross et al. (2019) took the routine catalogue's earthquakes as templates; Part B does the same for its hour, with the picks PhaseNet made in Session 5.
But are the two thresholds comparable? "8 × MAD" and "probability 0.5" are different rules, and a detector can always find more by lowering its threshold, at the price of more detections in noise. A fair comparison looks at both sides at once: for every threshold, the share of planted earthquakes found against the number of detections the same threshold makes on records with nothing planted. Tracing that out over all thresholds gives a curve, the receiver operating characteristic (ROC) of radar engineering; the better detector is the one whose curve lies higher and further left.
We need the second side: both methods on the same 37 s stretches with nothing planted. For template matching we keep the height of every correlation peak (in MAD units) for each of the 40 templates, and also for their reversed copies; for PhaseNet, every peak of its P probability.
quiet = REC[live, :, :NOISE]
quiet_f = bandpass(quiet)
hours_quiet = len(live) * NOISE / FS / 3600
peaks_tm, peaks_rev = [], []
for i in tpl_ids:
for out, tp in ((peaks_tm, template_at(i)), (peaks_rev, reverse(template_at(i)))):
for x in quiet_f:
c = scan(tp, x)
med = np.median(c); mad = np.median(np.abs(c - med))
z = (c - med) / mad
out.append(z[peaks(z, 3.0)]) # every peak above 3 MAD, 2 s apart
peaks_tm, peaks_rev = np.concatenate(peaks_tm), np.concatenate(peaks_rev)
prob_q = phasenet_prob(quiet.astype(np.float32))
peaks_pn = np.concatenate([pr[1][picks(pr[1], vmin=0.05)] for pr in prob_q])
def per_hour(heights, threshold, n_scans):
return (heights > threshold).sum() / n_scans / hours_quiet
k_mad, v_pn = np.arange(4, 16.5, 0.5), np.linspace(0.05, 0.95, 19)
fig, axs = plt.subplots(1, 2, figsize=(9, 3.4), sharey=True)
for ax, kind in zip(axs, (kinds[0], kinds[2])):
for snr, ls in ((2, "-"), (4, "--")):
t = trials[(trials.kind == kind) & (trials.snr == snr)]
ax.semilogx([max(per_hour(peaks_tm, k, len(tpl_ids)), 1e-2) for k in k_mad], [100 * (t.score_cc > k).mean() for k in k_mad],
ls, c="k", label=f"template matching, SNR {snr}")
ax.semilogx([max(per_hour(peaks_pn, v, 1), 1e-2) for v in v_pn], [100 * (t.score_pn >= v).mean() for v in v_pn],
ls, c="tab:green", label=f"PhaseNet, SNR {snr}")
ax.plot(max(per_hour(peaks_tm, 8, len(tpl_ids)), 1e-2), 100 * (t.score_cc > 8).mean(), "o", c="k", ms=4)
ax.plot(max(per_hour(peaks_pn, 0.5, 1), 1e-2), 100 * (t.score_pn >= 0.5).mean(), "o", c="tab:green", ms=4)
ax.set(title=kind, xlabel="detections per hour with nothing planted")
axs[0].set_ylabel("planted earthquakes found (%)")
axs[1].legend(fontsize=7, loc="center right")
fig.suptitle("At any threshold: template matching wins for its own earthquake, PhaseNet for the neighbour (dots: the rules used)", y=1.02)
plt.tight_layout()
plt.show()
print(f"{hours_quiet:.1f} h of records before catalogued earthquakes; detections per hour at the rules used:")
print(f" one template, 8 x MAD: {per_hour(peaks_tm, 8, len(tpl_ids)):.2f}; its reversed copy {per_hour(peaks_rev, 8, len(tpl_ids)):.2f}")
print(f" PhaseNet, P probability 0.5: {per_hour(peaks_pn, 0.5, 1):.1f}")
2.7 h of records before catalogued earthquakes; detections per hour at the rules used: one template, 8 x MAD: 0.12; its reversed copy 0.00 PhaseNet, P probability 0.5: 17.8
(The rate for template matching is that of one template; a bank of many templates detects in proportion more often. A rate of zero is drawn at 0.01.)
For the template's own earthquake, template matching's curve lies above PhaseNet's at every rate of detections in noise: there is no threshold at which PhaseNet catches up. For the typical neighbour it is the other way round: template matching finds some neighbours only once its threshold is so low that it detects tens to hundreds of times an hour, and at the same rate PhaseNet finds more. The improvement is therefore real, but it belongs to copies only — whatever thresholds are chosen.
Which finds more on these records? "Detections with nothing planted" are not all false alarms: these records are full of small earthquakes the catalogue does not list. Two tests, each with a null, tell us whether PhaseNet's picks are earthquakes. A real earthquake gives an S arrival a fraction of a second after P (the cell prints the S − P times of the ball at SQK); we count PhaseNet P picks followed by an S pick within that range, and the same count at random moments, which is the rate expected by chance. Then we compare with A.4's scan of all 138 templates on the same stretches: where a PhaseNet pick and a template detection fall within 0.5 s of each other, two methods that share nothing agree.
sp_ball = (S[with_s, SQK] - P[with_s, SQK]) / FS
lo, hi = np.quantile(sp_ball, [0.05, 0.95])
print(f"S - P at SQK for the ball: {lo:.2f}-{hi:.2f} s (5th-95th percentile)")
pn_p = [picks(pr[1]) for pr in prob_q]
pn_s = [picks(pr[2]) for pr in prob_q]
n_pn = sum(len(p) for p in pn_p)
s_after = lambda t, k: np.any((pn_s[k] - t >= lo * FS) & (pn_s[k] - t <= hi * FS))
with_s_after = sum(s_after(t, k) for k in range(len(live)) for t in pn_p[k])
chance = np.mean([s_after(t, k) for k in range(len(live)) for t in rng.integers(0, NOISE, 50)])
print(f"PhaseNet P picks: {n_pn} ({n_pn / hours_quiet:.1f} per hour); followed by an S pick at the ball's S - P: "
f"{with_s_after / n_pn:.0%}; at random moments: {chance:.1%}")
bank = det[(det.kind == "template") & (det["where"] == "before the catalogued earthquake")]
rec_of, at_of = np.divmod(bank["sample"].to_numpy() + 10, REC.shape[-1]) # minute, and P inside it
seen = np.isin(rec_of, live) & (at_of < NOISE) # inside what PhaseNet scanned
distinct = pd.DataFrame({"rec": rec_of[seen], "t": at_of[seen]}).sort_values(["rec", "t"])
distinct = distinct[(distinct.rec.diff() != 0) | (distinct.t.diff() > 100)] # one per event (1 s)
both = sum(np.any(np.abs(pn_p[np.searchsorted(live, r)] - a) <= 50) for r, a in zip(distinct.rec, distinct.t))
print(f"the 138-template bank (A.4) on the same stretches: {len(distinct)} distinct detections "
f"({len(distinct) / hours_quiet:.1f} per hour); also picked by PhaseNet: {both} ({both / len(distinct):.0%})")
S - P at SQK for the ball: 0.49-0.69 s (5th-95th percentile) PhaseNet P picks: 48 (17.8 per hour); followed by an S pick at the ball's S - P: 67%; at random moments: 0.1% the 138-template bank (A.4) on the same stretches: 26 distinct detections (9.6 per hour); also picked by PhaseNet: 11 (42%)
Most of PhaseNet's picks on these records are followed by an S pick at the right delay, far more often than chance, so most are probably small earthquakes missing from the catalogue rather than false alarms. The template bank finds fewer events in total, and part of what it finds PhaseNet picks as well; the rest are copies too weak for PhaseNet. So PhaseNet finds more events, and template matching adds events PhaseNet misses — the conclusion of the planted experiment, seen again on real background. "False alarm" is a statement about the referee as much as about the method; Part B makes the same comparison on a real hour, with the analysts' picks and the array as referees.
Round two: ten stations¶
Everything so far used one station. A real scan uses the whole network, and so does a real picker, through association: a pick on one station alone could be noise, so a catalogue declares an earthquake only when several stations pick arrivals that fit one source. A.5 showed why the network helps template matching; does it help both methods equally?
For 60 of the earthquakes the file holds the records of all 10 stations on a common clock: 37 s of noise, then the earthquake. We plant earthquakes into the 10-station noise as before, now with every station's seismogram of the planted earthquake, all scaled by one factor so that the stations keep their true relative amplitudes; the SNR is still the one at SQK, the second-closest station. Template matching now averages the correlation traces of every station and component, each shifted by its template's own arrival time — the moveout. For PhaseNet we count on how many stations it picks the planted P, and declare a detection when at least 3 do, a typical minimum for association. The fourth kind, an earthquake elsewhere, is left out because we have it on SQK only.
NET = pack["net"].astype(np.float64) * pack["net_scale"][..., None] # (60, 10 stations, E N Z, 4700)
NETf = bandpass(NET)
NOISE_NET = 3700 # the window starts at 37 s
def template_net(i):
"""(station, piece, component, offset): P on Z and S on E and N wherever the analyst picked them,
offsets counted from the window start (1 s before the origin), which is the moveout."""
out = []
for s in range(NS):
if P[i, s] >= 10:
out.append((s, Xf[i, s, 2, P[i, s] - 10:P[i, s] + 90], 2, P[i, s] - 10))
if S[i, s] >= 10 and S[i, s] + 90 <= X.shape[-1]:
out += [(s, Xf[i, s, c, S[i, s] - 10:S[i, s] + 90], c, S[i, s] - 10) for c in (0, 1)]
return out
def scan_net(template, xf, only=None):
"""Average over every station and component of the correlation, shifted by the moveout, so the
output is indexed by the window start of a matching earthquake. `only`: a list of stations."""
n = xf.shape[-1] - X.shape[-1]
return np.mean([ncc(t, xf[s, c])[off:off + n] for s, t, c, off in template if only is None or s in only], axis=0)
all_p = np.where((P >= 10).sum(1) >= 8)[0] # earthquakes picked on 8 or more stations
sim_net = pair_cc(Xf) # A.3's 10-station similarity
tpl_net = rng.choice(np.intersect1d(all_p, with_s), 25, replace=False)
noise_sqk = NETf[:, SQK, 2, :NOISE_NET].std(-1)
kinds_net, snrs_net = kinds[:3], [0.25, 0.5, 1, 2, 4, 8, 16]
rows, data = [], []
for i in tpl_net:
for kind in kinds_net:
j = (i if kind == kinds[0] else all_p[np.nanargmax(sim_net[i, all_p])] if kind == kinds[1]
else rng.choice(all_p[(all_p != i) & (sep[i, all_p] < 2)]))
for snr in snrs_net:
k = rng.integers(len(NET))
x = NET[k, :, :, :NOISE_NET].copy()
sig = X[j] - X[j].mean(-1, keepdims=True)
a = snr * noise_sqk[k] / np.abs(bandpass(sig[SQK])[2, P[j, SQK]:P[j, SQK] + 200]).max()
at = int(rng.integers(200, NOISE_NET - X.shape[-1]))
x[:, :, at:at + X.shape[-1]] += a * sig
rows.append(dict(template=i, planted=j, kind=kind, snr=snr, at=at))
data.append(x)
net_trials, data = pd.DataFrame(rows), np.stack(data)
xf_all = bandpass(data)
for name, only in (("template matching, SQK", [SQK]), ("template matching, 10 stations", None)):
hit = []
for r, xf in zip(net_trials.itertuples(), xf_all):
c = scan_net(template_net(r.template), xf, only)
med = np.median(c); mad = np.median(np.abs(c - med))
hit.append(c[r.at - 20:r.at + 21].max() > med + 8 * mad)
net_trials[name] = hit
prob = phasenet_prob(data.reshape(-1, 3, NOISE_NET).astype(np.float32)).reshape(len(net_trials), NS, 3, NOISE_NET)
net_trials["stations picking P"] = [sum(np.any(np.abs(picks(pr[s, 1]) - (r.at + P[r.planted, s])) <= 50)
for s in range(NS) if P[r.planted, s] >= 0)
for r, pr in zip(net_trials.itertuples(), prob)]
net_trials["PhaseNet, 3 or more stations"] = net_trials["stations picking P"] >= 3
cols = ["template matching, SQK", "template matching, 10 stations", "PhaseNet, 3 or more stations"]
recall_net = net_trials.groupby(["kind", "snr"])[cols].mean()
half_net = pd.DataFrame({c: {k: snr_half(recall_net.loc[k, c]) for k in kinds_net} for c in cols})
fig, axs = plt.subplots(1, 3, figsize=(10, 3.2), sharey=True)
for ax, kind in zip(axs, kinds_net):
for col, c in zip(cols, ("0.65", "k", "tab:green")):
ax.semilogx(snrs_net, 100 * recall_net.loc[kind, col], "o-", c=c, label=col)
gap_arrow(ax, half_net.loc[kind, cols[1]], half_net.loc[kind, cols[2]])
ax.set(title=kind, xlabel="signal-to-noise ratio at SQK", xticks=snrs_net, xticklabels=[f"{s:g}" for s in snrs_net])
no_minor(ax)
axs[0].set_ylabel("planted earthquakes found (%)")
fig.legend(*axs[0].get_legend_handles_labels(), loc="lower center", ncol=3, fontsize=8, bbox_to_anchor=(0.5, -0.08))
fig.suptitle("Ten stations: the gain grows for copies; the typical neighbour is still missed", y=1.03)
plt.tight_layout()
plt.show()
print("signal-to-noise ratio at SQK at which half are found:")
print(half_net.round(2).to_string())
peak = np.array([[np.abs(bandpass(X[j, s, 2])[max(P[j, s], 0):max(P[j, s], 0) + 200]).max() for s in range(NS)]
for j in net_trials.planted.unique()])
snr_rel = (peak / NETf[:, :, 2, :NOISE_NET].std(-1).mean(0)) / (peak[:, [SQK]] / NETf[:, SQK, 2, :NOISE_NET].std(-1).mean())
print("each station's SNR relative to SQK's (median over the planted earthquakes): " +
", ".join(f"{st} {v:.1f}" for st, v in zip(stations.station, np.median(snr_rel, 0))))
signal-to-noise ratio at SQK at which half are found:
template matching, SQK template matching, 10 stations PhaseNet, 3 or more stations
the template itself 1.29 0.48 5.96
its closest look-alike 5.66 2.59 5.19
a typical neighbour NaN NaN 6.29
each station's SNR relative to SQK's (median over the planted earthquakes): SB4 0.2, SQK 1.0, GDXB 0.8, DRK 1.0, AL3 0.9, HVC 0.6, MCL 0.5, FNF 0.6, HER 0.0, RTB 0.4
Three results, read from the curves and the table under them.
- Ten stations make template matching more sensitive to its own earthquake: every station adds a correlation trace that peaks at the same moment, while their noise averages down (A.5).
- They help with the closest look-alike too, but it still levels off well below 100 %, and they do not help at all with the typical neighbour. The stack lowers the noise, but it cannot raise the correlation above the two earthquakes' similarity: averaging ten stations of a poor match is still a poor match.
- PhaseNet needs about the same signal as on one station. Association asks each station to pick on its own, and several stations see the planted earthquake about as clearly as SQK does (the last line printed above), so three of them pick it at about the SNR at which SQK alone would. Association does not add sensitivity; it removes single-station false picks.
So at network scale the division of labour is sharper still: template matching reaches deeper for copies, and only for copies. (Each point of these curves rests on 25 planted earthquakes, and of the single-station curves on 40, so differences of about 10 % between neighbouring points are within chance.)
The false alarms on the 10-station noise, again with nothing planted:
quiet_net = NETf[:, :, :, :NOISE_NET]
fa = {"templates": 0, "reversed templates": 0}
for i in tpl_net:
tp = template_net(i)
for name, tt in (("templates", tp), ("reversed templates", [(s, t[::-1].copy(), c, o) for s, t, c, o in tp])):
for xf in quiet_net:
c = scan_net(tt, xf)
med = np.median(c); mad = np.median(np.abs(c - med))
fa[name] += len(peaks(c, med + 8 * mad))
hours_net = len(NET) * NOISE_NET / FS / 3600
pq = phasenet_prob(NET[:, :, :, :NOISE_NET].reshape(-1, 3, NOISE_NET).astype(np.float32)).reshape(len(NET), NS, 3, NOISE_NET)
def associate(pick_lists):
"""Events: 3 or more P picks (any stations) within 3 s."""
t_p, n, k = np.sort(np.concatenate(pick_lists)), 0, 0
while k < len(t_p):
group = t_p[(t_p >= t_p[k]) & (t_p < t_p[k] + 3 * FS)]
n, k = (n + 1, k + len(group)) if len(group) >= 3 else (n, k + 1)
return n
net_picks = [[picks(pq[e, s, 1]) for s in range(NS)] for e in range(len(NET))]
associated = sum(associate(p) for p in net_picks)
# chance: shift each station's picks by its own random amount (wrapping round), which keeps every
# station's rate but destroys any common moment
by_chance = np.mean([sum(associate([(p + rng.integers(NOISE_NET)) % NOISE_NET for p in ps]) for ps in net_picks)
for _ in range(20)])
print(f"{hours_net:.1f} h of 10-station records before catalogued earthquakes:")
print(f" template matching: {fa['templates'] / len(tpl_net) / hours_net:.2f} detections per template per hour; "
f"reversed templates {fa['reversed templates'] / len(tpl_net) / hours_net:.2f}")
print(f" PhaseNet: {associated / hours_net:.1f} associated events per hour (3 or more stations within 3 s); "
f"with each station's picks shifted at random: {by_chance / hours_net:.1f}")
0.6 h of 10-station records before catalogued earthquakes: template matching: 0.13 detections per template per hour; reversed templates 0.00 PhaseNet: 25.9 associated events per hour (3 or more stations within 3 s); with each station's picks shifted at random: 11.9
PhaseNet associates many events per hour in records that the catalogue calls quiet. About half that rate also arises by chance — the same picks with each station shifted at random still line up three at a time — so roughly half are coincidences of single-station picks, and the excess is probably small earthquakes the catalogue missed. A real associator also checks that the pick times fit one source location, which removes most chance groupings. Each template, by contrast, finds only its few copies, and its reversed version almost nothing. The two methods see different parts of the same seismicity.
What each method needs and gives¶
The table gathers what the two methods need and give, with costs measured on this computer.
per_hour_cc = time_cc / len(trials) / (NOISE / FS / 3600)
per_hour_pn = time_pn / len(trials) / (NOISE / FS / 3600)
print(f"one hour of one three-component station: template matching {1000 * per_hour_cc:.0f} ms per template "
f"(numpy on the processor, band-pass included), so {per_hour_cc * len(with_s):.0f} s for the 138-template bank of A.4; "
f"PhaseNet {1000 * per_hour_pn:.0f} ms on {DEVICE}, whatever the number of earthquakes")
one hour of one three-component station: template matching 95 ms per template (numpy on the processor, band-pass included), so 13 s for the 138-template bank of A.4; PhaseNet 73 ms on mps, whatever the number of earthquakes
| template matching | a trained picker (PhaseNet) | |
|---|---|---|
| needs | one recorded earthquake per template, and its picks | many labelled examples, once, to train |
| finds | copies of its templates, at a much lower SNR | any earthquake, if the arrival is clear enough |
| blind to | earthquakes unlike every template | arrivals too weak for its threshold |
| timing | relative to the template, to about a millisecond (A.6) | absolute picks, as good as the labels (Session 5) |
| cost | grows with the number of templates | fixed per hour of data |
| more stations | more sensitive: correlations add, noise averages down (A.5) | fewer false picks through association, about the same sensitivity |
| false alarms | measured with a null (A.4); coda of nearby earthquakes is not caught by it | measured against analysts and arrays (Part B) |
A.8 · Between the two: from one template to a learned basis¶
A template is blind to earthquakes that do not resemble it. Could we widen what it accepts without training a network? One classical answer is a subspace detector. Instead of one template, use what many earthquakes have in common.
Stack the template vectors of many earthquakes as the rows of a matrix — each row the P window on Z followed by the S windows on E and N, demeaned and of unit length. Its singular value decomposition (Sessions 1–3) gives orthonormal directions \(\mathbf{u}_1, \mathbf{u}_2, \dots\), ordered by how much of the rows' energy each explains. Keep the first \(d\) of them as the columns of a matrix \(U\). A stretch of record \(\mathbf{x}\) (demeaned) is then scored by the fraction of its length that lies in their span,
which for \(d = 1\) is just the correlation with the first direction — the shape most earthquakes of the ball share. Larger \(d\) admits more varied earthquakes. This is a learned representation: the directions are fitted to data, and \(U^\top\mathbf{x}\) is a \(d\)-number summary of a seismogram, an embedding. It is learned without labels, and it is linear.
We test it on the hardest case for template matching, the typical neighbour. For each planted earthquake we build \(U\) from all the other earthquakes of the ball, so the planted one is never in its own basis.
sp = int(np.median(S[with_s, SQK] - P[with_s, SQK])) # one S - P for the whole ball (samples)
def vec_at(x, p):
"""P on Z, then S on E and N at P + sp, 1 s each from 0.1 s before."""
return np.r_[x[2, p - 10:p + 90], x[0, p + sp - 10:p + sp + 90], x[1, p + sp - 10:p + sp + 90]]
basis_ids = usable[P[usable, SQK] + sp + 90 <= X.shape[-1]]
V = unit(np.stack([vec_at(Xf[i, SQK], P[i, SQK]) for i in basis_ids]))
def subspace_scan(U, xf):
"""||U^T x|| / ||x|| for every stretch of a filtered record, indexed by P - 0.1 s."""
n = xf.shape[-1] - sp - 100
sw = np.lib.stride_tricks.sliding_window_view
W = np.concatenate([sw(xf[2], 100)[:n], sw(xf[0, sp:], 100)[:n], sw(xf[1, sp:], 100)[:n]], axis=1)
W = W - W.mean(1, keepdims=True)
nw = np.linalg.norm(W, axis=1)
return np.linalg.norm(W @ U, axis=1) / np.where(nw > 0, nw, 1)
rows = []
for j in rng.choice(basis_ids, 40, replace=False):
_, sv, Vh = np.linalg.svd(V[basis_ids != j], full_matrices=False)
for d in (1, 3, 10):
for snr in (1, 2, 4, 8, 16):
k = rng.choice(live[live != j])
x = REC[k, :, :NOISE].copy()
sig = X[j, SQK][:, P[j, SQK] - 100:P[j, SQK] + 800]
sig = sig - sig.mean(-1, keepdims=True)
x[:, (at := int(rng.integers(300, NOISE - 900))):at + 900] += snr * z_noise[k] / np.abs(bandpass(sig)[2, 100:300]).max() * sig
c = subspace_scan(Vh[:d].T, bandpass(x))
med = np.median(c); mad = np.median(np.abs(c - med))
q = at + 100 - 10
rows.append(dict(d=d, snr=snr, hit=c[q - 20:q + 21].max() > med + 8 * mad))
energy = np.cumsum(sv ** 2) / np.sum(sv ** 2)
_, _, Vh_all = np.linalg.svd(V, full_matrices=False)
alarms = {}
for d in (1, 3, 10): # detections with nothing planted
n = 0
for x in quiet_f:
c = subspace_scan(Vh_all[:d].T, x)
med = np.median(c); mad = np.median(np.abs(c - med))
n += len(peaks(c, med + 8 * mad))
alarms[d] = n / hours_quiet
sub = pd.DataFrame(rows).groupby(["d", "snr"]).hit.mean().unstack()
print("share of the earthquakes' energy in the first 1, 3, 10 directions: " + ", ".join(f"{energy[d - 1]:.0%}" for d in (1, 3, 10)))
print("detections per hour with nothing planted: " + ", ".join(f"d = {d}: {v:.1f}" for d, v in alarms.items()) +
f"; one template: {per_hour(peaks_tm, 8, len(tpl_ids)):.2f}; PhaseNet: {per_hour(peaks_pn, 0.5, 1):.1f}")
print("typical neighbours found (%), by dimension d (rows) and SNR (columns):")
print((100 * sub).round(0).astype(int).to_string())
print("for comparison, at the same SNRs: one template",
(100 * recall.loc["a typical neighbour", "template matching"].loc[[1, 2, 4, 8, 16]]).round(0).astype(int).tolist(),
"and PhaseNet", (100 * recall.loc["a typical neighbour", "PhaseNet"].loc[[1, 2, 4, 8, 16]]).round(0).astype(int).tolist())
share of the earthquakes' energy in the first 1, 3, 10 directions: 8%, 21%, 49% detections per hour with nothing planted: d = 1: 24.8, d = 3: 9.6, d = 10: 3.3; one template: 0.12; PhaseNet: 17.8 typical neighbours found (%), by dimension d (rows) and SNR (columns): snr 1 2 4 8 16 d 1 12 28 38 42 42 3 8 28 28 55 48 10 2 10 30 48 52 for comparison, at the same SNRs: one template [0, 0, 0, 0, 0] and PhaseNet [0, 0, 10, 75, 98]
The subspace finds typical neighbours that a single template misses entirely, and the table shows its price: with nothing planted, one direction detects far more often than one template — more often than PhaseNet picks — while ten directions detect much less often for a similar share of neighbours. Either way the subspace finds fewer neighbours than PhaseNet at SNR 8 and above, and nothing much below SNR 2. The first few directions explain only a small share of the energy, because the waveforms of earthquakes even 1 km apart differ in detail; a linear basis cannot capture what they have in common without also admitting noise.
That is the gap a neural network fills. Its layers build nonlinear features: a convolution, a ReLU, another convolution, so that "an impulsive arrival followed by a slower, larger one on the horizontals" becomes something the network can respond to whatever the exact wiggles. Trained on labels it becomes a picker (Session 5); trained without labels, for example as an autoencoder that compresses a seismogram and rebuilds it, its compressed layer is a learned embedding in which similar seismograms lie close together, and clustering in it separates, for instance, waveforms recorded at different distances (Mousavi et al. 2019). Correlation is the special case of one fixed linear feature.
A.9 · The programs the field runs¶
- EQcorrscan (Chamberlain et al. 2018), matched-filter and subspace detection in Python, built on ObsPy: the program a student would use to run template matching for real.
- FMF, the Fast Matched Filter (Beaucé et al. 2018): the correlation of A.2 written for processors and graphics cards, for scans of thousands of templates.
- QTM, the Quake Template Matching scan of Southern California (Ross et al. 2019), which added about ten times more earthquakes to the catalogue: the state of the art.
- CCTorch, our group's code: correlation, MAD thresholds and subsample lags on graphics cards in
PyTorch.
nccin A.2 is distilled from it; Part B uses the same conventions. - The northern California correlation catalogue (Schaff & Waldhauser 2005), correlation differential times for two decades of earthquakes, measured as in A.6 and used for HypoDD relocations; Schaff & Waldhauser (2010) found that correlation detection lowers the detection threshold by about one magnitude unit at Parkfield.
Takeaways¶
Machine learning
- Normalised cross-correlation is the cosine of the angle between two demeaned vectors, so it measures shape and ignores amplitude.
- Sliding a template along a record is a convolution — a one-kernel convolutional network whose kernel is a copy of one earthquake — and the FFT and cumulative sums make it cheap.
- A detection threshold means nothing until its false-alarm rate is measured with a null; noise is heavier-tailed than a Gaussian, and the median absolute deviation is a robust scale for it.
- Averaging correlations over \(N\) channels shrinks the noise roughly by \(\sqrt{N}\) while a real copy stays, which is why network stacking works.
- The matched filter is the most sensitive detector for a known signal and blind to unknown ones; a trained network is general but needs a clearer signal. They complement each other.
- Two detectors are compared fairly only at equal rates of detection in noise, by tracing each one's share found against that rate over all its thresholds (the ROC curve), not at one threshold each.
- A subspace from the singular value decomposition is a linear learned embedding; a network's features are its nonlinear generalisation.
Seismology
- Two earthquakes have similar seismograms when they are close compared with a wavelength and share a mechanism, so lower frequencies stay similar over larger distances.
- Template matching finds only earthquakes near ones already known, which suits the tightly clustered seismicity of the Geysers.
- Correlation lags time arrivals to about a millisecond, far better than picks, which is why double-difference relocation uses them.
- Many detections that look like noise are small uncatalogued earthquakes; two independent methods agreeing is the evidence, and the coda of nearby earthquakes is a false alarm the usual null does not catch.
References¶
- Beaucé, E., Frank, W. B., & Romanenko, A. (2018). Fast Matched Filter (FMF): an efficient seismic matched-filter search for both CPU and GPU architectures. Seismol. Res. Lett., 89, 165–172. https://doi.org/10.1785/0220170181
- Chamberlain, C. J., Hopp, C. J., Boese, C. M., Warren-Smith, E., Chambers, D., Chu, S. X., Michailos, K., & Townend, J. (2018). EQcorrscan: repeating and near-repeating earthquake detection and analysis in Python. Seismol. Res. Lett., 89, 173–181. https://doi.org/10.1785/0220170151
- Geller, R. J., & Mueller, C. S. (1980). Four similar earthquakes in central California. Geophys. Res. Lett., 7, 821–824. https://doi.org/10.1029/GL007i010p00821
- Gibbons, S. J., & Ringdal, F. (2006). The detection of low magnitude seismic events using array-based waveform correlation. Geophys. J. Int., 165, 149–166. https://doi.org/10.1111/j.1365-246X.2006.02865.x
- Mousavi, S. M., Zhu, W., Ellsworth, W., & Beroza, G. (2019). Unsupervised clustering of seismic signals using deep convolutional autoencoders. IEEE Geosci. Remote Sens. Lett., 16, 1693–1697. https://doi.org/10.1109/LGRS.2019.2909218
- Poupinet, G., Ellsworth, W. L., & Frechet, J. (1984). Monitoring velocity variations in the crust using earthquake doublets: an application to the Calaveras Fault, California. J. Geophys. Res., 89, 5719–5731. https://doi.org/10.1029/JB089iB07p05719
- Ross, Z. E., Trugman, D. T., Hauksson, E., & Shearer, P. M. (2019). Searching for hidden earthquakes in Southern California. Science, 364, 767–771. https://doi.org/10.1126/science.aaw6888
- Schaff, D. P., & Waldhauser, F. (2005). Waveform cross-correlation-based differential travel-time measurements at the Northern California Seismic Network. Bull. Seismol. Soc. Am., 95, 2446–2461. https://doi.org/10.1785/0120040221
- Schaff, D. P., & Waldhauser, F. (2010). One magnitude unit reduction in detection threshold by cross correlation applied to Parkfield (California) and China seismicity. Bull. Seismol. Soc. Am., 100, 3224–3238. https://doi.org/10.1785/0120100042
- 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