{'fraction_non_alphanumeric': 0.03333455303889641, 'fraction_numerical': 0.014819422591386439, 'mean_word_length': 5.05718085106383, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Earthquake Phase Association using a Bayesian Gaussian Mixture Model

Earthquake phase association algorithms aggregate picked seismic phases from a network of seismometers into individual earthquakes and play an important role in earthquake monitoring. Dense seismic networks and improved phase picking methods produce massive earthquake phase data sets, particularly for earthquake swarms and aftershocks occurring closely in time and space, making phase association a challenging problem. We present a new association method, the Gaussian Mixture Model Association (GaMMA), that combines the Gaussian mixture model for phase measurements (both time and amplitude), with earthquake location, origin time, and magnitude estimation. We treat earthquake phase association as an unsupervised clustering problem in a probabilistic framework, where each earthquake corresponds to a cluster of P and S phases with hyperbolic moveout of arrival times and a decay of amplitude with distance. We use a multivariate Gaussian distribution to model the collection of phase picks for an event, the mean of which is given by the predicted arrival time and amplitude from the causative event. We carry out the pick assignment for each earthquake and determine earthquake parameters (i.e., earthquake location, origin time, and magnitude) under the maximum likelihood criterion using the Expectation-Maximization (EM) algorithm. The GaMMA method does not require the typical association steps of other algorithms, such as grid-search or supervised training. The results on both synthetic test and the 2019 Ridgecrest earthquake sequence show that GaMMA effectively associates phases from a temporally and spatially dense earthquake sequence while producing useful estimates of earthquake location and magnitude.


Introduction
Earthquake catalogs are fundamental products that are widely used in seismology to study and model various aspects of seismicity. Extensive efforts have been made to generate more complete catalogs with many more smaller earthquakes and more precise location and magnitude estimates. These high-resolution high-precision catalogs have the potential to reveal relationships among earthquakes and illuminate active structures that would otherwise remain hidden (Waldhauser & Schaff, 2008;Hauksson et al., 2012;Yoon et al., 2015;Ross, Trugman, et al., 2019;Park et al., 2020;Tan et al., 2021;Beroza et al., 2021). A standard earthquake-monitoring workflow from seismic waveforms to earthquake catalogs includes several tasks, including earthquake detection, phase picking (Allen, 1978), phase association (Yeck et al., 2019), earthquake location (Klein, 2002), and magnitude estimation (Richter, 1935).
The phase picking step detects seismic phases such as P-wave and S-wave phases at each seismic station. The phase association step aggregates these phases from multiple stations of a seismic network into separate groups associated with each earthquake. Earthquake location and magnitude are then estimated from the associated phase information, i.e., arrival time and amplitude. The resulting catalog can be further enhanced through template matching (Gibbons & Ringdal, 2006;Shelly et al., 2007;Peng & Zhao, 2009) or subspace projection (Harris & Dodge, 2011;Barrett & Beroza, 2014) that use the detected earthquakes as templates to re-scan the waveforms and detect small earthquakes with similar waveforms.
The phase picking step has been significantly improved by deep-learning-based pickers that learn from manual picks labeled by analysts to detect millions of phase picks from raw seismic waveforms (Ross et al., 2018;W. Zhu & Beroza, 2018;Mousavi et al., 2020). The rapidly growing volume of automatic picks and the ongoing growth of seismic networks makes developing effective phase association methods crucial. Phase association has not yet received the attention that has been devoted to other earthquake monitoring tasks. Association methods based on back-projection are most commonly used in classic earthquake monitoring systems, such as GLASS3 (Yeck et al., 2019), Earthworm (Friberg et al., 2010), and SeisComP3 (Weber et al., 2007). These approaches usually deploy a grid-search and back-project phase picks based on the expected moveout with distance. An earthquake and its initial location is declared based on the number of phase picks inside a spatial grid that are consistent with a candidate location. Although back-projectionbased association is robust and effective, its performance is limited for dense earthquake sequences when earthquakes occur so closely in time and space that interpreting the maxima resulting from back-projection becomes problematic. Studies have continued to focus on improving the grid-search and back-projection approach for different scenarios (Draelos et al., 2015;Gibbons et al., 2016;Zhang et al., 2019). Meanwhile, several new approaches have been proposed to solve the earthquake phase association using: graph theory (McBrearty, Gomberg, et al., 2019), the RANSAC algorithm (Woollam et al., 2020;L. Zhu et al., 2021), and deep learning (Ross, Yue, et al., 2019;Dickey et al., 2020). When combined with the rapid development of phase picking methods, better association methods have the potential to improve significantly the overall performance of earthquake monitoring pipelines.
We propose an association method based on a Bayesian Gaussian mixture model, which is an unsupervised machine learning method for clustering (Bishop, 2006) that has been widely used in different research fields, such as image processing (Permuter et al., 2006), speech recognition (Reynolds & Rose, 1995), and earthquake studies (Ross et al., 2020;Seydoux et al., 2020). Earthquake phase association can be treated as an unsupervised clustering problem, with groups of phase picks, in time and space, arising from a discrete set of earthquake origins. We combine the Gaussian mixture model with earthquake location, origin time, and magnitude estimation, so that the GaMMA method can cluster phase picks based on the physical constraints of arrival time moveout and amplitude decay with distance. Phase amplitude information is often neglected because it can be difficult to account for in conventional association methods; however, GaMMA is designed such that it can use phase arrival time, phase-type identification, and amplitude information for association while simultaneously estimating the underlying event source characteristics (i.e., location and magnitude). Moreover, GaMMA does not require extra association steps of grid-search or supervised training. These attributes make GaMMA an appealing approach to address the challenges arising in processing of large numbers of automatic picks in earthquake monitoring workflows.

Method
The objective of earthquake phase association is to post-process a large collection of phases picked on individual seismic stations and cluster them into groups of seismic phases originating from a same earthquake event, so that subsequent earthquake characterization tasks can be performed on individual events. Phase arrival times from the same earthquake follow a hyperbolic moveout that is determined by the hypocentral distance and the Earth model (i.e. seismic wave speed). This moveout allows association algorithms to distinguish between phases from different earthquakes. In this work, we extend the association problem by using both phase arrival time and phase amplitude information. On average, the phase amplitude scales with earthquake magnitude and decays with the hypocentral distance. Thus, the amplitude provides additional information to improve phase association. We formulate the association problem as follows: Given N seismic phases (x i , y i , z i , t i , a i ), i.e., arrival time t i and amplitude recorded a i at the i-th seismic station located at (x i , y i , z i ), we seek to group these phases into K earthquakes and estimate the underlying source parameters (x k , y k , z k , t k , m k ), i.e., location (x k , y k , z k ), origin time t k , and magnitude m k of the k-th earthquake.
We solve this association problem using the Gaussian mixture model (Permuter et al., 2006), which is an unsupervised clustering method that groups N data points into K clusters by maximizing the probability that these N data can be explained by a mixture of K Gaussian distributions.
We incorporate the physical constraints on phase arrival time and amplitude into a Gaussian mixture model to make it suitable for our association problem. Figure 1 illustrates how we model the Gaussian distributions to calculate the probability of the sequence of phases generated by the two causative earthquakes. The mathematical details of GaMMA are explained in the following sections.  algorithm to update iteratively: phase assignment, earthquake source parameters, and the mean and standard deviation of Gaussian distribution. This iteration converges to the correct phase assignment and earthquake source parameters to solve the association problem. Note that we use both phase arrival time and amplitude information to model the Gaussian distributions.

Bayesian Gaussian Mixture model
We cast the association problem in a probabilistic framework where we use a Gaussian mixture distribution to model the probability of each phase pick: where x i represents a phase pick including arrival time, phase type, and amplitude (t i , a i ) values at the i-the station. φ k is the mixture component coefficient of the k-th earthquake. N represents a Gaussian distribution, and µ k is the mean of the Gaussian distribution. µ k represents the theoretical phase arrival time and amplitude (t ik ,â ik ) on each i-th station, determined by the k-th earthquake. Λ k is the precision (inverse covariance) matrix of the Gaussian distribution. w i is the phase picking quality score between [0, 1]. n is the number of feature dimensions, which is 1 if only time information is used or 2 if both time and amplitude information are used. Based on the Gaussian mixture distribution, we can calculate the probability of a set of recorded phases (x 1 , x 2 , ..., x N ). We assume these observations (X) are independent and identically distributed (i.i.d.), then the log likelihood function is given by: We can find the assignment from N phase picks to K earthquakes, which is the goal of association, and the corresponding earthquake source parameters by maximizing the log likelihood of (4). A limitation of the Gaussian mixture model formulation is that we need to assume the number of underlying earthquakes K. To address this unknown, we implement the Bayesian Gaussian mixture model (Bishop, 2006), which uses variational inference to calculate approximate posterior distributions for the parameters of a Gaussian mixture distribution. Three conjugate priors are introduced. In particular, we use a Dirichlet prior for the mixture component coefficient, p(φ) = D(α 0 ) (Ferguson, 1973), which controls the concentration of mixture components; a Gaussian prior for the mean conditioned on the precision, p(µ k |Λ k ) = N (m 0 , β 0 Λ k ); and a Wishart prior for the precision, p(Λ k ) = W(W 0 , ν 0 ) (Wishart, 1928), which controls the estimation of covariance. The Bayesian model penalizes parameters that are away from the priors, which balances data fitting and model complexity. The mixture components (i.e., earthquakes) that do not contribute to the explaining the data (i.e., picks) will have approximately zero mixture coefficients so that we can choose a large number of components in the mixture model without over-fitting. In practice, we can initialize the space with many redundant earthquake hypocenters, and the Bayesian GMM suppresses unnecessary sources to infer a accurate number of earthquakes with associated picks.

Expectation-Maximization (EM) algorithm
We use the the Expectation-Maximization (EM) algorithm to solve the maximum likelihood estimation of p(x). To consider the physical constraints on phase arrival time and amplitude for association, we incorporate the estimate of earthquake location, origin time, and magnitude into the EM algorithm. We then iteratively update the assignments from picks to earthquakes in the E-step and optimize the earthquake parameters in the M-step: E-step: where γ ik is the probability that phase pick x i , i.e., arrival time t i and wave amplitude a i , is generated by the k-th earthquake.

M-step:
1. Effective number of picks assigned to the k-th earthquake: 2. Earthquake location, origin time, and magnitude of the k-th earthquake: 3. Theoretical travel time, amplitude, and statistics of residuals: where L is a loss function of the residuals between the picked phase arrival time t i and the theoretical arrival timet ik from the k-th earthquake. Minimization of the loss function l gives an estimate of the earthquake location and origin time (x k , y k , z k , t k ). m k is the magnitude of the k-th earthquake. F t represents the function used to calculate theoretical phase arrival timet ik . F a represents the function to calculate theoretical phase amplitudeâ ik based on earthquake magnitude m k , and F a represents the function to estimate earthquake magnitude using phase amplitude. d ik is the distance from the k-th earthquake to i-th seismic station. Here we decouple the optimization of earthquake magnitude (Equation (9)) from the optimization of earthquake location and time (Equation (8)). We use arrival times to constrain earthquake location and use phase amplitudes to constrain earthquake magnitude. Note that although we decouple the two optimizations, the precision matrix (Equation (11)) considers the correlation between arrival time and amplitude residuals. In this way, both arrival time and amplitude information are used for the association process of clustering picks among earthquakes (Equation (5)). For the Bayesian Gaussian mixture model, we add another stage in the M-step to update the posterior parameters: The E-step is modified as: where α = k α k and ψ is the digamma function (Abramowitz & Stegun, 1964). The explanation of these updating rules is detailed in Bishop (2006)'s textbook.

Earthquake location and magnitude estimation
The iteration of the EM algorithm both updates the clustering of picks based on earthquakes and optimizes the corresponding earthquake source parameters. We focus on efficient association rather than accuracy of earthquake source parameters, which can be realized once phases are properly associated, so we choose two basic approaches to estimate approximate earthquake locations and magnitudes. We optimize Equation (8) with a Huber loss function (Huber, 1992) as the target to reduce the effect of outliers: where the hyper-parameter δ is set to 1 second in this test. For this proof-of-concept study, we use a uniform velocity model to calculate the theoretical phase travel-time: We then solve the minimization of Equation (8) using the BFGS algorithm (Fletcher, 2013). Advanced earthquake location algorithms and complex velocity models can also be applied to solving Equation (8) but at a higher computational cost.
To estimate earthquake magnitude in Equation (9), we use a linear relationship between log phase amplitude and earthquake magnitude: Station correction terms can been added to consider site effects, i.e. site amplification factor (Münchmeyer et al., 2020). Based on the measured phase amplitude type, e.g., displacement, peak ground velocity, or peak ground acceleration, we can choose from among the Richter empirical magnitude relationship (Richter, 1935), the Richter simulation-based prediction (Al-Ismail et al., 2020), or a simplified ground motion prediction equations (Picozzi et al., 2018) for Equation (22).

Results
We demonstrate the performance of GaMMA first on a synthetic example and then on six days of data from the 2019 Ridgecrest, California earthquake sequence.

Synthetic test
We first created a synthetic experiment to demonstrate the association results of GaMMA. We generated a sequence of phases including both P-and S-phases from six earthquake events. To model the errors that exist in real data, we added a 0.5s random error in the phase arrival times and scaled the phase amplitude (peak ground velocity (PGV)) by a random factor between 0.3-3. We further added 30% false positive picks at random times. In total, 178 P-and S-phase picks from 40 stations were used for association (left panels of Figure 2). We used a simplified ground motion prediction equation of PGV from Picozzi et al.
(2018)'s work: log P GV = −2.175−1.68 log R+0.93M , where R is hypocentral distance and M is earthquake magnitude. The ground truth result is shown in the middle panels of Figure 2. The symbol size represents the relative size of phase amplitude and earthquake magnitude. We conducted two association tests using in the first test only the arrival time information (Figure 2(a)) and in the second test both the arrival time and amplitude information (Figure 2(b)). The same parameters and initialization were used for both cases. For each we initialized the earthquake locations at the center of research area and uniformly distributed the earthquake origin times. The association results are shown in the right panels of Figure 2. At least five of the six true earthquakes are successfully associated in both tests. However, the sixth event in the lower right corner of Figure 2(a, b)(ii) (marked in brown) is successfully associated only when amplitudes are used in conjunction with travel times (Figure 2(b)(iii)). Without amplitudes (Figure 2(a)(iii)), the result includes multiple incorrect event associations marked by colors pink and purple. This occurs because the Pphases of this event have arrival times that overlap with another distant event's moveout (the event marked in red) and could be mistakenly associated with this distant event. The extra information provided by amplitude adds the necessary extra constraint on distance, in addition to time, that allows the sixth event to be associated correctly. The test with amplitude information also correctly estimates the earthquake magnitudes. This synthetic experiment demonstrates that GaMMA benefits from using both the time and amplitude information in the association process.

Test on the 2019 Ridgecrest earthquake sequence
We next applied the GaMMA method to part of the 2019 Ridgecrest earthquake sequence to evaluate its performance on a real earthquake sequence. We focused on the initial six days of the sequence when a large number of earthquakes occurred and migrated from the southwest-striking fault to the northwest-striking fault. We applied the PhaseNet model (W. Zhu & Beroza, 2018) to extract picks from waveforms of "HH," "BH," "EH," and "HN" channels of 23 stations of the "CI" network within 1 degree of the location (-117.504W , 35.705N ). We measured the PGV value over a 8s window after the phase arrival time. We then associated the detected 651,994 P-picks and 686,291 S-picks using the GaMMA method. We used a uniform velocity model with v p = 6 km/s and v s = v p /1.75 for earthquake location estimation and a simple ground motion prediction equation as above for earthquake magnitude estimation. Because we used such a simple moveout behavior, we do not expect the earthquake locations to be highly accurate, but they are close enough for successful association, given the relatively short source-receiver distances involved. Figure 3 shows the statistics of the 34,791 associated earthquakes from 598,218 P-picks and 633,010 Spicks, leaving 53,776 P-picks and 53,281 S-picks unassociated. We also plotted the 9,873 earthquakes in the SCSN catalog (SCEDC, 2013) for comparison. Both the associated earthquake locations and magnitudes agree with the SCSN catalog (Figure 3(b) and (c)). Figure 4 shows an association example with a dense sequence of picks occurring during a 8-minute period. GaMMA associates 32 events during this period, while there are only 3 events in the SCSN catalog and 22 events in Ross, Idini, et al. (2019)'s template matching catalog. Figure 5 shows the residual distributions of the associated earthquake location and magnitude compared with the SCSN catalog. The statistics of mean, standard deviation (STD), and median absolute error (MAE) can be found in Table 1. The covariance matrix in Figure 5(d) shows that most of the associated earthquakes have small residuals of phase arrival time and amplitude, indicating that the phase picks match well with the theoretical values determined by the causative earthquakes found by association. These location and magnitude estimates can be further improved through the application of established earthquake location and magnitude algorithms once the picks have been associated.
We compared the catalog generated by GaMMA with three state-of-art catalogs (Ross, Idini, et al., 2019;Liu et al., 2020;Shelly, 2020). Table 2 shows the earthquake numbers in these catalogs during the same period. In each of these cases we assumed these catalogs as ground truth and analyzed whether the earthquakes they contain are also detected in GaMMA's catalog within a 5s window. Based on the recall rate, more than 95% of earthquakes in the catalogs of SCSN, Liu et al. (2020), and Shelly (2020) are successfully associated by GaMMA. The low precision and F1-score are due to the large number of new earthquakes associated by GaMMA. To verity whether these new earthquakes are reasonable, we compared the magnitude distributions of the four catalogs ( Figure 6). Most of these new earthquakes associated by GaMMA have a small magnitude and follow the Gutenberg-Richter magnitude-frequency relationship (Gutenberg, 1956), which suggests that they may be legitimate detections of real earthquakes. Figure 7 shows seismic waveforms of six newly detected events in Figure 4. We can see clear earthquake signals in these examples. Meanwhile, these signals are relatively weak and can only be detected at a few stations. Comprehensive comparison among these catalogs is a subject of future research.

Discussion
The Gaussian mixture model (GMM) is an effective and widely used unsupervised learning method for clustering. We have combined the GMM with earthquake location and magnitude estimation to develop a     Figure 3: Association results of the Ridgecrest dataset: (a) associated earthquake frequency, (b) associated earthquake magnitude, (c) associated earthquake location. Note that because we use a uniform velocity model during association, we do not expect the earthquake locations to be accurate, but they are close enough for effective association.  novel seismic phase association method (GaMMA). We treat earthquake phase association as a clustering problem that aims to cluster phases based on the causative earthquakes. There are several advantageous features of this unsupervised learning approach. First, GaMMA does not require grid-search or training commonly used in other association methods. Second, GaMMA can flexibly consider phase arrival time, amplitude, P/S type, and pick quality. While it is difficult for conventional association methods to consider phase amplitude information, GaMMA can easily include phase amplitude information both to improve association and to estimate earthquake magnitude. Finally, GaMMA optimizes the association result in a probabilistic framework and estimates the covariance of time and amplitude residuals. These advantages suggest that GaMMA is a promising approach for improved earthquake phase association, which is in turn an important component of improved earthquake monitoring. We note that there are several limitations of GaMMA that need to be considered. First, the time complexity of the Gaussian mixture model scales with O(K ·N ), where K is the number of clusters (earthquakes) and N is the number of samples (phase picks), so the computational cost could become prohibitive for a long earthquake sequence. In practice, however, it is straightforward and effective to segment a long sequence into relatively shorter windows to improve the association speed by processing data in parallel, as phases that are separated by a certain time interval (depending on the source-receiver distance) cannot come from a single earthquake. In this work, we used the DBSCAN algorithm (Schubert et al., 2017) to divide picks into sub-windows for association. The example in Figure 4 shows one sub-window from the seven hour sequence. Second, the clustering results of the Gaussian mixture model are affected by the initialization state. In this work, we used a simple strategy to initialize the earthquake locations uniformly in time and space. The number of initialized earthquakes is proportional to the ratio between the number of phase picks divided and the number of stations. This simple strategy worked well in the experiments described above. Improving initialization strategies has the potential to increase the association performance further. Third, GaMMA ensures that one pick is only assigned to one earthquake, while conventional back-projection methods may attribute one pick to several earthquakes; however, GaMMA does not consider station-based constraints, such as that only one pair of P-or S-picks from one station is assigned to each earthquake. McBrearty, Gomberg, et al. (2019) accounts for this constraint using a constrained ILP (integer linear programming) solution in their association method, and there is a strong correspondence between that technique and our technique. One potential solution that would allow introduction of the station-based constraint is to add a  normalizing scheme similar to Equation (5) over stations. Finally, GaMMA assumes that the residuals of pick time or amplitude follow a Gaussian distribution. This assumption is not accurate for false positive picks. Adding a mixture component of background uniform distribution to account false positive picks might be a helpful extension to be considered in future research (Melchior & Goulding, 2018).

Conclusions
We have developed a new association method based on a Bayesian Gaussian mixture model, GaMMA, which solves the phase association problem as an unsupervised clustering problem. To consider the physical constraints of phase arrival time and amplitude with earthquake location and magnitude, we incorporate optimization of earthquake location and magnitude into the Expectation-Maximization (EM) algorithm commonly used for solving the Gaussian mixture model. GaMMA, thus, can use both arrival time and amplitude information to cluster picks from the same earthquake and simultaneously estimates both the earthquake location and magnitude from each cluster of picks. The experiment results on both synthetic tests and the 2019 Ridgecrest earthquake sequence demonstrate the effectiveness of the GaMMA method in associating a dense sequence of P-and S-phase picks. GaMMA provides an unsupervised learning approach to solve the challenging phase association problem resulting from the increasingly wide applications of deeplearning-based phase pickers. The improved performance of GaMMA can associate more earthquakes from massive automatic phase picks, thus enriching earthquake catalogs and improving earthquake monitoring.