Geophys. J . Int.

(1994) 119, 857-868

Asymptotic theory for diffusive electromagnetic imaging

Jean Virieux,' Carlos Flores-Luna1,2and Dominique Gibert3
' Institut de Ge'odynamique, UNSA Rue Albert Einstein, 06560 Valbonne, France
'CICESE, Depto. Geojisica Aplicada KM 107 Carretera. Tijuana-Ensenada, Ensenada 22800, Baja California, Mexico
Geosciences Rennes, Rennes I Avenue GCnPral Leclerc, 35042 Rennes cidex, France

Accepted 1994 June 6. Received 1994 June 6; in original form 1993 August 4

Key words: electrical conductivity, electromagnetic surveys, inversion.

INTRODUCTION
Recovering conductivity information in the Earth is a
difficult task because the electromagnetic response is mainly
controlled by a diffusion process broadening the signal.
Many attempts have encountered limited success (Weidelt
1972, 1975; Barnett 1984; 'IXpp, Hohmann & Swift 1984;
Smith 1988; Tarits 1989) associated with the difficult task of
locating the conductivity anomaly responsible for the
electromagnetic signal observed at the free surface. In many
papers (Zhdanov & Frenkel 1983; Zhdanov 1988; Zhdanov
& Booker 1993; Zhdanov & Keller 1994), Zhdanov and his
coworkers have proposed a 'reverse' continuation of the
field in a manner similar to seismic-wave migration. Lee,
McMechan & Aiken (1987) and Levy, Oldenburg & Wang
(1988) have applied seismic techniques for layered media
with complex matrices for propagator leading to the
difficulty of strong variations of the amplitudes. In both
cases, the complete Green function is used and one may
think that the success of many seismic migration techniques
based on the separation of traveltime fitting and amplitude
adjustment has not yet been exploited.
Does a way exist to construct an asymptotic solution for
the diffusion electromagnetic problem which will be helpful
in an 'electromagnetic' imaging? It will be a more direct way
than to convert electromagnetic data into pseudo-seismic
data in order to use standard techniques of seismic imaging.
Recent papers by Lee, Liu & Morrison (1989), Gilbert &
Virieux (1991) and Lee & Xie (1993) have suggested this

latter strategy to increase the stability and accuracy of

electromagnetic imaging. In a recent paper, Nekut (1994)
has achieved this direct strategy for phase/traveltime
interpretation.
We shall show that, indeed, it is possible to construct a
solution which might be useful for least-squares inversion of
electromagnetic data. The inversion we propose comes
directly from recent works on efficient techniques to
construct seismic images (see Jin et al. 1992, for cited
references).
SOLUTION FOR A HOMOGENEOUS
MEDIUM
This section specifies our equations and the construction of a
solution for a homogeneous medium. This construction will
be such that extension to inhomogeneous media is possible.
In a 3-D medium with conductivity a, and permeability
p o without sources of displacement current, Maxwell
equations for the electric field E and the magnetic
susceptibility field H reduce to the following equations

_
aH --LVxE,
at

Po

aE
u~,E+ E - =
at

V x H,

where X denotes the cross product. For the transverse

electric mode (TE), the electric field E has one constant
direction which can be defined as the y-direction. The vector

857

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

SUMMARY
W e propose an asymptotic theory for diffusive electromagnetic imaging. Three steps
are required t o perform this imaging. (1) A high-frequency solution is first
constructed which mimics the one usually found in wave-propagation phenomena.
(2) This solution, valid for a smooth continuous description of the resistivity in the
medium, is used in a first-order Born approximation leading to a linear relation
between the resistivity perturbation of the subsurface and the perturbation of the
electric signal obtained at the free surface. (3) This linear relation is asymptotically
inverted by using an iterative quasi-Newtonian inversion based on a least-squares
criterion developed by Jin et al. (1992). Although the extension to smooth
heterogeneous reference medium is possible, we have only tested the inversion
scheme for homogeneous reference media as Zhdanov & Frenkel (1983) previously
did with another method.

858

J . Virieux et al.

E can be written (0, E,,, 0) and we assume invariance in the

y-direction. The magnetic susceptibility field H has only two
components ( H r ,0, H:). We found

(3)
which might be considered as an extension of the wave
equation with the extra diffusive term u,,aE,ldt. Although
transformation exists to convert diffusion into propagation
(Levy et al. 1988; Lee et al. 1989; Gibert & Virieux 1991),
deep differences arising from the diffusive extra term are
found which prevent a straightforward parallelism between
techniques used in propagation and in diffusion.
The separation between the T E mode and the TM mode
is not more valid at interfaces of zeroth-order discontinuity
for the conductivity (Wait 1981). Boundary conditions have
to be applied which couple the TM and TE modes. We do
not consider interfaces in the following and we shall
consider only a smooth continuous medium necessary for
our asymptotic imaging.
Going to the Fourier domain with the following sign
convention given by
(4)
and keeping the same notation for the Fourier transform of
the function E,, we deduce the Helmholtz equation:

again similar but different to the one deduced from wave

