Stability and uncertainty assessment of geoelectrical resistivity model parameters: a new hybrid metaheuristic algorithm and posterior probability density function approach

Estimating a reliable subsurface resistivity structure using conventional techniques is challenging due to the nonlinear nature of the inverse problems. The performance of the inversion techniques can be pretty ambiguous based on the optimal error, although traditional methods have proven to be quite effective. In this work, the impacts of the constraints accessible from a borehole are examined for further assessment and to enhance algorithm effectivity. The vPSOGWO strategy is a new approach that is based on a model search space without any prior information, and it describes the hybridization of particle swarm optimization (PSO) with the Grey Wolf Optimizer (GWO). To understand the efficiency and novelty of the algorithm, it has been validated on two different kinds of synthetic resistivity data with various sets of noise and, subsequently, applied to three field datasets of different geological terrains. The analyzed results suggest that the subsurface resistivity model shows considerable uncertainty. Thus, it is superior to examine the histograms and posterior probability density functions (PDFs) of such solutions to exemplify the global solution. A PDF with a 68.27 % confidence interval (CI) selects a region with a higher probability. Therefore, the inverted models are used to estimate the mean global solution and the most negligible uncertainties, where the mean global solution represents the best solution. Our vPSOGWO-inverted outcomes have been proven to be more accurate than classic PSO, GWO, and state-of-the-art variants of classic approaches. As a result, this novel method plays a vital role in vertical electrical sounding (VES) data inversion.

1 Introduction

The vertical electrical resistivity sounding (VES) technique is an economical and simple method that has been used to determine the layered parameters in a wide range of applications in the hydrogeological, groundwater, mineral, geothermal, hydrocarbon, engineering, and environmental fields, among others (Sen et al., 1993; Sharma, 2012; Panda et al., 2018). VES data interpretation is challenging due to its unstable, nonunique solution and algorithm sensitivity (Narayan et al., 1994; Oldenburg and Li, 1994; Singh et al., 2005, 2013). Therefore, many researchers have developed several inversion algorithms to improve accuracy and stability and to reduce uncertainty in the solutions. These inversion techniques are grouped into local and global optimization techniques. In the local inversion techniques, a logical initial guess is required to get the solution. This has led researchers to think about alternative methods via which a broad range of parameters can be established. Researchers have developed various metaheuristic optimization algorithms to solve various real-world problems. These algorithm types, inspired by natural phenomenon, include ant colony optimization (ACO; Colorni et al., 1991), the Bat Algorithm (Yang, 2010), biogeographically based optimization (Simon, 2008), differential evolution (DE; Storn and Price, 1997), the Firefly Algorithm (Yang, 2010), the Genetic Algorithm (GA; Whitley, 1994; Mitchell, 1998), the Gravitational Search Algorithm (GSA; Rashedi et al., 2009), the Grey Wolf Optimizer (GWO; Mirjalili et al., 2014), and particle swarm optimization (PSO; Kennedy and Eberhart, 1995). These optimization techniques aim to have an optimum solution and fast convergent rate to obtain global minima. However, unique characteristics, namely, exploration and exploitation, in global optimization algorithms persist. For example, PSO has a very high potential regarding exploitation, implying that the algorithm performs well with respect to a local search (Şenel et al., 2019), but it is inferior regarding exploration, which means that the algorithm has less ability with respect to establishing the starting position near global minima and, due to low exploration characteristics, it gets trapped at the local minima (Eiben and Schippers, 1998; Mirjalili and Hashim, 2010). Therefore, integrating two algorithms with the opposite characteristics is the best way to solve the exploration characteristics and exploitation characteristics and to provide a more accurate and reliable solution, compared with results obtained with an individual algorithm. Many authors have developed various hybrid metaheuristic algorithms such as PSOGA for fundamental function analysis, PSOACO for data mining, PSODE for global optimization using the standard function, and PSOGSA using the standard function (Esmin et al., 2013; Lai and Zhang, 2009; Rashedi et al., 2009).

This study focuses on a variable-weight hybrid algorithm, known as vPSOGWO (Şenel et al., 2019), that fuses the exploration ability of PSO with the exploration ability of GWO. In this algorithm, some random particles of PSO are replaced with new ones obtained from GWO. In prior work, the constant-weight hybrid technique of PSO and GWO, known as HPSOGWO, has been used by some authors for different applications, such as for single-area-unit commitment problems (Kamboj, 2015), mathematical problems (Singh and Singh, 2017), and benchmark functions and real-world issues (Şenel et al., 2019). However, to the best of our knowledge, none of these researchers have tested these methods on geophysical data inversion. Thus, in this study, the applicability of the vPSOGWO algorithm is demonstrated on synthetic data with noise, synthetic data without noise, and various field resistivity sounding data to estimate the resistivity distribution in a 1D Earth's subsurface model. This work also calculates the posterior probability density functions (PDFs) with a 68.27 % confidence interval (CI) and correlation matrix on all accepted models to determine the mean global model and uncertainty. As a result, we analyzed and compared the effectiveness of the proposed algorithms with classic PSO, GWO, and state-of-the-art variants of classic methods. Our analysis advocates for the fact that the vPSOGWO algorithm produces a more accurate and reliable model with excellent stability, the least model uncertainty, and the ability to successfully resist noise.

2 Forward modeling algorithm

With the help of the kernel function (Koefoed, 1979) and Schlumberger resistivity configuration (Fig. 1), the forward code is developed and synthetic resistivity datasets are produced from known parameters, such as the current electrode spacing and the number of geological multilayers of true resistivities and their corresponding thicknesses. The mathematical expression of the apparent resistivity is expressed as follows:

$$\begin{array}{}\text{(1)}& {\mathit{\rho }}_a(\mathit{s},\mathit{m})={\mathit{\rho }}_{\mathrm{1}}+{\mathit{s}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{1}}{\int }_{\mathrm{0}}^{\mathrm{\infty }}{T}_{\mathrm{1}}(\mathit{\lambda },\mathit{m}){\mathit{J}}_{\mathrm{1}}\left(\mathit{\lambda }\mathit{s}\right)\mathrm{d}\mathit{\lambda },\end{array}$$

