Relocate with HypoDD¶

HypoDD relocates earthquakes with the double-difference method (Waldhauser and Ellsworth, 2000). It inverts the differences of the travel times of event pairs at common stations. Errors of the velocity model largely cancel for close events, so HypoDD sharpens the relative locations of a cluster without station corrections.

qseek export hypodd writes a HypoDD project folder from a run: the picks, the stations, the velocity model and the control files for ph2dt and hypoDD, ready to run.

On the Campi Flegrei example, 20 May 2024, Qseek exports 424 detections with 5893 picks. ph2dt links 347 of them and hypoDD relocates 340 in about 2 s. On these 340 events, the median absolute double-difference residual of the catalog differential times falls from 73 ms at the Qseek locations to 54 ms at the HypoDD locations, and from 67 ms to 34 ms for P. The playground runs this export and the comparison for you.

Export a run¶

Export the detections of a finished run. Start the export in the directory of the search configuration, so that the velocity model file is found:

qseek export hypodd my-search/ my-search-hypodd/

Qseek writes these files:

The export selects the detections and picks:

• A pick needs a confidence of at least min_pick_confidence (default 0.3), and its residual to the modeled arrival must not exceed max_residual (default 1 s). The confidence is the pick weight in phase.dat, limited to 1: machine learning pickers give a probability from 0 to 1, the STA/LTA image function the peak of its image, which can exceed 1.
• A detection needs at least min_picks of these picks (default 6). max_rms drops detections with a larger residual RMS, min_distance_border detections close to the border of the search volume.
• station.dat lists the stations with exported picks.
• The travel times are the observed picks minus the origin time. Station corrections of the run are not applied; the double-difference method does not need them.

HypoDD needs unique integer event IDs and station labels of up to 7 characters. Qseek numbers the detections in time order and uses the station code as the label, or network and station code if the station code is not unique. hypodd_results.py maps the IDs in hypoDD.reloc back to the detections with event_ids.csv.

Run HypoDD¶

Build ph2dt and hypoDD from the HypoDD distribution (version 2.1), then run both in the project folder:

cd my-search-hypodd/
HYPODD_BIN=~/src/HypoDD/bin ./run.sh

HYPODD_BIN is the directory of the binaries; leave it out if they are in your PATH. ph2dt writes the catalog differential times dt.ct and the initial locations event.sel, hypoDD the relocations hypoDD.reloc and its log hypoDD.log.

Check the log before you use the relocations:

• Linked events: ph2dt lists the events it selected and the weakly linked events in ph2dt.log. Events without enough links to their neighbors are not relocated.
• Condition number: the LSQR solver should reach a condition number (CND in the iteration table) of about 40 to 80. Raise DAMP in hypoDD.inp if it is higher, lower it if it is lower. On Campi Flegrei, the default damping of 80 gives a CND of 40 to 49.
• Shifts: the mean shifts DX, DY, DZ should fall to the noise level of the data within the last iterations. The centroid shift OS should stay below the location uncertainty of the detections; HypoDD does not constrain the absolute position of a cluster well.

Results¶

After hypoDD, run.sh runs hypodd_results.py. It writes the relocated events with their Qseek detections to two files:

• hypodd_relocations.csv: one event per row, sorted by origin time.
• hypodd_relocations.yaml: the events as Pyrocko events, named by their origin time like the Qseek detections. Open them in Pyrocko Snuffler or load them with pyrocko.model.load_events. This file needs Pyrocko: run.sh runs the script with python3, set PYTHON to the Python of your Qseek installation, e.g. PYTHON=.venv/bin/python ./run.sh.

To load the CSV file in QGIS, add it as a delimited text layer with the geometry definition Well known text (WKT), the field WKT_geom and the CRS EPSG:4326.

After you change hypoDD.inp and run hypoDD by hand, convert the relocations again:

python3 hypodd_results.py

Cross-correlation¶

Waveform cross-correlation measures the differential times of close events more precisely than picks: for similar waveforms, to a fraction of a sample. With cross_correlation, the export correlates the waveforms of the detections and writes the differential times to dt.cc. hypoDD combines them with the catalog differential times of ph2dt (IDAT=3).

{
  "cross_correlation": {}
}

The export loads the waveforms with the waveform provider of the run, so start it in the directory of the search configuration. For each event pair, it correlates the P and S phases at the stations of both events:

• Event pairs: each event with up to max_neighbors nearest events closer than max_separation (default 20 events within 2 km). hypoDD skips pairs with events that ph2dt did not keep.
• Windows: window_p and window_s start before and end after the exported pick, or the modeled arrival at stations without a pick (modeled_arrivals). The P window ends before the S window starts, so at close stations it holds the P wave only. The window of the first event is the template; the window of the second event is longer by the maximum lag on both sides.
• Filter: one zero-phase Butterworth bandpass for all channels, bandpass (default 1 to 15 Hz), applied to the windows with a padding of three periods of the low corner. Channels with gaps in the padded windows are skipped.
• Correlation: the normalized correlation of the components of the phase, stacked: Z for P, the horizontals for S by default. The maximum is interpolated with a parabola to a fraction of a sample. Maxima at the maximum lag and below min_correlation (default 0.7) are rejected; the weight is the squared correlation coefficient.
• Differential times: the travel time differences of the matched windows, relative to the origin times in event.sel. The windows only select the waveforms, so a modeled arrival gives the same differential time as a pick. The origin time correction OTC in dt.cc is 0.

