Fast least-squares imaging with surface-related multiples: application to a North-Sea data set

Ning Tu and Felix J. Herrmann
University of British Columbia, Earth, Ocean and Atmospherical Sciences, Vancouver, British Columbia, Canada

Introduction

Marine seismic data usually contain strong surface-related multiples, which bounce more than once between the free surface (i.e., the air-water contact) and subsurface reflectors. Amongst these multiples, waves that bounce in the water column are often the most pronounced, which can lead to challenging short-period multiples. Acting on first impulse, surface-related multiples are typically removed during conventional seismic data processing. While these approaches, such as Surface-Related Multiple Elimination (SRME, Verschuur et al., 1992), have been successful in removing the imprint of multiples on the final image, they discard valuable information contained in these multiples whose wave-number contents and apparent aperture are enriched by the fact that multiples undergo multiple bounces. We illustrate this effect in Figure 1 . For this reason, we propose to embrace surface-related multiples in our imaging formulation so we can maximally leverage their complementary information.

The idea of using surface-related multiples by itself is not new. Berkhout (1993) presented with his \(\mathbf{W^-R^+W^+}\) model the theoretical foundation for later attempts to image multiples by modifying the source wavefield, either using the conventional cross-correlation imaging condition (i.e., RTM, Guitton, 2002; Schuster et al., 2004; Liu et al., 2011; Zuberi and Alkhalifah, 2013), or the deconvolutional imaging condition to control the coherent imaging artifacts associated with multiples using RTM (Muijs et al., 2007; N. D. Whitmore et al., 2010; Lu et al., 2014). As pointed out by Poole et al. (2010) and Tu et al. (2013b), even the latter approach may fail in the presence of complex geologies, as shown in Muijs et al. (2007), Figure (11b). Recent work by Lin et al. (2010), Verschuur (2011), and Wong et al. (2012) showed that these artifacts can be largely removed by casting the imaging-with-multiples problem as a proper iterative least-squares migration problem at the expense of expensive, and perhaps prohibitive in 3D, simulation costs. The main contributions of this work is fourfold, namely (i) we incorporate surface-related multiples in RTM via wave-equation solves w.r.t. areal sources, and thereby avoid expensive multiple predictions that plague SRME and Estimation of Primaries by Sparse Inversion (EPSI, Groenestijn and Verschuur, 2009; Lin and Herrmann, 2013); (ii) we estimate the source function on-the-fly during the iterations making use of multiples (Tu et al., 2013a), following ideas from variable projection and using the unique relation between primaries and multiples and the source (Verschuur et al., 1992); (iii) we remove cross-talk between multiples by promoting curvelet-domain sparsity in the model space; and finally (iv) we reduce the number of wave-equation solves to roughly the same number of a single RTM using all shots (Tu and Herrmann, 2015).

To demonstrate the possible advantages of our joint imaging (i.e., using both primaries and multiples), we apply this method (Tu and Herrmann, 2015) to a North-sea field data set. In particular, we are interested in specific benefits of imaging primaries and multiples together instead separately (e.g., in Guitton, 2002). After demonstrating the performance of the joint inversion approach on a stylized example, we apply it to a real North-sea data set.

Imaging with multiples & source estimation

For a given background-velocity model \(\mathbf{m}_0\), reverse-time migration (RTM) derives from the following linear relationship \[ \begin{equation} \delta\mathbf{d} = \nabla\mathbf{F}[\mathbf{m}_0;\mathbf{s}]\delta\mathbf{m} \label{eq:lmpri} \end{equation} \] between the observed upgoing primary wavefields collected in the vector \(\delta\mathbf{d}\), the unknown source wavefields collected in the vector \(\mathbf{s}\), and the unknown medium perturbation collected in \(\delta\mathbf{m}\). This linear relationship is based on the single-scattering approximation, embodied in the linearized Born scattering operator \(\nabla\mathbf{F}\).

During conventional imaging, surface-related multiples are removed from the observed data; the source is estimated, and a background-velocity model is built using (traveltime) tomography or migration-velocity analyses. In this paper, we assume the background-velocity model \(\mathbf{m}_0\) to be known while we estimate both the source signature \(\mathbf{s}\) and the image \(\delta\mathbf{m}\) from the total upgoing wavefield, which contains imprints of both the source signature and the surface-related multiples. Without loss of generality, we assume spatially impulsive sources with the same time-harmonic source signature \(\mathbf{s}_i=q(\omega)\) for all \(n_s\) sources \(i=1\cdots n_s\) with \(\omega\) the angular frequency.