where J₁ is the first-order Bessel function, λ represents the integration variables, s is the current electrode spacing, and m is the model.

Figure 1 Schlumberger array configuration for the three-layer case: C1 and C2, through which the current is injected, are current electrodes with spacing s, while P1 and P2 are potential electrodes with spacing b.

Figure 2 Three-layer synthetic data: (a) observed data (*) and best-fitted calculated apparent resistivity curve (>68.27 % PDF); (b) the 1D mean model (>68.27 % PDF) for the true model (black), vPSOGWO (red), GWO (blue), and PSO (green).

For each layer, the kernel's resistivity transform T_k has been determined by Pekeris (1940). The apparent resistivity, T_k(λ), is convoluted with linear filter theory to compute the following:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathit{T}}_k\left(\mathit{\lambda }\right)& ={\mathit{\rho }}_k\times \left({\mathit{T}}_{k+\mathrm{1}}\left(\mathit{\lambda }\right)+{\mathit{\rho }}_k\tanh \left(\mathit{\lambda }{\mathit{t}}_k\right)\right)/\\ & \left({\mathit{\rho }}_k+{\mathit{T}}_{k+\mathrm{1}}\left(\mathit{\lambda }\right)\tanh \left(\mathit{\lambda }{\mathit{t}}_k\right)\right).\end{array}\end{array}$$

3 Inverse modeling algorithm

The geophysical inverse problem can be formulated through a forward modeling operator/functional with the aim of achieving the geophysical model/solution that best illuminates the observed data. This operator integrates the geophysical problems and maps between the observed data y and the solution x as follows:

$$\begin{array}{}\text{(3)}& \mathit{y}=f\left(\mathit{x}\right).\end{array}$$

The inversion techniques minimize the cost functional/misfit functional, which is generally a degree of the relationship between the N number of observed data (y_o) and the calculated data (y_c). This misfit functional can be introduced here as a mean square error (MSE) and can be defined as follows:

$$\begin{array}{}\text{(4)}& \mathrm{MSE}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N{\left(y_{\mathrm{o}}-y_{\mathrm{c}}\right)}^{\mathrm{2}}.\end{array}$$

3.1 Particle swarm optimization

Particle swarm optimization (PSO) is based on the social behavior of animals, such as the schooling behavior of fish or the flocking behavior of birds (Kennedy and Eberhart, 1995). When birds go in search of food, they scatter randomly within a search space before they can determine the position of food. While searching for food, there is always a bird who is aware of the position of food, and they share this information with others. Using this method, each bird is referred to as a particle and is represented by geophysical solutions/models (i.e., here, a particle is a resistivity layer parameter). The capability/fitness of each swarm/bird is estimated between the N number of observed data (y_o), which measure the swarm and the food distance, and the computed data (y_c), which measure the swarm and the estimated position (resistivity layer parameter/solution) of the prey distance using Eq. (4).

Figure 3 Convergence curve for the best-fitted model parameters for the vPSOGWO algorithm.

Figure 4 Convergent curve, also known as the error versus iteration curve, for three-layer noiseless synthetic data.

The best position among particles with information about it is stored in memory for each iteration. The new velocity and position of the population pool are accepted if their possibility is large, otherwise they are rejected. In that case, the particles are randomly distributed in the search space in order to escape the local optima. The search continues until it gains a maximum possibility or reaches the maximum iteration. In the global search space, the position of each particle is updated by the following two mathematical equations:

$$\begin{array}{}\text{(5)}& {\displaystyle }\begin{array}{rl}{\mathit{v}}_i\left(t+\mathrm{1}\right)& ={\mathit{v}}_i\left(t\right)+c_{\mathrm{1}}\times \mathrm{rand}\left({\mathit{x}}_{\mathrm{p}}\left(t\right)-{\mathit{x}}_i\left(t\right)\right)\\ & +c_{\mathrm{2}}\times \mathrm{rand}\times \left({\mathit{x}}_g-{\mathit{x}}_i\left(t\right)\right),\end{array}\text{(6)}& {\displaystyle }{\mathit{x}}_i(t+\mathrm{1})={\mathit{x}}_i\left(t\right)+{\mathit{v}}_i(t+\mathrm{1}).\end{array}$$

Here, v_i represents the velocity of the ith particle with position x_i, x_p is the best position obtained by the ith particle, x_g is the best position, t is the number of the iteration, i represents the number of the model ($i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots },N)$, rand represent the random values with a range of [0, 1], and the coefficients c₁ and c₂ represent the optimization parameter. The disadvantage of the PSO algorithm is that, while directing particles to random positions, it has a small possibility of escaping the local minima.

3.2 The Grey Wolf Optimizer (GWO) algorithm

The GWO algorithm mimics the leadership hierarchy and hunting mechanics of grey wolves and is used to solve both standard and real-life problems. In the grey wolf community, animals are divided in four groups: (i) alpha animals, (ii) beta animals, (iii) delta animals, and (iv) omega animals. Alpha, beta, and delta animals are the fittest wolves, and they guide omega animals towards promising areas of the search space. The alpha is the pack leader and generally makes important and final decision for all of the wolves; thus, the alpha represents the fittest solution. The betas are subordinates and help the alphas in their decision-making. However, betas cannot force alphas into any decision; they can only order the lower wolves. The beta group takes orders from the alpha group, enforces orders with respect to the other groups, and sends feedback back to the alpha group. All of the groups are dominant with respect to the omega group. Nevertheless, the omega group is an important component of the pack during hunting, as they play the role of the scapegoat and are only allowed to eat at the end. If a wolf is not part of the alpha, beta, or omega group, they are known as delta and only summit to alpha and beta groups. In the GWO algorithm, the alpha group represents the best position, i.e., geophysical model/solution. In our case, the geophysical model is the resistivity layer parameters. The beta and delta groups are consecutive best solutions, and the omega group is the best solution that always follows the other groups. The capability/fitness of each wolf is estimated between the observed data (which measure the wolf and prey distance) and the computed data (which measure the wolf and the estimated position of the prey distance) using Eq. (4).

