Focal mechanism & moment tensor · Part B: the 2016 Mw 5.0, real and twin¶
For the same Geysers earthquakes, network moment tensors and first-motion mechanisms disagree — how much of that is the Earth, and how much is the inversion?
Part A built the methods — FPFIT, HASH, TDMT and CAP — on planted mechanisms. This notebook takes them to one earthquake, worked in depth: the Mw 5.0 of 14 December 2016, the largest in the catalogue window of the Geysers field, and the one earthquake here with both a first-motion mechanism and waveform moment tensors. They disagree by some sixty degrees, which makes it the sharpest form of the session's question.
Each data type is worked in three steps.
- The real records — what the methods return, and what question that raises.
- A twin — an earthquake we invent, recorded by the same stations along the same rays, with noise measured from the real data where it can be. Its truth is known, so we can measure what each method recovers, and change one setting at a time to see what matters.
- The real records again, read in the light of the twin: what the real answer means.
You can run it without running Part A first.
The questions this notebook answers¶
- What do the three published mechanisms of this earthquake say, and how are they labelled? (B.1)
- How good are its first motions, and what can a first-motion mechanism from readings like these be trusted to say? (B.2–B.4)
- What does a waveform inversion of its broadband records return, how does it depend on choices the data cannot make, and is its non-double-couple part real? (B.5–B.8)
- So how much of the sixty degrees is the Earth, and how much the inversion? (B.9)
import importlib.util
import sys
import time
import urllib.request
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# The course's modules sit beside the notebook on DataHub; Colab opens the notebook alone, so there
# they are fetched. `geysers_mech` holds FPFIT, HASH and CAP in full; Part A wrote out the core of each.
BASE = "https://raw.githubusercontent.com/AI4EPS/EPS207_Observational_Seismology/main/docs/notebooks"
for module in ("geysers_data.py", "geysers_func.py", "geysers_loc.py", "geysers_mech.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))
from geysers_func import sdr_to_mt, mt_to_vec, vec_to_mt, nodal_planes, kagan, mt_decompose
import geysers_mech as gm
rng = np.random.default_rng(207)
plt.rcParams.update({"figure.dpi": 120, "font.size": 9, "axes.grid": True,
"grid.alpha": 0.3, "axes.axisbelow": True})
def sdr(M, plane=0):
"""One nodal plane of a tensor as three rounded numbers, for printing."""
return tuple(round(float(v)) for v in nodal_planes(M)[plane])
def ray_rows(azimuth, takeoff):
"""One row of the design matrix per ray, so that the P amplitude is `row @ m` (Part A, A.1)."""
a, i = np.radians(azimuth), np.radians(takeoff)
g = np.stack([np.sin(i) * np.cos(a), np.sin(i) * np.sin(a), np.cos(i)], axis=-1)
return np.stack([g[..., 0] ** 2, g[..., 1] ** 2, g[..., 2] ** 2,
2 * g[..., 0] * g[..., 1], 2 * g[..., 0] * g[..., 2],
2 * g[..., 1] * g[..., 2]], axis=-1)
def dc_grid(step_strike, step_dip, step_rake):
"""Every double couple on a strike/dip/rake grid, each of unit size, so that |g @ m| <= 1."""
s, d, r = [a.ravel() for a in np.meshgrid(np.arange(0, 360, step_strike), np.arange(step_dip, 90.1, step_dip),
np.arange(-180, 180, step_rake), indexing="ij")]
V = np.array([mt_to_vec(sdr_to_mt(*x)) for x in zip(s, d, r)])
return np.stack([s, d, r], 1), V
def fpfit_simple(azimuth, takeoff, y, quality, V=None):
"""FPFIT's misfit on one 5/5/10-degree grid, with its tie-break (Part A, A.10)."""
V = V5 if V is None else V
g = ray_rows(azimuth, takeoff)
A = g @ V.T
w = gm.obs_weight(quality)[:, None] * np.sqrt(np.abs(A))
F = (w * (np.sign(A) != y[:, None])).sum(0) / w.sum(0)
tie = np.flatnonzero(F <= F.min() + 1e-9)
return tie[np.argmax(w.sum(0)[tie])]
def angle_trials(azimuth, takeoff, n=30, sigma=5.0):
"""n versions of the ray angles: the given ones, and n-1 perturbed by sigma degrees (HASH, A.11)."""
return (np.column_stack([azimuth] + [azimuth + rng.normal(0, sigma, len(azimuth)) for _ in range(n - 1)]),
np.column_stack([takeoff] + [takeoff + rng.normal(0, sigma, len(takeoff)) for _ in range(n - 1)]))
SDR5, V5 = dc_grid(5.0, 5.0, 10.0)
HASH_GRID = gm.hash_grid(5.0)
print("seed 207: every number below is reproducible")
seed 207: every number below is reproducible
The earthquake¶
B.1 · Three answers for one earthquake¶
On 14 December 2016 an Mw 5.0 struck the middle of the producing field — the largest event in the field window in the catalogue. Its USGS event page is nc72737985, and it carries three mechanisms from three methods: a first-motion mechanism from the Northern California network (FPFIT on the polarities), and two waveform moment tensors, one from Berkeley and one from the USGS, each fitted to long-period records independently of the other.
One detail on that page matters before any number does: each product has a status. Both waveform tensors are marked reviewed; the first-motion mechanism is marked preliminary. A preliminary product is the automatic solution, never checked by an analyst.
NC = 1e7 * np.array([[-2.496e16, -2.492e16, 3.796e15], # Berkeley's tensor, converted from
[-2.492e16, 3.829e16, 7.257e15], # (up, south, east) newton-metres to
[3.796e15, 7.257e15, -1.333e16]]) # (north, east, down) dyne-centimetres
US = 1e7 * np.array([[-2.84e16, -4.192e16, -7.88e15], # the USGS tensor, likewise
[-4.192e16, 5.021e16, 2.291e16],
[-7.88e15, 2.291e16, -2.18e16]])
FIRST_MOTION = sdr_to_mt(140.0, 75.0, -140.0) # the network's first-motion mechanism
for name, M in (("Berkeley, waveforms (reviewed)", NC), ("USGS, waveforms (reviewed)", US)):
d = mt_decompose(M)
print(f"{name:32s} Mw {d['Mw']:.2f} DC {d['dc']:3.0f} % CLVD {d['clvd']:3.0f} % plane {sdr(M)}")
print(f"{'network, first motions (prelim.)':32s} {'':34s} plane {sdr(FIRST_MOTION)}")
print(f"\nthe two waveform tensors differ from each other by {kagan(NC, US):.1f}°")
print(f"the first-motion mechanism differs from Berkeley's by {kagan(FIRST_MOTION, NC):.1f}° "
f"and from the USGS's by {kagan(FIRST_MOTION, US):.1f}°")
Berkeley, waveforms (reviewed) Mw 5.08 DC 49 % CLVD 51 % plane (335, 83, -165) USGS, waveforms (reviewed) Mw 5.21 DC 26 % CLVD 74 % plane (158, 84, 166) network, first motions (prelim.) plane (140, 75, -140) the two waveform tensors differ from each other by 13.5° the first-motion mechanism differs from Berkeley's by 61.6° and from the USGS's by 60.3°
Two notes on what the page reports. The magnitudes printed here are computed from each tensor's scalar moment with one convention; the agencies report Mw 5.01 (Berkeley) and Mww 5.1 (USGS) by their own, so compare the mechanisms, not these magnitudes. And the agencies computed their tensors for sources at 5 km (Berkeley) and 11.5 km (USGS), well below the 1.5 km at which the first arrivals place the hypocentre. (The page also carries a second, regional USGS tensor, Mwr 5.0; we use the two above.)
The two reviewed waveform tensors agree with each other about as well as two agencies ever do. The preliminary first-motion mechanism is some sixty degrees from both. Before asking what that means about the earthquake, look at the first motions themselves, with the tools of Part A.
First motions¶
B.2 · The real readings¶
The readings come from the course's polarity asset, with three stations left out — their readings disagree with the network's own published mechanisms most of the time, so they cannot be trusted — and the borehole channels' polarity convention already corrected. Run Part A's methods on them: FPFIT, as the network does, and HASH, with a 5° uncertainty on every take-off angle and azimuth.
from geysers_data import polarities
lab = polarities()
lab = lab[lab.polarity.isin(["U", "D"]) & ~lab.station.isin(["PSB", "SRB", "SSB"])].copy()
lab["y"] = np.where(lab.polarity.values == "U", 1, -1)
g5 = lab[lab.event_id == "nc72737985"].copy()
q5 = np.where(g5.remark.values == "i", 0, 2) # FPFIT's classes, as in A.10
f5 = gm.fpfit(g5.azimuth.values, g5.takeoff_angle.values, g5.y.values, q5)[0]
OURS_FM = sdr_to_mt(f5["strike"], f5["dip"], f5["rake"])
az5_mc = np.column_stack([g5.azimuth.values] + [g5.azimuth.values + rng.normal(0, 5, len(g5)) for _ in range(29)])
ih5_mc = np.column_stack([g5.takeoff_angle.values] + [g5.takeoff_angle.values + rng.normal(0, 5, len(g5)) for _ in range(29)])
h5 = gm.hash_mech(az5_mc, 180 - ih5_mc, g5.y.values.astype(float), (q5 > 0).astype(int), _grid=HASH_GRID)[0][0]
G5 = ray_rows(g5.azimuth.values, g5.takeoff_angle.values)
explained = lambda M: (np.sign(G5 @ mt_to_vec(M)) == g5.y.values).mean()
print(f"{len(g5)} usable first motions, {(g5.remark == 'i').sum()} impulsive, at "
f"{g5.distance_km.min():.1f}-{g5.distance_km.max():.1f} km")
print(f"our FPFIT returns {sdr(OURS_FM)}, {kagan(OURS_FM, FIRST_MOTION):.0f}° from the network's preliminary one")
print(f"HASH grades this event {h5['quality']}: probability {h5['prob']:.2f} (the share of its acceptable set within 45° "
f"of the answer), rms angle of the set about the answer {np.mean(h5['rms']):.0f}°")
print(f"\nshare of these first motions explained by: our FPFIT {explained(OURS_FM):.2f}, the network's "
f"{explained(FIRST_MOTION):.2f}, Berkeley's tensor {explained(NC):.2f}, the USGS's {explained(US):.2f}")
37 usable first motions, 21 impulsive, at 1.5-28.8 km our FPFIT returns (138, 78, -112), 29° from the network's preliminary one HASH grades this event D: probability 0.40 (the share of its acceptable set within 45° of the answer), rms angle of the set about the answer 40° share of these first motions explained by: our FPFIT 0.78, the network's 0.76, Berkeley's tensor 0.68, the USGS's 0.65
Our FPFIT, fed these readings, lands about thirty degrees from the network's preliminary answer. The polarity asset keeps only the stations within 30 km, though, and the network's FPFIT saw more: the course's phase asset keeps every polarity the network read for this earthquake, out to several hundred kilometres. Give our FPFIT exactly those, with the network's own pick classes and the borehole correction.
from geysers_data import phases
ph5 = phases(2016)
ph5 = ph5[(ph5.event_id == 72737985) & (ph5.phase_type == "P") & ph5.phase_polarity.isin(["U", "D"])].copy()
ph5["y"] = np.where(ph5.phase_polarity == "U", 1, -1)
ph5.loc[(ph5.network == "BG") & (ph5.instrument == "DP") & (ph5.component == "Z"), "y"] *= -1 # the borehole flip
ph5["q"] = np.rint(4 * (1 - ph5.phase_score)).astype(int) # pick class 0..4
f_net = gm.fpfit(ph5.azimuth.values, ph5.takeoff_angle.values, ph5.y.values, ph5.q.values,
(ph5.review_status == "automatic").values)[0]
M_net = sdr_to_mt(f_net["strike"], f_net["dip"], f_net["rake"])
print(f"{len(ph5)} polarities to {ph5.distance_km.max():.0f} km; pick classes {ph5.q.value_counts().sort_index().to_dict()}")
print(f"our FPFIT on them returns {sdr(M_net)}: {kagan(M_net, FIRST_MOTION):.1f}° from the network's preliminary mechanism")
186 polarities to 298 km; pick classes {0: 6, 1: 30, 2: 150}
our FPFIT on them returns (140, 75, -140): 0.0° from the network's preliminary mechanism
Given the network's readings, our FPFIT returns the network's mechanism exactly: the preliminary product is FPFIT on those readings, nothing more. And the readings are poor. Most of them are the network's third class of pick, which FPFIT itself assumes is wrong nearly a quarter of the time; the thirty-seven within 30 km are barely more than half impulsive; and on those HASH gives the result its lowest grade, with an acceptable set tens of degrees wide. That our FPFIT moves by thirty degrees when the distant readings are dropped says the same thing: these readings do not pin the mechanism down. The waveform tensors explain fewer of the near readings than the first-motion solutions do — but those solutions were fitted to them.
So what can a first-motion mechanism from readings like these be trusted to say? On real data there is no truth to check against. A twin supplies one.
B.3 · A twin of these first motions¶
The twin is this earthquake's own geometry — the rays to all 41 stations that recorded it, since in a twin every station is honest — with a mechanism we plant, and readings that go wrong the way real Geysers readings go wrong. That last part has to be measured, not assumed.
A loss should match how the data actually fail. The network has published a first-motion mechanism for hundreds of Geysers earthquakes, so we can ask how often a real reading disagrees with it — by how close the ray is to a nodal plane, and by the analyst's own note on the arrival: impulsive (a sharp onset) or emergent (one that grows out of the noise). We measure it here so that the twin is noisy in the way the Geysers readings are.
pub = lab[lab.strike.notna()].copy()
pub["q"] = np.where(pub.remark == "i", 0, 1) # 0 impulsive, 1 emergent
EVENTS = [g for _, g in pub.groupby("event_id") if len(g) >= 20 and (g.q == 0).sum() >= 10]
rec = []
for g in EVENTS:
A = ray_rows(g.azimuth.values, g.takeoff_angle.values) @ mt_to_vec(sdr_to_mt(g.strike.iloc[0], g.dip.iloc[0], g.rake.iloc[0]))
rec.append(pd.DataFrame(dict(q=g.q.values, a=np.abs(A), bad=(np.sign(A) != g.y.values).astype(float),
station=g.station.values)))
rec = pd.concat(rec)
EDGES = np.array([0, 0.1, 0.2, 0.3, 0.4, 0.55, 0.7, 0.85, 1.01])
rec["bin"] = np.digitize(rec.a, EDGES[1:-1])
RATE = rec.groupby(["q", "bin"]).bad.mean() # the measured error rate
EMERGENT = rec.groupby("bin").q.mean() # the measured share of emergent readings
STATION = rec.groupby("station").bad.mean() / rec.bad.mean() # each station's reliability, relative
print(f"{len(EVENTS)} earthquakes, {len(rec):,} readings; wrong overall {100 * rec.bad.mean():.1f} %")
print(f"impulsive: {100 * RATE[0].iloc[0]:.0f} % wrong on a nodal plane, {100 * RATE[0].iloc[-1]:.0f} % far from one; "
f"emergent: {100 * RATE[1].iloc[0]:.0f} % and {100 * RATE[1].iloc[-1]:.0f} %")
from sklearn.linear_model import LogisticRegression
mid = 0.5 * (EDGES[:-1] + EDGES[1:])
fig, ax = plt.subplots(figsize=(6, 3.6))
for q, col, lab_ in ((0, "#2c3e50", "impulsive"), (1, "#e67e22", "emergent")):
t = rec[rec.q == q]
r = t.groupby("bin").bad.agg(["mean", "size"])
se = np.sqrt(r["mean"] * (1 - r["mean"]) / r["size"])
ax.errorbar(mid[r.index], 100 * r["mean"], yerr=196 * se, fmt="o", color=col, ms=5, capsize=2,
label=f"{lab_} readings ({len(t):,})")
fit = LogisticRegression(C=1e4).fit(t[["a"]], t.bad)
xs = np.linspace(0, 1, 100)
ax.plot(xs, 100 * fit.predict_proba(pd.DataFrame(dict(a=xs)))[:, 1], color=col, lw=1.2)
ax.set(xlabel="|A| the published mechanism predicts at the station (0 = on a nodal plane)",
ylabel="readings that disagree (%)", ylim=(0, 55),
title="Real readings go wrong most often near a nodal plane,\nand emergent ones far more often (lines: logistic fits)")
ax.legend(fontsize=8, frameon=False)
plt.show()
410 earthquakes, 12,937 readings; wrong overall 18.1 % impulsive: 20 % wrong on a nodal plane, 9 % far from one; emergent: 44 % and 29 %
Both effects are strong. Near a nodal plane the true first swing is small, the noise is comparable to it, and the reading approaches a coin flip; an emergent arrival is harder to read at any distance from the plane. The smooth fall with distance from the plane is the shape logistic regression assumes (the lines), and it is the fact FPFIT's \(\sqrt{|A|}\) weight was designed around. (These rates are measured against the network's own mechanisms, which were fitted to these same readings, so they are if anything optimistic.)
Does FPFIT's weighting help, on readings like these?¶
Down-weighting emergent readings is a bet that they are the unreliable ones. Test the bet on a planted world: on this earthquake's rays, plant a mechanism, let small first swings be emergent more often than large ones, at the rates just measured, and vary how often an emergent reading is actually wrong, from never to 40 % of the time. Impulsive readings are kept wrong only 2 % of the time, a choice made so the emergent readings' effect stands alone. Compare counting wrong readings against FPFIT.
rays = polarities().query("event_id == 'nc72737985'").sort_values("station")
AZ, IH = rays.azimuth.to_numpy(), rays.takeoff_angle.to_numpy()
G = ray_rows(AZ, IH)
AV5 = G @ V5.T # rays x trials: the radiation
truths = [(rng.uniform(0, 360), np.degrees(np.arccos(rng.uniform(0, 1))), rng.uniform(-180, 180)) for _ in range(150)]
print(f"measured share of emergent readings: {EMERGENT.iloc[0]:.0%} nearest a nodal plane, {EMERGENT.iloc[-1]:.0%} farthest from one")
BET = {}
for p_em in (0.0, 0.1, 0.2, 0.3, 0.4):
err_count, err_fpfit = [], []
for tr in truths:
Mt = sdr_to_mt(*tr)
A = G @ mt_to_vec(Mt)
emergent = rng.random(len(A)) < EMERGENT.reindex(np.digitize(np.abs(A), EDGES[1:-1])).values
y = np.sign(A) * np.where(rng.random(len(A)) < np.where(emergent, p_em, 0.02), -1, 1)
wrong = (np.sign(AV5) != y[:, None]).sum(0)
err_count.append(kagan(sdr_to_mt(*SDR5[np.argmin(wrong)]), Mt))
err_fpfit.append(kagan(sdr_to_mt(*SDR5[fpfit_simple(AZ, IH, y, np.where(emergent, 2, 0))]), Mt))
BET[p_em] = (np.array(err_count), np.array(err_fpfit))
print(f"emergent readings wrong {100 * p_em:3.0f} %: median error, counting {np.median(err_count):5.1f}° "
f"FPFIT {np.median(err_fpfit):5.1f}°")
fig, ax = plt.subplots(figsize=(6, 3.6))
ps = sorted(BET)
for i, (lab_, col) in enumerate((("count the wrong readings", "#7f8c8d"), ("FPFIT (emergent weighted down, √|A|)", "#c0392b"))):
med = np.array([np.median(BET[p][i]) for p in ps])
lo = np.array([np.percentile(BET[p][i], 25) for p in ps]); hi = np.array([np.percentile(BET[p][i], 75) for p in ps])
ax.errorbar(np.array(ps) * 100 + (i - 0.5) * 1.2, med, yerr=[med - lo, hi - med], fmt="o-", color=col, capsize=2, lw=1.5, label=lab_)
ax.set(xlabel="how often an emergent reading is wrong (%) [impulsive: 2 %]", ylabel="Kagan angle to the truth (°)",
title="On this earthquake's rays, FPFIT's weights barely matter\n(median and quartiles of 150 planted mechanisms)")
ax.legend(fontsize=7, frameon=False, loc="upper left")
plt.show()
measured share of emergent readings: 19% nearest a nodal plane, 15% farthest from one
emergent readings wrong 0 %: median error, counting 7.9° FPFIT 8.8°
emergent readings wrong 10 %: median error, counting 8.9° FPFIT 8.6°
emergent readings wrong 20 %: median error, counting 9.5° FPFIT 11.3°
emergent readings wrong 30 %: median error, counting 11.9° FPFIT 10.4°
emergent readings wrong 40 %: median error, counting 11.0° FPFIT 9.8°
The bet barely matters here. At every error rate the two methods land within about three degrees of each other, now one ahead and now the other, which is the scatter of 150 planted mechanisms rather than a trend. The reason is in the printed share: at The Geysers most readings are impulsive, even near a nodal plane, so down-weighting the emergent few changes little on a network of this size. A weight is a claim about the noise, and it pays only when the claim is true and the readings it reweights are numerous enough to move the answer. In machine learning these are sample weights, and the same rule holds there.
Do the error bars contain the truth?¶
The test that matters for an error bar. Plant 60 mechanisms on this earthquake's rays, flip readings as the real ones flip, and make the take-off angles we believe wrong by 0, 5 or 10 degrees. For each, ask whether FPFIT's reported 90 % range of strike, dip and rake contains the truth, and whether HASH's acceptable set does; and sort HASH's answers by its own grade.
def fpfit_box_contains(f, M):
"""Is either nodal plane of M inside FPFIT's reported strike/dip/rake ranges?"""
for s, d, r_ in nodal_planes(M):
ds, dr = abs((s - f["strike"] + 180) % 360 - 180), abs((r_ - f["rake"] + 180) % 360 - 180)
if ds <= max(f["d_strike"], 5) and abs(d - f["dip"]) <= max(f["d_dip"], 5) and dr <= max(f["d_rake"], 10):
return True
return False
cover = []
for sig in (0.0, 5.0, 10.0):
for _ in range(60):
tr = (rng.uniform(0, 360), np.degrees(np.arccos(rng.uniform(0, 1))), rng.uniform(-180, 180))
Mt = sdr_to_mt(*tr)
A = G @ mt_to_vec(Mt)
emergent = rng.random(len(A)) < EMERGENT.reindex(np.digitize(np.abs(A), EDGES[1:-1])).values
q = emergent.astype(int)
p_wrong = np.array([RATE.get((qq, b), 0.1) for qq, b in zip(q, np.digitize(np.abs(A), EDGES[1:-1]))])
y = np.sign(A) * np.where(rng.random(len(A)) < p_wrong, -1, 1)
az_b, ih_b = AZ + rng.normal(0, sig, len(AZ)), IH + rng.normal(0, sig, len(IH)) # the angles we believe
f = gm.fpfit(az_b, ih_b, y, np.where(emergent, 2, 0))[0]
az_mc, ih_mc = angle_trials(az_b, ih_b)
sols, idx = gm.hash_mech(az_mc, 180 - ih_mc, y, q, _grid=HASH_GRID)
h = sols[0]
n_t, s_t = gm.sdr_vectors(*nodal_planes(Mt)[0])
cover.append(dict(sig=sig, fpfit_box=fpfit_box_contains(f, Mt),
hash_set=bool(gm.mech_rot(n_t, s_t, HASH_GRID[1][idx], HASH_GRID[0][idx])[0].min() <= 7.5),
hash_err=kagan(sdr_to_mt(*h["sdr"]), Mt), grade=h["quality"], rms=float(np.mean(h["rms"])),
share=len(idx) / len(HASH_GRID[0])))
cover = pd.DataFrame(cover)
print((100 * cover.groupby("sig")[["fpfit_box", "hash_set"]].mean()).round(0).rename(
columns=dict(fpfit_box="FPFIT's 90 % box contains the truth (%)", hash_set="HASH's set contains it (%)")).to_string())
print(f"HASH's set is wide: median rms spread {cover.rms.median():.0f}°, median {cover.share.median():.0%} of all mechanisms on the grid")
print(cover.groupby("grade").hash_err.agg(["size", "median"]).round(1).rename(
columns=dict(size="events", median="median error (°)")).to_string())
FPFIT's 90 % box contains the truth (%) HASH's set contains it (%)
sig
0.0 47.0 98.0
5.0 45.0 100.0
10.0 37.0 100.0
HASH's set is wide: median rms spread 32°, median 7% of all mechanisms on the grid
events median error (°)
grade
A 10 18.0
B 84 21.2
C 63 26.4
D 23 39.0
fig, axes = plt.subplots(1, 2, figsize=(10, 3.7))
g_ = cover.groupby("sig")[["fpfit_box", "hash_set"]].mean() * 100
w_ = 1.6
axes[0].bar(g_.index - w_ / 2, g_.fpfit_box, w_, color="#c0392b", label="FPFIT's 90 % box")
axes[0].bar(g_.index + w_ / 2, g_.hash_set, w_, color="#e67e22", label="HASH's acceptable set")
axes[0].axhline(90, color="k", ls="--", lw=1)
axes[0].text(-2.4, 91.5, "90 %", fontsize=7)
axes[0].set(xlabel="error in the take-off angles (°)", ylabel="contains the truth (%)", xticks=[0, 5, 10], ylim=(0, 110),
title="An error bar that ignores errors in G is too narrow")
axes[0].legend(fontsize=7, frameon=False, loc="lower left")
order = [q for q in "ABCD" if (cover.grade == q).any()]
data = [cover.loc[cover.grade == q, "hash_err"].values for q in order]
axes[1].boxplot(data, tick_labels=[f"{q}\n({len(v)})" for q, v in zip(order, data)], showfliers=False,
medianprops=dict(color="#e67e22", lw=2))
axes[1].set(xlabel="HASH quality grade (number of planted events)", ylabel="Kagan angle to the truth (°)",
title="HASH's grade tells you which answers to trust")
plt.show()
FPFIT's range is computed from the readings alone, so it knows nothing of the errors in the angles or of how often real readings are wrong near a nodal plane, and it contains the truth only a third to a half of the time against the 90 % it claims. HASH's set contains the truth essentially always, because the angle errors were part of how it was built; the price is a wide set, and it says so. And HASH's grade does what a confidence should do: the better the grade, the closer the answer tends to be to the truth, from grade A to grade D. In machine-learning terms that is selective prediction — knowing which of your answers to trust.
Neither program gets a better single mechanism from the extra work; the point is an honest error bar. An ensemble's first job is to say how sure it is.
B.4 · The real readings again¶
The twin says two things about first-motion answers from readings like these: FPFIT's reported range is too narrow to trust, and HASH's acceptable set is honest, with a grade that says which answers are good. This earthquake's HASH grade is the lowest. So the right question about its sixty-degree disagreement is not "which answer is right" but "do these readings rule the waveform answers out?" — is either waveform tensor inside the set of mechanisms these readings allow?
H5_SET = gm.hash_mech(az5_mc, 180 - ih5_mc, g5.y.values.astype(float), (q5 > 0).astype(int), _grid=HASH_GRID)[1]
def beyond_set(M):
"""How far M lies from the nearest member of this event's HASH acceptable set (0 = inside)."""
n, s_ = gm.sdr_vectors(*nodal_planes(M)[0])
return float(gm.mech_rot(n, s_, HASH_GRID[1][H5_SET], HASH_GRID[0][H5_SET])[0].min())
print(f"HASH's acceptable set for this event: {len(H5_SET):,} mechanisms ({len(H5_SET) / len(HASH_GRID[0]):.0%} of the grid); the nearest member is "
f"{beyond_set(FIRST_MOTION):.0f}° from the network's first-motion mechanism, "
f"{beyond_set(US):.0f}° from the USGS tensor and {beyond_set(NC):.0f}° from Berkeley's")
HASH's acceptable set for this event: 4,824 mechanisms (16% of the grid); the nearest member is 1° from the network's first-motion mechanism, 8° from the USGS tensor and 13° from Berkeley's
The readings weakly exclude the waveform answers — nothing stronger. By the standard of the twin above, where a mechanism counted as inside the set if a member lay within 7.5° of it, both waveform tensors sit just outside: the USGS's by a little, Berkeley's by a little more. The twin says the set almost never misses a truth, so that is evidence against the waveform tensors — but weak evidence. The set is wide, about a sixth of every possible mechanism; the tensors are not double couples (half or more of each is CLVD), so the double couple compared here is only their nearest one; and the twin assumed the readings and angles go wrong no worse than the measured rates. The first-motion solutions explain more of the readings than the waveform tensors do, but they were fitted to them.
So this earthquake is the clearest example of how the two data types differ: a first-motion mechanism is only as good as the first swing of the \(P\) wave at each station, while the waveform tensor reads two and a half minutes of long-period ground motion. The rest of this notebook inverts those waveforms ourselves, and asks how certain that answer is.
Waveforms¶
B.5 · The broadband records, and a first inversion¶
A seismogram is also \(d = Gm\)¶
A.1 used the \(P\) amplitude at one instant. A whole seismogram works the same way, and for the same reason: the equations of motion are linear, so the ground displacement at a station is a sum of the six moment-tensor components, each multiplying a waveform rather than a number.
The \(g_k(t)\) are Green's functions: the seismogram that a unit \(M_{xx}\), a unit \(M_{yy}\) and so on would each produce at that station, through the Earth between them. So the design matrix has one row per sample rather than per ray, and its columns are seismograms. Everything from A.2 to A.7 — the normal equations, the SVD, damping, the deviatoric constraint, the two error bars — carries over unchanged.
Computing \(g_k(t)\) for a layered Earth is a piece of numerical seismology in its own right; the Green's functions here were computed with Saikia's frequency–wavenumber program for gil7, the northern California velocity model Berkeley uses for its own moment tensors. What arrives is not six seismograms per station but ten fundamental Green's functions, from which the six columns are built by rotating to the station's azimuth — because the Earth model is the same in every direction, so only the azimuth changes between stations.
from geysers_data import green_functions, mw5_waveforms
Z = green_functions()
W = mw5_waveforms()
STATIONS = list(W["station"])
DIST, AZ = W["distance_km"], W["azimuth"]
print(f"Green's functions: {list(Z.files)[:4]} ... for "
f"{len(Z['depth_km'])} source depths x {len(Z['distance_km'])} distances, "
f"{Z['ZDD'].shape[2]:,} samples at {float(Z['dt']):.1f} s")
print(f"waveforms: {W['waveform'].shape} (station, component, sample) at "
f"{float(W['sampling_rate']):.0f} samples per second, in centimetres of displacement")
print(f"\n{'station':>8} {'distance':>10} {'azimuth':>9}")
for s, d_, a_ in zip(STATIONS, DIST, AZ):
print(f"{s:>8} {d_:9.1f} km {a_:8.1f}°")
Green's functions: ['depth_km', 'distance_km', 'station', 'dt'] ... for 10 source depths x 8 distances, 2,048 samples at 0.1 s
waveforms: (8, 3, 4400) (station, component, sample) at 4 samples per second, in centimetres of displacement
station distance azimuth
HOPS 27.6 km 313.6°
MNRC 35.2 km 79.6°
CVS 62.6 km 147.7°
MCCM 75.3 km 182.6°
GASB 93.1 km 6.6°
SUTB 101.9 km 63.4°
BKS 117.6 km 153.1°
ORV 141.5 km 54.5°
KEYS = ("TSS", "TDS", "RSS", "RDS", "RDD", "ZSS", "ZDS", "ZDD", "REX", "ZEX")
def mt_design(gf, azimuth):
"""The design matrix for one station: rows [T; R; Z] in time, columns the moment tensor.
This is the construction TDMT (Dreger, 2003) uses for Berkeley's own solutions. The column order is [Mxx, Myy, Mxy, Mxz, Myz, Mzz]. Every entry is a fundamental Green's
function times a trigonometric function of the azimuth: the source radiates a fixed pattern and
each station reads one slice of it.
"""
a = np.radians(azimuth)
s1, c1, s2, c2 = np.sin(a), np.cos(a), np.sin(2 * a), np.cos(2 * a)
g = {k: np.asarray(gf[k], float) for k in KEYS}
zero = np.zeros_like(g["TSS"])
T = np.stack([0.5 * s2 * g["TSS"], -0.5 * s2 * g["TSS"], -c2 * g["TSS"],
-s1 * g["TDS"], c1 * g["TDS"], zero], axis=1)
R = np.stack([g["RDD"] / 6 - 0.5 * c2 * g["RSS"] + g["REX"] / 3,
g["RDD"] / 6 + 0.5 * c2 * g["RSS"] + g["REX"] / 3,
-s2 * g["RSS"], c1 * g["RDS"], s1 * g["RDS"],
g["REX"] / 3 - g["RDD"] / 3], axis=1)
Zc = np.stack([g["ZDD"] / 6 - 0.5 * c2 * g["ZSS"] + g["ZEX"] / 3,
g["ZDD"] / 6 + 0.5 * c2 * g["ZSS"] + g["ZEX"] / 3,
-s2 * g["ZSS"], c1 * g["ZDS"], s1 * g["ZDS"],
g["ZEX"] / 3 - g["ZDD"] / 3], axis=1)
return np.concatenate([T, R, Zc], axis=0)
def mt_tensor(b, scale=-1.0e20):
"""The solution vector as a tensor in north-east-down, in dyne-centimetres.
`scale` carries two things: the Green's functions are for a 1e20 dyne-cm source, and the minus
one maps the convention they are computed in onto the one used everywhere else here.
"""
mxx, myy, mxy, mxz, myz, mzz = scale * np.asarray(b, float)
return np.array([[mxx, mxy, mxz], [mxy, myy, myz], [mxz, myz, mzz]])
# the same self-check as A.1, one level up: build a synthetic from a known tensor and recover it
gf_demo = {k: Z[k][6, 3] for k in KEYS} # 7 km depth, the fourth station
G_demo = mt_design(gf_demo, float(AZ[3]))
m_demo = np.array([0.4, -0.9, 0.3, 0.2, -0.5, 0.5])
b_demo = np.linalg.lstsq(G_demo, G_demo @ m_demo, rcond=None)[0]
print("recovering a known tensor from its own synthetic seismogram:",
f"largest error {np.abs(b_demo - m_demo).max():.2e}")
recovering a known tensor from its own synthetic seismogram: largest error 6.66e-16
Preparing the data and the Green's functions the same way¶
Two rules decide everything about this step.
Both sides get the identical filter. A moment tensor inversion compares a recorded seismogram with a synthetic one, and any filter applied to one must be applied to the other, or the comparison is between different things. Here that is a band-pass from 0.02 to 0.06 Hz — periods of 17 to 50 seconds. Long periods are used deliberately: at those wavelengths the Earth looks smooth enough for a one-dimensional model like gil7 to be right, while at 1 Hz it certainly is not.
The taper must not touch the signal. A record already begins 100 s before the origin time, so it is filtered and then cut. A Green's function begins at the origin time, so it is given 100 s of leading zeros first and the taper lands in those. Getting this wrong destroys the first half-minute of every trace, which is exactly the part that carries the arrival.
import obspy
FMIN, FMAX, NPTS, PAD = 0.02, 0.06, 150, 100.0
SR, T0 = float(W["sampling_rate"]), float(W["start_offset_s"])
GDT, DEPTHS = float(Z["dt"]), Z["depth_km"]
def band(x, dt, prepend=False):
"""Band-pass, resample to one sample per second, then drop the leading `PAD` seconds."""
y = np.asarray(x, float)
if prepend:
y = np.concatenate([np.zeros(int(round(PAD / dt))), y])
tr = obspy.Trace(y, header={"delta": dt})
tr.detrend("linear")
tr.taper(0.05)
tr.filter("bandpass", freqmin=FMIN, freqmax=FMAX, corners=3, zerophase=True)
tr.resample(1.0, no_filter=True)
return tr.data[int(round(PAD)):]
REC = {s: np.array([band(W["waveform"][i, j][int(round((-PAD - T0) * SR)):], 1.0 / SR)
for j in range(3)]) for i, s in enumerate(STATIONS)}
print(f"records prepared: {len(REC)} stations x 3 components, "
f"{REC[STATIONS[0]].shape[1]} samples at 1 s from the origin time")
records prepared: 8 stations x 3 components, 900 samples at 1 s from the origin time
One more practical step, and it is not a detail. A one-dimensional velocity model never has exactly
the right travel times, so the synthetic arrives a few seconds early or late. Every implementation
of this method allows a time shift per station, and how it is chosen matters. A shift chosen by
fitting each station on its own, with all six components free, is easily fooled, because such a free
fit can match almost any shift. So we start from that per-station guess and then re-pick every
station's shift against the joint fit of all stations, three times. B.6 tests this rule on a twin
where the true shifts are known. (Berkeley's code calls the shift Zcor and chooses it by
cross-correlation against the joint solution in the same spirit.)
NLAG = 30
def station_blocks(stations, depth, npts=NPTS):
"""Per station: the design matrix, and the record at every candidate shift (0..NLAG-1 s)."""
k = int(np.argmin(np.abs(DEPTHS - depth)))
out = []
for s in stations:
i = STATIONS.index(s)
G = mt_design({n: band(Z[n][k, i], GDT, prepend=True)[:npts] for n in KEYS}, float(AZ[i]))
D = np.stack([REC[s][:, z:z + npts].ravel() for z in range(NLAG)])
out.append((s, G, D, DIST[i]))
return out
def invert(stations, depth, repick=3, model="full"):
"""Weighted least squares with per-station shifts: each station's own best shift first, then the
shifts re-picked against the joint fit. Returns the solution vector, its variance reduction, the
stacked system and the shifts."""
blocks = station_blocks(stations, depth)
rmin = min(b[3] for b in blocks)
own = [int(np.argmax([1 - np.sum((d - G @ np.linalg.lstsq(G, d, rcond=None)[0]) ** 2) / np.sum(d ** 2)
for d in D])) for _, G, D, _ in blocks]
shifts = list(own)
for it in range(repick + 1):
G_all = np.vstack([G for _, G, _, _ in blocks])
d_all = np.concatenate([D[z] for (_, _, D, _), z in zip(blocks, shifts)])
w_all = np.concatenate([np.full(len(D[0]), r / rmin) for _, _, D, r in blocks])
if model == "deviatoric":
A = np.column_stack([G_all[:, 0] - G_all[:, 5], G_all[:, 1] - G_all[:, 5], G_all[:, 2], G_all[:, 3], G_all[:, 4]])
c = np.linalg.solve(A.T @ (A * w_all[:, None]), A.T @ (w_all * d_all))
b = np.array([c[0], c[1], c[2], c[3], c[4], -(c[0] + c[1])])
else:
b = np.linalg.solve(G_all.T @ (G_all * w_all[:, None]), G_all.T @ (w_all * d_all))
if it < repick:
shifts = [int(np.argmin(((D - G @ b) ** 2).sum(1))) for _, G, D, _ in blocks]
vr = 1 - np.sum(w_all * (d_all - G_all @ b) ** 2) / np.sum(w_all * d_all ** 2)
return b, float(vr), (G_all, d_all, w_all), dict(zip(stations, zip(own, shifts)))
A first inversion¶
The source depth is inside \(G\), not beside it: change the depth and every Green's function changes, so depth is searched — invert at each depth and keep the one that fits best. With the shift rule most codes begin from, each station's own fit:
first = []
for dep in DEPTHS:
b_, vr_, _, _ = invert(STATIONS, float(dep), repick=0)
first.append((float(dep), vr_, mt_tensor(b_)))
print(f"{'depth':>6} {'variance reduction':>20} {'Mw':>6} {'Kagan to Berkeley':>19}")
for dep, vr_, M_ in first:
print(f"{dep:5.0f} km {100*vr_:17.1f} % {mt_decompose(M_)['Mw']:6.2f} {kagan(M_, NC):17.1f}°")
fit_depth, fit_vr, M_fit = max(first, key=lambda t: t[1])
print(f"\nbest-fitting depth {fit_depth:.0f} km: variance reduction {100*fit_vr:.1f} %, Mw {mt_decompose(M_fit)['Mw']:.2f}, "
f"{kagan(M_fit, NC):.0f}° from Berkeley's tensor and {kagan(M_fit, US):.0f}° from the USGS's; "
f"the hypocentre is at {g5.event_depth_km.iloc[0]:.1f} km")
depth variance reduction Mw Kagan to Berkeley
1 km 42.2 % 5.02 21.6°
2 km 46.4 % 4.95 24.3°
3 km 63.2 % 5.10 28.2°
4 km 67.5 % 5.12 34.2°
5 km 69.3 % 5.10 38.2°
6 km 72.6 % 5.12 20.5°
7 km 78.0 % 5.13 15.5°
8 km 76.9 % 5.14 11.1°
9 km 78.0 % 5.16 11.7°
10 km 75.3 % 5.17 15.4°
best-fitting depth 7 km: variance reduction 78.0 %, Mw 5.13, 16° from Berkeley's tensor and 21° from the USGS's; the hypocentre is at 1.5 km
A good fit, a moment close to the agencies', and a mechanism about twenty degrees from theirs — a little more than the agencies are from each other, and far from the first-motion answer. But the answer rests on choices: a depth chosen because it fits best, deep below where the earthquake was located; a rule for the time shifts; a velocity model that is certainly wrong. And the tensor has a large non-double-couple part. Are those choices sound, and is that part real? On real data nothing can check. A twin can.
B.6 · A twin of these records: TDMT and CAP with a wrong velocity model¶
The twin is this earthquake's own recording geometry — the same eight stations at the same distances and azimuths — with a source we plant, at 7 km. The records are made not in gil7, the model the inversion assumes, but in Part A's three test Earths: for each station its own layered model, gil7 with the top 5 km slower by a random 0–15 % and the crust from 5 to 25 km changed by a random amount up to 5 %, drawn three times (Part A, A.12, tabulates them). The inversion still assumes gil7.
from geysers_data import true_earth_green
E = true_earth_green()
AZ_W, DIST_W, DT = AZ, DIST, GDT # the same stations, on the Green's functions' clock
K7 = int(np.argmin(np.abs(DEPTHS - 7.0)))
GF = [{k: Z[k][K7, i] for k in KEYS} for i in range(8)] # gil7: what the inversion assumes
TT = [gm.travel_times(dd, 7.0) for dd in DIST_W] # gil7 P and S times
TRUTH = (335.0, 83.0, -165.0)
M_true = sdr_to_mt(*TRUTH)
M0 = 10 ** (1.5 * 5.0 + 16.1) # dyne-cm, Mw 5.0
def design_vector(M_):
"""A north-east-down tensor in the Green's functions' order and convention."""
return -np.array([M_[0, 0], M_[1, 1], M_[0, 1], M_[0, 2], M_[1, 2], M_[2, 2]]) / 1e20
def synthetic(M_, earth=None, noise=0.02, seed=0):
"""Records of a source M_ at the 8 stations, made in gil7 (earth=None) or in a test Earth."""
r = np.random.default_rng(seed)
b_ = design_vector(M_)
out = []
for i in range(8):
g = GF[i] if earth is None else {k: E[k][earth, i] for k in KEYS}
tr = {c: gm.wave_design(g, AZ_W[i])[c] @ b_ for c in "TRZ"}
pk = max(np.abs(v).max() for v in tr.values())
out.append({c: v + r.normal(0, noise * pk, len(v)) for c, v in tr.items()})
return out
def cap_windows(data, shifts=True):
"""CAP's windows (Part A, A.13): Pnl 0.05-0.2 Hz and +-2 s, surface waves 0.02-0.1 Hz and +-4 s."""
return [gm.cap_station(GF[i], data[i], AZ_W[i], DIST_W[i], DT, band=((0.05, 0.2), (0.02, 0.1)),
tp=TT[i][0], ts=TT[i][1], pnl_shift=2.0 if shifts else 0.0,
surf_shift=4.0 if shifts else 0.0) for i in range(8)]
print(f"three test Earths: top 5 km slower by {100 * E['shallow_slowdown'].min():.0f}-"
f"{100 * E['shallow_slowdown'].max():.0f} % path by path")
three test Earths: top 5 km slower by 0-15 % path by path
Plant ten mechanisms — this earthquake's published one and nine random ones — at Mw 5.0 and 7 km, make records in gil7 itself and in each of the three test Earths, add 2 % noise, and invert with gil7 Green's functions five ways:
- TDMT, no shift: the whole 150 s at 0.02–0.06 Hz, least squares, all six components;
- TDMT, a shift per station from its own fit: slide each station's record to where a free six-component fit of that station alone matches best, then invert — a rule that is easy to write;
- TDMT, shifts re-picked against the joint fit: start from those, then re-pick every station's shift three times against the fit to all stations — the rule B.5 used;
- CAP, windows and shifts as above;
- CAP with its shifts turned off, to isolate what the shifts do.
def tdmt(data, rule="none", band=(0.02, 0.06), npts=150, max_shift=10, repick=3):
"""Whole-trace least squares at 1 sample/s, stations weighted by distance. `rule` sets the
per-station shift: "none"; "own", where the station's own free 6-component fit is best; or
"joint", starting there and re-picked `repick` times against the fit to all stations."""
step = int(round(1.0 / DT))
lags = range(-max_shift, max_shift + 1) if rule != "none" else [0]
Gs, Ds, ws = [], [], []
for i in range(8):
D = gm.wave_design(GF[i], AZ_W[i])
G_i = np.vstack([gm.bandpass(D[c], DT, *band)[::step][:npts] for c in "TRZ"])
rec = {c: gm.bandpass(np.asarray(data[i][c]), DT, *band)[::step] for c in "TRZ"}
Gs.append(G_i)
Ds.append(np.array([np.concatenate([np.roll(rec[c], -z)[:npts] for c in "TRZ"]) for z in lags]))
ws.append(DIST_W[i] / DIST_W.min())
# each station's own free fit at every lag, scored by its variance reduction, as `invert` does
pick = [int(np.argmax(1 - ((Dz - Dz @ np.linalg.pinv(G_i).T @ G_i.T) ** 2).sum(1) / (Dz ** 2).sum(1)))
for G_i, Dz in zip(Gs, Ds)]
G_all = np.vstack(Gs)
w_all = np.concatenate([np.full(len(G_i), w) for G_i, w in zip(Gs, ws)])
for it in range((repick if rule == "joint" else 0) + 1):
d_all = np.concatenate([Dz[k] for Dz, k in zip(Ds, pick)])
b = np.linalg.solve(G_all.T @ (G_all * w_all[:, None]), G_all.T @ (w_all * d_all))
pick = [int(np.argmin(((Dz - G_i @ b) ** 2).sum(1))) for G_i, Dz in zip(Gs, Ds)]
mxx, myy, mxy, mxz, myz, mzz = -1e20 * b
return np.array([[mxx, mxy, mxz], [mxy, myy, myz], [mxz, myz, mzz]])
def cap(data, shifts=True):
"""The published CAP: a 10-degree grid with the Mw line search, then 2 degrees about its best."""
st = cap_windows(data, shifts)
r = gm.cap_search(st, 5.0, 0.1)
s0 = r["sdr"]
fine = gm.dc_grid((s0[0] - 10, s0[0] + 10, 2), (max(s0[1] - 10, 1), min(s0[1] + 10, 90), 2), (s0[2] - 10, s0[2] + 10, 2))
r = gm.cap_search(st, r["mw"], 0.0, sdr=fine)
return sdr_to_mt(*r["sdr"])
planted = [TRUTH] + [(rng.uniform(0, 360), np.degrees(np.arccos(rng.uniform(0, 1))), rng.uniform(-180, 180)) for _ in range(9)]
rows = []
for earth in (None, 0, 1, 2):
for n, tr in enumerate(planted):
Mt = sdr_to_mt(*tr)
data = synthetic(Mt * M0, earth, seed=n)
rows.append(dict(earth="gil7 itself" if earth is None else f"test Earth {earth + 1}",
tdmt_own=kagan(tdmt(data, "own"), Mt), tdmt_joint=kagan(tdmt(data, "joint"), Mt),
tdmt=kagan(tdmt(data), Mt),
cap_off=kagan(cap(data, shifts=False), Mt), cap=kagan(cap(data), Mt)))
CMP = pd.DataFrame(rows)
LAB = dict(tdmt_own="TDMT, shift from each station's own fit", tdmt_joint="TDMT, shift re-picked against the joint fit",
tdmt="TDMT, no shift",
cap_off="CAP with shifts off", cap="CAP")
print(f"median Kagan angle to the truth (°), {len(planted)} planted mechanisms each:")
print(CMP.groupby("earth", sort=False)[list(LAB)].median().round(1).rename(columns=LAB).to_string())
wrong = CMP[CMP.earth != "gil7 itself"]
print(f"\nall three test Earths together ({len(wrong)} inversions each): " +
", ".join(f"{LAB[k]} {wrong[k].median():.1f}°" for k in LAB))
print(f"CAP beats TDMT without shifts in {np.mean(wrong.cap < wrong.tdmt):.0%} of them")
median Kagan angle to the truth (°), 10 planted mechanisms each:
TDMT, shift from each station's own fit TDMT, shift re-picked against the joint fit TDMT, no shift CAP with shifts off CAP
earth
gil7 itself 4.5 4.9 5.3 1.2 1.5
test Earth 1 9.3 4.3 8.1 8.5 6.8
test Earth 2 9.2 4.8 6.8 4.0 4.0
test Earth 3 10.1 5.2 10.4 9.1 2.7
all three test Earths together (30 inversions each): TDMT, shift from each station's own fit 10.1°, TDMT, shift re-picked against the joint fit 4.5°, TDMT, no shift 7.7°, CAP with shifts off 8.0°, CAP 3.6°
CAP beats TDMT without shifts in 77% of them
fig, ax = plt.subplots(figsize=(10, 3.9))
colors = dict(tdmt_own="#16a085", tdmt_joint="#16a085", tdmt="#16a085", cap_off="#2c3e50", cap="#2c3e50")
marks = dict(tdmt_own="o", tdmt_joint="^", tdmt="s", cap_off="s", cap="o")
for gi, (earth, sub) in enumerate(CMP.groupby("earth", sort=False)):
for mi, key in enumerate(LAB):
v = sub[key].values
ax.scatter(gi * 6 + mi + rng.uniform(-0.15, 0.15, len(v)), v, s=12, color=colors[key],
marker=marks[key], alpha=0.8)
ax.plot([gi * 6 + mi - 0.3, gi * 6 + mi + 0.3], [np.median(v)] * 2, color=colors[key], lw=2.5)
ax.set_xticks([g * 6 + 2 for g in range(4)], list(CMP.earth.unique()))
for key in LAB:
ax.scatter([], [], color=colors[key], marker=marks[key], label=LAB[key])
ax.set(ylabel="Kagan angle to the truth (°)", ylim=(0, 45),
title="Planted mechanisms (dots) and medians (bars): CAP's per-window shift absorbs a wrong velocity model")
ax.legend(fontsize=7, frameon=False, ncol=2, loc="upper right")
plt.show()
In gil7 itself no shift is needed, and every method lands within about five degrees of the truth. The test Earths separate them. Over the three test Earths together, CAP is the most accurate, and TDMT with its shifts re-picked against the joint fit is close behind — the rule B.5 used on the real records. With its shifts turned off CAP falls back to the level of TDMT without shifts, so what absorbs the wrong velocity model is the shifts, not the windows alone. And the plausible rule — slide each station to where its own free fit is best — is the worst over the test Earths, worse than no shift at all: a six-component fit to one station can match almost any shift, so the chosen shift follows the noise and the wrong velocity model together. A latent variable must be fitted against the model you are testing, not against a model free enough to fit anything.
In gil7 itself, part of CAP's edge is its restriction to double couples, which is right here because every planted source is one.
B.7 · In the twin: when does the data justify one more parameter?¶
A.6 fitted with and without the isotropic column and compared the misfit. The misfit, though, always improves when a parameter is added — a model with more freedom fits the data it was fitted to at least as well. The honest question is whether the extra parameter predicts data it was not fitted to. That is the bias–variance trade-off from the model's side: too few parameters and the model cannot represent the source (bias); too many and it fits the noise (variance), and which wins depends on how much data there is.
Three nested models: a double couple (strike, dip, rake and size — four numbers), a deviatoric tensor (five: the double couple plus a CLVD, a compensated linear-vector dipole, which is what a non-planar or multi-fault rupture can look like), and the full tensor (six: plus the isotropic part). The test: fit using \(k\) of the 8 stations, predict the waveforms of each station left out, and ask how often the richer model predicts it better. Do this for every choice of \(k\) stations, for \(k\) from 2 to 7, on two twins whose truth we set: a pure double couple, and a source with a 30 % CLVD part — each in the three test Earths.
Every station's contribution is a few small matrices — \(G_i^\top G_i\), and \(G_i^\top d_i\) and \(|d_i|^2\) at every time shift — so a fit on any subset is a sum, and a held-out score a formula. The time shifts are fitted against the model being tested, as B.6 taught.
from itertools import combinations
BAND, NPTS, PRE = (0.02, 0.06), 150, 12 # 1 sample per second; shifts within +-12 s
STEP = int(round(1.0 / DT))
G_ST = [np.vstack([gm.bandpass(gm.wave_design(GF[i], AZ_W[i])[c], DT, *BAND)[::STEP][:NPTS] for c in "TRZ"]) for i in range(8)]
WEIGHT = DIST_W / DIST_W.min()
T5 = np.array([[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1], [-1, -1, 0, 0, 0]], float)
DC_TRIALS = -gm.dc_tensors(gm.dc_grid((0, 355, 5), (5, 90, 5), (-180, 175, 5))) # unit double couples
def station_sums(records):
"""Per station: G'G, and G'd and |d|^2 at every lag; records start PRE s before the origin."""
A = np.stack([g.T @ g for g in G_ST])
D = [np.stack([np.concatenate([gm.bandpass(r[c], DT, *BAND)[::STEP][z:z + NPTS] for c in "TRZ"])
for z in range(2 * PRE + 1)]) for r in records]
return A, np.stack([d @ g for d, g in zip(D, G_ST)]), np.stack([(d ** 2).sum(1) for d in D])
def fit(A, c, e, idx, model):
"""Weighted least squares on the stations `idx`, the lags re-picked three times against the fit."""
lag = np.array([np.argmax(np.einsum("lk,kj,lj->l", c[i], np.linalg.pinv(A[i]), c[i]) / e[i]) for i in range(8)])
for _ in range(4):
SA, Sc = np.einsum("s,sij->ij", WEIGHT[idx], A[idx]), np.einsum("s,sk->k", WEIGHT[idx], c[idx, lag[idx]])
if model == "full":
b = np.linalg.solve(SA, Sc)
elif model == "deviatoric":
b = T5 @ np.linalg.solve(T5.T @ SA @ T5, T5.T @ Sc)
else:
num, den = DC_TRIALS @ Sc, np.einsum("tk,kj,tj->t", DC_TRIALS, SA, DC_TRIALS)
a = np.maximum(num / den, 0)
j = int(np.argmax(a * num))
b = DC_TRIALS[j] * a[j]
lag[idx] = np.argmax(2 * c[idx] @ b - e[idx], 1) # the residual at each lag, up to a constant
return b
def heldout(A, c, e, b, j):
"""Variance reduction on station j, at the shift that suits the prediction best."""
return np.max(1 - (e[j] - 2 * c[j] @ b + b @ A[j] @ b) / e[j])
def learning_curve(sums):
"""For k = 2..7 training stations: how often deviatoric beats DC, and full beats deviatoric,
on the held-out stations, over every subset of k stations."""
out = {}
for k in range(2, 8):
wins = {"deviatoric": [], "full": []}
for idx in combinations(range(8), k):
idx = np.array(idx)
b = {m: fit(*sums, idx, m) for m in ("double couple", "deviatoric", "full")}
for j in set(range(8)) - set(idx):
vr = {m: heldout(*sums, b[m], j) for m in b}
wins["deviatoric"].append(vr["deviatoric"] > vr["double couple"])
wins["full"].append(vr["full"] > vr["deviatoric"])
out[k] = {m: np.mean(v) for m, v in wins.items()} | {"n": len(wins["full"])}
return out
def clvd_source(frac):
"""The Mw 5.0's double couple with a CLVD about its T axis carrying `frac` of the norm."""
dc = M_true / np.sqrt(2)
w, v = np.linalg.eigh(dc)
cl = v @ np.diag([-1.0, -1.0, 2.0]) @ v.T / np.sqrt(6)
return (np.sqrt(1 - frac ** 2) * dc + frac * cl) * np.sqrt(2)
def padded(records):
"""Records of the twin, with PRE s of silence before the origin so that early shifts exist."""
return [{c: np.concatenate([np.zeros(int(PRE / DT)), v]) for c, v in r.items()} for r in records]
def averaged(curves):
"""Pool several noise realisations of one twin: the win rates weighted by their counts."""
return {k: {m: np.average([c[k][m] for c in curves], weights=[c[k]["n"] for c in curves]) for m in ("deviatoric", "full")}
| {"n": sum(c[k]["n"] for c in curves)} for k in curves[0]}
LC = {}
for frac in (0.0, 0.3):
for earth in (None, 0, 1, 2):
LC[(frac, earth)] = averaged([learning_curve(station_sums(padded(synthetic(clvd_source(frac) * M0, earth,
noise=0.03, seed=10 * (3 if earth is None else earth) + s_)))) for s_ in range(4)])
print(f"each rate over 4 noise realisations; with 7 fitting stations that is {LC[(0.0, 0)][7]['n']} held-out predictions")
print("deviatoric beats double couple on held-out stations (%), by number of fitting stations 2..7:")
for frac in (0.0, 0.3):
for earth in (None, 0, 1, 2):
print(f" {'pure DC ' if frac == 0 else '30 % CLVD'}, {'gil7 itself ' if earth is None else f'test Earth {earth + 1}'}: " +
" ".join(f"{100 * LC[(frac, earth)][k]['deviatoric']:3.0f}" for k in range(2, 8)))
print("full beats deviatoric (%):")
for frac in (0.0, 0.3):
print(f" {'pure DC ' if frac == 0 else '30 % CLVD'}: " + " ".join(
f"{100 * min(LC[(frac, e)][k]['full'] for e in range(3)):3.0f}-{100 * max(LC[(frac, e)][k]['full'] for e in range(3)):3.0f}" for k in range(2, 8)))
each rate over 4 noise realisations; with 7 fitting stations that is 32 held-out predictions deviatoric beats double couple on held-out stations (%), by number of fitting stations 2..7: pure DC , gil7 itself : 46 47 47 47 48 53 pure DC , test Earth 1: 46 42 43 46 54 66 pure DC , test Earth 2: 39 40 40 41 42 62 pure DC , test Earth 3: 47 52 56 56 52 53 30 % CLVD, gil7 itself : 63 67 73 76 77 78 30 % CLVD, test Earth 1: 45 53 58 63 68 69 30 % CLVD, test Earth 2: 48 55 59 66 67 66 30 % CLVD, test Earth 3: 60 69 73 78 80 84 full beats deviatoric (%): pure DC : 35- 40 34- 39 33- 44 35- 45 33- 52 41- 50 30 % CLVD: 35- 41 34- 40 32- 44 34- 46 33- 49 28- 53
ks = np.arange(2, 8)
fig, axes = plt.subplots(1, 2, figsize=(10.5, 3.8), sharey=True)
for ax, key, title in ((axes[0], "deviatoric", "deviatoric (5) against double couple (4)"),
(axes[1], "full", "full (6) against deviatoric (5)")):
for frac, col, lab in ((0.0, "#34495e", "planted: pure double couple"), (0.3, "#27ae60", "planted: 30 % CLVD")):
r = np.array([[100 * LC[(frac, e)][k][key] for k in ks] for e in range(3)])
ax.fill_between(ks, r.min(0), r.max(0), color=col, alpha=0.2, lw=0)
ax.plot(ks, np.median(r, 0), "o-", color=col, lw=1.5, ms=4, label=lab + ", 3 test Earths")
ax.plot(ks, [100 * LC[(frac, None)][k][key] for k in ks], "--", color=col, lw=1.3, label=lab + ", gil7 itself")
ax.axhline(50, color="0.4", ls="--", lw=1)
ax.set(xlabel="stations used for the fit", ylim=(0, 100), title=title)
axes[0].set_ylabel("held-out station predicted better\nby the richer model (%)")
axes[0].legend(fontsize=7, frameon=False, loc="lower left")
fig.suptitle("A real parameter wins more often as data accumulate; a wrong velocity model lets an unneeded one win sometimes", fontsize=10)
fig.tight_layout()
plt.show()
Start with the dashed lines, gil7 itself, where the only thing an extra parameter can fit is noise. For the pure double couple the fifth parameter wins about half the time at most — a coin flip at best: what it fits does not carry over to a station it has not seen. For the source with a real CLVD part it already wins more often than not with two stations, and wins more and more often as stations are added. That is the bias–variance trade-off deciding, from data, whether a parameter belongs.
Now the bands, the test Earths. The real CLVD still rises clearly. But the pure double couple's extra parameter now wins more than it should as stations are added: the three test Earths are wrong in a way shared by every path (the shallow crust is slower everywhere), and a CLVD term can partly stand in for that shared error, so it carries over to unseen stations too. A parameter that predicts better is not necessarily a real property of the source: it may be absorbing the model's error. Seismologists know this as spurious non-double-couple components from unmodelled structure, and it is why a claimed CLVD must beat what a wrong velocity model alone can produce — which is exactly what a twin measures. The sixth parameter, the isotropic part, which neither twin has, wins about half the time at most throughout (right panel): a coin flip at best.
B.8 draws the real records' curve.
B.8 · The real records again: choices made by held-out prediction¶
The source depth is inside \(G\), not beside it: change the depth and every Green's function changes. So depth cannot be solved for linearly; it is searched. And the shift rule — each station's own fit, or re-picked against the joint fit — is a second choice. Both are nuisance choices, and B.7 said how to make them: not by which fits the data best, since a more flexible choice always fits better, but by which predicts a station it was not fitted to. For every depth and both rules: invert on seven stations and predict the eighth, for each of the eight in turn — this is leave-one-out cross-validation — and pool the eight predictions' residuals.
def heldout_vr(depth, repick):
"""Leave-one-station-out cross-validation: fit seven stations, predict the eighth at the shift that
suits the prediction, and pool the eight residuals with the fit's own distance weights. Returns
the pooled held-out variance reduction and the worst single station's."""
blocks = station_blocks(STATIONS, depth)
rmin = min(b[3] for b in blocks)
res = dat = 0.0
worst = np.inf
for j, s in enumerate(STATIONS):
bj, _, _, _ = invert([x for x in STATIONS if x != s], depth, repick=repick)
_, Gj, Dj, rj = blocks[j]
r = ((Dj - Gj @ bj) ** 2).sum(1)
k = int(np.argmin(r))
res += rj / rmin * r[k]
dat += rj / rmin * (Dj[k] ** 2).sum()
worst = min(worst, 1 - r[k] / (Dj[k] ** 2).sum())
return 1 - res / dat, worst
rows = []
for dep in DEPTHS:
for rule, rp in (("own fit", 0), ("joint fit", 3)):
b, vr, _, _ = invert(STATIONS, float(dep), repick=rp)
M = mt_tensor(b)
ho, worst = heldout_vr(float(dep), rp)
rows.append(dict(depth=float(dep), rule=rule, repick=rp, vr=vr, heldout=ho, worst=worst,
Mw=mt_decompose(M)["Mw"], to_NC=kagan(M, NC), to_US=kagan(M, US), to_FM=kagan(M, OURS_FM), M=M))
SENS = pd.DataFrame(rows)
print(SENS[["depth", "rule", "vr", "heldout", "worst", "Mw", "to_NC", "to_US"]].rename(columns=dict(
vr="fit (VR)", heldout="held-out VR", worst="worst station", to_NC="° to Berkeley", to_US="° to USGS")).round(3).to_string(index=False))
near = SENS[SENS.heldout >= SENS.heldout.max() - 0.03]
print(f"choices within 0.03 of the best held-out prediction: depths {near.depth.min():.0f}-{near.depth.max():.0f} km, "
f"{near.to_NC.min():.0f}-{near.to_NC.max():.0f}° from Berkeley's tensor")
fit_best = SENS.loc[SENS.vr.idxmax()]
ho_best = SENS.loc[SENS.heldout.idxmax()]
print(f"\nbest fit: {fit_best.depth:.0f} km ({fit_best.rule}); best held-out prediction: {ho_best.depth:.0f} km ({ho_best.rule}); "
f"the hypocentre is at {g5.event_depth_km.iloc[0]:.1f} km")
print(f"across all {len(SENS)} choices: Mw {SENS.Mw.min():.2f}-{SENS.Mw.max():.2f}, {SENS.to_NC.min():.0f}-{SENS.to_NC.max():.0f}° "
f"from Berkeley's tensor, and never closer than {SENS.to_FM.min():.0f}° to our first-motion mechanism")
depth rule fit (VR) held-out VR worst station Mw ° to Berkeley ° to USGS 1.0 own fit 0.422 0.544 -1.004 5.024 21.551 32.399 1.0 joint fit 0.630 0.547 -0.608 4.962 37.590 34.196 2.0 own fit 0.464 0.565 -0.989 4.949 24.269 27.103 2.0 joint fit 0.651 0.515 -1.534 4.991 22.647 23.636 3.0 own fit 0.632 0.630 -1.510 5.100 28.164 24.755 3.0 joint fit 0.691 0.559 -2.959 5.029 21.009 23.797 4.0 own fit 0.675 0.561 -2.772 5.118 34.182 28.184 4.0 joint fit 0.665 0.477 -4.595 5.033 22.199 28.404 5.0 own fit 0.693 0.673 -0.192 5.101 38.154 34.174 5.0 joint fit 0.721 0.673 -0.158 5.032 22.283 14.455 6.0 own fit 0.726 0.684 -0.257 5.116 20.542 27.345 6.0 joint fit 0.755 0.690 -0.293 5.117 23.695 25.673 7.0 own fit 0.780 0.712 -0.776 5.131 15.506 21.197 7.0 joint fit 0.797 0.704 -0.925 5.141 32.477 29.773 8.0 own fit 0.769 0.710 -0.655 5.142 11.060 20.478 8.0 joint fit 0.799 0.702 -0.837 5.152 30.809 27.850 9.0 own fit 0.780 0.710 -0.723 5.164 11.696 16.364 9.0 joint fit 0.815 0.706 -0.665 5.171 41.560 35.471 10.0 own fit 0.753 0.695 -0.350 5.167 15.420 16.785 10.0 joint fit 0.809 0.693 -0.635 5.179 37.706 31.931 choices within 0.03 of the best held-out prediction: depths 6-10 km, 11-42° from Berkeley's tensor best fit: 9 km (joint fit); best held-out prediction: 7 km (own fit); the hypocentre is at 1.5 km across all 20 choices: Mw 4.95-5.18, 11-42° from Berkeley's tensor, and never closer than 56° to our first-motion mechanism
Four things in that table.
Held-out prediction rules out the shallowest sources, and cannot choose among the rest. The in-sample fit rises almost steadily with depth. Held-out prediction is clearly worse for sources at 1–4 km, and then flat: every choice from 6 to 10 km lies within 0.03 of the best. The data do not pick a depth within that range. They do agree with both agencies that the long-period waves are best explained by a source well below the 1.5 km of the hypocentre, which is located from the first arrivals and describes where the rupture began.
How held-out scores are combined matters. The held-out score here pools the eight left-out stations with the same distance weights the fit uses, so it scores unseen data with the fit's own loss. The worst-station column shows why that matters: at 3 and 4 km one station is predicted very badly. A median over the stations would ignore that station and prefer 3 km; the pooled score does not. Score held-out data with the loss you fitted, or say why not.
The two shift rules predict about equally well at the depths the data prefer. At the shallow depths the joint re-pick does worse, and its worst station is predicted worse still: re-picking every shift against a joint fit only helps when that fit is good.
The orientation is uncertain by tens of degrees. Among the choices held-out prediction cannot separate, the mechanism ranges over tens of degrees from Berkeley's tensor, more than the two agencies differ from each other. That is A.7's second error bar on real data: the uncertainty in the answer from what the Green's functions assume, not from noise. Yet no choice brings the waveform answer anywhere near the first-motion one.
The rest of this notebook uses the choice with the best held-out prediction, and keeps in mind that its neighbours are as good.
best_depth, best_rp = float(ho_best.depth), int(ho_best.repick)
b, vr, (G, d, w), shifts = invert(STATIONS, best_depth, repick=best_rp)
M_ours = mt_tensor(b)
dec = mt_decompose(M_ours)
print(f"source at {best_depth:.0f} km, shifts from the {ho_best.rule}, {FMIN}-{FMAX} Hz, {NPTS} s of record at 1 s")
print(f"variance reduction {100*vr:.1f} %; Mw {dec['Mw']:.2f} DC {dec['dc']:.0f} % CLVD {dec['clvd']:.0f} % ISO {dec['iso']:.0f} %")
print(f"nodal planes {sdr(M_ours, 0)} and {sdr(M_ours, 1)}")
print(f"Kagan angle to Berkeley's: {kagan(M_ours, NC):.1f}°, to the USGS's: {kagan(M_ours, US):.1f}° "
f"(which is {kagan(NC, US):.1f}° from Berkeley's)")
source at 7 km, shifts from the own fit, 0.02-0.06 Hz, 150 s of record at 1 s variance reduction 78.0 %; Mw 5.13 DC 32 % CLVD 58 % ISO 10 % nodal planes (340, 79, -179) and (250, 89, -11) Kagan angle to Berkeley's: 15.5°, to the USGS's: 21.2° (which is 13.5° from Berkeley's)
syn = G @ b
t = np.arange(NPTS)
fig, axes = plt.subplots(len(STATIONS), 3, figsize=(10.5, 1.05 * len(STATIONS)), sharex=True)
for i, s in enumerate(STATIONS):
for j, comp in enumerate(("transverse", "radial", "vertical")):
ax = axes[i, j]
o = i * 3 * NPTS + j * NPTS
ax.plot(t, d[o:o + NPTS], color="0.15", lw=1.0)
ax.plot(t, syn[o:o + NPTS], color="#c0392b", lw=1.0)
ax.set(yticks=[], xlim=(0, NPTS))
ax.grid(False)
for sp in ("top", "right", "left"):
ax.spines[sp].set_visible(False)
if i == 0:
ax.set_title(comp, fontsize=9)
if j == 0:
ax.set_ylabel(f"{s}\n{DIST[i]:.0f} km", rotation=0, ha="right", va="center", fontsize=8)
for ax in axes[-1]:
ax.set_xlabel("seconds after the origin time (each station at its own shift)", fontsize=8)
fig.suptitle(f"Recorded (black) and the fitted tensor's prediction (red) — variance reduction {100*vr:.0f} %", fontsize=10)
fig.tight_layout()
plt.show()
The same records through CAP¶
B.6 found CAP the more robust method when the Earth model is wrong. Run the published CAP on the same eight records: Pnl and surface-wave windows, each pair slid within a few seconds, a grid over double couples and a line search in magnitude. As gCAP does, each record is first aligned on its observed \(P\) arrival where the network picked one, so the window shifts only have to absorb what the model gets wrong after that. And as for TDMT, records and Green's functions are filtered with 100 s of lead-in, so the filter's start-up lands before the signal.
CAP filters its two kinds of window separately, and B.6 used the bands that suit it on paper: Pnl at 0.05–0.2 Hz, where the body waves are strong, and surface waves at 0.02–0.1 Hz. Here each kind of window is run at that band and at TDMT's 0.02–0.06 Hz, in all four combinations, to see how much the answer depends on the choice.
from scipy.signal import detrend, resample_poly
from geysers_data import phases
ph_all = phases(2016)
ph_p = ph_all[(ph_all.event_id == 72737985) & (ph_all.phase_type == "P") & (ph_all.review_status == "manual")]
T_PICK = {r.station: (pd.Timestamp(r.phase_time) - pd.Timestamp(str(W["origin_time"])[:26])).total_seconds()
for r in ph_p.itertuples() if r.station in STATIONS}
X10 = resample_poly(W["waveform"], 5, 2, axis=2) # 4 -> 10 samples per second
LEAD = 100.0 # seconds of lead-in before the origin
def cap_real(depth, bands):
"""CAP on the real records at one depth and pair of bands: the best double couple and magnitude,
and each window's variance reduction at that answer (stations x 5)."""
k = int(np.argmin(np.abs(DEPTHS - depth)))
st, lags = [], []
for i, s in enumerate(STATIONS):
tp, ts = gm.travel_times(DIST[i], DEPTHS[k])
lag = T_PICK.get(s, tp) - tp # observed minus predicted P
i0 = int(round((-T0 + lag - LEAD) * 10))
seg = detrend(X10[i, :, i0:i0 + 2048 + int(LEAD * 10)], axis=1)
seg[:, :100] *= np.hanning(200)[:100] # taper the lead-in's first 10 s
gf = {n: np.r_[np.zeros(int(LEAD * 10)), Z[n][k, i]] for n in KEYS}
st.append(gm.cap_station(gf, dict(zip("TRZ", seg)), float(AZ[i]), float(DIST[i]), GDT,
band=bands, tp=tp + LEAD, ts=ts + LEAD))
lags.append(lag)
r = gm.cap_search(st, 5.0, 0.1) # coarse grid, magnitude line search
s0 = r["sdr"]
r = gm.cap_search(st, r["mw"], 0.0, sdr=gm.dc_grid((s0[0] - 10, s0[0] + 10, 2), (max(s0[1] - 10, 1), min(s0[1] + 10, 90), 2),
(s0[2] - 10, s0[2] + 10, 2)))
rec2 = np.array([[w["rec2"] for w in x] for x in st])
return r, 1 - r["window_err"] / rec2, lags
print("P aligned on the network's pick at " + ", ".join(f"{s} ({T_PICK[s] - gm.travel_times(DIST[i], best_depth)[0]:+.1f} s)"
for i, s in enumerate(STATIONS) if s in T_PICK))
CAP_BANDS = {("0.05-0.2", "0.02-0.1"): ((0.05, 0.2), (0.02, 0.1)),
("0.02-0.06", "0.02-0.1"): ((FMIN, FMAX), (0.02, 0.1)),
("0.05-0.2", "0.02-0.06"): ((0.05, 0.2), (FMIN, FMAX)),
("0.02-0.06", "0.02-0.06"): ((FMIN, FMAX), (FMIN, FMAX))}
cap_runs = {}
print(f"{'Pnl band (Hz)':>14s} {'surface band (Hz)':>18s} Mw ° to Berkeley ° to USGS ° to our TDMT")
for key, bands in CAP_BANDS.items():
r, _, _ = cap_real(best_depth, bands)
M_ = sdr_to_mt(*r["sdr"])
cap_runs[key] = M_
print(f"{key[0]:>14s} {key[1]:>18s} {r['mw']:5.2f} {kagan(M_, NC):12.1f} {kagan(M_, US):9.1f} {kagan(M_, M_ours):12.1f}")
P aligned on the network's pick at HOPS (+0.2 s), MNRC (+0.7 s), CVS (+0.2 s), MCCM (+1.0 s), SUTB (+0.0 s), BKS (+0.9 s), ORV (+0.2 s) Pnl band (Hz) surface band (Hz) Mw ° to Berkeley ° to USGS ° to our TDMT
0.05-0.2 0.02-0.1 4.83 42.2 37.8 31.4
0.02-0.06 0.02-0.1 4.86 46.7 42.0 35.8
0.05-0.2 0.02-0.06 4.87 39.1 33.8 29.8
0.02-0.06 0.02-0.06 5.02 21.1 33.5 27.4
The band moves CAP's answer by twenty to twenty-five degrees. With both kinds of window at 0.02–0.06 Hz, the band TDMT used, CAP lands about twenty degrees from Berkeley's tensor, within the spread of the table above; with either kind of window reaching higher frequencies it lands about forty degrees from both agencies.
No held-out test can choose the band: records filtered differently are different data, and their misfits cannot be compared. What the comparison settles is narrower. Given the same records in the same band as TDMT, CAP — a different method, restricted to double couples, with its own windows and shifts — lands within the spread that TDMT's own nuisance choices produce. B.6's twin could not have warned about the band: its Earths were layered, so the higher frequencies were modelled as well as the long periods, whereas in the real Earth shorter periods are more sensitive to the three-dimensional structure that gil7 leaves out. The rest of this notebook uses CAP at 0.02–0.06 Hz.
M_cap = cap_runs[("0.02-0.06", "0.02-0.06")]
What is resolved, and is the CLVD part real?¶
The tensor is one answer with six numbers in it, and they are not equally believable. First the conditioning, and A.6's test of the isotropic column; then B.7's test of whether the extra parameters are real. (The singular values are those of the weighted design matrix of all eight stations, as in A.3.)
sv = np.linalg.svd(G * np.sqrt(w)[:, None], compute_uv=False)
print("singular values, divided by the largest:", np.round(sv / sv[0], 3), f" condition number {sv[0]/sv[-1]:.1f}")
b_dev, vr_dev, _, _ = invert(STATIONS, best_depth, repick=best_rp, model="deviatoric")
M_dev = mt_tensor(b_dev)
print(f"\nsix free components : variance reduction {100*vr:.1f} %, isotropic part {dec['iso']:.0f} %")
print(f"deviatoric (five) : variance reduction {100*vr_dev:.1f} %, isotropic part 0 % by construction")
print(f" the extra freedom buys {100*(vr - vr_dev):.1f} points of variance and moves the mechanism {kagan(M_ours, M_dev):.1f}°")
singular values, divided by the largest: [1. 0.648 0.423 0.342 0.299 0.229] condition number 4.4 six free components : variance reduction 78.0 %, isotropic part 10 % deviatoric (five) : variance reduction 77.8 %, isotropic part 0 % by construction the extra freedom buys 0.2 points of variance and moves the mechanism 0.6°
A condition number near 4 is small: at these periods the eight stations determine all six components, so whether the isotropic part is real is not a question of resolution. Nor can the misfit settle it: the sixth column buys a fraction of a point of variance reduction, but more freedom always fits at least as well. B.7's test can. Fit using \(k\) of the eight stations, predict each station left out, and count how often the richer model predicts it better — for every subset of \(k\) stations. Each station contributes a few small matrices, so every fit is a sum.
from itertools import combinations
blocks = station_blocks(STATIONS, best_depth)
A_st = np.stack([Gb.T @ Gb for _, Gb, _, _ in blocks])
c_st = np.stack([Db @ Gb for _, Gb, Db, _ in blocks]) # G'd at every shift
e_st = np.stack([(Db ** 2).sum(1) for _, _, Db, _ in blocks])
W_ST = DIST / DIST.min()
T5 = np.array([[1, 0, 0, 0, 0], [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1], [-1, -1, 0, 0, 0]], float)
DC_TRIALS = -gm.dc_tensors(gm.dc_grid((0, 355, 5), (5, 90, 5), (-180, 175, 5)))
def fit_subset(idx, model):
"""B.7's fit on the stations `idx`: shifts from each station's own fit, re-picked against the joint fit."""
lag = np.array([np.argmax(np.einsum("lk,kj,lj->l", c_st[i], np.linalg.pinv(A_st[i]), c_st[i]) / e_st[i]) for i in range(8)])
for _ in range(4):
SA = np.einsum("s,sij->ij", W_ST[idx], A_st[idx]); Sc = np.einsum("s,sk->k", W_ST[idx], c_st[idx, lag[idx]])
if model == "full":
b_ = np.linalg.solve(SA, Sc)
elif model == "deviatoric":
b_ = T5 @ np.linalg.solve(T5.T @ SA @ T5, T5.T @ Sc)
else:
num, den = DC_TRIALS @ Sc, np.einsum("tk,kj,tj->t", DC_TRIALS, SA, DC_TRIALS)
a = np.maximum(num / den, 0); j = int(np.argmax(a * num)); b_ = DC_TRIALS[j] * a[j]
lag[idx] = np.argmax(2 * c_st[idx] @ b_ - e_st[idx], 1)
return b_
def heldout(b_, j):
return np.max(1 - (e_st[j] - 2 * c_st[j] @ b_ + b_ @ A_st[j] @ b_) / e_st[j])
LCR = {}
for k in range(2, 8):
wins = {"deviatoric": [], "full": []}
for idx in combinations(range(8), k):
idx = np.array(idx)
fits_ = {m: fit_subset(idx, m) for m in ("double couple", "deviatoric", "full")}
for j in set(range(8)) - set(idx):
v = {m: heldout(fits_[m], j) for m in fits_}
wins["deviatoric"].append(v["deviatoric"] > v["double couple"]); wins["full"].append(v["full"] > v["deviatoric"])
LCR[k] = {m: 100 * np.mean(x) for m, x in wins.items()} | {"n": len(wins["full"])}
print("held-out station predicted better (%), by number of fitting stations:")
print(f"{'fitting stations':>32s}" + "".join(f"{k:>6d}" for k in LCR))
print(f"{'deviatoric beats double couple':>32s}" + "".join(f"{LCR[k]['deviatoric']:6.0f}" for k in LCR))
print(f"{'full beats deviatoric':>32s}" + "".join(f"{LCR[k]['full']:6.0f}" for k in LCR))
print(f"{'held-out predictions':>32s}" + "".join(f"{LCR[k]['n']:6d}" for k in LCR))
fig, ax = plt.subplots(figsize=(6.4, 3.9))
ks = list(LCR)
ax.plot(ks, [LCR[k]["deviatoric"] for k in ks], "o-", color="#27ae60", lw=2.4, label="deviatoric (5) beats double couple (4)")
ax.plot(ks, [LCR[k]["full"] for k in ks], "s--", color="#d35400", lw=1.8, label="full (6) beats deviatoric (5)")
ax.axhline(50, color="0.4", ls="--", lw=1)
for k in ks:
ax.text(k, 3, f"n={LCR[k]['n']}", fontsize=6, ha="center", color="0.4")
ax.set(xlabel="stations used for the fit", ylabel="held-out station predicted better (%)", ylim=(0, 100),
title="The real Mw 5.0: the fifth parameter earns its place, the sixth does not")
ax.legend(fontsize=7, frameon=False, loc="upper left")
plt.show()
held-out station predicted better (%), by number of fitting stations:
fitting stations 2 3 4 5 6 7
deviatoric beats double couple 52 59 59 66 73 88
full beats deviatoric 40 41 44 48 52 50
held-out predictions 168 280 280 168 56 8
Read it against B.7. The isotropic part (the sixth parameter) wins held-out prediction about half the time at most, whatever the number of stations — what B.7's twins did, none of which had an isotropic part. Report it as not resolved.
The CLVD part (the fifth parameter) is different: it wins about half the time with two stations, and more and more often as stations are added — three times in four with six. With five or six stations that is above anything the test Earths produced for a pure double couple in B.7, and in the range of the twins with a real CLVD part. (With seven stations there are only eight held-out predictions, too few to read.) So the data support a non-double-couple part of this size, with the caveat B.7 made precise: a velocity model wrong in the same way on every path can mimic part of it. The Geysers is a geothermal field where such sources are expected and have been reported (Guilhem et al., 2014; Boyd et al., 2015); this is the evidence, and this is its error bar. Each rate in B.7 and here is a proportion of a few dozen to a few hundred predictions, uncertain by several points, so read the trend across the number of stations, not any single value.
Exercise. Invert with only the four nearest stations, then with only the four most distant, and print for each the variance reduction, Mw, the double-couple share and the Kagan angle to the eight-station answer and to Berkeley's. Which subset reproduces the published tensor better, and would the variance reduction have told you in advance?
order = [s for _, s in sorted(zip(DIST, STATIONS))]
for label, subset in (("the four nearest", order[:4]), ("the four most distant", order[-4:])):
bs, vrs, _, _ = invert(subset, best_depth, repick=best_rp)
Ms = mt_tensor(bs)
ds_ = mt_decompose(Ms)
print(f"{label:22s} {', '.join(subset):28s} VR {100*vrs:5.1f} % Mw {ds_['Mw']:.2f} "
f"DC {ds_['dc']:3.0f} % {kagan(Ms, M_ours):5.1f}° from all eight, {kagan(Ms, NC):5.1f}° from Berkeley")
the four nearest HOPS, MNRC, CVS, MCCM VR 78.5 % Mw 5.16 DC 47 % 26.1° from all eight, 34.2° from Berkeley
the four most distant GASB, SUTB, BKS, ORV VR 82.6 % Mw 5.13 DC 72 % 19.6° from all eight, 24.8° from Berkeley
Four stations give the same moment either way, but not the same mechanism: each half of the network moves the answer by twenty to twenty-five degrees from the eight-station one, as much as the choice of depth does. The distant four fit their records a little better and land a little closer to Berkeley's tensor. But the two variance reductions measure fits to different records, so on their own they could not have said which tensor is better; only the comparison with the full network and the agencies' does. What the exercise shows most plainly is how much the answer depends on which stations it rests on.
B.9 · The answer¶
Put everything on one focal sphere, with this earthquake's first motions plotted on top.
def focal_xy(azimuth, takeoff):
az, ih = np.asarray(azimuth, float) % 360, np.asarray(takeoff, float)
up = ih > 90
az = np.where(up, az + 180, az) % 360
ih = np.where(up, 180 - ih, ih)
r = np.sqrt(2) * np.sin(np.radians(ih) / 2)
return r * np.sin(np.radians(az)), r * np.cos(np.radians(az))
def plane_xy(strike, dip, n=181):
f, d_ = np.radians(strike), np.radians(dip)
along = np.array([np.cos(f), np.sin(f), 0.0])
down = np.array([-np.cos(d_) * np.sin(f), np.cos(d_) * np.cos(f), np.sin(d_)])
tt = np.linspace(0, 2 * np.pi, n)[:, None]
v = np.cos(tt) * along + np.sin(tt) * down
v = np.where(v[:, [2]] < 0, -v, v)
return focal_xy(np.degrees(np.arctan2(v[:, 1], v[:, 0])), np.degrees(np.arccos(np.clip(v[:, 2], -1, 1))))
def focal_sphere(ax, M, color="#2c3e50", lw=1.2, ls="-"):
for strike, dip, _ in nodal_planes(M):
px, py = plane_xy(strike, dip)
cut = np.where(np.abs(np.diff(px)) + np.abs(np.diff(py)) > 0.5)[0] + 1
ax.plot(np.insert(px, cut, np.nan), np.insert(py, cut, np.nan), color=color, lw=lw, ls=ls)
ax.add_patch(plt.Circle((0, 0), 1, fill=False, color="0.4", lw=1))
ax.set(xlim=(-1.15, 1.15), ylim=(-1.15, 1.15), xticks=[], yticks=[])
ax.set_aspect("equal"); ax.grid(False)
for sp in ax.spines.values():
sp.set_visible(False)
fig, axes = plt.subplots(1, 3, figsize=(11.2, 4.4))
x5, y5 = focal_xy(g5.azimuth.values, g5.takeoff_angle.values)
panels = [("first motions", f"network (preliminary) and our FPFIT: {kagan(FIRST_MOTION, OURS_FM):.0f}° apart",
[(FIRST_MOTION, "#2c3e50", 2.4, "-", "network (preliminary)"), (OURS_FM, "#c0392b", 1.4, (0, (5, 2)), "our FPFIT")]),
("waveforms", f"ours {kagan(NC, M_ours):.0f}° (TDMT) and {kagan(NC, M_cap):.0f}° (CAP) from Berkeley's",
[(NC, "#2c3e50", 2.4, "-", "Berkeley"), (US, "#7f8c8d", 1.6, (0, (1, 2)), "USGS"),
(M_ours, "#c0392b", 1.4, (0, (5, 2)), "our TDMT"), (M_cap, "#8e44ad", 1.4, (0, (3, 1, 1, 1)), "our CAP")]),
("the two data types", f"{kagan(OURS_FM, M_ours):.0f}° apart",
[(OURS_FM, "#c0392b", 2.2, "-", "from first motions"), (M_ours, "#2980b9", 2.2, "-", "from waveforms")])]
for ax, (title, note, items) in zip(axes, panels):
for M, c, lw, ls, label in items:
focal_sphere(ax, M, color=c, lw=lw, ls=ls)
ax.plot([], [], color=c, lw=lw, ls=ls, label=label)
ax.scatter(x5[g5.y > 0], y5[g5.y > 0], s=24, c="#2c3e50", zorder=3)
ax.scatter(x5[g5.y < 0], y5[g5.y < 0], s=24, facecolors="white", edgecolors="#2c3e50", zorder=3)
ax.set_title(f"{title}\n{note}", fontsize=9)
ax.legend(fontsize=7, loc="lower center", bbox_to_anchor=(0.5, -0.2), ncols=2, frameon=False, handlelength=1.8)
fig.suptitle(f"The Geysers Mw 5.0 of 14 December 2016, with its {len(g5)} usable first motions", fontsize=10, y=1.02)
fig.tight_layout()
plt.show()
print(f"first motions, network vs ours (near) : {kagan(FIRST_MOTION, OURS_FM):5.1f}°")
print(f"waveforms, Berkeley vs the USGS : {kagan(NC, US):5.1f}°")
print(f"waveforms, Berkeley vs our TDMT / CAP : {kagan(NC, M_ours):5.1f}° / {kagan(NC, M_cap):5.1f}°")
print(f"first motions vs waveforms, published : {kagan(FIRST_MOTION, NC):5.1f}°")
print(f"first motions vs waveforms, ours : {kagan(OURS_FM, M_ours):5.1f}°")
first motions, network vs ours (near) : 28.8° waveforms, Berkeley vs the USGS : 13.5° waveforms, Berkeley vs our TDMT / CAP : 15.5° / 21.1° first motions vs waveforms, published : 61.6° first motions vs waveforms, ours : 75.8°
Within each data type the implementations agree about as well as the data allow: our FPFIT reproduces the network's preliminary mechanism from the network's readings (B.2), and moves by thirty degrees on the near readings alone; our two waveform inversions lie within the spread of defensible choices B.8 measured, and so does the gap between the two agencies' reviewed tensors. Across the two data types the answers are far apart — farther than any choice of depth, shift rule, band or method moves either one.
How much of that is the Earth, and how much the inversion? The notebook cannot split the sixty degrees into exact shares, but it can say where each part lives. The inversion's choices — depth, shifts, band, stations — move the waveform answer by up to a few tens of degrees, because eight long-period stations do not pin them down. The first-motion answer carries the width of HASH's acceptable set, tens of degrees, because this earthquake's first motions are poor, and B.2 showed that the preliminary mechanism itself moves thirty degrees with the choice of readings. Neither spread alone reaches sixty degrees, and the first-motion readings only weakly exclude the waveform tensors (B.4). Even with good data the two need not agree exactly — they sample different parts of the rupture through different parts of the Earth, and Guilhem et al. (2014), comparing them for M ~ 3 earthquakes in this field, call them complementary but different. The machine-learning version of the lesson is the one this session kept returning to: an answer is only as good as its worst-understood input, and an honest error bar is what tells you which input that is.
What to take away¶
Machine learning¶
- A sample weight is a claim about the noise, and it pays only when the claim is true and the reweighted readings are numerous; on this earthquake's rays FPFIT's down-weighting of emergent readings barely changes the answer.
- An error bar computed from the data alone (FPFIT's range) contained the planted truth about a third of the time; one that also perturbs the inputs (HASH's set) contained it almost always, and its grade sorted good answers from bad.
- A nuisance variable such as a time shift must be fitted against the model being tested: in the twin, re-picking the shifts against the joint fit nearly matched CAP, while shifts from each station's free fit did worse than no shift at all.
- An extra parameter belongs in a model only if it predicts held-out data more often as data are added, and a twin with a wrong velocity model shows how often a spurious parameter passes that test.
- Nuisance choices such as depth and shift rule are made by held-out prediction, never by the in-sample misfit, which favours the most flexible choice; and the held-out score must use the loss that was fitted, since a median over stations hid a badly predicted one and chose a different depth.
Seismology¶
- The Mw 5.0's preliminary first-motion mechanism is FPFIT applied to the network's mostly low-quality readings; HASH grades the first-motion answer D, and the readings only weakly exclude the reviewed waveform tensors.
- On the real records CAP's answer moves by twenty to twenty-five degrees with the frequency band; at TDMT's long-period band it lands within the spread of TDMT's own choices.
- Held-out prediction rules out a source at 1–4 km but cannot choose between 6 and 10 km; like both agencies, it places the long-period source well below the 1.5 km hypocentre.
- The Mw 5.0's moment is firm and its waveform mechanism moves by tens of degrees with depth and shift choices; held-out prediction supports its CLVD part but not an isotropic part.
- Most of the sixty degrees between the first-motion and waveform answers comes from the poor first motions; the inversion's choices add a spread of tens of degrees on the waveform side.
References¶
- Boyd, O. S., Dreger, D. S., Lai, V. H., & Gritto, R. (2015). A systematic analysis of seismic moment tensor at The Geysers geothermal field, California. BSSA, 105(6), 2969–2986. https://doi.org/10.1785/0120140285
- Dreger, D. S. (2003). TDMT_INV: time domain seismic moment tensor INVersion. In International Handbook of Earthquake and Engineering Seismology (International Geophysics 81), 1627. https://doi.org/10.1016/s0074-6142(03)80290-5
- Gilbert, F. (1971). Excitation of the normal modes of the Earth by earthquake sources. Geophysical Journal International, 22(2), 223–226. https://doi.org/10.1111/j.1365-246X.1971.tb03593.x
- Guilhem, A., Hutchings, L., Dreger, D. S., & Johnson, L. R. (2014). Moment tensor inversions of M ~ 3 earthquakes in the Geysers geothermal fields, California. JGR: Solid Earth, 119, 2121–2137. https://doi.org/10.1002/2013JB010271
- Hardebeck, J. L., & Shearer, P. M. (2002). A new method for determining first-motion focal mechanisms. BSSA, 92(6), 2264–2276. https://doi.org/10.1785/0120010200
- Li, J., Zhu, W., Biondi, E., & Zhan, Z. (2023). Earthquake focal mechanisms with distributed acoustic sensing. Nature Communications, 14, 4181. https://doi.org/10.1038/s41467-023-39639-3
- Reasenberg, P., & Oppenheimer, D. (1985). FPFIT, FPPLOT and FPPAGE: Fortran computer programs for calculating and displaying earthquake fault-plane solutions. USGS Open-File Report 85-739. https://doi.org/10.3133/ofr85739
- Saikia, C. K. (1994). Modified frequency-wavenumber algorithm for regional seismograms using Filon's quadrature. Geophysical Journal International, 118(1), 142–158. https://doi.org/10.1111/j.1365-246X.1994.tb04680.x
- Skoumal, R. J., Hardebeck, J. L., & Shearer, P. M. (2024). SKHASH: a Python package for computing earthquake focal mechanisms. Seismological Research Letters, 95(4), 2519–2526. https://doi.org/10.1785/0220230329
- Zhu, L., & Helmberger, D. V. (1996). Advancement in source estimation techniques using broadband regional seismograms. BSSA, 86(5), 1634–1641. https://doi.org/10.1785/bssa0860051634
- Zhu, L., & Ben-Zion, Y. (2013). Parametrization of general seismic potency and moment tensors for source inversion of seismic waveform data. Geophysical Journal International, 194(2), 839–843. https://doi.org/10.1093/gji/ggt137