Marchenko Green’s Function Retrieval in Layered Elastic Media from Two-Sided Reflection and Transmission Data

: By solving a Marchenko equation, Green’s functions at an arbitrary (inner) depth level inside an unknown elastic layered medium can be retrieved from single-sided reﬂection data, which are collected at the top of the medium. To date, it has only been possible to obtain an exact solution if the medium obeyed stringent monotonicity conditions and if all forward-scattered (non-converted and converted) transmissions between the acquisition level and the inner depth level were known a priori. We introduce an alternative Marchenko equation by revising the window operators that are applied in its derivation. We also introduce an auxiliary equation for transmission data, which are collected at the bottom of the medium, and a coupled equation, which is based on both reﬂection and transmission data. We show that the joint system of the Marchenko equation, the auxiliary equation and the coupled equation can be succesfully inverted when broadband reﬂection and transmission data are available. This results in a novel methodology for elastodynamic Green’s function retrieval from two-sided data. Apart from these data, our approach requires P - and S -wave transmission times between the inner depth level and the top of the medium, as well as two angle-dependent amplitude scaling factors, which can be estimated from the data by enforcing energy conservation.


Introduction
Inversion of the Marchenko equation has proven to be an effective tool for the retrieval of Green's functions in an unknown acoustic medium from single-sided reflection data [1,2].For an introduction to this subject, the numerical implementation of the Marchenko equation, field data applications and recent developments, see [3][4][5][6][7][8][9][10][11], respectively.An equivalent (Marchenko) equation has also been derived for wave propagation in elastic media [12][13][14].Inversion of this equation requires a priori knowledge of all forward-scattered (non-converted and converted) waveforms [15].Moreover, a unique solution can only be obtained if the medium obeys stringent monotonicity conditions [16], which are often not met in realistic scenarios.Once a solution to the Marchenko equation is found, it can be used for various purposes, such as wavefield retrieval inside an unknown medium [17], the imaging of elastic medium properties [18] or the suppression of multiple undesired reflections in reflection data [19].
In this paper, we show that the conditions for elastodynamic Green's function retrieval are significantly better when an elastic volume can be accessed from two sides, as is the case in particular laboratory experiments [20], non-destructive testing [21,22], brain imaging [23,24], transcranial ultrasound focusing [25,26], transcranial photoacoustics [27,28] and when using auxiliary downhole receivers in seismic data acquisition [29,30].Although the underlying representations of our work could be extended to account for lateral variations [19], the presence of a free surface [31,32] and intrinsic attenuation [33,34], we restrict ourselves to a layered lossless medium for simplicity.
In Section 2, we derive a system of forward equations that relate the (multi-component) focusing function at a specified focal depth z I to observed reflection and transmission data, which are to be acquired at depths of z U < z I and z L > z I .In Section 3, we show how the system can be inverted for the focusing function and two unknown amplitude scaling factors, α and β, which are related to transmission losses of (non-converted) Pand S-waves, respectively, between depth levels z I and z U .Apart from the recorded (reflection and transmission) data, our scheme requires two direct arrivals, which are represented by pulses of unit amplitude, delayed with the (non-converted) Pand S-wave travel times from z I to z U .In this way, we can apply exact (data-driven) Marchenko redatuming of two-sided data in a layered elastic medium, which is the main contribution of this paper.The retrieved focusing functions can be transformed into Green's functions as if there were virtual Por S-wave sources at z I , which could eventually be used for imaging and inversion of the elastic medium's properties.We close the paper with a discussion in Section 4.

Forward Equations
After providing some preliminaries in Section 2.1, we discuss the causality cones of multi-component Green's functions and focusing functions in Section 2.2.We propose novel window operators for Green's functions (based on non-converted P-wave travel times) and focusing functions (based on non-converted S-wave travel times).With the help of these operators, we derive (reflection-based) Marchenko equations in Section 2.3, (transmission-based) auxiliary equations in Section 2.4 and (transmission-and reflectionbased) coupled equations in Section 2.5.In Section 2.6, we take these equations together, leading to a joint system.Finally, we present a relation to convert focusing functions into Green's functions in Section 2.7.

