Focusing waves in an unknown medium without wavefield decomposition

: The Gel’fand-Levitan equation, the Gopinath-Sondhi equation, and the Marchenko equation are developed for one-dimensional inverse scattering problems. Recently, a version of the Marchenko equation based on waveﬁeld decomposition has been introduced for focusing waves in multi dimensions. However, waveﬁeld decomposition is a limitation when waves propagate horizontally at the focusing level. Here, the Marchenko equation for focusing without waveﬁeld decomposition is derived, and by iteratively solving the Marchenko equation, the Green’s function for an arbitrary location in the medium is retrieved from the scattered waves recorded on a closed receiver array and an estimate of the direct-wave without waveﬁeld decomposition. V C 2021 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) (http://creativecommons.org/licenses/by/4.0/) .


Introduction
Inverse scattering (Chadan and Sabatier, 1989;Colton and Kress, 1998;Gladwell, 1993) uses scattered waves to determine the scattering properties of a medium. Burridge (1980) shows that the Gel'fand-Levitan equation and the Gopinath-Sondhi equation have the same structure as the Marchenko equation, and shows that the Marchenko equation can be used for medium reconstruction (Burridge, 1980;Newton, 1980a). The solution of the one-dimensional (1D) Marchenko equation is an exact integral equation to make the connection between the scattered data and the scatterer potential. Rose (2001Rose ( , 2002 defines focusing as finding an incident wave that becomes a delta function at a prescribed focus location and time inside the medium. He shows that this incident wave follows from the scattered data and uses the Marchenko equation for 1D inverse scattering problems. Broggini and Snieder (2012) utilize Rose's approach and introduce a scheme in one dimension to retrieve the Green's function containing single-scattered and multiply scattered waves of the inhomogeneous medium. Wapenaar et al. (2012) show the virtual source creation in two dimensions using the recorded data but the proposed method excludes horizontally propagating energy at the virtual source level. Wapenaar et al. (2013Wapenaar et al. ( , 2014 derive the three-dimensional Marchenko equation for wavefield focusing and, therefore, for the Green's function retrieval; however, their solution requires up/down decomposition of the wavefield, which also excludes horizontally propagating energy at the focusing level. This is a limitation when the medium has steeply dipping structures because the horizontally scattered waves and refracted waves cannot be fully represented with the up/down decomposition. Recently, there have been several studies to address the limitation of the Marchenko method due to the up/down separation of the Marchenko equation. Kiraz et al. (2020) show wavefield focusing for an arbitrary point inside an unknown highly scattering inhomogeneous medium using the data acquired on a closed boundary. Diekmann and Vasconcelos (2021) and Wapenaar et al. (2021) present alternative approaches to Green's function retrieval without up/down decomposition each with their own pros and cons.
The Marchenko schemes proposed in one dimension provide an exact solution for focusing, and for Green's function retrieval in the medium. Green's function retrieval is of importance for imaging applications in many fields. The ability to focus waves opens up applications ranging from scattering kidney stones to performing imaging, monitoring in seismology, and non-destructive testing.
In this paper, we propose a two-dimensional (2D) Marchenko equation for focusing waves in a highly scattering inhomogeneous medium. We show that by iteratively solving the Marchenko equation, the Green's function for an arbitrary point in a a) Author to whom correspondence should be addressed. b) ORCID: 0000-0003-1445-0857. c) ORCID: 0000-0002-1620-8282. strongly scattering inhomogeneous medium can be retrieved without wavefield decomposition at the focal point. As opposed to the current Marchenko algorithms that use single-sided acquisition methods, we include waves propagating in all directions at the focal point using the contributions from a closed array. Our scheme is an extension of the 1D Marchenko algorithms proposed by Newton (1980a), Rose (2001Rose ( , 2002, and Broggini and Snieder (2012) into two dimensions.

Theory
Consider the acoustic wave equation where q is density, x is the angular frequency, f is the source term, p is pressure, and c is the velocity. We use the acoustic wave equation for a constant velocity and variable density in the numerical examples in this paper. We define the Green's function Gðx; x s ; tÞ as the solution to the wave equation LG ¼ dðx À x s ÞdðtÞ, with the differential operator L ¼ qr Á ðq À1 rÞ À c À2 @ 2 =@t 2 . Here, x s is the source location and the Green's function is the response to a source at x s recorded at the receiver location x. We use the following convention for the Fourier transform: f ðtÞ ¼ ð1=2pÞ Ð FðxÞ expðÀixtÞdx, where i is the imaginary unit. In the frequency domain Gðx; x s ; xÞ satisfies LG ¼ dðx À x s Þ, with the differential operator L ¼ qr Á ðq À1 rÞ þ x 2 =c 2 .
We describe an iterative solution to focus a wavefield in the medium to a pre-defined location at t ¼ 0 when injected into the medium. Our solution requires the direct-wave information modeled in the homogeneous medium (when q and c are constant) for a source at the focusing point x s . This is the known Green's function G 0 ðx; x s ; tÞ in a homogeneous medium. Sending this direct-wave back into the inhomogeneous medium from a circular receiver array with the radius R in a time-reversed order creates a focus at the focal point at t ¼ 0; however, in addition to the focal spot, other waves are present around the focusing point, and Fig. 1(a) shows the snapshot at t ¼ 0 of the time-reversed direct-wave injection into the heterogeneous medium shown in Fig. 2 (about which the details will be provided in Sec. 2.1). This shows that emitting the time-reversed direct-wave into an inhomogeneous medium does not restrict the focused field to the focusing point, and our goal is to remove the waves at other locations than the focusing point in Fig. 1(a). Figure 1(b) shows the snapshot at t ¼ 0 of the wavefield injection obtained by the iterative algorithm we propose. As shown in Fig.  1(b), our algorithm creates a wavefield that focuses to the pre-defined focal point, which acts as a virtual source, and suppresses other waves at t ¼ 0. We obtain our focusing wavefield by only using the direct-wave information modeled in the homogeneous medium (when q and c are constant) and the recorded scattering response. Unlike the conventional Marchenko methods, our method does not require the decomposition of the Marchenko equation to achieve focusing. In Sec. 2.1, we discuss the iterative Marchenko equation we propose for wavefield focusing and show how to obtain better focuses in the medium than one can achieve with the direct waves only.

Iterative scheme and the Marchenko equation
We define the ingoing wavefield, U in ðn 0 ; tÞ, and outgoing wavefield, U out ðn; tÞ, wheren 0 andn denote the locations on the circle with radius R; they are related via the scattering response Aðn;n 0 ; tÞ of the inhomogeneous medium. Following Broggini and Snieder (2012), Rose (2001Rose ( , 2002, Wapenaar et al. (2013), and Wapenaar et al. (2014), we design a wavefield that becomes a delta function at the focus location with an iterative scheme that relates the ingoing wave U in k to the outgoing wave U out k at iteration k as The ingoing and outgoing waves and the scattering operator are defined on a circle with radius R, but for brevity, we omit the parameter dependence on R in Eq.
(2). The iterative scheme starts with injecting a delta function into the medium and the ingoing wavefield for the first iteration gives where t d ðn 0 Þ is the arrival time of the direct waves that propagates from the focusing point to the pointn 0 on the circle. Following Broggini and Snieder (2012), the purpose of the iterative scheme is to reconstruct a wavefield that after interacting with the heterogeneities in the medium collapses onto a delta function at the focusing point at t ¼ 0. We create a symmetric field in time for Àt d ðn 0 Þ < t < t d ðn 0 Þ. We later show that the symmetry in time leads to focusing. To achieve the symmetry for the iterative scheme, we define the ingoing wavefield as U in k ðn 0 ; sÞ ¼ U in 0 ðn 0 ; sÞ À Hðn 0 ; sÞU out kÀ1 ðn; ÀsÞ; where Hðn 0 ; sÞ is a window function and defined as Hðn 0 ; sÞ ¼ 1 when Àt d ðn 0 Þ < s < t d ðn 0 Þ, and otherwise Hðn 0 ; sÞ ¼ 0.
When the iterative scheme converges (hence when U out k ¼ U out kÀ1 ), the iteration number can be dropped. Inserting Eq. (4) into Eq. (2) then gives U out ðn; tÞ ¼ þ ð Aðn;n 0 ; t À sÞU in 0 ðn 0 ; sÞdsdn 0 À þ ð t d Àt d Aðn;n 0 ; t À sÞU out ðn 0 ; ÀsÞ dsdn 0 ; with t d ¼ t d À , where we introduce as a small positive constant to exclude the direct-wave at t d . If we define K ¼ ÀU out and substitute this into Eq. (5) using Eq. (3), we obtain Aðn;n 0 ; t À sÞKðn 0 ; ÀsÞ dsdn 0 ¼ 0: Burridge (1980) shows that the 1D Marchenko equation, Gel'fand-Levitan equation, and the Gopinath-Sondhi equations of inverse scattering can be written in symbolic notation as shows the time interval, R is the recorded data, and K is the function we solve for. Eq. (6) has the same structure as the equations derived by Burridge (1980) and, therefore, gives a 2D Marchenko equation without using up/down decomposition. Eq. (6) also has a similar relation with the equations derived by Newton (1980bNewton ( , 1981Newton ( , 1982 using the scattering data in multi-dimensional media.

Numerical example and Green's function retrieval
We illustrate our method with a 2D numerical example. Figure 2 shows the source and receiver geometry of a 2D acoustic medium. The red asterisk in Fig. 2 denotes the virtual source location and the blue line represents a circle on which 400 equidistant sources and receivers are placed. The virtual source location, x s ¼ ðx; zÞ, is at x ¼ 4 cm and z ¼ 0.8 cm. The medium has a constant background velocity and density, c 0 ¼ 2 km/s and q 0 ¼ 2 g/cm 3 , respectively. Figure 2 also shows four different elliptical-shaped scatterers located in the medium with densities q 1 ¼ 4.5 g/cm 3 , q 2 ¼ 5 g/cm 3 , q 3 ¼ 7.5 g/cm 3 , Fig. 2. Geometry of the 2D model. Sources and receivers (400 each) are located on the blue circle, and the red asterisk shows the virtual source location x s . The elliptical scatterers have contrasting densities which are given on the right-hand side. ARTICLE asa.scitation.org/journal/jel and q 4 ¼ 6 g/cm 3 , respectively. We use finite-difference modeling with absorbing boundaries and the source wavelet is a Ricker wavelet (Ricker, 1953) with a central frequency of 2 MHz. A challenge of the used geometry is that the focusing point is located inside one of the scatterers, which has a reflection coefficient of about 40% at the boundaries. As a result, the source generates strong reverberations within the scatterer (see supplementary material, 1 movie 1).
The ingoing wavefield in the finite-difference modeling can be implemented by either changing the finitedifference stencil at the circular array, or by using the equivalent sources f (in equation form) in the acoustic wave equation (1) to produce the desired ingoing wavefield. We use the equivalent sources in the acoustic wave equation (1) for the finite-difference implementation where the equivalent sources are given by the normal derivative of the ingoing wavefield (see supplementary material 1 ). To solve the Marchenko equation iteratively, we start with U in 0 ðn 0 ; tÞ ¼ U d ðn 0 ; ÀtÞ where U d is the time-reversed direct-wave in the homogeneous background medium. We send the ingoing wave U in 0 from the receiver array into the medium and use the outgoing wave recorded at the array in Eq. (4) to determine the ingoing wave for the next iteration. We use seven iterations to get close to convergence but more iterations might be needed for more complicated media where velocity and density are varying.
We next inject the wavefield obtained by the iterative solution on the boundary. Figure 3(a) shows the total wavefield, U total ðn 0 ; tÞ ¼ U in ðn 0 ; tÞ þ U out ðn 0 ; tÞ, recorded on the boundary for the 7th iteration, which consists of the superposition of the ingoing and outgoing wavefield. The wavefield in Fig. 3(a) is symmetric in time for Àt d ðn 0 Þ < t < t d ðn 0 Þ (approximately between À5 and 5 ls). If we take the difference between the total wavefield in Fig. 3(a) and its time-reversed version, i.e., U total ðn 0 ; tÞ À U total ðn 0 ; ÀtÞ, all events in the interval Àt d ðn 0 Þ < t < t d ðn 0 Þ vanish as shown in Fig. 3(b). A small amount of energy remains in Fig. 3(b) for Àt d ðn 0 Þ < t < t d ðn 0 Þ, this is due to numerical inaccuracies in our solution of the Marchenko equation. Since U total ðn 0 ; tÞ À U total ðn 0 ; ÀtÞ is anti-symmetric in time, it vanishes for t ¼ 0, also after injecting it into the medium. Hence, we diagnose the focusing by showing the time derivative ð@=@tÞðU total ðn 0 ; tÞ À U total ðn 0 ; ÀtÞÞ, injected into the medium. Figure 3(b) shows that for positive times, the wavefield U total ðn 0 ; tÞ À U total ðn 0 ; ÀtÞ vanishes at the receivers for t < t d ðn 0 Þ. If we consider this wavefield at t ¼ 0, the direct waves radiated at t ¼ 0 from x s arrive at a receiver location Rðn 0 Þ at t d ðn 0 Þ. Suppose that waves would radiate at t ¼ 0 from a point x 6 ¼ x s . For some receivers, those waves would arrive at a time t < t d ðn 0 Þ; however, as shown in Fig. 3(b), no waves arrive at time t < t d ðn 0 Þ. This means that waves do not radiate from any point x 6 ¼ x s at t ¼ 0. Therefore, the time-derivative of the wavefield U total ðn 0 ; tÞ À U total ðn 0 ; ÀtÞ, injected into the medium, is only non-zero at t ¼ 0 at the point x s , and the wavefield focuses at t ¼ 0 at the virtual source location (also see supplementary material 1 ).
We let pðx; tÞ denote the total wavefield in the interior that is associated with the wavefield U total ðn 0 ; tÞ on the boundary, and pðx; ÀtÞ denote the time-reversed version of this wavefield. The homogeneous Green's function [G h ðx; x s ; tÞ ¼ Gðx; x s ; tÞ À Gðx; x s ; ÀtÞ] (Oristaglio, 1989), for the virtual source location x s and the receiver location x is, up to a multiplicative constant, obtained from (also see supplementary material 1 ) G h ðx; x s ; tÞ ¼ pðx; tÞ À pðx; ÀtÞ: If we want to focus a wavefield at the virtual source location where there is no actual source located, we must have a nonzero incident wavefield. The causal and acausal Green's functions satisfy the inhomogeneous acoustic wave equation, but the homogeneous Green's function G h satisfies the homogeneous wave equation (Oristaglio, 1989). Equation (7), therefore, retrieves the Green's function for t > 0 for the virtual source location x s . Unlike other (interferometric) Green's function retrieval methods (Campillo and Paul, 2003;Duroux et al., 2010;Roux et al., 2004;Sabra et al., 2005;Schuster, 2009;Snieder and Larose, 2013;Wapenaar et al., 2005;Weaver and Lobkis, 2001), no physical receiver is required at the position of the virtual source, and unlike other Marchenko methods (Wapenaar et al., 2013(Wapenaar et al., , 2014, we do not rely on an up/down decomposition of the wavefield. When one applies the Marchenko algorithm to two points in the interior, one obtains the  Green's function for these two points recorded on the boundary. Using interferometric techniques, these Green's functions can be used to reconstruct the Green's function for waves propagating between two points in the interior (Brackenhoff et al., 2019;Singh and Snieder, 2017). Figure 3(c) shows the Green's function obtained from Eq. (7) with x taken at the boundary (blue lines), superimposed on the directly modeled Green's function (red lines). For clarity, the traces have been multiplied by exp(2t) to emphasize the scattered waves. The latest arrival time for the single-scattered waves for our geometry is about 18 ls. All waves arriving after 18 ls, therefore, are multiply scattered waves. For earlier times, the Green's function consists of a combination of single-scattered waves and multiply scattered waves. As a result of our iterative solution, we retrieve the directwave and the scattered waves. Figure 4 shows normalized vertical cross-sections of the wavefield at t ¼ 0 taken from Figs. 1(a) and 1(b) for x ¼ 4 cm. The red trace denotes the cross section of Fig. 1(a) and the blue trace denotes the cross section of Fig. 1(b). The snapshots (see Fig. 1) and the cross-sections (see Fig. 4) show that the reconstructed Green's function creates a focus only around the focusing point and cancels other arrivals around the focusing point to a large extent, whereas the results one can achieve with using only direct waves contain other arrivals that distort the focusing.

Conclusion
We derive the 2D Marchenko equation for wavefield focusing and Green's function retrieval for an arbitrary point in an unknown highly scattering inhomogeneous medium with a closed receiver array. We successfully retrieve the Green's function for a pre-defined location and the comparison to the directly modeled Green's function is found to be excellent [see Fig. 3(c)]. The cross-sections in Fig. 4 show that we can create better focusing in the medium than one can achieve with the direct waves only. Our retrieved Green's function contains both the single-and multiply scattered waves of the heterogeneous medium model. Because we use a constant background velocity model, our method requires the direct-wave information modeled only in the homogeneous medium (when q is constant), and the recorded scattering response Aðn;n 0 ; tÞ to solve the Marchenko equation iteratively like other multi-dimensional Marchenko methods proposed earlier (Wapenaar et al., 2013(Wapenaar et al., , 2014; however, it does not require wavefield decomposition. We show that after the convergence, we retrieve the Green's function for any desired location in the medium without relying on prior information about the scatterers in the medium and wavefield decomposition to solve the Marchenko equation. The Marchenko equation we propose forms the basis for imaging the interior of a medium inside a closed array without up/down decomposition and makes the Marchenko methods more appropriate for imaging steeply dipping structures.