Figure 5 (a) Histogram and (b) posterior PDF of all 10 000 solutions corresponding to the output of each run for the three-layer synthetic Earth model.

Figure 6 Correlation plot between model parameters (off-diagonal) and the posterior PDF curve (diagonal) for those models whose PDF exceeds 68.27 % of the confidence interval (CI).

Figure 7 Four-layer synthetic data: (a) observed data (*) and the best-fitted calculated apparent resistivity curve (>68.27 % PDF); (b) the 1D mean model (>68.27 % PDF) for the true model (black), vPSOGWO (red), GWO (blue), and PSO (green).

Figure 8 Convergent curve, also known as the error versus iteration curve, for four-layer noiseless synthetic resistivity sounding data.

Figure 9 Histogram of the logarithmic mean square error for vPSOGWO, GWO, and PSO over 10 000 models. The x axis of the three histograms represents the misfit error corresponding to 10 000 models.

Figure 10 (a) Histogram and (b) posterior PDF of all 10 000 solutions corresponding to the output of each run for four-layer synthetic resistivity sounding data.

Figure 11 Correlation plot between the model parameters (off-diagonal) and posterior PDF curve (diagonal) for models whose PDF exceeds 68.27 % of the confidence interval (CI).

Hunting in the grey wolf community has been divided into three groups: searching for prey, encircling the prey, and attacking the prey. The encircling nature of the wolves is defined by the following equations:

$$\begin{array}{}\text{(7)}& {\displaystyle }d=\left|c\times \left(t\right)-{\mathit{x}}_i\left(t\right)\right|,\text{(8)}& {\displaystyle }{\mathit{x}}_i\left(t+\mathrm{1}\right)={\mathit{x}}_{\mathrm{p}}\left(t\right)-a\times d.\end{array}$$

Here, x_p is the prey position, x_i is the grey wolves' positions, and a and c are the vectors mathematically formulated as follows:

$$\begin{array}{}\text{(9)}& {\displaystyle }a=a_{\mathrm{1}}\times \left(\mathrm{2}\times \mathrm{rand}-\mathrm{1}\right),\text{(10)}& {\displaystyle }c=\mathrm{2}\times \mathrm{rand}.\end{array}$$

Here, $a_{\mathrm{1}}=\mathrm{2}\times \left(\mathrm{1}-t/l\right)$ and varies from 2 to 0 in decreasing order with increasing iteration (t), l represents the maximum iteration, and rand is the random number in the range of [0, 1].

In the grey wolf community, the alpha group leads, the beta and the delta groups search for the prey location, and the omega group follows the other groups. Therefore, the alpha group gives the best solution, while the respective second- and third-best solutions are provided by the beta and the delta groups. Thus, the remaining wolves, i.e., the omega group, follow the best solution (the other wolf groups) to obtain the best location. This is mathematically equated as follows:

$$\begin{array}{}\text{(11)}& d_{\mathit{\alpha },\mathit{\beta },\mathit{\delta }}=\mathrm{|}{\mathit{c}}_{\mathrm{1},\mathrm{2},\mathrm{3}}\times {\mathit{x}}_{\mathit{\alpha },\mathit{\beta },\mathit{\delta }}-\mathit{x}\mathrm{|}.\end{array}$$

The best location/position for alpha, beta, and delta wolves in each iteration is given by x_α, x_β, and x_δ, respectively:

$$\begin{array}{}\text{(12)}& {\mathit{x}}_{\mathrm{1},\mathrm{2},\mathrm{3}}=\mathrm{|}{\mathit{x}}_{\mathit{\alpha },\mathit{\beta },\mathit{\delta }}-{\mathit{a}}_{\mathrm{1},\mathrm{2},\mathrm{3}}\times {\mathit{d}}_{\mathit{\alpha },\mathit{\beta },\mathit{\delta }}\mathrm{|}.\end{array}$$

Here, x_p(t+1) describes the updated position of the prey in the (t+1) iteration and is obtained from the mean position of three best wolves in the population; thus,

$$\begin{array}{}\text{(13)}& {\mathit{x}}_{\mathrm{p}}\left(t+\mathrm{1}\right)=\left({\mathit{x}}_{\mathrm{1}}+{\mathit{x}}_{\mathrm{2}}+{\mathit{x}}_{\mathrm{3}}\right)/\mathrm{3}.\end{array}$$

The values of a are utilized by wolves who force the search to move away from the prey. When a≥1, hunting is abandoned in order to find a better solution; in contrast, when a<1, the wolves are forced to attack the prey. In Eq. (9), a varies in the range of $\left[-\mathrm{2}{a}_{\mathrm{1}},\mathrm{2}{a}_{\mathrm{1}}\right]$.

Figure 12 Three-layer field data over Mount Turner, northern Queensland, Australia: (a) observed data (*) and the best-fitted calculated apparent resistivity curve (>68.27 % PDF); (b) the 1D mean model (>68.27 % PDF) for the true model (black), vPSOGWO (red), GWO (blue), and PSO (green).

Figure 13 Five-layer field data: (a) observed data (*) and the best-fitted calculated apparent resistivity curve (>68.27 % PDF); (b) the 1D mean model (>68.27 % PDF) for the true model (black), vPSOGWO (red), GWO (blue), and PSO (green).

Figure 14 Six-layer field data over Keshiari (Kharagpur), India: (a) observed data (*) and the best-fitted calculated apparent resistivity curve (>68.27 % PDF); (b) the 1D mean model (>68.27 % PDF) for the true model (black), vPSOGWO (red), GWO (blue), and PSO (green).

