An introduction to Marchenko methods for imaging

INTRODUCTION

The aim of seismic imaging is to map unknown heterogeneities in the earth’s subsurface, given a wavefield measured on or close to the earth’s surface. An approximate seismic wave speed model, usually called a velocity model, is required to map the subsurface accurately. This model provides a basic level of understanding about how seismic waves propagate through the subsurface and allows seismic information measured at the surface to be mapped to approximately correct subsurface locations. Much effort goes into estimating the velocity model using migration velocity analysis (Yilmaz, 2001; Sava and Biondi, 2004), traveltime tomography (Stork, 1992; Jones, 2010), and full-waveform inversion (Tarantola, 1984; Pratt et al., 1998; Virieux and Operto, 2009), but it is always imperfect. In particular, it is usually far smoother than the true earth and therefore is not kinematically accurate; in other words, it does not map waves that reflect from abrupt interfaces to their true subsurface positions. Even when these errors are sufficiently small that the image produced is correctly positioned, the inaccuracies usually cause other additional artifacts to be superimposed on the final image. Artifacts that are often most troublesome are those created by recorded seismic waves that reflect more than once in the subsurface, called multiples. This tutorial explains a set of methods that account for such waves so that these artifacts do not occur.

Marchenko methods (Rose, 2001; Broggini et al., 2012) are data-driven methods that use measured surface seismic data and an approximate velocity model to calculate the signal that would have been recorded at the surface if an impulsive, frequency band-limited source had fired at each chosen subsurface image point — including multiples. The estimated signals are called (frequency band-limited) Green’s functions, and they are exactly the information needed for accurate subsurface imaging (Behura et al., 2014; Wapenaar et al., 2014), seismic redatuming (Wapenaar et al., 2014), or identifying and removing multiples (Meles et al., 2015, 2016).

The name Marchenko comes from the author of the original work on inverse scattering (Marchenko, 1955), who devised methods to estimate Green’s functions in the field of quantum mechanics in one dimension (Snieder [2015] provides more information about this application). More recently, a solution to the so-called Marchenko equations was formulated for geophysical applications that allow 2D and 3D media to be imaged under certain approximations (Wapenaar et al., 2013; Lomas and Curtis, 2017).

This paper presents an intuitive introduction to the Marchenko method and its applications. The aim is not to introduce fundamentally new concepts, but to provide an easily accessible guide to some of the key concepts and methods that already exist. Additionally, a Marchenko MATLAB code and a relevant data set accompany this paper: The code is well-commented, easily editable, and adaptable for 2D seismic problems. It is constructed so as to give readers further insight into the workflow used to calculate Green’s functions using Marchenko methods and to allow them to experiment and gain comfort with the methods — rather than being geared toward computationally efficient construction of large seismic images. Nevertheless, all 2D examples presented herein were constructed using this code; hence, it is perfectly sufficient to be used to teach and learn about Marchenko methods and to process small data sets. Other codes exist in the public domain (e.g., Thorbecke et al., 2017), but our code is designed specifically for user experimentation, and so it is written in a more intuitively accessible (higher level) programming language. It is therefore also ideal as an aid for teaching about Marchenko methods in masters or professional development courses.

Marchenko methods are simplest, most intuitive, and most accurate for 1D problems, so the first section of this paper introduces Green’s functions estimation for a simple 1D medium in which full wavefields can be displayed and understood. We then introduce the reader to 2D examples, the accompanying MATLAB code, and the application of Marchenko imaging. Throughout this paper, multiple data sets are used: All are constructed in acoustic media and exclude free surface multiples because such data allow the simplest and most studied form of Marchenko methods to be applied. However, theory exists for Marchenko methods using elastic data (da Costa Filho et al., 2014, 2015; Wapenaar, 2014; Wapenaar and Slob, 2014) and data containing surface-related multiples (Singh et al., 2015, 2017). There are also limited examples of applications to real data (Ravasi et al., 2016; Jia et al., 2018; Wapenaar et al., 2018). For an entirely nonmathematical introduction to Marchenko methods, we refer readers to van der Neut et al. (2015c), or for an introduction to 1D Marchenko methods, see Cui et al. (2018); for a thorough introduction to the more sophisticated mathematical aspects of the Marchenko methods, see Slob et al. (2014b), Wapenaar et al. (2014), or van der Neut et al. (2015b). Our tutorial fills the niche between these studies by introducing the concepts, mathematics, and computational machinery in an accessible way, with a code designed to facilitate experimentation and education.

