Inversion of residual gravity anomalies using tuned PSO

Many kinds of particle swarm optimization (PSO) techniques are now available and various efforts have been made to solve linear and non-linear problems as well as onedimensional and multi-dimensional problems of geophysical data. Particle swarm optimization is a metaheuristic optimization method that requires intelligent guesswork and a suitable selection of controlling parameters (i.e. inertia weight and acceleration coefficient) for better convergence at global minima. The proposed technique, tuned PSO, is an improved technique of PSO, in which efforts have been made to choose the controlling parameters, and these parameters have been selected after analysing the responses of various possible exercises using synthetic gravity anomalies over various geological sources. The applicability and efficacy of the proposed method is tested and validated using synthetic gravity anomalies over various source geometries. Finally, tuned PSO is applied over field residual gravity anomalies of two different geological terrains to find the model parameters, namely amplitude coefficient factor (A), shape factor (q) and depth (z). The analysed results have been compared with published results obtained by different methods that show a significantly excellent agreement with real model parameters. The results also show that the proposed approach is not only superior to the other methods but also that the strategy has enhanced the exploration capability of the proposed method. Thus tuned PSO is an efficient and more robust technique to achieve an optimal solution with minimal error.


Introduction
The gravity method is based on the measurement of gravity anomalies caused by the density variation due to source anomalies.The gravity method has been used in a wide range of applications as a reconnaissance method for oil exploration and as a secondary method for mineral exploration, to find out the approximate geometry of the source anomalies, bedrock depths and shapes of the earth.The interpretation of geophysical data involves solving an inverse problem; many techniques have been developed to invert the geophysical data to estimate the model parameters.These methods can be broadly categorized into two groups: (1) local search techniques (e.g.steepest descent method, conjugate gradient method, ridge regression, Levenberg-Marquardt method) and (2) global search techniques (e.g.simulated annealing, genetic algorithms, particle swarm optimization, ant colony optimization).Local search techniques are simple and require a very good initial presumption -one that is close enough to the true model for a successful convergence.On the other hand, global search methods may provide an acceptable solution but are computationally time intensive.Several local and global inversion techniques have been developed to interpret gravity anomalies (Thanassoulas et al., 1987;Shamsipour et al., 2012;Montesinos et al., 2005;Qiu, 2009;Toushmalani, 2013).However, PSO has been successfully applied to many fields, such as model construction, biomedical images, electromagnetic optimization, hydrological problems, etc. (Cedeno and Agrafiotis, 2003;Wachowiak et al., 2004;Boeringer and Werner, 2004;Kumar and Reddy, 2007;Eberhart and Shi, 2001;El-Kaliouby and Al-Garni, 2009) but in the geophysical field PSO has a limited number of applications (Alvarez et al., 2006;Shaw and Srivastava, 2007).
Published by Copernicus Publications on behalf of the European Geosciences Union.R. Roshan and U. K. Singh: Inversion of residual gravity anomalies using tuned PSO In this paper, improved particle swarm optimization, known as tuned PSO, has been discussed.This PSO method is a global optimization technique, has artificial intelligence, preserves its past experience to avoid trapping in local minima and converges at global minima.This technique has a very good exploration capability (searching capability) and global convergence capability (ability of finding the optimal solution).These capabilities are modulated by inertia weight (w) and acceleration coefficients.Therefore, the selection of suitable controlling parameters play a vital role in its performance, enabling it to achieve an optimal solution.In traditional PSO, controlling parameters (i.e.w, c 1 and c 2 ) are generally taken as 0.9, 1.4 and 1.4 for reasonable results (Kennedy, 1999;Das et al., 2008).However, controlling parameters in PSO analyses are usually data dependent; hence no unique architecture can be generalized for every application.In this paper PSO parameters (i.e.w, c 1 and c 2 ) are tuned in several exercises to find the best tuning of the controlling parameters.They are demonstrated using synthetic gravity anomalies over various geometrical bodies and their efficacy is compared.On the basis of performance the method is finally applied to field gravity anomalies to compute the essential model parameters such as amplitude coefficient factor (A), shape factor (q) and depth (z) and horizontal location(x 0 ) of the source geometry.Thus, the derived method provides a more accurate, consistent and stable solution than conventional PSO and other methods.

