Gravity Data Imaging Using Local Wavenumber-Based Algorithm: Sustainable Development Cases Studies

A fast effective inversion algorithm is proposed herein to interpret gravitational responses caused by mineralized/ore sources (sphere, vertical and horizontal cylinders). The algorithm relies on local wavenumber and correlation imaging techniques. The correlation factor (R) between the local wavenumber of observed gravitational field and that of computed field was calculated, and the maximum Rmax was considered to correspond to the best true model (parameters). The proposed algorithm was applied to two theoretical examples, including an example contaminated with regional background and another multisource example. Besides, the proposed approach was used on three different real field cases for mining/ore investigation from Canada and Cuba. From the results obtained from the theoretical and real examples and by comparing the results with drilling and literature information, it was concluded that the method is effective, is applicable even for more than one source, is accurate, and does not necessitate any prior knowledge of the source shape.

Several techniques have been created to analyze gravity data from different structures by approximating these different source structures to simple geometric shapes (e.g., spheres and cylinders) (Mohan et al., 1986;Salem et al., 2004;Essa 2013;Singh & Biswas, 2016) and then estimating the various parameters of the sources (e.g., amplitude factor, depth, horizontal position). Among these techniques are graphical ones that rely on characteristic points, nomograms, and master curves (Siegel et al., 1957;Grant & West, 1965;Nettleton, 1976;Prakasa Rao et al., 1986;Reynolds, 1997;Kara & Kanli, 2005;Essa, 2007Essa, , 2012; Euler and Wernerdeconvolution (Werner, 1953;Thompson, 1982;Stavrev, 1997;Zhang et al., 2000;Ghosh, 2016), derivatives method (Saad, 2006;Essa, 2007;Ekinci et al., 2013), least squares approach (Gupta, 1983;Abdelrahman & Sharafeldin, 1995;Essa, 2011Essa, , 2014Abdelrahman & Essa, 2015) and 2D, 2.5D, and 3D modeling (Chai & Hinze, 1988;Pinto & Casas, 1996;Zhang et al., 2001;Eshaghzadeh & Hajian, 2018). Most of these previous methods have difficulties in inversion, which include the necessity of prior information (geological information), dependence of few specific points and accuracy in regional-residual anomaly separation, and ignoring other points along a profile, resulting in ill-posed and ambiguity issues (Zhdanov, 2002;Tarantola, 2005;Essa & Munschy, 2019). These different problems lead to high uncertainty in predicted source parameters and estimated anomaly or model shift (Nettleton, 1976; Table 1. Definitions of A, m, and q for source geometry shown in Figure 1; r represents the radius (m) of source body, q represents density contrast (gm/cm 3 ), and k represents the universal gravitational constant (6.67 9 10 11 m 3 /kg s 2 ) Source geometry A q m Vertical cylinder pkqr 2 0.5 0 Horizontal cylinder 2pkqr 2 1 1 Sphere 4 3 pkqr 3 1.5 1 these techniques include genetic algorithm (Amjadi & Naji, 2013;Di Maio et al., 2016Kaftan, 2017), simulated annealing (Sen & Stoffa, 2013;Biswas et al., 2014Biswas et al., , 2017Biswas, 2016), particleswarm-optimization (Toushmalani, 2013;Singh & Biswas, 2016;Essa & Munschy, 2019), and the differential evolution technique (Wu et al., 2014;Ekinci et al., 2016;Balkaya et al., 2017). An efficient imaging algorithm was developed in this study to fully interpret gravitational data caused by various underground structures (sphere, horizontal cylinder, and vertical cylinder). This approach relies upon calculating the correlation factor (R) parameter between the local wavenumber of observed gravitational field and that of the computed field. The model that corresponds to the highest R parameter was interpreted as the best model. This approach can be used in different applications, such as mineral and ore exploration, as it determines the different structure parameters, which include amplitude factor (A), depth (h), body origin (w), and shape factor (q), and it does not necessitate any prior knowledge of the source shape. Moreover, the method can be applied to estimate multisource parameters. To confirm the efficacy and the applicability of the suggested approach, this method was used to invert gravity data of three various theoretical examples with and without 20% Gaussian noise levels and three field examples from Canada and Cuba.

METHODOLOGY
The gravitational anomaly (g) caused by different geometrical structures (vertical cylinder, sphere, and horizontal cylinder) at horizontal position (x l , z) along a profile ( Fig. 1) Asfahani & Tlas, 2008;Essa, 2014;Essa et al., 2020) can be calculated as: where x o is the location of the structure in horizontal direction, z o is the depth of the structure, q is the shape factor, m is the exponent constant varied according to the value of q (i.e., for q = 0.5, m = 0; for q = 1 and 1.5, m = 1, as presented in Table 1), and A is the amplitude factor that relies upon the type of structure (Table 1). The observed local wavenumber (K Obs ) can be formulated as (Ma et al., 2017): where Substituting Eq. 3 into Eq. 2, K Obs can be calculated as: @g @z @g @x 2 4 3 5 ¼ 1 AS 2 @ 2 g @x@z : @g @x ! À @ 2 g @x 2 : @g @z where AS represents the amplitude of the analytic signal (Nabighian, 1972): By taking the horizontal ð @g @x Þ and vertical ð @g @z Þ derivatives of Eq. 1 analytically, substituting them into Eq. 4, and abbreviating the resultant equation, the calculated local wavenumber (K Cal ) is: For K Obs and K Cal , the correlation factor (s) can be obtained as (Ma et al., 2017): With Eq. 7, the R between K Obs and K Cal is calculated, and the maximum R corresponds to real body parameters (Ma et al., 2017). Figure 2 shows a flowchart of the sequence procedures of the proposed algorithm. Once the best optimal parameters are chosen from the search space based on the maximum R, the 2D image of the R of the preferred source (i.e., shape factor q and its related m) can be constructed as a function of the depth (m) of the subsurface. The solid black dot on the imaging section represents the correct location for depth and location.

THEORETICAL EXAMPLES
This section applies the proposed technique to three different noisy and noise-free synthetic examples to examine the efficiency and applicability of the proposed approach in interpreting gravity anomalies.

Model 1
A pure gravity profile caused by horizontal cylinder was computed using the following parameters: A = 150 mGal.m, z o = 4 m, q = 1, x o = 51 m, and profile length of 100 m (Fig. 3a). Interpretation began with the calculation of horizontal and vertical derivatives of the observed data ( Fig. 3b), and then, K Obs was calculated using Eq. 4 (Fig. 3c). Next, R was calculated by applying Eq. 7 (Fig. 3d) using different values of q (Table 2). Figure 3d shows the maximum R = 1 (black circle) located at x o = 51 m, z o = 4 m, and q = 1 (Table 3), indicating that the suggested method is highly efficient. The proposed method was used to estimate the inverted parameters (Table 3), and the errors of the different parameters account to 0%.  Table 3. True and recovered model parameters for Model 1: noise-free example presented in Figure 2 Model parameters True Recovered To investigate the stability and the performance of the proposed approach to noisy data, a 20% Gaussian random noise was added to the previous gravity anomaly (Fig. 4a). First, the vertical and horizontal derivatives of the noisy data were calculated (Fig. 4b), and then, K Obs was calculated applying Eq. 4 (Fig. 4c). Equation 7 was then applied to calculate R (Fig. 4d). Figure 4d shows the maximum R = 0.9061 (black circle) located at x = 50 m, z = 4.7 m, and q = 1 (Table 4), indicating that the suggested technique can be applied with high performance to noisy data. The proposed approach was used to estimate the inverted parameters (Table 4), and the errors of A, z o , x o , and q were 1.69, 17.5, 1.96, and 0%, respectively.

Model 2
In some real-world cases, multiple structures or bodies can cause gravity anomalies. Therefore, this study aimed to explore the applicability of the suggested approach in the presence of multisource models. A composite gravity anomaly was computed using a horizontal cylinder model (with the following parameters: A = 120 mGal.m, z o = 3 m, x o = 30  Table 4. True and recovered model parameters for Model 1: the noisy synthetic example presented in Figure 3 Model parameters True Recovered  m, and q = 1) and a sphere model (with the following parameters: A = 550 mGal.m 2 , z o = 5 m, x o = 80 m, and q = 1.5) (Fig. 5a). Interpretation began by calculating the vertical and horizontal derivatives of the observed composite profile (Fig. 5b). The K Obs was calculated using Eq. 4 (Fig. 5c), and the R was calculated using Eq. 7 (Fig. 5d). Figure 5d shows the maximum R of the first body to be 0.74 (first black circle) located at x o = 30 m, z o = 3.4 m, and q = 1 and the maximum R of the second body to be 0.70 (second black circle) located at x o = 80 m, z o = 5.2 m, and q = 1.5 (Table 5). The errors of A, z o , x o , and q were 2.06, 13.33, 0, and 0%, respectively, for the horizontal cylinder structure, while the errors of A, z o , x o , and q were 2.66, 4, 0, and 0%, respectively, for the spherical structure.
To use the suggested approach in the presence of noisy data, a 20% Gaussian random noise was added to the previous composite gravity anomaly (Fig. 6a). First, the vertical and horizontal derivatives of the noisy composite profile were calculated (Fig. 6b), the K Obs was computed (Fig. 6c), and R was calculated (Fig. 6d). Figure 6d shows the maximum R of the first body to be 0.58 (first black circle) located at x o = 30 m, z o = 3.7 m, and q = 1 and the maximum R of the second body to be 0.56 (second black circle) located at x o = 80 m and z o = 5.6 m, and q = 1.5 (Table 6). The errors of A, z o , x o , and q were 5.21%, 23.33%, 0%, and 0%, respectively, for the horizontal cylinder structure, while the errors of A, z o , x o , and q were 7.58, 12, 0, and 0%, respectively, for the spherical structure.
Thus, it can be concluded that the proposed method is highly applicable in multisource cases.

Model 3
The effectiveness of the proposed technique in the presence of regional anomaly was studied. Regional anomaly (first order) was computed, and a composite gravity anomaly consisting of vertical cylinder was computed using the following parameters: A = 190 mGal.m, z o = 5 m, x o = 51 m, and q = 0.5. Moreover, a 10% Gaussian random noise was added (Fig. 7a). Using the same procedure mentioned above, the vertical and horizontal derivatives of the observed composite profile were computed (Fig. 7b). The K Obs (Fig. 7c) and the R were then calculated (Fig. 7d). Figure 7d shows the maximum R to be 0.40 (black circle) located at x o = 56 m, z o = 5.9 m, and q = 0.5 (Table 7). The errors of A, z o , x o , and q were 9.56, 18, 9.80, and 0%, respectively. The results so far show the effectiveness of the proposed technique even if the data were contaminated with regional background.

FIELD EXAMPLES
To test the performance of the proposed approach, it was applied to three different real cases, one from Cuba and two from Canada.

Mobrun Anomaly, Qué bec, Canada
The Mobrun area is situated about 34 km northeast of Rouyn-Noranda (Barrett et al., 1992) (Fig. 8). It is underlain by the Renault Formation, which is lined with Destor rocks from north to south. The Renault Formation consists of large fragmental rhyolitic and andesitic units (Barrett et al., 1992). The footwall of the Mobrun area is composed of Copper Hill rhyolite unit, which is underlain by an andesitic-rhyolite and andesite sequence followed by a level of felsic pyroclastic rocks, while the hanging wall is composed of massive rhyolite that similarly correlates to the footwall stratigraphy of the main complex along the strike to the west. Brecciated rhyolite is present immediately beneath the main complex (Caumartin & Caill, 1990;Barrett et al., 1992). The Mobrun deposit is made up of two enormous sulfide lens complexes. The main complex is made up of five ore bodies. Figure 9 shows the residual gravity map over the Mobrun ore. A profile was extracted from the residual gravity map, and this profile was taken perpendicular to the massive sulfide body (Grant & West, 1965;Sivakumar Sinha, & Ram Babu, 1985;Roy et al., 2000;Biswas, 2015) (Fig. 10a). The gravity profile was 230 m long, and it was sampled at 2 m intervals. The gravity profile was subjected to the local wavenumber approach by computing the Figure 8. Geological map of the Noranda property, Canada (after Barrett et al., 1992). Open squares represent currently producing mines; closed ones represent past mines. horizontal and vertical derivatives (Fig. 10b), followed by computing the K Obs (Fig. 10c). The R was then calculated (Fig. 10d) using different values of q (Table 8). Figure 10d shows the maximum R to be 0.99 (black circle) located at x o = 2 m, z o = 46 m, and q = 1 (Table 9). Table 9 summarizes the predicted structure parameters, indicating that the body is a horizontal cylinder. Figure 10a shows a comparison between the observed and the inverted gravity profiles, which have an excellent coincidence. Figure 10a also shows that the present approach had less RMS error between the observed and calculated anomaly than was obtained by other techniques (Essa et al., 2020). Table 9 presents the comparison between the inverted parameters of the proposed method and those of the various methods in the literature.

Camaguey Chromite Anomaly, Cuba
The Camaguey chromite deposits formed at the intersection of serpentinized dunite and peridotite rocks and feldspathic rocks (Davis et al., 1957) (Fig. 11). Davis et al. (1957) collected chromite gravity data throughout a dipping sheet-shaped chromite-bearing ore deposit in the Camaguey District of Cuba as part of a United States Geological Survey exploration expedition. Figure 12 illustrates the gravity map over the Camaguey chromite deposits. A profile was extracted from the residual gravity map and was taken perpendicular to the strike of the chromite deposit (Fig. 13a). The gravity profile length is 89 m and was sampled at 2 m intervals. The local wavenumber approach was applied to the gravity profile by Figure 9. The Mobrun anomaly, Canada. Residual gravity anomaly map (Ekinci et al., 2016). AB denotes the gravity anomaly profile (Fig. 10) subjected to interpretation.
computing the horizontal and vertical derivatives (Fig. 13b), followed by computing the K Obs (Fig. 13c). After that, the R was calculated (Fig. 13d) using different values of q (Table 10). Figure 13d shows the maximum R to be 0.97 (black circle) located at x o = À 2 m, z o = 15 m, and q = 1 (Table 11). Table 11 presents the predicted parameters, indicating that the body is a horizontal cylinder. Figure 13a shows the comparison between the observed and the inverted gravity profiles, which are closely matched. In addition, Figure 13a also shows that the present approach had less RMS error between the observed and calculated anomaly than was obtained by other techniques (Essa et al., 2020). The depth obtained from the proposed method matched well with the drilling information (Fig. 14). Table 10 presents the comparison between the inverted parameters of the proposed method and those from different techniques in the literature.

Faro Gravity Anomaly, Yukon, Canada
Faro Mine lies 15 km north of Faro in the south-central Yukon Territory of Northern Canada.   It is located in the northwest area of the Faro Complex (Tang, 2011) (Fig. 15). The geography of the property is dictated by the Yukon Plateau and the adjacent mountains of Anvil Range, with altitudes above 1800 m. The regional geology of the Faro Mine can be illustrated by the Anvil district geology (Tang, 2011) (Fig. 16). Figure 17a shows the residual gravity profile above the Faro mine. Figure 17b displays the geological cross section of the profile. The profile length was 805 m and was sampled at 12 m intervals (Fig. 18a). The local wavenumber approach was applied to the gravity profile by computing the horizontal and vertical derivatives (Fig. 18b), followed by computing the K Obs (Fig. 18c). After that, the R was calculated (Fig. 18d) using different val-ues of q (Table 12). Figure 18d shows the maximum R to be 0.99 (black circle) located at x o = 378 m, z o = 290 m, and q = 1.5 (Table 13). Table 13 presents the inverted parameters, indicating that the body is a sphere. Figure 18a shows the comparison between the observed and the inverted gravity profiles, with relatively good matching. In addition, Figure 18a also shows that the present approach had less RMS error between the observed and calculated anomaly than was obtained by other techniques (Essa et al., 2020). Table 13 presents the comparison between the inverted parameters of the proposed method and those from different techniques in the literature. Figure 12. The Camaguey anomaly, Cuba. Residual gravity anomaly map (Davis et al., 1957;Roy, 2001).
AB denotes the gravity anomaly profile (Fig. 13) subjected to interpretation.

CONCLUSIONS
In this paper, an efficient inversion imaging algorithm was applied to model gravity data caused by different sources (sphere, vertical cylinder, and horizontal cylinder). The algorithm demonstrated here can be used in mineral/ore exploration because it can predict different structure parameters, including amplitude factor (A), depth (z o ), body origin (x o ), and shape factor (q), with high accuracy and with no need of a priori information. The suggested algorithm applies the correlation factor (R) between the local wavenumber of the observed gravitational field and that of the computed field. The results show that the maximum R reflects the best-estimated model. In addition, the proposed approach represents an imaging algorithm which offers good and fast (only taking several seconds of) imaging for the subsurface depth and location of hidden anomalous sources. The efficiency, accuracy, and stability of the proposed algorithm were tested by applying it to three different synthetic cases,   Figure 14. The Camaguey anomaly, Cuba. Drilling information; F, H-G, B-A, and J-C represent labels of drill holes (redrawn from Figure 4 of Davis et al., 1957).  including a pure example, a noisy example contaminated with different percentages (10% and 20%) of Gaussian random noise (i.e., normal distribution with mean = 0 and standard deviation = 1), and an example with regional background effects. The applicability of the proposed algorithm was checked on three real mineral/ore exploration examples from Canada and Cuba. The resultant models for the real cases correlated very well with drilling data and with different results published in the literature. Finally, from the present study, the proposed algorithm is suitable for mineral/ore deposits exploration, and it is recommended for application to geothermal investigation as well.