Rapid Estimation of Earthquake Magnitude and Source Parameters Using Genetic Algorithms

: To minimize the impacts of large losses and optimize the emergency response when a large earthquake occurs, an accurate early warning of an earthquake or tsunami is crucial. One important parameter that can provide an accurate early warning is the earthquake’s magnitude. This study proposes a method for estimating the magnitude, and some of the source parameters, of an earthquake using genetic algorithms (GAs). In this study, GAs were used to perform an inversion of Okada’s model from earthquake displacement data. In the ﬁrst stage of the experiment, the GA was used to inverse the displacement calculated from the forward calculation in Okada’s model. The best performance of the GA was obtained by tuning the hyperparameters to obtain the most functional conﬁguration. In the second stage, the inversion method was tested on GPS time series data from the 2011 Tohoku Oki earthquake. The earthquake’s displacement was ﬁrst estimated from GPS time series data using a detection and estimation formula from previous research to calculate the permanent displacement value. The proposed method can estimate an earthquake’s magnitude and four source parameters (i.e., length, width, rake, and slip) close to the real values with reasonable accuracy.


Introduction
Earthquakes with large magnitudes always result in large losses and can even trigger a subsequent disaster such as a tsunami. To provide an early warning for tsunamis and facilitate emergency responses to earthquakes, one important parameter that needs to be known early and as accurately as possible is the earthquake's magnitude [1]. Tsunami early warning systems in many countries such as Indonesia use the estimated magnitude of the earthquake as the main reference for providing tsunami early warnings [2,3]. If the earthquake's magnitude is estimated at 6.5 or greater, the system will issue a warning that the earthquake has the potential to trigger a tsunami [3]. The tsunami potential of an earthquake must be decided within 3-5 min. However, the magnitude estimation performed by conventional estimation methods is far from accurate in such a short period of time. In some previous cases of large earthquakes that have triggered tsunamis, the early warning was inaccurate because the magnitude was underestimated at crucial times [4][5][6][7].
In this study, the proposed rapid estimation method using GAs was tested on the 2011 Tohoku Oki earthquake data to measure its ability to estimate the earthquake's source parameters and magnitude.

Materials and Methods
This study aimed to estimate earthquakes' magnitudes and source parameters from GPS time series data recorded during an earthquake event. Here, the parameters were estimated using a GA to inverse the elastic dislocation equation based on Okada's model. The hyperparameters of the GA were tuned to obtain the best configuration settings that yielded the best parameter estimation. The selected configuration was then tested on real displacement time series data from the GPS stations affected by the 2011 Tohoku Oki earthquake. GPS time series data from all stations were processed to estimate the permanent displacement value in both the horizontal and vertical components. The estimated permanent displacement was then inversed using the GA based on Okada's model to obtain the source parameters.

Inversion of Okada's Model
A fault is a large fracture in the Earth's crust that causes two sides of rock blocks to become displaced relative to each other. An earthquake occurs due to a sudden release of energy that causes the fault to slip suddenly. The plane of the fracture where the fault slips is called the fault plane, and is represented by parameters such as strike, dip, depth, and rake [29]. Figure 1a illustrates a fault plane. A slip on the fault plane will cause some displacements on the seafloor of the Earth, which is called surface deformation. Okada's elastic dislocation model, proposed by Okada in 1985, features a set of analytical solutions for surface displacements in a homogeneous half-space. This model translates the effect of a fault-plane slip into the ground displacements represented in u x , u y , u z [30], as illustrated in Figure 1b.
Okada's model is widely used for co-seismic displacement analysis. The displacement field u i (x 1 , x 2 , x 3 ) due to a dislocation ∆u j (ξ 1 , ξ 2 , ξ 3 ) across a surface Σ in an isotropic medium is given by the following formula [31]: A fault in an earthquake can be modeled as a collection of rectangular patches, where each patch has the following parameters [32]: a.
The latitude and longitude of the fault plane's centroid point; b.
The length and width of the fault field (in meters or kilometers); c. The depth (in meters or kilometers); d.
In simple terms, Okada's model formulates these nine parameters to calculate the coseismic displacement fault values, tilts, and compression of an earthquake. This process can be expressed as a function of the nine parameters, which results in the following set of three-dimensional displacements [32]: Okada(E, N, depth, strike, dip, length, width, rake, slip) = [uE, uN, uZ] (2) where E and N are the coordinates of the observation at the geographic reference (east and north) relative to the centroid of the fault plane, and [uN, uE, uZ] is a set of displacement values at the observation points in the directions of north (N), east (E), and vertical (Z). Elastic dislocation modeling, in particular Okada's model, requires the assumption of the fault plane's centroid. In actual conditions, the centroid of an earthquake's fault plane can be assumed to be the hypocenter location obtained from the seismic observations. For the inversion of Okada's model, a calculation in the inverse direction of the model is used to find some or all of the earthquake's source parameters in the model function above from the displacement values at the observation points.

