Velocity-independent Marchenko focusing in time- and depth-imaging domains for media with mild lateral heterogeneity

The Marchenko method retrieves Green’s functions between the acquisition surface and any arbitrary point in the medium. The process generally involves solving an inversion starting with an initial focusing function, e.g., a direct-wave Green’s function from the desired subsurface position, typically obtained using an approximate velocity model. We have formulated the Marchenko method in the time-imaging domain. In that domain, we recognize that the traveltime of the direct-wave Green’s function is related to the Cheop’s traveltime pyramid commonly used in time-domain processing, which in turn can be readily obtained from the local slopes of the common-midpoint gathers. This observation allows us to substitute the velocity-model-based initial focusing operator with that from a data-driven slope estimation process. Moreover, we found that working in the time-imaging domain allows for the specification of the desired subsurface position in terms of vertical time, which is connected to the Cartesian depth position via the time-to-depth conversion. Our results suggest that the prior velocity model is only required when specifying the position in depth, but this requirement can be circumvented by making use of the time-imaging domain within its usual assumptions (e.g., mild lateral heterogeneity). Provided that those assumptions are satisfied, the estimated Green’s functions from the proposed method have comparable quality to those obtained with the knowledge of a prior velocity model.


INTRODUCTION
Green's functions between the surface and any subsurface point are the main ingredient in seismic redatuming and imaging. The Marchenko method is a framework to obtain such information using solely the reflection data at the surface and an initial estimate of the direct-wave Green's function from the desired subsurface position da Costa Filho et al., 2014;Slob et al., 2014;Wapenaar et al., 2014aWapenaar et al., , 2014bWapenaar et al., , 2017Ravasi, 2017;Singh et al., 2017). Given a prior approximate (smooth) velocity model of the subsurface in Cartesian coordinates, one can specify a subsurface position and obtain the direct-wave Green's function from this position to the surface using conventional forward extrapolation (Wapenaar et al., 2014b;Thorbecke et al., 2017). An alternative strategy involves a separate inversion for the direct-wave Green's function from the common focal point technology based on the same starting velocity model (Berkhout, 1997;Thorbecke, 1997). Therefore, much like conventional seismic migration, a caveat to the current Marchenko imaging implementation is the requirement of a priori velocity knowledge.
Conventional seismic imaging can be accomplished in either the time or the depth domain. The former generally performs with higher computational efficiency but becomes less accurate than the latter when dealing with geologically complex areas such as subsalt regions (Yilmaz, 2001). The shortcomings of timeimaging methods are largely due to the following (Fomel, 2013(Fomel, , 2014 1) Approximate direct-wave Green's functions are used for imaging, which typically depend on hyperbolic or nonhyperbolic traveltime approximations. 2) Each time-domain image point is associated with its own approximate effective velocity under the assumption of straightray geometry relative to the surface. 3) When lateral heterogeneity is present, the final images are generated in distorted coordinates defined by image rays (Hubral, 1977) as shown in Figure 1.
However, in areas with moderately complex geology where such assumptions are approximately valid, we can turn these limitations into advantages. In particular, recent research on time-domain imaging has led to an alternative data-driven time-imaging workflow for improved efficiency and accuracy with local event slopes from common-midpoint (CMP) gathers instead of velocity (Fomel, 2007). This development leads to an opportunity of relying on these velocity-independent data-driven techniques from time-domain imaging to estimate the direct-wave Green's functions that can be used in a new form of the Marchenko method.
In this paper, we first study the Marchenko method in the timeimaging domain and establish relationships between the focusing functions obtained from the Marchenko methods in the timeand depth-imaging domains. Making use of the slope-based time-domain processing workflow, we subsequently propose a scheme to obtain the direct-wave Green's function from the desired subsurface position on a reflector to the surface. We show that the newly estimated direct-wave Green's function can be used in the Marchenko method and leads to comparable results to those that rely on the prior knowledge of a velocity model. We rely on two synthetic models with and without lateral heterogeneity and a field data set from the North Sea to demonstrate the proposed approach.

Reciprocity theorems on curvilinear surfaces
The key components to deriving the single-sided Green's function representations for the Marchenko method are the one-way reciprocity theorems of convolution and correlation type (Wapenaar and Grimbergen, 1996;Wapenaar et al., 2014a;van der Neut et al., 2015). Assuming that the image rays are well defined with no caustics, we first recognize that a constant depth level in the Cartesian coordinates generally corresponds to a curved level in the timeimaging domain and vice versa ( Figure 2). In other words, the current Marchenko method has already been implemented with respect to a curved level in the time-imaging domain. To show that a converse relationship exists, we need to find the single-sided representations for a curvilinear level in depth that corresponds to a constant time surface.
Because the time-imaging domain is defined by image rays (Figure 1), it represents a special curvilinear coordinate system of the semiorthogonal type (Sava and Fomel, 2005) due to the orthogonality between the ray direction and the wavefront (Figure 2). In such curvilinear systems, we denote coordinates along the curved boundary surfaces as ξ ¼ ðξ 1 ; ξ 2 ; ξ 3 Þ and assume that there exists a one-to-one reversible map between the ray coordinates ξ and the physical Cartesian coordinates. We use ξ in our subsequent derivation to ensure its generality and will specify the meaning to each time-domain coordinate ξ i in a later section.
The one-way reciprocity theorems for semiorthogonal curvilinear systems in the Fourier ω domain are given by (Frijlink and Wapenaar, 2010): where equations 1 and 2 represent the one-way reciprocity theorems of convolution and correlation type, respectively. Here, p denotes the flux-normalized wavefields in the frequency domain decomposed into upgoing ð−Þ and downgoing ðþÞ constituents with respect to ξ 3 at the acquisition surface S a and the focusing surface S f . The superscript Ã denotes the complex conjugation. The integrations are done along the coordinates ξ a and ξ f at surfaces S a and S f , which no longer need to represent constant-depth surfaces. Similar to the current derivation of the Marchenko method, the subscripts A and B denote the two acoustic statesthe truncated medium and the true medium, respectively ( Figure 3). Moreover, the considered volume between S a and S f is assumed to have equal medium parameters in both states and is source free.

Depth coordinates Time coordinates
Surface of variable depth Surface of constant time = image wavefront Figure 2. The relationship between the variable-depth surface in Cartesian coordinates and the constant time surface in the time-imaging (image-ray) coordinates. The solid arrow lines denote image rays that originate from the acquisition surface S a at right angles. The focusing surface of S f is chosen as the surface of equal imageray traveltime (image wavefront) that gets mapped to a flat surface in the time-imaging domain.

Image ray
Depth coordinates Time coordinates Figure 1. The relationship between depth-and time-imaging coordinates. An example image ray originating from x a with an orthogonal slowness vector to the surface is shown. Every point along this ray is mapped to the same lateral distance in the time domain with a different traveltime t a .

Marchenko equations in image-ray coordinates
We can clearly observe that equations 1 and 2 are similar to the one-way reciprocity theorems for the case of flat datum levels except that the integrations are now done over curvilinear surfaces. Therefore, we can follow a similar procedure as in the previous works and derive the single-sided Green's function representations (Wapenaar et al., 2014a(Wapenaar et al., , 2014bvan der Neut et al., 2015). We first consider two acoustic states shown in Figure 3 and specify the pertaining parameters in equations 1 and 2 as follows: representing a point source for a downgoing wavefield just above S a .
ωÞ representing the reflected response in the truncated medium from ξ a to ξ A a .
2) State A at S f : ωÞ representing the transmission response from the surface S a to S f .
representing a point source for a downgoing wavefield just above S a .
ωÞ representing the total reflected response in the true medium from ξ a to ξ B a .

