A statistical methodology for classifying earthquake detections and for earthquake parameter estimation in smartphone-based earthquake early warning systems

Smartphone-based earthquake early warning systems (EEWSs) are emerging as a complementary solution to classic EEWSs based on expensive scientific-grade instruments. Smartphone-based systems, however, are characterized by a highly dynamic network geometry and by noisy measurements. Thus, there is a need to control the probability of false alarms and the probability of missed detection. This study proposes a statistical methodology to address this challenge and to jointly estimate in near real-time earthquake parameters like epicenter and depth. The methodology is based on a parametric statistical model, on hypothesis testing and on Monte Carlo simulation. The methodology is tested using data obtained from the Earthquake Network (EQN), a citizen science initiative that implements a global smartphone-based EEWS. It is discovered that, when the probability to miss an earthquake is fixed at 1%, the probability of false alarm is 0.8%, proving that EQN is a robust smartphone-based EEW system.

/fams. . of monitoring smartphones is small. Both false alarms and missed detections may undermine people's trust in the EQN. In the pivotal study by Finazzi and Fassò [15], a statistical methodology is developed for identifying in real-time earthquake occurrence. The study, however, does not take into account the spatial dimension of the smartphone network, making the detection algorithm prone to false alarms. Moreover, the methodology does not allow to estimate important earthquake parameters such as epicenter and depth. In Finazzi et al. [16], instead, the EQN detection capabilities are modeled within a probabilistic framework. It is discovered that the EQN missed some relatively strong earthquakes that were supposed to be detected by the smartphone network. These considerations and findings suggest that there is room to improve EQN's methods and algorithms.
This study proposes a statistical methodology for 1) controlling the probability of false alarms, 2) controlling the probability of missed detection, 3) classifying a detection between true and false earthquake, and 4) estimating earthquake epicenter and depth (if the detection is classified as a true earthquake).
The methodology is based on a statistical parametric model, statistical hypothesis testing, and Monte Carlo simulation. Contrary to model-less approaches (see for instance [3]), the methodology exploits the fact that the spatio-temporal dynamic of seismic waves is well-known. This information is retained by the statistical model, and it helps to both classify the EQN detection and to estimate the earthquake parameters.
Due to the peculiarity of the specific application, real-time is a constraint. Ideally, classification and earthquake parameter estimation should not exceed 1 or 2 s of computing time.
The smartphone-based EQN is used to test the statistical methodology, which is then applied to some true and false EQN detections.

