GEOPHYSICS, VOL. 61, NO. 2 (MARCH-APRIL 1996); P. 394-408, 18 FIGS.

3-D inversion of magnetic data

Yaoguo Li* and Douglas W. Oldenburg*

inherent nonuniqueness. By Gauss’ theorem, if the field distri-

ABSTRACT bution is known only on a bounding surface, there are infinitely
many equivalent source distributions inside the boundary that
We present a method for inverting surface magnetic
can produce the known field. Any magnetic field measured on
data to recover 3-D susceptibility models. To allow the
the surface of the earth can be reproduced by an infinitesimally
maximum flexibility for the model to represent geologi-
thin zone of magnetic dipoles beneath the surface. From a
cally realistic structures, we discretize the 3-D model
mathematical perspective, this means there is no depth reso-
region into a set of rectangular cells, each having a
lution inherent in magnetic field data. A second source for
constant susceptibility. The number of cells is generally
far greater than the number of the data available, and nonuniqueness is the fact that magnetic observations are finite
thus we solve an underdetermined problem. Solutions in number and are inaccurate. If there exists one model that
are obtained by minimizing a global objective function reproduces the data, there are other models that will repro-
composed of the model objective function and data duce the data to the same degree of accuracy. The severity of
misfit. The algorithm can incorporate a priori informa- the nonuniqueness problem for magnetic data is illustrated in
tion into the model objective function by using one or Figures l-3. (The gray scale in all figures indicates suscepti-
more appropriate weighting functions. The model for bility in SI units for model sections and magnetic data in nT for
inversion can be either susceptibility or its logarithm. If data plots.) A 3-D dipping prism of uniform susceptibility in
susceptibility is chosen, a positivity constraint is imposed Figure 1 produces the surface magnetic field shown in Figure 2,
to reduce the nonuniqueness and to maintain physical which consists of 441 data. Slices of a 3-D susceptibility model
realizability. Our algorithm assumes that there is no that adequately reproduces the 441 data are shown in Figure 3.
remanent magnetization and that the magnetic data are That result, however, bears little resemblance to the true
produced by induced magnetization only. All minimiza- model. Susceptibility is concentrated near the surface and
tions are carried out with a subspace approach where displays zones of negative values. This mathematical model
only a small number of search vectors is used at each solution provides little information about the true structure
iteration. This obviates the need to solve a large system that is useful.
of equations directly, and hence earth models with many Faced with this extreme nonuniqueness, previous authors
cells can be solved on a deskside workstation. The have mainly taken two approaches in the inversion of magnetic
algorithm is tested on synthetic examples and on a field data. The first is parametric inversion, where the parameters of
data set. a few geometrically simple bodies are sought in a nonlinear
inversion and values are found by solving an overdetermined
problem. This methodology is suited for anomalies known to
INTRODUCTION be generated by simple causative bodies, but it requires a great
deal of a priori knowledge about the source expressed in the
Magnetic surveying has been used widely over the years, form of an initial parameterization, an initial guess for param-
resulting in a great amount of data with enormous area1 eter values, and limits on the susceptibility allowed (e.g.,
coverage. Magnetic data have been used for mapping geolog- Bhattacharyya, 1980; Zeyen and Pous, 1991). Nonuniqueness
ical structures, especially in the reconnaissance stage of explo- is not generally an issue because only a small subset of possible
ration, but when used in detailed prospecting, robust and models is considered due to the restrictive nature of the
efficient inversion algorithms must be used. However, a prin- inversion algorithm. A related, but unique, approach in Wang
cipal difficulty with the inversion of the potential data is the and Hansen (1990) assumes polyhedronal causative bodies and

Presented at the 63rd Annual International Meeting, Society of Exploration Geophysicists. Manuscript received by the Editor May 2, 1994; revised
manuscript received June 29, 1995.
*UBC-Geophysical Inversion Facility, Dept. of Geophysics and Astronomy, University of British Columbia, 129-2219 Main Mall, Vancouver, BC
V6T 1Z4, Canada.
© 1996 Society of Exploration Geophysicists. All rights reserved.

394
 3-D Inversion of Magnetic Data 395
inverts for the position of the vertices of these bodies using the able information. Last and Kubik (1983) choose to minimize
spectrum of the magnetic data. The method is general in the total volume of the causative body so that the final model
principle but has difficulties both in constructing the causative is compact and structurally simple. Guillen and Menichetti
bodies from the recovered vertices and in obtaining the (1984) minimize the moment of inertia of the causative body
susceptibility distribution. with respect to the center of gravity or an axis passing through
In the second approach to inverting magnetic data, the earth it. Their inversion result is guided by the estimate of the central
is divided into a large number of cells of fixed size but of depth and dip of the causative body. These approaches have
unknown susceptibility. Nonuniqueness of solution is recog- merit but they are not flexible enough to handle problems we
nized and the algorithm produces a single model by minimizing are concerned with. This is especially true of methods that
an objective function of the model subject to fitting the data. attempt to collapse the anomalous susceptibility into a single
Green (1975) minimizes a weighted model norm with respect body; such a solution is rarely an adequate representation of
to a reference model, and this allows the interpreter to guide geologic structure.
the inversion by varying the weighting according to the avail- In our inversion approach, we first make a decision about
the variable in which the interpretation is to be made, that is,
whether susceptibility, log susceptibility, or some function of
susceptibility is sought. Next, we form a multicomponent
objective function that has the flexibility to generate different
types of models. The form of this objective function is such that
it can correct for the undesirable aspects of the mathematically
acceptable model in Figure 3, namely-the concentration of
susceptibility near the surface, the excessive structure, and the
existence of negative susceptibilities. Our objective function
incorporates an optional reference model so that the con-
structed model is close to that. It penalizes roughness in three
spatial directions, and it has a depth weighting designed to
distribute the susceptibility with depth. Additional 3-D weight-
ing functions in the objective function can be used to incorpo-
rate further information about the model. Such information
might be available from other geophysical surveys, geological
data, or the interpreter’s qualitative or quantitative under-
standing of the geologic structure and its relation to the
magnetic susceptibility. These 3-D weighting functions can also
be used to answer questions about the existence of suscepti-
bility features found from previous inversions. Negative sus-
ceptibilities are prevented by making a transformation of

