Data-driven retrieval of primary plane-wave responses

Seismic images provided by reverse time migration can be contaminated by artefacts associated with the migration of multiples. Multiples can corrupt seismic images, producing both false positives, i.e. by focusing energy at unphysical interfaces, and false negatives, i.e. by destructively interfering with primaries. Multiple prediction / primary synthesis methods are usually designed to operate on point source gathers, and can therefore be computationally demanding when large problems are considered. A computationally attractive scheme that operates on plane-wave datasets is derived by adapting a data-driven point source gathers method, based on convolutions and cross-correlations of the reflection response with itself, to include plane-wave concepts. As a result, the presented algorithm allows fully data-driven synthesis of primary reflections associated with plane-wave source responses. Once primary plane-wave responses are estimated, they are used for multiple-free imaging via plane-wave reverse time migration. Numerical tests of increasing complexity demonstrate the potential of the proposed algorithm to produce multiple-free images from only a small number of plane-wave datasets.


Introduction
Most standard processing steps, e.g. velocity analysis (Yilmaz (2001)) and reverse time migration (Whitmore (1983); McMechan (1983); Zhu et al. (1998); Gray et al. (2001); Mulder and Plessix (2004)), are based on linear (Born) approximations, for which multiply scattered waves represent a source of coherent noise. When linearized methods are employed, multiples should then be suppressed to avoid concomitant artefacts. Free-surface multiples particularly affect seismic images resulting from marine data (Wiggins (1988)), and many algorithms have been designed to attenuate the presence of free-surface multiples (for a comprehensive review see Dragoset et al. (2010)). On the other hand, internal multiples strongly contaminate both land (Kelamis et al. (2006)) and marine data (van Borselen (2002)). Fewer techniques have been designed to estimate and remove internal multiples. The seminal method by Jakubowicz 1 (1998) uses combinations of three observed reflections to predict and remove internal multiples. However, this scheme requires prior information about reflections to allow proper multiple prediction and removal. On the other hand, applications of inverse scattering methods (Weglein et al. (1997)) can predict all orders of internal multiple reflections with approximate amplitudes in one step without model information (ten Kroode (2002); Löer et al. (2016); Zhang et al. (2019a)).
Multiple-related artefacts can also be dealt with via Marchenko methods. Marchenko redatuming estimates Green's functions between arbitrary locations inside a medium and real receivers located at the surface ; Wapenaar et al. (2012Wapenaar et al. ( , 2014; da Costa Filho et al. (2014)). In Marchenko redatuming, Green's functions are estimated using reciprocity theorems involving so called 'focusing functions', i.e. wavefields which achieve focusing properties in the subsurface ). In contrast to seismic interferometry, Marchenko redatuming requires an estimate of the direct wave from the virtual sources to the surface receivers, only one sided illumination of the medium and no physical receivers at the position of the virtual sources ; Wapenaar et al. (2014)). Focusing functions and redatumed Green's functions can provide multiple-free images directly ; Wapenaar et al. (2014)). Moreover, combining Marchenko methods and convolutional interferometry allows estimating internal multiples in the data at the surface (Meles et al. (2015); da Costa Filho et al. (2017b)). Other applications of the Marchenko method include microseismic source localization (Behura and Snieder (2013); van der Neut et al. (2017); Brackenhoff et al. (2019)), inversion ; van der Neut and Fokkema (2018)), homogeneous Green's functions retrieval (Reinicke and Wapenaar (2019); Wapenaar et al. (2018)) and various wavefield focusing techniques (Meles et al. (2019)). Despite its requirements on the quality of the reflection data, and more specifically its frequency content, the Marchenko scheme has already been successfully applied to a number of field datasets (Ravasi et al. (2016); van der Neut et al. (2015b); Jia et al. (2018);da Costa Filho et al. (2017a); Staring et al. (2018); Zhang and Slob (2020b)). Further developments have also shown how a successful Marchenko redatuming can be achieved either via correct deconvolution of the source wavelet from the measured data or by including wavelet information in the Marchenko equations (Ravasi (2017); Slob and Wapenaar (2017); Becker et al. (2018)). Recent advances in Marchenko methods led to revised derivations which resulted in fully data driven demultiple / primary synthesis algorithms (van der Neut and Wapenaar (2016); Zhang et al. (2019b); Zhang and Slob (2019). Different from standard Marchenko applications, in these revised derivations the focusing functions are projected to the surface, thus leading to the retrieval of specific properties of reflections responses in the data at the surface (i.e., internal multiples/primaries) instead of redatumed Green's functions. We refer to the class of applications introduced by van der Neut and Wapenaar (2016) and Zhang et al. (2019b) as to 'data domain Marchenko methods'.
Inspired by work on areal-source methods for primaries (Rietveld et al. (1992)), Marchenko redatuming and imaging schemes were recently adapted to include plane-wave concepts (Meles et al. (2018)). Here, we follow a similar approach and extend the applications of data domain Marchenko methods, originally derived for point sources, to plane-wave sources. The benefit of using plane-wave data for imaging, i.e. an overall reduction in the data volume and the possibility to get subsurface images by migrating fewer plane-wave gathers than shot gathers (Dai and Schuster (2013); Wang et al. (2018); Stoffa et al. (2006); Schultz and Claerbout (1978)) is then combined with a fully data-driven demultiple scheme.

Data domain Marchenko method
In this section we briefly summarize the primary reflections retrieval algorithm recently proposed by Zhang et al. (2019b) and in sections 3.2 and 3.3 discuss how it can be extended to include plane-wave concepts. First, we briefly introduce the definitions and properties of the so-called Marchenko focusing functions, upon which the work on projected focusing functions is based. Following standard notation, we indicate time as t and the position vector as x = (x H , z), where z stands for depth and x H for the horizontal coordinates (x, y). An acoustically transparent acquisition boundary ∂D 0 is defined at z 0 = 0 and points in ∂D 0 are denoted as x 0 = (x H , z 0 ). Similarly, points along an arbitrary horizontal depth level ∂D i are indicated as x i = (x F , z i ), where z i indicates the depth of ∂D i and x F denotes the horizontal coordinates of a focal point at this depth. Note that boundaries ∂D 0 and ∂D i in 2D and 3D are lines and planes, respectively (for a comprehensive analysis of generalized Marchenko concepts in 2D and 3D, see Wapenaar et al. (2018)). The focusing function f 1 (x 0 , x i , t) is the solution of the source-free wave equation in a truncated medium which focuses at the focal point x i . We define the truncated medium as being identical to the physical medium between ∂D i and ∂D 0 , and reflection-free elsewhere ). The focusing function f 1 (x 0 , x i , t) is decomposed into down-and up-going components, indicated by f + 1 (x 0 , x i , t) and f − 1 (x 0 , x i , t), respectively. The down-going component of the focusing function, f + 1 (x 0 , x i , t), is the inverse of the transmission response T (x i , x 0 , t ) of the above mentioned truncated medium, i.e. where and T (x i , x 0 , t ) can be decomposed into direct and coda components, indicated by d and m subscripts, respectively: and Using source-receiver reciprocity, Eq. 1 can be generalized as: where δ(x H −x H ) is now a two-dimensional delta function along ∂D 0 . The up-going component of the focusing function, f − 1 (x 0 , x i , t), is by definition the reflection response of the truncated medium to f + 1 (x 0 , x i ), and it is equivalent to: where R(x 0 , x 0 , t) is the impulse reflection response (with the source ignited at time t = 0 to allow standard Marchenko derivations) at the surface of the physical medium, with x 0 , x 0 denoting receiver/source locations. This relationship is valid for −t d + ε < t < t d + ε, where t d is the one-way traveltime from a surface point x 0 to x i and ε is a small positive value accounting for the finite bandwidth of the data. Note that, unlike for the original Marchenko scheme, we have chosen an asymmetric time interval, following Zhang et al. (2019b). For this time interval, the coda of the down-going focusing function, namely f + 1m (x 0 , x i , t), satisfies the following relationship: Next we project the focusing functions to the surface. The projected focusing functions v − and v + m are then introduced as: and where the variable z i indicates that these functions depend on the depth level along which standard Marchenko focusing functions are defined. Note that differently than in previous literature (van der Neut and Wapenaar (2016); Zhang et al. (2019b)) we now make explicit the dependence of v − and v + m on z i (Zhang and Slob (2020a)). By convolving and integrating in space along ∂D i both sides of Eqs. 5 and 6 with T d as indicated in Eq. 4, we obtain: and for ε < t < t 2 + ε, where for convenience we have replaced the dependence on z i by the new variable t 2 = t 2 (x 0 , x 0 , z i ) corresponding to the two-way traveltime from a surface point x 0 to the specular reflection at a (hypothetical) interface at level z i and back to the surface point x 0 . Different from previous literature on this subject, we make all the relevant variables in v − and v + m explicit, by considering also t 2 . Note that for t < ε and t > t 2 + ε both v − and v + m are zero, which is why the integrals on the right-hand side are evaluated only for the time interval < t < t 2 + . Using the time-domain formalism introduced in van der Neut et al. (2015a) we rewrite Eqs. 9 and 10 as: and where R indicates a convolution integral operator of the measured data R with any wavefield, the superscript indicates time-reversal and Θ t 2 +ε ε is a muting operator removing values outside of the interval (ε, t 2 + ε).
Terms in Eq. 11 are rearranged using Eq. 12 to get: which, under standard convergence conditions (Fokkema and van den Berg (1993)), is solved by: This procedure allows to retrieve v − (x 0 , x 0 , t, t 2 ), whose last event, when its two-way travel time t is equal to t 2 (x 0 , x 0 , z i ) is a transmission loss compensated primary reflection in R(x 0 , x 0 , t) (Zhang et al. (2019b)). In practice, the transmission loss compensated primary is obtained by computing v − via Eq. 14 for all values t 2 (i.e., by considering the corresponding windowing operator Θ t 2 +ε ε ), and by storing results in a new, parallel dataset at t = t 2 . Similarly to other Marchenko schemes, in practical applications only a few terms of the series in Eq. 14 need to be computed to achieve proper convergence ). Moreover, following Zhang and Staring (2018), instead of computing t 2 as the space-and model-dependent two-way traveltime via a chosen depth level z i , we can evaluate Eq. 14 for all possible constant values t 2 (to include values large enough to allow waves to reach the bottom of the model and come back to the surface) and store results at t =t 2 . In this way the (transmission-compensated) primary reflection response in R(x 0 , x 0 , t) is then fully retrieved.

