PDF VersionMarkdown Version

Least-squares extended imaging with surface-related multiples

R. Kumar\(^1\), N. Tu\(^{1}\), T. van Leeuwen\(^{2}\), and F. J. Herrmann\(^1\)
\(^{1}\) Seismic Laboratory for Imaging and Modeling (SLIM), University of British Columbia, Canada
\(^{1}\) Mathematical Institute, Utrecht University, The Netherlands


Common image gathers are used in building velocity models, inverting for anisotropy parameters, and analyzing reservoir attributes. Typically, only primary reflections are used to form image gathers as multiples can cause artifacts that interfere with the events of interest. However, it has been shown that multiples can actually provide extra illumination of the subsurface, especially for delineating the near-surface features. In this paper, we aim to form common image gathers directly from the data with surface related multiples by applying concepts that have been used to successfully deal with surface-related multiples in imaging. We achieve this by effectively inverting an extended migration operator. This results in extended images with better near-surface illumination that are free of artifacts that can hamper velocity analysis. In addition, being able to generate extended images directly from the total data avoids the need for (time-consuming) pre-processing. Synthetic examples on a layered model show that the proposed formulation is promising.


An extended image is a multi-dimensional correlation of source and receiver wavefields, as a function of subsurface offsets (see Sava and Vasconcelos, 2011, and references therein for a recent overview). There are many different applications in which extended images are used extensively. A prime example is the construction of angle-domain common-image gathers (ADCIGs) for migration velocity analysis (Biondi and Symes, 2004, Shen and Symes, 2008; Symes, 2008; Sava and Vasconcelos, 2011; Kumar et al., 2013, 2014). Extended images can also be used to study rock properties and fluid indicators by estimating the amplitudes of reflected waves as a function of incident angle at the interface. One of the current limitations is that the extended imaging conditions cannot handle primaries and multiple reflections simultaneously. An obvious way to deal with this is to separate the data into primaries and multiples first, and form separate (extended) images from them (Lu et al., 2014). However, the separation of primary and multiple wavefields is challenging, computationally expensive, and may damage the primary signal. Therefore, we would like to use primaries and multiples simultaneously when forming image gathers. To achieve this, we follow earlier work on joint imaging of primaries and surface-related multiples by (Tu et al., 2013) and (Tu and Herrmann, accepted Jan 2015) and adapt it to extended imaging. The way that multiples are used in this imaging scheme is remniscent of the EPSI (Estimation of Primaries via Sparse Inversion) approach, proposed by (Groenestijn and Verschuur, 2009). This approach maps surface-related multiples back to primaries via multi-dimensional deconvolutions. In this abstract, instead of mapping the data to the primary wavefield and subsequently mapping it to an extended image, we aim to map the total up-going wavefield to the extended image directly. Effectively, this produces a true-amplitude extended image that does not suffer from the typical coherent artifacts that are common when using correlation-based multiple-to-primary mapping. We show the potential benefits of the proposed formulation on a two-layer synthetic velocity model, using the full two-way wave-equation to produce the extended image.

Extended imaging with free-surface multiples

