An analysis of acquisition-related subsampling effects on Marchenko focusing, redatuming, and primary estimation

By means of Marchenko methods, the full-waveform Green’s function to a virtual source in the subsurface can be retrieved from single-sided reflection data. Given a smooth velocity model (Wapenaar et al., 2014), this is accomplished by solving a set of coupled Marchenko equations by direct least-squares (LSQR) inversion (van der Neut et al., 2015a) or iterative substitution (van der Neut et al., 2015b). The method, in theory, retrieves internal multiples with correct amplitudes. One of the challenges of this method is that the multidimensional convolution operation of the reflection response with the focusing function requires a regularly sampled, dense data set with colocated sources and receivers. Similar multidimensional convolution operations are also used — and fairly well-understood — in surface-related multiple elimination (SRME) (Verschuur et al., 1992; Lopez and Verschuur, 2015). Here, we look into the effects of regular subsampling in the Marchenko method specifically — showing similarities and differences to the effects seen in SRME applications. Haindl et al. (2018) formulate the Marchenko equations with sparse inversion to alleviate the 2D irregular source subsampling effect. This study observed spatial gaps in the output redatumed Green’s function and the focusing function. Ravasi (2017) also studies the effect of irregular subsampling in sources in the Rayleigh-Marchenko method with LSQR inversion and observed artifacts, but no spatial gaps in the retrieved signals, thereby showing the Rayleigh approach’s flexibility in handling the source subsampling. Sripanich and Vasconcelos (2019) provide a comprehensive study on the effects of a limited aperture on Marchenko focusing functions and their radiation behavior of focusing. In synthetic and field data studies, respectively, Jia et al. (2019) and Staring and Wapenaar (2020) present 3D Marchenko redatuming, in which they apply data interpolation to fill the gaps caused by subsampling. More recently, in a companion study to this paper, regular and irregular subsampling effects on 3D redatuming are discussed by Ravasi and Vasconcelos (2020), in the context of high-performance computing. To account for the imperfect source sampling, Wapenaar and van IJsseldijk (2020) develop a new discrete representation containing point-spread functions that can be multidimensionally deconvolved to deblur the focusing functions and Green’s functions. Van IJsseldijk and Wapenaar (2020) then alter the iterative Marchenko scheme to integrate this new representation and test it on synthetic data. All of the preceding studies alleviate the subsampling effects using different approaches, including sparse inversion, LSQR inversion, and/or data interpolation. Yet, there is still a need to build a basic understanding of these effects’ origin and patterns, which could not only inform acquisition sampling choices in the future but also crucially help in the design and quality control of preprocessing routines tailored toward Marchenko methods.

Although the original form of the Marchenko system is designed to perform redatuming, i.e., retrieve subsurface wavefields, the framework can also be used for surface-to-surface internal multiple prediction — or more specifically — to predict primary-only data. To that end, an alternative, surface-projected form of the Marchenko equations presented by van der Neut and Wapenaar (2016) was recently used to eliminate multiples directly from reflection data without requiring model information (Zhang and Staring, 2018). In a comparative study, Zhang et al. (2019a) illustrate the benefits and shortcomings of Marchenko-based data-driven scheme for the internal multiple reflection elimination over the theory of inverse scattering series. The former method was further improved to compensate the transmission loss (Zhang et al., 2019b). In this context, i.e., Marchenko-based primary estimation (MPE), the effects of acquisition sampling remain largely unaddressed. Very recently, this projected scheme was also used in the augmented Marchenko method to account for short-period internal multiples in media with a horizontally layered overburden by Elison et al. (2020).

In recent years, Marchenko-based methods have been developed in two related, yet distinct, contexts: depth-domain redatuming and data-domain primary estimation. In this paper, we choose to look into Marchenko approaches in both of these contexts because they typically pertain to different applications — namely, redatuming is tied to full-wavefield imaging at depth, whereas primary estimation is usually a preprocessing step prior to model building and imaging. To this end, building on the aforementioned studies that look into acquisition aspects of the Marchenko approaches in specific contexts, we intend to provide a basic understanding of acquisition sampling effects in focusing, redatuming, and primary estimation.

