Reconstructing the details of subsurface structures deep beneath complex overburden structures, such as subsalt, remains a challenge for seismic imaging. Over the past few years, the Marchenko redatuming approach has proven to reliably retrieve full-wavefield information in the presence of complex overburden effects. When used for redatuming, current practical Marchenko schemes cannot make use of a priori subsurface models with sharp contrasts because of their requirements regarding initial focusing functions, which for sufficiently complex media can result in redatumed fields with significant waveform inaccuracies. Using a scattering framework, we evaluate an alternative form of the Marchenko representation that aims at retrieving only the unknown perturbations to focusing functions and redatumed fields. From this framework, we have developed a two-step practical focusing-based redatuming scheme that first solves an inverse problem for the background focusing functions, which are then used to estimate the perturbations to focusing functions and redatumed fields. In our scheme, initial focusing functions are significantly different from previous approaches because they contain complex waveforms encoding the full transmission response of the a priori model. Our goal is the handling of not only highly complex media but also realistic data — band-limited, unevenly sampled, free-surface-multiple contaminated data. To that end, we combine the versatility of Rayleigh-Marchenko redatuming with our scattering-based scheme allowing an extended version of the method able to handle single-sided band-limited multicomponent data. This scattering-Rayleigh-Marchenko strategy accurately retrieves wavefields while requiring minimum preprocessing of the data. In support of the new methods, we evaluate a comprehensive set of numerical tests using a complex 2D subsalt model. Our numerical results indicate that the scattering approaches retrieve accurate redatumed fields that appropriately account for the complexity of the a priori model. We find that the improvements in wavefield retrieval translate into measurable improvements in our subsalt images.


Accurate estimation of the stratigraphy and properties of complex subsurface geology has long represented a challenge for seismic imaging, especially in the presence of salt (or basalt) formations in the overburden (Leveille et al., 2011). This is largely due to uneven illumination of the target area arising from the complex propagation in the overburden as well as the presence of strong reverberations in the recorded surface data (Jones and Davison, 2014). On the one hand, advances in seismic acquisition techniques, such as the adoption of wide-azimuth surveys, have played a major role in providing more even subsurface illumination, which in turn can help to reduce the ill-posed nature of the associated imaging problem. On the other hand, it is only by incorporating better physics in the imaging process that the full potential of the recorded seismic data can be exploited. This is for example justified by the uplift in the quality and focusing of reflectors provided by high-end, wave-equation-based migration methods such as reverse time migration (Baysal et al., 1983; McMechan, 1983; Farmer et al., 2009), data- or image-domain least-squares migration (Nemeth et al., 1999; Fletcher et al., 2016; Arasanipalai et al., 2019), and full-waveform inversion (Tarantola, 1984; Virieux and Operto, 2009; Esser et al., 2016; Wang et al., 2019). Moreover, a variety of approaches, emerged under the pursuit of imaging with multiples, have also shown their effectiveness in such settings; these methods turn multiple reverberations originating from high impedance contrasts around salt and basalt structures into useful signal that can complement the illumination of primaries (Malcolm et al., 2008; Liu et al., 2015).

Alongside model-driven approaches that aim at iteratively reconstructing the velocity (and/or reflectivity) model of the subsurface, from which a more accurate propagation of the underlying wavefield can be numerically modeled, a great deal of research has been recently devoted to target-oriented methods. By targeting a specific area of the subsurface, these methods promise to provide reservoir-characterization-ready models of the subsurface (da Costa et al., 2019; Cui et al., 2020; Guo and Alkhalifah, 2020). The reduced computational domain and size of the model to invert for can allow for the inclusions of higher frequencies and even more complex physics (e.g., elastic, attenuation) in the modeling operator used to describe the data. Nevertheless, the quality of the final image and/or inverted properties is strongly dependent on the ability to accurately redatum surface data at the target level of interest. In this regard, novel wavefield extrapolation techniques that go beyond the single-scattering (i.e., Born) approximation open up new ways to create such target data.

Marchenko redatuming (Wapenaar et al., 2014b; van der Neut et al., 2015b) has recently been demonstrated to be a reliable method for retrieving full-wavefield subsurface responses in a variety of geologic settings, which can accommodate moderately complex overburden geology. The Marchenko method makes use of the recorded reflection response and an estimate of the direct transmission response in the medium to estimate so-called focusing functions that serve as operators to focus energy at specific points in the subsurface. Once the redatuming step has been performed, the resulting subsurface wavefields can be used to create artifact-free images of the subsurface (Wapenaar et al., 2014b) or to estimate local properties (Cui et al., 2020) by naturally including all orders of multiples present in seismic reflection data. From this perspective, Marchenko redatuming represents an ideal platform for creating a subsurface wavefield suitable for target-oriented imaging and inversion beneath complex structures such as salt and other environments such as subbasalt or other highly complex tectonic settings.

To that end, however, the original Marchenko scheme has strict requirements with regard to the reflection response; namely, it requires large aperture acquisition geometries, it should represent data from a surface-multiple-free survey, and it should have a flat broadband frequency content. When it comes to practical implementations of the method, the requirement for an accurate deconvolution of the source wavelet, removal of surface-related multiples (SRMs), and source/receiver colocation hinders its application to real data sets. Moreover, treating SRMs as noise that must be removed from the data is not only technically challenging and time consuming but it may also not be beneficial in complex media. Surface multiples can carry complementary information compared with primaries because, when compared with surface-multiple-free data, they are exposed to longer propagation times and possibly different propagation paths in the subsurface as a consequence of the natural interactions of finite-aperture data with a free surface in a highly complex medium. Furthermore, on the subject of data preprocessing for Marchenko redatuming, source deconvolution with an erroneously estimated source wavelet leads to strong coherent artifacts in local images (Mildner et al., 2019b). In many ways, the success of single-sided Green’s functions retrieval depends on the availability of accurate focusing functions. Dukalski et al. (2019) describe how to correct for the effect of short-period multiples in the Marchenko framework, by adding energy conservation constraints to the focusing problem. To relax on the originally strict acquisition requirements, Rayleigh-Marchenko redatuming (R-Marchenko) was introduced by Ravasi (2017) and Slob and Wapenaar (2017). The main advantage of R-Marchenko is the use of a band-limited operator defined in terms of vertical particle velocity rather than a broadband reflection response, as is the case in standard Marchenko implementations. Such consideration significantly reduces the need for most preprocessing steps.

More importantly in the context of this work, previous versions of the Marchenko scheme could not accommodate highly complex media such as subsalt — particularly in the presence of an a priori model containing sharp contrasts. This is because currently available Marchenko schemes explicitly require an estimate of a single direct arrival from a chosen focal point to the survey surface. In most approaches so far, this estimation is either based on an initial smooth velocity model (no sharp contrasts) or by selecting/windowing first arrivals obtained by modeling, thus not being able to fully use complete waveform information contained in the modeling through complex velocity models with high (or low) acoustic impedance inclusions, such as in the case of salt bodies. Despite initial reports of successful attempts imaging complex subsalt structures (Jia et al., 2018; Staring et al., 2018; Staring and Wapenaar, 2020), Vasconcelos et al. (2014) show that the accuracy of conventional Marchenko redatuming can be considerably lower when dealing with highly complex media. To better reproduce seismic subsurface responses in such conditions, Vasconcelos and Sripanich (2019) propose a modified Marchenko scheme that can incorporate available information from previously estimated, complex subsurface models by introducing an auxiliary Marchenko system of equations for a reference model as a constraint to the solution in the true medium.

