Multiple inversions of Rayleigh wave dispersion curve for geotechnical site characterization using particle swarm optimization dan genetic algorithm

: The inversion of the Rayleigh wave dispersion curve is a crucial step in the multi-channel analysis of surface waves (MASW) method, used to obtain the shear wave velocity ( Vs ) profile. The nonlinear and multimodal nature of the dispersion curve makes a global optimization approach, such as particle swarm optimization (PSO) and genetic algorithm (GA), the optimal choice for inversion. This study aims to compare the performance of multiple inversions of PSO (MI-PSO) and multiple inversions of GA (MI-GA) in solving the inversion problem of the Rayleigh wave dispersion curve. The test results indicate that the utilized MI-PSO outperforms MI-GA in terms of computational time and accuracy of the obtained model


Introduction
Geotechnical site characterization using the inversion of Rayleigh wave data has developed in recent years.This is because the inversion of Rayleigh wave data is able to provide a profile of shear wave velocity (Vs) versus depth with good accuracy, where these two parameters are important parameters in geotechnical studies such as research on liquefaction, shallow cavities, subway routes, and bedrock layers (Karray, 2010).In the multichannel analysis of surface waves (MASW) method, Rayleigh wave dispersion curve inversion is the last step in data processing (Park et al., 2007;Le Ngal et al., 2019).
Inversion of the Rayleigh wave dispersion curve, like many geophysical problems, is nonlinear (Dal Moro et al., 2007).As a result, the linearized inversion approach (Song et al., 2012), a locally based optimization method for getting shear wave profiles (Vs) versus depth, has difficulty obtaining suitable estimating models.This occurs because the accuracy of linear inversion results is primarily controlled by the a priori information and initial model provided.As a result, global optimization is a better technique for overcoming the non-linear nature of the Rayleigh wave dispersion curve inversion.Global optimization methods that are often used in geophysical data inversion are genetic algorithms (La Hamimu et al., 2011;Safani et al., 2013;Rubayin et al., 2019), simulated annealing (Beaty et al., 2002;Lu et al., 2016;Wang et al., 2021), and particle swarm optimization (Kennedy and Eberhart, 1995;Raha et al., 2019;Oyeyemi et al., 2023).
Even though global optimization methods are more suitable for solving non-linear problems, the accuracy of the inversion results produced by each global optimization method may be different.So, this research looks at how well the particle swarm optimization (PSO) and genetic algorithm (GA) methods work for inversion.In each of these optimization approaches, multiple inversions are carried out, and then the parameters of the inversion results (i.e., shear wave velocity and depth) obtained in each single inversion are weighted.

Multiple Inversions with Particle Swarm Optimization (MI-PSO)
The MI-PSO algorithm we propose in this research refers to a single inversion of PSO by Raha (2019).The stages of the MI-PSO algorithm are as follows: 1. Determine the upper limit (  ) and lower limit (  ) of the shear wave velocity (Vs) and layer thickness (d) as the particles to be estimated.The combination of boundaries between the two particles is arranged in the following matrices: as the upper bound and   = [  1 …   +1   1 …    ] as the lower bound (where subscript D indicates the number of finite-thickness layers).In the initial settings, the number of particles/population M, the maximum number of inversions Imax, and the maximum number of iterations kmax are also determined; 2. Loop the inversion I from 1 to Imax; 3. Initialize the initial population   randomly with an upper bound   and a lower bound   .Initialize the initial velocity of all particles with v0 = 0. Calculate the misfit of each particle.Initialize the best position of each particle (p i ) from the initial position of each particle.Selection of the global best position (g) of the particle with the best misfit obtained 4. Loop the iteration I from 1 to k max , where at each iteration the following steps are carried out: • Update particle velocity using Equation (1),   ( + 1) =   () + 1(  −   ())  + 2( −   ())  (1) and the conditions for the maximum velocity value (v max ) utilize the following equation: (2) where   is the velocity component of the i th particle, t and t+1 denote two consecutive iterations of the algorithm,   is the "personal best" of the particle, which is the best solution obtained during optimization by the individual, and  is the "global best," which is the overall best solution obtained by the group.The acceleration coefficients c1 and c2 are known as "cognitive coefficient" and "social coefficient," respectively, and are set to the values c1 = 2 and c2 = 2, respectively.  and   are two diagonal matrices of random numbers from a uniform distribution between [0, 1].w is the inertial weight to balance the local search and global search capabilities of the PSO.In this study, the inertia weight (w) is 1.The constant γ is the speed limitation coefficient (© = 0.05).• Update the position of each particle using Equation (3): ( + 1) =   () +   ( + 1) (3) where   is the position of the i th particle.

Multiple Inversions with Genetic Algorithm (MI-GA)
John Holland was the first to use GA (Goldberg, 1989).It is a random search algorithm that was created by mimicking the processes of natural selection and genetics.GA works in string structures, which are like biological structures that develop over time by following the rules of "survival of the fittest" through the random exchange of information.In each generation, new strings are made by taking parts of the previous generation and putting them together (Roetzel et al., 2020).This algorithm is divided into four steps that employ random probabilities (La Hamimu, 2011;Rubaiyn et al., 2019) (2020).We develop multiple inversions using genetic algorithm (MI-GA) with the following procedure: 1. Determine the initial values of the model solution • Carry out mutations with probability Pm; 5.The model solution with the smallest misfit value obtained in the inversion is the global model solution (g) at t = kmax, with the misfit being   ; 6. Taking the average of the obtained model solutions as the final model solution.

