General representations for wavefield modeling and inversion in geophysics

General representations for wavefield modeling

and inversion in geophysics

Kees Wapenaar

1

ABSTRACT

Acoustic, electromagnetic, elastodynamic, poroelastic, and electroseismic waves are all governed by a unified ma-trix-vector wave equation. The matrices in this equation obey the same symmetry properties for each of these wave phe-nomena. This implies that the wave vectors for each of these phenomena obey the same reciprocity theorems. By substi-tuting Green’s matrices in these reciprocity theorems, unified wavefield representations are obtained. Analogous to the well-known acoustic wavefield representations, these unified representations find applications in geophysical modeling, migration, inversion, multiple elimination, and interfero-metry.

INTRODUCTION

Wavefield representations play an important role in forward and inverse geophysical problems, such as modeling, migration, inver-sion, multiple elimination, and, recently, interferometry. Various au-thors have derived acoustic and elastodynamic wavefield represen-tations by substituting Green’s functions into the Rayleigh and Ray-leigh-Betti reciprocity theorems, respectively共Morse and Feshbach, 1953; de Hoop, 1958; Gangi, 1970, 2000; Aki and Richards, 1980; Fokkema and van den Berg, 1993兲. In this paper, we follow a similar approach for a general matrix-vector wave equation that governs acoustic, electromagnetic, elastodynamic, poroelastic, or elec-troseismic wave propagation. First, we derive general convolution and correlation reciprocity theorems for this wave equation, supple-mented with boundary conditions for imperfectly coupled interfac-es. Next, we introduce a Green’s matrix as the point-source solution of the general wave equation. By substituting this Green’s matrix into the reciprocity theorems, we obtain general convolution- and correlation-type representations. We conclude this paper by briefly

discussing a number of applications of these general representations in seismic modeling, migration, inversion, multiple elimination, and interferometry.

MATRIX-VECTOR WAVE EQUATION

Diffusion, flow, and wave phenomena can each be captured by the differential matrix-vector equation,

ADtu + Bu + Dxu = s 共1兲

共Wapenaar and Fokkema, 2004兲, where u = u共x,t兲 is a vector con-taining space- and time-dependent field quantities, s = s共x,t兲 is a source vector, A = A共x兲 and B = B共x兲 are matrices containing space-dependent material parameters, and Dxis a matrix containing the spatial differential operators⳵1,⳵2, and⳵3. Dtdenotes the material

time derivative, defined as Dt=⳵t+ v0·⵱ =⳵t+vk0⳵k, where⳵tis

the time derivative in the reference frame and v0_{= v}0共x兲 the space-dependent flow velocity of the material;vk0denotes the kth

compo-nent of v0_{. Einstein’s summation convention applies to repeated} sub-scripts; lower-case Latin subscripts共except t兲 run from 1 to 3. In ex-ploration geophysics, we consider nonmoving media; hence, from here, onward we replace Dtby⳵t.

For acoustic wave propagation in an attenuating fluid, the vectors and matrices in equation 1 are defined by

u =

冢

p v₁ v2 v₃

冣

, s =

冢

q f₁ f2 f₃

冣

, 共2兲

with p = p共x,t兲 the acoustic pressure, vi=vi共x,t兲 the particle

veloci-ty, q = q共x,t兲 the volume injection rate, fi= fi共x,t兲 the external

force;

Manuscript received by the Editor December 16, 2006; revised manuscript received March 28, 2007; published online August 23, 2007. 1_{Delft University of Technology, Department of Geotechnology, Delft, the Netherlands. Email: c.p.a.wapenaar@tudelft.nl}

© 2007 Society of Exploration Geophysicists. All rights reserved.

A =

冢

␬

0 0 0 0

␳

0 0 0 0

␳

0 0 0 0

␳

冣

, B =

冢

bp 0 0 0 0 bv 0 0 0 0 bv 0 0 0 0 bv

冣

, 共3兲

with␬ = ␬共x兲 the compressibility, ␳ = ␳共x兲 the mass density, bp

= bp_{共x兲 and b}v_{= b}v_{共x兲 the loss terms, and}

Dx=

冢

0

⳵

1

⳵

2

⳵

3

⳵

1 0 0 0

⳵

2 0 0 0

⳵

3 0 0 0

冣

. 共4兲

For electromagnetic diffusion and/or wave propagation in matter, we have

u =

冉

E

H

冊

, s =

冉

− Je

− Jm

冊

, 共5兲

with E = E共x,t兲 and H = H共x,t兲 the electric and magnetic field vec-tors, Je_{= J}e共x,t兲 and Jm_{= J}m共x,t兲 the external electric and magnetic

current density vectors;

A =

冉

⑀

O

O

␮

冊

, B =

冉

␴

e

O

O

␴

m

冊

, 共6兲

with⑀ = ⑀共x兲 and ␮ = ␮共x兲 the permittivity and permeability ten-sors,␴e_=␴e_{共x兲 and␴}m_=␴m_{共x兲 the electric and magnetic}

conduc-tivity tensors, O the null-matrix, and

Dx=

冉

O D₀T D0 O

冊

, D0=

冢

0 −

⳵

3

⳵

2

⳵

3 0 −

⳵

1 −

⳵

2

⳵

1 0

冣

共7兲 共superscript T denotes matrix transposition only; it does not denote operator transposition兲.

For elastodynamic wave propagation in a solid, we have uT

=共vT_,−␶ 1 T_,−␶ 2 T_,−␶ 3 T兲 共with v and␶

ithe particle velocity and

trac-tion vectors兲, sT_=共fT_,h 1 T_,h 2 T_,h 3 T_{兲 共with f and h}

ithe external force and

deformation rate vectors兲, and A, B, and Dxare 12⫻12 matrices 共de Hoop and de Hoop, 2000; Wapenaar and Fokkema, 2004兲.

For electroseismic wave propagation in a saturated porous solid 共Pride, 1994兲, we have uT₌共ET_,HT_,兵vs其T_,−␶

1 T_,−␶ 2 T_,−␶ 3 T_,wT_{, p}f兲

共with superscripts s and f referring to the solid and fluid phase, re-spectively兲, where E and H are the average electric and magnetic field vectors, vs_and␶

ithe solid particle velocity and bulk traction

vectors, w =␸共vf_{− v}s兲 the filtration velocity 共with␸ the porosity兲,

and pf _{the fluid pressure. Moreover, s}T_{=共− 兵J}e_其T_,

−兵Jm其T_,fT_,h 1 T_,h 2 T_,h 3

T_,兵ff其T_,hf兲, where Je_{and J}m_{are again the external}

electric and magnetic current density vectors, f and ff_{the external}

forces on the bulk and on the fluid, and hiand hfthe modified

exter-nal deformation rates for the bulk and the fluid. Fiexter-nally, for this situa-tion A, B, and Dxare 22⫻22 matrices. Omitting E, H, Je, and Jm from u and s gives the field and source vectors for poroelastic waves 共Biot, 1956兲 governed by 16⫻16 matrices A, B, and Dx. On the oth-er hand, omitting w, pf_{, f}f_{and h}f_{and reorganizing B results in the}

electrokinetic equations for a piezoelectric system共Auld, 1973兲, with 18⫻18 matrices A, B and Dx.

In all cases, matrices A共x兲 and B共x兲 can be replaced by convolu-tional operators A共x,t兲ⴱ and B共x,t兲ⴱ to account for more general at-tenuation mechanisms.

We define the temporal Fourier transform of a space- and time-de-pendent quantity p共x,t兲 as

pˆ共x,

␻

兲 =

冕

−⬁ ⬁

exp共− j

␻

t兲p共x,t兲dt, 共8兲

where j is the imaginary unit and␻ the angular frequency. Applying the Fourier transform to all terms in matrix-vector equation 1, with Dtreplaced by⳵tand A共x兲 and B共x兲 replaced by convolutional

