Hunger Games Search for the elucidation of gravity anomalies with application to geothermal energy investigations and volcanic activity studies

1 Introduction

The gravity method in geophysics, a traditional, passive, non-invasive, and cost-efficient method, measuring the variation in the Earth’s gravitational field caused by the difference of rock mass properties at a certain position, has been applied in a large spectrum of projects, including ore/mineral exploration [1,2,3,4,5], archaeology investigation [6,7,8,9], radioactive waste management [10], weapons and unexploded ordnances detection [11,12], hydrocarbon/reservoir (oil and gas) recognition [13,14,15], geothermal energy and volcanic activity study [16,17,18], underground cavity or tunnel identification [19,20,21], geotechnical and engineering application [22,23], groundwater monitoring [24,25], landfill mapping [26,27], and subsurface structure imaging [28,29,30]. It obtains the same quality as other constantly applied methods, which can remotely sense potential targets in terms of capturing their features and variations [31,32,33].

Especially, in geothermal energy and volcanic activity studies, the gravity method enables insights into the geometrical information of the causative sources, which are crucial for targeting the potential geothermal reservoirs [34] and clarify the geodynamic process based on the spatial and temporal volcanic variations, which have got up from earthquakes and interior movements, respectively [15,35]. Intrinsically, the goal of delineating the measured gravity response is to estimate the model parameters (amplitude coefficients, depths, and spatial locations) of buried sources from a precarious process due to the non-unique and ill-posed mathematical nature of the gravity inversion problem [36,37,38,39], namely, that numerous configurations are capable of reproducing observations within a given error tolerance. Such gravity structures can be considered as spherical, cylindrical, dike-like, fault-like structures, and sheets [18,40,41]. For these structures, versatile quantitative techniques have been proposed and implemented for better interpretation performances, which can be categorized as follows:

1. Conventional and non-conventional methods, which rely on the characteristic points, nomograms, window and depth curves, matching curves, various transformations, gradient-based minimization techniques, moving average, Euler deconvolution, depth to extreme point, and fair function [42,43,44,45,46,47,48,49,50,51,52].

2. Remarkably, nature-inspired derivative-free metaheuristic algorithms have been effectively applied with proper modifications, including particle swarm optimization (PSO), genetic algorithm, differential evolution, differential search, simulated annealing, ant colony optimization, gravitational search algorithm, bat algorithm (BA), manta-ray foraging optimization, and modified barnacles mating optimization [53,54,55,56,57,58,59,60,61,62,63,64,65,66].

However, generally, the correctness and accuracy of the results obtained by the interpretation techniques aforementioned rely on the precision by which the anomaly of the source structure(s) is extracted from the entire measured information. The merits of these algorithms have been verified in various ways even though the so-called derivative-based local optimizers require a priori information. Their success largely depends on the choice of the initial guess, and therefore, they cannot entirely explore the model space and traditional global optimizers may be captured by massive local minima existed within the model space and may not be exploitable enough to approximate the global optimum when dealing with some problems.

Notably, in the late 1990s, researchers proposed the joint inversion technique [67] to better address the above-mentioned non-unique and ill-posed mathematical nature of geophysical methods [68,69]. This approach efficiently reduces the ambiguity and non-uniqueness inherent in independent inversion procedures, leading to higher stability and accuracy in parameter estimation [70,71,72,73].

The effectiveness of Hunger Games Search (HGS) in interpreting the geomagnetic data caused by ore bodies was first validated in the study by Ai et al. [74]. Ai et al. [75] further verified the effectiveness of HGS for solving self-potential data inversion problems. We are thus motivated to further implement this novel global optimizer, HGS, for the delineation of gravity anomalies generated by different type of causative sources with application to geothermal energy identifications and volcanic activity studies. The superiority of HGS was verified in three theoretical cases of different circumstances. Practically, HGS was also experimented on two field cases from India and Japan. Prior to the inversion experiments, we implemented modal analyses to map the error topography of the objective function predetermined to understand the resolvability characteristics of model parameters and to investigate their effects. In applications, we successfully eliminated the regional background effect with the second moving average (SMA) technique. Some studies assessing post-inversion uncertainty were conducted to comprehend the consistency and credibility of the model solutions. Finally, the results of HGS applications were thoroughly, systematically, and unbiasedly compared with those of the well-established PSO algorithm regarding the trade-off of balancing the exploration and exploitation phases in metaheuristics. Previous geophysical, geological, and borehole studies were further gathered to generate coherent and convincing evaluations toward the two global optimizers involved.

2 Methods

2.1 Forward modeling

The general expression of gravity anomalies g(x _i ) for a spherical structure, an infinite horizontal cylinder, or a semi-infinite vertical cylinder, at any given point x _i along the principal profile on the free surface in a Cartesian coordinate system (Figure 1), is given by previous studies [18,76,77] as

(1) g ( x i ) = k z 0 [ ( x i − x 0 ) 2 + z 0 2 ] q , i = 1 , 2 , 3 … M ,

k = 4 3 π G Δ σ R 3 , 2 π G Δ σ R 2 , and π G Δ σ R 2 z 0 for q = 1.5 , 1 , and 0.5 , respectively .

where q stands for a shape factor (dimensionless), which determines the shape of the buried geological body; q = 1.5 (sphere), 1 (horizontal cylinder), and 0.5 (vertical cylinder); Δσ means the density contrast (g/cm³) between the causative source and the host rock; G denotes the universal gravitational constant (m³ kg⁻¹ s⁻²); x ₀ defines the horizontal coordinate of the geometrical center of the buried causative structure (m); k represents the amplitude factor (mGal m^2q−1); and R is the radius of the buried object (m).

Figure 1

Geometries and parameters of (a) a sphere, (b) a horizontally placed infinite cylinder, and (c) a vertically placed semi-infinite cylinder.

3 HGS optimizer

The HGS global optimization algorithm mathematically simulates the rules used by the animals during hunger, such as instinctive behaviors, preferences, and decisions. Yang et al. [79] clearly demonstrated the superiority of the HGS algorithm over the existing well-known global optimizers using 23 challenging benchmark functions. Ai et al. [74] further showed the effectiveness of HGS in the inversion of geomagnetic data. The mathematical background of the tailored HGS optimizer within this study can be briefly summarized as follows:

Communicating and preying characters of the predators are formulated by ref. [79]

where r 1 and r 2 are the numbers randomized between 0 and 1, rand(∼) produces a normally distributed random number, t is the current iteration, W 1 → and W 2 → are the weights of hunger that change for each parameter, X ( t ) → represents the location of the parameter set (R, Δσ, x ₀, z ₀, q) in the model space, X b → denotes the location of the best parameter set in the current iteration, and l is a predefined number between (0, 1), which should be tuned for the optimization, and

In equation (6), F ( i ) is the fitness value of each parameter, N represents the population dimension, BF is the best fitness in the current iteration, sech ( x ) = 2 / e x + e − x , R → = 2 × u × rand − u , and u = 2 ( 1 − ( t /Max_iter ) ) . Max_iter is the maximum iteration. The hunger weights in equation (5) are defined as follows [79]:

where r 3 , r 4 , and r 5 are the randomized numbers between 0 and 1, and hungry ( i ) is given by

where

and

where r 6 and r are also randomized numbers between 0 and 1; F ( i ) is the parameter fitness value; BF and WF are the best and worst fitness values, respectively; UL and LL are the upper and lower search limits. All the randomized numbers are regenerated in each iteration. Notably, LH (basically varies within [0, 1,000]) is also a predefined number that should be tuned before the application. In all global optimizers, algorithm-based control parameters dictate the performance in terms of accuracy, convergence speed, and stability. The best values of these parameters may vary depending on the mathematical nature of the handling problem. Therefore, before the data inversion studies performed via HGS, parameter tuning test of l and LH must be carried out [53,78,80].

4 Objective function and stopping criteria

Employing HGS for interpretation of gravity anomalies associated with geothermal investigations and volcanic activities involves treating the model parameter set (R, Δσ, x ₀, z ₀, q) as search agents. These search agents are guided by numerically imitated hunger activities to iteratively converge. We terminated the iteration process when the predefined maximum iteration is reached. The objective function used to be minimized is the root-mean-square-error (RMSE) as follows:

where g _obs is the observed data, g _cal denotes the data calculated from the model parameters estimated, and M is the data size.

5 Modal analysis

Pre- and post-inversion uncertainty appraisal analyses are vital for understanding the ambiguity and credibility of the estimated model solutions due to the well-acknowledged non-uniqueness and ill-posedness of geophysical inverse problems [58,81,82]. We tend to investigate the modal type of the determined objective function by means of mapping and categorizing the error topographies within each parameter pairs so that the resolvability and possible uncertainties between model parameters can be acknowledged. Without the loss of generality, the objective function within the modal analysis is constructed as follows:

In equation (13), g p = 3 ( x i ) is the simulated gravity anomaly generated by a multi-structure model consisting of a sphere, a horizontally placed infinite cylinder, and a vertically placed semi-infinite cylinder, where its true parameter set is presented in Table 1 (the profile length is 120 km, and the station interval is 3 km). We add 50% perturbations to the original parameter set (R, Δσ, x ₀, z ₀, q) to generate the model space (search range in Table 1). Figure 2 exhibits the distributions of log₁₀(RMSE). The topography maps are normalized to [0, 1] for better visualization. The middle locations of the presented maps (black symbols) locate the global minimum. The distribution maps of R–x ₀, x ₀–Δσ, z ₀–Δσ, and z ₀–x ₀ pairs clearly reveal that the global minima values are enclosed with almost rounded contours and local ones are barely situated, which exactly shows the uni-modal feature [83] and illustrates that these pairs can be assessed individually by implementing an efficient optimizer. Nevertheless, various smooth plates, elongated valleys, or basins with different flat bottoms representing the lowest fitness region displaying the multi-modal characteristics [84] are observed in other maps (e.g., R–Δσ, R–z ₀, R–q, q–x ₀, q–Δσ, and q–z ₀ in Figure 2), which indicate the corresponding parameter pairs hold different sensitivities, making the problem complicated, decreasing the resolvability characteristics, and increasing the uncertainties in the solutions retrieved. We can thus conclude that the objective function yields the composite modality [85] since the uni-modality and multi-modality co-exist regarding all the model parameter pairs. Error topographies of the constructed objective function show the characteristics of non-linearity, high dimensionality, and various shapes in paired model spaces, which means that estimations of the model parameters are challenging. However, these difficulties can be overcome if an inversion algorithm has sufficient ability to approach the global minimum as close as possible without compromising its robustness in order not to be captured by the local minima.

Figure 2

Modal analysis (parameter set is (R, Δσ, x ₀, z ₀, q)) of the defined objective function (50% fluctuations are added to the original parameter set to generate the model space; distribution of log₁₀(RMSE) in paired model parameter spaces (a–j) are illustrated; notably, these topographies are normalized to better demonstrating the variation characteristics).

6 Parameter-tuning studies for HGS and PSO

The effectiveness of global optimizers depends to a large extent on finely tuned control parameters because there is no algorithm-based control parameter value that suits all types of inverse problems [53,74,86]. Hence, a crucial step of effectively applying a metaheuristic is to find the best user-defined control parameter(s). As aforementioned, l and LH control the exploitation and exploration stages of HGS. To investigate the effect of tuning l and LH of HGS, a noise-free residual gravity anomaly derived from the previous multi-structure model was inverted (model parameters are still kept in Table 1). 50% oscillations of the true parameter set (R, Δσ, x ₀, z ₀, q) were adopted again, and the profile length (120 km) and station interval (3 km) were remained. HGS was independently performed 30 times with N = 80 and Max_iter = 140 in order to suppress the stochastic nature of metaheuristics so that more objective evaluation can be achieved. As l and LH of HGS vary successively, the mean values obtained are displayed in Figure 3a. Correspondingly, Figure 3b demonstrates the variation of the standard deviations of the calculated RMSE values, illustrating the inversion uncertainties. Curves in Figure 3 are normalized into [0,1] for better visualization. Different performances, which are controlled by various combinations of l and LH and very similar variation behaviors (Figure 3a and b), indicate that ill-determined l and LH can intensively reduce the performance of HGS in terms of affecting the solution correctness and the stableness of approximating the true parameter values. Moreover, l is a dominant factor in controlling the performance of HGS comparing to LH. The excellent performance of HGS is intensively corrupted with the amplification of l, however, mildly contaminated with the increment of LH. We finally determined the best combination of LH and l (LH = 10, l = 0.01).