Extension to horizontal plane-wave data
In this paper, following a similar approach to what was recently proposed to extend Marchenko redatuming from point-source to horizontal plane-wave concepts (Meles et al. (2018)), we consider integral representations of the projected focusing functions v − and v + m . More precisely, we first define new projected focusing functions V − (x 0 , t, t 2 ) and V + m (x 0 , t, t 2 ) as: where T 2 = T 2 (x 0 , z i ) is the two-way traveltime of a horizontal plane-wave propagating down from the surface to the specular reflection at a (hypothetical) interface at level z i and back to the surface point x 0 . We then integrate Eqs. 9 and 10 along ∂D 0 to obtain: and for ε < t < T 2 + ε and where R PW (x 0 , t) ≡ ∂D 0 dx 0 R(x 0 , x 0 , t) is by definition the horizontal plane-wave source response of the medium (i.e., the source emits a vertically downward propagating plane wave). Using again the time-domain formalism we can therefore rewrite Eqs. 17 and 18 as: and and therefore: which is solved by: This procedure allows to retrieve V − (x 0 , t, T 2 ), whose last event, when its two-way travel time t is equal to T 2 (x 0 , z i ), is a transmission loss compensated primary reflection in R PW (x 0 , t).
Instead of computing T 2 as the space-and model-dependent two-way traveltime via a chosen depth level z i , we can evaluate Eq. 22 for constant valuesT 2 . By computing Eq. 22 for all possible constant valuesT 2 and storing results at t =T 2 , the (transmission-compensated) primary reflection response in R PW (x 0 , t) is then fully retrieved. Note that in practical applications, the integrals along ∂D 0 in Eqs. 15-18 and in the definition of R P W are replaced by summations over source locations.