THE MARCHENKO METHOD

There are multiple applications of Marchenko methods (imaging, redatuming, constructing primaries, and multiple removal), but they all have the same foundation, namely, Green’s function estimation. Green’s functions are the waves that arrive at a receiver location due to the firing of a spatiotemporally impulsive source. We represent these Green’s functions as $G(x_0,x_i,t)$, where $x_0$ is the location of a receiver on the recording surface, $x_i$ is a source point in the subsurface, and $t$ represents the time domain. In this syntax, each term is a signal with two locations ($x$): The second always denotes the source location, and the first is the receiver location. The Marchenko method estimates Green’s functions between an arbitrarily chosen image point (or an artificial or virtual source) within the subsurface, and any point within the surface acquisition array (Figure 1).

The most basic form of Green’s function estimation is to assume or estimate an initial approximate velocity model and estimate Green’s functions $G(x_0,x_i,t)$ using either ray propagation or wavefield calculation through that model. This is standard practice in reverse time migration (RTM) (e.g., Baysal et al., 1983). Marchenko methods provide a workflow to estimate Green’s functions but decomposed into two constituent parts: The first part consists of all waves that are upgoing ($-$) at the image point in the earth’s subsurface, whereas the second part consists of the downgoing ($+$) waves. This includes components of the wavefield that have undergone multiple reflections, so-called multiples. In other words, two Green’s functions can be constructed from each subsurface image point, which are recorded at the surface: The first $G^{-}$ contains signals that start at the image point as a source wavefield propagating upward and the second $G^{+}$ contains signals that initially propagate downward from the image point and are reflected back up to the surface.

For simplicity, let us begin with the 1D Marchenko method. Two pieces of information are needed to calculate the decomposed Green’s functions $G^{+}$ and $G^{-}$. The first is the reflectivity from a point source at the surface measured by a point receiver at the surface, denoted by $R(x_0,x_0,t)$; in the real world, this is an idealized version of a 1D surface seismic reflection data after surface-related multiple removal. The second is an estimate of the direct (nonreflected) wave arrival $T_d(x_i,x_0,t)$ between the surface source and an image point. The decomposed Green’s functions $G^{+/-}$ between $x_0$ and $x_i$ are related to the reflectivity $R$ through additional terms $f^{+}$ and $f^{-}$, which are called focusing functions and are the subject of the next section:

Equations 1 and 2 are defined in the time domain. Symbol $-t$ (and also later in this paper, superscript *) denotes the time reversal of the signal that precedes it. This is accomplished if we flip the positive and negative time axis of the initial signal, an operation that corresponds to complex conjugation in the frequency domain. The symbol $\otimes$ represents a time-domain convolution, which is equivalent to multiplication in the frequency domain.

It is worth noting that equations 1 and 2 differ from those given in most of the existing literature on Marchenko methods. To aid intuition, we have created a virtual source at the image point inside the subsurface rather than a virtual receiver (the latter is more common). Comparing these cases, the direction of wave propagation is reversed: $G^{+}(x_0,x_i,t)=G^{-}(x_i,x_0,t)$. However, the property of source-receiver reciprocity states that these are equivalent (identical signals). We will continue to use the virtual source syntax for the remainder of this paper.

Focusing functions

Focusing functions are keys for understanding Marchenko methods. Imagine throwing a stone into a still pond on a windless day: Ripples diverge from the location of impact, propagating as waves across the water surface. Let us imagine that these ripples are recorded on some closed boundary of receivers that surrounds the impact point. If we waited until all of the energy had settled, we could then use the receivers as sources to inject the recorded wavefield back into the pond. If we do this in time-reversed order (inject the last wave recorded at each receiver first), the original ripples will be recreated, but this time, they would converge inward rather than propagating outward (Cassereau and Fink, 1992). They will all eventually refocus at the impact point, then diverge outward again, creating another wavefield that can be recorded at the receiver boundary. In this thought experiment, the (time-reversed) wavefield injected on the boundary is called a focusing function: It defines exactly which waves we should inject to focus the ingoing energy at the impact point.

