GaMMA: Gaussian Mixture Model Associator

1. Install

pip install git+https://github.com/wayneweiqiang/GaMMA.git

The implementation is based on the Gaussian mixture models in scikit-learn

2. Related papers

• Zhu, W., McBrearty, I. W., Mousavi, S. M., Ellsworth, W. L., & Beroza, G. C. (2022). Earthquake phase association using a Bayesian Gaussian mixture model. Journal of Geophysical Research: Solid Earth, 127(5), e2021JB023249.

Generating automatic picks:

• Zhu, W., & Beroza, G. C. (2019). PhaseNet: a deep-neural-network-based seismic arrival-time picking method. Geophysical Journal International, 216(1), 261-273.

Adding 1D velocity model support:

• Xi, Z., Wei, S. S., Zhu, W., Beroza, G. C., Jie, Y., & Saloor, N. (2024). Deep learning for deep earthquakes: insights from OBS observations of the Tonga subduction zone. Geophysical Journal International, 238(2), 1073-1088.
• Jie, Y., Wei, S. S., Zhu, W., Freymueller, J. T., & Elliott, J. (2025). Background Seismicity and Aftershocks of the 2020-2021 Large Earthquakes at the Alaska Peninsula Revealed by a Deep-learning-based Catalog. Authorea Preprints.

3. Examples

• Inputs (the same names as ADLoc since 1.8):
  • picks: columns station_id, phase_time, phase_type ("P" or "S"), phase_score, phase_amplitude (m/s).
  • stations: columns station_id and the coordinates named in dims (default ["x_km", "y_km", "z_km"], z_km positive down).
  • The frame follows proj, as in ADLoc. With proj (a pyproj Proj), every station goes into that one frame and events also carry x_km, y_km in it. Without proj, one run can cover any area, up to the whole globe: picks cluster on the globe, each segment is projected around its own stations (azimuthal equidistant, as ADLoc does per event), and its events are located within 300 km of those stations. dims (the stations' own projected columns, as before 1.13) is deprecated but still works.
  • Events come back with the coordinates named in dims, plus depth_km. Names from before 1.8 (id, timestamp, type, prob, amp, x(km), min_picks_per_eq, ...) still work and are translated with one warning.
• Hyperparameters:
  • xlim_km, ylim_km, zlim_km: the region (km) where events can be: the locator's bounds, and zlim_km also sets the depth where each event starts. bfgs_bounds, when given, overrides the bounds.
  • use_amplitude (default = True): If using amplitude information.
  • vel (default = {"P": 6.0, "S": 6.0 / 1.73}): velocity for P and S phases (uniform velocity model).
  • eikonal (optional): When specified, the 1D velocity model overrides the uniform velocity model:

  {"vel": {"Z": [depths],
  "P": [p_velocities],
  "S": [s_velocities]},
  "h": grid_spacing,
  "zlim_km": z_range}  # optional: the tables' depth range; default the events' zlim_km

  The tables reach across the farthest event-station pair (the region and the stations together); xlim_km and ylim_km here, when given, set that range instead. Set the tables' zlim_km deeper than the events' (e.g. (0, 60), below the Moho) so far stations' first arrivals (Pn) are modelled while events stay within the region; a range that does not cover the events' zlim_km is widened to cover it, with a warning. Since 1.10.0 the 1D model is read at absolute depth, as in ADLoc and ADTomo: one travel-time table per station depth (stations within 100 m share one), so a station 2 km above sea level sees the same layers as one at sea level, 2 km further up. Event depths (z_km, depth_km) are below sea level; before 1.10.0 the model hung from each station, and depths were measured from the station surface (about the stations' mean elevation shallower).
  • segment (default = "dbscan"): how to cut a long sequence of picks into segments, which speeds up association a lot. "dbscan": the hierarchical DBSCAN cut below (the default again since 1.12.0). "link": groups of picks that one event could produce (stations within 100 km whose time differences fit the distance between them), each at most segment_max_picks (default = 200) picks, split further in time (the default from 1.7.0 to 1.11). None: all picks in one segment. It replaces use_dbscan, which still works with a warning (True → "dbscan", False → None).
  • dbscan_eps (default = 10.0s): The maximum time between two picks for one to be considered as a neighbor of the other. See details in DBSCAN
  • dbscan_min_samples (default = 3): The number of samples in a neighborhood for a point to be considered as a core point. See details in DBSCAN
  • oversampling_factor (default = 10): The initial number of clusters is determined by (Number of picks)/(Number of stations)/(Inital points) * (oversampling factor).
  • initial_points (default=[1,1,1] for (x, y, z) directions): Initial earthquake locations (cluster centers). For a large area over 10 degrees, more initial points are helpful, such as [2,2,1].
  • covariance_prior (default = 5 × (1.5², 0.5²) = (11.25, 1.25)) and degrees_of_freedom_prior (default = 5): the Wishart prior on each event's time (s²) and amplitude (log10²) residuals, as in scikit-learn: degrees_of_freedom_prior pseudo-picks with spreads of 1.5 s and 0.5. A larger time spread merges events; a smaller one splits them.
  • Filtering low quality association: each event is relocated on its own picks; picks with residuals above the limits below are dropped and the event relocated again, until its picks stop changing. The event is kept if enough picks remain.
    • min_picks: Minimum picks for associated earthquakes. We can also specify minimum P or S picks:
      • min_p_picks: Minimum P-picks for associated earthquakes.
      • min_s_picks: Minimum S-picks for associated earthquakes.
    • max_residual_time (default = 1.0): Max phase time residual (s)
    • max_residual_amplitude (default = 1.0): Max phase amplitude residual (log10 scale)
    • The old names max_sigma11 and max_sigma22 still work with a warning, as max_residual_time and max_residual_amplitude; max_sigma12 is ignored.
    • min_station_support (default = 0.3; 0 keeps every event): drops events that picks lined up by chance. A real event is picked by the stations nearest to it; a chance event is a few scattered picks with silent stations in between. Support = picking stations observed / expected among an event's nearest 8 reliable stations: about 1 for a real event, about 0 by chance. Each station's reliability is measured on the run's own large events (15+ picks), so stations that are often silent (offline, noisy, skipped by analysts) count little. 0.3 asks roughly 3 of the nearest 8 reliable stations to pick; events picked at more than 8 stations need proportionally less (0.15 at 16). Weak events (5–8 picks) are often dropped too, since silence cannot tell them from chance well; a larger value (e.g. 0.4) removes more false events and a few more real ones.

Note the association speed is controlled by dbscan_eps and oversampling_factor. Larger values are preferred, but at the expense of a slower association speed.

• Synthetic Example

• Example using PhaseNet picks

See details in the notebook: example_phasenet.ipynb

• Example using Seisbench

See details in the notebook: example_seisbench.ipynb

More examples can be found in the earthquake detection workflow -- QuakeFlow