oper-ators A共x,t兲ⴱ and B共x,t兲ⴱ, we obtain

j

␻

Aˆ uˆ + Bˆuˆ + Dxuˆ = sˆ, 共9兲

where, apart from the field and source vectors uˆ共x,␻兲 and sˆ共x,␻兲, the material parameter matrices Aˆ 共x,␻兲 and Bˆ共x,␻兲 are in their general form frequency-dependent and complex-valued. Note that j␻Aˆ and Bˆ could be combined into one material-parameter matrix. However, we prefer to keep these terms separated to acknowledge the different character of the parameters in Aˆ and Bˆ.

For each situation, matrices Aˆ , Bˆ, and Dxobey the symmetry rela-tions KAˆ K = Aˆ = AˆT_{, 共10兲} KBˆ K = BˆT, 共11兲 KDxK = − Dx= − Dx T , 共12兲

where K is a real-valued diagonal matrix, obeying the property K = K−1_{. For example, for the acoustic and electromagnetic situations} discussed above, we have

K =

冢

1 0 0 0 0 − 1 0 0 0 0 − 1 0 0 0 0 − 1

冣

and K =

冉

− I O O I

冊

, 共13兲

respectively, with I the identity matrix.

MATRIX-VECTOR BOUNDARY CONDITION

At any position in space where the medium parameters in matrices Aˆ and Bˆ are discontinuous, the wavefield quantities in vector uˆ should obey boundary conditions. The same is true at a fracture with imperfect coupling. For this situation, the wavefield quantities in uˆ may exhibit a finite jump. We call both types of medium singularities interfaces. In the following, we consider the most general case for which the medium parameters are different at both sides of the inter-face and the media are not in perfect contact with each other.

Consider an interface with normal vector n =共n1,n2,n3兲Tbetween two materials共see Figure 1兲. In linearized form, the boundary condi-tions at an imperfect interface can be formulated in the space-fre-quency domain as

关Muˆ兴 = − j

␻

Yˆ 具Muˆ典 共14兲

condi-tions, Yˆ = Yˆ共x,␻兲 is a matrix containing the boundary parameters, and关·兴 and 具·典 represent the jump and the average across the inter-face, respectively, as stated by

关pˆ共x,

␻

兲兴 = lim

h↓0共pˆ共x + hn,

␻

兲 − pˆ共x − hn,

␻

兲兲, 共15兲

具pˆ共x,

␻

兲典 = lim

h↓0共pˆ共x + hn,

␻

兲 + pˆ共x − hn,

␻

兲兲/2, 共16兲

where x is chosen at the interface.

For acoustic waves, the matrices in equation 14 are defined as

M =

冉

1 0 0 0

0 n1 n2 n3

冊

and Yˆ =

冉

0

␳

ˆ

b

␬

ˆb ₀

冊

. 共17兲

The superscript b denotes that␳ˆb₌␳ˆb共x,␻兲 and␬ˆb₌␬ˆb共x,␻兲 are boundary parameters. The dimension of each boundary parameter is equal to the dimension of the corresponding volumetric parameter, multiplied by meter. For vanishing␳ˆb_and␬ˆb_{, equation 14, with uˆ, M,} and Yˆ defined in equations 2 and 17, reduces to the standard bound-ary conditions for perfectly coupled fluids, i.e.,关pˆ兴 = 0 and 关vˆini兴

= 0. When␳ˆb_and␬ˆb_{are nonzero, 1/j}␻␳ˆb_{is the hydraulic boundary} permeability,␬ˆb_{the boundary compressibility, and 1/}␬ˆb_{= K}ˆb _the boundary stiffness. Note that the dimension of the boundary stiffness Kˆb_{is that of stiffness per meter共i.e., Pa/m兲. Therefore, Kˆ}b_{is also} called the specific boundary stiffness.

For elastodynamic waves, the matrices in equation 14 are defined as

M =

冉

I O O O

O n₁I n₂I n₃I

冊

and Yˆ =

冉

O Sˆb

␳

ˆb _O

冊

, 共18兲

where␳ˆb_{and Sˆ}b_{are the boundary density and compliance tensors,} respectively. When␳ˆb_{= O}共which is usually a good approximation兲 and n =共0,0,1兲T_{共i.e., the interface is horizontal兲, equation 14, with}

M and Yˆ defined in equation 18, reduces to the linear slip model of Schoenberg共1980兲 when Sˆb_{is diagonal and real-valued, to the} ex-tended linear slip model of Pyrak-Nolte et al.共1990兲 when Sˆb_is diag-onal and complex-valued, or to the general boundary model 共includ-ing shear-induced conversion兲 of Nakagawa et al. 共2000兲 when the nondiagonal elements of Sˆb_{are also nonzero. Liu et al.}共1995, 2000兲 relate the parameters in the compliance tensor to the details of the microstructure of the interface.

For electromagnetic waves, equation 14 is a generalization of the Kaufman and Keller共1983兲 conductive interface model. For po-roelastic waves, it is a generalization of the Gurevich and Schoen-berg共1999兲 permeable interface model; when Yˆ vanishes, it reduces to the Deresiewicz and Skalak共1963兲 open-pore boundary condition

for the perfectly coupled porous solids. For electroseismic waves, equation 14 combines the boundary conditions for electromagnetic and poroelastic waves.

In all cases Yˆ obeys the symmetry relations

YˆTN = − NYˆ and Yˆ†J = JYˆ* 共19兲

共Wapenaar et al., 2004兲, where superscript ⴱ denotes complex conju-gation and † complex conjuconju-gation and transposition. For example, for the acoustic and elastodynamic situations discussed above, we have N =

冉

0 1 − 1 0

冊

, J =

冉

0 1 1 0

冊

共20兲 and N =

冉

O I − I O

冊

, J =

冉

O I I O

冊

, 共21兲 respectively.

CONVOLUTION-TYPE RECIPROCITY THEOREM

In general, a reciprocity theorem interrelates two independent states in one and the same domain共de Hoop, 1966; Fokkema and van de Berg, 1993兲. Here, we derive a reciprocity theorem for the general wave vector uˆ described in the previous sections. We introduce two independent states共i.e., wavefields, medium parameters, boundary parameters, and source functions兲 that will be distinguished by the subscripts A and B,共see Table 1兲. In the frequency domain, each of these states obeys the general matrix-vector-wave equation 9. We consider the interaction quantity uˆATKDxuˆB−共DxuˆA兲TKuˆB. Using

wave equation 9 as well as symmetry relations 10 and 11 for both states, we obtain

uˆA T

KDxuˆB−共DxuˆA兲TKuˆB= uˆA T KsˆB− sˆA T KuˆB − uˆA T K兵j

␻

共AˆB− AˆA兲 +共BˆB− BˆA兲其uˆB. 共22兲

This is the local form of the convolution-type matrix-vector reci-procity theorem. We call this convolution type because the products in the frequency domain共uˆATKsˆBetc.兲 correspond to convolutions in

the time domain. Next, we consider an arbitrary spatial domainD with boundary⳵D and outward-pointing normal vector n 共see Figure 2兲. Note that⳵D does not necessarily coincide with a physical bound-ary. For the moment, we assume that the medium parameters in both

n

Figure 1. Interface between two media with imperfect coupling.

Table 1. States for the unified reciprocity theorems.

State A State B

Wavefields uˆA共x,␻兲 uˆB共x,␻兲

Medium parameters _兵Aˆ

A,BˆA其共x,␻兲 兵AˆB,BˆB其共x,␻兲

Boundary parameters _Yˆ_{A共x,␻兲 Y}ˆ_B共x,␻兲 Source functions sˆA共x,␻兲 sˆB共x,␻兲

states are continuous inD. We integrate both sides of equation 22 over this domain and apply Gauss’ theorem in matrix-vector form 共equation A-5, seeAppendix A兲 to the left-hand side. This yields

冖

⳵D uˆA T KNxuˆBd2x =

