Classical acoustic wave-field representations consist of volume and boundary integrals, of which the integrands contain specific combinations of Green's functions, source distributions, and wave fields. Using a unified matrix-vector wave equation for different wave phenomena, these representations can be reformulated in terms of Green's matrices, source vectors, and wave-field vectors. The matrix-vector formalism also allows the formulation of representations in which propagator matrices replace the Green's matrices. These propagator matrices, in turn, can be expressed in terms of Marchenko-type focusing functions. An advantage of the representations with propagator matrices and focusing functions is that the boundary integrals in these representations are limited to a single open boundary. This makes these representations a suitable basis for developing advanced inverse scattering, imaging and monitoring methods for wave fields acquired on a single boundary.

The aim of this paper is to give a systematic treatment of different types of wave-field representations (with Green's functions, propagator matrices, and Marchenko-type focusing functions), discuss their mutual relations, and indicate some new applications.

  • Representations with Green's functions. A Green's function is the response of a medium to an impulsive point source. It is named after George Green, who, in a privately published essay,1 introduced the use of impulse responses in field representations.2 Wave-field representations with Green's functions have been formulated, among others, for optics,3 acoustics,4,5 elastodynamics,6–9 and electromagnetics.10–12 They find numerous applications in forward modeling problems,13–15 inverse source problems,12,16 inverse scattering problems,5,17–19 imaging,20–25 time-reversal acoustics,26 and Green's function retrieval from ambient noise.27–29 

  • Representations with propagator matrices. In elastodynamic wave theory, a matrix formalism has been introduced to describe the propagation of waves in laterally invariant layered media.30,31 This formalism was refined by Gilbert and Backus,32 who coined the name propagator matrix. In essence, a propagator matrix “propagates” a wave field (represented as a vectorial quantity) from one plane in space to another. Using perturbation theory, the connection between wave-field representations with Green's functions and the propagator matrix formalism was discussed.33 The propagator matrix has been extended for laterally varying isotropic and anisotropic layered media.34,35 Propagation invariants for laterally varying layered media have been introduced, and the propagator matrix concept has been proposed for the modeling of reflection and transmission responses.36–39 The propagator matrix has also been used in a seismic imaging method, which accounts for multiple scattering in a model-driven way.40 

  • Representations with Marchenko-type focusing functions. Building on a one-dimensional (1D) acoustic autofocusing method, it has been shown that the wave field inside a laterally invariant layered medium can be retrieved with the Marchenko method from the single-sided reflection response at the surface of the medium.41–43 This concept was extended to a three-dimensional (3D) Marchenko wave-field retrieval method for laterally varying media.44 Central in the 3D Marchenko method are wave-field representations containing focusing functions. These representations have found applications in imaging methods45–47 and inverse source problems,48 which account for multiple scattering in a data-driven way.

The setup of this paper is as follows. In Sec. II, we briefly review the matrix-vector wave equation for laterally varying media,34,35 generalized for different wave phenomena, and we briefly discuss the concept of the Green's matrix, propagator matrix, and Marchenko-type focusing function. The symmetry properties of the matrix-vector wave equation allow the formulation of unified matrix-vector wave-field reciprocity theorems,49,50 which are reviewed in Sec. III. These reciprocity theorems form the basis for a systematic treatment of the different types of wave-field representations mentioned above. Traditionally, a wave-field representation is obtained by replacing one of the states in a reciprocity theorem by a Green's state. In Sec. IV, we follow this approach for the matrix-vector reciprocity theorems. By replacing one of the wave-field vectors by the Green's matrix (and the source vector by a unit source matrix), we obtain wave-field representations with Green's matrices. Analogous to this, in Sec. V, we replace one of the wave-field vectors in the reciprocity theorems by the propagator matrix and, thus, obtain wave-field representations with propagator matrices. In Sec. VI, we discuss a mixed form, obtained by replacing one of the wave-field vectors by the Green's matrix and the other by the propagator matrix. In Sec. VII, we discuss the relation between the propagator matrix and Marchenko-type focusing functions and use this relation to derive wave-field representations with focusing functions. We end with conclusions in Sec. VIII.

We use a unified matrix-vector wave equation as the basis for the derivations in this paper. In the space-frequency domain, it has the following form:32–36 

3qAq=d,
(1)

with

q=(q1q2),d=(d1d2),A=(A11A12A21A22).
(2)

Here, q(x,ω) is a space- and frequency-dependent N×1 wave-field vector, where x denotes the Cartesian coordinate vector (x1,x2,x3) (with the positive x3-axis pointing downward) and ω is the angular frequency. The N/2×1 sub-vectors, q1(x,ω) and q2(x,ω), contain wave-field quantities, which are specified for different wave phenomena in Table I. Operator 3 stands for the differential operator /x3. Matrix A(x,ω) is an N × N operator matrix; it contains the space- and frequency-dependent anisotropic medium parameters and horizontal differential operators 1 and 2. The definitions of this operator matrix for different wave phenomena can be found in many of the references mentioned in Sec. I. The N×1 vector d(x,ω) contains the space- and frequency-dependent source functions. A comprehensive overview of the operator matrices and source vectors for the wave phenomena considered in Table I (with some minor modifications) is given in Ref. 51.

TABLE I.

The wave-field sub-vectors q1(x,ω) and q2(x,ω) for the different wave phenomena. For acoustic waves in stationary fluids, p and v3 stand for the acoustic pressure and the vertical component of the particle velocity, respectively. For electromagnetic waves, Eα and Hα (α=1,2) are the horizontal components of the electric and magnetic field strength, respectively. For elastodynamic waves, vk and τk3 (k = 1,2,3) are the particle velocity and traction components, respectively. The same quantities appear in the vectors for the other wave phenomena, where, for poroelastodynamic and seismoelectric waves, the quantities are averaged in the bulk, fluid, or solid as indicated by the superscripts b, f, and s, respectively. Finally, ϕ denotes the porosity.

Nq1q2
Acoustic p v3 
Electromagnetic E0=(E1E2) H0=(H2H1) 
Elastodynamic v=(v1v2v3) τ3=(τ13τ23τ33) 
Poroelastodynamic (vsϕ(v3fv3s)) (τ3bpf) 
Piezoelectric 10 (vH0) (τ3E0) 
Seismoelectric 12 (vsϕ(v3fv3s)H0) (τ3bpfE0) 
Nq1q2
Acoustic p v3 
Electromagnetic E0=(E1E2) H0=(H2H1) 
Elastodynamic v=(v1v2v3) τ3=(τ13τ23τ33) 
Poroelastodynamic (vsϕ(v3fv3s)) (τ3bpf) 
Piezoelectric 10 (vH0) (τ3E0) 
Seismoelectric 12 (vsϕ(v3fv3s)H0) (τ3bpfE0) 

For the wave phenomena considered in Table I, the operator matrix A obeys the symmetry properties

AtN=NA,
(3)
AK=KA¯,
(4)
A*J=JA¯,
(5)

with

N=(OIIO),K=(OIIO),J=(IOOI),
(6)

where O and I are the N/2×N/2 zero and identity matrices, respectively. The superscript t denotes transposition (meaning that the matrix is transposed and the operators in the matrix are also transposed with 1t=1 and 2t=2), The “*” denotes complex conjugation, and “” denotes transposition and complex conjugation. In general, the medium parameters in A are complex valued and frequency dependent, accounting for losses. The bar above a quantity means that this quantity is defined in the adjoint medium. Hence, if A is defined in a lossy medium, then A¯ is defined in an effectual medium and vice versa52 (a wave propagating through an effectual medium gains energy). For lossless media, the bar can be dropped. For the wave phenomena considered in Table I, the power-flux density j in the x3-direction is related to the sub-vectors q1 and q2 according to

j=14(q1q2+q2q1).
(7)

As a special case, we consider acoustic waves in an inhomogeneous stationary medium with complex-valued and frequency-dependent compressibility κ(x,ω) and mass density ρkl(x,ω). The latter is defined as a tensor to account for effective anisotropy, for example, as a result of fine layering at the micro scale.53 The mass density tensor is symmetric, that is, ρkl(x,ω)=ρlk(x,ω). We introduce the inverse of the mass density tensor, the specific volume tensor ϑkl(x,ω), via ϑklρlm=δkm. Einstein's summation convention holds for repeated subscripts (unless otherwise noted); Latin subscripts run from 1 to 3 and Greek subscripts run from 1 to 2. For the acoustic situation, the vectors and matrix in Eq. (1) are given by

q=(pv3),d=(ϑ331ϑ3lfl1iωα(bαβfβ)+q),
(8)
A=(ϑ331ϑ3ββiωϑ331iωκ1iωαbαββαϑα3ϑ331),
(9)

with

bαβ=ϑαβϑα3ϑ331ϑ3β,
(10)

where i is the imaginary unit and q(x,ω) and fl(x,ω) are sources in terms of the volume-injection rate density and external force density, respectively.

Finally, for an isotropic medium, using ϑkl=ρ1δkl (where ρ is the scalar mass density), we obtain for the source vector and operator matrix54–58 

d=(f31iωα(1ρfα)+q),A=(0iωρiωκ1iωα1ρα0).
(11)

In the space-time domain, a Green's function is the response to an impulsive point source with the impulse defined as δ(t). The Fourier transform of δ(t) equals one; hence, in the space-frequency domain, the Green's function is the response to a point source with unit amplitude for all frequencies. We introduce the N × N Green's matrix G(x,xA,ω) for an unbounded, arbitrary inhomogeneous, anisotropic medium as the solution of

3GAG=Iδ(xxA),
(12)

where xA=(x1,A,x2,A,x3,A) defines the position of the point source, and I is an N × N identity matrix. Here, I has a size that is different from that in Eq. (6). For simplicity, we use one notation for differently sized identity matrices (the size always follows from the context). Also, for the zero matrix O, we use a single notation for differently sized matrices. Similar to the operator matrix A, the Green's matrix is partitioned as

G(x,xA,ω)=(G11G12G21G22)(x,xA,ω).
(13)

Equation (12) does not have a unique solution. To specify a unique solution, we demand that the time-domain Green's function G(x,xA,t) is causal, hence,

G(x,xA,t<0)=O.
(14)

This condition implies that G is outward propagating for |xxA|.

The simplest representation involving the Green's matrix is obtained when q and G reside in the same medium throughout space and both are outward propagating for |xxA|. Whereas q(x,ω) is the response to a source distribution d(x,ω) [Eq. (1)], G(x,xA,ω) is the response to a point source Iδ(xxA) for an arbitray source position xA [Eq. (12)]. Because Eq. (1) and (12) are linear, a representation for q(x,ω) follows by applying the superposition principle, according to

q(x,ω)=3G(x,xA,ω)d(xA,ω)d3xA,
(15)

where is the set of real numbers. This representation is a special case of the more general representations with the Green's matrices, which are derived in a more formal way in Sec. IV.

Next, we discuss the 2 × 2 acoustic Green's matrix as a special case of the N × N Green's matrix. For this situation, G is partitioned as

G(x,xA,ω)=(Gp,fGp,qGv,fGv,q)(x,xA,ω).
(16)

Here, the first superscript (p or v) refers to the observed wave-field quantity at x (the acoustic pressure or vertical component of the particle velocity), whereas the second superscript (f or q) refers to the source type at xA (the vertical component of the force or volume-injection rate). The unit of a specific element of the Green's matrix is the ratio of the units of the observed wave-field quantity at x and the source quantity at xA. For example, [Gp,f(x,xA,ω)]=[p]/[f]=m2. For an isotropic medium, all elements can be expressed in terms of the upper-right element as

Gv,q(x,xA,ω)=1iωρ(x,ω)3Gp,q(x,xA,ω),
(17)
Gp,f(x,xA,ω)=1iωρ(xA,ω)3,AGp,q(x,xA,ω),
(18)
Gv,f(x,xA,ω)=1iωρ(x,ω)(3Gp,f(x,xA,ω)δ(xxA)).
(19)

Here, 3,A stands for the differentiation with respect to the source coordinate x3,A. Equations (17) and (19) follow directly from Eqs. (11), (12), and (16). Equation (18) follows from Eq. (17) and a source-receiver reciprocity relation, which is derived in Sec. IV A.

We illustrate Gp,q, decomposed into plane waves, for a horizontally layered lossless isotropic medium. To this end, we first define the spatial Fourier transform of a space- and frequency-dependent function u(x,ω) along the horizontal coordinates xH=(x1,x2) according to

ũ(s,x3,ω)=2exp{iωsxH}u(xH,x3,ω)d2xH,
(20)

where s=(s1,s2), and s1 and s2 are the horizontal slownesses. This transform decomposes function u(x,ω) at a given depth level x3 into monochromatic plane wave components. Next, we define the inverse temporal Fourier transform per slowness value s as

u(s,x3,τ)=1π0ũ(s,x3,ω)exp{iωτ}dω,
(21)

where denotes the real part and τ is the intercept time.59 We apply these transforms to the Green's function Gp,q(x,xA,ω) and source function δ(xxA), choosing xA=(0,0,x3,0) and setting s2=0 (the field is cylindrically symmetric in the considered horizontally layered isotropic medium). We, thus, obtain Gp,q(s1,x3,x3,0,τ), which is the plane wave response (as a function of x3 and τ) to a source function δ(x3x3,0)δ(τ).

An example of a horizontally layered medium is shown in Fig. 1(a). The propagation velocities in the layers are indicated by cn (with c=1/κρ). The depth level of the source is chosen as x3,0=0 m; see Fig. 1(b). The vertical green line in Fig. 1 is the line τ = 0, left of which the field is zero due to the causality condition [similar to that in Eq. (14)]. Figure 1(b) further shows the numerically modelled Green's function Gp,q(s1,x3,x3,0,τ) as a function of x3 and τ for a single horizontal slowness s1=1/3500 s/m (the red and blue arrows indicate the downgoing and upgoing waves, respectively). This is the wave field that would be measured by a series of acoustic pressure receivers, vertically above and below the volume-injection rate source at x3,0=0 m. Each trace shows Gp,q(s1,x3,x3,0,τ) for a specific depth x3 as a function of τ (actually, at each depth, the Green's function has been convolved with a time-symmetric wavelet with a central frequency of 50 Hz to get a nicer display).