In addition, two types of Marchenko methods are addressed in this paper because they require different inputs and temporal-windowing approaches. Albeit relying on the same general operations of multidimensional convolution and correlation, these two types of Marchenko approaches can behave differently in the presence of subsampling. There are two main reasons for this: (1) Marchenko redatuming does not perform time stepping and uses the full wavefield that has been retrieved by focusing for a given redatuming level, and primary estimation applies time stepping and just selects the samples at that time after applying the Marchenko method at a given time step and (2) the windowing operators themselves in Marchenko redatuming and primary estimation are different, leading to differences in the physical contributions of the convolution processes in either of the applications. We observe that subsampled input data cause artifacts for both methods and may further lead to divergence for the primary estimation, depending on the degree of subsampling and the algorithm used to solve the system, i.e., Neumann series expansion or LSQR inversion. Because we wish to analyze subsampling at varying levels of aliasing, we maintain a fixed bandwidth and only vary the sampling. In practice, sampling and aliasing are relative to the actual data bandwidth and, as a general rule, adequately sampling the integrals requires satisfying the general spatial Nyquist criteria for the integrands — similarly to well-known applications such as SRME or evaluating Kirchhoff-type integrals.

Our study of the subsampling effects is structured in a simple, yet general, framework that could inform acquisition design or help in tuning and controlling data processing. With this goal in mind, we limit ourselves to 2D regular subsampling in just one dimension (source or receiver). We believe that the insight of this case shares similarities with the other cases, e.g., irregular subsampling in both dimensions. Because three dimensions and data interpolation are related, yet complex, topics, they are left for future research and therefore are not included in this study. First, we look at subsampling numerical examples in two dimensions to mimic subsampling caused by realistic acquisition design. Then, we analyze the results of the aforementioned two Marchenko methods in the context of subsampling without data interpolation, also by looking at specific yet simple measures to alleviate the effects of subsampling, such as invoking spatial reciprocity.

In all of the current forms of Marchenko methods — whether in the context of redatuming or primary estimation — the theory clearly outlines that monopole (e.g., pressure-to-pressure) and dipole (e.g., particle velocity-to-pressure) fields are not only required but must be used in specific source and receiver configurations (Wapenaar et al., 2014; van der Neut et al., 2015b). When it comes to field data, e.g., in offshore acquisition settings, although monopole fields are easily available from streamer and ocean-bottom acquisitions, dipole fields are usually only available on the receiver side in the form of acceleration or velocity measurements from multicomponent sensors. Thus, in this paper, we show the effects of subsampling in the presence of dipole sources — which are generally not available in practice — as well as dipole receivers. We do this to provide insights on the importance and effects of dipole integration choices that may impact choices on how to use multicomponent sensor data for Marchenko applications or how to preprocess 1C data to be used in Marchenko schemes.

The layout of this paper is as follows. We begin by presenting the Marchenko integral representations, while paying particular attention to the discrete form of the terms within the iterative solution of the coupled Marchenko system, followed by a brief review of the primary estimation method (Zhang and Staring, 2018). With a layered benchmark model, we then illustrate our observations on the subsampling effects on Marchenko focusing for multiple scenarios, including the joint effects of subsampling versus the dipole source/receiver integration variable choice and the effects of focusing depth jointly with subsampling. Then, a complex subsalt model is used to illustrate subsampling effects together with limited aperture on Marchenko redatuming in a more realistic scenario. Subsequently, we show with numerical examples for both Marchenko systems that integrating on the monopole dimension instead of the (theoretically correct) dipole dimension actually yields negligible errors compared with the true results for media with mild lateral heterogeneity. Finally, for the primary estimation, we compare the joint effects of different subsampling schemes versus different approaches to solve the Marchenko system with a layered model and provide an operator norm analysis for the iterative substitution method.

One key point to note, here as well as throughout this paper, is that we refer to monopole and dipole sources and receivers so as to relate the numerical results to the fields in the preceding equations. However, in real life, practical acquisition, dipole fields are not necessarily always available. For example, in an offshore acquisition setting, sources are generally of the monopole type only — this means that to obtain dipole sources for calculations one would need to apply, e.g., local obliquity factor corrections. On the receiver side, however, dipole observations correspond to either direct particle velocity or acceleration measurements, which are generally present in ocean-bottom systems and available in some multicomponent streamer systems.