FIG. 2. The total field anomaly produced by the slab model in

FIG. 1. Slices through a 3-D magnetic susceptibility model Figure 1. The inducing field has direction I = 75° and D = 25°
composed of a dipping slab in a nonsusceptible half-space. The and a strength of 50 000 nT. Uncorrelated Gaussian noise, with
slab is buried at a depth of 50 m and extends to 400-m depth a standard deviation of 2% of the datum magnitude plus 1 nT,
at a dip angle of 45° . The gray scale indicates the value of is added to the data. The gray scale indicates the magnetic
magnetic susceptibility in SI units. anomaly in nT.
396 Li and Oldenburg

variables and solving a nonlinear inverse problem. The numer- by inverting a field data set over a copper-gold porphyry
ical solution for the inversion is accomplished by dividing the deposit and a subsequent discussion.
earth into a large number of cells so that relatively complex
geologic bodies can be constructed. The computational diffi- INVERSION METHODOLOGY
culties often encountered in solving large matrix systems are
avoided by working explicitly with a generalized subspace
Each magnetic anomaly datum observed above the surface
algorithm.
can be evaluated by calculating the projection of the anoma-
The paper begins by outlining our inversion methodology
lous magnetic field onto a given direction. Let the source
and empirically estimating parameters for the depth weighting
region be divided into a set of rectangular cells by an orthog-
based upon synthetic inversion of single 3-D prisms. Data from
onal 3-D mesh and assume a constant magnetic susceptibility
two synthetic models are then inverted. The paper concludes
value within each cell. Further we assume that there is no
remanent magnetization and that the demagnetization effect is
negligible. Thus only the induced magnetization is considered.
This magnetization is uniform within each cell and is given by
the product of the susceptibility and the inducing geomagnetic
field The magnetic anomaly at a location on, or above, the
surface is related to the subsurface susceptibility by a linear
relationship

(1)
where d = is the data vector and
is the susceptibility in the cells. The matrix has as
elements which quantify the contribution of a unit suscep-
tibility in thejth cell to the ith datum. Closed form solutions for
were first presented in Bhattacharyya (1964) and later
simplified in Rao and Babu (1991) into a form more suitable
for fast computer implementation. The function is the
projection onto a given direction of the magnetic field that is
produced by a rectangular cell, so equation (1) is valid for
computing different magnetic anomalies. For example, a pro-
jection onto the vertical direction gives the vertical magnetic
anomaly while a projection onto the ambient geomagnetic field
direction yields the total magnetic anomaly. Thus, the method
presented here can be used to invert different types of mag-
netic data and in the following, we simply refer to them as the
magnetic data with the understanding that it is direction
specific.
Our inverse problem is formulated as an optimization
problem where an objective function of the model is minimized
subject to the constraints in equation (1). For magnetic
inversion, the first question that arises concerns definition of
the “model.” Two possible choices are and but any
function can, in principle, be used. In general, we prefer
to invert for since the field anomaly is directly proportional
to the susceptibility that varies on a linear scale. But depending
upon the expected dynamic range of susceptibility and the
physical interpretation attached to its value or variation, it may
be that is more desirable. To accommodate this, we
introduce the generic symbol m for the model with the
understanding that it might be or any monotonic
function Having defined a model, we next construct an
objective function, which when minimized, produces a model
that is geophysically interpretable. The details of the objective
function are problem dependent, but generally we need the
FIG. 3. The susceptibility model constructed by minimizing flexibility to be close to a reference model m 0 and also require
subject to fitting the data in Figure 2. As a mathematical that the model be relatively smooth in three spatial directions.
solution, this model provides little, if any, information about Here we adopt a right-handed Cartesian coordinate system
the subsurface susceptibility distribution. It effectively illus-
trates the nonuniqueness inherent to the inversion of static with x positive north and z positive down. Let the model
magnetic field data. objective function be
 3-D Inversion of Magnetic Data 397

extra information is incorporated, the inversion derives a

