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

Authors: Khalid S. Essa
Journal: Journal of Geological Research
 (2012)

---

## Abstract

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 formfz=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.


---

## Body

## 1. 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-gmax  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 gmax  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–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 equationz=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–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.


## 2. 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 Figure1 and given in Abdelrahman et al. [17] as(1)g(xi,z,q)=Azm(xi2+z2)q,
where(2)A={43πGσR3,2πGσR2,πGσR2,m={1,1,0,q={32forasphere1forahorizontalcylinder12foraverticalcylinderR≪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). xi is the position coordinate, σ is the density contrast, G is the universal gravitational constant, and R is the radius.Residual gravity anomalies and schematic diagrams for various simple geometrical structures: (a) vertical cylinder, (b) horizontal cylinder, and (c) sphere.
(a)(b)(c)At the origin (xi=0), (1) gives the following relationship:(3)g(0)=Az2q-m.Using (1), we obtain the following normalized equation at xi=±N and xi=±M where N=1,2,3,… and M=1,2,3,…(4)g(N)g(0)=(z2N2+z2)q,g(M)g(0)=(z2M2+z2)q.LetF=(g(N)/g(0)) and T=(g(M)/g(0)) then from (4), we get:(5)z=e([(ln⁡F/ln⁡T)*(ln⁡(z2/(M2+z2)))]+ln⁡(N2+z2))/2,M≠N.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:(6)zf=f(zj),
where zj is the initial depth, and zf is the revised depth; zf will be used as the zj for the next iteration. The iteration stops when |zf-zj|≤e, 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.Once, the depth (z) is known, the shape factor (q) can be estimated from the following form:(7)qc=ln⁡Fln⁡(zc2/(N2+zc2)),
where zc is the estimated depth. Finally, knowing the shape factor (q), the amplitude coefficient (A) can be estimated from the following form:(8)Ac=g(0)zc2qc-m,
where qc is the estimated shape factor.For eachN 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.


## 3. Synthetic Examples

We computed three different residual gravity anomalies, each consisting of the effect of local structure (semiinfinite vertical cylinder, horizontal cylinder, and sphere). The model equations representing the models are(9)g1(xi)=100(xi2+32)0.5,g2(xi)=1200(xi2+42),g3(xi)=2500(xi2+52)1.5.The three gravity anomalies are shown in Figure1. 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: (10)grand(xi)=g(xi)[1+(RND(i)-0.5)*0.1],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.
Using synthetic dataUsing data with 10% random errorsM (km)Computed depth(km)Computed shape factorComputed amplitude factor(mGal)Computed depth(km)Computed shape factorComputed amplitude factor(mGal)2.003.000.50100.002.970.4794.333.003.000.50100.002.900.54104.014.003.000.50100.002.420.3865.535.003.000.50100.002.970.4794.406.003.000.50100.003.310.58134.607.003.000.50100.003.150.53112.768.003.000.50100.003.260.56126.739.003.000.50100.003.190.54118.4110.003.000.50100.003.060.50102.89Average (km)3.000.50100.003.030.51105.96Percent of error0.000.000.001.002.005.96Table 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.
M (km)Using synthetic dataUsing data with 10% random errorsComputed depth (km)Computed shape factorComputed amplitude factor (mGal)Computed depth (km)Computed shape factorComputed amplitude factor (mGal)2.004.001.00300.004.661.19350.253.004.001.00300.004.181.12349.384.004.001.00300.004.151.11327.145.004.001.00300.004.001.03331.936.004.001.00300.004.051.06362.687.004.001.00300.003.780.93336.278.004.001.00300.004.001.03329.879.004.001.00300.004.281.01318.8510.004.001.00300.004.491.11350.93Average (km)4.001.00300.004.171.07339.69Percent of error0.000.000.004.257.0013.23Table 3
Numerical results of the present method applied to the sphere synthetic example (q=1.5, z=5 km, and A=500 mGal * km2; profile length = 20 km; sampling interval = 1 km) without and with 10% random noise.
Using synthetic dataUsing data with 10% random errorsM (km)Computed depth (km)Computed shape factorComputed amplitude factor (mGal)Computed depth (km)Computed shape factorComputed amplitude factor (mGal)2.005.001.50500.002.791.39334.343.005.001.50500.004.481.17146.714.005.001.50500.004.991.44407.925.005.001.50500.003.001.20271.766.005.001.50500.005.071.49493.737.005.001.50500.005.071.49496.488.005.001.50500.005.301.62844.659.005.001.50500.005.251.59750.9810.005.001.50500.005.161.54602.08Average (km)5.001.50500.004.561.44483.18Percent of error0.000.000.008.804.663.36wheregrand(xi) is the contaminated anomaly value at xi, 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.Table1 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.


### 3.1. 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 * km2; station separation interval = 1 km) buried at different depths. The computed gravity anomaly g(xi)was contaminated with random errors with a noise level of 10 mGal using the following equation: (11)Δgrand(xi)=g(xi)+10(RND(i)-0.5),whereΔgrand(xi) is the contaminated anomaly value at xi.Following the interpretation method, (5), (7), and (8) were used to determine depth, shape factor, and amplitude coefficient, respectively. The percentage of error in model parameters was plotted against the model depth for comparison. Numerical results are shown in Figure 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.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.


### 3.2. Application to Composite Anomalies

The composite gravity anomaly in mGal, shown in Figure3, 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:(12)Δg(xi)=3750(xi2+25)(2Dhorizontalcylinder)+300(xi2+36)0.5(3Dverticalcylinder).Figure 3
Noisy composite gravity anomaly consisting of the combined effects of an intermediate structure (horizontal cylinder withA=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].In Figure3, 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(xi) is also contaminated with a 5% random error. Following the same interpretation method, the result is shown in Table 4. The 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.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.
Using synthetic dataUsing data with 5% random errorsM (km)Computed depth (km)Computed shape factorComputed amplitude factor (mGal)Computed depth (km)Computed shape factorComputed amplitude factor (mGal)3.007.291.542252.785.380.76393.014.005.871.011812.756.911.222686.235.005.300.82590.165.930.91706.476.004.940.72435.485.730.85561.927.004.670.64359.255.460.78423.998.004.430.58211.875.020.67280.069.004.200.52179.184.500.54187.7110.004.000.47156.334.180.47152.74Average (km)5.090.79749.725.390.78674.02


### 3.3. 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 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.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.
(a)(b)


### 3.4. The Sensitivity of the Removing Regional Field

The composite gravity anomaly (in mGal), shown in Figure5, which consists of a local structure (sphere with z=3 units and A=750 mGal * units2; profile length = 20 units) and a 2nd-order polynomial of regional structure. The model equation is(13)Δg(xi)=2250(xi2+9)1.5(Asphericalmodel)-0.02xi2+2xi+10(2nd-orderregionalfield).Figure 5
Composite gravity anomaly of a buried sphere model (q=1.5, A=750 mGal * units2 and z=3 units; profile length = 20 units) and second-order regional.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–16], the residual separated and inverted the residual anomaly (Table 5).Table 5
Numerical results of the present method applied to composite gravity anomaly (in mGal), which consists of a local structure (sphere withz=3 units and A=750 mGal * units2; profile length = 20 units) and a 2nd-order polynomial of regional structure.
M (units)Composite synthetic gravity data invertedUsing graphical separation techniques, remaining residual data invertedUsing analytical method, the remaining residual data invertedz(units)qA (mGal * unit2)z (units)qA (mGal * unit2)z(units)qA (mGal * unit2)3.002.290.73135.652.540.99230.873.001.50750.004.002.040.61109.292.420.92195.123.001.50750.005.001.810.5194.782.300.85166.593.001.50750.006.001.590.4387.682.190.78145.153.001.50750.007.001.400.3785.452.080.72129.563.001.50750.008.001.220.3286.621.980.67118.343.001.50750.009.001.070.2790.421.890.63110.233.001.50750.0010.000.940.2496.491.810.59104.353.001.50750.00Average (unit)1.550.4398.302.150.77150.023.001.50750.00


## 3.1. 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 * km2; station separation interval = 1 km) buried at different depths. The computed gravity anomaly g(xi)was contaminated with random errors with a noise level of 10 mGal using the following equation: (11)Δgrand(xi)=g(xi)+10(RND(i)-0.5),whereΔgrand(xi) is the contaminated anomaly value at xi.Following the interpretation method, (5), (7), and (8) were used to determine depth, shape factor, and amplitude coefficient, respectively. The percentage of error in model parameters was plotted against the model depth for comparison. Numerical results are shown in Figure 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.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.


## 3.2. Application to Composite Anomalies

The composite gravity anomaly in mGal, shown in Figure3, 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:(12)Δg(xi)=3750(xi2+25)(2Dhorizontalcylinder)+300(xi2+36)0.5(3Dverticalcylinder).Figure 3
Noisy composite gravity anomaly consisting of the combined effects of an intermediate structure (horizontal cylinder withA=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].In Figure3, 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(xi) is also contaminated with a 5% random error. Following the same interpretation method, the result is shown in Table 4. The 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.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.
Using synthetic dataUsing data with 5% random errorsM (km)Computed depth (km)Computed shape factorComputed amplitude factor (mGal)Computed depth (km)Computed shape factorComputed amplitude factor (mGal)3.007.291.542252.785.380.76393.014.005.871.011812.756.911.222686.235.005.300.82590.165.930.91706.476.004.940.72435.485.730.85561.927.004.670.64359.255.460.78423.998.004.430.58211.875.020.67280.069.004.200.52179.184.500.54187.7110.004.000.47156.334.180.47152.74Average (km)5.090.79749.725.390.78674.02


## 3.3. 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 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.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.
(a)(b)


## 3.4. The Sensitivity of the Removing Regional Field

The composite gravity anomaly (in mGal), shown in Figure5, which consists of a local structure (sphere with z=3 units and A=750 mGal * units2; profile length = 20 units) and a 2nd-order polynomial of regional structure. The model equation is(13)Δg(xi)=2250(xi2+9)1.5(Asphericalmodel)-0.02xi2+2xi+10(2nd-orderregionalfield).Figure 5
Composite gravity anomaly of a buried sphere model (q=1.5, A=750 mGal * units2 and z=3 units; profile length = 20 units) and second-order regional.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–16], the residual separated and inverted the residual anomaly (Table 5).Table 5
Numerical results of the present method applied to composite gravity anomaly (in mGal), which consists of a local structure (sphere withz=3 units and A=750 mGal * units2; profile length = 20 units) and a 2nd-order polynomial of regional structure.
M (units)Composite synthetic gravity data invertedUsing graphical separation techniques, remaining residual data invertedUsing analytical method, the remaining residual data invertedz(units)qA (mGal * unit2)z (units)qA (mGal * unit2)z(units)qA (mGal * unit2)3.002.290.73135.652.540.99230.873.001.50750.004.002.040.61109.292.420.92195.123.001.50750.005.001.810.5194.782.300.85166.593.001.50750.006.001.590.4387.682.190.78145.153.001.50750.007.001.400.3785.452.080.72129.563.001.50750.008.001.220.3286.621.980.67118.343.001.50750.009.001.070.2790.421.890.63110.233.001.50750.0010.000.940.2496.491.810.59104.353.001.50750.00Average (unit)1.550.4398.302.150.77150.023.001.50750.00