In this paper, we first build on the work of Vasconcelos and Sripanich (2019) and present a scattering-based Marchenko framework designed to handle highly complex media, while making full use of state-of-the-art velocity models with sharp contrasts. To account for band-limited data in the presence of free-surface multiples and with realistic acquisition geometries, we combine the versatility of R-Marchenko with the scattering-based Marchenko’s ability to handle redatuming in complex media. This enables us to derive a new scattering-based Rayleigh-Marchenko method (SR-Marchenko) that shows a data-robust, superior performance to previously developed approaches in highly complex media. The new strategy naturally accounts for the compatibility between real and modeled data in the perturbation operator, which in turn allows for rather direct application on field multicomponent data sets. To evaluate the effectiveness of the new method, we present a comprehensive set of numerical tests using a complex 2D subsalt model. The focusing and Green’s functions produced by the SR-Marchenko, R-Marchenko, and the original Marchenko schemes are compared against each other as well as against a numerically modeled pressure field. Such wavefields are further used as input for a step of source redatuming via multidimensional deconvolution (MDD) followed by localized imaging by means of reverse time migration (RTM). At all stages of redatuming, we observe that the wavefields retrieved from our new approach are superior in terms of event continuity and amplitude fidelity, suggesting that wavefield redatuming and imaging in complex media by means of Marchenko equations can greatly benefit from further information on the underlying subsurface media, i.e, migration velocity models.


In this section, we review the main aspects involved in the traditional implementations of single-sided wavefield focusing and extrapolation, followed by an extension of Marchenko redatuming accounting for scattering effects in connection with an a priori reference model. Contrary to the current practice of approximating the initial focusing function by a time-reversed direct arrival, this strategy offers the possibility of estimating a reference-medium inverse transmission response, accounting for the effects of geometric spreading and the kinematics of primaries as well as internal multiples in the given model. This more accurate focusing approach positively impacts the subsequent redatuming step, which in turn benefits the resulting imaging process. However, estimating the full transmission response is sensitive to potential inaccuracies in the reference model. Such errors may lead to a time shift in the retrieved focusing functions as well as the arising of artifacts, which ultimately influence the quality of the retrieved Green’s functions (Thorbecke et al., 2013; Broggini et al., 2014). Further analysis of this subject and its implications for redatuming are addressed in the “Discussion” section.

Marchenko focusing and redatuming

One-way reciprocity theorems (Wapenaar and Grimbergen, 1996) provide the basis to estimate Green’s functions, g(xr,xf;t), from a virtual source, xf, at a depth level Di to a set of receivers xr located along the surface D0 using a single-sided focusing operator conventionally called the focusing function f(xr,xf;t) (Wapenaar et al., 2014a). Such a wavefield is defined in a reference configuration (i.e., a truncated medium) and collapses at zero time at a given depth level when injected from the surface. The original Marchenko framework is restricted to a lossless medium, in this case, supporting only acoustic waves, and it does not consider evanescent wave modes at either the acquisition or focusing levels (Wapenaar et al., 2004). According to Wapenaar et al. (2014b), the Marchenko integral equations relate the focusing function in the truncated medium with the true medium’s Green’s function, using the measured reflection response as the kernel of a multidimensional convolution (MDC). In theory, the reflection response R(xr,xr;t) is a full-bandwidth field that does not include surface-related multiple (SRM) effects and is, in general, assumed to be known. This procedure, in the frequency domain, is described by
where the frequency dependence has been omitted for notational brevity. Here, * represents the complex conjugation, i.e., time reversal in the time domain. A typical Marchenko setup consists of a set of sources xr and receivers xr collocated on a transparent surface D0. Similarly, the virtual source location xf lays on an arbitrary focusing depth level Di. The required two-way response from xf to xr is given in terms of its up- and downgoing components g and g+ according to g=g+g+. However, the up- and downgoing focusing functions f± are defined in a truncated version of the medium designed to be reflection-free outside the region between levels D0 and Di.
From a practical perspective, sources and receivers are positioned at discrete locations. Therefore, it is useful to discretize the Marchenko system of equations 1 and 2 and recast them in a compact matrix-vector notation. In the discrete Marchenko framework (van der Neut et al., 2015b), the reflection response acts as a multidimensional filter on the focusing operators returning the desired Green’s functions, such that
denotes the discrete version of the Marchenko system. The Marchenko-operator matrix is composed of an identity matrix I together with matrix operators R and R*, which apply multidimensional space-time convolution and correlation — note that this operator form applies to time- and frequency-domain versions of the system, although, of course, the numerical structure of the operators themselves differ. As per the focusing f± and Green’s g± functions, they are concatenated in block-vector form (i.e., array stacking) in time space to form the vectors of the system. Given that equations 1 and 2 represent an underdetermined system of two equations with four unknown functions (g, g+, f, and f+), additional constraints need to be introduced to solve such a system. In the particular case of relatively smooth media with moderately curved interbeds and limited acquisition aperture, it is plausible to assume a causality argument stating that focusing and Green’s functions lay on different space-time regions. Provided that this condition holds, g and f do not mutually interfere with one another (Wapenaar et al., 2014b). To numerically enforce causality, we implement a window matrix Ψ designed to be symmetric in time such that it removes all causal and acausal events within the time interval td(xr,xf)ttd(xr,xf), i.e., Ψg=0 and Ψg+=0. However, the downgoing focusing function can be seen as the composition of a direct wave followed by a coda containing all different orders of scattering that arrive at later times f+=fd++fm+. The action of the window matrix on f+ only removes its direct wave, i.e., Ψf+=fm+, and it does not have any effect on its upgoing counterpart Ψf=f. In light of these arguments, the number of unknowns in equations 1 and 2 is reduced once Ψ is applied to both sides of the equations (van der Neut et al., 2015b). Finally, a simplified matrix-vector representation for the discrete case follows after rearranging the remaining terms:

Evaluation of equation 4 requires the direct part of the downgoing focusing function (fd+) to be known. In principle, if the direct part of the transmission response in the truncated reference medium Td is available, fd+ can be obtained from the focusing condition, i=Tdfd+, by direct inversion (van der Neut et al., 2015b). Here, vector i represents a band-limited delta function at the focal point xf at time zero. From a practical standpoint, instead of inverting for the initial focusing function, one common choice is to approximate it by taking the time reversal of the modeled transmitted direct Green’s function, i.e., fd+gd*. This approximation demands some caution because it mostly accounts for kinematics and does not properly account for more subtle amplitude effects such as geometric spreading, turning waves, or attenuation (Wapenaar et al., 2014a).

Redatuming of scattered fields