3.3 Variable-weight hybrid PSOGWO (vPSOGWO)

Despite usefulness of the PSO technique with respect to achieving successful results in real-world problems, it tends to fall into the local minima, causing the solution to move away from global minima. This tendency for deterioration within the local minima is stopped by the explorative ability of the GWO algorithm. Therefore, the variable-weight hybrid PSOGWO, known as vPSOGWO, fuses the exploitation potential of PSO with the exploration potential of GWO to overcome each other's discrepancy via the implementation of varying weight. Due to the involvement of two distinct variants running together to solve the problem, this hybrid vPSOGWO is called a coevolutionary hybrid algorithm. The encircling behavior of each wolf is updated by the following:

$$\begin{array}{}\text{(14)}& {\displaystyle }{\mathit{d}}_{\mathit{\alpha },\mathit{\beta },\mathit{\delta }}=\mathrm{|}{\mathit{c}}_{\mathrm{1},\mathrm{2},\mathrm{3}}\times {\mathit{x}}_{\mathit{\alpha },\mathit{\beta },\mathit{\delta }}-w\times \mathit{x}\mathrm{|},\text{(15)}& {\displaystyle }\text{where},w=w_{\max }-\left(w_{\max }-w_{\min }\right)\times t/l.\end{array}$$

Here, w_max=0.9 and w_min=0.2 are found to be more appropriate after tuning for our study.

The best location/position (geophysical model) for alpha, beta, and delta wolves in each iteration is given by x_α, x_β, and x_δ, respectively.

$$\begin{array}{}\text{(16)}& {\mathit{x}}_{\mathrm{1},\mathrm{2},\mathrm{3}}=\mathrm{|}{\mathit{x}}_{\mathit{\alpha },\mathit{\beta },\mathit{\delta }}-{\mathit{a}}_{\mathrm{1},\mathrm{2},\mathrm{3}}\times {\mathit{d}}_{\mathit{\alpha },\mathit{\beta },\mathit{\delta }}\mathrm{|},\end{array}$$

where

$$\begin{array}{}\text{(17)}& {\displaystyle }{a}_{\mathrm{1},\mathrm{2},\mathrm{3}}=a_{\mathrm{1}}\times \left(\mathrm{2}\times \mathrm{rand}-\mathrm{1}\right),\text{(18)}& {\displaystyle }{c}_{\mathrm{1},\mathrm{2},\mathrm{3}}=\mathrm{0.5}\left(\text{chosen after tuning}\right),\text{(19)}& {\displaystyle }{a}_{\mathrm{1}}=\mathrm{2}\times \left(\mathrm{1}-t/l\right).\end{array}$$

The updated velocity and position of vPSOGWO are as follows:

$$\begin{array}{}\text{(20)}& {\displaystyle }\begin{array}{rl}{\mathit{v}}_i\left(t+\mathrm{1}\right)& =w\times {\mathit{v}}_i\left(t\right)+c_{\mathrm{1}}\times \mathrm{rand}\times \left({\mathit{x}}_{\mathrm{1}}-{\mathit{x}}_i\left(t\right)\right)\\ & +c_{\mathrm{2}}\times \mathrm{rand}\times \left({\mathit{x}}_{\mathrm{2}}-{\mathit{x}}_i\left(t\right)\right)+c_{\mathrm{3}}\times \mathrm{rand}\\ & \times \left({\mathit{x}}_{\mathrm{3}}-{\mathit{x}}_i\left(t\right)\right),\end{array}\text{(21)}& {\displaystyle }{\mathit{x}}_i(t+\mathrm{1})={\mathit{x}}_i\left(t\right)+{\mathit{v}}_i(t+\mathrm{1}).\end{array}$$

Here, the value of 1.5 is found to be more suitable for each of the coefficients (c₁, c₂, and c₃) after tuning the parameters in the present study (Roshan and Singh, 2017).

Algorithm 1 The vPSOGWO algorithm.

Max_Iter: maximum iterations set

Pop_no: population size

Para: number of parameters

Fitness = infinite: already set

Lb and Ub: set lower bound (Lb) and upper bound (Ub) for different parameters

Initialize particles randomly

Procedure

for l=1 to Max_Iter do

for i=1 to Pop_no do

for j=1 to Para do

check the Lb and Ub for randomly created particles

end for

end for

for i=1 to Pop_no do

Calculate the fitness form cost function

Update the wolves' fitness and position

end for

Update a1, a, c, and w using Eqs. (15) and (17)–(19)

for i=1 to Pop_no do

for j=1 to Para do

Update position of x₁, x₂, and x₃ using Eqs. (14) and (16)

Update best particle velocity and position using Eqs. (20)–(21)

end for

end for

end for

Table 1 Optimization of the mean model result for three-layer synthetic resistivity sounding data.

Table 2 Correlation matrix using a 68.27 % PDF limit for three-layer synthetic resistivity sounding data.

4 Statistical simulation for global model and uncertainty estimation

The proposed algorithms yield good-fitting models, but the evaluation of a global solution requires numerous techniques. This is noteworthy with respect to selecting the region of solution/model search space, where we find enormous solutions. The methods for selecting the region of model space were chosen to envisage the global solution and reduce the uncertainty in the ultimate solution (Mosegaard and Tarantola, 1995; Sen and Stoffa, 1996). Thus, many solutions and the associated error estimated were kept in memory for consequent statistical measurements. Therefore, 10⁸ solutions were generated for each algorithm using the logarithmic mean square error, and every computed response corresponding to each model fits well with the observed response. However, the model parameters obtained, which lie within the search range in multidimensional space, may differ from each other. Hence, the mean model from the model parameters is defined as follows (Ross, 2009):

$$\begin{array}{}\text{(22)}& {\widehat{\mathit{m}}}_i={\displaystyle \frac{\mathrm{1}}{M}}{\sum }_{j=\mathrm{1}}^M{\mathit{m}}_{i,j},\end{array}$$