Given monochromatic source and receiver wavefields as data matrices \(U\) and \(V\), where each column represents the wavefield for a single source, the extended image is given by \[ \begin{equation*} E = VU^*, \end{equation*} \] where \(^*\) denotes the conjugate-transpose. An element \(e_{ij}\) of resulting matrix \(E\) captures the interaction between gridpoints \(i\) and \(j\). All conventional extended images can be extracted from this matrix, for example by selecting only elements with lateral interaction. We take the wavefields \(U\) and \(V\) to be the solutions of a time-harmonic wave equation \[ \begin{equation*} HU = P_s^TQ, \quad \mbox{and} \quad H^*V = P_r^TD, \end{equation*} \] where \(H\) is a discretization of the Helmholtz operator \((\omega^2 m + \nabla^2)\), \(m\) is the squared slowness, the matrix \(Q\) represents the sources, \(D\) is the data matrix and the matrices \(P_s, P_r\) sample the wavefield at the source and receiver positions (and hence, their transpose injects the sources and receivers into the grid). We can now write the extended image as \[ \begin{equation*} E = H^{-*}P_r^TDQ^*P_sH^{-*}. \label{eq:extimage} \end{equation*} \] Note here that \(H^{-*}\) represents an adjoint wave-equation solve, which constitutes the main computational cost in forming the extended image. However, if \(D\) contains primaries as well as surface related multiples this extended image will contain artifacts. In order to incorporate multiples in imaging, derived from the well-known SRME relation that the source \(Q\) should be replaced by an areal source \(Q-D\), which is effectively the total downgoing wavefield at the surface: \[ \begin{equation*} E = H^{-*}P_r^TD(Q^*-D^*)P_sH^{-*}. \end{equation*} \] The rationale behind this is that the multi-dimensional correlation \(DD^*\) is an approximation of the surface-related multiples, such that \(DQ^* - DD^*\) contains only primaries. However, such a correlation-based prediction is rarely accurate and the resulting extended image will contain artifacts. To counter this, we follow and invert the above relation for \(E\). To formalize this idea, we introduce an extended Born demigration operator which maps the extended image to the total upgoing wavefield \[ \begin{equation*} \mathcal{F}(E) = P_rH^{-1}EH^{-1}P_s^T(Q-D). \end{equation*} \] Evaluation of this extended demigrtion operator would require \(n_s + n_r\) wave-equation solves, where \(n_s,n_r\) are the number of source and receivers and \(\mathcal{O}(n^2)\) memory, where \(n\) is the number of subsurface gridpoints to store the extended image \(E\). This is prohibitive, even in 2D. To adress this, we follow and work with a subsampled extended image \(\widetilde{E}\) which is obtained by extracting \(l\) columns from \(E\). Effectively, we are computing the extended image only on \(l\) subsurface points, thus creating so-called image-point gathers (CIPs). Formally, we achieve this by multiplying \(E\) by the tall matrix \(W = [{\bf{w}}_1, \ldots, {\bf{w}}_{l}]\), where \({\bf{w}}_i = [0, \ldots, 0, 1, 0,\ldots, 0]\). Common image gathers (CIGs) can be constructed efficiently in a similar manner. The extended demigration operator is now given by \[ \begin{equation*} \mathcal{F}(\widetilde{E}) = P_rH^{-1}EWW^TH^{-1}P_s^T(Q-D), \end{equation*} \] and its evaluation requires \(2l\) wave-equation solves and \(\mathcal{O}(nl)\) memory. Typically, \(l\) can be one or two orders of magnitude less than \(n_s,n_r\) and \(n\) so this constitutes a significant saving. Finally, the least-squares estimation problem for \(\widetilde{E}\) is given by \[ \begin{equation} \label{ls} \underset{\widetilde{E}}{\text{minimize}} \quad {\frac{1}{2}\|\mathcal{F}(\widetilde{E}) - D\|_F^2}, \end{equation} \] where \(\|.\|^2_F\) is the Frobenius norm of the matrix (i.e., the sum of the squared entries). In the above expression we assume that the source signatures are known, however, in the proposed framework these can be estimated as well (A. Aravkin et al., 2013). Another possible extension is the use of source encoding to further reduce the computational cost (Tu and Herrmann, 2012).


To test the proposed algorithm for computing common-image gathers, we use a synthetic two-layer velocity model (with a grid sampling of 10 m) as shown in Figure 1(a). The sources and receivers are located on the surface with a spatial sampling of 10m. We use a 2D finite-difference frequency-domain code to generate the synthetic data sets. The source signature is a Ricker wavelet with a peak frequency 10 Hz. We solve \(\ref{ls}\) and reverse-time migration (Tu and Herrmann, accepted Jan 2015) using a few iterations of LSQR, a standard iterative linear least-squares solver. Figure 1(b) shows the background velocity model used to form least-squares reverse-time migration (RTM) images and common-image gathers (CIGs). This velocity model is kinematically correct, so we expect to see reflectors at the correct depth-position using primary reflections only in least-squares RTM images and focused CIGs. We construct common-image gathers (CIGs) at \(x = 1000\)m for all \(z\). Figures 3(a), 4(a) show the RTM image and CIG when we only use primary reflection data and \(Q\) as the source function. As expected, the least-squares migration and image point gather is fully focused. Next, we use the total up-going wavefield as data, but still use \(Q\) as the source function. Since the multiple reflections are not properly dealt with, we can see the focused ghost event generated by the multiple wavefields in Figures 3(b), 4(b). As mentioned before, the presence of such ghost events can mislead the velocity analysis process. Finally, we form the RTM image and image point gather using the total up-going wavefield along with the areal source function \(Q - D\). We can see in Figures 3(c), 4(c) that incorporation of areal sources in the least-squares RTM and extended images remove the ghost reflector.