Focusing functions used in Marchenko methods are intuitively similar to those above. The only conceptual differences are that the source point in this case is the subsurface source in Figure 1, and that the receiver boundary is at the earth’s surface and so is only on one side of the source point. Downgoing focusing functions are related to wavefields that, if injected at the earth’s surface, would focus (collapse all of their energy to a point) at a specific location in the subsurface (here, the location of any chosen virtual source or image point). However, in the case of focusing in the subsurface, this only occurs in an idealized (truncated) model of the earth’s subsurface structure, which is homogeneous below the depth of that point, but which has the true earth’s structure above that depth.

Focusing means that there is a time at which the waves at a certain depth only exist at one specific image point — everywhere else at that depth the wavefield is zero. The function $f^{+}$ is the wavefield that we would have to inject at the surface (at point $x_0$) in order for the wavefield to focus at the image point; hence, this wavefield is downgoing at the surface. Function $f^{-}$ is the wavefield that we would record at the surface as we inject $f^{+}$ in the truncated model; hence, $f^{-}$ is upgoing at the surface. Both wavefields are shifted along the time axis such that the focus occurs at time zero.

Figure 2 includes a standard representation of a focusing function (Slob et al., 2014b) for a simple 1D subsurface model that consists of layers with varying density and a constant velocity. First, Figure 2a shows the wavefield that develops in space and time when a simple impulsive source (convolved with a Ricker wavelet) is injected at the surface at time $t=0$ (note that time is on the horizontal axis). This consists of a direct wave (the first continuous, linear wave on the left) and a set of (singly and multiply) reflecting waves. At a particular image point in depth (e.g., $1400\mathrm{m}$ — indicated by an arrow in Figure 2), multiple waves arrive; hence, there is not a focus of energy. To create such a focus, additional energy must be injected to cancel all but the direct wave at that point.

The focusing function is the signal at the surface ($\mathrm{depth}=0\mathrm{m}$) in Figure 2b, shown in Figure 2c. The downgoing component $f^{+}$ is the signal injected at the surface to create the focus at $\mathrm{depth}=1400\mathrm{m}$, shown by the circle in Figure 2b. The upgoing component $f^{-}$ is the reflected response observed at the surface from this injected signal. It can be seen that three pulses of energy ($\alpha$, $\beta$, and $\gamma$) are injected at $x_0$ in addition to the initial pulse $\delta$ to cancel out the reflected components of the wavefield observed in Figure 2a. These three signals together with $\delta$ make up the complete downgoing focusing function $f^{+}$ and all of the up-coming waves at $\mathrm{depth}=0\mathrm{m}$ comprise $f^{-}$.

Iterative solution

The Marchenko method works by first calculating the focusing functions and then using equations 1 and 2 to estimate the Green’s functions. Although the relationships between focusing functions and Green’s functions in those equations are relatively simple, they do not explain how one can calculate focusing functions. Several methods have been proposed to do this (Broggini et al., 2014; van der Neut et al., 2015a, 2015b); here, we present the method of Wapenaar et al. (2014) because it can be understood most intuitively.

In a 1D system, we assume that we know the reflectivity at the surface $R(x_0,x_0,t)$ as well as the direct arrival between the surface and the image point $T_d(x_i,x_0,t)$ — which identifies the chosen image point $x_i$. The first step in estimating the focusing functions is to set

Equation 3 inverts the direct arrival, commonly approximated as performing a time reversal (switching the time axis of $T_d$ and setting the signal at positive times to zero: $T_d(x_i,x_0,t{)}^{-1}\approx T_d(x_i,x_0,-t)$). The result is used as a first approximation for $f^{+}$ denoted $f_0^{+}$ (Wapenaar et al., 2014) and forms the component $\delta$ from Figure 2b. This makes intuitive sense: If we time reverse the direct wave between the image point and the surface, it will propagate back to its source point (the image point) and create a pulse of energy there at zero time, just as in the example of ripples on the pond. Unfortunately, though, as it propagates back into the subsurface, some of its energy will scatter or reflect from heterogeneities in the earth, creating a more complex part of the wavefield that will disrupt the focus. Marchenko methods design energy to inject to destructively interfere with these scattered waves, reducing them to zero amplitude.

The estimate for $f_0^{+}$ can then be used to estimate $f_0^{-}$:

Within the square brackets, equation 4 convolves the initial estimate of $f_0^{+}$ with the reflectivity, which is equivalent to injecting $f_0^{+}$ into the earth and recording the result at the surface. Again, this is equivalent to injecting the time-reversed wavefield in the pond and recording the reflecting waves on the source boundary. The additional term in this equation $\theta$ is a focusing-location-dependent window, which removes all energy that arrives at times greater than or equal to the direct arrival and is symmetric in time. It may appear counterintuitive to apply a window that removes all energy at these times because this is the data that we are ultimately trying to estimate in the Green’s functions. However, this stage of the Marchenko method estimates the focusing functions, and these functions only exist at times before the direct arrival and after the time-reversed direct arrival. Outside this window is where the Green’s function exists, but its accuracy is dependent on the accuracy of the focusing functions (by equations 1 and 2). Furthermore, this in itself is an approximation because we are assuming that the Green’s function and focusing functions can be separated by a windowing operator in the space-time domain, which is not always the case (e.g., when the focusing location is on or near a subsurface interface). Nevertheless, we work with these approximations and now iterate to a solution as follows.

In Figure 2b, three wave packets are labeled ($\alpha$, $\beta$, and $\gamma$): These are injected at the surface in addition to the initial impulsive source $\delta$ that is used to obtain the reflectivity in Figure 2a. These additional wave packets make up $M_k^{+}$, which is the coda (later part) of $f_k^{+}$, where $k$ is the number of iterations:

As a demonstration of how the focusing functions are estimated using equations 3 and 4 (and equations 6 and 7 below), Figures 3 and 4 show a series of raypath diagrams that explain their various traveltime relationships. These figures use a similar display format to that presented by van der Neut et al. (2015b). In Figure 3, the three primary reflections from the reflectivity are depicted individually (middle column) and convolved with the inverted direct arrival $f_0^{+}$ (left column) from equation 3. The main point of Figure 3 is to show that the results of this convolution are a series of nonphysical signals, each of which is contained within the pass window of $\theta$ and make up $f_0^{-}$. They are not physical because each signal on the right is made up of combinations of energy that has positive (solid) and negative (dashed) traveltimes, which survive the windowing operation in equation 3. Despite being nonphysical at this stage, they can be used to construct the focusing functions by progressing them to the next iteration:

To estimate $M_k^{+}$ using equation 6, we start with the estimate $f_{k-1}^{-}$ produced by the previous iteration (or the initial iteration $f_0^{-}$), time reverse it, and then convolve it with the reflectivity. The same windowing operator $\theta$ as above is then applied to isolate the focusing functions of interest.

A schematic of the first iteration ($k=1$) to calculate $M_1^{+}$ is shown in Figure 4. The columns on the left of Figure 4 show a subset of the (time reversed) columns on the right of Figure 3, a subset was selected (rows 2 and 3 from Figure 3 only) because these are the only components that contribute to $M_1^{+}$. This column is convolved across rows with primary reflections from $R$; again, we have only included a subset of these reflections because these are the components that contribute to $M_1^{+}$. The solutions to this convolution step are shown in the third column, and the time reversal of this result is given in the final column.

The final column represents $M_1^{+}$ and is composed of three signals each of which is made up of the time-reversed direct arrival (the left column in Figure 3) plus an additional time lag. This additional time lag is equal to the two-way traveltime through one or more subsurface layers. This traveltime information is what the Marchenko method requires to accurately account for internal multiples. To demonstrate this we have included the traveltimes $\alpha =-0.24\mathrm{s}$, $\beta =-0.16\mathrm{s}$, and $\gamma =0.16\mathrm{s}$ from Figure 2b and 2c in the column depicting $M_1^{+}$ in Figure 4. This shows that at $\mathrm{depth}=0\mathrm{m}$, we have calculated the traveltimes of the signals that we need to inject into the subsurface to destructively interfere with the multiply scattered components of the reflectivity so as to cancel them out above the focus point. We have therefore demonstrated that the additional components of the focusing functions that are used in Figure 2b to remove multiples from the seismic reflection data can be formed by convolving the data with itself (and with the direct wave estimate). Although the traveltimes of $\alpha$, $\beta$, and $\gamma$ are correct, the amplitudes of the energy constructed in iteration 1 do not cause perfect cancellation of the internal multiples; the amplitudes are corrected in subsequent iterations.

The waves in $M_k^{+}$ can then be injected into the subsurface and what would return to the surface can be calculated (the convolution with $R$); the results are windowed with $\theta$, and whatever remains is added to the estimate of $f_0^{-}$:

Equations 6 and 7 are iterated until the solutions for $M_k^{+}$ and $f_k^{-}$ in consecutive iterations have converged to stable values. When the solutions have converged, the total downgoing focusing function $f_k^{+}$ can be constructed by summing the inverted direct arrival $f_0^{+}$ and $M_k^{+}$ using equation 5.

It is worth noting that within equations 4–7, the quantity from the previous iteration is convolved with the reflectivity $R$, which in practice always contains a source term (no matter whether real or modeled data are used). To avoid iteratively convolving multiple source terms together, the source wavelet was initially deconvolved from the reflectivity. If the reflectivity was not deconvolved, the effective source wavelet would change with each iteration, so consecutive iterations would be inconsistent. An alternative approach that avoids deconvolution is illustrated in the 2D example below.

Green’s function estimation

Once calculated, the focusing functions can be used to estimate directionally decomposed Green’s functions ($G^{+/-}$) using equations 1 and 2; summing those two signals gives the total Green’s function $G(x_0,x_i,t)=G^{-}(x_0,x_i,t)+G^{+}(x_0,x_i,t)$. This Green’s function represents the signal that would have been measured at the surface if there had been a source at the image point (or vice versa, by source-receiver reciprocity).

Because this experiment is synthetic, we can test the accuracy of the Green’s function estimates by comparing them with the true Green’s function — one that is obtained by modeling an actual source at the image point, as shown in Figure 5. The final panels of Figure 5a and 5b show the estimated and true Green’s functions for two different subsurface image points. For both examples, they match in time, and amplitudes are correct for all arrivals. In the first two panels of Figure 5a and 5b, the Green’s functions are shown in decomposed form as obtained directly from equations 1 and 2: We observe well-separated events in the upgoing and downgoing components. It can be seen that there are no measured downgoing arrivals in Figure 5a, which is to be expected given that downgoing waves from a virtual source at the image point would have to be reflected back upward to be recorded at the surface; there are no interfaces below the image point in Figure 5a (as shown by the reflector locations in Figure 2), so no such reflection can occur. By contrast, the image point in Figure 5b lies between reflectors; therefore, the downgoing and upgoing signals contain arrivals.

Marchenko methods are only able to construct events that follow raypaths for which the energy was recorded in the original reflectivity $R$. Therefore, in the examples presented in Figure 5, we have only plotted the estimated Green’s function to a maximum time of $1.4\mathrm{s}$. This equates to times preceding the recording time of the reflectivity $R$ minus the traveltime of the direct arrival $T_d$. This shows that it is important to have sufficiently long recording times for Marchenko methods to be effective.

So far, we have discussed a simple, 1D example that demonstrates the methodology clearly and in which the solutions are essentially perfect. The next section extends the examples to higher dimensions, where the results are more prone to errors.

MARCHENKO METHODS IN HIGHER DIMENSIONS

The example above illustrated for 1D problems that all of the information required to determine Green’s functions between a subsurface image point and the surface is contained within just two signals $R(x_0,x_0,t)$ and $T_d(x_i,x_0,t)$. In two dimensions, this is not the case because the reflectivity from one surface source (${\mathbf{x}}_0^{\prime }$) is measured by multiple receivers (${\mathbf{x}}_0$). It is worth noting that our notation must now change to account for the extra spatial coordinate, where $\mathbf{x}=(x_1,x_2)$. Nevertheless, in two or three dimensions, concepts similar to those in Figures 3 and 4 hold. Indeed, although for 1D problems, those diagrams have axes of depth and time; they apply with similar geometries (but incorrect angles of transmission and reflection) to 2D problems if the horizontal time axis is replaced with the horizontal space axis. Each arrival in the desired Green’s function at any particular angle is constructed by the interference of other specific arrivals at other particular receiver and source combinations in the reflectivity.

Rather than requiring that we selected specific arrivals to convolve to construct each arrival in the Green’s function, Marchenko methods sum (integrate) over all possible sources along the acquisition array (boundary) and rely on destructive interference to cancel out unwanted energy. A similar cancellation occurs in RTM (Kaelin and Guitton, 2006) and in seismic interferometry (van Manen et al., 2005, 2006; Wapenaar and Fokkema, 2006). This only works accurately in two dimensions if the reflectivity is of the correct form, and in practice, this means that a vertical spatial derivative (often called a dipole) source or receiver needs to be created or measured. In the 2D examples presented in this paper, the reflectivity is from a pressure (monopole) source measured by a vertical particle velocity (dipole) receiver.

