A Fast Interpretation Method for Inverse Modeling of Residual Gravity Anomalies Caused by Simple Geometry

An inversion technique using a fast method is developed to estimate, successively, the depth, the shape factor, and the amplitude coefficient of a buried structure using residual gravity anomalies. By defining the anomaly value at the origin and the anomaly value at different points on the profile, the problem of depth estimation is transformed into a problem of solving a nonlinear equation of the form 𝑓(𝑧)=0. Knowing the depth, the shape factor can be estimated and finally the amplitude coefficient can be estimated. This technique is applicable for a class of geometrically simple anomalous bodies, including the semiinfinite vertical cylinder, the infinitely long horizontal cylinder, and the sphere. The efficiency of this technique is demonstrated with gravity anomaly due to a theoretical model, in each case with and without random errors. Finally, the applicability is illustrated using the residual gravity anomaly of Mobrun ore body, situated near Noranda, QC, Canada. The interpreted depth and the other model parameters are in good agreement with the known actual values.


Introduction
Inversion of gravity data is nonunique in the sense that the observed gravity anomalies in the plane of observation can be explained by a variety of density distributions.One way to solve this ambiguity is to assign a suitable geometry to the anomalous body with a known density followed by inversion of gravity anomalies [1].Although simple models may not be geologically realistic, they are usually are sufficient to analyze sources of many isolated anomalies [2].The interpretation of such an anomaly aims essentially to estimate the body parameters such as shape, depth, and radius.Several graphical and numerical methods have been developed for interpreting gravity anomalies caused by simple bodies.The simplest formula used to approximate depth to a causative body from the residual gravity data is the half-g max rule [3,4].However, the drawback with this approach is that it is highly subjective and can lead to large errors [5].Gupta [6] presented a numerical approach to determine the depth to cylindrical and spherical models from the residual gravity data.Abdelrahman [7] argued that inserting the maximum gravity value g max as a known parameter in Gupta's formulation may lead to large error in the calculation of depth in the existence of noise.Recently, several computer-based methods of inverting gravity data to determine model parameters have been presented with various levels of success [8][9][10].A simple method proposed by Essa [11] is used to determine the depth and shape factor of simple shapes from residual gravity anomalies along the profile.Another automatic method, the least-squares method, was proposed by Asfahani and Tlas [12], by which the depth and amplitude coefficient can be determined.The principal difficulty with the inversion methods is the inherent nonuniqueness of the solution [13].Therefore, there is still a need for an interpretation technique that is robust, rapid, and can provide parameters of the bodies in field situations.
In this paper, an inversion technique based on nonlinear equation z = f (z) to analyze gravity anomalies due to simple structures.The inversion technique simultaneously estimates the depth, the nature of the source (shape factor), and the amplitude coefficient of the buried structures.The accuracy of the result obtained by this procedure depends upon the accuracy to which the residual anomaly can be separated from the Bouguer anomaly.In most cases, graphical methods [5] or standard numerical methods [14][15][16] can be used to separate the residual gravity anomaly attributable to the buried structure from the Bouguer data.Also, the accuracy of the result of the present method depends on the extent to which the source body conforms to one of the assumed geometries.The methodology is illustrated with theoretical models, in each case with and without random errors, and tested by the gravity anomaly of the Mobrun ore body, situated in the mining district of Noranda, QC, Canada.