Our formulation for imaging with multiples derives from three important modifications to conventional RTM where the source signature \(\mathbf{s}\) is assumed to be known—i.e., \[ \begin{equation*} \delta \mathbf{m}_\mathrm{rtm}=\nabla\mathbf{F}^\ast[\mathbf{m}_0;\mathbf{s}(q)]\delta\mathbf{d}, \end{equation*} \] where the image is obtained by applying the adjoint of \(\nabla\mathbf{F}\) to the de-multipled data—namely,

1. we replace the above spatially impulsive sources by areal sources that contains the downgoing data—i.e., we have \[ \begin{equation} \tilde{\mathbf{s}}(q)\mapsto \begin{bmatrix}\mathbf{s}(q)-\delta\tilde{\mathbf{d}}\end{bmatrix}. \label{eq:as} \end{equation} \] Note that now the observed data \(\delta\tilde{\mathbf{d}}\) also contains all upgoing surface-related multiples. Contrary to conventional imaging, the source wavefield is now made up of the downgoing point-source wavefield plus the downgoing data (i.e., the sign-reversed upgoing data as we approximate the surface reflectivity by \(-1\));

2. instead of carrying out migration with the areal source described earlier, we solve the curvelet-domain (denoted by the matrix \(\mathbf{C}\)) sparsity-promoting optimization program \[ \begin{equation} \begin{aligned} & \mathrm{minimize}_{\delta\mathbf{x}}\quad \|\delta\mathbf{x}\|_1 \\ & \mathrm{subject\ to} \quad\|\delta\tilde{\mathbf{d}}-\nabla\mathbf{F}[\mathbf{m}_0;\tilde{\mathbf{s}}(q)]\mathbf{C}^*\delta\mathbf{x}\|_{\mathrm{2}}^{\mathrm{2}} \leq \sigma^2, \end{aligned} \label{eq:l1inv} \end{equation} \] within some user-specified tolerance \(\sigma\). Similarly to Robust EPSI (REPSI, Lin and Herrmann, 2013), this sparsity-promoting program inverts the surface-related multiples but, instead of estimating the primaries, it produces a sparse curvelet representation for the multi-free least-squares image \(\delta\mathbf{x}_{\mathrm{LS}}\) with a give source estimate \(q\). The image is then readily derived by taking the inverse curvelet transform via \(\delta\mathbf{m}_{\mathrm{LS}}=\mathbf{C}^\ast\delta\mathbf{x}_{\mathrm{LS}}\);

3. we estimate the source function on-the-fly after the \(k^{th}\) iteration of the above optimization program—i.e, \[ \begin{equation} \mathrm{minimize}_{q}\quad\|\delta\tilde{\mathbf{d}}-\nabla\mathbf{F}[\mathbf{m}_0;\tilde{\mathbf{s}}(q)]\mathbf{C}^*\delta\mathbf{x}^k\|_{\mathrm{2}}^{\mathrm{2}}. \label{eq:src} \end{equation} \] As we will demonstrate, the proposed formulation produces high-resolution, high-fidelity, artifact-free images as well as accurate estimates for the source signature. The imprint of the source and surface-related multiples are removed by combining the linearized forward model with SRME. The latter relates the total upgoing wavefield to the surface-free Green’s function acting on the total downgoing wavefield. Compared to SRME our approach has the advantage of inverting, instead of removing, the surface-related multiples and thereby preserving multiple-scattering information. In addition, the proposed formalism uses wave-equation solvers to carry out the multi-dimensional convolutions implicitly (Tu and Herrmann, 2015), which can lead to considerable computational savings. However, as REPSI, our approach is based on an iterative inversion procedure that requires multiple linearized Born scatterings and migrations, which can be prohibitively expensive for large models. Before providing details how to reduce the computational cost, let us first describe common difficulties one encounters when imaging with surface-related multiples.

Challenges during imaging with multiples

Data preprocessing