. EQN's detection algorithm
Before formalizing the classification and the earthquake parameter estimation problems, it is useful to detail the output of the earthquake detection algorithm currently implemented by the EQN [15]. For any given area of radius 30 km, the algorithm compares the number of triggering smartphones in the last 10 s with the number of active smartphones. A triggering smartphone is a smartphone that detected an acceleration above a threshold, while an active smartphone is a smartphone known to monitor earthquakes. If the ratio between triggering smartphones and active smartphones exceeds a threshold, an earthquake is claimed to be detected. The output of the detection algorithm consists of the detection location and the list of the triggering smartphones (triggers for short), which are identified by their spatial coordinates (latitude and longitude) and the triggering time.
. Problem formalization An earthquake detection made by an EQN is defined in terms of k j > 0 triggers, where j is the index of the generic detection. In general, k j is not a constant, meaning that each detection is characterized by a different number of triggers. Each trigger is described by the feature vector as follows: where t i ∈ R is the triggering time, while lat i , lon i ∈ S 2 are the smartphone coordinates, with S 2 being the sphere embedded in R 3 .
is the data point, and the feature Let Y = {−1, 1} ∋ y be the label space. For each earthquake detection, y = 1 if the detection is false while y = −1 if the detection is related to a true earthquake.
The aim is to learn a hypothesis map h : X → Y such that y ≈ h (X) for any data point X (i.e., for any future EQN detection). The map h is highly non-linear since the information content of X is determined by the spatio-temporal dynamics of the seismic waves and spatial distribution of the smartphones at the time of the earthquake.
A statistical parametric model f : X → is adopted to understand if X is generated by a true earthquake. The unknown model parameter vector is θ ∈ = R s , with s ≪ k j as the vector size. The hypothesis map is then h (X) = g f (X) = g (θ ). Note that s is constant, and it does not depend on the dimension of X.
When dealing with EEW systems, it is required to control two parameters: the probability α of missed detections (true earthquakes which are not detected by the system) and the probability β of false detections (detections which are not related to any occurred earthquake). It is thus reasonable to adopt a 0/1 loss function as follows: L X, y , g = 1 if y =ŷ 0 else, and to learn a g that minimized the Bayes risk g = arg min g E L X, y , g .
As discussed by Jung [17], solving (Equation 1) requires knowing the joint probability distribution p X, y . Instead, we rely on the fact that it is relatively easy to simulate EQN detections under different smartphone geometries and different earthquake parameters. This induces a variability on X and on the number of triggers k j . Assuming to have a data set D = X (1) , y (1) , ..., X (m) , y (m) and that D is a representative sample of p X, y , we define the empirical risk as follows: and g is learned from the following minimization problem: Note that solving (Equation 2) is equivalent to solvê . /fams. .
where it is made explicit that the probabilities of missed and false detections depend on g.
From an EEW perspective, the solution provided by Equation (3) is not necessarily the best. In some contexts, a missed detection has a larger negative impact than a false detection, while in other contexts, it is the opposite. In this case, one probability is fixed to the desired level, and the other probability is minimized. Two other minimization problems for learning g are the following: . Statistical parametric model and classification In this section, we propose a statistical parametric model for the generic data point X. The observed triggering time for a smartphone sensing an earthquake is modeled as where t * i is the expected triggering time, while with as the distance between the hypocentre and the smartphone location, v is the seismic wave speed, and t O ∈ R is the earthquake origin time.
In Equation (8), D i,E is the distance between the epicenter (lat E , lon E ) ∈ S 2 and the smartphone location, d E ∈ [0, 500] is the earthquake depth, and R is the earth radius (6, 371 km). Here, it is assumed that all smartphones either detect the primary seismic wave (v = 7.8 km/s) or they all detect the secondary wave (v = 4.5 km/s). This assumption is justified by the fact that earthquake detection is based on smartphones within a radius of 30 km, which is a relatively small area.
The role of the random component ǫ i is to model the difference between the expected and the observed triggering time. This difference is mainly due to the smartphone detection delay and a seismic wave velocity that may differ from the expected value.
Equations (6)(7)(8) fully define the statistical model f and the model parameter vector is θ

. . Model estimation
Model estimation is based on the maximum likelihood method. For a generic EQN detection, the log-likelihood function based on the joint probability distribution of The t i are assumed to be independent. This assumption is realistic because smartphones do not share a common clock, detection delays are independent, and the detection by each smartphone is influenced by local factors (e.g., where the smartphone is located, at which floor of the building, and the accelerometer sensitivity).
Maximum likelihood estimates of lat E , lon E , d E , and t O are given by arg min The solution of Equation (10) cannot be obtained in a closed form due to the non-linearity of Equation (8) hence, estimates are obtained via numerical optimization using the BFGS Quasi-Newton method [18]. As usual, to avoid local minima, the numerical optimization algorithm is run multiple times starting from random initial values for lat E , lon E , d E , and t O . The minimization in Equation (10) is possible because for any "proposed" values of the model parameters, t * i can be computed using Equations (7), (8) and then compared with the observed t i .
At convergence, the BFGS quasi-network method also returns the Hessian matrix. Since maximum likelihood estimates for model parameters are obtained from a minimization problem, the Hessian is equivalent to the observed Fisher information matrix. The variance-covariance matrix of the three parameters is then the inverse of the Hessian matrix from which standard errors are easily computed.
Finally, the maximum likelihood estimate of the variance is as follows: whereˆ t i = t i −t * i is computed after replacing in Equations (7) and in Equation (8) the maximum likelihood estimates of latitude, longitude, and depth, whileμ is the mean of theˆ t i .

