Deep learning · Part B: picking one hour of the Geysers

← Open this notebook in Colab. Part A, Deep learning: the methods: open in Colab.

Run a published picker on one hour of continuous Geysers data from 131 stations. Do its picks agree with the analysts', what does the catalogue miss, and does the sampling rate matter?

Part A built a U-Net picker and trained it on labelled event windows. Here we take the network the field actually uses, PhaseNet (Zhu & Beroza 2019), with its published weights, and run it the way it is run in practice: on continuous records, hour after hour, on every station, with nobody telling it where the earthquakes are. There are no labels for most of what it finds, so we judge it two ways: against the picks the Northern California Seismic Network's analysts made for the same hour, and against itself, by asking whether the picks on 131 stations line up into earthquakes.

You can run it without running Part A first. A GPU makes it faster (Colab: Runtime → Change runtime type → GPU), but every step also runs on a CPU.

The questions this notebook answers

1. Which data do we have for one hour, and how are they read from two public archives? (B.1)
2. How does PhaseNet turn a continuous three-component record into picks? (B.2)
3. What does one hour of picks on 131 stations look like? (B.3)
4. Do PhaseNet's picks agree with the analysts' picks? (B.4)
5. How many earthquakes does the catalogue miss in this hour? (B.5)
6. PhaseNet was trained on 100 Hz data. What should we do with data recorded at 200 or 500 Hz? (B.6)
7. Where does the picker go wrong? (B.7)

B.1 · One hour, two kinds of network

The hour. We work on one hour: from 23:30 on 19 October 2022 to 00:30 on 20 October (UTC; it was 16:30–17:30 in California). The first half is ordinary. Just before midnight an M 3.3 starts a burst — another M 3.3 within a minute, an M 3.2 three and a half minutes later, and dozens of smaller earthquakes after them (the list is printed below). A quiet half and a busy half in one hour let us see what the picker does on noise and what it does when earthquakes overlap.

The permanent network. The final project's station table (gd.sites()) lists every station inside the 40 km study circle. Most belong to the BG network, the borehole geophones Calpine runs over the field, recorded at 500 samples per second (500 Hz); the others are the USGS (NC) and Berkeley (BK) short-period, broadband and strong-motion stations, recorded at 100 Hz. A station in the table is not always a station with data: we keep those that were open on that day and have files in the NCEDC public archive. Their continuous records are stored as one file per channel per day, on Amazon S3.