FIG. 1.

(Color online) The (a) horizontally layered medium and (b) Green's function Gp,q(s1,x3,x3,0,τ) (fixed s1), convolved with a wavelet, are shown.

FIG. 1.

(Color online) The (a) horizontally layered medium and (b) Green's function Gp,q(s1,x3,x3,0,τ) (fixed s1), convolved with a wavelet, are shown.

Close modal

The horizontal slowness s1 is related to the propagation angle αn in layer n via sinαn=s1cn, hence, in layer 1 (with c1=1600 m/s), the propagation angle is α1=27.2°. In the thin layer, with c3=3600 m/s, we obtain sinα3=1.03, meaning that α3 is complex valued. This implies that the wave is evanescent in this layer. Because the layer is thin, the wave tunnels through the layer and continues with a lower amplitude as a downgoing wave in the lower half-space.

For a lossless medium, a homogeneous Green's function is the superposition of a Green's function and its complex conjugate (or, in the time domain, its time-reversed version). The superposition is chosen such that the source terms of the two functions cancel each other, hence, a homogeneous Green's function obeys a wave equation without a source term.19,20 Here, we extend this concept for the matrix-vector wave equation for a medium with losses.

Let G(x,xA,ω) be, again, the outward propagating solution of Eq. (12) for an unbounded, arbitrary inhomogeneous, anisotropic medium. We introduce the Green's matrix of the adjoint medium, G¯(x,xA,ω), as the outward propagating solution of

3G¯A¯G¯=Iδ(xxA).
(22)

Pre- and post-multiplying all terms by J and, subsequently, using Eq. (5) and JJ=I gives

3JG¯JA*JG¯J=Iδ(xxA).
(23)

Subtracting the complex conjugate of all terms in Eq. (23) from the corresponding terms in Eq. (12) yields

3GhAGh=O,
(24)

with

Gh(x,xA,ω)=G(x,xA,ω)JG¯*(x,xA,ω)J.
(25)

Because Gh(x,xA,ω) obeys a wave equation without a source term, we call it the homogeneous Green's matrix. According to Eqs. (6), (13), and (25), it is partitioned as

Gh(x,xA,ω)=({G11G¯11*}{G12+G¯12*}{G21+G¯21*}{G22G¯22*})(x,xA,ω).
(26)

We introduce the N × N propagator matrix W(x,xA,ω) for an arbitrary inhomogeneous anisotropic medium as the solution of the unified matrix-vector wave equation (1) but without the source vector d. Hence,

3WAW=O.
(27)

Similar to the operator matrix A and Green's matrix G, the propagator matrix is partitioned as

W(x,xA,ω)=(W11W12W21W22)(x,xA,ω).
(28)

Equation (27) does not have a unique solution. To specify a unique solution, we impose the boundary condition

W(x,xA,ω)|x3=x3,A=Iδ(xHxH,A),
(29)

where xH,A denotes the horizontal coordinates of xA, hence, xH,A=(x1,A,x2,A). Because Eq. (27) is first order in 3, a single boundary condition suffices. Note that W(x,xA,ω) only depends on the medium parameters between the depth levels x3,A and x3. This is different from the Green's matrix G(x,xA,ω), which, for an arbitrary inhomogeneous medium, depends on the entire medium [this is easily understood from the illustration in Fig. 1(b)].

The simplest representation involving the propagator matrix is obtained when q and W reside in the same medium in the region between x3,A and x3 and both have no sources in this region. Whereas q(x,ω) obeys no boundary conditions in this region, W(x,xA,ω) collapses to Iδ(xHxH,A) at the depth level x3,A [Eq. (29)]. Applying the superposition principle again yields

q(x,ω)=DAW(x,xA,ω)q(xA,ω)d2xA,
(30)

where DA is the horizontal boundary defined as x3=x3,A.32–36 Note that W(x,xA,ω) propagates the field vector q from the depth level x3,A to x3, hence, the name “propagator matrix.” This representation is a special case of the more general representations with the propagator matrices, which are derived in a more formal way in Sec. V.

When we replace q(xA,ω) by a Green's matrix G(xA,xB,ω) with xB outside the region between x3,A and x3, we obtain

G(x,xB,ω)=DAW(x,xA,ω)G(xA,xB,ω)d2xA.
(31)

This is the simplest relation between the Green's matrix and propagator matrix. It is a special case of the more general representations with the Green's matrices and propagator matrices, which are derived in a more formal way in Sec. VI A. Alternatively, we may replace G(x,xB,ω) by the homogeneous Green's matrix Gh(x,xB,ω), where xB may be located anywhere because the homogeneous Green's matrix has no source at xB, hence,

Gh(x,xB,ω)=DAW(x,xA,ω)Gh(xA,xB,ω)d2xA.
(32)

This relation will be derived in a more formal way in Sec. VI B.

Next, we discuss the 2 × 2 acoustic propagator matrix as a special case of the N × N propagator matrix. For this situation, W is partitioned as

W(x,xA,ω)=(Wp,pWp,vWv,pWv,v)(x,xA,ω).
(33)

The first and second superscripts refer to the wave-field quantities at x and xA, respectively (where the superscript p, again, stands for the acoustic pressure and v is the vertical component of the particle velocity). The unit of a specific element of the propagator matrix is the ratio of the units of the wave-field quantities at x and xA. For example, [Wp,v(x,xA,ω)]=[p]/[v]=kgs1m2. For an isotropic medium, the elements can be expressed in terms of the upper-right element as

Wv,v(x,xA,ω)=1iωρ(x,ω)3Wp,v(x,xA,ω),
(34)
Wp,p(x,xA,ω)=1iωρ(xA,ω)3,AWp,v(x,xA,ω),
(35)
Wv,p(x,xA,ω)=1iωρ(x,ω)3Wp,p(x,xA,ω).
(36)

Equations (34) and (36) follow directly from Eqs. (11), (27), and (33). Equation (35) follows from Eq. (34) and a source-receiver reciprocity relation, which is derived in Sec. V A.

We illustrate the elements Wp,p and Wp,v, decomposed into plane waves, for the horizontally layered medium of Fig. 1(a). We choose, again, xA=(0,0,x3,0) and set s2=0. Usually, the propagator matrix is considered in the frequency domain, but to facilitate the comparison with the Green's function in Fig. 1(b), we consider the time-domain functions Wp,p(s1,x3,x3,0,τ) and Wp,v(s1,x3,x3,0,τ) [obtained via the transforms of Eqs. (20) and (21)], once more, for a single horizontal slowness s1=1/3500 s/m; see Figs. 2(a) and 3(a). At x3=x3,0 (the horizontal green lines in Figs. 2(a) and 3(a)), the boundary conditions for these functions are Wp,p(s1,x3,0,x3,0,τ)=δ(τ) and Wp,v(s1,x3,0,x3,0,τ)=0, respectively [this follows from applying the transforms of Eqs. (20) and (21) to Eq. (29), with xH,A=(0,0)]. Following these functions along the depth coordinate, starting at x3=x3,0, we observe an acausal upgoing wave and a causal downgoing wave until we reach the first interface. Here, both events split into upgoing and downgoing waves below the interface. The waves tunnel through the high-velocity layer, split again, and continue with higher amplitudes in the next layer. This illustrates that the evanescent waves may lead to unstable behaviour of the propagator matrix and should be handled with care. Note that Wp,p(s1,x3,x3,0,τ) and Wp,v(s1,x3,x3,0,τ) are symmetric and asymmetric in time, respectively. This is best seen in Figs. 2(b) and 3(b), which show the last traces of both elements. The same holds for elements Wv,v and Wv,p, which are not shown.

FIG. 2.

(Color online) The (a) propagator element Wp,p(s1,x3,x3,0,τ) (fixed s1), convolved with a wavelet, and (b) last trace of (a) are shown.

FIG. 2.

(Color online) The (a) propagator element Wp,p(s1,x3,x3,0,τ) (fixed s1), convolved with a wavelet, and (b) last trace of (a) are shown.

Close modal
FIG. 3.

(Color online) The (a) propagator element Wp,v(s1,x3,x3,0,τ) (fixed s1), convolved with a wavelet, and (b) last trace of (a) are shown.

FIG. 3.

(Color online) The (a) propagator element Wp,v(s1,x3,x3,0,τ) (fixed s1), convolved with a wavelet, and (b) last trace of (a) are shown.

Close modal

For an arbitrary inhomogeneous anisotropic medium, the time-domain version of boundary condition (29) reads

W(x,xA,t)|x3=x3,A=Iδ(xHxH,A)δ(t).
(37)

This boundary condition has a form similar to the focusing condition for the focusing functions appearing in the multidimensional Marchenko method.44 In Sec. VII A, we discuss the general relations between the propagator matrix and Marchenko-type focusing functions. Here, we present a short preview of these relations by considering the acoustic propagator matrix in the horizontally layered lossless isotropic medium of Fig. 1(a). We combine the elements Wp,p and Wp,v, decomposed into plane waves (Figs. 2 and 3), as60 

Fp(s1,x3,x3,0,τ)=Wp,p(s1,x3,x3,0,τ)s3,0ρ0Wp,v(s1,x3,x3,0,τ),
(38)

where the vertical slowness s3,0 is defined as

s3,0=1/c02s12,
(39)

where c0 and ρ0 are the propagation velocity and mass density of the upper half-space x3x3,0, respectively. Due to the symmetric form of Wp,p and the asymmetric form of Wp,v, half of the events double in amplitude and the other half of the events cancel. The result is shown in Fig. 4. For the interpretation of this focusing function, we start at the bottom of Fig. 4(a). The blue arrows indicate the upgoing waves, which are tuned such that after the interaction with the tunneling waves in the thin layer and downgoing wave just above the thin layer (indicated by a red arrow), they continue as a single upgoing wave, which finally focuses at the depth level x3,0 as a temporal delta function δ(τ) and continues as an upgoing wave into the homogeneous upper half-space.

FIG. 4.