equation. At high frequency, the diffusive term in u,w can
be neglected compared with
and we are going back
strictly speaking to the wave equation used for instance in
high-frequency georadar imaging (Davis, Lytle & Laine
1979; Radcliff & Balanis 1981). Here we are interested in a
lower frequency range where the term in
can be
neglected instead, and eq. (5) reduces to

A generating solution of this equation is a complex plane

whose Fourier transform is 2 K , ( ' / - i o s r ) .

The
function K , is the modified Bessel function of zero order.
Taking the approximation of the function KO for large
arguments, an asymptotic solution can be written for the
2-D case as follows

where r is the distance between the line source (0,O) and the
receiver (x, z ) and H ( t ) is the Heaviside function. The time
4
the
at which the solution is maximum is t = ~ o u , , r 2 /while
mean square distance of the diffusion from the source at
time t is 4t/p,uo (Carslaw & Jaeger 1978, p. 256). In
practical applications, the penetration depth of diffusive
electromagnetic signal from the surface of the Earth ranges
from a few kilometres up to many hundreds of kilometres
for commonly used frequencies (Cagniard 1953).
The solution (8) must be analysed carefully when
compared with the exact solution (7). We must underline
that we have only constructed an approximate asymptotic
solution in the 2-D homogeneous case. Moreover a
frequency-dependent scaling factor 11appears in
the solution. This factor came for the spatial extension of
the source. One can find that a similar factor 1 1 6 appears
in the 2-D asymptotic solution of the wave propagation
corresponding to the Hilbert transform. As for the wave
propagation, this factor will be retained during our
construction of an asymptotic solution for a smoothly
varying medium. Then, it will be introduced at the end of
the procedure in the final solution to take into account the
geometry of the problem as done for wave propagation
(Bleistein 1984).
Obviously, this solution has a specific damping and
oscillating factor we hope to observe for a smooth
inhomogeneous medium. Following the basic idea leading to
the high-frequency approximation for the wave-equation
solution (Aki & Richards 1980), we are going to construct
an asymptotic solution for the diffusion equation.
We must underline that similar exact solutions can be
constructed for 1-D and 3-D geometries. These solutions in
the Fourier domain display the same specific damping and
oscillating factor without any approximation which was not
the case for 2-D geometry. Consequently, the asymptotic
solution we are going to construct is the exact solution for
the homogeneous 1-D and 3-D cases, while it is only an
approximation for the homogeneous 2-D cases. This means
that the 2-D geometry is the worst situation to be considered
for asymptotic diffraction tomographic reconstruction. The
reconstruction method we propose will use exact Green
functions when the 1-D or 3-D background media are
homogeneous. This validates the use of the asymptotic

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

For the transverse magnetic mode (TM), simidar equations

are obtained. This system of equations shows similarity but
also differences with the system of equations for SH motion
in an elastic medium. The intrinsic difference comes from
the term uoEv which introduces the diffusion. The system
(2) can be transformed into the telegraph equation

wave exp(k.x) with - i w u , p o = k Z . The selection of the

square root of -i is such that the plane-wave amplitude
should decrease with distance.
Let us briefly recall what are solutions in a homogeneous
medium for an impulsive signal -4n8(r). Let us consider a
line source. The solution will be

Diffusive electromagnetic imaging

solution and our numerical attempt to estimate the accuracy
of 2-D geometrical reconstructions.

the origin, the diffusive and propagation solutions are,

respectively,

ASYMPTOTIC SOLUTION F O R A SMOOTH

MEDIUM

I,\

Let us assume now that the medium has a smooth variation

of the conductivity u(x). We look for a solution in the
frequency domain which behaves as
(9)

(Vr)' = u(x)p,

+ 2VA, - V Z = 0
VZAk-,- AkVZr+ 2VAk

-(lek)

in a 3-D medium. These solutions have a very simple

geometrical interpretation. The wave solution is the delta
source signal shifted by the traveltime T and scaled by A,,
(Fig. 1). The diffusive solution is more complex to analyse:
the delta source signal is transformed into a localized
damped function maximum at r2/6 (Fig. 2) which is similar
to a phase shift although the diffusive energy stays around
the source. The geometrical effect of the 2-D medium will
add a tail to both signals, while reducing down to the 1-D
geometry again deforms the signal. The geometry of
medium influences the propagation and diffusion of the
signal as summarized in Fig. 1 and Fig. 2. The diffusion
solution (11) adds an extra difficulty because it exhibits a
time decrease of l/@ which is converted in a l / f decrease
in the 2-D case and l / f i i n the 1-D case. In our numerical
illustrations, we shall only use the 2-D geometry for both
sources and medium.
The smearing of the diffusive solution is the feature
preventing any simple reconstruction of the conductivity of
the medium. Our chance is that this smearing is controlled
by the pseudo-phase function in the asymptotic solution.
This chance is also a difficulty: the pseudo-phase might be
split into a damping and a phase and the damping be as fast
as the oscillation of the phase, preventing any arguments
based on the rapid oscillation of the phase term as we shall
see for the Hessian estimation in the inversion procedure.
The eikonal equation [Vz(x)12= o(x)p, allows ray tracing
( b r v e n y , Molotkov & PSenCii 1977) with a 'pseudoslowness' vector defined by p = V z at each point sampled
along the ray. The product u(x)p, can be compacted into
one single notation as u2(x) where u might be called
pseudo-slowness. For the eikonal, ray-tracing equations are
deduced as soon as a sampling parameter is selected along
rays. Let us define s as the curvilinear coordinate on the ray.

(10)

V r = 0 k > = 1.

WAVE
lD

The first two equations retain our attention because solving

them will give the pseudo-phase r(x) and the amplitude
A,(x) used for the zero-order term of the solution. These are
identical to both the eikonal and the transport equations
used for ray tracing of seismic waves and permit us to
compute r(x) and A,(x) with any standard ray-tracing
program (see LambarC 1992, for an example of efficient ray
tracing used for imaging). The intrinsic property of the
diffusion is contained in the ansatz we have taken. The
signal has its amplitude split in two parts: geometrical
spreading A, and a 'bulk' attenuation contained in the
pseudo-phase function. These two effects contribute in a
different manner to the electromagnetic diffusion amplitude.
Let us compare more precisely the zero-order term of eq.
(9) with its propagative equivalent for which the traveltime
of the wave is taken as T. For a Dirac source at

Impulsive source

rJ Heaviside step

[Hilbert function

Figure 1. Propagation of the impulsive signal from a source along a

ray. Note the difference of shape due to the medium geometry.

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

This ansatz for the solution is the one used in the ray theory
except that the io term for wave propagation has been
replaced by -Gfor diffusive transport. It has been
previously proposed by Tikhonov (1965) on a purely
theoretical basis. This term has no obvious property in the
Fourier domain as the translation property of the
exponential terms. Nevertheless, we have already constructed an analytical solution in the previous section with
this factor. The function 7 will be called the pseudo-phase
function following previous works on the link between
diffusion and propagation equations (Lee et al. 1989; Gibert
& Virieux 1991). The dimension of the pseudo-phase is the
square root of time.
For the wave-propagation equation, many other ansatz
have been suggested to take into account specific
geometries: ray theory singularities at caustics (Ludwig
1966), evanescent waves (Choudhary & Felsen 1973). A
seismological classification has been given by Chapman
(1985), while a review of different ansatz has been given by
Borovikov & Kinber (1974). These different ansatz show the
interest in developing specific asymptotic solutions for each
problem at hand.
Inserting the ansatz (9) in eq. (6) which is still valid for a
continuous inhomogeneous conductivity u ( x ) , we deduce
series of equations in powers of G:

A,V27

859

860

J. Virieux et al.

DIFFUSION
1D

Let us underline again that this solution is only

approximate in a 2-D homogeneous medium while the
corresponding solution for 1-D and 3-D homogeneous
media are exact solutions. Considering asymptotic solutions
is a valid assumption, especially when the background
medium is nearly homogeneous.

BORN APPROXIMATION

with wave propagation arising from the medium geometry.

u(x)iwp,,E,,(r,

r,)

+ VZE,,(r,r,) = -4zS(r

r,)

(16)

is expanded into
The transport equation which is the second one we are
interested in gives the evolution of the amplitude term as we
move along the ray. One finds an estimation of the
amplitude at s from s

a,(x)iopo AEy+ Vz AE,, = -4n[iwp,Aa(x)EY/4z],

(17)

the solution of which can be written as a convolution of the

Green function E;, solution of the left-hand side of eq. (17),
and the source term i w p o Aa(x)E,/4a leading to

AE,(r, r5, w )
where J is the Jacobian used in the definition of an
elementary surface perpendicular to the ray parameterized
by two coordinates y , and y z associated with the curvilinear
coordinate s.
The electrical field will be in the frequency domain

where S ( w ) is the Fourier transform of the time function of

the source and # is the intensity of the electromagnetic
source. The pseudo-phase z is equal to r
m = r u(r).
In order to estimate the intensity of eq. (13), we might
look at a canonical problem where the high-frequency
solution is known. Fortunately, a solution given by eq. (8)
exists for a homogeneous medium which has been
constructed in the previous section. Because the Jacobian is
equal to r in a 2-D homogeneous medium, we deduce the
intensity for an isotropic source at the origin

r, x , w ) Au(x)Ey(x,r,, w ) dx2,
where the domain of integration A is over diffracting points
x. The first-order Born approximation is obtained by
replacing E, with E; in the integral, leading to a linear
functional between Au and AE,,:

AE,(r, rs, 0)
(19)

valid for a smoothly varying medium by identification with

the solution (8). The asymptotic solution in a 2-D smooth
varying medium for a rs source is

where we have again introduced the specific frequencydependent factor for the 2-D geometry already mentioned in
the previous paragraph.

Assuming now that solutions have the asymptotic form in

the smooth reference 2-D medium

we obtain for a given couple of (source, receiver) the simple

relation

AEy(r, rs, 0)

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

Figure 2. Diffusion of the impulsive signal from a source along a

pseudo-ray. Note the less dramatic difference of shape compared

Let us now consider a perturbation of the conductivity.

Because in this paper we assume an invariance along the
direction y , we consider a 2-D reference medium with a
smoothly varying conductivity ao(x). For a point source (a
line source in a 3-D medium with one invariant direction) at
the source position r,, the Green function at position r is
denoted E;(r, rJ. We want to study media slightly different
from this smooth reference medium by a conductivity
perturbation given by a ( x ) = u o ( x )+ Au(x). We assume that
the magnetic permeability po remains constant. The Green
function E,, will be split into the known Green function E:
and the perturbation AE,,. The equation

Diffusive electromagnetic imaging

861

where

come from the product of the Green function connecting the

source and the diffracting point and the Green function
connecting the diffracting point and the receiver.
The asymptotic form has simplified the propagation of the
field because ray tracing is symmetrical: the asymptotic
Green function from the diffracting point up to the receiver
is equivalent to the one from the receiver down to the
diffracting point.
Going back to the time domain is possible with a very
simple integral expression which is the Fourier inverse of eq.
J
(21):

T(r, y

5 rS)2

-1
2

e-r(r.x.rs)2/41

H ( t ) dx.

(23)

The kernel of the expression (23) in time has a maximum

amplitude at t,,, = T2(r)(1/2 + l/G) with a sharper
decrease for time lower than t,,, (exponential term of time)
than for time higher than t,,, (inverse power of time). As a
consequence, main contributions to the integral are located
on a shell which can be made finite by numerical cut-off of
the kernel. We must underline that the situation is slightly
more complex than for 2-D seismic imaging: the
contribution to the seismogram comes from an isochrone
surface along which the image is concentrated and defined
by a t equal to traveltime (Bernard & Madariaga 1984;
Spudich & Frazer 1984). In diffusive electromagnetic
imaging, we have generalized the so-called isochrone
property to a less local but still finite shell.

-3

180

5
u

;a -5
Q)

g o
a
a

-6
0

&

-7

-8
- 11

3
L

log(frequency)

-180
-1

,i/
2
log(frequency)

I
5

Figure 3. Comparison between the exact solution (solid line) and the asymptotic Born approximation (dashed line) for a line diffractor. The
difference arises essentially from the asymptotic approximation. Both amplitude and phases are plotted.

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

Expressions (21) and (23) are valid in a 2-D or a 3-D

medium with an invariance in the y-direction for the
definition of the T E mode. Expressions (21) in the
frequency domain and (23) in the time domain are our basic
linearized forward modelling for which we shall construct an
inverse operator. The technique has been previously
proposed by Jin et al. (1992) based on the work performed
by Beylkin (1985) and Beylkin & Burridge (1990) among
others on asymptotic Radon transforms.
We check the Born approximation with two simple
examples. In the first one, we shall consider a line diffractor
of conductivity aI= 0.1 S m- embedded into a homogeneous medium of conductivity uo= 10-S m- I . We compute
the Born approximation for an incident isotropic line source
of unitary strength diffracted by a line located lOOm just
below it for frequencies between 10 - Hz and 10 Hz. Fig. 3
shows the exact solution compared with the Born asymptotic
solution. We must say that the difference in the scattered
field arises essentially from the difference between the
complete 2-D solution and its high-frequency approximation. As the frequency becomes higher, the discrepancy in
amplitude and in phase becomes negligible. In 2-D and 3-D
geometries, we should have a perfect agreement because the
embedding medium is homogeneous.
In the second example, we compare the solution for a
cylinder diffractor at a depth of 100m in a homogeneous
medium of conductivity 0.01 S m- excited by a line source.
We still find a good agreement while we are summing the
expression (23) inside the cylinder to recover the diffracted
field (Fig. 4). When the conductivity contrast becomes more
important, the disagreement between the complete solution
and the Born approximation increases at higher frequencies
(Fig. 4a). As the cylinder becomes wider, the agreement
decreases when the radius increases (Fig. 4b). The Born
approximation is valid at relatively high frequencies for

862

J . Virieux et al.

(4

180
Slm

-1
----

I
I

2
log(frequency)

-1

2
log(frequency)

(bl

180

/'

50 m

--

d)

z
s

-5
-1

2
log(frequency)

-180
-1

2
log(frequency)

Figure 4. (a) 'Comparison between the exact solution and the Born approximation for a cylinder diffractor of radius 20m. Notice the
disagreement introduced by the Born approximation when the contrast of conductivity increases. For clarity, only the phase response for the
0.1 S m-' case is included. (b) Comparison between the exact solution and the Born approximation for a cylinder diffractor of conductivity
0.1 S m-'. Notice the disagreement introduced by the Born approximation when the size of the cylinder increases. For clarity, only the phase
response for the 50 m radius case is included.

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

a
-180

-5
-1

Diffusive electromagnetic imaging

small perturbations of the conductivity and for small sizes of
diffractor objects. We shall see in our inversion the effects of
these disagreements. These effects will be the main
limitation of our reconstruction technique.
We must remember that the main misfit arises from our
transformation of the 2-D Green function into its
high-frequency approximation. For 3-D and 1-D homogeneous media, the Green function agrees with its
high-frequency approximation as we have checked. Therefore, we d o expect better results with the Born
approximation for 3-D structures.

863

Inversion by a least-squares method

The definition of the misfit function requires an explicit
definition of the inner product in data space:

(AE, 1 AE:)y =

dw AE,*(r,ra, w )

rp,r

X Q(r, xo, r,,

AWr, rs, 0 )

0)

(25)

where * denotes the complex conjugate. The sum and the

integral extend over the data space 9. The covariance
matrix 0 is diagonal by construction with elements:

ASYMPTOTIC I N V E R S I O N

AEy = G Au

(24)

where G : A - + 9 is the integral operator of the two-way

Green functions (see eqs 21 and 23).
The solution of the linearized inverse problem consists of
finding the inverse of operator G. We look for an
approximate solution to the inversion of eq. (24) by the
optimization method which minimizes a misfit function
between observed and calculated electric signals (see
Tarantola 1987, for a description of inverse theory). We
adopt the least-squares norm Y2of the sum of the square of
the difference between observed and predicted signals at
each frequency. We shall see that this criterion leads to a
quasi-Newtonian iterative method as shown by Jin et al.
(1992).

Source

Free Surface

Receiver

where p(r, %, r,) = Vt(r, x,,, rs), is the gradient of the

two-way pseudo-phase function r (Fig. 5).
The particular form of the covariance matrix Q is
introduced in order to correct for the decrease with distance
and for the spectral content of the Green function. One can
see in the definition of covariance matrix 0 only the
influence of the geometrical spreading Aoand not the effect
of the diffusion. The covariance matrix depends on q,the
coordinate of the point at which the model will be inverted.
While the choice of a diagonal matrix is rather standard in
the optimization method, Q is usually independent of where
we are in the model space. Our choice borrowed from the
work done by Jin et al. (1992) for seismic profiles is in fact a
preconditioning applied to the gradient of the misfit
function. The additional complex form of 0 comes from the
specific diffusive tail.
The covariance matrix 0 upgrades weak late signals:
instabilities might arise when noise exists in the late part of
the time signal or when the reference medium leads to
strong defocusing in a time window of the signal where
energetic pulses are observed. Summing over the data
acquisition geometry reduces these incoherent instabilities
considerably. If not enough, a priori information defining
the maximum amplitude for perturbation of conductivity
stabilizes the procedure.
For the inverse problem we also need a definition of the.
inner product between two functions A u and A d in the
model space A:

(Au I Au)A=

i,Au*(x) Au(x) dx.

The misfit function is defined by

S(Au) = 1/2(AE; - G AU I AE,- G A u ) ~ ,

\ /
Diffractor poin

(28)

where AE; are the observed data and GAu are the
synthetic ones estimated through eq. (24). We formulate the
inverse problem as
find A c : min [~(Au)],
ACT

(29)

and we obtain the classical system of normal equations:

GG A u = Gi AE,
Figure 5. Geometry of pseudo-rays around the diffractor point.

(30)

for all x where Gt is the adjoint operator to the forward

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

In order to pose the inverse problem properly we must

define both model and data spaces, and the operators
between these two spaces. The model space A(x) is the
space of all possible perturbations of the conductivity A u at
each point x of the medium. The data space consists of all
electric perturbations of the T E mode for many sources rs
and receivers r at different frequencies. For simplicity of the
analytical developments, we shall consider Fourier transformed data, but at the end we shall return to the time
domain to take advantage of the finite domain of A from
where significant contributions to the electric signal come.
This signal is defined on the domain 9 ( r , r-, w ) . The
linearized forward problem can be written in a compact
form

J . Virieux et al.

864

integral operator G, and defined by the classical relationship

(AE, G Acr),/

= (Gt A E v

1 ACT),#,

(31)

Let us underline that the choice of the factor 0 has

modified the t,,, compared with the one obtained for the
forward linearized expression.

while the generalized inverse C g

is equal to (GtG)-'Gt.

Asymptotic expression of the Hessian

Estimation of the gradient
The adjoint operator G' might be expressed by a kernel YC
with the following expression

The operator H - ' in eq. (34) is the formal inverse of G'G

and is very difficult to estimate for many inverse problems
(see Jin et al. 1992, for a short discussion and Tarantola
1987, for an extensive review). In our approach, an
asymptotic estimation is possible. We found from eqs (21)
and (32) that

Yt(r, x, r,, w )
-1
2

= -poA(r, x, rs)Q*(r, x,), rs, w ) G *

G*r('.'Js)

(33)

which reduces to

The formal solution of the system (30) is

Am(%)

= H-'(xo, x ) Y " ( x ) ,

(34)

where

y O ( x )= Gt(r, x , r,) AE;"(r,

rs, w )