GAs
A GA is an optimization algorithm based on the principle of biological scientific selection. GAs are widely used in linear and nonlinear optimization problems [15,33,34]. In this study, the inversion of Okada's model is also carried out through an optimization approach using GAs.
Based on Equation (2), the inversion process takes displacement data as the input to obtain the nine parameters of Okada's model. However, this study only estimates four parameters through inversion: the length, width, rake, and slip parameters. The other parameters of E, N, depth, dip, and strike are predetermined parameters based on the earthquake fault model used in the experiment. GAs, as heuristic search algorithms, seek the best solutions for the estimated parameters.
In this study, the estimation results of the four variables along with the predetermined parameters are then used as input for the forward calculation of Okada's model to obtain the displacement estimation values. There are two displacement values involved in this inversion process: actual displacement and estimated displacement. The actual displacement is calculated from GPS time series data using the permanent displacement formula, while the estimated displacement is obtained from the estimation using Okada's model.
The generic GA is shown in Figure 2. The algorithm starts by initializing a population that consists of individuals representing solution candidates for the problem to be solved. The fitness of each candidate solution for the problem is then evaluated using an objective function called the fitness function. Two individuals with highest fitness score are then selected for the next parent candidate to produce new offspring through a crossover phase. The next step is the mutation phase to maintain or increase the genetic diversity of the population. This process will continue until the fitness value stabilizes and the algorithm has converged to an optimal solution. The following paragraphs provide more detailed descriptions about the implementation of each step in this study.

Population Initialization
The first step in this phase is to define individual representations as candidate solutions for the source parameter estimation problem. As mentioned previously, in this study the GA only estimates four parameters through the inversion, i.e., length, width, rake, and slip. Hence, in this study, each individual (I) consists of four chromosomes, I = [c 1 , c 2 , c 3 , c 4 ], where c i is the chromosome representing the candidate solution for the ith estimated variable, i.e., length, width, rake, and slip, respectively.
Each chromosome is represented as a number of bits. Therefore, the number of bits parameter (n_bits) is tuned in this experiment. The number of bits will determine the number of possible values generated for each variable. Tuning of the number of bits parameter is described in the scenarios and settings section in Section 3.1. Each chromosome has a predetermined range of values that depends on the estimated variable it represents. According to the range of values for the Okada model parameters previously described in Section 2.1 and the parameter values of fault models tested, the search value for the length parameter ranges from 25 to 750 km, the width parameter from 10 to 300 km, the rake parameter from 60 • to 120 • , and the slip parameter from 0.1 to 25 m. The parameters of length, width, and rake are specified as integer data type, while the slip parameter is of real data type. The fault models tested in this study are described in detail in Section 3.1.1.
In addition to individual representation, another hyperparameter that needs to be determined in this phase is the size of the population, which denotes the number of individuals generated in one population. The size of the population is also included in the hyperparameter tuning of the experiment.