. . EQN detection classification
Among all elements of θ , the parameter that carries information about how the EQN detection should be classified is σ 2 ǫ . Indeed,σ 2 ǫ tends to be small when the earthquake is true (and triggering times follow the seismic wave dynamic) whileσ 2 ǫ tends to be large when the detection is not related to an earthquake event. This implies that g (θ) reduces to g σ 2 ǫ . In this study, g is chosen to be a statistical hypothesis test on σ 2 ǫ . The system of hypothesis is given by The null hypothesis is rejected when the variance is higher than expected, namely, when smartphone triggering times do not follow the propagation law of the primary or secondary seismic wave. As customary in the statistical hypothesis testing, the probability α is fixed, and it represents the probability to reject the null hypothesis when it is actually true (namely, it is the probability to miss a true earthquake).
The test statistic is as follows: which, under the null hypothesis, is distributed as a chi-square with k − 4 degrees of freedom (df ), where 4 is the number of estimated parameters in Equation (10). The null hypothesis is rejected if T > q (1−α),df , whereT is obtained replacing σ 2 ǫ withσ 2 ǫ in Equation (13), while q (1−α),df is the (1 − α)-quantile of a chi-square distribution with df degrees of freedom, usually called the critical value. In practice, an EQN detection is a true earthquake unless data bring enough evidence that the detection is actually false.
Since we do not know which seismic wave is detected by the smartphones, two models f are estimated: one with v = 7.8 km/s and another with v = 4.5 km/s in Equation (7). This brings to two estimated values for σ 2 ǫ and two hypothesis tests are implemented. The detection is classified as a false earthquake if the null hypothesis is rejected under both tests; otherwise, the earthquake is classified as true. It is worth noting that the statistical hypothesis test is equivalent to a linear map. Indeed, setting then g = w ′ φ, and the earthquake detection classification is based on the following rule: Finally, δ is obtained by solving the problem Algorithm 1 summarizes the steps for classifying an EQN detection and for estimating the earthquake parameters in case the detection is classified as a true earthquake.

. Simulation study
The minimization problem in Equation (17) has no closedform solution. For this reason, we implement a Monte Carlo simulation that aims to simulate a data set D and to minimize Equation (17).
A total of 1,000 true EQN detections and 1,000 false EQN detections are simulated considering the true locations of 1,000 smartphones of the EQN in Lima (Peru).
The probability of missed detection is fixed to α = 0.01 while δ is made varying from 0.1 to 1.5 with step 0.1. For each value of δ, β (δ) is computed by estimating the model f and by implementing the hypothesis test (Equation 13) overall data points X (j) in D. Finally,δ is the value of δ that minimizes β (δ).

. . Simulation of true detections
For simulating a true earthquake, the following aspects are taken into account: the earthquake epicenter and depth, the arrival time of the seismic wave at the smartphone locations, the earthquake detectability by the smartphone, and the error on the triggering time. Finally, we account for the fact that smartphones may detect events unrelated to the earthquake.   the error on the triggering time is simulated from a zero mean normal distribution with variance σ 2 ε = 1.67. Such variance guarantees that the 1st and the 99th percentiles of the error distribution are around −3 and 3 s, respectively, which are realistic values for an error on the triggering time.
Of the remaining 30% of smartphones which do not trigger, 6% are made triggering at random with a triggering time uniformly .
/fams. .  generated in the range [0, 12] s. This implies that when the earthquake is detected by the EQN detection algorithm, the list of triggering smartphones may include triggers unrelated to the earthquake dynamic.
Once the list of triggering smartphones is defined and sorted by triggering time, the EQN detection algorithm is applied to the list. The algorithm stops when the detection condition is satisfied, and the sub-list of triggers that concurred with the earthquake detection is given as the output. Figure 1 shows an example of a simulated true earthquake. Two separated regions can be visually identified, one with triggering smartphones (those that concurred with the detection) and another with non-triggering smartphones not yet reached by the seismic waves.