4)
State B at S f : ωÞ representing the downgoing part at S f of the Green's function between S a and S f .
ωÞ representing the upgoing part at S f of the Green's function between S a and S f . The superscripts in ξ A a and ξ B a are used to distinguish between the coordinates at S a of state A and B. The ξ a in the argument of R and R A is used to denote a general position on S a . From this point onward, we will also omit the argument ω for conciseness. The above substitution leads to To retrieve the Green's function from the two above equations, we define the focusing function f þ 1 as the inverse of transmission T in the truncated medium (state A) as follows: where ξ A f denotes the desired focusing position on the surface S f in the truncated medium (state A). As a result, f þ 1 ðξ A a ; ξ A f Þ represents a purpose-built downgoing field that travels from ξ A a on the surface S a to focus at ξ A f on S f in the truncated medium. We apply the focusing function (equation 5) to equations 3 and 4, and we integrate over S a , which leads to the desired one-sided integral representations: where we define the reflected response f − 1 of the truncated medium to f þ 1 as In other words, f − 1 ðξ a ; ξ A f Þ represents the upgoing reflected response of the truncated medium at ξ a on the surface S a to the propagation of f þ 1 ðξ A a ; ξ A f Þ that is designed to focus at ξ A f on S f . From equations 6-8, we can see how different ξ variables in the single-sided Green's function representations have emerged. For example, g − ðξ A f ; ξ B a Þ corresponds to the Green's function from ξ B a in the true medium (state B) to the focusing position ξ A f defined in state A according to equation 5. For simplicity, we recast equations 6 and 7 by switching the notation according to ξ B a ¼ ξ a , ξ A f ¼ ξ f and changing the dummy integrating variable to ξ 0 a . This leads to the following final system of equations in the frequency domain:

Reflection-free
State A State B Figure 3. The two acoustic states A and B used in the derivation of the system of Marchenko equations (equation 9). The setup is similar to that for the conventional Marchenko derivation in the Cartesian coordinates, except that the surfaces S a and S f are flat in some choice of semiorthogonal curvilinear coordinates, which are the image-ray coordinates in our case. State A is identical to state B above S f but has a reflection-free region (homogeneous) below. The reflection and transmission responses (operators) in state A are denoted by R A and T, respectively. State B is the true medium, where the region below S f can be heterogeneous and cause upward scattering. The term R denotes the total reflection response related to the originally recorded data at the surface S a . The desired down-and upgoing Green's functions at S f are denoted by g þ and g − , which can be obtained from solving the Marchenko system.
which has the same integral form as the original single-sided Green's function representations in Cartesian depth coordinates, except for the integration over the curvilinear boundary S a . The variables ξ a and ξ 0 a are defined on the surface S a , whereas ξ f is defined at the focusing level S f . In case that the acquisition surface S a is flat A similar reflection response R to the constant-depth case can be used, and we can recast the fully curvilinear system in equation 9 as which has a similar form to the Cartesian-space one-sided representations of the conventional Marchenko approach (e.g., Wapenaar et al., 2014b;van der Neut et al., 2015). We do, however, point out two key points that differentiate this system from its original counterpart: • In equation 10, the focusing functions and Green's functions retain the conventional Cartesian x a as their (data-related) surface argument, but the focal-point argument is ξ f which here lies on the curvilinear depth datum S f connected to a fixed time-domain datum ( Figure 2; more on this point below).
• As a result, the AE superscripts in this theoretical construct correspond to propagation directions normal to the curvilinear surface local to ξ f , meaning that the ideal g AE differs in radiation from their conventional Marchenko counterparts at some fixed ξ f ¼ x f . In practice, whether this radiation distinction expresses itself on the retrieved fields is entirely dependent on how one chooses to solve the focusing problem numerically.
Though these are notable differences, it is also very convenient that in the integral kernels of equation 10, we see the conventional Cartesian reflection response Rðx a ; x 0 a Þ. Because of it, this final system of coupled Marchenko equations allows us to solve for the time-domain-coordinate focusing functions from the usual, datarelated, reflection response R, by applying causality constraints to equation 10 as described by Wapenaar et al. (2014b) and van der Neut et al. (2015).
Because equation 10 is similar to that in the case of a constantdepth focusing level, we can argue that the form of the corresponding Marchenko equations remains the same as long as there exists a transformation between the Cartesian coordinates and some semiorthogonal curvilinear coordinates, whose level curve matches the desired curvilinear datum level. In the case of the time-imaging domain, the coordinate transformation is defined by the mapping of image rays for the time-to-depth conversion (Cameron et al., 2007;Iversen and Tygel, 2008;Sripanich and Fomel, 2018). In this paper, we define the image-ray coordinates as ξ 1 ¼ x 0 , ξ 2 ¼ y 0 , and ξ 3 ¼ t 0 (Figure 1). The first two coordinates x 0 and y 0 define the escape location of the image rays at the acquisition surface S a and t 0 is their one-way traveltime. The curved datum level at depth then corresponds to an image-wavefront surface tied to some constant time t 0 (Figure 2). Given the same focusing functions at the acquisition surface, the focusing position defined in the time-imaging domain (x 0 , y 0 , and t 0 ) can be translated to its corresponding Cartesian position through the same mapping. Equipped with these results, we can proceed with making use of efficient time-domain techniques to solve the Marchenko equations and obtain focusing functions associated with some specified position in the time-imaging domain.