(Color online) The (a) focusing function Fp(s1,x3,x3,0,τ) (fixed s1), convolved with a wavelet, and (b) last trace of (a) are shown.

FIG. 4.

(Color online) The (a) focusing function Fp(s1,x3,x3,0,τ) (fixed s1), convolved with a wavelet, and (b) last trace of (a) are shown.

Close modal

Note that, although we interpret the focusing function here in terms of upgoing and downgoing waves, it is derived from the propagator matrix (which does not rely on up/down decomposition), vertical slowness s3,0, and the mass density ρ0 of the upper half-space. Moreover, the focusing function is defined in the actual medium rather than in a truncated version of the actual medium as is usually the case for the Marchenko-type focusing functions.44 Hence, the only assumption is that a real-valued vertical slowness s3,0 exists in the upper half-space. Other than that, the focusing function Fp(s1,x3,x3,0,τ), as defined in Eq. (38), does not require a truncated medium, does not rely on up/down decomposition inside the medium, and does (at least in principle) not break down when the waves become evanescent inside the medium. These properties also hold for the more general version of the Marchenko-type focusing functions defined in Sec. VII.

As the basis for the derivation of the general representations with Green's matrices (Sec. IV), propagator matrices (Sec. V), or a combination thereof (Sec. VI), we introduce unified matrix-vector wave-field reciprocity theorems. In general, a wave-field reciprocity theorem interrelates two wave states (sources, wave fields, and medium parameters) in the same spatial domain.12 The reciprocity theorems have been formulated for acoustic,4 electromagnetic,61 elastodynamic,62,63 poroelastodynamic,64 piezoelectric,65,66 and seismoelectric waves.67 The matrix-vector equation, discussed in Sec. II, allows a unified formulation of the reciprocity theorems for these different wave phenomena,49,50 which extends the theory of the propagation invariants.36–39 Using the symmetry properties of operator matrix A, formulated in Eqs. (3) and (4), the matrix-vector reciprocity theorems can be derived49,50

D(dAtNqB+qAtNdB)d3x=D0DMqAtNqBn3d2x+DqAtN(AAAB)qBd3x
(40)

and

D(dAKqB+qAKdB)d3x=D0DMqAKqBn3d2x+DqAK(A¯AAB)qBd3x.
(41)

Here, D denotes a domain enclosed by two infinite horizontal boundaries D0 and DM at the depth levels x3,0 and x3,M with outward pointing normals n3=1 and n3=1, respectively; see Fig. 5. The subscripts A and B refer to two independent states. These theorems hold for lossless media and media with losses.51 Equation (40) is a convolution-type reciprocity theorem as products like qAtNqB in the frequency domain correspond to convolutions in the time domain (similar to those in Ref. 63). Equation (41) is a correlation-type reciprocity theorem because products such as qAKqB in the frequency domain correspond to correlations in the time domain (as in Ref. 18).

FIG. 5.

The configuration for the matrix-vector reciprocity theorems, Eqs. (40) and (41), and for the representations with Green's matrices is depicted.

FIG. 5.

The configuration for the matrix-vector reciprocity theorems, Eqs. (40) and (41), and for the representations with Green's matrices is depicted.

Close modal

A special case is obtained when the sources, wave fields, and medium parameters are identical in both states. We may then drop the subscripts A and B, and Eq. (41) simplifies to

D14(dKq+qKd)d3x=D0DM14qKqn3d2x+D14qK(A¯A)qd3x.
(42)

Because 14qKq=14(q1q2+q2q1)=j [see Eq. (7)], Eq. (42) formulates the unified power balance. The term on the left-hand side is the power generated by the sources in D. The first term on the right-hand side is the power flux through the boundary D0DM (i.e., the power leaving the domain D), and the second term on the right-hand side is the dissipated power in D.

A wave-field representation is obtained by replacing one of the states in a reciprocity theorem by a Green's state.6–9 In this section, we follow this approach to derive wave-field representations with Green's matrices from the matrix-vector reciprocity theorems discussed in Sec. III.

Before we derive wave-field representations, we first derive a symmetry property of the Green's matrix. To this end, we replace both of the wave-field vectors qA and qB in the reciprocity theorem (40) by the Green's matrices G(x,xA,ω) and G(x,xB,ω), respectively. Accordingly, we replace the source vectors dA and dB by Iδ(xxA) and Iδ(xxB), respectively, where xA and xB denote the source positions; see Fig. 5. Furthermore, we replace D by 3 such that the boundary integral in Eq. (40) vanishes (the Sommerfeld radiation condition). Both of the Green's matrices are defined in the same medium, hence, AA=AB. This implies that the second integral on the right-hand side of Eq. (40) also vanishes. From the remaining integral, we, thus, obtain the following symmetry property of the Green's matrix:

Gt(xB,xA,ω)N=NG(xA,xB,ω).
(43)

The Green's matrix on the left-hand side is the response to a source at xA, observed by a receiver at xB. Similarly, the Green's matrix on the right-hand side is the response to a source at xB, observed by a receiver at xA. Hence, Eq. (43) is a unified source-receiver reciprocity relation.

Using this relation and JN=NJ, for the homogeneous Green's matrix defined in Eq. (25), we find

Ght(xB,xA,ω)N=NGh(xA,xB,ω).
(44)

We derive a representation of the convolution type for the actual wave-field vector q(x,ω), emitted by the actual source distribution d(x,ω) in the actual medium; the operator matrix in the actual medium is defined as A(x,ω). We let state B in the reciprocity theorem (40) be this actual state, hence, we drop the subscript B from qB, dB, and AB. For state A, we choose the Green's state. Therefore, we replace qA(x,ω) in the reciprocity theorem (40) by G(x,xA,ω) and dA(x,ω) by Iδ(xxA). We keep the subscript A in AA(x,ω) to account for the fact that, in general, this operator matrix is defined in a medium that may be different from the actual medium. We, thus, obtain

χD(xA)Nq(xA,ω)=DGt(x,xA,ω)Nd(x,ω)d3x+D0DMGt(x,xA,ω)Nq(x,ω)n3d2x+DGt(x,xA,ω)N{AAA}q(x,ω)d3x,
(45)

where χD(x) is the characteristic function for the domain D, which is defined as