Forward modelling for generating the synthetic gravity anomalies
A general expression of gravity anomaly caused by a sphere, an infinite long horizontal cylinder and a semi-infinite vertical cylinder have been used for generating the gravity anomalies in a forward problem that is given in Eq. (1) (Abdelrahman et al., 1989) as follows: where 3 2 for a sphere, 1 for a horizontal cylinder, 1 2 for a vertical cylinder; R z.

,
where A, q and z represent amplitude coefficient factor, shape factor and depth and x i , σ , G and R are the position coordinate, density contrast, universal gravitational con-stant and radius of geometrical bodies.The unit of amplitude coefficient factor (A) changes as the geometry of the model changes, depending on shape factor (q) and exponent m.To avoid the confusion of using a unit of amplitude coefficient factor we use expression A (mGal km 2q−m ).From the expression it is clear that the unit of amplitude coefficient factor (A) for a spherical shape is mGal km 2 and in case of a cylinder mGal km.mGal is the traditional unit of gravity anomaly, mainly used in gravity surveys.The SI unit of Gal is 0.01 m s −2 while the shape factor is dimensionless.Two types of synthetic gravity data have been created over spherical and vertical cylindrical geometrical models using forward modelling of the Eq. ( 1).The value of parameters for the spherical model, i.e. amplitude coefficient factor (A), shape factor and depth have been taken as 600 mGal km 2 , 1.5 and 5.0 km respectively.Similarly, the values of parameters for the vertical cylindrical model (A = 200 mGal km, q = 0.5 and z = 3.0 km) have been selected.The shape factor approaches zero as the structure becomes a near-horizontal bed and approaches 1.5 as the structure becomes a perfect sphere (point mass).As in the formulae x i is the position coordinate; at the origin x i = 0 Eq. ( 1) then becomes (2) 10 % white Gaussian noise is added to synthetic gravity anomalies using the following equation.
g noisy (x) = awgn(g(x), 0.1), where awgn is additive white Gaussian noise, g(x) represents gravity anomaly value in which noise is to be added and 0.1 is the magnitude of noise at 10 %.
3 Tuned particle swarm optimization (tuned PSO) Tuned particle swarm optimization (tuned PSO) is an improved particle swarm optimization (PSO) method, named after the fine tuning of its learning parameters.The PSO technique is an evolutionary computational technique inspired by the social behaviour of the particles (Eberhart and Kennedy, 1995).Each particle as a potential solution of the problem knows its best values (P best ) and its position.Moreover, each particle knows the best value in the group (G best ) among the P best .All of the best values are based on the objective function (Q) so that each problem can be solved.Each particle tries to modify its position in the current velocity.The velocity of each particle can be updated using the following equations (Santos, 2010): where v k i is the velocity of the ith particle at the kth iteration, x k i represents the current position of the ith particle at the kth iteration, rand( ) is a random number in the range of 0, 1. c 1 and c 2 are constants known as the cognitive coefficient and the social coefficient respectively and ω is an inertia weight.
The objective function is calculated by the following equation (Santos, 2010).
where N is the number of iterations, and v o i and v c i are observed and calculated gravity anomalies measured at point p(x i ).

Selection of learning parameter for tuned PSO modelling
In this paper, a judicious selection of the parameters (i.e.ω, c 1 , and c 2 ) has been discussed for controlling the convergence behaviour of the tuned-PSO-based algorithm.The settings of these parameters determine how it optimizes the search space.These algorithms with a suitable selection of parameters become more powerful for their practical applications.

Inertia weight
Inertia weight (ω) controls the momentum of the particle (Eberhart and Shi, 2001;Eberhart and Kennedy, 1995).Here, two kinds of source geometry are adopted to evaluate more suitable ranges of parameters in the tuned PSO.For tuning the inertia weight, 0.4, 0.7 and 0.9, have been taken for two different acceleration coefficients at 1.4 and 2.0.From Fig. 1 and Table 1, it is clear that the best convergence is performed by the algorithm at inertia weight 0.7.This value of inertia weight produces a high convergence rate with a smaller number of iterations than the other values.

The maximum velocity
Maximum velocity (V max ) describes the maximum number of changes of position coordinates that can take place during each iteration.The concept of V max was introduced to avoid explosion and divergence (Das et al., 2008).Generally, a full search range is set as V max for the particles' positions.For example, if a particle has position vector x = (x 1 , x 2 , x 3 ) and if −n ≤ x i ≤ n where i = 1, 2, 3 then the maximum velocity becomes 2n.In our case, the range of particle position −15 ≤ x i ≤ 15 is more appropriate, so we have taken 30 as the value of V max .