As we mentioned earlier, our formulation derives from the SRME relation, which requires dense fixed-spread sampling with co-located sources and receivers. For marine towed-streamer acquisition, we fill in the missing half of data using reciprocity, and interpolate the missing (cross-line/near-offset) traces as in standard SRME. To obtain a reasonable estimate for this total upgoing wavefield, the direct waves and the downgoing wavefield at the surface (i.e., receiver ghosts) need to be removed, and the upgoing wavefield needs to be extrapolated from the receiver level to the free surface (Verschuur et al., 1992). The source-side ghosts, on the other hand, should be left intact so that we mimic a dipole source as required by the theory (Verschuur et al., 1992; Tu and Herrmann, 2015).

The need to invert

Contrary to previous attempts to mitigate the imprint of surface-related multiples, the proposed formulation calls, as in REPSI, for a sparsity-promoting inversion of the linear forward modeling operator that ties perturbations in the velocity model to data that contain surface-related multiples. Failure to invert this operator leads to images that suffer from the imprint of the source and, perhaps more importantly, from artifacts related to coherent interferences between the areal sources and the observed data due to the presence of multiples.

To explain the physical origin of these artifacts, we included in Figure 2 a cartoon of the ray paths of primaries associated with conventional RTM (Figure 2a) and RTM with areal sources (Figure 2b). In this plots, the star indicates the seismic source and the triangles indicate the receivers. Moreover, the solid arrows depict the ray paths of the primaries \(\mathbf{P}_0\), and the \(1^{st}\) to \(3^{rd}\)-order surface-related multiples \(\mathbf{M}_1, \cdots, \mathbf{M}_3\). We use solid circles to delineate imaging points that correspond to true reflectors, and dashed circles to delineate imaging points that produce phantom reflectors. As we can see from the illustration in Figure 2a, conventional RTM produces erroneous events that do not correspond to physical reflectors when the input data contains surface-related multiples. Only primaries contribute to imaging the true reflectors. (We ignore internal multiples for now.) Conversely, when we include the total down-going wavefield as an areal source, the RTM image contains a mix of proper and erroneous contributions from interactions between multiples present in the data and the areal sources. Neighboring — e.g., \(\mathbf{P}_0\) and \(\mathbf{M}_1\), and \(\mathbf{M}_1\) and \(\mathbf{M}_2\), etc. — pairs of multiples image to physical reflectors while the others create coherent interferences, which are similar to artifacts when carrying out conventional migration on data that contain surface-related multiples.

To illustrate interferences between surface-related multiples on the source and receiver sides, we include a stylized synthetic example consisting of a homogeneous water body, which also serves as the background model, and a single shallow perturbation as depicted in Figure 3a . In this model, we generate linearized primary-only data, by applying equation \(\ref{eq:lmpri}\) with spatially impulsive sources, and data with surface-related multiples by injecting the areal sources given by equation \(\ref{eq:as}\) . The output of conventional RTM with the primary-only data is of reasonable quality, as shown in Figure 3b, demonstrating RTM’s ability to reproduce subsurface structures from multiple-free data. Unfortunately, the same observation does not apply when we image data including surface-related multiples (juxtapose Figures 3b and 3c). In that case, the image contains erroneous multiple-related artifacts. Unfortunately, including the total downgoing wavefield as areal sources during RTM is insufficient as we can see the result still contains interferences that can lead to wrong interpretations (juxtapose Figures 3b and 3d).

As we will demonstrate below, these artifacts can be removed by properly inverting the multiple-prediction operator, now disguised as areal sources in the Born scattering operator, with curvelet-domain sparsity promotion as outlined above.

Imaging & source estimation with multiples

As the stylized example in the previous section clearly illustrates, using areal sources (in equation \(\ref{eq:as}\)) alone is obviously inadequate to image the surface-related multiples. Furthermore, RTM also needs prior knowledge of the temporal source signature of the spatially impulsive sources, which is typically unknown unfortunately. However, when we invert the linearized relationship with sparsity promotion, we obtain far superior results for the above example as shown in Figure 4 where we either know the source signature beforehand, or we estimate this source signature on the fly.

