Seismic and electromagnetic interferometry: Retrieval of the earth's reflection response using crosscorrelation

Interferometry –

Retrieval of the Earth’s Reflection

Response Using Crosscorrelation

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus Prof. dr. ir. J.T. Fokkema, voorzitter van het College voor Promoties,

in het openbaar te verdedigen,

op maandag 10 december 2007 om 15.00 uur door

Deyan DRAGANOV

Samenstelling promotiecommissie:

Rector Magnificus voorzitter

Prof. dr. ir. C.P.A. Wapenaar, Technische Universiteit Delft, promotor Prof. dr. ir. G.E.W. Bauer, Technische Universiteit Delft

Prof. dr. W.A. Mulder, Technische Universiteit Delft Prof. dr. ir. A. Gisolf, Technische Universiteit Delft Prof. dr. ir. A.J. Berkhout, Technische Universiteit Delft Drs. A. Verdel, Shell International Exploration and

Production

Dr. A. Curtis, University of Edinburgh

Jan Thorbecke heeft als begeleider in belangrijke mate aan de totstandko-ming van het proefschrift bijgedragen.

ISBN-978-90-9022528-9

Copyright c _{2007, by D. Draganov, Delft University of Technology, Delft, The} Nether-lands.

All rights reserved. No part of this publication may be reproduced, stored in a retrieval system or transmitted in any form or by any means, electronic, mechanical, photocopying, recording or otherwise, without the prior written permission of the author.

SUPPORT

The research reported in this thesis was financially supported by the Netherlands Research Centre for Integrated Solid Earth Science ISES and by the Dutch Technology Foundation STW, applied science division of NWO, and the Technology Program of the Ministry of Economic Affairs (grant DTN.4915).

Cover: The amplitude spectrum of the transmission response from Figure 3.2(a) as it is seen in the temporal frequency domain.

Table of Contents iii

1 Introduction 1

1.1 Why do we want to retrieve the reflection response? . . . 1 1.2 What are Seismic and Electromagnetic Interferometry? . . . 2 1.3 Outline of this thesis . . . 4 2 Theory of Seismic and Electromagnetic Interferometry 7 2.1 SI and EMI relations based on a two-way wavefield reciprocity 11 2.1.1 Reciprocity theorems . . . 11 2.1.2 SI relations for impulsive sources in an acoustic medium 16 2.1.3 SI relation for impulsive sources in an acoustic medium

with a free surface . . . 22 2.1.4 SI relation for transient sources in acoustic medium . 22 2.1.5 SI relation for noise sources in acoustic medium . . 23 2.1.6 SI relations for impulsive sources in an elastic medium 24 2.1.7 SI relation for impulsive sources in an elastic medium

2.1.11 EMI relations for transient sources . . . 41 2.1.12 EMI relations for noise sources . . . 42 2.2 SI relations based on a one-way wavefield reciprocity . . . . 43 2.2.1 Reciprocity theorems . . . 43 2.2.2 SI relations for impulsive sources in an acoustic

me-dium . . . 48 2.2.3 SI relation for transient sources in an acoustic medium 51 2.2.4 SI relation for noise sources in an acoustic medium . 51 2.2.5 SI relations for a direct migration of noise recordings 52 2.2.6 SI relations for impulsive sources in an elastic medium 54 2.2.7 SI relation for transient sources in an elastic medium 55 2.2.8 SI relation for noise sources in an elastic medium . . 56 3 Retrieval of the reflection response: numerical modelling results 59 3.1 Retrieval of the acoustic reflection response . . . 61

3.1.1 Retrieval of the acoustic reflection response from re-cordings from transient subsurface sources . . . 61 3.1.2 Retrieval of the acoustic reflection response from

re-cordings from white-noise subsurface sources . . . . 77 3.2 Retrieval of the elastodynamic reflection response . . . 85

3.2.1 Retrieval of the elastodynamic reflection response

from recordings from transient subsurface sources . 85 3.2.2 Retrieval of the elastodynamic reflection response

from recordings from white-noise subsurface sources 92 3.3 Retrieval of the electromagnetic reflection response . . . 96 4 Retrieval of the seismic reflection response: application to

laboratory and field data 101

4.1 Application to laboratory data from transient sources . . . . 102 4.2 Application to field data from background-noise sources . . 111

5 Conclusions 125

Bibliography 129

Samenvatting 141

Acknowledgements 145

Introduction

1.1

Why do we want to retrieve the reflection

res-ponse?

a subsurface structure, is even larger as the imaging of the deep subsurface structures requires very powerful sources.

The above-mentioned problems imply that certain areas of the Earth are not accessible for imaging or that we should look for alternative methods of ob-taining seismic reflection recordings. Such an alternative is (passive) Seis-mic Interferometry, the subject of this thesis. It can be applied at all scales - laboratory, exploration, regional, and global (for the later see Ruigrok [2006]).

1.2

What are Seismic and Electromagnetic

Interfe-rometry?

Schuster[2001] and Schuster et al.[2004] introduced the concept of Inter-ferometric Imaging for the process of constructing a migrated seismic image from crosscorrelated results. Based on stationary-phase arguments, the au-thors showed that the method is applicable not only to passively acquired seismic noise data from subsurface noise sources, but also to conventional reflection data shot and recorded at the surface. Using a one-way reciprocity theorem of the correlation type,Wapenaar et al.[2002,2004b] proved Claer-bout’s conjecture for any 3D inhomogeneous acoustic or elastic medium for both impulsive and white-noise sources, where the sources could be at depth or at the surface. Berkhout and Verschuur [2003,2006] proposed a scheme for turning multiple reflections into primaries based on a weighted corre-lation process. Slob et al. [2006a,b] showed how one can retrieve magnetic recordings from crosscorrelation or crossconvolution of electro-magnetic fields. They named this process Electroelectro-magnetic Interferometry (EMI). Later, the concept of retrieval of new responses through crosscorre-lation was extended to any type of field [Wapenaar et al.,2006].

arrivals between different seismological stations in Mexico. Derode et al. [2003] showed by physical reasoning how the time-reversal principle can be applied to retrieve the Green’s function between two points. This result is very closely connected to Green’s function retrieval through crosscorrelation since a crosscorrelation is nothing but a convolution of a signal with a second signal which is reversed in time.Snieder[2004,2006] showed that the main contribution to the retrieved Green’s function comes from sources around the stationary point. Curtis et al.[2006] showed the connection between the different methods in the different fields in a comprehensive way.

Because of the independent developments of the methods in the different fields, different names appeared for the same process of retrieval of a sig-nal between two points from crosscorrelation or time reversal. Examples are: acoustic daylight imaging, Green’s function retrieval, interferometric imaging, virtual source method, and passive seismics. To avoid confusion, the term ”Seismic Interferometry” was adopted for a special issue of Geo-physics (July-August, 2006) to describe the application of these methods to seismic exploration and seismology. Later, the term ”Electromagnetic In-terferometry” was adopted to describe the retrieval of electromagnetic fields for electromagnetic exploration.

1.3

Outline of this thesis

Chapter 2 gives the derivations of relations for SI and EMI. These derivations are based on two-way and one-way wavefield reciprocity relations of the correlation type.

Chapter 3 shows results from numerical modelling examples that serve to il-lustrate the SI and EMI relations and to investigate the applicability of SI and EMI relations to practical situations. Using finite-difference acoustic and finite-element elastic modelling codes, recordings at the surface were mod-elled from deterministic as well as noise sources in the subsurface. These recordings were then crosscorrelated in accordance with the derived SI and EMI relations to retrieve recordings at the surface as if from sources at the surface.

Theory of Seismic and

Electromagnetic

Interferometry

The reciprocity theorems can describe relations between two-way wavefields or between one-way wavefields. Such relations are called two-way wavefield reciprocity theorems and one-way wavefield reciprocity theorems, respec-tively.

In section2.1, followingWapenaar and Fokkema[2006],Slob et al.[2007], andSlob and Wapenaar[2007], derivations are presented of the SI and EMI relations using two-way wavefield reciprocity theorems. Section2.2, follow-ing Wapenaar et al. [2004b], is dedicated to derivations of the SI relations using one-way wavefield reciprocity theorems. Even though the end results look very similar, there are some basic differences in the relations in the two sections, which make them useful in different practical applications.

The SI and EMI relations are derived under the assumption of a lossless medium. Of course in practice there are always losses in the medium. As shown bySlob et al.[2007];Slob and Wapenaar[2007], in case of losses in the medium, the derived SI and EMI relations can still be used, but later ar-rivals, i.e., primaries and multiples with long travel paths, may not be recon-structed. Furthermore, as shown in Section3.1.1, big losses in the medium may give rise to ghost events.

In both sections the derivations of the SI and EMI relations follow the same basic workflow:

• defining the system of two equations describing the wave propagation

in a matrix-vector from;

• choosing an interaction quantity and deriving from it a reciprocity

the-orem;

• for each of the cases (acoustic, elastic, or electromagnetic) deriving the

relationship for seismic (electromagnetic) interferometry in the case of – impulsive sources,

– transient sources, – and noise sources.