Extension to dipping plane-wave data
In standard Marchenko derivations it is assumed that point sources are fired at t = 0 ; Zhang et al. (2019b)). Since dipping plane-waves are associated with many sources excited at different times, we cannot expect standard algorithms, such as that in Eq. 22, to 6 predict primaries when delayed source gathers are considered. To illustrate how to proceed when dipping plane-waves are taken into account, we first consider the obvious corresponding projected focusing functions: and where p is a ray parameter vector and T 2 = T 2 (x 0 , p, z i ) is the two-way traveltime of a planewave with ray parameter p, propagating down from the surface to the specular reflection at a (hypothetical) interface at level z i and back to the surface point x 0 . Substituting Eqs. 9 and 10 into Eqs. 23 and 24, and indicating the reflection response associated with a dipping plane-wave source characterized by ray parameter vector p as R DW (x 0 , p, t) ≡ ∂D 0 dx 0 R(x 0 , x 0 , t−p·x H ), we obtain: and for ε + p · x H < t < T 2 + ε. The relationship between V − (x 0 , p, t, T 2 ) and V + m (x 0 , p, t, T 2 ), using again the time-domain formalism, is then established by: and Combining Eqs. 27 and 28 together we finally get: which is solved by: (30) This procedure allows to retrieve V − (x 0 , p, t, T 2 ), whose last event, when its two-way travel time t is equal to T 2 (x 0 , p, z i ), is a transmission loss compensated primary reflection in R DW (x 0 , p, t). Note that, in principle, the muting operators in Eq. 30, similarly to those in Eqs. 14 and 22, are space-and model-dependent. However, in analogy to the previous cases, the upper boundary of the muting operators in Eq. 30 can be taken parallel to the lower one (see Fig. 1), thus exhibiting a space-dependent but model-independent shape, i.e. T 2 (x 0 , p, z i )+ε ≈ ε+T 2 +p·x H for a generic constant valueT 2 . By computing Eq. 30 for all possible constant valuesT 2 and storing results at t =T 2 + p · x H , the (transmission-compensated) primary reflection response in R DW (x 0 , p, t) is then fully retrieved. The performance of the algorithm in Eq. 30 is assessed in the following numerical examples.