(35)

is the gradient of S ( A a ) at A a = O and H is the Hessian

giving the information of the curvature of the misfit
function. Explicitly the gradient is given by

lp(,.,x,

),

112

-*

r(r.x.rs)

AE;"(r,

ra,w )

(36)

in the frequency domain. Expressing AE;hS(r, rsrw ) as the

Fourier transform of AEyhs(r, rs, t ) leads to the following
expression:

The term under the integral over frequency has a damping

term which prevents any arguments of equally balanced
contributions for the whole frequency spectrum as it is for
seismic inversion. This is a fundamental difference that we
have underlined previously but we can argue that the main
contribution to the integral is when the phase of the
exponential is zero because the diffusive damping term is a
slowly decaying function. If the background medium is
sufficiently smooth, this occurs when x is close to %. With
local expansion of the pseudo-phase and the amplitude

(37)
We have used the complex conjugate operation giving
6*
= fi and the transformation of w into w ' = - 0 . By
doing so, we are able to write the gradient as a single
integral in the time domain
dofie-

x {cos

[mp

\/?;r(r.XO.rs'e-~~'(~-~))

( x - X")]

- sin -\/;siz p (x - % ) I }

The time-domain integration in the asymptotic inversion

kernel is around t,,, = z(r, x , rJZ/2 which corresponds to a
shell in the diffraction domain A. We have found the
equivalent of the isochrone summation for the asymptotic
seismic migration related to the principle of coinciding time
of Claerbout (1985). Apart from an integration over a time
interval instead of a single value for seismic migration, we
have an expression that is manageable for practical
implementation.

which has significant values when x is close to x.,

assume the more drastic approximation,

H ( x , %) - H(x0, X

0 ) W - %),

(42)
We shall
(43)

justified by the exponential decrease away from the position

%. The Hessian can now be estimated as

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

where the reconstructed function is obtained at position x

assuming a covariance Q depending on % . F r o m eq. (31)
and from the forward problem (21),

DifSusive electromagnetic imaging

The final and precise shape of the operator H is mainly
controlled by the data acquisition geometry and eq. (44)
should be considered only as an estimation for an iterative
inversion method. The nearly diagonal structure of the
operator H makes the computation of the inverse easier and
leads to

M%)
= 1/H(%, xg)

Yo(%).

(45)

Iterative quasi-Newton inversion method

Let us note h the approximation of the Hessian H . The
quasi-Newton solution of the inverse problem (29) is
Au(x)"+' = Au(x)"

+K ' y "

(46)

where y" is the gradient of the misfit function S(Au)

calculated around the value of A u at the nth iteration:
yn = G + ( A E O ~ SG A ~ " ) .

(47)

As shown by Jin et al. (1992), the iterative method

converges to the limit

A test with the analytical complete solution

lim ha" = (G-gG)-'G-gAEohS

n--r-

which shows that the iterative method corrects for any bias
in the Hessian estimation.
A self-consistent test
Before performing an inversion for a complete solution, we
want to test the internal coherence of the linear inversion we
propose. We have computed the Born approximation for a
cylinder of radius 20 m at a depth of 100 m. The conductivity

-200 m

of the cylinder is 0.1 S m-' embedded in a homogeneous

medium of conductivity 0.01 S m-I. We consider a single line
source right above the cylinder and five receivers at the free
surface distributed by steps of 50m on both sides of the
source.
Figure 6 shows both the true cylinder and the recovered
conductivity contrast. The recovered maximum value is only
0.01 Sm-' while we should have found a value around
0.09Sm-'. The wider extension of the recovered image
explains this low maximum value. With the blurred image of
Fig. 6, the predicted signals can explain most parts of the
synthetic signals. Unfortunately, this is a drawback of the
diffusion phenomenon which implies an integration over a
pseudo-isochrone shell to recover the image. We find the
typical 'smile' associated with the data acquisition geometry.
During iterations, the misfit function decreases from 130 at
the first iteration down to 66 at the 50th iteration.
If we consider other sources translated on the horizontal
axis, we are able to stack different pictures of the medium
and to improve our resolution. The cylinder shape is better
resolved as expected in Fig. 7 with a noticeable reduction of
the 'smile'. The misfit function goes down to 44 when
normalized by the number of sources, lower than the misfit
function for a single source. The extension of the recovered
image still biases strongly the maximum amplitude of the
conductivity contrast. We do expect better results when
considering other geometries than the surface-to-surface
geometry.

The previous test has been performed using Born

computation as synthetic data. What happens when
considering the exact solution of diffusion by a cylinder? We
have already compared true and Born solutions in the
frequency domain in a previous paragraph. From these
analytical solutions in the frequency domain, we compute
the solution in the time domain for the same geometry by
fast Fourier transform. Then, we perform asymptotic
inversions similar to the ones previously computed with
Born synthetic solutions.

200 m
Om

200 m
Figure 6. Conductivity contrast recovered after 50 iterations for one source at the free surface above the cylinder. Inverted data are the Born
solution.

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

In order to check the final gradient expression (38) and

the final Hessian expression (44)o we have computed the
reconstructed conductivity at a single diffracting point %
from the scattered field computed by the Born approximation. We have deduced the exact conductivity contrast
between the homogeneous referenceJ medium and the
scattering point. This analytical checking makes us confident
in our expressions for imaging diffusive fields.

865

866

J. Virieux et al.

-200 m

200 m
Om

Figure 7. Conductivity contrast recovered after 50 iterations for five sources distributed at the free surface on both sides of the cylinder.
Inverted data are the Born solution.

For inversion of the exact solution by a scattering

cylinder, we find a deeper image with the typical smile
associated with one source geometry (see Fig. 8). The
maximum amplitude of conductivity contrast is as high as
0.03 S m-. The amplitude difference between the Born
solution and the complete solution in the time domain easily
explains this amplification (see Fig. 9). The slight shift of the
maximum diffusion pulse between these two solutions
explains also why the image is deeper than the true cylinder.
The reduction of the misfit function is from 475 at the first
iteration to 166 at the 50th iteration. Fig. 9 gives an example
of the residual signal left unexplained by the analytical
inversion.
Performing the inversion using more sources removes the
smile geometry as shown in Fig. 10 but does not improve
the position of the strongest conductivity amplitude which is
still deeper than the true position. Changing the geometry
structure by adding data recorded inside wells might
improve this position problem already noticed with the Born
approximation and amplified when using the exact solution.
We see in this final synthetic example the difficulty

-200 m

inherent in diffusion phenomena which blurs the image by

the spatial extension of the so-called pseudo-isochrone shell
for image recovering. This distortion of the image arises also
from the Born approximation which, at intermediate
frequencies, shows relatively poor agreement. Starting from
a better reference medium with the already included
low-frequency content of the conductivity might improve the
picture, because, for the relatively higher spatial frequencies, the associated diffusion tail in time will be easier
to handle.
DISCUSSION A N D CONCLUSION
We have developed an analytical inversion for diffusive
electromagnetism. This formalism draws its features from
the seismic inversion approach. We extend the isochrone
line concept of wave propagation to a pseudo-isochrone
shell concept for diffusion. By doing so, we were able to
construct an inversion scheme with explicit formula for the
gradient and the Hessian of the misfit function. Because
these formula are analytical, they are insensitive to noise

200 m
Om

200 m
Figure 8. Conductivity contrast recovered after 50 iterations for-onesource at the free surface above the cylinder. Inverted solution is the exact
analytical solution.

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

200 m

Diffusive electromagnetic imaging

861

z
1.0 -

0.5 -

0.0000

0.0001

0.0002

0.0003

Time (sec)
Figure 9. Synthetic electric signal (triangles) recorded at the receiver above the cylinder as well as the predicted signal (circles) computed for
the cylinder image. The residual is also shown (diamonds), as well as the synthetic Born solution used in the forward modelling (dots).

which are simply not focused back into the medium when
incoherent as already shown for seismic data (Lambar6 et al.
1992).
The use of the Born approximation has limitations when
we try to recover a spatially extended object on a
homogeneous background. The possible extension to a
smooth inhomogeneous reference medium will improve the
position of heterogeneities by requiring only a fit of the
relatively high-frequency part of the electric signal when the
smooth background velocity has already been obtained. In

-200 m

that sense, this reconstruction technique is only a partial one

in terms of spatial resolution.
In any case, we believe that, because the main diffusion
pulse which is the one fitted in our present applications
would have been already contained in our initial model, the
spatially high-frequency content of the image will be better
resolved. Of course, we might analyse the effect of noise for
this inversion scheme as well as we might analyse
performances on real data. This will be the purpose of
object work.

200 m

Om

200 m
Figure 10. Conductivity contrast recovered after 50 iterations for five sources distributed at the free surface on both sides of the cylinder. The
inverted solution is the exact analytical solution.

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

-0.5

868

J. Virieux et al.

ACKNOWLEDGMENTS
This work was partly supported by the CNRS-INSU through
the Tomographic group and MEN-DRED through Jeune
Equipe RuaDE. We thank two anonymous reviewers for
their helpful comments.

REFERENCES

Downloaded from http://gji.oxfordjournals.org/ at Texas A&M University Evans Library on June 4, 2014

Aki, K. & Richards, P., 1980. Quantitative Seismology: Theory and

Methods, W. H. Freeman and Co, San Francisco, CA.
Barnett, C.T., 1984. Simple inversion of time-domain electromagnetic data, Geophysics, 49, 925-933.
Bernard, P. & Madariaga, R., 1984. A new asymptotic method for
the modelling of near-field accelerograms: Bull. seism. SOC.
Am., 74,539-559.
Beylkin, G., 1985. Imaging of discontinuities in the inverse
scattering problem by inversion of a causal generalized Radon
transfrom, J. math. Phys., 26, 99-108.
Beylkin, G. & Burridge, R., 1990. Linearized inverse scattering
problems in acoustics and elasticity, Wave Motion, 12, 15-52.
Bleistein, N., 1984. Mathematical Methods for Wave Phenomena,
Academic Press, Inc., Orlando, FL.
Borovikov, V.A. & Kinber, B.Y.E., 1974. Some problems in the
asymptotic theory of diffraction, Proc. IEEE, 62, 1416-1437.
Cagniard, L., 1953. Basic theory of the magneto-telluric method of
geophysical prospecting, Geophysics, 18, 605-635.
Carslaw, H.S. & Jaeger, J.C., 1978. Conduction of Heat in Solids,
Oxford University Press, Oxford.
b r v e n f , V., Molotkov, LA. & PSenEii, I., 1977. Ray Method in
Seismology, Charles University Press, Prague.
Chapman, C.H., 1985. Ray theory and its extension: WKBJ and
Maslov seismograms, J. Geophys., 58, 27-43.
Choudhary, S. & Felsen, L.B., 1973. Asymptotic theory for
inhomogeneous waves, IEEE Trans. Antennas Propag., AP-21,
827-842.
Claerbout, J.F., 1985. Imaging the Earths Interior, Blackwell,
Oxford.
Davis, D.T., Lytle, R.J. & Laine, E.F., 1979. Use of high-frequency
electromagnetic waves for mapping an in situ coal gasification
burn front, In situ, 3, 95-119.
Gibert, D. & Virieux, J., 1991. Electromagnetic imaging and
simulated annealing, J. geophys. Res., 96, 8057-8067.
Jin, S., Madariaga, R., Virieux, J. & Lambart, G., 1992.
Two-dimensional asymptotic iterative elastic inversion,
Geophys. J. Int., 108,575-588.
Lambare, G., Virieux, J., Madariaga, R. & Jin, S., 1992. Iterative
asymptotic inversion: applications for one parameter,
Geophysics, 57, 1138-1 154.
Lee, K.H. & Xie, G., 1993. A new approach to imaging with

low-frequency electromagnetic fields, Geophysics, 58,780-796.

Lee, S . , McMechan, G.A. & Aiken, C.L., 1987. Phase-field imaging:
the electromagnetic equivalent of seismic migration,
Geophysics, 52, 1678-693.
Lee, K.H., Liu, G. & Morrison, H.F., 1989. A new approach to
modeling the electromagnetic response of conductive media,
Geophysics, 54, 1180-1192.
Levy, S., Oldenburg, D.W. & Wang, J., 1988. Subsurface imaging
using magnetotelluric data, Geophysics, 53, 104-117.
Ludwig, D., 1966. Uniform asymptotic expansion at a caustic,
Commun. Pure appl. Math., 29, 215-250.
Nekut, A.G., 1994. Electromagnetic ray-trace tomography,
Geophysics, 59, 371-377.
Radcliff, R.D. & Balanis, C.A., 1981.Electromagnetic geophysical
imaging incorporating refraction and reflection, IEEE Trans.
Antennas Propag., AP-29,288-292.
Smith, J.T., 1988. Rapid inversion of multi-dimensional magnetotelluric data, PhD thesis, University of Washington.
Spudich, P. & Frazer, L.N., 1984. Use of ray theory to calculate
high-frequency radiation from earthquake sources having
spatially variable rupture velocity and stress drop, Bull. seism.
SOC. A m . , 74,2061-2082.
Tarantola, A., 1987. Inverse Problem Theory, Elsevier Science
Publishers BV, Amsterdam.
Tarits, P., 1989. Contribution des sondages Clectromagnitiques
profonds
Iktude du manteau supCrieur terrestre, ThPse
dEtat, UniversitC Paris 7.
Tikhonov, A.N., 1965. Mathematical basis of the theory of
electromagnetic soundings, Zh. Vychisl. Mat. Mat. Fiz., 5,
545-548.
Tripp, A.C., Hohmann, G.W. & Swift, C.M., 1984. Twodimensional resistivity inversion, Geophysics, 49, 1708-1717.
Wait, J.R., 1981. Wave Propagation Theory, Pergamon Press,
Oxford.
Weidelt, P., 1972. The inverse problem of geomagnetic induction, J.
Geophys., 38,257-289.
Weidelt, P., 1975, Inversion of two-dimensional conductivity
structures, Phys. Earth planet. Inter., 10, 282-291.
Zhdanov, M.S., 1988. Integral Transforms in Geophysics, p. 361,
Springer Verlag, Berlin.
Zhdanov, M.S. & Frenkel, M.A., 1983. The solution of the inverse
problems on the basis of the analytical continuation of the
transient electromagnetic field in the reverse time., J.
Geomagn. Geoelect., 35, 747-765.
Zhdanov, M.S. & Booker, J.R., 1993. Underground imaging by
electromagnetic migration, SEG Expanded Abstracts, pp.
81-82, 63rd Annual Meeting, SOC. Explor. Geophys.,
Washington, DC.
Zhdanov, M.S. & Keller, G., 1994. The Geoelectrical Methods in
Geophysical Exploration, p. 873, Elsevier, Amsterdam,
London, New York, Tokyo.