The relations presented in this chapter are given in the frequency domain. A time-Fourier transform of a time- and space-dependent quantity is defined asuˆ (x, ω) =R∞

t=0u(x, t) e

−jωt_{dt, where x = (x}

thesis, Einstein’s summation convention applies to repeated lowercase Latin subscripts from 1 to 3, unless stated otherwise.

The list below describes in short all the symbols used in the this chapter in the order they appear in the text.

Symbol Meaning

D Domain

∂D Surface bounding the domain_D

n_{= nm}i_{m Outward-pointing vector normal to ∂D; and the vectors} i_{m denote the base vectors of a Cartesian reference frame in a}

three-dimensional Euclidian space in the mth_-direction.

x_{= xm}i_{m Cartesian coordinates vector (m)} ρ Volume density of mass (kg/m3)

κ Compressibility (1/Pa=m2/N)

spqkl Compliance tensor (1/Pa)

cijkl Stiffness tensor (P a)

εkr Electric permittivity (s/(Ω m))

µji Magnetic permeability (Ω s/m)

j Imaginary unit

ω Angular frequency (rad)

u Wavefield vector

s Source vector

A Material properties matrix D Differentiation operator matrix

p Acoustic pressure (Pa)

vi Particle velocity vector (m/s)

τij stress tensor (Pa)

Er Electric field vector (V/m)

Hi Magnetic field vector (A/m)

q Source distribution in terms of density of volume injection rate (1/s)

hpq Source distribution in terms of external deformation-rate density (1/s)

J_ke External source volume densities of electric currents (A/m2)

Jm

j External source volume densities of magnetic currents (V/m2)

A, B Two states of the domain D

K Operator matrix_{= diag(−1, 1)} T _{Vector or matrix transposition}

δik Kronecker delta function

∆ijkl = 1₂(δikδjl+ δjkδil)

ǫkmi The antisymmetric (Levi-Civita) tensor of rank three

Nx Operator matrix= antidiag(ni, ni)

† _{Vector or matrix complex conjugation and transposition} ∗ _{Complex conjugation}

G Green’s function (1/m)

δ Dirac Delta function

ℜ Real part of a parameter

c Propagation velocity (m/s)

α ray angle (degr)

s Wavefield source

S Source power spectrum

F Shaping filter

N Source noise spectrum

cP Propagation velocity for P-waves (m/s)

cS Propagation velocity for S-waves (m/s)

ζij Inverse of the magnetic permeability (m/(Ω s))

ξij Inverse of the electric permittivity (Ω m/s)

P Flux-normalized one-way wavefield quantity vector B One-way pseudo-differential operator matrix S Flux-normalized one-way source quantity vector

L Composition operator L−1 Decomposition operator

P+(−) Downgoing(upgoing) one-way wavefield

S+(−) _{Downgoing(upgoing) one-way source field}

Φ+(−) _{Downgoing(upgoing) P-wave potential}

Ψ+(−) _{Downgoing(upgoing) S-wave potential}

N Operator matrix_{= antidiag(1, −1)} J Operator matrix_{= diag(1, −1)}

R Flux-normalized reflection response of a medium

T Flux-normalized transmission response of a medium

r− _{Reflection coefficient of the free surface}

W+(−) Forward extrapolation operator

R Flux-normalized reflection response matrix T Flux-normalized transmission response matrix r− Matrix with reflection coefficient of the free surface I Identity matrix

F Shaping filter matrix

NNN Noise spectra matrix

2.1

SI and EMI relations based on a two-way

wave-field reciprocity

2.1.1 Reciprocity theorems

Consider a 3D inhomogeneous, lossless medium with a bounded domain_D in it. Three different cases are looked at for the medium - acoustic, elas-tic, or electromagnetic. In the acoustic case, the medium is described with the material properties mass density ρ(x) and compressibility κ (x); in the

However, because the medium is taken lossless, the electric and magnetic conductivities are set to zero. Domain_{D is enclosed by a boundary ∂D with} an outward-pointing normal vector n. For the medium, inside as well as outside_{D, in the frequency domain the system of two equations defining the} wave propagation can be written in matrix-vector form as

jωA(x)uˆ(x, ω) +D ˆu(x, ω) = ˆs(x, ω) . (2.1)

In the above equation, the variableuˆ(x, ω) represents the wave fields, ˆs(x, ω)

- the source vector, A(x) - the frequency-independent material properties,

andD is an operator matrix representing differentiation with respect to the three spatial dimensions. Table 2.1 shows the expansion of these variables in the acoustic, elastic, and electromagnetic case.

To obtain reciprocity relations, two independent states_{A and B are} consid-ered in domain _{D . State A is described by the wave quantity} uˆA(x, ω), by the source quantity ˆsA(x, ω) and by the material matrixA(x). State B is characterized by the same quantities, but the corresponding subscripts_A are replaced by subscripts_{B. Further, to arrive at a reciprocity relation} be-tween the wavefields in the two states, wavefield interactions are considered in domainD

One wavefield interaction can be described by a local convolution between the wavefields in the two states [de Hoop,1995] and can be written as

ˆ uT_A(x, ω)KD ←−+D  ˆ uB(x, ω) = −  D ˆuA(x, ω) T K ˆuB(x, ω) +uˆTA(x, ω)K  D ˆuB(x, ω)  . (2.2)

In the above equation, the superscript T denotes matrix transposition, D

←−

means that the differentiation operatorDacts on the wave quantity left of it, andK_{= diag (−1, 1). Further, to obtain the right-hand side of the equation,} the properties KT _{= K, KD}

= −DK, andDT _{= Dwere used. The latter} property, in the acoustic case, is easily verified by inspecting Table2.1. To verify this property in the elastic case, the relation∆ijkl = ∆klji should be used. In the electromagnetic case, the propertyDT _{=Dis verified by using} the fact that_−ǫkmi= ǫimk. A similar remark holds for ǫjmr.

Parameter Acoustic case Elastic case Electromagnetic case ˆ u(x, ω) pˆ(x, ω) ˆ vi(x, ω) ! ˆ vi(x, ω) −ˆτkl(x, ω) ! ˆ_{Er(x, ω)} ˆ Hi(x, ω) ! ˆs_{(x, ω)} qˆ(x, ω) ˆ fi(x, ω) ! ˆ fi(x, ω) ˆ hpq(x, ω) ! − ˆJe k(x, ω) − ˆJ_jm(x, ω) ! A_(x) κ(x) 0 0 ρ(x) ! ρ_(x) 0 0 spqkl(x) ! εkr(x) 0 0 µji(x) ! D 0 ∂i ∂i 0 ! 0 ∆ijkl∂j ∆pqji∂j 0 ! 0 _−ǫkmi∂m ǫjmr∂m 0 !

Table 2.1: Expansion of the variables in equation 2.1 in the acoustic, elastic, or electromagnetic case. In the first row, p stands for acoustic pressure,ˆ ˆvi denotes

the particle velocity vector, τˆkl stands for the stress tensor, ˆEr describes the

elec-tric field vector, and ˆHi describes the magnetic field vector. In the second row,qˆ

represents the source distribution in terms of volume injection-rate density, ˆfi

rep-resents the source distribution in terms of external volume force density, ˆhpqstands

for the external deformation-rate density, and ˆJe

k and ˆJjm stand for the external

source volume densities of electric and magnetic currents, respectively. In the third row, κ denotes the compressibility, ρ denotes the volume density of mass, spqkl

de-scribes the compliance, εkrdescribes the electric permittivity, and µjidenotes the

magnetic permeability. In the fourth row,∆ijkl = 1₂(δikδjl+ δjkδil), where δabis

the Kronecker delta function, while ǫkmi is the antisymmetric (Levi-Civita) tensor

of rank three, where ǫkmi = 1, for kmi = {123, 231, 312}, and ǫkmi = −1, for

equation2.2and using the propertyAT A(x)K=KAA(x), results in ˆ uTA(x, ω)K D_←−uˆB(x, ω) +uˆTA(x, ω)KD ˆuB(x, ω) = ˆ uT_A(x, ω)K ˆsB(x, ω) − ˆsTA(x, ω)K ˆuB(x, ω) + jωuˆT_A(x, ω)K  AA(x) −AB(x)  ˆ uB(x, ω) . (2.3)

Equation2.3represents the local form of the two-way wavefield convolution-type reciprocity theorem [Rayleigh, 1878; de Hoop, 1988; Fokkema and van den Berg, 1993] because multiplication of two quantities in the fre-quency domain corresponds to a convolution in the time domain.

Integration of equation2.3over the domain_{D and application of Gauss’} di-vergence theorem to the left-hand side of equation2.3yields (after swapping the places of the two sides)

Z D  ˆ uT_A(x, ω)K ˆsB(x, ω) −ˆsTA(x, ω)K ˆuB(x, ω) +jωuˆTA(x, ω)K  AA(x) −AB(x)  ˆ uB(x, ω)  d3x₌ I ∂D ˆ uT_A(x, ω)KNxuˆB(x, ω) d2x, (2.4)

where in the acoustic case Nx = antidiag (ni, ni), in the elastic case it is Nx = antidiag (∆ijklnj,∆pqjinj), and in the electromagnetic case it is

Nx = antidiag (−ǫkminm, ǫjmrnm). Equation 2.4 is known as the global form of the two-way wavefield convolution-type reciprocity theorem. When the material quantities are assumed to be the same for both states, i.e., AA(x) =AB(x), equation2.4can be rewritten as

Z D  ˆ uT_A(x, ω)K ˆsB(x, ω) −ˆsAT (x, ω)K ˆuB(x, ω) d3x= I ∂D ˆ uT_A(x, ω)KNxuˆB(x, ω) d2x. (2.5)

as ˆ u†_A(x, ω)D ←−+D  ˆ uB(x, ω) =  D ˆuA(x, ω) † ˆ uB(x, ω) +uˆ†A(x, ω)  D ˆuB(x, ω)  , (2.6)

where the superscript†means complex conjugation and transposition, hence ˆ

u†_A(x, ω) is the time-reversed variant ofuˆT

A(x, ω). Then, after substitution of equation2.1 into the right-hand side of equation 2.6 and making use of the propertyA†_A(x) = AA(x), which follows from the invariance to time-reversal of the wave equation, one arrives at

ˆ u†_A(x, ω)D ←−uˆB(x, ω) +uˆ†A(x, ω)D ˆuB(x, ω) = ˆ u†_A(x, ω)ˆsB(x, ω) + ˆs † A(x, ω)uˆB(x, ω) + jωuˆ†_A(x, ω)  AA(x) −AB(x)  ˆ uB(x, ω) . (2.7)

Equation2.7represents the local form of the two-way wavefield correlation-type reciprocity theorem [Bojarski, 1983; de Hoop, 1988; Wapenaar and Berkhout, 1989; Fokkema and van den Berg, 1993] since a multiplication of a quantity with a complex-conjugate of another quantity in the frequency domain corresponds to a correlation in the time domain.

Integration of equation 2.7 over the domain _{D and application of Gauss’} divergence theorem to the left-hand side of equation2.7gives

Z D n ˆ u†_A(x, ω)ˆsB(x, ω) + ˆs†A(x, ω)uˆB(x, ω) +jωuˆ†_A(x, ω)  AA(x) −AB(x)  ˆ uB(x, ω)  d3x₌ I ∂D ˆ u†_A(x, ω)NxuˆB(x, ω) d2x. (2.8)

For the case when the material properties in both states are the same, equa-tion2.8can be rewritten as

Z D n ˆ u†_A(x, ω)ˆsB(x, ω) + ˆs † A(x, ω)uˆB(x, ω) o d3x₌ I ∂D ˆ u†_A(x, ω)NxuˆB(x, ω) d2x. (2.9)

The the global forms of the reciprocity theorems2.4and2.8(or2.5and2.9, respectively) can be used to obtain the wavefields inside domain_{D only from} measurements of the wavefields along its boundary ∂_{D. In this sense, these} theorems are used, for example, to derive integral equations for forward and inverse two-way wavefield extrapolation [Wapenaar and Berkhout,1989]. In the following subsections, the reciprocity theorems2.5 and 2.9are used to derive source-receiver reciprocity relations followed by relations for seismic and electromagnetic interferometry.

2.1.2 SI relations for impulsive sources in an acoustic medium For an acoustic medium, as specified in the acoustic-case column of Table 2.1, the reciprocity relations2.5and2.9can be rewritten as

To make use of the reciprocity theorems, certain choices should be made for the wavefields and the source fields in the states _{A and B. Taking only} impulsive sources of volume injection-rate densityqˆ(x, ω) at points xAand

x_{B in D for statesA and B, respectively, one can define the recorded}

pres-surepˆ(x, ω) at x in both states as the observed impulse response (Green’s

function ˆG x, xA(B), ω
) (see Table 2.2). The expressions for the particle velocities observed at x for states _{A and B are obtained from the} substi-tution of the respective quantities for qˆ(x, ω), ˆfi(x, ω), and ˆp(x, ω) from Table2.2into the wave equation2.1in case of an acoustic medium.

Parameter State_A State_B

ˆ q(x, ω) δ_{(x − x}A) δ(x − xB) ˆ fi(x, ω) 0 0 ˆ p(x, ω) Gˆ(x, xA, ω) Gˆ(x, xB, ω) ˆ vi(x, ω) _jωρ(−1_x)∂iGˆ(x, xA, ω) _jωρ(−1_x)∂iGˆ(x, xB, ω)

Table 2.2: Choices for wavefield and source parameters for the acoustical states_A andB to be used in the two-way wavefield acoustic reciprocity theorems.

equation2.12 will go to zero as it will be of order O(∆−1_{) [Fokkema and} van den Berg, 1993]. If for this ∂_{D the right-hand side vanishes, then,} be-cause it is independent of the choice of the boundary, it should vanish for any choice of ∂_{D as long as the points x}A and xB are inside it. This results in the acoustic source-receiver reciprocity relation

ˆ

G(xB, xA, ω) = ˆG(xA, xB, ω) . (2.13) This relation shows the well-known property that interchanging the place of source and receiver will not change the measurement result.

Let one now substitute the quantities from the middle and the right columns of Table2.2 into their corresponding places in the correlation-type acoustic reciprocity theorem 2.11. Making again use of the sifting property of the delta function gives

sources on ∂_{D and observed Green’s functions ∂}iG xˆ A(B), x, ω
ni due to dipole sources on ∂_D.

Equation2.15is exact and valid for any lossless inhomogeneous (inside, as well as outside ∂_{D) acoustic medium. It represents the basis for acoustic SI.} Van Manen et al.[2005] developed an effective modelling scheme based on a relation very similar to equation 2.15, where the difference comes from using another source definition for the Green’s function.

Even though equation 2.15 is very useful for modelling or laboratory in-vestigations, there are several difficulties in its application to field data, like for example in petroleum exploration. One difficulty is the need to evaluate two correlation products. To overcome this, let one assume that in a small region (relative to the dominant wavelength) around ∂_{D the medium} param-eters change smoothly. The Green’s functions observed in domain_{D can be} written as a sum of waves that propagate initially inward from the sources on ∂_{D (superscript in) and of waves that propagate initially outward from} the sources (superscript out):

ˆ

G xA(B), x, ω
= ˆGin xA(B), x, ω
+ ˆGout xA(B), x, ω
. (2.16)

or 2ℜn ˆG(xA, xB, ω) o = I ∂D −1 jωρ(x)n ˆG in∗ (xA, x, ω) ∂iGˆin(xB, x, ω) −∂iGˆin(xA, x, ω) ∗ ˆ Gin(xB, x, ω) + ˆGout ∗ (xA, x, ω) ∂iGˆout(xB, x, ω) −∂iGˆout(xA, x, ω) ∗ ˆ Gout(xB, x, ω) o nidx2 + I ∂D −1 jωρ(x)n ˆG in∗ (xA, x, ω) ∂iGˆout(xB, x, ω) −∂iGˆin(xA, x, ω) ∗ ˆ Gout(xB, x, ω) + ˆGout ∗ (xA, x, ω) ∂iGˆin(xB, x, ω) −∂iGˆout(xA, x, ω) ∗ ˆ Gin(xB, x, ω) o nid2x. (2.18) When the dominant wavelengths of the fields are small compared to the sizes of the inhomogeneities (high-frequency regime), then for each constituent

ˆ

Ginw xA(B), x, ω
and ˆGoutw xA(B), x, ω
direct wave, scattered wave, etc. -of the Green’s function, one has

∂iGˆinw xA(B), x, ω
ni = −j ω c(x)|cos αw(x)| ˆG in w xA(B), x, ω
(2.19a) ∂iGˆoutw xA(B), x, ω
ni = +j ω c(x)|cos αw(x)| ˆG out w xA(B), x, ω
, (2.19b) where c(x) is the local propagation velocity and αw(x) is the local angle between the pertinent ray and the normal on ∂_{D. Note that in the above} equation, no summation takes place over w. The main contribution to the integrals in equation2.18 comes from regions on ∂_{D around the stationary} points [Schuster et al., 2004; Snieder, 2004; Wapenaar et al., 2004a]. For such regions,_{|cos α}w(x)| for ˆGw(xA, x, ω) and ˆGw(xB, x, ω) will be iden-tical. This means that in equation 2.18 the second integral will disappear, while the terms under the first integral will add together. Hence,

Because the inward- and outward-propagating waves at the points x on ∂_D cannot be measured separately, equation2.16 for states_{A and B is used to} rewrite equation2.20as 2ℜn ˆG(xA, xB, ω) o + ghosts = I ∂D 2 jωρ(x)  ∂iGˆ(xA, x, ω) ∗ ˆ G(xB, x, ω) nid2x, (2.21) where ghosts= I ∂D 2 jωρ(x) n ∂iGˆin(xA, x, ω) ∗ ˆ Gout(xB, x, ω) +∂iGˆout(xA, x, ω) ∗ ˆ Gin(xB, x, ω) o nid2x. (2.22) The left-hand side of equation2.21 contains not only the retrieved Green’s function, but also spurious events (ghosts) that will appear as a result from the crosscorrelations in the right-hand side (see Figure 3.12). These ghosts result from the integration over correlation products of inward-propagating waves in one of the states with outward-propagating waves in the other state. When the sources are randomly distributed in space, i.e., the surface ∂_{D is} sufficiently irregular, the correlation products will not interfere coherently and the ghost terms will be strongly weakened and can be ignored, while the retrieved Green’s function will correctly contain all scattered fields from inside and outside ∂_{D (see Figure}3.15andDraganov et al.[2003]).

If, on the other hand, the medium at and outside ∂_{D is assumed to be} homo-geneous with propagation velocity c and mass density ρ, then ghosts = 0.

Furthermore, if ∂_{D is taken to be a sphere with a sufficiently large radius,} then all rays are normal to ∂_{D and consequently α = 0. Then, from equation} 2.21one finally obtains

2ℜn ˆG(xA, xB, ω) o ≈ _ρc2 I ∂D ˆ G∗(xA, x, ω) ˆG(xB, x, ω) d2x. (2.23) The approximation in equation 2.23 involves only amplitude errors; when

∂D differs significantly from a sphere, these errors may be significant.

canceled completely when equation2.15is used, may here give rise to arte-facts. Nevertheless, the application of relation 2.23 will correctly retrieve the phases of all arrivals. That is why relation2.23is considered acceptable for SI in practical applications.

2.1.3 SI relation for impulsive sources in an acoustic medium with a free surface

Let the medium _{D represent part of the Earth and be bounded by a free} surface. In that case, the boundary ∂_{D can be written as consisting of two} parts: ∂_D0, coinciding with the free surface, and ∂Dm, representing the part of the boundary in the subsurface. The integral on the right-hand side of equation2.14can be split into an integral over ∂_Dmand an integral over ∂D0. At the free surface, the Green’s function ˆG x, xA(B), ω
, which represents observed pressure, will be zero and as a result the integral over ∂_D0 will be zero too. This means that the right-hand side of equation 2.14 needs to be evaluated only over ∂Dm. Following the same steps as above, one can conclude that in the case when a free surface is present, to retrieve the Green’s function in the left-hand side of equation2.23, one needs to evaluate the right-hand side only over ∂Dm, i.e., one needs to have sources only in the subsurface. Intuitively, one can look at the free surface as a mirror that reflects the subsurface sources and creates an illusion of sources distributed over a closed boundary.

2.1.4 SI relation for transient sources in acoustic medium

In the previous subsection, the sources along the boundary ∂_{D were taken} to be impulsive. In practice, one will have to deal with band-limited sources characterized by a wavelet sˆ(x, ω). Then, the observed wavefields at the

points xAand xB in the two states are

ˆ

Substitution of equation 2.24 for states _{A and B into the SI relation} 2.23 gives 2ℜn ˆG(xA, xB, ω)o ˆS0(ω) ≈ 2 ρc I ∂D ˆ F(x, ω) ˆpobs∗(xA, x, ω) ˆpobs(xB, x, ω) d2x, (2.25) where ˆ F (x, ω) = Sˆ0(ω) ˆ S(x, ω) (2.26)

is a shaping filter with ˆS(x, ω) = ˆs∗(x, ω) ˆs(x, ω) the power spectrum of

the sources along ∂_{D and ˆ}S0(ω) an average, arbitrarily chosen power spec-trum. The SI relation 2.25 can be used when separate recordings can be taken at xA and xB from each of the transient sources on ∂D. The shaping filter ˆF (x, ω) compensates for the different power spectra of the sources.

This requires that these power spectra are known.

2.1.5 SI relation for noise sources in acoustic medium

It is not always possible to record separate responses from the sources on

∂D, as required for the application of SI relation 2.25. Often, the sources will be overlapping in time or one may not even know exactly when the sources were active. In such cases, there is an alternative solution if one may assume the sources to be mutually temporally uncorrelated with equal power spectrum ˆS(ω) and to be acting simultaneously. This means that at

the observation points xAand xB, the recorded wavefields will be

ˆ pobs(xA, ω) = I ∂D ˆ G(xA, x, ω) ˆN(x, ω) d2x, (2.27a) ˆ pobs(xB, ω) = I ∂D ˆ G(xB, x′, ω) ˆN(x′, ω) d2x′, (2.27b) where ˆN(x, ω) and ˆN(x′_{, ω) are the spectra of the sources at x and x}′_, respectively. Because the sources are assumed uncorrelated,

_ˆ

where_{h•i stands for spatial ensemble average. In such a case, the correlation} of the observed wavefields at xA and xBcan be written as

ˆpobs∗(xA, ω) ˆpobs(xB, ω) = I ∂D ˆ G∗(xA, x, ω) ˆG(xB, x, ω) ˆS(ω) d2x. (2.29) Comparing equation2.29to equation2.23, it can be concluded that

2ℜn ˆG(xA, xB, ω)o ˆS(ω) ≈

2 ρc ˆp

obs∗

(xA, ω) ˆpobs(xB, ω) . (2.30) The above relation is very easy to apply in practice – one only needs to install seismic receivers, record the noise responses, and then crosscorrelate the recorded responses following equation 2.30. A disadvantage is that no compensation can be applied for the different spectra of the different sources. As the sources are assumed uncorrelated, they can be called ”noise sources” and the recorded wavefields can be called ”noise” in the sense that one may not be able to see any clear arrivals. Nevertheless, what is recorded from such sources are propagating fields. The advantage of using relation2.30is that one does not need to have any a priori information about the sources. In practice, though, it will be difficult to find simultaneously acting noise sources. It will be most likely that the noise sources will be acting separately in time or will be partly overlapping in time. For this reason, to ensure the fulfillment of the assumption that the sources are uncorrelated in time, the noise recordings should be very long. In such a case, the ensemble average is then approximated by averaging over different time windows.

2.1.6 SI relations for impulsive sources in an elastic medium For an elastic medium, making use of the elastic-case column in Table2.1, the reciprocity relations2.5and2.9can be rewritten as

Z

D

n

−ˆτij,A(x, ω) ˆhij,B(x, ω) − ˆvi,A(x, ω) ˆfi,B(x, ω)

+ˆhij,A(x, ω) ˆτij,B(x, ω) + ˆfi,A(x, ω) ˆvi,B(x, ω)

o d3x₌ I

∂D{ˆv

and

Z

D

n

−ˆτij,A∗ (x, ω) ˆhij,B(x, ω) + ˆv∗i,A(x, ω) ˆfi,B(x, ω)

−ˆh∗ij,A(x, ω) ˆτij,B(x, ω) + ˆfi,A∗ (x, ω) ˆvi,B(x, ω)

o d3x₌ I

∂D−ˆv ∗

i,A(x, ω) ˆτij,B(x, ω) − ˆτij,A∗ (x, ω) ˆvi,B(x, ω) njd2x, (2.32) respectively.

Just like in the acoustic case, certain choices are made about the wavefields and the source fields in the states _{A and B. Here, the choice is for} impul-sive point sources of external volume force density ˆfi(x, ω) at points xA and xB in D for states A and B, respectively, while the external deforma-tion rate density ˆhij(x, ω) is taken to be zero. Then, the particle velocity

ˆ

vi(x, ω) at x in both states can be defined as the observed impulse responses (Green’s functions ˆGv,f_i,p(q) x_{, xA(B)}, ω
, see Table_{2.3). The expressions for}

the stresses observed at x for states_{A and B are obtained from the} substi-tution of the corresponding quantities for ˆfi(x, ω), ˆhij(x, ω), and ˆvi(x, ω) from Table2.3into the matrix-vector equation2.1for an elastic medium and the usage of the inverse of the compliance, which is the stiffness tensor cijkl obeying cijklsklmn= sijklcklmn = ₂1 (δimδjn+ δinδjm).

Substitution of the expressions in the middle and right columns of Table2.3 for their corresponding quantities in the convolution-type reciprocity theo-rem 2.31 and the use of the sifting property of the delta function, results in − ˆGv,f_q,p(xB, xA, ω) + ˆGv,fp,q(xA, xB, ω) = I ∂D n ˆ_Gv,f i,p (x, xA, ω) ˆGτ,fij,q(x, xB, ω) − ˆGτ,f_ij,p(x, xA, ω) ˆGv,fi,q (x, xB, ω) o njd2x. (2.33) The right-hand side of equation2.33represents a surface integral over prod-ucts of causal-in-time Green’s functions. If the surface of integration ∂D is

boundary, it should vanish for any choice of ∂_{D as long as the points x}Aand

x_{Bare inside it. This results in the elastodynamic source-receiver reciprocity}

relation

ˆ

Gv,fq,p (xB, xA, ω) = ˆGv,fp,q(xA, xB, ω) . (2.34) If, though, the choice is for impulsive point sources of external deformation rate density at the point xB in D for state B, while for state A the choice is for a source of external volume force density at xA, see Table 2.4, then the substitution of the parameters from Table2.4into the reciprocity relation 2.31followed by the same reasoning as above, will produce

ˆ

Gτ,fqr,p(xB, xA, ω) = ˆGv,hp,qr(xA, xB, ω) . (2.35) Let one now substitute the quantities from the middle and the right columns of Table 2.3 into their corresponding places in the correlation-type elasto-dynamic reciprocity relation2.32. Making use of the sifting property of the

Parameter State_A State_B

ˆ fi(x, ω) δ(x − xA) δip δ(x − xB) δiq ˆhij(x, ω) 0 0 ˆ vi(x, ω) Gˆv,fi,p (x, xA, ω) Gˆv,fi,q (x, xB, ω) ˆ τij(x, ω) _jω1 cijkl(x) ∂lGˆv,f_k,p(x, xA, ω) _jω1 cijkl(x) ∂lGˆv,f_k,q(x, xB, ω) , ˆGτ,f_ij,p(x, xA, ω) , ˆGτ,fij,q(x, xB, ω)

Table 2.3:Choice of wavefield and source parameters for the elastodynamic states

A and B, to be used in two-way wavefield elastodynamic reciprocity theorems. The

delta function gives n ˆ_Gv,f q,p(xB, xA, ω) o∗ + ˆGv,f_p,q(xA, xB, ω) = − I ∂D n ˆ_Gv,f i,p (x, xA, ω) ∗ ˆ Gτ,f_ij,q(x, xB, ω) + ˆGτ,f_ij,p(x, xA, ω) ∗ ˆ Gv,f_i,q (x, xB, ω) o njd2x. (2.36) Further application of the source-receiver reciprocities2.34and2.35to both sides of equation2.36results in

2ℜn ˆGv,f_p,q(xA, xB, ω) o = − I ∂D n ˆ_Gv,f p,i (xA, x, ω) ∗ ˆ Gv,h_q,ij(xB, x, ω) + ˆGv,h_p,ij(xA, x, ω) ∗ ˆ Gv,f_q,i (xB, x, ω) o njd2x. (2.37) In the time domain, the products on the right-hand side of the above equa-tion correspond to crosscorrelaequa-tions, while the left-hand side corresponds to a Green’s function and its time-reversed version. In this way, relation

Parameter State_A State_B

ˆ fi(x, ω) δ(x − xA) δip 0 ˆhij(x, ω) 0 δ(x − xB) δiqδjr ˆ vi(x, ω) Gˆv,fi,p (x, xA, ω) _jωρ1 ∆ijkl∂jGˆτ,hkl,qr(x, xB, ω) , ˆGv,h_i,qr(x, xB, ω) ˆ τij(x, ω) _jω1 cijkl(x) ∂lGˆv,f_k,p(x, xA, ω) Gˆτ,hij,qr(x, xB, ω) , ˆGτ,f_ij,p(x, xA, ω)

2.37provides the possibility to retrieve the complete Green’s function and its time-reversed version between the points xAand xB by constructive and de-structive interference (summation) of crosscorrelations of observed Green’s functions ˆGv,f_p,i x_{A(B), x, ω}
, due to impulsive force sources on ∂D, and

ob-served Green’s functions ˆGv,hq,ij xA(B), x, ω
nj due to impulsive deformation sources on ∂_{D. Equation}2.37is exact and valid for any lossless inhomoge-neous (inside, as well as outside ∂_{D) elastic medium.} Van Manen et al. [2006] developed an effective modelling scheme based on a relation very similar to equation2.37.

But just as in the acoustic case, in the application of this relation one faces the two difficulties of having to perform two correlations separately and needing the responses of two different types of sources. To overcome these problems, it is assumed that in a small region around ∂_{D the medium is homogeneous} and isotropic. Because of this, the observed wavefields can be written as sim-ple sum of observed P- and S-waves. Thus, the Green’s functions observed at x due to sources at xA and xB can be represented in terms of Green’s function potentials for observed P- and S-waves [Wapenaar and Berkhout, 1989;Aki and Richards,2002]:

ˆ

Gv,f_i,p(q) x_{, xA(B), ω
=} −1

jωρ n

∂iGˆφ,f_0,p(q) x, xA(B), ω
+ ǫijk∂jGˆφ,f_k,p(q) x, xA(B), ω

o

with ∂kGˆφ,f_k,p(q) x, xA(B), ω
= 0 (2.38) where the superscript φ means that the observed quantity is a P-wave when the left subscript is0, or an S-wave when the left subscript is k. The Green’s

function potentials obey the Helmholtz equations

∂i∂iGˆφ,f_0,p(q) x, xA(B), ω
+ ω2 c2_PGˆ φ,f 0,p(q) x, xA(B), ω
= 0 (2.39a) ∂i∂iGˆ φ,f k,p(q) x, xA(B), ω
+ ω2 c2 S ˆ Gφ,f_k,p(q) x_{, xA(B)}, ω
= 0 , _(2.39b)

terms of the observed Green’s functions as ˆ Gφ,f_0,p(q) x_{, xA(B), ω}
= −ρc 2 P jω ∂iGˆ v,f i,p(q) x, xA(B), ω
(2.40a) ˆ Gφ,f_k,p(q) x_{, xA(B)}, ω
= ρc 2 S jω ǫkji∂jGˆ v,f i,p(q) x, xA(B), ω
. (2.40b) Furthermore, at the boundary ∂_{D the Green’s function potentials can be} di-vided into inward- and outward-propagating parts:

ˆ Gφ,f_K,p(q) x_{, xA(B)}, ω
= ˆ_Gφ,f,in K,p(q) x, xA(B), ω
+ ˆG φ,f,out K,p(q) x, xA(B), ω
, (2.41) where K = 0, 1, 2, 3. Expressions for the Green’s functions ˆGτ,f_ij,p(x, xA, ω) and ˆGτ,f_ij,q(x, xB, ω) in terms of the Green’s function potentials can be written as ˆ Gτ,f_ij,p(q) x_{, xA(B)}, ω
= −1 (jω)2ρ n cijkl∂l∂kGˆφ,f_0,p(q) x, xA(B), ω
+cijkl∂lǫkmn∂mGˆφ,f_n,p(q) x, xA(B), ω
o . (2.42)

Substitution of equation 2.38 and equation2.42 in conjunction with equa-tion 2.41 for states _{A and B into the right-hand side of equation} 2.36 re-sults in products between inward- and outward-propagating P-waves and S-waves.Wapenaar and Berkhout[1989] show that for a horizontal boundary, products between only P-waves and products between only S-waves remain, while the cross-products between P- and S-waves cancel each other. They show further that only products between waves traveling in opposite direc-tions will remain. I.e., for a horizontal boundary ∂_{D, the right-hand side of} equation2.36can be written as

(note that the complex conjugation reverses the propagation direction, hence, all terms are indeed products of waves traveling in the opposite directions). The above result is not exact because the evanescent waves that propagate parallel to the horizontal boundary were neglected.

The integral will now be analyzed for an arbitrary closed surface ∂_{D. In the} high-frequency regime, the integral on the right-hand side of equation2.36 receives its main contributions from stationary points along the boundary

∂D. At such stationary points a local coordinate system is chosen, in which

the local x3 axis is in the direction of the outward-pointing normal on ∂D. For each of these local coordinate systems, one obtains a result similar to the result in2.43. When all the stationary points are accounted for, the total result would include also waves that propagate in the horizontal direction in the global coordinate system, including evanescent waves. As a result of this reasoning, going back to the global coordinate system, equation2.36can be written as 2ℜn ˆGv,f_p,q(xA, xB, ω) o = 2 jωρ I ∂D n ∂jGˆφ,f,in0,p (x, xA, ω) ∗ ˆ Gφ,f,in_0,q (x, xB, ω) +∂jGˆφ,f,out0,p (x, xA, ω) ∗ ˆ Gφ,f,out_0,q (x, xB, ω) +∂jGˆφ,f,ink,p (x, xA, ω) ∗ ˆ Gφ,f,in_k,q (x, xB, ω) +∂jGˆφ,f,outk,p (x, xA, ω) ∗ ˆ Gφ,f,out_k,q (x, xB, ω) o njd2x. (2.44) As the medium in a small region around ∂_{D is assumed homogeneous and} isotropic, the Green’s functions observed at xA and xB due to sources at

x can be represented in terms of Green’s function potentials from P- and

S-wave sources [Wapenaar and Berkhout,1989]:

Then, ˆ Gv,φ_p(q),0 x_{A(B), x, ω}
= −ρc 2 P jω ∂iGˆ v,f

p(q),i xA(B), x, ω
, (2.46a)

ˆ Gv,φ_p(q),k x_{A(B), x, ω
=} ρc 2 S jω ǫkji∂jGˆ v,f p(q),i xA(B), x, ω
. (2.46b) The use of reciprocity relations2.34together with equations2.40 and2.46, results in the reciprocity relation for the Green’s function potentials

ˆ

Gφ,f,in_K,p(q) x_{, xA(B)}, ω
= ˆ_Gv,φ,out

p(q),K xA(B), x, ω
, (2.47a)

ˆ

Gφ,f,out_K,p(q) x_{, xA(B)}, ω
= ˆ_Gv,φ,in

p(q),K xA(B), x, ω
, (2.47b) where K = 0, 1, 2, 3. Thus, relation2.44can be rewritten as

crosscorrelations on the right-hand side (see Figure3.30for examples). When the sources are randomly distributed in space, i.e., when the surface ∂_{D is} sufficiently irregular, the ghost correlation products will not interfere coher-ently and the ghost terms will be strongly weakened and can be ignored. At the same time, the retrieved Green’s function will contain all scattered fields from inside and outside ∂_{D (see Figure} 3.32). If the medium outside ∂_D is taken to be homogeneous with propagation velocity cP and cS for P- and S-waves, respectively, and mass density ρ, the outward-propagating part of the Green’s function potentials will not come back and then ghosts= 0.

In the high-frequency regime, the normal derivatives of each constituent of the Green’s function potentials can be written as

∂jGˆv,φ,in_p(q),K,w xA(B), x, ω
nj = −j ω cK |cos αK,w(x)| ˆG v,φ,in p(q),K,w xA(B), x, ω
(2.51a) ∂jGˆ v,φ,out p(q),K,w xA(B), x, ω
nj = +j ω cK |cos αK,w(x)| ˆG v,φ,out p(q),K,w xA(B), x, ω
, (2.51b) where cK = cP for K = 0 and cK = cS for K = 1, 2, 3. Further, if ∂D is taken to be a sphere with a sufficiently large radius, then all rays are normal to ∂_{D and consequently α = 0. Consequently, equation}2.49can finally be expressed as 2ℜn ˆGv,f_p,q(xA, xB, ω) o ≈ 2 ρcK I ∂D  ˆ_Gv,φ p,K(xA, x, ω) ∗ ˆ Gv,φ_q,K(xB, x, ω) d2x. (2.52)

The approximation in equation 2.52 involves only amplitude errors; when

∂D differs significantly from a sphere these errors may be significant.

2.1.7 SI relation for impulsive sources in an elastic medium with a free surface

Just as in the acoustic case, the medium_{D can be taken to be bounded (at} the top) by a free surface. Then, ∂_{D = ∂D}m ∪ ∂D0, where ∂D0 is the part coinciding with the free surface. The integral on the right-hand side of equation2.36needs then to be evaluated only over ∂_Dmbecause the traction

ˆ

Gτ,f_ij,p(q) x_{, xA(B), ω}
is zero on ∂D_{0. Following the same steps as above, one}

can conclude that in the case when a free surface is present, to retrieve the Green’s function in the left-hand side of equation2.52, one needs to evaluate the right-hand side only over sources in the subsurface.

2.1.8 SI relation for transient sources in an elastic medium In practice, the sources along the boundary ∂_{D will not be impulsive, as in} the previous subsection, but will be band-limited. Let the P- and S-wave sources be characterized by a wavelet with a spectrumsˆK_{(x, ω) for K =}

0, 1, 2, 3. Then the observed wavefields at the points xA and xB for the two states are ˆ vobs_p(q),K x_A(B), x, ω
= ˆ_Gv,φ p(q),k xA(B), x, ω
ˆs K (x, ω) . (2.53) Substitution of equation2.53for the two points into the SI relation2.52gives

2ℜn ˆGv,fp,q(xA, xB, ω)o ˆS0(ω) ≈ 2 ρcK I ∂D ˆ FK(x, ω) ˆvp,Kobs∗(xA, x, ω) ˆvobsq,K(xB, x, ω) d2x, (2.54) where ˆ FK(x, ω) = Sˆ0(ω) ˆ SK_{(x, ω)} (2.55) is a shaping filter with ˆSK(x, ω) = ˆsK∗(x, ω) ˆsK(x, ω) the power spectrum

2.1.9 SI relation for noise sources in an elastic medium To make use of the SI relation2.54, one needs separate recordings from each source-type at each source position on ∂_{D. Such information will be} avail-able only in special cases or after some preprocessing. When such recordings are not available, an alternative solution is to assume the sources to be act-ing simultaneously and to be mutually temporally uncorrelated with equal power spectrum cP

cKSˆ(ω). This means that at the observation points xAand

x_{B the recorded wavefields will be}

ˆ v_pobs(xA, ω) = I ∂D ˆ Gv,φ_p,K(xA, x, ω) ˆNK(x, ω) d2x, (2.56a) ˆ v_qobs(xB, ω) = I ∂D ˆ Gv,φ_q,L(xB, x′, ω) ˆNL(x′, ω) d2x′, (2.56b) where ˆNK(x, ω) and ˆNL(x′, ω) are the spectra of the different source types at x and x′, respectively. Because the sources are assumed uncorrelated for any K _{6= L and x 6= x}′_{, ˆ} N_K∗ (x, ω) ˆNL(x′, ω) = cP cKδKLδ(x − x ′_{) ˆS(ω) .} (2.57) Then, the correlation of the observed wavefields at xAand xB can be written as ˆvobs∗ p (xA, ω) ˆvqobs(xB, ω) = cP cK I ∂D  ˆ_Gv,φ p,K(xA, x, ω) ∗ ˆ Gv,φ_q,K(xB, x, ω) ˆS(ω) d2x. (2.58) Comparing equation2.58to equation2.52, it can be concluded that

2ℜn ˆGv,f_p,q (xA, xB, ω)o ˆS(ω) ≈

2 ρcP

ˆvobs∗

2.1.10 EMI relations for impulsive sources

In the electromagnetic case, making use of the right-most column in Table 2.1, the reciprocity relations 2.5 and 2.9 can be rewritten as [Slob et al., 2007] Z D n ˆ_{Er,A(x, ω) ˆJ}e r,B(x, ω) − ˆHj,A(x, ω) ˆJj,Bm (x, ω) − ˆJr,Ae (x, ω) ˆEr,B(x, ω) + ˆJj,Am (x, ω) ˆHj,B(x, ω) o d3x₌ I ∂D nmǫmrj ˆHj,A(x, ω) ˆEr,B(x, ω) − ˆEr,A(x, ω) ˆHj,B(x, ω)  d2x (2.60) and − Z D n ˆ_E∗ r,A(x, ω) ˆJ e r,B(x, ω) + ˆH ∗ j,A(x, ω) ˆJ m j,B(x, ω)

+ ˆJr,Ae∗ (x, ω) ˆEr,B(x, ω) + ˆJj,Am∗(x, ω) ˆHj,B(x, ω)

o d3x₌ I ∂D nmǫmrj ˆHj,A∗ (x, ω) ˆEr,B(x, ω) + ˆEr,A∗ (x, ω) ˆHj,B(x, ω)  d2x_, (2.61) respectively.

A choice is made for impulsive point source of external volume density of the electric current ˆJre(x, ω) at points xAand xBinD for statesA and B, respec-tively, while the external volume density of the magnetic current ˆJm

j (x, ω) is taken to be zero. In such a way, the electric and the magnetic field vectors can be expressed in terms of Green’s functions, see Table2.5.

Substituting the expressions in the middle and right columns of Table2.5for their corresponding quantities in the convolution-type reciprocity theorem 2.60, making use of the sifting property of the delta function, and applying the same arguments as the ones directly after equation2.33, yields

ˆ

GE,J_q,pe(xB, xA, ω) = ˆGE,J

e

p,q (xA, xB, ω) . (2.62) If, instead, a choice is made in state_{B for a source of external volume density} of the magnetic current ˆJm

as above would produce

ˆ

GH,J_q,p e(xB, xA, ω) = − ˆGE,J

m

p,q (xA, xB, ω) . (2.63) When the quantities from the middle and the right columns of Table2.5are substituted into their corresponding places in the correlation-type electro-magnetic reciprocity relation2.61, followed by the application of the sifting property of the delta function, results in

n ˆ_GE,Je q,p (xB, xA, ω) o∗ + ˆGE,J_p,q e(xA, xB, ω) = − I ∂D nmǫmrjn ˆGH,J e j,p (x, xA, ω) ∗ ˆ GE,J_r,q e(x, xB, ω) + ˆGE,Jr,p e(x, xA, ω) ∗ ˆ GH,J_j,q e(x, xB, ω) o d2x_{. (2.64)}

If this is followed by the application of the reciprocity relations 2.62 and

Parameter State_A State_B

ˆ J_ke_{(x, ω)} δ_{(x − x}A) δkp δ(x − xB) δkq ˆ J_jm(x, ω) 0 0 ˆ Er(x, ω) GˆE,J e r,p (x, xA, ω) GˆE,J e r,q (x, xB, ω) ˆ Hi(x, ω) −ζ_jωijǫjmr∂mGˆE,J e r,p (x, xA, ω) −ζ_jωijǫjmr∂mGˆE,J e r,q (x, xB, ω) , ˆGH,J_i,p e(x, xA, ω) , ˆGH,J e i,q (x, xB, ω)

Table 2.5:Choice of wavefield and source parameters for the electromagnetic states

A and B to be used in the two-way wavefield electromagnetic reciprocity theorems.

The superscripts in the Green’s function notation represent the receiver quantity (electric E or magnetic H field vector) and source quantity (electric Je or mag-netic Jm current), respectively, while the subscripts represent the components of the observed quantity (k, i) and the source components (p, q), respectively. ζij is

2.63, then the result is 2ℜn ˆGE,J_p,qe(xA, xB, ω) o = I ∂D nmǫmrjn ˆGE,J m p,j (xA, x, ω) ∗ ˆ GE,J_q,r e(xB, x, ω) + ˆGE,Jp,re(xA, x, ω) ∗ ˆ GE,Jq,j m(xB, x, ω) o d2x_{. (2.65)}

The above equation shows, that to retrieve the electric field at point xAdue to an electric current source at point xB, one needs to sum the crosscorrelations of observed electric field at xA and xB due to impulsive sources of electric and magnetic currents at points x along ∂_{D. Equation}2.65is exact and valid for any inhomogeneous (inside, as well as outside ∂_{D) lossless medium.} Just like in the acoustic and the elastic case, some simplifications of equation 2.65 are needed to make it easy for practical applications. The first simpli-fication stems from the desire to measure at xA and xB only the response from one source type. To achieve this, using Table 2.5, the right-hand side of equation2.64is rewritten only in terms of observed electrical fields, while

Parameter State_A State_B

ˆ J_ke(x, ω) δ_{(x − x}A) δkp 0 ˆ J_jm_{(x, ω)} 0 δ_{(x − x}B) δjq ˆ Er(x, ω) GˆE,J e r,p (x, xA, ω) ξ_jωrkǫkmi∂mGˆH,J m i,q (x, xB, ω) , ˆGE,Jr,q m(x, xB, ω) ˆ Hi(x, ω) −ζ_jωijǫjmr∂mGˆE,J e r,p (x, xA, ω) GˆH,J m i,q (x, xB, ω) , ˆGH,J_i,p e(x, xA, ω)

Table 2.6: Alternative choice of wavefield and source parameters for the elec-tromagnetic states _{A and B to be used in the two-way wavefield electromagnetic} convolution-type reciprocity theorem. ξrkis the inverse of the electric permittivity

the source-receiver reciprocity has been used on the left-hand side: 2ℜn ˆGE,J_p,qe(xA, xB, ω) o = − I ∂D nmǫmrj  −ζ_jωjiǫinl∂nGˆE,J e l,p (x, xA, ω) ∗ ˆ GE,J_r,q e(x, xB, ω) + ˆGE,J_r,p e(x, xA, ω) ∗ −ζ_ji jω ǫinl∂nGˆ E,Je l,q (x, xB, ω)  d2x_{. (2.66)}

Assuming that the medium in a small region around ∂_{D is homogeneous and} isotropic, the inverse of the magnetic permeability can be taken outside the integral (and switching back to magnetic permeability):

2ℜn ˆGE,J_p,qe(xA, xB, ω) o = − _jωµ1 I ∂D nmǫmrj n ǫjnl∂nGˆE,J e l,p (x, xA, ω) ∗ ˆ GE,J_r,q e(x, xB, ω) − ˆGE,J_r,p e(x, xA, ω) ∗ ǫjnl∂nGˆE,J e l,q (x, xB, ω) o d2x_{. (2.67)}

Using the fact that ǫmrjǫjnl = δmnδrl− δmlδrn, equation2.67becomes

2ℜn ˆGE,Jp,qe(xA, xB, ω) o = − _jωµ1 I ∂D n nm∂mGˆE,J e r,p (x, xA, ω) ∗ ˆ GE,Jr,q e(x, xB, ω) − ˆGE,J_r,p e(x, xA, ω) ∗ nm∂mGˆE,J e r,q (x, xB, ω) o d2x − _jωµ1 I ∂D n nm ˆGE,J e r,p (x, xA, ω) ∗ ∂rGˆE,J e m,q (x, xB, ω) −nm  ∂rGˆE,J e m,p (x, xA, ω) ∗ ˆ GE,J_r,q e(x, xB, ω) o d2x_{. (2.68)}

Because the electric permittivity is constant in a small region around ∂_{D, the} electric field is divergence free, i.e., ∂rGˆE,J

e

that one can write 2ℜn ˆGE,J_p,qe(xA, xB, ω) o = − _jωµ1 I ∂D n nm∂mGˆE,J e r,p (x, xA, ω) ∗ ˆ GE,J_r,q e(x, xB, ω) − ˆGE,J_r,p e(x, xA, ω) ∗ nm∂mGˆE,J e r,q (x, xB, ω) o d2x − _jωµ1 I ∂D nm∂rn ˆGE,J e r,p (x, xA, ω) ∗ ˆ GE,J_m,qe(x, xB, ω) − ˆGE,Jm,pe(x, xA, ω) ∗ ˆ GE,Jr,q e(x, xB, ω) o d2x_{. (2.69)}

The second integral on the right-hand side of the above equation has a van-ishing contribution. This can be shown by decomposing the partial derivative in its normal and tangential components:

∂r = nr(nj∂j) − ǫmijniǫjklnk∂l. (2.70) Making use of this relation, the second integral on the right-hand side of equation2.69can be written as

1 jωµ I ∂D nm∂rn ˆGE,J e r,p (x, xA, ω) ∗ ˆ GE,J_m,qe(x, xB, ω) − ˆGE,J_m,pe(x, xA, ω) ∗ ˆ GE,J_r,q e(x, xB, ω) o d2x₌ 1 jωµ I ∂D nmnr(nj∂j)n ˆGE,J e r,p (x, xA, ω) ∗ ˆ GE,J_m,qe(x, xB, ω) − ˆGE,Jm,pe(x, xA, ω) ∗ ˆ GE,Jr,q e(x, xB, ω) o d2x − _jωµ1 I ∂D nmǫmijniǫjklnk∂ln ˆGE,J e r,p (x, xA, ω) ∗ ˆ GE,Jm,qe(x, xB, ω) − ˆGE,J_m,pe(x, xA, ω) ∗ ˆ GE,J_r,q e(x, xB, ω) o d2x_{. (2.71)}

The first integral on the right-hand side of equation2.71is zero, which can be verified by directly expanding the notation after substituting the subscripts

is split into two parts ∂_D1 and ∂D2, which share a common bounding curve

L; the unit tangent vectors τ and −τ of L point in opposite directions. Then, 1 jωµ I ∂D nm∂rn ˆGE,J e r,p (x, xA, ω) ∗ ˆ GE,J_m,qe(x, xB, ω) − ˆGE,J_m,pe(x, xA, ω) ∗ ˆ GE,J_r,q e(x, xB, ω) o d2x₌ −_jωµ1 Z ∂D1 nmǫrijniǫjklnk∂ln ˆGE,J e r,p (x, xA, ω) ∗ ˆ GE,J_m,qe(x, xB, ω) − ˆGE,J_m,pe(x, xA, ω) ∗ ˆ GE,J_r,q e(x, xB, ω) o d2x −_jωµ1 Z ∂D2 nmǫrijniǫjklnk∂ln ˆGE,J e r,p (x, xA, ω) ∗ ˆ GE,J_m,qe(x, xB, ω) − ˆGE,Jm,pe(x, xA, ω) ∗ ˆ GE,Jr,q e(x, xB, ω) o d2x_{. (2.72)}

Application of Stoke’s curl theorem to the right-hand side of the above equa-tion produces 1 jωµ I ∂D nm∂rn ˆGE,J e r,p (x, xA, ω) ∗ ˆ GE,J_m,qe(x, xB, ω) − ˆGE,J_m,pe(x, xA, ω) ∗ ˆ GE,J_r,q e(x, xB, ω) o d2x₌ − nm 1 jωµ I L ǫmijniτjn ˆGE,J e r,p (x, xA, ω) ∗ ˆ GE,J_m,qe(x, xB, ω) − ˆGE,Jm,pe(x, xA, ω) ∗ ˆ GE,Jr,q e(x, xB, ω) o dx + 1 jωµ I L nmǫmijniτjn ˆGE,J e r,p (x, xA, ω) ∗ ˆ GE,Jm,qe(x, xB, ω) − ˆGE,J_m,pe(x, xA, ω) ∗ ˆ GE,J_r,q e(x, xB, ω) o dx, (2.73)

for each of the line integrals. In this way, equation2.69finally becomes 2ℜn ˆGE,J_p,qe(xA, xB, ω) o = − _jωµ1 I ∂D n nm∂mGˆE,J e p,r (xA, x, ω) ∗ ˆ GE,J_q,r e(xB, x, ω) − ˆGE,J_p,re(xA, x, ω) ∗ nm∂mGˆE,J e q,r (xB, x, ω) o d2x_{, (2.74)}

where the reciprocity relation2.62was used.

Following a similar analysis path as the one for equations2.42till2.51, i.e., taking the medium outside ∂_{D to be homogeneous and isotropic, making a} high-frequency approximation, and assuming ∂_{D to be a sphere with} suffi-ciently large radius, one can finally write

2ℜn ˆGE,J_p,qe(xA, xB, ω) o ≈ − _cµ2 I ∂D  ˆ_GE,Je p,r (xA, x, ω) ∗ ˆ GE,J_q,r e(xB, x, ω) d2x, (2.75) where c = _εµ1 1 2

is the propagation velocity outside ∂_{D. The} approxi-mation in the above equation refers only to amplitude errors, but when ∂_D differs significantly from a sphere these errors might be significant. Never-theless, the application of relation2.75 will retrieve correctly the phases of all arrivals. That is why relation 2.75is considered acceptable for EMI for practical applications.

2.1.11 EMI relations for transient sources

In the previous subsection, the sources along the boundary ∂_{D were assumed} impulsive. In practice, the sources will rather be band-limited. Consider electric current sources characterized by a wavelet with a spectrumˆsj(x, ω)

for j = 1, 2, 3. In this case, at the points xA and xB the observed electric fields for both states would be

ˆ

E_p(q),jobs x_A(B), x, ω
= ˆ_GE,Je

Substituting equation2.76into the EMI relation2.75results in 2ℜn ˆGE,J_p,qe(xA, xB, ω)o ˆS0(ω) ≈ − _cµ2 I ∂D ˆ Fj(x, ω) ˆE_p,jobs(xA, x, ω) ∗ ˆ E_q,jobs(xB, x, ω) d2x. (2.77) In the above equation,

ˆ

Fj(x, ω) = Sˆ0(ω) ˆ

Sj_{(x, ω)} (2.78) is a shaping filter with ˆSj(x, ω) the power spectrum of the electric current

sources along ∂_{D and ˆ}S0(ω) an average, arbitrary chosen power spectrum. The EMI relation2.77can be applied in practice if separate recordings can be made from each of the transient sources on ∂_{D. The shaping filter ˆ}Fj(x, ω)

compensates for the known different power spectra of the sources. 2.1.12 EMI relations for noise sources

It is not always possible to have separate recordings from the electrical cur-rent sources along ∂_{D. To overcome this problem, alternatively one can} assume the sources to be mutually uncorrelated with an equal power spec-trum ˆS(ω). For simultaneously acting sources, at the observation points xA and xB the recorded electric fields will be

ˆ

Epobs(xA, ω) =

I

∂D

ˆ

GE,Jp,i e(xA, x, ω) ˆNi(x, ω) d2x (2.79a)

ˆ E_qobs(xB, ω) = I ∂D ˆ GE,J_q,j e(xB, x′, ω) ˆNj(x′, ω) d2x′, (2.79b) where ˆNi(x, ω) and ˆNj(x′, ω) represent the spectra of the electric source at

x and x′_{, respectively. Because the sources are assumed uncorrelated for any}

i_{6= j and x 6= x}′, it follows that

Using the above property, the correlation of the observed electric fields at xA and xB is ˆ E_pobs∗(xA, ω) ˆEqobs(xB, ω) = I ∂D 2 cµ ˆG E,Je p,j (xA, x, ω) ∗ ˆ GE,Jq,j e(xB, x, ω) ˆS(ω) d2x. (2.81) Comparing equation2.81to equation2.75, it can be written that

2ℜn ˆGE,J_p,q e(xA, xB, ω)o ˆS(ω) ≈ −

ˆ

E_pobs∗(xA, ω) ˆEqobs(xB, ω), (2.82) which means that the crosscorrelation of the observed electromagnetic noise at two points will retrieve the electromagnetic Green’s function between these two points convolved with the power spectrum of the noise.

2.2

SI relations based on a one-way wavefield

re-ciprocity

2.2.1 Reciprocity theorems

Consider a 3D inhomogeneous, lossless medium. In this medium, one choo-ses a domain _{D with mass density ρ (x) and compressibility κ (x), if the} medium is acoustic, or with compliance sijkl(x) in case of an elastic medium; domain _{D contains no discontinuities in the medium parameters. This} do-main is bounded by ∂_{D with an outward-pointing normal vector n. The} boundary consists of two parallel planes, which stretch to infinity and are horizontal. The top plane is ∂_D0 and the bottom one is ∂Dm. The x3 axes points downwards.

Choosing the preferred direction of propagation to be the vertical direction (along the x3 axis), the one-way wavefield equation can be written in the space-frequency domain in matrix-vector form as

∂3Pˆ (x, ω) −Bˆ (x, ω)Pˆ (x, ω) =Sˆ (x, ω) . (2.83) In the above equation,Pˆ (x, ω) represents the flux-normalized one-way wave

quantities andSˆ(x, ω) represents the flux-normalized one-way source

Parameter Acoustic medium Elastic medium ˆ P(x, ω) ˆP +_{(x, ω)} ˆ P−_{(x, ω)} !           ˆ Φ+_{(x, ω)} ˆ Ψ+ x1(x, ω) ˆ Ψ+ x2(x, ω) ˆ Φ−_{(x, ω)} ˆ Ψ− x1(x, ω) ˆ Ψ− x2(x, ω)           ˆ S(x, ω) ˆS +_{(x, ω)} ˆ S−(x, ω) !           ˆ S_Φ+(x, ω) ˆ S_Ψx+ 1(x, ω) ˆ S_Ψx+ 2(x, ω) ˆ S_Φ−(x, ω) ˆ S_Ψx− ₁(x, ω) ˆ S_Ψx− 2(x, ω)          

Table 2.7:Expansion of the variables in equation2.83in case of an acoustic and an elastic medium. The superscripts_{+ and − denote downgoing and upgoing} propaga-tion, respectively, ˆP±stands for one-way wavefield, ˆS±- for one-way source field,

ˆ

Φ is the P-wave potential , ˆΨx1and ˆΨx2are the S-wave potentials with polarizations

in the x2x3- and x1x3-plane, respectively, ˆSΦ is the source P-wave potential, and

ˆ

S_Ψx1 and ˆSΨx2 are the source S-wave potentials with x2x3- and x1x3-polarization,

and an elastic medium. Bˆ (x, ω) is a one-way pseudo-differential operator

matrix defined as ˆ

B_{(x, ω) = −jω}Λˆ_{(x, ω) + ˆ}Θ_{(x, ω) , (2.84)}

where the first term_−jωΛˆ_{(x, ω) accounts for the vertical propagation of the}

one-way wavefield, while the operator matrix

ˆ

Θ_{(x, ω) = −L}ˆ−1_{(x, ω) ∂3L}ˆ_{(x, ω) (2.85)}

accounts for the vertical scattering of the one-way wavefield due to vari-ations of the medium parameters in the vertical direction [Corones et al., 1983;Fishman et al.,1987;Wapenaar and Berkhout, 1989;de Hoop, 1992; Wapenaar and Grimbergen, 1996;Haines and de Hoop, 1996]. In equation 2.85,Lˆ(x, ω) is a composition operator turning the one-way wavefields into

two-way wavefields and Lˆ−1(x, ω) is a decomposition operator turning the

two-way wavefields into one-way wavefields.

Two independent states _{A and B are taken in the domain D. State A is} described in the space-frequency domain by the one-way wavefield quantity

ˆ

PA(x, ω), by the one-way source quantity SˆA(x, ω), and by the operator matrix BˆA(x, ω). State B is characterized by the same variables, but the corresponding subscripts_{A are substituted by the subscript B.}

The one-way interaction quantity describing the local convolution of the wavefields in the two states, which is the counterpart of the two-way acoustic interaction quantity2.2, is ∂3 n ˆ PT_A(x, ω)N ˆPB(x, ω) o =n∂3PˆTA(x, ω) o N ˆPB(x, ω) +PˆTA(x, ω)N n ∂3PˆB(x, ω) o , (2.86) whereN_{= antidiag (1, −1).}

Equation2.87represents the local form of the one-way wavefield convolution-type reciprocity theorem [Wapenaar and Grimbergen,1996].

Integration of equation 2.87 over the domain _{D and application of Gauss’} divergence theorem to the left-hand side of equation2.87gives

Z D n ˆ PAT (x, ω)BˆTA(x, ω)N ˆPB(x, ω) +PˆTA(x, ω)N ˆBB(x, ω)PˆB(x, ω) +PˆT_A(x, ω)N ˆSB(x, ω) +SˆTA(x, ω)N ˆPB(x, ω) o d3x₌ I ∂D n ˆ PT_A(x, ω)N ˆPB(x, ω) o n3d2x. (2.88) Equation2.88represents the global form of the one-way wavefield convolu-tion-type reciprocity theorem.

Making use of the symmetry property ˆ

BTA(x, ω)N= −N ˆBA(x, ω) (2.89) and assuming the medium parameters in both states to be the same, which means that the one-way pseudo-differential operator matrices for the two states are equal BˆA(x, ω) = BˆB(x, ω), then the first two terms in the left-hand side of equation2.88cancel and the convolution-type reciprocity rela-tion simplifies to Z D n ˆ PAT (x, ω)N ˆSB(x, ω) +SˆTA(x, ω)N ˆPB(x, ω) o d3x₌ I ∂D n ˆ PT_A(x, ω)N ˆPB(x, ω) o n3d2x. (2.90) Instead of equation2.86, as a starting point can be taken the interaction quan-tity describing the local correlation of the wavefields in the two state:

∂3 n ˆ P†_A(x, ω)J ˆPB(x, ω) o =n∂3Pˆ†A(x, ω) o J ˆPB(x, ω) +Pˆ†_A(x, ω)Jn∂3PˆB(x, ω) o , (2.91)