χD(x)={1for xinsideD,12for xonD0DM,0for xoutsideD.
(46)

Using the symmetry property of the Green's matrix, formulated by Eq. (43), we obtain

χD(xA)q(xA,ω)=DG(xA,x,ω)d(x,ω)d3xD0DMG(xA,x,ω)q(x,ω)n3d2xDG(xA,x,ω){AAA}q(x,ω)d3x.
(47)

This is the unified wave-field representation of the convolution type with the Green's matrix. The left-hand side is the wave-field vector q at a specific point xA, multiplied with the characteristic function. According to the right-hand side, this field consists of a contribution from the source distribution in D (the first integral), a contribution from the wave field on the boundary of D (the second integral), and a contribution caused by the contrasts between the operator matrices in the Green's and the actual state (the third integral).

Note that Eq. (15) follows as a special case of Eq. (47) if we choose the same medium parameters in both states and replace D by 3 [except that the roles of x and xA are interchanged because we used the symmetry property of Eq. (43) in the derivation of Eq. (47)].

Next, we consider another special case. Again, we choose the same medium parameters in both states, but this time we replace D by the entire half-space below D0, which we choose to be source-free in the actual state. Assuming xA lies in the lower half-space, we obtain

q(xA,ω)=D0G(xA,x,ω)q(x,ω)d2xfor x3,A>x3,0.
(48)

Using Eqs. (8), (16), and (18) for the isotropic acoustic situation, for the upper element of q(xA,ω), Eq. (48) yields

p(xA,ω)=D0(1iωρ(x,ω){3Gp,q(xA,x,ω)}p(x,ω)+Gp,q(xA,x,ω)v3(x,ω))d2x,
(49)

which is the well-known acoustic Kirchhoff-Helmholtz integral. Hence, the representation of Eq. (48) is the unified Kirchhoff-Helmholtz integral. It can be used for forward wave-field extrapolation of the wave-field vector q(x,ω) from D0 to any point xA below D0. Note that the Green's matrix G(xA,x,ω) depends on the medium parameters of the entire half-space below D0. In Sec. V B, we derive a relation similar to Eq. (48) but with the Green's matrix replaced by the propagator matrix, which depends only on the medium parameters between x3,0 and x3,A.

Finally, we replace q(x,ω) by G(x,xB,ω) in Eq. (48). Because we assumed that the lower half-space is source-free for q, we choose xB in the upper half-space. We, thus, obtain

G(xA,xB,ω)=D0G(xA,x,ω)G(x,xB,ω)d2xfor x3,A>x3,0>x3,B.
(50)

This expression shows that the Green's matrix can be composed from a cascade of Green's matrices, assuming a specific order of the depth levels at which the Green's sources and receivers are situated. In Sec. V B, we derive a similar relation for propagator matrices for an arbitrary order of depth levels.

Equation (50) can be seen as a generalization of the representation that underlies the Green's function retrieval by cross-convolution [where G(xA,xB,ω) is the unknown68] or by multidimensional deconvolution [where G(xA,x,ω) is the unknown69].

We derive a representation of the correlation type for the actual wave-field vector q(x,ω), emitted by the actual source distribution d(x,ω) in the actual medium. As in Sec. IV B, we let the state B be this actual state, and state A be the Green's state. Hence, making the same substitutions as in Sec. IV B, this time in the reciprocity theorem (41), using Eq. (43) and N1K=J, yields

χD(xA)q(xA,ω)=DJG*(xA,x,ω)Jd(x,ω)d3xD0DMJG*(xA,x,ω)Jq(x,ω)n3d2xDJG*(xA,x,ω)J{A¯AA}q(x,ω)d3x.
(51)

This is the unified wave-field representation of the correlation type with the Green's matrix.

As a special case, we derive a representation for the homogeneous Green's matrix Gh. To this end, for state B, we choose the Green's matrix in the actual medium, hence, we replace q(x,ω) by G(x,xB,ω) and d(x,ω) by Iδ(xxB). For state A, we replace the Green's matrix by that in the adjoint of the actual medium, therefore, we replace G(xA,x,ω) by G¯(xA,x,ω) and AA by A¯. With this choice, the contrast operator A¯AA=A¯¯A vanishes. Making these substitutions in Eq. (51), taking xA and xB both inside of D, and using Eq. (25), we obtain

Gh(xA,xB,ω)=D0DMJG¯*(xA,x,ω)JG(x,xB,ω)n3d2xforx3,M>x3,{A,B}>x3,0.
(52)

This is the unified homogeneous Green's matrix representation. It finds applications in optical, acoustic, and seismic holography,20,25 inverse source problems,16 inverse scattering methods,19 and Green's function retrieval by cross correlation.27–29,70,71 A disadvantage is that the integral is taken along the two boundaries D0 and DM, whereas in many practical situations, the measurements are only available on a single boundary. Using the propagator matrix, in Sec, VI B, we present a single-sided unified homogeneous Green's function representation.

Finally, using Eqs. (16)–(18) for the isotropic acoustic situation, for the upper-right element of Gh(xA,xB,ω), Eq. (52) yields

Ghp,q(xA,xB,ω)=1iωD0DM(1ρ¯*(x,ω){3G¯p,q(xA,x,ω)}*Gp,q(x,xB,ω)1ρ(x,ω){G¯p,q(xA,x,ω)}*3Gp,q(x,xB,ω))n3d2x,
(53)

where

Ghp,q(xA,xB,ω)=Gp,q(xA,xB,ω)+{G¯p,q(xA,xB,ω)}*,
(54)

according to Eqs. (16) and (26). For a lossless constant density medium, using Gp,q(x,xB,ω)=iωρG(x,xB,ω), Eq. (53) is the well-known scalar homogeneous Green's function representation,19,20 which is applied to the configuration of Fig. 5.

We follow an approach similar to that used in Sec. IV to derive wave-field representations. However, this time, we replace one of the states in the reciprocity theorems by a propagator state (instead of a Green's state). Hence, this leads to the wave-field representations with the propagator matrices.

We start with deriving the symmetry properties of the propagator matrix. To this end, we replace the wave-field vectors qA and qB in the reciprocity theorem (40) by the propagator matrices W(x,xA,ω) and W(x,xB,ω), respectively. We define horizontal boundaries DA and DB, containing the points xA and xB, respectively; see Fig. 5. We replace D0 and DM in Eq. (40) by these boundaries (and D by the region enclosed by these boundaries). Note that the boundary condition (29) implies W(x,xA,ω)=Iδ(xHxH,A) for x at DA and W(x,xB,ω)=Iδ(xHxH,B) for x at DB. The propagators obey the wave equation (27) without sources, hence, we set dA and dB to O. This implies that the integral on the left-hand side of Eq. (40) vanishes. Both propagator matrices are defined in the same medium, hence, AA=AB. This implies that the second integral on the right-hand side vanishes. Evaluating the remaining integral along the boundary DADB yields (irrespective of the arrangement of DA and DB)

Wt(xB,xA,ω)N=NW(xA,xB,ω).
(55)

This is the first unified symmetry relation for the propagator matrix.

A second symmetry relation can be derived from the reciprocity theorem (41). To this end, we replace qA and qB by W¯(x,xA,ω) and W(x,xB,ω), respectively. We replace D0DM in Eq. (41) by DADB (and D by the enclosed region). We set dA and dB to O, therefore, the integral on the left-hand side vanishes. The propagator matrices are defined in mutually adjoint media, hence, AA=A¯B. This implies that the second integral on the right-hand side of Eq. (41) vanishes. From the remaining integral, we, thus, obtain

W¯(xB,xA,ω)K=KW(xA,xB,ω).
(56)

From Eqs. (55) and (56), using KN1=J, we find that

W¯*(xB,xA,ω)J=JW(xB,xA,ω).
(57)

Note that in this last equation, xB and xA appear in the same order on the left- and right-hand sides.

We derive a representation of the convolution type for the actual wave-field vector q(x,ω), emitted by the actual source distribution d(x,ω) in the actual medium; the operator matrix in the actual medium is defined as A(x,ω). We let state B in reciprocity theorem (40) be this actual state, hence, we drop the subscript B from qB, dB, and AB. For state A, we choose the propagator state. Therefore, we replace qA(x,ω) in the reciprocity theorem (40) by W(x,xA,ω) and dA(x,ω) by O. We keep the subscript A in AA(x,ω) to account for the fact that, in general, this operator matrix is defined in a medium that may be different from the actual medium. We replace D0DM in the reciprocity theorem (40) by D0DA, where DA is the boundary containing xA. Here and in the following, DA is below D0, thus, x3,A>x3,0. The domain enclosed by these boundaries is called DA; see Fig. 6. Unlike in the classical, decomposition-based derivations of the Marchenko method,44,72 the domain DA does not define a truncated medium; in general, the medium below DA is inhomogeneous. Applying the mentioned substitutions in Eq. (40), using the boundary condition (29) and symmetry relation (55) (with xB replaced by x), yields

q(xA,ω)=DAW(xA,x,ω)d(x,ω)d3x+D0W(xA,x,ω)q(x,ω)d2xDAW(xA,x,ω){AAA}q(x,ω)d3x.
(58)

This is the unified wave-field representation of the convolution type with the propagator matrix. Note the analogy with the representation of the convolution type with the Green's matrix, Eq. (47). An important difference is that the boundary integral in Eq. (58) is single sided.

FIG. 6.

The configuration for the representations with the propagator matrices is shown.

FIG. 6.

The configuration for the representations with the propagator matrices is shown.

Close modal

We consider a special case by choosing the same medium parameters in DA in both states and choosing DA to be source free. We obtain

q(xA,ω)=D0W(xA,x,ω)q(x,ω)d2x.
(59)

This is the special case that was already presented in Eq. (30) [except that the roles of x and xA are interchanged because we used the symmetry property of Eq. (55) in the derivation of Eq. (58)].

Note the analogy of Eq. (59) with Eq. (48), which contains a Green's matrix instead of the propagator matrix. The propagator matrix in Eq. (59) depends only on the medium parameters between x3,0 and x3,A. Equation (59) is used in Sec. VII as the basis for deriving the Marchenko-type representations.

Next, we again consider Eq. (58) in which we replace q(x,ω) by W(x,xB,ω) and, hence, d(x,ω) by O and A by AA. This implies that the first and third integrals on the right-hand side vanish. Moreover, we replace D0 by DC at an arbitrarily chosen depth level x3,C; see Fig. 6. This yields

W(xA,xB,ω)=DCW(xA,x,ω)W(x,xB,ω)d2x.
(60)

This expression shows that the propagator matrix can be composed from a cascade of propagator matrices.32–36 It is similar to Eq. (50) with the cascade of the Green's matrices, but unlike in Eq. (50), where x3,A>x3,0>x3,B, the arrangement of x3,A,x3,B, and x3,C in Eq. (60) is arbitrary (because there are no sources in both states).

Finally, we replace xB in Eq. (60) by xA=(xH,A,x3,A). Applying the boundary condition (29) to the left-hand side, we obtain

Iδ(xH,AxH,A)=DCW(xA,x,ω)W(x,xA,ω)d2x.
(61)

This equation defines W(xA,x,ω) as the inverse of W(x,xA,ω).32–36 

We aim to derive a representation of the correlation type for the actual wave-field vector q(x,ω), emitted by the actual source distribution d(x,ω) in the actual medium. Similar to Sec. V B, we let state B be this actual state. For state A, we choose the adjoint propagator state. Hence, we replace qA(x,ω),AA(x,ω), and dA(x,ω) by W¯(x,xA,ω),A¯A(x,ω), and O, respectively. Furthermore, we again replace D0DM by D0DA. Applying these substitutions in the reciprocity theorem (41), using the boundary condition (29) and symmetry relation (56) (with xB replaced by x), yields, once more, Eq. (58). Therefore, due to the additional symmetry property of the propagator matrix [Eq. (56)], the correlation-type reciprocity theorem does not lead to a new representation.

We follow an approach, which is similar to that in Secs. IV and V, to derive wave-field representations. However, this time, we replace one of the states in the reciprocity theorem of the convolution type by a Green's state and the other state by a propagator state. This leads to wave-field representations with Green's matrices and propagator matrices.

In the reciprocity theorem (40), we choose the propagator state for state A and the Green's state for state B. Hence, in state A, we replace qA(x,ω) and dA(x,ω) by W(x,xA,ω) and O, respectively. In state B, we replace qB(x,ω) and dB(x,ω) by G(x,xB,ω) and Iδ(xxB), respectively. The medium parameters in states A and B may be different, thus, we keep the subscripts in AA(x,ω) and AB(x,ω). Last but not least, we replace D0DM by D0DA and D by the enclosed region DA; see Fig. 6. With these substitutions, using boundary condition (29) and symmetry relation (55), the reciprocity theorem (40) yields

G(xA,xB,ω)χDA(xB)W(xA,xB,ω)=D0W(xA,x,ω)G(x,xB,ω)d2xDAW(xA,x,ω){AAAB}G(x,xB,ω)d3x,
(62)

where χDA(x) is the characteristic function for domain DA, defined similarly as χD(x) in Eq. (46). Equation (62) is a representation for GW (when xB is inside DA) in terms of integrals containing G and W. When G and W are defined in the same medium, this representation simplifies to

G(xA,xB,ω)χDA(xB)W(xA,xB,ω)=D0W(xA,x,ω)G(x,xB,ω)d2x.
(63)

Finally, when xB (the position of the source of the Green's matrix) lies outside of DA (i.e., above D0 or below DA), then the latter expression simplifies to

G(xA,xB,ω)=D0W(xA,x,ω)G(x,xB,ω)d2xforx3,B<x3,0orx3,B>x3,A.
(64)

This is a single-sided representation of the Green's matrix, which uses the propagator matrix. This special case was already presented in Eq. (31) (except that the roles of x and xA are interchanged as we used the symmetry property of Eq. (55) in the derivation of Eq. (64)].

Note the analogy of Eq. (64) with Eq. (50), which contains a Green's matrix instead of the propagator matrix. The propagator matrix in Eq. (64) depends only on the medium parameters between x3,0 and x3,A. Moreover, unlike the representation of Eq. (50), which only holds for x3,B<x3,0, Eq. (64) also holds for x3,B>x3,A.

We follow the same procedure as in Sec. VI A, except that in state B, we choose the homogeneous Green's matrix; hence, we replace qB(x,ω) and dB(x,ω) in the reciprocity theorem (40) by Gh(x,xB,ω) and O, respectively. We, thus, obtain

Gh(xA,xB,ω)=D0W(xA,x,ω)Gh(x,xB,ω)d2xDAW(xA,x,ω){AAAB}Gh(x,xB,ω)d3x.
(65)

When Gh and W are defined in the same medium, this representation simplifies to

Gh(xA,xB,ω)=D0W(xA,x,ω)Gh(x,xB,ω)d2x.
(66)

This is a single-sided representation of the homogeneous Green's matrix, which uses the propagator matrix. This special case was already presented in Eq. (32) [except that the roles of x and xA are interchanged as we used the symmetry property of Eq. (55) in the derivation of Eq. (66)]. Note that unlike the Green's matrix representation of Eq. (64), the homogeneous Green's matrix representation of Eq. (66) has no restrictions for the position of xB (as there is no source at xB).

Equation (66) is the single-sided counterpart of the unified homogeneous Green's matrix representation of Eq. (52). Whereas in Eq. (52) the integral is taken along two boundaries D0 and DM, in Eq. (66) the integral is taken only along D0. This is an important advantage for the practical situations in which the measurements are often available only on a single boundary. In Sec. VI C, we indicate how Eq. (66) can be used in a process called source and receiver redatuming.

Finally, using Eqs. (16), (17), (26), (33), and (35) for the isotropic acoustic situation, assuming real-valued ρ(x), Eq. (66) yields for the upper-right element of Gh(xA,xB,ω)

Ghp,q(xA,xB,ω)=1iωD01ρ(x)({3Wp,v(xA,x,ω)}Ghp,q(x,xB,ω)Wp,v(xA,x,ω)3Ghp,q(x,xB,ω))d2x,
(67)

with Ghp,q(xA,xB,ω) as defined in Eq. (54).

Redatuming is the process of moving sources and/or receivers from the acquisition boundary to positions inside of the medium.22,57,73,74 Traditionally, this is done with the Kirchhoff-Helmholtz integrals with the Green's functions defined in a smooth background medium, thus, ignoring multiple scattering. Here, we derive a unified redatuming method, which accounts for multiple scattering, using the single-sided homogeneous Green's matrix representation of Sec. VI B as the starting point.

Let xS and xR denote the source and receiver coordinates, respectively, at the acquisition boundary D0. Then the Green's matrix G(xR,xS,ω) represents the medium's response, measured with the sources and receivers at the surface; see Fig. 7(a). Assuming that the medium is lossless, the homogeneous Green's matrix Gh(xR,xS,ω) is obtained from this Green's matrix via Eq. (25).

FIG. 7.

The illustration of source and receiver redatuming. All of the responses are represented by simple rays, but in reality these are multi-component wave fields, including primaries, multiples, converted, refracted, and evanescent waves. (a) The response G(xR,xS,ω) at the surface is shown. The homogeneous Green's matrix Gh(xR,xS,ω) is obtained from this response with Eq. (25). (b) The source redatuming is depicted. Using Eq. (69), the homogeneous Green's matrix Gh(xR,xB,ω) is obtained for a virtual source at xB. (c) The receiver redatuming is shown. Using Eq. (70), the homogeneous Green's matrix Gh(xA,xB,ω) is obtained for a virtual receiver at xA.

FIG. 7.

The illustration of source and receiver redatuming. All of the responses are represented by simple rays, but in reality these are multi-component wave fields, including primaries, multiples, converted, refracted, and evanescent waves. (a) The response G(xR,xS,ω) at the surface is shown. The homogeneous Green's matrix Gh(xR,xS,ω) is obtained from this response with Eq. (25). (b) The source redatuming is depicted. Using Eq. (69), the homogeneous Green's matrix Gh(xR,xB,ω) is obtained for a virtual source at xB. (c) The receiver redatuming is shown. Using Eq. (70), the homogeneous Green's matrix Gh(xA,xB,ω) is obtained for a virtual receiver at xA.

Close modal

First, we discuss a method for source redatuming. We rename some variables in Eq. (66) (xBxR,xAxB,xxS) and transpose all matrices, which gives

Ght(xB,xR,ω)=D0Ght(xS,xR,ω)Wt(xB,xS,ω)d2xS.
(68)

Using Eqs. (44) and (55), we obtain

Gh(xR,xB,ω)=D0Gh(xR,xS,ω)W(xS,xB,ω)d2xS.
(69)

The latter expression describes the redatuming of the sources from xS at the acquisition boundary to a virtual-source position xB inside of the medium; see Fig. 7(b).

Next, we discuss receiver redatuming. We replace x by xR in Eq. (66), which gives

Gh(xA,xB,ω)=D0W(xA,xR,ω)Gh(xR,xB,ω)d2xR.
(70)

This expression describes the redatuming of the receivers from xR at the acquisition boundary to a virtual-receiver position xA inside of the medium; see Fig. 7(c). The Green's matrix under the integral is the output of the source redatuming method, which is described by Eq. (69).

Combining Eqs. (69) and (70) gives the following expression for the combined source and receiver redatuming:

Gh(xA,xB,ω)=D0D0W(xA,xR,ω)Gh(xR,xS,ω)W(xS,xB,ω)d2xSd2xR.
(71)

This expression redatums the actual sources and receivers from the acquisition boundary D0 to the virtual sources and receivers inside of the medium. Because it takes multiple scattering into account, it generalizes the classical source and receiver redatuming.22,57,73,74 Equation (71) resembles expressions for source-receiver interferometry75,76 but with the integrals along a closed boundary replaced by the integrals along the open acquisition boundary D0 and two of the Green's functions replaced by the propagator matrices. When the medium is known, these matrices can be numerically modelled. Alternatively, the propagator matrices can be expressed in terms of focusing functions which, at least for the acoustic situation, can be retrieved via the Marchenko method from the reflection response at the acquisition boundary. In Sec. VII, we introduce the relations between the propagator matrix and focusing functions (Sec. VII A) and derive Marchenko-type representations, which relate the focusing functions to the reflection response at the acquisition boundary (Sec. VII B).

In Sec. II E, we indicated that there is a relation between the propagator matrix and the Marchenko-type focusing functions. Here, we derive this relation for the propagator matrix of the unified matrix-vector wave equation, assuming that the medium is lossless. Our starting point is Eq. (59), which is repeated here for convenience,

q(xA,ω)=D0W(xA,x,ω)q(x,ω)d2x,
(72)

with the boundary condition

W(xA,x,ω)|x3,A=x3,0=Iδ(xH,AxH)
(73)

for x at D0. From here on, we assume that the half-space above D0 (including D0) is homogeneous and isotropic, the medium below D0 is arbitrary inhomogeneous, anisotropic, and source-free, and xA is chosen at or below D0 (hence, x3,Ax3,0).

Without loss of generality, we can decompose q(x,ω) in the upper half-space into downgoing and upgoing plane waves. To this end, we defined the spatial Fourier transform ũ(s,x3,ω) of a space- and frequency-dependent function u(x,ω) in Eq. (20). For a function of two space variables, u(xA,x,ω), we define the spatial Fourier transform along the horizontal components of the second space variable as

ũ(xA,s,x3,ω)=2u(xA,xH,x3,ω)exp{iωsxH}d2xH.
(74)

Using these definitions and Parseval's theorem, we rewrite Eq. (72) as

q(xA,ω)=ω24π22W̃(xA,s,x3,0,ω)q̃(s,x3,0,ω)d2s,
(75)

with the boundary condition

W̃(xA,s,x3,0,ω)|x3,A=x3,0=Iexp{iωsxH,A}.
(76)

In the upper half-space, we relate the wave-field vector q̃ to the downgoing and upgoing wave-field vectors p̃+ and p̃, respectively, via

q̃=(q̃1q̃2)=(L̃1+L̃1L̃2+L̃2)(p̃+p̃)forx3x3,0
(77)

(Refs. 55–57 and 77). Next, we renormalize the downgoing and upgoing wave-field vectors. To this end, we define the downgoing and upgoing parts of q̃1 as

q̃1±=L̃1±p̃±
(78)

and rewrite Eq. (77) as

q̃=(q̃1q̃2)=(IID̃1+D̃1)D̃1(q̃1+q̃1)b̃1forx3x3,0,
(79)

with

D̃1±=L̃2±(L̃1±)1.
(80)

In  Appendix A, we give explicit expressions for D̃1± for the acoustic, electromagnetic, and elastodynamic situation. The substitution of Eq. (79) into Eq. (75) gives

q(xA,ω)=ω24π22Ỹ1(xA,s,x3,0,ω)b̃1(s,x3,0,ω)d2s,
(81)

with

Ỹ1(xA,s,x3,0,ω)=W̃(xA,s,x3,0,ω)D̃1(s).
(82)

Let F̃1(xA,s,x3,0,ω) be the upper-right N/2×N/2 matrix of the block matrix Ỹ1(xA,s,x3,0,ω). According to Eqs. (28), (82), and the definition of D̃1(s) in Eq. (79), it is defined as

F̃1(xA,s,x3,0,ω)=(W̃11+W̃12D̃1)(xA,s,x3,0,ω).
(83)

This is a generalization of the definition of the acoustic focusing function for a horizontally layered medium, which is defined in Eq. (38). From Eqs. (28) and (76), it follows that F̃1(xA,s,x3,0,ω) obeys the boundary condition

F̃1(xA,s,x3,0,ω)|x3,A=x3,0=Iexp{iωsxH,A},
(84)

or, applying an inverse spatial and temporal Fourier transform,

F1(xA,x,t)|x3,A=x3,0=Iδ(xH,AxH)δ(t),
(85)

for x at D0. Hence, F1(xA,x,t) is indeed a focusing function (where xA is a variable and x is the focal point at D0).

To analyse the upper-left N/2×N/2 matrix of the block matrix Ỹ1(xA,s,x3,0,ω), we first establish some symmetry properties. The symmetry of W̃ follows from the Fourier transform of Eq. (57) for the lossless situation, hence,

W̃(xA,s,x3,0,ω)=JW̃*(xA,s,x3,0,ω)J.
(86)

For the sub-matrices of W̃, this implies that

W̃αβ(xA,s,x3,0,ω)=JααW̃αβ*(xA,s,x3,0,ω)Jββ,
(87)

(no summation convention) where J11=J22=I. The symmetry of D̃1± follows from its explicit definitions in  Appendix A. In all of the cases, the following symmetry holds:

D̃1+(s)=D̃1(s).
(88)

When we ignore the evanescent waves at and above D0, then D̃1±(s) is real valued; see  Appendix A. Hence, given that F̃1(xA,s,x3,0,ω) is the upper-right N/2×N/2 matrix of Ỹ1(xA,s,x3,0,ω), we derive from Eqs. (28), (82), (87), and (88) and the structure of D̃1 indicated in Eq. (79) that the upper-left N/2×N/2 matrix is equal to F̃1*(xA,s,x3,0,ω). Using this in Eq. (81), for the upper N/2×1 vector of q(xA,ω), we obtain

q1(xA,ω)=ω24π22F̃1*(xA,s,x3,0,ω)q̃1+(s,x3,0,ω)d2s+ω24π22F̃1(xA,s,x3,0,ω)q̃1(s,x3,0,ω)d2s.
(89)

By applying Parseval's theorem again, we obtain

q1(xA,ω)=D0F1*(xA,x,ω)q1+(x,ω)d2x+D0F1(xA,x,ω)q1(x,ω)d2x.
(90)

This relation has previously been derived via another route for the situations of acoustic waves (q1=p) and elastodynamic waves (q1=v),78 but without explicitly defining the focusing function F1(xA,x,ω). Here, we have an explicit expression for the Fourier transform of this focusing function [Eq. (83)]. In Sec. VII B, we use Eq. (90) as the basis for deriving the Marchenko-type Green's matrix representations.

Next, we derive a representation similar to Eq. (90) for q2(xA,ω). To this end, we define the downgoing and upgoing parts of q̃2 as

q̃2±=L̃2±p̃±
(91)

and rewrite Eq. (77), analogous to Eq. (79), as

q̃=(q̃1q̃2)=(D̃2+D̃2II)D̃2(q̃2+q̃2)b̃2forx3x3,0,
(92)

with

D̃2±=L̃1±(L̃2±)1=(D̃1±)1.
(93)

The substitution of Eq. (92) into Eq. (75) gives

q(xA,ω)=ω24π22Ỹ2(xA,s,x3,0,ω)b̃2(s,x3,0,ω)d2s,
(94)

with

Ỹ2(xA,s,x3,0,ω)=W̃(xA,s,x3,0,ω)D̃2(s).
(95)

Let F̃2(xA,s,x3,0,ω) be the lower-right N/2×N/2 matrix of the block matrix Ỹ2(xA,s,x3,0,ω). According to Eqs. (28), (95), and the definition of D̃2(s) in Eq. (92), it is defined as

F̃2(xA,s,x3,0,ω)=(W̃21D̃2+W̃22)(xA,s,x3,0,ω).
(96)

F̃2 is a focusing function with focusing properties similar to those of F̃1, expressed in Eqs. (84) and (85).

To analyse the lower-left N/2×N/2 matrix of the block matrix Ỹ2(xA,s,x3,0,ω), from Eqs. (88) and (93), we first derive

D̃2+(s)=D̃2(s).
(97)

When we ignore the evanescent waves at and above D0, then D̃2±(s) is real valued. Hence, given that F̃2(xA,s,x3,0,ω) is the lower-right N/2×N/2 matrix of Ỹ2(xA,s,x3,0,ω), we derive from Eqs. (28), (87), (95), and (97) and the structure of D̃2 indicated in Eq. (92) that the lower-left N/2×N/2 matrix is equal to F̃2*(xA,s,x3,0,ω). Using this in Eq. (94), we obtain for the lower N/2×1 vector of q(xA,ω)

q2(xA,ω)=ω24π22F̃2*(xA,s,x3,0,ω)q̃2+(s,x3,0,ω)d2s+ω24π22F̃2(xA,s,x3,0,ω)q̃2(s,x3,0,ω)d2s.
(98)

Applying Parseval's theorem again, we obtain

q2(xA,ω)=D0F2*(xA,x,ω)q2+(x,ω)d2x+D0F2(xA,x,ω)q2(x,ω)d2x.
(99)

From Eqs. (78), (80), and (91), we obtain q̃2±=D̃1±q̃1±. Substituting this into Eq. (98), we obtain a representation for q2(xA,ω) in terms of q̃1+ and q̃1. Combining this with Eq. (89) into a single equation, we obtain Eq. (81) with [using Eq. (88)]

Ỹ1(xA,s,x3,0,ω)=(F̃1*(xA,s,x3,0,ω)F̃1(xA,s,x3,0,ω)F̃2*(xA,s,x3,0,ω)D̃1(s)F̃2(xA,s,x3,0,ω)D̃1(s)).
(100)

Once Ỹ1 is known, the propagator matrix W̃ follows from inverting Eq. (82), according to

W̃(xA,s,x3,0,ω)=Ỹ1(xA,s,x3,0,ω){D̃1(s)}1.
(101)

Equations (100) and (101), with D̃1(s) defined in Eq. (79), express the propagator matrix explicitly in terms of the focusing functions.

For the acoustic situation, using Eq. (A16), Eqs. (83) and (96) become

F̃p(xA,s,x3,0,ω)=(W̃p,ps3,0ρ0W̃p,v)(xA,s,x3,0,ω),
(102)
F̃v(xA,s,x3,0,ω)=(ρ0s3,0W̃v,p+W̃v,v)(xA,s,x3,0,ω),
(103)

with s3 defined as in Eq. (A15). For the acoustic situation, Eq. (100) yields

Ỹ1(xA,s,x3,0,ω)=({F̃p(xA,s,x3,0,ω)}*F̃p(xA,s,x3,0,ω)s3,0ρ0{F̃v(xA,s,x3,0,ω)}*s3,0ρ0F̃v(xA,s,x3,0,ω)).
(104)

Using this in Eq. (101), together with

{D̃1(s)}1=(12ρ02s3,012ρ02s3,0),
(105)

we obtain an explicit expression for W̃(xA,s,x3,0,ω) in terms of the acoustic focusing functions. Transforming this to the space-frequency domain, replacing s3,0 by (1/ω)H1 (where H1 is the square-root of the Helmholtz operator ω2/c02+αα in the homogeneous upper half-space54,56,57) yields60 

W(xA,x,ω)=({Fp(xA,x,ω)}iωρ0H11(x,ω){Fp(xA,x,ω)}1iωρ0H1(x,ω){Fv(xA,x,ω)}{Fv(xA,x,ω)})
(106)

for x at D0, where denotes the imaginary part.

We use Eqs. (90) and (99) as the starting point for deriving two Marchenko-type Green's matrix representations. We follow a similar procedure as in Ref. 78 (Appendix B), generalized for the different wave phenomena considered in this paper. We replace the N/2×1 vector q1(x,ω) by a modified version Γ12(x,xS,ω) of the N/2×N/2 Green's matrix G12(x,xS,ω). Here, G12(x,xS,ω) stands for the q1-type field observed at x in response to a unit d2-type source at xS. We choose xS=(xH,S,x3,S) in the upper half-space at a vanishing distance ϵ above D0, hence, x3,S=x3,0ϵ. Our aim is to modify this Green's matrix such that for x at D0 (i.e., just below the source), its downgoing part simplifies to

Γ12+(x,xS,ω)|x3=x3,0=Iδ(xHxH,S).
(107)

First, we derive the properties of the downgoing part of G12(x,xS,ω) for x at D0. To this end, consider the inverse of Eq. (79), which reads

(q̃1+q̃1)=((Δ̃1)1D̃1(Δ̃1)1(Δ̃1)1D̃1+(Δ̃1)1)(D̃1)1(q̃1q̃2),
(108)

for x3x3,0 with

Δ̃1=D̃1+D̃1.
(109)

The upper-right N/2×N/2 matrix (Δ̃1)1 transforms q̃2 into a downgoing field vector q̃1+. In a similar way, this matrix transforms a unit d2-type source in a homogeneous half-space into the downgoing part of the Green's matrix G̃12 just below this source, according to

limx3x3,SG̃12+(s,x3,0,x3,S,ω)={Δ̃1(s)}1.
(110)

The explicit expressions for (Δ̃1)1 are given in  Appendix A. In Eq. (110), the source is located at (0,x3,S). Next, we consider G12(x,xS,ω) for a laterally shifted source position (xH,S,x3,S). Applying a spatial Fourier transform along the horizontal source coordinate xH,S, using Eq. (74) with xH replaced by xH,S, yields G̃12(x,s,x3,S,ω). For the downgoing part just below the source, we obtain a phase-shifted version of the Green's matrix of Eq. (110), according to

limx3x3,SG̃12+(x,s,x3,S,ω)={Δ̃1(s)}1exp{iωsxH}.
(111)

This suggests that the modified Green's matrix be defined as

Γ̃12(x,s,x3,S,ω)=G̃12(x,s,x3,S,ω)Δ̃1(s),
(112)

such that

limx3x3,SΓ̃12+(x,s,x3,S,ω)=Iexp{iωsxH}.
(113)

Transforming this back to the space domain yields, indeed, Eq. (107). Next, we define the reflection response R12(x,xS,ω) of the medium below D0 as the upgoing part of the modified Green's matrix Γ12(x,xS,ω) for x at D0, hence,

R12(x,xS,ω)=Γ12(x,xS,ω).
(114)

Substituting q1(xA,ω)=Γ12(xA,xS,ω) and q1±(x,ω)=Γ12±(x,xS,ω) into Eq. (90), using Eqs. (107) and (114), gives

Γ12(xA,xS,ω)=D0F1(xA,x,ω)R12(x,xS,ω)d2x+F1*(xA,xS,ω)
(115)

for x3,Ax3,0. This is the first Marchenko-type representation. In a similar way, we derive a second Marchenko-type representation for a modified version Γ22(x,xS,ω) of the Green's matrix G22(x,xS,ω). We derive the properties of the downgoing part of G22(x,xS,ω) for x at D0. Consider the inverse of Eq. (92), which reads

(q̃2+q̃2)=((Δ̃2)1(Δ̃2)1D̃2(Δ̃2)1(Δ̃2)1D̃2+)(D̃2)1(q̃1q̃2),
(116)

for x3x3,0 with

Δ̃2=D̃2+D̃2.
(117)

The upper-right N/2×N/2 matrix (Δ̃2)1D̃2 transforms q̃2 into a downgoing field vector q̃2+. In a similar way, this matrix transforms a unit d2-type source in a homogeneous half-space into the downgoing part of the Green's matrix G̃22 just below this source, according to

limx3x3,SG̃22+(s,x3,0,x3,S,ω)={Δ̃2(s)}1D̃2(s).
(118)

The explicit expressions for (Δ̃2)1D̃2 are given in  Appendix A. Steps similar to those below Eq. (110) lead to

limx3x3,SG̃22+(x,s,x3,S,ω)={Δ̃2(s)}1D̃2(s)exp{iωsxH}.
(119)

This suggests that the modified Green's matrix be defined as

Γ̃22(x,s,x3,S,ω)=G̃22(x,s,x3,S,ω){D̃2(s)}1Δ̃2(s),
(120)

such that

Γ22+(x,xS,ω)|x3=x3,0=Iδ(xHxH,S).
(121)

We define the reflection response R22(x,xS,ω) as

R22(x,xS,ω)=Γ22(x,xS,ω)
(122)

for x at D0. Substituting q2(xA,ω)=Γ22(xA,xS,ω) and q2±(x,ω)=Γ22±(x,xS,ω) into Eq. (99), using Eqs. (121) and (122), gives

Γ22(xA,xS,ω)=D0F2(xA,x,ω)R22(x,xS,ω)d2x+F2*(xA,xS,ω)
(123)

for x3,Ax3,0. This is the second Marchenko-type representation.

For the acoustic case, using Eqs. (A16) and (A17), Eqs. (112) and (120) become

Γ̃12(x,s,x3,S,ω)=2s3,0ρ0G̃p,q(x,s,x3,S,ω)=2iωρ03,SG̃p,q(x,s,x3,S,ω),
(124)
Γ̃22(x,s,x3,S,ω)=2G̃v,q(x,s,x3,S,ω),
(125)

or in the space-frequency domain [using Eq. (18)],

Γ12(x,xS,ω)=2Gp,f(x,xS,ω),
(126)
Γ22(x,xS,ω)=2Gv,q(x,xS,ω).
(127)

For this situation, R12(x,xS,ω) and R22(x,xS,ω) in the representations of Eqs. (115) and (123) are the upgoing parts of 2Gp,f(x,xS,ω) and 2Gv,q(x,xS,ω), respectively, for x at D0. Note that according to Eq. (43), Gv,q(x,xS,ω)=Gp,f(xS,x,ω).

In their general form, the representations of Eqs. (115) and (123) are generalizations of the previously derived representations for the 3D Marchenko method for acoustic72,79,80 and elastodynamic wave fields81,82 (but note that the subscripts 1 and 2 of the focusing functions have a different meaning than in those papers; here, they refer to the wave-field components q1 and q2). In most of the previous work on the 3D Marchenko method, one of the underlying assumptions is that the wave field can be decomposed into downgoing and upgoing waves in the interior of the medium. Only recently, several authors proposed to avoid the decomposition inside the medium.78,83,84 The representations discussed in this section expand on this. In these representations, decomposition into downgoing and upgoing waves and negligence of the evanescent waves occurs only in the upper half-space. However, inside the medium, no wave-field decomposition takes place. Moreover, the evanescent waves inside the medium (for example, in high-velocity layers) are accounted for by the representations of Eqs. (115) and (123). These representations, transformed to the time domain, form the basis for the development of Marchenko schemes, aiming at resolving the focusing functions F1 and F2 from the reflection responses R12 and R22, respectively, at the acquisition boundary D0. Such schemes have been successfully developed for precritical acoustic data.45–47,85–88 For more complex situations, research on retrieving the focusing functions from the reflection responses is ongoing (for example, Ref. 89). Once the focusing functions are found, they can be used to define the propagator matrix via Eqs. (100) and (101), which can, subsequently, be used in the representations of Secs. V and VI.

We have discussed different types of wave-field representation in a systematic way. Classical wave-field representations contain Green's functions. By starting with a unified matrix-vector wave equation, we have formulated wave-field representations with Green's matrices, analogous to the classical representations. For example, the classical Kirchhoff-Helmholtz integral follows as a special case of the unified representation with the Green's matrix. Another special case is the classical homogeneous Green's function representation.

Using the same matrix-vector formalism, we formulated wave-field representations with propagator matrices. Unlike a Green's matrix, a propagator matrix depends only on the medium parameters between the two depth levels for which this matrix is defined. The representations with the propagator matrices have a similar form as those with the Green's matrices. An important difference is that the boundary integrals in the representations with the propagator matrices are single sided.

We also formulated representations containing a mix of Green's matrices and propagator matrices. A special case is the single-sided homogeneous Green's function representation as the counterpart of the classical closed-boundary homogeneous Green's function representation.

We have shown that the propagator matrix is related to Marchenko-type focusing functions. We have used this relation to reformulate the representations with the propagator matrix into representations with focusing functions. For the acoustic situation, these focusing functions can be retrieved from the single-sided reflection response of the medium by applying the Marchenko method. For more complex situations, research on retrieving these focusing functions from the reflection response is ongoing. Once the focusing functions are known, they can be used to construct the propagator matrix. Subsequently, the propagator matrix can be used in the representations to obtain the wave field inside of the medium which, in turn, can be used, for example, for imaging or monitoring. Unlike earlier imaging methods, which use the propagator matrix, this is a data-driven approach (because the propagator matrix is retrieved from the reflection response) and, hence, it does not require a detailed model of the medium.

The author thanks Roel Snieder, Sjoerd de Ridder, and Evert Slob for fruitful discussions and two anonymous reviewers for their positive evaluations of the original manuscript and useful suggestions for further improvements. This research is funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant No. 742703).

We derive the explicit expressions for D̃1±,(Δ̃1)1, and (Δ̃2)1D̃2 in the homogeneous isotropic upper half-space x3x3,0 for the acoustic, electromagnetic, and elastodynamic situations.

1. Unified expressions in the horizontal slowness domain

For the lossless homogeneous isotropic upper half-space x3x3,0, we transform the unified matrix-vector wave equation (1) to the horizontal slowness domain using the spatial Fourier transform of Eq. (20). This yields

3q̃Ãq̃=d̃,
(A1)

where

q̃=(q̃1q̃2),d̃=(d̃1d̃2),Ã=(Ã11Ã12Ã21Ã22).
(A2)

The symmetry property, defined by Eq. (3), transforms to

{Ã(s,ω)}tN=NÃ(s,ω).
(A3)

We decompose matrix à as55,57,77

Ã=L̃Λ̃L̃1,
(A4)

where

L̃=(L̃1+L̃1L̃2+L̃2),Λ̃=(Λ̃+OOΛ̃).
(A5)

With a proper scaling of the columns of matrix L̃, these matrices obey the following symmetry properties:

{L̃(s)}tN=N{L̃(s)}1,
(A6)
{Λ̃(s,ω)}tN=NΛ̃(s,ω).
(A7)

The substitution of Eq. (A4) into Eq. (A1) and premultiplying all terms with L̃1 yields

3p̃Λ̃p̃=L̃1d̃,
(A8)

where

p̃=L̃1q̃=(p̃+p̃),
(A9)

where p̃+ and p̃ are wave-field vectors containing the downgoing (+) and upgoing (-) waves, respectively. Following Eqs. (80), (93), (109), and (117), we define

D̃1±(s)=L̃2±(L̃1±)1,{Δ̃1(s)}1=(D̃1+D̃1)1,
(A10)
D̃2±(s)=L̃1±(L̃2±)1,{Δ̃2(s)}1=(D̃2+D̃2)1.
(A11)
2. Acoustic wave equation

For the acoustic wave equation, the 1 × 1 sub-matrices of à are

Ã12=iωρ,Ã21=iω(κsr2/ρ),Ã11=Ã22=0,
(A12)

where the radial slowness sr is defined as

sr2=sαsα=s12+s22.
(A13)

Here, κ and ρ are the compressibility and mass density of the homogeneous upper half-space, respectively. For convenience, here and in the remainder of this appendix, we do not use the subscripts 0 to indicate the parameters of the upper half-space. The 1 × 1 sub-matrices of L̃ and Λ̃ are

L̃1±=(ρ2s3)1/2,L̃2±=±(s32ρ)1/2,Λ̃±=±iωs3,
(A14)

where the vertical slowness s3 is defined as

s3={1/c2sr2forsr21/c2,isr21/c2forsr2>1/c2,
(A15)

where the propagation velocity c is defined as c=1/κρ. The two expressions in Eq. (A15) represent the propagating and evanescent waves. On substitution of Eq. (A14) into Eqs. (A10) and (A11), we obtain

D̃1±(s)=±s3ρ,{Δ̃1(s)}1=ρ2s3,
(A16)
D̃2±(s)=±ρs3,{Δ̃2(s)}1D̃2(s)=12.
(A17)
3. Electromagnetic wave equation

For the electromagnetic wave equation, the 2 × 2 sub-matrices of à are

Ã12=iω(μs12εs1s2εs1s2εμs22ε),Ã21=iω(εs22μs1s2μs1s2μεs12μ),
(A18)
Ã11=Ã22=O.
(A19)

Here, ε and μ are the permittivity and permeability of the upper half-space, respectively. The 2 × 2 sub-matrices of L̃ and Λ̃ are

L̃1±=i2sr(s1(s3ε)1/2s2(μs3)1/2s2(s3ε)1/2s1(μs3)1/2),
(A20)
L̃2±=±i2sr(s1(εs3)1/2s2(s3μ)1/2s2(εs3)1/2s1(s3μ)1/2),
(A21)
Λ̃±=±iω(s300s3),
(A22)

where sr and s3 are defined in Eqs. (A13) and (A15), respectively, but this time with the propagation velocity c defined as c=1/εμ. On substitution of Eqs. (A20) and (A21) into Eqs. (A10) and (A11), we obtain

D̃1±(s)=±1μs3(c2s22s1s2s1s2c2s12),
(A23)
{Δ̃1(s)}1=12εs3(c2s12s1s2s1s2c2s22),
(A24)
D̃2±(s)=±1εs3(c2s12s1s2s1s2c2s22),
(A25)
(Δ̃2)1D̃2=(120012).
(A26)
4. Elastodynamic wave equation

For the elastodynamic wave equation, the 3 × 3 sub-matrices of à are

Ã11=iω(00s100s2λλ+2μs1λλ+2μs20),
(A27)
Ã12=iω(1μ0001μ0001λ+2μ),
(A28)
Ã21=iω(ρν1s12μs22(ν2+μ)s1s20(ν2+μ)s1s2ρμs12ν1s22000ρ),
(A29)
Ã22=(Ã11)t,
(A30)

where

ν1=4μ(λ+μλ+2μ),ν2=2μ(λλ+2μ),
(A31)

where λ and μ are the Lamé parameters and ρ is the mass density of the upper half-space. The 3 × 3 sub-matrices of L̃ and Λ̃ are

L̃1±=1(2ρ)1/2(s1(s3P)1/2s1(s3S)1/2srs2cSsr(s3S)1/2s2(s3P)1/2s2(s3S)1/2srs1cSsr(s3S)1/2±(s3P)1/2±sr(s3S)1/20),
(A32)
L̃2±=(ρ2)1/2cS2(±2s1(s3P)1/2s1(cS22sr2)sr(s3S)1/2s2(s3S)1/2cSsr±2s2(s3P)1/2s2(cS22sr2)sr(s3S)1/2±s1(s3S)1/2cSsr(cS22sr2)(s3P)1/22sr(s3S)1/20),
(A33)
Λ̃±=±iω(s3P000s3S000s3S),
(A34)

where sr is defined as in Eq. (A13), and the vertical slownesses s3P and s3S are defined as

s3P,S={1/cP,S2sr2forsr21/cP,S2,isr21/cP,S2forsr2>1/cP,S2.
(A35)

Here, cP and cS are the P- and S-wave velocities of the upper half-space, respectively, defined as cP=(λ+2μ)/ρ and cS=μ/ρ. On substitution of Eqs. (A32) and (A33) into Eqs. (A10) and (A11), we obtain

D̃1±(s)=ρcS2s3Ps3S+sr2(±((cS2s22)s3P+s22s3S)±s1s2(s3Ps3S)s1(cS22S)±s1s2(s3Ps3S)±((cS2s12)s3P+s12s3S)s2(cS22S)s1(cS22S)s2(cS22S)±s3ScS2),
(A36)
{Δ̃1(s)}1=12ρ(s12s3P+(1cS2s12)1s3S(1s3P1s3S)s1s20(1s3P1s3S)s1s2s22s3P+(1cS2s22)1s3S000s3P+sr2s3S),
(A37)
{Δ̃2(s)}1D̃2(s)=(120s1(ScS212)s3S012s2(ScS212)s3Ss1(ScS212)s3Ps2(ScS212)s3P12),
(A38)

where S=s3Ps3S+sr2.

1.
G.
Green
,
An Essay on the Application of Mathematical Analysis to the Theories of Electricity and Magnetism
(
Nottingham
,
1828
). [Reprinted in three parts in J. Reine Angew. Math.
39
(
1
),
73
89
(
1850
); 44(4), 356–374 (1852); and 47(3), 161–221 (1854).]
2.
L.
Challis
and
F.
Sheard
, “
The Green of Green functions
,”
Phys. Today
56
(
12
),
41
46
(
2003
).
3.
M.
Born
and
E.
Wolf
,
Principles of Optics
(
Pergamon
,
London
,
1965
), Chap. 8.
4.
J. W. S.
Rayleigh
,
The Theory of Sound. Volume II
(
Dover
,
New York
,
1878
), Chap. 14, reprinted 1945.
5.
N.
Bleistein
,
Mathematical Methods for Wave Phenomena
(
Academic
,
Orlando
,
1984
), Chap. 6.
6.
L.
Knopoff
, “
Diffraction of elastic waves
,”
J. Acoust. Soc. Am.
28
,
217
229
(
1956
).
7.
A. T.
de Hoop
, “
Representation theorems for the displacement in an elastic solid and their applications to elastodynamic diffraction theory
,” Ph.D. thesis,
Delft University of Technology
,
Delft
,
1958
.
8.
A. F.
Gangi
, “
A derivation of the seismic representation theorem using seismic reciprocity
,”
J. Geophys. Res.
75
,
2088
2095
, (
1970
).
9.
Y. H.
Pao
and
V.
Varatharajulu
, “
Huygens' principle, radiation conditions, and integral formulations for the scattering of elastic waves
,”
J. Acoust. Soc. Am.
59
,
1361
1371
(
1976
).
10.
J. A.
Kong
,
Electromagnetic Wave Theory
(
John Wiley & Sons
,
Hoboken
,
1986
), Chap. 5.
11.
C.
Altman
and
K.
Suchy
,
Reciprocity, Spatial Mapping and Time Reversal in Electromagnetics
(
Kluwer
,
Dordrecht
,
1991
), Chap 3.
12.
A. T.
de Hoop
,
Handbook of Radiation and Scattering of Waves
(
Academic
,
London
,
1995
), Chap. 7, 15, and 28.
13.
F. J.
Hilterman
, “
Three-dimensional seismic modeling
,”
Geophysics
35
,
1020
1037
(
1970
).
14.
L. N.
Frazer
and
M. K.
Sen
, “
Kirchhoff-Helmholtz reflection seismograms in a laterally inhomogeneous multi-layered elastic medium, I, Theory
,”
Geophys. J. R. Astron. Soc.
80
,
121
147
(
1985
).
15.
M.
Mansuripur
, “
A tutorial on the classical theories of electromagnetic scattering and diffraction
,”
Nanophotonics
10
,
315
342
(
2021
).
16.
R. P.
Porter
and
A. J.
Devaney
, “
Holography and the inverse source problem
,”
J. Opt. Soc. Am.
72
,
327
330
(
1982
).
17.
A. J.
Devaney
, “
A filtered backpropagation algorithm for diffraction tomography
,”
Ultrasonic Imaging
4
,
336
350
(
1982
).
18.
N. N.
Bojarski
, “
Generalized reaction principles and reciprocity theorems for the wave equations, and the relationship between the time-advanced and time-retarded fields
,”
J. Acoust. Soc. Am.
74
,
281
285
(
1983
).
19.
M. L.
Oristaglio
, “
An inverse scattering formula that uses all the data
,”
Inverse Probl.
5
,
1097
1105
(
1989
).
20.
R. P.
Porter
, “
Diffraction-limited, scalar image formation with holograms of arbitrary shape
,”
J. Opt. Soc. Am.
60
,
1051
1059
(
1970
).
21.
W. A.
Schneider
, “
Integral formulation for migration in two and three dimensions
,”
Geophysics
43
,
49
76
(
1978
).
22.
A. J.
Berkhout
,
Seismic Migration. Imaging of Acoustic Energy by Wave Field Extrapolation. A. Theoretical Aspects
(
Elsevier
,
Amsterdam
,
1982
), Chap. 7.
23.
J. D.
Maynard
,
E. G.
Williams
, and
Y.
Lee
, “
Nearfield acoustic holography: I. Theory of generalized holography and the development of NAH
,”
J. Acoust. Soc. Am.
78
,
1395
1413
(
1985
).
24.
C.
Esmersoy
and
M.
Oristaglio
, “
Reverse-time wave-field extrapolation, imaging, and inversion
,”
Geophysics
53
,
920
931
(
1988
).
25.
C.
Lindsey
and
D. C.
Braun
, “
Principles of seismic holography for diagnostics of the shallow subphotosphere
,”
Astrophys. J. Suppl. Ser.
155
,
209
225
(
2004
).
26.
M.
Fink
and
C.
Prada
, “
Acoustic time-reversal mirrors
,”
Inverse Probl.
17
,
R1
R38
(
2001
).
27.
A.
Derode
,
E.
Larose
,
M.
Tanter
,
J.
de Rosny
,
A.
Tourin
,
M.
Campillo
, and
M.
Fink
, “
Recovering the Green's function from field-field correlations in an open scattering medium (L)
,”
J. Acoust. Soc. Am.
113
,
2973
2976
(
2003
).
28.
K.
Wapenaar
, “
Synthesis of an inhomogeneous medium from its acoustic transmission response
,”
Geophysics
68
,
1756
1759
(
2003
).
29.
R. L.
Weaver
and
O. I.
Lobkis
, “
Diffuse fields in open systems and the emergence of the Green's function (L)
,”
J. Acoust. Soc. Am.
116
,
2731
2734
(
2004
).
30.
W. T.
Thomson
, “
Transmission of elastic waves through a stratified solid medium
,”
J. Appl. Phys.
21
,
89
93
(
1950
).
31.
N. A.
Haskell
, “
The dispersion of surface waves on multilayered media
,”
Bull. Seismol. Soc. Am.
43
,
17
34
(
1953
).
32.
F.
Gilbert
and
G. E.
Backus
, “
Propagator matrices in elastic wave and vibration problems
,”
Geophysics
31
,
326
332
(
1966
).
33.
B. L. N.
Kennett
, “
The connection between elastodynamic representation theorems and propagator matrices
,”
Bull. Seismol. Soc. Am.
62
,
973
983
(
1972
).
34.
B. L. N.
Kennett
, “
Seismic waves in laterally inhomogeneous media
,”
Geophys. J. R. Astron. Soc.
27
,
301
325
(
1972
).
35.
J. H.
Woodhouse
, “
Surface waves in a laterally varying layered structure
,”
Geophys. J. R. Astron. Soc.
37
,
461
490
(
1974
).
36.
A. J.
Haines
, “
Multi-source, multi-receiver synthetic seismograms for laterally heterogeneous media using F-K domain propagators
,”
Geophys. J. Int.
95
,
237
260
(
1988
).
37.
B. L. N.
Kennett
,
K.
Koketsu
, and
A. J.
Haines
, “
Propagation invariants, reflection and transmission in anisotropic, laterally heterogeneous media
,”
Geophys. J. Int.
103
,
95
101
(
1990
).
38.
K.
Koketsu
,
B. L. N.
Kennett
, and
H.
Takenaka
, “
2-D reflectivity method and synthetic seismograms for irregularly layered structures—II. Invariant embedding approach
,”
Geophys. J. Int.
105
,
119
130
(
1991
).
39.
H.
Takenaka
,
B. L. N.
Kennett
, and
K.
Koketsu
, “
The integral operator representation of propagation invariants for elastic waves in irregularly layered media
,”
Wave Motion
17
,
299
317
(
1993
).
40.
C. P. A.
Wapenaar
,
N. A.
Kinneging
, and
A. J.
Berkhout
, “
Principle of prestack migration based on the full elastic two-way wave equation
,”
Geophysics
52
,
151
173
(
1987
).
41.
J. H.
Rose
, “ 
‘Single-sided’ focusing of the time-dependent Schrödinger equation
,”
Phys. Rev. A
65
,
012707
(
2001
).
42.
J. H.
Rose
, “ 
‘Single-sided’ autofocusing of sound in layered materials
,”
Inverse Probl.
18
,
1923
1934
(
2002
).
43.
F.
Broggini
and
R.
Snieder
, “
Connection of scattering principles: A visual and mathematical tour
,”
Eur. J. Phys.
33
,
593
613
(
2012
).
44.
K.
Wapenaar
,
J.
Thorbecke
,
J.
van der Neut
,
F.
Broggini
,
E.
Slob
, and
R.
Snieder
, “
Green's function retrieval from reflection data, in absence of a receiver at the virtual source position
,”
J. Acoust. Soc. Am.
135
,
2847
2861
(
2014
).
45.
M.
Ravasi
,
I.
Vasconcelos
,
A.
Kritski
,
A.
Curtis
,
C. A.
da Costa Filho
, and
G. A.
Meles
, “
Target-oriented Marchenko imaging of a North Sea field
,”
Geophys. J. Int.
205
,
99
104
(
2016
).
46.
M.
Staring
,
R.
Pereira
,
H.
Douma
,
J.
van der Neut
, and
K.
Wapenaar
, “
Source-receiver Marchenko redatuming on field data using an adaptive double-focusing method
,”
Geophysics
83
,
S579
S590
(
2018
).
47.
X.
Jia
,
A.
Guitton
, and
R.
Snieder
, “
A practical implementation of subsalt Marchenko imaging with a Gulf of Mexico data set
,”
Geophysics
83
,
S409
S419
(
2018
).
48.
J.
Van der Neut
,
J. L.
Johnson
,
K.
van Wijk
,
S.
Singh
,
E.
Slob
, and
K.
Wapenaar
, “
A Marchenko equation for acoustic inverse source problems
,”
J. Acoust. Soc. Am.
141
,
4332
4346
(
2017
).
49.
C. P. A.
Wapenaar
, “
Reciprocity theorems for two-way and one-way wave vectors: A comparison
,”
J. Acoust. Soc. Am.
100
,
3508
3518
(
1996
).
50.
A. J.
Haines
and
M. V.
de Hoop
, “
An invariant imbedding analysis of general wave scattering problems
,”
J. Math. Phys.
37
,
3854
3881
(
1996
).
51.
K.
Wapenaar
, “
Unified matrix-vector wave equation, reciprocity and representations
,”
Geophys. J. Int.
216
,
560
583
(
2019
).
52.
A. T.
de Hoop
, “
Time-domain reciprocity theorems for acoustic wave fields in fluids with relaxation
,”
J. Acoust. Soc. Am.
84
,
1877
1882
(
1988
).
53.
M.
Schoenberg
and
P. N.
Sen
, “
Properties of a periodically stratified acoustic half-space and its relation to a Biot fluid
,”
J. Acoust. Soc. Am.
73
,
61
67
(
1983
).
54.
J. P.
Corones
, “
Bremmer series that correct parabolic approximations
,”
J. Math. Anal. Appl.
50
,
361
372
(
1975
).
55.
B.
Ursin
, “
Review of elastic and electromagnetic wave propagation in horizontally layered media
,”
Geophysics
48
,
1063
1081
(
1983
).
56.
L.
Fishman
and
J. J.
McCoy
, “
Derivation and application of extended parabolic wave theories. I. The factorized Helmholtz equation
,”
J. Math. Phys.
25
,
285
296
(
1984
).
57.
C. P. A.
Wapenaar
and
A. J.
Berkhout
,
Elastic Wave Field Extrapolation
(
Elsevier
,
Amsterdam
,
1989
), Chap. 3, 4, 11, and 12.
58.
M. V.
de Hoop
, “
Generalization of the Bremmer coupling series
,”
J. Math. Phys.
37
,
3246
3282
(
1996
).
59.
P. L.
Stoffa
,
Tau-p—A Plane Wave Approach to the Analysis of Seismic Data
(
Kluwer Academic
,
Dordrecht
,
1989
), Chap. 2.
60.
K.
Wapenaar
and
S.
de Ridder
, “
On the relation between the propagator matrix and the Marchenko focusing function
,”
Geophysics
87
(
2
),
A7
A11
(
2021
).
61.
H. A.
Lorentz
, “
The theorem of Poynting concerning the energy in the electromagnetic field and two general propositions concerning the propagation of light
,”
Versl. Afd. Natuurkd. K. Akad. Wet.
4
,
176
187
(
1895
).
62.
L.
Knopoff
and
A. F.
Gangi
, “
Seismic reciprocity
,”
Geophysics
24
,
681
691
(
1959
).
63.
A. T.
de Hoop
, “
An elastodynamic reciprocity theorem for linear, viscoelastic media
,”
Appl. Sci. Res.
16
,
39
45
(
1966
).
64.
E. G.
Flekkøy
and
S. R.
Pride
, “
Reciprocity and cross coupling of two-phase flow in porous media from Onsager theory
,”
Phys. Rev. E
60
,
4130
4137
(
1999
).
65.
B. A.
Auld
, “
General electromechanical reciprocity relations applied to the calculation of elastic wave scattering coefficients
,”
Wave Motion
1
,
3
10
(
1979
).
66.
J. D.
Achenbach
,
Reciprocity in Elastodynamics
(
Cambridge University Press
,
Cambridge
,
2003
), Chap. 14.
67.
S. R.
Pride
and
M. W.
Haartsen
, “
Electroseismic wave properties
,”
J. Acoust. Soc. Am.
100
,
1301
1315
(
1996
).
68.
E.
Slob
and
K.
Wapenaar
, “
Electromagnetic Green's functions retrieval by cross-correlation and cross-convolution in media with losses
,”
Geophys. Res. Lett.
34
,
L05307
, (
2007
).
69.
K.
Wapenaar
and
J.
van der Neut
, “
A representation for Green's function retrieval by multidimensional deconvolution
,”
J. Acoust. Soc. Am.
128
,
EL366
EL371
(
2010
).
70.
K.
Wapenaar
,
E.
Slob
, and
R.
Snieder
, “
Unified Green's function retrieval by cross correlation
,”
Phys. Rev. Lett.
97
,
234301
(
2006
).
71.
R.
Snieder
,
K.
Wapenaar
, and
U.
Wegler
, “
Unified Green's function retrieval by cross-correlation; connection with energy principles
,”
Phys. Rev. E
75
,
036103
(
2007
).
72.
E.
Slob
,
K.
Wapenaar
,
F.
Broggini
, and
R.
Snieder
, “
Seismic reflector imaging using internal multiples with Marchenko-type equations
,”
Geophysics
79
,
S63
S76
(
2014
).
73.
J. R.
Berryhill
, “
Wave-equation datuming before stack
,”
Geophysics
49
,
2064
2066
(
1984
).
74.
K.
Hokstad
, “
Multicomponent Kirchhoff migration
,”
Geophysics
65
,
861
873
(
2000
).
75.
A.
Curtis
and
D.
Halliday
, “
Source-receiver wavefield interferometry
,”
Phys. Rev. E
81
,
046601
(
2010
).
76.
D.
Halliday
,
A.
Curtis
, and
K.
Wapenaar
, “
Generalized PP + PS = SS from seismic interferometry
,”
Geophys. J. Int.
189
,
1015
1024
(
2012
).
77.
B. L. N.
Kennett
,
N. J.
Kerry
, and
J. H.
Woodhouse
, “
Symmetries in the reflection and transmission of elastic waves
,”
Geophys. J. R. Astron. Soc.
52
,
215
229
(
1978
).
78.
K.
Wapenaar
,
R.
Snieder
,
S.
de Ridder
, and
E.
Slob
, “
Green's function representations for Marchenko imaging without up/down decomposition
,”
Geophys. J. Int.
227
,
184
203
(
2021
).
79.
K.
Wapenaar
,
F.
Broggini
,
E.
Slob
, and
R.
Snieder
, “
Three-dimensional single-sided Marchenko inverse scattering, data-driven focusing, Green's function retrieval, and their mutual relations
,”
Phys. Rev. Lett.
110
,
084301
(
2013
).
80.
J.
Van der Neut
,
I.
Vasconcelos
, and
K.
Wapenaar
, “
On Green's function retrieval by iterative substitution of the coupled Marchenko equations
,”
Geophys. J. Int.
203
,
792
813
(
2015
).
81.
K.
Wapenaar
and
E.
Slob
, “
On the Marchenko equation for multicomponent single-sided reflection data
,”
Geophys. J. Int.
199
,
1367
1371
(
2014
).
82.
C. A.
da Costa Filho
,
M.
Ravasi
,
A.
Curtis
, and
G. A.
Meles
, “
Elastodynamic Green's function retrieval through single-sided Marchenko inverse scattering
,”
Phys. Rev. E
90
,
063201
(
2014
).
83.
M. S. R.
Kiraz
,
R.
Snieder
, and
K.
Wapenaar
, “
Focusing waves in an unknown medium without wavefield decomposition
,”
JASA Express Lett.
1
,
055602
(
2021
).
84.
L.
Diekmann
and
I.
Vasconcelos
, “
Focusing and Green's function retrieval in three-dimensional inverse scattering revisited: A single-sided Marchenko integral for the full wave field
,”
Phys. Rev. Res.
3
,
013206
(
2021
).
85.
J.
Brackenhoff
,
J.
Thorbecke
, and
K.
Wapenaar
, “
Virtual sources and receivers in the real Earth: Considerations for practical applications
,”
J. Geophys. Res.
124
,
11802
11821
, (
2019
).
86.
R.
Pereira
,
M.
Ramzy
,
P.
Griscenco
,
B.
Huard
,
H.
Huang
,
L.
Cypriano
, and
A.
Khalil
, “
Internal multiple attenuation for OBN data with overburden/target separation
,” in
SEG, Expanded Abstracts
(
2019
), pp.
4520
4524
.
87.
L.
Zhang
and
E.
Slob
, “
A field data example of Marchenko multiple elimination
,”
Geophysics
85
,
S65
S70
(
2020
).
88.
M.
Staring
and
K.
Wapenaar
, “
Three-dimensional Marchenko internal multiple attenuation on narrow azimuth streamer data of the Santos Basin, Brazil
,”
Geophys. Prosp.
68
,
1864
1877
(
2020
).
89.
C.
Reinicke
,
M.
Dukalski
, and
K.
Wapenaar
, “
Comparison of monotonicity challenges encountered by the inverse scattering series and the Marchenko demultiple method for elastic waves
,”
Geophysics
85
,
Q11
Q26
(
2020
).