Practical Inversion of 3D Time Domain Data Using EH3DT

1
Practical Inversion of 3D Time Domain Data Using
EH3DT The San Nicolás deposit, Zacatecas,
Mexico Scott Napier, Doug Oldenburg and Eldad
Haber
The Geophysical Inversion Facility University of
British Columbia
The Geophysical Inversion Facility University of
British Columbia
Introduction
Evaluating preliminary results
Evaluating final results
gradient algorithm. The inversion algorithm
is described in Haber, Oldenburg and Shekhtman
(2006). The system of equations for forward
modelling can be abbreviated to A(m)u q where m
is the conductivity model of the earth, A(m) is
the operation of Maxwell's equations on that
model, u contains the fields and q contains
sources and boundary conditions. The predicted
data are sampled from the forward modeled fields
in u by an operator Q. The inverse problem is
defined as a minimization of the functional F in
equation 2. This minimization is achieved by the
Gauss-Newton method. In the interest of
computational efficiency the linear system for
the step to solve equation 2 is pre-conditioned
with a matrix composed of the L-BFGS vectors of
the Quasi-Newton method.
The interpretation of surface time domain data
has been hampered over the years because of the
inability to invert data to recover a 3D
conductivity model. Interpretations have often
been limited to simplistic plate models, 2D
inversions, or combining the results of 1D
inversions. For many geologic situations these
approaches have proved to be insufficient. The
purpose of this poster is to show
(1)
Preliminary inversions of the San Nicolás
deposit located a conductor at the correct
lateral location but the depth was far too
shallow and the conductivity much lower than
expected. Additionally the overburden layer was
not present and undesirable artifacts existed in
the recovered model. The artifacts appear to
be closely related to the geometry of the
transmitters and receivers. Figure 7 shows a
typical initial inversion result. This effect is
a prime example of the non-uniqueness of the 3D
EM inversion problem. The San Nicolás UTEM
survey has only a few source geometries. This
combined with the extremely high variation in
sensitivity of cells throughout the model led to
these deficiencies in the recovered model. A
method of mitigating these problems is required
that time domain electromagnetic (TDEM) data
collected in a field survey can be inverted to
recover a 3D distribution of electrical
conductivity. That the inversion of TDEM data can
yield substantial benefits for interpreters and
that the process of 3D inversion of TDEM data is
practical and can become routine in exploration
work.
Figure 11 The data misfit as a function of
iteration number for the final inversion
Figure 10 shows the data fit for the large loop
transmitter and Figure 12 shows both north-south
and east-west cross sections through the inverted
result The results are excellent. The recovered
model shows an overburden layer, and a deep
conductor. They are distinct and their
conductivities are similar to the conductivities
surveyed from hand samples of drill core. Note
the excellent correspondence in both lateral
location and depth between drilling and the
inverted conductivity sections. The use of
multiple transmitters along with sensitivity
weighting helps to mitigate the non-uniqueness
problems encountered in early inversions.
Additionally the use of a conductive layer in the
reference model based on the overburden improves
the inversion results vastly, resolving the
overburden and target accurately.
Figure 7 Preliminary inversion result showing
loop geometry (red), receiver locations (black)
and an outline of the sulphide deposit (red).
Note the conductivity recovered is too shallow,
and artifacts on the surface appear closely
related to survey geometry
Figure 1 Map showing the location of the San
Nicolás deposit in the state of Zacatecas, Mexico
(2)
Geologic background
A practical inversion procedure
Unit
Density Susceptibility Resistivity
Chargeability
Unit
Density Susceptibility Resistivity
Chargeability
-
3
(g/cc) (S.I. x10
-
3
) (ohm
-
m) (msec)
(g/cc) (S.I. x10
) (ohm
-
m) (msec)
Inversion of 3D Time domain is a complex process.
For this reason we have developed a simple
procedure in order to assist with the process.
Following the steps outlined here, along with
careful attention to detail will help to achieve
successful inversions. In the following we will
discuss the various steps in this inversion
process within the context of the San Nicolás
deposit and show how following this simple
procedure can lead to excellent inversion
results.
Qal
2 0
-
10
50 5
Qal
2.0 0
-
10
50 5
Sensitivity weighting
-
Tv
2.3 0
-
5(20)
20
-
30
10
-
30
Tv
2.3 0
5(20)
20
-
30
10
-
30
150 20
-
40
Mst./Lst
.
2.4 0
150 20
-
40
Mst./Lst
.
2.4 0
Mafic Vol.
2.7 0
-
5
80 30
-
50
Mafic Vol.
2.7 0
-
5
80 30
-
50
Mafic/Int
Vol.
2.7 0
-
5
80 30
-
50
Mafic/Int
Vol.
2.7 0
-
5
80 30
-
50
3.5 10
20 200
3.5 10
20 200
Sulphide
Sulphide
We use the method of sensitivity weighting to
address the non-uniqueness issues discussed in
the previous section. We follow the formulation
of Chen et al (2002) based on the current density
in a representative background halfspace. Under
Qtz
. Rhyolite
2.4 0
-
10
100 10
-
20
Qtz
. Rhyolite
2.4 0
-
10
100 10
-
20
Graphitic
Mst
.
2.4 0
-
5
100 30
-
70
Graphitic
Mst
.
2.4 0
-
5
100 30
-
70
Figure 2 Left - showing a cross section through
the deposit at 400m S. Right - a table of
physical properties for geologic units in the
survey area.
Figure 10 a) predicted data. b) observed data
for the inversion result shown in Figure 8. c)
the decay curve shown in is for the data location
indicated by the star in a)
some simplifying assumptions the approximate UTEM
sensitivity, J can calculated analytically. Using
this approximation we can calculate a sensitivity
weight for each cell in the model via equation
3. By scaling this quantity to an
appropriate value using a parameter ? and
applying a threshold to small values we can
achieve a sensible value for a weighting whose
elements compose the matrix Wm from equation 2.
Figure 8 shows the results of a synthetic
inversion with and without sensitivity weighting.
The sensitivity weighted inversion has the
correct conductivity, much better depth
San Nicolás is a Copper-Zinc volcanogenic
massive sulphide deposit located in central
Mexico in the state of Zacatecas. The deposit is
a continuous but geometrically complex body of
sulphides which is covered by 175-250 m of
variable composition overburden. The local
geology is also somewhat complex and contains
numerous sedimentary and volcanic units. An
east-west cross section through the deposit at
450 m south, as interpreted from a drilling
program is shown in Figure 2 . Major geologic
units are labeled.
(3)
Figure 4 A simplified procedure for 3D inversion
Pre-inversion considerations
Geophysical background
The San Nicolás deposit has a long history of
being a test case for the application of
geophysics to the exploration of massive sulphide
deposits. The deposit has been investigated in
the past using a number of geophysical surveys.
Additionally the deposit has been thoroughly
drilled and the recovered drill core has been
logged for physical properties. As a result, both
the geology and physical properties of the
deposit are reasonably well understood. The
sulphide deposit presents an electrical
conductivity contrast with most of the geologic
units in the area, however the main overburden
unit has a conductivity close to that found in
the sulphide. This makes EM interpretation at San
Nicolás challenging. A UTEM survey was
conducted over the San Nicolás deposit in
December of 1998. Figure 3 shows the location of
both receivers and transmitters. The UTEM time
domain method uses a loop transmitter with a
saw-tooth type current waveform. Receivers
typically measure dB/dt. For the 3D inversion of
the San Nicolás survey we consider only z
component data as we found other components
measured tended to be extremely noise affected.
It is necessary to
The process of discretization requires an
estimate of a representative background model.
The background model determines the way in which
both the transmitter waveform in time and the
spatial domain are discretized.
Figure 8 Synthetic inversion using SN large loop
geometry (Figure 3, blue). a) Cross section
through block model consisting of 2 ?m block
buried at 200 m depth b) Inverted section without
sensitivity weighting. c) Inverted model with
sensitivity weight
Figure 12 a) EW cross sections through inverted
model at 450 m south. b) NS cross section through
inverted model at 1380 m west. c) EW cross
section though conductivity model from drilling
at 450m south. d) NS cross section through
conductivity model from drilling at 1380 m west.
resolution and attenuated artifacts at the
surface. The synthetic inversion is also useful,
as it illustrates the depth of investigation of
the UTEM survey. It is clear that the survey is
not sensitive to the conductor below about 400 m
depth.
Conclusions
Figure 5 a) discretization of the current
waveform in time. 38 time steps. Fields are
sampled and compared to UTEM data at the times
marked in blue b) large forward modeling mesh,
330,000 cells. c) inversion mesh for final
results, 240,000 cells. Preliminary inversions
used a much coarser inversion mesh of 90,000
cells
Final inversion
The example has clearly shown that 3D time domain
inversion can provide exceptional results in
areas where traditional interpretation of 3D EM
data becomes difficult if not impossible. The
result is an achievement unmatched by any other
EM interpretation at the site and we have clearly
demonstrated the considerable value of 3D TDEM
inversion at San Nicolás . In a more general
sense and perhaps much more importantly we have
shown that, though 3D time domain inversion is
complicated, with the use of a simple procedure
and careful attention to detail the process can
become commonplace and routine.
convert the UTEM data into a format suitable for
3D inversion. This requires an intimate
understanding of the data, its normalizations,
and its orientation conventions. This tends to an
overlooked aspect of inversion but the importance
of thoroughly understanding the data can not be
underestimated
Given a background model it is then possible
develop a discretization in time and in space.
The key to discretization is that it should
ensure that the accuracy of the forward modelling
is such that errors due to modelling are well
below the expected response to the target. At
UBC-GIF we have developed a set of rules of
thumb to assist with discretization, in both
space and time. Because of the complexity of
the inversion computation and the requirement
that we satisfy boundary conditions, most
practical surveys will require more cells than is
practical given todays computing resources. For
this reason we employ a method of reducing the
volume. Here the background model is used to
calculate a corrective source term on a smaller
mesh that reproduces the results
The final inversion used 3 transmitters and data
from 9 separate time channels to give 3523 data.
The time discretization had 38 time steps and the
inversion mesh had 240,000 cells. The actual
discretizations are shown in Figures 4a and 4c
respectively. Convergence took 33 iterations over
9 values of the regularization parameter ?. A
refined reference model was also used in this
inversion which consisted of a conductive, 20 ?m,
layer over a halfspace. The layer was 90m thick.
This corresponds to average overburden thickness
and conductivity. This type of a priori
information is precisely of the sort that should
be available to a geophysicist in an exploration
setting.
Acknowledgments
We would like to thank Teck Cominco Ltd. for
providing, and SJ Geophysics Ltd. for collecting
the UTEM data. We would like to thank Nigel
Phillips for his extensive work done on the San
Nicolás deposit as a part of his Master's
thesis.
for the large forward modelling mesh. This allows
the inversion to take place on a vastly reduced
domain. The data were divided into separate time
channels. and standard deviations were assigned
to the data based on a percentage of
Figure 3 Map showing the survey geometry for
the UTEM data to be inverted. The Tx loops are
shown in red, blue and black. The data from each
loop is shown in the corresponding colour.
The inversion algorithm
References
The details of the forward modeling
techniques used in EH3DT are outlined in Haber et
al (2004) and the references therein. Beginning
with Maxwell's equations in equation 1, initial
conditions E0, H0 and boundary conditions, n H
0, the domain is discretized in both time and
space. The initial conditions are propagated
through time using a Backward Differentiation
Formula (BDF). Given the fields at time tn, this
gives the fields everywhere inside the
discretized volume at time tn1. The electric and
magnetic fields are decomposed into potentials
and the resulting system of linear equations are
solved with a conjugate
Chen, J. , Oldenburg, D.W. , Haber, E., and
Bishop, J. , 2002, The down-hole magnetometric
resistivty method A 3-D inversion scheme,
unpublished work
Haber E., Oldenburg D. W., and Shekhtman R. 2006
Inversion of Time domain 3D Electromagnetic data,
accepted GJI
Figure 6 Comparison of predicted data from EH3DT
with an independent 1D code. This ensures that
the discretization used provides accurate
results.
Haber, E., Ascher, U., and Oldenburg, D. W.,
2004, Inversion of 3D electromagnetic data in
frequency and time using an inexact all-at-once
approach Geophysics, 69, 12161228.
the datum and a floor value. In consideration of
known modelling errors along with an analysis of
measurement errors the percentage used was higher
for early channels and lower for late time
channels
Figure 9 20 ?m iso-surface of inverted result
showing overburden layer and deposit
Phillips, N, 2001, Geophysical inversion in an
integrated exploration environment Examples from
the San Nicolás deposit Ph.D. thesis,
University of British Columbia.
2
1
3
4