Retrieving the Elastodynamic Green’s Function of an Arbitrary Inhomogeneous Medium
by Cross Correlation
Kees Wapenaar
*
Department of Geotechnology, Delft University of Technology, P.O. Box 5028, 2600 GA Delft, The Netherlands
(Received 18 August 2004; published 16 December 2004)
A correlation-type reciprocity theorem is used to show that the elastodynamic Green’s function of
any inhomogeneous medium (random or deterministic) can be retrieved from the cross correlation of
two recordings of a wave field at different receiver locations at the free surface. Unlike in other
derivations, which apply to diffuse wave fields in random media or irregular finite bodies, no
assumptions are made about the diffusivity of the wave field. In a second version, it is assumed that
the wave field is diffuse due to many uncorrelated sources inside the medium.
DOI: 10.1103/PhysRevLett.93.254301 PACS numbers: 43.20.+g, 43.60.+d, 91.30.–f
Recently it has been shown by various authors that the
Green’s function of a random medium or an irregular
finite body can be obtained by cross correlating the re-
cordings of a diffuse wave field at two receiver positions
[1–5]. The resulting Green’s function is the wave field that
would be observed at one of these receiver positions if
there were an impulsive source at the other. This theoreti-
cal result has been successfully demonstrated with ultra-
sonic measurements [4–7] and seismic surface waves
[8,9]. The accuracy of the reconstructed Green’s function
depends on the amount of disorder of the medium pa-
rameters and the duration of the signal. Ideally the cross
correlations should be done in the equipartitioned regime
(where the net energy flux is equal to zero), which takes
place after sufficiently long multiple scattering of the
wave field between the heterogeneities in the disordered
medium [7].
An initially independent line of research deals with the
reconstruction of the seismic reflection response of a
deterministic medium from passive recordings of the
transmission response. Already in 1968 Claerbout
showed that the autocorrelation of the transmission re-
sponse of a horizontally layered earth yields the super-
position of the reflection response and its time-reversed
version [10]. The source in the subsurface may be a
transient or a noise signal; in both cases the source
signature in the reconstructed reflection response is the
autocorrelation of the source signal in the subsurface.
This method has been applied to microearthquake data
[11]. The derivation in [10] was strictly one dimensional.
Later Claerbout conjectured for the 3D situation that ‘‘by
cross-correlating noise traces recorded at two locations on
the surface, we can construct the wave field that would be
recorded at one of the locations if there was a source at
the other’’ [12]. Although it was not explicitly stated, this
conjecture applies to deterministic media: in exploration
seismology the earth is usually considered to be built up
of geological layers with smoothly varying properties,
separated by well-defined curved interfaces and faults
which act as the main reflectors; scattering due to disor-
der of the parameters within the geological layers is
generally considered a second order effect. Numerical
modeling studies have been carried out to confirm
Claerbout’s conjecture [13]. These modeling studies
showed that ‘‘longer time series, and a white spatial
distribution of random noise events would be necessary
for the conjecture to work in practice.’’ The cross-
correlation approach has been applied successfully to
helioseismic data [12,14]. Recently Claerbout’s conjec-
ture has been proven by the author [15,16]. The proof
also explains the empirical observations of the numerical
modeling studies.
In this Letter we derive a relation between the elasto-
dynamic Green’s function and the cross correlation of
observed wave fields that holds at the free surface of
random as well as deterministic media. The approach is
quite different from that in [1–7], which holds only for
diffuse wave fields. It is also different from the derivation
in [15,16], which is based on coupled one-way wave
equations for acoustic down-going and up-going waves,
assumes a certain degree of smoothness of the medium
parameters, and ignores evanescent waves.
The basis for our derivation is a reciprocity theorem,
which relates two independent elastodynamic states
(wave fields and sources) in one and the same medium
[17]. One can distinguish between reciprocity theorems of
the convolution type and of the correlation type [18].
Correlation-type reciprocity theorems contain correla-
tions between the wave fields and sources in both states.
Since it is our aim to retrieve the Green’s function from
the cross correlation of observed wave fields, the
correlation-type reciprocity theorem is a natural choice.
Let an elastodynamic wave field be characterized by
the space- and time-dependent particle velocity v
i
x;t
and stress tensor
ij
x;t. Here x x
1
;x
2
;x
3
and t de-
note the Cartesian coordinate vector and time, respec-
tively. Subscripts i and j take on the values 1, 2, and 3. We
use Einstein’s summation convention for repeated lower-
PRL 93, 254301 (2004)
PHYSICAL REVIEW LETTERS
week ending
17 DECEMBER 2004
0031-9007=04=93(25)=254301(4)$22.50 254301-1 2004 The American Physical Society
case subscripts. We define the temporal Fourier transform
of a space- and time-dependent quantity as
^
v
i
x;!
R
expj!tv
i
x;tdt, where j is the imaginary unit and
! the angular frequency. In the space-frequency domain
the particle velocity and stress tensor in an inhomoge-
neous anisotropic lossless medium obey the equation of
motion j!
^
v
i
@
j
^
ij
^
f
i
and the stress-strain relation
j!s
ijkl
^
kl
@
j
^
v
i
@
i
^
v
j
=2
^
h
ij
, where @
j
is the par-
tial derivative in the x
j
direction, x the mass density of
the medium, s
ijkl
x its compliance,
^
f
i
x;! the external
volume force, and
^
h
ij
x;! the external deformation rate.
We consider the ‘‘interaction quantity’’
^
v
i;A
^
ij;B
^
ij;A
^
v
i;B
, where the asterisk denotes complex conjuga-
tion and subscripts A and B are used to distinguish the two
independent states. Note that the product
^
v
i;A
^
ij;B
is the
Fourier transform of the cross correlation of v
i;A
x;t and
ij;B
x;t; a similar remark holds for ^
ij;A
^
v
i;B
. The elas-
todynamic reciprocity theorem of the correlation type is
obtained by applying the operator @
j
to the interaction
quantity, substituting the equation of motion and the
stress-strain relation for states A and B, using the sym-
metry relations ^
ij
^
ji
and s
ijkl
s
klij
, integrating the
result over a spatial domain D with boundary @D and
outward pointing normal vector n n
1
;n
2
;n
3
, and ap-
plying the theorem of Gauss [17]. This gives
I
@D
f
^
v
i;A
^
t
i;B
^
t
i;A
^
v
i;B
gd
2
x
Z
D
f
^
h
ij;A
^
ij;B
^
f
i;A
^
v
i;B
^
ij;A
^
h
ij;B
^
v
i;A
^
f
i;B
gd
3
x;
(1)
with the traction
^
t
i
at the boundary @D defined as
^
t
i
^
ij
n
j
. We apply this theorem to the wave fields in an
arbitrary inhomogeneous anisotropic lossless medium,
bounded by a free surface (Fig. 1). We choose @D such
that it consists of a part of the free surface, denoted by
@D
0
, and an arbitrarily shaped surface @D
1
inside the
medium. Furthermore, we assume that the sources
^
h
ij;A
,
^
f
i;A
,
^
h
ij;B
,and
^
f
i;B
in D are zero. Hence, Eq. (1) becomes
Z
@D
0
f
^
v
i;A
^
t
i;B
^
t
i;A
^
v
i;B
gd
2
x
Z
@D
1
f
^
v
i;A
^
t
i;B
^
t
i;A
^
v
i;B
gd
2
x: (2)
We introduce sources in terms of boundary conditions at
the free surface @D
0
. This is possible since at a free
surface the traction is zero everywhere, except at those
positions where a source traction is applied. In state A we
apply a source traction in the x
p
direction at x
A
2 @D
0
with source function st. Hence, in the space-frequency
domain the traction for x 2 @D
0
is given by
^
t
i;A
x;!
x x
A
ip
^
s!, where x x
A
is a 2D Dirac delta
function and
ip
is the Kronecker delta function. The
observed particle velocity at x due to the source at x
A
is
expressed as
^
v
i;A
x;!
^
G
v;t
i;p
x; x
A
;!
^
s!, where x can
be anywhere at the free surface or in the medium.
^
G
v;t
i;p
x; x
A
;! represents a Green’s function, with the fol-
lowing notation convention: the two coordinate vectors
between the brackets (here x and x
A
) represent the ob-
servation point and the source point, respectively; the
superscripts (here v and t) represent the observed quantity
(velocity) and the source quantity (traction), respectively;
the subscripts (here i and p) represent the component of
the observed quantity and the source quantity, respec-
tively. Similarly, in state B we apply a source traction in
the x
q
direction at x
B
2 @D
0
with the same source func-
tion st. Hence,
^
t
i;B
x;!x x
B
iq
^
s! for x 2
@D
0
and
^
v
i;B
x;!
^
G
v;t
i;q
x; x
B
;!
^
s! for any x.
Substituting these expressions in the left-hand side of
Eq. (2), and using reciprocity of the Green’s function
[i.e.,
^
G
v;t
q;p
x
B
; x
A
;!
^
G
v;t
p;q
x
A
; x
B
;!], yields
Z
@D
0
f
^
v
i;A
^
t
i;B
^
t
i;A
^
v
i;B
gd
2
x f
^
G
v;t
q;p
x
B
;x
A
;!g
^
G
v;t
p;q
x
A
;x
B
;!j
^
s!j
2
2Ref
^
G
v;t
p;q
x
A
;x
B
;!gj
^
s!j
2
;
(3)
where Refg denotes the real part. In order to evaluate the
right-hand side of Eq. (2), we note that in states A and B
the velocities and tractions for x 2 @D
1
, due to the trac-
tion sources at x
A
and x
B
, can be expressed as
^
v
i;A
x;!
^
G
v;t
i;p
x; x
A
;!
^
s!
^
G
v;f
p;i
x
A
; x;!
^
s!;
^
t
i;A
x;!
^
G
t;t
i;p
x; x
A
;!
^
s!
^
G
v;h
p;i
x
A
; x;!
^
s!;
^
v
i;B
x;!
^
G
v;t
i;q
x; x
B
;!
^
s!
^
G
v;f
q;i
x
B
; x;!
^
s!;
^
t
i;B
x;!
^
G
t;t
i;q
x; x
B
;!
^
s!
^
G
v;h
q;i
x
B
; x;!
^
s!:
Note that the second Green’s function in each of these
equations is the reciprocal of the first and can therefore be
interpreted as an observation at the free surface (at x
A
or
x
B
), due to a source at x 2 @D
1
.In
^
G
v;f
p;i
x
A
; x;! the
observed quantity is velocity and the source is a volume
n
A
x
B
x
()ρ x
()
ijkl
s x
I∂
D
I∂
D
1
0
n
A
x
B
x
()ρ x
()
ijkl
s x
I∂
D
I∂
D
I∂
D
I∂
D
1
0
FIG. 1. Inhomogeneous anisotropic lossless medium,
bounded by a free surface.
PRL 93, 254301 (2004)
PHYSICAL REVIEW LETTERS
week ending
17 DECEMBER 2004
254301-2
force; in
^
G
v;h
p;i
x
A
; x;! the observed quantity is again velocity but the source is a specific type of deformation rate source
(actually
^
G
v;h
p;i
n
j
^
G
v;h
p;ij
, where the source for
^
G
v;h
p;ij
is a deformation rate tensor and n
j
is the j component of the normal
at @D
1
). Substituting these expressions in the right-hand side of Eq. (2) (and using Eq. (3) for the left-hand side) we
obtain
2Ref
^
G
v;t
p;q
x
A
; x
B
;!g
Z
@D
1
f
^
G
v;f
p;i
x
A
; x;!g
^
G
v;h
q;i
x
B
; x;!f
^
G
v;h
p;i
x
A
; x;!g
^
G
v;f
q;i
x
B
; x;!d
2
x: (4)
The left-hand side of this equation represents the Fourier
transform of the superposition of the time domain
Green’s function G
v;t
p;q
x
A
; x
B
;t and its time-reversed ver-
sion at the free surface. Since the Green’s function is
causal, it does not interfere with its time-reversed version
(except at t 0), so it is recovered by multiplying the
two-sided result with the Heaviside step function Ht.
Note that this Green’s function represents the velocity in
the x
p
direction at x
A
due to a source traction in the x
q
direction at x
B
. The terms under the integral on the right-
hand side represent cross correlations of particle veloc-
ities in the x
p
and x
q
directions at x
A
and x
B
, respectively,
due to volume force and deformation rate sources at x on
an arbitrarily shaped surface @D
1
in the medium.
Equation (4) is exact, so in theory it is possible to retrieve
the exact Green’s function, including the coda, of any
inhomogeneous anisotropic lossless medium (random or
deterministic) from cross correlations of observed parti-
cle velocities at the free surface. Unlike in the approach
for diffuse wave fields [1–7], which requires only one or a
few sources, in Eq. (4) it is assumed that volume force and
deformation rate sources are available everywhere on @D
1
and that the response of each source is measured sepa-
rately. In general not all sources are equally important.
When the aim is to retrieve surface waves, the main
contribution comes from sources at @D
1
close to the
free surface. On the other hand, in order to retrieve the
reflection response, sources at @D
1
below x
A
and x
B
give
the main contribution.
In the following we make some approximations which
circumvent several of the assumptions mentioned above.
When the medium outside @D
1
is homogeneous and
source-free, Eq. (4) may be approximated by
2Ref
^
G
v;t
p;q
x
A
; x
B
;!g
Z
@D
1
f
^
G
v;
p;k
x
A
; x;!g
^
G
v;
q;k
x
B
; x;!d
2
x; (5)
which follows from an analysis similar as in [19]. The
Green’s functions in the right-hand side represent again
velocities at x
A
and x
B
due to sources at x 2 @D
1
.The
superscript denotes that these sources are quasi P-wave
sources (for k 1) and quasi S-wave sources with differ-
ent polarizations (for k 2; 3). The accuracy of Eq. (5)
depends on the curvature of @D
1
(when @D
1
would be
planar, the only approximation would be that evanescent
waves are neglected). Next we assume mutually uncorre-
lated noise sources
^
N
k
x;! and
^
N
l
x
0
;! for any x and
x
0
at @D
1
, obeying the relation h
^
N
k
x;!
^
N
l
x
0
;!i
kl
x x
0
^
S!, where hi denotes a spatial ensemble
average and
^
S! the power spectrum of the noise (which
is assumed to be the same for all sources). Inserting this
relation in the right-hand side of Eq. (5) yields
2Ref
^
G
v;t
p;q
x
A
;x
B
;!g
^
S! hf
^
v
obs
p
x
A
;!g
^
v
obs
q
x
B
;!i;
(6)
where
^
v
obs
p
x
A
;!
Z
@D
1
^
G
v;
p;k
x
A
; x;!
^
N
k
x;!d
2
x (7)
and
^
v
obs
q
x
B
;!
Z
@D
1
^
G
v;
q;l
x
B
; x
0
;!
^
N
l
x
0
;!d
2
x
0
: (8)
Here
^
v
obs
p
x
A
;! and
^
v
obs
q
x
B
;! are the observed particle
velocities at x
A
and x
B
at the free surface due to a
distribution of noise sources at an arbitrarily shaped
surface @D
1
inside the medium. The average in Eq. (6)
is taken over different realizations of the source distribu-
tion. In the time domain Eq. (6) becomes
Z
1
1
fG
v;t
p;q
x
A
; x
B
; t
0
G
v;t
p;q
x
A
; x
B
;t
0
gSt t
0
dt
0
Z
1
1
v
obs
p
x
A
;t t
0
v
obs
q
x
B
;t
0
dt
0
: (9)
According to this equation, the cross correlation of the
observed particle velocities at x
A
and x
B
yields the elas-
todynamic Green’s function between x
A
and x
B
, con-
volved with the autocorrelation of the noise sources.
The advantage of Eq. (9) over Eq. (4) is that no separate
measurements of the responses of all sources at @D
1
are
required; these responses can be measured simulta-
neously, according to Eqs. (7) and (8). Note that by omit-
PRL 93, 254301 (2004)
PHYSICAL REVIEW LETTERS
week ending
17 DECEMBER 2004
254301-3
ting the subscripts p, q, k,andl in Eqs. (7)–(9) we obtain
an expression for the acoustic Green’s function at the free
surface of an inhomogeneous lossless fluid.
Equation (9) has a striking resemblance to the result of
[1–7]. In Eq. (9) there are no specific assumptions made
about the medium, but there should be many mutually
uncorrelated sources present, which leads again to a dif-
fuse wave field.
An application of Eq. (9) is passive seismic imaging. In
this application v
obs
p
x
A
;t and v
obs
q
x
B
;t represent pas-
sive recordings by geophones at the earth’s surface of
noise generated by natural sources in the subsurface. By
placing many geophones at the surface, the reflection
response can be reconstructed for many x
A
’s and x
B
’s.
By downward extrapolating these responses into the sub-
surface one can subsequently form an image of the sub-
surface [20–24]. This procedure makes no assumptions
with respect to the diffusivity of the wave field. Suppose
that the medium is fully deterministic and that there is
only one noise source present in the subsurface. Although
for this situation it is not possible to reconstruct the
Green’s function as an intermediate result, the passive
imaging method still maps the primary reflection re-
sponse (i.e., the ballistic wave) to its correct scattering
origin in depth as long as the specular reflection point at
the free surface lies within the array of geophones. This
has been shown with stationary phase analysis [20 –22]
and confirmed with numerical modeling studies [24].
In conclusion, we have shown that the elastodynamic
Green’s function of any inhomogeneous medium (random
or deterministic) can be retrieved from the cross correla-
tion of two recordings of a wave field at different receiver
locations. Unlike in the derivations in [1–7], which apply
to diffuse wave fields in random media or finite bodies,
we have made no assumptions about the distribution of
the medium parameters. Using a reciprocity theorem of
the correlation type, we derived an exact representation
for the Green’s function in terms of an integral of cross
correlations of observed wave fields at two points at the
free surface [Eq. (4)]; the integral is along an arbitrarily
shaped surface inside the medium, which contains the
sources of the observed wave fields. Next we assumed that
these sources are mutually uncorrelated noise sources,
which led to an expression [Eq. (9)] analogous to that in
[1–7]. Finally, we indicated that for passive imaging of
the ballistic wave the wave field need not be diffusive at
all: a single source in a deterministic medium suffices, as
long as the receiver array covers the specular reflection
point at the free surface [20 –22].
This work is supported by the Netherlands Research
Centre for Integrated Solid Earth Science (ISES). I thank
my colleagues Deyan Draganov, Evert Slob, and Jesper
Spetzler for carefully reading the manuscript.
*Electronic address: C.P.A.Wapenaar@CiTG.TUDelft.NL
[1] R. L. Weaver and O. I. Lobkis, Phys. Rev. Lett. 87, 134301
(2001).
[2] B. A. van Tiggelen, Phys. Rev. Lett. 91, 243904 (2003).
[3] R. K. Snieder, Phys. Rev. E 69, 046610 (2004).
[4] O. I. Lobkis and R. L. Weaver, J. Acoust. Soc. Am. 110,
3011 (2001).
[5] R. L. Weaver and O. I. Lobkis, Ultrasonics 40, 435 (2002).
[6] A. Derode, E. Larose, M. Tanter, J. de Rosny, A. Tourin,
M. Campillo, and M. Fink, J. Acoust. Soc. Am. 113, 2973
(2003).
[7] A. E. Malcolm, J. A. Scales, and B. A. van Tiggelen, Phys.
Rev. E 70, 015601 (2004).
[8] M. Campillo and A. Paul, Science 299, 547 (2003).
[9] N. M. Shapiro and M. Campillo, Geophys. Res. Lett. 31,
L07614 (2004).
[10] J. F. Claerbout, Geophysics 33, 264 (1968).
[11] M. R. Daneshvar, C. S. Clay, and M. K. Savage,
Geophysics 60, 1178 (1995).
[12] J. Rickett and J. Claerbout, Leading Edge 18, 957 (1999).
[13] J. Rickett and J. F. Claerbout, Stanford Exploration
Project Technical Report, 1996.
[14] T. L. Duvall, S. M. Jefferies, J.W. Harvey, and M. A.
Pomerantz, Nature (London) 362, 430 (1993).
[15] C. P. A. Wapenaar, D. Draganov, J.W. Thorbecke, and J.T.
Fokkema, in Expanded Abstracts of the 72nd Annual
Meeting of the Society of Exploration Geophysicists
(Society of Exploration Geophysicists, Tulsa, OK, 2002).
[16] C. P. A. Wapenaar, Geophysics 68, 1756 (2003).
[17] A. T. de Hoop, Appl. Sci. Res. 16 , 39 (1966).
[18] N. N. Bojarski, J. Acoust. Soc. Am. 74, 281 (1983).
[19] C. P. A. Wapenaar and G. C. Haime
´
, Geophys. Prospect.
38, 23 (1990).
[20] G. T. Schuster, in Extended Abstracts of the 63rd Annual
Meeting of the European Association of Geoscientists
and Engineers (European Association of Geoscientists
and Engineers, Houten, The Netherlands, 2001).
[21] G.T. Schuster, J. Yu, J. Sheng, and J. Rickett, Geophys. J.
Int. 157, 838 (2004).
[22] C. P. A. Wapenaar, D. Draganov, J. van der Neut, and J.W.
Thorbecke, in Expanded Abstracts of the 74th Annual
Meeting of the Society of Exploration Geophysicists
(Society of Exploration Geophysicists, Tulsa, OK, 2004).
[23] B. Artman, D. Draganov, C. P. A. Wapenaar, and B.
Biondi, in Extended Abstracts of the 66th Annual
Meeting of the European Association of Geoscientists
and Engineers (European Association of Geoscientists
and Engineers, Houten, The Netherlands, 2004).
[24] D. Draganov, C. P. A. Wapenaar, B. Artman, and B.
Biondi, in Expanded Abstracts of the 74th Annual
Meeting of the Society of Exploration Geophysicists
(Ref. [22]).
PRL 93, 254301 (2004)
PHYSICAL REVIEW LETTERS
week ending
17 DECEMBER 2004
254301-4