model that not only fits the data, but more importantly, also
has a likelihood of representing the earth. From the viewpoint
of magnetic inversion, such an approach allows one to con-
struct a most-likely earth model that uses all available infor-
mation, and it can also be used to explore the nonuniqueness.
These two aspects form the foundation of a responsible
interpretation.
The kernels (values of for the surface magnetic data
decay with depth. It is for this reason that an inversion that
minimizes dv subject to fitting the
data will generate a susceptibility that is concentrated near the
surface. To counteract the geometric decay of the kernels and
(2)
to distribute susceptibility with depth, we introduce a weighting
of the form into and optionally
where functions ws, wx, wy , and wz are spatially dependent include it in The values of and z 0 are investigated in the
weighting functions while and are coefficients following section, but their choice essentially allows equal
that affect the relative importance of different components in chance for cells at different depths to be nonzero.
the objective function. Here, w ( z ) is a depth weighting The next step in setting up the inversion is to define a misfit
function. It is convenient to write equation (2) as measure. Here we use the 2-norm measure
refers to the first term in equation (2)
and refers collectively to the remaining three terms that (3)
involve variation of the model in three spatial directions. and we assume that the contaminating noise on the data is
The objective function in equation (2) has the flexibility of independent and Gaussian with zero mean. Specifying to
constructing many different models. The reference model m 0 be a diagonal matrix whose ith element is where is the
may be a general background model that is estimated from standard deviation of the ith datum, makes a chi-squared
previous investigations, or it could be the zero model. The variable distributed with N degrees of freedom. Accordingly
reference model would generally be included in but can be
= N provides a target misfit for the inversion.
removed if desired from any of the remaining terms. Often we The inverse problem is solved by finding a model m that
are more confident in specifying the value of the model at a minimizes and misfits the data by a predetermined amount.
particular point than in supplying an estimate of the gradient. This is accomplished by minimizing
The relative closeness of the final model to the reference where is our target misfit and is a Lagrangian
model at any location is controlled by the function w,. For multiplier. To perform a numerical solution, we first discretize
example, if the interpreter has high confidence in the reference the objective function in equation (2) using a finite-difference
model at a particular region, he or she can specify ws to have approximation according to the mesh defining the susceptibil-
increased amplitude there compared to other regions of the ity model. This yields

(4)

model. The weighting functions wx, wy , and wz can be designed where m and m0 are M- length vectors. The individual matrices
to enhance or attenuate structures in various regions in the are calculated straightforwardly once the
model domain. If geology suggests a rapid transition zone in model mesh and the weighting functions ws, wx, wy, wz, and
the model, then a decreased penalty for variation can be put w ( z ) are defined (see Appendix). The cumulative matrix
there, and the constructed model will exhibit higher gradients is then formed. For our formulation, the matrix is
provided that this feature does not contradict the data. There- never computed explicitly but we shall use it to derive our final
fore, the reference model and four 3-D weighting functions equations.
allow for the incorporation into the inversion of additional The inverse problem is solved by minimizing with an
information other than the magnetic data. The additional appropriate minimization technique. To reduce computation
information can be from previous knowledge about the sus- and to invoke positivity, we use a subspace methodology. In its
ceptibility, from other geophysical surveys, or from the inter- general form, the subspace technique allows the model param-
preter’s qualitative or quantitative understanding about the eter to be both positive and negative, and thus to ensure
geologic structure and its relation to susceptibility. When this positive susceptibility, we may need to invoke a transformation
398 Li and Oldenburg

of variables. Whether or not the transformation is required of that achieves the smallest misfit is taken. The search is
depends upon the relationship between mi and If usually accomplished by solving equation (10) a number of
so that interpretations are carried out in the logarithmic times using different values. Once the optimum value of is
domain, then no further transformation is necessary since found, the system is solved again to obtain the coefficients
will be positive irrespective of the sign of m i . However, if and the model perturbation. This iterative process is continued
and is a positive function, then a until the final expected data misfit is achieved and the model
transformation is required. All possibilities can be handled by objective function undergoes no significant decrease with
introducing a new parameter p, such that where successive iterations. Subspace vectors v i are generated mainly
f ( p ) is a monotonic function whose inverse and first-order from the gradients of the data and model objective functions.
derivative exist. This mapping is then incorporated directly into The data are grouped to form subobjective functions of misfit,
the subspace minimization process. and a steepest descent vector corresponding to each subobjec-
Let p ( n ) denote the parameter vector at the nth iteration and tive function is used as a subspace vector. Partitioning of the
denote the sought perturbation. Performing a Taylor ex- data can be formed by grouping data that are spatially close, or
pansion of the perturbed model objective function about the by grouping data such that each group has approximately the
point p ( n ) yields same contribution to the total data misfit. Both approaches
have worked well. The model objective function is partitioned
(5) and the gradient vector associated with each of the four
components in the model objective function provides addi-
where is a diagonal matrix with elements tional subspace vectors. In addition, a constant vector is always
included, and the selected subspace vectors are orthonormal-
(6) ized before being used in the search. More details on the
implementation of the subspace method for the linear inverse
problem can be found in Oldenburg and Li (1994).
A similar Taylor expansion applied to the misfit objective
The final item of practical importance is the specification of
functional yields
the mapping needed to ensure positivity of susceptibility. The
(7) positivity is required since we are dealing only with induced
magnetization, and the presence of negative susceptibility is
At each iteration we desire a perturbation that minimizes negligible in practical geophysical applications. Although our
equation (4) subject to generating a data misfit of formalism permits the minimization of m = the two most
where is the target misfit at the nth iteration. In the common situations are m = ln and m = When m =
subspace technique we represent the perturbation as ln we set p = m and hence the matrix in equation (10)
is the identity matrix. If m = we use the two-stage mapping
proposed in Oldenburg and Li (1994). It is composed of an
(8) exponential segment and a straight line. The two segments are
joined together such that the mapping and its first derivative
are both continuous. The mapping is given by
where the M-length vectors v i ( i = 1, q ) are as yet arbitrary.
Writing the objective function to be minimized in terms of the
coefficients yields
(11)

where p = p 1 is the transition point between exponential and

linear segments, and is selected to be small enough such
(9) that susceptibilities smaller than are not significantly differ-
Differentiating with respect to the coefficients a yields the final ent from zero when the final interpretation is carried out. Here,
equations and hence p l are chosen so that the ratio
does not exceed about two orders of magnitude. This prevents
the elements Fii from becoming too disparate. We note that
the ith row of is multiplied by Fii, and if this value is too
small, the ith row of is essentially annihilated and there will
be no possibility of adjusting the value of the ith cell. However,
if the ratio is too small, the flexibility in the mapping will be
restricted and this affects the convergence rate of the
algorithm. In the limit that the nonlinear mapping
(10) degenerates into a linear truncation and the inversion will not
converge. However, between the above two extremes, there is
We note that the matrix is q x q and therefore the system of a wide range of values for the ratio that can yield a good
equations is easily solved if q is small. At each iteration, we mapping. Based upon numerical experiments (Oldenburg and
search for a value of that yields the target misfit for that Li, 1994), we have chosen a value of 50.0 for this ratio for the
iteration. If the target misfit cannot be reached, then the value examples throughout this paper.
 3-D Inversion of Magnetic Data 399
DEPTH WEIGHTING

It is well known that static magnetic data have no inherent The decay factor e-az causes the constructed model m,(z) to
depth resolution. For instance, when minimizing have structure concentrating toward the region of small z in the
structures tend to concentrate near the surface classic model construction that minimizes since the
regardless of the true depth of the causative bodies. In terms of model will be a linear combination of the kernels, i.e.,
model construction, this is a direct manifestation of the nature
of the kernels whose amplitudes rapidly diminish with depth.
The tendency to put structure at the surface can be overcome
(13)
by introducing a depth weighting to counteract this natural
decay. Intuitively, a weighting that approximately compensates
for the decay gives cells at different depths equal probability to This is shown in Figure 4a and 4b for two different models.
enter into the solution with a nonzero susceptibility. Before These models are constructed from five data ( i = 0,4) to which
proceeding with the details of the weighting function for noise has been added. It is apparent that the constructed
magnetic inversion, we illustrate the necessity, and effective- model is shifted toward small z where the amplitude of kernels
ness, of such a weighting function using a simple 1-D problem. is relatively large. One way to counteract the bias is to seek a
Consider a set of data d = ( d 1, , . . . , d N ) T generated from the solution in model space that is spanned by the nondecaying
equation portion of the kernels, in this case just the cosine functions.
The desired model would have the form

(12)
(14)
where the kernels are

FIG. 4. A 1-D example showing the use of a weighting function in the inversion procedures to counteract the natural decay in the
kernel function. In all panels the dashed line shows the true model. Panels (a) and (b) show, for the two different true models,
respectively, the model constructed using the original kernel functions with the decaying factor e-az. Notice the shift of the
recovered model towards the small z region. Panels (c) and (d) show the weighted models recovered by applying a weighting function
w ( z ) = e-az/2. They are better representations of the true model.
400 Li and Oldenburg

where are coefficients. Free from the influence from the in both directions, and 2% Gaussian noise is then added. The
decay factor, a model constructed from this set of basic observation is assumed to be 1 m above the surface and the
functions should have a better chance of having significantly inducing field has I = 75°, D = 25°. The region directly beneath
high values at depth. the data grid is taken as the model domain and discretized into
We accomplish this by finding an appropriate weighting 4000 cells (20 cells in each horizontal direction and 10 along
function w ( z ). We first rewrite the data equation as depth) of 50 m on a side.
Given the stated data parameters and model discretization,
the estimated value of z 0 in the depth weighting function is
( 1 5 ) 25 m. Figure 5 shows the comparison of the kernel beneath a
datum point and the function w 2( z ). This weighting function is
where are the weighted kernels and m w ( z ) is the used to invert surface data caused by the susceptible prism, and
weighted model. Then the inverse problem is solved by mini- the results of minimizing are shown in Figure 6. Each
mizing and the solution is given by panel in the figure is the cross-section through the center of the
model obtained by inverting the data set produced by a cube at
a different depth. They are rather good recoveries in terms of
(16) source depth, which is indicated by the superimposed outline
of the true body in each section.
In the above analysis we have established a practical way for
Dividing by the weighting function and substituting in
estimating an appropriate depth weighting function that dis-
yields
tributes the susceptibility more uniformly with depth. The
weighting is valid when the model objective function consists
only of In general, we like to include a penalty against
roughness and thereby produce a model that is smooth. To
(17) incorporate the above weighting scheme in the spatial varia-
tions, we make the following argument. Since minimizing
This equation can be made identical to equation (14) by
tends to provide a reasonable depth distribution, we wish only
choosing Carrying out the weighted inversion
to improve the model’s smoothness while maintaining the
for the above two data sets produces models shown in
depth characteristic. A conceptually consistent approach
Figures 4c and 4d. They are much better representations of
would be to apply the roughness measures to the weighted
true models.
model. We form a generic model objective function
This methodology is then applied to the inversion of surface
magnetic data by finding the appropriate weighting function
that counteracts the depth decay of the data kernels. There is
no distinct separable factor defining the decay in the kernel,
therefore we resort to an empirical estimate. Since the decay
rate depends upon the observation height as well as the size
and aspect ratios of the cells making up the 3-D model, such
estimates are expected to be problem dependent. Numerical
experiments indicate that the function of the form ( z + z 0)-3
closely approximates the kernel’s decay directly under the
observation point, given a correctly chosen value of z 0. This is
consistent with the fact that, to first order, a cubic-shaped cell
acts like a dipole source whose magnetic field decays by inverse
distance cubed. The value of z 0 can be obtained by matching
-3
the function ( z + z 0) with the kernel function beneath the
observation point. Thus, a reasonable candidate for the depth
weighting function is given by

(18)

The susceptibility model constructed by minimizing a model

objective function consisting of only i.e.,

(19)

subject to fitting the data should place the recovered anomaly FIG. 5. Comparison of the kernel function (solid) directly
at approximately the depth of the causative body. This hypoth- beneath the observation point with the estimated curve
(dashed) given by w 2( z ) = ( z + z 0)-3 with z 0 = 25 m. The
esis is tested by inverting surface data produced by a suscep- source cell is a cube of 50 m on a side. Here, z denotes the
tible cubic body at three different, depths. The cube is 200 m on depth to the center of the cell. Both curves are normalized for
a side. Data are calculated over a 21 X 21 grid of 50-m spacing comparison.
 3-D Inversion of Magnetic Data 401

it is observed that straightforward inclusion of the depth weight-

(20) ing derived above into the 3-D weighting function in the form of
can yield reasonable results.

where the depth weighting is applied inside the derivatives of PRACTICAL ASPECTS OF DATA PREPARATION
the roughness components and the reference model m 0 can be
removed from any term if desired. This type of depth weighting The data used in the inversion are the residual data obtained
has proven to work satisfactorily on a number of synthetic by subtracting a regional field from the initial observation. The
examples and is the default choice in our algorithm. The inversion algorithm has been developed under the assumptions
examples to be presented in the following sections all use this that the surface magnetic anomaly is produced by the induced
depth weighting function. magnetization only and that there are no remanent magneti-
Before proceeding further, we remark that the above zation or demagnetization effects present. Incorrect removal of
weighting represents only one possibility. One could poten- regional field, or any deviation from the above assumptions, is
tially design a different weighting by incorporating the depth expected to cause a deterioration in the inversion results.
weighting in the usual 3-D weighting functions ws, wx, wy, wz. Furthermore, the susceptibility distribution is mathematically
Such an approach applies the depth weighting outside the deriv- represented by a piece-wise constant function defined on a
ative operators directly. However, the decay rate of the depth user-specified grid of cells. Magnetic sources, however, have a
weighting for each component will be different, and it is difficult to wide range of physical sizes. In some cases, source dimensions
establish a consistent rule for the choice of the different weight- will be significantly smaller than the size of cells in the
ings. In addition, the extra set of parameters required by such a mathematical model. If measurements are taken close to such
weighting scheme introduces more subjectivity into the inversion a source, the resulting anomaly will have a width that is
process. We have not explored this approach in detail; however, significantly smaller than that produced by a single cell in the
mathematical model and this may produce artifacts. We ame-
liorate this problem by inverting data that have been upward
continued to a height approximately equal to the width of the
surface cells in the model. We arrive at this conclusion from a
numerical experiment. We first generate the magnetic field
from a small localized surface source that is assumed to be a
cube of width At each height h above the surface, a
one-parameter inverse problem is carried out to find a uniform
susceptibility of a large surface cube that has a width of L and
shares a common horizontal center with the small cube. If HL
is the field of the large cell that best reproduces then the
misfit functional,

(21)

FIG. 6. Cross-sections through the center of the recovered

model for a cube at a central depth of 150,200, and 250 m. The
cube is 200 m on a side. The inversion uses the weighting FIG. 7. The misfit between magnetic field as a result of a small
function derived from the kernel decay estimated in Figure 5. cubic source and the field as a result of a larger cubic model
The true position of the cube is outlined in each cross-section. cell having a best fitting susceptibility. The numbers indicate
As the true source depth increases and, as a result, the the ratio of the cell width. The misfit is plotted as a function of
high-frequency content in the data decreases, the recovered the observation height normalized by the width of the model
model becomes increasingly smooth and attains a smaller cell. Note that the misfit decreases rapidly until the height is
amplitude. However, the depth of the recovered model is close approximately equal to the width of the model cell, and that it
to the true value. changes slowly thereafter.
402 Li and Oldenburg

can be computed, where is the surface area of the data map. The inversion uses 54 subspace vectors and achieves the
Figure 7 shows the misfit function r ( h ) for trial values of = expected misfit in 13 iterations. The recovered model is shown
0.1, 0.2, 0.4. We note that r ( h ) decreases rapidly until h L, in Figure 9. It is smoother, has a slightly lower amplitude than
and that it changes slowly thereafter. Since the above misfit the model in Figure 8, and it recovers the essential features of
analysis is a worst case scenario because the contaminating the true model such as the depth and dip angle.
body is located at the surface, the suggestion of upward It is observed, in this example and in other synthetic and
continuing the data to a height approximately equal to the field test examples, that minimizing either the first term in the
width of surface cells may be somewhat conservative, and model objective function in equation (20), or using all
inversionists may want to vary this. However, in many field four terms, generates models that are reasonable representa-
surveys, magnetically susceptible small bodies exist close to the tions of the true structure. In the absence of prior information,
surface and hence upward continuing the data prior to inver- both models can provide useful information about the subsur-
sion is prudent. face susceptibility distribution. However, the model minimiz-
ing can be obtained at less computational cost. Further-
SYNTHETIC EXAMPLES

As the first example, we invert the total field anomaly data

given in the Introduction. The model consists of a 3-D dipping
slab buried in a nonsusceptible half-space (slab model).
Figure 1 shows three slices through the slab model. The
susceptibility of the slab is 0.06 (SI unit). Under an inducing
field with a strength of 50 000 nT and a direction at I = 75° and
D = 25°, the slab model produces the surface total magnetic
anomaly shown in Figure 2, which consists of 441 data over a
21 X 21 grid of 50-m spacing. The data have independent
Gaussian noise added whose standard deviation is equal to 2%
of the accurate datum magnitude plus 1 nT. We invert these
441 noise-contaminated data to recover the susceptibility of an
earth model parametrized by 4000 cells of 50 m on a side (20
cells in each horizontal direction and 10 in depth).
The data are partitioned into 49 groups to provide 49 search
vectors for the subspace algorithm. In addition, each compo-
nent in the model objective function provides one basis vector,
and a constant vector is included. For the depth weighting, the
value of z 0 is estimated as 25 m. The additional 3-D weightings
in the objective function are all set to unity. The reference
susceptibility model is set to zero. For the nonlinear mapping,
we choose = 0.0002 and = 0.01.
First, we invert the data by minimizing an objective function
composed only of the and using m = as the model
parameter. A total of 51 subspace vectors are used at each
iteration. The inversion reaches the expected misfit in 13
iterations but a few extra iterations are performed in an
attempt to further reduce the value of the model objective
function while keeping the misfit at the target value. By
iteration 18, the objective function is decreasing by less than
1% per iteration, and the process is terminated. The con-
structed susceptibility model is shown in Figure 8 and can be
compared with the true model in Figure 1. The tabular shape
of the anomaly and its dipping structure are clear, and the
depth extent is reasonably recovered. The amplitude of the
recovered model is slightly higher than the true value, but the
dip angle inferred from the recovered model is close to the true
value. We point out that the model sections should be plotted
using gray shading for each cell to reflect the piece-wise
constant nature of the model. However, when the model has
only a small number of cells in each spatial direction, the
structural trends are more readily shown when contours are
used. For this reason, we have contoured the model sections. FIG . 8. Model obtained from inverting the data shown in
Figure 2 by minimizing only which has the depth weight-
Next, the same data are inverted using a model objective ing applied. This is to be compared with the true model in
function that includes penalty terms on spatial roughness, Figure 1. The major features in the true model, such as dip
The depth weighting is applied to all terms, as in equation (20). angle and depth exte nt, are evident in the recovered model.
 3-D Inversion of Magnetic Data 403
more, the depth weighting in this case is rather well supported as an initial model. The available prior information can be
by mathematical analysis whereas it is an argued extension for incorporated into the second inversion by forming a reference
the three roughness components. Therefore, a reasonable model and 3-D weighting functions, ws, wx, wy, wz.
approach to inverting field data might be a two-step process. We now invert the same data by using m = ln as the
The data can be inverted first by minimizing and the model. It is not possible to incorporate a zero susceptibility as
resultant model may be used in the interpretation as a prelim- the reference model, so we minimize an objective function
inary result. If there are interesting features present and if one consisting of with the reference model removed. The same
desires to refine the model by incorporating prior information depth weighting is applied to all terms of Since = em,
to enhance or attenuate the structural complexity in different the positivity of the susceptibility is ensured without invoking
regions, a second inversion can be carried out using an the transformation of variables. The result is shown in
objective function consisting of both and The model Figure 10a. This is a cross-section at x = 500 m and plotted on
obtained by minimizing can then be used in this inversion a logarithmic scale in accordance with the model used in the
inversion. The inverted susceptibility shows the presence of the
dipping anomaly as a broad region of high susceptibility.
However, the interpretation based upon such a model can be
complicated by the variations of susceptibility that are small
and have little effect on the surface data. We have replotted
the cross-section on a linear scale in Figure 10b and the
anomalous region is now delineated more clearly. Its top
portion indicates the tabular body and defines the depth to the

FIG. 10. The model obtained from inverting the data shown in
Figure 2 by using m = ln as the model and minimizing
FIG. 9. The model derived from inverting the slab model data with the reference model removed. The inverted logarithmic
in Figure 2 by minimizing the model objective function having susceptibility in cross-section at x = 500 m is shown in (a) and
both and The same depth weighting is used. This it is replotted on a linear scale in (b). As a comparison, the
model appears to be smoother and has a smaller amplitude result obtained by using m = and the same objective
than that in Figure 8. function is shown in (c).
404 Li and Oldenburg
top and dipping angle. The anomaly terminates at a shallower scale as demonstrated here, and since the magnetic data are
depth than the true model and has a nearly horizontal exten- linearly related to the susceptibility, we generally prefer to work
sion to the left. As an exact comparison, Figure 10c is the with the susceptibility K as the model in the inversion.
susceptibility model obtained by minimizing but using As the second example we invert the total field anomaly data .
m = as model and invoking the positivity. This is a smoother produced by a slightly more complicated model and with two
model and exhibits more gradual changes in the susceptibility. different inducing field directions. The true model is shown in
It has a slightly deeper extent than the model in Figure 10b. With Figure 11 in the same format as before. It is a dipping slab
the exception of details toward the bottom, however, both models having its top and bottom portions offset to simulate the result
provide almost the same information about the anomalous sus- of a normal faulting. The faulted slab strikes north. The data
ceptibility region. It might be concluded that inversion using from this model, when the inducing field has a direction of I =
either linear or logarithmic susceptibility is viable for practical 45° and D = 45°, are shown in Figure 12. Again Gaussian noise
applications. However, we note that the presentation in has been added to the data. The inversion minimizes an
Figure 10b is inconsistent with the model used in the inversion. objective function consisting of and that have the
Since the inverted susceptibility is easier to interpret on a linear same depth weighting and nonlinear mapping as used to
produce the results in Figure 9. Figure 13 displays the recov-
ered model in three slices. It shows two distinct anomalous
regions of susceptibility that correspond to those in the true
model. The dipping structure is evident from the top block. On
plan view, the strike direction and the strike length of the
anomaly are also well recovered.
When the inducing field direction is I = 0° and D = 45°, the
surface anomaly with added Gaussian noise is that shown in
Figure 14. Carrying out the inversion using an identical model
objective function generates the model shown in Figure 15. It
is similar to the model shown in Figure 13, which is recovered
under an inducing field at 45° inclination. Again, the two
separate blocks, the dipping direction, and the length and
direction of the strike, are all reasonably recovered. This is a
positive result in that, although the surface anomalies have
very different expressions under different inducing field direc-
tions, the inversion algorithm is able to consistently recover the
source structure. Moreover, the algorithm had no difficulty in
inverting data generated from an inducing field having zero
inclination; such data often pose problems in interpretations
that include a reduction to pole.
We emphasize that positivity has played a pivotal role in all
the inversions. Magnetic data generally have regions of nega-

FIG . 11. The second synthetic test example. The top and FIG . 12. The surface total field anomaly produced by the
bottom portions of the anomalous susceptibility are offset to faulted slab in Figure 11, under an inducing field at I = 45° and
simulate a norm fault structure. It also has a large strike length D. = 45°. Uncorrelated Gaussian noise is again added to the
in the north direction. data.
 3-D Inversion of Magnetic Data 405
tive values that result from dipping bodies or inclined inducing FIELD EXAMPLE
field, or both. Without positivity, the constructed susceptibility
is often negative and the dipping bodies appear more vertical. As the final example, we invert field data taken over a
Recovery of correct dip and, to some extent, depth to the top copper-gold porphyry deposit at Mt. Milligan in central British
of the anomalous body, are often the result of invoking Columbia. The host rocks for the deposit are early Mesozoic
positivity. Once the positivity is imposed, it is no longer true volcanic and sedimentary rocks and contain intrusive monzo-
that an equivalent stratum that reproduces the data exists at nitic rocks that have accessory magnetite. Porphyry-style alter-
any depth. Therefore, cells of anomalous susceptibility cannot ation and copper-gold mineralization are contemporaneous
be placed arbitrarily close to the surface, and no equivalent with the intrusive events. The copper and gold are known to be
source can be constructed with negative susceptibilities. This concentrated in the potassic alteration assemblage, which is
restricts the class of admissible models and, consequently, mainly around the contact of the monzonite intrusions and
reduces the nonuniqueness. may extend outward and into fractured volcanic rocks. Among
other minerals, magnetite is one of the strong indicators of the
potassic alteration. Ground magnetic data are acquired in the
region at 12.5-m spacing along lines in the east direction and
spaced 50 m apart. Our study of the data set has focused on a
1.2 km x 1 km area, which covers a large monzonite body
known as the MBX stock and contains a reasonably isolated set
of magnetic anomalies. Fairly detailed information about the
geology is available through a major drilling program, but no
susceptibility logs were available.
Magnetic data from a larger area were first upward contin-
ued to 20 m. A regional field was then defined and removed
from the upward continued data. The continuation operation
suppresses the noise in the data and also facilitates the
discretization of the topographic surface for the model so that
all observation points remain above the discretized surface.
Although the original data were collected at 12.5-m spacing,
we use the data at 25-m spacing. This yields 1029 data points at
varying elevations. Figure 16 shows the data contoured accord-
ing to their horizontal locations. The direction of the inducing
field is I = 75° and D = 25.73°. Several major magnetic highs
are observed in the map. However, the influence of anomalies
adjacent to the map is also visible along the edges. We choose
a model domain that is horizontally larger than the data area,
coincides at the top with the highest point on the topographic
surface, and extends to 450-m depth. The model is discretized
horizontally at a 25-m interval beneath the area of data. In the

FIG. 13. The susceptibility model recovered from the data

shown in Figure 12. It is seen that both the top and bottom FIG. 14. The surface total field anomaly produced by the
block of the true model are recovered and the strike direction faulted slab in Figure 11 under an inducing field at I = 0° and
and length are well defined. D = 45°. Uncorrelated Gaussian noise is added to the data.
406 Li and Oldenburg

vertical direction, the first 100 m is divided at a 12.5-m interval hundred subspace vectors generated by dividing the data map
so that the surface can be adequately discretized onto the into small subareas are used in the inversion. We use a
model mesh. Below the depth of 100 m, an interval of 25 m is nonlinear mapping with = 0.0002 and = 0.02. The
used. This results in a mesh with 52 x 44 x 22 cells. Once the recovered model is shown in Figure 17 as one plan-section and
mesh is defined, the topography is discretized onto it. The three cross-sections. From the plan-section, two concentrated
43 428 cells below this surface define the susceptibility model, susceptibility highs are observed in the central region. Sur-
and the inverse problem is therefore formalized by inverting rounding them are three linear anomalies trending northeast.
1029 data to recover the susceptibilities in those cells. The In the cross-sections, the major anomalies are seen at moder-
depth weighting is referenced to the top of the model domain. ate depths but there is considerable variation in the depth to
Each datum is assumed to have an error whose standard the top. There are also smaller anomalies extending to the
deviation is equal to 5% of its magnitude plus 10 nT. The error surface. In general, there are more detailed structures near the
estimate includes not only the repeatability of the instrument surface and the model becomes increasingly smooth at greater
reading but also the geological noise and errors introduced by depths. As required by the objective function, there is no
the inaccurate recording position and by separating the anom- excessive structure associated with each unit of high suscepti-
alous field from the initial total field measurements. One bility region. Comparison with drill logs indicates that the
recovered magnetic susceptibility highs are mostly associated
with the monzonite intrusions and with faults or fracture zones.
Figure 18 compares the recovered susceptibility model with the
geology (Cam DeLong, personal communication) in the cross-
section at x = 600 m. The large susceptibility high is spatially
well-correlated with the MBX stock and reflects the initial
magnetite content in the intrusion. Two smaller susceptibility
highs are present east of the stock. The high at y = 650 m
coincides with the boundary of stock and porous trachytic units
while the high at y = 900 m coincides with the upper portion
of the Rainbow dyke. These are locations of the most intensive
potassic alterations and the susceptibility highs are indicative.
of the magnetite produced by the alteration process. Over all,
this is a rather encouraging result.

CONCLUSION

We have developed an algorithm to invert surface magnetic

data for general 3-D susceptibility distributions. Although we
have illustrated the algorithm using examples on the scale
pertinent to mining applications, the method is general and
applicable to problems on different scales ranging from envi-
ronmental to regional investigations. To overcome the inher-
ent nonuniqueness, we obtain the solution by minimizing a

FIG. 16. The extracted total field anomaly from ground mag-
FIG. 15. The susceptibility model recovered from the data netic data at Mt. Milligan Copper-gold porphyry deposit. The
shown in Figure 14. This model is similar to that shown in data are contoured according to their horizontal locations in
Figure 13. this map, although thev are at different elevations.
 3-D Inversion of Magnetic Data 407
specific objective function of the model. Our model objective
function has the ability to incorporate prior information into the
inversion via a reference model and 3-D weighting functions. A
crucial feature of the objective function is a depth weighting
function that counteracts the natural decay of the kernel func-
tions. The parameters of the depth weighting depend upon the
discretization of the model but are easily calculated. The minimi-
zation is carried out using a subspace technique that reduces the
computational effort and allows the positivity constraint of sus-
ceptibility to be incorporated. Both susceptibility and logarithmic
susceptibility can potentially be used as the model in the inver-
sion. Since the data are linearly related to susceptibility, and since
usually absolute values of susceptibility are required for interpre-
tation rather than relative values, especially in regions of very low
FIG. 18. Comparison of the recovered susceptibility model in a
susceptibility, we have generally chosen to work with susceptibil- cross-section ( x = 600) with the geology for the Mt. Milligan
ity. To suppress the noise from small magnetic bodies near the deposit. The susceptibility high within the MBX stock reflects
surface, we recommend in general that the data be upward the initial magnetite in the intrusive while the susceptibility
continued to a height comparable with the width of the surface highs near the Rainbow dyke are related to the magnetite
produced by potassic alteration.
cell before inversion.

Applications of our inversion to synthetic data sets have

produced models representative of the true structures and
demonstrated the ability of the algorithm to construct consis-
tent models at different magnetic latitudes. Inversion of field
data has produced a susceptibility model that is consistent with
the known geology and mineralization information. These
results represent an encouraging conclusion: although the
inversion of magnetic data seems impossibly nonunique when
one has a large number of cells, inversions using a properly
designed model objective function can produce susceptibility
distributions that yield meaningful geologic information.

ACKNOWLEDGMENTS

This work was supported by an NSERC IOR grant and an

industry consortium “Joint and Cooperative Inversion of Geo-
physical and Geological Data.” Participating companies are
Placer Dome, BHP Minerals, Noranda Exploration, Cominco
Exploration, Falconbridge, INCO Exploration & Technical
Services, H u d s o n B a y E x p l o r a t i o n a n d D e v e l o p m e n t ,
Kennecott Exploration Company, Newmont Gold Company,
Western Mining Corporation, and CRA Exploration Pty.

REFERENCES

Bhattacharyya, B. K., 1964, Magnetic anomalies due to prism-shaped

bodies with arbitrary magnetization: Geophysics, 29, 5 17-53 1.
1980, A generalized multibody model for inversion of magnetic
anomalies: Geophysics, 45, 255-270.
Green, W. R., 1975, Inversion of gravity profiles by use of a Backus-
Gilbert approach: Geophysics, 40, 763-772.
Guillen, A., and Menichetti, V., 1984, Gravity and magnetic inversion
with minimization of a specific functional: Geophysics, 49, 1354-
1360.
Last, B. J., and Kubik, K., 1983, Compact gravity inversion: Geophys-
ics, 48, 713-721.
Oldenburg, D. W., and Li, Y., 1994, Subspace linear inversion meth-
ods, Inverse Problems, 10, 915-935.
Rao, D. B., and Babu, N. R., 1991, A rapid methods for three-
dimensional modeling of magnetic anomalies: Geophysics, 56,
1729-1737.
FIG. 17. The recovered susceptibility model shown in one Wang, X., and Hansen, R. O., 1990, Inversion for magnetic anomalies
plan-section and three cross-sections. The plan-section is at the of arbitrary three-dimensional bodies: Geophysics, 55, 1321-1326.
depth of 150 m and the three cross-sections are at x = 600,500, Zeyen, H., and Pous, J., 1991, A new 3-D inversion algorithm for
and 400 m, respectively. magnetic total field anomalies: Geophys. J. Int., 104, 583-591.
408 Li and Oldenburg

APPENDIX
MODEL OBJECTIVE FUNCTION

Our inversion method uses a model objective function of the Each component matrix can be written as the product of three
form individual matrices and one oefficient. That is,

(A-3)

where are diagonal matrices representing the spatially

dependent 3-D weighting functions, are the finite-difference
operators for each component, and is a diagonal matrix
representing the discretized form of depth weighting function

The elements of are given by They are defined

(A-1) over each cell for and over each interface between
adjacent cells in the respective directions for and
The numerical evaluation of this functional is carried out by has elements on its diagonal, where Ax, Ay,
introducing the model mesh and evaluating all terms using a and are the cell width. The matrix has two elements
finite-difference approximation. The discretized model objec- in each row, where is the distance between
tive function has the form the centers of cells adjacent in the x - direction. Similarly,
and have elements and respec-
tively, where and are the distances between centers of
adjacent cells in the y- and z - directions. Once the mesh is
defined and all weighting functions, ws, wx, wy , wz, and w ( z )
are chosen, equation (A-3) is evaluated straightforwardly
(A-2) and is formed.