It is important to note here that although redatuming requires the numerical solutions only once per redatuming location, MPE requires equation 18 to be solved for every time sample of the preceding iteration — for primary prediction for a single chosen shot gather. Hereafter, we rely on this algorithm to analyze the subsampling effects in MPE.

For the purpose of benchmarking, we use the same constant-wavespeed (2000 m/s), variable-density layered model (Figure ²) from van der Neut et al. (2015b), who use it to provide a comprehensive study of the integral contributions to the Marchenko scheme. We choose this model such that the reader can refer to individual events and understand how they are generated — van der Neut et al. (2015b) provide a detailed analysis of how each event is reconstructed by redatuming in relation to the events present on the surface data, via an interferometric interpretation of the Marchenko integrals acting in the practical iterative scheme. Back to our case, Figure ² and ² shows the regular sampling masks for the source and receiver dimensions, which are used in all examples of Marchenko focusing and redatuming in this study. The full reflection data have 201 colocated sources and receivers with a spacing interval of $10m$⁠. The acquisition spread lies within the range of $x=(−1000m,1000m)$ at the acquisition surface with the time sampling rate $dt=2.5ms$⁠. Two focal points are chosen for this model: one at $(XF=0m,ZF=2670m)$ and the other $(XF=0m,ZF=1125m)$⁠, as shown in Figure ². Our study in this section will focus on the deeper focal point as discussed by van der Neut et al. (2015b), with the shallower one used later as a comparison. In our examples, all reflection responses are modeled by a spectral-element (SEM) wave solver from the Salvus package (Afanasiev et al., 2019) with monopole sources and dipole receivers.

Because of the symmetry of the chosen model, the reflection responses are the same for colocated shot and receiver gathers, which allows us to freely integrate over either dimension. Figure ³ shows the focusing functions in the time domain obtained by integrating on different chosen variables, with source- or receiver-domain subsampled data. Here, we show results after a fixed number of iterations (13 in this case) of the Marchenko scheme — the number of iterations is chosen arbitrarily to be more than sufficient for the scheme to achieve convergence, which in ideal sampling conditions occurs in a mere few (e.g., 3–5) iterations. Figure ³ shows the benchmark focusing function obtained without subsampling. Figure ³ displays the same focusing function obtained by subsampling in either the source or receiver and integrating on the same corresponding variable. Because of the subsampling, aliasing artifacts appear at far offsets but both focusing solutions in Figure ³ are fully sampled in space, i.e., the surface-coordinate dimensions. In comparison, Figure ³ displays clear spatial gaps obtained by subsampling in the source or receiver but integrating on the other variable, together with the far-offset artifacts shown in Figure ³. A comparison of traces at far offset $xr=910m$ in the three focusing functions of Figure ³–³ displayed in Figure ³ shows that artifacts appear in both integration scenarios. From the trace plot, we can also see that the artifacts in Figure ³ seem more unstable than those in Figure ³ due to poor sampling near the stationary-point contributions — these stationary-point regions are discussed by van der Neut et al. (2015b). Figure ³–³ shows the magnified plots of Figure ³–³ around the circles, respectively.

When integration is conducted over the subsampled dimensions (Figure ³), the summation process in every element of the output upgoing focusing function $f^1−$ in equations 5 and 9 becomes increasingly inaccurate as a result of the increased subsampling interval, but they remain well sampled in space because the row space of equations 5 and 9 is well sampled. This property carries over in the subsequent convolution operations needed in the iterative scheme or in any other iterative linear solver (e.g., LSQR), leading to error propagation at later iterations. As a consequence of the negligible scaling in the source-receiver reciprocity in terms of the horizontally independent obliquity factor present in the monopole-dipole configuration (in this case because the medium is laterally homogeneous), subsampling and integrating on the source or receiver dimension are equivalent. When the subsampling and integration are performed over different dimensions, the summation process in every nonzero element of the output upgoing focusing function $f^1−$ in equations 5 and 9 will be accurate because the column space of equations 5 and 9 is well sampled. However, subsampling will result in spatial gaps in the output focusing function because equations 5 and 9 are now subsampled in the row space. At later iterations, the row-space gaps will remain and inaccuracies arise in the nonzero elements of the output focusing functions because the input focusing functions gaps carry over into the summation process, as shown in Figure ³. A diagram of all four cases in integration versus subsampling is shown in Figure ⁴.

