Application of particle swarm optimization for gravity inversion of 2.5-D sedimentary basins using variable density contrast

Particle swarm optimization (PSO) is a global optimization technique that works similarly to swarms of birds searching for food. A MATLAB code in the PSO algorithm has been developed to estimate the depth to the bottom of a 2.5-D sedimentary basin and coefficients of regional background from observed gravity anomalies. The density contrast within the source is assumed to vary parabolically with depth. Initially, the PSO algorithm is applied on synthetic data with and without some Gaussian noise, and its validity is tested by calculating the depth of the Gediz Graben, western Anatolia, and the Godavari sub-basin, India. The Gediz Graben consists of Neogen sediments, and the metamorphic complex forms the basement of the graben. A thick uninterrupted sequence of Permian–Triassic and partly Jurassic and Cretaceous sediments forms the Godavari sub-basin. The PSO results are better correlated with results obtained by the Marquardt method and borehole information.


Introduction
The gravity method is a natural source method which helps in locating masses of greater or lesser density than the surrounding formations.It is used as a reconnaissance survey in hydrocarbon exploration.Sedimentary basins, which are characterized by negative gravity anomalies, are the location of almost all of the world's hydrocarbon reserves.Interpretation of gravity data is a mathematical process of trying to optimize the parameters of the initial model in order to get a good match to the observed data.Interpretation of gravity data is always associated with the ambiguity problem.Ambiguity in gravity anomalies can be over-come by assigning a mathematical geometry to the anomalycausing body with a known density contrast (Rama Rao and Murthy, 1978).Bott (1960) used stacked prism model to describe the cross-section of a sedimentary basin, whereas Talwani (1959) used the polygonal model to describe source geometry.The parabolic density function is used to remove the complications associated with exponential (Cordell, 1973), cubic (Garcia-Abdeslem, 2005) and quadratic (Gallardo-Delgado et al., 2003) density functions.The Marquardt inversion through residual gravity anomalies delineates the structure of a sedimentary basin by estimating regional background (Chakravarthi and Sundararajan, 2007;Marquardt, 1963).Several authors have developed 2-D/2.5-Dlocal optimization techniques over the 2-D/2.5-Dsedimentary basin (Annecchione et al., 2001;Barbosa et al., 1999;Bhattacharya and Navolio, 1975;Gadirov et. al, 2016;Litinsky, 1989;Morgan and Grant, 1963;Murthy et al., 1988;Murthyan and Rao, 1989;Rao, 1994;Won and Bavis, 1987) to interpret gravity anomalies with constant density function.In many publications over 3-D gravity field computation with an approximation of geological bodies by 3-D polygonal horizontal prism has been applied (Eppelbaum and Khesin, 2004;Khesin et al. 1996).Rao (1990) used a quadratic density function, which is comparatively reliable, to analyse gravity anomalies over basins having a limited thickness, whereas Chakravarthi and Rao (1993) carried out modelling and inversion of gravity anomalies with a parabolic density function.
Particle swarm optimization (PSO) is a robust stochastic optimization technique based on the movement and intelligence of swarms, which was developed by Kennedy and Eberhart (1995).PSO applies the concept of social optimization in problem solving in various fields.In this paper, a MATLAB code based on PSO is developed to interpret the gravity anomalies of 2.5-D sedimentary basins, where the density varies parabolically with depth.PSO-analysed results are consistent with and more accurate than other techniques and also agree significantly well with borehole information.

Theory
In Bott's approach, a sedimentary basin is approximated by a series of vertical prisms.The gravity anomaly g b at any point on the profile AA as shown in Fig. 1.
The gravity effect of lth prism is given as (2) The parabolic density function at any depth w is given by Finally, after integration of Chakravarthi and Sundararajan ( 2006), Eq. (2) becomes , where Here, N is the number of observations, G is the universal gravitational constant, C and D are coefficients of regional background, c is half width of the prism, z 1 and z 2 are depths to the top and bottom of the basin, 2L is strike length of the prism, a is the offset of the profile from the centre of the prism and ρ 0 and α are constants of the parabolic density function at depth z.
Since the profile AA does not pass through the centres of each prism, Eq. ( 4) has to be calculated twice by putting L−a and L + a for L and taking the average.The initial depths of the basin calculated using observed anomaly g 0 are given as The profile AA entirely covers the lateral dimensions of the sedimentary basin; therefore the depth of the basin on either side of the profile become zero.So,