With cross_correlation, the default iterations follow Table 1 of the HypoDD user guide: 10 iterations with down-weighted cross-correlation data, so the catalog data restore the large-scale picture, then 15 iterations in which the cross-correlation data dominate for event pairs closer than 2 km, at last closer than 500 m. Check the RMSCC and CC columns of the iteration table in hypoDD.log: the residuals of the cross-correlation data should fall to a few milliseconds, while most data stay in use.

Choose the settings for your data. The defaults are a starting point for local seismicity recorded at about 100 Hz; at lower sampling rates, the windows hold fewer samples and the lag is less precise:

• The bandpass should hold the energy of the smallest events above the noise, below 90% of the Nyquist frequency.
• A window should hold the phase and its first oscillations, not the coda.
• The maximum lag must exceed the error of the arrival times of both events, but a large lag lets the correlation jump by a period of the dominant frequency.

On Campi Flegrei, the export correlates 1391 event pairs with 8458 differential times in about 30 s. On 340 common events, the median absolute double-difference residual of the cross-correlation times is 54 ms at the plain Qseek locations, 48 ms with station corrections (SSST), 47 ms after hypoDD with catalog data only and 21 ms after hypoDD with both data types. The cross-correlation residuals do not depend on the picks, so they also compare Qseek locations with HypoDD on independent data.

The export logs how many traces it dropped: without data covering the windows and the filter padding, e.g. at gaps or at the start and end of the archive, or with a Nyquist frequency below the low corner of the bandpass. If no event pair has min_observations differential times, the export warns and writes catalog differential times only.

The filtered waveforms of the events stay in a cache of cache_size (default 2 GB); events that do not fit are loaded again.

Velocity model¶

The export takes the 1D velocity model of the ray tracer of the P phase: the Pyrocko Cake or fast marching model, written as layers with their P velocity and Vp/Vs ratio (IMOD=1). A constant velocity becomes HypoDD's straight-ray model (IMOD=5). 3D models are not exported.

HypoDD needs layers of constant velocity. Gradient layers are split into layers of at most max_layer_thickness (default 500 m) down to the bottom of the search volume. Each layer gets the harmonic mean velocity of its depth range, which keeps the vertical travel time. HypoDD allows 30 layers; raise max_layer_thickness if the split model needs more, or use a model with fewer layers.

Depths in HypoDD are in kilometers below sea level, like the depths of Qseek. HypoDD places the top of the model at the elevation of each station, so the velocity of the first layer applies from the station down to the second layer. The top of the first layer is written as 1 km above sea level: for an event at the top of the first layer, hypoDD 2.1 reads the velocity at the source outside its velocity model, which can make the inversion fail with NaN.

Sea level is the top of the model in HypoDD:

• Shallow events: detections above sea level are set to 0 km depth. Events that hypoDD moves above sea level are air-quakes, also when they are below the stations. By default they stay at their depth of the previous iteration; remove_airquakes removes them instead. On Campi Flegrei, hypoDD relocates 340 events when it keeps the air-quakes and 323 when it removes them.
• Stations below sea level: in a layered model, hypoDD moves borehole and ocean-bottom stations below sea level up to 0 m elevation, and the export warns about them. Only the constant velocity model (IMOD=5) keeps them below sea level.

Settings¶

Change the selection and the parameters of ph2dt and hypoDD with a JSON file. Unknown fields are errors, so a misspelled setting does not fall back to its default:

qseek export hypodd my-search/ my-search-hypodd/ --config hypodd.json

{
  "min_picks": 6,
  "max_rms": null,
  "min_distance_border": 0.0,
  "min_pick_confidence": 0.3,
  "max_residual": 1.0,
  "max_layer_thickness": 500.0,
  "ph2dt": {
    "min_weight": 0.0,
    "max_distance": null,
    "max_separation": 5000.0,
    "max_neighbors": 20,
    "min_links": 8,
    "min_observations": 8,
    "max_observations": null
  },
  "hypodd": {
    "max_distance": null,
    "min_links": 8,
    "initial_locations": "catalog",
    "solver": "LSQR",
    "remove_airquakes": false,
    "iterations": [
      {
        "n_iterations": 5,
        "weight_p": 1.0,
        "weight_s": 0.5,
        "max_residual": null,
        "max_separation": null,
        "weight_cc_p": -999.0,
        "weight_cc_s": -999.0,
        "max_residual_cc": null,
        "max_separation_cc": null,
        "damping": 80.0
      },
      {
        "n_iterations": 5,
        "weight_p": 1.0,
        "weight_s": 0.5,
        "max_residual": 6.0,
        "max_separation": 4000.0,
        "weight_cc_p": -999.0,
        "weight_cc_s": -999.0,
        "max_residual_cc": null,
        "max_separation_cc": null,
        "damping": 80.0
      },
      {
        "n_iterations": 5,
        "weight_p": 1.0,
        "weight_s": 0.5,
        "max_residual": 4.0,
        "max_separation": 2000.0,
        "weight_cc_p": -999.0,
        "weight_cc_s": -999.0,
        "max_residual_cc": null,
        "max_separation_cc": null,
        "damping": 80.0
      }
    ]
  },
  "cross_correlation": null
}