Figure 3

Analyzing the effect of tuning l and LH of HGS: (a) variation behavior of the mean RMSE values obtained via the increment of l and LH, (b) inverted uncertainties (STD values) of HGS via the variation of l and LH. The mean(RMSE) values and the STD(RMSE) values are normalized into [0,1] to better visualize the variation features.

In order to make unbiased comparative studies, we performed tuning studies for the well-known PSO further. Inertia weight (w) and two social factors (c ₁ and c ₂) are the algorithm-based control parameters of PSO. We used the same N and Max_iter numbers as in the HGS-tuning studies. Nine parameter sets suggested in the literature were taken into consideration in the applications. The results of 30 independent runs of PSO are given in Table 2. The noise-less gravity response was generated by the introduced multi-source model. It is clearly seen that the best statistical results for the handling inversion problem were obtained with the use of the parameter set suggested by Carlisle and Dozier [87].

7 Synthetic data tests

The HGS method was applied to three synthetic tests with different properties to investigate its merits of inversion. The first test involves a residual gravity anomaly produced by the introduced multi-structure model. To generate the second test, we added a regional background (third-order polynomial) together with a local interference due to a near-surface spherical source to the residual anomaly simulated. In the last synthetic test, the second case anomaly was further contaminated with a certain degree of normally distributed random noise. The superiority of HGS using the parameters LH = 10 and l = 0.01 was testified by comparing it with the results of the PSO algorithm (w = 0.729, c ₁ = 2.041, and c ₂ = 0.948). Furthermore, we investigated the validity of implementing the SMA approach to suppress the regional information.

7.1 Synthetic test 1

As discussed earlier, more generally, measured gravity anomalies are induced by multiple sources. We thus tested and compared the performance of HGS and PSO using anomalies derived from the multi-structure model:

(14) g p = 3 ( x i ) π G = 4 3 × 0.5 × 12 , 000 3 × 24 , 000 [ ( x i + 50 , 000 ) 2 + 24 , 000 2 ] 3 / 2 + 4 × 6 , 000 2 × 30 , 000 [ ( x i − 40 , 000 ) 2 + 30 , 000 2 ] + − 0.5 × 9 , 000 2 [ ( x i + 1 , 000 ) 2 + 12 , 000 2 ] 1 / 2 .

The profile length and the station interval remained 120 and 3 km, respectively. Yellow symbols with a black line through the center in Figure 4 illustrate the simulated gravity anomaly, where deep gravity source (vertical cylinder) produces a broad negative anomaly curve (blue symbols with a black line through the center), while shallow causative structures (sphere and horizontal cylinder) generate narrow and sharp positive anomaly curves (black and green symbols with a black line through the center). For the inversion process, 30 independent runs were carried out with 80 search agents (N = 80) and 140 iterations for both HGS and PSO. Notably, we obtained the same number of model solutions as the number of independent runs; we thus followed the idea of averaging solutions within a specific range of the minimum [74,75] to generate the final output of both optimizers. We initially sorted the results from 30 independent runs in an ascending order according to the RMSE values and subsequently selected the first two solution sets to compute the mean output and associated uncertainties as the inversion results.

Figure 4

(a) Synthetic residual gravity anomaly (yellow symbols with a black line through the center) of the multi-source model superposed on a regional background (pink symbols with a black line through the center) and local interference (dark blue symbols with a black line through the center). Red symbols with a black line through the center in (a) represent the gravity anomaly corrupted by undesired regional and local information. Modeled subsurface structure is depicted in (b).

Table 3 shows the model space and the solutions obtained. The inversion responses of the HGS and PSO methods are compared with the residual gravity anomaly in Figure 5a. Figure 5b and c shows the behaviors of the mean values and standard deviations of calculated RMSE results against iterations. Figure 5d shows the yielded relative errors (REs) of the model parameters obtained via HGS and PSO. The RE between the true and predicted parameter sets is determined by

(15) RE = ∣ S true − S pre ∣ ∣ S true ∣ ,

where S _pre and S _true are the estimated parameter set and the true parameter set, respectively. From results (Figure 5 and Table 3), the model parameters are well recuperated by HGS, and estimated anomaly of HGS is in good agreement with the simulated one. However, PSO algorithm produces unrecovered model parameters and larger fitting errors. In addition, HGS converged less than 60 iterations (indicated by a blue arrow in Figure 5b); however, PSO took more than 100 iterations to reach the convergence plateau (indicated by a black arrow in Figure 5b). These findings clearly show that the PSO algorithm in this test failed to balance the exploration and exploitation stage properly to escape the capture of massive local minima and efficiently approximate the global minimum.

Table 3

Inverted results of the residual and contaminated gravity anomaly induced by multiple sources using HGS and PSO without and with the implementation of the SMA method