The corresponding focusing functions in Figure ³–³ are transformed to the frequency-wavenumber (⁠$f−k$⁠) domain, as shown in Figure ⁵, ⁵, and ⁵, separately. Normalized by the maximum value of Figure ⁵, ⁵, and ⁵ shows the differences between the benchmark response (Figure ⁵) and those in Figure ⁵ and ⁵, respectively. By comparing all of the panels in Figure ⁵, we see that the aliasing effects caused by subsampling and integrating on the same dimension are relatively weaker — compared with subsampling and integrating on different dimensions in which artifacts are generally stronger.

Figure ⁶ illustrates how the spatially dependent artifacts at the far offset come into being in Figure ³, as a result of the column-space subsampling in equation 5. Figure ⁶, ⁶, and ⁶ shows the convolution gathers (i.e., the integrands) of $Rf1d+$ at different receiver positions, where the gaps represent the missing sources. Their summation (properly scaled by the subsampling interval), corresponding to different elements in the vector of the output $f^1−(xr;xF)$⁠, is shown by the blue traces in Figure ⁶, ⁶, and ⁶, respectively. For reference, we overlay the corresponding traces of the focusing function without source subsampling. As can be seen from Figure ⁶, the focusing function traces obtained by the first convolutional operation with subsampling are very close to the reference ones at the near offset because the traveltime varies slowly with respect to $xs$⁠. Toward larger offsets, the rate of traveltime variations in the integrand increases, leading to steeper slopes in the integrand and thus aliasing-induced artifacts.

Figure ⁷–⁷ shows the errors caused by subsampling and integrating on the same dimension from the first to the fourth iterations. As can be seen, most of the error comes from the second iteration and then the error decreases in magnitude with increasing iteration numbers — this is expected assuming that the iterative scheme is applied within its convergence criteria of the spectral radius (the maximum magnitude eigenvalue) of the operator, and as such our observations should hold for any media that satisfy them (Dukalski and de Vos, 2018). This is the main reason that we focus our analysis here on the artifacts arising at the leading iteration of the Marchenko redatuming scheme.

Because most data sets have fixed source and receiver apertures, the behavior of focusing functions and resulting redatumed fields changes with the focal-point depth. Hence, we complement our observations by studying the effects of focusing depth together with subsampling. Figure ⁸ and ⁸ shows the focusing functions calculated by integrating on source without and with subsampling at a shallower focusing depth $(XF=0m,ZF=1125m)$⁠. Comparing Figures ³ and ⁸, we see that the artifacts at the far offset brought by the same subsampling appear to be stronger for the shallow focusing depth than those for the deep focusing depth. This effect can be explained by comparing the traveltime slopes of the two input focusing functions as shown in Figure ⁸ and ⁸ because the shallower depth focusing function in Figure ⁸ has steeper slopes at far offsets, the integrand events at far offsets will be more aliased than those from a deeper focal point (Figure ⁶), given a fixed subsampling rate.

The Marchenko methods have a requirement of monopoles and dipoles, either on the source or receiver side, as we show in equation 3. Because of the symmetrical geometry in previous examples, the source-receiver transpose yields the same reflection record; therefore, integrating on the dipole or the monopole dimension yields the same result in the Marchenko method for laterally homogeneous media, as shown in Figure ³. However, this may no longer be true for laterally inhomogeneous media because, e.g., the angles of rays coming in and out at a fixed source-receiver pair may not be equal anymore. In practice, real dipoles are typically only available on the receiver side, for example, through particle velocity or acceleration measurements. On the source side, they are generally unavailable on offshore experiments because the source is monopole in its nature. The common practice in industry when one needs a dipole on the source side is to apply model-based corrections that would require knowledge of the local wavenumber and impedance at the source locations. Here, we show with numerical examples that, for smoothly laterally varying 2D media, the error caused by integrating over the monopole dimension is generally negligible — so long as the near-surface properties in the vicinity of sources and receivers do not have pronounced variations. In the case in which the near-surface properties vary rapidly laterally, in addition to obliquity factor, differences in local impedance at source and receiver locations must be taken into account.