where i represents the layer parameters, j is the number of models, and m_i,j is the jth model with i number of parameters.

All algorithms are executed for 10 000 runs with 1000 iterations to obtain the best model parameters. It is important to mention that multiple runs are crucial in vPSOGWO, as 1000 weightage points lay between the inertial weights of 0.9 and 0.2, such that each weightage point yields a fitted model in a run. As a result, 10 000 runs provide 10 000 chances for each weightage point to fetch the best-fitted model.

Table 3 Stability test for three-layer synthetic resistivity sounding data using different search ranges.

Table 4 Optimization mean model result for three-layer synthetic resistivity sounding data with 10 % noise.

Therefore, the posterior covariance matrices are defined using the following equation (Ross, 2009):

$$\begin{array}{}\text{(23)}& \mathrm{Cov}\left({\mathit{m}}_{i,k}\right)={\displaystyle \frac{\mathrm{1}}{M-\mathrm{1}}}{\sum }_{j=\mathrm{1}}^M\left({\mathit{m}}_{i,j}-{\widehat{\mathit{m}}}_i\right)\times \left(m_{k,j}-{\widehat{\mathit{m}}}_k\right).\end{array}$$

Posterior correlation matrices are described using the following equation:

$$\begin{array}{}\text{(24)}& \mathrm{Corr}\left({\mathit{m}}_{i,k}\right)=\mathrm{Cov}\left({\mathit{m}}_{i,k}\right)/\sqrt{\mathrm{Cov}\left({\mathit{m}}_{i,i}\right)\times \mathrm{Cov}\left({\mathit{m}}_{k,k}\right)}.\end{array}$$

Here, i and k lie between 1 and total number of parameters.

The square roots of the diagonal elements of the covariance matrix define the uncertainty in the solution, and the correlation matrix gives a rough idea about the relation between the model parameters. If the parameters do not provide a global solution, the apparent resistivity curve corresponding to the mean model will not adequate to the observed value. The posterior correlation matrix corresponding to the indigenous solution will not yield an actual correlation between the parameters obtained via linear regression. For further analysis, posterior PDFs and histograms are calculated over all accepted models. The 1D posterior PDF for various parameters with mean ${\widehat{\mathit{m}}}_i$ and standard deviation σ_i is given as follows (Ross, 2009):

$$\begin{array}{}\text{(25)}& p\left({\mathit{y}}_i,{\widehat{\mathit{m}}}_i,{\mathit{\sigma }}_i\right)=\left(\mathrm{1}/{\mathit{\sigma }}_i\sqrt{\mathrm{2}\mathit{\pi }}\right)\times \exp (-{\left({\mathit{y}}_i-{\widehat{\mathit{m}}}_i\right)}^{\mathrm{2}}/\mathrm{2}{\mathit{\sigma }}_i^{\mathrm{2}}),\end{array}$$

where y is the solution/model parameter's output store after 10 000 runs of an algorithm and i is between 1 and the number of model parameters.

Different techniques are based on the posterior PDF to obtain the global solution. One of these techniques is to pick the model parameters with the highest probability values. Another method based on the PDF is to normalize (0 to 1) each model parameter by its respective highest probability value. The best model is considered to have the highest sum of normalized probability values (Sharma, 2012). Furthermore, the best model can also be determined by taking the mean of each parameter with probably more significance than the threshold probability. However, these techniques fail to provide the global model.

Therefore, proceeding with a new approach to the study, we introduce a confidence interval (CI) more significant than 68.27 % as a benchmark for all model parameters. According to the empirical rule, 68.27 % of the data lie within 1 standard deviation of the mean (Ross, 2009). Thus, the model parameters below a CI of 68.27 % are discarded, and the remaining parameters are used to determine the mean solution and uncertainty. This means that the model represents the global solution with less uncertainty.

Table 5 Correlation matrix using a 68.27 % PDF limit for three-layer synthetic resistivity sounding data with 10 % noise.

Table 6 Optimization of the mean model result for four-layer synthetic resistivity sounding data.

5 Computational information

The code used in this work was developed in MATLAB R2019a on a Windows 10 platform with the following configuration: an HP Z240 Tower Workstation and an Intel Xeon E3-1225 v6 3.30 GHz CPU with 32.0 GB of RAM and a 64-bit operating system (OS). However, global optimization is a time-consuming process, as it requires many forwarding calculations to obtain the best-fitted result.

6 Results and discussion

The applicability of the new algorithm (vPSOGWO), GWO, and PSO has been assessed by inverting several cases of synthetic and field data extracted from different geological terrains (Dixon and Doherty, 1977; Panda et al., 2018). Both synthetic and field datasets were computed and optimized using the developed algorithms, maintaining a population size of 10 and 1000 iterations for 10 000 runs, leading each algorithm to analyze 10⁸ models. We discuss the inverted results of the algorithms with respect to their application in a few example synthetic and field cases in the following sections.

6.1 Example 1: synthetic data – three-layer case

Initially, to access the applicability and efficacy of the proposed algorithms, a synthetic apparent resistivity sounding dataset measured with a Schlumberger array was generated considering a three-layer Earth model sandwich with a highly resistant layer of 500.0 Ωm and a thickness of 150.0 m between two low-resistance layers of 8.0 and 5.0 Ωm, respectively. The synthetic data were computed in the MATLAB environment, as shown in Fig. 2a using asterisks. Figure 2a shows the three-layer synthetic data with the best-fitted calculated apparent resistivity curve (>68.27 % PDF), while Fig. 2b presents the 1D mean model (>68.27 % PDF) for the true model (black), vPSOGWO (red), GWO (blue), and PSO (green).

Table 7 Correlation matrix using a 68.27 % PDF limit for four-layer synthetic resistivity sounding data.

Table 8 Optimization of the mean model result for four-layer synthetic resistivity sounding data with 10 % noise.