Without the implementation of the SMA method
Parameter set	Ground truth	Search spaces	PSO results	RE	HGS results	RE
Inversion of gravity anomaly induced by multiple structures (a sphere, a horizontally placed infinite cylinder, and a vertically placed semi-infinite cylinder)
R (m)	12,000/6,000/9,000	6,000–18,000/3,000–9,000/4,500–13,500	13029.9 ± 3524.1/6978.9 ± 1322.5/8553.3 ± 1003.7	0.1925 ± 0.2042/0.2233 ± 0.1169/0.0761 ± 0.0861	12,000 ± 0/6,000 ± 0/9,000 ± 0	0 ± 0/0 ± 0/0 ± 0
Δσ (g/cm3)	0.5/2/−0.5	0.25–0.75/1–3/−0.75 to 0.25	0.6139 ± 0.1702/2.2046 ± 0.5246/−0.5413 ± 0.1598	0.3378 ± 0.1506/0.2331 ± 0.0546/0.2151 ± 0.2076	0.5 ± 0/2 ± 0/−0.5 ± 0	0 ± 0/0 ± 0/0 ± 0
x 0 (m)	−50,000/40,000/−1,000	−75,000 to −25,000/20,000–60,000/−1,500 to −500	−52700.3 ± 1343.6/41517.8 ± 2116.6/−1102.5 ± 156.1	0.0540 ± 0.0268/0.0432 ± 0.0463/0.1288 ± 0.1235	−50,000 ± 0/40,000 ± 0/−1,000 ± 0	0 ± 0/0 ± 0/0 ± 0
z 0 (m)	24,000/30,000/12,000	12,000–36,000/15,000–45,000/6,000–18,000	18800.6 ± 3968.6/27947.5 ± 3630.4/7515.1 ± 643.3	0.2166 ± 0.1653/0.0942 ± 0.0913/0.3737 ± 0.0536	24,000 ± 0/30,000 ± 0/12,000 ± 0	0 ± 0/0 ± 0/0 ± 0
q	1.5/1/0.5	0.75–2.25/0.5–1.5/0.25–0.75	1.5478 ± 0.0321/1.0261 ± 0.0170/0.5284 ± 0.0304	0.0319 ± 0.0214/0.0261 ± 0.0170/0.0568 ± 0.0609	1.5 ± 0/1 ± 0/0.5 ± 0	0 ± 0/0 ± 0/0 ± 0
RMSE between the residual gravity anomaly and the inverted residual responses			3.86 ± 0.65		0 ± 0
Inversion of gravity anomaly induced by multiple structures further corrupted by an undesired regional and local effect
R (m)			12726.7 ± 4822.9/6707.2 ± 1832.3/6211.8 ± 354.06	0.2727 ± 0.2354/0.1963 ± 0.2371/0.3097 ± 0.0393	15411.5 ± 682.7/7716.6 ± 1110.3/5326.6 ± 267.4	0.2842 ± 0.0568/0.2861 ± 0.1850/0.4081 ± 0.0297
Δσ (g/cm3)			0.5609 ± 0.08/1.9632 ± 0.1422/−0.4636 ± 0.1304	0.1499 ± 0.1243/0.0574 ± 0.0247/0.1805 ± 0.1647	0.6063 ± 0.1219/2.4644 ± 0.3034/−0.5855 ± 0.0695	0.2560 ± 0.1705/0.2322 ± 0.1517/0.1710 ± 0.1391
x 0 (m)			−44736.2 ± 12572.02/56893.7 ± 2024.1/−636.5 ± 101.2	0.2031 ± 0.1339/0.4223 ± 0.0506/0.3634 ± 0.1012	−25602.2 ± 537.7/55050.9 ± 1013.6/−1424.08 ± 106.6	0.4879 ± 0.0107/0.3762 ± 0.0253/0.4240 ± 0.1066
z 0 (m)			16424.8 ± 3876.2/36068.2 ± 6216.004/13563.4 ± 1719.5	0.3156 ± 0.1615/0.2022 ± 0.2072/0.1348 ± 0.1367	14499.7 ± 2895.4/33703.05 ± 1278.7/16511.8 ± 1263.2	0.3958 ± 0.1206/0.1234 ± 0.0426/0.3759 ± 0.1052
q			2.0753 ± 0.1959/0.9698 ± 0.0292/0.4482 ± 0.0264	0.38354 ± 0.1306/0.0325 ± 0.025/0.1034 ± 0.0528	1.6315 ± 0.0243/0.9965 ± 0.0089/0.4347 ± 0.0072	0.0877 ± 0.0162/0.0057 ± 0.0068/0.1304 ± 0.0144
Mean RMSE between the gravity anomalies estimated and the inverted			14.66 ± 0.11		13.003 ± 0.14
RMSE between the residual gravity anomaly and the inverted residual responses			40.88 ± 0.37		41.79 ± 0.39

With the implementation of the SMA method (the calculated mean results using different s values are given)
Parameter set	Ground truth	Search spaces	PSO results	RE	HGS results	RE
Inversion of gravity anomaly induced by multiple structures further corrupted by an undesired regional and local effect
R (m)	12,000/6,000/9,000	6,000–18,000/3,000–9,000/4,500–13,500	12028.7 ± 1790.6/5675.1 ± 726.5/9106.3 ± 1453.6	0.2708 ± 0.1008/0.2549 ± 0.1214/0.2085 ± 0.0635	12476.9 ± 787.4/6008.8 ± 426.9/8900.1 ± 448.3	0.0736 ± 0.0676/0.0685 ± 0.0585/0.0589 ± 0.0855
Δσ (g/cm3)	0.5/2/−0.5	0.25–0.75/1–3/−0.75 to 0.25	0.4711 ± 0.0457/1.9340 ± 0.1667/−0.5257 ± 0.0996	0.1950 ± 0.0851/0.1660 ± 0.1093/0.2556 ± 0.0747	0.5064 ± 0.0411/1.9676 ± 0.1204/−0.4901 ± 0.0340	0.0805 ± 0.0741/0.0660 ± 0.0505/0.0561 ± 0.0609
x 0 (m)	−50,000/40,000/−1,000	−75,000 to −25,000/20,000–60,000/−1,500 to −500	−50850.1 ± 6652.4/42583.6 ± 4581.7/−1036.3 ± 134.8	0.1287 ± 0.0784/0.1841 ± 0.0404/0.1787 ± 0.0837	−49856.4 ± 287.6/40756.06 ± 1577.9/−1056.8 ± 43.6	0.0070 ± 0.0134/0.0259 ± 0.0361/0.0673 ± 0.0552
z 0 (m)	24,000/30,000/12,000	12,000–36,000/15,000–45,000/6,000–18,000	24307.7 ± 1536.3/30698.7 ± 3766.9/12462.4 ± 962.8	0.1908 ± 0.0890/0.1576 ± 0.0886/0.0880 ± 0.0421	24285.1 ± 2393.7/31395.7 ± 1350.7/11904.9 ± 237.5	0.0669 ± 0.0691/0.0791 ± 0.0651/0.0097 ± 0.0188
q	1.5/1/0.5	0.75–2.25/0.5–1.5/0.25–0.75	1.6770 ± 0.1290/1.1212 ± 0.0710/0.4936 ± 0.0290	0.1262 ± 0.0760/0.1374 ± 0.0455/0.0635 ± 0.0205	1.5103 ± 0.0271/1.0002 ± 0.0076/0.4968 ± 0.0091	0.0131 ± 0.0147/0.0199 ± 0.0203/0.0120 ± 0.020
Mean RMSE between the calculated SMA gravity anomalies and the inverted ones			0.49 ± 0.01		0.39 ± 0.01
RMSE between the residual gravity anomaly and the inverted residual responses			59.03 ± 22.28		26.62 ± 23.23
Inversion of gravity anomaly induced by multiple structures corrupted by an undesired regional effect and local interference. 40% of random noise is further introduced
R (m)			11609.4 ± 1288.1/5719.5 ± 615.9/695.1 ± 695.1	0.2920 ± 0.0491/0.3068 ± 0.0580/0.2959 ± 0.0413	12288.1 ± 443.03/6015.1 ± 598.3/8988.4 ± 538.6	0.1013 ± 0.0556/0.1295 ± 0.0689/0.0852 ± 0.0470
Δσ (g/cm3)			0.5472 ± 0.0680/2.0310 ± 0.2356/−0.5074 ± 0.4809	0.3163 ± 0.0596/0.2589 ± 0.03198/0.2875 ± 0.0233	0.5118 ± 0.0395/1.9544 ± 0.1069/−0.4880 ± 0.0113	0.1162 ± 0.0497/0.1431 ± 0.0550/0.0914 ± 0.0339
x 0 (m)			−50642.4 ± 4546.4/40135.5 ± 2454.5/−1153.6 ± 131.9	0.2757 ± 0.0347/0.2737 ± 0.0422/0.3290 ± 0.0351	−49360.4 ± 2822.8/41346.6 ± 1210.5/−1001.4 ± 52.01	0.0569 ± 0.0870/0.1053 ± 0.0482/0.1196 ± 0.0694
z 0 (m)			25059.8 ± 2237.9/31389.5 ± 2296.2/11534.3 ± 1818.8	0.2681 ± 0.0252/0.2698 ± 0.0562/0.1909 ± 0.1118	24614.9 ± 1879.8/31697.2 ± 1950.8/11392.5 ± 1564.4	0.0979 ± 0.0568/0.1420 ± 0.0658/0.0718 ± 0.1189
q			1.7980 ± 0.0738/1.2115 ± 0.0596/0.5214 ± 0.0545	0.2035 ± 0.0543/0.2169 ± 0.0527/0.1179 ± 0.0643	1.4947 ± 0.0186/1.0646 ± 0.0718/0.5049 ± 0.0196	0.0166 ± 0.0119/0.0742 ± 0.0663/0.0324 ± 0.0344
Mean RMSE between the calculated SMA gravity anomalies and the inverted ones			0.65 ± 0.02		0.49 ± 0.01
RMSE between the residual gravity anomaly and the inverted residual responses			64.64 ± 29.56		38.12 ± 21.21