The 2D model used in our test is shown in Figure ⁹. The full reflection data have 201 colocated sources (monopole) and receivers (dipole) with spacing intervals of $10m$ lying within the range of $x=(3000m,5000m)$ at the acquisition surface with the time sampling rate of $dt=2.4ms$⁠. The focal point is located at $(XF=4000m,ZF=1000m)$⁠. The coupled Marchenko system is solved with 13 iterations. Figure ¹⁰–¹⁰ shows the Green’s function obtained from SEM modeling, Marchenko redatuming by integrating on the dipole dimension and the monopole dimension, respectively. These Green’s functions correspond to waves generated from a monopole source sitting at the focal point and recorded at receivers at the acquisition surface, or vice versa, due to the source-receiver reciprocity. The focusing functions calculated with dipole and monopole integration as well as their difference are shown in Figure ¹¹–¹¹, respectively. For primary estimation, Figure ¹² shows the reference reflection data with the source position $Xs=4000m$⁠. Figure ¹²–¹² shows the result from dipole and monopole integration as well as their difference, respectively. Note that our observations are entirely consistent with well-established experience in SRME: amplitude effects — those related to the incidence-angle-dependent obliquity factors — are second order with respect to kinematic errors that result from sampling deficiencies.

To illustrate how acquisition sampling can impact results in a more realistic scenario, we study the effect of a limited aperture, together with source integration and subsampling in the redatumed Green’s function with a salt model, which has a constant density and varying wavespeed (Figure ¹³). This model has also been used by Vargas and Vasconcelos (2020) in the context of scattering-based focusing. We choose this model such that the effects of density-related impedance changes do not play a role, which allows us to focus on the kinematics and reconstruction of the events — in any case. Because sampling-related inaccuracies are controlled by the event slope (kinematics) and bandwidth, we expect that tests on a variable-density version of this model would yield the same conclusions as those we present later. The full reflection data have 201 colocated sources and receivers with spacing intervals of $40m$ in the subsalt model. The acquisition spread lies within the range of $x=(4000m,12,000m)$ at the acquisition surface with the time sampling rate of $dt=4ms$⁠. The focal point is located $(XF=8130m,ZF=4400m)$⁠, as shown in Figure ¹³. The same subsampling scheme is used as shown in Figure ². Figure ¹⁴ shows the Green’s functions obtained with the full aperture. As observed in Figure ¹⁴, subsampling by a ratio of eight brings strong artifacts for the complex salt model. A comparison between Figure ¹⁴ and ¹⁴ illustrates that the acquisition with smaller aperture also increases aperture-induced artifacts in the Green’s function. Figure ¹⁴ shows that a smaller aperture and subsampling compound into more complicated artifacts. Here, all of the Green’s functions are calculated by solving the coupled Marchenko equations by inversion with the LSQR solver (Paige and Saunders, 1982; Ravasi and Vasconcelos, 2020).

Completing our analysis of subsampling effects on redatuming, here we study the receiver subsampling effects on MPE in a four-layer model with colocated 201 sources and receivers, as shown in Figure ¹⁵. Figure ¹⁶ shows the reference reflection data with the source location $Xs=0m$⁠, where four strong primaries and the ensuing relatively weak multiples can be observed. Figure ¹⁶ shows the primary prediction result — using the algorithm presented in Algorithm 1 — without subsampling, where only four primaries are visible and the multiples are suppressed significantly. Figure ¹⁶–¹⁶ shows the results of primary estimation with three different subsampling schemes (evenly zeroing out 66%, 75%, and 80%), calculated with 20 iterative substitutions, respectively. Here, as previously, the choice of 20 iterations is arbitrary, chosen to be sufficiently higher than our observed number of iterations to achieve convergence on fully sampled data. Figure ¹⁶ shows the result using the same subsampling scheme as shown in Figure ¹⁶ but calculated with LSQR instead of explicit Neumann iterations. As expected, the strength of the artifacts increases with increasingly subsampled data. The errors in the first three iterations for the 75% subsampling scheme are given in Figure ¹⁷–¹⁷, which show that the major error comes from the first iteration and can decrease in subsequent iterations.