## 4. 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: (14)μ=∑i=1N[g(xi)-gc(xi)]2N,whereg(xi) is the observed gravity values, and gc(xi) 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.


### 4.1. Mobrun Sulfide Body

The residual gravity anomaly profile (Figure6) 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 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 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.34 m, 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).Table 6
Numerical results of the present method applied to the Mobrun field example, Canada (best-fit in bold).
NM (m)Depthz(m)Shape factorqAmplitude coefficientA (mGal)Standard errorμ(mGal)33.5067.0025.680.6155.750.4933.50100.5026.540.6455.480.3433.50134.0028.660.7155.510.1367.00100.5028.640.6856.070.1667.00134.0033.340.7859.100.02100.50134.0042.360.9372.380.21−33.50−67.0024.890.5258.711.01−33.50−100.5029.530.6457.140.16−33.50−134.0030.740.6857.450.08−67.00−100.5041.300.8468.430.14−67.00−134.0040.140.8266.460.11−100.50−134.0038.150.7963.860.06Average (m)32.500.7260.530.06Table 7
Comparative results of Mobrun field example, Canada.
UsingGrant and West [20] methodUsing Roy et al. [21] methodUsing present methodz(m)3029.4433.3q—0.770.78A (mGal)——59.1Figure 6
The measured gravity response (black circles) over a Mobrun sulfide ore body in Noranda, QC, Canada, and the predicted response (white circles) computed from the present inversion method.