Figure1(a)True velocity model. (b) Background velocity model used in extended imaging.
Figure2(a)Primary reflections only. (b) Total up-going data (primaries and multiples).

Figure3Least-squares reverse-time migration. (a) Primary data only. (b) Total up-going data is used but areal source is not used. We can clearly see the ghost reflector when we do not use the areal source. (c) Total up-going data is used along with the areal source. Incorporation of areal source remove the ghost reflector.
Figure4Least-squares common-image gather extracted along \(x = 1000\)m and for all \(z\). (a) Primary data only. (b) Total up-going data is used but areal source is not used. (c) Total up-going data is used along with the areal source. We can clearly see the effect of ghost reflectors (in middle) which can be remove (in right) via incorporating the areal source.


We have proposed a methodology for forming extended image gathers from data containing both primaries and multiples. The approach is based on earlier developments in simultaneous least-squares imaging of primaries and multiples, and matrix-probing methods to compute the extended image in a computationally efficient manner. The result is a method that creates a true-amplitude extended image gather directly from data containing surface related multiples, at a computational cost that is roughly equivalent to a least-squares migration. Numerical results on a simple synthetic model are promising and encourage further research.


This work was financially supported in part by the Natural Sciences and Engineering Research Council of Canada Discovery Grant (RGPIN 261641-06) and the 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, Statoil, Sub Salt Solutions, Total SA, WesternGeco, Woodside.

Aravkin, A., Leeuwen, T. van, and Tu, N., 2013, Sparse seismic imaging using variable projection: In Acoustics, speech and signal processing (iCASSP), 2013 iEEE international conference on (pp. 2065–2069). doi:10.1109/ICASSP.2013.6638017

Biondi, B., and Symes, W. W., 2004, Angle-domain common-image gathers for migration velocity analysis by wavefield-continuation imaging: Geophysics, 69, 1283. Retrieved from http://link.aip.org/link/GPYSA7/v69/i5/p1283/s1\&Agg=doi

Groenestijn, G. J. A. van, and Verschuur, D. J., 2009, Estimating primaries by sparse inversion and application to near-offset data reconstruction: Geophysics, 74, A23–A28. doi:10.1190/1.3111115

Kumar, R., Leeuwen, T. van, and Herrmann, F. J., 2013, AVA analysis and geological dip estimation via two-way wave-equation based extended images: SEG technical program expanded abstracts. doi:10.1190/segam2013-1348.1

Kumar, R., Leeuwen, T. van, and Herrmann, F. J., 2014, Extended images in action: Efficient WEMVA via randomized probing: EAGE. doi:10.3997/2214-4609.20141492

Lu, S., Whitmore, D., Valenciano, A., and Chemingui, N., 2014, Enhanced subsurface illumination from separated wavefield imaging: First Break, 32, 87–92.

Sava, P., and Vasconcelos, I., 2011, Extended imaging conditions for wave-equation migration: Geophysical Prospecting, 59, 35–55. doi:10.1111/j.1365-2478.2010.00888.x

Shen, P., and Symes, W. W., 2008, Automatic velocity analysis via shot profile migration: Geophysics, 73, VE49–VE59. Retrieved from http://dx.doi.org/10.1190/1.2972021

Symes, W. W., 2008, Migration velocity analysis and waveform inversion: Geophysical Prospecting, 56, 765–790.

Tu, N., and Herrmann, F. J., 2012, Least-squares migration of full wavefield with source encoding: EAGE. Retrieved from https://www.slim.eos.ubc.ca/Publications/Public/Conferences/EAGE/2012/tu2012EAGElsm/tu2012EAGElsm.pdf

Tu, N., and Herrmann, F. J., accepted Jan 2015, Fast imaging with surface-related multiples by sparse inversion: Geophysical Journal International. Retrieved from https://www.slim.eos.ubc.ca/Publications/Private/Submitted/2014/tu2014fis/tu2014fis.pdf

Tu, N., Aravkin, A. Y., Leeuwen, T. van, and Herrmann, F. J., 2013, Fast least-squares migration with multiples and source estimation: EAGE. doi:10.3997/2214-4609.20130727