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 form . 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.
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- 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 ββ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 equation 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 Figure 1 and given in Abdelrahman et al. [17] as whereIn (1), is the depth, 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 () approaches zero as the structure becomes a nearly horizontal bed, and approaches 1.5 as the structure becomes a perfect sphere (point mass). is the position coordinate, Ο is the density contrast, is the universal gravitational constant, and is the radius.
(a)
(b)
(c)
At the origin (), (1) gives the following relationship:
Using (1), we obtain the following normalized equation at and where and
Let and then from (4), we get:
Equation (5) can be solved for using the standard methods for solving nonlinear equations (e.g., [19]), and its iteration form can be expressed as: where is the initial depth, and is the revised depth; will be used as the for the next iteration. The iteration stops when , where is a small predetermined real number close to zero. The source depth is determined by solving one nonlinear equation in . Any initial guess for works well because there is always one global minimum. Theoretically, two different values of and are enough to determine the depth. In practice, more than two values of and are preferable because of the presence of noise in the data.
Once, the depth () is known, the shape factor () can be estimated from the following form: where is the estimated depth. Finally, knowing the shape factor (q), the amplitude coefficient (A) can be estimated from the following form: where is the estimated shape factor.
For each and value, we compute the values of the model parameters (, , and ) from (5), (7), and (8), respectively. Theoretically, the anomaly values at the origin and any two and distances are just enough to determine the model parameters. However, in practice, it is recommended to use all possible combinations of and 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 , , and . 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
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 and 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 () is the contaminated anomaly value at , 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.
3.1. Effect of Random Noise
We compute a gravity anomaly due to a sphere model (profile length = 40βkm, , and βmGalβ*βkm2; station separation interval = 1βkm) buried at different depths. The computed gravity anomaly was contaminated with random errors with a noise level of 10 mGal using the following equation:
where is the contaminated anomaly value at .
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.
3.2. 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 units, βmGalβ*βunit, and a station separation of 1 depth unit) and an interference from neighboring rocks (a 3D semiinfinite vertical cylinder with units, and 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 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 (, units; β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.
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 () 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 (, units, and βmGal; profile length = 20 units, units, and 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 , , and 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.
(a)
(b)
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 units and βmGalβ*βunits2; 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β16], the residual separated and inverted the residual anomaly (Table 5).
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:
where is the observed gravity values, and 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 (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 and values. Then, we computed the standard error (ΞΌ) between the observed values and the values computed from estimated parameters , , and for each and value. The results are shown in Table 6 for the cases of and 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 βm and βm. The best-fit model parameters are βm, , and 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 () and the shape factor of an infinite horizontal cylinder (). 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).
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 ( and ). 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.
Appendix
At the origin (), (A.1) gives the following relationship:
Using (A.1), we obtain the following normalized equation at and , where and
Let and by taking algorithm to both sides:
by dividing (A.4a) and (A.4b), we get
From (A.4c), we get
By taking an exponential to both sides, we get
Equation (A.5) can be solved for using the standard methods for solving nonlinear equations (e.g., [19]), and its iteration form can be expressed as
Acknowledgments
The author wishes to express his sincere thanks to Professor Steven Forman, the editor, Professor Michel Chouteau, and the reviewer for their excellent suggestions, keen interests and thorough review that improved the paper. Many thanks to Professor El-Sayed M. Abdelrahman, Geophysics Department, Faculty of Science, Cairo University, for his constant help and encouragement.