The search limit for novel inversions techniques (vPSOGWO, GWO, and PSO) is carefully chosen, as shown in Table 1. Each algorithm, including vPSOGWO, runs 10 000 times to perform statistical analysis and determine the global mean model with the least uncertainty. Figure 3 shows the convergence curve of the resistivity layer parameters using vPSOGWO. We found no changes in the convergence pattern after 590 iterations, and layer parameters became stable. The convergence curves, in terms of error versus iterations, for three existing algorithms are shown in Fig. 4. It is observed that vPSOGWO, GWO, and PSO converged at 590, 950, and 380 iterations and have mean square errors of $\mathrm{1.586}\times {\mathrm{10}}^{-\mathrm{8}}$, $\mathrm{5.238}\times {\mathrm{10}}^{-\mathrm{8}}$, and $\mathrm{5.792}\times {\mathrm{10}}^{-\mathrm{8}}$, respectively, whereas ridge regression has an error of 0.633.

The 10 000 models inverted are used to find the posterior PDF and histogram for each parameter. As shown in Fig. 5b, the peak of the posterior PDF is roughly close to the actual model parameter. The histogram shown in Fig. 5a suggests that ρ₂ and h₂ have a broader range. This represents the equivalence problem associated with the resistive layer, as the uncertainty in each algorithm was found to be large considering all of the accepted models. Therefore, selecting the models with a posterior PDF with a CI greater than 68.27 % reduces the uncertainty in the model, increases the resolution of a solution, and helps estimate the best mean model close to the actual model (Table 1). Table 1 shows the model parameters and uncertainty for the proposed algorithms.

Here, two approaches are used to present the mean solution with its uncertainty estimation: (1) the mean solution for all accepted best-fitted solutions obtained from 10 000 runs for all three algorithms is used or (2) the mean model calculated from the solution with a posterior PDF (for which values are greater than 68.27 % CI from all accepted solution parameters) is used.

Here, we observed that the second-layer parameters for PSO and GWO are too diverted from actual values (higher uncertainty) due to their inability to balance exploitation and exploration properties. In contrast, the hybrid vPSOGWO algorithm provides more accurate results and falls within the uncertainty ranges (Table 1). Therefore, the hybrid algorithm has a more balanced nature with respect to exploitation and exploration than PSO and GWO. As shown in Table 2, the posterior correlation matrix illustrates that the first-layer resistivity reveals a feeble correlation with other associated parameters. However, there is a negative correlation found between ρ₂ and h₂ (both parameters have a trade-off relationship). In contrast, a positive correlation is observed between ρ₂ and h₁ (i.e., resistivity of the second layer increases with increasing the thickness of the first layer and vice versa). Similarly, a correlation can also be seen between third-layer resistivity and second-layer thickness, but it is inverse in nature.

Figure 6 represents the correlation plot between model parameters (off-diagonal) and the posterior PDF curve (diagonal) for models with a CI greater than 68.27 % for all parameters. No significant error differences are found between the observed and calculated apparent resistivity data for all three algorithms (Fig. 2a). However, the error difference in the 1D model and the result for the mean model with a CI of 68.27 % are presented in Fig. 2b and Table 1, respectively.

To check the stability of the parameter, the hybrid algorithm is tested with three different search spaces, as shown in Table 3. Consequently, it estimates the mean model and uncertainty for 100 runs. Table 3 suggests that using a broader search space does not cause the result to divert too much from the actual model. In this example, the computational time required for one run with 30 data points was 1.54, 1.49, and 1.48 s for vPSOGWO, GWO, and PSO, respectively.

The proposed optimization is also performed using the same synthetic data with 10 % Gaussian noise and keeping the search range shown in Table 1. The same procedure is applied to determine the mean model from all best-fitted solutions and solutions with a posterior PDF with a CI greater than 68.27 % for parameters of all of the solutions (Table 4). Although 10 % noise is added, the result obtained from the mean model for a posterior PDF of 68.27 % for the hybrid algorithm is not much different compared to actual values. At the same time, the error was observed to slightly increase: $\mathrm{1.309}\times {\mathrm{10}}^{-\mathrm{5}}$, $\mathrm{1.313}\times {\mathrm{10}}^{-\mathrm{5}}$, and $\mathrm{1.327}\times {\mathrm{10}}^{-\mathrm{5}}$ for vPSOGWO, GWO, and PSO, respectively. Table 5 depicts the correlation matrix of vPSOGWO, which clearly shows interdependence with values of 0.3315 and −0.7879 for the first- and second-layer parameters. Similarly, we can also determine the relation between second-layer resistivity and first-layer thickness (0.6142), third-layer resistivity, and the second-layer thickness (−0.7618). Hence, the result obtained using the proposed technique shows good agreement with the actual model values.

6.2 Example 2: synthetic data – four-layer case

A four-layer Earth model with a thin, relatively low-resistance layer (24.0 Ωm) sandwiched between two high-resistance layers (840.0 Ωm and 8400.0 Ωm) is considered for the demonstration of the proposed algorithms. Table 6 illustrates the actual model for synthetic data, the search range, and the inverted results. vPSOGWO, GWO, and PSO converge at iterations 590, 674, and 750 and with associated errors of $\mathrm{3.624}\times {\mathrm{10}}^{-\mathrm{8}}$, $\mathrm{1.370}\times {\mathrm{10}}^{-\mathrm{8}}$, and $\mathrm{2.097}\times {\mathrm{10}}^{-\mathrm{7}}$, respectively, as shown in Fig. 8, whereas the error estimated using the ridge regression method is 0.383. Instead of higher error in vPSOGWO compared with GWO, it can be observed that the error scale for the vPSOGWO algorithm is narrower than the other two algorithms, which is an essential factor for determining the mean model (Fig. 9). Hence, the mean model is affected by the error scale, as shown in Fig. 9.