The Method
The general vertical component of the gravity anomaly expression produced by a sphere (3D), an infinite long horizontal cylinder (2D), and a semiinfinite vertical cylinder (3D) is shown in Figure 1 and given in Abdelrahman et al. [17] as where for a sphere 1 for a horizontal cylinder 1 2 for a vertical cylinder R z. ( In (1), z is the depth, q is the shape factor, for example, the shape factors for the semiinfinite vertical cylinder (3D) (the gravity response in case of the semiinfinite vertical cylinder is only applicable when the radius of the cylinder is much smaller than the distance from observation position to the top of the cylinder.This is called the "thin vertical rod approximation."For the general case of the right vertical cylinder, the gravity response is much more complicated), horizontal cylinder (2D), and sphere (3D) are 0.5, 1.0, and 1.5, respectively.Also, the shape factor for the finite vertical cylinder is approximately 1 [18].The shape factor (q) approaches zero as the structure becomes a nearly horizontal bed, and approaches 1.5 as the structure becomes a perfect sphere (point mass).x i is the position coordinate, σ is the density contrast, G is the universal gravitational constant, and R is the radius.At the origin (x i = 0), (1) gives the following relationship: (3) Using (1), we obtain the following normalized equation at x i = ± N and x i = ± M where N = 1, 2, 3, . . .and M = 1, 2, 3, . . .
Let F = (g(N)/g(0)) and T = (g(M)/g(0)) then from (4), we get: z = e ([(ln F/ ln T) * (ln(z 2 /(M 2 +z 2 )))]+ln(N 2 +z 2 ))/2 , M / = N. (5) Equation ( 5) can be solved for z using the standard methods for solving nonlinear equations (e.g., [19]), and its iteration form can be expressed as: where z j is the initial depth, and z f is the revised depth; z f will be used as the z j for the next iteration.The iteration stops when Once, the depth (z) is known, the shape factor (q) can be estimated from the following form: where z c is the estimated depth.Finally, knowing the shape factor (q), the amplitude coefficient (A) can be estimated from the following form: where q c is the estimated shape factor.For each N and M value, we compute the values of the model parameters (z, q, and A) from ( 5), (7), and (8), respectively.Theoretically, the anomaly values at the origin and any two N and M distances are just enough to determine the model parameters.However, in practice, it is recommended to use all possible combinations of N and M values to determine the most appropriate source parameters solutions from all gravity data.We then measure the goodness of fit between the observed and computed gravity data for each set of solutions.The simplest way to compare two gravity profiles is to compute the standard error (μ) between the observed values and the values computed from estimated values of z, q, and A. The model parameters which give the least root mean sum squares differences are the best.In this way, we can select the best-fit source parameters solutions from all gravity data.

Synthetic Examples
We computed three different residual gravity anomalies, each consisting of the effect of local structure (semiinfinite Horizontal distance (km) Residual gravity anomaly (mGal) Residual gravity anomaly  vertical cylinder, horizontal cylinder, and sphere).The model equations representing the models are The three gravity anomalies are shown in Figure 1.Equations ( 5), (7), and (8) were applied to the residual anomaly profiles, yielding the model parameters: the depth, the shape factor, and the amplitude coefficient, respectively, solutions for all possible N and M points.The computed model parameters for the three models (a semiinfinite vertical cylinder model, a horizontal cylinder model and a sphere model) are summarized in Tables 1-3, respectively.In order to examine the influence of the noise on this approach, a 10% random noise has been added to the synthetic data using the following expression: where g rand (x i ) is the contaminated anomaly value at x i , and RND (i) is a pseudorandom number whose range is (0,1).The interval of the pseudorandom number is an open interval, that is, it does not include the extremes 0 and 1.
Table 1 shows the model parameters (z, q, and A) in the case of using a semiinfinite vertical cylinder model results are the same when using synthetic data.The depth is within 1%, the shape factor is within 2% and the amplitude coefficient is within 5.9%.Table 2 shows the model parameters in the case of using a horizontal cylinder model results are the same when using synthetic data.After adding 10% random errors in the synthetic data, the depth is within 4.2%, the shape factor is within 7%, and the amplitude coefficient is within 13.2%.Table 3 shows the model parameters in the case of using a sphere model results are the same when using synthetic data.After adding 10% random errors in the synthetic data, the depth is within 8.8%, the shape factor is within 4.6%, and the amplitude coefficient is within 3.3%.In all cases examined, the exact values of the depth (z), the shape factor (q), and the amplitude coefficient (A) were obtained when using synthetic data without random errors.However, in studying the error response of the present method, synthetic examples contaminated with 10% random errors were considered.Good results are obtained by using the present algorithm-particularly for shape and depth estimation, which is a primary concern in gravity prospecting and other geophysical work.

Effect of Random
Noise.We compute a gravity anomaly due to a sphere model (profile length = 40 km, q = 1.5, and A = 5000 mGal * km 2 ; station separation interval = 1 km) buried at different depths.The computed gravity anomaly g(x i ) was contaminated with random errors with a noise level of 10 mGal using the following equation: where Δg rand (x i ) is the contaminated anomaly value at x i .Following the interpretation method, ( 5), (7), and (8) were used to determine depth, shape factor, and amplitude  2. We verified numerically that the depth, shape factor, and amplitude coefficient are within 3.5%, 2%, and 11.9%, respectively.Good results are obtained by using the present algorithm because our technique is robust in the presence of noise.

Application to Composite
Anomalies.The composite gravity anomaly in mGal, shown in Figure 3, consists of the combined effects of an intermediate structure of interest (a 2D horizontal cylinder with z = 5 units, A = 750 mGal * unit, and a station separation of 1 depth unit) and an interference from neighboring rocks (a 3D semiinfinite vertical cylinder with z = 6 units, and A = 300 mGal; a station separation of 1 depth unit) The anomaly was computed using the following expression: In Figure 3, anomaly 1 is the anomaly due to the intermediate structure of interest, and anomaly 2 is the anomalies due to the interference from neighboring structures represented by a horizontal cylinder and a vertical cylinder, respectively.The composite gravity anomaly Δg(x i ) is also contaminated with a 5% random error.Following the same interpretation method, the result is shown in Table 4.The Table 3: Numerical results of the present method applied to the sphere synthetic example (q = 1.5, z = 5 km, and A = 500 mGal * km 2 ; profile length = 20 km; sampling interval = 1 km) without and with 10% random noise.

Using synthetic data
Using data with 10% random errors result is generally in very good agreement with the model parameters (q = 1, z = 5 units; A = 750 mGal * unit).Good results are obtained by using the present algorithm for model parameters which are of primary concern in gravity prospecting and other geophysical work.It is also emphasized that the method works well when dealing with gravity data having interference from neighboring anomalies and noise.

The Effect of the Offset in the Origin Point of the Gravity
Profile.When interpreting real gravity data, inaccurate selection of the origin point of the gravity profile can lead to errors in estimating the gravity parameters.In order to examine this effect, we have introduced some successive errors (δx) of 0, ±0.2, ±0.4, . .., ±1.5 units to the horizontal coordinate x of (1) of the gravity forward modeling calculations.The corresponding gravity dataset of a semiinfinite vertical cylinder model (q = 0.5, z = 5 units, and A = 1200 mGal; profile length = 20 units, N = 2 units, and M = 3 units) has been inverted following the same procedures.Maximum absolute errors in the gravity inverse parameters (z, q, and A) are found to be 53%, 106%, and 180%, respectively (Figure 4(a)).This figure shows that the algorithm proposed has recovered almost the true values of z, q, and A for each anomalous body at zero offset (δx = 0 unit).They also show that as δx increases, the absolute error in z, q, and A in general increases too, and the sign of the error in these parameters can change.In order to examine the accuracy and stability of the introduced algorithm on the combined effects of inaccurate origin selection and noise contamination, 5% random noise has been added to the various gravity dataset described in the previous paragraph, and then inverted.The maximum absolute errors in the gravity parameters (z, q, and A) obtained from inversion for the vertical cylinder models are found to be 71%, 99%, and 108%, respectively (Figure 4(b)).The analysis introduced demonstrates that Offset (in units) Model parameters: q = 0.5, A = 1200 mGal.Profile length = 20 units z = 5 units, N = 2 units, and M = 3 units Error in model parameters (%) Offset (in units) Model parameters: q = 0. the present method is stable and can estimate the gravity parameters with a reasonable accuracy depending upon the embedded noise level and the offset value in the origin point.

3.4.
The Sensitivity of the Removing Regional Field.The composite gravity anomaly (in mGal), shown in Figure 5, which consists of a local structure (sphere with z = 3 units and A = 750 mGal * units 2 ; profile length = 20 units) and a 2nd-order polynomial of regional structure.The model equation is Using a separation technique to remove the effect of the regional structure (graphical method; [5]) and reestimate the model parameters (z, q, and A) (Table 5).Table 5 shows that the results are still incorrect because the field is still contaminated by remaining regional field.By using different separation methods [14][15][16], the residual separated and inverted the residual anomaly (Table 5).

Field Example
The fast algorithm has been adapted for interpreting residual gravity anomalies related to three different types of structures, for example, a sphere, a vertical cylinder, and a horizontal cylinder.The standard error (μ) is used in this paper as statistical preference criterion in order to compare the observed and calculated values.This μ is given by the following mathematical relationship: where g(x i ) is the observed gravity values, and g c (x i ) is the calculated gravity values.A residual gravity field anomaly taken from Canada has been interpreted using the proposed technique in order to examine its applicability and stability.

Mobrun Sulfide Body.
The residual gravity anomaly profile (Figure 6) over the Mobrun sulfide body in Noranda, QC, Canada (after [20]) was digitized at an interval of 33.5 m.The method was applied to the anomaly profile using a sampling interval of 33.5 m to determine the model parameters of the buried structure using all successful combinations of N and M values.Then, we computed the standard error (μ) between the observed values and the values computed from estimated parameters z, q, and A for each N and M value.The results are shown in  of the optimum set (0.02 mGal).The optimum set is given at N = 67 m and M = 134 m.The best-fit model parameters are z = 33.34m, q = 0.72, and A = 59.1 mGal (Figure 6).This suggests that the shape of the ore body resembles a semiinfinite vertical cylinder, probably with a large radius.This is because the shape factor computed by the present method (0.72) is located between the shape factor of a perfect semiinfinite vertical cylinder (q = 0.5) and the shape factor of an infinite horizontal cylinder (q = 1.0).It is evident from the field example that our method gives good insight from gravity data of a short profile length concerning the nature of the source body.This is because the geologic situation is not complicated.The present method may not be applied to real data in complex geologic situations to obtain reliable or detailed information about the different shallow sources.This is true because each gravity measurement determines, at the station location, the sum of all effects from the surface downward.In complex geologic situations, the gravity profile is seldom a simple picture of a single isolated disturbance but always is a combination of two or more anomalies of shallow origin and very broad anomalies of regional nature, which may have their origin below the section within which the geologic interest lies.The shape and the depth to the top of the ore body obtained by the present method agree very well with those obtained from Roy et al. [21] and drilling information [20] (Table 7).

Conclusion
The problem of determining the appropriate depth, shape factor, and amplitude coefficient of a buried structure from

Figure 1 :
Figure 1: Residual gravity anomalies and schematic diagrams for various simple geometrical structures: (a) vertical cylinder, (b) horizontal cylinder, and (c) sphere.

A = 5000 mGal * km 2 Figure 2 :
Figure 2: The effect of the depth of burial of a sphere model on the inverted parameters for synthetic data contaminated by 10% random noise.

Figure 3 :
Figure 3: Noisy composite gravity anomaly consisting of the combined effects of an intermediate structure (horizontal cylinder with A = 750 mGal * unit and z = 5 units) [anomaly 1] and interference from neighboring structure (semiinfinite vertical cylinder with A = 300 mGal and z = 6 units) [anomaly 2].

5 ,AFigure 4 :
Figure 4: The effect of the offset in the origin point of a semiinfinite vertical cylinder model: (a) on the inverted model parameters for noise free synthetic data and (b) on the inverted model parameters for synthetic data contaminated by 5% random noise.
where e is a small predetermined real number close to zero.The source depth is determined by solving one nonlinear equation in z.Any initial guess for z works well because there is always one global minimum.Theoretically, two different values of N and M are enough to determine the depth.In practice, more than two values of N and M are preferable because of the presence of noise in the data.

Table 1 :
Numerical results of the present method applied to the semiinfinite vertical cylinder synthetic example (q = 0.5, z = 3 km, and A = 100 mGal; profile length = 20 km; sampling interval = 1 km) without and with 10% random noise.

Table 2 :
Numerical results of the present method applied to the horizontal cylinder synthetic example (q = 1.0, z = 4 km, and A = 300 mGal * km; profile length = 20 km; sampling interval = 1 km) without and with 10% random noise.

Table 4 :
Numerical results of the present method applied to the horizontal cylinder synthetic example (q = 1, z = 5 km, and A = 750 mGal * km; profile length = 40 km; sampling interval = 1 km) without and with 5% random noise.

Table 6
for the cases of N and M values where the μ difference between the modelled and observed data is less than 1.1 mGal.Also we computed the set of mean values.The set of mean values of the model parameters is rejected because it has a larger μ value (0.06 mGal) than the μ-value of the optimum set (0.02 mGal).Also we computed the set of mean values.The set of mean values of the model parameters is rejected because it has a larger σ value (0.055 mGal) than the μ-value

Table 5 :
Numerical results of the present method applied to composite gravity anomaly (in mGal), which consists of a local structure (sphere with z = 3 units and A = 750 mGal * units 2 ; profile length = 20 units) and a 2nd-order polynomial of regional structure.

Table 6 :
Numerical results of the present method applied to the Mobrun field example, Canada (best-fit in bold).

Table 7 :
Comparative results of Mobrun field example, Canada.