We can make the following observations: (i) sparsity-promoting inversion of data with multiples, but without using the areal sources, leads to a high resolution image that contains artifacts related to the multiples (see Figure 4a); (ii) when provided with the correct source wavelet, carrying out the inversion with areal sources leads to a high-resolution artifact-free image (see Figure 4b); (iii) the same applies when we invert for the wavelet on the fly (see Figure 4c ); (iv) source estimation produces accurate estimates for the source signature (see Figure 4d).

The above observations are truly remarkable because we have, following the procedure outlined above, been able to solve three problems in one go, namely (i) we successfully jointly map primaries and surface-related multiples to true reflectors so we maximally benefit from the higher Signal-to-Noise Ratio (SNR) of the primaries and the extra illumination of the multiples; (ii) we estimate the unknown source function accurately while removing its imprint; and (iii) we properly scale the contributions from primaries and multiples that both contribute to the image. While the results for this stylized example are certainly encouraging, we will need to first remedy the computational costs associated with this procedure, which relies on many iterations involving computationally expensive migrations and demigrations.

Fast randomized iterative inversion

Least-squares migrations, especially with sparsity-promoting, are in most cases computationally prohibitively expensive because they rely on expensive iterations, the simulation cost of which scales linearly with the number of (monochromatic) source experiments. To reduce these costs, we borrow ideas from Herrmann and Li (2012) and limit the total number of wave-equation solves by working with randomized subsets of frequencies and source experiments during each iteration. As described by the authors (Tu and Herrmann, 2015), this type of randomized procedure—where we regularly draw new randomized experiments for each subproblem to expedite the convergence—leads to significant cost savings bringing the required number of wave-equation solves of our iterative procedure down to roughly equal a single RTM with all sources. Remember that by virtue of using areal sources, our imaging procedure also avoids expensive iterative multi-dimensional convolutional operations required by SRME and EPSI. Instead, these expensive operations are carried out by the wave-equation solver (Tu and Herrmann, 2015), leading to what we believe a computationally viable approach that can be applied to 2D field data sets as shown below, and in the future, to 3D field data.

North-sea case study

To evaluate the performance of the proposed method, we consider a 2D relatively shallow (\(90 \mathrm{m}\) water depth) marine data set from the North Sea, acquired by PGS with dual-sensor streamers. Because this type of acquisition gives us access to both the pressure and the vertical particle-velocity wavefields, this data set is ideal for us to accurately separate the observed wavefield into its up- and down-going constituents. To get the best results, we apply the following preprocessing steps: near-offset data interpolation using the Radon transform; up-down wavefield decomposition to remove the receiver ghost and removal of the direct waves. These processing steps, in combination with application of source-receiver reciprocity to mimic a fixed spread acquisition, gives us an estimate for the downgoing wavefield with 401 co-located sources and receivers with 12.5m grid spacing that serves as input to our imaging algorithm. Only the first \(4092\,\mathrm{ms}\) of the data record are imaged up to 48Hz. We derive the \(3\times 5\mathrm{km}\) background-velocity model with 6.25m grid spacing from the stacking velocity. To assess the quality of our imaging results, we also ran REPSI to get independent estimates for the primaries and surface-related multiples. The least-squares (LS) image (we will use this term for images obtained by solving the above sparsity-promoting program) from the primary-only data is included in Figure 5a for reference. The image of the total upgoing data using the proposed method is shown in Figure 5b . For both images, we estimate the source wavelet on the fly. On the computational side, both images are obtained with roughly the same wave-equation solves as a single RTM with all the data (in the joint-imaging case, the source estimation step necessitates some extra data modeling that results in 50% more simulation cost, we in return reduce the number of iterations by one-third for Figure 5b).

Joint versus separate inversions

With the proposed method, we can either image primaries and multiples separately, or image the total upgoing wavefields in one step via joint inversion with on-the-fly source estimation. Compared to imaging with primaries and multiples separately, this joint approach has the advantage that the images from these two components are properly scaled as we optimize over the source wavelet to fit both primaries and multiples. As a result, we avoid the source-scaling ambiguity, which hampers imaging with primaries only (Tu et al., 2013a). In addition, we produce one combined imaging result that maximally uses the information residing in both primaries and multiples.