冕

D 兵uˆA T KsˆB − sˆA T KuˆB其d3x −

冕

D uˆA T K兵j

␻

共AˆB− AˆA兲 +共BˆB− BˆA兲其uˆBd3x. 共23兲

Here, Nxis a matrix containing the components of normal vector n, organized in the same way as matrix Dx. For example, for the acous-tic situation Nxis defined as

Nx=

冢

0 n₁ n₂ n₃ n1 0 0 0 n2 0 0 0 n3 0 0 0

冣

. 共24兲

For all situations, Nxobeys the symmetry relation

KNxK = − Nx = − Nx T

, 共25兲

analogous to equation 12. Equation 23 is the unified convolution-type reciprocity theorem for a domainD in which the medium pa-rameters are continuous. It interrelates the wavefield quantities 共con-tained in uˆAand uˆB兲, the medium parameters 共contained in AˆA, AˆB,

BˆA, and BˆB兲 and the source functions 共contained in sˆAand sˆB兲 of states

A and B. The left-hand side is a boundary integral that contains a spe-cific combination of the wavefield quantities of states A and B at the boundary of domainD. The first integral on the right-hand side inter-relates the wavefield quantities and the source functions inD. The second integral on the right-hand side contains the differences be-tween the medium parameters in both states; obviously, this integral vanishes when the medium parameters in both states are identical.

We now extend the reciprocity-theorem equation 23 for the situa-tion in whichD contains internal interfaces共or fractures兲 with imper-fect coupling. To this end, we subdivideD into M continuous re-gions, according toD = D1艛D2· · · ·艛DM, see Figure 3. RegionDm

is enclosed by surface⳵Dmwith outward-pointing normal vector nm.

The boundaries between these regions represent the imperfect inter-nal interfaces. Note that each interinter-nal interface is part of two surfaces ⳵Dm, with opposite-pointing normal vectors nm, see Figure 3.

Since the medium parameters in regionDmare continuous,

reci-procity-theorem equation 23 applies to each of these regions. Sum-ming both sides of this equation over m again yields equation 23 for the total domainD, with an extra integral over the internal interfaces on the left-hand side as stated by

冕

⳵Dint 共共uˆA T KNxuˆB兲1+ 共uˆA T KNxuˆB兲2兲d2x, 共26兲

where⳵Dintconstitutes the total of all internal interfaces; the sub-scripts 1 and 2 denote the two sides of the internal interfaces. Using the general boundary condition 14 for imperfect interfaces and the first of the symmetry relations in equation 19, this internal interface integral can be rewritten as

冕

⳵Dint uˆA T MTN共I − ZˆA −1 ZˆB兲MuˆBd2x 共27兲

共Wapenaar et al., 2004兲 with

Zˆ = 共I + j

␻

Yˆ /2兲−1_{共I − j}

_␻

_{Yˆ /2兲. 共28兲}

Zˆ obeys the symmetry relation

ZˆTN = NZˆ−1, 共29兲

which follows from equations 19 and 28. For small Yˆ , equation 27 simplifies to j

␻

冕

⳵Dint uˆA T MTN共YˆB− YˆA兲MuˆBd2x. 共30兲

In the integrals in equations 27 and 30, uˆA, uˆB, and M are all chosen at

the same side of the interfaces共but which side is arbitrary兲.

Adding the internal interface integral of equation 27 to the left-hand side of equation 23, we obtain

冖

⳵D uˆA T KNxuˆBd2x +

冕

⳵Dint uˆA T MT_{N共I − Zˆ} A −1_Zˆ B兲MuˆBd2x =

冕

D 兵uˆA T KsˆB− sˆA T KuˆB其d3x −

冕

D uˆA T K兵j

␻

共AˆB− AˆA兲 +共BˆB− BˆA兲其uˆBd3x. 共31兲 n x1 x3 x2 ∂

Figure 2. Configuration for the reciprocity theorems.

n3 n1 n1 n1 n2 n2 n2 _n 3 n3 ∂ 3 ∂ 2 ∂ 1 3 1 2

We conclude this section by specifying equation 31 for the acoustic situation, assuming small boundary parameters. Upon substitution of equations 2, 3, 13, 17, 20, and 24, we obtain

冖

⳵D共pˆ

Avˆi,B− vˆi,ApˆB兲nid2x + j

␻

冕

⳵Dint

兵共

␬

ˆ_Bb −

␬

ˆ_Ab兲pˆApˆB

−共

␳

ˆ_Bb −

␳

ˆ_Ab兲vˆ_i,Anivˆj,Bnj其d2x

=

冕

D

兵pˆAqˆB −vˆi,Afˆi,B− qˆApˆB+ fˆi,Avˆi,B其d3x

− j

␻

冕

D 兵共

␬

ˆB−

␬

ˆA兲pˆApˆB− 共

␳

ˆB−

␳

ˆA兲vˆi,Avˆi,B其d3x −

冕

D 兵共bˆB p − bˆA p_兲pˆ ApˆB−共bˆBv − bˆvA兲vˆi,Avˆi,B其d3x, 共32兲

which has the familiar form known from, e.g., Fokkema and van den Berg共1993兲 and de Hoop 共1995兲, but with an extra integral over the internal interfaces.

CORRELATION-TYPE RECIPROCITY THEOREM

Porter共1970兲 and Bojarski 共1983兲 formulated reciprocity rems with back-propagating wavefields. Here, we extend these theo-rems for the general wave vector uˆ. We consider the interaction quantity uˆA†DxuˆB+共DxuˆA兲†uˆB. Since superscript † denotes

transpo-sition and complex conjugation, uˆA†represents a back-propagating

wavefield. Using wave equation 9 as well as the symmetry relations in equations 10 and 11 for both states, we obtain

uˆA†DxuˆB +共DxuˆA兲†uˆB = uˆA†sˆB + sˆA†uˆB − uˆA†兵j

␻

共AˆB− AˆA*兲

+共BˆB+ BˆA

†_兲其uˆ

B. 共33兲

This is the local form of the correlation-type matrix-vector reciproc-ity theorem. We say correlation type because the products in the fre-quency domain共uˆA†sˆBetc.兲 correspond to correlations in the time

do-main. Next, we consider again domainD with boundary⳵D and out-ward-pointing normal vector n, see Figure 2. For the moment, we as-sume that the medium parameters in both states are continuous inD. We integrate both sides of equation 33 over this domain and apply the theorem of Gauss in matrix-vector form共equation A-4兲 to the left-hand side. This yields

冖

⳵DuˆA † NxuˆBd2x =

冕

D 兵uˆA † sˆB + sˆA † uˆB其d3x −

冕

D uˆA †_兵j

_␻

_共Aˆ B− AˆA *_兲 +共BˆB+ BˆA†兲其uˆBd3x. 共34兲

Equation 34 is the unified correlation-type reciprocity theorem for a domainD in which the medium parameters are continuous.

We now extend this reciprocity theorem for the situation in which D contains internal interfaces with imperfect coupling, see Figure 3. Since the medium parameters in regionDmare continuous,

reciproc-ity theorem 34 applies to each of these regions. Summing both sides of this equation over m again yields equation 34 for the total domain

D, with an extra integral over the internal interfaces on the left-hand side, as stated by

冕

⳵Dint 共共uˆA † NxuˆB兲1+ 共uˆA † NxuˆB兲2兲d2x. 共35兲

Using the general boundary condition 14 for imperfect interfaces and the second of the symmetry relations in equation 19, this internal interface integral can be rewritten as

冕

⳵Dint

uˆA

†

MTJ共I − 共ZˆA

⬘

兲−1ZˆB兲MuˆBd2x 共36兲

共Wapenaar et al., 2004兲, with

Zˆ

⬘

=共I + j

␻

Yˆ*/2兲−1共I − j

␻

Yˆ*/2兲. 共37兲

Note that Zˆ_⬘obeys the symmetry relation