Estimating focusing and Green’s functions from single-sided Marchenko redatuming has proven to be challenging in highly complex media. Although the performance of Marchenko redatuming compared to conventional reverse time extrapolation has been shown to be superior, several reconstructed seismic events are missing for subsurface characterized by sharp impedance contrasts, boulders, dipping layers, or strong diffractors (Vasconcelos et al., 2014). A major drawback of this approach is that some difficulties arise when trying to reproduce relative and absolute amplitudes, as well as retrieving all possible scattered arrivals. However, in such environments, greater emphasis is generally placed on building a seismic velocity model because an accurate knowledge of the velocity field becomes even more important for successful migration results. For example, high-end migration models that can capture valuable information regarding the salt’s sharp structure, provide additional constraints that may be used for purposes other than the migration itself. Provided that such models are available, our goal is to take full advantage of them to retrieve full-wave scattered fields for redatuming and imaging purposes. To that end, here we extend the Marchenko framework to reconstruct perturbed focusing functions with respect to a background medium, and as a result, scattered Green’s functions are also retrieved.

Our derivation starts by considering two independent Marchenko integral systems in the matrix-vector form (equation 3), the former describing the relationship between focusing and broadband seismic responses in the true medium, and the latter providing the connections between corresponding fields in a given reference (migration) model,
where R represents the real medium’s reflection response at the transparent surface and R0 represents the corresponding reflectivity operator in the reference (e.g., migration velocity) model. Again, * in the superscript denotes the complex conjugation, that is, the reverse time version of the corresponding operator. For the sake of notation, both systems of equation 5 can be written in a more compact form as g=Mf and g0=M0f0, where M and M0 are the real- and reference-medium Marchenko operators, respectively. Given the known background model, the quantities g0 and R0 are calculated by forward modeling, so that one may have access to f0 by inversion. Then, we introduce the perturbation quantities δg=gg0, and δf=ff0, which represent the scattered-wavefield effects in g and f as a result of the differences between the real and reference media. By subtracting the Marchenko system in the reference from that in the actual medium, an extended Marchenko system for scattered fields is derived:
As can be noted, this new representation δg=Mδf+δMf0 describes the relationship between the unknown Green’s function perturbation δg, in terms of R, R0, and f0, where one identifies the perturbation operator δM=MM0 as the fundamental element reflecting discrepancies between the two media. We referred to such formulation as scattering Marchenko or S-Marchenko. In real applications, care must be taken into account when considering the proper scaling between the two reflection responses. Solving for δf+ and δf is done by restricting the space-time support to the interval td(xr,xf)ttd(xr,xf), where the windowing operator Ψ separates causal from acausal events. This operator admits all events before the first arrival in δg±, whereas it eliminates any other event mapping outside the window (Wapenaar et al., 2014b). Introducing this causality constraint allows us to recast equation 6 into a system that relates focusing-function perturbations δf only to the focusing functions in the background f0, i.e.,
A simplified form of equation 7 reads δMf0=Mδf, where M and δM take into account the influence of the windowing function Ψ acting on R. Finally, the aforementioned equation is solved as the constrained optimization problem,
minδfMδfδMf022s.t.  g0=M0f0.
From a practical perspective, numerical implementation is done by following a two-step, alternated optimization process. First, full-waveform forward modeling in the reference model is implemented to reconstruct R0 and g0, followed by inversion of the full (i.e., without windowing) Marchenko system in the reference model for an estimate of f0. Then, the true-medium response R together with the modeled reference-model response R0, and f0 from the previous step, are used to invert for δf.


The idea of reproducing full-wavefield responses from virtual sources at depth using focusing operators as presented in the previous section requires several key preprocessing steps that may introduce inaccuracies or even hinder the implementation of Marchenko redatuming in real data sets (Ravasi et al., 2016; da Costa Filho et al., 2017; Jia et al., 2018; Mildner et al., 2019a). The original approach demands knowledge of the acquisition wavelet and dense spatial distribution of sources and receivers. In effect, the surface reflection response needs to be accurately deconvolved by the source wavelet. As Mildner et al. (2019b) show, an erroneous removal of the source fingerprint leads to artifacts in local images based on redatumed data. In addition, the formalism presented above does not consider the effects of the free surface; as a result, SRM elimination (SRME) (Verschuur et al., 1992; van Borselen et al., 1996; Amundsen, 2001) is a requirement. Although most field data sets require 3D SRME for reasonable multiple suppression (Dragoset et al., 2010), it is also well-known that this process may introduce distortions to target signals, thus affecting data quality before redatuming. Alternatively, in the context of Marchenko redatuming, a natural way to circumvent the heavy toll of preprocessing consists of transforming the up- and downgoing Green’s functions (equations 1 and 2), with a band-limited aerial source via MDC (equation 9). Such a source can seamlessly account for the effects of all orders of internal and SRM enclosed in the downgoing vertical particle velocity measured at the receiver side (Ravasi, 2017).

Rayleigh Marchenko

Provided the availability of dual-sensor data, one can use the recorded pressure p(xr,xs;t) to decompose the vertical particle velocity vz(xr,xs;t) into its up- and downgoing constituents vz±(xr,xs;t) (Wapenaar, 1998). In marine seismic data studies, and in particular for ocean-bottom-cable or -node surveys, such fields are emitted from sources at a different level Ds above the receiver line and below a free surface (see Figure 1). Ravasi (2017) uses vz+ as the aerial source at the receiver level, which then propagates into the medium with the help of the Green’s functions, a process that grants access to the pressure field that would have been measured at Ds as if there was a source at the focal point xf. This convolutional operation is described in equation 9. What is also interesting to note is that vz+(xr,xs;t) and vz(xr,xs;t) are related to one another through the reflection response R(xr,xr) at the receiver level, as indicated by the following Rayleigh integral representation (equation 10):
The advantage of considering vertical particle velocities is that they naturally contain the effects of SRM while preserving the band-limited character of observed seismic fields. Contrary to the data required by the reflection operator R (which demand accurate absolute-value scaling), vz± fields are naturally scaled with respect to each other; therefore, estimation of an absolute scaling factor is not required in practical implementations (van der Neut et al., 2015c). Introducing equations 1 and 2 into equation 9 with the help of equation 10 eliminates the contribution of the reflection response R(xr,xr), yielding
The set of equations 11 and 12 represent the Rayleigh-Marchenko formulation for single-sided redatuming problems (Ravasi, 2017), hereafter referred to as R-Marchenko. An alternative derivation is found directly from one-way reciprocity theorems of the convolution and correlation type with interaction quantity z{PA+Vz,B+PAVz,B+} (Wapenaar and Grimbergen, 1996; Slob and Wapenaar, 2017). Referring back to the discrete formulation of equations 1 and 2, the R-Marchenko system can also be expressed in matrix-vector notation by defining Vz± as multidimensional discrete filters applied to the focusing functions f± to reconstruct, p±, the up- and downgoing subsurface pressure wavefields at the focusing level Di. Without loss of generality, it is written as
Solving the system of equation 13 relies on a causality argument along the lines of the one stated in the previous section. A new window Θ is built to remove all causal and acausal events within the time interval td(xs,xf)ttd(xs,xf), i.e., Θp=0 and Θp+=0. Application of this muting function to both sides of equations 11 and 12 allows the focusing functions to be retrieved from the following system:
where fm+ and f are evaluated in terms of an initial focusing function typically replaced by the direct part of a seismic response obtained from a given smooth model, fd+gd*. The focusing process for Marchenko and R-Marchenko schemes depends on the correct estimation of gd*. Depending on the overburden complexity in terms of fine-layer impedance variations, additional amplitude and phase corrections may be accounted for in a second focusing step as a way to compensate for the effects of short-period multiples (Dukalski et al., 2019).