Particle swarm optimization
PSO uses a number of particles (solutions) that constitute a swarm moving around in the search space looking for the best solution.Each particle adjusts its "flying" according to its own flying experience as well as the flying experience of other particles.Each particle keeps track of its coordinates in the solution space which is associated with the best solution (fitness) that has been achieved so far by that particle.This value is called personal best, P best.Another best value that is tracked by PSO is the best value obtained so far by any particle in the neighborhood of that particle.This value is called global best, Gbest.The basic idea of PSO lies in accelerating each particle towards its P best and the Gbest locations with a random weighted acceleration at each time step (Mohapatra and Das, 2013).
where k is the number of iterations; t is the particle number; V k t is the velocity of the tth particle at k iterations; X k t is the position of tth particle at k iterations; P best t is the best position of individual tth particle (local best position); Gbest is the best position of all particles (global best position); rand 1 and rand 2 are the independent uniformly random numbers in the range [0, 1]; c 1 and c 2 are the positive learning factor which controls the maximum step length; and w is the inertial weight factor that controls the speed of the particles.Equation ( 7) gives the updated velocity based on the current velocity, current position, local, best position and global best position.This process is repeated until the desired result is obtained.The schematic diagram/flow chart of the PSO algorithm is shown in Fig. 2.

Examples
The MATLAB code based on PSO is applied to a synthetic model of a sedimentary basin and real field data sets from Gediz Graben, western Anatolia, and the Godavari sub-basin, India.

Synthetic Example
We have used a synthetic gravity anomaly of 45 × 10 3 m length at a 1.0 × 10 3 m station interval.In Bott's algorithm, the prism will be of equal width, 1.0 × 10 3 m, but with different strike lengths.Here parabolic density function is used with the constants ρ 0 = −0.65×10 3 Kg m −3 and α = 0.04 × 10 3 Kg m −3 per 1000 m.The profile does not bisect the strike lengths of prism, and so offset distance of the profile from the centre of each prism is mentioned in the code.We have added an interference term, Ax + B, with A = −0.007mgal per 1000 m and B = −10 mgal for the regional background.The required result is found at 15 iterations with a root mean square error (RMSR) of 2.9369e−006 from the Marquardt algorithm.
We have used the same synthetic gravity anomaly for the PSO algorithm.The Fig. 3 shows the learning process of P best and Gbest in terms of error and iterations.The best result is found with 57 iterations and 50 models (Fig. 3).So it is seen that after 57 iterations and 50 models, the calculated anomalies match the synthetic anomaly and estimated depths coincide with the actual structure where the RMSE is 2.8383e − 004.Gaussian noises of 5 and 10 % are added to the synthetic data to perceive the robustness of the PSO algorithm.PSO does not find the true depths, but it gives values close to the true depths.The upper part of Fig. 4 shows the synthetic and PSO-calculated gravity anomalies of a synthetic model of a 2.5-D sedimentary basin, and the lower part shows the inferred depth structure obtained from the PSO and Marquardt algorithm for synthetic data.Figures 5 and 6 show the synthetic data with 5 and 10 % Gaussian noises and calculated gravity anomalies obtained from the PSO algorithm, and inferred depth structure obtained by the PSO and Marquardt algorithm.
3 Field example