Zˆ†J = J共Zˆ

⬘

兲−1, 共38兲

which follows from equations 19 and 37. For small Yˆ , equation 36 simplifies to

j

␻

冕

⳵Dint

uˆA†MTJ共YˆB− YˆA*兲MuˆBd2x. 共39兲

Adding the internal interface integral of equation 36 to the left-hand side of equation 34, we obtain

冖

⳵D uˆA†NxuˆBd2x +

冕

⳵Dint uˆA†M T J共I − 共ZˆA

⬘

兲−1ZˆB兲MuˆBd2x =

冕

D 兵uˆA†sˆB + sˆA†uˆB其d3x −

冕

D uˆA†兵j

␻

共AˆB− AˆA*兲 +共BˆB+ BˆA †_兲其uˆ Bd3x. 共40兲

Note that when the medium and boundary parameters, sources, and wavefields are identical in both states, this reciprocity theorem reduces共omitting the subscripts A and B兲 to

2R

冕

D uˆ†sˆd3x =

冖

⳵D uˆ†Nxuˆd2x +

冕

D

uˆ†兵− 2

␻

I共Aˆ兲 + Bˆ + Bˆ†其uˆd3x

+

冕

⳵Dint

uˆ†_MT_{J共I − 共Zˆ}

_⬘

_兲−1_{Zˆ 兲Muˆd}2_{x, 共41兲}

We conclude this section by specifying equation 40 for the acous-tic situation, assuming small boundary parameters. Upon substitu-tion of equasubstitu-tions 2, 3, 17, 20, and 24, we obtain

冖

⳵D共pˆA * vˆi,B+ vˆi,A * pˆB兲nid2x + j

␻

冕

⳵Dint 兵共

␬

ˆ_Bb −

␬

ˆ_Ab*兲pˆ_A*pˆB +共

␳

ˆ_Bb −

␳

ˆb*_A兲vˆ_i,A* nivˆj,Bnj其d2x =

冕

D 兵pˆA * qˆB +vˆ_i,A* fˆi,B+ qˆA * pˆB+ fˆ_i,A* vˆi,B其d3x − j

␻

冕

D 兵共

␬

ˆB−

␬

ˆA*兲pˆA*pˆB+共

␳

ˆB−

␳

ˆA*兲vˆi,A* vˆi,B其d3x −

冕

D 兵共bˆB p + bˆA p*_兲pˆ A * pˆB+共bˆB v + bˆA v*_兲vˆ i,A * vˆi,B其d3x, 共42兲 which has the familiar form known from, e.g., Fokkema and van den Berg共1993兲 and de Hoop 共1995兲, but with an extra integral over the internal interfaces.

GREEN’S MATRIX

The wavefield vector uˆ共x,␻兲 and the source vector sˆ共x,␻兲 are L ⫻1 vectors, where the value of L depends on the type of wavefield considered. A Green’s function is defined as the wavefield that would be obtained if the source were an impulsive point source␦共x − x⬘兲␦共t兲, or, in the frequency domain, a point source␦共x − x⬘兲 with unit spectrum. Since the source vector sˆ contains L different source functions, we may define L different Green’s wavefield vectors. We define the lth Green’s wavefield vector共with 1ⱕlⱕL兲 as the causal solution of general wave equation 9 with boundary condition 14, with source vector sˆ共x,␻兲 replaced by il␦共x − x⬘兲, where ilis the L

⫻1 unit vector 共0, ¯,1, ¯,0兲T_{, with ‘1’ on the lth position. Hence,}

in the space-frequency domain the lth Green’s wave vector obeys the relations

j

␻

Aˆ gˆl+ Bˆ gˆl+ Dxgˆl= il

␦

共x − x

⬘

兲 共43兲

and

关Mgˆl兴 = − j

␻

Yˆ 具Mgˆl典, 共44兲

where gˆl= gˆl共x,x⬘,␻兲 is the lth L⫻1 Green’s wave vector observed

at x, due to a point source of the lth type at x_⬘. Due to their causal be-haviour in the time domain, the components of these Green’s vectors obey the Kramers-Kronig relations.

Equations 43 and 44 each represent L matrix-vector equations for the L Green’s wave vectors gˆl, with 1ⱕlⱕL. For example, for the

acoustic situation共L = 4兲, equation 43 reads

冢

␩

ˆ

⳵

₁

⳵

₂

⳵

₃

⳵

1

␨

ˆ 0 0

⳵

2 0

␨

ˆ 0

⳵

3 0 0

␨

ˆ

冣

冢

Gˆp,q_共x,x

_⬘

_,

_␻

_兲 Gˆ₁v,q共x,x

⬘

,

␻

兲 Gˆ₂v,q共x,x

⬘

,

␻

兲 Gˆ₃v,q共x,x

⬘

,

␻

兲

冣

=

冢

␦

共x − x

⬘

兲 0 0 0

冣

, 共45兲 for l = 1共with␩ˆ = j␻␬ˆ + bˆp_and␨ˆ = j␻␳ˆ + bˆv兲,

冢

␩

ˆ

⳵

₁

⳵

₂

⳵

₃

⳵

1

␨

ˆ 0 0

⳵

2 0

␨

ˆ 0

⳵

3 0 0

␨

ˆ

冣

冢

Gˆ,1p,f共x,x

⬘

,

␻

兲 Gˆ_1,1v,f共x,x

⬘

,

␻

兲 Gˆ_2,1v,f共x,x

⬘

,

␻

兲 Gˆ_3,1v,f共x,x

⬘

,

␻

兲

冣

=

冢

0

␦

共x − x

⬘

兲 0 0

冣

, 共46兲 for l = 2, etc. The superscripts of the Green’s functions refer to the type of observed wavefield at x and the source type at x_⬘, respective-ly; the subscripts denote the different components. We now combine the L Green’s vectors into a Green’s matrix and the L source vectors into a source matrix,

共gˆ1. . . gˆl. . . gˆL兲共x,x

⬘

,

␻

兲 = Gˆ共x,x

⬘

,

␻

兲, 共47兲

共i1. . . il. . . iL兲

␦

共x − x

⬘

兲 = I

␦

共x − x

⬘

兲, 共48兲

where Gˆ 共x,x⬘,␻兲 is the L⫻L Green’s wavefield matrix and I is the L⫻L identity matrix. With this notation, equations 43 and 44 for l = 1 . . . L can be combined into

j

␻

Aˆ Gˆ + BˆGˆ + DxGˆ = I

␦

共x − x

⬘

兲 共49兲

and

关MGˆ兴 = − j

␻

Yˆ 具MGˆ典, 共50兲

respectively. These are the general wave equation and boundary con-dition for the Green’s matrix Gˆ 共x,x_⬘,␻兲.

Note that the wave vector uˆ共x,␻兲 and the Green’s matrix Gˆ 共x,x_⬘,␻兲 obey the same linear wave equation with the same linear boundary conditions, but with different source functions, sˆ共x,␻兲 and I␦共x − x⬘兲, respectively. Hence, we may apply the superposition principle to express the wave vector as

uˆ共x,

␻

兲 =

冕

Ds

Gˆ 共x,x

⬘

,

␻

兲sˆ共x

⬘

,

␻

兲d3x

⬘

, 共51兲 whereDsis the domain occupied by the source distribution. In the

following, we derive representations for the Green’s matrix Gˆ 共x,x_⬘,␻兲, and, implicitly via equation 51, for the wave vector uˆ共x,␻兲.

CONVOLUTION-TYPE REPRESENTATION

We consider again the piecewise continuous domain D with boundary⳵D and outward-pointing normal vector n 共Figure 3兲. We assume that the boundaries between the different regions are imper-fect interfaces; the combination of all internal interfaces is represent-ed by⳵Dint.

B. The medium parameters in this state are the actual parameters, and, consequently, the Green’s matrix in state B is defined in the ac-tual medium, again see Table 2.