SLOPE-BASED TIME-DOMAIN PROCESSING
The process of time-domain imaging can be conceptually summarized as shown in Figure 4 (Fowler, 1997). The recorded CMP data, which are dependent on midpoint m, half-offset h, and time t are first migrated to zero-offset (x n , t n ) through a combination of normal and dip moveout operations. Poststack time migration subsequently maps the result to (x 0 , t 0 ) to correct subsurface reflectionpoint positions. The entire process constitutes the time-imaging routine, and it is equivalent to the prestack time migration process stacking along the Cheops traveltime pyramid (Claerbout, 1996). Fomel (2007) shows that under the regular assumptions of hyperbolic traveltime and straight-ray geometry of time imaging, prestack time migration (mapping) can be done with local event slopes as follows: where t is the two-way reflection traveltime. The terms p m ¼ ∂t∕∂m and p h ¼ ∂t∕∂h are estimated local event slopes from the CMP gathers in the midpoint and half-offset directions, respectively. From Figure 4, we can deduce that the traveltime of direct-wave Green's function from the same subsurface position is represented by the same value of t 0 in this mapping (migration) process. Consequently, the desired traveltime of the direct-wave Green's function from that location is a contour of time t 0 of the one-way traveltime map (equation 11). Therefore, we can summarize the steps to obtain slope-based direct-wave Green's function and solve the Marchenko system as follows: 1) Given the CMP gathers, measure the slopes of primaries using methods such as plane-wave destruction (Fomel, 2002), which iteratively solves a partial differential equation that governs a local plane wave. 2) Generate the traveltime t 0 ðm; h; tÞ and distance x 0 ðm; h; tÞ maps according to equations 11 and 12. 3) Remap the t 0 data according to the x 0 data to correctly position the traveltime across the midpoint coordinate and obtain t 0 ðx 0 ; h; tÞ. 4) Specify a desired focusing position in (x 0 , t 0 ) and obtain the traveltime of the direct-wave Green's function from the corresponding contour of t 0 ðx 0 ; h; tÞ. 5) Run a hyperbolic regression using picked t 0 to ensure a smooth hyperbolic curve on discretized grids and obtain a best-fitted velocity. 6) Generate the hyperbolic traveltime curve from the obtained velocity and convolve it with a zero-phase wavelet. 7) Scale the amplitudes with appropriate geometric spreading factor for straight-ray geometry to obtain an approximate directwave Green's function. 8) Use the time-reversed direct-wave Green's function as the initial focusing function and solve the time-constrained Marchenko system that follows from equation 10 (Wapenaar et al., 2014b;van der Neut et al., 2015).