Table 9 Correlation matrix using a 68.27 % PDF limit for four-layer synthetic resistivity sounding data with 10 % noise.

Table 10 Optimization of the mean model result for three-layer field resistivity sounding data.

Table 11 Correlation matrix using a 68.27 % PDF limit for three-layer field resistivity sounding data.

Table 12 Optimization of the mean model result for five-layer field resistivity sounding data.

The symbol “–” in the table denotes no information. VFSA stands for very fast simulated annealing.

To reduce uncertainty and increase the resolution of the model, model parameters containing a posterior PDF with a CI greater than 68.27 % are selected. In Table 6, the true model lies within the uncertainty range of the hybrid vPSOGWO, whereas GWO and PSO have failed to keep the true model within its uncertainty range in the second-, third-, and fourth-layer parameters. In the case of ridge regression, the uncertainty level of the model parameters is too high. For example, in the case of the third layer, both resistivity and thickness have an uncertainty approx. 44 times higher than the actual value.

The inverted 10 000 models are also computed in this example to find the posterior PDFs and histograms for each parameter. The peak of the posterior PDF is roughly close to the actual solution, as shown in the histograms in Fig. 10a and b, which reveal that ρ₂ and h₂ have a broader range that signifies the equivalence problem associated with the resistive layer. The uncertainty in each algorithm is found to be large considering all of the accepted models. However, picking the models with a greater posterior PDF CI than 68.27 % reduces the uncertainty in the model and increases the resolution of a solution.

The correlation plot between model parameters (off-diagonal) and the posterior PDF curve (diagonal) for models with a CI greater than 68.27 % for all parameters is shown in Fig. 11. There are also no significant error differences between the computed and observed apparent resistivity data for all three optimization algorithms.

The correlation matrix of a four-layer model of synthetic resistivity data is shown in Table 7. It illustrates that the first-layer parameters are correlated by a correlation matrix of 0.7383. A strong negative correlation was found between the second-layer parameters (−0.9798), and the third-layer parameters are strongly correlated with each other by a positive correlation matrix of 0.9983. Figure 7a shows the fitness between four-layer synthetic (*) and computed apparent resistivity data obtained for vPSOGWO, GWO, and PSO. The difference in fitness curves for all three optimization techniques cannot be determined, as the observed error is significantly negligible. However, the error difference can be observed in the 1D resistivity–depth models obtained from the mean model with a PDF CI of 68.27 %, as shown in Fig. 7b. Table 6 shows the mean model with a posterior PDF with a CI greater than 68.27 % for all accepted parameters in the four-layer Earth model case. In this example, the computational time required for one run with 27 data points for vPSOGWO, GWO, and PSO is 1.94, 1.84, and 1.85 s (PSO), respectively.

The optimization techniques are also executed using the same four-layer model of synthetic data with 10 % Gaussian noise and keeping the search range shown in Table 6. The same procedure is applied to determine the mean model from all of the best-fitted models and models of a posterior PDF with a CI greater than 68.27 % for all model parameters presented in Table 8. Although 10 % noise is added, the result obtained from the mean model for the posterior PDF CI of 68.27 % for the hybrid algorithm is not much different from actual values. At the same time, the experimental error is $\mathrm{3.831}\times {\mathrm{10}}^{-\mathrm{4}}$, $\mathrm{3.831}\times {\mathrm{10}}^{-\mathrm{4}}$, and $\mathrm{3.870}\times {\mathrm{10}}^{-\mathrm{4}}$ for vPSOGWO, GWO, and PSO, respectively.

Table 9 illustrates the correlation matrix of the hybrid algorithm, which clearly describes interdependence, as seen by values of 0.7644, −0.9665, and 0.9980 for the first- and second-, and third-layer parameters, respectively. Similarly, we also find a relation between second-layer resistivity and first-layer thickness (0.3605) and the resistivity of the fourth layer and thickness of the third layer (0.0549). This means that the outcome from the advocated proposed is very close to the real model.

6.3 Example 3: field data – three-layer case

For the first field example, we used a one three-layer case of vertical electrical resistivity sounding data measured with a Schlumberger array over Mount Turner, northern Queensland, Australia, interpreted by Dixon and Doherty (1977, Fig. 2a), as shown in Fig. 12a. After selecting a suitable search range, three novel algorithms, namely, vPSOGWO, GWO, and PSO, are executed to reconstruct the model interpreted by Dixon and Doherty (1977). The search range and the comparison of the proposed algorithms with the previous result (Dixon and Doherty, 1977) are presented in Table 10. Our results (for a 68.27 % CI) are closed to the development given by Dixon and Doherty (1977). The convergent error for the best-fitted model in vPSOGWO is $\mathrm{3.681}\times {\mathrm{10}}^{-\mathrm{4}}$, GWO is $\mathrm{3.697}\times {\mathrm{10}}^{-\mathrm{4}}$, and PSO is $\mathrm{3.682}\times {\mathrm{10}}^{-\mathrm{4}}$.

Table 11 presents a correlation matrix that shows a negative correlation between the first-layer parameters and a positive correlation between the second-layer parameters. A positive correlation is also observed between ρ₃ and h₂, which maintain the same model data. Figure 12a shows the apparent resistivity curve, while the 1D model obtained from the mean model with a 68.27 % CI result is shown in Fig. 12b. In this example, the computational time required for one run with 14 data points is 0.90, 0.83, and 0.78 s for vPSOGWO, GWO, and PSO respectively.

6.4 Example 4: field data – five-layer case