The one-way focusing functions in equation 4 can be found by solving the coupled Marchenko equations either by Neumann-series iterative substitution (e.g., Wapenaar et al., 2014b; van der Neut et al., 2015b) or by direct inversion (iterative or otherwise, e.g., van der Neut et al., 2015a; Ravasi, 2017). Note that equation 4 is of the form y=[IA]x, a Fredholm integral equation of the second kind, which is amenable to solutions in terms of Neumann series expansion. A detailed analysis of this kind of solution is offered by Dukalski and de Vos (2017). It is shown that for strong scattering regimes, and in the presence of SRMs, convergence is not necessarily guaranteed. More broadly speaking, the convergence of such series depends on whether or not the spectral radius of operator A is bounded by unity. Unlike the standard Marchenko equations, the Rayleigh Marchenko system is a Fredholm integral of the first kind; in consequence, it cannot be expanded into a Neumann series. Without additional constraints or conditions, it is not possible to solve the R-Marchenko for the focusing functions using iterative substitution; thus, direct inversion must be used instead. Either equation 4 or 14 can be treated as a linear inverse problem (van der Neut et al., 2015a) in which iterative solvers (e.g., least-squares with QR decomposition [LSQR] — Paige and Saunders, 1982) prove to converge to a solution whose residuals decay monotonically, making such solvers an attractive option for Marchenko problems. Furthermore, an advantage of using direct inversion lies in the fact that data do not necessarily need to be complete. Indeed, even in cases in which it is irregularly sampled or polluted with noise, one may impose constraints in such a way that a satisfactory solution can still be estimated (Haindl et al., 2021). In such cases, experience has shown that sparsity-promoting inversion (Hennenfent and Herrmann, 2008; Beck and Teboulle, 2009; van den Berg and Friedlander, 2009) may be even more suitable than least-squares inversion — provided that a suitable sparsifying basis transformation is known and available (Haindl et al., 2021).

Scattering Rayleigh Marchenko

Building on the ideas of Vasconcelos and Sripanich (2019), we present an extended formulation of the previously described R-Marchenko method. This framework aims at reconstructing perturbed focusing functions with respect to a background medium from band-limited data while preserving the effects of the SRM while also allowing for substantial flexibility in terms of acquisition geometries. As a result, the corresponding scattered pressure fields at depth are retrieved. Even though the R-Marchenko scheme relaxes some of the acquisition and preprocessing requirements, it still suffers from similar difficulties that traditional Marchenko endures when it comes to wavefield extrapolation in complex media. In subsurface characterized by geologic intrusions exhibiting high impedance contrast, as is the case of diapiric traps, even a fairly accurate estimation of a direct wave as a proxy for the initial focusing function in Marchenko may not be sufficient to guarantee its success. Once again, we rely on the availability of high-end migration models containing information on the salt’s sharp structure and contrast magnitudes to define two independent Rayleigh Marchenko integral systems in the matrix-vector form (equation 13), one describing the relationship between focusing and pressure fields in the true medium and the other being a reference (e.g., migration) model,
Here, up- and downgoing pressure fields p± are given in terms of focusing functions f± through the real medium’s vertical particle velocity measured at the surface Vz±. Similarly, in the background model, V0,z± relates focusing functions f0± with pressure fields p0±. For the sake of notation, both systems in equations 15 can be written in a more compact form as p=Vf and p0=V0f0, respectively. Because we assume the background model to be known, quantities p0± and V0,z± can be estimated by forward modeling; as a result, one may have access to f0 by inversion of the reference system of equations 15 or 5 when the broadband fields g0± and R0 are modeled instead. Then, we introduce the perturbation quantities δp=pp0 and δf=ff0. They represent the scattered-wavefield effects in p and f as a result of the differences between the real and reference media. In light of this definition, any feature in the data that cannot be correctly explained by the migration model should appear in the perturbations. In terms of focusing functions and redatuming operators, the perturbed pressure fields arise after combining the Marchenko system in the actual medium, p=Vf, with the one in the reference, p0=V0f0. Therefore, an alternative Marchenko system for scattered fields reads as follows:
with δp=Vδf+δVf0. We refer to this new representation as scattering-based Rayleigh-Marchenko or SR-Marchenko. It describes the relationship between the unknown pressure field perturbations δp, in terms of f0 and those in δf via vertical particle velocity propagator V. A key element of our system in equation 16 is the perturbation operator δV=VV0, which captures the differences between data from the reference model and field data, i.e., the true medium. When it comes to applications to standard surveys, it is important that the data in the operators V and V0 are properly scaled relative to one another. We ensure operator compatibility by transplanting the measured downgoing vertical particle velocity vz+ into the reference medium at the receiver level Vz+V0,z+, and we used it as the aerial source. The seismic response of the reference model is then given by the one-way Rayleigh integral representation (equation 10) through a numerically modeled reflection response R0. Then, the following conditions relate the wavefield propagation in the true and reference medium: V0,z=R0Vz+, V0,z+=Vz+. The fields in the reflection response R0 do not consider multiples related to reflectors above the level D0, while introducing additional information from the background model into the upgoing vertical particle velocity v0,z. This argument in the SR-Marchenko system leads to a further simplification of the perturbation operator, as δVz+=0 and δVz=Vz+R0Vz+ because all medium perturbations occur below D0, so field perturbations only pertain to the upgoing fields.
Taking into consideration that the system of equation 16 is underdetermined, additional restrictions need to be imposed to reduce the number of unknowns. Once more, we invoke the same causality arguments introduced in the previous section and apply a muting window to reformulate equation 16 into a system that only relates focusing functions in the background, f0, with their perturbations counterparts, δf, i.e.,
which for simplicity, it can also be expressed in compact format as δVf0=Vδf, with V and δV encoding the influence of the windowing function Θ acting on V and δV. Unlike equation 4, equation 17 is a Fredholm integral of the first kind and does not admit a solution in terms of iterative substitution (van der Neut et al., 2015b). Finally, similar to equation 7, the aforementioned equation is solved as a two-step optimization problem:
minδfVδfδVf022s.t.  g0=M0f0.

From an operational point of view, numerical implementation is done first by constructing R0 and g0 using full-waveform forward modeling in the reference model, followed by inversion of the reference Marchenko system in equation 5 to find f0. Second, we solve the system of equation 17 for δf using the true-medium responses, Vz±, the perturbation operators δVz±, and f0 from the previous step. Finally, the perturbation in the pressure field δp is given according to equation 16, and the total redatumed field is then found as p=p0+δp.