Note: Parameters given like “12,000/6,000/9,000” represent the parameter of the sphere, the horizontal cylinder, and the vertical cylinder, respectively.

Figure 5

Inverted results of model 1, (a) fittings between the simulated gravity anomaly (pink symbols) and the inverted responses of HGS and PSO (blue and black symbols with solid lines through the center), (b) the variation behavior of the mean RMSE via iterations, (c) the variation behavior of the standard deviation of RMSE via iterations, and (d) calculated REs of the inverted model parameters using HGS and PSO for the three causative sources (error bars are given in (a) and (d)).

7.2 Synthetic test 2

To further understand the performance of HGS, we introduced the regional and local effects to the previous residual anomaly to generate the second test using the following expression:

(16) g gen 2 ( x i ) = g p = 3 ( x i ) + g re ( x i ) + g local ( x i ) ,

g re ( x i ) = 0.5 × 10 − 13 x i 3 + 10 − 10 x i 2 + 10 − 10 x i + 10 , g local ( x i ) = 4 3 × 0.5 × 3 , 000 3 × 9 , 000 [ ( x i − 50 , 000 ) 2 + 9 , 000 2 ] 3 / 2 .

where g gen 2 ( x i ) denotes the contaminated result, g _re(x _i ) is the undesired regional background (third-order polynomial), and g _local(x _i ) represents the local perturbation generated by a near-surface spherical interference.

Red symbols with a black line through the center in Figure 4 represent the corrupted gravity anomaly g gen 2 ( x i ) . HGS and PSO were implemented again with the same settings and without the involvement of the SMA technique. Model space consistent with the first experiment was used. After each iteration, the mean, the standard deviation of RMSE values, and the RE of estimated parameters were calculated to investigate the performance of approaching the global minimum (Figure 6). Likewise, Table 3 keeps the quantitative results of estimations. It is clear that the mean RMSE value increased significantly considering the test 2. Moreover, the deviations between the calculated and the true model parameters of three causative bodies also increased. We thus conclude that both HGS and PSO exhibited undesirable performance with the imposed corruption. Next, the SMA method is implemented to reduce the regional effects.

Figure 6

Inverted results of the second test without the implementation of the SMA method: (a) fittings between the synthetic residual gravity anomaly (pink symbols) and the inverted responses of HGS and PSO (blue and black symbols with solid lines through their center), (b) the variation behavior of the mean RMSE via iterations, (c) the variation behavior of the standard deviation of RMSE via iterations, and (d) calculated RE values of inverted model parameters using HGS and PSO for three causative sources (error bars are also given in (a) and (d)).

Figure 7 illustrates the calculated SMA gravity anomalies (red lines) for various s values (0.7, 1.3, 1.6, 1.9, 2.2, and 2.5 × dx, dx equals the station interval). Blue and black lines represent the inverted responses using HGS and PSO, respectively. They are generally in very good agreement with the red lines except the anomalies generated by the interfering spherical source situated at the right corner. Correspondingly, inverted gravity anomalies of HGS and PSO comparing with the synthetic residual one obtained RMSE of two scenarios (RMSE between the pure residual gravity anomaly and the inverted residual responses of HGS and PSO; RMSE between the calculated SMA gravity anomalies of the second test and the estimated results of HGS and PSO) changing with different s values, and the RE of retrieved parameters are shown in Figure 8. Table 3 gives the detailed inversion results further.

Figure 7

Calculated SMA gravity anomalies (red solid lines) and estimated SMA gravity responses using HGS and PSO (blue and black solid lines; error bars are shown) for different s values (0.7, 1.3, 1.6, 1.9, 2.2, and 2.5 × dx, dx = 3 km).

Figure 8