The swarm size
It is quite a common practice in the PSO literature to limit the range of the number of particles.Van den Bergh and Engel- brecht ( 2001) have shown that, although there is a slight improvement of the optimal value with increasing swarm size, a larger swarm increases the number of function evaluations to converge to an error limit.However, Eberhart and Shi illustrated that the population size has hardly any effect on the performance of the PSO method.Therefore, in this paper, population size is set at 100.

The acceleration coefficients
The acceleration coefficients are the learning coefficients which provide stability for the exploration of the particle.
There are two kinds of acceleration coefficients: (i) the cognitive coefficient c 1 , which contributes towards the self exploration of a particle and (ii) the social coefficient c 2 , which contributes towards the motion of the particles in a global direction.To find the best tuning of learning parameters, various values of c 1 , c 2 (i.e.c 1 = c 2 = 1.0, 1.2, 1.4, 1.6, 1.8, and 2.0) and inertia weights (i.e.0.4, 0.7 and 0.9) are taken, and various exercises have been made using the two different geometrical bodies by fixing each of the inertia weights (Table 1).The results were analysed and it was found that the values of c 1 and c 2 (i.e.c 1 = c 2 = 1.4) are the best-tuned acceleration coefficients for our case.These values of acceleration coefficients have been used to invert the gravity anomalies, which provides significant improvement and produces optimal solutions for the geological bodies.

Discussion and results
Initially the set of the controlling parameters (w, c 1 , c 2 ) were taken as (0.4, 1.2, 1.2), (0.4, 1.4, 1.4), (0.4, 1.6, 1.6), (0.4, 1.8, 1.8) and (0.4, 2.0, 2.0).The rms error was observed during the examinations of the sensitivity of the parameters and the applicability of the proposed algorithm using each set of pa-Table 1. Performance of the acceleration coefficients c 1 and c 2 using the synthetic gravity anomalies over spherical and vertical cylindrical geometrical bodies.

Gravity data
Weighting rms error Synthetic w = 0.4 0.004899 0.002899 0.00014 0.000907 0.000853 0.000861 spherical body w = 0.7 0.002532 0.000118 0.000013 0.000087 0.000187 0.000247 w = 0.9 0.005215 0.000118 0.000063 0.000379 0.000167 0.002695 Synthetic w = 0.4 0.004892 0.003231 0.000327 0.000835 0.000704 0.000932 vertical w = 0.7 0.001913 0.000318 0.000011 0.000065 0.000207 0.000511 Cylindrical body w = 0.9 0.003259 0.000551 0.000189 0.000183 0.001265 0.002747 rameters.Similarly, the same observations were made after replacing the inertia weight of 0.4 with 0.7 and 0.9, keeping the same values of acceleration coefficients (c 1 and c 2 ) as shown in Table 1.We analysed Table 1 and suggest that the value of the rms error is at a minimum for the set of parameters (0.7, 1.4, 1.4) in comparison to other sets of parameters.In addition, the algorithm using the sets of parameters (0.7, 1.4, 1.4) has a lower number of local minima and faster convergence than the other (Fig. 1).The proposed technique using tuned parameters was demonstrated on synthetic and field residual gravity anomalies to find the amplitude coefficient factor (A), shape factor (q) and depth (z) of various geometrical bodies of different geological terrains.

