  The association of seismic wave arrivals with causative earthquakes becomes
progressively more challenging as arrival detection methods become more
sensitive, and particularly when earthquake rates are high. For instance,
seismic waves arriving across a monitoring network from several sources may
overlap in time, false arrivals may be detected, and some arrivals may be of
unknown phase (e.g., P- or S-waves). We propose an automated method to
associate arrivals with earthquake sources and obtain source locations
applicable to such situations. To do so we use a pattern detection metric based
on the principle of backprojection to reveal candidate sources, followed by
graph-theory-based clustering and an integer linear optimization routine to
associate arrivals with the minimum number of sources necessary to explain the
data. This method solves for all sources and phase assignments simultaneously,
rather than in a sequential greedy procedure as is common in other association
routines. We demonstrate our method on both synthetic and real data from the
Integrated Plate Boundary Observatory Chile (IPOC) seismic network of northern
Chile. For the synthetic tests we report results for cases with varying
complexity, including rates of 500 earthquakes/day and 500 false
arrivals/station/day, for which we measure true positive detection accuracy of
> 95%. For the real data we develop a new catalog between January 1, 2010 -
December 31, 2017 containing 817,548 earthquakes, with detection rates on
average 279 earthquakes/day, and a magnitude-of-completion of ~M1.8. A subset
of detections are identified as sources related to quarry and industrial site
activity, and we also detect thousands of foreshocks and aftershocks of the
April 1, 2014 Mw 8.2 Iquique earthquake. During the highest rates of aftershock
activity, > 600 earthquakes/day are detected in the vicinity of the Iquique
earthquake rupture zone.