The dense array. From 2018 to 2024 the California Energy Commission (CEC) ran a temporary array of 91 three-component geophones, about 500 m apart, over a 5 × 5 km patch of the field, recorded at 200 Hz. Its records are stored as one file per channel per hour, on Google Cloud Storage (gs://quakeflow_waveform/Geysers_CEC, public). The files name the nodes S00–S95; the coordinate table names them ST001–ST100. The two match number for number (S12 is ST012), except that S00 is ST100.

67 project stations open on the day, inside 40 km

Which of them have files that day? The archive answers by listing a folder. We keep one instrument per station, preferring a velocity sensor (broadband HH, borehole DP, short-period EH) to an accelerometer (HN, CN), because PhaseNet was trained on velocity records.

57 of them have files on both days: {'DP': 32, 'EH': 13, 'CN': 5, 'HH': 4, 'HN': 3}
open but no files: ['1216', '68035', '68328', '68676', 'ADS2', 'ADSP', 'COB', 'GALB', 'HOPS', 'N005']

74 of the 91 CEC nodes recorded in the hour (453 hourly files)

The catalogue and the analysts' picks. NCEDC serves its catalogue through an FDSN event service, and asked for includearrivals, it returns every pick that located each earthquake. A pick carries an evaluation_mode: manual means an analyst made or reviewed it, automatic means the network's real-time picker made it and nobody changed it.

67 catalogue events in the hour, M -0.1–3.3; 1747 picks: {('P', 'automatic'): 632, ('P', 'manual'): 1054, ('S', 'manual'): 61}
1,686 P picks and 61 S picks: one S for every 28 P
the earthquakes of M ≥ 3: M 3.3 at 23:59:04, M 3.3 at 23:59:50, M 3.2 at 00:02:34

The catalogue made 1,686 P picks for this hour and 61 S picks — one S for every 28 P. Its real-time picker makes only P picks; the S picks are the analysts'. An S pick is harder to make by hand, and on many of the network's stations only the vertical component is recorded, on which S is weak.

One order for every station. Every figure below lists stations in one order: by their distance along the line A–A′ of Session 4, Clustering, which runs north-west to south-east along the long axis of the field's seismicity. Rows close together are stations at about the same position along that axis (a station off to the side of the line is placed at the point of the line nearest to it). The 74 dense-array nodes, which sit in the middle of the line, fill most of the middle rows, mixed with the BG boreholes around them; the figures shade their rows.

131 stations in all: {'CEC node': 74, 'BG': 37, 'short-period': 11, 'broadband': 7, 'strong-motion': 2}

The map shows the two networks at two scales. On the left, the permanent stations spread over the 40 km circle; on the right, the dense array sits directly over the earthquakes of the hour.

B.2 · PhaseNet on a continuous record

What PhaseNet is. PhaseNet (Zhu & Beroza 2019) is the U-Net of Part A, trained on "more than seven million waveform samples extracted from over thirty years of earthquake recordings" with the P and S arrival times that NCEDC analysts picked. Every training record was sampled at 100 Hz and cut into 30 s windows of 3,001 samples, and each arrival was labelled with a Gaussian of standard deviation 0.1 s, as in Part A. Its weights are published, so we can run exactly the network the paper trained.

Reading the records. First, two readers, one for each archive. A CEC node's hour is two hourly files per component (23:00 on the 19th and 00:00 on the 20th), fetched whole. A network station's hour sits inside two day files of about 50 MB each, of which we need half an hour from each. A miniSEED file is a sequence of fixed 4,096-byte records, and each record's header states the time of its first sample, so we can find the records that cover our hour by bisection — read the header of the middle record, keep the half that contains the hour, repeat — and then download only those, with an HTTP Range request. A few dozen small requests replace downloading the whole day.

miniSEED record length in the NCEDC archive: 4,096 bytes

From a record to the network's input. Each step below copies what PhaseNet's own code does with a continuous record before the network sees it.

1. Merge. A day of data can arrive in several pieces; each piece has its mean removed and gaps are filled with zeros.
2. High-pass at 1 Hz, with a causal filter (it uses only past samples), so that nothing is smeared to before an onset.
3. Choose the sampling rate. PhaseNet never sees a time axis, only samples, so the rate we give it changes what it sees. We keep the recorded rate of the dense array (200 Hz) and of the 100 Hz stations, and lower the BG boreholes from 500 to 250 Hz; B.6 shows why.
4. Standardise. Each component has a running mean subtracted and is divided by a running standard deviation, both measured in windows of 3,000 samples centred every 1,500 samples and interpolated in between. The training windows were standardised as a whole; a long record cannot be, because one large earthquake would make every small one in the hour vanish.

PhaseNet's training windows: 3,001 samples (30 s at 100 Hz); standardisation of a long record: windows of 3,000 samples, one every 1,500

The network. It is Part A's UNet2018 with one difference that matters only for the published weights. The weights were trained in TensorFlow, whose "same" padding for a convolution of stride 4 and kernel 7 adds one sample on the left and two on the right, where Part A's version adds three on each side. With the wrong padding, every stage of the encoder is shifted by part of a sample; the shifts add up through the network, and the published weights then give wrong probabilities. The three padding rules below make the class reproduce PhaseNet.

The network contains only convolutions, so it accepts a record of any length. We give it the whole hour in one pass rather than cutting it into 30 s windows: no window edges, and no seams to stitch.

PhaseNet: 268,443 parameters, on mps

From probabilities to picks. The network returns, for every sample, the probability that it is noise, a P arrival or an S arrival. A pick is a local maximum of the P (or S) probability that is at least 0.3 and is the largest value within half a second on either side — the rule of Part A (PhaseNet's code writes it as 50 samples, which is half a second at 100 Hz). The threshold trades missed arrivals against false ones, as Part A showed; 0.3 is PhaseNet's default. Where we need fewer false picks (the figures, the detections of B.5), we keep only picks of 0.5 or more.

a pick: a local maximum of at least 0.3, the largest within 0.5 s either side

One node, one earthquake. We run the hour of node S56 through these steps and look at the M 1.3 of 00:18:38, whose epicentre is 0.5 km from the node — a typical size for the hour. The three components go in; two probability curves come out; their peaks are the picks.

rates given to PhaseNet: {'CEC': 200.0, 'BG': 250.0}, every other station 100 Hz; S56 is 0.5 km from the epicentre

P pick 0.39 s after the origin time, probability 0.98
S pick 0.80 s after the origin time, probability 0.72

M 1.34 at 00:18:37

The P probability rises to almost 1 at the first motion, and the S probability peaks where the horizontal components grow. The S peak is lower and wider than the P peak. Read it as the network rating S as harder to place than P — an analyst would agree — but not as a calibrated chance: Part A showed that a probability from a network need not mean what it says.

B.3 · One hour on 131 stations

Now every station, at the rates chosen in B.2. Reading the hour takes a minute or two the first time (most of it is the dense array); after that it comes from the disk. Picking the hour on all stations takes seconds on a GPU and a few minutes on a CPU.

read 131 stations (598 MB) in 33 s; no data in the hour: ['BG.FFA', 'NC.NLB']; 129 stations with data

picked in 81 s: 27,714 picks, {'P': 13993, 'S': 13721}

Each dot below is one pick with a probability of at least 0.5, placed at its time and at its station's row; rows follow A–A′, and the shaded rows are dense-array nodes. An earthquake appears as a column: it reaches every station within a few seconds, so its picks stand almost on one vertical line. The three boxes are drawn again underneath, stretched out in time.

Three things stand out.

• The burst. For the first half hour there are a few columns; from 23:59 on there is hardly a gap. The catalogue lists 67 earthquakes in this hour, and the top panel suggests that PhaseNet sees many more. B.5 counts them.
• The moveout. In each zoom the P picks (blue) form a curve: the stations nearest the source pick first, and the farther ones, at both ends of A–A′, later. The S picks (red) follow the same shape, later and spread wider, because S travels at a little more than half the speed of P.
• Picks that belong to no column. Isolated dots, some on a few rows that pick much more often than their neighbours. They are either very small earthquakes seen by one or two sensitive stations, or noise. B.6 and B.7 look at them.

B.4 · Does PhaseNet agree with the analysts?

The catalogue's picks are on the permanent stations only (the dense array was never part of the routine network). For each of them we look for a PhaseNet pick of the same phase on the same station within 0.5 s. The share that finds one is the recall: the fraction of the analysts' picks that PhaseNet also made. The time difference, PhaseNet minus catalogue, is the residual. We report its median, which says whether PhaseNet picks systematically early or late, and its median absolute deviation (MAD), which says how far a typical pick strays; both ignore the few matches far out in the tails, which the 0.5 s window lets in.

phase     picks   n  recall  median_ms  mad_ms  mean_ms  std_ms
    P automatic 627    0.83       -4.0     6.0    -4.91   48.24
    P    manual 657    0.90        0.0    10.0     3.50   61.11
    S    manual  60    0.97       17.0    26.0    30.52   86.19

PhaseNet made 9,921 picks on these stations; the catalogue 1,344

PhaseNet finds most of the analysts' picks, and where it does, a typical P pick lies within about a hundredth of a second of theirs (the MAD). Three details are worth a look.

• PhaseNet agrees more closely with the network's automatic P picks than with the manual ones, and sits a few milliseconds earlier than the automatic picker.
• Against the manual P picks the histogram has a tail to the right: on some stations PhaseNet picks later than the analyst — perhaps where the first motion is weak and the analyst caught it while the network did not.
• S is harder: the residuals are wider and PhaseNet picks S slightly later than the analysts did. There are only 61 manual S picks in the hour, so the S statistics are thin.

How does this compare with the paper? On its test set of NCEDC records, Zhu & Beroza (2019, Table 1) report, for residuals below 0.5 s, a mean of 2 ms and a standard deviation of 52 ms for P, and 3 ms and 83 ms for S. The mean_ms and std_ms columns above are the same two statistics for this hour, and they are close to the paper's: PhaseNet does about as well on the Geysers hour as on its own test set. The standard deviation is much larger than the MAD because it is pulled out by the few residuals far in the tails.

The last line of the printout is the larger story. PhaseNet made many times more picks on these stations than the catalogue holds. Some of those belong to earthquakes the catalogue does contain but at stations the analysts did not use; many belong to earthquakes the catalogue does not contain; some are false. The analysts' picks cannot tell those apart — they only test recall. B.5 uses the dense array to sort them.

B.5 · How many earthquakes does the catalogue miss?

A pick on one station could be anything. Picks on many stations at once mean an earthquake, because noise at stations hundreds of metres apart does not line up in time. The dense array lets us use this directly. An earthquake under the array reaches its nodes within about a second (the zooms of B.3), so we call it a detection when at least 15 nodes make a P pick (probability ≥ 0.5) within one second of each other. A run of such moments with no pause of 2 s or more counts as one detection.

This is the simplest form of association, the step that groups picks into earthquakes; Session 2, Location & relocation, located earthquakes from picks that were already grouped. It says nothing yet about where each earthquake is; it only counts them.

a detection: at least 15 nodes with P probability ≥ 0.5 within 1 s; a pause of 2 s ends it
matched when the detection is within 1.5 s of origin time + 0.7 s
135 detections in the hour; 60 match a catalogue earthquake, 75 do not
of the 67 catalogue earthquakes, 61 were detected; 8 of the 67 have no magnitude in the catalogue
unmatched detections: 77% come within a minute after a catalogue earthquake, 96% after the first M ≥ 3
catalogue earthquakes not detected (magnitude, distance from the array centre):
                       time  mag   km
2022-10-19T23:59:46.850000Z 1.43 2.14
2022-10-20T00:00:56.190000Z  NaN 0.44
2022-10-20T00:04:14.210000Z  NaN 1.07
2022-10-20T00:04:26.660000Z  NaN 1.29
2022-10-20T00:16:04.240000Z 0.85 0.89
2022-10-20T00:18:52.390000Z 0.75 1.00

The figure counts both as the hour goes on: the catalogue's earthquakes and the array's detections.

Before the burst, the array and the catalogue both count only a few earthquakes. From the first M 3.3 (dashed lines mark the three M ≥ 3) the array's count pulls away. Most of the unmatched detections come during the burst, within a minute of a catalogued earthquake (printed above) — the time when earthquakes follow each other within seconds and their records overlap.

Are the extra detections real earthquakes? Below are four of them, picked at random: the vertical records of every third dense-array node, in A–A′ order, with PhaseNet's picks. Each trace is scaled by its own noise level before the event, so that what rises above the noise is visible.

In each panel, P picks (blue) line up across many nodes and S picks (red) follow a few tenths of a second later on the same nodes — the signature of a small earthquake close beneath the array. In the first panel the arrivals are weak and only part of the array picks them; in the other three they are clear on most nodes. These four, at least, are earthquakes the catalogue does not contain.

Two cautions before counting all the unmatched detections as missed earthquakes. The matching uses the origin times only, so a catalogued earthquake whose P reaches the array unusually late, or whose S triggers a second detection after a pause, can be counted twice. And the catalogue is the routine network's — tens of stations, BG boreholes among them — which has to decide what is an earthquake with far fewer stations close to each source than the array has.

B.6 · Which sampling rate?

PhaseNet was trained on records sampled at 100 Hz, but it has no notion of seconds: it sees a sequence of numbers. Give it a record sampled at 500 Hz and every wave is spread over five times as many samples as at 100 Hz — a 25 Hz wave, 4 samples per cycle at 100 Hz, takes 20 samples, which to the network looks like a 5 Hz wave did in training. The usual practice is therefore to resample everything to 100 Hz. That has a cost: to resample to 100 Hz without aliasing, everything above about 40 Hz has to be filtered out first.

Does that cost matter here? The Geysers earthquakes are small, and small earthquakes radiate their energy at high frequencies. Below is the spectrum of the P wave of the M 1.3 from B.2 at the BG borehole station SQK, recorded at 500 Hz, next to the noise just before it.

at SQK, S arrives 0.39 s after P; the P window is 0.37 s long

SQK, M 1.34: 17% of the P-wave energy above 2 Hz lies above 40 Hz, where the signal is 11 times the noise

A good part of this small earthquake's P-wave energy lies above 40 Hz, well above the noise (the share is printed above). A picker given the 100 Hz version of this record never sees it. This is one earthquake at one station, but small earthquakes in general radiate more of their energy at high frequencies than large ones do. So should the BG boreholes keep the 500 Hz they were recorded at? We run them three ways: resampled to 100 Hz, resampled to 250 Hz, and left at 500 Hz.

We need two measures, because a picker can fail in two ways: it can miss arrivals, and it can invent them.

1. Recall of the analysts' picks (B.4): how many of the arrivals we know about does it find?
2. Share of picks that follow an array detection. A real arrival at a station comes from an earthquake; most earthquakes in this hour were detected by the dense array (B.5), and their waves reach a station within a few seconds of that detection (at most the station's distance divided by 4 km/s, plus 3 s for S). A pick that follows no detection is either from an earthquake the array did not see or is false. Random times give the share expected by chance.

 rate  recall_P  recall_S  picks  n_follow  n_none  follow  chance
  100      0.77      0.96   4936      4399     537    0.89    0.16
  250      0.93      1.00   9027      6418    2609    0.71    0.15
  500      0.90      1.00   9630      5741    3889    0.60    0.15

The figure puts the two measures side by side: recall of the analysts' P picks on the left, the share of picks that follow an array detection on the right (the dotted line is what random times would give).

The rate changes two things the network sees: how many samples a wave is spread over, and how many seconds the 3,000-sample standardisation window spans. (The pick rule is set in seconds, so it is the same at every rate.) The table and the figure show what that does.

From 100 to 250 Hz, recall rises sharply. At 100 Hz the boreholes miss many more of the analysts' picks: the energy above 40 Hz was part of what made those arrivals pickable. The picks that follow a detection (n_follow) increase — more arrivals of real earthquakes are found — and so do the picks that follow none (n_none), so the share that follows a detection falls.

From 250 to 500 Hz, recall does not improve. Compare n_follow and n_none between the two rows: going to 500 Hz buys no more of the picks that follow a detection, while the picks that follow none keep growing. At 500 Hz every wave is spread over twice as many samples again, and short bursts of noise begin to look to the network like the arrivals it learned.

This is the trade-off between recall and precision that Part A met with the probability threshold, now set by a different knob. The rate is a hyperparameter, and choosing it took two measures: on recall alone, 250 and 500 Hz look alike; on the share of picks after a detection alone, 100 Hz looks best. Together they point to 250 Hz for the boreholes, which is what B.3 used — a judgement, since how many false picks are worth how many found arrivals depends on what the picks are for. The dense array keeps its recorded 200 Hz and the other stations their 100 Hz. The general lesson: a test against labels measures what the labels contain — here, recall — and says nothing about the picks they do not contain.

Try it: raise the NC stations, recorded at 100 Hz, to 200 Hz with pick_station(row, 200.0). Raising the rate adds no high frequencies — there are none to add — but spreads every wave over twice as many samples. Measure both quantities, and decide whether you would keep it.

B.7 · Where does the picker go wrong?

Every station is run the same way, but stations are not alike. We count, for each station, its picks and the share of them that follow an array detection (B.6), and list the stations at both ends.

stations with fewer than 20 picks in the hour (km: distance from the dense array):
net   sta          kind   km  picks
 NC   GSS  short-period 22.4    0.0
CEC   S45      CEC node  0.4   18.0
 NC   GSG  short-period  9.8    0.0
 NC  GWKB     broadband 37.4   11.0
 NP  1744 strong-motion 23.6   10.0
 BG   FFA            BG  8.6    0.0
 NC   GSX     broadband 24.8   12.0
 NC   NLB  short-period 23.7    0.0
 CE 68208 strong-motion 36.4   19.0
 NC  NMTB     broadband 31.2   11.0

the five stations with the most picks that follow no detection:
net sta     kind  picks  follow
 BG BUC       BG  623.0    0.26
CEC S38 CEC node  398.0    0.59
 BG SB4       BG  413.0    0.62
 BG MCL       BG  366.0    0.62
CEC S09 CEC node  362.0    0.65

The figure shows every station's picks in the hour, in A–A′ order: the part that follows an array detection and the part that does not.

Three kinds of failure show up.

• Stations without signal. A few stations give no picks at all. For two of them the archive holds files for the day but no data for this hour (B.3). Others are tens of kilometres from the array (the km column), where the small earthquakes of the hour are too weak to record. One, NC.GSG, records nothing above its noise even for the M 3.2 at about 10 km: a dead channel. A picker cannot tell a dead station from a quiet one; a count of picks per station, read with its distance, can.
• Stations that pick too often. A handful of stations — some BG boreholes, some dense-array nodes — make many picks that follow no detection. Some may be real earthquakes too small or too far for the array to detect; on a station that picks several times a minute through the quiet first half hour (the BUC row in the B.3 figure), most must be noise. A practical pipeline checks each station's pick rate and the share of its picks that associate with earthquakes, and looks at the records of the outliers.
• S on one component. Most NC short-period stations record only the vertical component. S arrives mostly on the horizontals, so on these stations PhaseNet sees one third of its usual input (the missing components are zeros, and the standardisation scales the live one up to make up for them). Compare their S picks with the three-component stations:

                P      S  S_per_P
components                       
1           318.0  273.0     0.86
3           131.0  112.0     0.85
the analysts' S picks on these stations, by number of components: {3: 10, 1: 5}
components
1    0.86
3    0.87

PhaseNet makes about as many S picks per P pick on the vertical-only stations as on the three-component ones. Whether those S picks are right, the analysts can hardly tell us: they made only five S picks on these stations in the hour. One check needs no labels: an S pick should come a plausible time after a P pick on the same station (here between 0.1 and 8 s; Part A used the same physics). The last lines give that share for the two kinds of station, and it is the same. By this check the S picks on vertical-only stations look as plausible as the others — but the check cannot tell a correctly timed S from one that is a few tenths of a second off. Part A showed that the catalogue's S labels almost vanish on vertical-only stations, and PhaseNet was trained on that catalogue, so it has probably seen few examples of S on a vertical component alone. An S pick from such a station deserves less trust than its probability suggests.

What we found

The questions at the top, answered.

1. The data. One hour on 129 stations with data: 55 permanent stations within 40 km (BG boreholes at 500 Hz, other networks at 100 Hz), read by byte range from the NCEDC archive, and 74 dense-array nodes at 200 Hz from Google Cloud Storage.
2. The picker. PhaseNet is a U-Net that turns three components into P and S probabilities; on a continuous record it needs a running standardisation and is applied to the whole record at once.
3. One hour. Earthquakes appear as columns of picks across the stations, with P and S moveout along A–A′, and a quiet half-hour turns into a burst.
4. Agreement with the analysts. PhaseNet finds most of the analysts' picks (a share of 0.83 to 0.97) and agrees with them to about 10 ms for P; S is wider.
5. The catalogue. The dense array detects about twice as many earthquakes as the catalogue holds in this hour; the four we looked at are clear earthquakes, and most of the extra detections come during the burst.
6. The sampling rate. Lowering the boreholes' rate to 100 Hz loses real high-frequency arrivals; keeping their full 500 Hz adds picks that follow no detection without finding more. We chose 250 Hz with two measures, recall and the share of picks that follow a detection.
7. Failures. Stations without signal, a few noisy stations, and S picks on vertical-only stations that the analysts' picks cannot check.

Takeaways

Machine learning

• A network trained on one sampling rate sees samples, not seconds, so the rate we give it at inference is a hyperparameter that changes what it sees.
• A test against labels measures recall only; false picks need a second test that does not rely on the labels, such as agreement across many stations.
• A fully convolutional network can be applied to a record of any length in one pass, with a running standardisation standing in for the per-window standardisation it was trained with.

Seismology

• A deep-learning picker run on continuous data reproduces the analysts' P picks to about a hundredth of a second and makes far more S picks than a routine catalogue contains.
• In this hour, about half of the earthquakes the dense array detected are not in the routine catalogue, most of them during the burst, when records of successive earthquakes overlap.
• Small, close earthquakes can carry a sizeable share of their energy above 40 Hz, so data recorded at high rates should not be resampled to 100 Hz before picking without checking what is lost.

References

• Zhu, W., & Beroza, G. C. (2019). PhaseNet: a deep-neural-network-based seismic arrival-time picking method. Geophysical Journal International, 216(1), 261. https://doi.org/10.1093/gji/ggy423
• NCEDC (2014). Northern California Earthquake Data Center. UC Berkeley Seismological Laboratory. Dataset. https://doi.org/10.7932/NCEDC