Misfit and Similarity Index
The observed Rayleigh wave dispersion curves are inverted to determine the optimal combination of shear wave velocities (Vs) and the corresponding thicknesses (d) for each subsurface layer.This is done by minimizing the discrepancy between the measured and theoretical dispersion curves.The stiffness matrix method developed by Kausel and Rosset (1981) and Olafsdottir (2019) is used in the forward modeling method to calculate the theoretical Rayleigh wave dispersion curve.
The following is the objective function used to compute the misfit between the observed and theoretical dispersion curves: where    () and    () are the observed and calculated (or theoretical) Rayleigh wave phase velocities as a function of frequency (f), respectively, and N is the number of measured frequencies.
The accuracy of the algorithm used is checked through the similarity index (SI) value, which is expressed by the following equation: where    is the model parameter resulting from the inversion,    is the actual model parameter, and M is the number of model parameters.

RESULTS AND DISCUSSIONS
The MI-PSO and MI-GA were used to test the inversion process using two sets of synthetic data that mimicked the properties of subsurface layers that are often found in the field.The two subsurface layering models are classified into the low velocity layer (LVL) case and the complex case, which is a combination of low velocity layer (LVL) and high velocity layer (HVL).

LVL Case
The LVL case is characterized by the presence of a soft layer between the hard layers.The softness and hardness of a rock layer are represented by the value of the shear wave velocity (Vs) at a certain depth range (d).In this study, the soft rock layer is represented by a Vs value < 200 m/s.
The parameters of the earth's layer system that influence the Rayleigh wave phase velocity consist of primary wave velocity (Vp), shear wave velocity (Vs), layer thickness (d), and density (ρ).Of these four parameters, shear wave velocity (Vs) and layer thickness (d) are the two that most dominantly influence the Rayleigh wave phase velocity (Safani, 2007).Therefore, Vs and d are two parameters estimated in the inversion process, while the other two parameters are assumed to be known through a priori information.
Table 1 shows the model parameters and inversion search space for the LVL case.This model consists of five layers, including a half-space.The soft layer is in the second layer with a value of Vs = 150 m/s and is flanked by two harder layers with a value of Vs = 200 m/s for each layer.As an initial stage, these four parameters are used as input variables to calculate the Rayleigh wave phase velocity, which is then used as measured data in the inversion process.
The MI-PSO and MI-GA begin with determining the model search space, which includes the upper and lower limit values of shear wave velocity (Vs) and layer thickness (d) (Table 1).Both algorithms' inversion processes were repeated ten times, with 400 iterations for each inversion.The number of populations generated in the search space is 30.Inversion using these two global optimization methodologies not only assesses the accuracy of the inversion results produced by each approach, but also the efficiency of the time required to perform each inversion.This inversion process utilizes a laptop with the following specifications: an AMD Ryzen 3 4300U processor (2.7 GHz), 8192 MB of DDR4 memory, and Windows 10 HSL 64-bit.The results obtained for each inversion are shown in Table 2 for the MI-PSO and Table 3 for the MI-GA.The time required to complete each inversion using the MI-PSO varies, with the fastest time being 1330 s and the longest time being 1358 s.Meanwhile, the MI-GA completes the fastest calculation time of 1400 s, and the longest time is 1463 s.This shows that the time-consuming MI-PSO algorithm is more efficient than MI-GA.Inspection of the resulting misfit shows that the MI-PSO shows a smaller misfit in all inversions compared to the MI-GA.The small misfit in the MI-PSO is directly proportional to the accuracy of the resulting model solution.The average of the ten model solutions using the MI-PSO (green color) for the LVL case is more accurate than the model solution using the MI-GA (red color) (Figure 1). Figure 1 depicts the findings of the average model solution from these two methods, which include a comparison of the dispersion curves (Figure 1a) and a comparison of the average model profiles from each algorithm (Figure 1b).Table 4 details the average value of the model solution provided by the MI-PSO and MI-GA approaches, as well as the relative errors and similarity index.The relative error values in Table 4 explain the greater accuracy of the average model with the MI-PSO compared to the average model with the MI-GA, where the relative error for practically all parameters with the MI-PSO method is significantly smaller than with the MI-GA.Multiple inversions of Rayleigh wave dispersion curve for… Aside from that, the similarity index of MI-PSO method is higher than MI-GA, at 98.01% versus 91.95%.