Numerical Examples
We explore the potential of the proposed scheme for the retrieval of plane-wave source primary reflections with numerical examples involving increasingly complex 2D models. Evaluation of the series in Eq. 22 requires computation of the operators R and R and of the plane-wave reflection response R PW (x 0 , t). The reflection responses in R and R need to be recorded with wide band and properly sampled (according to Nyquist criterion in space and time) co-located sources and receivers placed at the surface of the model. In the following numerical examples, source-receiver sampling is set to 10m, while gathers R PW (x 0 , t) are computed with a 20 Hz Ricker Wavelet. All data used here are simulated with a Finite Difference Time Domain solver (Thorbecke et al. (2017)).
For our first numerical experiment we consider a 2D model with gently dipping interfaces 8 Figure 2: (a) Velocity and (b) density models used in the first numerical experiment.  (see Fig. 2). The recording surface is reflection free. The dataset associated with a horizontal plane-wave source fired at the surface of this model is shown in Fig. 3(a). Notwithstanding the geometrical simplicity of the model, due to the strong impedance variations the data are contaminated with many internal multiples, as indicated by the red arrows. We then apply to this dataset the method as described in section 3.2. More precisely, we compute V − via Eq. 22 for all valuesT 2 , and by storing results at t =T 2 we build a parallel dataset, which theoretically only involves primaries. Note that the algorithm is fully data driven, and no model information or any human intervention (e.g., picking) is involved in the process. For this dataset we only computed the first 20 terms of the series in Eq. 22. The result of the procedure is shown in Fig.  3(b). We then image both datasets in Fig. 3 via standard plane-wave reverse time migration (based on the zero lag of the cross-correlation between the source and receiver wavefields, Claerbout (1985)) using a smoothed version of the true velocity distribution in Fig. 2 and constant density. Migration results are shown in Fig. 4. When the full dataset is migrated, internal multiples contaminate the image as shown in Fig. 4(a), producing many false positive artefacts (indicated by red arrows). The image is much cleaner when the dataset in 3(b) is migrated. Each interface is properly recovered, as demonstrated by a comparison between Figs. 2 and 4(b). Green arrows in 4(b) point at physical interfaces which are invisible in Fig.  4(a), where they are attenuated by interfering multiple-related artefacts. Black arrows point at physical interfaces only partially resolved. The relatively poor performances in imaging dipping interfaces is not due to residual internal multiples, but to the intrinsic limitations of horizontal plane-wave imaging. However, note that only one demultipled plane-wave response and a single migration were required to produce the multiple-free image in Fig. 4(b). We conclude that for gently dipping models horizontal plane-wave datasets are sufficient to produce satisfactory results.
In the second example (Fig. 5) we consider a more challenging model with critical features for any Marchenko method, i.e. the presence of thin layers, diffractors and dipping layers ; Zhang et al. (2019b); Dukalski et al. (2019)). We initially follow the same imaging strategy as for the first example. We first compute the dataset associated with a horizontal plane-wave source fired at the surface of the model shown in Fig. 5(a). Also for this dataset we only computed the first 20 terms of the series in Eq. 22. Given the complexity of the model, many events, primaries as well as internal multiples (red arrows) cross each-other, especially in the lower part of the plane-wave gather. Picking specific events in the gather in Fig. 6(a) would be challenging. However, as discussed above, our method does not involve any human intervention, and by applying the same scheme as for the first model we retrieve the dataset shown in Fig. 6(b), where primaries otherwise overshadowed by interfering multiples are clearly visible (green arrows). We then migrate datasets in Fig. 6(a) and (b), and show in Fig. 7(a) and (b) the corresponding images. Large portions of the image in Fig. 7(a) associated with the dataset in Fig. 6(a) are dominated by noise due to the presence of internal multiples (red arrows). On the other hand, the image in Fig. 7(b), which is associated with the estimated primaries in Fig. 6(d), is much cleaner, with fewer artefacts (red arrows) contaminating limited domains of the image. Note that relatively poor imaging performances of dipping interfaces (black arrows in Fig. 7(b)) are not necessarily associated with shortcomings of the discussed demultiple method but with the intrinsic limitation of what can be illuminated by a single plane-wave experiment. For this specific model we then decide to process and migrate also dipping plane-wave data. In total we then consider 10 additional datasets, uniformly ranging from −25 • to 25 • (as discussed in section 3.3, the angle of the plane-wave is implemented by adding time delays to the shot positions on the horizontal array). Representative dipping planewave data are shown in Fig. 6(b,c), next to the corresponding processed gathers (in Fig. 6(e,f)). Red and green arrows point again at internal multiples and recovered primaries, respectively. We finally consider aggregate plane-wave migrated images. By migrating a total of 11 fulldata gathers, the image in Fig. 7(c) is obtained. While thanks to the better illumination the improvement over the image in Fig. 7(a) is clear, some of the key features of the final result are still misleading (red arrows point at false positives associated with the migration of internal multiples). A significantly better result is obtained when the 11 processed gathers are imaged and stacked (Fig. 7(d)). The dipping features poorly visible in (b) are now properly resolved. This example shows that the proposed method can successfully process dipping planewave datasets and therefore benefit from the corresponding improved illumination. Residual artefacts in the migrated image indicated by the red arrow in Fig. 7(d) are likely due to the presence of thin layers, diffractors and dipping layers that are known to be critical in Marchenko methods.   Fig. 6(a). Red arrows point at artefacts related to internal multiples. (b) Standard plane-wave reverse time migration of the dataset in Fig. 6(b). Black arrows indicate dipping interfaces that are only partially recovered due to the poor illumination provided by a single plane-wave experiment. Note that these interfaces are also not properly imaged in (a). (c) Aggregate RTM of 11 plane wave full datasets (uniformly ranging from −25 • to 25 • ). Red arrows point at artefacts related to internal multiples. (d) Aggregate Standard RTM of synthesized primaries. Green and red arrows indicate interfaces barely visible in (a) and minor residual artefacts, respectively. Differences in amplitude between images in (a,c) and (b,d) are due to multiples removal and transmission loss compensation (see Fig. 6). 14