To substantiate this claim, we conduct a series of experiments. First, we compare images from primaries only (Figure 6a); from multiples only (Figure 6b), and from primaries and multiples jointly (Figure 6c). All results are obtained with our sparsity-promoting inversion method using their respective source wavefields (see Discussions in Tu and Herrmann, 2015). When we compare the images obtained from primaries and multiples only, we find that in general the image of the primaries alone is, as expected, cleaner. For example, the structures indicated by the red arrows are better resolved with the primaries, and so are the two sections marked inside the blue rectangles. However, while scaled differently, the image from the multiples alone reveals most of the major reflectors (cf. Figures 6a and 6b) and even leads to an improved image for the reflector indicated by the grey arrows in Figure 6b .

Comparing our results obtained by inverting primaries- and multiples-only data with the joint inversion suggests that inverting the total data may lead to superior results that benefit from the relatively noise-free primaries and extra illumination from the multiples. As the results in Figure 6c indicate, our joint inversion strikes a good compromise and produces a single image that combines the benefits of both primaries and multiples.

Suppressing imaging artifacts from multiples

To further substantiate our claim that the proposed inversion scheme can remove the coherent artifacts associated with multiples, we study the detrimental effects that surface-related multiples can have when handled incorrectly—i.e., when included but not inverted. For this purpose, we again draw your attention to imaging with multiples only by comparing the result of a single RTM (Figure 7a) with the result obtained with sparsity-promoting inversion (Figure 7b). In both cases, we use areal sources defined in equation \(\ref{eq:as}\) . Even though using areal sources leads to true reflectors such as the sea floor (c.f. the solid circles in Figure 2b), it also result in interferences (c.f. the dashed circles in Figure 2b), which rings across the section (Figure 7a). The presence of these artifacts exposes a fundamental shortcoming of RTM’s correlation-based imaging condition, which we overcome largely by sparsity-promoting inversion. As we can be readily observe from Figure 7b and the zoomed in sections A and B included in Figure 8a and 8b, the strong artifacts associated with the multiples in Figure 7a are mostly suppressed in Figure 7b . These results demonstrate that the energy of surface-related multiples can be mapped to true reflectors via a sparsity-promoting inversion procedure using the downgoing receiver data as areal sources. The results also support our claim that inverting primaries and multiples jointly allows us to benefit from the multiples without the need to separate primaries and multiples upfront, which is an computationally expensive proposition.

Conclusions

Because surface-related multiples bounce more than once between the free surface and the subsurface reflectors, they generate coherent artifacts when not treated properly during imaging—even in situations where conventional sources are replaced by areal sources that include the downgoing data. However, if we combine these areal sources with a proper sparsity-promoting inversion procedure—during which the deconvolved image and source signature are estimated—we arrive at a formulation that maximally leverages primary and multiple energy via a joint-inversion procedure. During this inversion, the adverse coherent artifacts arising from multiples are turned into coherent images, and we can benefit from the multiples that are better sampled than primaries as their bounce points act as secondary sources. Aside from mapping multiple-related artifacts to the correct reflectors, our method also avoids the expensive combination of wave-equation based surface-related multiple elimination and subsequent reverse-time migration of primaries and/or multiples. Instead, our method proposes to image the primaries and multiples jointly during an iterative sparsity-promoting inversion procedure. By adaptively estimating the source wavelet, we achieve an optimal balance between the contributions from primaries and multiples. By using randomized subsampling techniques, we limit the number of required wave-equation solves to roughly equal a single reverse-time migration with all data. Since the modeling of the multiple is done by the wave-equation solver, no explicit multiple prediction and removal are needed. Application of our inversion framework on a real marine field dataset clearly demonstrates the advantage of inverting primaries and multiples jointly compared to estimating and imaging these wave constituents separately.

Acknowledgements

The author would like to thank PGS for providing the data set and granting us the permission to publish the imaging results. We also thank Dr. Shaoping Lu and Steve Kelly from PGS for beneficial discussions, and Prof. Eric Verschuur for preprocessing the data set and providing the stacking velocity.

This work was financially supported in part by the Natural Sciences and Engineering Research Council of Canada Collaborative Research and Development Grant DNOISE II (CDRP J 375142-08). This research was carried out as part of the SINBAD II project with support from the following organizations: BG Group, BGP, Chevron, ConocoPhillips, CGG, ION GXT, Petrobras, PGS, WesternGeco, and Woodside.