Figure ¹⁶ shows that the iterative substitution method fails to converge for this example when the subsampling degree increases to 80%. Based on the theory of Neumann series expansion, in Figure ¹⁸, we examine the operator norm of $A=Wϵt2−ϵR⋆Wϵt2−ϵR$ in equation 18 for different subsampling schemes as a function of time. In general, the operator norm first increases and then plateaus with time $t2$⁠. This occurs as more finite-energy arrivals are included by the windowing operator $Wϵt2−ϵ$ as the time increases, but the contribution of an increasing number of arrivals is offset by the decreasing energy of high-order multiples that arrive late (Zhang and Staring, 2018). As shown in Figure ¹⁸, the operator norm with full data for primary estimation with 20 iterations lies below that of the Marchenko redatuming for the chosen focal depth (0.8094), which satisfy the convergence condition $||A||<1$⁠. It is important to note that $||A||<1$ guarantees convergence, whereas for $||A||>1$⁠, it may diverge not guaranteed, but true for our example — for details on Neumann series convergence, we refer readers to Dukalski and de Vos (2018). Subsampling evenly by 66% and 75% brings the operator norm to approximately 1.2 and 2.1, respectively, where the iterations diverge slowly but lie within an acceptable range (see Figure ¹⁶ and ¹⁶). Subsampling evenly by 80% further pushes the operator norm to approximately 3.1, resulting in a much faster divergence as shown in Figure ¹⁶. Here, we point out again that this operator norm analysis does not carry the same implications when the problem is solved without relying on explicitly evaluating the Neumann-series-based solution.

Our study with the four-layer model shows that for Marchenko redatuming the subsampling effects jointly depend on (1) the choice of integration variable and (2) the subsampling dimension in the multidimensional convolutions within the Marchenko scheme. Because of source-receiver reciprocity, the freedom to choose an appropriate integration variable can be advantageous in cases in which only one of the dimensions (source or receiver) is subsampled. For the same reason, when the sampling of reflection data is symmetric with respect to the source and receiver dimensions, integrating on either dimension yields the same result, regardless of subsampling. When the data are regularly subsampled on one dimension, integrating on the other dimension yields inaccurate focusing results, with focusing function gaps at surface locations corresponding to the missing data. However, in this same sampling scenario, integrating on the subsampled dimension yields focusing functions with artifacts but no spatial gaps, with their strength affected by the chosen focusing depth. The analysis of the integrand gather shows that the mechanism behind lies at the integration within the multidimensional convolution between the reflection response and the focusing function formulated in the form of a matrix-vector multiplication in the frequency domain. Naturally, we expect that subsampling in sources and receivers will compound and cause stronger artifacts together with gaps, compared with subsampling in a single dimension. In addition, limiting acquisition aperture will also bring artifacts (Sripanich and Vasconcelos, 2019), which can be compounded by the subsampling effects, as we show with our salt model example.

Furthermore, we show with a model of mild lateral heterogeneity that conducting integration over the monopole dimension instead of the dipole brings negligible errors to both Marchenko systems investigated here, which may facilitate our choice of integration variable in combination with the subsampling scheme. This also reveals that the Marchenko methods are not highly sensitive to relatively small errors contained within the reflection data. Therefore, it implies that in practice interpolating the gaps of the input reflection data should be preferred over postprocessing the output focusing functions or Green’s functions, as shown by Jia et al. (2019) and Staring and Wapenaar (2020) for 3D redatuming.

One important take away from our results with different combinations of dipole and monopole field quantities is that, although theory strictly requires the use of dipole sources — which are not practically available — it is entirely feasible to rely on dipole data instead, e.g., those acquired with multicomponent sensors. Although replacing dipole sources by receivers may violate spatial reciprocity in heterogeneous media, we believe that, in most offshore conditions, the associated errors are likely negligible. In the absence of any dipole data, e.g., the case of conventional streamer data, we recommend that a monopole-to-dipole filter be applied in the best-sampled domain, likely to be the receiver domain in streamer-based data. In the case of ocean-bottom nodes, which are multicomponent but often highly sparse, it may be better to apply dipole corrections on the more densely sampled source domain and perform convolutions over sources.

