A Fortran 77 computer code for damped least-squares inversion of Slingram electromagnetic anomalies over thin tabular conductors

doi:10.1016/j.cageo.2004.02.003

Computers & Geosciences

Volume 30, Issue 6, July 2004, Pages 591-599

https://doi.org/10.1016/j.cageo.2004.02.003 Get rights and content

Abstract

A FORTRAN 77 computer code is presented that permits the inversion of Slingram electromagnetic anomalies to an optimal conductor model. Damped least-squares inversion algorithm is used to estimate the anomalous body parameters, e.g. depth, dip and surface projection point of the target. Iteration progress is controlled by maximum relative error value and iteration continued until a tolerance value was satisfied, while the modification of Marquardt's parameter is controlled by sum of the squared errors value. In order to form the Jacobian matrix, the partial derivatives of theoretical anomaly expression with respect to the parameters being optimised are calculated by numerical differentiation by using first-order forward finite differences.

A theoretical and two field anomalies are inserted to test the accuracy and applicability of the present inversion program. Inversion of the field data indicated that depth and the surface projection point parameters of the conductor are estimated correctly, however, considerable discrepancies appeared on the estimated dip angles. It is therefore concluded that the most important factor resulting in the misfit between observed and calculated data is due to the fact that the theory used for computing Slingram anomalies is valid for only thin conductors and this assumption might have caused incorrect dip estimates in the case of wide conductors.

Introduction

Slingram (so-called Horizontal Loop Electromagnetic, HLEM) is a maximum coupled system which uses horizontal transmitter and receiver coil pairs kept at a fixed distance apart (usually 30, 60 or 90 m). The receiver system normally measures both in-phase and quadrature components of the secondary field as a percentage of the primary field intensity. The system is commonly used in exploring for conductive ore bodies and for groundwater exploration in fractured zones (Palacky et al., 1981; McNeill, 1990).

In electromagnetic methods, calculation of the response of a conductive target with arbitrary shape is more complicated than that in other geophysical methods. The structures for which an analytical solution can be obtained are rather limited, e.g. a thin perfectly conductive vertical or dipping dike, sphere and disk-shaped conductors (Grant and West, 1965). Duckworth et al. (1991) suggested a method for quantitative depth estimates by transforming Slingram anomalies to a form which is free from the effect of the coil separation, however, the interpretation of Slingram data is generally achieved using type curves or Argand diagrams (Lowrie and West, 1965; Nair et al., 1968; Parasnis, 1971; Hanneson and West, 1984). For this purpose, Ketola and Puranen (1967) offered a detailed model catalogue for the type curves for Slingram measurements.

As an adjunct to conventional analysis performed using type curves, we offer a FORTRAN 77 code named SLINV.FOR for the quantitative interpretation of Slingram anomalies. The program computes the parameters of a dipping conductive dike type mineral deposit from Slingram measurements using the damped least-squares inversion algorithm suggested by Marquardt (1963). The parameters obtained are burial depth (z), dip angle (α) measured from vertical and the surface projection point (X₀) of the vein (Fig. 1).

