Imaging with surface-related multiples to overcome large acquisition gaps

To get the best result for seismic imaging using primary reflections, data with densely-spaced sources and receivers are ideally preferred. However, dense acquisition can sometimes be hindered by various obstacles, like platforms or complex topography. Such areas with large data gaps may deter exploration or monitoring, as conventional imaging strategies would either provide poor seismic images or turn out to be very expensive. Surface-related multiples travel along different paths compared to primaries, illuminating a wider subsurface area and hence making them valuable in case of data with large gaps. We propose different strategies of using surface-related multiples to get around the problem of imaging in the case of a large data gap. Conventional least-squares imaging methods that incorporate surface-related multiples do so by re-injecting the measured wavefield in the forward-modelling process, which makes it still sensitive to missing data. We introduce a ‘non-linear’ inversion approach in which the surface multiples are modelled from the original source field. This makes the method less dependent on the receiver geometry, therefore, effectively exploiting the information from surface multiples in cases of limited illumination. However, such an approach is sensitive to the knowledge of the source properties. Therefore, we propose a ‘hybrid’ method that combines the non-linear imaging method with the conventional ‘linear’ multiple imaging method, which further improves our imaging result. We test the methods on numerical as well as field data. The results indicate substantial removal of artefacts in the image derived from linear imaging methods due to incomplete data, by exploiting the surface multiples to a maximum extent.


Introduction
Seismic data can sometimes end up with large gaps due to several reasons such as bad topography, lack of access, infrastructure (e.g. an obstructing platform), legal or environmental considerations or failing equipment. If we use conventional seismic imaging algorithms that rely on primary reflection data only, migration artefacts become unavoidable (Hou & Symes 2017). A common practice in such scenarios is re-acquiring the data specifically targeting the missed zone. For instance, to illuminate areas under physical obstructions such as platforms, undershooting has become a standard practice (Hill 1986;OGP 2011). However, the large offsets between the sources and receivers take a toll on the vertical resolution of the data and as a result more efforts and expense are put in to acquiring high-resolution data. For example, Games et al. (2015) describe a high-resolution seismic undershooting survey to specifically address the loss of vertical resolution in the case of large offset data in conventional two-streamer undershooting around physical obstructions. This calls for alternative imaging strategies that make the most out of the available data to mitigate this problem.
Although in the past multiples have been deemed as noise, it is well acknowledged that multiples provide better seismic illumination, which can help in overcoming the limitations of the employed acquisition geometry (Berkhout & Verschuur 1994;Guitton 2002;Muijs et al. 2007;Liu et al. 2011;Lu et al. 2013;Zuberi & Alkhalifah 2013;Zhang & Schuster 2014).
Free-surface multiples travel along different propagation paths compared to primaries. This property is useful as they not only illuminate areas that are unreachable by the primaries but also contribute to the vertical resolution because of improved angular illumination . Therefore, in cases with limited illumination, migration methods that include surface-related multiples have an advantage (Tu et al. 2011;Makinen & Rickett 2016). Surface-related multiples can be migrated by considering the re-injected total measured response R ∩ ⃗ P − as the illuminating source wavefield (Verschuur & Berkhout 2011). This helps in migrating the multiples by establishing a linear relationship between the multiples and their sub-event primaries. Since the measured data are used both as incident and reflected field, calibration of source wavelets is no longer necessary (Davydenko & Verschuur 2017). Surface-related multiples have also been included in migrating the total measured response by considering the original source wavefield ⃗ S + and the re-injected total measured response R ∩ ⃗ P − as the total illuminating source wavefield ⃗ Q + (Berkhout 2014b;Lu et al. 2018). This helps in migrating both the primaries and the multiples simultaneously, thereby eliminating the need for multiple removal for migration. However, since these methods rely on re-injection of recorded data as incident wavefield, it is assumed that the complete downgoing wavefield has been recorded. Thus, for the area under large gaps, current methods have been incapable of exploiting the power of multiple scattering.
In this paper, we present a 'non-linear' imaging scheme that compensates for unrecorded data by synthetically modelling all the surface-related multiples, starting from the original source wavefield. In this way the dependence on receiver geometry becomes weaker, allowing us to overcome large acquisition gaps. This can be viewed as the extension of the least-squares migration method for primaries (Nemeth et al. 1999) or by using multiples (Brown & Guitton 2005;Zhang & Schuster 2014;Liu et al. 2016) to overcome incomplete data issues specifically by focusing on exploiting surface multiples in various ways. Furthermore, it relates to the model-based multiple imaging method proposed by Jiang et al. (2007), except with a closed-loop approach. We also discuss a hybrid method that combines the existing methods along with our proposed non-linear imaging scheme in order to exploit the benefits of both the methods. We demonstrate this methodology on finitedifference modelled data as well as field data from the Vøring basin.
Note that in this paper we are only concerned with the use of surface-related multiples in imaging. However, the forward model that we describe can be extended to include internal multiples (Berkhout 2014a;Davydenko & Verschuur 2017). In that case, internal multiples from primaries and internal multiples from the surface-related multiples are included as well. Such extension will be limited to internal multiples that travel in a down-up-down mode (so with three, five, seven, etc. subsurface bounces). Zuberi & Alkhalifah (2014) and Aldawood et al. (2015) describe how imaging could be extended to include internal multiples from so-called duplex waves (i.e. with two bounces in the subsurface). Marchenko-based methods for internal multiple imaging Slob et al. 2014;da Costa Filho et al. 2015) or internal multiple removal (Zhang & Slob 2019) are usually applied after surface multiple removal, so they might be less adequate in this case.