Discussion
In section 3.2 we have extended a recently proposed primary synthesis method, devised for point source gathers, to horizontal plane-wave source data. The new scheme still needs full point-source data as input, but its output is a horizontal plane-wave response. The method is based on integration of point-source responses over the acquisition surface (e.g., Eqs. 9 and 10), which allows the derivation of relationships associated with plane-wave sources (e.g., Eqs. 17 and 18). Both the point-source and plane-wave primary synthesis methods are totally data driven, and both are implemented by inversion of the same family of linear operators, i.e.: Each operator is defined by a different value ofT 2 . In previous literature that underlies this contribution, an integration over the focusing surface was used to adapt Greens' functions redatuming methods to virtual plane-wave redatuming (Meles et al. (2018)). While conceptually similar, there is a subtle yet very important difference between the methods discussed here and previous methods on virtual plane-wave redatuming. Whereas in any Marchenko redatuming scheme (e.g. for point or plane virtual sources) a different, model dependent, window operator for each point or plane is required, as focusing is achieved in the subsurface, the window operators discussed here are the same for each input data, as the focusing operators are projected to the surface. Since the operators in Eq. 31 are linear and do not depend on the specific gather they are applied to, any linear combination of point-source data can be processed at once, provided that all the corresponding sources are fired at the same time (see section 3.2 for more details). The proposed method can then be used, without any modification, to blended-source data as well as to individual point sources and horizontal plane-wave gathers. This is shown in Fig. 8, where the algorithm is applied to a dataset associated with 5 sources with different spectra fired at the same time ( Fig. 8(a)). Application of the proposed scheme results in the gather shown in Fig. 8(b). A nearly identical result (relative difference smaller than 0.1%) is achieved when the method is applied to each single point source gather separately, after which the corresponding results are summed together. In section 3.3 we extended the primary synthesis method for dipping plane-wave source data, which helps to improve the illumination of dipping interfaces in the subsurface.

Conclusions
We have shown that recent advances in data domain Marchenko methods can be extended to incorporate plane-wave source concepts. More specifically, we have discussed how to retrieve estimates of the primary responses to a plane-wave source. The retrieved primaries can then be used via standard reverse time migration to produce images free of artefacts related to internal multiples. Whereas previous data domain Marchenko methods are applied to point source gathers and therefore tend to be rather expensive for large datasets, the proposed method provides good imaging results by only involving a small number of primary synthesis steps and the corresponding plane-wave reverse time migration. The plane-wave source primary synthesis algorithm discussed in this paper could then be used as an initial and inexpensive processing step, potentially guiding more expensive target imaging techniques. In this paper we have only discussed 2D examples and internal multiples, but an obvious extension would be allowing surface source primary synthesis in 3D problems as well as incorporating free surface multiples. Finally, applications of data domain Marchenko methods to field data have already been performed. Future work will then focus on applying plane-wave primary synthesis methods to field data too.