Distances are in meters, as everywhere in Qseek; the export converts them to the kilometers of HypoDD. The ph2dt defaults follow the HypoDD user guide for a dense local network: event pairs up to 5 km apart and at least 8 differential times per pair. The three default iteration sets weight P twice as strongly as S, then remove outliers beyond 6 and 4 standard deviations and limit the pair separation to 4 and 2 km.

--force replaces an existing export directory only after the export succeeded. The control files are plain text with comments. You can also edit ph2dt.inp and hypoDD.inp in the project folder and run HypoDD again.

min_picks: PositiveInt = 6

max_rms: PositiveFloat | None = None

min_distance_border: float = 0.0

min_pick_confidence: float = 0.3

max_residual: PositiveFloat = 1.0

max_layer_thickness: PositiveFloat = 500.0

ph2dt: Ph2DTSettings

hypodd: HypoDDSettings

cross_correlation: CrossCorrelation | None = None

set_cc_iterations() -> None

correlate(search: Search, events: list[tuple[int, EventDetection, datetime]], event_picks: list[list[tuple[NSL, float, float, str]]], max_distance: float, stations: dict[NSL, tuple[float, float, float]]) -> dict[tuple[int, int], list[DifferentialTime]]

write_phases(file: Path, events: list[tuple[int, EventDetection, datetime]], event_picks: list[list[tuple[NSL, float, float, str]]], labels: dict[NSL, str]) -> dict[str, int]

get_velocity_model(search: Search) -> tuple[list[HypoDDLayer], int]

get_subclasses() -> tuple[type[Exporter], ...]

min_weight: float = 0.0

max_distance: PositiveFloat | None = None

max_separation: PositiveFloat = 5000.0

max_neighbors: PositiveInt = 20

min_links: PositiveInt = 8

min_observations: PositiveInt = 8

max_observations: PositiveInt | None = None

max_distance: PositiveFloat | None = None

min_links: PositiveInt = 8

initial_locations: Literal['catalog', 'centroid'] = 'catalog'

solver: Literal['LSQR', 'SVD'] = 'LSQR'

remove_airquakes: bool = False

iterations: list[IterationSet]

n_iterations: PositiveInt = 5

weight_p: float = 1.0

weight_s: float = 0.5

max_residual: PositiveFloat | None = None

max_separation: PositiveFloat | None = None

weight_cc_p: float = UNUSED

weight_cc_s: float = UNUSED

max_residual_cc: PositiveFloat | None = None

max_separation_cc: PositiveFloat | None = None

damping: PositiveFloat = 80.0

bandpass: tuple[PositiveFloat, PositiveFloat] = (1.0, 15.0)

window_p: PhaseWindow = PhaseWindow(seconds_before=0.1, seconds_after=0.5, max_lag=0.2, components='Z')

window_s: PhaseWindow = PhaseWindow(seconds_before=0.2, seconds_after=1.0, max_lag=0.3, components='NE12')

min_correlation: float = 0.7

max_separation: PositiveFloat = 2000.0

max_neighbors: PositiveInt = 20

min_observations: PositiveInt = 4

modeled_arrivals: bool = True

channels: list[str] | None = None

cache_size: ByteSize = ByteSize(2 * 1024 ** 3)

n_parallel: PositiveInt = 8

padding: float

select_pairs(events: list[CorrelationEvent]) -> list[tuple[int, int]]

correlate(events: list[CorrelationEvent], waveform_provider: WaveformProvider) -> dict[tuple[int, int], list[DifferentialTime]]

log_stats(stats: Counter[str], n_pairs: int) -> None

load_waveforms(event: CorrelationEvent, waveform_provider: WaveformProvider, stats: Counter[str] | None = None) -> list[Trace]

filter_waveforms(traces: list[Trace], spans: dict[NSL, tuple[float, float]], stats: Counter[str] | None = None) -> list[Trace]

correlate_pair(event_1: CorrelationEvent, event_2: CorrelationEvent, traces_1: list[Trace], traces_2: list[Trace]) -> list[DifferentialTime]

seconds_before: NonNegativeFloat

seconds_after: PositiveFloat

max_lag: PositiveFloat

components: str