Consider the following property of the delta function

冕

D

␦

共x − x

⬘

兲u共x兲d3_{x =}

_␹

D共x

⬘

兲u共x

⬘

兲, 共52兲

where␹D共x⬘兲 is the characteristic function for domain D, defined as

␹

D共x

⬘

兲 =

冦

1 for x

⬘

苸 D 1 2 for x

⬘

苸

⳵

D 0 for x

⬘

苸 R3\兵D 艛

⳵

D其.

冧

共53兲

Upon substitution of the states of Table 2 into the convolution-type reciprocity theorem共equation 31兲, using this property of the delta function, we obtain

␹

D共x

⬙

兲KG¯ˆT共x

⬙

,x

⬘

,

␻

兲K −

␹

D共x

⬘

兲Gˆ共x

⬘

,x

⬙

,

␻

兲 =

冖

⳵DKG ¯ˆT_共x,x

_⬘

,

␻

兲KNxGˆ 共x,x

⬙

,

␻

兲d2x +

冕

D KG¯ˆT共x,x

⬘

,

␻

兲K⌬Hˆ共x,

␻

兲Gˆ共x,x

⬙

,

␻

兲d3x +

冕

⳵Dint KG¯ˆT共x,x

⬘

,

␻

兲K⌬Hˆb共x,

␻

兲Gˆ共x,x

⬙

,

␻

兲d2x, 共54兲 with the contrast functions⌬Hˆ and ⌬Hˆb_{defined as}

⌬Hˆ = j

␻

共Aˆ − A¯ˆ 兲 + 共Bˆ − B¯ˆ 兲, 共55兲 ⌬Hˆb_{= KM}T

N兵I − Z¯ˆ−1Zˆ 其M. 共56兲

Equation 54 is the general convolution-type representation of the Green’s matrix. Applications are discussed in a later section. Here, we derive a reciprocity relation for the Green’s matrix. To this end we replace the background parameters in state A by the actual pa-rameters; hence, the last two integrals on the right-hand side of equa-tion 54 vanish. Then, we replaceD by R3_{, to make the characteristic} functions in the left-hand side of equation 54 both equal 1. Finally, we assume that outside some sphere with finite radius, the medium is homogeneous, isotropic, and nonporous, which implies that the first integral on the right-hand side vanishes as well共Sommerfeld radia-tion condiradia-tions, Born and Wolf, 1965; Pao and Varatharajulu, 1976; de Hoop, 1995兲. This leaves

KGˆT共x

⬙

,x

⬘

,

␻

兲K = Gˆ共x

⬘

,x

⬙

,

␻

兲. 共57兲 Of course a similar relation holds for the Green’s matrix in the back-ground medium. Equation 57 formulates source-receiver reciprocity for a piecewise continuous medium with imperfect interfaces. For example, for the acoustic situation, we have

冢

Gˆp,q − Gˆ1v,q − Gˆ2v,q − Gˆ3v,q − Gˆ_,1p,f Gˆ_1,1v,f Gˆ_2,1v,f Gˆ_3,1v,f − Gˆ_,2p,f Gˆ_1,2v,f Gˆ_2,2v,f Gˆ_3,2v,f − Gˆ_,3p,f Gˆ_1,3v,f Gˆ_2,3v,f Gˆ_3,3v,f

冣

共x

⬙

,x

⬘

,

␻

兲 =

冢

Gˆp,q Gˆ_,1p,f Gˆ_,2p,f Gˆ_,3p,f Gˆ1v,q Gˆ1,1v,f Gˆ1,2v,f Gˆ1,3v,f Gˆ₂v,q Gˆ_2,1v,f Gˆ_2,2v,f Gˆ_2,3v,f Gˆ₃v,q Gˆ_3,1v,f Gˆ_3,2v,f Gˆ_3,3v,f

冣

共x

⬘

,x

⬙

,

␻

兲. 共58兲 CORRELATION-TYPE REPRESENTATION

For the derivation of a correlation-type representation of the Green’s matrix, we substitute the states of Table 2 into the reciproci-ty theorem of correlation-reciproci-type共equation 40兲, which gives

␹

D共x

⬙

兲G¯ˆ†共x

⬙

,x

⬘

,

␻

兲 +

␹

D共x

⬘

兲Gˆ共x

⬘

,x

⬙

,

␻

兲 =

冖

⳵D G¯ˆ†共x,x

⬘

,

␻

兲NxGˆ 共x,x

⬙

,

␻

兲d2x +

冕

D G¯ˆ†共x,x

⬘

,

␻

兲⌬Hˆ共x,

␻

兲Gˆ共x,x

⬙

,

␻

兲d3x +

冕

⳵Dint G¯ˆ†共x,x

⬘

,

␻

兲⌬Hˆb共x,

␻

兲Gˆ共x,x

⬙

,

␻

兲d2x, 共59兲

with the contrast functions⌬Hˆ and ⌬Hˆb_{now defined as}

⌬Hˆ = j

␻

共Aˆ − A¯ˆ*_{兲 + 共Bˆ + B¯}ˆ†_{兲, 共60兲}

⌬Hˆb_{= M}T

J兵I − 共Z¯ˆ

⬘

兲−1Zˆ 其M. 共61兲

Equation 59 is the general correlation-type representation of the Green’s matrix.

APPLICATIONS

Here we discuss a number of applications of the general convolu-tion-type and correlaconvolu-tion-type representations. This overview is not exhaustive but serves as an illustration.

Table 2. Green’s states for the unified representations.

Forward wavefield extrapolation

The most straightforward application of the convolution-type rep-resentation共equation 54兲 is forward wavefield extrapolation. Con-sider the configuration of Figure 4. The boundary⳵D consists of an acquisition surface⳵D1and a hemisphere⳵D0in the upper half-space with its midpoint at x_⬘. The source domainDsis sited below the

ac-quisition surface⳵D1. When we let the radius of the hemisphere go to infinity and assume that beyond some finite radius the medium is ho-mogeneous, isotropic, and nonporous, the contribution of the bound-ary integral over⳵D0vanishes. Assuming the contrasts⌬Hˆ and ⌬Hˆb defined by equations 55 and 56 are negligible inD, we obtain from equation 54共using equation 57兲

Gˆ 共x

⬘

,x

⬙

,

␻

兲 = −

冕

⳵D1

G¯ˆ 共x

⬘

,x,

␻

兲NxGˆ 共x,x

⬙

,

␻

兲d2x. 共62兲

Multiplying both sides by sˆ共x⬙,␻兲 and integrating over the source domainDs, we obtain共using equation 51兲

uˆ共x

⬘

,

␻

兲 = −

冕

⳵D1

G¯ˆ 共x

⬘

,x,

␻

兲Nxuˆ共x,

␻

兲d2x. 共63兲

This expression formulates forward extrapolation of the wavefield uˆ共x,␻兲 at acquisition surface ⳵D1, because of sources below this sur-face, to any point x_⬘above this surface. By substituting the vectors and matrices for the acoustic situation, we obtain

pˆ共x

⬘

,

␻

兲 = −

冕

⳵D1 兵G¯ˆp,q_共x

_⬘

_,x,

_␻

_兲vˆ j共x,

␻

兲 + G¯ˆ_,jp,f共x

⬘

,x,

␻

兲pˆ共x,

␻

兲其njd2x, 共64兲 vˆi共x

⬘

,

␻

兲 = −

冕

⳵D1 兵G¯ˆ i v,q_共x

_⬘

,x,

␻

兲vˆj共x,

␻

兲 + G¯ˆ_i,jv,f共x

⬘

,x,

␻

兲pˆ共x,

␻

兲其njd2x, 共65兲

which are the well-known Kirchhoff-Helmholtz integrals共Morse and Feshbach, 1953; Berkhout, 1985兲, with many applications in