Forward modeling theory presented by Wesley (1958) was used in the least-squares algorithm. Since the theory is based on the perfectly conductive thin conductors, Slingram anomalies are purely in-phase and do not have a quadrature component. The theory has also been widely used by various researchers (i.e. Grant and West, 1965; Ketola and Puranen, 1967; Duckworth et al., 1993). Duckworth and Krebes (1995) proposed the expression of normalized secondary field at the receiver coil as a percentage of the secondary field using the notation provided by Wesley (1958) as $H(x)=100 L^{3} π tan^{−1} (a/L) L^{3} − π L^{3} + a p^{2} L^{2} + L^{3} π π 2q^{3} + tan^{−1} (b/q) q^{3} + b q^{2} p^{2} + L^{3} (q^{2} −3L^{2}) sin^{2} α π π 2q^{5} + tan^{−1} (b/q) q^{5} + b q^{4} p^{2} + 2L^{3} b sin^{2} α πq^{2} p^{4} [(x_{t} +x_{r})^{2} −L^{2} cos^{2} α]+ L^{3} 4πr_{r} r_{t} p^{2} [a+b cos 2α−c sin 2α]+ L^{3} sin α πp^{4} r_{r} [((x_{t} +x_{r})+L cos α)(d sin α+a cos α)]− L^{3} sin α πp^{4} r_{t} [((x_{t} +x_{r})−L cos α)(d sin α−a cos α)]+ L^{3} 2πr_{r} r_{t} p^{4} [(bd^{2} −ac^{2}) sin^{2} α+ abc sin 2α−ab(a+b) cos^{2} α],$ where L is the coil separation and α is the dip angle of the conductor (Fig. 1). Other notations used for simplifying Eq. (1) are as follows: $x_{r} =(x−X_{0}) cos α+z sin α+ L 2 cos α,x_{t} =(x−X_{0}) cos α+z sin α− L 2 cos α,z_{r} =z cos α−(x−X_{0}) sin α− L 2 sin α,z_{t} =z cos α−(x−x_{0}) sin α+ L 2 sin α,r_{r} = x_{r}^{2} +z_{r}^{2},r_{t} = x_{t}^{2} +z_{t}^{2},$ $a=2 r_{r} r_{t} cos φ_{r} −φ_{t} 2,b=2 r_{r} r_{t} cos φ_{r} +φ_{t} −3π 2,$ $c=2 r_{r} r_{t} cos φ_{r} +φ_{t} 2,d=2 r_{r} r_{t} cos φ_{r} −φ_{t} −3π 2,$ $p= L^{2} +a^{2},q= (x_{t} +x_{r})^{2} +(L sin α)^{2},$ $z_{r} <0 ⇒φ_{r} =− π 2 − tan^{−1} x_{r} z_{r},z_{r} >0 ⇒φ_{r} = π 2 − tan^{−1} x_{r} z_{r},z_{t} <0 ⇒φ_{t} =− π 2 − tan^{−1} x_{t} z_{t},z_{t} >0 ⇒φ_{t} = π 2 − tan^{−1} x_{t} z_{t},$ where x is the midpoint distance of transmitter and receiver coils. According to the forward solution theory of Wesley (1958) in Eq. (1), the dip angle of the conductor, α, should be in the range of 0–90°. If the dip is in the other direction (e.g. less than 0°), the direction of the profile should be reversed for the inversion process.

Section snippets

Inversion method

In inverse modeling, a geometrical model is chosen with an initial guess of the body parameters and then the process is iteratively advanced until a satisfactory fit is obtained between observed and calculated anomalies. Several methods, such as the gradient method, ridge regression, the Gauss method and singular value decomposition have been used by other authors to determine target parameters automatically, especially for the inversion of potential field data arising from anomalous bodies

Calculation of partial derivatives