In the previous section, we introduced a set of formulas to estimate Green’s functions using Marchenko methods. These formulas are extended to two dimensions by changing the 1D convolutions to multidimensional convolutions and integrating across all sources on the acquisition boundary (we denote this boundary as $\partial {\mathbb{D}}_0$). For example, the 2D versions of equations 1 and 2 for source redatuming are

To implement equations 8 and 9, the focusing functions need to be available between the focusing location and the surface sources and receivers (e.g., $f^{+}({\mathbf{x}}_0,{\mathbf{x}}_i,t)$ and $f^{+}({\mathbf{x}}_0^{\prime },{\mathbf{x}}_i,t$)). In practice, the 2D formulation of these functions requires interchangeability between the two. We therefore impose the condition that the source array and receiver array are colocated ${\mathbf{x}}_0^{\prime }={\mathbf{x}}_0$.

A further consideration for implementation of the 2D Marchenko method is the direct arrival estimate $T_d$ and windowing function $\theta$, which now need to be estimated in two dimensions, as they were above in one dimension. These functions are now calculated between a single focusing location and multiple surface sources/receivers. This increases the complexity of these functions and the potential for errors in estimating them; nevertheless, if this is done following the same workflow as 1D Marchenko methods, the relationships discussed in the previous sections still hold.

Green’s function estimation in two dimensions

In the following 2D example, the source wavelet in the reflectivity is not the same as that in the focusing functions. We use the $20\mathrm{Hz}$ Ricker wavelet shown in Figure 6b for the focusing function, whereas we use a flat spectrum wavelet shown in Figure 6a as the source wavelet for our reflectivity. The flat spectrum wavelet is defined in the frequency domain so as to have an amplitude of one over the range of frequencies of interest (the frequencies contained within the Ricker wavelet), as demonstrated in Figure 6c. Using this formulation removes the need for deconvolution because it ensures that the frequency content of the updated focusing functions does not change between iterations of the Marchenko method (in each iteration, they are multiplied by a source wavelet that does not change the shape of the current wavelet within the frequency band of interest [Thorbecke et al., 2017]).

Figure 7 shows a subsurface model that is used to demonstrate 2D Marchenko methods. This subsurface model has variable density and velocity. There are 188 symmetrically spread sources and receivers colocated at 16 m intervals along the surface of the model (depth = 0 m). A point is also marked in the subsurface, ${\mathbf{x}}_i=(1000\mathrm{m},800\mathrm{m}$), which identifies a chosen image point for Marchenko Green’s function calculation.

Two models have been used to create the two input data sets required for the Marchenko method. The reflectivity from surface sources measured at the surface receivers and excluding free-surface multiples has been created using the true model given in Figure 7a. This represents surface seismic reflection data after free-surface multiples have been removed. An estimate of the direct arrival between the surface and the image point has been created using the smooth model in Figure 7b. In practice, we do not have access to the true model; hence, we have used a smoothed version of the true model for our direct arrival calculations, assuming that in practice some initial or reference velocity model will be available. The direct arrival signal can be modeled with finite-difference solvers (Figure 7c), or approximated using eikonal solvers to find the traveltime at which a scaled source wavelet can be assumed to arrive. In this example, for accuracy, both data sets were created using finite-difference solutions to the acoustic wave equation (in Figure 7b, a source was fired at the image point and recorded along the surface receiver array, giving Figure 7c). See van der Neut and Wapenaar (2016) for an alternative solution to solving the Marchenko method if an estimated velocity model is not available.

By iterating the 2D form of equations 1–7 (Wapenaar et al., 2014), the focusing functions and Green’s functions are obtained, and the final Green’s function estimates are shown in Figure 8b. For simplicity, we have not included every component of the estimated Marchenko Green’s function (e.g., focusing functions) — for these, we refer readers to the accompanying code package within which these figures are included.

To test the accuracy of the Marchenko solution, the calculated Green’s functions in Figure 8b are plotted beside the true solutions in Figure 8a. The latter panel shows a solution computed in the true model in Figure 7a by firing a source at the image point. The estimated signal shows a good match, with negligible errors visible. The errors that are present can be attributed to limited boundary coverage by the acquisition array, errors in the finite-difference solution, and windowing artifacts.

