# Imaging strategies using focusing functions with applications to a North Sea field

Imaging strategies using focusing functions with applications to a North Sea field Summary Seismic methods are used in a wide variety of contexts to investigate subsurface Earth structures, and to explore and monitor resources and waste-storage reservoirs in the upper ∼100 km of the Earth’s subsurface. Reverse-time migration (RTM) is one widely used seismic method which constructs high-frequency images of subsurface structures. Unfortunately, RTM has certain disadvantages shared with other conventional single-scattering-based methods, such as not being able to correctly migrate multiply scattered arrivals. In principle, the recently developed Marchenko methods can be used to migrate all orders of multiples correctly. In practice however, using Marchenko methods are costlier to compute than RTM—for a single imaging location, the cost of performing the Marchenko method is several times that of standard RTM, and performing RTM itself requires dedicated use of some of the largest computers in the world for individual data sets. A different imaging strategy is therefore required. We propose a new set of imaging methods which use so-called focusing functions to obtain images with few artifacts from multiply scattered waves, while greatly reducing the number of points across the image at which the Marchenko method need be applied. Focusing functions are outputs of the Marchenko scheme: they are solutions of wave equations that focus in time and space at particular surface or subsurface locations. However, they are mathematical rather than physical entities, being defined only in reference media that equal to the true Earth above their focusing depths but are homogeneous below. Here, we use these focusing functions as virtual source/receiver surface seismic surveys, the upgoing focusing function being the virtual received wavefield that is created when the downgoing focusing function acts as a spatially distributed source. These source/receiver wavefields are used in three imaging schemes: one allows specific individual reflectors to be selected and imaged. The other two schemes provide either targeted or complete images with distinct advantages over current RTM methods, such as fewer artifacts and artifacts that occur in different locations. The latter property allows the recently published ‘combined imaging’ method to remove almost all artifacts. We show several examples to demonstrate the methods: acoustic 1-D and 2-D synthetic examples, and a 2-D line from an ocean bottom cable field data set. We discuss an extension to elastic media, which is illustrated by a 1.5-D elastic synthetic example. Interferometry, Body waves, Controlled source seismology, Wave propagation, Wave scattering and diffraction 1 INTRODUCTION Seismic methods using actively fired sources are used in a wide variety of contexts to explore the upper ∼100 km of the Earth’s subsurface. This is important both beneath land and offshore, in order to identify and monitor Earth resources (ores, minerals, hydrocarbons and water reservoirs), potential storage reservoirs in host rock for waste or by-products (nuclear waste or CO2), suitable bedrock for foundations of buildings and other infrastructure, and for examining deeper structures such as the Earth’s lower crust and upper mantle. The family of seismic imaging methods called migration are an integral part of applied seismology, and are used more broadly in medical imaging and non-destructive evaluation of components and structures. Within geophysics, they construct high-resolution images of subsurface structures from recorded seismic survey data. In a standard seismic data acquisition, seismic waves are created by surface sources (e.g. air guns and vibroseis), which propagate through the subsurface and are recorded at an array of receivers (e.g. hydrophones and geophones). Using state-of-the-art methods, images are usually constructed from standard seismic data using so-called reverse-time migration (RTM, Baysal et al.1983; McMechan 1983; Chang 1986). In RTM, a synthetic source mimicking the real source is propagated computationally through a mainly smooth model which approximates the subsurface seismic velocities with estimates that can be obtained from other methods such as velocity analysis, tomography or waveform inversion (Jones 2010). The received wavefield is computationally backpropagated (propagated backwards in time) through this same smooth model. At points in space where these two wavefields coincide over significant periods of time, it is likely that a subsurface scatterer exists. This is because at such locations in the physical seismic experiment, the source field converts into the receiver-side (scattered) field, hence the fact that source and receiver fields match. A common approach is therefore to calculate the zero-time lag cross-correlation (a measure of similarity) between these two numerically estimated wavefields at every point in our computational grid; the set of such values constitutes the final image. In common with most other seismic processing and migration methods, RTM relies on the single-scattering (or first-order Born) assumption (Oristaglio 1989). Seismic data always contain multiply scattered waves (herein called multiples), and these create spurious structures in the final image since within RTM they cannot be backpropagated from the receivers to each correct scattering location without already having knowledge of the other scattering locations. Multiples must therefore be removed from seismic data prior to migration. Several commonly used methods exist to suppress free surface-related and internal multiples from pre-stack data, as described for example in Verschuur (1992), Weglein et al. (1997), Ziolkowski et al. (1999) and Amundsen (2001). Other methods alter the principles of single-scattering-based imaging to account for multiples during the migration process, and may even use them to enhance the final image (Schuster et al.2004; Berkhout & Verschuur 2006; Muijs et al.2007; Malcolm et al.2009; Berkhout 2012). One particular recently developed redatuming method, the Marchenko method (Wapenaar et al.2013), can be used to obtain images of the subsurface devoid of most multiple-related artifacts. It works by solving a so-called Marchenko equation that relates surface reflection data, transmission estimates (waves traveling directly from the surface to each imaging location in the subsurface), and the critical wavefields for constructing images which are the surface-to-subsurface Green’s functions and focusing functions. The Green’s function is the full-wavefield response at a subsurface location due to impulsive sources at the surface, including all orders of scattering. RTM approximates these Green’s functions in the propagation step described above by omitting all scattering since the propagation model is always dominantly smooth. Green’s functions estimated using Marchenko methods, on the other hand, include all orders of scattering. These Green’s functions could therefore be used in place of the propagation steps in RTM. The other important component of Marchenko methods are focusing functions. These are solutions of the wave equation which propagate through a reference medium in such a way that they reach a chosen focusing depth as a localized delta function in time and space at each focus point. The reference medium in which a focusing function is defined contains all true reflectors above its focusing depth but is homogeneous below that depth. The focused field therefore then diverges as a downgoing field through the homogeneous part of the reference medium, and the focusing functions are mathematical entities that cannot be realized in the true Earth. Nevertheless since they account for the true Earth down to the depth of focus, they clearly contain information about the Earth, a fact that we exploit below. Marchenko methods are part of a larger class of high-frequency redatuming methods, which include classical redatuming (Berryhill 1979, 1984; Kinneging et al.1989), and interferometric redatuming (Curtis et al.2006; Bakulin & Calvert 2006; Schuster & Zhou 2006; Curtis & Halliday 2010; Halliday & Curtis 2010). Accurate low-frequency models of the subsurface (in particular, accurate one-way traveltimes) are required for a many high-frequency redatuming methods (excluding, most notably, seismic interferometry), and these have been traditionally been supplied by inversions for the full subsurface model. However, there have been important advances in lower frequency redatuming such as target-oriented inversion (Tang & Biondi 2011), localized full-waveform inversion (Yuan et al.2017) and box tomography (Masson & Romanowicz 2017), which may be beneficial for target-oriented applications by providing more detailed models of subsurface target areas. Marchenko methods have recently been used for a variety of purposes, including internal multiple attenuation (Meles et al.2015; da Costa Filho et al.2017), target-oriented imaging (Wapenaar et al.2014; Ravasi et al.2016), overburden elimination (van der Neut & Wapenaar 2016), synthesizing primaries (Meles et al.2016, 2017) and virtual subsurface wavefield estimation (Wapenaar et al.2016). It has also been extensively used for imaging (e.g. Wapenaar et al.2014; Singh et al.2015; da Costa Filho et al.2015). So far, most Marchenko imaging methods have required computing the solution to the Marchenko equation at each imaging point—a costly endeavor. Alternatively, they have been used to invert for a redatumed reflection response devoid of overburden interactions which is the used in conventional imaging (Wapenaar et al.2014; Ravasi et al.2016; Staring et al.2017). Herein, we show how Marchenko focusing functions computed at far fewer locations can be used to obtain significantly better images than RTM. This is achieved by substituting source and receiver wavefields in conventional RTM by focusing functions. We identify three different source/receiver wavefield pairs which give rise to three distinct imaging methods, all of which account for multiple scattering in the Earth’s subsurface. Finally, we propose a post-imaging filter based on da Costa Filho & Curtis (2016) which can be used to enhance the results of two of these methods. We compare all methods to conventional RTM in 1-D and 2-D acoustic synthetic media, as well as in a test involving real data from a 2-D line from an ocean bottom cable (OBC) field survey in the North Sea. Finally, we discuss extensions to elastic media using a 1.5-D elastic synthetic medium in a worked example. 2 WAVE THEORY The idea of using focusing operators to better image a medium is not new. It has been used in several different contexts, including acoustic time-reversal focusing experiments (Fink 1992), as well as migration and inversion of seismic data (Thorbecke 1997; Berkhout & Verschuur 2005). The main practical constraint imposed in geophysics is that data are usually only collect on one side of a medium (the Earth’s surface) in contrast to medical imaging which the subject can be entirely surrounded by sources and receivers. Rose (2001) published the first method to obtain focusing wavefields using single-sided data in a medium with several scatterers. This was generalized to provide 3-D focusing functions in 3-D acoustic media through the work of Wapenaar et al. (2013), and to 3-D elastic media by da Costa Filho et al. (2014) and Wapenaar & Slob (2014). Similarly to Wapenaar et al. (2014), in our notation, the Marchenko method provides Green’s functions G(xF, x0, t) as well as focusing functions denoted f1(x0, xF, t) where $$\mathbf {x}_0 = (\hat{\mathbf {x}}_0, z_0)$$ is a location at the surface z = z0, $$\mathbf {x}_F = (\hat{\mathbf {x}}_F, z_F)$$ is the focusing location at depth z = zF and t is time. Here, $$\hat{\mathbf {x}}$$ represents the horizontal coordinates of x. A focusing function is defined in a reference medium which is homogeneous below zF, and is equal to the true (unknown) medium between the surface and the focusing depth, as shown in Fig. 1. They have the characteristic that at zero time, the wavefield observed at the depth of zF focuses: that is, it becomes a delta function at $$\hat{\mathbf {x}}_F$$. These functions are commonly decomposed with regard to their directions of propagation related to their first coordinate (on the surface at which they are emitted or received). Therefore, with a focus at xF, $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$ refers to the downgoing focusing field departing from x0, and $$f_1^-(\mathbf {x}_0, \mathbf {x}_F, t)$$ is the upgoing field which then arrives at x0. The superscript + is inherited from the convention that our z-axis is positive downwards. Similarly, G − (xF, x0, t) and G + (xF, x0, t) are defined as up- and downgoing Green’s functions, respectively, where now the superscripts + and − refer to the direction of propagation at the subsurface focus point xF. Figure 1. View largeDownload slide Reference medium in which focusing functions $$f_1^\pm (\mathbf {x}_0, \mathbf {x}_F)$$ are defined. Figure 1. View largeDownload slide Reference medium in which focusing functions $$f_1^\pm (\mathbf {x}_0, \mathbf {x}_F)$$ are defined. Focusing fields along a horizontal line at depth have been used for imaging structures from below (Wapenaar et al.2014; Ravasi et al.2016). This shows that focusing functions contain information about the reflection properties of the medium. With an alternative interpretation of these fields, we now show other ways in which they can be used for imaging. The focusing functions are related to the surface reflection data $$R(\mathbf {x}_0^{\prime \prime }, \mathbf {x}_0, t)$$ and to the subsurface Green’s function $$G^-(\mathbf {x}_F, \mathbf {x}_0^{\prime \prime }, t)$$ by the following equation:   \begin{eqnarray} G^-(\mathbf {x}_F, \mathbf {x}_0^{\prime \prime }, t) &=& - f_1^-(\mathbf {x}_0^{\prime \prime }, \mathbf {x}_F, t) \nonumber\\ &&+ \int \limits _{\partial \mathbb {D}_0}\int \limits _{-\infty }^\infty R(\mathbf {x}_0^{\prime \prime }, \mathbf {x}_0, t-\tau ) f^+_1(\mathbf {x}_0, \mathbf {x}_F, \tau )\,\mathrm{d}{\tau }\,\mathrm{d}^2{\mathbf {x}_0},\nonumber\\ \end{eqnarray} (1)where $$R(\mathbf {x}_0^{\prime \prime }, \mathbf {x}_0, t)$$ corresponds to the pressure of the reflected field recorded at $$\mathbf {x}_0^{\prime \prime }$$ created by vertical force sources at x0, multiplied by −2 (Broggini et al.2014). Eq. (1) is obtained by considering two wavestates A and B: wavestate A corresponds to Green’s functions which exist in the true medium, while the focusing functions correspond to wavestate B and exist in a reference medium which is truncated below xF (Fig. 1). For general wavestates which propagate through media that are identical between two arbitrary boundaries $$\partial \mathbb {D}_0$$ and $$\partial \mathbb {D}_F$$, Wapenaar et al. (2014) derive a one-way acoustic reciprocity theorem for pressure normalized wavefields in the angular frequency domain $$p_A^\pm (\mathbf {x}, \omega )$$ and $$p_B^\pm (\mathbf {x}, \omega )$$:   \begin{eqnarray} &&{-\int \limits _{\partial \mathbb {D}_0} \frac{1}{\rho }\lbrace p_A^+ (\partial _z p_B^-) + p_A^- (\partial _z p_B^+)\rbrace \,\mathrm{d}^2{\mathbf {x}_0} }\nonumber\\ &&= \int \limits _{\partial \mathbb {D}_F} \frac{1}{\rho }\lbrace (\partial _z p_A^+)p_B^- + (\partial _z p_A^-) p_B^+\rbrace \,\mathrm{d}^2{\mathbf {x}_F}, \end{eqnarray} (2)where arguments of $$p_A^\pm$$ and $$p_B^\pm$$ have been omitted, and pressure normalized fields are those whose sum equals that of the pressure, that is p+ + p− = p. Substituting appropriate boundary conditions in eq. (2) yields eq. (1) (see appendix A of Wapenaar et al.2014). By choosing different wavestates, as long as they are defined in media which coincide between two arbitrary boundaries $$\partial \mathbb {D}_0$$ and $$\partial \mathbb {D}_F$$, we may use the same one-way acoustic reciprocity theorem in eq. (2) to obtain novel relationships. We use this freedom to consider wavestate A to refer to focusing wavefields in the reference medium as before, but now consider wavestate B to be the Green’s functions $$\overline{G}$$ also in the reference medium, as opposed to the true medium which was assumed when deriving eq. (1). More explicitly, we substitute $$p_A^\pm (\mathbf {x}, \omega ) = f_1^\pm (\mathbf {x}, \mathbf {x}_F^{\prime }, \omega )$$ and $$p_B^\pm (\mathbf {x}, \omega ) = \overline{G}^\pm (\mathbf {x}, \mathbf {x}_0^{\prime }, \omega )$$ in eq. (2) to obtain   \begin{eqnarray} &&{-\!\int \limits _{\partial \mathbb {D}_0} \!\frac{1}{\rho }\lbrace f_1^+(\mathbf {x}_0, \mathbf {x}_F^{\prime }) [\partial _z \overline{G}^-(\mathbf {x}_0,\mathbf {x}_0^{\prime })] \!+\! f_1^-(\mathbf {x}_0, \mathbf {x}_F^{\prime }) [\partial _z \overline{G}^+(\mathbf {x}_0,\mathbf {x}_0^{\prime })]\rbrace \mathrm{d}^2{\mathbf {x}_0}}\nonumber \\ && = \int \limits _{\partial \mathbb {D}_F} \frac{1}{\rho }\lbrace [\partial _z f_1^+(\mathbf {x}_F, \mathbf {x}_F^{\prime })] \overline{G}^-(\mathbf {x}_F, \mathbf {x}_0^{\prime }) \!+\! [\partial _z f_1^-(\mathbf {x}_F, \mathbf {x}_F^{\prime })]\nonumber\\ &&\quad\times \overline{G}^+(\mathbf {x}_F, \mathbf {x}_0^{\prime })\rbrace \,\mathrm{d}^2{\mathbf {x}_F}, \end{eqnarray} (3)where the frequency dependence has been suppressed for notational brevity. A property of the reference medium used for both wavestates is that no reflections come from above $$\partial \mathbb {D}_0$$ or from below $$\partial \mathbb {D}_F$$. This condition is such that the vertical derivative $$\partial _z\overline{G}^+(\mathbf {x}_0, \mathbf {x}_0^{\prime })$$ vanishes on $$\partial \mathbb {D}_0$$, and both $$\overline{G}^-(\mathbf {x}_F, \mathbf {x}_0^{\prime })$$ and $$\partial _z f_1^-(\mathbf {x}_F, \mathbf {x}_F^{\prime })$$ also vanish on $$\partial \mathbb {D}_F$$. In addition, following previous Marchenko work, we assume that $$\partial _z \overline{G}^+(\mathbf {x}_0, \mathbf {x}_0^{\prime }) = -\frac{1}{2} \iota \omega \rho (\mathbf {x}_0^{\prime })\delta (z_0 -z_0^{\prime })$$ at $$\partial \mathbb {D}_0$$, that is, Green’s functions have a vertical force point-source mechanism (and that is the only downgoing field at the surface since there are no reflectors above that level). Applying these conditions to eq. (3) yields   $$\int \limits _{\partial \mathbb {D}_0} \overline{R}(\mathbf {x}_0^{\prime }, \mathbf {x}_0)f_1^+(\mathbf {x}_0, \mathbf {x}_F^{\prime }) \,\mathrm{d}^2{\mathbf {x}_0} = f_1^-(\mathbf {x}_0^{\prime }, \mathbf {x}_F^{\prime }),$$ (4)where $$\overline{R}(\mathbf {x}_0^{\prime }, \mathbf {x}_0) \equiv \frac{-2}{\iota \omega \rho (\mathbf {x}_0)}\partial _z \overline{G}^-(\mathbf {x}_0,\mathbf {x}_0^{\prime })$$. If we define the volume encompassed between depths ε above and below $$\partial \mathbb {D}_0$$ as V (Fig. 1), we may rewrite eq. (4) as   \begin{eqnarray} \int \limits _V \overline{R}(\mathbf {x}_0^{\prime }, \mathbf {x}_0)f_1^+(\mathbf {x}_0, \mathbf {x}_F^{\prime })\delta (z-z_0) \mathrm{d}^3{\mathbf {x}_0} \!=\! \int \limits _{-\varepsilon }^{\varepsilon } f_1^-(\mathbf {x}_0^{\prime }, \mathbf {x}_F^{\prime })\delta (z-z_0)\mathrm{d}{z},\nonumber\\ \end{eqnarray} (5)where z0 is the depth coordinate of x0. Considering that similarly to R, $$\overline{R}$$ may be viewed as a scaled dipole-to-monopole surface response in the reference medium multiplied by −2, it is effectively a Green’s function. As such, the left-hand side of eq. (5) is a volume integral between a Green’s function and a source function, namely $$f_1^+(\mathbf {x}_0, \mathbf {x}_F^{\prime })\delta (z-z_0)$$. A direct consequence of this is that a line source of $$f_1^+$$ produces an $$f_1^-$$ field localized along $$\partial \mathbb {D}_0$$. Simply put, $$f_1^-$$ is the received reference medium response to source field $$f_1^+$$, where both the action of the source field and the recorded response are at the surface. The Marchenko method is then not only useful for producing redatumed Green’s functions, but also focusing functions which are essentially virtual data acquisitions (surveys) that take place in an imagined medium which contains the true medium above the focusing depth, that is, the reference medium in Fig. 1. Here, we explore some of interesting properties of these virtual surveys which are derived from the focusing functions, and propose imaging methods which exploit these virtual acquisitions, in particular an imaging method that targets specific reflectors. 3 IMAGING METHODS In RTM, each seismic point source is injected into a computational model of the medium that is a best-possible estimate of the true earth subsurface, and is propagated forwards in time to simulate the downgoing field from the source throughout the subsurface. The received wavefield is injected into the same model and propagated backwards in time to estimate what this field was like throughout the subsurface before it was recorded. Conveniently, this second step can be applied using a standard forward propagation simulation of the time-reversed input data because the wave equation is symmetric in time (forward and backward propagations are identical mathematical operations). The two wavefields are then combined to construct an image, the particular manner in which the fields are combined being called the imaging condition. When imaging with a virtual survey, the same process can be employed. If we use $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$ as the source field, it must be injected along the whole surface source array simultaneously (as opposed to being a point source at a specific location as in conventional RTM). In this case, the receiver array measures $$f_1^-(\mathbf {x}_0, \mathbf {x}_F, t)$$ as opposed to the recorded surface data, hence $$f_1^-$$ must be injected along the receiver array and backpropagated in time. The nature of the focusing source/receiver pair will have two effects on the image produced using a single virtual focus point xF. First, it will be more localized than an outward-spreading field from a point source as the new source focuses towards the focus point. While the outward-spreading wavefield will tend to reach all subsurface locations given sufficient time to propagate, the focusing source will be approximately contained in the triangle (or downwards-pointing cone, in 3-D) whose base is the source array on it and whose tip is the focus point. Second, since it is defined in the reference medium, it cannot be used to image anything below xF. In light of this, in order to image the medium using these virtual surveys, one has to compute the image for several focus points in order to illuminate the whole subsurface, and these foci must be located below the deepest reflector of interest. The velocity model used to image these virtual surveys is essentially the same as used in conventional imaging, with the difference that they need extend only down to the virtual source depth. The nature of the focusing functions has been explored and exemplified in Slob et al. (2014) and Wapenaar et al. (2014). They show that $$f_1^+(\mathbf {x}, \mathbf {x}_F,t)$$ contains a time-reversed direct wave $$f_{1,d}^+$$ which scatters as it propagates through the reference medium. A coda of multiply scattered waves, $$f_{1,m}^+$$, destructively interferes with the scattering of $$f_{1,d}^+$$ in such way that the only arrival reaching the focusing depth is the time-reversed direct wave as a delta function localized at xF. Note that as $$f_{1,d}^+$$ propagates through the medium it will reflect off every interface in the reference medium; therefore $$f_1^-$$ will contain, amongst other events, primaries resulting from source field $$f_{1,d}^+$$. Above we suggested to migrate source/receiver pair $$f_1^+$$/$$f_1^-$$ from several focus points. However, knowing that $$f_1^-$$ contains a primary from the direct wave $$f_{1,d}^+$$ suggests that results may improve when using $$f^+_{1,d}$$ instead of $$f_1^+$$ as source-side wavefield. This is because RTM assumes that fields consist only of primaries which can be downwards propagated using direct waves, an assumption which is valid when using $$f_{1,d}^+$$, but not when using $$f_1^+$$. Indeed, it has been shown that using a single-event source field (e.g. a direct wave) in the imaging condition eliminates some crosstalk between source-side wavefields and unrelated events in the receiver-side wavefield (Wapenaar et al.1987; da Costa Filho et al.2015). This is corroborated by other recent work: focusing functions have the property of focusing not only in space (at a single subsurface location), but also in time (nothing arrives at the focusing depth before or after the focusing pulse). Meles et al. (2017) have shown that a direct consequence of focusing is that only primaries are present in $$f_1^-$$, at least in horizontally layered media. All multiples are canceled by destructive interference with $$f_{1,m}^+$$ and any residual multiple would destroy the focus. Some of these primaries are created by $$f_{1,d}^+$$, and some of them are created by $$f_{1,m}^+$$. Meles et al. (2017) also showed that in 1-D media the last-arriving event in $$f_1^-$$ is a primary created by $$f_{1,d}^+$$. Therefore, for any given focus point, we may generate an upgoing focusing gather $$f_1^-(\mathbf {x}_0, \mathbf {x}_F, t)$$ such that the last event of this gather will be a primary, $$f_{1,p}^-(\mathbf {x}_0, \mathbf {x}_F, t)$$, originating from a known source, $$f_{1,d}^+(\mathbf {x}_0, \mathbf {x}_F, t)$$. This allows us to migrate only primaries by considering source/receiver pairs composed of $$f_{1,d}^+$$ and $$f_{1,p}^-$$. The primary recovered is the one reflected from the closest reflector immediately above the focusing location, and thus can be used to image that reflector individually. If all reflectors are to be imaged, subsurface focus points must therefore be located below every reflector. We thus derive three new methods to image the subsurface, each consisting or an RTM-type operation using different source/receiver fields. First, the pair $$f_1^+/f_1^-$$ can be used to image the medium above the focus points. Second, we can restrict the source field to the direct wave so as to use only the energy from $$f_{1,d}^+/f_1^-$$ which reduces crosstalk artifacts in the image. The problem with the second method is that some of the primary energy in $$f_1^-$$ was created by the coda $$f_{1,m}^+$$ which is now not included in the source field, hence errors in amplitude may arise and some crosstalk artifacts remain. The third method removes all crosstalk by using pair $$f_{1,d}^+/f_{1,p}^-$$, whose receiver-side field consists of a single primary reflection from the first reflector above the focus point. The third method has the property that it images a single reflector, allowing targeted images to be produced, but has the disadvantage that the last arriving energy in $$f_1^-$$ must be identified and isolated prior to imaging. Another possibility to remove crosstalk present in the first two methods is to use the post-imaging filter proposed by da Costa Filho & Curtis (2016). In that work, two images IC and IF are combined to create a new image using the following expression   $$I_{\rm H}(\mathbf {x}_I) = \frac{E[I_{\rm C}](\mathbf {x}_I)\, I_{\rm F}(\mathbf {x}_I) + E[I_{\rm F}](\mathbf {x}_I)\, I_{\rm C}(\mathbf {x}_I)}{E[I_{\rm C}](\mathbf {x}_I)+E[I_{\rm F}](\mathbf {x}_I)}$$ (6)where xI is the image point, and E[ · ] is the non-negative Hilbert envelope applied to each vertical image trace individually. The combined image IH has the property that events which appear in both constituent images remain in IH, while events which appear in only one or the other of IC or IF are attenuated. Below, we apply the combined images method using as constituent images the conventional RTM image (IC), and a focusing image (IF). In both images, true reflectors appear at the same locations while artifacts may appear in different spatial locations. The combined images can then be used to attenuate those artifacts which appear at different locations. 4 NUMERICAL EXAMPLES We now present several examples to demonstrate and explore these methods. The first two examples are in acoustic media: a 1-D synthetic model with variable velocity, and a 2-D synthetic model with vertical and horizontal density variations. We also apply acoustic methods to data from a 2-D line from an OBC field acquisition in the North Sea. Finally, in Section 5, we provide a simple worked elastic example. 4.1 1-D synthetic First, we present a 1-D synthetic example with focusing functions that have been calculated analytically. The medium and a focusing wavefield are shown in Fig. 2. Events associated with the downwards pointing arrows are injected as part of $$f_1^+(0\,\mathrm{km}, 4.5\,\mathrm{km}, t)$$, with the earliest (leftmost) event being $$f_{1,d}^+$$. The other injected components, $$f_{1,m}^+$$, are seen to contribute no energy to the focus; they are present to destroy events that would otherwise reflect downwards and destroy the focus at depth. It is also seen that the last (rightmost) event arriving at the surface as part of $$f_1^-(0\,\mathrm{km}, 4.5\,\mathrm{km}, t)$$ is a sum of primaries, one of which originates from $$f_{1,d}^+$$. Figure 2. View largeDownload slide Synthetic 1-D model overlain by a focusing wavefield. The black circle denotes the focusing point at depth z = zf and time t = 0, and the dashed line denotes the focusing depth. Arrows pointing downwards depict events in $$f_1^+$$ and those pointing upwards depict events in $$f_1^-$$. Note that since there is only one spatial dimension, the horizontal axis is time. Figure 2. View largeDownload slide Synthetic 1-D model overlain by a focusing wavefield. The black circle denotes the focusing point at depth z = zf and time t = 0, and the dashed line denotes the focusing depth. Arrows pointing downwards depict events in $$f_1^+$$ and those pointing upwards depict events in $$f_1^-$$. Note that since there is only one spatial dimension, the horizontal axis is time. Next, we migrate source/receiver pairs $$f_1^+$$ and $$f_1^-$$, and $$f_{1,d}^+$$ and $$f_1^-$$ using a smooth propagation model, and compare it to conventional RTM, which we will from now on simply refer to as RTM. These images (traces in 1-D examples) are shown, respectively, in Figs 3(a)–(c). All figures have imaged the true reflectors exactly the same way. These appear at incorrect depths (0.62, 1.35 and 3.24 km) because the migration velocity model used a constant velocity of 1 km s−1 in order to eliminate any backscattering which might occur in the imaging step. They also appear with different absolute and relative amplitudes: in the first two images the receiver wavefield is $$f_1^-$$ which has several contributions to its primaries when compared to RTM. Moreover, the first two images have different source-time functions, which again result in different amplitudes. However, this effect is minor: taking the maximum energy of the third reflector reveals peaks at 10.303 × 10−2, 10.304 × 10−2 and 8.964 × 10−2 for Figs 3(a)–(c), respectively. Figure 3. View largeDownload slide Images from sources/receiver wavefields composed from (a) $$f_1^+$$ and $$f_1^-$$, and (b) $$f_{1,d}^+$$ and $$f_1^-$$. (c) Conventional RTM image. (d) Image obtained by combining (a) and (c) using eq. (6). Figure 3. View largeDownload slide Images from sources/receiver wavefields composed from (a) $$f_1^+$$ and $$f_1^-$$, and (b) $$f_{1,d}^+$$ and $$f_1^-$$. (c) Conventional RTM image. (d) Image obtained by combining (a) and (c) using eq. (6). Fig. 3(a) shows two artifacts, one which results from the interaction between $$f_{1,d}^+$$ and an unrelated event in $$f_1^-$$. This same artefact is seen in Fig. 3(b), albeit with smaller amplitude. Fig. 3(a) also contains a second, smaller artefact which results from the interaction between $$f_{1,m}^+$$ and an unrelated event in $$f_1^-$$. This crosstalk artefact is absent from Fig. 3(b) as expected from the discussion above. In this simple model, we do not show an individual reflector imaged using the last-arriving event from this focus point (the third method above), as it suffices to see that Figs 3(a) and (b) show the last reflector is correctly recovered. Therefore, the last-arriving event in $$f_1^-$$ has the traveltime of the desired primary and correctly images the deepest reflector. However, using RTM the last imaged reflector is spurious, and we see a total of four artifacts in the area imaged (two between the second and third reflectors, and two after the third reflector). These RTM artifacts generally appear at different locations than those in the focusing-related images. Hence, the combination of a focusing image with the RTM image using the method of da Costa Filho & Curtis (2016) is likely to produce an image devoid of most artifacts. Indeed, the result shown in Fig. 3(d), which exhibits the image obtained by combining images in panels (Figs 3a and c), shows a virtually artefact-free image. 4.2 2-D synthetic We now use a 2-D synthetic acoustic model with a constant velocity of 2.4 km s−1 and densities shown in Fig. 4. Varying only densities allows straight ray paths to be assumed, so that events observed in the data and focusing functions can be interpreted more easily. However, this is not a limitation of the new methods as Marchenko focusing functions can also be found for media with varying velocities, including highly complex media (Behura et al.2014; Vasconcelos et al.2015; Jia et al.2017). Figure 4. View largeDownload slide Acquisition and imaging geometry for 2-D synthetic model. Stars and triangles represent the 201 colocated surface sources and receivers. White dots represent the locations where focusing functions $$f_1^\pm$$ were computed. Figure 4. View largeDownload slide Acquisition and imaging geometry for 2-D synthetic model. Stars and triangles represent the 201 colocated surface sources and receivers. White dots represent the locations where focusing functions $$f_1^\pm$$ were computed. We choose several depth levels, each containing 151 focusing locations at which we compute focusing functions as shown by the white dots in Fig. 4. The bottom-most depth level at 1.6 km is used for imaging with pairs $$f_1^+$$ and $$f_1^-$$, and $$f_{1,d}^+$$ and $$f_1^-$$, displayed, respectively, in Figs 5(a) and (b); in this case, their Marchenko reference medium has the same reflectors as the full medium since the focusing level is beneath all reflectors. Figure 5. View largeDownload slide Images using source/receiver wavefields (a) $$f_1^+$$ and $$f_1^-$$, (b) $$f_{1,d}^+$$ and $$f_1^-$$ and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. (c) Conventional RTM image. (d) Combined image using images in panels (a) and (c). (e) Combined image using images in panels (b) and (c). Images in panels (a)–(e) use only virtual points in the bottom-most line of focus points, while (f) uses all virtual points shown in Fig. 4. Arrows indicate spurious reflectors. All panels are clipped at 10 per cent of their maximum amplitude. Figure 5. View largeDownload slide Images using source/receiver wavefields (a) $$f_1^+$$ and $$f_1^-$$, (b) $$f_{1,d}^+$$ and $$f_1^-$$ and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. (c) Conventional RTM image. (d) Combined image using images in panels (a) and (c). (e) Combined image using images in panels (b) and (c). Images in panels (a)–(e) use only virtual points in the bottom-most line of focus points, while (f) uses all virtual points shown in Fig. 4. Arrows indicate spurious reflectors. All panels are clipped at 10 per cent of their maximum amplitude. The focusing functions at other depth levels are used to image with last-arriving primaries only, that is, using the respective functions $$f_{1,d}^+$$ and $$f_{1,p}^-$$ for each depth level. The summation of these six images (each of which contains one reflector corresponding to a single depth level) is seen in Fig. 5(f). In order to use this method to recover all reflectors using only primaries, focus points are necessary below each interface. In this synthetic example, we choose a set of depths which we know are guaranteed to recover all reflections. In field data, this is not guaranteed. In addition, since this method requires picking, we do not believe that it is useful in areas with little knowledge of the imaging target. This method is much more appropriate for areas with a priori knowledge of the region wherein the imaging target is located, at which we wish to construct a more focused image. This is explored in the following field data example. As expected, Fig. 5(a) contains several artifacts resulting from crosstalk as indicated by the white arrows. These are reduced in Fig. 5(b), which features fewer spurious interfaces and artifacts with smaller amplitudes. Fig. 5(c) shows a standard RTM image, which has more artifacts in deeper sections compared to Figs 5(a) and (b). Since these artifacts appear in different locations, we combined the focusing-derived images with the standard RTM image to generate the combined images in Figs 5(d) and (e). These images are virtually free from artifacts, and also have fewer acquisition imprints (i.e. low-frequency artifacts common in RTM which can be seen mostly around the top reflector) than their constituent images. In addition, the combined image in Fig. 5(e) has slightly fewer artifacts than that in Fig. 5(d), consistent with the fact that fewer artifacts appear in Fig. 5(b) compared to Fig. 5(a). Fig. 5(f) also shows near perfect reconstruction of primaries. No coherent spurious interfaces are imaged, and all true reflectors have been recovered. However, this performance comes at the cost of picking the last event in $$f_1^-$$ (in this case performed automatically), as well as an increased computational cost due to the increased number of focus points. 4.3 Field data The methods were applied to an OBC field data set from the seabed above Volve oil field, located in the gas/condensate-rich Sleipner area 200 km off the North Sea coast of Norway (Fig. 6). Prior to applying the Marchenko method, the data were pre-processed with spatial interpolation, up/down wavefield separation, redatuming by multidimensional deconvolution, among others; the pre-processing is described by Ravasi et al. (2015, 2016), who detail the steps required to ensure the assumptions of the Marchenko method are met. Figure 6. View largeDownload slide Geographic location of the Volve field. Figure 6. View largeDownload slide Geographic location of the Volve field. Fig. 7(a) shows the velocity model used for the Marchenko method and for imaging, and Fig. 7(b) shows a magnification of the portion of the model in which we are interested. In both figures, the black dots show the 411 locations where focusing functions $$f_1^\pm$$ were computed using the pre-processed data. Figure 7. View largeDownload slide (a) Acquisition and imaging geometry. Triangles represent the 235 OBC receivers and stars represent the sources redatumed to the same receiver locations. The black box represents the imaged area, which is magnified in (b). Black dots represent locations where focusing functions were computed and the larger red dot shows the focusing location for functions in Fig. 8. The black arrow indicates a particularly strong contrast in the velocity model. Figure 7. View largeDownload slide (a) Acquisition and imaging geometry. Triangles represent the 235 OBC receivers and stars represent the sources redatumed to the same receiver locations. The black box represents the imaged area, which is magnified in (b). Black dots represent locations where focusing functions were computed and the larger red dot shows the focusing location for functions in Fig. 8. The black arrow indicates a particularly strong contrast in the velocity model. In Fig. 8, we show an example of a focusing function from the location marked by the large red dot in Fig. 7. In Fig. 8(c), we magnify a portion of the upgoing field shown in Fig. 8(b). Our picked final event in this gather is shown in Fig. 8(d). While it is clear from Fig. 8(c) that there is energy below the picked event, we verified that most of that energy is not coherent across all focusing locations (e.g. coherent energy must appear in every $$f_1^-$$ gather, and must at no point be crossed by the window used in the Marchenko iteration). Since the strong contrast in the velocity model (black arrow in Fig. 7b) appears across the entire model, we expected to see its associated primary along all focusing locations. The identified event shown in Fig. 8(d) was picked manually for a regularly spaced subset of the focusing locations, and the rest were picked automatically by applying interpolated muting windows to the focusing functions. The imaging examples below show that this procedure recovered the correct primary. Figure 8. View largeDownload slide (a) $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$ and (b) $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$, for $$\mathbf {x}_F = (6.7\,\mathrm{km}, 3.3\,\mathrm{km})$$ (large red dot in Fig. 7). The area enclosed by the black box in (b) is magnified in (c), and the result from windowing its last arrival is shown in (d). All panels are clipped at 50 per cent of the maximum amplitude of (b). Figure 8. View largeDownload slide (a) $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$ and (b) $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$, for $$\mathbf {x}_F = (6.7\,\mathrm{km}, 3.3\,\mathrm{km})$$ (large red dot in Fig. 7). The area enclosed by the black box in (b) is magnified in (c), and the result from windowing its last arrival is shown in (d). All panels are clipped at 50 per cent of the maximum amplitude of (b). Fig. 9 shows the migrated images of the section of interest. Figs 9(a) and (b) are obtained, respectively, by using source/receiver wavefield pairs $$f_1^+$$ and $$f_1^-$$, and $$f_{1,d}^+$$ and $$f_1^-$$. These are virtually the same, with no notable differences other than in illumination and amplitudes. This happens because the coda of the downgoing focusing function $$f_{1,m}^+$$ has relatively low amplitude, and therefore contributes much less to the image than the direct part $$f_{1,d}^+$$, as evidenced in Fig. 8(a). The section of interest of the RTM image is shown in Fig. 9(c). The RTM image was provided from previous processing and fewer imaging points, spaced 8 m apart rather than 4 m as used in Figs 9(a) and (b). This may cause a decrease in resolution. In addition, the RTM image shows severely disrupted reflectors near the centre of the image when compared to Figs 9(a) and (b). These discontinuities may be caused by either multiples or backscattering which cross the reflector; such effects are not expected in our method where multiples do not play a role and where backscattering effects are diminished in the reference medium. However, the RTM image shows more continuity towards the edges of the image. The focusing images, which have low illumination near lateral edges, make identifying continuous reflectors harder at those locations. Figure 9. View largeDownload slide Images using source/receiver wavefields (a) $$f_1^+$$ and $$f_1^-$$, (b) $$f_{1,d}^+$$ and $$f_1^-$$ and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. (c) Conventional RTM image. (d) Combined image using images in panels (a) and (c). (e) Combined image using images in panels (b) and (c). Arrows indicate reflectors which appear better imaged in (a), (b) or (f) than in (c). Figure 9. View largeDownload slide Images using source/receiver wavefields (a) $$f_1^+$$ and $$f_1^-$$, (b) $$f_{1,d}^+$$ and $$f_1^-$$ and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. (c) Conventional RTM image. (d) Combined image using images in panels (a) and (c). (e) Combined image using images in panels (b) and (c). Arrows indicate reflectors which appear better imaged in (a), (b) or (f) than in (c). We show the combined images in Figs 9(d) and (e). In order to combine them, we interpolated the RTM image onto the finer grid used in the focusing images which may cause a decrease in resolution in the combined images compared to panels (a) and (b). The combined images also show generally fewer acquisition imprints. However, they are severely harmed by the discontinuities which appear in RTM, leading to a worse interpretability than their constituent non-RTM images. Nevertheless, continuous reflectors are much clearer in the combined image than in either of the constituent images. Fig. 9(f) shows the migration of the primary picked from $$f_1^-$$ at each point along the focus point line. An example of these gathers was shown in Fig. 8. We expected the primary to be related to the strong velocity contrast in the model indicated in Fig. 7. Indeed, Fig. 9(f) shows that we have individually imaged that strong contrast in the model, which also appears prominently in Figs 9(a) and (b). 5 DISCUSSION The results shown above demonstrate the advantages and limitations of each imaging method. The first method uses $$f_1^+$$ and $$f_1^-$$ as source/receiver wavefield pairs to create the image and provides images that are inherently different from those obtained using standard RTM, provided that accurate focusing is achieved in the Marchenko method. Using 1-D synthetic data where the focusing was near-perfect, this first method and standard RTM show, as expected, artifacts in different locations. Also in the 2-D synthetic model, artifacts mostly appear in different locations, suggesting that the focusing achieved was also acceptable. For field data, though hard to see from the individual images themselves, the combined images (discussed below) also show that artifacts from RTM appear in a different locations than in our first method. Nevertheless, our first method exhibits crosstalk between $$f_{1}^+$$ and unrelated events in $$f_{1}^-$$. In order to reduce these, we tested substituting $$f_1^+$$ for its direct wave component $$f_{1,d}^+$$. We note that this is not the only possible approach to reduce artifacts, with deconvolution of receiver wavefields by source wavefields being another option. Deconvolution, however, is often unstable. Our approach, on the other hand, is stable and reliably produces images with fewer artifacts, as shown in the synthetic models. However, for the field data it showed no significant improvement over our previous method. This is because in this particular case the coda components are weak as evidenced in Fig. 8(a), so the direct wave also dominates the downgoing field when we use the full $$f_1^+$$ as the source field. Nevertheless, when compared to RTM, this method exhibits improvements in reflector continuity, likely because of the effects of backscattering which are not accounted for in RTM. Similarities between the images obtained by our methods and RTM were exploited by using the combined imaging method of da Costa Filho & Curtis (2016). Our synthetic results show near perfect recovery of all reflectors with minimal interference from multiple or acquisition-related artifacts. Using field data, the combined images also reduced artifacts, but the low quality of the RTM image and the lack of continuity near the edges of the focusing images degraded the resulting combined image. Our final imaging method uses the windowed $$f_1^-$$ functions, denoted $$f_{1,p}^-$$. These retain only their last events which are primaries related to the reflector immediately above the focus point. This method allows the imaging of individual reflectors without interference (crosstalk) from other parts of the wavefield. However, if an image of the entire medium is desired, each reflector must be imaged separately which in turn requires more computational power to calculate focusing functions at an increased number of focus points. Nevertheless, it has been shown to provide clean images of individual reflectors, even in complex field data. While picking for the 2-D synthetic data was done automatically, picking consistent arrivals in the field data was performed manually to ensure that other events in $$f_1^-$$ where not picked in error by the automatic picking routine. This increases processing time significantly, and hence this method may be most applicable for target-oriented applications for which the extra time may be worthwhile if it results in a better image of a reservoir. The methods shown for acoustic media have a straightforward extension to elastic media, with certain caveats which we highlight below using the synthetic 1.5-D model in Fig. 10(a). da Costa Filho et al. (2014), Wapenaar (2014) and Wapenaar & Slob (2014) have shown how focusing functions may be computed in elastic media. Our first and second methods, which do not require windowing, may then be used to compute elastic images by migrating source/receiver pairs composed of down- and upgoing elastic focusing functions. The resulting SS wave images are shown in Fig. 10(d) using fields $$f_1^+$$ and $$f_1^-$$, and Fig. 10(b) using $$f_{1,d}^+$$ and $$f_1^-$$. In this simple model, there is no crosstalk between $$f_{1,m}^+$$ and events in $$f_1^-$$ in Fig. 10(d) or between $$f_{1,d}^+$$ and events in $$f_1^-$$ in Fig. 10(b), but there are S-to-P conversion-related artifacts. They compare well to RTM (Fig. 10c), and combining our image in Fig. 10(b) with RTM provides the virtually artefact-free image shown in Fig. 10(e). Figure 10. View largeDownload slide (a) 1.5-D synthetic density model used to generate elastic data. Colocated sources and receivers are shown as stars and triangles, and focus locations are shown as dots. P and S velocities are constant at 2.5 and 1.3 km s−1, respectively. SS images are constructed using source/receiver wavefields (b) $$f_{1,d}^+$$ and $$f_1^-$$, (d) $$f_{1}^+$$ and $$f_1^-$$, and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. RTM is shown in (c), and its combined image with (b) is shown in (e). White arrows indicate artifacts. Images are clipped at 5 per cent of their maximum amplitude. Figure 10. View largeDownload slide (a) 1.5-D synthetic density model used to generate elastic data. Colocated sources and receivers are shown as stars and triangles, and focus locations are shown as dots. P and S velocities are constant at 2.5 and 1.3 km s−1, respectively. SS images are constructed using source/receiver wavefields (b) $$f_{1,d}^+$$ and $$f_1^-$$, (d) $$f_{1}^+$$ and $$f_1^-$$, and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. RTM is shown in (c), and its combined image with (b) is shown in (e). White arrows indicate artifacts. Images are clipped at 5 per cent of their maximum amplitude. The third method which requires windowing $$f_1^-$$ fields cannot be extended as easily to elastic media. In particular, the assumption that the last arriving wave is a pure-mode primary used for acoustic media does not hold for all $$f_1^-$$ gathers. A focus originating from a compressional time-reversed direct wave will create not only P waves, but also conversions which will arrive at or after the last pure-mode P primary in the upgoing focusing functions: these have certainly reflected only once, but they also forward-scatter through conversion, making them non-compliant with the first-order Born assumption and harder to migrate. However, the assumption is valid for S-wave foci: the last event arriving in the upgoing field created by an S direct wave will also be a pure-mode S wave. Any conversion into P will have traveled faster than the pure S wave. As such, the third method of primaries can be extended directly to pure-mode S waves in elastic media. This is verified by Fig. 10(d), which uses SS data to show the recovery of only the reflector directly above the focus points, without any conversion-related artifacts. These three methods and their application to several models provides significant insights into several areas of active research in seismics. In standard RTM, one assumes that the point-source function only creates primaries. This is known as the linearization, or first-order Born approximation, whereby the receiver wavefield created by the source function is linearly related to the medium parameters. In order for this approach to be valid, data used in standard RTM must first have all multiples removed. In our first method, the source function $$f_1^+$$ creates all events in $$f_1^-$$. Indeed, since all events in $$f_1^-$$ are primaries, this imaging method is truly linear in the sense that events in $$f_1^-$$ are linearly related to the source function ($$f_1^+$$) by the reference medium parameters. This is achieved without any multiple removal. The second method, while not respecting wave propagation fully ($$f_{1,d}^+$$ only generates some of the primary energy in $$f_{1}^-$$), has been shown to be a useful approximation in terms of imaging as the approximation made in the physics trades-off particularly with fewer artifacts caused by the migration method itself. The third method also creates clean images in practice, and can be seen as a truly linear method which uses $$f_{1,d}^+$$ as source and only one primary in $$f_{1}^-$$ as the received wavefield. Finally, we have shown new examples of how the combined imaging conditions reduce artifacts in pairs of images, including a successful application to field data. So-called Marchenko imaging requires Green’s and focusing functions to be computed at all imaging locations. This requires the initial focusing field to be computed from each source/receiver location to the subsurface. In addition, it requires that we solve the Marchenko equation separately at all subsurface imaging locations. The cost of obtaining the image is then applying an imaging condition to these fields, the cost of which is much lower than the aforementioned steps. If the initial focusing field is computed with finite differences and stored for all sources and subsurface locations, the cost of computing these fields is the same as computing the source RTM wavefields. The solution to the Marchenko equation relies on computing cross-correlations and convolutions with the reflection data, a process which supplants the receiver-side backpropagation in RTM. This step is usually significantly more costly than RTM, and one cross-correlation/convolution is taken to be equivalent to receiver-side backpropagation (for all source locations). Therefore, Marchenko imaging methods are generally considered to be several times more costly than conventional RTM (Slob et al.2014). Focusing functions, on the other hand, allow structures to be targeted with a small number of focus points. This means that the initial focusing fields and Marcheko solutions need to be computed at fewer locations. For example, in the 2-D synthetic example we computed these fields at 151 focusing locations, while a Marchenko image would require 301151 = 751 × 401 of these computations corresponding to each image location in the model. After focusing functions have been computed, they need to be used in the forward and backwards extrapolation steps like RTM, which dominates the cost. In the 2-D synthetic, we performed 151 of these computations for our focusing images versus 201 for standard RTM; in the field data, these numbers were 411 and 235, respectively. Therefore, whether the focusing images or RTM will be more computationally intensive depends on the application. At the cost of performing both the standard RTM and our focusing RTM, the combined images can deliver very clean images. As such, these methods can be used in a target oriented way to monitor the evolution of producing reservoirs over time. Moreover, the first two methods can be used in place of demultiple plus RTM or Marchenko redatuming plus RTM, especially for geologies which generate many multiples potentially driving down the total processing cost. 6 CONCLUSION We have presented novel methods to image subsurface structures using specially crafted virtual acquisitions. These acquisitions are obtained by the outputs of the Marchenko method, and are composed of source/receiver wavefields given by the down- and upgoing focusing functions, respectively. Additionally, we present a method based on these acquisitions which can image individual reflectors: this method involves windowing upgoing focusing functions based on generally applicable properties of the traveltimes of events in focusing functions, making it feasible to perform on field data. We apply these methods to field data from the North Sea, to acoustic 1-D and 2-D synthetic data, and to elastic 1.5-D synthetic data. We discuss the advantages and limitations of the methods, as well as implications they might have for other areas of seismic processing. Acknowledgements The authors thank the Edinburgh Interferometry Project sponsors (Schlumberger Gould Research, Statoil and Total) for supporting this research. We thank Statoil ASA, the Statoil Volve team and the Volve license partners ExxonMobil E&P Norway and Bayerngas Norge for release of the data. The first author also thanks CAPES for funding. The authors also thank the Editor, Lapo Boschi, as well as reviewer Filippo Broggini and an anonymous reviewer for their contributions. REFERENCES Amundsen L., 2001. Elimination of free-surface related multiples without need of the source wavelet, Geophysics , 66( 1), 327– 341. https://doi.org/10.1190/1.1444912 Google Scholar CrossRef Search ADS   Bakulin A., Calvert R., 2006. The virtual source method: theory and case study, Geophysics , 71( 4), SI139– SI150. https://doi.org/10.1190/1.2216190 Google Scholar CrossRef Search ADS   Baysal E., Kosloff D.D., Sherwood J. W.C., 1983. Reverse time migration, Geophysics , 48( 11), 1514– 1524. https://doi.org/10.1190/1.1441434 Google Scholar CrossRef Search ADS   Behura J., Wapenaar K., Snieder R., 2014. Autofocus imaging: image reconstruction based on inverse scattering theory, Geophysics , 79( 3), A19– A26. https://doi.org/10.1190/geo2013-0398.1 Google Scholar CrossRef Search ADS   Berkhout A.J., 2012. Combining full wavefield migration and full waveform inversion, a glance into the future of seismic imaging, Geophysics , 77( 2), S43– S50. https://doi.org/10.1190/geo2011-0148.1 Google Scholar CrossRef Search ADS   Berkhout A.J., Verschuur D.J., 2005. Removal of internal multiples with the common-focus-point (CFP) approach: Part 1—explanation of the theory, Geophysics , 70( 3), V45– V60. https://doi.org/10.1190/1.1925753 Google Scholar CrossRef Search ADS   Berkhout A.J., Verschuur D.J., 2006. Imaging of multiple reflections, Geophysics , 71( 4), SI209– SI220. https://doi.org/10.1190/1.2215359 Google Scholar CrossRef Search ADS   Berryhill J.R., 1979. Wave-equation datuming, Geophysics , 44( 8), 1329– 1344. https://doi.org/10.1190/1.1441010 Google Scholar CrossRef Search ADS   Berryhill J.R., 1984. Wave-equation datuming before stack, Geophysics , 49( 11), 2064– 2066. https://doi.org/10.1190/1.1441620 Google Scholar CrossRef Search ADS   Broggini F., Wapenaar K., van der Neut J., Snieder R., 2014. Data-driven Green’s function retrieval and application to imaging with multidimensional deconvolution, J. geophys. Res. , 119( 1), 425– 441. https://doi.org/10.1002/2013JB010544 Google Scholar CrossRef Search ADS   Chang W.-F., 1986. Reverse-time migration of offset vertical seismic profiling data using the excitation-time imaging condition, Geophysics , 51( 1), 67, doi:10.1190/1.1442041. https://doi.org/10.1190/1.1442041 Google Scholar CrossRef Search ADS   Curtis A., Halliday D., 2010. Source-receiver wave field interferometry, Phys. Rev. E , 81( 4), 46601. https://doi.org/10.1103/PhysRevE.81.046601 Google Scholar CrossRef Search ADS   Curtis A., Gerstoft P., Sato H., Snieder R., Wapenaar K., 2006. Seismic interferometry—turning noise into signal, Leading Edge , 25( 9), 1082– 1092. https://doi.org/10.1190/1.2349814 Google Scholar CrossRef Search ADS   da Costa Filho C.A., Curtis A., 2016. Attenuating multiple-related imaging artifacts using combined imaging conditions, Geophysics , 81( 6), S469– S475. https://doi.org/10.1190/geo2016-0113.1 Google Scholar CrossRef Search ADS   da Costa Filho C.A., Ravasi M., Curtis A., Meles G.A., 2014. Elastodynamic Green’s function retrieval through single-sided Marchenko inverse scattering, Phys. Rev. E , 90( 6), 063201. https://doi.org/10.1103/PhysRevE.90.063201 Google Scholar CrossRef Search ADS   da Costa Filho C.A., Ravasi M., Curtis A., 2015. Elastic P- and S-wave autofocus imaging with primaries and internal multiples, Geophysics , 80( 5), S187– S202. https://doi.org/10.1190/geo2014-0512.1 Google Scholar CrossRef Search ADS   da Costa Filho C.A., Meles G.A., Curtis A., 2017. Elastic internal multiple analysis and attenuation using Marchenko and interferometric methods, Geophysics , 82( 2), Q1– Q12. https://doi.org/10.1190/geo2016-0162.1 Google Scholar CrossRef Search ADS   Fink M., 1992. Time reversal of ultrasonic fields. I. Basic principles., IEEE Trans. Ultrason. Ferroelectr. Freq. Control , 39( 5), 555– 66. https://doi.org/10.1109/58.156174 Google Scholar CrossRef Search ADS PubMed  Halliday D., Curtis A., 2010. An interferometric theory of source-receiver scattering and imaging, Geophysics , 75( 6), SA95– SA103. https://doi.org/10.1190/1.3486453 Google Scholar CrossRef Search ADS   Jia X., Guitton A., Singh S., Snieder R., 2017. Subsalt Marchenko imaging: a Gulf of Mexico example, in SEG Technical Program Expanded Abstracts 2017 , pp. 5588– 5592, Society of Exploration Geophysicists. Jones I.F., 2010. An Introduction To: Velocity Model Building , EAGE Publications bv. Google Scholar CrossRef Search ADS   Kinneging N.A., Budejicky V., Wapenaar C. P.A., Berkhout A.J., 1989. Efficient 2D and 3D shot record redatuming, Geophys. Prospect. , 37( 5), 493– 530. https://doi.org/10.1111/j.1365-2478.1989.tb02220.x Google Scholar CrossRef Search ADS   Malcolm A.E., Ursin B., de Hoop M.V., 2009. Seismic imaging and illumination with internal multiples, Geophys. J. Int. , 176( 3), 847– 864. https://doi.org/10.1111/j.1365-246X.2008.03992.x Google Scholar CrossRef Search ADS   Masson Y., Romanowicz B., 2017. Box tomography: localized imaging of remote targets buried in an unknown medium, a step forward for understanding key structures in the deep Earth, Geophys. J. Int. , 211( 1), 141– 163. https://doi.org/10.1093/gji/ggx141 Google Scholar CrossRef Search ADS   McMechan G.A., 1983. Migration by extrapolation of time–dependent boundary values, Geophys. Prospect. , 31( 3), 413– 420. https://doi.org/10.1111/j.1365-2478.1983.tb01060.x Google Scholar CrossRef Search ADS   Meles G.A., Löer K., Ravasi M., Curtis A., da Costa Filho C.A., 2015. Internal multiple prediction and removal using Marchenko autofocusing and seismic interferometry, Geophysics , 80( 1), A7– A11. https://doi.org/10.1190/geo2014-0408.1 Google Scholar CrossRef Search ADS   Meles G.A., Wapenaar K., Curtis A., 2016. Reconstructing the primary reflections in seismic data by Marchenko redatuming and convolutional interferometry, Geophysics , 81( 2), Q15– Q26. https://doi.org/10.1190/geo2015-0377.1 Google Scholar CrossRef Search ADS   Meles G.A., da Costa Filho C.A., Curtis A., 2017. Estimating and imaging with primaries constructed from single-sided focusing wavefields, in 79th EAGE Conference and Exhibition , doi:10.3997/2214-4609.201701056. Muijs R., Robertsson J. O.A., Holliger K., 2007. Prestack depth migration of primary and surface-related multiple reflections: Part I - Imaging, Geophysics , 72, S59, doi:10.1190/1.2422796. https://doi.org/10.1190/1.2422796 Google Scholar CrossRef Search ADS   Oristaglio M.L., 1989. An inverse scattering formula that uses all the data, Inverse Probl. , 5( 6), 1097– 1105. https://doi.org/10.1088/0266-5611/5/6/015 Google Scholar CrossRef Search ADS   Ravasi M., Vasconcelos I., Kritski A., Curtis A., da Costa Filho C.A., Meles G.A., 2015. Marchenko imaging of Volve field, North Sea, in 77th EAGE Meeting and Technical Exhibition 2015 , doi:10.3997/2214-4609.201412938. Ravasi M., Vasconcelos I., Kritski A., Curtis A., da Costa Filho C.A., Meles G.A., 2016. Target-oriented Marchenko imaging of a North Sea field, Geophys. J. Int. , 205( 1), 99– 104. https://doi.org/10.1093/gji/ggv528 Google Scholar CrossRef Search ADS   Rose J.H., 2001. “Single-sided” focusing of the time-dependent Schrödinger equation, Phys. Rev. A , 65( 1), 0127071. https://doi.org/10.1103/PhysRevA.65.012707 Google Scholar CrossRef Search ADS   Schuster G.T., Zhou M., 2006. A theoretical overview of model-based and correlation-based redatuming methods, Geophysics , 71( 4), SI103– SI110. https://doi.org/10.1190/1.2208967 Google Scholar CrossRef Search ADS   Schuster G.T., Yu J., Sheng J., Rickett J., 2004. Interferometric/daylight seismic imaging, Geophys. J. Int. , 157( 2), 838– 852. https://doi.org/10.1111/j.1365-246X.2004.02251.x Google Scholar CrossRef Search ADS   Singh S., Snieder R., Behura J., van der Neut J., 2015. Marchenko imaging: imaging with primaries, internal multiples, and free-surface multiples, Geophysics , 80( 5), S164– S174. https://doi.org/10.1190/geo2014-0494.1 Google Scholar CrossRef Search ADS   Slob E., Wapenaar K., Broggini F., Snieder R., 2014. Seismic reflector imaging using internal multiples with Marchenko-type equations, Geophysics , 79( 2), S63– S76. https://doi.org/10.1190/geo2013-0095.1 Google Scholar CrossRef Search ADS   Staring M., Pereira R., Douma H., van der Neut J., Wapenaar C., 2017. Adaptive double-focusing method for source-receiver Marchenko redatuming on field data, in SEG Technical Program Expanded Abstracts 2017 , pp. 4808– 4812, Society of Exploration Geophysicists. Tang Y., Biondi B., 2011. Target-oriented wavefield tomography using synthesized Born data, Geophysics , 76( 5), WB191– WB207. https://doi.org/10.1190/geo2010-0383.1 Google Scholar CrossRef Search ADS   Thorbecke J., 1997. Common focus point technology, PhD thesis , Delft University of Technology. van der Neut J., Wapenaar K., 2016. Adaptive overburden elimination with the multidimensional Marchenko equation, Geophysics , 81( 5), T265– T284. https://doi.org/10.1190/geo2016-0024.1 Google Scholar CrossRef Search ADS   Vasconcelos I., Wapenaar K., van der Neut J., Thomson C., Ravasi M., 2015. Using inverse transmission matrices for Marchenko redatuming in highly complex media, in SEG Technical Program Expanded Abstracts 2015 , pp. 5081– 5086, Society of Exploration Geophysicists. Verschuur D.J., 1992. Adaptive surface-related multiple elimination, Geophysics , 57( 9), 1166– 1177. https://doi.org/10.1190/1.1443330 Google Scholar CrossRef Search ADS   Wapenaar C. P.A., Kinneging N., Berkhout A., 1987. Principle of prestack migration based on the full elastic two-way wave equation, Geophysics , 52( 2), 151– 173. https://doi.org/10.1190/1.1442291 Google Scholar CrossRef Search ADS   Wapenaar K., 2014. Single-sided Marchenko focusing of compressional and shear waves, Phys. Rev. E , 90( 6), 063202. https://doi.org/10.1103/PhysRevE.90.063202 Google Scholar CrossRef Search ADS   Wapenaar K., Slob E., 2014. On the Marchenko equation for multicomponent single-sided reflection data, Geophys. J. Int. , 199( 3), 1367– 1371. https://doi.org/10.1093/gji/ggu313 Google Scholar CrossRef Search ADS   Wapenaar K., Broggini F., Slob E., Snieder R., 2013. Three-dimensional single-sided Marchenko inverse scattering, data-driven focusing, Green’s function retrieval, and their mutual relations, Phys. Rev. Lett. , 110( 8), 084301. https://doi.org/10.1103/PhysRevLett.110.084301 Google Scholar CrossRef Search ADS PubMed  Wapenaar K., Thorbecke J., van der Neut J., Broggini F., Slob E., Snieder R., 2014. Marchenko imaging, Geophysics , 79( 3), WA39– WA57. https://doi.org/10.1190/geo2013-0302.1 Google Scholar CrossRef Search ADS   Wapenaar K., van der Neut J., Slob E., 2016. Unified double- and single-sided homogeneous Green’s function representations, Proc. R. Soc. A: Math. Phys. Eng. Sci. , 472( 2190), 20160162. https://doi.org/10.1098/rspa.2016.0162 Google Scholar CrossRef Search ADS   Weglein A.B., Gasparotto F.A., Carvalho P.M., Stolt R.H., 1997. An inverse-scattering series method for attenuating multiples in seismic reflection data, Geophysics , 62( 6), 1975– 1989. https://doi.org/10.1190/1.1444298 Google Scholar CrossRef Search ADS   Yuan S., Fuji N., Singh S., Borisov D., 2017. Localized time-lapse elastic waveform inversion using wavefield injection and extrapolation: 2-D parametric studies, Geophys. J. Int. , 209( 3), 1699– 1717. https://doi.org/10.1093/gji/ggx118 Google Scholar CrossRef Search ADS   Ziolkowski A., Taylor D.B., Johnston R. G.K., 1999. Marine seismic wavefield measurement to remove sea-surface multiples, Geophys. Prospect. , 47( 6), 841– 870. https://doi.org/10.1046/j.1365-2478.1999.00165.x Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Imaging strategies using focusing functions with applications to a North Sea field

