The joint application of a metaheuristic algorithm and a Bayesian statistics approach for uncertainty and stability assessment of nonlinear magnetotelluric data
In this paper, we have developed three algorithms, namely hybrid weighted particle swarm optimization (wPSO) with the gravitational search algorithm (GSA), known as wPSOGSA; GSA; and PSO in MATLAB to interpret one-dimensional magnetotelluric (MT) data for some corrupted and non-corrupted synthetic data, as well as two examples of MT field data over different geological terrains: (i) geothermally rich area, island of Milos, Greece, and (ii) southern Scotland due to the occurrence of a significantly high electrical conductivity anomaly under crust and upper mantle, extending from the Midland Valley across the Southern Uplands into northern England. Even though the fact that many models provide a good fit in a large predefined search space, specific models do not fit well. As a result, we used a Bayesian statistical technique to construct and assess the posterior probability density function (PDF) rather than picking the global model based on the lowest misfit error. The study proceeds using a 68.27 % confidence interval for selecting a region where the PDF is more prevalent to estimate the mean model which is more accurate and close to the true model. For illustration, correlation matrices show a significant relationship among layer parameters. The findings indicate that wPSOGSA is less sensitive to model parameters and produces more stable and reliable results with the least uncertainty in the model, compatible with existing borehole samples. Furthermore, the present methods resolve two additional geologically significant layers, one highly conductive (less than 1.0 Ωm) and another resistive (300.0 Ωm), over the island of Milos, Greece, characterized by alluvium and volcanic deposits, respectively, as corroborated by borehole stratigraphy.
The magnetotelluric (MT) method is a natural source electromagnetic method that explores various natural resources, namely hydrocarbon, minerals, geothermal prospects, groundwater, and metalliferous ores, etc. (Nabighian and Asten, 2002; Simpson and Bahr, 2005). Due to its instability, non-unique solution, and algorithm sensitivity, MT data interpretation is thought-provoking. Many researchers have attempted and developed various inversion algorithms to interpret and improve the model accuracy, convergence speed, and stability and reduce the uncertainty of the solutions (Kirkpatrick et al., 1983; Constable et al., 1987; Rodi and Mackie, 2001; Li et al., 2018; Zhang et al., 2019; Khishe and Mosavi, 2020). There are mainly two categories of the inversion algorithm: first, the local optimization methods, namely the conjugate gradient, Levenberg–Marquardt/ridge regression, the Gauss-Newton algorithm, the steepest descent, and Occam's inversion, which require a good initial guess (Shaw and Srivastava, 2007; Wen et al., 2019; Roy and Kumar, 2021). Another category is global optimization techniques (i.e., ant colony optimization, genetic algorithm, particle swarm optimization, gravitational search algorithm, and simulated annealing, etc.), which do not require an initial guess. Many researchers have carried out numerous metaheuristic optimization algorithms to invert MT data (Dosso and Oldenburg, 1991; Pérez-Flores and Schultz, 2002; Miecznik et al., 2003; Sen and Stoffa, 2013). These algorithms are inspired by natural phenomena and have various geophysical applications, including particle swarm optimization (Kennedy and Eberhart, 1995; Essa et al., 2023), the genetic algorithm (Whitley, 1994), Bat algorithm (Yang, 2010a; Essa and Diab, 2023), differential evolution (Storn and Price, 1997), biogeographically based optimization (Simon, 2008), the firefly algorithm (Yang, 2010b), grey wolf optimizer (Mirjalili et al., 2014), ant colony optimization (Colorni et al., 1991), the gravitational search algorithm (Rashedi et al., 2009), and the novel barnacles mating optimizer algorithm (Ai et al., 2022).
However, unique characteristics, namely exploration and exploitation, persist in any global optimization algorithm. For example, the particle swarm optimization (PSO) algorithm has a very high potential for exploitation, which implies that the algorithm performs well in a local search but is inferior in exploration (S̨enel et al., 2019). This suggests that the algorithm has a limited capacity to estimate the best model in an extensive search range. Because of low exploration characteristics, it gets trapped at the local minima (Mirjalili and Hashim, 2010). So, integrating the two algorithms with opposite characteristics is the best way to balance exploration and exploitation characteristics to achieve better solutions than the results obtained from an individual algorithm.
Here, we utilized weighted hybrid PSOGSA (wPSOGSA), a new global optimization method that takes into account the algorithm based on natural behavior seen in birds, fish, and insects, known as PSO, and gravity-based Newton's law (with high exploration capability), known as the gravity search algorithm (GSA). Researchers interested in artificial intelligence and developing effective optimization algorithms for comparative analysis of different metaheuristic algorithms (Pace et al., 2022) have been drawn to notable characteristics in such social behavior. The wPSOGSA, PSO, and GSA are used to estimate the resistivity distribution of 1D multi-layered Earth models using synthetic (noise-free and noisy) data for three- and four-layer cases, taken from Shaw and Srivastava (2007) and Xiong et al. (2018), respectively, and field MT sounding data for six- and four-layer cases, taken from Hutton et al. (1989) and Jones and Hutton (1979), respectively.
Furthermore, numerous (here 10 000) models that fit well are optimized to obtain the mean model, which is preceded by calculating the posterior probability density function (PDF) based on Bayesian concepts using all accepted models to find the optimal mean solution with the least uncertainty, as well as a correlation matrix to determine the relationships among the layer parameters. Thus, the research reveals that the wPSOGSA method may be utilized to provide a more accurate and reliable model with superior stability, a quick rate of convergence, and the least amount of model uncertainty.
2.1 Synthetic and field data
Different MT datasets are utilized to evaluate the proposed wPSOGSA's effectiveness, sensitivity, stability, and robustness in outlining the genuine subsurface structure. These datasets are noise-free and Gaussian noise synthetic data produced for several geological formations, and two MT field data have been optimized for analysis.
To demonstrate and evaluate the robustness of the present algorithms, we have generated a synthetic MT apparent resistivity and apparent phase dataset without noise and with noise levels (10 % and 20 % noise) considering a three-layer typical continental crustal model with a total thickness of 33 000 m (i.e., 33.0 km) and a resistivity of the upper crust of 30 000.0 Ωm with 15 000 m (i.e., 15 km) thickness (high resistive layer) and a resistivity of the middle crust of 5000.0 Ωm with 18 000 m (i.e., 18.0 km) thickness (reasonable low resistive layer) underlain by 1000 Ωm (low resistive) half space taken from Shaw and Srivastava (2007).
For the second example of the synthetic data, a typical four-layer HK-type of Earth model taken from Xiong et al. (2018) is generated by forward modeling equations for the demonstration of the wPSOGSA, PSO, and GSA, and their performance is compared with improved differential evolution (IDE) results obtained by Xiong et al. (2018).
We utilized the first example of field data taken from Hutton et al. (1989) on the island of Milos, Greece. Milos is a part of the South Aegean active volcanic arc, an example of an emergent volcanic edifice (Stewart and McPhie, 2006) formed by monogenetic effusive and explosive magmatism pulses. Milos is the world's biggest exporter of bentonite, and it also has a diverse variety of metalliferous and non-metalliferous mineral reserves. It is a conserved on-land laboratory for studying shallow underwater hydrothermal ore-forming processes. The accompanying shallow subsurface hydrothermal venting fields have had significantly less attention. Dawes (1986) used magnetotelluric data to assess the resistivity structure of the geothermal area on Milos' west side. With around 3.0 km spacing, 37 MT probes in the bandwidth of 100–0.01 Hz and 12 investigations in the bandwidth of 0.01–0.0001 Hz were installed, along with various profiles that were perpendicular to the Zephyria graben in the W–E direction, as well as along the graben in the S–N direction (Hutton et al., 1989).
Another field example of MT data from Jones and Hutton (1979) was picked to illustrate our technique from Newcastleton (55.196∘ N, 2.796∘ W, in geographic coordinates), Southern Uplands of Scotland. By the Southern Uplands fault, the Southern Uplands are isolated from the Midland Valley. The bulk of the Southern Uplands comprises Silurian–lower Paleozoic sedimentary deposits such as greywackes and shales that originated in the Iapetus Ocean during the late Neoproterozoic and early Paleozoic geologic eras. These rocks emerged from the seafloor as an accretionary wedge during the Caledonian orogeny. The majority of the rocks are coarse greywacke, a kind of sandstone that has been poorly metamorphosed and contains angular quartz, feldspar, and small rock fragments. The Midland Valley and northern England, on the other hand, are known for their thick Carboniferous layers, which are used to measure coal.
2.2 Forward modeling – magnetotelluric formulation for 1-D Earth
The ability to formulate an effective inversion method requires a thorough understanding of the forward modeling technique for the issue of interest. Factors like frequency range, actual resistivity, and layer thickness are used to create datasets of synthetic MT apparent resistivity, ρa(ω), and apparent phase, φa(ω). The electromagnetic impedance (Z) for layered structures is described in terms of an orthogonal electric field, magnetic field, wavenumber (k), reflection coefficient (R), and exponent factor (τf) with angular frequency (ω) as (Ward and Hohmann, 1988)
where the wavenumber , Ex and Ey are components of the electric field, and Hx and Hy are the components of the magnetic field.
If displacement currents are not taken into account, Eq. (1) becomes
Noisy impedance is calculated by the following equation:
If the angle between impedance phase with Ex is 45∘, then the resistivity (ρ) in half space of impedance Z(ω) and time period (T) can be written as
Thus, the apparent resistivity and apparent phase are defined (Cagniard, 1953; Ward and Hohmann, 1988) as follows:
where the exponent factor , the induction parameter , h is the layer thickness, μ0 is the magnetic permeability for free space, Z∗ is the complex conjugate of impedance, and “rand” is used for generating a random number between 0 and +1.
2.3 Global optimization technique
The techniques that we have used for the joint modeling of metaheuristic global optimization, namely PSO, GSA, and wPSOGSA in Step 1 and posterior Bayesian probability density function technique in Step 2, obtain the global model by utilizing the synthetic data generated by using forward modeling and field MT apparent resistivity and phase curves, depicted in the schematic diagram (Fig. 1).
2.4 Optimization and error estimation
In the present study, we have implemented a new innovative global optimization technique known as wPSOGSA, in which swarm particles and mass particles provide the best particle, i.e., the best model. The best model is chosen based on the fitness of the particles, and the cost function or objective function is used to estimate this fitness. Thus the magnetotelluric (MT) inverse problem can be formulated through the forward modeling operator, f(x), to achieve the resistivity model, which illuminates the observed data ρ and φ. This operator combines the problem of physics and inverts the observed apparent resistivity and phase data to the resistivity–depth model, x, as
The cost function (fitness of the particle) is a mathematical relation between observed and calculated data, and the root mean square error (RMSE) is defined as follows:
where N is the total observed data points, ρ and φ are the observed apparent resistivity and phase, and ρC and φC are the computed apparent resistivity and phase data.
2.5 Particle swarm optimization
The particle swarm optimization (PSO) technique is a widespread evolutionary optimization approach for determining the optimal global solution to a nonlinear inverse problem (Kennedy and Eberhart, 1995). This technique is analogous to the particle's natural behavior in search of food with the help of collaborative support from the population represented by geophysical resistivity solutions/models (known as particles) in a swarming group. The best model/position obtained among the particles so far is stored for each iteration, which helps in search of the global best solution, defined by the fitness of each particle estimated using Eq. (8). The velocity and location of the kth particle at the tth iteration are defined as
where w is the inertia weight set between 0 and 1; c1 and c2 are a personal learning coefficient and a global learning coefficient, respectively; vk(t) is the velocity of the kth particle at tth iteration; rand is used for a random number between 0 and 1; xp is the present best solution, xg is the global best solution, and xk(t) is the position of the kth particle at the tth iteration. Particles change their position at each iteration to approach an optimum solution. The first, second, and third terms in Eq. (9) represent exploratory ability, private thought, and particle collaboration, respectively.
2.6 Gravitational search algorithm
The gravitational search algorithm (GSA) is a metaheuristic algorithm based on Newton's gravitational law (Rashedi et al., 2009), which states that mass particles attract each other with a gravitational force that is directly proportional to the product of their masses and inversely proportional to the square of the distance between them. It signifies that massive particles (here, particle represents the resistivity layer model/solution) attract the neighboring lighter particles. Similar to PSO, the gravitational search optimizer works with a population of particles known as mass particles in the universe. Thus the best model, solution, or particle is achieved among the mass particles. The best model is defined by each particle's capability (i.e., the fitness), calculated using Eq. (8). The initialization of their position in the search spaces is given by
where N and D are the number of particles and models, i.e., the dimension of the model, and “up”, and “down” are the upper and lower limit of the search range, respectively.
During execution time, the gravitational acting force on agent kth from agent jth at a specific time (t) is defined as
where Ma,j and Mp,k are the active and passive gravitational masses for particle j and k, respectively; xj(t) is the position of the particle j at a time t for various parameters; Rk,j (t) is the Euclidian distance between two particles; and ε is a small constant.
Here, the gravitational constant G(t) at a specific time t is defined as in Kunche et al. (2015), and the acceleration of the kth agent at the tth iteration for models is ack(t), defined as
where Fk,j(t) the gravitational acting force on agent k from agent j, and Mk(t) is the mass of the object at a specific time (t).
where α, G0, “iter”, and “maxiter” are descending coefficients, starting value of gravitational constant, current iteration, and maximum iterations, respectively.
The following equations are used to update the particle's velocity and location:
All the particles are randomly placed in the search range using Eq. (11), and then the particle's velocity is initialized. Meanwhile, the gravitational constant, total forces, and acceleration are computed, and the locations are updated. The end criterion, the misfit error (i.e. 10−9), is taken in our study.
2.7 Weighted hybrid PSOGSA (wPSOGSA)
The weighted hybrid of PSO and GSA, known as wPSOGSA, integrates two essential characteristics, exploration (i.e., the ability of an algorithm to search the whole range of a given parameter) and exploitation (i.e., the ability to converge the solution nearest to the best solution) of the global optimization algorithm that increases its efficiency and converges the objective function to achieve global minima. The velocity and location of the particles updated in the wPSOGSA are illustrated in the schematic diagram (see Fig. 2).
The wPSOGSA combines the characteristic of social thinking of PSO and the searching capability of GSA; thus, the particle's velocity is defined as
where vk(t) is the velocity of the particle k at iteration t, w is the weight function (i.e., the constant which helps to control the momentum of the algorithm to perform optimization properly), ack(t) is the acceleration of agent k, xg is the best solution, and rand is a random number that lies between 0 and 1. At each iteration, particles updated their location to achieve the best solution defined as
The algorithm starts by randomly initializing the velocity, mass, and acceleration of the particles. The cost function is evaluated for all particles for specified iterations to get the most optimal solution, and inverted results are updated at each iteration. Equations (12), (17), and (18) are used to update the gravitational force, velocity, and location of particles after initialization. However, the velocity and position stop updating their values when the algorithm converge, and the least error of the cost function is reached.
2.8 Bayesian probability density function
In a Bayesian framework, the probability distribution of the model parameters (known as posterior probability distribution) is computed using given observed data and models obtained from inversion. The posterior for a model is calculated using Bayes' theorem and previous model space information. Individual model parameter ranges are incorporated into the prior knowledge. The two fundamental stages in the Bayesian statistics method are the representation of previous knowledge as a probability density function and calculation of the likelihood functional derived from data misfit (Tarantola and Valette, 1982). Specific characteristics, such as the best fitting model, mean model, and correlation matrix, may be determined from posterior distribution of models. According to Bayes' theorem,
As a result, our priori distribution function (f(μ)) for the parameter, xu; mean priori information, M; and t2, the mean uncertainty, is defined as
and the likelihood function is
Hence the posterior density function calculated for a parameter (xu) using mean (μ) and variance (σ2) defined (Lynch, 2007) as
The posterior Bayesian PDF is calculated from accepted models within a set of parameters, as shown below:
where P(X|E) is the posterior probability distribution of the parameter (X) given the evidence (E), P(X) is the prior information of (X), and L(E|X) is the likelihood function of X.
After the application of the PDF, the study further proceeds by choosing a confidence interval (CI) of 68.27 % that is based on the empirical rule, known as the 68–95–99.7 rule (Ross, 2009). The model parameters below 68.27 % CI are discarded, and the remaining parameters are used for determining the mean model and uncertainty. Thus, the mean model (Pj) is calculated using the best models having a PDF within a 68.27 % CI, defined in the following equation:
Here accepted models are used to calculate the correlation matrix (i.e., correlation among model parameters lie between −1 and 1) using the following equation (Tarantola, 2005):
Here, N is the total number of models; d is used for the number of the layer parameters; and Pj,k is the jth model parameter of the kth model, where l and j both vary from 1 to d (number of layer parameters). CovP(l,j) is the covariance matrix between model parameters l and j, Pl,k is the lth model parameter of the kth model, and CorP(l,j) is the correlation matrix between model parameters l and j.
The effectiveness, sensitivity, stability, and robustness of the proposed wPSOGSA in identifying the authentic subsurface structure are evaluated using various MT datasets. These datasets consist of synthetic data with no noise and Gaussian noise, which simulate different geological formations. Additionally, two MT field datasets have been optimized for analysis.
3.1 Application to synthetic MT data – three-layer case
The first example of synthetic MT data was executed for 10 000 runs, keeping the same lower and upper bounds as given in Table 1 and iterations to 1000. Figure 3 shows (a) the observed apparent resistivity with the computed data; (b) the observed apparent phase with the computed data; (c) the 1D inverted model by wPSOGSA (red color), GSA (green color), and PSO (blue color) with a true model (black color); and (d) the relation between the misfit and iterations for the noise-free synthetic data.
The misfit curve as shown in Fig. 3d is gradually decreasing with increasing iterations and becomes constant where the algorithm converges. The PSO, GSA, and wPSOGSA converge at iterations 492, 35, and 316, with associated errors of , , and and associated computational times of 27.06, 1.75, and 3.35 s, respectively. Thus, the curves show that wPSOGSA converges at the least RMSE. In contrast, PSO, GSA, and wPSOGSA using 10 % noisy synthetic data converge at 102, 88, and 358 iterations, with associated errors of 0.00435, 0.00439, and 0.00426 and associated computational times of 5.61, 4.40, and 3.80 s, respectively.
The 20 % noisy synthetic MT data were executed for 10 000 runs, keeping the same lower bound, upper bound, and iterations. The well-fitted inverted MT response (see Fig. 4) was as follows: (a) the corrupted synthetic and calculated apparent resistivity data, (b) the corrupted synthetic and calculated apparent phase data, (c) the inverted 1D depth model, and (d) convergence response in terms of misfit error versus iterations. We analyzed Fig. 4d and found that the PSO, GSA, and wPSOGSA using noisy synthetic data converge at iterations 236, 7, and 73, with associated errors of 0.0394, 0.0408, and 0.0393, respectively.
3.1.1 Bayesian analysis and uncertainty in model parameters
Two methods are used to estimate a mean solution and uncertainty: one method is the mean solution for all accepted best-fitted solutions acquired from 10 000 runs for all three global optimization techniques; another method is the model derived from all approved solutions using the posterior Bayesian PDF within 1 standard deviation. To get the global best solutions in our study, we incorporated the posterior PDF based on the Bayesian approach to enhance the efficacy of the inverted model and minimize the uncertainty in the model. The process for obtaining the mean solution proceeds by selecting an initial threshold error, which is essential because the smaller the threshold value, the more significant the number of models with less uncertainty in the model parameters (Sharma, 2012). Thus, a more considerable threshold gives a lower number of models with enormous uncertainty in the model parameter (Sen and Stoffa, 1996; Sharma, 2012). The study proceeds by calculating the PDF for each parameter value using Eq. (22). In order to select values of each parameter that having a higher posterior PDF, a 68.27 % CI is used. The mean model obtained from selected model parameters is near to the actual model.
Figure 5 shows the output of the posterior Bayesian PDF, which selects model parameters with less error. The straight lines (dashed lines) present the actual value of the respective layer parameters. The first-layer thickness, second-layer thickness, and first-layer resistivity have higher uncertainties, i.e., 61.25 m, 51.47 m, and 210.61 Ωm, respectively, whereas the second-layer resistivity and third-layer resistivity have lower uncertainty, i.e., 17.71 and 0.03 Ωm, respectively.
Table 1 shows the inverted layer parameters using wPSOGSA, GSA, and PSO for noise-free and noisy synthetic MT data based on the posterior Bayesian PDF, as well as the actual model and the search range. In addition, layered properties of synthetic data corrupted with 10 % and 20 % random noise are compared and statistically analyzed. Our findings, as shown in Table 1, were compared to those obtained using the genetic algorithm (GA), ridge regression (RR), and PSO by Shaw and Srivastava (2007), which consistently outperform GA and RR, which is closer to the genuine model.
The mean values of the accepted model parameters are 30 243.42±471.26, 5007.04±39.59, 1000.02±0.064, 14 969.33±136.82, and 18 029.76±114.90, and their associated amounts of uncertainty are 1.5 %, 0.78 %, 0.0064 %, 0.91 %, and 0.63 %. On the basis of a low posterior PDF and high uncertainty, we have taken (ρ1) and (h1) for the exercise to show the models are not biased to the selected models.
As well as being based on the histograms (see Fig. 6), the posterior PDF and uncertainty of the inverted layer parameters resistivity (ρ1) and thickness (h1) for the three-layered synthetic MT data have been taken to depict the global solution using presented algorithm. Here we prepared the cross-plots of ρ1 versus h1 using (a) wPSOGSA, (b) PSO, and (c) GSA, showing all accepted models (red circle), selected models with a misfit error of less than a threshold error of 10−4 (magenta circle), models of a PDF greater than 95 % (blue circle), models of a PDF greater than 75 % (green circle), models of a PDF greater than 68.27 % (yellow circle), and mean models, i.e., model parameters which have a PDF greater than 68.27 % (black asterisk), as shown in Fig. 7. It is noticed that all inverted results give the global solution, which has a good agreement with the true model, whereas wPSOGSA gives more accurate results than the other two algorithms, PSO and GSA, as shown in Table 2.
3.1.2 Sensitivity, correlation matrix, and model parameters
The accepted models, which have a posterior PDF value within 68.27 % CI, are used to calculate the correlation matrix. This correlation matrix gives the relationship among model parameters. Thus, the lower correlation value gives weak relation among the parameters and vice versa. The correlation matrix of PSO, GSA, and wPSOGSA was examined on one set of synthetic data, as shown in Figs. 8, 9, and 10, demonstrating the sensitivity among inverted model parameters. The value of correlation matrix, 1.0, indicates that the two parameters are strongly correlated.
Figure 8 shows that first-layer resistivity is correlated highly positively with a first-layer thickness (0.97) and second-layer thickness (0.98), while the second-layer resistivity (−0.99) and third-layer resistivity (−0.81) are substantially negative connected. Second-layer resistivity is correlated with the third-layer resistivity (0.87), which has a significant positive relationship, while second-layer resistivity has a significant negative correlation with the first-layer thickness (−0.99) and the second-layer thickness (−1.00). First-layer thickness (−0.92) and second-layer thickness (−0.90) are very negatively associated with third-layer resistivity, while first-layer thickness is extremely positively correlated with second-layer thickness (0.99).
Figure 9 indicates that first-layer resistivity is highly associated with second-layer thickness (1.00) and weakly with second-layer resistivity (−1.00), third-layer resistivity (−1.00), and first-layer thickness (−1.00). Second-layer resistivity (−1) is highly linked with a second-layer thickness (−1.00), while third-layer resistivity (1.00) and first-layer thickness are strongly correlated (1.00). Third-layer resistivity has a highly positive correlation with first-layer thickness (1.00) and a strong negative correlation with second-layer thickness (−1.00), whereas first-layer thickness has a significant negative correlation with second-layer thickness (−1.00).
Figure 10 shows the correlation matrix of wPSOGSA. The analyses reveal that the first-layer resistivity is strongly negative with the second-layer resistivity, substantially negative (−0.92) with the third-layer resistivity, weakly positive (0.30) with the first-layer thickness, and considerably (0.63) with the second-layer thickness. Second-layer resistivity is slightly positive (0.31) when compared to third-layer resistivity (0.43) but substantially negative when compared to first-layer thickness. Third-layer resistivity has a slightly negative correlation (−0.23) with first-layer thickness but a moderately negative correlation (−0.71) with second-layer thickness, and first-layer thickness has a negative correlation (−0.71). Thus, the conclusion can be made that the layer parameters are independent of others, so changing one will have no effect on the other, compared to the result obtained via PSO and GSA.
3.1.3 Stability analysis
We used two different search ranges for stability evaluation of the proposed wPSOGSA and executed the algorithms over three layers of synthetic MT data, one of which is expanded, and the other is contracted by 10 % of the initial search range. We infer from three layers of synthetic data that results fluctuate by approximately 3 % from the true value when the search range is changed. This variation is about 10 % on average for synthetic data corrupted with 30 % random noise, as shown in Table 2.
3.2 Application to synthetic MT data – four-layer case
For the second example of the synthetic data, a typical four-layer HK-type of Earth model to analyze the performance of the present algorithm with improved differential evolution (IDE) results obtained by Xiong et al. (2018). Analysis over noisy synthetic data is done by corrupting synthetic data with 10 % and 20 % Gaussian random noise to mimic the real field data because different types of noises influence apparent resistivity data. Following that, all three optimization methods are run using the noisy synthetic data. As the misfit error increases with the noise in the data, the Bayesian PDF of 68.27 % CI is calculated with respect to the threshold misfit error of 0.01, and thus the mean model is calculated.
Enormous uncertainty is shown in the inverted results; hence, we calculated the mean model for 68.27 % CI using the posterior Bayesian PDF to reduce the uncertainty and produce the global best solution. The optimized results obtained from the posterior PDF and the true model are shown in Table 3. Figure 11 illustrates that the inverted responses for PSO, GSA, and wPSOGSA are well-fitting as follows (a) observed and calculated apparent resistivity data, (b) observed and calculated apparent phase data, (c) 1-D depth model, and (d) convergence response of present algorithms. We have estimated the layer parameters for synthetic data corrupted with 20 % random noise for comparative analysis and found that PSO, GSA, and wPSOGSA converge at iterations 96, 556, and 187, with associated errors of 3.69, 4.04, and 3.69, respectively.
Additionally, the synthetic data corrupted with 10 % random noise are also used for the execution of inversion, keeping the search range, a number of particles, and iterations the same as before; it was observed that PSO, GSA, and wPSOGSA converge at iterations 151, 2, and 250, with associated errors of 1.7609, 1.95, and 1.76, respectively. The posterior Bayesian PDF for threshold data with 68.27 % CI is calculated similarly to a three-layer case to minimize the uncertainty in inverted results.
For the stability evaluation of presented algorithms over four layers of synthetic MT data, similar to the three-layer case, we used two different search ranges and executed the algorithms for 1000 iterations. The method exhibits good results with four layers of synthetic data and reveals minimal variation for noise-free data. For 30 % contaminated data, the variation is approximately 10 % and 12 % in Case 1 and Case 2, respectively. The outputs do not change much across runs and provide consistent results, as shown in Table 4.
3.3 Application to field MT data – island of Milos, Greece
In one-dimensional MT data for site G5 near borehole M2 (Hutton et al., 1989), as shown in Fig. 12, the apparent resistivity and phase values are inverted using wPSOGSA, PSO, and GSA, keeping the same set of controlling parameters as for noisy synthetic data, such as the swarm size, inertia weight (w), personal learning coefficient (c1), global learning coefficient (c2), descending coefficient (α), and initial value of the universal gravitational constant (G0).
Figure 13 shows the calculated data and model parameters as (a) a match between observed and computed apparent resistivity data and (b) a match between observed and computed apparent phase data; (c) the 1D inverted model; and (d) the convergence response of wPSOGSA (red color), GSA (green color), and PSO (blue color), along with the true model (black color). Figure 13c depicts alluvium deposits with a resistivity of 1.0 Ωm with 15 m thickness as the top layer, and volcanic deposits with a resistivity of 300 Ωm and 10 m thickness lie beneath the alluvium deposits. A very high conducting layer of resistivity less than 1.0 Ωm is estimated, equivalent to the green lahar under the high resistivity volcanic deposits. The next layer below, with higher resistivity, corresponds to the crystalline foundation. In the geothermal zone's depths, the resistivity drops again. The resistivity in the depth range of about 1000 m, which is similar to earlier studies, was explored, and the findings of the proposed algorithm were discovered to be in good agreement with the model developed by Dawes in Hutton et al. (1989).
Figure 13d reveals that the algorithms converge at iterations 218, 1, and 425, with corresponding errors of 0.0494, 0.0518, and 0.0493 for PSO, GSA, and wPSOGSA, respectively. The hybrid algorithm has the least error between observed and computed data. The algorithms are executed for 1000 iterations and 10 000 models, and findings are compared with available stratigraphy. The result is derived using the Monte Carlo technique by Hutton et al. (1989). After examining our optimized effects from Fig. 13 and Table 5, hybrid wPSOGSA outperformed PSO and GSA.
3.3.1 Bayesian analysis and uncertainty in model parameters
A posterior Bayesian method determines the global model and related uncertainty. Figure 14 shows another uncertainty study that examined the six-layered resistivity model over the geothermal field, island of Milos, Greece, and found that the peak values of the posterior PDF for all model parameters are much closer to the actual value of the layer parameters, providing less uncertainty. We have analyzed the wPSOGSA inverted results from Fig. 14 and Table 5 and found that the first, second, third, fourth, fifth, and sixth layers' resistivity with uncertainty in associated layer parameters is 1.23±0.49, 297.61±53.43, 0.55±0.02, 2.41±0.16, 14.18±1.76, and 99.92±0.37 Ωm. Similarly, the associated thicknesses with uncertainty are 14.51±1.35, 9.85±1.35, 127.39±6.01, 823.01±7.57, and 2750.88±63.07 m. Thus, the analysis suggests the lower uncertainties in each layer's parameters except resistivity of the first and second layers.
Table 5 compares optimized results obtained from all three presented algorithms based on the posterior Bayesian PDF under 68.27 % CI condition. However, the 1D depth model inverted from wPSOGSA shows good agreement with the available borehole M-2 (Hutton et al., 1989). As a result, the hybrid algorithm is functioning better, and the findings are encouraging.
3.3.2 Sensitivity, correlation matrix, and model parameters
Here, a similar study of the correlation matrix is carried out for field example from the island of Milos, Greece, using all accepted models, which have posterior PDF values within 68.27 % CI. The correlation matrix of PSO, GSA, and wPSOGSA was examined over the field MT data as shown in Figs. 15, 16, and 17, demonstrating the sensitivity among inverted model parameters, and found an almost similar correlation among the layer parameters for three-layer synthetic study. From correlation analyses, we noticed that the values are showing moderate and weak correlation among parameters in the wPSOGSA case, indicating that wPSOGSA is linearly independent of layer parameters. This indicates that the parameter is less affected by other layer parameters and resistivity curves, whereas the correlation among layer parameters for field data using GSA and PSO is either strongly positive or strongly negative, which shows that the parameters are dependent on each other. Thus a change in one parameter affects the other, and the apparent resistivity curve is also very much involved.
3.4 Application to field MT data – Newcastleton, Southern Uplands, Scotland
Another field example of MT data was picked to illustrate our technique from Newcastleton (55.196∘ N, 2.796∘ W, in geographic coordinates), Southern Uplands of Scotland. The Southern Uplands are isolated from the Midland Valley by the Southern Uplands fault. The location of the MT site and the geology of the study area are shown in Fig. 18.
During 9 d, in the frequency range of 0.1 to 0.0001 Hz, the variations of the magnetic and telluric fields concerning the time at four sites along a line perpendicular to the anomaly's strike were recorded, keeping a high signal-to-noise ratio, where the anisotropy ratios were very near to 1, and the skew factor is less than 0.1 for the majority of periods. Due to low anisotropy ratios and the skew factor, the resistivity distribution under this location is one-dimensional (Jones and Hutton, 1979). Here one set of MT data is inverted using PSO, GSA, and wPSOGSA to obtain the best fitting apparent resistivity curve, apparent phase curve, and 1D depth model, as shown in Fig. 19a, b, and c, respectively. Figure 19 shows a realistic one-dimensional resistivity variation, with a phase response ranging from 60∘ at 100 s to 35∘ at 1000 s, which can only be obtained by establishing a conducting zone at lower crustal–upper mantle levels (Jones and Hutton, 1979).
The execution time for wPSOGSA (33 s) is the lowest as compared to GSA (34 s) and PSO (53 s). The convergence iterations of PSO, GSA, and wPSOGSA are 79, 101, and 65, and its associated misfit errors are 3.79, 4.72, and 3.70, respectively.
The inverted MT model is illustrated in Fig. 19c, which depicts two low conductive zones at a depth of 21 and 400 km. The first conductive layer (70 Ωm) with a thickness of 28 km is underlain by a high resistive top layer of thickness of 21 km, and the second very high conductive layer (less than 1.0 Ωm) at a depth of 400 km is underlain by a high resistive layer (550 Ωm) of thickness 351 km. Thus, the last layer of a very high conductive zone (i.e., resistivity less than 1.0 Ωm) as a lower crust–upper mantle conductor at a depth of 400 km is estimated. At 400 m depths, a conducting zone meets both the amplitude and phase long period responses. This explanation is directly equivalent to accepted models derived from Monte Carlo models for the structure underlying the Southern Uplands.
The analysis on the two synthetic MT datasets shows that the proposed algorithm works very well and provides encouraging well-fitted calculated apparent resistivity and apparent phase data with the associated observed data. Also from the study, it is noted that the proposed algorithm is less sensitive to the search range and the constraints used in this algorithm. And the comparison of the model parameters shows that the output from the proposed algorithm is very precise in comparison to the true model and faster than the individual algorithms and the other algorithms used in the previous papers.
A research study conducted by Hutton et al. (1989) compared various techniques, including Parker H+, Dawes (a combined Monte Carlo–Hedgehog approach developed by Dawes), Jupp–Vozoff, Fischer, and Parker D+ inversions, to interpret the subsurface geology of the island of Milos, Greece. Among these techniques, the Dawes algorithm proved to be the most effective in identifying a reservoir interface at the borehole's depth. Consequently, these methods were utilized to compile 1-D models along several traverses. All three models, Parker H+, Parker D+, and Dawes, indicate the presence of a resistivity boundary at a depth of approximately 1000 m (Hutton et al., 1989), where the geothermal reservoir has been detected. The research findings through the proposed algorithm reveal that the resistivity of the crystalline basement beneath the geothermal site is abnormally low (<20 Ωm) in the uppermost portion and remains below 50 Ωm at depths of at least 10 km.
Figure 19c presents the inverted MT model of the Southern Uplands, Scotland, displaying two zones of low conductivity at depths of 21 and 400 km. Contrary to the findings of Jain and Wilson (1967), there is strong evidence suggesting that the conducting zone (70 Ωm) beneath the Southern Uplands exists at a depth exceeding 20 km. Furthermore, there is a second layer of extremely high conductivity (representing lower crust–upper mantle of less than 1.0 Ωm) at a depth of 400 km, which is underlain by a highly resistive layer (550 Ωm) spanning a thickness of 351 km. At a depth of 400 m, this conducting zone shows a direct alignment with both the amplitude and phase responses of long period measurements, which can be noticed in the model derived from Monte Carlo simulations of the structure underlying the Southern Uplands (see Fig. 19). The results from the proposed algorithm as well as from PSO and GSA (see Fig. 19) demonstrate the presence of a highly conductive layer at depths exceeding 20 and 400 km, corroborating the findings of Jones and Hutton (1979).
The study presented wPSOGSA along with PSO and GSA and evaluated their efficacy and applicability to the MT data. These algorithms narrate the appraisal of 1D resistivity models from apparent resistivity, apparent phase, and frequency datasets. So, synthetic and field MT data from various geological terrains were used to demonstrate the relevance of these methods, which are further carried out by applying multiple runs, generating a large number of models that fit the apparent resistivity and apparent phase curves. Furthermore, these best-fitting models within a specified range are then chosen for statistical analysis. The statistical analysis includes the posterior PDF based on the Bayesian approach with a 68.27 % CI, the correlation matrix, and stability analysis in order to understand the accuracy of the mean model and its uncertainty. However, the solution from the posterior PDF based on the Bayesian of wPSOGSA is better than GSA and PSO, explaining the reliability of the proposed inversion algorithm. In general, conventional techniques can effectively resolve the model in random noise, but they can miscarry in methodical error or inappropriate models. Also, the performance of the proposed algorithms on field datasets has been analyzed based on the mean model, uncertainty, correlation, and stability of layered Earth models, and it was found that the results obtained from wPSOGSA are reliable, stable, and more accurate than the available results, which are well adapted to borehole lithology.
The datasets used for the present study and analysis have been taken from published papers, e.g. Shaw and Srivastava (2007), Xiong et al. (2018), Hutton et al. (1989), Jones and Hutton (1979).
M: conceptualization of the study, methodology, computer code, analysis, and drafting of the manuscript. KS: methodology, computer code, analysis, and drafting of the manuscript. UKS: supervision, suggestions, and editing.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors would like to express their gratitude to the IIT(ISM), Dhanbad, for providing a pleasant environment to pursue this study and support for the research. We also express our gratitude to the chief editor, associate editor, and anonymous reviewers, whose suggestions and comments enabled us to better understand the issues and considerably improve our manuscript.
This paper was edited by Norbert Marwan and reviewed by two anonymous referees.
Ai, H., Essa, K. S., Ekinci, Y. L., Balkaya, Ç., Li, H., and Géraud, Y.: Magnetic anomaly inversion through the novel barnacles mating optimization algorithm, Sci. Rep., 12, 22578, https://doi.org/10.1038/s41598-022-26265-0, 2022.
Cagniard, L.: Basic theory of the magneto-telluric method of geophysical prospecting, Geophysics, 18, 605—635, https://doi.org/10.1190/1.1437915, 1953.
Colorni, A., Dorigo, M., and Maniezzo, V.: Distributed Optimization by Ant Colonies, Proceedings of the First European Conference on Artificial Life, Paris, France, 134–142 pp., 1991.
Constable, S. C., Parker, R. L., and Constable, C. G.: Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data, Geophysics, 52, 289–300, https://doi.org/10.1190/1.1442303, 1987.
Dawes, G. J. K.: Magnetotelluric feasibility study: Island of Milos, Greece, Luxembourg, Edinburgh Univ. (UK). Dept. of Geophysics, Luxembourg, Report Number EUR-10674, Reference Number: ERA-13-007410, EDB-88-008365, 1986.
Dosso, S. E. and Oldenburg, D. W.: Magnetotelluric appraisal using simulated annealing, Geophys. J. Int., 106, 379–385, https://doi.org/10.1111/j.1365-246X.1991.tb03899.x, 1991.
Essa, K. S. and Diab, Z. E.: Gravity data inversion applying a metaheuristic Bat algorithm for various ore and mineral models, J. Geodyn., 155, 101953, https://doi.org/10.1016/j.jog.2022.101953, 2023.
Essa, K. S., Abo-Ezz, E. R., Géraud, Y., and Diraison, M.: A successful inversion of magnetic anomalies related to 2D dyke-models by a particle swarm scheme, J. Earth Syst. Sci., 132, 65, https://doi.org/10.1007/s12040-023-02075-4, 2023.
Hutton, V. R. S., Galanopoulos, D., Dawes, G. J. K., and Pickup, G. E.: A high resolution magnetotelluric survey of the Milos geothermal prospect, Geothermics, 18, 521–532, https://doi.org/10.1016/0375-6505(89)90054-0, 1989.
Jain, S. and Wilson, C. D. V.: Magneto-Telluric Investigations in the Irish Sea and Southern Scotland, Geophys. J. Int., 12, 165–180, https://doi.org/10.1111/j.1365-246X.1967.tb03113.x, 1967.
Jones, A. G. and Hutton, R.: A multi-station magnetotelluric study in southern Scotland – I. Fieldwork, data analysis and results, Geophys. J. Int., 56, 329–349, https://doi.org/10.1111/j.1365-246X.1979.tb00168.x, 1979.
Kennedy, J. and Eberhart, R.: Particle swarm optimization, in: Proceedings of ICNN'95 – International Conference on Neural Networks, 1942–1948 vol.4, https://doi.org/10.1109/ICNN.1995.488968, 1995.
Khishe, M. and Mosavi, M. R.: Chimp optimization algorithm, Expert Syst. Appl., 149, 113338, https://doi.org/10.1016/j.eswa.2020.113338, 2020.
Kirkpatrick, S., Gelatt C., D., and Vecchi M. P.: Optimization by Simulated Annealing, Science, 220, 671–680, https://doi.org/10.1126/science.220.4598.671, 1983.
Kunche, P., Sasi Bhushan Rao, G., Reddy, K. V. V. S., and Uma Maheswari, R.: A new approach to dual channel speech enhancement based on hybrid PSOGSA, Int. J. Speech Technol., 18, 45–56, https://doi.org/10.1007/s10772-014-9245-5, 2015.
Leslie, A. G., Millward, D., Pharaoh, T., Monaghan, A. A., Arsenikos, S., and Quinn, M.: Tectonic synthesis and contextual setting for the Central North Sea and adjacent onshore areas, 21CXRM Palaeozoic Project, https://nora.nerc.ac.uk/id/eprint/516757/1/21CXRM_Tectonic_synthesis_Leslieetal_CR_15_125N_Finalv2.pdf (last access: 24 March 2016), 2015.
Li, S.-Y., Wang, S.-M., Wang, P.-F., Su, X.-L., Zhang, X.-S., and Dong, Z.-H.: An improved grey wolf optimizer algorithm for the inversion of geoelectrical data, Acta Geophys., 66, 607–621, https://doi.org/10.1007/s11600-018-0148-8, 2018.
Lynch, S. M.: Introduction to applied Bayesian statistics and estimation for social scientists, Springer, New York, https://doi.org/10.1007/978-0-387-71265-9, 2007.
Miecznik, J., Wojdyła, M., and Danek, T.: Application of nonlinear methods to inversion of 1D magnetotelluric sounding data based on very fast simulated annealing, Acta Geophys. Pol., 51, 307–322, 2003.
Mirjalili, S. and Hashim, S. Z. M.: A new hybrid PSOGSA algorithm for function optimization, in: 2010 International Conference on Computer and Information Application, 374–377, https://doi.org/10.1109/ICCIA.2010.6141614, 2010.
Mirjalili, S., Mirjalili, S. M., and Lewis, A.: Grey Wolf Optimizer, Adv. Eng. Softw., 69, 46–61, https://doi.org/10.1016/j.advengsoft.2013.12.007, 2014.
Nabighian, M. N. and Asten, M. W.: Metalliferous mining geophysics – State of the art in the last decade of the 20th century and the beginning of the new millennium, Geophysics, 67, 964–978, https://doi.org/10.1190/1.1484538, 2002.
Pace, F., Raftogianni, A., and Godio, A.: A Comparative Analysis of Three Computational-Intelligence Metaheuristic Methods for the Optimization of TDEM Data, Pure Appl. Geophys., 179, 3727–3749, https://doi.org/10.1007/s00024-022-03166-x, 2022.
Pérez-Flores, M. A. and Schultz, A.: Application of 2-D inversion with genetic algorithms to magnetotelluric data from geothermal areas, Earth Planets Space, 54, 607–616, https://doi.org/10.1186/BF03353049, 2002.
Rashedi, E., Nezamabadi-pour, H., and Saryazdi, S.: GSA: A Gravitational Search Algorithm, Inf. Sci., 179, 2232–2248, https://doi.org/10.1016/j.ins.2009.03.004, 2009.
Rodi, W. and Mackie, R. L.: Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion, Geophysics, 66, 174–187, https://doi.org/10.1190/1.1444893, 2001.
Ross, S.: Probability and statistics for engineers and scientists, Elsevier, New Delhi, https://www.sciencedirect.com/book/9780123704832/introduction-to-probability-and-statistics-for-engineers-and-scientists (last access: 2014), 2009.
Roy, A. and Kumar, T. S.: Gravity inversion of 2D fault having variable density contrast using particle swarm optimization, Geophys. Prospect., 69, 1358–1374, https://doi.org/10.1111/1365-2478.13094, 2021.
Sen, M. K. and Stoffa, P. L.: Bayesian inference, Gibbs' sampler and uncertainty estimation in geophysical inversion1, Geophys. Prospect., 44, 313–350, https://doi.org/10.1111/j.1365-2478.1996.tb00152.x, 1996.
Sen, M. K. and Stoffa, P. L.: Global Optimization Methods in Geophysical Inversion, Cambridge University Press, Cambridge, https://doi.org/10.1017/CBO9780511997570, 2013.
S̨enel, F. A., Gökçe, F., Yüksel, A. S., and Yiğit, T.: A novel hybrid PSO–GWO algorithm for optimization problems, Eng. Comput., 35, 1359–1373, https://doi.org/10.1007/s00366-018-0668-5, 2019.
Sharma, S. P.: VFSARES – a very fast simulated annealing FORTRAN program for interpretation of 1-D DC resistivity sounding data from various electrode arrays, Comput. Geosci., 42, 177–188, https://doi.org/10.1016/j.cageo.2011.08.029, 2012.
Shaw, R. and Srivastava, S.: Particle Swarm Optimization: A new tool to invert geophysical data, Geophysics, 72, F75–F83, https://doi.org/10.1190/1.2432481, 2007.
Simon, D.: Biogeography-Based Optimization, IEEE Trans. Evol. Comput., 12, 702–713, https://doi.org/10.1109/TEVC.2008.919004, 2008.
Simpson, F. and Bahr, K.: Practical Magnetotellurics, Cambridge University Press, https://doi.org/10.1017/CBO9780511614095, 2005.
Stewart, A. L. and McPhie, J.: Facies architecture and Late Pliocene – Pleistocene evolution of a felsic volcanic island, Milos, Greece, Bull. Volcanol., 68, 703–726, https://doi.org/10.1007/s00445-005-0045-2, 2006.
Storn, R. and Price, K.: Differential Evolution – A Simple and Efficient Heuristic for global Optimization over Continuous Spaces, J. Glob. Optim., 11, 341–359, https://doi.org/10.1023/A:1008202821328, 1997.
Tarantola, A.: Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics, https://doi.org/10.1137/1.9780898717921, 2005.
Tarantola, A. and Valette, B.: Generalized nonlinear inverse problems solved using the least squares criterion, Rev. Geophys., 20, 219–232, https://doi.org/10.1029/RG020i002p00219, 1982.
Ward, S. H. and Hohmann, G. W.: 4. Electromagnetic Theory for Geophysical Applications, in: Electromagnetic Methods in Applied Geophysics: Volume 1, Theory, Society of Exploration Geophysicists, 130–311, https://doi.org/10.1190/1.9781560802631.ch4, 1988.
Wen, L., Cheng, J., Li, F., Zhao, J., Shi, Z., and Zhang, H.: Global optimization of controlled source audio-frequency magnetotelluric data with an improved artificial bee colony algorithm, J. Appl. Geophys., 170, 103845, https://doi.org/10.1016/j.jappgeo.2019.103845, 2019.
Whitley, D.: A genetic algorithm tutorial, Stat. Comput., 4, 65–85, https://doi.org/10.1007/BF00175354, 1994.
Xiong, J., Liu, C., Chen, Y., and Zhang, S.: A non-linear geophysical inversion algorithm for the mt data based on improved differential evolution, Eng. Lett., 26, 161–170, 2018.
Yang, X.-S.: A New Metaheuristic Bat-Inspired Algorithm, in: Nature Inspired Cooperative Strategies for Optimization (NICSO 2010), edited by: González, J. R., Pelta, D. A., Cruz, C., Terrazas, G., and Krasnogor, N., Springer Berlin Heidelberg, Berlin, Heidelberg, 65–74, https://doi.org/10.1007/978-3-642-12538-6_6, 2010a.
Yang, X.-S.: Firefly algorithm, stochastic test functions and design optimisation, Int J Bio Inspired Comput, 2, 78–84, https://doi.org/10.48550/arxiv.1003.1409, 2010b.
Zhang, Z., Ding, S., and Jia, W.: A hybrid optimization algorithm based on cuckoo search and differential evolution for solving constrained engineering problems, Eng. Appl. Artif. Intell., 85, 254–268, https://doi.org/10.1016/j.engappai.2019.06.017, 2019.