Obtained results of the second test with the implementation of the SMA method: (a) fittings between the synthetic residual gravity anomaly (pink symbols) and the responses obtained via HGS and PSO (blue and black symbols with solid lines through their center), (b) the variation behavior of the calculated RMSE between the uncontaminated residual gravity anomaly and the obtained residual responses of HGS and PSO for different s values, (c) the variation behaviors of the computed RMSE between the calculated SMA gravity anomalies and the estimated SMA gravity responses obtained via HGS and PSO for different s values, and (d) REs of the calculated mean estimation of the yielded model parameters using HGS and PSO for three causative sources from the SMA gravity anomalies for different s values (error bars are also given in (a)–(d)).

Considering all the information stored in Figures 7 and 8 and Table 3, it is clear that the SMA technique proved useful in eliminating the regional effect from the composite gravity anomaly, and therefore, both algorithms yielded relatively agreeable results. However, examining the results in detail, one can see that the mathematical nature of simulating the special hunger behaviors of animals yielded lower data misfits, more accurate model parameter values, and higher inversion stability than PSO regardless of the high complexity and ill-posedness of the objective function, the number of causative sources, and s value of the SMA routine.

7.3 Synthetic test 3

Finally, to further inspect the attainability of both algorithms, 40% of normally distributed random noise is added to test 2 using the following equation:

(17) g ˜ gen 2 ( x i ) = g gen 2 ( x i ) + 0 .4 × min ( ∣ g gen 2 ( x ) ∣ ) × [ rand 1 ( M ) − rand 2 ( M ) ] ,

where min() returns the minimum value of the input, and rand₁() and rand₂() return an array containing pseudorandom values between [0, 1] of a given size. Using similar processing procedures, the SMA results of the noise-corrupted gravity anomalies and estimated SMA gravity responses using HGS and PSO (blue and black solid lines) for different s values are exposed in Figure 9. For simplicity, we directly give the inversion results detained in Table 3. According to the results presented in Figure 9 and Table 3, they reveal that the added 40% random noise only produces minor corruption to the calculated SMA gravity anomalies (red lines in Figure 9). Moreover, the noise effect will be diminished gradually as the s value increases.

Figure 9

SMA results of noise-corrupted gravity anomalies (red solid lines) and estimated SMA gravity responses using HGS and PSO (blue and black solid lines; error bars are also shown) for different s values (0.7, 1.3, 1.6, 1.9, 2.2, and 2.5 × dx, dx = 3 km).

Therefore, the performance of HGS and PSO preserves in a good manner. Unsurprisingly, HGS still yields better performance than PSO in terms of fitting errors, inversion accuracy of model parameters, and stability of optimizing. Additionally, a novel denoising approach termed modified non-local means is tailored for eradicating the noise component when dealing with potential noisy data with high magnitudes [96].

7.4 Post-inversion uncertainty appraisal study

As mentioned in the modal analysis section, considered inverse problem here is unstable and error-prone because of the complexity of escaping the capture of massive local minima while effectively exploiting the global minimum for optimizing. Hence, post-inversion uncertainty assessment is essential for validating the credibility of the model solutions obtained. To carry out this analysis, calculated model parameters obtained from 30 independent runs of HGS and PSO were selected and sorted in an ascending order of their misfit values. Note that the solutions of the first synthetic test were used. We assumed An as a variable, which determines the number of sorted parameters that were used to compute the final results by calculating the mean responses. Figure 10 explicitly exhibits that HGS obtains strong coherency and low sensitivity with the increment of An. However, the performance of PSO is significantly contaminated with a larger An (there’s a 36 mGal gap between the obtained mean response of An = 2 and An = 27). Notably, this phenomenon correlates with the aforementioned high inversion uncertainties for inversion of PSO, which further verifies that PSO is incapable of supporting robust outcomes in this presented inverse problem. Moreover, Figure 10 explains why we sorted and partly selected the model solutions from 30 independent implementations in applications, i.e., to compare the performance of the involved two global methods in their best situation and yield more objective evaluations.

Figure 10

Post-inversion uncertainty appraisal analysis: yielded mean responses of HGS and PSO varying with An.

8 Application to geothermal investigation and volcanic activity studies with post-inversion uncertainty appraisal analyses

Two available field data from India and Japan for geothermal investigation and volcanic activity study were analyzed to further examine the applicability of the introduced methodology. Selected profile data were processed as residual gravity anomalies and delineated by Essa and Diab [18] using BA (previous geophysical estimations are presented in Table 4). Therefore, we also initially treated these digitized data sets as residual data and applied HGS algorithm (LH = 10, l = 0.01) to perform further inversions with the comparison of the PSO method (w = 0.729, c ₁ = 2.041, and c ₂ = 0.948). Both algorithms were independently run 30 times with 80 search agents and 140 iterations. Subsequently, we exclusively implemented HGS to the calculated SMA gravity anomalies to understand whether these anomalies are contaminated with unwanted information. Finally, the solutions obtained were interpreted in an integrated manner with known geological settings, geophysical estimations, and/or available borehole validations.

8.1 The Phenaimata gravity anomaly, Gujarat, India

As depicted in Figure 11a, the Phenaimata igneous complex locates northern the Narmada River and belongs to the Narmada–Tapti tectonic zone (NTTZ). The Narmada rift appears to constrain the positioning of this plug-like entity. This intrusive, located 23 km southwest of Chhota Udaipur on the left bank of the Heran River (Figure 11b), is a bimodal igneous complex, with volcanic and plutonic rock types representing tholeiitic and alkaline magmatism. Notably, Phenaimata plug is a differentiated igneous complex, which is comprised of basalt, layered gabbro, diorite, nepheline-syenite, lamprophyres, and granophyres [97,98]. The Phenaimata complex is comprised of 2/3 basalts and 1/3 alkaline plutonic series. According to the study by Hari et al. [98], the Phenaimata igneous complex yielded suitable environment to produce orthopyroxene gabbro. Based on the studies of the petrological modeling, magma accumulation in the crust, contemporaneous assimilation, and fractional crystallization are the main reasons for the origination of the gabbroic rocks. Besides, reverse magnetization has been found in olivine gabbro and syenite rocks, which implied emplacement toward the end of the major Deccan volcanism period in the Phenaimata igneous complex [99,100].

Figure 11