We illustrate the practical implementation of our scattering-based redatuming schemes, S-Marchenko and SR-Marchenko, and compare their performance against conventional Marchenko and R-Marchenko, respectively. With the idea of evaluating the effectiveness of these methods, we designed a complex model describing different orders of scattering, which exhibits a high impedance contrast overburden mimicking the effects of an inhomogeneous (dirty) salt body (Figure 2). Beneath the complex overburden, this model is characterized by laterally heterogeneous sediments, interbed discontinuities, and sharp unconformities. The medium expands 16.26 km in the horizontal direction and extends up to 8.0 km in depth. In addition, the reference model that we use as a supplement in the following tests resembles the output of conventional high-end velocity model building, where the salt and sea bottom are represented by sharp discontinuities, whereas the sediment wavespeed is kept smooth and the salt wavespeed is constant. We first demonstrate wavefield reconstruction at virtual receivers located on a target level inside the medium, followed by interferometric redatuming. The former is carried through MDD, implemented to build local extended images (EIs). Finally, we migrate the EIs and discuss the main consequences of using accurate redatumed full-wavefields accounting for internal and free multiple reflections.

Wavefield extrapolation beneath complex overburden

The acquisition setup of our numerical experiments consists of 201 receivers regularly distributed every 40 m. To simulate the full reflection response R in the true medium, we use a finite-difference acoustic solver and sequentially inject 201 sources on the model surface. It is worth mentioning that the shot gathers are recorded for 8 s, allowing the observation of deeper events, and that their maximum aperture is 4 km. Similarly, a numerically modeled reflection response R0 in the migration-velocity model is implemented. For data modeling, a broadband impulse source with a flat spectrum in the range of 1–80 Hz is used. It should be noted that the relative source amplitude is preserved to ensure survey compatibility when constructing the perturbation operator δR. The extrapolation algorithm starts by modeling the full transmission response in the background model associated with a focal point at 4.4 km depth in the middle of the physical domain at 8.13 km using a 20 Hz Ricker wavelet (Figure 3d). For the original Marchenko redatuming scheme, the coda part of this transmission is typically discarded and only the direct part is preserved to find the initial focusing function, fd+, either by inversion of i=g0,dfd+, or approximated by time reversal fd+g0,d*. Contrary to that strategy, in S-Marchenko, we retain the full wavefield to solve for the true focusing function f0 (Figure 3a) in this model and we show that the coda part of g0 carries essential information for redatuming. Because no windowing constrains intervene in this step, internal multiples arriving after the direct wave contribute to the focusing process, as suggested by events appearing outside of the domain that would have been muted otherwise by a windowing operator. Now that f0 is available (Figure 3a), focusing function perturbations are estimated by our second inversion step. Figure 3b shows that f contains additional events explained by the reflection data in the true medium that add further information to the total focusing function, f in Figure 3c. As observed in Figure 3e, the scattered Green’s function δg retrieves reflection events coming from the deeper parts of the medium and internal multiples arriving after 3.5 s that do not appear in the reference Green’s function (g0 in Figure 3d). Comparing the information in the S-Marchenko-based focusing function (Figure 4a) to that in the standard approach (Figure 4b) exposes event mapping outside the range td(xr,xf)ttd(xr,xf) violating the causality argument invoked in defining the windowing operator used in Marchenko schemes. In practice, our S-Marchenko scheme can only update δf within the td(xr,xf)ttd(xr,xf) interval, but the δf update and the full-time range of f0 contribute to enhanced retrieval of δg — we provide further thoughts on this in the “Discussion” section. Furthermore, Figure 4c displays a substantial trace mismatch among focusing operators. At this point, we verify the retrieved Green’s function accuracy by comparing it with a directly modeled benchmark (Figure 5a). The recovered responses, corresponding to the same selected focal point, are depicted in Figure 5b for the S-Marchenko, and Figure 5c for the standard Marchenko. Although to some extent the fields retrieved by the conventional scheme in Figure 5c match those in Figure 5a, upon closer inspection, we observe substantial amplitude inaccuracies and many missing events. Likewise, examining the waveforms of the traces extracted from a receiver located at 8.6 km reveals a mismatch even for the early reflections. Such behavior is more pronounced for events arriving at later times (Figure 5d). However, our scattering-based equivalent (Figure 5b) achieves a much better fit especially for the earlier arrivals — not only the amplitudes but also the kinematics of the whole gather is substantially better retrieved. In this case, trace evaluation against the benchmark (Figure 5d) supports the hypothesis that including more accurate transmitted wavefields (implicit in our f0 fields), by means of a priori model information, can be critical for waveform fidelity in Marchenko-redatumed fields for media with this level of complexity.

Having exemplified the performance of Marchenko versus S-Marchenko, we move to our second experiment. In this case, the efficiency of the R-Marchenko and SR-Marchenko is subject to evaluation. Note that opposite to the previous case study, we anticipate that those Rayleigh-Marchenko versions could be more suitable candidates for real seismic surveys, considering that they require fewer preprocessing steps while still handling band-limited data containing internal and free-surface multiples (Ravasi, 2017). In addition, the R-Marchenko schemes can accommodate multiple acquisition geometries, as well as surveys with blended/simultaneous sources. To test these schemes, we extend our model with a 200 m water column to mimic an ocean-bottom survey geometry. On the source side, pressure and vertical particle velocity data for 201 sources located at a depth of 10 m are generated using a finite-difference vector-acoustic solver. Our synthetic ocean-bottom setup consists of 201 receivers regularly distributed every 40 m over a line at a depth of 200 m (Figure 2). Data are then decomposed into their up- and downgoing components (Wapenaar, 1998) followed by convolution with a 10 Hz Ricker wavelet to impose the band-limited character of the data to be used in the Rayleigh-Marchenko operator Vz±. Intentionally, this survey configuration allows us to use the full reflection response R0 from the previous experiment to generate the upgoing vertical particle velocity in the reference model V0,z through MDC with its downgoing counterpart assumed to be observed in the true medium, V0,z+. Because the medium is assumed to be unknown only below the receiver level (i.e., the sea bottom), using Vz+V0,z+ ensures proper operator relative scaling while simplifying the SR-Marchenko operator. In this paper, we purposefully include only results without free-surface effects in our Rayleigh-Marchenko analysis because those effects are far more dominant than those related to internal multiples and salt-body interactions, whereas the overall behavior of the R- and SR-Marchenko methods show analogous patterns as those seen in the examples below — they are representative of these two methods with or without free-surface effects. For a detailed comparison of R-Marchenko versus the original scheme in the presence of free-surface multiples, we refer the reader to Ravasi (2017).