Fitness Calculation
The optimization problem in the source parameter estimation in this study is to find, for the four estimated variables, the optimal solutions that result in the minimum error value between the estimated displacement and the actual displacement. In other words, the closer the estimated displacement value is to the actual displacement value, the more optimal the model parameter estimation value is. Therefore, the objective function used in this study is the minimize function, and the minimized parameter is the error value between the estimated displacement and the actual displacement. The smaller the error value is, the more individuals that will be generated that will fit the optimal solution.
Error value representations used in the objective functions in this study are the mean square error (mse), square error (se), and the number of elements in the displacement matrix with a specific error value. Two objective functions were tested in the experiment to find the optimal estimation as follows: where e dis is the displacement error between the estimated displacement and actual or referred displacement; dis est i is the estimated displacement result from element data i; dis re f i represents the referenced displacement value for element data i, N is the number of data points; and n (e dis j >1) represents the number of elements for which the displacement errors are >1. The estimated displacement (dist est ) is obtained from the forward calculation of Okada's model based on Equation (2) as follows: where E, N, depth, strike, and dip are the predetermined parameters based on the earthquake's fault model used in the calculation as described previously, while v 1 , v 2 , v 3 , v 4 are the estimated variables of the source parameter length, width, rake, and slip, respectively. GA generates the estimation values for the estimated variables. The more optimal the estimation value for the variables, the more minimal the displacement error value will be. In other words, the estimated displacement (dist est ) values are getting closer to the reference displacement (dist re f ) values. Both objective functions in Equations (3) and (4) include the square error between estimated displacement and actual displacement, and the number of elements in the displacement matrix for which the displacement errors are >1. However, the second objective function also includes the mean square error of the displacements. Hence, the second objective function attempts to capture the central value of the displacement errors, in addition to accumulating the error values. One of the objective functions will be selected as the recommended optimal function based on the experimental results.

Parent Selection
The purpose of selecting parent candidates is to obtain better individuals or solution candidates than the previous generation. This study uses the tournament selection method to perform this selection. The selected parent candidate is the fittest individual among k randomly chosen individuals in the population, with k representing a predetermined hyperparameter in the tournament selection method. Individuals with the smallest error values will have the highest probability to be selected.

Crossover
The crossover method used in this study is one point crossover, which designates only one point to split the bit strings of the selected parents. Hence, each pair of selected parents will create two children. The splitting point is randomly selected in the experiment. The bit string representation of each created child refers to the bit string of its parents. The left part of the splitting point in the first child's chromosomes refer to the first parent's bit strings, and the right part of the splitting point refers to the second parent's bit strings. The second child is created by inverting the process.
The crossover process is controlled by a parameter called crossover rate (r_cross), which denotes the probability of the crossover to be applied. Since the crossover concept is something that naturally occurs in real life, the crossover rate is generally set high.

Mutation
A mutation is used to maintain genetic variation within the population. Therefore, mutation is an important step to avoid a local minimum during the search for the optimum solution. Mutation involves bit flipping of the created child. The mutation rate denotes the probability of the mutation, which is typically low and set to 1

L bit
where L bit is the length of the bit strings representing the individual [35]. The formula for the mutation rate (r_mut) is as follows [35]: where n_bits denotes the number of bits in a chromosome, and n_chrom is the number of chromosomes in an individual.
The set values for all GA hyperparameters used and tuned in this study are represented in the scenarios and setting subsection in Section 3. The inversion of Okada's model using the GA will be referred to as the Okada-GA inversion model in this paper.

Moment Magnitude (Mw) Estimation
Moment is related to the total energy released during an earthquake event. Moment is proportional to the slip on the fault multiplied by the fault area that slips. The formula is shown below [36]: where M o is moment seismicity, µ is the rigidity parameter or shear modulus of the crust (approximately 3 × 10 10 N/m), L is the fault length (m), W is the fault width (m), and l s is the slip length (m). Moment seismicity can be converted into an earthquake's magnitude standard called moment magnitude using the following formula [36]:

Permanent Displacement Detection and Estimation
This study uses the permanent displacement calculation formula proposed in [5] to detect and estimate the permanent displacement that occurs at certain observation points. This formula uses a characteristic function, D, based on the parameters of the short-term average and long-term average as follows [5]: where STA(t) and LTA(t) are the short-term and long-term averages of a time series in a certain t, expressed as follows [5]: where x i is the norm of the horizontal components at t = i, α and β are the suitable time window lengths based on the noise properties of the real-time kinematic GPS time series (this study assumes α = 60s and β = 600s), and p i is the weighting parameter [5].
For the standard deviation, the generic formula is as follows: where {x 1 , x 2 , . . . , x i } are the observed values of the sample items, x is the mean of the data observed, and N is number of data points observed. D is a key value for detecting permanent displacement. Permanent displacement is said to occur and be detected when the value of D passes a certain threshold value. The threshold used is expressed in K, which is defined as follows [5]: where SD(D) is the standard deviation of D.
To detect permanent displacement, the formula uses only the norm of the horizontal components to minimize the chance of false detection. The horizontal component has a higher signal-to-noise ratio than the vertical data [5]. However, to estimate the displacement value, the vertical components are also used in the estimation formula.
The algorithm used to detect and estimate permanent displacement is shown in Figure 3. The estimation process begins when the system detects the occurrence of permanent displacement. The process is completed when no more movement is detected at the observation point and D has reached its maximum value. The permanent displacement value is calculated as the difference between the observation point coordinates before and after the earthquake event occurs [5].

Results
Two experimental stages were involved in this study, namely, tuning of the GA hyperparameters in the Okada-GA inversion model and testing of the model using the observation data. The scenarios and settings of each stage along with its results are outlined below.

GA Hyperparameter Tuning Experiments for Okada-GA Inversion Model
The purpose of the first-stage experiments in this study was to test the ability of the GAs to carry out inversion of Okada's model to estimate four earthquake source parameters (i.e., length, width, rake, and slip) and to obtain the hyperparameter setting that resulted in the best performance of the Okada-GA inversion model.

Experimental Scenarios and Settings
In simple terms, the task of the GA in the Okada-GA inversion model was to convert the earthquake displacement value into the estimated parameters (i.e., length, width, rake, and slip) based on Okada's inversion model. In this stage of the experiment, the displacement value used for the inversion was obtained from the forward calculation of Okada's model using a predetermined earthquake fault model. The displacement was referred to as Okada's displacement. As Okada's model involves observation points in its forward calculation, this stage of the experiment used 737 GPS stations related to the 2011 Tohoku Oki earthquake as the observation points. Hence, the fault model used in the experiment was also the fault model of the 2011 Tohoku Oki earthquake event. To obtain valid results for the GA performance when performing the inversion of Okada's model, more than one fault model was used for the experiment. The four predetermined earthquake fault models used in this stage of the experiment are shown in Table 1. The first model was based on prior knowledge about the 2011 Tohoku Oki earthquake, the second model referred to a study published formally by the Tectonics Observatory California Institute of Technology [37], and the third and fourth models were based on models published by the Geospatial Information Authority of Japan (GSI) [38]. Figure 4 shows the fault model and distribution of observation points in the experiment.
To obtain a reliable estimation and represent the actual conditions of the data, some random error values were added to Okada's displacement; these values were about 2-5 mm for each displacement value in the horizontal component. Both the forward and the inversion calculations of Okada's model in this experiment involved the 737 GPS stations as observation points. As noted earlier, all parameters of each fault model presented in Table 1 were used as input for the forward calculation of Okada's model to generate Okada's displacement. In this stage of experiment, the Okada's displacement of each fault model was inversed using the GA to obtain four source parameters (i.e., the length, width, rake, and slip of the fault plane) by searching for the optimum solution that yielded the minimum fitness function. The range of search values for each estimated parameter and its data type were previously described in Section 2.2.1. From the estimated parameters obtained using the GA, the estimated displacement was calculated using Okada's model. In this study, the estimated displacement was referred to as GA displacement.
GA performance strongly relies on the hyperparameter settings. Therefore, this study also included the tuning process in its experiments. In this stage, the hyperparameters of the GA were tuned to obtain the optimal estimation results. These hyperparameters included the number of generations (n_gen), size of the population (n_pop), number of chromosome bits (n_bits), crossover rate (r_cross), and mutation rate (r_mut). The previously mentioned tournament selection hyperparameter k was set as a fixed parameter in the experiment: k = 4. Each hyperparameter combination generated in the tuning process was tested on the four fault models using the two objective functions described previously. At the end of the tuning process, the optimal set of hyperparameters along with the optimal objective function were selected for each fault model.
The values set for each hyperparameter tuned in the experiments are outlined in Table 2. According to the number of values of each hyperparameter, 108 hyperparameter combinations were generated for each combination of fault model and objective function. Hence, a total number of 864 experiments were run in the hyperparameter tuning process, with different hyperparameter settings for each. The Okada-GA inversion model was run for each experiment. All experiments in this study were run using Python programming on a computer with the following hardware specification: 16 GB RAM DDR4 and 4.2 GHz processor.

