Seismic Phase Association Based on the Maximum Likelihood Method

Phase association is a process that links seismic phases triggered at the stations of a seismic network to declare the occurrence of earthquakes. During phase association, a set of phases from different stations is examined to determine the common origin of phases within a specific region, predominantly on the basis of a grid search and the sum of observations. The association of seismic phases in local earthquake monitoring systems or earthquake early warning systems is often disturbed not only by transient noises, but also by large regional or teleseismic events. To mitigate this disturbance, we developed a seismic phase association method, binder_max, which uses the maximum likelihood method to associate seismic phases. The method is based on the framework of binder_ew, the phase associator of Earthworm, but it uses a likelihood distribution of the arrival information instead of stacking arrival information. Applying binder_max to data from seismic networks of South Korea and Ohio, United States, we found a significant improvement in the robustness of the method against large regional or teleseismic events compared to binder_ew. Our results indicate that binder_max can associate seismic phases of local earthquakes to the same degree as binder_ew as well as can avoid many of the false associations that have limited binder_ew.


INTRODUCTION
When seismic phases are triggered at the stations belonging to a seismic network, they are examined to determine whether a set of phases originated from a common location. Although the concept of association is rather simple, seismic phase association is an intricate task in real-time automated earthquake monitoring. This is because there can be various origins of seismic signals, such as small to large, shallow to deep, or near to distant earthquakes. In addition, various transient noises can be falsely detected as seismic phases and may be misidentified. In real-time earthquake monitoring systems, phase association is predominantly based on a grid search algorithm (e.g., Earthworm, Antelope, SeisComp3, and GLASS3). Detected arrivals are back-projected and stacked for each grid that the region of interest was discretized into. Then, the association stack is examined to investigate whether the count of back-projected arrivals on each grid exceeds a certain threshold (Johnson et al., 1997;Yeck et al., 2019).
This approach has also been used for locating earthquakes (Zhou, 1994;Font et al., 2004;Theunissen et al., 2012). The equal differential time (EDT) surface is defined as the collection of spatial grids that satisfy the time difference tolerance between P-wave arrivals at pairs of stations. The spatial grid traversed by the most EDT surfaces represents the hypocenter of the earthquake. Within this context, Lomax (2005) used a Gaussian distribution of the EDT surface to investigate earthquake locations. The search of the maximum likelihood together with the Gaussian distribution simply becomes a classical least squares problem and is subject to strong bias in the presence of outliers (Tarantola and Valette, 1982a;Tarantola and Valette, 1982b). Therefore, Lomax (2005) summed a Gaussian distribution of the EDT surface rather than multiplied, which leads to a highly complex EDT-likelihood function (Lomax, 2005). Tarantola and Valette (1982a) showed that the use of a heavytailed probability function can manage the inclusion of outliers in constructing the likelihood function through the multiplication of probability functions. Consequently, Sheen (2015) introduced a robust maximum-likelihood earthquake location method (MAXEL) that employs the Student's t-distribution, which decreases more slowly than the Gaussian function, and constructs the likelihood function by multiplying the Student's t-distribution of the EDT surfaces to locate an earthquake. Results from many previous studies indicate that the MAXEL method demonstrates good results for locating events from a small number of triggers, even with the presence of outliers (Sheen, 2015;Sheen et al., 2016;Chen et al., 2019).
The MAXEL method was used on the Earthworm platform to develop eqmaxel, a fully functional Earthworm module (Lisowski et al., 2017). This module uses a list of phases obtained from eqassemble, which gathers pick information using binder_ew, an association module of Earthworm. Results from previous studies indicate that eqmaxel can filter out both regional and teleseismic events that otherwise would have been falsely located inside the network using the traditional Earthworm hypoinverse solution (Chen et al., 2019;Lisowski et al., 2017). However, it was found that eqmaxel can still be affected by outliers that were falsely associated by binder_ew (Lisowski et al., 2017). After an event is falsely initiated from a small number of phases including an outlier, the associator repeatedly places phases that match with a false location into the list of phases that are used in eqmaxel. Accordingly, we modified the original association module in Earthworm, binder_ew, into a new binder, binder_max, based on the maximum likelihood location algorithm of MAXEL.