For redatuming purposes, we select a focal point at position xf=(4.4,6.0) km and compute g0 similar to the previous experiment. Figure 6a6c depicts focusing functions from our band-limited ocean-bottom experiment. As in our previous study, the pressure response from the reference model (Figure 6d) offers access to crucial waveform information that is otherwise underestimated or neglected. As above, the scattered pressure wavefield (Figure 6e) reveals waveform features that play a major role in the substantial improvement of the total retrieved pressure field (Figure 6e). In the interest of assessing the relative accuracy in both redatuming scenarios, we analyze the full focusing functions (Figure 7) because they largely influence the redatuming performance. Whereas in the SR-Marchenko, a full-waveform transmitted field is encoded in the f0 used to retrieve f (Figure 7b), in the original R-Marchenko, an initial focusing estimate extracted from the directly transmitted first-arrival in the reference medium is used instead. When comparing the focusing response from the proposed scattering-based method (Figure 7a) with the standard Rayleigh-based focusing operator (Figure 7b), we observe significant waveform discrepancies between one another. First, the focusing function in the SR-Marchenko expands beyond the limits in the muting window due to the more complex f0, whereas that in the R-Marchenko lies within the domain imposed by causality. Second, events missing in Figure 7b appear in Figure 7a. Such arrivals are fundamental to ensure wavefield reconstruction and cancellation during the focusing process for retrieving p — ultimately resulting in better quality retrieved pressure fields. Keeping in mind the introduced windowing constraints in our current SR-Marchenko scheme, the information outside the windowing function is inherited from the solution of f0 in the reference model — but this information plays a role in the improved retrieval of events in δf. A final remark can be made when comparing the absolute and relative amplitudes of events retrieved by both methods (Figure 7c). There, we observe discrepancies between the two approaches that are noticeably significant — these differences are responsible for the absolute and offset-dependent amplitude differences between common events present in the p fields resulting from either approach.

Finally, we take a closer look at the recovered pressure fields and compare the performance of our SR-Marchenko approach relative to the R-Marchenko scheme. In Figure 8, pressure fields retrieved by the R- and SR-Marchenko schemes are compared with that obtained by forward finite-difference modeling. A general inspection of Figure 8b and 8c leads to the conclusion that SR-Marchenko consistently provides a robust estimation of the pressure fields considerably closer to the benchmark. Not only it is kinematically more accurate, but also the amplitude is better resolved on most of the retrieved events, contrary to those from R-Marchenko. Furthermore, a close examination of the retrieved fields reveals events missing in Figure 8c that are recovered in SR-Marchenko. Likewise, pronounced artifacts observed in Figure 8c appear suppressed in Figure 8b. Following from our analysis of the focusing fields above, we attribute such improvement to the relevant refinements added by the perturbed pressure wavefields to the background-related ones, which subsequently enhances the final pressure response. Extracting a seismic trace interval for detailed comparison (Figure 8d) indicates an overall better waveform fit by the SR-Marchenko in contrast with the considerable misalignment in the R-Marchenko reconstruction.

Target-oriented structural imaging

In the next examples, we design a virtual survey, using the different flavors of Marchenko-based extrapolated wavefields, to recreate the data that would have been recorded at a datum level by an actual seismic experiment, as if the sources and receivers were located just above the datum. Such data sets, in connection with traditional imaging conditions, are referred to as extended images (EIs) (Vasconcelos et al., 2010; Vasconcelos and Rickett, 2013; Ravasi et al., 2014). From that perspective, the media below the target level backscatters an upgoing pressure field p when illuminated from above with an incident downgoing pressure field P+. This interaction is described by the convolutional process p=P+R, where R is the local reflection response at the focusing depth and assumes a reflection-free medium above that level. An estimation of the local reflection response, with illumination from above the target, is obtained by implementing MDD (Wapenaar et al., 2008, 2011; Vasconcelos and Rickett, 2013); in Appendix  A, we describe our approach to preconditioned, time-domain MDD in more detail.

We start by reconstructing up- and downgoing pressure fields — by means of our Marchenko schemes — from the surface to an array of 151 virtual points spanning 3 km from 6.0 to 9.0 km horizontally, at a depth of 4.4 km. The virtual survey is shown in Figure 2 enclosed by the dashed target box beneath the salt body. Figure 9 shows the local reflection responses corresponding to a virtual source located at position xf=(4.4,7.5) km where different extrapolated wavefields are used in the deconvolution process. Figure 9b and 9c corresponds to the virtual shot gathers resulting from using broadband S-Marchenko and Marchenko-based redatumed fields, respectively. We identify multiple coherent reflection events when evaluating against a synthetic gather modeled by finite differences, using only the target-medium properties (Figure 9a). Although here we display only one of the 151 common-source EIs retrieved by MDD, the behavior displayed in these gathers is representative of the overall MDD results. In spite of the implicit numerical instabilities introduced by the MDD problem, a considerably more consistent reconstruction is achieved when relying on scattering-based Marchenko fields contrary to those from the original scheme. Although some events are underestimated in the conventional method, Figure 9c, a stronger contribution for some of the events arriving after 1.5 s is appreciated in Figure 9b. Similarly, we present EIs obtained from band-limited redatumed data (Figure 9e and 9f). Contrary to the previous case, redatuming is done with the R-Marchenko and SR-Marchenko methods as input for the MDD. Figure 9e shows the result for SR-Marchenko, whereas Figure 9f shows the result for R-Marchenko. In these EI gathers, given the relatively simple geology beneath the salt, we do not expect to have strong interfering or discontinuous events. However, as observed in Figure 9e, the scattering-based method presents a more stable result than that of the conventional R-Marchenko when compared to the benchmark (Figure 9d). The prominent event discontinuities seen in the R-Marchenko and original-scheme EIs result from poorly reconstructed events in their corresponding downgoing fields — these do not capture the correct waveform behavior as a function of the subsurface offset in subsalt areas, leading to illumination gaps in the retrieved EIs. However, S- and SR-Marchenko downgoing fields contain more arrivals with improved waveform fidelity, resulting in EI gathers with better subsurface-offset amplitude compensation that are considerably closer to the benchmark responses.

To assess how the differences in EI gathers may translate into imaging, we conduct target-oriented RTM on the virtual shot gathers predicted by the different extrapolation methods. We compute 151 single-shot images and then stack them together to find the target-medium (box in Figure 2) migrated images (Figure 10a10c). This is in line with targeted Marchenko imaging approaches presented in previous studies (e.g., Wapenaar et al., 2014b; Ravasi, 2017). In this paper, for the sake of brevity, we refrain from showing comparisons between the Marchenko-based target-oriented images and the surface-data-based migrated images of the target area. Overall, all of our Marchenko-based target images are superior to their surface-based RTM counterparts, with improvements comparable to those seen in previous studies (e.g., Wapenaar et al., 2014b; Ravasi, 2017), and with our S- and SR-Marchenko-based images being considerably better than those from surface-data-based RTM.

Here, we focus our analysis on comparing the imaging performance of our approaches versus existing Marchenko schemes, strictly in the context of target-oriented imaging. As such, Figure 10c corresponds to the conventional Marchenko local reflection response (Figure 9c), Figure 10b corresponds to the S-Marchenko result (Figure 9b), and Figure 10a corresponds to the RTM of the target-medium benchmark (Figure 9a). Our results show that the scattering-based image (Figure 10b) yields better fidelity to the benchmark in contrast to the traditional Marchenko-based image. In this case, the gain in resolution due to the better handling of complex salt-influenced wavefield is noticeable, especially for the large offsets as observed toward the image edges. In addition, the details of stratigraphic features are also better resolved. As for the Marchenko-based image, the overall incorrectly extrapolated up- and downgoing wavefields may result in structures that appear to be shallower than the true interfaces. This likely results from the improper handling of the direct-arrival phase distortions due to the salt body interactions. Regarding band-limited structural images corresponding to the Rayleigh-based Marchenko fields (Figure 10a10c), a reliable reconstruction of the main reflectors can be observed in the SR-Marchenko image (Figure 10b). Unlike that of R-Marchenko (Figure 10c), in which additional artifacts are spotted for the deeper areas at approximately 7.0 km, the SR-based image exhibits a closer representation of the subsurface as in the reconstruction in Figure 10a. Moreover, substantial differences in image quality are visible in the near surface (approximately 5.4 km) as a consequence of an erroneous phase and amplitude estimation of the early arrivals in the redatumed wavefields (Figure 8d).