In order to constitute the Jacobian matrix, we need to have partial derivatives of the analytical expression of Slingram anomaly with respect to the parameters being estimated. Because the analytical solution of the simple geometrical model on which the inversion will be performed is already known, the expressions for the partial derivatives are usually obtained by analytical differentiation in the inversion process. In the present study, however, the partial derivatives with respect to depth (z

Description of the program

The main program SLINV.FOR contains 4 functions and 12 subroutines, each indicated in italics in Fig. 3. The program reads coil separation (L), initial estimates of the parameters (p_i), a tolerance value and observed anomaly values (F_{X_i}^obs) with corresponding distances of observation points (x_i) from input data file SLINPUT.DAT and prints out the final iteration number (ITER), initial estimates and optimised parameters (p_i), observed and calculated data (F_{X_i}^obs and F_{X_i}^cal) with differences

Application to theoretical data

In order to test the efficiency and accuracy of the present inversion program, a theoretical anomaly is produced for a coil separation of L=60 m (200 ft) with the actual body parameters z=20 m, α=45° and X₀=150 m using Eq. (1). Fig. 4 shows the theoretical data and its inversion result. The final solution obtained at the 7th iteration with an RMS error of 0.0028 shows satisfactory fit between theoretical and optimised curves.

Application to field data

Two Slingram anomalies obtained over a dipping graphitic shale from

Results and discussion

The method presented may provide a useful basis for the quantitative interpretation of Slingram anomalies. It is evident from the results given in Fig. 5, Fig. 6 that the most reliable parameters provided by the method are estimates of depth (z) and surface projection point (X₀) of the conductor. The largest difference, however, between the results of Duckworth et al. (1991) and the present study was observed in dip estimates of the conductor.

Considering the assumptions of the theory used for

Acknowledgements

We thank Dr. Venkata Raju and one anonymous reviewer for their helpful and constructive comments on the earlier version of this paper.

References (18)

K Duckworth et al.
A Quickbasic program for the computation of the response of a thin tabular perfect conductor to a two coil horizontal coplanar electromagnetic prospecting system
Computers & Geosciences
(1995)
M.R Nair et al.
Experimental studies on the electromagnetic response of tilted half-planes to a horizontal loop prospecting system
Geoexploration
(1968)
C Thanassoulas et al.
Gravity inversion of a fault by Marquardt's method
Computers & Geosciences
(1987)
R.D Atchuta et al.
Inversion of gravity and magnetic anomalies over some bodies of simple geometric shape
Pure and Applied Geophysics
(1985)
K Duckworth et al.
A method for obtaining depth estimates from the geometry of Slingram profiles
Geophysics
(1991)
K Duckworth et al.
A coincident coil frequency-domain electromagnetic prospecting system
Canadian Journal of Exploration Geophysics
(1993)
F.S Grant et al.
Interpretation Theory in Applied Geophysics
(1965)
J.E Hanneson et al.
The horizontal loop electromagnetic response of a thin plate in a conductive earthPart II—computational results and examples
Geophysics
(1984)
Ketola, N., Puranen, M., 1967. Type curves for the interpretation of Slingram (horizontal loop) anomalies over tabular...

There are more references available in the full text version of this article.

Cited by (5)

Parameter estimation by Differential Search Algorithm from horizontal loop electromagnetic (HLEM) data
2018, Journal of Applied Geophysics
Citation Excerpt :
The so-called algorithm can also take into account the variation of conductor dip measured from the vertical and coil separation for the normalized secondary field at the location of the receiver coil. Dondurur and Sarı (2004) used the algorithm in a damped least-squares scheme to estimate the parameters of a dipping conductive dike-type mineral deposit from Slingram measurements. The detailed expression for the forward modeling and its notations can be found in both studies and Appendix A.
We present an efficient inversion tool for parameter estimation from horizontal loop electromagnetic (HLEM) data using Differential Search Algorithm (DSA) which is a swarm-intelligence-based metaheuristic proposed recently. The depth, dip, and origin of a thin subsurface conductor causing the anomaly are the parameters estimated by the HLEM method commonly known as Slingram. The applicability of the developed scheme was firstly tested on two synthetically generated anomalies with and without noise content. Two control parameters affecting the convergence characteristic to the solution of the algorithm were tuned for the so-called anomalies including one and two conductive bodies, respectively. Tuned control parameters yielded more successful statistical results compared to widely used parameter couples in DSA applications. Two field anomalies measured over a dipping graphitic shale from Northern Australia were then considered, and the algorithm provided the depth estimations being in good agreement with those of previous studies and drilling information. Furthermore, the efficiency and reliability of the results obtained were investigated via probability density function. Considering the results obtained, we can conclude that DSA characterized by the simple algorithmic structure is an efficient and promising metaheuristic for the other relatively low-dimensional geophysical inverse problems. Finally, the researchers after being familiar with the content of developed scheme displaying an easy to use and flexible characteristic can easily modify and expand it for their scientific optimization problems.
MRQP: A windows-based mixed-language program for mineral resource quantitative prediction
2008, Computers and Geosciences
Dynamic mechanisms of the post-seismic deformation following large events: Case study of the 1999 Chi-Chi earthquake in Taiwan of China
2009, Science in China, Series D: Earth Sciences
Inversion of viscous properties of crust and mantle from the GPS temporal series measurements
2006, Acta Geophysica Sinica
Depth estimates for slingram electromagnetic anomalies from dipping sheet-like bodies by the normalized full gradient method
2005, Pure and Applied Geophysics

^☆: Code available from server at http://www.iamg.org./CGEditor/index.htm

View full text

A Fortran 77 computer code for damped least-squares inversion of Slingram electromagnetic anomalies over thin tabular conductors☆

Abstract

Introduction

Section snippets

Inversion method

Calculation of partial derivatives

Description of the program

Application to theoretical data

Application to field data

Results and discussion

Acknowledgements

Computers & Geosciences

Geoexploration

Computers & Geosciences

Inversion of gravity and magnetic anomalies over some bodies of simple geometric shape

Pure and Applied Geophysics

A method for obtaining depth estimates from the geometry of Slingram profiles

Geophysics

A coincident coil frequency-domain electromagnetic prospecting system

Canadian Journal of Exploration Geophysics

Interpretation Theory in Applied Geophysics

The horizontal loop electromagnetic response of a thin plate in a conductive earthPart II—computational results and examples

Geophysics