(a) Simplified geological and tectonic map of the Deccan traps (green area) and adjacent area showing locations of some important igneous complexes in the NW and central India. Gray color rectangles indicate the present study area covering Phenaimata igneous complex [101,102], and (b) geology and physiography of the Phenaimata igneous intrusive and its vicinities. B–B′ represents the profile considered for gravity elucidation using the introduced metaheuristic [98,101].

Figure 12 gives the Bouguer gravity anomaly map in the Phenaimata igneous exposure. It shows an elliptical-shaped anomaly closure to the north of this igneous complex, with the principal axis of the anomalies aligned roughly in the ENE–WSW direction, which aligns with the Narmada–Tapti tectonic zone’s orientation. A north–south B-B′ 30 km long profile was taken approximately perpendicular to the closure anomaly of Figure 12 using a 0.8 km sampling interval. The B–B′ residual gravity anomaly profile (Figure 13a) was obtained after discarding the regional information from the raw gravity anomaly by the study by Essa and Diab [18]. Based on the inversion result of the study by Essa and Diab [18], the intrusive body of the Phenaimata igneous complex anomaly was well approximated by a spherical source (Table 4). In our case, initially, we applied HGS and PSO to the anomaly and considered it as a residual response induced by a spherical source again. Wide parameter-searching ranges were assigned (Table 5). Besides, it is worth noting that inverted parameter sets using HGS and PSO here were still calculated in the close vicinity of the low misfit regions (An = 2) of the objective function in order to compute the final result and compare the performance of HGS and PSO in their best manner. The next field case was processed in the same manner. Figure 13a shows the comparison between the field anomaly (red symbols with a black line through the center) and the reconstructed mean responses of HGS and PSO (blue and black symbols with solid lines through their center). As Figure 13a depicts, the assumed spherical source approximated the buried source successfully and HGS yielded lower fitting errors (the calculated RMSE of HGS is 1.65 mGal, which is smaller than the 4.22 mGal of PSO). Retrieved parameter sets (R, Δσ, x ₀, z ₀, q) are given in Table 5 in detail (PSO obtains higher standard deviations, and it thus yields lower stability). Figure 13b illustrates the yielded mean responses of HGS and PSO varying with An. It can be clearly observed that the performance of PSO decreases with a larger An (there is an 8.34 mGal gap between the obtained mean response of An = 2 and An = 27). This phenomenon well correlates with the high standard deviations for the inversion of PSO. However, only minor fluctuations occur in HGS (smaller gap 3.18 mGal between the mean response of An = 2 and An = 27), which verifies that HGS is generally insensitive to larger An values, which further confirms the robustness of HGS and guarantees its credibility of retrieving model parameters. Finally, Table 5 also gives the calculated mean model parameters of HGS using the SMA method for s = 0.5, 0.75, and 1 × dx (dx = 0.8 km), which are in good agreement with the yielded results without the implementation of the SMA method and previous geophysical estimation in Table 4 [18]. It thus validates that the B–B′ profile gravity anomaly is barely corrupted by the regional field effect. Furthermore, according to the study by Singh et al. [101], the density contrast (Δσ) is 0.21 g/cm³ and the depth to the top of the subsurface structure is about 2.76 km, which correlates with the inversion results (Table 5) of HGS relatively good.

Figure 12

Obtained Bouguer gravity anomaly of Phenaimata igneous complex, Gujarat, India (symbol + represents the measurement stations). B–B′ shows the profile considered for gravity elucidation [101].

Figure 13

Inverted results of the profile B–B′, Phenaimata igneous complex, Gujarat, India: (a) fittings between the field gravity anomaly (red symbols with a black line through the center) and the inverted mean responses of HGS and PSO (blue and black symbols with solid lines through their center), and (b) yielded mean responses of HGS and PSO varying with An.

Table 5

Interpretation results using HGS and PSO without and with the implementation of the SMA method

Parameter set	Search spaces	PSO results	HGS results (without SMA)	HGS results (with SMA; the mean output of various s values)
Profile B–B′ of the Phenaimata gravity anomaly, Gujarat, India
R (m)	1,500–5,100	1501.4 ± 14.1	1804.2 ± 13.3	1967.2 ± 12.2
Δσ (g/cm3)	0.1258–0.4280	0.1378 ± 0.0148	0.2574 ± 0.0290	0.2095 ± 0.0020
x 0 (m)	500–1,700	997.6 ± 38.8	1059.03 ± 0.07	1066.9 ± 2.6
z 0 (m)	2,650–9,010	3629.7 ± 25.1	5399.9 ± 0.2	5392.7 ± 16.9
q	0.75–2.55	1.2716 ± 0.0062	1.3133 ± 0.0004	1.3149 ± 0.0004
RMSE (mGal) between the observed gravity anomaly and the inverted residual responses via PSO and HGS are 4.22 ± 0.01 and 1.65 ± 0.15, respectively.
RMSE (mGal) between the calculated SMA anomaly and the SMA result obtained via HGS is 0.11 ± 0 (the mean output of various s values)
Profile A–B of the Hohi volcanic gravity anomaly, Kyushu, Japan
R (m)	3,000–10,200	3000.2 ± 80.04	3060.1 ± 71.6	2993.1 ± 20.5
Δσ (g/cm3)	−2.8535 to −0.8392	−0.8393 ± 0.0001	−0.9645 ± 0.0959	−1.0312 ± 0.0084
x 0 (m)	−5,000 to 5,000	−228.8 ± 13.7	−264.5 ± 3.2	−264.3 ± 1.5
z 0 (m)	1,950–6,630	2677.09 ± 0.01	3799.9 ± 0.01	3760.2 ± 5.2
q	0.315–1.071	0.5219 ± 0.0001	0.5168 ± 0.0058	0.5163 ± 0.0008
RMSE (mGal) between the observed gravity anomaly and the inverted residual responses via PSO and HGS is 3.79 ± 0.06 and 1.13 ± 0.09, respectively
RMSE (mGal) between the calculated SMA anomaly and the SMA result obtained via HGS is 0.04 ± 0 (the mean output of various s values)

8.2 Hohi volcanic gravity anomaly, Kyushu, Japan