We selected another field example using a vertical electrical resistivity sounding dataset from Keshiari (Kharagpur), West Bengal, India, as a five-layer Earth subsurface model to determine the aquifer zone (Panda et al., 2018, their Fig. 3). This area is covered with different geological units, such as laterite, clay, and sand, and laterite material restricts the aquifer recharge process and is the most problematic area with respect to groundwater potential. We inverted these data for a five-layer Earth structure parameter using the vPSOGWO, GWO, and PSO inversion algorithms. The results are shown in Table 12 for the available models, borehole samples, and the search space for vPSOGWO, GWO, and PSO. The computed apparent resistivity curve for all three algorithms (–) and the field data (indicated by asterisks) are shown in Fig. 13a. Their error differences are significant (Fig. 13a, Table 12). The inverted 1D layered model using all algorithms obtained from the mean model with a 68.27 % CI is shown in Fig. 13b. In this example, the computational time for one run with 28 data points is 2.55, 2.43, and 2.45 s for vPSOGWO, GWO, and PSO, respectively.

Table 13 Correlation matrix using a 68.27 % PDF limit for five-layer field resistivity sounding data.

The result obtained from the mean solution of all accepted solutions and solutions with a PDF CI greater than 68.27 % aimed at all parameters using the developed techniques is presented in Table 12. The final mean models are comparable with lithological data from a 100 m deep tube well near VES6. The convergent error for vPSOGWO, GWO, and PSO is $\mathrm{4.498}\times {\mathrm{10}}^{-\mathrm{4}}$, $\mathrm{4.541}\times {\mathrm{10}}^{-\mathrm{4}}$, and $\mathrm{4.566}\times {\mathrm{10}}^{-\mathrm{4}}$, respectively, whereas the error for very fast simulated annealing (VFSA) obtained by Panda et al. (2018) is $\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{2}}$. The correlation matrix clearly shows a strong relation between the parameters of the first layer (−0.9736), the second layer (0.8434), and the third layer (−0.9907) as well as a moderate relation between the parameters of the fourth layer (0.5653). We have noticed a moderate interdependence of ρ₃ with h₂ and of ρ₅ with h₄, which can be seen in Table 13.

Table 14 Optimization of the mean model result for six-layer field resistivity sounding data.

6.5 Example 5: field data – six-layer case

As a six-layer case study, we again applied the vPSOGWO, GWO, and PSO algorithms to invert field apparent resistivity data extracted near a borehole in Apulia, South Italy, for hydrogeological purposes (Sen et al., 1993). The search range was taken from Sen et al. (1993), but the fourth-layer and upper-bound thickness of the fifth layer increase by 50 m, as shown in Table 14. The reproduced field data (*) and inverted field data (–) are shown in Fig. 14a. The misfit error obtained is $\mathrm{2.830}\times {\mathrm{10}}^{-\mathrm{4}}$, $\mathrm{3.243}\times {\mathrm{10}}^{-\mathrm{4}}$, and $\mathrm{3.133}\times {\mathrm{10}}^{-\mathrm{4}}$ for vPSOGWO, GWO, and PSO, respectively, whereas the error using simulating annealing (SA) in Sen et al. (1993) is 0.017. Table 14 also includes the mean model for a CI of 100 % and 68.27 % using the proposed algorithms and previously published literature. It is observed that few parameters obtained fall within the uncertainty of corresponding parameters of vPSOGWO. The vPSOGWO-inverted results provide higher similarity with the borehole information than the results from SA (Sen et al., 1993). The interdependence between the layer parameters can be seen from the correlation matrix shown in Table 15. A strong correlation among parameters of the first layer (0.8211), the second layer (−0.9327), and the third layer (0.9766) is shown by the correlation matrix, which is comparable to the correlation matrix that has been presented in Sen et al. (1933, their Table 13). A moderate correlation between fourth-layer (−0.5246) and fifth-layer parameters (0.4486) is also observed. It is also notable that there is a sensible relation between sixth-layer resistivity and fifth-layer thickness, keeping the same model data.

The error differences between computed data and observed data are significant, as shown in Fig. 14a and Table 12. The inverted 1D layered models obtained from the mean model with a CI of 68.27 % are shown in Fig. 14b. In this example, the computational time for one run with 28 data points is 3.58, 3.44, and 3.45 s for vPSOGWO, GWO, and PSO, respectively. The inverted results from vPSOGWO, GWO, and PSO have been shown along with the borehole data and published results (Sen et al., 1993) in Table 14. In addition to published findings, it is notable that the outcomes obtained from the hybrid algorithm (vPSOGWO) satisfy the borehole information very well compared with the other algorithms.

Table 15 Correlation matrix using a 68.27 % PDF limit for six-layer field resistivity sounding data.

7 Conclusion

We conducted an assessment of three metaheuristic algorithms, namely PSO, GWO, and vPSOGWO, to determine their effectiveness and suitability with respect to solving geoelectrical inverse issues. These challenges include the estimation of 1D resistivity models based on geoelectrical resistivity sounding data. The relevance of these algorithms was validated using synthetic and field resistivity sounding data signifying the different kinds of Earth subsurface stratigraphy. An enormous solution (100 000 000 from 10 000 runs) was assessed. Subsequently, the best-fitted solutions were chosen within a predefined value for statistical measurements. The statistical study, including a posterior PDF with a 68.27 % CI, a mean solution, a posterior solution correlation matrix, and a covariance matrix using search space, was carried out to refine the solutions and obtain the global mean solution with the least uncertainty. These statistical simulations yield essential information as to the reliability of an inversion algorithm. In general, conventional techniques can be quite effective in resolving the model with respect to random noise but can fail with respect to systematic error and inappropriate models. Our investigation with the application of the developed algorithm, including the statistical simulation for different multilayer resistivity parameters, resulted in a quantitative appraisal of uncertainty in the derived model parameters. We observed that the output of the hybrid algorithm in terms of the mean model or error might be similar to either PSO or GWO (attributed to the exploration characteristics of GWO and the exploitation characteristics of PSO). The vPSOGWO, GWO, and PSO algorithm performance has been analyzed based on the uncertainty and stability and on the mean model of the layered Earth structure. We found that vPSOGWO gives much closer results than the results inverted from other two algorithms and than conventional methods and that these results are consistently better than the previously published results and are correlated well with borehole information.