, Volume 213 (1) – Apr 1, 2018
13 pages

/lp/ou_press/imaging-strategies-using-focusing-functions-with-applications-to-a-BI42XClyEt
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx562
Publisher site
See Article on Publisher Site

### Abstract

Summary Seismic methods are used in a wide variety of contexts to investigate subsurface Earth structures, and to explore and monitor resources and waste-storage reservoirs in the upper ∼100 km of the Earth’s subsurface. Reverse-time migration (RTM) is one widely used seismic method which constructs high-frequency images of subsurface structures. Unfortunately, RTM has certain disadvantages shared with other conventional single-scattering-based methods, such as not being able to correctly migrate multiply scattered arrivals. In principle, the recently developed Marchenko methods can be used to migrate all orders of multiples correctly. In practice however, using Marchenko methods are costlier to compute than RTM—for a single imaging location, the cost of performing the Marchenko method is several times that of standard RTM, and performing RTM itself requires dedicated use of some of the largest computers in the world for individual data sets. A different imaging strategy is therefore required. We propose a new set of imaging methods which use so-called focusing functions to obtain images with few artifacts from multiply scattered waves, while greatly reducing the number of points across the image at which the Marchenko method need be applied. Focusing functions are outputs of the Marchenko scheme: they are solutions of wave equations that focus in time and space at particular surface or subsurface locations. However, they are mathematical rather than physical entities, being defined only in reference media that equal to the true Earth above their focusing depths but are homogeneous below. Here, we use these focusing functions as virtual source/receiver surface seismic surveys, the upgoing focusing function being the virtual received wavefield that is created when the downgoing focusing function acts as a spatially distributed source. These source/receiver wavefields are used in three imaging schemes: one allows specific individual reflectors to be selected and imaged. The other two schemes provide either targeted or complete images with distinct advantages over current RTM methods, such as fewer artifacts and artifacts that occur in different locations. The latter property allows the recently published ‘combined imaging’ method to remove almost all artifacts. We show several examples to demonstrate the methods: acoustic 1-D and 2-D synthetic examples, and a 2-D line from an ocean bottom cable field data set. We discuss an extension to elastic media, which is illustrated by a 1.5-D elastic synthetic example. Interferometry, Body waves, Controlled source seismology, Wave propagation, Wave scattering and diffraction 1 INTRODUCTION Seismic methods using actively fired sources are used in a wide variety of contexts to explore the upper ∼100 km of the Earth’s subsurface. This is important both beneath land and offshore, in order to identify and monitor Earth resources (ores, minerals, hydrocarbons and water reservoirs), potential storage reservoirs in host rock for waste or by-products (nuclear waste or CO2), suitable bedrock for foundations of buildings and other infrastructure, and for examining deeper structures such as the Earth’s lower crust and upper mantle. The family of seismic imaging methods called migration are an integral part of applied seismology, and are used more broadly in medical imaging and non-destructive evaluation of components and structures. Within geophysics, they construct high-resolution images of subsurface structures from recorded seismic survey data. In a standard seismic data acquisition, seismic waves are created by surface sources (e.g. air guns and vibroseis), which propagate through the subsurface and are recorded at an array of receivers (e.g. hydrophones and geophones). Using state-of-the-art methods, images are usually constructed from standard seismic data using so-called reverse-time migration (RTM, Baysal et al.1983; McMechan 1983; Chang 1986). In RTM, a synthetic source mimicking the real source is propagated computationally through a mainly smooth model which approximates the subsurface seismic velocities with estimates that can be obtained from other methods such as velocity analysis, tomography or waveform inversion (Jones 2010). The received wavefield is computationally backpropagated (propagated backwards in time) through this same smooth model. At points in space where these two wavefields coincide over significant periods of time, it is likely that a subsurface scatterer exists. This is because at such locations in the physical seismic experiment, the source field converts into the receiver-side (scattered) field, hence the fact that source and receiver fields match. A common approach is therefore to calculate the zero-time lag cross-correlation (a measure of similarity) between these two numerically estimated wavefields at every point in our computational grid; the set of such values constitutes the final image. In common with most other seismic processing and migration methods, RTM relies on the single-scattering (or first-order Born) assumption (Oristaglio 1989). Seismic data always contain multiply scattered waves (herein called multiples), and these create spurious structures in the final image since within RTM they cannot be backpropagated from the receivers to each correct scattering location without already having knowledge of the other scattering locations. Multiples must therefore be removed from seismic data prior to migration. Several commonly used methods exist to suppress free surface-related and internal multiples from pre-stack data, as described for example in Verschuur (1992), Weglein et al. (1997), Ziolkowski et al. (1999) and Amundsen (2001). Other methods alter the principles of single-scattering-based imaging to account for multiples during the migration process, and may even use them to enhance the final image (Schuster et al.2004; Berkhout & Verschuur 2006; Muijs et al.2007; Malcolm et al.2009; Berkhout 2012). One particular recently developed redatuming method, the Marchenko method (Wapenaar et al.2013), can be used to obtain images of the subsurface devoid of most multiple-related artifacts. It works by solving a so-called Marchenko equation that relates surface reflection data, transmission estimates (waves traveling directly from the surface to each imaging location in the subsurface), and the critical wavefields for constructing images which are the surface-to-subsurface Green’s functions and focusing functions. The Green’s function is the full-wavefield response at a subsurface location due to impulsive sources at the surface, including all orders of scattering. RTM approximates these Green’s functions in the propagation step described above by omitting all scattering since the propagation model is always dominantly smooth. Green’s functions estimated using Marchenko methods, on the other hand, include all orders of scattering. These Green’s functions could therefore be used in place of the propagation steps in RTM. The other important component of Marchenko methods are focusing functions. These are solutions of the wave equation which propagate through a reference medium in such a way that they reach a chosen focusing depth as a localized delta function in time and space at each focus point. The reference medium in which a focusing function is defined contains all true reflectors above its focusing depth but is homogeneous below that depth. The focused field therefore then diverges as a downgoing field through the homogeneous part of the reference medium, and the focusing functions are mathematical entities that cannot be realized in the true Earth. Nevertheless since they account for the true Earth down to the depth of focus, they clearly contain information about the Earth, a fact that we exploit below. Marchenko methods are part of a larger class of high-frequency redatuming methods, which include classical redatuming (Berryhill 1979, 1984; Kinneging et al.1989), and interferometric redatuming (Curtis et al.2006; Bakulin & Calvert 2006; Schuster & Zhou 2006; Curtis & Halliday 2010; Halliday & Curtis 2010). Accurate low-frequency models of the subsurface (in particular, accurate one-way traveltimes) are required for a many high-frequency redatuming methods (excluding, most notably, seismic interferometry), and these have been traditionally been supplied by inversions for the full subsurface model. However, there have been important advances in lower frequency redatuming such as target-oriented inversion (Tang & Biondi 2011), localized full-waveform inversion (Yuan et al.2017) and box tomography (Masson & Romanowicz 2017), which may be beneficial for target-oriented applications by providing more detailed models of subsurface target areas. Marchenko methods have recently been used for a variety of purposes, including internal multiple attenuation (Meles et al.2015; da Costa Filho et al.2017), target-oriented imaging (Wapenaar et al.2014; Ravasi et al.2016), overburden elimination (van der Neut & Wapenaar 2016), synthesizing primaries (Meles et al.2016, 2017) and virtual subsurface wavefield estimation (Wapenaar et al.2016). It has also been extensively used for imaging (e.g. Wapenaar et al.2014; Singh et al.2015; da Costa Filho et al.2015). So far, most Marchenko imaging methods have required computing the solution to the Marchenko equation at each imaging point—a costly endeavor. Alternatively, they have been used to invert for a redatumed reflection response devoid of overburden interactions which is the used in conventional imaging (Wapenaar et al.2014; Ravasi et al.2016; Staring et al.2017). Herein, we show how Marchenko focusing functions computed at far fewer locations can be used to obtain significantly better images than RTM. This is achieved by substituting source and receiver wavefields in conventional RTM by focusing functions. We identify three different source/receiver wavefield pairs which give rise to three distinct imaging methods, all of which account for multiple scattering in the Earth’s subsurface. Finally, we propose a post-imaging filter based on da Costa Filho & Curtis (2016) which can be used to enhance the results of two of these methods. We compare all methods to conventional RTM in 1-D and 2-D acoustic synthetic media, as well as in a test involving real data from a 2-D line from an ocean bottom cable (OBC) field survey in the North Sea. Finally, we discuss extensions to elastic media using a 1.5-D elastic synthetic medium in a worked example. 2 WAVE THEORY The idea of using focusing operators to better image a medium is not new. It has been used in several different contexts, including acoustic time-reversal focusing experiments (Fink 1992), as well as migration and inversion of seismic data (Thorbecke 1997; Berkhout & Verschuur 2005). The main practical constraint imposed in geophysics is that data are usually only collect on one side of a medium (the Earth’s surface) in contrast to medical imaging which the subject can be entirely surrounded by sources and receivers. Rose (2001) published the first method to obtain focusing wavefields using single-sided data in a medium with several scatterers. This was generalized to provide 3-D focusing functions in 3-D acoustic media through the work of Wapenaar et al. (2013), and to 3-D elastic media by da Costa Filho et al. (2014) and Wapenaar & Slob (2014). Similarly to Wapenaar et al. (2014), in our notation, the Marchenko method provides Green’s functions G(xF, x0, t) as well as focusing functions denoted f1(x0, xF, t) where $$\mathbf {x}_0 = (\hat{\mathbf {x}}_0, z_0)$$ is a location at the surface z = z0, $$\mathbf {x}_F = (\hat{\mathbf {x}}_F, z_F)$$ is the focusing location at depth z = zF and t is time. Here, $$\hat{\mathbf {x}}$$ represents the horizontal coordinates of x. A focusing function is defined in a reference medium which is homogeneous below zF, and is equal to the true (unknown) medium between the surface and the focusing depth, as shown in Fig. 1. They have the characteristic that at zero time, the wavefield observed at the depth of zF focuses: that is, it becomes a delta function at $$\hat{\mathbf {x}}_F$$. These functions are commonly decomposed with regard to their directions of propagation related to their first coordinate (on the surface at which they are emitted or received). Therefore, with a focus at xF, $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$ refers to the downgoing focusing field departing from x0, and $$f_1^-(\mathbf {x}_0, \mathbf {x}_F, t)$$ is the upgoing field which then arrives at x0. The superscript + is inherited from the convention that our z-axis is positive downwards. Similarly, G − (xF, x0, t) and G + (xF, x0, t) are defined as up- and downgoing Green’s functions, respectively, where now the superscripts + and − refer to the direction of propagation at the subsurface focus point xF. Figure 1. View largeDownload slide Reference medium in which focusing functions $$f_1^\pm (\mathbf {x}_0, \mathbf {x}_F)$$ are defined. Figure 1. View largeDownload slide Reference medium in which focusing functions $$f_1^\pm (\mathbf {x}_0, \mathbf {x}_F)$$ are defined. Focusing fields along a horizontal line at depth have been used for imaging structures from below (Wapenaar et al.2014; Ravasi et al.2016). This shows that focusing functions contain information about the reflection properties of the medium. With an alternative interpretation of these fields, we now show other ways in which they can be used for imaging. The focusing functions are related to the surface reflection data $$R(\mathbf {x}_0^{\prime \prime }, \mathbf {x}_0, t)$$ and to the subsurface Green’s function $$G^-(\mathbf {x}_F, \mathbf {x}_0^{\prime \prime }, t)$$ by the following equation:   \begin{eqnarray} G^-(\mathbf {x}_F, \mathbf {x}_0^{\prime \prime }, t) &=& - f_1^-(\mathbf {x}_0^{\prime \prime }, \mathbf {x}_F, t) \nonumber\\ &&+ \int \limits _{\partial \mathbb {D}_0}\int \limits _{-\infty }^\infty R(\mathbf {x}_0^{\prime \prime }, \mathbf {x}_0, t-\tau ) f^+_1(\mathbf {x}_0, \mathbf {x}_F, \tau )\,\mathrm{d}{\tau }\,\mathrm{d}^2{\mathbf {x}_0},\nonumber\\ \end{eqnarray} (1)where $$R(\mathbf {x}_0^{\prime \prime }, \mathbf {x}_0, t)$$ corresponds to the pressure of the reflected field recorded at $$\mathbf {x}_0^{\prime \prime }$$ created by vertical force sources at x0, multiplied by −2 (Broggini et al.2014). Eq. (1) is obtained by considering two wavestates A and B: wavestate A corresponds to Green’s functions which exist in the true medium, while the focusing functions correspond to wavestate B and exist in a reference medium which is truncated below xF (Fig. 1). For general wavestates which propagate through media that are identical between two arbitrary boundaries $$\partial \mathbb {D}_0$$ and $$\partial \mathbb {D}_F$$, Wapenaar et al. (2014) derive a one-way acoustic reciprocity theorem for pressure normalized wavefields in the angular frequency domain $$p_A^\pm (\mathbf {x}, \omega )$$ and $$p_B^\pm (\mathbf {x}, \omega )$$:   \begin{eqnarray} &&{-\int \limits _{\partial \mathbb {D}_0} \frac{1}{\rho }\lbrace p_A^+ (\partial _z p_B^-) + p_A^- (\partial _z p_B^+)\rbrace \,\mathrm{d}^2{\mathbf {x}_0} }\nonumber\\ &&= \int \limits _{\partial \mathbb {D}_F} \frac{1}{\rho }\lbrace (\partial _z p_A^+)p_B^- + (\partial _z p_A^-) p_B^+\rbrace \,\mathrm{d}^2{\mathbf {x}_F}, \end{eqnarray} (2)where arguments of $$p_A^\pm$$ and $$p_B^\pm$$ have been omitted, and pressure normalized fields are those whose sum equals that of the pressure, that is p+ + p− = p. Substituting appropriate boundary conditions in eq. (2) yields eq. (1) (see appendix A of Wapenaar et al.2014). By choosing different wavestates, as long as they are defined in media which coincide between two arbitrary boundaries $$\partial \mathbb {D}_0$$ and $$\partial \mathbb {D}_F$$, we may use the same one-way acoustic reciprocity theorem in eq. (2) to obtain novel relationships. We use this freedom to consider wavestate A to refer to focusing wavefields in the reference medium as before, but now consider wavestate B to be the Green’s functions $$\overline{G}$$ also in the reference medium, as opposed to the true medium which was assumed when deriving eq. (1). More explicitly, we substitute $$p_A^\pm (\mathbf {x}, \omega ) = f_1^\pm (\mathbf {x}, \mathbf {x}_F^{\prime }, \omega )$$ and $$p_B^\pm (\mathbf {x}, \omega ) = \overline{G}^\pm (\mathbf {x}, \mathbf {x}_0^{\prime }, \omega )$$ in eq. (2) to obtain   \begin{eqnarray} &&{-\!\int \limits _{\partial \mathbb {D}_0} \!\frac{1}{\rho }\lbrace f_1^+(\mathbf {x}_0, \mathbf {x}_F^{\prime }) [\partial _z \overline{G}^-(\mathbf {x}_0,\mathbf {x}_0^{\prime })] \!+\! f_1^-(\mathbf {x}_0, \mathbf {x}_F^{\prime }) [\partial _z \overline{G}^+(\mathbf {x}_0,\mathbf {x}_0^{\prime })]\rbrace \mathrm{d}^2{\mathbf {x}_0}}\nonumber \\ && = \int \limits _{\partial \mathbb {D}_F} \frac{1}{\rho }\lbrace [\partial _z f_1^+(\mathbf {x}_F, \mathbf {x}_F^{\prime })] \overline{G}^-(\mathbf {x}_F, \mathbf {x}_0^{\prime }) \!+\! [\partial _z f_1^-(\mathbf {x}_F, \mathbf {x}_F^{\prime })]\nonumber\\ &&\quad\times \overline{G}^+(\mathbf {x}_F, \mathbf {x}_0^{\prime })\rbrace \,\mathrm{d}^2{\mathbf {x}_F}, \end{eqnarray} (3)where the frequency dependence has been suppressed for notational brevity. A property of the reference medium used for both wavestates is that no reflections come from above $$\partial \mathbb {D}_0$$ or from below $$\partial \mathbb {D}_F$$. This condition is such that the vertical derivative $$\partial _z\overline{G}^+(\mathbf {x}_0, \mathbf {x}_0^{\prime })$$ vanishes on $$\partial \mathbb {D}_0$$, and both $$\overline{G}^-(\mathbf {x}_F, \mathbf {x}_0^{\prime })$$ and $$\partial _z f_1^-(\mathbf {x}_F, \mathbf {x}_F^{\prime })$$ also vanish on $$\partial \mathbb {D}_F$$. In addition, following previous Marchenko work, we assume that $$\partial _z \overline{G}^+(\mathbf {x}_0, \mathbf {x}_0^{\prime }) = -\frac{1}{2} \iota \omega \rho (\mathbf {x}_0^{\prime })\delta (z_0 -z_0^{\prime })$$ at $$\partial \mathbb {D}_0$$, that is, Green’s functions have a vertical force point-source mechanism (and that is the only downgoing field at the surface since there are no reflectors above that level). Applying these conditions to eq. (3) yields   $$\int \limits _{\partial \mathbb {D}_0} \overline{R}(\mathbf {x}_0^{\prime }, \mathbf {x}_0)f_1^+(\mathbf {x}_0, \mathbf {x}_F^{\prime }) \,\mathrm{d}^2{\mathbf {x}_0} = f_1^-(\mathbf {x}_0^{\prime }, \mathbf {x}_F^{\prime }),$$ (4)where $$\overline{R}(\mathbf {x}_0^{\prime }, \mathbf {x}_0) \equiv \frac{-2}{\iota \omega \rho (\mathbf {x}_0)}\partial _z \overline{G}^-(\mathbf {x}_0,\mathbf {x}_0^{\prime })$$. If we define the volume encompassed between depths ε above and below $$\partial \mathbb {D}_0$$ as V (Fig. 1), we may rewrite eq. (4) as   \begin{eqnarray} \int \limits _V \overline{R}(\mathbf {x}_0^{\prime }, \mathbf {x}_0)f_1^+(\mathbf {x}_0, \mathbf {x}_F^{\prime })\delta (z-z_0) \mathrm{d}^3{\mathbf {x}_0} \!=\! \int \limits _{-\varepsilon }^{\varepsilon } f_1^-(\mathbf {x}_0^{\prime }, \mathbf {x}_F^{\prime })\delta (z-z_0)\mathrm{d}{z},\nonumber\\ \end{eqnarray} (5)where z0 is the depth coordinate of x0. Considering that similarly to R, $$\overline{R}$$ may be viewed as a scaled dipole-to-monopole surface response in the reference medium multiplied by −2, it is effectively a Green’s function. As such, the left-hand side of eq. (5) is a volume integral between a Green’s function and a source function, namely $$f_1^+(\mathbf {x}_0, \mathbf {x}_F^{\prime })\delta (z-z_0)$$. A direct consequence of this is that a line source of $$f_1^+$$ produces an $$f_1^-$$ field localized along $$\partial \mathbb {D}_0$$. Simply put, $$f_1^-$$ is the received reference medium response to source field $$f_1^+$$, where both the action of the source field and the recorded response are at the surface. The Marchenko method is then not only useful for producing redatumed Green’s functions, but also focusing functions which are essentially virtual data acquisitions (surveys) that take place in an imagined medium which contains the true medium above the focusing depth, that is, the reference medium in Fig. 1. Here, we explore some of interesting properties of these virtual surveys which are derived from the focusing functions, and propose imaging methods which exploit these virtual acquisitions, in particular an imaging method that targets specific reflectors. 3 IMAGING METHODS In RTM, each seismic point source is injected into a computational model of the medium that is a best-possible estimate of the true earth subsurface, and is propagated forwards in time to simulate the downgoing field from the source throughout the subsurface. The received wavefield is injected into the same model and propagated backwards in time to estimate what this field was like throughout the subsurface before it was recorded. Conveniently, this second step can be applied using a standard forward propagation simulation of the time-reversed input data because the wave equation is symmetric in time (forward and backward propagations are identical mathematical operations). The two wavefields are then combined to construct an image, the particular manner in which the fields are combined being called the imaging condition. When imaging with a virtual survey, the same process can be employed. If we use $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$ as the source field, it must be injected along the whole surface source array simultaneously (as opposed to being a point source at a specific location as in conventional RTM). In this case, the receiver array measures $$f_1^-(\mathbf {x}_0, \mathbf {x}_F, t)$$ as opposed to the recorded surface data, hence $$f_1^-$$ must be injected along the receiver array and backpropagated in time. The nature of the focusing source/receiver pair will have two effects on the image produced using a single virtual focus point xF. First, it will be more localized than an outward-spreading field from a point source as the new source focuses towards the focus point. While the outward-spreading wavefield will tend to reach all subsurface locations given sufficient time to propagate, the focusing source will be approximately contained in the triangle (or downwards-pointing cone, in 3-D) whose base is the source array on it and whose tip is the focus point. Second, since it is defined in the reference medium, it cannot be used to image anything below xF. In light of this, in order to image the medium using these virtual surveys, one has to compute the image for several focus points in order to illuminate the whole subsurface, and these foci must be located below the deepest reflector of interest. The velocity model used to image these virtual surveys is essentially the same as used in conventional imaging, with the difference that they need extend only down to the virtual source depth. The nature of the focusing functions has been explored and exemplified in Slob et al. (2014) and Wapenaar et al. (2014). They show that $$f_1^+(\mathbf {x}, \mathbf {x}_F,t)$$ contains a time-reversed direct wave $$f_{1,d}^+$$ which scatters as it propagates through the reference medium. A coda of multiply scattered waves, $$f_{1,m}^+$$, destructively interferes with the scattering of $$f_{1,d}^+$$ in such way that the only arrival reaching the focusing depth is the time-reversed direct wave as a delta function localized at xF. Note that as $$f_{1,d}^+$$ propagates through the medium it will reflect off every interface in the reference medium; therefore $$f_1^-$$ will contain, amongst other events, primaries resulting from source field $$f_{1,d}^+$$. Above we suggested to migrate source/receiver pair $$f_1^+$$/$$f_1^-$$ from several focus points. However, knowing that $$f_1^-$$ contains a primary from the direct wave $$f_{1,d}^+$$ suggests that results may improve when using $$f^+_{1,d}$$ instead of $$f_1^+$$ as source-side wavefield. This is because RTM assumes that fields consist only of primaries which can be downwards propagated using direct waves, an assumption which is valid when using $$f_{1,d}^+$$, but not when using $$f_1^+$$. Indeed, it has been shown that using a single-event source field (e.g. a direct wave) in the imaging condition eliminates some crosstalk between source-side wavefields and unrelated events in the receiver-side wavefield (Wapenaar et al.1987; da Costa Filho et al.2015). This is corroborated by other recent work: focusing functions have the property of focusing not only in space (at a single subsurface location), but also in time (nothing arrives at the focusing depth before or after the focusing pulse). Meles et al. (2017) have shown that a direct consequence of focusing is that only primaries are present in $$f_1^-$$, at least in horizontally layered media. All multiples are canceled by destructive interference with $$f_{1,m}^+$$ and any residual multiple would destroy the focus. Some of these primaries are created by $$f_{1,d}^+$$, and some of them are created by $$f_{1,m}^+$$. Meles et al. (2017) also showed that in 1-D media the last-arriving event in $$f_1^-$$ is a primary created by $$f_{1,d}^+$$. Therefore, for any given focus point, we may generate an upgoing focusing gather $$f_1^-(\mathbf {x}_0, \mathbf {x}_F, t)$$ such that the last event of this gather will be a primary, $$f_{1,p}^-(\mathbf {x}_0, \mathbf {x}_F, t)$$, originating from a known source, $$f_{1,d}^+(\mathbf {x}_0, \mathbf {x}_F, t)$$. This allows us to migrate only primaries by considering source/receiver pairs composed of $$f_{1,d}^+$$ and $$f_{1,p}^-$$. The primary recovered is the one reflected from the closest reflector immediately above the focusing location, and thus can be used to image that reflector individually. If all reflectors are to be imaged, subsurface focus points must therefore be located below every reflector. We thus derive three new methods to image the subsurface, each consisting or an RTM-type operation using different source/receiver fields. First, the pair $$f_1^+/f_1^-$$ can be used to image the medium above the focus points. Second, we can restrict the source field to the direct wave so as to use only the energy from $$f_{1,d}^+/f_1^-$$ which reduces crosstalk artifacts in the image. The problem with the second method is that some of the primary energy in $$f_1^-$$ was created by the coda $$f_{1,m}^+$$ which is now not included in the source field, hence errors in amplitude may arise and some crosstalk artifacts remain. The third method removes all crosstalk by using pair $$f_{1,d}^+/f_{1,p}^-$$, whose receiver-side field consists of a single primary reflection from the first reflector above the focus point. The third method has the property that it images a single reflector, allowing targeted images to be produced, but has the disadvantage that the last arriving energy in $$f_1^-$$ must be identified and isolated prior to imaging. Another possibility to remove crosstalk present in the first two methods is to use the post-imaging filter proposed by da Costa Filho & Curtis (2016). In that work, two images IC and IF are combined to create a new image using the following expression   $$I_{\rm H}(\mathbf {x}_I) = \frac{E[I_{\rm C}](\mathbf {x}_I)\, I_{\rm F}(\mathbf {x}_I) + E[I_{\rm F}](\mathbf {x}_I)\, I_{\rm C}(\mathbf {x}_I)}{E[I_{\rm C}](\mathbf {x}_I)+E[I_{\rm F}](\mathbf {x}_I)}$$ (6)where xI is the image point, and E[ · ] is the non-negative Hilbert envelope applied to each vertical image trace individually. The combined image IH has the property that events which appear in both constituent images remain in IH, while events which appear in only one or the other of IC or IF are attenuated. Below, we apply the combined images method using as constituent images the conventional RTM image (IC), and a focusing image (IF). In both images, true reflectors appear at the same locations while artifacts may appear in different spatial locations. The combined images can then be used to attenuate those artifacts which appear at different locations. 4 NUMERICAL EXAMPLES We now present several examples to demonstrate and explore these methods. The first two examples are in acoustic media: a 1-D synthetic model with variable velocity, and a 2-D synthetic model with vertical and horizontal density variations. We also apply acoustic methods to data from a 2-D line from an OBC field acquisition in the North Sea. Finally, in Section 5, we provide a simple worked elastic example. 4.1 1-D synthetic First, we present a 1-D synthetic example with focusing functions that have been calculated analytically. The medium and a focusing wavefield are shown in Fig. 2. Events associated with the downwards pointing arrows are injected as part of $$f_1^+(0\,\mathrm{km}, 4.5\,\mathrm{km}, t)$$, with the earliest (leftmost) event being $$f_{1,d}^+$$. The other injected components, $$f_{1,m}^+$$, are seen to contribute no energy to the focus; they are present to destroy events that would otherwise reflect downwards and destroy the focus at depth. It is also seen that the last (rightmost) event arriving at the surface as part of $$f_1^-(0\,\mathrm{km}, 4.5\,\mathrm{km}, t)$$ is a sum of primaries, one of which originates from $$f_{1,d}^+$$. Figure 2. View largeDownload slide Synthetic 1-D model overlain by a focusing wavefield. The black circle denotes the focusing point at depth z = zf and time t = 0, and the dashed line denotes the focusing depth. Arrows pointing downwards depict events in $$f_1^+$$ and those pointing upwards depict events in $$f_1^-$$. Note that since there is only one spatial dimension, the horizontal axis is time. Figure 2. View largeDownload slide Synthetic 1-D model overlain by a focusing wavefield. The black circle denotes the focusing point at depth z = zf and time t = 0, and the dashed line denotes the focusing depth. Arrows pointing downwards depict events in $$f_1^+$$ and those pointing upwards depict events in $$f_1^-$$. Note that since there is only one spatial dimension, the horizontal axis is time. Next, we migrate source/receiver pairs $$f_1^+$$ and $$f_1^-$$, and $$f_{1,d}^+$$ and $$f_1^-$$ using a smooth propagation model, and compare it to conventional RTM, which we will from now on simply refer to as RTM. These images (traces in 1-D examples) are shown, respectively, in Figs 3(a)–(c). All figures have imaged the true reflectors exactly the same way. These appear at incorrect depths (0.62, 1.35 and 3.24 km) because the migration velocity model used a constant velocity of 1 km s−1 in order to eliminate any backscattering which might occur in the imaging step. They also appear with different absolute and relative amplitudes: in the first two images the receiver wavefield is $$f_1^-$$ which has several contributions to its primaries when compared to RTM. Moreover, the first two images have different source-time functions, which again result in different amplitudes. However, this effect is minor: taking the maximum energy of the third reflector reveals peaks at 10.303 × 10−2, 10.304 × 10−2 and 8.964 × 10−2 for Figs 3(a)–(c), respectively. Figure 3. View largeDownload slide Images from sources/receiver wavefields composed from (a) $$f_1^+$$ and $$f_1^-$$, and (b) $$f_{1,d}^+$$ and $$f_1^-$$. (c) Conventional RTM image. (d) Image obtained by combining (a) and (c) using eq. (6). Figure 3. View largeDownload slide Images from sources/receiver wavefields composed from (a) $$f_1^+$$ and $$f_1^-$$, and (b) $$f_{1,d}^+$$ and $$f_1^-$$. (c) Conventional RTM image. (d) Image obtained by combining (a) and (c) using eq. (6). Fig. 3(a) shows two artifacts, one which results from the interaction between $$f_{1,d}^+$$ and an unrelated event in $$f_1^-$$. This same artefact is seen in Fig. 3(b), albeit with smaller amplitude. Fig. 3(a) also contains a second, smaller artefact which results from the interaction between $$f_{1,m}^+$$ and an unrelated event in $$f_1^-$$. This crosstalk artefact is absent from Fig. 3(b) as expected from the discussion above. In this simple model, we do not show an individual reflector imaged using the last-arriving event from this focus point (the third method above), as it suffices to see that Figs 3(a) and (b) show the last reflector is correctly recovered. Therefore, the last-arriving event in $$f_1^-$$ has the traveltime of the desired primary and correctly images the deepest reflector. However, using RTM the last imaged reflector is spurious, and we see a total of four artifacts in the area imaged (two between the second and third reflectors, and two after the third reflector). These RTM artifacts generally appear at different locations than those in the focusing-related images. Hence, the combination of a focusing image with the RTM image using the method of da Costa Filho & Curtis (2016) is likely to produce an image devoid of most artifacts. Indeed, the result shown in Fig. 3(d), which exhibits the image obtained by combining images in panels (Figs 3a and c), shows a virtually artefact-free image. 4.2 2-D synthetic We now use a 2-D synthetic acoustic model with a constant velocity of 2.4 km s−1 and densities shown in Fig. 4. Varying only densities allows straight ray paths to be assumed, so that events observed in the data and focusing functions can be interpreted more easily. However, this is not a limitation of the new methods as Marchenko focusing functions can also be found for media with varying velocities, including highly complex media (Behura et al.2014; Vasconcelos et al.2015; Jia et al.2017). Figure 4. View largeDownload slide Acquisition and imaging geometry for 2-D synthetic model. Stars and triangles represent the 201 colocated surface sources and receivers. White dots represent the locations where focusing functions $$f_1^\pm$$ were computed. Figure 4. View largeDownload slide Acquisition and imaging geometry for 2-D synthetic model. Stars and triangles represent the 201 colocated surface sources and receivers. White dots represent the locations where focusing functions $$f_1^\pm$$ were computed. We choose several depth levels, each containing 151 focusing locations at which we compute focusing functions as shown by the white dots in Fig. 4. The bottom-most depth level at 1.6 km is used for imaging with pairs $$f_1^+$$ and $$f_1^-$$, and $$f_{1,d}^+$$ and $$f_1^-$$, displayed, respectively, in Figs 5(a) and (b); in this case, their Marchenko reference medium has the same reflectors as the full medium since the focusing level is beneath all reflectors. Figure 5. View largeDownload slide Images using source/receiver wavefields (a) $$f_1^+$$ and $$f_1^-$$, (b) $$f_{1,d}^+$$ and $$f_1^-$$ and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. (c) Conventional RTM image. (d) Combined image using images in panels (a) and (c). (e) Combined image using images in panels (b) and (c). Images in panels (a)–(e) use only virtual points in the bottom-most line of focus points, while (f) uses all virtual points shown in Fig. 4. Arrows indicate spurious reflectors. All panels are clipped at 10 per cent of their maximum amplitude. Figure 5. View largeDownload slide Images using source/receiver wavefields (a) $$f_1^+$$ and $$f_1^-$$, (b) $$f_{1,d}^+$$ and $$f_1^-$$ and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. (c) Conventional RTM image. (d) Combined image using images in panels (a) and (c). (e) Combined image using images in panels (b) and (c). Images in panels (a)–(e) use only virtual points in the bottom-most line of focus points, while (f) uses all virtual points shown in Fig. 4. Arrows indicate spurious reflectors. All panels are clipped at 10 per cent of their maximum amplitude. The focusing functions at other depth levels are used to image with last-arriving primaries only, that is, using the respective functions $$f_{1,d}^+$$ and $$f_{1,p}^-$$ for each depth level. The summation of these six images (each of which contains one reflector corresponding to a single depth level) is seen in Fig. 5(f). In order to use this method to recover all reflectors using only primaries, focus points are necessary below each interface. In this synthetic example, we choose a set of depths which we know are guaranteed to recover all reflections. In field data, this is not guaranteed. In addition, since this method requires picking, we do not believe that it is useful in areas with little knowledge of the imaging target. This method is much more appropriate for areas with a priori knowledge of the region wherein the imaging target is located, at which we wish to construct a more focused image. This is explored in the following field data example. As expected, Fig. 5(a) contains several artifacts resulting from crosstalk as indicated by the white arrows. These are reduced in Fig. 5(b), which features fewer spurious interfaces and artifacts with smaller amplitudes. Fig. 5(c) shows a standard RTM image, which has more artifacts in deeper sections compared to Figs 5(a) and (b). Since these artifacts appear in different locations, we combined the focusing-derived images with the standard RTM image to generate the combined images in Figs 5(d) and (e). These images are virtually free from artifacts, and also have fewer acquisition imprints (i.e. low-frequency artifacts common in RTM which can be seen mostly around the top reflector) than their constituent images. In addition, the combined image in Fig. 5(e) has slightly fewer artifacts than that in Fig. 5(d), consistent with the fact that fewer artifacts appear in Fig. 5(b) compared to Fig. 5(a). Fig. 5(f) also shows near perfect reconstruction of primaries. No coherent spurious interfaces are imaged, and all true reflectors have been recovered. However, this performance comes at the cost of picking the last event in $$f_1^-$$ (in this case performed automatically), as well as an increased computational cost due to the increased number of focus points. 4.3 Field data The methods were applied to an OBC field data set from the seabed above Volve oil field, located in the gas/condensate-rich Sleipner area 200 km off the North Sea coast of Norway (Fig. 6). Prior to applying the Marchenko method, the data were pre-processed with spatial interpolation, up/down wavefield separation, redatuming by multidimensional deconvolution, among others; the pre-processing is described by Ravasi et al. (2015, 2016), who detail the steps required to ensure the assumptions of the Marchenko method are met. Figure 6. View largeDownload slide Geographic location of the Volve field. Figure 6. View largeDownload slide Geographic location of the Volve field. Fig. 7(a) shows the velocity model used for the Marchenko method and for imaging, and Fig. 7(b) shows a magnification of the portion of the model in which we are interested. In both figures, the black dots show the 411 locations where focusing functions $$f_1^\pm$$ were computed using the pre-processed data. Figure 7. View largeDownload slide (a) Acquisition and imaging geometry. Triangles represent the 235 OBC receivers and stars represent the sources redatumed to the same receiver locations. The black box represents the imaged area, which is magnified in (b). Black dots represent locations where focusing functions were computed and the larger red dot shows the focusing location for functions in Fig. 8. The black arrow indicates a particularly strong contrast in the velocity model. Figure 7. View largeDownload slide (a) Acquisition and imaging geometry. Triangles represent the 235 OBC receivers and stars represent the sources redatumed to the same receiver locations. The black box represents the imaged area, which is magnified in (b). Black dots represent locations where focusing functions were computed and the larger red dot shows the focusing location for functions in Fig. 8. The black arrow indicates a particularly strong contrast in the velocity model. In Fig. 8, we show an example of a focusing function from the location marked by the large red dot in Fig. 7. In Fig. 8(c), we magnify a portion of the upgoing field shown in Fig. 8(b). Our picked final event in this gather is shown in Fig. 8(d). While it is clear from Fig. 8(c) that there is energy below the picked event, we verified that most of that energy is not coherent across all focusing locations (e.g. coherent energy must appear in every $$f_1^-$$ gather, and must at no point be crossed by the window used in the Marchenko iteration). Since the strong contrast in the velocity model (black arrow in Fig. 7b) appears across the entire model, we expected to see its associated primary along all focusing locations. The identified event shown in Fig. 8(d) was picked manually for a regularly spaced subset of the focusing locations, and the rest were picked automatically by applying interpolated muting windows to the focusing functions. The imaging examples below show that this procedure recovered the correct primary. Figure 8. View largeDownload slide (a) $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$ and (b) $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$, for $$\mathbf {x}_F = (6.7\,\mathrm{km}, 3.3\,\mathrm{km})$$ (large red dot in Fig. 7). The area enclosed by the black box in (b) is magnified in (c), and the result from windowing its last arrival is shown in (d). All panels are clipped at 50 per cent of the maximum amplitude of (b). Figure 8. View largeDownload slide (a) $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$ and (b) $$f_1^+(\mathbf {x}_0, \mathbf {x}_F, t)$$, for $$\mathbf {x}_F = (6.7\,\mathrm{km}, 3.3\,\mathrm{km})$$ (large red dot in Fig. 7). The area enclosed by the black box in (b) is magnified in (c), and the result from windowing its last arrival is shown in (d). All panels are clipped at 50 per cent of the maximum amplitude of (b). Fig. 9 shows the migrated images of the section of interest. Figs 9(a) and (b) are obtained, respectively, by using source/receiver wavefield pairs $$f_1^+$$ and $$f_1^-$$, and $$f_{1,d}^+$$ and $$f_1^-$$. These are virtually the same, with no notable differences other than in illumination and amplitudes. This happens because the coda of the downgoing focusing function $$f_{1,m}^+$$ has relatively low amplitude, and therefore contributes much less to the image than the direct part $$f_{1,d}^+$$, as evidenced in Fig. 8(a). The section of interest of the RTM image is shown in Fig. 9(c). The RTM image was provided from previous processing and fewer imaging points, spaced 8 m apart rather than 4 m as used in Figs 9(a) and (b). This may cause a decrease in resolution. In addition, the RTM image shows severely disrupted reflectors near the centre of the image when compared to Figs 9(a) and (b). These discontinuities may be caused by either multiples or backscattering which cross the reflector; such effects are not expected in our method where multiples do not play a role and where backscattering effects are diminished in the reference medium. However, the RTM image shows more continuity towards the edges of the image. The focusing images, which have low illumination near lateral edges, make identifying continuous reflectors harder at those locations. Figure 9. View largeDownload slide Images using source/receiver wavefields (a) $$f_1^+$$ and $$f_1^-$$, (b) $$f_{1,d}^+$$ and $$f_1^-$$ and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. (c) Conventional RTM image. (d) Combined image using images in panels (a) and (c). (e) Combined image using images in panels (b) and (c). Arrows indicate reflectors which appear better imaged in (a), (b) or (f) than in (c). Figure 9. View largeDownload slide Images using source/receiver wavefields (a) $$f_1^+$$ and $$f_1^-$$, (b) $$f_{1,d}^+$$ and $$f_1^-$$ and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. (c) Conventional RTM image. (d) Combined image using images in panels (a) and (c). (e) Combined image using images in panels (b) and (c). Arrows indicate reflectors which appear better imaged in (a), (b) or (f) than in (c). We show the combined images in Figs 9(d) and (e). In order to combine them, we interpolated the RTM image onto the finer grid used in the focusing images which may cause a decrease in resolution in the combined images compared to panels (a) and (b). The combined images also show generally fewer acquisition imprints. However, they are severely harmed by the discontinuities which appear in RTM, leading to a worse interpretability than their constituent non-RTM images. Nevertheless, continuous reflectors are much clearer in the combined image than in either of the constituent images. Fig. 9(f) shows the migration of the primary picked from $$f_1^-$$ at each point along the focus point line. An example of these gathers was shown in Fig. 8. We expected the primary to be related to the strong velocity contrast in the model indicated in Fig. 7. Indeed, Fig. 9(f) shows that we have individually imaged that strong contrast in the model, which also appears prominently in Figs 9(a) and (b). 5 DISCUSSION The results shown above demonstrate the advantages and limitations of each imaging method. The first method uses $$f_1^+$$ and $$f_1^-$$ as source/receiver wavefield pairs to create the image and provides images that are inherently different from those obtained using standard RTM, provided that accurate focusing is achieved in the Marchenko method. Using 1-D synthetic data where the focusing was near-perfect, this first method and standard RTM show, as expected, artifacts in different locations. Also in the 2-D synthetic model, artifacts mostly appear in different locations, suggesting that the focusing achieved was also acceptable. For field data, though hard to see from the individual images themselves, the combined images (discussed below) also show that artifacts from RTM appear in a different locations than in our first method. Nevertheless, our first method exhibits crosstalk between $$f_{1}^+$$ and unrelated events in $$f_{1}^-$$. In order to reduce these, we tested substituting $$f_1^+$$ for its direct wave component $$f_{1,d}^+$$. We note that this is not the only possible approach to reduce artifacts, with deconvolution of receiver wavefields by source wavefields being another option. Deconvolution, however, is often unstable. Our approach, on the other hand, is stable and reliably produces images with fewer artifacts, as shown in the synthetic models. However, for the field data it showed no significant improvement over our previous method. This is because in this particular case the coda components are weak as evidenced in Fig. 8(a), so the direct wave also dominates the downgoing field when we use the full $$f_1^+$$ as the source field. Nevertheless, when compared to RTM, this method exhibits improvements in reflector continuity, likely because of the effects of backscattering which are not accounted for in RTM. Similarities between the images obtained by our methods and RTM were exploited by using the combined imaging method of da Costa Filho & Curtis (2016). Our synthetic results show near perfect recovery of all reflectors with minimal interference from multiple or acquisition-related artifacts. Using field data, the combined images also reduced artifacts, but the low quality of the RTM image and the lack of continuity near the edges of the focusing images degraded the resulting combined image. Our final imaging method uses the windowed $$f_1^-$$ functions, denoted $$f_{1,p}^-$$. These retain only their last events which are primaries related to the reflector immediately above the focus point. This method allows the imaging of individual reflectors without interference (crosstalk) from other parts of the wavefield. However, if an image of the entire medium is desired, each reflector must be imaged separately which in turn requires more computational power to calculate focusing functions at an increased number of focus points. Nevertheless, it has been shown to provide clean images of individual reflectors, even in complex field data. While picking for the 2-D synthetic data was done automatically, picking consistent arrivals in the field data was performed manually to ensure that other events in $$f_1^-$$ where not picked in error by the automatic picking routine. This increases processing time significantly, and hence this method may be most applicable for target-oriented applications for which the extra time may be worthwhile if it results in a better image of a reservoir. The methods shown for acoustic media have a straightforward extension to elastic media, with certain caveats which we highlight below using the synthetic 1.5-D model in Fig. 10(a). da Costa Filho et al. (2014), Wapenaar (2014) and Wapenaar & Slob (2014) have shown how focusing functions may be computed in elastic media. Our first and second methods, which do not require windowing, may then be used to compute elastic images by migrating source/receiver pairs composed of down- and upgoing elastic focusing functions. The resulting SS wave images are shown in Fig. 10(d) using fields $$f_1^+$$ and $$f_1^-$$, and Fig. 10(b) using $$f_{1,d}^+$$ and $$f_1^-$$. In this simple model, there is no crosstalk between $$f_{1,m}^+$$ and events in $$f_1^-$$ in Fig. 10(d) or between $$f_{1,d}^+$$ and events in $$f_1^-$$ in Fig. 10(b), but there are S-to-P conversion-related artifacts. They compare well to RTM (Fig. 10c), and combining our image in Fig. 10(b) with RTM provides the virtually artefact-free image shown in Fig. 10(e). Figure 10. View largeDownload slide (a) 1.5-D synthetic density model used to generate elastic data. Colocated sources and receivers are shown as stars and triangles, and focus locations are shown as dots. P and S velocities are constant at 2.5 and 1.3 km s−1, respectively. SS images are constructed using source/receiver wavefields (b) $$f_{1,d}^+$$ and $$f_1^-$$, (d) $$f_{1}^+$$ and $$f_1^-$$, and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. RTM is shown in (c), and its combined image with (b) is shown in (e). White arrows indicate artifacts. Images are clipped at 5 per cent of their maximum amplitude. Figure 10. View largeDownload slide (a) 1.5-D synthetic density model used to generate elastic data. Colocated sources and receivers are shown as stars and triangles, and focus locations are shown as dots. P and S velocities are constant at 2.5 and 1.3 km s−1, respectively. SS images are constructed using source/receiver wavefields (b) $$f_{1,d}^+$$ and $$f_1^-$$, (d) $$f_{1}^+$$ and $$f_1^-$$, and (f) $$f_{1,d}^+$$ and $$f_{1,p}^-$$. RTM is shown in (c), and its combined image with (b) is shown in (e). White arrows indicate artifacts. Images are clipped at 5 per cent of their maximum amplitude. The third method which requires windowing $$f_1^-$$ fields cannot be extended as easily to elastic media. In particular, the assumption that the last arriving wave is a pure-mode primary used for acoustic media does not hold for all $$f_1^-$$ gathers. A focus originating from a compressional time-reversed direct wave will create not only P waves, but also conversions which will arrive at or after the last pure-mode P primary in the upgoing focusing functions: these have certainly reflected only once, but they also forward-scatter through conversion, making them non-compliant with the first-order Born assumption and harder to migrate. However, the assumption is valid for S-wave foci: the last event arriving in the upgoing field created by an S direct wave will also be a pure-mode S wave. Any conversion into P will have traveled faster than the pure S wave. As such, the third method of primaries can be extended directly to pure-mode S waves in elastic media. This is verified by Fig. 10(d), which uses SS data to show the recovery of only the reflector directly above the focus points, without any conversion-related artifacts. These three methods and their application to several models provides significant insights into several areas of active research in seismics. In standard RTM, one assumes that the point-source function only creates primaries. This is known as the linearization, or first-order Born approximation, whereby the receiver wavefield created by the source function is linearly related to the medium parameters. In order for this approach to be valid, data used in standard RTM must first have all multiples removed. In our first method, the source function $$f_1^+$$ creates all events in $$f_1^-$$. Indeed, since all events in $$f_1^-$$ are primaries, this imaging method is truly linear in the sense that events in $$f_1^-$$ are linearly related to the source function ($$f_1^+$$) by the reference medium parameters. This is achieved without any multiple removal. The second method, while not respecting wave propagation fully ($$f_{1,d}^+$$ only generates some of the primary energy in $$f_{1}^-$$), has been shown to be a useful approximation in terms of imaging as the approximation made in the physics trades-off particularly with fewer artifacts caused by the migration method itself. The third method also creates clean images in practice, and can be seen as a truly linear method which uses $$f_{1,d}^+$$ as source and only one primary in $$f_{1}^-$$ as the received wavefield. Finally, we have shown new examples of how the combined imaging conditions reduce artifacts in pairs of images, including a successful application to field data. So-called Marchenko imaging requires Green’s and focusing functions to be computed at all imaging locations. This requires the initial focusing field to be computed from each source/receiver location to the subsurface. In addition, it requires that we solve the Marchenko equation separately at all subsurface imaging locations. The cost of obtaining the image is then applying an imaging condition to these fields, the cost of which is much lower than the aforementioned steps. If the initial focusing field is computed with finite differences and stored for all sources and subsurface locations, the cost of computing these fields is the same as computing the source RTM wavefields. The solution to the Marchenko equation relies on computing cross-correlations and convolutions with the reflection data, a process which supplants the receiver-side backpropagation in RTM. This step is usually significantly more costly than RTM, and one cross-correlation/convolution is taken to be equivalent to receiver-side backpropagation (for all source locations). Therefore, Marchenko imaging methods are generally considered to be several times more costly than conventional RTM (Slob et al.2014). Focusing functions, on the other hand, allow structures to be targeted with a small number of focus points. This means that the initial focusing fields and Marcheko solutions need to be computed at fewer locations. For example, in the 2-D synthetic example we computed these fields at 151 focusing locations, while a Marchenko image would require 301151 = 751 × 401 of these computations corresponding to each image location in the model. After focusing functions have been computed, they need to be used in the forward and backwards extrapolation steps like RTM, which dominates the cost. In the 2-D synthetic, we performed 151 of these computations for our focusing images versus 201 for standard RTM; in the field data, these numbers were 411 and 235, respectively. Therefore, whether the focusing images or RTM will be more computationally intensive depends on the application. At the cost of performing both the standard RTM and our focusing RTM, the combined images can deliver very clean images. As such, these methods can be used in a target oriented way to monitor the evolution of producing reservoirs over time. Moreover, the first two methods can be used in place of demultiple plus RTM or Marchenko redatuming plus RTM, especially for geologies which generate many multiples potentially driving down the total processing cost. 6 CONCLUSION We have presented novel methods to image subsurface structures using specially crafted virtual acquisitions. These acquisitions are obtained by the outputs of the Marchenko method, and are composed of source/receiver wavefields given by the down- and upgoing focusing functions, respectively. Additionally, we present a method based on these acquisitions which can image individual reflectors: this method involves windowing upgoing focusing functions based on generally applicable properties of the traveltimes of events in focusing functions, making it feasible to perform on field data. We apply these methods to field data from the North Sea, to acoustic 1-D and 2-D synthetic data, and to elastic 1.5-D synthetic data. We discuss the advantages and limitations of the methods, as well as implications they might have for other areas of seismic processing. Acknowledgements The authors thank the Edinburgh Interferometry Project sponsors (Schlumberger Gould Research, Statoil and Total) for supporting this research. We thank Statoil ASA, the Statoil Volve team and the Volve license partners ExxonMobil E&P Norway and Bayerngas Norge for release of the data. The first author also thanks CAPES for funding. The authors also thank the Editor, Lapo Boschi, as well as reviewer Filippo Broggini and an anonymous reviewer for their contributions. REFERENCES Amundsen L., 2001. Elimination of free-surface related multiples without need of the source wavelet, Geophysics , 66( 1), 327– 341. https://doi.org/10.1190/1.1444912 Google Scholar CrossRef Search ADS   Bakulin A., Calvert R., 2006. The virtual source method: theory and case study, Geophysics , 71( 4), SI139– SI150. https://doi.org/10.1190/1.2216190 Google Scholar CrossRef Search ADS   Baysal E., Kosloff D.D., Sherwood J. W.C., 1983. Reverse time migration, Geophysics , 48( 11), 1514– 1524. https://doi.org/10.1190/1.1441434 Google Scholar CrossRef Search ADS   Behura J., Wapenaar K., Snieder R., 2014. Autofocus imaging: image reconstruction based on inverse scattering theory, Geophysics , 79( 3), A19– A26. https://doi.org/10.1190/geo2013-0398.1 Google Scholar CrossRef Search ADS   Berkhout A.J., 2012. Combining full wavefield migration and full waveform inversion, a glance into the future of seismic imaging, Geophysics , 77( 2), S43– S50. https://doi.org/10.1190/geo2011-0148.1 Google Scholar CrossRef Search ADS   Berkhout A.J., Verschuur D.J., 2005. Removal of internal multiples with the common-focus-point (CFP) approach: Part 1—explanation of the theory, Geophysics , 70( 3), V45– V60. https://doi.org/10.1190/1.1925753 Google Scholar CrossRef Search ADS   Berkhout A.J., Verschuur D.J., 2006. Imaging of multiple reflections, Geophysics , 71( 4), SI209– SI220. https://doi.org/10.1190/1.2215359 Google Scholar CrossRef Search ADS   Berryhill J.R., 1979. Wave-equation datuming, Geophysics , 44( 8), 1329– 1344. https://doi.org/10.1190/1.1441010 Google Scholar CrossRef Search ADS   Berryhill J.R., 1984. Wave-equation datuming before stack, Geophysics , 49( 11), 2064– 2066. https://doi.org/10.1190/1.1441620 Google Scholar CrossRef Search ADS   Broggini F., Wapenaar K., van der Neut J., Snieder R., 2014. Data-driven Green’s function retrieval and application to imaging with multidimensional deconvolution, J. geophys. Res. , 119( 1), 425– 441. https://doi.org/10.1002/2013JB010544 Google Scholar CrossRef Search ADS   Chang W.-F., 1986. Reverse-time migration of offset vertical seismic profiling data using the excitation-time imaging condition, Geophysics , 51( 1), 67, doi:10.1190/1.1442041. https://doi.org/10.1190/1.1442041 Google Scholar CrossRef Search ADS   Curtis A., Halliday D., 2010. Source-receiver wave field interferometry, Phys. Rev. E , 81( 4), 46601. https://doi.org/10.1103/PhysRevE.81.046601 Google Scholar CrossRef Search ADS   Curtis A., Gerstoft P., Sato H., Snieder R., Wapenaar K., 2006. Seismic interferometry—turning noise into signal, Leading Edge , 25( 9), 1082– 1092. https://doi.org/10.1190/1.2349814 Google Scholar CrossRef Search ADS   da Costa Filho C.A., Curtis A., 2016. Attenuating multiple-related imaging artifacts using combined imaging conditions, Geophysics , 81( 6), S469– S475. https://doi.org/10.1190/geo2016-0113.1 Google Scholar CrossRef Search ADS   da Costa Filho C.A., Ravasi M., Curtis A., Meles G.A., 2014. Elastodynamic Green’s function retrieval through single-sided Marchenko inverse scattering, Phys. Rev. E , 90( 6), 063201. https://doi.org/10.1103/PhysRevE.90.063201 Google Scholar CrossRef Search ADS   da Costa Filho C.A., Ravasi M., Curtis A., 2015. Elastic P- and S-wave autofocus imaging with primaries and internal multiples, Geophysics , 80( 5), S187– S202. https://doi.org/10.1190/geo2014-0512.1 Google Scholar CrossRef Search ADS   da Costa Filho C.A., Meles G.A., Curtis A., 2017. Elastic internal multiple analysis and attenuation using Marchenko and interferometric methods, Geophysics , 82( 2), Q1– Q12. https://doi.org/10.1190/geo2016-0162.1 Google Scholar CrossRef Search ADS   Fink M., 1992. Time reversal of ultrasonic fields. I. Basic principles., IEEE Trans. Ultrason. Ferroelectr. Freq. Control , 39( 5), 555– 66. https://doi.org/10.1109/58.156174 Google Scholar CrossRef Search ADS PubMed  Halliday D., Curtis A., 2010. An interferometric theory of source-receiver scattering and imaging, Geophysics , 75( 6), SA95– SA103. https://doi.org/10.1190/1.3486453 Google Scholar CrossRef Search ADS   Jia X., Guitton A., Singh S., Snieder R., 2017. Subsalt Marchenko imaging: a Gulf of Mexico example, in SEG Technical Program Expanded Abstracts 2017 , pp. 5588– 5592, Society of Exploration Geophysicists. Jones I.F., 2010. An Introduction To: Velocity Model Building , EAGE Publications bv. Google Scholar CrossRef Search ADS   Kinneging N.A., Budejicky V., Wapenaar C. P.A., Berkhout A.J., 1989. Efficient 2D and 3D shot record redatuming, Geophys. Prospect. , 37( 5), 493– 530. https://doi.org/10.1111/j.1365-2478.1989.tb02220.x Google Scholar CrossRef Search ADS   Malcolm A.E., Ursin B., de Hoop M.V., 2009. Seismic imaging and illumination with internal multiples, Geophys. J. Int. , 176( 3), 847– 864. https://doi.org/10.1111/j.1365-246X.2008.03992.x Google Scholar CrossRef Search ADS   Masson Y., Romanowicz B., 2017. Box tomography: localized imaging of remote targets buried in an unknown medium, a step forward for understanding key structures in the deep Earth, Geophys. J. Int. , 211( 1), 141– 163. https://doi.org/10.1093/gji/ggx141 Google Scholar CrossRef Search ADS   McMechan G.A., 1983. Migration by extrapolation of time–dependent boundary values, Geophys. Prospect. , 31( 3), 413– 420. https://doi.org/10.1111/j.1365-2478.1983.tb01060.x Google Scholar CrossRef Search ADS   Meles G.A., Löer K., Ravasi M., Curtis A., da Costa Filho C.A., 2015. Internal multiple prediction and removal using Marchenko autofocusing and seismic interferometry, Geophysics , 80( 1), A7– A11. https://doi.org/10.1190/geo2014-0408.1 Google Scholar CrossRef Search ADS   Meles G.A., Wapenaar K., Curtis A., 2016. Reconstructing the primary reflections in seismic data by Marchenko redatuming and convolutional interferometry, Geophysics , 81( 2), Q15– Q26. https://doi.org/10.1190/geo2015-0377.1 Google Scholar CrossRef Search ADS   Meles G.A., da Costa Filho C.A., Curtis A., 2017. Estimating and imaging with primaries constructed from single-sided focusing wavefields, in 79th EAGE Conference and Exhibition , doi:10.3997/2214-4609.201701056. Muijs R., Robertsson J. O.A., Holliger K., 2007. Prestack depth migration of primary and surface-related multiple reflections: Part I - Imaging, Geophysics , 72, S59, doi:10.1190/1.2422796. https://doi.org/10.1190/1.2422796 Google Scholar CrossRef Search ADS   Oristaglio M.L., 1989. An inverse scattering formula that uses all the data, Inverse Probl. , 5( 6), 1097– 1105. https://doi.org/10.1088/0266-5611/5/6/015 Google Scholar CrossRef Search ADS   Ravasi M., Vasconcelos I., Kritski A., Curtis A., da Costa Filho C.A., Meles G.A., 2015. Marchenko imaging of Volve field, North Sea, in 77th EAGE Meeting and Technical Exhibition 2015 , doi:10.3997/2214-4609.201412938. Ravasi M., Vasconcelos I., Kritski A., Curtis A., da Costa Filho C.A., Meles G.A., 2016. Target-oriented Marchenko imaging of a North Sea field, Geophys. J. Int. , 205( 1), 99– 104. https://doi.org/10.1093/gji/ggv528 Google Scholar CrossRef Search ADS   Rose J.H., 2001. “Single-sided” focusing of the time-dependent Schrödinger equation, Phys. Rev. A , 65( 1), 0127071. https://doi.org/10.1103/PhysRevA.65.012707 Google Scholar CrossRef Search ADS   Schuster G.T., Zhou M., 2006. A theoretical overview of model-based and correlation-based redatuming methods, Geophysics , 71( 4), SI103– SI110. https://doi.org/10.1190/1.2208967 Google Scholar CrossRef Search ADS   Schuster G.T., Yu J., Sheng J., Rickett J., 2004. Interferometric/daylight seismic imaging, Geophys. J. Int. , 157( 2), 838– 852. https://doi.org/10.1111/j.1365-246X.2004.02251.x Google Scholar CrossRef Search ADS   Singh S., Snieder R., Behura J., van der Neut J., 2015. Marchenko imaging: imaging with primaries, internal multiples, and free-surface multiples, Geophysics , 80( 5), S164– S174. https://doi.org/10.1190/geo2014-0494.1 Google Scholar CrossRef Search ADS   Slob E., Wapenaar K., Broggini F., Snieder R., 2014. Seismic reflector imaging using internal multiples with Marchenko-type equations, Geophysics , 79( 2), S63– S76. https://doi.org/10.1190/geo2013-0095.1 Google Scholar CrossRef Search ADS   Staring M., Pereira R., Douma H., van der Neut J., Wapenaar C., 2017. Adaptive double-focusing method for source-receiver Marchenko redatuming on field data, in SEG Technical Program Expanded Abstracts 2017 , pp. 4808– 4812, Society of Exploration Geophysicists. Tang Y., Biondi B., 2011. Target-oriented wavefield tomography using synthesized Born data, Geophysics , 76( 5), WB191– WB207. https://doi.org/10.1190/geo2010-0383.1 Google Scholar CrossRef Search ADS   Thorbecke J., 1997. Common focus point technology, PhD thesis , Delft University of Technology. van der Neut J., Wapenaar K., 2016. Adaptive overburden elimination with the multidimensional Marchenko equation, Geophysics , 81( 5), T265– T284. https://doi.org/10.1190/geo2016-0024.1 Google Scholar CrossRef Search ADS   Vasconcelos I., Wapenaar K., van der Neut J., Thomson C., Ravasi M., 2015. Using inverse transmission matrices for Marchenko redatuming in highly complex media, in SEG Technical Program Expanded Abstracts 2015 , pp. 5081– 5086, Society of Exploration Geophysicists. Verschuur D.J., 1992. Adaptive surface-related multiple elimination, Geophysics , 57( 9), 1166– 1177. https://doi.org/10.1190/1.1443330 Google Scholar CrossRef Search ADS   Wapenaar C. P.A., Kinneging N., Berkhout A., 1987. Principle of prestack migration based on the full elastic two-way wave equation, Geophysics , 52( 2), 151– 173. https://doi.org/10.1190/1.1442291 Google Scholar CrossRef Search ADS   Wapenaar K., 2014. Single-sided Marchenko focusing of compressional and shear waves, Phys. Rev. E , 90( 6), 063202. https://doi.org/10.1103/PhysRevE.90.063202 Google Scholar CrossRef Search ADS   Wapenaar K., Slob E., 2014. On the Marchenko equation for multicomponent single-sided reflection data, Geophys. J. Int. , 199( 3), 1367– 1371. https://doi.org/10.1093/gji/ggu313 Google Scholar CrossRef Search ADS   Wapenaar K., Broggini F., Slob E., Snieder R., 2013. Three-dimensional single-sided Marchenko inverse scattering, data-driven focusing, Green’s function retrieval, and their mutual relations, Phys. Rev. Lett. , 110( 8), 084301. https://doi.org/10.1103/PhysRevLett.110.084301 Google Scholar CrossRef Search ADS PubMed  Wapenaar K., Thorbecke J., van der Neut J., Broggini F., Slob E., Snieder R., 2014. Marchenko imaging, Geophysics , 79( 3), WA39– WA57. https://doi.org/10.1190/geo2013-0302.1 Google Scholar CrossRef Search ADS   Wapenaar K., van der Neut J., Slob E., 2016. Unified double- and single-sided homogeneous Green’s function representations, Proc. R. Soc. A: Math. Phys. Eng. Sci. , 472( 2190), 20160162. https://doi.org/10.1098/rspa.2016.0162 Google Scholar CrossRef Search ADS   Weglein A.B., Gasparotto F.A., Carvalho P.M., Stolt R.H., 1997. An inverse-scattering series method for attenuating multiples in seismic reflection data, Geophysics , 62( 6), 1975– 1989. https://doi.org/10.1190/1.1444298 Google Scholar CrossRef Search ADS   Yuan S., Fuji N., Singh S., Borisov D., 2017. Localized time-lapse elastic waveform inversion using wavefield injection and extrapolation: 2-D parametric studies, Geophys. J. Int. , 209( 3), 1699– 1717. https://doi.org/10.1093/gji/ggx118 Google Scholar CrossRef Search ADS   Ziolkowski A., Taylor D.B., Johnston R. G.K., 1999. Marine seismic wavefield measurement to remove sea-surface multiples, Geophys. Prospect. , 47( 6), 841– 870. https://doi.org/10.1046/j.1365-2478.1999.00165.x Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Unlimited reading Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere. ### Stay up to date Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates. ### Organize your research It’s easy to organize your research with our built-in tools. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. ### Monthly Plan • Read unlimited articles • Personalized recommendations • No expiration • Print 20 pages per month • 20% off on PDF purchases • Organize your research • Get updates on your journals and topic searches$49/month

14-day Free Trial

Best Deal — 39% off

### Annual Plan

• All the features of the Professional Plan, but for 39% off!
• Billed annually
• No expiration
• For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588$360/year

billed annually

14-day Free Trial