SEISMIC PHASE ASSOCIATION
The phase associator of Earthworm, binder_ew, was originally designed by Johnson et al. (1997) and further developed by several other researchers. Since its early development, a large number of hyperparameters have been introduced to solve ad hoc problems which required careful tuning and reviewing of the results. Therefore, the level of tuning could affect the performance of the associator. In this study, for the sake of simplicity, only picks triggered on the vertical component were used and considered as P-arrivals. Notably, binder_ew also nucleates new events based on P-arrivals and uses S-arrivals only in a later step, namely, the "update procedure." Both associators obtain picks from a shared memory, PICK_RING, and attempts to identify the occurrence of the smallest number of earthquakes from a list of picked times. The fundamental procedure for the association of phases consists of two steps. When a new pick is passed into the ring, possible associations with active events are investigated. If this fails, a grid search is conducted. Picks from the list that are not associated with an event are back-projected onto each grid. The back-projection of binder_ew is similar to that utilized in the EDT approach proposed by Zhou (1994). For a perfectly determined hypocenter, the difference between the arrival times t i and t j picked at two stations i and j should be equal to the difference between the computed travel times T i and T j from a theoretical hypocenter to each station. Because of picking errors and random noises, such as velocity uncertainty in the subsurface, the hypocenter would be on a thick spatial surface where the residual Δ between the differential arrival and travel times is less than the tolerance error. The residual Δ ij for two stations i and j is calculated for each grid as follows: A set of EDT surfaces can be constructed relative to the latest pick t 0 . If the number of picks with a residual smaller than the selected threshold is greater than or equal to four, a new event is initiated and the picks are associated with that event. Then, the L1 locator in binder_ew is used to locate the initial hypocenter.
In this study, the concept of MAXEL was introduced into binder_ew for the initiation of an event. The MAXEL method is based on maximum-likelihood estimation with the EDTs of P-arrivals and Student's t distribution. We briefly introduce the MAXEL method; however, refer to Sheen (2015) for further details.
The likelihood distribution can be defined as the sum of a set of N likelihood functions generated by an N-1 independent probability function from N arrivals (Sheen, 2015). Instead of simply stacking back-projected picks and counting the number of picks, we calculated the likelihood distribution at each grid as follows: where P is the probability density function represented by Student's t-distribution for the EDT and the function of the residual in Eq. 1. Student's t-distribution has a probability density function given by in which Γ is the gamma function. The parameter ] is the number of degrees of freedom and set to be one in this study (Sheen, 2015). In addition, we used a jackknife resampling technique to discriminate outliers from the picks (Sheen, 2015). If the normalized likelihood value at the maximum location is greater than 0.4 (Sheen et al., 2016), the location represents the hypocenter and the origin time of the event can be estimated. Then, the residual between the pick time at station i and the theoretical arrival time, which is the sum of the estimated origin time and the travel time from the grid of a hypocenter candidate to the station, can be calculated. If the number of picks with a residual smaller than the threshold is greater than or equal to four, a new event is initiated and the L1 locator, which is identical to that of binder_ew, locates the initial hypocenter.
The initiated event may draw more picks from the ring or may be removed during the update procedure in the associator. As new picks are triggered and sent to PICK_RING, some are associated with a previously initiated event. When the residual of a new pick for an existing event is less than the threshold, the pick is associated with the event, and the location can be updated by the L1 locator. The event can be removed if the root mean square error estimated from the L1 locator is too large, which allows the associated picks to be released back to the ring.