Preliminaries
Let (x, y, z) be an Euclidean coordinate system with the z-axis pointing downwards, whereas t denotes time.We consider a layered lossless isotropic elastic medium, which is characterized by P-wave velocity c P (z), S-wave velocity c S (z) and mass density ρ(z).Let z U and z L be two depth levels, which are located above and below all heterogeneties in the medium, respectively (hence, constant medium properties are assumed above z U and below z L ).Elastodynamic wave propagation is considered in the (x, z)-plane, where wavefields are assumed to be constant in the y-direction.All wavefields are decomposed into fluxnormalized up-and downgoing P-, Svand Sh-components in the (p, τ)-domain [35-37], where p is the rayparameter and τ is the intercept time.The Sh-components are decoupled from the Pand Sv-components and will not be considered in this paper (for notational convenience, component Sv will be referred to as S).
In Figure 1, we show a layered elastic medium, which will be used throughout our paper as a running example.For convenience, we have chosen the vertical dimension of the model to be 1 m.However, all quantities can be rescaled to fit a particular application in, e.g., ultrasound or seismology applications.More information on the design of the medium and the parameters that are used for modeling are provided in Appendix A. Our objective is to retrieve the Green's responses at z U and z L to a Green's source at a specified depth level z I (see the dashed magenta line in Figure 1) from recorded (reflection and transmission) data.These Green's functions are represented by the following matrix: Here, G −+ U and G −− U are the upgoing (indicated by the first superscript −) Green's functions at z U from a down-and upwards radiating (indicated by the second superscript + or −) virtual source at z I , respectively.These quantities contain distinguished PP-, PS-, SPand SS-components and are organized as For the representations of Green's functions, we make use of so-called focusing functions [13], which are represented by the matrix Here, F − U and F + U are the up-and downgoing focusing functions (organized akin to Equation (2) at z U .These functions are defined in a fictitious medium where the halfspace below z I is homogeneous.They focus 'from above' at z I and continue as a downgoing wavefield below this depth level (see [15] for details).Similarly, F + L and F − L are the downand upgoing focusing functions at z L .These functions are defined in a medium where the halfspace above z I is homogeneous.These functions focus 'from below' at z I and continue as an upgoing wavefield above this depth level.Finally, Z is an operator that reverses the signs of p and τ.For example, applying this operator to F + U yields Our objective is to retrieve the focusing functions and Green's functions from reflection and transmission data, to be acquired at z U and z L .Let R U,PP , R U,PS , R U,SP and R U,SS be the PP-, PS-, SPand SS-reflection responses at z U (for their definitions, see Appendix B.1).Based on these recordings, we can construct an operator R U that convolves a wavefield with the reflection response.When applied to F + U , this multidimensional convolution is defined as Similar operators R L , T LU and T UL can be constructed for convolution with the reflection response at z L , the transmission response from z U to z L and the transmission response from z L to z U , respectively.Apart from R U , R L , T LU , T UL and Z, we make use of two window operators that will be defined in the following section.