## 4.1. Mobrun Sulfide Body

The residual gravity anomaly profile (Figure6) 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 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 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.34 m, 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).Table 6
Numerical results of the present method applied to the Mobrun field example, Canada (best-fit in bold).
NM (m)Depthz(m)Shape factorqAmplitude coefficientA (mGal)Standard errorμ(mGal)33.5067.0025.680.6155.750.4933.50100.5026.540.6455.480.3433.50134.0028.660.7155.510.1367.00100.5028.640.6856.070.1667.00134.0033.340.7859.100.02100.50134.0042.360.9372.380.21−33.50−67.0024.890.5258.711.01−33.50−100.5029.530.6457.140.16−33.50−134.0030.740.6857.450.08−67.00−100.5041.300.8468.430.14−67.00−134.0040.140.8266.460.11−100.50−134.0038.150.7963.860.06Average (m)32.500.7260.530.06Table 7
Comparative results of Mobrun field example, Canada.
UsingGrant and West [20] methodUsing Roy et al. [21] methodUsing present methodz(m)3029.4433.3q—0.770.78A (mGal)——59.1Figure 6
The measured gravity response (black circles) over a Mobrun sulfide ore body in Noranda, QC, Canada, and the predicted response (white circles) computed from the present inversion method.


## 5. Conclusion

The problem of determining the appropriate depth, shape factor, and amplitude coefficient of a buried structure from the residual gravity data of a short or a long profile length can be solved using the present method. A simple and rapid inversion approach is formulated to use the anomaly values at the origin and two pairs of measured data points (±N and ±M). The repetition of the method using all possible combinations of such pairs of measured points will lead to the best-fitting model. This happens when these two pairs of points contain the least amount of noise in the entire set of measured data. It is also emphasized that the calculated gravity anomaly of a set of mean values of the model parameters obtained by the present method does not necessarily guarantee it matches the observed anomaly values when the data contain measurement errors. The advantages of this method over previous graphical and numerical techniques used to interpret gravity data are (1) all the three model parameters can be obtained from all observed data, (2) the method is automatic, and (3) the method works well even when gravity data was noisy. Moreover, the advantage of the present method over the least-squares method is that the method does not require computation of analytical or numerical derivatives with respect to the model parameters. Also, the disadvantages of this technique are very difficult to interpret more complicated structures and the accuracy of the results depending upon the removing unwanted field (regional). Finally, in view of the above facts, we envisage the application of this method in solving various problems related to potential field data interpretation in the future.



---

*Source: 327037-2012-03-22.xml*