EXAMPLES
In this section, we first apply the proposed workflow to two cases of synthetic horizontally layered and laterally heterogeneous media. The primary-only CMP gathers to be used for slope estimation are generated with Kirchhoff modeling based on an accurate two-point ray tracing scheme (Sripanich and Fomel, 2014), whereas the full data set with all orders of multiples used in the construction of the reflection operator R in equation 10 is obtained from finite-difference modeling. We subsequently show a field-data example from the North Sea previously studied by Szydlik et al. (2007) and Ravasi et al. (2015Ravasi et al. ( , 2016. In practice, we note that the required CMP gathers for slope estimation can be obtained from the results of prior demultiple processes using, for example, Radon-based filtering of multiples. These demultipled gathers are the same inputs that are generally used in prestack time/depth migrations.

Horizontally layered model
We first consider a horizontally layered model shown in Figure 5 and look at the Green's function from (0 m, 1000 m) on the third reflector. The input CMP corresponding to this position is shown in Figure 6a. Its estimated dip and the associated one-way traveltime map t 0 ðx 0 ; h; tÞ generated according to the proposed workflow are shown in Figure 6b and 6c, respectively. The contour of  Velocity-independent Marchenko Q61 phase wavelet and multiplied by the weight to properly account for geometric spreading. For simplicity, we express this weight along the desired traveltime contour as a panel (Figure 6d). The weighting can then be easily accomplished by a point-wise multiplication to the panel of the convolved traveltime contour. This weight is determined from the ratio 1∕ ffiffi ffi r p in a 2D medium (or 1∕r in 3D), where r is the distance from the specified point to different positions on the surface in the time-imaging domain. The use of this ratio is permissible due to the straight-ray and effective velocity assumptions used in the time-imaging domain.
The true Green's function from forward modeling is shown in Figure 7a overlain by the traveltime prediction (the dotted-dashed magenta line) of the direct wave using the proposed slope-based workflow. We use this prediction to generate the initial focusing function and solve the Marchenko system (equation 10). The result from the usual velocity-based workflow (van der Neut et al., 2015), with the initial focusing function generated by numerical wave propagation in a background velocity model, is shown in Figure 7b in comparison with that from the proposed slope-based method in Figure 7c. A comparison between the central traces of the true and estimated Green's functions from Figure 7 is shown in Figure 8. We observe that the result from the slope-based workflow is slightly inferior due to the assumption of zero-phase wavelet, which leads to small errors in the predicted coda. That aside, the results are comparable in quality indicating the validity of the proposed method. We also investigate another potential choice of wavelet that corresponds to the characteristic phase behavior of 2D far-field Green's functions in the "Discussion" section.

Laterally heterogeneous model
Next, we turn to a laterally heterogeneous model with lateral heterogeneity (Figure 9). In this example, the image rays are no longer vertical and the focusing positions in the time-and depthimaging domains are related through the mapping defined by image rays. We consider the Green's function from (−35 m, 1000 m) on the third reflector because this is the position at which the image ray originating from (0 m, 0 m) will pass through. Figure 10 shows the inputs and intermediate products of the proposed slope-based workflow. Note that there is a nonzero dip in the midpoint direction (Figure 10c). In this example, the contour of t 0 ¼ 0.4165 s is chosen. The final comparison of Green's functions is shown in Figure 11, and a comparison of central traces is shown in Figure 12  again are comparable in quality, and similar conclusions can be drawn. We emphasize that the focusing position is defined in terms of (x 0 , t 0 ) as opposed to (x, z) in the usual Marchenko workflow. To confirm that the specified focusing (x 0 , t 0 ) translate to the Cartesian (x, z) according to the image ray mapping, we back propagate the computed focusing function, and the result is shown in Figure 13. The dashed line is vertical, whereas the solid line is the image ray originating from (0 m, 0 m). We can observe that the response is now positioned along the image ray and is at the reflector that we chose to compute our slope for the traveltime prediction in the first place.

Field-data example from the North Sea
For the last example, we test our proposed approach against an ocean-bottom cable data set from the North Sea. The reference velocity model in Cartesian coordinates of this data set is shown in Figure 14. The receiver line placed at the sea bottom of 92 m depth spans a length of 6 km (from 3147 to 9098 m) and contains 235 sensors with an approximate interval of 25 m. The shot interval is 50 m along the 12 km sail line across the extent of the model. We specify two focusing points at a shallow position (6000 m, 1306 m) and at a deep position (6000 m, 2167 m), where we will compare the results of our approach with those from the velocity-based workflow.
To facilitate the comparison between results from our proposed approach and those from the previous velocity-based workflow, we first convert the reference velocity model (Figure 14) from depth to time according to a 1D-medium assumption. This process is approximately valid in this case due to the nearly lateral homogeneous characteristic of the model close to the desired focusing positions around the midpoint 6000 m. A comparison between the velocity models before and after this process is shown in Figure 15. We will use the velocity in Figure 15b to obtain a CMP gather for our slopebased workflow that in principle, should allow a direct comparison between our results and those from the velocity-based workflow in depth.
An example raw CMP gather of this data set at the location 6000 m is shown in Figure 16. To obtain an accurate estimation of slopes that corresponds to primary reflections, we apply a demultiple process with the use of a hyperbolic Radon transform based on the converted reference velocity in Figure 15b. This process amounts to transforming the input raw CMP gather using the hyperbolic Radon transform with given reference velocity, selective windowing of nearly flat events, and inverse transforming these selected data. This requirement of primary-only gathers in our pro-  Figure 11. A comparison in the laterally heterogeneous model of (a) the true Green's function, (b) estimated ones using the true velocity to generate initial focusing functions, and (c) that from the proposed slope-based workflow. The dotteddashed magenta line in (a) denotes the estimated traveltime of the direct wave using local slopes.

a) b) c)
Velocity-independent Marchenko Q63 posed approach substitutes the previous need of a depth velocity model in the conventional Marchenko workflow. As described in the previous section, the slope of primary-only gathers is used to obtain the initial focusing function, but the reflection operator R in equation 10 still makes use of the full data set containing all or-ders of multiples. We will discuss more about this point in the next section. The resulting demultipled gather at a similar location is shown in Figure 17 together with the estimated slope and chosen t 0 contours at 0.6 and 1.0 s associated with the two focusing positions. Figure 18 shows the direct-wave Green's functions obtained from the conventional workflow with wave extrapolation in the reference velocity model (Figure 14) overlain by the traveltime prediction (the dotted-dashed magenta line) using the proposed slopebased workflow. We can observe a good agreement between the results from both approaches, which corroborates our proposition. In our experiment, we observe that the processing step on hyperbolic regression plays a crucial role in controlling the final quality of the traveltime predictions. The number of samples of picked traveltime contours used for the regression must be chosen in accordance with the offset range and the quality of measured slopes at a large offset. Here, we choose to run a regression in the offset range of −350 to 350 and −600 to 600 for the shallow and deep focusing positions, respectively. Using the obtained traveltime predictions, we follow the proposed workflow and obtain the final estimates of f þ    for the shallow focusing position shown in Figure 19. Their counterparts from the conventional velocity-based workflow are shown in Figure 20. The results appear comparable in quality to similar prominent events present and no apparent artifacts. To further compare the results, we investigate the difference between Figures 19 and 20 as shown in Figure 21. We can observe that the only difference in estimated f þ 1 and g þ corresponds to the estimate of the direct-wave Green's function. This is not surprising because we arbitrarily choose a zero-phase wavelet as an approximate waveform in our workflow. However, the estimated f − 1 and g − appear quite different between the results from the two methods. We emphasize that this observation does not invalidate our slope-based workflow because in the Marchenko system (equation 10), as discussed by van der Neut et al. (2015), all of the physical events in the upgoing Green's function g − are retrieved with the correct amplitudes at the first iteration by the action of the reflection operator R on the initial focusing function (the time-reversed direct-wave Green's function), whereas later iterations amount to eliminating artifacts. However, physical events in the downgoing Green's function g þ are updated throughout all the iterations. Therefore, the fact that we do not observe any difference in the g þ panel apart from the direct-wave Green's function, but do observe some small difference throughout the g − panelwhile using the same R in both experimentsis due entirely to the difference in our choice of the direct-wave Green's function (i.e., the zero-phase waveform). The absence of any prominent artifacts between both results also suggests that our Marchenko solver is functioning properly. This find- ing is also supported by a comparison of the central traces of the estimated f − 1 and g − as shown in Figure 22. There, we observe that the coda events captured by the velocity-based and the proposed workflows are similar, except for minor amplitude differences and slight shifts due to the use of different direct-wave Green's functions. These relatively small mismatches lead to the errors shown in the difference panels in Figure 21b and 21d. Similar results for the deep focusing position are shown in Figures 23, 24, 25, and 26, from which the same conclusions can be drawn.
Finally, it is important to emphasize that in practice, we do not have access to such a reference velocity model ( Figure 14) and its corresponding velocity in time (Figures 15b and 27a) at the early stage of processing. Given raw CMP gathers, we would have to follow the conventional preprocessing procedure (e.g., surface-consistent correction, deghosting, and demultiple) and obtain primaryonly gathers before implementing the proposed workflow. These processes may also involve a choice of background velocity trend that requires quality control. Therefore, a more realistic use of this data set can be achieved by first picking a background velocity trend from the velocity analysis (Figure 27b). We subsequently use this picked velocity in our simple demultiple process with the hyperbolic Radon transform. Figure 28 shows the resulting gather together with the estimated slope and chosen t 0 contours associated with the two focusing positions. The final traveltime predictions after hyperbolic regression are shown in Figure 29 indicating a good performance. In this example, these predictions lead to almost similar results to those in Figures 19 and 23.

