- Email: [email protected]
An iterative least-squares approach to depth determination from residual magnetic anomalies due to thin dikes
WPLlECl GEOPHYSICS ELSEVIER
Journal of Applied Geophysics 34 (1996) 213-220
An iterative least-squares approach to depth determination from residual magnetic anomalies due to thin dikes E.M. Abdelrahman,
S.M. Sharafeldin
Geophysics Department, Faculty of Science, Cairo University, Giza, Egypt
Received 30 March 1995; accepted 8 August 1995
Abstract We have developed a least-squares approach to depth determination from residual magnetic anomalies due to thin dikes. By defining the zero-anomaly distance x0 and the anomaly value T( 0) at the origin on the profile, the problem of depth determination from magnetic data has been transformed into the problem of finding a solution of a non-linear equation. Procdures are also formulated to estimate the amplitude coefficient and the index parameter. The method is applied to synthetic data with and without random noise. The obtained depth agrees with the model depth within f 4%. The percentage error approaches zero as the sampling interval decreases to 0.25 depth units. Error in the depth parameter due to errors in x,, and T( 0) was also studied through imposing + l-IO% errors on these parameters on a synthetic profile due to a thin dike. The maximum error in depth is f 12%.
1. Introduction
The aim of magnetic surveys is to locate rocks or minerals having unusu.al magnetic properties, which reveal themselves as anomalies in the intensity of the Earth’s magnetic field. Magnetic anomalies are often used in a qualitative way to assist regional geologic interpretations. However, sometimes an individual magnetic anomaly is found which stands out so clearly that it can be separated from the regional background and the neighboring interferences, and which is so simple in appearance that it can be considered as caused by a single magnetized body. In this case, quantitative methods of interpretation can be used to determine the parameters of the magnetized body by assuming a model with simple geometry. The model is considered realistic if the form and magnitude of the calculated magnetic effects are close to the observed residual magnetic anomalies. 0926-9851/96/$15.00 0 1996 Elsevier Science B.V. All rights reserved SSDZO926-9851(95)00017-8
The thin dike model is widely used in magnetic interpretations to find the depth, and the magnitude and direction of magnetization, of intrusive igneous rocks in the forms of dikes and veins which very frequently carry magnetic minerals of economic interest. Several graphical methods have been developed for interpreting the residual magnetic anomalies due to thin dikes (e.g. Parker Gay, 1963; Grant and West, 1965; Stanley, 1977; Prakasa Rao et al., 1986). These methods use characteristic points, distances, curves, and nomograms for interpretation. The drawback with most graphical methods is that they are subject to error because of the use of only a few points and distances or because of the matching with standardized curves, which leads to errors in the derived depths. On the other hand, several numerical methods have been described in the geophysical literature. Clearly, the most used of these methods are the Werner deconvolution (Werner, 1953; Hartman et al., 1971; Jain,
E.M. Abdelrahman,
214
SM. Sharafeldin
/Journal
1976) and the Euler deconvolution (Hood, 1965; Slack et al., 1967; Thompson, 1982; Reid et al., 1990). In these methods, the depth determination problem is transformed into a problem of finding a solution of a system of linear equations. The methods are sensitive both to errors in anomaly amplitude resolution and in determination of vertical and horizontal gradients which are highly sensitive to noise (Steenland, 1968). Moreover, one important characteristic of the Werner and Euler deconvolution methods is that any four to seven data points of the anomaly profile are sufficient to produce an estimate of the thin dike parameters. By using different sets of points, different parameters may be obtained. Under ideal conditions, all estimates should of course be consistent. In practice, however, due to noise, this is seldom the case (Silva, 1989). In this paper, a least-squares approach is presented to estimate the depth of a buried structure from residual magnetic anomalies. Using a previously published formula, the depth determination problem has been parameterized and transformed into a problem of finding a solution of a non-linear equation of the formf( z) = 0. The solution to this non-linear equation is obtained by minimizing a function in the least-squares sense. The method has been applied to synthetic data with and without random noise, and is tested on two field examples from Arizona and Brazil.
of Applied Geophysics
34 (1996) 213-220
P(x.2)
Fig. 1. Cross-sectional
view of the thin dike model (2-D).
total, vertical and horizontal fields were given by Parker Gay (1963, table 1, p. 169). At the origin (xi = 0), Eq. 1 gives the following relationship:
K=To
(2)
cos 8
where T(0)is the anomaly value at the origin (Fig. 2). Setting Eq. 1 to zero, we obtain the following equation: tan 8=-z
(3) x0
2. The method Let us consider a thin dike which extends to infinity in both strike directions and downdip (2-D). The requirements of “thinness” is satisfied only by making the depth much greater than the thickness of the dike. Following Parker Gay (1963) and Atchuta Rao et al. ( 1980), the general expression, T,for the magnetic anomaly either in total, vertical, or horizontal field at a point P along the x-axis (Fig. 1) of an arbitrary magnetized thin dike is given by: T(
Xi,z,
0) =
si;*e 1:*ys ‘,
zXi
Model parameters z= 3 units.
X (in arbitrary units)
m
i= 1, 2, 3, . . .. N
-40
(1) where z is the depth to the top of the body, xi is a discrete point along x where the observed anomaly is located, K is the amplitude coefficient, and 8 is the index parameter. The values of K and 8 for the anomalies in the
-60 Fig. 2. A typical total magnetic anomaly profile over a thin dike. The zero-anomaly distance, x0, the anomaly value at the origin T( 0) , the positions of maximum vahx (M) , and minimum value (m) , x, and x,,,, respectively, are illustrated.
KM. Abdelrahmun, S.M. Sharafeldin / Journal of Applied Geophysics 34 (1996) 213-220
where x0 is the zero anomaly distance (Fig. 2). Substituting Eqs. 2 and 3 in Eq. 1, we obtain the following non-linear equation in z:
T(O)W(w)
T(Xi,Z) =
(4)
X0 where:
w(xivz)=
z2(xo --Xi> (,.f
+
The unknown minimizing: 4(Z)
= 2
LL(Xi)
z2)
depth z in Eq. 4 can be obtained by
- T( O> w(xi,Z)
lx012
(5)
i=l
where L(xi) denotes the observed magnetic anomaly at Xi. Setting the derivative of &z) with respect to z to zero leads to the following equation after simplification:
i=l
Z=
(6)
~~(0)~(XO--Xi)'/XO(~+Z2)3 i=l
Eq. 6 can be solved for z using the standard methods for solving non-linear equations. Here, it is solved by an iterative method (Demidovich and Maron, 1973; Press et al., 1986). The iterative form of Eq. 6 is given as:
zf=f
(7)
where zj is the initial depth and z, is the revised depth. zf will be used as Zj for the next iteration. The iteration stops when 1zf- Zj 1I q, where q is a small predetermined real number close to zero. The source body depth is determined by solving one non-linear equation in z. Any initial guess for z works well because there is al ways one global minimum. The experience with the minimization technique for two or more unknowns is that it always produces good results from synthetic data with or without random noise. In the case of field data, good results may only be obtained when using very good initial guesses on the model parameters. The optimization problem for the depth parameter is highly non-linear; increasing the number of parameters to be solved simultaneously also
215
increases the dimensionality of the energy surface, thereby greatly increasing the probability of the optimization stalling in a local minimum on that surface. Thus, common sense dictates that the non-linear optimization should be restricted to as few parameters as is consistent with obtaining useful results. This is why we propose a solution for only one unknown, 2. Once z is known, the index parameter, 0, can be determined from Eq. 3. From 0, the effective angle of magnetization and the dip of the dike can be estimated. The value of 8 thus obtained lies between -90 and 90”, although in reality it can be any value between 0 and - 360” (Parker Gay, 1963, figs. 3-6). The correct value can be obtained from examination of the field profile. Positions of the dominant positive and the dominant negative determined the actual value by following the rules given in Table 1. Finally, knowing 0, the amplitude coefficient, K, can be determined from Eq. 2. To this point, we have assumed knowledge of the axes of the magnetic profile so that T( 0) and x0 can be determined. Otherwise, T( 0) and x0 can be determined using methods described by Stanley ( 1977). As illustrated in Fig. 2, the line M-m intersects the anomaly profile at xi = 0. The baseline of the anomaly profile lies a distance M - T( 0) above the minimum m. It can also be determined by using the fact that x0 is located midway between x, and x,. Stanley’s (1977) methods work not only for the total intensity anomaly but also for the anomalies in vertical and horizontal fields.
3. Synthetic examples
3.1. EfSect of random noise Numerical results for various test cases are shown in Tables 2 and 3. In each case, a search algorithm was used to determine x0 based on linear interpolations between observed values (Davis, 1973). The search Table 1 Criteria for determining
the actual value of 0
Position and sign of dominant anomaly
Actual value of 0
Dominant positive on the left of origin Dominant positive on the right of origin Dominant negative on the left or right of origin
B” e-360 0- 180
216
E.M. Abdelrahman, S.M. Sharafeldin /Journal of Applied Geophysics 34 (1996) 213-220
algorithm was needed to find x,, because the true x0 value might be changed due to adding the random error to the observed values. It was verified numerically that Eq. 6 gives accurate values of z when synthetic data with or without random noise are analyzed. The depth obtained is within + 4% (Table 2). It was also shown that Eqs. 2 and 3 give reliable values of K and 8 (Table 2). The amplitude coefficient and the index parameter obtained are within + 2%. The percentage of error in model parameters approaches zero as the sampling interval decreases to 0.25 depth units (Table 3). As the sampling interval decreases within the profile length, a more accurate x0 is obtained and z in Eq. 6 approaches the actual depth of the dike. This illustrates that our method is less sensitive to errors in the residual anomaly than the Werner and Euler deconvolution methods and therefore noisy data will still give reliable depth estimates.
3.2. Effect of error in x0 and T(0) In studying the error response of the least-squares method, synthetic examples of a dike (profile length = 80 units, z = 3 units, 8 = - 60”, K = 100 units, and sampling interval= 1 unit) were considered in which errors of k 1, k 2, f 3, . . . . + 10% were imposed in both x0 and T( 0). Following the proposed method, values of the three model parameters (z, K, and 13)were computed and the percentage of errors in the model parameters were mapped (Figs. 3, 4 and 5). When x0 and T(0) both have errors of equal magnitude and of the same signs simultaneously, the interpreted z and K will not differ much from the true values (Figs. 3 and 4). In the case of the index parameter, 0, the interpreted values will vary slightly from the true values (Fig. 5). The maximum error in the latter case is f 4%. On the other hand, when x0 and T(0) both possess errors of equal magnitude and opposite signs simultaneously, the interpreted z and K values will vary widely from the actual values. The maximum error in z and K is
Fig. 3. A map showing error response in depth estimates. Abscissa: % error imposed in x,,. Ordinate: % error imposed in T( 0). C.I. = 1%.
E.M. Abdelrahman, SM. Sharafeldin / Journal of Applied Geophysics 34 (1996) 213-220
Table 2 Synthetic examples for depth determination units, and sampling interval = 1 unit
(computed
from Eq. 6). Other model parameters
217
are: 0 = - 60”, K= 100 units, profile length = 80
Model depth (units)
Using synthetic data
% error in Z
% error in 8
% error in K
Using data with 10% random noise
% error in Z
% error in 6
% error in K
2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7
2.06 2.58 3.04 3.50 4.04 4.54 5.01 5.52 6.03 6.52 7.00
3.47 3.35 1.58 0.17 1.22 0.97 0.30 0.45 0.60 0.36 0.07
-1.08 -1.04 - 0.48 - 0.05 - 0.37 - 0.29 - 0.09 -0.13 -0.17 -0.10 - 0.02
- 1.39 - 1.85 -0.88 -0.10 - 0.67 -0.53 -0.16 - 0.25 -0.32 -0.19 - 0.04
2.07 2.59 3.04 3.51 4.04 4.56 5.00 5.53 6.04 6.50 6.99
3.80 3.60 1.64 0.55 1.24 1.44 0.00 0.68 0.67 0.11 -0.01
- 1.03 -0.76 -0.44 0.10 - 0.63 -0.13 - 0.20 - 0.03 -0.13 -0.12 - 0.07
- 1.83 -1.36 -0.83 0.19 - 1.14 - 0.25 - 0.37 -0.06 - 0.24 -0.23 -0.13
+ 12 and + 8%, respectively (Figs. 3 and 4). The computed index parameter, however, will not vary much from the true values (Fig. 5). Finally, when either x0
or T( 0) is noise free, the percentage error in the model parameters is smaller than the error imposed on the other.
Fig. 4. A map showing error response in amplitude coefficient estimates. Abscissa: C.I. = 1%.
% error imposed in x0. Ordinate: % error imposed in T(0).
218
E.M. Abdelrahman, S.M. Sharafeldin / Journal of Applied Geophysics 34 (1996) 213-220
Table 3 Synthetic examples showing the effect of sampling interval on model parameter 19= - 60”, K= 100 units, and profile length is 80 units
determination.
Using synthetic data
The actual model parameters
are: Z= 3 units,
Data with 10% random noise
Sampling interval (units)
No. of data points
% error in Z
% error in 0
% error in K
% error in Z
% error in 19
% error in K
1 0.5 0.25 0.125
81 161 321 641
1.58 0.59 0.03 0.01
- 0.48 -0.18 -0.01 -0.01
0.88 0.33 -0.02 -0.01
1.71 0.30 0.13 0.24
-0.71 -0.10 0.02 0.10
- 1.27 -0.30 0.02 0.03
4. Field examples
depth to the top of the body (64 m) . This profile of 750 m length was digitized at an interval of 25 m. The values of T(0) and x0 used in the present method are 359 nT and 50 m, respectively (Fig. 6). The model parameters determined are: z= 66 m, fI= - 53” and K = 596.5 nT. The depth agrees very well with the depth of 64 m obtained from drilling. The depth and the index parameter obtained by Parker Gay ( 1963) using standardized curves were 70 m and - 50”, respec-
4. I. Pima copper mine anomaly Fig. 6 shows a vertical magnetic anomaly from the Pima copper mine, Arizona (Parker Gay, 1963, fig. 10, p. 198), which represents the anomaly due to a thin dike. Drilling information showed that the mineralized zone to be 11 m thick which is much less than the actual
-8
-6
-4 IB
Fig. 5. A map showing error response in index parameter C.I. = 1%.
OF
-2
ERF
estimates. Abscissa:
% error imposed in x,,. Ordinate:
% error imposed in T(0).
LM Abdelrahman, SM. Sharafeldin / Journal of Applied Geophysics 34 (1996) 213-220
219
2ol nT
tively. The computed anomaly profile using the parameters obtained by our method also agrees with the observed profile (Fig. 6). The model misfits some measurable values because the field data contain measurement errors.
(1 &ma -
1.54
m)
x
4.2. Pamaiba dike anomaly Fig. 7 gives a total Imagnetic anomaly above a Mesozoic diabase dike intruded into Paleozoic sediments from the Parnaiba basin, Brazil (Silva, 1989, fig. 10). The depth to the outcropping dike (sensor height) is 1.9 m (Fig. 7). The :zero crossings were determined using Stanley’s ( 1977) method. This anomaly profile of 24.64 m length was digitized at an interval of 1.54 m. The values of T( 0) and x,, used in the least-squares method are -49 nT and - 5.39 m, respectively (Fig. 7). The model parameters obtained are: z = 3.5 m, 8= 33.3” and K= - 58.6 nT. Silva ( 1989) applied a method to the same magnetic data employing Mfitting technique and obtained a depth of about 3.5 m. The depth to the top is over estimated by the two methods. This is reasonable because the upper part of the dike was weathered and the magnetite present was oxidized, loosing most of its magnetic property (Silva, 1989). Moreover, Fig. 7 shows that there is a lack of “thinness”. However, because many of the characteristics of the thick dikle anomalies are similar to those for thin dikes (Parker Gay, 1963), our method can be applied not only to magnetic anomalies due to thin I600
I
. *. . .
Observed Computed
anomaly anomaly
(1 unil
25m)
00
200
.11* -16
-4
-200
I
--
cnorn0ly Comqu!cd cnomoly
Observed
-50
1
Fig. 7. Total magnetic anomaly (above) over an outcropping dike (below) in the Pamaiba basin, Brazil (Silva, 1989). The base line and the zero crossings are determined using Stanley’s (1977) method.
nT -
I~~~I~S1,II” -12 -6
_
=
dikes but also to those due to thick dikes to obtain reliable depth estimates. Finally, these two field examples emphasize one of the principle advantages of the least-squares methods: a reliable depth can be obtained in spite of irregularities in the anomaly curves (Figs. 6 and 7) that would more seriously affect methods of depth estimation based on only a few isolated points and distances of the curves.
X
I 4
0
12
16
j
Fig. 6. Vertical magnetic anomaly over the Pima copper deposit in Arizona. The base line and the zero crossings shown are the same as those given by Parker Gay ( 1963).
5. Discussion and conclusions The problem of determining the depth of a buried thin dike from the residual magnetic anomaly has been
220
E.M. Abdelrahman,
S.M. Sharafeldin / Journal of Applied Geophysics 34 (1996) 213-220
transformed into the problem of finding a solution of a non-linear equation. The presented method is very simple to apply. The advantages of the present method over previous graphical techniques, which used only a few points and distances or standardized curves, are ( 1) all observed values can be used, and (2) the method is automatic. The advantages of the present technique over the Werner and Euler deconvolution methods are ( 1) the depth estimates are always consistent, (2) measured or calculated horizontal and vertical gradients are not required and therefore noisy data will still give reliable depth estimates, particularly when a complete anomaly profile is used, and (3) the method is less sensitive to errors in the residual anomaly. Very good results can be obtained by decreasing the sampling interval. Finally, the technique can be extended to anomalies in the vertical and horizontal fields due to spheres and to anomalies in total, vertical, and horizontal fields due to long horizontal cylinders.
Acknowledgements The authors would like to thank the Editors and three capable reviewers for their excellent suggestions, valuable comments on the manuscripts, and improvement of this work. Thanks are also due to to Prof. Dr. Y.E. Abdelhady, Head of the Geophysics Department, Cairo University, for his constant encouragement.
References Atchuta Rao, D., Ram Babu, H.V. and Sankar Narayan, P.V., 1980. Relationship of magnetic anomaliei due to subsurface features
and the interpretation of sloping contacts. Geophysics, 45: 3236. Davis, J.G., 1973. Statistics and Data Analysis in Geology. Wiley, Chichester. Demidovich, B.P. and Maron, LA., 1973. Computational Mathematics. Mir, Moscow, 691 pp. Grant, F.S. and West, G.F., 1965. Interpretation Theory in Applied Geophysics. McGraw-Hill, New York, NY, 584 pp. Hartman, R.R., Teskey, D.J. and Friedberg, J.L., 1971. A system for rapid digital aeromagnetic interpretation. Geophysics, 36: 891918. Hood, P., 1965. Gradient measurements in aeromagnetic surveying. Geophysics, 30: 891-902. Jain, S., 1976. An automatic method of direct interpretation of magnetic profiles. Geophysics, 41: 53 1-541. Parker Gay, S., 1963. Standard curves for interpretation of magnetic anomalies over long tabular bodies. Geophysics, 28: 161-200. Prakasa Rao, T.K.S., Subrahmanyam, M. and Srikrishna Murthy, A., 1986. Nomograms for the direct interpretation of magnetic anomalies due to long horizontal cylinders. Geophysics, 51: 21562159. Press, W.H., Flannery, B.P., Teukolsky, S.A. and Vetterling, W.T., 1986. Numerical Recipes, The Art of Scientific Computing. Cambridge Univ. Press, Cambridge. Reid, A.B., Allsop, J.M., Gmser, H., Millett, A.J. and Somerton, I.W., 1990. Magnetic interpretation in three dimensions using Euler deconvolution. Geophysics, 55: 80-91. Silva, J.B.C., 1989. Transformation of nonlinear problems into linear ones applied to the magnetic field of a two-dimensional prism. Geophysics, 54: 114-121. Slack, H.A., Lynch, V.M. and Langan, L., 1967. The geomagnetic gradiometer. Geophysics, 32: 877-892. Stanley, J.M., 1977. Simplified magnetic interpretation of the geologic contact and thin dike. Geophysics, 42: 1236-1240. Steenland, N.C., 1968. Discussion on “The geomagnetic gradiometer” by H.A. Slack, V.M. Lynch and L. Langan (Geophysics, October 1967, 877-892). Geophysics, 33: 681-684. Thompson, D.T., 1982, EULDPH-a new technique for making computer-assisted depth estimates from magnetic data. Geophysics, 47: 31-37. Werner, S., 1953. Interpretation of magnetic anomalies at sheet-like bodies. Sver. Geol. Under. Ser. C. Arsb., 43: N:6.
Journal of Applied Geophysics 34 (1996) 213-220
An iterative least-squares approach to depth determination from residual magnetic anomalies due to thin dikes E.M. Abdelrahman,
S.M. Sharafeldin
Geophysics Department, Faculty of Science, Cairo University, Giza, Egypt
Received 30 March 1995; accepted 8 August 1995
Abstract We have developed a least-squares approach to depth determination from residual magnetic anomalies due to thin dikes. By defining the zero-anomaly distance x0 and the anomaly value T( 0) at the origin on the profile, the problem of depth determination from magnetic data has been transformed into the problem of finding a solution of a non-linear equation. Procdures are also formulated to estimate the amplitude coefficient and the index parameter. The method is applied to synthetic data with and without random noise. The obtained depth agrees with the model depth within f 4%. The percentage error approaches zero as the sampling interval decreases to 0.25 depth units. Error in the depth parameter due to errors in x,, and T( 0) was also studied through imposing + l-IO% errors on these parameters on a synthetic profile due to a thin dike. The maximum error in depth is f 12%.
1. Introduction
The aim of magnetic surveys is to locate rocks or minerals having unusu.al magnetic properties, which reveal themselves as anomalies in the intensity of the Earth’s magnetic field. Magnetic anomalies are often used in a qualitative way to assist regional geologic interpretations. However, sometimes an individual magnetic anomaly is found which stands out so clearly that it can be separated from the regional background and the neighboring interferences, and which is so simple in appearance that it can be considered as caused by a single magnetized body. In this case, quantitative methods of interpretation can be used to determine the parameters of the magnetized body by assuming a model with simple geometry. The model is considered realistic if the form and magnitude of the calculated magnetic effects are close to the observed residual magnetic anomalies. 0926-9851/96/$15.00 0 1996 Elsevier Science B.V. All rights reserved SSDZO926-9851(95)00017-8
The thin dike model is widely used in magnetic interpretations to find the depth, and the magnitude and direction of magnetization, of intrusive igneous rocks in the forms of dikes and veins which very frequently carry magnetic minerals of economic interest. Several graphical methods have been developed for interpreting the residual magnetic anomalies due to thin dikes (e.g. Parker Gay, 1963; Grant and West, 1965; Stanley, 1977; Prakasa Rao et al., 1986). These methods use characteristic points, distances, curves, and nomograms for interpretation. The drawback with most graphical methods is that they are subject to error because of the use of only a few points and distances or because of the matching with standardized curves, which leads to errors in the derived depths. On the other hand, several numerical methods have been described in the geophysical literature. Clearly, the most used of these methods are the Werner deconvolution (Werner, 1953; Hartman et al., 1971; Jain,
E.M. Abdelrahman,
214
SM. Sharafeldin
/Journal
1976) and the Euler deconvolution (Hood, 1965; Slack et al., 1967; Thompson, 1982; Reid et al., 1990). In these methods, the depth determination problem is transformed into a problem of finding a solution of a system of linear equations. The methods are sensitive both to errors in anomaly amplitude resolution and in determination of vertical and horizontal gradients which are highly sensitive to noise (Steenland, 1968). Moreover, one important characteristic of the Werner and Euler deconvolution methods is that any four to seven data points of the anomaly profile are sufficient to produce an estimate of the thin dike parameters. By using different sets of points, different parameters may be obtained. Under ideal conditions, all estimates should of course be consistent. In practice, however, due to noise, this is seldom the case (Silva, 1989). In this paper, a least-squares approach is presented to estimate the depth of a buried structure from residual magnetic anomalies. Using a previously published formula, the depth determination problem has been parameterized and transformed into a problem of finding a solution of a non-linear equation of the formf( z) = 0. The solution to this non-linear equation is obtained by minimizing a function in the least-squares sense. The method has been applied to synthetic data with and without random noise, and is tested on two field examples from Arizona and Brazil.
of Applied Geophysics
34 (1996) 213-220
P(x.2)
Fig. 1. Cross-sectional
view of the thin dike model (2-D).
total, vertical and horizontal fields were given by Parker Gay (1963, table 1, p. 169). At the origin (xi = 0), Eq. 1 gives the following relationship:
K=To
(2)
cos 8
where T(0)is the anomaly value at the origin (Fig. 2). Setting Eq. 1 to zero, we obtain the following equation: tan 8=-z
(3) x0
2. The method Let us consider a thin dike which extends to infinity in both strike directions and downdip (2-D). The requirements of “thinness” is satisfied only by making the depth much greater than the thickness of the dike. Following Parker Gay (1963) and Atchuta Rao et al. ( 1980), the general expression, T,for the magnetic anomaly either in total, vertical, or horizontal field at a point P along the x-axis (Fig. 1) of an arbitrary magnetized thin dike is given by: T(
Xi,z,
0) =
si;*e 1:*ys ‘,
zXi
Model parameters z= 3 units.
X (in arbitrary units)
m
i= 1, 2, 3, . . .. N
-40
(1) where z is the depth to the top of the body, xi is a discrete point along x where the observed anomaly is located, K is the amplitude coefficient, and 8 is the index parameter. The values of K and 8 for the anomalies in the
-60 Fig. 2. A typical total magnetic anomaly profile over a thin dike. The zero-anomaly distance, x0, the anomaly value at the origin T( 0) , the positions of maximum vahx (M) , and minimum value (m) , x, and x,,,, respectively, are illustrated.
KM. Abdelrahmun, S.M. Sharafeldin / Journal of Applied Geophysics 34 (1996) 213-220
where x0 is the zero anomaly distance (Fig. 2). Substituting Eqs. 2 and 3 in Eq. 1, we obtain the following non-linear equation in z:
T(O)W(w)
T(Xi,Z) =
(4)
X0 where:
w(xivz)=
z2(xo --Xi> (,.f
+
The unknown minimizing: 4(Z)
= 2
LL(Xi)
z2)
depth z in Eq. 4 can be obtained by
- T( O> w(xi,Z)
lx012
(5)
i=l
where L(xi) denotes the observed magnetic anomaly at Xi. Setting the derivative of &z) with respect to z to zero leads to the following equation after simplification:
i=l
Z=
(6)
~~(0)~(XO--Xi)'/XO(~+Z2)3 i=l
Eq. 6 can be solved for z using the standard methods for solving non-linear equations. Here, it is solved by an iterative method (Demidovich and Maron, 1973; Press et al., 1986). The iterative form of Eq. 6 is given as:
zf=f
(7)
where zj is the initial depth and z, is the revised depth. zf will be used as Zj for the next iteration. The iteration stops when 1zf- Zj 1I q, where q is a small predetermined real number close to zero. The source body depth is determined by solving one non-linear equation in z. Any initial guess for z works well because there is al ways one global minimum. The experience with the minimization technique for two or more unknowns is that it always produces good results from synthetic data with or without random noise. In the case of field data, good results may only be obtained when using very good initial guesses on the model parameters. The optimization problem for the depth parameter is highly non-linear; increasing the number of parameters to be solved simultaneously also
215
increases the dimensionality of the energy surface, thereby greatly increasing the probability of the optimization stalling in a local minimum on that surface. Thus, common sense dictates that the non-linear optimization should be restricted to as few parameters as is consistent with obtaining useful results. This is why we propose a solution for only one unknown, 2. Once z is known, the index parameter, 0, can be determined from Eq. 3. From 0, the effective angle of magnetization and the dip of the dike can be estimated. The value of 8 thus obtained lies between -90 and 90”, although in reality it can be any value between 0 and - 360” (Parker Gay, 1963, figs. 3-6). The correct value can be obtained from examination of the field profile. Positions of the dominant positive and the dominant negative determined the actual value by following the rules given in Table 1. Finally, knowing 0, the amplitude coefficient, K, can be determined from Eq. 2. To this point, we have assumed knowledge of the axes of the magnetic profile so that T( 0) and x0 can be determined. Otherwise, T( 0) and x0 can be determined using methods described by Stanley ( 1977). As illustrated in Fig. 2, the line M-m intersects the anomaly profile at xi = 0. The baseline of the anomaly profile lies a distance M - T( 0) above the minimum m. It can also be determined by using the fact that x0 is located midway between x, and x,. Stanley’s (1977) methods work not only for the total intensity anomaly but also for the anomalies in vertical and horizontal fields.
3. Synthetic examples
3.1. EfSect of random noise Numerical results for various test cases are shown in Tables 2 and 3. In each case, a search algorithm was used to determine x0 based on linear interpolations between observed values (Davis, 1973). The search Table 1 Criteria for determining
the actual value of 0
Position and sign of dominant anomaly
Actual value of 0
Dominant positive on the left of origin Dominant positive on the right of origin Dominant negative on the left or right of origin
B” e-360 0- 180
216
E.M. Abdelrahman, S.M. Sharafeldin /Journal of Applied Geophysics 34 (1996) 213-220
algorithm was needed to find x,, because the true x0 value might be changed due to adding the random error to the observed values. It was verified numerically that Eq. 6 gives accurate values of z when synthetic data with or without random noise are analyzed. The depth obtained is within + 4% (Table 2). It was also shown that Eqs. 2 and 3 give reliable values of K and 8 (Table 2). The amplitude coefficient and the index parameter obtained are within + 2%. The percentage of error in model parameters approaches zero as the sampling interval decreases to 0.25 depth units (Table 3). As the sampling interval decreases within the profile length, a more accurate x0 is obtained and z in Eq. 6 approaches the actual depth of the dike. This illustrates that our method is less sensitive to errors in the residual anomaly than the Werner and Euler deconvolution methods and therefore noisy data will still give reliable depth estimates.
3.2. Effect of error in x0 and T(0) In studying the error response of the least-squares method, synthetic examples of a dike (profile length = 80 units, z = 3 units, 8 = - 60”, K = 100 units, and sampling interval= 1 unit) were considered in which errors of k 1, k 2, f 3, . . . . + 10% were imposed in both x0 and T( 0). Following the proposed method, values of the three model parameters (z, K, and 13)were computed and the percentage of errors in the model parameters were mapped (Figs. 3, 4 and 5). When x0 and T(0) both have errors of equal magnitude and of the same signs simultaneously, the interpreted z and K will not differ much from the true values (Figs. 3 and 4). In the case of the index parameter, 0, the interpreted values will vary slightly from the true values (Fig. 5). The maximum error in the latter case is f 4%. On the other hand, when x0 and T(0) both possess errors of equal magnitude and opposite signs simultaneously, the interpreted z and K values will vary widely from the actual values. The maximum error in z and K is
Fig. 3. A map showing error response in depth estimates. Abscissa: % error imposed in x,,. Ordinate: % error imposed in T( 0). C.I. = 1%.
E.M. Abdelrahman, SM. Sharafeldin / Journal of Applied Geophysics 34 (1996) 213-220
Table 2 Synthetic examples for depth determination units, and sampling interval = 1 unit
(computed
from Eq. 6). Other model parameters
217
are: 0 = - 60”, K= 100 units, profile length = 80
Model depth (units)
Using synthetic data
% error in Z
% error in 8
% error in K
Using data with 10% random noise
% error in Z
% error in 6
% error in K
2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7
2.06 2.58 3.04 3.50 4.04 4.54 5.01 5.52 6.03 6.52 7.00
3.47 3.35 1.58 0.17 1.22 0.97 0.30 0.45 0.60 0.36 0.07
-1.08 -1.04 - 0.48 - 0.05 - 0.37 - 0.29 - 0.09 -0.13 -0.17 -0.10 - 0.02
- 1.39 - 1.85 -0.88 -0.10 - 0.67 -0.53 -0.16 - 0.25 -0.32 -0.19 - 0.04
2.07 2.59 3.04 3.51 4.04 4.56 5.00 5.53 6.04 6.50 6.99
3.80 3.60 1.64 0.55 1.24 1.44 0.00 0.68 0.67 0.11 -0.01
- 1.03 -0.76 -0.44 0.10 - 0.63 -0.13 - 0.20 - 0.03 -0.13 -0.12 - 0.07
- 1.83 -1.36 -0.83 0.19 - 1.14 - 0.25 - 0.37 -0.06 - 0.24 -0.23 -0.13
+ 12 and + 8%, respectively (Figs. 3 and 4). The computed index parameter, however, will not vary much from the true values (Fig. 5). Finally, when either x0
or T( 0) is noise free, the percentage error in the model parameters is smaller than the error imposed on the other.
Fig. 4. A map showing error response in amplitude coefficient estimates. Abscissa: C.I. = 1%.
% error imposed in x0. Ordinate: % error imposed in T(0).
218
E.M. Abdelrahman, S.M. Sharafeldin / Journal of Applied Geophysics 34 (1996) 213-220
Table 3 Synthetic examples showing the effect of sampling interval on model parameter 19= - 60”, K= 100 units, and profile length is 80 units
determination.
Using synthetic data
The actual model parameters
are: Z= 3 units,
Data with 10% random noise
Sampling interval (units)
No. of data points
% error in Z
% error in 0
% error in K
% error in Z
% error in 19
% error in K
1 0.5 0.25 0.125
81 161 321 641
1.58 0.59 0.03 0.01
- 0.48 -0.18 -0.01 -0.01
0.88 0.33 -0.02 -0.01
1.71 0.30 0.13 0.24
-0.71 -0.10 0.02 0.10
- 1.27 -0.30 0.02 0.03
4. Field examples
depth to the top of the body (64 m) . This profile of 750 m length was digitized at an interval of 25 m. The values of T(0) and x0 used in the present method are 359 nT and 50 m, respectively (Fig. 6). The model parameters determined are: z= 66 m, fI= - 53” and K = 596.5 nT. The depth agrees very well with the depth of 64 m obtained from drilling. The depth and the index parameter obtained by Parker Gay ( 1963) using standardized curves were 70 m and - 50”, respec-
4. I. Pima copper mine anomaly Fig. 6 shows a vertical magnetic anomaly from the Pima copper mine, Arizona (Parker Gay, 1963, fig. 10, p. 198), which represents the anomaly due to a thin dike. Drilling information showed that the mineralized zone to be 11 m thick which is much less than the actual
-8
-6
-4 IB
Fig. 5. A map showing error response in index parameter C.I. = 1%.
OF
-2
ERF
estimates. Abscissa:
% error imposed in x,,. Ordinate:
% error imposed in T(0).
LM Abdelrahman, SM. Sharafeldin / Journal of Applied Geophysics 34 (1996) 213-220
219
2ol nT
tively. The computed anomaly profile using the parameters obtained by our method also agrees with the observed profile (Fig. 6). The model misfits some measurable values because the field data contain measurement errors.
(1 &ma -
1.54
m)
x
4.2. Pamaiba dike anomaly Fig. 7 gives a total Imagnetic anomaly above a Mesozoic diabase dike intruded into Paleozoic sediments from the Parnaiba basin, Brazil (Silva, 1989, fig. 10). The depth to the outcropping dike (sensor height) is 1.9 m (Fig. 7). The :zero crossings were determined using Stanley’s ( 1977) method. This anomaly profile of 24.64 m length was digitized at an interval of 1.54 m. The values of T( 0) and x,, used in the least-squares method are -49 nT and - 5.39 m, respectively (Fig. 7). The model parameters obtained are: z = 3.5 m, 8= 33.3” and K= - 58.6 nT. Silva ( 1989) applied a method to the same magnetic data employing Mfitting technique and obtained a depth of about 3.5 m. The depth to the top is over estimated by the two methods. This is reasonable because the upper part of the dike was weathered and the magnetite present was oxidized, loosing most of its magnetic property (Silva, 1989). Moreover, Fig. 7 shows that there is a lack of “thinness”. However, because many of the characteristics of the thick dikle anomalies are similar to those for thin dikes (Parker Gay, 1963), our method can be applied not only to magnetic anomalies due to thin I600
I
. *. . .
Observed Computed
anomaly anomaly
(1 unil
25m)
00
200
.11* -16
-4
-200
I
--
cnorn0ly Comqu!cd cnomoly
Observed
-50
1
Fig. 7. Total magnetic anomaly (above) over an outcropping dike (below) in the Pamaiba basin, Brazil (Silva, 1989). The base line and the zero crossings are determined using Stanley’s (1977) method.
nT -
I~~~I~S1,II” -12 -6
_
=
dikes but also to those due to thick dikes to obtain reliable depth estimates. Finally, these two field examples emphasize one of the principle advantages of the least-squares methods: a reliable depth can be obtained in spite of irregularities in the anomaly curves (Figs. 6 and 7) that would more seriously affect methods of depth estimation based on only a few isolated points and distances of the curves.
X
I 4
0
12
16
j
Fig. 6. Vertical magnetic anomaly over the Pima copper deposit in Arizona. The base line and the zero crossings shown are the same as those given by Parker Gay ( 1963).
5. Discussion and conclusions The problem of determining the depth of a buried thin dike from the residual magnetic anomaly has been
220
E.M. Abdelrahman,
S.M. Sharafeldin / Journal of Applied Geophysics 34 (1996) 213-220
transformed into the problem of finding a solution of a non-linear equation. The presented method is very simple to apply. The advantages of the present method over previous graphical techniques, which used only a few points and distances or standardized curves, are ( 1) all observed values can be used, and (2) the method is automatic. The advantages of the present technique over the Werner and Euler deconvolution methods are ( 1) the depth estimates are always consistent, (2) measured or calculated horizontal and vertical gradients are not required and therefore noisy data will still give reliable depth estimates, particularly when a complete anomaly profile is used, and (3) the method is less sensitive to errors in the residual anomaly. Very good results can be obtained by decreasing the sampling interval. Finally, the technique can be extended to anomalies in the vertical and horizontal fields due to spheres and to anomalies in total, vertical, and horizontal fields due to long horizontal cylinders.
Acknowledgements The authors would like to thank the Editors and three capable reviewers for their excellent suggestions, valuable comments on the manuscripts, and improvement of this work. Thanks are also due to to Prof. Dr. Y.E. Abdelhady, Head of the Geophysics Department, Cairo University, for his constant encouragement.
References Atchuta Rao, D., Ram Babu, H.V. and Sankar Narayan, P.V., 1980. Relationship of magnetic anomaliei due to subsurface features
and the interpretation of sloping contacts. Geophysics, 45: 3236. Davis, J.G., 1973. Statistics and Data Analysis in Geology. Wiley, Chichester. Demidovich, B.P. and Maron, LA., 1973. Computational Mathematics. Mir, Moscow, 691 pp. Grant, F.S. and West, G.F., 1965. Interpretation Theory in Applied Geophysics. McGraw-Hill, New York, NY, 584 pp. Hartman, R.R., Teskey, D.J. and Friedberg, J.L., 1971. A system for rapid digital aeromagnetic interpretation. Geophysics, 36: 891918. Hood, P., 1965. Gradient measurements in aeromagnetic surveying. Geophysics, 30: 891-902. Jain, S., 1976. An automatic method of direct interpretation of magnetic profiles. Geophysics, 41: 53 1-541. Parker Gay, S., 1963. Standard curves for interpretation of magnetic anomalies over long tabular bodies. Geophysics, 28: 161-200. Prakasa Rao, T.K.S., Subrahmanyam, M. and Srikrishna Murthy, A., 1986. Nomograms for the direct interpretation of magnetic anomalies due to long horizontal cylinders. Geophysics, 51: 21562159. Press, W.H., Flannery, B.P., Teukolsky, S.A. and Vetterling, W.T., 1986. Numerical Recipes, The Art of Scientific Computing. Cambridge Univ. Press, Cambridge. Reid, A.B., Allsop, J.M., Gmser, H., Millett, A.J. and Somerton, I.W., 1990. Magnetic interpretation in three dimensions using Euler deconvolution. Geophysics, 55: 80-91. Silva, J.B.C., 1989. Transformation of nonlinear problems into linear ones applied to the magnetic field of a two-dimensional prism. Geophysics, 54: 114-121. Slack, H.A., Lynch, V.M. and Langan, L., 1967. The geomagnetic gradiometer. Geophysics, 32: 877-892. Stanley, J.M., 1977. Simplified magnetic interpretation of the geologic contact and thin dike. Geophysics, 42: 1236-1240. Steenland, N.C., 1968. Discussion on “The geomagnetic gradiometer” by H.A. Slack, V.M. Lynch and L. Langan (Geophysics, October 1967, 877-892). Geophysics, 33: 681-684. Thompson, D.T., 1982, EULDPH-a new technique for making computer-assisted depth estimates from magnetic data. Geophysics, 47: 31-37. Werner, S., 1953. Interpretation of magnetic anomalies at sheet-like bodies. Sver. Geol. Under. Ser. C. Arsb., 43: N:6.