Hyperparameter Tuning Results
Based on the experimental results obtained from the hyperparameter tuning process, the GA performance in the Okada-GA inversion model was measured using the following parameters: (1). The displacement error value, which was defined as the difference between GA displacement and Okada's displacement. This displacement error value was represented by using the RMSE of all the displacement values generated at each observation point.
The error value was referred to as the error_disp parameter; (2). The estimated magnitude Mw, which was generated by the GA estimation compared with the reference magnitude value generated from the Mw calculation based on the length, width, and slip parameter values in the forward calculation of Okada's model.
The best estimation result was the estimation with the minimal displacement error value and magnitude value error for all conducted experiments. Table 3 shows the best estimation result for each fault model tested in this stage of the experiment, compared with their original values as presented earlier in Table 1.   The hyperparameter combination settings shown in Table 3 that resulted in the best estimation performance were then selected to be used in the next stage of the experiment. The settings here included the number of generations, size of the population, number of chromosome bits, crossover rate, and mutation rate for each fault model. Table 4 shows the best configuration settings recommended by the tuning process for each fault model used. In the following experiments, the selected objective function and hyperparameter setting shown in Table 4 were used to test the Okada-GA inversion model on real observation data, i.e., the GPS timeseries data from the 2011 Tohoku Oki and 2018 Lombok earthquake events.

Testing the Okada-GA Inversion Model Using the Observation Data
The testing in this step aimed to measure the ability of the Okada-GA inversion model to estimate source parameters from the real data of an earthquake. The estimated parameters were the same as those used in the previous stage (i.e., length, width, rake, and slip), in addition to the magnitude Mw of the earthquake. In this stage, the GA with the best hyperparameter settings and the selected objective function was tested on the real displacement time series data from two earthquake events, i.e., the 2011 Tohoku Oki and 2018 Lombok earthquake events.
For all experiments carried out in this stage, the difference between the estimated parameters and reference parameters was expressed as the RMSE using the following standard formula: where x i est is the ith estimated value, x i re f is the ith reference value, and N is the number of data points.
The results of the testing for each data event are described in the following sections.