Causality Cones of Green's Functions and Focusing Functions
In moderately inhomogeneous acoustic media, the Green's function and the focusing function are separated in time, except for a single overlapping event, which is commonly referred to as the direct wave.In the derivation of the acoustic Marchenko equation, this observation is exploited by truncating wavefields either before [2,11] or after [38,39] the direct wave.In elastic media, there can be a multitude of overlapping events, making the situation significantly more cumbersome [16].To illustrate this problem, we show the (symmetrized) causality cones of multi-component Green's functions and focusing functions in Figure 2. In particular, we refer the reader to the orange areas, where the Green's functions and focusing functions may overlap.Because of this potential overlap, we have designed two distinct time window operators: one for Green's functions, which is based on the (non-converted) direct P-wave travel time τ P d and one for focusing functions, which is based on the (non-converted) direct S-wave travel time τ S d .First, we discuss the window operator for Green's functions, which is based on the travel time τ P d of the (non-converted) P-wave, propagating from z I outwards.In Figure 2a, we can see that the Green's function and its time-reversed counterpart vanish in the interval L that reside at the boundary of the interval [−τ P d , τ P d ] (corresponding to the direct non-converted P-wave transmissions).Now, we can partition the Green's function that we defined earlier in Equation (1) as where Ld are referred to as the Green's function codasand O is a zero matrix.We design a window matrix Θ [P] that removes all data outside the interval [−τ P d , τ P d ].When we apply this matrix to the Green's function in Equation ( 6), it follows that In this formulation, Θ [P] U and Θ

[P]
L are operators that remove all data outside the intervals [−τ P Ud , τ P Ud ] and [−τ P Ld , τ P Ld ], respectively.Here, τ P Ud and τ P Ld are the direct P-wave travel times for propagation from z I to z U and z L , respectively.Ud and ±τ S Ld of the direct (non-converted) S-wave transmissions.All wavefields are stricly zero in the gray areas (whereas they may be non-zero in the yellow and orange areas).The areas where the Green's functions and focusing functions can overlap are indicated in orange.
We proceed with the window operator for focusing functions, which is based on the travel time τ S d of the (non-converted) direct S-wave, propagating from z I outwards.As illustrated in Figure 2b, the focusing functions and their time-reversed counterparts vanish L that reside at the boundary of the interval (−τ S d , τ S d ) (corresponding to the direct non-converted S-wave transmissions).Now, we may partition the focusing function that we defined earlier in Equation (3) as where Ld are referred to as the focusing function codas.We design a window matrix Θ (S) that removes all data outside the interval (−τ S d , τ S d ).During the inversion that will be applied later in this paper, we wish to restrict F m to the interval (−τ S d , τ S d ).To enforce this in practice, we replace F m in Equation (8) with In this formulation, Θ U and Θ (S) L are operators that remove all data outside the intervals (−τ S Ud , τ S Ud ) and (−τ S Ld , τ S Ld ), respectively.Here, τ S Ud and τ S Ld are the direct S-wave travel times for propagation from z I to z U and z L , respectively.

Marchenko Equations
In Appendix B.1, we derive the following system of Green's function representations that are based on reflection data: where I is a 2 × 2 identity matrix.When we apply the operator Θ [P] to both sides of this equation, it follows, with the help of Equation ( 7), that Next, we may substitute Equation ( 9) and rewrite the result as Here, we have used the fact that We refer to Equation (12) as a system of reflection-based Marchenko equations, which could be inverted for the unknown components of the focusing function F m .The blockdiagonal structure of this system reveals that the Marchenko equations at z U and z L are decoupled.Operator Θ (S) restricts the unknown focusing function , this leads to an underdetermined system of equations, which cannot be unconditionally inverted.We illustrate this by plotting the singular values of A Mar in Figure 3a (red curve) for data from the model that we presented above in Figure 1.In this case, A Mar contains 2192 columns but only 928 independent rows.To increase the rank of this matrix, we propose to use auxiliary transmission data, for which we derive a similar system of equations in the following section.

Auxiliary Equations
In Appendix B.2, we derive the following system of Green's function representations that are based on transmission data: When we apply operator Θ [P] to both sides of this equation, it follows, with the help of Equation (7), that When we substitute Equation ( 9), we find eventually that We refer to 15 as a system of auxiliary equations, which can be interpreted as a transmission-based inverse problem for F m .In Figure 3a (green curve), we show that the governing matrix A Aux of this problem is rank-deficient (at least for the medium in Figure 1).Nevertheless, this matrix can provide complementary information to A Mar , as illustrated in Figure 3b (orange curve).In this case, we have concatenated the rows of A Mar and A Aux , leading to a matrix of rank > 928.However, the rank is still far below 2192, which is the number of unknowns for this problem.In the next section, we show how we can improve on this by coupling reflection and transmission data, leading to yet another system of equations.