Review of migration schemes that use surface-related multiples
We use the discrete matrix notation as given by Berkhout (1982) to describe the forward modelling in the frequency domain (for a given frequency). The monochromatic wavefield for one source experiment measured at the depth level z m is indicated as ⃗ P(z m ). Wavefield operators that connect wavefields between two levels z m and z n are written as a matrix. For example, the matrix W + (z m , z n ) denotes a monochromatic one-way forward wavefield propagation operator that extrapolates wavefield from depth level z n to z m . The superscripts + and − describe downward and upward propagation, respectively.
In order to include surface-related multiples in a migration method, Verschuur & Berkhout (2011) made some modifications to the incident wavefield. Instead of using only the source wavefield ⃗ S + (z 0 ) generated at the surface, the total downgoing wavefield at the surface ⃗ Q + (z 0 ), comprising both the original source wavefield ⃗ S + (z 0 ) and the downward reflected total wavefield, is used as the incident wavefield: Here ⃗ P − (z 0 ) refers to the total 'recorded' wavefield at the acquisition surface and R ∩ (z 0 , z 0 ) represents the downward reflectivity operator at the surface z 0 , which is usually taken as R ∩ (z 0 , z 0 ) = −I (Verschuur & Berkhout 2011). (We 743 use the upgoing wavefield at the surface after receiver deghosting.) The downgoing wavefield ⃗ P + (z m , z 0 ) at depth level z m is given by forward extrapolating ⃗ Q + (z 0 ) to z m : while the upgoing wavefield ⃗ P − (z m , z 0 ) at depth level z m comes from inverse-extrapolating the total recorded wavefield ⃗ P − (z 0 ) to z m : Similar to imaging using primary data, every event now (primaries and surface-related multiples) is accounted for in the primary event and primary subevent of the multiple, therefore leading to a 'linear' relationship between the forward modelled data and reflectivity [if the effect of internal multiples is excluded-note that if full wavefield migration (Berkhout 2014b) is used, internal multiples will be included]. Reflectivity can now be derived approximately using one of the traditional imaging principles [e.g. by the correlation of ⃗ P + (z m , z 0 ) and ⃗ P − (z m , z 0 )], thereby migrating the primaries and multiples simultaneously (Berkhout & Verschuur 1994;Verschuur & Berkhout 2011;Zuberi & Alkhalifah 2013;Weglein 2015). Since conventional imaging conditions are not equipped to handle complex re-injected wavefields, it is replaced by a least-squares minimisation process with an iterative inversion method (i.e. a closed-loop approach) (Davydenko & Verschuur 2012;Verschuur & Berkhout 2015). Figure 1a depicts the flowchart of the closed-loop process. In each iteration for a given velocity model and a starting reflectivity model, we compare our forward modelled data with the measured data to obtain the residual (P obs − P mod ). The field data is re-injected in the forward modelling to explain the multiples. The iterative minimisation method takes care of the cross-talks generated with every progressing iteration. This produces an image that explains the multiples properly; given the forward model includes sub-event primaries. The cost function J i of the iterative minimisation process is given by (4) For a given iteration i, the current reflectivity operator matrix R ∪ i is used to model seismic data using the propagation operators W (Berkhout 2014a). (We refer to 'R ∪ i ' as 'R i ' from here on, as we do not consider downward reflectivity in our modelling, except at the free surface.) The scattering is generated by the reflectivity model while the propagation effects are handled by the velocity model. In these methodologies, since we assume the reflection parameter to be angle-independent, R i is taken to be a diagonal matrix. The residual data is then used to update the reflectivity, thereby also suppressing possible cross-talk (Berkhout 2014b). In the case of multiples-imaging, the method can overcome the limitations of coarse source sampling (Lu et al. 2014;Davydenko & Verschuur 2017). However, as the recorded wavefield becomes the incident wavefield for imaging, the performance of the method is dependent on good receiver coverage and density.