FIGURE
EQN triggers for the earthquake occurred on October close to Genoa (Italy). The diameter of circles is proportional to the triggering time. The number of triggering smartphones is n = 21, and 99% confidence intervals is presented in brackets. Real earthquake parameters are taken from the website of the European-Mediterranean Seismological Centre (https://www.emsc-csem.org).

. . Simulation of false detections
To simulate a false detection, we assume that smartphones trigger at random with a triggering time that does not follow the law of seismic wave propagation. Only 30% of the smartphones are made triggering, and the triggering time is uniformly sampled in the range [0, 12] s. Figure 2 shows an example of a simulated false EQN detection. Contrary to true earthquakes, no specific spatial pattern on the triggers is observed.

. . Simulation results
The minimization of Equation (17) is attained whenδ = 0.6 and β is found to be equal to 0.008 (conditionally on α = 0.01). Figure 3 shows the empirical distributions ofσ 2 ǫ for both true and false simulated EQN detections. Although the detection classification is based on the hypothesis test (and not directly onσ 2 ǫ ), the overlapping between distributions suggests that classification errors are possible.
A by-product of detection classification is the estimate of the earthquake parameters. Figure 4 shows the box plots of errors .
/fams. . on earthquake epicenter and depth. Both errors have a median of around 18 km, suggesting that along with the detection classification (true/false), the model output can be exploited to provide preliminary estimates of the earthquake parameters.

. Real data example
The methodology developed in this study is applied to true and false detections made by the EQN. As a true earthquake, the event occurred near Genova (Italy) on 4 October 2022 at 21:41:10.5 UTC is considered. Figure 5 depicts the triggering smartphones (n = 21), while estimation and classification results are reported in Table 1 for v = 7.8 and v = 4.5 km/s, respectively.
For both seismic wave velocities, we can observe that latitude and longitude are accurately estimated, while the error in depth is not negligible. Nonetheless, the true values are within the 99% confidence intervals evaluated from the standard errors on the model parameters. In addition, the earthquake is classified as true under both velocities since both observed test statistics are lower than the test critical value. This happens because triggers are close to the epicenter, and primary and secondary seismic waves are nearly concurrent.
The estimation and classification results were obtained in less than 1 s using an Intel(R) Core(TM) i7-9750H CPU @2.60GHz, suggesting that the approach can be adopted for real-time applications. Figure 6 shows the n = 108 triggers of a false detection occurred near Acapulco (Mexico) on 25 September 2022, at 09:55:45 UTC. In this case, the computed test statistics are 1039.7 and 1026.0 for v = 4.5 and 7.8 km/s, respectively, while the critical value is 141.62. H 0 is rejected in both cases and the detection is claimed as false. In this particular case, the detection was caused by a strong lightning bolt. The speed of sound, however, is around 0.3 km/s, a value much smaller than the speed of primary and secondary seismic waves.

. Discussion
The methodology developed in this study allows to classify detections made by smartphone-based earthquake early warning systems between true (related to a real earthquake) and false. This is done analyzing the information content of the smartphone triggers that contributed to the detection.
With respect to classic classification problems, the data point describing the triggers has a varying dimension which depends on the smartphone network geometry. The proposed solution is based on two steps. First, a statistical parametric model is used to convert the data point into a parameter vector with a fixed (and small) dimension. Second, a hypothesis test is implemented for classification.
While we do not claim our choices of f and g to be optimal, both steps are based on well-established statistical methods. With respect to the specific choice of g, it is worth discussing that a simpler alternative is the linear map g * = δ ′ φ, with δ = (δ, 1) ′ and φ = 1, −σ 2 ǫ ′ . In this case, the classification is based on the more intuitive comparisonσ 2 ǫ ⋚ δ. This simpler solution, however, does not take into account neither the actual number of triggers for the specific detection (10 or 1,000 makes a difference in the uncertainty ofσ 2 ǫ ) nor the fact that the distribution of σ 2 ǫ is known under the null hypothesis (that the detection is related to a true earthquake). Using hypothesis testing, we are thus able to retain a part of the information which is lost when X is synthesized with θ .

. Conclusion
Classification and earthquake parameter estimation are performed in near real time, making the statistical methodology suitable to be implemented in operational systems. On the contrary, the methodology does not fully exploit the information available on the EQN system. Specifically, the modeling is only on the triggering smartphones, while the active nontriggering smartphones are ignored. Knowing, at the EQN detection time, which smartphones have not (yet) triggered may better constraint epicenter and depth, thus improving their estimates.
In addition, for an EEWS like EQN that works globally, it would be important to study if the data set D generated by the Monte Carlo simulation is a representative sample of p X, y . If not, the observed α and β probabilities might deviate from the expected ones.
Finally, a limit of the approach proposed by this study is that the statistical methodology is applied downstream of EQN detections. Ideally, the detection, the classification, and the earthquake parameter estimation problems should be jointly addressed in a unified approach. In this regard, the vast literature on wireless sensor networks may help propose a solution under the real-time constraint.
These open problems, along with the estimation of the earthquake magnitude, will be the focus of future works.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions
FF: conceptualization, writing-review, and editing. FM: investigation, methodology, validation, and writing-original draft preparation. All authors contributed to the article and approved the submitted version.

Funding
This article was funded by the European Union's Horizon 2020 Research and Innovation Program under grant agreement RISE No. 821115.