Gediz Graben, western Anatolia
The first field case study of the interpretation of gravity anomalies has been taken from Gediz Graben, western Anatolia.The PSO technique has been applied using 29 vertical prisms, each with equal width of 2.0 × 10 3 m but with dif-  ferent strike lengths, whereas Chakravarthi and Sundararajan (2007) used the same prism and interpreted gravity anomaly by the Marquardt algorithm using a parabolic density function whose constants are ρ 0 = −1.407× 10 3 and α = 2.26935 × 10 3 Kg m −3 per 1000 m.
We have used a similar number of prisms in PSO to improve the results.So with 65 iterations and 50 models, we achieve a good fit between observed and PSO-analysed gravity anomalies with a RMSE of 0.0083.The maximum thickness of the graben is inferred as 1.87×10 3 m, which matches well with 1.8×10 3 m as estimated by Sari and Salk (2002), as compared to 1.64×10 3 m obtained by Chakravarthi and Sundararajan (2007).The upper part of Fig. 7 shows the observed and PSO-calculated gravity anomalies over Gediz Graben, western Anatolia, and the lower part show the inferred depth structure obtained from the PSO and Marquardt algorithm.

Godavari sub-basin
The Godavari sub-basin is one of the major basins of the Pranhita-Godavari valleys (Rao, 1982), whose approximate strike length is 220 × 10 3 m.The gravity profile is taken for  study from the residual Bouguer gravity anomaly map of the Godavari sub-basin as shown in Fig. 8.We have used 29 vertical prisms, each with equal widths of 2.0 × 10 3 m but with different strike lengths for sedimentary basin modelling.The constants of parabolic density functions used for the Godavari sub-basin are ρ 0 = −0.5×10 3 and α = 0.1518259× 10 3 Kg m −3 per 1000 m (Chakravarthi and Sundararajan, 2004).So with 71 iterations and 45 models, we achieve a good fit between observed and PSO-analysed gravity anomalies.The RMSE is 0.0099.The maximum depth of the basin, obtained from PSO, is 4.09 × 10 3 m, which is quite close to the borehole information (Agarwal, 1995).Chakravarthi and Sundararajan (2005) obtained a maximum depth of 4.0× 10 3 m, whereas Ramanamurty and Parthasarathy (1988) suggested 4.5 × 10 3 m as the thickness of the basin.The up-

Conclusions
Particle swarm optimization in the MATLAB environment has been developed to estimate the model parameters of a 2.5-D sedimentary basin where the density contrast varies parabolically with depth.We have implemented the PSO algorithm on synthetic data with and without Gaussian noise and two field data sets.An observation has made that PSO is affected by some levels of noise, but estimated depths are close to the true depths.The results obtained from PSO using synthetic and field gravity anomalies are well correlated with borehole samples and provide more geological viability.Despite its long computation time, PSO is very simple to implement and is controlled by only one operator.

Figure 1 .
Figure 1.The 2.5-D sedimentary basin and its approximation by juxtaposing prisms.

Figure 2 .
Figure 2. The detail schematic diagram/flow chart of PSO techniques.

Figure 3 .
Figure 3. Iteration verses RMSE of P best and Gbest using the PSO technique through synthetic gravity anomaly.

Figure 4 .
Figure 4. Synthetic and calculated gravity anomalies with parabolic density function due to a synthetic model of a 2.5-D sedimentary basin, obtained from the PSO algorithm, and inferred depth structure obtained from the PSO and Marquardt algorithm.

Figure 5 .
Figure 5. Synthetic data with 5 % Gaussian noise and calculated gravity anomalies obtained from the PSO algorithm, and inferred depth structure obtained from the PSO and Marquardt algorithm.

Figure 6 .
Figure 6.Synthetic data with 10 % Gaussian noise and calculated gravity anomalies obtained from the PSO algorithm, and inferred depth structure obtained from the PSO and Marquardt algorithm.

Figure 7 .
Figure 7. Observed and calculated gravity anomalies with parabolic density function, Gediz Graben, western Anatolia, obtained from the PSO algorithm, and inferred structure obtained from the PSO and Marquardt algorithm.

Figure 8 .
Figure 8. Residual Bouguer gravity anomaly map of the Godavari sub-basin (modified after Chakravarthi and Sundararajan, 2005) and gravity anomaly profile taken for study.

Figure 9 .
Figure 9. Observed and calculated residual Bouguer gravity anomalies of parabolic density function of the Godavari sub-basin obtained from the PSO algorithm, and inferred depth structure from the PSO and Marquardt algorithm.