The 2011 Tohoku Oki Earthquake Event
The time series data used for the testing were 1 Hz real-time kinematic-GPS data recorded at 737 observation points during this major earthquake. Here, the displacement of the 2011 Tohoku Oki earthquake was first calculated using the permanent displacement detection and estimation formula [5] described earlier. The displacement result from the formula using the real GPS time series data is referred to, in the next section, as the referenced permanent displacement. From the permanent displacement detection and calculation process, information about the permanent displacement detection and estimation time needed by each GPS station was obtained. Figure 6 shows the number of GPS stations involved in the calculation of the permanent displacement completed at a certain elapsed time t. The estimated permanent displacement was then inversed using the GA to obtain the estimated source parameters. The permanent displacement results were used as input for the Okada-GA inversion model. The fault centroid, depth, strike, and dip parameters used for the inversion were the same parameters of the fault models used in the previous stage of the experiment. Comparison between the GA displacements and the referenced permanent displacement can be seen in Figure 7. The blue vectors are GA displacements, and the black vectors are the referenced permanent displacements of the 2011 Tohoku Oki earthquake event. The red pentagram indicates the centroid of each fault model.  Table 5. Table 5 also presents the error_disp between the GA displacement and referenced permanent displacement obtained from the permanent displacement formula [5] using the real time series data of 2011 Tohoku Oki earthquake event. The error_disp is illustrated earlier in Figure 7. From the estimated results of the length, width, rake, and slip of the parameters, the Mw value was estimated using Equations (7) and (8). Tables 3 and 5 show the estimated results of Mw for the first and second stages of the experiment, respectively.
In this stage of the experiment, in addition to the displacement results to be used as input into the GA, information about the elapsed time required by each GPS station was also retrieved to complete the permanent displacement formula as shown in Figure 6. Here, t pd is the elapsed time required to complete the calculation of the estimated permanent displacement, and t GA is the GA computational time required to obtain the parameter estimation results and Mw, calculated from the time of the permanent displacement is input into the GA. Hence, the total elapsed time required by the Okada-GA inversion model proposed in this study to complete the source parameters and magnitude Mw estimation (t est ) was as follows: Table 6 shows the GA computational time required (t GA ) to obtain the estimated source parameters and the estimated Mw of the 2011 Tohoku Oki earthquake for each fault model tested in the experiment using the previously mentioned hardware specification. According to the data shown in Figure 6, Equation (14), and the t GA data presented in Table 6, the progress of the Mw estimation could be observed per elapsed time, as shown in Figure 8

The 2018 Lombok Earthquake Event
The Okada-GA inversion model was also tested on the 2018 Lombok earthquake's time series data. Due to limited observation data, the time series data used for testing were 1 Hz real-time kinematic-GPS data recorded at one observation point during this major earthquake. As in the testing experiment for the 2011 Tohoku Oki earthquake described previously, the displacement of the 2018 Lombok earthquake was first calculated. The displacement results were then used as input data for the Okada-GA inversion model.
In this testing experiment, the predetermined parameters used for the inversion, i.e., hypocenter location, strike, and dip, referred to the focal mechanism information published by the International Seismological Centre (ISC) [39], as can be seen in Table 7. The estimation results of the parameters length, width, rake, and slip were 73 km, 73 km, 60 • , and 0.171 m, respectively. Figure 9 shows the comparison between the GA displacements and the referenced permanent displacement of the 2018 Lombok earthquake case. The GA displacement and the referenced displacement are represented by the blue vector and the black vector, respectively. The red pentagram indicates the epicenter of the 2018 Lombok earthquake. From the source parameters estimation results, the Mw estimation was calculated using Equations (7) and (8). The estimated Mw result was 6.898, with the t pd = 64 s and the t GA = 28.263 s, using the hardware specification previously described. Hence, the total elapsed time required by the Okada-GA inversion model proposed in this study to complete the source parameters and magnitude Mw estimation (t est ) of the 2018 Lombok earthquake event was 92.263 s.