In Figure 8, all of the amplitudes have been scaled to values between 1 and $-1$. This has been done for comparison purposes because the Marchenko methods implemented in this paper cannot estimate true absolute amplitudes of Green’s functions. The primary reason for this is that the direct arrival was approximated at the start of the Marchenko method as $f_0^{+}({\mathbf{x}}_0,{\mathbf{x}}_i,t)\approx T_d({\mathbf{x}}_i,{\mathbf{x}}_0,-t)$: We do not know the amplitude of the true inverse, so it is impossible to estimate a Marchenko solution with the true absolute amplitude because the initial focusing function estimate is implicit in the final solution (see equation 5).

MARCHENKO CODE PACKAGE

Accompanying this paper is a set of well-commented MATLAB codes for 2D Green’s function estimation. The first of these (CODE_1) is the code used to create Figure 8, and running it without edits should produce a version of that figure. For this code, the inputs are precomputed and the variables already set; below we discuss the operation of this. However, because that code is inflexible (the image point is fixed), we have included a second code for user experimentation, which we introduce at the end of this section. By choosing an appropriate setup, users should be able to use the second code to reproduce similar outputs of the first — a good exercise for learning and teaching purposes.

Data

Four data sets accompany this MATLAB code, as summarized in Table 1. All of these data sets are stored in the time domain, although for computational efficiency, most of the code operations are performed in the frequency domain.

The direct arrival, windowing function, and true signals are all common shot gathers from a source at the image point with identical dimensions: $i=2001$ and $j=188$, where $i$ is the number of time samples and $j$ is the number of surface receivers. In this example, the sampling interval was $0.002\mathrm{s}$ and the recording time was from $-2$ to $2\mathrm{s}$. The fourth data set is the reflectivity, which has an extra dimension $k$ to account for multiple shot locations, with $i=2001$, $j=188$, and $k=188$.

Codes

The code called ICCR_marchenko.m, follows the same workflow introduced earlier in this paper. Algorithm 1 shows the corresponding pseudocode. The equations referred to in Algorithm 1 are the 1D versions given herein, but the code uses the equivalent 2D versions given in Wapenaar et al. (2014). We have not defined a measure of convergence within this code; instead, a desired number of iterations is input by the operator, the default being five.

One additional feature in the MATLAB code has not been discussed above: A spatial taper is applied during each summation (integration) of sources along the acquisition boundary. Its purpose is to account for the limited acquisition coverage: It ensures cancellation of edge effects from the extremities of the acquisition array because these would otherwise create spurious energy in the solutions (if the boundary was infinite, as assumed in 2D or 3D Marchenko theory, this would not be required). The taper takes a half cosine shape at either end of the array, and the number of points to be tapered at either end of the array can be varied by the user (the default is 20% of the number of receivers).

A second code in the accompanying package (CODE_2), operates using the same fundamentals as the code introduced above, but it has been implemented as matrix multiplications for computational efficiency and changed to allow alternative virtual source locations to be used. To the latter end, the direct arrival and window are no longer precomputed; instead, we have included a function that computes these using an eikonal solution (ICCR_marchenko_eik.mat). The input seismic data remains the same, but there is no longer a true solution available for comparison. This code is set up so that the input data can be changed; to do this, a measured or calculated reflectivity data set in the correct form will be required as well as an eikonal solution through an estimated velocity model.

MARCHENKO IMAGING

Marchenko imaging creates images of the subsurface using the estimated Green’s functions. This works in a similar way to conventional imaging algorithms such as RTM in the sense that the similarity of two signals is tested by the use of an imaging condition (Claerbout, 1971), and the result is used to identify subsurface inhomogeneities.

Marchenko imaging should in theory be able to be used to calculate more accurate images than standard methods because we have access to more accurate Green’s functions $G^{+/-}$, which account for multiply reflected waves. The imaging condition applied for our implementation of Marchenko imaging is a (zero-lag) crosscorrelation between the downgoing Green’s function $G^{+}$ and the direct arrival estimate $T_d$, as proposed by da Costa Filho et al. (2015):

For comparison purposes, we have also implemented an alternative imaging condition to approximate standard methods:

where $G_0^{+}\approx R\otimes T_d^{*}$, which is equivalent to a back-propagated wavefield used in RTM methods (for further details, see da Costa Filho and Curtis, 2016).

The two signals in equations 10 and 11 should be most similar when the image point ${\mathbf{x}}_i$ is on a reflector, so the crosscorrelation will produce maxima at those points. We note that as with conventional RTM, there are several alternative imaging conditions that could be applied. Comparisons and discussion of these is beyond the scope of this paper, but more detail is given in da Costa Filho et al. (2015) and Singh and Snieder (2017b).

We imaged the model in Figure 7a using equations 10 and 11, with results shown in Figure 9. Image points were selected on a 4 m grid inside the imaged area. Directionally decomposed Green’s functions were calculated at each of these points. The only difference between these Green’s functions and those calculated in Figure 8 is in the direct arrival input, which here was constructed by placing a source wavelet at the traveltime calculated using an eikonal solver, rather than calculating $T_d$ using finite-difference methods; this saves on computational cost because eikonal solvers require relatively few operations compared with finite-difference methods. This example therefore also illustrates that this approximate method of modeling the direct field can be sufficient for some imaging applications.

The image calculated using equation 10 at each image point is shown in Figure 9a. All interfaces in the subsurface are identified with few artifacts in the solution, despite the presence of internal multiples in the surface-acquisition data that is input to the Marchenko method. The clarity of the image is mainly due to the accuracy of the Green’s function estimates $G^{+}$, which are calculated by Marchenko methods and used by equation 10. For comparison, the solution in Figure 9b is calculated using equation 11 with the approximate $G_0^{+}$ shows clear artifacts due to internal multiples.

We have included a final code (CODE_3) in the software package, which can be used to implement Marchenko imaging using the methods discussed above. There is an increased computational cost associated with implementing this because a Green’s function now needs to be calculated at each and every imaging point. Inside the code, the image can be target-oriented by defining the image point spacing and a limited or targeted subsurface location.

DISCUSSION

The aim of this tutorial is to provide beginners to Marchenko methods with an accessible way to understand the topic and to allow them to begin to experiment with the methods on synthetic examples using easily understandable and editable MATLAB code. Marchenko methods as introduced in this paper use surface seismic data and an estimate of the subsurface velocities to estimate Green’s functions between a subsurface image point and the surface. These estimated signals have a steadily increasing number of applications, including subsurface imaging and seismic redatuming, but also for removing multiples (Meles et al., 2015), constructing primaries (Meles et al., 2016), and calculating Green’s functions, where the source and receiver are inside the subsurface (Wapenaar et al., 2016; Singh and Snieder, 2017a). Marchenko methods also have applications outside of seismology, for example, for imaging using ground-penetrating radar (Slob et al., 2014a) or for medical imaging (van der Neut et al., 2017).

The results are promising and offer a data-driven method that improves on current imaging methods, in particular by correctly predicting the arrival of multiply reflected waves at image points. Given the novelty of these methods, there are still aspects that are poorly understood — areas of ongoing research. They include exploring how to apply Marchenko in real, dissipative media with seismic attenuation, the effects on Green’s function estimates of velocity model errors and corresponding poor estimates of direct arrivals (in time and amplitude), the effect of various types of noise in the reflectivity field, and the cost and effort of scaling 2D Marchenko methods to three dimensions.

The accompanying MATLAB codes for the estimation of acoustic Green’s functions using 2D Marchenko methods also comes with an accompanying data set, which can be used to reconstruct Figure 8. The codes can easily be adapted for users to include their own data sets. The data will need to be formatted correctly, the details for which are discussed in the earlier sections of this tutorial and in comments within the code. Data sets of a similar simplicity should achieve similarly positive results.

CONCLUSION

In this paper, we have introduced Marchenko methods, a set of novel, data-driven techniques that can be applied to seismic redatuming and imaging problems. We have shown that these methods can accurately estimate directionally decomposed Green’s functions from virtual subsurface source locations to surface receiver locations in one and two dimensions, and this includes the multiply scattered components of the Green’s functions. However, all of the methods we have presented are based on synthetic seismic data sets; extending these methods to more realistic data sets and examples is an area of active research, and with the accompanying MATLAB code, readers have the tools to contribute to this.

DATA AND MATERIALS AVAILABILITY

The data and code used for this paper can be downloaded from http://software.seg.org/2019/0002.