For a more detailed analysis, we select multiple close-ups (indicated by the boxes in Figures 10 and 11) at different locations (Figures 12 and 13). For the first row in Figure 12 matching the green box in Figure 10, the most prominent differences between the S-Marchenko and Marchenko RTM images are indeed in the shallow section of the target medium; still, noticeable differences between the two images also exist at deeper locations. In box 1 (the top row), we see that the SR-Marchenko image yields a consistently more reliable representation of the complex seismic facies present in the benchmark image, which could impact potential interpretations in an exploration or reservoir characterization context. The next row corresponds to the magnified sections associated with box 2 across the multiple images. Even though most reflectors are present in both cases, a better continuity in the S-Marchenko image (Figure 12b) exposes a remarkable enhancement. We note that this is not the case for the reflectors reconstructed by the Marchenko-based image, which in fact, displays a prominent image artifact in the left-center portion of Figure 12c. An additional observation concerns image reconstruction at far offsets: The last row in Figure 12, box 3 in Figure 10, displays a much better resolved reflector below 7.0 km for the S-Marchenko case, not only due to the accurate prediction of discontinuities but also because it can erroneously be interpreted as a mild dipping layer in the traditional Marchenko image. Figure 13 is the band-limited counterpart of Figure 12 for the same magnified areas using the SR-Marchenko (Figure 13b) and R-Marchenko (Figure 13c) schemes contrasted against a 20 Hz RTM benchmark (Figure 13a). In general terms, similar improvements to the ones described for the S-Marchenko case are observed. Based on the previous examples, we expect that in the case of data containing SRMs, SR- and R-Marchenko imaging would result in images like those in Figures 12 and 13, whereas conventional Marchenko imaging becomes significantly degraded by the presence of free-surface effects (Ravasi, 2017).


In this paper, we present two novel redatuming scattering-based schemes for wavefield redatuming with the Marchenko framework whose objective is to account for highly complex overburden media such as subsalt or subbasalt. With numerical experiments in a complex salt model, not only do we expose issues with the previous schemes, but we also convincingly show that our scattering approaches result in fields and images with considerably superior waveform fidelity. At the point of writing this study, these are the best scattering-based estimates we have to offer; however, we point out that they are approximate and likely could further improve. One key reason for this is our use of the windowing operator in the scattering schemes. As it stands, the windowing is a necessary constraint that allows for the elimination of Green’s function perturbations from the Marchenko representations, e.g., as shown in equations 6 and 16, leading to the focusing-only systems that can numerically be solved. In practice, however, we indeed observe that the reference focusing functions f0 contain events not only after the “direct arrival” — as is assumed by conventional Marchenko schemes — but also before it, as shown in Figures 3 and 6. This, in fact, is physically expected because in such complex media, waveform triplications/caustics are common, so not only are there multiple events corresponding to the “direct” wave but also other waves may arrive before some of the later paths of the ballistic/direct field. Here, we note that previous iterations of the Marchenko scheme (Wapenaar et al., 2014a) are shown to handle triplications in the retrieved coda, whereas those in the direct wave need to be properly included in the estimated initial direct wave. In our case, (1) caustics are handled on the direct/ballistic arrivals, (2) including multiple arrivals of the direct field in f0, with (3) the correct dynamic/amplitude behavior properly accounted for because f0 is obtained by inversion. As a consequence of this, the behavior we see in f0 should, in principle, also apply to f and thus to δf. However, because we rely on the same windowing operator — and thus the same causality arguments — of conventional Marchenko schemes, our practical approach only retrieves updates to δf that lie within the causality window. Any waves contributing to δf that would lie outside the causality window are currently not retrieved by our practical schemes. In the tests that we conducted so far, even without estimating δf contributions lying outside the causality window, we obtain redatuming results with considerably better accuracy than those produced by the conventional approach. But that is not to say that these missing contributions do not matter because there are still inaccuracies and missing events in our reconstructions (e.g., as shown in Figures 5 and 8). Building a better understanding of the role of focusing function wave components outside of the causality window, together with developing approaches that can estimate them, is the subject of ongoing research.

Continuing with the discussion on causality arguments and windowing, we point out that there is also an important distinction between the reflection operator and Rayleigh-based Marchenko constructs. When defining causality constraints for the representations such as in equations 3 and 6, the windowing operator has two purposes: The first purpose is to annihilate g from the system, and the second is to isolate fm+ and f. Mathematically, this is so because the representations have a form, e.g., gf=Rf+, in which the left side contains a superposition of the Green’s and focusing functions. This is not the case for Rayleigh-Marchenko representations. There, the representations have a form, e.g., p=Vzf++Vz+f, in which there is no superposition of pressure fields and focusing functions on the left side. As a result, the only role of causality constraints in Rayleigh-Marchenko is to eliminate the p wavefields from the system. This seemingly innocent difference between the two types of representations implies that windowing operators can in fact be different in practical implementations of either scheme. For example, this in principle allows for greater flexibility in the context of R-Marchenko windowing where different, time-asymmetric operators could be applied separately for up- and downgoing pressure field components. We have begun to look into the implications of this observation, and we will report on the outcome in a future publication.

Given that the scattering schemes we present here make explicit use of a priori models with sharp contrasts included, it is natural to question the sensitivity of our approaches to errors in these models. First, it is important to note that, even for the conventional Marchenko redatuming approach, in the case of complex media (e.g., subsalt) one requires the direct-arrival estimate to be as kinematically accurate as possible. In previous Marchenko studies in subsalt environments — which included field data (Jia et al., 2018; Staring et al., 2018; Staring and Wapenaar, 2020) — complex salt models were in fact used for the estimation of direct arrivals for the initial focusing functions. This is done because, in the presence of large contrasts such as salt, model smoothing can lead to significant phase distortions that affect the direct-wave kinematics. This is not only true for Marchenko redatuming, but it is generally a well-known fact in imaging practice, e.g., in the case of RTM in subsalt environments. Therefore, in that sense, here our approaches make better use of the same type of high-end models already in use for the purpose of general depth imaging as well as Marchenko redatuming. However, model errors do have an effect in Marchenko redatuming. In the context of the original Marchenko framework, previous studies (Thorbecke et al., 2013; Broggini et al., 2014) have shown that velocity model errors lead to image distortions similar to those in conventional migration, but that the handling of internal multiples and their artifacts is rather robust with respect to relatively large velocity errors. In our case, because the benefits arising from the scattering schemes are rooted on complex wave events in f0 — and those depend on the sharp impedance contrast in the model — we expect our methods to have greater sensitivity to model parameters than the conventional Marchenko, and in particular to the geometry of sharp interfaces. In general, however, we expect that current industry-standard models built to maximize migration image quality — i.e., to optimize focusing in subsurface angle gathers — in complex environments to be sufficiently accurate for our redatuming purposes. Verifying this claim, as well as building a more thorough understanding of the role of model accuracy in our approaches is beyond the scope of this paper, but it is the subject of current study.