Coupled Equations
It is observed that the Green's function matrix G can be eliminated from the reflectionand transmission-based representations by subtracting Equation ( 13) from Equation (10).This leads to Here, we divided by a factor 2 to achieve a better amplitude balance with the matrices that were derived in the previous sections.After the substitution of Equation ( 9), we obtain We refer to 17 as a system of coupled equations, which can be interpreted as another inverse problem for F m .Although we have increased the number of rows signficantly (up to 16,384 in our running example) by not applying the window operator Θ [P] , the matrix A Cou is still rank-deficient, as shown in Figure 3a (blue curve).However, adding the rows of either A Mar or A Aux to the rows of A Cou results in a full-rank matrix, as illustrated by the magenta and cyan curves in Figure 3b.An intuitive understanding of this observation is that the subtraction of Equation ( 13) from Equation (10) (which was required for the construction of A Cou ) has reduced the row space of our system matrix, which can be compensated for by adding complementary rows from either the Marchenko or auxiliary system.

Joint System of Equations
Although the concatenation of matrix A Cou and either A Aux or A Mar seems sufficient by itself to construct a full-rank matrix, we choose to merge all three matrices, leading to the overall system As indicated by the black curve in Figure 3b, matrix A has full rank, and hence can be inverted.When we apply singular-value decomposition A = UΣV t and define the pseudo-inverse as A ‡ = VΣ ‡ U t (where Σ ‡ contains the reciprocals of all non-zero singular values), we may now write F m = A ‡ B. Akin to the acoustic Marchenko problem, a range of alternative solvers can be used to compute the pseudo-inverse [40,41].

Construction of Green's Functions from Focusing Functions
Either (the reflection-based) Equation ( 10) or (the transmission-based) Equation ( 13) can be used to convert focusing functions into Green's functions.Alternatively, we may take the average of both approaches, leading to We use this result later in this paper to construct Green's functions from (retrieved) focusing functions.

Inversion
In order to construct matrix B in Equation ( 18), we require a priori knowledge of four direct arrivals: Ud and F −SS Ld .In Section 3.1, we show that these arrivals can be expressed in terms of two travel times, τ P Ud (for P-wave transmission from z I to z U ) and τ S Ud (for S-wave transmission from z I to z U ), as well as two amplitude scaling factors, α and β.In Section 3.2, we present a procedure to estimate these scaling factors.In Section 3.3, we apply this procedure to retrieve focusing functions and Green's functions from numerical data.

