Edinburgh Research Explorer
Modeling of Wave Propagation in Inhomogeneous Media
Citation for published version:
van Manen, D-J, Robertsson , JOA & Curtis, A 2005, 'Modeling of Wave Propagation in Inhomogeneous
Media', Physical Review Letters, vol. 94, no. 16, 164301, pp. 1-4.
https://doi.org/10.1103/PhysRevLett.94.164301
Digital Object Identifier (DOI):
10.1103/PhysRevLett.94.164301
Link:
Link to publication record in Edinburgh Research Explorer
Document Version:
Publisher's PDF, also known as Version of record
Published In:
Physical Review Letters
Publisher Rights Statement:
Published in Physical Review Letters by the American Physical Society (2005)
General rights
Copyright for the publications made accessible via the Edinburgh Research Explorer is retained by the author(s)
and / or other copyright owners and it is a condition of accessing these publications that users recognise and
abide by the legal requirements associated with these rights.
Take down policy
The University of Edinburgh has made every reasonable effort to ensure that Edinburgh Research Explorer
content complies with UK legislation. If you believe that the public display of this file breaches copyright please
contact openaccess@ed.ac.uk providing details, and we will remove access to the work immediately and
investigate your claim.
Download date: 09. Aug. 2022
Modeling of Wave Propagation in Inhomogeneous Media
Dirk-Jan van Manen
*
and Johan O. A. Robertsson
WesternGeco Oslo Technology Centre, Solbraveien 23, 1383 Asker, Norway
Andrew Curtis
School of GeoSciences, University of Edinburgh, Grant Institute, West Mains Road, Edinburgh EH9 3JW, United Kingdom
(Received 8 February 2005; published 27 April 2005)
We present a methodology providing a new perspective on modeling and inversion of wave propagation
satisfying time-reversal invariance and reciprocity in generally inhomogeneous media. The approach
relies on a representation theorem of the wave equation to express the Green function between points in
the interior as an integral over the response in those points due to sources on a surface surrounding the
medium. Following a predictable initial computational effort, Green’s functions between arbitrary points
in the medium can be computed as needed using a simple cross-correlation algorithm.
DOI: 10.1103/PhysRevLett.94.164301 PACS numbers: 43.20.+g, 43.60.+d
Introduction.—Many applications in diverse fields such
as communications analysis, waveform inversion, imaging,
survey and experimental design, and industrial design
require a large number of modeled solutions of the wave
equation in different media. The most complete methods of
solution, such as finite differences (FD), which accurately
model all high-order interactions between scatterers in a
medium, typically become prohibitively expensive for
realistically complete descriptions of complex media and
the geometries of sources and receivers and hence for
solving realistic problems based on the wave equation.
Here we show that the key to breaking this apparent
paradigm lies in a basic reciprocity argument in combina-
tion with recent theoretical advances in the fields of time-
reversed acoustics [1] and seismic interferometry [2–5].
In time-reversed acoustics, invariance of the wave equa-
tion for time reversal can be exploited to focus a wave field
through a highly scattering medium on an original source
point [6]. Cassereau and Fink [7,8] realized that the acous-
tic representation theorem [9] can be used to time reverse
a wave field in a volume by creating secondary sources
(monopole and dipole) on a surface surrounding the me-
dium such that the boundary conditions correspond to the
time-reversed components of a wave field measured there.
These secondary sources give rise to the backpropagating,
time-reversed wave field inside the medium that collapses
onto itself at the original source location. Note that since
there is no source term absorbing the converging wave
field, the size of the focal spot is limited to half a (domi-
nant) wavelength in accordance with diffraction theory [7].
The diffraction limit was overcome experimentally by
de Rosny and Fink [10] by introducing the concept of an
‘‘acoustic sink.’’
In interferometry, waves recorded at two receiver loca-
tions are correlated to find the Green function between the
locations. Interferometry has been successfully applied to
helioseismology [11], ultrasonics [2], and exploration seis-
mics [3]. Recently it was shown that there exists a close
link between the time-reversed acoustics and interferom-
etry disciplines when Derode et al. [1] analyzed the emer-
gence of the Green function from field-field correlations in
an open scattering medium in terms of time-reversal sym-
metry. The Green function can be recovered as long as the
sources in the medium are distributed forming a perfect
time-reversal device. A rigorous proof for the general case
of an arbitrary inhomogeneous elastic medium was pre-
sented by Wapenaar [4].
Theory and method.—Our starting point is the acous-
tic wave equation in the space-frequency domain:
@
i
1
@
i
p!
2
=p f, where p px;! denotes
the pressure field at location x and frequency !, x,
and x denote the mass density and incompressibility,
respectively, and f fx;! is a source term, denoting the
change of volume injection rate density over time. Now
consider two states A and B that could occur in the same
medium independently: @
i
1
@
i
p
A
!
2
=p
A
f
A
and @
i
1
@
i
p
B
!
2
=p
B
f
B
. The acoustic repre-
sentation theorem can be derived by multiplying the equa-
tion for the first state by p
B
x;! and the equation for the
second state by p
A
x;!, subtracting and integrating the
results over a volume V, applying Gauss’s theorem to
convert the volume integral to a surface integral and iden-
tifying state A with a mathematical state [i.e., a state
involving (analytic) Green’s functions rather than mea-
sured quantities [9]]: f
A
xx x
0
and p
A
x
Gx; x
0
, where x denotes the Dirac delta distribution
and Gx; x
0
Green’s function due to a source at x
0
.
Following a reciprocity argument, interchanging the coor-
dinates x $ x
0
and dropping the superscripts for state B,
this procedure yields
px
Z
V
Gx; x
0
fx
0
dV
0
Z
S
1
x
0
r
0
Gx; x
0
px
0
Gx; x
0
r
0
px
0
ndS
0
; (1)
where r
0
Gx; x
0
denotes the gradient of the Green func-
PRL 94, 164301 (2005)
PHYSICAL REVIEW LETTERS
week ending
29 APRIL 2005
0031-9007=05=94(16)=164301(4)$23.00 164301-1 2005 The American Physical Society
tion with respect to primed coordinates and n denotes the
normal to the boundary. Thus, the wave field can be
computed everywhere inside the volume V once the source
f x
0
inside the volume and both the wave field px
0
and
its gradient r
0
px
0
on the surrounding surface S are
known. To time reverse a wave field in a volume V, the
wave field p and its gradient r
0
p, measured at the surface S
in a first step have to be time reversed on the surface such
that the time-reversed pressure field p
tr
x radiated from
the boundary can be written
p
tr
x
Z
S
1
x
0
r
0
Gx; x
0
p
x
0
Gx; x
0
r
0
p
x
0
ndS
0
; (2)
where an asterisk denotes complex conjugation, and we
have ignored the volume integral which corresponds to the
acoustic sink [10]. Note that Eq. (2) can be used to compute
the time-reversed wave field (including all high-order in-
teractions) at any location, not just at an original source
location. Now, assume that the wave field px
0
was due to
a point source at location x
1
and that we have px
0
Gx
0
; x
1
. By measuring the wave field in a second location
x
2
, the Green function and its time reverse between the
source point x
1
and the second point x
2
are observed [1,4]:
G
x
2
; x
1
Gx
2
; x
1
Z
S
1
x
0
r
0
Gx
2
; x
0
G
x
0
; x
1
Gx
2
; x
0
r
0
G
x
0
; x
1
ndS
0
;
(3)
where the negative forward Green function Gx
2
; x
1
arises from the missing acoustic sink [7,10]. Using reci-
procity, we can rewrite Eq. (3) so that it involves only
sources on the boundary enclosing the medium:
G
x
2
; x
1
Gx
2
; x
1
Z
S
1
x
0
r
0
Gx
2
; x
0
G
x
1
; x
0
Gx
2
; x
0
r
0
G
x
1
; x
0
ndS
0
:
(4)
Thus, the Green function between two points x
1
and x
2
can
be calculated once the Green functions between the enclos-
ing boundary and these points are known.
A highly efficient two-stage modeling strategy follows
from Eq. (4): first, the Green function terms G and r
0
G are
calculated from boundary locations to internal points in a
conventional forward modeling phase; in a second inter-
correlation phase, the integral is calculated requiring only
cross correlations and numerical integration. Since the
computational cost of typical forward modeling algorithms
(e.g., FD) does not significantly depend on the number of
receiver locations but mainly on the number of source
locations, efficiency and flexibility are achieved by storing
the wave field modeled for each of the boundary sources in
as many points as possible throughout the medium. To
calculate the Green function between two points the re-
cordings in the first point due to the dipole sources on the
boundary are cross correlated with the recordings in the
second point due to the monopole sources, and vice versa.
The resulting cross correlations are subtracted and numeri-
cally integrated over the boundary of source locations.
Unprecedented flexibility follows from the fact that the
Green function can be calculated between all pairs of
points that were defined up front and stored in the initial
modeling phase. Thus, we calculate a partial modeling
solution that is common to all Green functions, then a
bespoke component for each Green function.
Results.—Our method is illustrated using a FD imple-
mentation of the two-dimensional acoustic wave equation
for a typical modeling scenario in an exploration seismic
setting. In Fig. 1 the compressional wave velocity in a
4:6 4:6kmrepresentative region of an Earth model often
used to benchmark marine seismic imaging algorithms
[12] is shown. Note the high velocity (4500 m=s) salt
body on the right. In black, two points of interest (offset
1 km) are shown. The dotted line denotes the boundary
with N
S
912 source locations distributed with a density
consistent with the local spatial Nyquist frequency.
Outgoing (i.e., radiation or absorbing) boundary conditions
[13] are applied right outside the surface enclosing the
points of interest to truncate the computational domain.
Forward simulations were carried out for each of the
912 source locations on the boundary and the waveforms
stored at 90 000 points distributed throughout the model.
Note that because of the cross symmetry of the terms in the
integrand in Eq. (4), no sources are required along inter-
faces with homogeneous boundary conditions (e.g., the
Earth’s free surface). Depending on the particular wave
equation (scalar or vector), several forward simulations
may have to be carried out for each source location. In
the acoustic example two data sets are required: with
X−coordinate (m)
Depth (m)
1.9 1.95 2 2.05 2.1 2.15 2.2 2.25 2.3
x 10
4
0
500
1000
1500
2000
2500
3000
3500
4000
P−wave velocity (m/s)
1600
1800
2000
2200
2400
2600
2800
FIG. 1 (color). 2D acoustic marine seismic model (compres-
sional velocity) [12]. The color scale is clipped to display weak
velocity contrasts (velocity of salt is 4500 m=s).
PRL 94, 164301 (2005)
PHYSICAL REVIEW LETTERS
week ending
29 APRIL 2005
164301-2
monopole and dipole sources, respectively. However, when
the surface surrounding the medium has outgoing bound-
ary conditions, the wave field and its gradient (traction) are
directly related [14]. Hence, the normal derivatives can be
calculated from the wave field itself without additional
modeling. Figure 2 shows the integrand of Eq. (4), in-
versely weighted by boundary source density, for each
source location x
0
and for the two points x
1
and x
2
in
Fig. 1. In Fig. 3, the resulting Green function between
the points in Fig. 1 computed using Eq. (4), and a reference
trace computed by direct FD modeling, are shown in red
and blue, respectively. The signal at negative times corre-
sponds to the waves flowing back in time and the opposite
direction past the second point. The four insets show the
excellent match between the reference trace and the new
method in detail.
Interestingly, the time series in Fig. 2 bear little resem-
blance to the final Green function in Fig. 3. Equation (4)
sums these signals along the horizontal axis and hence
relies on the delicate constructive and destructive interfer-
ence of time-reversed waves backpropagating through the
medium, recombining and undoing the scattering at every
discontinuity to produce the Green function.
In Fig. 2, each column represents the set of all waves
traveling from point x
1
to a single boundary source, corre-
lated with the Green function from that boundary source
to x
2
. Some of the waves traveling from x
1
to this bound-
ary source may pass through x
2
before being recorded and
therefore have the remainder of their path in common with
waves emitted from x
2
in the same direction (or wave
number). The travel times associated with such identical
parts of the path are eliminated in the cross correlation and
the remaining part corresponds to an event in the Green
function from x
1
to x
2
. Similarly, some waves emitted from
x
2
may travel to the boundary source location via x
1
and
have a common section of path between x
1
and the bound-
ary source. Again travel time on the common section will
be eliminated and give rise to the same event in the Green
function from x
1
to x
2
at negative (acausal) times. Note
that the vector wave numbers involved for positive and
negative times are in general not parallel since they are
related to the background structure of the whole model
(one or the other may not exist for the same boundary
source). Hence, waves at positive and negative times are
Boundary source no.
Time (s)
100 200 300 400 500 600 700 800 900
−1
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
Pressure (kPa)
−0.1
−0.08
−0.06
−0.04
−0.02
0
0.02
0.04
0.06
0.08
0.1
FIG. 2 (color). Green’s function intercorrelation gather
(weighted) for the two points shown in Fig. 1. The low corre-
lation amplitude for boundary sources 250–310 corresponds to
the shadow of the salt body.
−5 0 5
−1.5
−1
−0.5
0
0.5
1
1.5
Time (s)
)aPk
( eruss
e
rP
0.5 1
−1
−0.5
0
0.5
1
2.5 3
−0.4
−0.2
0
0.2
0.4
−1 −0.5
−1
−0.5
0
0.5
1
−3 −2.5
−0.4
−0.2
0
0.2
0.4
0.4 0.6 0.8 1
0.5
1
2
4
8
Time (s)
(
gni
caps
ecruos yradnuoB λ
nim
)
A
B
FIG. 3 (color). (a) Waveform com-
puted by summation of the weighted
intercorrelation gather shown in Fig. 2
(red) compared to a conventional FD
computation (blue). The insets show par-
ticular events in the time series.
(b) Waveform computed after successive
subsampling of the intercorrelation
gather.
PRL 94, 164301 (2005)
PHYSICAL REVIEW LETTERS
week ending
29 APRIL 2005
164301-3
reconstructed differently. All cross correlations involv-
ing energy that does not pass through x
2
are eliminated
by destructive interference by the summation of the col-
umns [15].
The new method is particularly attractive in applications
where Green’s functions are desired between a large num-
ber of points interior to a medium, but where there are
few common sources or receiver points. No other exist-
ing method could offer full waveforms at comparable com-
putational cost. The method also offers great flexibility
where the exact interior points are not known in advance
since Green’s functions can be computed on an ‘‘as
needed’’ basis from Green’s functions between points on
the surrounding surface and its interior. We have shown
how the latter Green functions constitute a common com-
ponent of all Green’s functions in the medium through
Eq. (4). In the example above, this common component
is stored compressed by a factor of 50 compared to ex-
plicitly storing all desired Green’s functions between pairs
of interior points.
Whereas traditional approximate modeling methods
typically impose restrictions with respect to the degree of
heterogeneity in the medium of propagation or neglect
high-order scattering, the new time-reversal modeling
methodology allows us instead to compromise on the noise
level while maintaining high-order scattering and full het-
erogeneity in the medium. Recent experimental and theo-
retical work indicates that time-reversed imaging is robust
with respect to perturbations in the boundary conditions
[1,16]. For cases where the wave propagation is heavily
dominated by multiple scattering even a single source may
be sufficient to refocus the essential parts of a time-
reversed signal [17]. Also for more deterministic models,
such as the one in the example, it is possible to substan-
tially reduce the number of sources and still recover the
essential parts of the signal. In Fig. 3(b) we show the part of
the signal corresponding to the inset in the upper right
corner of Fig. 3(a) as we reduce the number of sources
around the boundary. Even for as few as
1
16
of the original
number of sources we are able to reproduce amplitude and
phase of an arrival of interest fairly accurately, but with an
increased noise level. Clearly, the required number of
sources will depend on the application. Our numerical
experiments thus confirm the robustness of the methodol-
ogy with respect to variations in parameters such as loca-
tion and discretization of integration surfaces.
We also experimented with exciting the boundary
sources simultaneously by encoding the source signals
using pseudonoise sequences [18] and with simultaneous
sources distributed randomly in the medium [1] as two
alternative ways to reduce the number of sources. There is
a well-known limit to the quality of separation of such
sequences of a given length when emitted simultaneously
[19]. Insufficient separation of sequences again is manifest
in an increased noise level in the final Green functions. In
all cases, the limits of separation caused relatively high
noise levels compared to the equivalent FD effort using the
method described above.
Thus, we have shown how recent insight into the rela-
tionship between Green’s theorem and time reversal can be
extended to the modeling of wave propagation by invoking
reciprocity. We expect that this may significantly change
the way we approach modeling and inversion of the wave
equation in future.
We would like to thank Roel Snieder for alerting us to
possible redundancy in the boundary sources, which was
tested using successive downsampling.
*Electronic address: DManen@oslo.westerngeco.slb.com
[1] A. Derode, E. Larose, M. Tanter, J. de Rosny, A. Tourin,
M. Campillo, and M. Fink, J. Acoust. Soc. Am. 113, 2973
(2003).
[2] R. L. Weaver and O. I. Lobkis, Phys. Rev. Lett. 87, 134301
(2001).
[3] K. Wapenaar, D. Draganov, J. Thorbecke, and J. Fokkema,
Geophys. J. Int. 156, 179 (2004).
[4] K. Wapenaar, Phys. Rev. Lett. 93, 254301 (2004).
[5] G. T. Schuster, in Extended Abstracts of the 63rd Annual
International Meeting of the European Association of
Geoscientists and Engineers (European Association of
Geoscientists and Engineers, Houten, The Netherlands,
2001), p. A-32.
[6] A. Derode, P. Roux, and M. Fink, Phys. Rev. Lett. 75,
4206 (1995).
[7] D. Cassereau and M. Fink, IEEE Trans. Ultrason.
Ferroelectr. Freq. Control 39, 579 (1992).
[8] D. Cassereau and M. Fink, J. Acoust. Soc. Am. 94, 2373
(1993).
[9] K. Wapenaar and J. Fokkema, J. Appl. Mech. 71, 145
(2004).
[10] J. de Rosny and M. Fink, Phys. Rev. Lett. 89, 124301
(2002).
[11] J. E. Rickett and J. F. Claerbout, Sol. Phys. 192, 203
(2000).
[12] D. Stoughton, J. Stefani, and S. Michell, in Expanded
Abstracts of the 71st Annual International Meeting of the
Society of Exploration Geophysicists (Society of
Exploration Geophysicists, Tulsa, OK, 2001), p. 1269.
[13] R. Clayton and B. Engquist, Bull. Seismol. Soc. Am. 67,
1529 (1977).
[14] E. Holvik, Ph.D. thesis, NTNU Trondheim, Norway, 2003.
[15] R. Snieder, Phys. Rev. E 69, 46610 (2004).
[16] R. Snieder and J. A. Scales, Phys. Rev. E 58, 5668 (1998).
[17] C. Draeger and M. Fink, J. Acoust. Soc. Am. 105, 611
(1999).
[18] P. Fan and M. Darnell, Sequence Design for Com-
munications Applications (Research Studies Press,
Hertfordshire, UK, 2003).
[19] L. R. Welch, IEEE Trans. Inf. Theory 20, 397 (1974).
PRL 94, 164301 (2005)
PHYSICAL REVIEW LETTERS
week ending
29 APRIL 2005
164301-4
![FIG. 1 (color). 2D acoustic marine seismic model (compressional velocity) [12]. The color scale is clipped to display weak velocity contrasts (velocity of salt is 4500 m=s).](/figures/fig-1-color-2d-acoustic-marine-seismic-model-compressional-3oz5slhx.webp)