seismic modeling共Frazer and Sen, 1985; Hill and Wuenschel, 1985; Wenzel et al., 1990; Druzhinin et al., 1998兲. Equation 63 is the gener-alization of the Kirchhoff-Helmholtz integral for any of the wave phenomena considered in this paper.

Inverse wavefield extrapolation

An expression for inverse wavefield extrapolation follows in a similar way from the correlation-type representation in equation 59. Consider the configuration of Figure 5. The boundary⳵D now con-sists of an acquisition surface⳵D0, a horizontal surface⳵D1between x_⬘and the source domainDs, and a cylindrical surface⳵Dcylwith a vertical axis through x_⬘共Figure 5 is a side-view of this configura-tion兲. When we let the radius of this cylindrical surface go to infinity, the contribution of the boundary integral over⳵Dcylvanishes for body waves. The boundary integral over⳵D1contains an evanescent wave contribution and a contribution proportional to the square of the reflection coefficients of the interfaces in domainD共Wapenaar and Berkhout, 1989兲. Ignoring these contributions and assuming that the medium and interfaces inD are lossless and the contrasts⌬Hˆ and⌬Hˆb_{defined by equations 60 and 61 are negligible inD, we} ob-tain from equation 59共using equation 57兲

Gˆ 共x

⬘

,x

⬙

,

␻

兲 = −

冕

⳵D0

KG¯ˆ*共x

⬘

,x,

␻

兲KNxGˆ 共x,x

⬙

,

␻

兲d2x,

共66兲 or, using equation 51,

uˆ共x

⬘

,

␻

兲 = −

冕

⳵D0

KG¯ˆ*共x

⬘

,x,

␻

兲KNxuˆ共x,

␻

兲d2x. 共67兲

This expression formulates inverse extrapolation of the wavefield uˆ共x,␻兲 at acquisition surface ⳵D0, because of sources below this sur-face, to any point x_⬘between this surface and the sources. It is a gen-eralization of the Kirchhoff-Helmholtz integral for inverse wave-field extrapolation 共Schneider, 1978; Berkhout, 1985; Bleistein, 1987; Tygel et al., 2000兲, with applications in seismic migration. Un-like in the derivation of equation 63 for forward extrapolation, we as-sumed that the medium and interfaces inD are lossless. When the medium and/or interfaces inD are not lossless, we can choose the reference parameters as follows: A¯ˆ = Aˆ*_{, B}¯ˆ = − Bˆ†_{, and Y}¯ˆ = Yˆ*_{. As} a consequence, the contrast functions⌬Hˆ and ⌬Hˆb_{defined by} equa-tions 60 and 61 are zero, so equaequa-tions 66 and 67 remain valid. How-ever, this choice of reference parameters implies that the Green’s function G¯ˆ 共x⬘,x,␻兲 propagating through the reference medium is

∂ 1 S r x' x'' x n G(x', x, )−ˆ ∞ ω Ĝ(x', x'', )ω Ĝ(x, x'', )ω ∂ 0

Figure 4. Configuration for forward wavefield extrapolation.

∂ 1 S r x' x'' x n G*(x', x, )−ˆ ∞ ω Ĝ(x', x'', )ω Ĝ(x, x'', )ω ∂ 0 ∂ cyl

growing exponentially共to compensate for the exponential decay in the wavefield uˆ共x,␻兲兲. Hence, the implementation of the inverse ex-trapolation integral for attenuating media should be done with ut-most care to avoid instabilities共Mittet et al., 1995; Zhang and Wap-enaar, 2002兲.

Boundary integral representation (perfect interfaces)

We derive a representation for the scattered wavefield above a per-fectly coupled interface. Consider the configuration of Figure 6, in which⳵D1represents an interface. Assuming the contrasts⌬Hˆ and ⌬Hˆb_{defined by equations 55 and 56 are negligible inD, we obtain} from equation 54共using equation 57兲

Gˆ 共x

⬘

,x

⬙

,

␻

兲 = G¯ˆ 共x

⬘

,x

⬙

,

␻

兲 −

冕

⳵D1

G¯ˆ 共x

⬘

,x,

␻

兲NxGˆ 共x,x

⬙

,

␻

兲d2x. 共68兲

We define the Green’s matrix as a superposition of an incident and a scattered contribution, as stated by

Gˆ 共x,x

⬙

,

␻

兲 = Gˆinc_共x,x

_⬙

_,

_␻

_{兲 + Gˆ}sct_共x,x

_⬙

_,

_␻

_{兲, 共69兲}

where Gˆinc共x,x⬙_,␻兲 = G¯ˆ 共x,x⬙_,␻兲. Equation 68 remains valid when outsideD, i.e., in the lower half-space, the actual and reference me-dium parameters are different. We choose the reference parameters in the lower half-space such that they are continuous across⳵D1, and homogeneous, isotropic, and nonporous beyond some finite domain in the lower half-space. Consequently,

冕

⳵D1 G¯ˆ 共x

⬘

,x,

␻

兲NxGˆinc共x,x

⬙

,

␻

兲d2x = O; 共70兲 hence, Gˆsct共x

⬘

,x

⬙

,

␻

兲 = −

冕

⳵D1 G¯ˆ 共x

⬘

,x,

␻

兲NxGˆsct共x,x

⬙

,

␻

兲d2x, 共71兲

or, using equation 51,

uˆsct共x

⬘

,

␻

兲 = −

冕

⳵D1

G¯ˆ 共x

⬘

,x,

␻

兲Nxuˆsct共x,

␻

兲d2x. 共72兲

Note the analogy with equation 63 for forward wavefield extrapola-tion. The main difference is that uˆ in equation 63 is the upgoing wavefield because of sources below⳵D1, whereas uˆsctin equation 72 is the upgoing scattered wavefield at interface⳵D1because of sourc-es above this interface. This scattered wavefield can be exprsourc-essed in terms of a reflection operator acting on the incident wavefield at⳵D1. For example, for the acoustic situation, it can be written as Nxuˆsct共x,␻兲⬇− R共x,␣兲KNxuˆinc共x,␻兲, where R共x,␣兲 is the local an-gle-dependent reflection coefficient. This is a generalization of what is commonly known as the Kirchhoff approximation共Bleistein, 1984兲.

Boundary integral representation (imperfect interfaces)

The derivation for the scattered wavefield above an imperfect in-terface is somewhat different. The inin-terface is now represented by ⳵Dint, whereas⳵D is a sphere with infinite radius. From equation 54, we thus obtain共using equation 57兲

Gˆ 共x

⬘

,x

⬙

,

␻

兲 = G¯ˆ 共x

⬘

,x

⬙

,

␻

兲 −

冕

⳵Dint

G¯ˆ 共x

⬘

,x,

␻

兲⌬Hˆb共x,

␻

兲Gˆ共x,x

⬙

,

␻

兲d2x,

共73兲 with⌬Hˆb_{defined by equation 56. Various choices are possible for} the reference medium. Let us choose a reference medium that is identical to the actual medium, except that it has an interface with perfect coupling, i.e., Y¯ˆ = O, and hence Z¯ˆ = I. For ⌬Hˆb_{, we thus} ob-tain共assuming small Yˆ兲

⌬Hˆb_{= KM}T

N兵I − Zˆ其M ⬇ j

␻

KMTNYˆ M. 共74兲

Moreover, for this choice, the reference Green’s function G¯ˆ 共x⬘,x_⬙,␻兲 in equation 73 is equal to the actual Green’s function Gˆ 共x⬘,x_⬙,␻兲 in equation 68.

Equation 73 is an integral equation of the second kind for Gˆ 共x_⬘,x_⬙,␻兲. It can be solved iteratively, according to

兵Gˆ共x

⬘

,x

⬙

,

␻

兲其共k兲= G¯ˆ 共x

⬘

,x

⬙

,

␻

兲 −

冕

⳵Dint G¯ˆ 共x

⬘

,x,

␻

兲⌬Hˆb_共x,

_␻

_兲 ⫻兵Gˆ共x,x

⬙

,

␻

兲其共k−1兲d2x, 共75兲 for kⱖ1, with 兵Gˆ共x

⬘

,x

⬙

,

␻

兲其共0兲= G¯ˆ 共x

⬘

,x

⬙

,

␻

兲. 共76兲

Volume integral representation

The derivation for a volume integral representation is similar to that for the imperfect interfaces, but instead of the contrast function ⌬Hˆb_{at the internal interfaces, we consider the contrast function}⌬Hˆ

inD, see Figure 7. For⳵D we take again a sphere with infinite radius. Hence, Gˆ 共x

⬘

,x

⬙

,

␻

兲 = G¯ˆ 共x

⬘

,x

⬙

,

␻

兲 −

冕

D G¯ˆ 共x

⬘

,x,

␻

兲⌬Hˆ共x,

␻

兲Gˆ共x,x

⬙

,

␻

兲d3x, 共77兲 with⌬Hˆ defined by equation 55. The iterative solution of this inte-gral equation is given by

兵Gˆ共x

⬘

,x

⬙

,

␻

兲其共k兲= G¯ˆ 共x

⬘

,x

⬙

,

␻

兲 −

冕

D G¯ˆ 共x

⬘

,x,

␻

兲⌬Hˆ共x,

␻

兲 ⫻兵Gˆ共x,x

⬙

,

␻

兲其共k−1兲d3x, 共78兲 for kⱖ1, with 兵Gˆ共x

⬘

,x

⬙

,

␻

兲其共0兲= G¯ˆ 共x

⬘

,x

⬙

,

␻

兲. 共79兲 For k = 1, equation 78 is the Born approximation, which is frequent-ly used as a representation of primary data in modeling and inversion 共Cohen and Bleistein, 1979; Raz, 1981; Bleistein and Cohen, 1982; Tarantola, 1984; Miller et al., 1987; Wu and Toksöz, 1987; Orista-glio, 1989兲. For k⬎1, equation 78 represents a Neumann series ex-pansion, which can be used for modeling primaries as well as

inter-nal multiples. For a discussion on the convergence aspects, see Fokkema and van den Berg共1993兲. Applications for the prediction of internal multiples in nonlinear inversion are discussed by, e.g., Snieder共1990兲, Ten Kroode 共2002兲 and Weglein et al. 共2003兲.

Surface-related multiple prediction and elimination (convolution approach)

Surface-related multiple prediction and elimination was intro-duced by Berkhout共1985兲 and Verschuur et al. 共1992兲 and it was based on reciprocity theory by Fokkema and van den Berg共1993兲 and van Borselen et al.共1996兲. The latter approach is generalized for the wave phenomena discussed in this paper as follows. Let⳵D con-sist of the acquisition surface⳵D0and a hemisphere⳵D1with infinite radius in the lower half-space共we assume that, beyond some finite radius, the medium in the lower half-space is homogeneous, isotro-pic, and nonporous兲. The Green’s matrix Gˆ共x⬘,x_⬙,␻兲 in this configu-ration共with x⬘and x_⬙at⳵D0兲 represents the actual data, including the multiples related to⳵D0, see Figure 8a. In the half-space below⳵D0, the reference medium is specified as identical to the actual medium. In the upper half-space, the reference parameters are homogeneous, isotropic, and nonporous, and specified as continuous across⳵D0. Hence, the Green’s matrix G¯ˆ 共x_⬘,x_⬙,␻兲 in the reference medium rep-resents the data without surface-related multiples, see Figure 8b. The relation between the two Green’s matrices follows from equation 54 and is given by

G¯ˆ 共x

⬘

,x

⬙

,

␻

兲 − Gˆ共x

⬘

,x

⬙

,

␻

兲 =

冕

⳵D0

G¯ˆ 共x

⬘

,x,

␻

兲NxGˆ 共x,x

⬙

,

␻

兲d2x. 共80兲

This expression can be used as the basis for modeling as well as elim-ination of surface-related multiples. For modeling applications, we assume that G¯ˆ 共x_⬘,x_⬙,␻兲, the response without surface-related multi-ples, is known. Then Gˆ 共x⬘,x_⬙,␻兲, the response with surface-related multiples, can be found by solving equation 80 iteratively, according to 兵Gˆ共x

⬘

,x

⬙

,

␻

兲其共k兲= G¯ˆ 共x

⬘

,x

⬙

,

␻

兲 −

冕

⳵D0 G¯ˆ 共x

⬘

,x,

␻

兲Nx兵Gˆ共x,x

⬙

,

␻

兲其共k−1兲d2x, 共81兲 for kⱖ1. The initial estimate is given by the response without multi-ples; hence,

兵Gˆ共x

⬘

,x

⬙

,

␻

兲其共0兲= G¯ˆ 共x

⬘

,x

⬙

,

␻

兲. 共82兲 When equation 80 is used for multiple elimination, then Gˆ 共x⬘,x_⬙,␻兲 is the known response, and equation 80 is solved iteratively, accord-ing to

Figure 7. Configuration for volume integral representation.

a)