The Hohi volcanic zone lies in central Kyushu, southwest Japan. It is a 45 km wide, 70 km long volcano-tectonic depression with Plio-Pleistocene volcanic materials widely dispersed in an E–W orientation. These volcanic rocks contain limited amount of volcanic-clastic material and are mostly comprised of dacitic pyroclastic-flow deposits and andesitic lava flows. From around 5 Ma to the present, the Hohi volcanic zone plays an important role in disclosing that volcanic activity and subsidence alternated under a regional extensional stress field, with the depression mainly remunerated by filling of an equal amount of volcanic material [103,104]. Based on the evidence of drill-hole, geochronologic, and gravity data, the buried Shishimuta caldera was identified beneath post-caldera lava domes and lacustrine deposits in the center of the Hohi volcanic zone. The Yabakei pyroclastic flow, erupted 1.0 Ma ago with a bulk volume of 110 km³, is the source of the Shishimuta caldera in the center of the Hohi volcanic zone (Figure 14). The Shishimuta caldera is an 8 km wide, 3 km deep breccia-filled funnel depression with a V-shaped negative Bouguer gravity anomaly of up to 36 mGal (Figure 15). The andesitic breccia fills the caldera and has a relatively low density; it was likely produced by the fragmentation of disturbed ceiling rock during the intense Yabakei eruption and subsequent collapse [105]. In this work, a NE–SW Bouguer gravity anomaly profile A–B with a length of 13.5 km and data sampling interval of 0.4 km (Figure 16a) was taken to be further interpreted using HGS and PSO. This anomaly was also studied by Essa and Diab [18] using BA (Table 4). Based on the model parameters retrieved, Essa and Diab [18] suggested that the source of the Shishimuta caldera anomaly is a vertical cylinder-like model. Likewise, initially, the anomaly was elucidated as a residual response induced by a vertical cylinder again within this work. Parameter-searching ranges of PSO and HGS can still be observed in Table 5. Calculated mean responses of HGS and PSO (blue and black symbols with solid lines through their center) are compared with the acquired one (red symbols with a black line through the center) in Figure 16a. Coherently, as Figure 16a depicts, the assumed vertical cylinder model approximated the buried source well and HGS yielded lower fitting errors again (calculated mean RMSE of HGS is 1.13 mGal, which is smaller than the 3.79 mGal of PSO). Detailed inversion results can be found in Table 5 (PSO still yields higher standard deviations regarding all the model parameters). Moreover, according to Figure 16b, much higher robustness for the inversion process of HGS can still be recognized. Finally, based on the delineated result from the SMA anomalies for s = 0.5, 0.75, and 1 × dx (dx = 0.4 km) in Table 5, it shows that the profile of Hohi volcanic gravity anomaly is also barely contaminated with the regional field effect. Furthermore, drill holes outside Shishimuta caldera (boreholes with dotted circles; a–k in Figure 14) encounter Pre-Tertiary rocks at depths of 2.0–2.5 km [106,107]. Drill holes up to 3 km deep within the caldera (bore-holes with opened circles; y and z in Figure 14), on the other hand, penetrate Plio-Pleistocene volcanic rocks instead of Pre-Tertiary materials [108]. The introduced HGS interpreted the depth to Pre-Tertiary basement rocks of the caldera about 3.8 km (deeper than the result 2.677 km of PSO), which agrees very well with the drilling (Figure 17) and depth modeling information to gravity basement of the caldera ([103,104] and [109]). Previous studies ([103,104] and [109]) introduced a gravity model that defined the V-shaped depression with depth to the Pre-Tertiary basement up to 3.8 km below sea level.

Figure 14

Detailed geologic map of the Kusu Basin and Shishimuta area of the Hohi volcanic zone, Kyushu, Japan [18,109]. Dotted circles represent the locations of drill holes, and opened circles refer to the location of 3,000∼m deep drill holes with spot coring.

Figure 15

Bouguer gravity map of the central part of Hohi volcanic zone, Kyushu, Japan (contour interval is 2 mGal). The dotted line correlates with the location of the Shishimuta caldera (taken from the study by Essa and Diab [18], and they modified after Komazawa and Kamata [109]). Lines A–B show the target profile considered for interpretation using involved algorithms.

Figure 16

Inverted results of the profile A–B, Hohi volcanic zone, Kyushu, Japan: (a) fittings between the observed anomaly (red symbols with a black line) and the inverted mean responses of HGS and PSO (blue and black symbols with solid lines) and (b) yielded mean responses of HGS and PSO varying with An.

Figure 17

Previous interpretive result of N–S profile of the Shishimuta caldera [18,103].

Another thing that needs to be clarified is that, as we can tell the differences from the theoretical tests and the field applications, PSO exhibited more erroneous performance in terms of larger fitting errors, higher inversion uncertainties, and lower robustness regarding the gravity responses that were generated by a multi-structure model. This is simply because the model parameters of multiple structures will enhance the non-uniqueness, non-linearity, and ill-posedness of the inversion process. However, HGS preserves its superiority of retrieving model parameters regardless of the complexity of the subsurface geological–geophysical settings.

9 Conclusions

Retrieving model parameters describing various subsurface structures from a precarious process due to the well-acknowledged non-uniqueness and ill-posedness of gravity inversion problems is the intrinsic aim in gravity data elucidation. We introduced a novel derivative-free nature-inspired global optimizer termed HGS and validated its superiority in determining model-describing parameters from gravity anomalies in detail. Prior to the implementation of HGS, modal analyses using error topography maps of model parameter pairs showed the non-linear and high-dimensional mathematical nature of the objective function tailored. To take precautions against these limitations, parameter-tuning studies were performed. With the tuned control parameters, HGS was applied to synthetically generated gravity anomalies and subsequently experimented on two field data examples from India and Japan for geothermal energy exploration and volcanic activity study. In addition, the effectiveness of implementing the SMA scheme to suppress the regional background was examined and certified with regard to the contaminated cases. PSO, the most commonly used global optimizer in geophysical inversions, was also unbiasedly applied to perform comparative experiments. Post-inversion uncertainty appraisal analyses were further implemented to understand the reliability of the model parameters obtained.

In general, simulated and field anomalies were more satisfactorily inverted with the HGS algorithm. The inverted outcomes of the HGS algorithm were well in line with the available geologic, geophysical, and/or drilling information. We conclude that the HGS optimizer is more efficient than PSO in balancing global exploration and the local exploitation processes in this type of geophysical inverse problems. This novel metaheuristic can, therefore, be considered as a promising alternative in geothermal reservoir identifications and volcanic activity studies.