Application to synthetic gravity anomalies
Two geometrical models, i.e. sphere and vertical cylinder, have been considered for testing the applicability and efficacy of tuned PSO.The synthetic gravity anomalies over the above-considered models are generated from Eq. ( 1).In addition, other data sets (noisy synthetic gravity anomalies) are also generated with 10 % Gaussian noise to perceive the efficacy of the proposed algorithm.In each case, the gravity profile length is 51 km and the data points are kept at equal intervals of 1 km.The proposed tuned PSO algorithm has been applied to the above-mentioned synthetic data sets.The optimized results obtained by tuned PSO for synthetic data without noise have been shown in Table 1a and Table 2a.
Similarly, the results for synthetic data with noise have been shown in Table 1b and Table 2b.
Figure 1 shows the iteration versus error.It suggests that the error is at a minimum and has a lower number of local minima at values of controlling parameters c 1 = 1.4,c 2 = 1.4 and w = 0.7.It means that the tuned PSO technique minimizes the number of local minima for solving the geophysical non-linear inverse problems.The synthetic gravity anomaly without noise and computed gravity anomaly by tuned PSO are shown in Figs.2a and 3a respectively.Similarly, the synthetic gravity anomaly with noise and computed gravity anomaly by tuned PSO are shown in Figs.2b and 3b.Figures 2a and 3a show that the calculated gravity anomaly curves by tuned PSO are matched well with the synthetic gravity anomaly curves for spherical and vertical cylindrical models respectively.The behaviour of p best and g best are shown in Fig. 4 and suggests that the error for g best decreases more rapidly with a high convergence rate.
Tables 2a and 3a show the values of the rms error using the synthetic data without noise.Also, Tables 2b and 3b show the rms error using the synthetic data with noise.The analysis of the tables reflects that the rms error is comparatively higher in the case of synthetic data with noise.However, the horizontal location (x 0 ) is a substantially stable parameter and varies on a small scale.Table 6 shows the computation time taken by the algorithm to find the solution for given number of iterations.The number of models is 100 in the entire analysis of synthetic gravity anomalies.It is clear from the table that the algorithm is powerful and takes less computation time to produce optimal results.

Mobrun sulfide body near Rouyn-Noranda, Canada
The Mobrun polymetallic deposit near Rouyn-Noranda comprises two complexes of massive lenses within mainly felsic volcanic rocks of the Archean Blake River Group (Barrett et al., 1992).Host volcanic rocks of mainly sulfide ore bodies are mostly massive, breciated, and tuffaceous rhyolites.The Mobrun ore body is located at a shallow depth; the top of the body is at a depth of approximately 17 m and extends to 175 m (Aghajani et al., 2009).Tuned PSO in MATLAB has been applied to field residual gravity anomalies.The anomaly profile length of 268 m has been taken from the Mobrun sulfide body, Noranda, Canada (Nettleton, 1976;Essa, 2012).It is seen from Fig. 5 that both anomaly curves, i.e. analysed from tuned PSO and observed gravity anomalies, are significantly well correlated with the optimal rms error of 0.0271 %.The results in terms of model parameters (amplitude coefficient factor, shape factor and depth) over the Mobrun ore body, analysed using the tuned PSO method, can seen in Table 4a.This table provides the optimum results obtained from tuned PSO with a 0.0271 % error.This agrees well with the results obtained from other methods.The calculated value of the shape factor, q, is 0.77 (Table 4a).This value over the Mobrun sulfide ore body reflects the shape of a semi-infinite vertical cylindrical geological body present at a depth of 30 m.Since the shape factor computed by the proposed method (q = 0.77) lies between the shape factor of a perfect semi-infinite vertical cylinder, i.e. q = 0.5, and the shape factor of an infinite horizontal cylinder, i.e. q = 1.0, it can be seen from Table 4b that the values of amplitude coefficient factor, shape factor and depth correspond to 60.0, 0.77 and 30 are more accurate than the results analysed by various authors.
Table 3.(a) Optimized model parameters, converged iteration and rms error in the inversion of synthetic gravity anomaly over a vertical cylindrical source model and (b) optimized parameters, converged iteration and rms error in the inversion of synthetic gravity anomaly with 10 % white Gaussian noise over a same-source model from tuned PSO.
(a) Optimized parameters, converged iteration and rms error in the inversion of synthetic gravity anomaly over a vertical cylindrical source model.Z (km) A (mGal km) q g 0 (mGal) x 0 (km) Iteration rms error (%)  (Nettleton, 1962).In this paper, tuned PSO in MATLAB environment has also been applied to another field case study.The gravity anomaly of Louga area, west coast of Senegal, western Africa (Essa, 2014) has been chosen for tuned PSO analysis as shown in Fig. 6.It has a profile length of 32 km.The results in terms of the model parameters (amplitude coefficient factor, shape factor and depth) over the Louga anomaly analysed from tuned PSO method can be seen in Table 5a.It is seen from Fig. 6 that both gravity anomalies curves analysed with tuned PSO and observed gravity anomalies are extremely well correlated with the optimal rms error of 0.0271 %.The optimum results of the model parameters amplitude coefficient factor (A), shape factor (q) and depth (z) are 545.30mGal km, 0.53 and 4.92 km respectively.They show significantly good agreement with the results obtained by various authors as shown in Table 5b.The tuned-PSOanalysed value of the shape factor confirms that the shape of the causative body is a semi-infinite vertical cylinder present at a depth of about 4.92 km.