One of our main goals with this work, particularly through the scattering-Rayleigh scheme, is to enable Marchenko-driven redatuming and imaging applications in complex media that also handle band-limited data with realistic acquisition configurations. In a companion paper, Ravasi and Vasconcelos (2021) present an open-source-based implementation of Marchenko operators for high-performance computing (HPC) environments — cluster or cloud-based — and demonstrate how such a framework enables research and development of approaches such as ours for 3D seismic data. The Pylops-based (Ravasi and Vasconcelos, 2020) codes that we developed and used in this study are already deployable under the HPC framework by Ravasi and Vasconcelos (2021); as such, our approach is ready for validation with 3D field seismic data. At the time of writing this paper, we have recently been given access to a 3D ocean-bottom, permanent-reservoir-monitoring (PRM) data set from the North Sea, which we will use as a testbed for validation and further research of the approaches herein.

As a final discussion point, we highlight that the improvements in waveform fidelity of redatumed fields in complex media brought by our scattering-based schemes bring new opportunities in the broader context of target-oriented inversion and monitoring. In the examples that we provide here, Marchenko-redatumed fields are used for estimation of local-reflection EIs, and they are subsequently used for local imaging by means of depth migration — this is inline with, e.g., the Marchenko imaging steps as described by Wapenaar et al. (2014b). Although this is one option for targeted imaging using Marchenko fields, it is by no means the only option. For example, when retrieved with sufficient waveform fidelity, Marchenko fields can be used in local full-waveform inversion approaches (Vasconcelos et al., 2017; Elison et al., 2018; Cui et al., 2020). Alternatively, Marchenko fields may be integrated in other inversion schemes, such as in an inverse-scattering-series-based scheme (e.g., Guo and Alkhalifah, 2020) or directly into wave-equation solutions as a constraint (van Leeuwen and Herrmann, 2015; Diekmann and Vasconcelos, 2021). These being some of the uses of Marchenko fields in inverse schemes, they attest to the potential of such fields in contributing to target-oriented inversion.

Finally, although in this paper we focus on targeted imaging in complex media, the scattering formalism and approaches that we present are equally suitable for time-lapse problems — in which the reference medium and data are no longer dictated by a background model, but are instead given in the form of the baseline model and observed data, respectively. The application of our scattering schemes to monitoring will be a focus in our upcoming field data studies in the context of PRM.


Although the recent rise of redatuming and imaging approaches based on the Marchenko framework showed great promise for tackling the challenges of imaging in highly complex media, such as subsalt or subbasalt, previous Marchenko schemes impose limitations that can severely hinder accurate wavefield retrieval in such media. We present scattering-based approaches to Marchenko redatuming that can handle highly complex media, by taking advantage of high-end subsurface models typically built for migration. Our approaches take wavefield redatuming accuracy to a level of waveform fidelity in redatuming far beyond that of conventional imaging, resulting in EI gathers and target-oriented images that are consequently highly accurate — thus realizing the potential of the Marchenko framework for imaging of complex subsurface environments. To facilitate the deployment of our scheme to band-limited data, with realistic acquisition and free-surface conditions, we combine our scattering approach with the Rayleigh-Marchenko scheme into a new method that can produce accurate wavefields and images from offshore or ocean-bottom 4C data sets. We support our new approaches with numerical experiments using a highly complex, dirty-salt 2D model. Given our robust theoretical and numerical constructs, we expect that our approaches will be useful for imaging and monitoring in challenging subsurface geologic conditions.


We are grateful to H. Peng, L. Diekmann, A. Tataris, and T. van Leeuwen for insightful discussions that contributed to the work in this paper. I. Vasconcelos and M. Ravasi are thankful to J. van der Neut, whose input contributed to the early development of our ideas on scattering-based Marchenko. We thank the sponsors of the Utrecht Consortium for Subsurface Imaging for their financial support.



EIs (Vasconcelos et al., 2010) R(xf,xf;ω) are virtual reflectivity gathers at a target level Di formed by the interaction of decomposed wavefields P+(xs,xf;ω) and P(xs,xf;ω), connecting the focal point xf to a set of sources (or receivers) xs at the surface level. In the frequency domain, the MDC process can be discretized as follows (Wapenaar et al., 2011; van der Neut and Herrmann, 2012):
The associated normal equations are then written as C=ΓR (e.g., Vasconcelos and Rickett, 2013), where the point spread function Γ is a blurring operator acting on R determining the crosscorrelation gather C,
Here, it is important to note the use of the convention that discrete wavefield operators have sources on the row space and receivers on columns — this is opposite to the commonly used convention in MDD (e.g., Vasconcelos et al., 2010). The reason for this choice is that our numerical iterative solvers (e.g., LSQR), under the Pylops framework (Ravasi and Vasconcelos, 2020, 2021), are designed to solve Ax=y for x, assuming A is a left-hand operator; numerically this means that we choose our row/column space convention so that the forward-operator fields are left-hand operators.
A common approach to deconvolution in the frequency domain is to compute an inverse matrix for each frequency, W(ω)=(Γ(ω)+ϵI)1, acting on the crosscorrelation gather, where ϵ aims at stabilizing the inversion and can be interpreted as a Tikhonov regularization term. In contrast, a time-domain representation of the very same problem requires the definition of an operator P^+=F1P+F, which acts on a vector and performs a step of MDC as described by Ravasi and Vasconcelos (2020). Here, F and F1 are the forward and inverse space-time (2D) Fourier transform operators, respectively. Although an explicit definition in computer memory of the equivalent matrix of such an operator is prohibitive, many iterative solvers only require the resulting forward and adjoint operations of the convolutional matrix operator on a given data vector — not an explicitly defined version of the operator in matrix form. With that in mind, using object-oriented representations of linear operators (Ravasi and Vasconcelos, 2020) is an attractive alternative to solve large-scale optimization problems in the time domain. Our implementation of MDD in the time domain allows for the introduction of two preconditioners and can be written in the following way:
where W is a muting operator imposing frequency-wavenumber constraints and Θ is a muting operator in the time-space domain enforcing causality on the reflection response.

Because C is a blurred version of R, it encompasses most of its data-space features in terms of the frequency-wavenumber support that can realistically be retrieved by MDD on a shot-by-shot basis. Consequently, the fk-domain spectrum of the crosscorrelation gathers are used to design a natural preconditioner W constraining the inversion within the expected frequency-wavenumber support of C. In the time domain, Θ has two roles: (1) It ensures the reflectivity responses from MDD are strictly nonnegative, and (2) it enforces wave-equation causality by annihilating contributions arriving before direct waves (based on the a priori migration model). Our experience shows that such preconditioning is essential to the retrieval of numerically robust, physically reliable, EI/virtual source gathers by time-domain MDD.

Biographies and photographs of the authors are not available.

Freely available online through the SEG open-access option.