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.
BROADBAND WAVEFIELD EXTRAPOLATION
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
Evaluation of equation 4 requires the direct part of the downgoing focusing function () to be known. In principle, if the direct part of the transmission response in the truncated reference medium is available, can be obtained from the focusing condition, , by direct inversion (van der Neut et al., 2015b). Here, vector represents a band-limited delta function at the focal point 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., . 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.
BAND-LIMITED WAVEFIELD EXTRAPOLATION
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).
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 , 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 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
From an operational point of view, numerical implementation is done first by constructing and using full-waveform forward modeling in the reference model, followed by inversion of the reference Marchenko system in equation 5 to find . Second, we solve the system of equation 17 for using the true-medium responses, , the perturbation operators , and from the previous step. Finally, the perturbation in the pressure field is given according to equation 16, and the total redatumed field is then found as .
SUBSALT SYNTHETIC EXAMPLES
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 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 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 . 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, , either by inversion of , or approximated by time reversal . Contrary to that strategy, in S-Marchenko, we retain the full wavefield to solve for the true focusing function (Figure 3a) in this model and we show that the coda part of 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 is available (Figure 3a), focusing function perturbations are estimated by our second inversion step. Figure 3b shows that contains additional events explained by the reflection data in the true medium that add further information to the total focusing function, in Figure 3c. As observed in Figure 3e, the scattered Green’s function 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 ( 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 violating the causality argument invoked in defining the windowing operator used in Marchenko schemes. In practice, our S-Marchenko scheme can only update within the interval, but the update and the full-time range of contribute to enhanced retrieval of — 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 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 . Intentionally, this survey configuration allows us to use the full reflection response from the previous experiment to generate the upgoing vertical particle velocity in the reference model through MDC with its downgoing counterpart assumed to be observed in the true medium, . Because the medium is assumed to be unknown only below the receiver level (i.e., the sea bottom), using 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 km and compute similar to the previous experiment. Figure 6a–6c 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 used to retrieve (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 , 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 — 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 in the reference model — but this information plays a role in the improved retrieval of events in . 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 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 when illuminated from above with an incident downgoing pressure field . This interaction is described by the convolutional process , where 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 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 10a–10c). 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 10a–10c), 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 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 , with (3) the correct dynamic/amplitude behavior properly accounted for because is obtained by inversion. As a consequence of this, the behavior we see in should, in principle, also apply to and thus to . 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 that lie within the causality window. Any waves contributing to 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 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 from the system, and the second is to isolate and . Mathematically, this is so because the representations have a form, e.g., , 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., , 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 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 — 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.
LOCAL EIs BY TIME-DOMAIN INTERFEROMETRIC REDATUMING
Because is a blurred version of , 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 -domain spectrum of the crosscorrelation gathers are used to design a natural preconditioner constraining the inversion within the expected frequency-wavenumber support of . In the time domain, has two roles: It ensures the reflectivity responses from MDD are strictly nonnegative, and 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.