Conclusions
In this paper, various synthetic gravity anomalies and field gravity anomalies have been adopted to evaluate the applicability and efficacy of tuned PSO algorithms and also to determine the suitable range of settings for the learning parameters (i.e.inertia weight and acceleration coefficients).On the basis of the performance, a novel algorithm PSO with suitable learning parameters has been implemented for gravity anomalies of assuming models such as spheres and vertical cylinders.This technique has been tested and demonstrated on synthetic gravity anomalies with and without Gaussian noise and finally applied to field residual gravity anomalies over the Mobrun sulfide ore body, Noranda, QC, Canada and the Louga anomaly on the western coast of Senegal, western Africa.This technique provides robust and plausible results even in the presence of noise and is consistent with the results obtained from other classical methods.Thus, tuned PSO is a  powerful tool that improves the results of classical PSO and other techniques significantly with less time and optimal error.

Data availability
The main aim of our work deals with the applicability of the proposed algorithm tuned PSO in the geophysical potential field.Initially tuned PSO is demonstrated on synthetic data created by Eq. ( 1) in the text (Abdelrahman et al., 1989) and details about synthetic data were given in Sect.4.1.Finally, proposed algorithms were applied over two kinds of field data, taken as follows: (i) residual gravity anomaly over Mobrun field area of profile length 268 m was digitized at interval of 8.4 m from Essa ( 2012) and (ii) the Louga anomaly, west coast of Senegal, western Africa, of profile length 32 km was digitized at interval of 1.0 km from Essa (2014).
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Iteration versus rms error plot at different acceleration coefficients and inertia weights.

Figure 2 .
Figure 2. (a) Synthetic gravity anomaly versus tuned PSO calculated gravity anomaly over a spherical model and (b) synthetic gravity anomaly versus tuned PSO calculated gravity anomaly over the same model with 10 % white Gaussian noise.

Figure 3 .
Figure 3. (a) Synthetic gravity anomaly versus tuned PSO calculated gravity anomaly over a vertical cylindrical model, (b) synthetic gravity anomaly versus tuned PSO calculated gravity anomaly over the same model with 10 % white Gaussian noise.

Figure 4 .
Figure 4. Iteration versus rms error of tuned PSO showing p best and g best over synthetic gravity anomalies.

Figure 5 .
Figure 5. Observed field gravity anomaly versus tuned PSO calculated gravity anomaly over Mobrun sulfide ore body, Canada.

Figure 6 .
Figure 6.Observed field gravity anomaly versus tuned PSO calculated gravity anomaly over the western Senegal anomaly, Louga area, western Africa.

Table 2 .
(a) Optimized model parameters, converged iteration and rms error in the inversion of synthetic gravity anomaly over a spherical source model and (b) optimized parameters, converged iteration and rms error in the inversion of synthetic gravity anomaly with 10 % white Gaussian noise over a same-source model from tuned PSO.(a) Optimized parameters, converged iteration and rms error in the inversion of synthetic gravity anomaly over a spherical source model.Optimized parameters, converged iteration and rms error in the inversion of synthetic gravity anomaly with 10 % white Gaussian noise over a spherical source model.

Table 4 .
Optimized parameters, converged iteration and rms error in the inversion of synthetic gravity anomaly with 10 % white Gaussian noise over a vertical cylindrical source model.(a)Analysed results and parameters (A, z and q) used to invert the gravity anomaly over Mobrun sulfide ore body and (b) comparative results over Mobrun field, Canada from various methods and tuned PSO.(a) Optimized parameters, converged iteration and rms error in the inversion of field gravity anomaly over Mobrun sulfide ore body.

Table 5 .
(a) Analysed results and parameters (A, z and q) used to invert the gravity anomaly over the western Senegal anomaly, Louga area, western Africa and (b) comparative results over same area from various methods and tuned PSO.(a) Optimized parameters, converged iteration and rms error in the inversion of field gravity anomaly over western Senegal (Louga area) anomaly.

Table 6 .
Synthetic source models, iterations and computation time (in seconds).