Incorporating surface multiples in a non-linear inversion method
Addressing the limitations of the linear inversion method in the case of missing data, we introduce a non-linear inversion scheme. In this method, all surface-related multiples are included incrementally in the incident wavefield using the modelled upgoing wavefield (Berkhout & Verschuur 2014) leading to a non-linear relationship between reflectivity and the forward modelled data. Since this method omits the reinjection of measured data, it is less affected by missing data. Due to the non-linear relationship, we also expect higher sensitivity to changes in the reflectivity model, as the errors in estimated reflectivities are amplified in higher-order multiples.
In this method we perform forward modelling with ⃗ S + (z 0 ) as our initial incident wavefield. In subsequent iterations, we modify the downgoing wavefield [ ⃗ Q + (z 0 )] i such that it includes the modelled surface multiples, based on reinjection of the modelled upgoing wavefield ⃗ P − m (z 0 , z 0 ) from the previous iteration. Therefore [ ⃗ Q + (z 0 )] i is given by where i refers to the iteration number. After making suitable changes for non-linear inversion in a closed-loop approach, we then derive the reflectivity image using the same minimisation method mentioned in equation (4). Figure 1b depicts the flowchart for this method. A small change is made in the way multiples are modelled in the non-linear closed-loop method from the linear closed-loop method (figure 1a), as can be seen in the flowchart in figure 1b. The simulated data instead of the measured data is being re-injected to get the multiples. Table 1 compares the different wavefields used for imaging in the methods described above.

Hybrid approach
While the non-linear method performs better in the case of incomplete sampling, it comes with its own shortcomings. As the relationship between the forward modelled data and reflectivity is non-linear, this method is naturally more sensitive to errors. This means that, although potentially a higher resolution can be obtained, the inversion process can be trapped in a local minimum. Since the non-linear method requires modelling the data starting from a source wavefield, the image is quite sensitive to the errors in the source wavelet. In contrast, the linear imaging method does not require knowledge of the wavelet because measured data is used as both incident and reflected field, thereby making it more robust. As the linear and non-linear methods complement each other on their drawbacks and benefits, a hybrid approach that combines the two methods by successively applying the methods, i.e. using a result from one method and putting it as an input for the other, will give the best solution in cases of missing data. The suggested strategy is as follows (Nath & Verschuur 2017).
(i) As the non-linear inversion method creates a better modelled data in the case of acquisition gaps, we first start with the non-linear inversion process on the data with gaps. (ii) We then use the modelled shot record gather from the non-linear inversion method to fill in the gaps of the measured seismic data. (iii) As the linear inversion method is more robust, we use this infilled data in the linear inversion method to obtain an intermediate subsurface image.
(iv) Finally, as the non-linear inversion method explains the missing data better, we feed this reflectivity image received from the previous step as a jump start to the non-linear inversion method along with the original measured data. This provides the final image.