Complex Case
The second model is a more complicated model that involves simultaneous presence of both LVL and HVL.Table 5 contains a tabular summary of the complex model and the search space for the Rayleigh wave dispersion data inversion procedure.This model profile consists of five layers, including half-space.The low velocity layer (LVL) is in the second and fourth layers, with shear wave velocity values of 150 m/s and 180 m/s, respectively.These two LVLs flank a HVL in the third layer with Vs = 240 m/s.In contrast to the Vs value, the Vp value increases with increasing layer depth, while the density (ρ) is set constant for all layers.Just like in the previous model, in this complex case, the four model parameters are used to generate the Rayleigh wave phase velocity, which is then used as "measurement data" in the inversion process.
The inversion process was carried out 10 times, with 400 iterations for each inversion.In this synthetic test, the search space parameters for Vs and d pairs for each layer are determined based on the true model values.However, in the case of field data, the boundaries of this search space are usually estimated from the corresponding observed phase velocity and wavelength.The number of model populations is the same as in the previous LVL model, namely 30 populations.
The solution model, time consumption, and misfit resulting from each inversion are presented in Table 6 for the MI-PSO and Table 7 for MI-GA.The minimum computing time spent to carry out one inversion process with the MI-PSO is 1254 s, and 1703 s is the maximum computing time.Meanwhile, in the MI-GA, the minimum computing time spent on one inversion process is 1515 s, and the maximum time is 2618 s.Misfit calculations for each single inversion using the MI-PSO and MI-GA show acceptable values.However, in general, inversion with MI-PSO provides smaller misfit values than MI-GA.The lowest misfit using the MI-PSO is 0.04%, and the largest misfit is 0.32%.Meanwhile, the MI-GA gives the lowest misfit of 0.76% and the largest misfit of 1.63%.This illustrates that the MI-PSO provides more accurate results than the MI-GA.So, even in this complex model, the MI-PSO algorithm still shows its superiority over MI-GA.The average model solution of all ten inversions of the complex case for both optimization methods is presented in Figure 2. The comparison of the observed dispersion curve and the calculated dispersion curve for the MI-PSO shows better agreement than that for the MI-GA (Figure 2a).This has implications for the final average model (i.e.Vs versus d profile), where the average model from MI-PSO (green thick line) matches the true model better than the average model from MI-GA (dashed red line) (Figure 2b).
The average model from the MI-PSO and MI-GA is presented in more detail in Table 8.The model parameter values (Vs and d) for each layer produced by the MI-PSO are closer to the true model and are supported by a smaller relative error compared to the MI-GA.The most striking discrepancy is seen in the thickness of the first and fourth layers estimated using the MI-GA, where the thickness of each layer is quite far from the true thicknesses.This discrepancy is confirmed by the relative errors of the two layers, which are quite high, namely 40.89% for the thickness of the first layer and 15.56% for the thickness of the second layer.The better SI value of the MI-PSO method (96.15%) compared to MI-GA (87.89%) again emphasizes the superiority of the MI-PSO.
Table 7. Model solutions for each inversion with the MI-GA for the complex case

CONCLUSION
The accuracy of the two developed algorithms, namely MI-PSO and MI-GA, was successfully tested using synthetic Rayleigh wave phase velocity data for two different cases.The test results show that both the MI-PSO and MI-GA algorithms are good at estimating the subsurface shear wave velocity (Vs) and layer thickness (d) for both LVL and complex cases.But MI-PSO outperforms MI-GA in solving the inverse problem of Rayleigh wave phase velocity, both in terms of computing time and the correctness of the final model solution.

Figure 1 .
Figure 1.Inversion results of Rayleigh wave dispersion curve for LVL case.

Figure
Figure 1b demonstrates how the two algorithms that produced shear wave velocity (Vs) versus depth (d) profiles can reconstruct the profile of the true model.However, the final model solution obtained with the MI-PSO provides better results compared to MI-GA.The LVL in the second layer can be estimated very well with the MI-PSO.
Evaluate the position in each iteration: if   () <   , update   =   , if   () < . , update g =   (). 5.The model solution obtained in the I-th inversion is the global model solution (g) at t = kmax, with the misfit being. .6. Calculate the average of the model solutions obtained as the final model solution.
1 …    ] and lower bound (  ) as   = [  1 …   +1   1 …   Pm = 0.01, learning rate η = 0.1, number of iterations kmax, and number of inversions Imax; 2. Loop the inversion I from 1 to Imax; 3. Initialize a random initial population with  1 =  0 and  2 =   + .* (  −   ),   =   + .* (  −   ), where rand is a random vector with the same dimensions.Then encode each population into a binary string; 4. Loop the iteration k from 1 to kmax; • Decode each binary string of particles into a decimal number, evaluating the misfit of each particle, and taking the value of the model solution with the minimum misfit; • Converts misfit function values to fit values; • Reproduce the generation model solution around the best model solution, and then encode; • Randomly pair binary strings with row swaps; • Perform crossovers with probability Pc; ], population size M, number of bits representing each variable   = 10, crossover probability Pc = 1, mutation probability

Table 1 .
LVL model parameters and inversion search space

Table 2 .
Model solutions for each inversion with the MI-PSO for the LVL case

Table 3 .
Model solutions for each inversion with the MI-GA for the LVL case

Table 4 .
Average model with MI-PSO and MI-GA for the LVL case

Table 5 .
Model parameters and inversion search space for the complex case

Table 6 .
Model solutions for each inversion with the MI-PSO for the complex case