EARTHWORM IMPLEMENTATION
Earthworm (Johnson et al., 1995) is an automatic earthquake monitoring system, which has modules for data acquisition, transport, and processing. The modules are independent programs that can communicate with each other using share memory regions, called message rings, which are a first-in firstout data structure. A message ring can contain seismic waveforms and information about picks and events. The configuration of the Earthworm platform can vary depending on what modules and rings are in use and can be easily customized (Olivieri and FIGURE 1 | Distribution of earthquakes and seismic stations (cyan triangles) used in this study. Circles represent earthquakes greater than or equal to magnitude 2.0, as reported by the Korea Meteorological Administration (KMA). Empty and red circles represent the events that were and were not associated by binder_max, respectively, while yellow squares indicate the epicenters of earthquakes associated by binder_max, which are linked with the epicenters located by the KMA. The blue solid line indicates the grid search area.
FIGURE 2 | Epicenters of regional earthquakes and associated results. Epicenters of the events associated by (A) binder_ew (circles) and (B) binder_max (squares). Red stars represent the epicenters located by the U.S. Geological Survey's National Earthquake Information Center (NEIC) or the Japanese National Research Institute for Earth Science and Disaster Prevention (NIED). Numbers in symbols represent the number of phases used for associations. An orange square indicates the events that were removed during the update procedure of the association.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 699281 Clinton, 2012). The configurations of the platform are depicted in Figure 1 of Quintiliani and Pintore (2013), Figure 3 of Chen et al. (2015), and Figure 2 of Bono et al. (2021). For a basic operation, data acquisition, picker, associator, and location modules would be required. The data acquisition module obtains seismic data from seismic instruments, other earthquake monitoring systems, or recorded waveforms, and inputs the data into a shared memory region, WAVE_RING. After retrieving the waveform from WAVE_RING, the picker module determines arrival times and then, feeds the pick information to another shared memory region named PICK_RING. The associator module pulls the pick information out from PICK_RING, and transports the event information to third shared memory region named HYPO_RING.
The Earthworm module tankplayer can play back historical events by sending the waveform to WAVE_RING; the module putpick can simulate the picker module by feeding the pick information to PICK_RING. The new association module binder_max will be distributed with the Earthworm platform in the next release in 2021. Therefore, users can run the earthquake monitoring procedure multiple times with the same waveform or the same picks. More details can be found at http://www.earthwormcentral.org.

Application to South Korea
The earthquake early warning system of South Korea began operation in 2015 . Since then, the system has successfully issued public warnings for earthquakes greater than magnitude 5.0. However, the system occasionally struggles with regional and teleseismic events, although no false warnings have been issued to date. For example, the May 30, 2015, Bonin Islands Mw 7.9 Earthquake that occurred at a depth of 664 km was split and falsely identified as four events by the system .
To demonstrate the capability of binder_max, we used picks from a pilot monitoring system, which uses the Earthworm platform of the Korea Institute of Geoscience and Mineral Resources (KIGAM), from March to December 2020. Data from 103 seismic stations were used in this study (Figure 1). The KMA uses seismic data from over 300 contributing stations for monitoring local earthquakes, whereas our pilot system used only a subset of the stations in South Korea. For the 10-months operational period of the pilot system, KMA cataloged 60 events with a magnitude greater than or equal to 2.0 but only 50 events were associated with binder_max. Note that 10 unassociated events occurred outside the seismic network. Moreover, binder_ew associated 50 events; however, two of the 10 unassociated events of binder_ew, which also occurred outside the network, differed from those of binder_max. However, both methods successfully associated and located all earthquakes that occurred within the network. Notably, the grid size for binder_max was 75 × 86 × 2, while that of binder_ew was 201 × 230 × 10 for the region of interest (Figure 1).
During this process, several regional and teleseismic earthquakes were misidentified by both associators (Supplementary Table S1). We found that binder_ew associated 32 events that were not in the catalog of KMA with 15 or more phases, and generated several split events with similar origin times, whereas binder_max only associated 11 events for the same conditions. Two out of 32 events were associated by both associators but were not listed in any worldwide earthquake catalogs. Visual inspection of seismograms of the two events revealed that these were small local events that occurred outside the network.
Results from our investigation of the association of binder_ew and binder_max for regional events occurring outside the search grid indicate that events were poorly located by both associators, owing to the search grid being too small to identify true hypocenters ( Figure 2). However, our results clearly indicate binder_max to be more robust than binder_ew for regional events. Ten events were associated with binder_ew, while only six events were associated with binder_max, of which one was removed during the update procedure. In addition, it is important to note that the numbers of associated picks for each event related to binder_max were much smaller than those related to binder_ew, and remaining picks were not associated with any other events by binder_max ( Figure 2B).
Of the 32 events, 20 regional or teleseismic earthquakes were incorrectly split into several events by binder_ew (Figure 3). The first example is the April 18, 2020, Bonin Island Mw 6.6 Earthquake occurred at a depth of 453 km and triggered a large number of seismic stations in South Korea. Similarly, as shown in Sheen et al. (2017), 16 events were independently associated with binder_ew, and although 13 of them were removed during the update procedure, three events remained with up to 28 phases. In contrast, binder_max initiated only one event that was associated with seven phases. The second example is the August 21, 2020, Banda Sea Mw 6.9 Earthquake at a depth of 627 km. This earthquake was also split into four events by binder_ew, and one event remained with 16 phases. However, although the same picks were used for the association, binder_max did not initiate an event, which is the desired result for a local earthquake monitoring system. Therefore, our results indicate that binder_max is robust enough to successfully process picks triggered by regional and teleseismic earthquakes (Figure 3). Only two of the 20 regional or teleseismic earthquakes that were incorrectly split by binder_ew, were split into two events by binder_max, and only one of each remained. In addition, nine of the 20 earthquakes were not initiated by binder_max.

Application to Ohio, United States
The state of Ohio monitors earthquakes using a seismic network known as OhioSeis (Hansen and Ruff, 2003). This network also uses the Earthworm system for base processing but has issues with false triggers from regional events and large mining blasts that occur in adjoining states. These false triggers due to noises or anthropogenic activities were often incorrectly located by binder_ew inside the network.
Association results for local and teleseismic earthquakes were obtained from the OhioSeis (Figure 4). While the region of grid Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 699281 search is the same for both associators, the number of search grids for binder_max was lower by about a factor of 35 than that for binder_ew. Although the January 22, 2021, M 2.4 Earthquake was associated and well located by both methods, the May 24, 2017, M 3.4 Earthquake was successfully associated by binder_max but split into two events by binder_ew. The location of the results from binder_ew also differed from the location determined by the NEIC. The November 24, 2019, M 6.3 Alaska Earthquake was located inside the network by binder_ew but binder_max located this at the boundary of the search grid. Moreover, the January 11, 2021, M 5.6 Indonesia Earthquake was not initiated by binder_max whereas binder_ew located this inside the network. As demonstrated in the examples from South Korea (Figure 3), binder_max was found to be more prone to not associate picks triggered by regional or teleseismic events, which is also the case for the OhioSeis. When this occurs, the earthquakes would be located at the boundary of the region of the grid search.
FIGURE 3 | Epicenters of associated events for regional and teleseismic earthquakes including association results for (A) the April 18, 2020 Bonin Island Mw 6.6 Earthquake at a depth of 453 km and (B) the August 21, 2020 Banda Sea Mw 6.9 Earthquake at a depth of 627 km. Red stars represent the epicenters located by the U.S. Geological Survey's National Earthquake Information Center (NEIC), while circles and squares represent the epicenters of events associated by binder_ew, and binder_max, respectively. Numbers in symbols represent the number of phases for the association, while numbers in parentheses denote the number of associated events (cyan) and removed events (orange), respectively.
FIGURE 4 | Association results from OhioSeis. Picks triggered by two local and two teleseismic events were associated by binder_ew (circles) and binder_max (squares). Stars represent the epicenters of the earthquakes determined by the U.S. Geological Survey's National Earthquake Information Center (NEIC). The region of grid search is represented by the blue solid line. Green and yellow symbols represent the epicenters of local events, while red and orange symbols depict the results for the November 24, 2019, M 6.3 Alaska Earthquake and the January 11, 2021, M 5.6 Indonesia Earthquake, respectively.

DISCUSSION
As described above, binder_max, which uses the maximum likelihood method, can associate better than binder_ew, the original phase associator of the Earthworm platform. The association is too complicated for optimal tuning and features numerous configurable parameters that control how to solve the problem (Dietz, 2002). Therefore, the improvement observed when using binder_max could be attributed to several factors. However, an optimal tuning of the associator would be beyond the scope of this study, which aims at highlighting the benefits of the maximum likelihood method to associate seismic phases. As shown in Figure 4A, binder_ew initiated several events due to the picks triggered by the April 18, 2020, Bonin Island Mw 6.6 Earthquake. It is shown that binder_ew initiated one of them close to the southeast corner of the grid search area but removed the event after associating 24 picks during the update procedure. This event was associated with five picks by binder_ew but also checked to be initiated with the same earlier four picks by binder_max. To illustrate the difference between the implementations of the concept of the EDT surface, the association grids of binder_max and binder_ew for this case are presented in Figure 5. On the one hand, binder_max calculates the likelihood values at each grid and determines the location of an event at the grid with the maximum likelihood value. On the other hand, binder_ew counts the number of the EDT surface that intersects a grid, and locates an event at the center of the grids with the maximum counts because there can be numerous grids intersected by the maximum number of EDT surfaces as shown in Figure 5. The true hypocenter of these picks is outside the grid search area in the southeast direction; thus, the likelihood values increase towards the hypocenter. However, the maximum count of intersected EDT surfaces is limited up to the number of picks and the grids with the maximum are elongated towards the true hypocenter, which makes it difficult to determine whether the hypocenter is inside the search area or outside it. Considering this difference between binder_max and binder_ew, we added a procedure to check if the initial location of an event is on the edge of the grid boundary. If so, the event is discarded.
It is evident that the computational cost for binder_max at each grid is more expensive than for binder_ew. However, as shown in Figure 5A, smooth likelihood distribution enables the use of a coarse grid for association and accordingly, binder_max can correctly associate picks with coarser grids than binder_ew. Although we do not provide the details here, the computation time for the association of binder_max, was less than a second, which would not decrease the performance of earthquake monitoring systems.
Deep learning-based automatic seismic phase pickers have shown remarkable improvement on the processing speed and the accuracy in both picking P and S phases, even on the cases with low SNR (e.g., Ross et al., 2018;Zhou and Beroza, 2019;Mousavi et al., 2020;Liao et al., 2021). Therefore, deep learning pickers can be used for routine seismic data processing and cataloging of events with local to regional seismic network (Baker et al., 2021;Walter et al., 2021). Although deep learning approach to automatic seismic phase association (Ross et al., 2019) was introduced, seismic phase association is still challenging (Mousavi et al., 2020;Baker et al., 2021;Walter et al., 2021). Based on this study, it is expected that binder_max will be extremely useful for real-time earthquake monitoring, earthquake early warning, and the association of large seismic picks generated by deep learning-based phase pickers.  Figure 3A are used for the association. Seismic stations triggered are represented by yellow triangles. Note that binder_max discarded this group of the picks because the location of the maximum likelihood lies on the edge of the grid boundary; in contrast, binder_ew initiated this as an event because the center of the grids with the maximum count lies inside the grid.

CONCLUSION
We developed a seismic phase associator, binder_max, based on the framework of the Earthworm phase associator binder_ew in conjunction with Sheen's (2015) maximum likelihood method, MAXEL. Instead of counting the number of picks satisfying a specific threshold at each grid, we used a likelihood distribution with Student's t distribution from the EDTs of P-arrivals, using a grid with a maximum of likelihood value to determine the initiation of an event.
Picks obtained from a pilot monitoring system of the KIGAM, South Korea were used to investigate the performance of binder_max. During the period from March to December 2020, all earthquakes with a magnitude greater than or equal to 2.0 occurring within the seismic network were associated and well located by both associators. However, our results clearly indicate that binder_max effectively processed picks from regional or teleseismic earthquakes, which were falsely associated and split into several events by binder_ew. The advantages of using binder_max are confirmed by the association results from the Ohio seismic network.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
D-HS implemented and applied the method and wrote the first draft of the manuscript. D-HS and PF organized the data. PF revised the manuscript.