Numerical example 1: FWMod data ('inverse crime')
Here, we illustrate the different methods using a synthetic model ( 2a). For numerical example 1, a full wavefield synthetic data is generated using full-wavefield modelling (Berkhout 2014a, FWMod) on this 2-D model, with the restriction of angle-independent reflectivity in both modelling and imaging. Using such an 'inverse crime' approach we can investigate the maximum expected performance of the proposed methodology. Figure 2c shows one such shot record with a gap in the data. Figure 3a shows the closed-loop image obtained from primaries-only data. The effects in the image due to missing receivers is clearly visible as the reflectors and shallow diffractors around the gap are not imaged properly. Figure 3b shows the imaging result using the primaries and multiples in a closed-loop linear inversion method mentioned earlier. We do see some improvements in this image compared to figure 3a, however the effect due to the gap is still evident. Figure 3c illustrates the imaging result from the proposed non-linear inversion method. The deeper reflectors as well as the shallow scatterers are better imaged here. It can be seen in figure 4a that the data are not modelled well near the gap, causing a poor minimisation of residual data at far offsets, as shown in figure 4b. The modelled data from the non-linear method has improved around the area with missing data as well as at the far offsets of the first-order multiples ( therefore, we see an improvement in the residual data as well (see figure 4d).
The modelled shot record gather from the non-linear imaging method (figure 5b) is used to fill the gaps in the measured seismic data, as shown in figure 5c. These shot record data are then used in our linear imaging method, which is more robust to get a reflectivity image. Figure 5d shows this intermediate imaging result of the hybrid approach. Since the non-linear inversion method better explains the missing data, in the final step of the hybrid approach, we feed this reflectivity image (figure 5d) to the non-linear inversion method as a jump-start along with the original measured data to give the result shown in figure 5e. Despite a very large gap, the method manages to fill in the image with the help of modelled as well as measured multiples. Figure 4e shows the final modelled data we get from the hybrid approach. The complementary residual data from this hybrid approach, displayed in figure 4f, shows that the modelling process performs successfully in explaining the data around the gap and also for the larger offsets. much better results than the other two approaches. The shallow diffractors are now imaged properly, the reflector peaks are also higher and sharper, as shown in figure 6c, and the residuals are reduced compared to the other methods mentioned previously.

Numerical example 2: finite difference data
In this example, synthetic data are generated using acoustic finite difference modelling on the model shown in figure 7. Figures 7a and 7b are the corresponding velocity and density models used to generate the data. Figure 7c is the resulting true reflectivity image. We use the four shots at the surface at 600 m, 900 m, 4500 m and 4800 m. Receivers are spread throughout the surface except for from x = 2400 m to x = 3000 m. Figure 8a shows the imaging result using the primariesonly data. The effects in the image due to missing receivers is clearly visible. Figure 8b shows the imaging result using primaries and multiples via the closed-loop linear inversion method mentioned earlier. Upon comparison with figure 8a, the benefits of using the free-surface multiples are already visible here. The improvements in the image can be noticed not only at the shallowest reflector but also in the deeper reflectors. Figure  be seen better imaged around the area with the receiver gap (marked with arrows). Figure 8d shows the final imaging result obtained via the hybrid method. A considerable improvement is seen upon comparison of the image in figure 8d with figures 8a or 8b. The modelled data from the linear imaging method (see figure 9a) are not modelled well around the missing receivers; the corresponding residual data (see figure 9b) also shows the same. The modelled and the residual data from the non-linear method show quite an improvement (figures 9c and 9d). Improvements can also be seen in the residual data from the hybrid method (figure 9f) compared to the one from the linear method, indicating better modelling of data. However, note that the improvements are less prominent as in the inverse crime example 1. This is because the forward-modelling algorithm used in migration (FWMod) is different from the modelling algorithm used to produce the data (finite difference algorithm), making it more realistic. This suggests that a better handling of angledependent reflectivity is necessary, as in this implementation of FWMod the reflectivity is considered angle-independent, while the input data has all amplitude variation with offset (AVO) effects.

Field data example
In this section, we compare the proposed methods on data from the Vøring basin in the Norwegian sea. We select a sub-set over 5800 m of recording and use four shot records with sources at 1600 m, 2100 m, 3600 m and 4100 m. Receivers are spread along the surface at a spacing of 25 m except for from x = 2600 m to x = 3600 m, yielding a 1-km receiver gap. Figure 10 shows one of these shot records. Note that we used the fully sampled data to construct split-spread shot records and near-offset interpolation to make initially completely sampled gathers (Kabir & Verschuur 1995). Such preprocessing was required for application of the Estimation of Primaries by Sparse Inversion (EPSI); see also Van Groenestijn & Verschuur (2009). The source wavelet in our example was derived using the EPSI process. For this experiment, the four selected shots could represent an ocean bottom node-type geometry, although sources and receivers are at the surface. Figure 11a shows the imaging result using primaries-only imaging. The effects in the image due to missing receivers can be seen below the gap area, marked with arrows. Figure 11b shows the imaging result using the closed-loop linear inversion method (using primaries and multiples). The contribution of multiples in filling the image gaps and improving the resolution have been marked with arrows. Figure 11c illustrates the imaging result from the nonlinear inversion method. It can be seen that the reflectors are better-imaged around the indicated area. Figure 11d shows the final imaging result of the hybrid approach. The reflectors are visibly better-imaged not only around the area with gap, but also at the outer locations. Figures 12a and 12c  Finally, in figure 13 we display the same shot record as figure 10 with the reconstructed data in the gap. Although the reconstruction for the surface multiples (starting at t = 4 seconds) is somewhat weaker, we see a proper reconstruction of the primaries. Note that this information originated from the multiples.

Discussion
Including multiples in the imaging method comes with a lot of benefits but meanwhile also introduces several problems unique to multiples. A closed-loop imaging method becomes necessary to minimise the resulting cross-talk from multiples (Davydenko & Verschuur 2017;Lu et al. 2018). Figure 14 shows the comparison of objective function values as a function of iterations for imaging methods using primaries-only, primaries and multiples in the linear method, and primaries and multiples in the non-linear method from Example 1. This demonstrates that including multiples will require more iterations to reach convergence than the method with primaries-only data. Since the non-linear method simultaneously builds the downgoing wavefield along with migration in every iteration, it requires a few more iterations to catch up with the linear method using primaries and multiples. Time (s) Figure 10. Field data example: measured shot record data for the source at 2100 m.
In this paper, we demonstrate the methods only in 2-D cases. Further merit of this method can be seen if we test the method in more realistic situations, i.e. using 3-D acquisition.
In our method, we exclusively utilise the surface-related multiples to image the area under a gap. Imaging with internal multiples (Malcolm et al. 2009;Davydenko & Verschuur 2013;Zuberi & Alkhalifah 2013 could also give way to illuminating such areas using a similar principle of iterative modelling. However, internal multiples are usually weaker than surface multiples, so exploiting their extra illumination is challenging. In our examples, we have considered reflectivity to be a scalar quantity. Since we will be using the non-linear imaging method in cases of gaps, introducing angle-dependent reflectivity (extending the parameters) may lead to over-fitting. Hence, extending the reflectivity parameter to account for AVO effects may require certain constraints in the reflectivity matrix R.

Time (s)
Receiver location (m) Figure 13. Field data example: modelled shot record data for the source at 2100 m after putting an infill in the gap obtained from the non-linear inversion process.
The use of the hybrid method still requires an accurate estimation of the source wavefield. For our field data example, this wavefield could be estimated from the surface multiples via the EPSI method. Also, it can be extracted from calibrating the primaries-only image with the linear imaging from the multiples, as described by Davydenko et al. (2015). More research is required for a stable and robust source wavefield estimation in this context. Primary wavefields are stronger than their associated multiples, and therefore they dictate the imaging result in methods that include both the primaries and surface multiples. In case of velocity errors in the model, linear and non-linear methods that use both primaries and surface multiples are dominated by the primaries in the imaging result. Figure 15 shows imaging results using the same data of example 1 with a 4% increase in the velocity model used for imaging. Figures 15a, 15b and 15c show the imaging result with sparse sources using the primaries-only method, using the linear primaries and multiples method and the nonlinear primaries and multiples method, respectively. The dominance of primaries in the imaging result can be clearly seen, albeit with a very small improvement over methods that use multiples. To investigate the effect of velocity errors on primaries and multiples separately, figure 16a shows the imaging result with the wrong velocity model using the primaries-only method with a dense coverage of sources and receivers. Figure 16b shows the imaging result using the multiples-only imaging method with sparse sources.  Figure 15. Imaging with 4% velocity error from example 1 using a sparse source configuration indicated with the red dots. Imaging result using closedloop imaging process (a) using primaries-only; (b) using primaries and surface multiples in a linear imaging method; (c) using primaries and surface multiples in the proposed non-linear inversion method. 756 Figure 16. Imaging with 4% velocity error from example 1. Imaging result using a closed-loop imaging process (a) using primaries from a densely sampled source configuration; (b) using multiples only from a sparse source configuration with a linear multiples-only imaging method. This image shows a high resemblance to figure 16a, despite the sparse sources. In spite of intuition, this seems to show that in practice, multiples are less sensitive to imaging with velocity errors in cases of acquisition restrictions. One of the reasons may be that the multiples tend to propagate under smaller angles than primaries for a similar offset. While a larger offset is ideal for velocity estimation methods, multiples can contribute in another key way; by providing a better and robust image. This helps in avoiding any local minima. This aspect can be strategically utilised for better illumination and velocity estimation in methods like joint migration inversion ) and wavefield tomography (Soubaras & Gratacos 2018).

Conclusions
We proposed different strategies for imaging with multiples in case of large acquisition gaps. We developed a non-linear imaging method in which the forward modelled data (including surface-related multiples) are modelled recursively using the original source wavefield. We choose this method as it is less dependent on receiver density compared to imaging methods that include surface-related multiples via re-injecting the measured data. The non-linear method, however, is more sensitive to velocity errors, source wavefield knowledge and forward-modelling imperfections compared to the linear imaging methods.
We tested these methods on 2-D synthetic data as well as field data; we see that it performs well in infilling the gap in the data even if these gaps are large. Since the linear imaging method-via measured data re-injection-has some advantages over the non-linear method, a hybrid approach using both methodologies further improves the result while tackling the limitations of the linear and non-linear methods. The merits of the hybrid method have been demonstrated on both the synthetic and field data.