Discussion
The estimation results of the first stage of the experiment in Table 3 showed that error_disp was, on average, 3.784 × 10 −3 m for all the models or 3.572 × 10 −3 m for the first model as the minimum error and 4.339 × 10 −3 m for the maximum error value. These values indicated that the Okada-GA inversion model generated estimated displacements close to the referenced displacements. The estimation results for each parameter in each model were also close to the reference values in all fault models tested. This meant that the Okada-GA inversion model proposed in this study could provide relatively accurate estimation results. Figure 5 also showed that the GA displacements almost coincided with Okada's displacements.
However, the estimation results in the second stage of the experiment were not as accurate as the first stage results. The average of error_disp in the second stage of the experiment (0.105 m) was larger than error_disp generated in the first experiment (3.78 × 10 −3 m). In the first stage, the referenced displacement values were Okada's displacements, which were fixed values generated from the forward calculation of Okada's model. Meanwhile, in the second stage, the reference values were the real permanent displacement values estimated from the permanent displacement formula [5], which presented significant errors compared with Okada's displacements.
According to the calculations using all the observation points, the average of error_disp in the RMSE between the real permanent displacements and Okada's displacements was 0.083 m for the N components and 0.148 m for the E components. For the vertical components, the RMSE of error_disp in the Z components was 0.098 m.
The Okada-GA inversion model proposed in this study used Okada's model as the main component of its objective function. Hence, the smaller the error_disp was between the real permanent displacement and Okada's displacement, the smaller the RMSE was between GA displacement and the reference value, which was Okada's displacement. Therefore, when the Okada-GA inversion model estimated the displacement values with real permanent displacement as its reference values, it was reasonable for the error difference to reach 0.105 m, but not when the Okada-GA inversion model estimated the displacement values using Okada's displacements as its reference values.
However, the estimated result of each source parameter (length, width, rake, and slip) shown in Table 5 was reasonable according to the 2011 Tohoku Oki earthquake's source model described and analyzed in previous research [40][41][42]. Some fault plane models analyzed in previous research ranged from 600 to 812 km in length and from 109 to 240 km in width [41,42], while the results estimated in this study used 668 km for the length parameter and 234 km for the width parameter in the third fault model tested. The range of estimated width parameter values that resulted from this study was between 153 to 234 km. The slip estimation value ranged between 25 and 60 m in most previous research [40,41], while this study estimated slips in the range of 10 to 24.598 km.
The Mw estimation results in both the first and second stages of the experiment were also close to the reference values in each model, both in the 2011 Tohoku Oki and 2018 Lombok earthquake events. Furthermore, Figure 8 showed that the proposed method provided a relatively stable estimation value from early observations (elapsed time t = 173.81 s)that is, as soon as a permanent displacement was detected and estimated. In the testing experiment on the 2011 Tohoku Oki earthquake's data, the first Mw estimations generated using each fault model in testing were 8.8, 8.823, 9.14, and 9.067 for fault model 1, 2, 3, and 4, respectively. In the 2018 Lombok earthquake testing, the Mw estimation generated based on the single observation point was 6.898 which was close to the referenced Mw, 6.9. Unlike the estimation results in previous studies, where the magnitude generally started with a low value [5,7,9,[11][12][13], the first generated magnitude value that was estimated in this study was already close to the reference value and then increased or decreased slowly until it stabilized.
According to the data shown in Figures 6 and 8, the first Mw estimation generated in elapsed time t = 173.81 s was generated from the observation data of a single GPS station. The estimation value was relatively stable compared with the Mw estimations at the following elapsed times generated from multiple GPS stations. In the case of the 2018 Lombok earthquake, due to the limited observation data, the Mw estimation was also generated from a single GPS station. However, the estimation result was close to the referenced value of Mw. It meant that in the proposed method, the number of GPS stations used as observation points did not significantly affect the values of the estimated parameter results.
The two main processes in the proposed model were permanent displacement estimation based on GPS time series data and the inversion of Okada's model via GA to convert the displacement data into estimated source parameters. Hence, the proposed model has the potential to be applied to almost any earthquake case, as long as the earthquake causes surface deformation and significant permanent displacement that can be estimated using GPS time series data. However, the proposed model still requires some predetermined source parameters to perform the inversion, such as strike and dip parameters. Therefore, for different earthquake types such as an outer rise event, there may be constraints that need to be studied further in order for the model to estimate properly.

Conclusions
One crucial parameter must be estimated carefully when an earthquake occurs: the earthquake's magnitude. Underestimating an earthquake's magnitude will result in an inaccurate early warning and emergency response. This paper proposed a method using a GA to estimate four source parameters (i.e., length, width, rake, and slip). The proposed method can also estimate the moment magnitude with greater reliability and accuracy than previous methods, at least for the earthquake case analyzed in this study: the 2011 Japan Tohoku Oki earthquake. However, this method currently only estimates four source parameters and needs prior information to carry out the inversion. Hence, the model proposed in this study should be improved in future work.

Conflicts of Interest:
The authors declare no conflict of interest.