b)

兵G¯ˆ 共x

⬘

,x

⬙

,

␻

兲其共k兲= Gˆ 共x

⬘

,x

⬙

,

␻

兲 +

冕

⳵D0

兵G¯ˆ 共x

⬘

,x,

␻

兲其共k−1兲NxGˆ 共x,x

⬙

,

␻

兲d2x,

共83兲 for kⱖ1. This time, the initial estimate is given by the response with multiples; hence,

兵G¯ˆ 共x

⬘

,x

⬙

,

␻

兲其共0兲= Gˆ 共x

⬘

,x

⬙

,

␻

兲. 共84兲 The product under the integral in equation 83 represents a convolu-tion process, producing high-order multiples from primaries and lower-order multiples, which, after addition to the first term, com-pensate the multiples in Gˆ 共x⬘,x_⬙,␻兲.

Surface-related multiple prediction and elimination (correlation approach)

Schuster共2001兲 and Berkhout and Verschuur 共2003兲 suggest an alternative to the convolution-based multiple prediction and elimi-nation approach, based on correlations. For the configuration dis-cussed above, assuming in addition that the medium is lossless, we obtain from equation 59

KG¯ˆ*共x

⬘

,x

⬙

,

␻

兲K + Gˆ共x

⬘

,x

⬙

,

␻

兲 =

冕

⳵D0 KG¯ˆ*共x

⬘

,x,

␻

兲KNxGˆ 共x,x

⬙

,

␻

兲d2x +

冕

⳵D1 KG¯ˆ*共x

⬘

,x,

␻

兲KNxGˆ 共x,x

⬙

,

␻

兲d2x. 共85兲

The products under the integrals in equation 85 represent a correla-tion process, producing primaries and low-order multiples from higher-order multiples. Unlike in the convolution representation, the integral along⳵D1in equation 85 does not vanish when the medium in the lower half-space is homogeneous, isotropic, and nonporous beyond some finite radius. On the other hand, it vanishes due to scat-tering loss when the medium in the lower half-space is sufficiently inhomogeneous共Wapenaar, 2006兲.

Interferometry (correlation approach)

Seismic interferometry deals with the generation of new seismic responses by cross-correlating wavefield measurements at different receiver positions 共Claerbout, 1968; Weaver and Lobkis, 2001; Schuster, 2001; Wapenaar et al., 2002; Campillo and Paul, 2003; Derode et al., 2003; Schuster et al., 2004; Sabra et al., 2005; Draga-nov et al., 2007兲. The measurements take place in the actual medium, so, the basic expression for interferometry is obtained by taking the reference state equal to the actual state in the representation of the correlation-type, equation 59. Using the symmetry properties of Gˆ , Aˆ , Bˆ, Zˆ_⬘, and Nx, this yields

␹

D共x

⬙

兲Gˆ共x

⬘

,x

⬙

,

␻

兲 +

␹

D共x

⬘