In the context of focusing and redatuming, it is important to point out that the behavior of subsampling-related errors in Marchenko operations is analogous to those occurring in SRME. In SRME literature and practice (Dragoset et al., 2006, 2010; Verschuur, 2013), acquisition-related subsampling has been studied by many authors in many contexts — 2D, 3D, synthetic, and field data. Throughout the literature in multiple elimination, subsampling effects have been extensively investigated, together with several approaches suggested to condition data for multiple suppression — including interpolation/regularization techniques. For example, in closed-loop SRME (Lopez and Verschuur, 2015), one may use a focal transformation and a sparse norm regularization to remove the undersampling noise. Although our study begins to highlight the connections between these methods, we expect that more about Marchenko can be learned from what has been studied in the context of SRME — e.g., by taking inspiration from data conditioning approaches used in that context. However, appropriately translating SRME lessons into Marchenko practice will require further testing and research, particularly in the context of field data examples.

Although here we base our discussions on 2D examples, the arguments that we build regarding the integration over monopole fields are general, and as such in principle apply to 3D media, i.e., having source and receiver integrals each act on a plane of source/receiver coordinates. We point out that the work of Ravasi and Vasconcelos (2020) contains examples of subsampling 3D Marchenko operators — indeed in agreement with this study. Their example does further highlight that although aliasing effects in three dimensions follows the same patterns as we present here, they can be lower in magnitude when either source or receiver integration is subsampled in only one of the two spatial coordinates. This observation is likely to depend on the degree of subsampling and on the degree of complexity — due to the subsurface heterogeneity — present in the 3D wavefields that make up the reflection operator.

A study of subsampling effects on the primary estimation based on the projected Marchenko system shows a similar pattern of artifacts that appear in the original coupled Marchenko system, which can be explained by the common multidimensional convolution and correlation between reflection data and focusing functions. As expected, the artifacts increase with the degree of subsampling and may cause divergence issue with iterative substitution when subsampling to a certain extent, which is illustrated by the analysis of the time-dependent operator norm lying in the core of the iteration. This problem can be mitigated using more robust iterative linear solvers, such as LSQR — but this would likely come at a greater computational cost because such solvers tend to convergence more slowly than the Neumann series, thus requiring a larger number of iterations (i.e., of operator evaluations; Dukalski and de Vos, 2018). In our case, LSQR yields consistent results, as opposed to Neumann-series-based iterative substitution, when presented with data with high degrees of subsampling. When the data are well sampled, the Neumann-series iterative scheme and the LSQR solver perform equally well. As in using the coupled Marchenko system for redatuming, the errors caused by subsampling decrease rapidly with the number of iterations — provided that the Neumann series remains stable.

By relying on representative numerical cases to mimic acquisition scenarios over varying subsurface complexity, we provide a study of the acquisition-related subsampling effects on the 2D Marchenko focusing, redatuming (original coupled Marchenko system), and primary estimation (projected Marchenko system) to shed light on the underlying physical and numerical mechanisms. For Marchenko focusing and redatuming, we show that the artifacts brought by subsampling are rooted in the integration within the multidimensional convolution between the reflection response and the focusing function, which also applies to primary estimation because of how similar the underlying operations are. These subsampling effects depend jointly on the integration variable choice and the subsampled dimension (the dipole source or receiver). Artifacts appear without spatial gaps when these two are identical, and spatial gaps are present with artifacts for different dimensions. Our study also shows that integrating on the monopole dimension also gives very close approximation to the results of both Marchenko systems for smoothly laterally varying media. In this regard, future studies are expected regarding the sampling requirement for convergence, the interpolation accuracy, etc. The artifacts in redatumed Greens’ functions depend on the reflection data sampling and focusing depths. In imaging, these artifacts will very likely cause small errors, but for full-waveform inversion, time lapse, and reservoir characterization, these could lead to larger errors.

Our study focuses only on regularly subsampling in one dimension of the input reflection data, i.e., source or receiver dimension. Meanwhile, the conclusions drawn from our study are limited to the 2D case. For three dimensions, the irregular subsampling effects are likely more complicated, due to the 2D integration over the source and receiver dimension, particularly in the case of highly complex media. Our results can inform acquisition design, but more importantly they tell us how to handle single and multicomponent data for Marchenko applications, while providing insight that can help in ensuring the quality of data conditioning approaches for Marchenko applications.

We thank all of the sponsors of the Utrecht Consortium for Subsurface Imaging for funding this project. We appreciate the open source MATLAB code Spot for providing a linear-operator toolbox that allows us to compute the operator norm and the Salvus software for SEM-based modeling.

Data associated with this research are available and can be obtained by contacting the corresponding author.

Biographies and photographs of the authors are not available.