Initialization
The direct arrival G −−PP Ud that is required for the construction of G −− Ud can be expressed in terms of a delayed unit pulse δ(τ − τ P Ud ) and an amplitude scaling factor, which we parameterize strategically as −α 1 2 for some (yet-unknown) α.This leads us to obtain The direct arrival, G ++PP Ld , that is required for the construction of G ++ Ld can be related to (with E P Ud as defined in Equation (20).This leads to Similarly, the direct wave F +SS Ud that is required for the construction of F + Ud can be expressed in terms of a time-advanced unit pulse δ(τ + τ S Ud ) and an amplitude scaling factor, which we parameterize strategically as β 1 2 for some (yet-unknown) β.This leads to The direct arrival F −SS Ld that is required for the construction of F − Ld can be related to Here, H +SS LUd is the first event of the SS-component of the inverse transmission response H + LU from z U to z L , which can be obtained via the inversion of as defined in Equation ( 22).This leads to We assume that the travel times τ P d and τ S d are known.The amplitude scaling factors α and β can be estimated from the data, as we discuss in the following section.

Estimation of Amplitude Scaling Factors α and β
In Appendix C, we show that the focusing function matrix F can be written as an explicit function of τ, α and β, according to In this formulation, k(τ) and l(τ) represent 2 × 1 vectors, which can be explicitly computed for each value of τ from the recorded data (see Appendix C for details).To estimate the amplitude scaling factors α and β, we make use of a relation for energy conservation [42,43], which has been used earlier for the estimation of amplitude scaling factors in acoustic media [44].This criterion can be written as When we substitute the expressions for F ± U (τ, α, β) from ( 24) into (25), subtract IS(τ) on both sides and evaluate the result at τ = 0, we find four expressions for α and β (i.e., the four entries of the 2 × 2 matrix equation).The first of these expressions (corresponding to the first diagonal entry of the matrix) is independent of β.We multiply this expression by α and define the left-hand side of the result as h P U (α).We find that In a similar way, the last of our four expressions (corresponding to the last diagonal entry of the matrix) is independent of α.We multiply this expression by β and define the left-hand side of the result as h S U (β).This leads to Two more expressions can be obtained by enforcing the energy conservation of the focusing function F L .We find, akin to Equation (25), that When we substitute the expressions for F ± L (τ, α, β) from Equation ( 24) into this result and repeat the abovementioned steps, we arrive at expressions for h P L (α) and h S L (β) (which are equivalent to Equations ( 26) and ( 27) with the subscript U replaced by L).The scaling factors α and β could be found by evaluating the roots of h P U (α), h P L (α), h S U (β) and h S L (β).However, we have chosen an alternative approach based on minimizing the cost functions and In Figure 4, we show both cost functions as computed from the numerical data of our running example.We can use Matlab's fminbnd routine to minimize these functions, yielding the estimates α ≈ 0.4716 and β ≈ 0.7946 (their true values being 0.4694 and 0.7934, respectively).Cost functions (a) J P (α) and (b) J S (β) for our numerical data.The dashed magenta lines denote the minima, 0.4716 (for α) and 0.7946 (for β), that were found using Matlab's fminbnd routine.The solid cyan lines denote the exact values 0.4694 (for α) and 0.7934 (for β), as extracted from the reference data.

Results
Now that α and β are resolved, the focusing function can be computed for our running example with the help of Equation (24).In Figures 5-8 we compare the PP-, SP-, PSand SS-components of the retrieved focusing functions with the results of direct modeling.We observe that all events have been recovered well, where the most significant differences (which can hardly be observed in the figure) can be attributed to the (small) errors in our estimates of α and β.Next, we compute the Green's functions from the retrieved focusing functions with the help of Equation (19).In Figures 9-12, we compare the PP-, SP-, PSand SS-components of the retrieved Green's functions with the results of direct modeling.Once more, we report an acceptable match, where the main differences (which can hardly be observed in the figures) can be attributed to errors in our estimates of α and β.

Discussion
Our methodology requires knowledge of the intercept times τ P Ud and τ S Ud , whereas the focal depth z I may be unknown.Hence, we can effectively retrieve Green's functions at a desired (intercept) time, even in the absence of velocity information [2].A velocity model is only required to convert intercept times into depths, akin to acoustic Marchenko imaging [45].An important observation in this context is that the construction of a virtual P-wave source is intrinsically decoupled from the construction of a virtual S-wave source in our formalism and both tasks may even be processed independently.Consequently, we may choose τ P Ud and τ S Ud at mutually different focal depths without affecting the accuracy of our results.Hence, we may conclude that neither c P (z) and c S (z), nor the ratio c P c S (z), is intrinsically required for the application of our methodology.
For our numerical simulations, we have designed a medium such that the arrival times of all waveforms coincide with exact time samples (see Appendix A).In this way, we could avoid problems related to discretization.In practical applications, data are recorded within a finite frequency band only, posing limitations to our resolution, especially in the presence of thin layers [2].It has been shown previously that some of these limitations can be overcome by enforcing energy conservation and minimum-phase conditions in the single-sided Marchenko equation [16,42,43].Similar strategies may be applied to the system of equations that we presented in this paper.
Although our methodology has been derived for a layered lossless medium with homogeneous halfspaces above z U and below z L , it could potentially be applied to a broader range of problems.Mild lateral variations of the medium's properties may be tolerable, akin to elastodynamic Marchenko imaging of single-sided data [12,18].The effects of dissipation might be incorporated by computing all correlation-based reflection and transmission operators in an effectual medium, akin to the equivalent two-sided acoustic problem [33,34].Heterogeneities above z U might be accounted for by convolving our representations with areal sources that take interactions with this part of the medium into account, akin to Rayleigh-Marchenko redatuming [46,47].A similar strategy may allow us to account for heterogeneities below z L .
The invertibility of matrix A in our formalism is likely to depend on the medium's properties, the available bandwidth and the wavefield components that can be emitted and recorded in practice.For transcranial applications, we may modify the theory [48] or apply redatuming [49] to account for spherical arrays.Moreover, it might be necessary to place our transducers in a water layer, which is impenetrable for S-waves.In terms of matrix algebra, such a configuration induces a projection of our multi-component wavefields to a reduced domain of P-waves only [50].It seems plausible that such a projection would affect the invertibility of A, but this remains to be investigated.

Conclusions
We have revised the window operators in the elastodynamic Marchenko equation.This leads to a system of equations that is intrinsically rank-deficient and hence cannot be solved without additional constraints.To overcome this issue, we have introduced an auxiliary equation (based on transmission data) and a coupled equation (based on reflection and transmission data).By concatenating these equations, we can construct a joint system for two-sided data that is invertible.Apart from the reflection and transmission data, this approach requires the direct (non-converted) Pand S-wave transmission times from the focal level to the upper acquisition array and two amplitude scaling factors, which can be retrieved by enforcing energy conservation.This leads to a methodology for velocityindependent true-amplitude Green's function retrieval in a lossless layered isotropic elastic medium from two-sided data.The methodology could potentially be extended to account for mild lateral variations, dissipation, anisotropy and heterogeneities above the upper acquisition array, as well as below the lower acquisition array.
where we have dropped subscript I (denoting the focal depth) for notational convenience.After taking the inverse Fourier transform (as defined in Equation (A3) and applying operator Z to the second and fourth rows, we obtain Equation (10), presented in the main text.

Appendix B.2. Transmission-Based Representations
We can derive a system of equivalent equations that are based on transmission data.First, we set m = U and n = I.We choose the actual medium in state A and place a source where we have dropped the subscript I once again for notational convenience.After taking the inverse Fourier transform (as defined in Equation (A3) and applying operator Z to the first and third row, we obtain Equation (13), presented in the main text.

Appendix C. Expression for F as a Function of τ, α and β
In this Appendix, we write F explicitly as a function of τ, α and β.We start with the substitution of Equations ( 20)- (23) into the definition of B Mar , which is given in the left-hand side of Equation (12).We write the result as where M I Aux and M I I Aux have been used to identify the constructed matrices.Finally, we may substitute Equations ( 22) and ( 23) into the definition of B Cou , which is given in the left-hand side of Equation ( 17).This yields where M I Mar and M I I Mar have been used to identify the constructed matrices.Next, we substitute Equations (A16)-(A18) into (18) and apply the pseudo-inverse of A to both sides of the result.This eventually leads to (A19) The direct focusing function matrix F d can be written in a similar form.This is achieved by substituting Equations ( 22) and ( 23) into the definition of F d , which is given in the right-hand side of Equation (8).The result can strategically be written as Here, we have indicated the arguments of all matrices for convenience to emphasize that we have expressed F explicitly as a function of τ, α and β.Note that K(τ) and L(τ) are independent ofthe scaling factors and hence can be computed from the recorded data, τ P d and τ S d .Finally, we can express our result as Equation (24), presented in the main text, by renaming the quantities that constitute matrices K(τ) and L(τ).

Figure 1 .
Figure 1.Example of a layered elastic medium with (a) P-wave velocity c P (in m•s −1 ), (b) S-wave velocity c S (in m•s −1 ) and (c) density ρ (in kg•m −3 ) as a function of depth z (in m).Above z U = 0 m and below z L = 1 m, the medium is homogeneous.The dashed magenta line indicates the focusing depth z I = 0.5 m, where a virtual source is to be constructed.

Figure 2 .
Figure 2. Symmetrized causality cone of (a) a Green's function (G or Z G) and (b) a focusing function for the medium in Figure 1 at p = 0.2 ms•m −1 , with the source/focal depth at z I = 0.5 m (indicated by the magenta dashed line; the black dashed lines indicate layer boundaries).The blue lines denote the travel times ±τ P Ud and ±τ P Ld of the direct (non-converted) P-wave transmissions.The red lines denote the travel times ±τ SUd and ±τ S Ld of the direct (non-converted) S-wave transmissions.All wavefields are stricly zero in the gray areas (whereas they may be non-zero in the yellow and orange areas).The areas where the Green's functions and focusing functions can overlap are indicated in orange.

Figure 3 .
Figure 3. (a) Singular values of the matrices A Mar (size: 928 × 2192), A Aux (size: 928 × 2192) and A Cou (size: 16,384 ×2192) for the model in Figure 1 at p = 0.2 ms•m −1 .(b) Singular values after concatenating various combinations of the matrices in (a).All curves have been normalized with respect to the highest singular value.The dots indicate the lowest singular values in the matrices.

Figure 4 .
Figure 4. Cost functions (a) J P (α) and (b) J S (β) for our numerical data.The dashed magenta lines denote the minima, 0.4716 (for α) and 0.7946 (for β), that were found using Matlab's fminbnd routine.The solid cyan lines denote the exact values 0.4694 (for α) and 0.7934 (for β), as extracted from the reference data.

Figure 5 .+ dτ 2 (Figure 6 .
Figure 5. PP-component of the retrieved focusing functions: (a) F − U , (b) F + U , (c) F + L and (d) F − L .The solid black traces were computed via direct modeling.The dashed orange traces were retrieved by means of our methodology.The blue lines have been drawn at ± τ P d + dτ 2 (where dτ = 2 µs denotes the intercept time sampling) to visualize the interval [−τ P d , τ P d ].The red lines have been drawn at ± τ S d − dτ 2 to visualize the interval (τ S d , τ S d ).

Figure 7 .
Figure 7. PS-component of the retrieved focusing functions (organized as in Figure5).

Figure 8 .
Figure8.SS-component of the retrieved focusing functions (organized as in Figure5).

Figure 9 .Figure 10 .
Figure 9. PP-component of the retrieved Green's functions: (a) G −+ U , (b) G −− U , (c) G +− L and (d) G ++ L .The solid black traces were computed via direct modeling.The dashed orange traces were retrieved using our methodology.The blue and red lines have been drawn at τ P d + dτ 2 and τ S d − dτ 2 , respectively.
Mar and M I I Mar have been used to identify the constructed matrices.In a similar way, we can substitute Equations (20)-(23) into the definition of B Aux , which is given in the left-hand side of Equation(15).This leads to

2 +
α, β) = D I (τ) + A ‡ (τ)M I (τ) D I I (τ) + A ‡ (τ)M I I (τ) . Here, T +PPLUd is the first event of the PP-component of the recorded transmission response T + LU from z U to z L .We assume that this event can be isolated from the transmission data by means of a time gate.Next, we define operator T P d for the 1D convolution of any signal with T +PP LUd .The direct Green's function at the lower level G ++PP Ud, which can be expressed in our notation as Z E P Ud α − 1 2 is the autocorrelation of the source signal s(τ), as defined in Appendix A. We assume that H +SS LUd can be isolated from H +SS