兲Gˆ†共x

⬙

,x

⬘

,

␻

兲 = −

冖

⳵D Gˆ 共x

⬘

,x,

␻

兲NxGˆ†共x

⬙

,x,

␻

兲d2x +

冕

D Gˆ 共x

⬘

,x,

␻

兲⌬Hˆ共x,

␻

兲Gˆ†_共x

_⬙

_,x,

_␻

_兲d3_x +

冕

⳵Dint Gˆ 共x

⬘

,x,

␻

兲K⌬Hˆb*共x,

␻

兲KGˆ†共x

⬙

,x,

␻

兲d2x, 共86兲 with ⌬Hˆ = − 2

␻

I共Aˆ兲 + Bˆ + Bˆ†, 共87兲 ⌬Hˆb_{= M}T_{兵J − Zˆ}†_{JZˆ 其M. 共88兲}

Equation 86 is a general representation of the Green’s matrix be-tween x_⬙and x_⬘in terms of cross-correlations of observed fields at x_⬙ and x_⬘because of sources at x on the boundary⳵D, on the internal imperfect interfaces⳵Dint, as well as in the domainD. The inverse Fourier transform of the left-hand side is ␹D共x⬙兲G共x⬘,x⬙,t兲 +␹D共x⬘兲GT共x⬙,x⬘,− t兲, from which G共x⬘,x⬙,t兲 is obtained by taking the causal part共assuming x⬙is located inD兲. When the medium and interfaces are lossless, it suffices to have sources on⳵D only, see Fig-ure 9. Note that⳵D is not necessarily a closed surface: When the me-dium is sufficiently inhomogeneous ⳵D can be an open surface 共Wapenaar, 2006兲. On the other hand, when the medium is dissipa-tive throughoutD and the radius of⳵D is sufficiently large, the boundary integral vanishes and sources are required throughoutD 共Snieder, 2006; Snieder et al., 2007兲.

The application of equation 86 in its current form requires inde-pendent measurements of the impulse responses of different types of sources at all x involved in the integrals. The right-hand side can be modified into a direct cross-correlation共i.e., without the integrals兲 of diffuse field observations at x_⬙and x_⬘, the diffusivity being caused by a distribution of uncorrelated noise sources, either on⳵D 共for loss-less media兲 or in D 共for dissipative media兲 共Wapenaar et al., 2006兲.

Equation 86 has also important applications in efficient modeling and inversion共van Manen et al., 2005, 2006兲. As mentioned above, for the lossless situation, only the boundary integral over⳵D needs to be evaluated. Hence, by modeling the responses of a distribution of

x' x x'' n Ĝ(x', x'', ) ω Ĝ(x', x, ) ω Ĝ(x'', x, ) ω ∂

sources on the 2D boundary, equation 86 allows us to determine the responses of all possible sources in the 3D volume enclosed by the boundary. This is very useful, for example, in nonlinear inversion, where the Green’s functions between all possible pairs of points in a volume are needed共see e.g., Weglein et al., 2003兲.

Interferometry (convolution approach)

When the dissipation of the medium is significant, interferometry according to the correlation approach requires a distribution of sources throughout the medium. As an alternative, Slob and Wap-enaar共2007兲 and Slob et al. 共2007兲 propose a convolution approach to interferometry. Taking the reference state equal to the real state in the convolution-type representation of equation 54 and using the symmetry property of Gˆ gives

兵

␹

D共x

⬙

兲 −

␹

D共x

⬘

兲其Gˆ共x

⬘

,x

⬙

,

␻

兲

=

冖

⳵D

Gˆ 共x

⬘

,x,

␻

兲NxKGˆT共x

⬙

,x,

␻

兲Kd2x. 共89兲

This is a representation of the Green’s matrix between x_⬙and x_⬘in terms of cross-convolutions of observed fields at x_⬙and x_⬘due to sources at x on the boundary⳵D only. Note that one of the observa-tion points should be inside this boundary and the other outside, see Figure 10共otherwise, the left-hand side of equation 89 vanishes兲. There are no restrictions with respect to the losses in the medium. The application of equation 89 requires independent measurements of the impulse responses of different types of sources at all x苸⳵D; a modification for uncorrelated noise sources is not possible for the convolution approach.

CONCLUSIONS

Starting with a unified matrix-vector-form wave equation and boundary conditions for acoustic, electromagnetic, elastodynamic, poroelastic, and electroseismic waves, we derived general convolu-tion- and correlaconvolu-tion-type wavefield representations. We discussed applications including forward and inverse wavefield extrapolation, boundary integral representations for perfect and imperfect interfac-es, volume integral representations共the Born approximation and the Neumann series expansion兲, multiple elimination, and seismic inter-ferometry, the latter two both in terms of convolutions and correla-tions. Each of these applications is a generalization of the well-es-tablished acoustic representations for any of the wave phenomena governed by the unified wave equation.

ACKNOWLEDGMENTS

This work is supported by the Netherlands Research Centre for In-tegrated Solid Earth Science. I thank the editors and reviewers 共Roald van Borselen and Enru Liu兲 for their constructive comments.

APPENDIX A

THE DIVERGENCE THEOREM OF GAUSS IN MATRIX-VECTOR FORM

For a scalar field a共x兲, the divergence theorem of Gauss reads

冕

D

⳵

ia共x兲d3x =

冖

⳵D

a共x兲nid2x. 共A-1兲

Here, we modify this theorem for the differential operator matrix Dx appearing in equations 1 and 9. Let DIJdenote the operator in row I

and column J of matrix Dx. The symmetry of Dx共equation 12兲 im-plies DIJ= DJI. We define a matrix Nx, which contains the compo-nents of the normal vector n, organized in the same way as matrix Dx. Hence, NIJ= NJI, where NIJdenotes the element in row I and column

J of matrix Nx. If we replace the scalar field a共x兲 by aI共x兲bJ共x兲, we

may generalize equation A-1 to

冕

D

DIJ兵aI共x兲bJ共x兲其d3x =

冖

⳵D

aI共x兲bJ共x兲NIJd2x, 共A-2兲

where the summation convention applies to repeated capital Latin subscripts, which may run from 1 to 4, 6, 12, 16, 18, or 22, depending on the choice of operator Dx. Applying the product rule for differen-tiation and using the symmetry property DIJ= DJI, we obtain for the

integrand in the left-hand side of equation A-2,

DIJ共aIbJ兲 = aIDIJbJ+ 共DJIaI兲bJ,

= aTDxb +共Dxa兲Tb, 共A-3兲

where a and b are vector functions, containing the scalar functions aI共x兲 and bJ共x兲, respectively. Rewriting the integrand in the

right-hand side of equation A-2 in a similar way, we obtain the divergence theorem of Gauss in matrix-vector form

冕

D

兵aT

Dxb +共Dxa兲Tb其d3x =

冖

⳵D

aTNxbd2x. 共A-4兲

Finally, we consider a variant of this equation. We replace a by Ka, where K is the real-valued diagonal matrix introduced in equations 10–12, obeying the property K = K−1_{. Using equation 12, we thus} obtain

冕

D 兵aT KDxb −共Dxa兲TKb其d3x =

冖

⳵Da T KNxbd2x. 共A-5兲 REFERENCES

Aki, K., and P. G. Richards, 1980, Quantitative seismology, vol. I: W. H. Freeman & Company.

Auld, B. A., 1973, Acoustic fields and waves in solids: Wiley-Interscience. Berkhout, A. J., 1985, Seismic migration: Imaging of acoustic energy by

wave field extrapolation: Elsevier Science Publ. Co., Inc.

Berkhout, A. J., and D. J. Verschuur, 2003, Transformation of multiples into primary reflections: 73rd Annual International Meeting, SEG, Expanded Abstracts, 1925–1928. x' x x'' n Ĝ(x', x'', )ω Ĝ(x', x, )ω Ĝ(x'', x, )ω ∂