DISCUSSION
We emphasize first that using the proposed method, a prior knowledge on a smooth velocity model is not required and a forward wavefield extrapolation is no longer needed to generate the direct-wave Green's function (the time-reversed initial focusing function). By using the knowledge of local slopes, we show how to obtain the direct-wave Green's function directly from the primary-only CMP gathers under the usual assumptions of time-domain imaging. The original recorded gathers are still required to construct the reflection operator R used in the Marchenko system.
In view of time-imaging assumptions, we emphasize that in the presence of strong lateral heterogeneity, image rays that define the time-imaging map may cross-cut each other and lead to caustics. In this study, we restrict ourselves to cases in which the lateral heterogeneity is mild, and the image rays are uniquely defined, with no complications from caustics. Future investigations are required to properly formulate the Marchenko equations appropriate for time-domain imaging under such circumstances. Moreover, we rely on the hyperbolic traveltime approximation, which is exact for plane reflectors in a constant-velocity (homogeneous isotropic) or elliptically anisotropic overburden and approximately valid for small-offset data in other cases. For example, we can already Velocity-independent Marchenko Q67 observe the decrease in accuracy in Figure 11a, where the estimated traveltime denoted by the dotted-dashed magenta line deviates slightly at larger distances. In such cases, other (more accurate) traveltime approximations can, in principle, be used and additional measurements (more parameters) in addition to local slopes may be needed (Khoshnavaz et al., 2016;Stovas and Fomel, 2016). The topic of traveltime approximation is one of the primary research subjects in time-domain imaging, and there are many choices of traveltime approximation with different ranges of accuracy and applicability (Tsvankin and Thomsen, 1994;Alkhalifah and Tsvankin, 1995;Gelchinsky et al., 1999;Jäger et al., 2001;Fomel and Stovas, 2010;Santos et al., 2011;Blias, 2013;Dell et al., 2013;Farra et al., 2016;Ravve and Koren, 2017;Sripanich et al., 2017).
Apart from the limitations that come with those time-imaging assumptions, we observe that the choice of wavelet also plays a notable role in controlling the final quality of the results from the Marchenko method and in ensuring consistent results from both setups in the time-and depth-imaging domains. In the current approach, the amplitude information is only partially handled by weighting according to the geometric spreading factor corresponding to the straight-ray assumption from the focusing position.
An alternative is to also consider the behavior of the far-field Green's function in a 2D constant-velocity medium. This involves an additional 45°phase shift and a division by ffiffiffi ffi ω p . The results from our proposed slope-based workflow with this choice of wavelet in the two synthetic models are shown in Figure 30. Similar to the case of a zero-phase wavelet, we can still observe minor differences between the exact Green's functions and the estimated ones. Finding a dynamically appropriate waveform with nonzero phase from slopes using, for example, wavelet estimation methods and far-field signatures is the subject of future research. In a recent notable example, the augmented Marchenko approach by Dukalski et al. (2019) relies on energy conservation for the focusing fields to yield amplitude and phase corrections to the initial focusing functionsthe approach currently applies only to laterally invariant media, but the framework could probably be extended to more complex scenarios. In our case, we note that our current choice of a zero-phase wavelet will still lead to a focus at the desired position, but the focused response will contain a slightly different amplitude signature than that computed from the velocity-based workflow, which intrinsically involves a nonzero-phase wavelet. Nonetheless, as shown by our experiments, the Marchenko method can still provide appropriate coda waveforms for focusing functions and Green's functions while accommodating various wavelet choices for the initial focusing functions used as an input to the system. In a companion paper (Sripanich and Vasconcelos, 2019), we discuss additional insights on Marchenko focusing, by considering the method in the time-imaging domain and analyzing the amplitude and directionality of the focused responses with respect to the surface data aperture.
Accurate knowledge of local slopes also plays a crucial role in the traveltime approximation of the direct-wave Green's functions. In this study, we use the dip estimator based on the plane-wave destruction filter described by Fomel (2002), which requires the computational cost of OðN d N f N i Þ, where N d is the number of data points, N f is the length of the filter, and N i is the total number of a) b) that is approximately laterally homogeneous (1D) and (b) picked velocity from velocity analysis, which is done following the conventional processing procedure.
Velocity-independent Marchenko Q69 iterations. Several other methods exist for this purpose (Harlan et al., 1984;Marfurt, 2006;Hale, 2007;Schleicher et al., 2009;Chen et al., 2013) and should be chosen on a case-by-case basis depending on the desired level of accuracy and efficiency. Moreover, in our scheme, it is crucial that only the slopes of target primary reflection events are used for traveltime prediction. This particular requirement is similar to most processing/imaging techniques that work with primary reflections. In view of the Marchenko method, this replaces the need of an approximate depth velocity model, which is relatively much harder to obtain and requires costly tomographic or full-waveform inversion. Provided the primary-only gathers, our method can readily provide estimates of the initial focusing functions based on data-driven slope measurements without directly computing for subsurface velocity. However, if the gathers still contain residual multiples, the traveltime predictions in equations 11 and 12 are expected to be erroneous because both equations were derived using the geometry of primary reflections (Fomel, 2007). For example, in our horizontally layered model, if the raw gather with primaries and multiples is used in the slope estimation process, the results would look like those shown in Figure 31. In comparison to those in Figure 6, we observe that the local slopes of primaries and multiples are estimated, leading to inaccurate traveltime estimates. One should expect that such errors will be higher in more complex models, in which primaries and multiples may overlap and lead to conflicting dips. An expeditious multiple removal based on local slopes (Figure 31b) using simple velocity filtering exists (Cooke et al., 2009) and may be incorporated in the workflow.
We note that the recent development of the plane-wave Marchenko method (Meles et al., 2018) can also potentially benefit from a consideration based on the time-imaging domain. In this scheme, a different focusing condition is proposed to solve for which represents a summation of conventional focusing functions along some focusing surface S f at depth z f . Using the one-way reciprocity theorem for curvilinear coordinates (equations 1 and 2) and considering the time-imaging domain, one can follow a similar procedure and which represents a summation of focusing functions along some focusing surface S f (image wavefront) at vertical time t f . The use of time-imaging tools (such as local slopes) in this context requires further investigation, but we note that the focusing position defined by t f in this case, again, does not require the knowledge of velocity as t f is directly related to the coordinates of the recorded CMP data. Finally, we highlight that this work can also be deemed complementary to that of van der Neut and Wapenaar (2016). In their study, it was shown that the single-sided Green's function representations in the Cartesian coordinates could be modified by redatuming with a choice of direct-wave Green's function related to some fictitious reflector. By doing so, revised single-sided representations can be developed, and the initial focusing function (strictly the inverse of the direct-wave Green's function) for the corresponding Marchenko method becomes a band-limited delta function. The direct-wave Green's function that was used to write these revised representations in the first place does not need to be specified explicitly, conveniently obviating the requirement for a smooth velocity model. Van der Neut and Wapenaar (2016) subsequently show that in this redatumed domain, it is possible to separate and eliminate the effects of primary and multiple reflections from the overburden above the chosen fictitious reflector, and an effective internal multiple elimination technique can be developed. However, this regime is not applicable to Green's function retrieval because the knowledge of the direct-wave Green's function is never explicitly computed. On the other hand, in this study, we approach the Marchenko method from the perspective of time-domain imaging and show that it is possible to obtain an approximate direct-wave Green's function straightforwardly from primary-only CMP gathers using local slopes without needing any velocity model. Therefore, it is conceptually possible to adopt the method of van der Neut and Wapenaar (2016) or its variant (Zhang and Staring, 2018) to obtain primaryonly gathers needed by our proposed approach, which then leads directly into our velocity-independent Marchenko-based focusing and redatuming workflow.

CONCLUSION
In this paper, we formulate a new form of the Marchenko system in the time-imaging domain defined by image rays. We show that the resulting Marchenko equations have the same form as in the previous case of a constant-depth datum level, except that the focal-point coordinates and radiation directions are now determined along curved space boundaries of constant time (image-wavefront surfaces) that define the time-imaging domain. A priori knowledge of subsurface velocity is no longer needed to generate the initial focusing function and the Marchenko scheme can be carried out by making use of local event slopes measured directly from recorded CMP gathers. The focusing positions are now defined in terms of image-ray coordinates related to the surface location and vertical traveltime, easily obtained from the data without any prior knowledge of the subsurface model. The resulting focusing functions and Green's functions obtained using the proposed method are of comparable quality to those conventionally obtained with prior knowledge of the subsurface velocity.

DATA AND MATERIALS AVAILABILITY
Data associated with this research are available and can be obtained by contacting the corresponding author. a) b) c) d) Figure 31. Inputs and products from the slope-based workflow in the horizontally layered model: (a) raw CMP gather that contains both the primaries and multiples, (b) estimated dip with a plane-wave destruction filter, (c) one-way traveltime map t 0 ðx 0 ; h; tÞ, and (d) the difference between the estimated one-way time in the presence of multiple (c) and that from the primary-only gather in Figure 6c. We can observe an error in traveltime predictions when the gather does not contain only primary reflections.