3D Marchenko internal multiple attenuation on narrow azimuth streamer data of the Santos Basin, Brazil
33D Marchenko internal multipleattenuation on narrow azimuth streamerdata of the Santos Basin, Brazil
Myrna Staring and Kees Wapenaar
Abstract
In recent years, a variety of Marchenko methods for the attenuation ofinternal multiples has been developed. These methods have been exten-sively tested on 2D synthetic data and applied to 2D field data, but onlylittle is known about their behaviour on 3D synthetic data and 3D fielddata. Particularly, it is not known whether Marchenko methods are suffi-ciently robust for sparse acquisition geometries that are found in practice.Therefore, we start by performing a series of synthetic tests to identifythe key acquisition parameters and limitations that affect the result of 3DMarchenko internal multiple prediction and subtraction using an adaptivedouble-focusing method. Based on these tests, we define an interpolationstrategy and use it for the field data application.Starting from a wide azimuth dense grid of sources and receivers, aseries of decimation tests is performed until a narrow azimuth streamergeometry remains. We evaluate the effect of the removal of sail lines, nearoffsets, far offsets and outer cables on the result of the adaptive double-focusing method. These tests show that our method is most sensitiveto the limited aperture in the crossline direction and the sail line spacingwhen applying it to synthetic narrow azimuth streamer data. The sail linespacing can be interpolated, but the aperture in the crossline direction isa limitation of the acquisition.Next, we apply the adaptive Marchenko double-focusing method to thenarrow azimuth streamer field data from the Santos Basin, Brazil. Internalmultiples are predicted and adaptively subtracted, thereby improving thegeological interpretation of the target area. These results imply that ouradaptive double-focusing method is sufficiently robust for the applicationto 3D field data, although the key acquisition parameters and limitationswill naturally differ in other geological settings and for other types ofacquisition.
The Santos Basin in Brazil is known for its oil-bearing carbonate reservoirs belowa highly reflective stratified salt layer (see figure 1). The salt layer generatesstrong internal multiples that pose a problem for seismic imaging (Cyprianoet al., 2015). Most imaging methods assume that the recorded wavefield was1 a r X i v : . [ phy s i c s . g e o - ph ] A p r nly reflected once and thus incorrectly interpret internal multiples as primariesfrom deeper reflectors. As a result, these methods create ghost reflectors that donot exist in reality. These ghost reflectors can interfere with the real reflectors inthe target area and thereby corrupt the image. Therefore, we wish to attenuateinternal multiples in order to obtain a reliable image of the target area.The attenuation of internal multiples is a challenge. Various methods havebeen proposed, ranging from filtering methods (e.g. Hampson, 1986; Foster andMosher, 1992; Zhou and Greenhalgh, 1994) that transform the reflection re-sponse to an alternative domain in which the primaries and the internal multi-ples separate, to wave equation-based methods that aim to predict the internalmultiples by convolving and correlating the reflection response with itself (e.g.Jakubowicz, 1998; Weglein et al., 1997). The application of filtering methodsis often challenging in settings with a complex overburden, since there is usu-ally no distinct difference in properties between the primary reflections fromthe target area and the (strong) internal multiples generated in an overbur-den. Also, the application of wave equation-based methods is not undisputed.Some wave equation-based methods require the manual identification of inter-nal multiple generators, thereby introducing bias and the risk of not correctlycapturing all internal multiple generators in the process. In addition, some waveequation-based methods predict internal multiples with incorrect amplitudes oruse a layer stripping approach that results in error accumulation from shallowto deep reflectors. In order to attenuate internal multiples in a complex settingsuch as the Santos Basin, an alternative method is needed.Marchenko methods (Ware and Aki, 1969; Broggini et al., 2012; Wapenaaret al., 2014) are data-driven and wave-equation-based methods that do not havethese drawbacks. These methods have the ability to consider the entire over-burden as a whole, instead of having to identify all individual internal multiplegenerators separately. In addition, Marchenko methods allow us to retrieveGreen’s functions including primaries as well as all orders of internal multi-ples at any desired depth level without having to resolve overlying layers first.When writing the retrieval of Green’s functions using the coupled Marchenkoequations as a Neumann series, Marchenko methods can be used for the pre-diction of internal multiples (van der Neut et al., 2015b). These predictionsin principle have the correct amplitude and phase. However, minor amplitudeand phase differences are usually present when applying the method to fielddata due to imperfect acquisition or preprocessing. A mild adaptive filter canbe used to correct for these minor differences. We previously reported on thesuccessful application of an adaptive Marchenko method (the adaptive double-focusing method) to 2D synthetic data and a 2D line of streamer data of theSantos Basin in Brazil (Staring et al., 2018b). Internal multiples were predictedand adaptively subtracted from the target area, which improved the geologicalinterpretation. In addition, we found that the adaptive double-focusing methodwas relatively robust for a sparse acquisition geometry in 2D and suitable for theapplication to large data volumes. In the hope that these properties also holdin 3D, we use this adaptive Marchenko method for the prediction and adap-tive subtraction of internal multiples from 3D narrow azimuth streamer data2cquired in the Santos Basin.The extension from 2D to 3D Marchenko methods may seem trivial in the-ory, but it is not the case in practice. Some aspects are similar, such as the datapreparation requirements that include noise suppression, signature deconvolu-tion, deghosting and the removal of surface-related multiples. However, aspectsrelated to the sampling of the acquired data are different. In addition to the in-line direction in 2D, there is in 3D also a crossline direction that typically has alimited aperture and less densely spaced sources and receivers. Also, streamersusually do not record responses at negative offsets, near offsets and far offsets inthe inline direction. A thorough understanding of the effect of these acquisitionlimitations on the result of Marchenko internal multiple attenuation would al-low us to estimate whether the application to any particular dataset is feasible.In addition, it would aid us in defining an interpolation strategy. Eventhoughsome researchers already applied a Marchenko method to 3D field data (Staringet al., 2018a; Pereira et al., 2018), they did not address the acquisition require-ments and limitations of 3D Marchenko methods in detail. The objective of thispaper is to gain a better understanding of the key acquisition parameters andlimitations that affect the application of the adaptive double-focusing methodto 3D data.In this paper, we first revise the theory of the adaptive Marchenko double-focusing method. Second, we perform a series of 3D synthetic tests to study theeffect of the acquisition parameters on the result of internal multiple predictionand adaptive subtraction using this method. Starting from a grid spacing of 50m (inline direction) by 75 m (crossline direction) co-located sources and receiverswith positive and negative offsets, near offsets, far offsets and a crossline apertureof 1.8 km (figure 2a), we step-by-step decimate the acquisition down to a narrowazimuth streamer geometry on which our 3D field data were acquired (figure 2b).Based on these tests, we identify the key limiting acquisition parameters and usethese to design an interpolation strategy for the field data application. Next, wetest the proposed interpolation strategy on 3D synthetic data. Finally, we applythe adaptive double-focusing method to 3D narrow azimuth streamer data. Inthe following discussion and conclusion section, we evaluate the performance ofthe adaptive double-focusing method. The adaptive Marchenko double-focusing method requires a preprocessed re-flection response R ( x R , x S , t ) acquired on a sufficiently dense grid of sources x S and receivers x R at the acquisition surface ∂ D . A smooth velocity model ofthe subsurface is needed to obtain the direct wave of the downgoing focusingfunction “ f +0 . The direct wave is obtained by modeling and time-reversing theresponse from sources at the redatuming level ∂ D i to receivers at the acquisitionsurface ∂ D using finite-difference modeling or an Eikonal solver (see figure 3).3he “ · symbol indicates an user-specified wavelet that is convolved with themodeled wavefield. The direct downgoing focusing function “ f +0 initiates the it-erative scheme that solves the coupled Marchenko equations. If the overburdenwere homogeneous, this initial wavefield would be sufficient to create a focusat the desired focal point at ∂ D i . Otherwise, a coda for the downgoing focus-ing function has to be retrieved using the following series (van der Neut et al.,2015a): “ f + ( x S , x (cid:48) F , t ) = ∞ (cid:88) i =0 “ f + i ( x S , x (cid:48) F , t ) = ∞ (cid:88) i =0 { Θ R (cid:63) Θ R} i “ f +0 ( x S , x (cid:48) F , t ) , (1)where i is the iteration number. Symbol x (cid:48) F denotes focal points at theredatuming level ∂ D i that become virtual sources. Operators R and R (cid:63) performa multidimensional convolution or correlation of the reflection response R withthe wavefield that it acts upon. Window functions Θ are tapered Heavisidestep functions that separate the causal and the acausal wavefields (i.e. Green’sfunctions and focusing functions) in time. See appendix A of Staring et al.(2018b) for details on the design of window function Θ . The first update of thecoda of the downgoing focusing function “ f +1 already contains many of the correctevents to compensate for the inhomogeneous overburden, but with incorrectamplitude. Higher-order estimates ( i = 2 , , , etc. ) are needed to obtain thecorrect amplitude.Using the downgoing focusing function “ f + , we can also retrieve the receiverredatumed upgoing Green’s function: “ G − ( x F , x S , t ) = Ψ R “ f + ( x S , x F , t )= ∞ (cid:88) j =0 “ G − j ( x F , x S , t ) = Ψ R ∞ (cid:88) j =0 { Θ R (cid:63) Θ R} j “ f +0 ( x S , x F , t ) , (2)where mute Ψ = I − Θ now selects the causal wavefield and symbol x F represents focal points at the redatuming level ∂ D i that become virtual receivers.The iteration number is given by j . Initial estimate “ G − is the standard receiver-redatumed upgoing Green’s function at x F . The first update “ G − contains a first-order estimate of the receiver-side internal multiples generated in the overburdenwith incorrect amplitude. Next updates ( “ G − , “ G − , etc.) contain higher-orderestimates that are necessary to obtain the correct amplitude. An additionalstep is needed to also remove source-side and source-and-receiver-side internalmultiples generated by the overburden.The retrieval of the upgoing Green’s function “ G − with a grid of sources atthe acquisition surface ∂ D and a grid of virtual receivers at the redatuming level ∂ D i is a single-focusing step. By creating double-focusing we also remove otherinternal multiples generated by the overburden. To this end, we convolve theupgoing Green’s function “ G − at virtual receivers with the downgoing focusing4unction “ f + at virtual sources (Wapenaar et al., 2016; van der Neut et al., 2018;Staring et al., 2018b): ““ G − + ( x F , x (cid:48) F , t ) = (cid:90) ∂ D “ G − ( x F , x S , t ) ∗ “ f + ( x S , x (cid:48) F , t ) d x S . (3)By applying this for many positions x F and x (cid:48) F at redatuming level ∂ D i , agrid of downward radiating virtual sources and virtual receivers that measurethe upgoing wavefield is created. The result is a redatumed Green’s function ““ G − + in the physical medium. Internal multiples generated by the overburden(figure 4a) have been removed, but later arriving internal multiples generatedby interactions between the target area and the overburden (figure 4b) andinternal multiples generated by the target area (figure 4c) remain. Accordingto Cypriano et al. (2015), the main internal multiples that contaminate theimage of the target area in the Santos Basin are generated between the waterbottom and the top of salt (see figure 1). By using double-focusing, we removethese internal multiples, while leaving some internal multiples below the targetarea behind. Note that one user-specified wavelet “ · has to be deconvolvedfrom the redatumed Green’s function ““ G − + . Also note that the integral over theacquisition surface ∂ D allows us to parallelize the implementation of the double-focusing method per pair of focal points, which makes this method particularlysuitable for the application to large 3D data volumes.Using equations 1 and 2, we can write equation 3 as a series: ““ G − + ( x F , x (cid:48) F , t ) = ∞ (cid:88) j =0 ∞ (cid:88) i =0 (cid:90) ∂ D “ G − j ( x F , x S , t ) ∗ “ f + i ( x S , x (cid:48) F , t ) d x S ≈ (cid:90) ∂ D “ G − ( x F , x S , t ) ∗ “ f +0 ( x S , x (cid:48) F , t ) d x S − (cid:90) ∂ D − “ G − ( x F , x S , t ) ∗ “ f +0 ( x S , x (cid:48) F , t ) d x S − (cid:90) ∂ D − “ G − ( x F , x S , t ) ∗ “ f +1 ( x S , x (cid:48) F , t ) d x S − (cid:90) ∂ D − “ G − ( x F , x S , t ) ∗ “ f +1 ( x S , x (cid:48) F , t ) d x S − ... (4)The first term “ G − ∗ “ f +0 is the standard source and receiver redatumed Green’sfunction including primaries and internal multiples. The second term − “ G − ∗ “ f +0 contains first-order predictions of receiver-side internal multiples generated bythe overburden, while the third term − “ G − ∗ “ f +1 contains first-order predictions ofsource-side internal multiples generated by the overburden and the fourth term − “ G − ∗ “ f +1 contains first-order predictions of source-and-receiver-side internal5ultiples generated by the overburden. Subsequent terms contain higher-orderestimates of the predicted internal multiples that are needed to obtain the cor-rect amplitude. Note that Marchenko methods in principle do not rely on anadaptive filter to accurately attenuate internal multiples. However, instead ofretrieving all terms in the series in equation 4 by correlating and convolvingthe data with itself many times (see equation 2), we propose to retrieve only afew updates and use an adaptive filter as a substitute for higher-order ampli-tude corrections. Also, an adaptive filter can correct for any minor amplitudeand phase differences that are present in the internal multiple predictions dueto imperfections in the data acquisition or preprocessing. Iteration numbers i and j that are needed to obtain predictions of all overburden internal multi-ples depend on the geological setting. In our case, new internal multiples werenot predicted beyond the third term in equation 4, so we only use the terms − “ G − ∗ “ f +0 and − “ G − ∗ “ f +1 for the prediction of internal multiples in this par-ticular setting. These predictions are treated as individual internal multiplepredictions, which are orthogonalized to the data prior to simultaneous adap-tive subtraction. We have chosen for an adaptive filter in the curvelet domain(Herrmann et al., 2008; Wu and Hung, 2015), since it can distinguish betweenprimaries and internal multiples in space, time and dip. Naturally, care has tobe taken not to subtract the primary reflections with the internal multiples. We perform a series of 3D synthetic tests to identify the key acquisition param-eters that affect the result of the adaptive double-focusing method. In order togenerate synthetic data that represent the geological contrasts in the area asrealistically as possible, we use a velocity model (see a 2D slice in figure 1) anda density model that are obtained from an acoustic inversion of field data basedon the original seismic image and migration velocity. The grid size of thesemodels is 18.75 m by 18.75 m by 10 m. Co-located sources and receivers arepositioned with a spacing of 50 m in the inline direction and a spacing of 75 min the crossline direction (figure 2a), thereby simulating an inline spacing of 50m and a sail line spacing of 75 m. The inline aperture is 20 km (offsets from -10km to 10 km) and the crossline aperture is 1.8 km. An acoustic finite-differencealgorithm is used to model data up to 30 Hz, such that the dominant wavelengthat the receivers is 50 m. The recording time is 8.5 s. Also, we generate an ini-tial focusing function “ f +0 in the smooth velocity model using an Eikonal solver.Geometrical spreading is part of the simulation. In addition, we convolve theresponse with an Ormsby wavelet with tapers at the low and the high ends.Starting from 24 lines of data modeled on this dense acquisition grid, westep-by-step decimate down to a realistic streamer acquisition geometry with acable spacing of 150 m, a sail line spacing of 450 m and a cable length of 6 km.The inline source and receiver spacing remain 50 m. Inline offsets range from250 m to 6250 m and the crossline aperture is 0.75 km (figure 2b). Throughoutthe decimation tests, we use the Marchenko double-focusing method to redatum6o a grid of co-located virtual sources and virtual receivers below the overburdenwith a spacing of 25 m by 37.5 m. We have chosen the redatuming level to bejust above the base of salt. The base of salt is the top of our reservoir and istherefore part of the target area. The main internal multiple generators in thisgeological setting, the water bottom and the top of salt, are part of the over-burden. Internal multiple predictions are obtained by convolving the individualupdates of the wavefields “ G − j and “ f + i , which are subsequently subtracted fromthe data using a 3D curvelet filter (Wu and Hung, 2015). We orthogonalize thepredictions and the data before subtraction, but do not use a global least squaresfilter for pre-conditioning. Parameters that need to be set are the number ofscales in the transform, the number of angles in the transform, the window sizeand some sparsity parameters that control the inversion. We extensively testdifferent filter settings and obtain the best results (the least damage of the pri-mary reflections) using 7 scales, 8 angles and tapered windows of 768 ms by 256traces. These settings are used for all synthetic examples shown here. First, we apply the adaptive double-focusing method to synthetic data generatedon the dense grid in figure 2a. Figure 5 shows redatumed common source gathersbefore and after internal multiple prediction and subtraction ( “ G − ∗ “ f +0 and ““ G − + from equation 4). The common source gathers are from a virtual source inthe middle of the grid of focal points, as indicated by the red star in figure2a. A difference is visible, especially in the yellow ellipses and along the yellowlines. It seems that conflicting seismic events were resolved, resulting in a bettercontinuity of the primary events.Next, we deconvolve an user-specified wavelet “ · and migrate the result.Figure 6 shows RTM images of the reflection response “ R at the acquisitionsurface (note that this image was truncated at the base of salt for comparison),the redatumed Green’s function including primaries and internal multiples G − ∗ “ f +0 and the redatumed Green’s function “ G − + after internal multiple predictionand subtraction. Figures 6a and 6b are comparable, thereby demonstrating thatsource-receiver redatuming was correctly performed (according to the standardprimary approach). A comparison between figure 6b and figure 6c shows adistinct difference. Internal multiples indicated by the yellow curved stripes infigures 6a and 6b are no longer visible in figure 6c. Also, the yellow ellipsesindicate areas where the removal of internal multiples is clearly visible. Overall,the continuity of the reflectors has improved. Below the vertical RTM imagesare depth slices of the 3D RTM volume at 5900 m depth, where the internalmultiples are present in figures 6a and 6b, but have been attenuated in figure 6c.Based on these results, we conclude that adaptive double-focusing performs wellin terms of redatuming and predicting and subtracting internal multiples whenapplying it to the initial dense acquisition geometry. Note that the image infigure 6c appears to have a lower frequency content compared to the images infigures 6a and 6b. It seems that the removed internal multiple reflections tend to7ave a high frequency, possibly due to the generation mechanism in the stratifiedsalt. The resulting images in other papers on internal multiple attenuation inthe Santos Basin (e.g. Griffiths et al., 2011; Cypriano et al., 2015) confirm thisobservation.In the following, we continue with synthetic data without negative offsets anduse source-receiver reciprocity for reconstruction before applying the adaptivedouble-focusing method. Since a sail line spacing of 75 m is not realistic, we study the effect of coarsersail line spacings on the result of our adaptive double-focusing method. Startingfrom the result in figure 6c with 75 m sail line spacing (here figure 7a), wecompare RTM images showing the result of adaptive double-focusing when usinga sail line spacing of 150 m (figure 7b), 300 m (figure 7c) and finally 450 m (figure7d). The result obtained from data with a sail line spacing of 150 m looksvery similar to the result obtained with 75 m sail line spacing, there are onlysome minor amplitude differences indicated by the arrows. A more significantdifference becomes visible when decimating from 150 m sail line spacing to 300m sail line spacing. Some internal multiples at the top of the image are no longerpredicted and subtracted, probably because the traces that are necessary for thereconstruction of these multiples are missing. The realistic scenario of 450 msail line spacing shows more internal multiples that could not be predicted andsubtracted, now in the deeper part of the image as well. The depth slices confirmthese observations: there is little difference between the results obtained with75 m and 150 m sail line spacing, but internal multiple attenuation becomes lesseffective when moving to a sail line spacing of 300 m and 450 m. Based on thesetests, we conclude that the sail line spacing is a key acquisition parameter thataffects our adaptive double-focusing method. Ideally, interpolation from 450 msail line spacing to 150 m sail line spacing would be applied prior to the fielddata application in this geological setting. We remark that although we expectthat the sail line spacing will also be a key acquisition parameter that affectsthe result of our adaptive double-focusing methods in other geological settings,the exact spacing at which the result is still acceptable will be different in everysetting. In the following synthetic tests, we continue with a sail line spacing of150 m, thereby assuming that the interpolation from 450 m sail line spacing to150 m sail line spacing can be carried out correctly.
The responses at near offsets are typically not recorded by streamers, so westudy the effect of removing the first 250 m of inline offsets. Figure 8 shows acomparison of RTM images with and without near offsets, both after internalmultiple prediction and subtraction. The removal of the near offsets deterioratesthe result somewhat in terms of a few remnant internal multiples (at the ellipseand at the arrows), but not as much as expected. The depth of the first reflector8nfluences how much the near offset responses contribute to the image. Sincethe water bottom in this setting is very deep (see figure 1), most reflectionsoriginating from this depth would have simply not been recorded by the first250 m of receivers. Even though the near offset responses do not have a largeeffect on the result of internal multiple prediction and removal in this very deepmarine setting, we will interpolate the field data for the missing offsets in orderto predict as many internal multiples as accurately as possible.
Next, we assume that we could correctly reconstruct the responses at near offsetsand we study the effect of removing the far offsets (the inline offsets 6250-10000m). Figure 9 shows the RTM images of the result of internal multiple predictionand subtraction using adaptive double-focusing, where figure 9a shows the resultwhen including far offsets in the reflection response and figure 9b shows the resultwhen excluding far offsets from the reflection response. Only minor differencesare visible in the result, mostly in terms of amplitudes. Surprisingly, the faroffsets seem to have little impact on the result of adaptive double-focusing,similar to the near offsets. Verschuur (2013) reports that missing offsets havea particularly large effect on multiple prediction methods in a shallow watersetting. Since we are in a very deep marine setting, missing offsets seem to onlyhave a minor effect.
Lastly, we study the effect of removing the outer cables. Instead of a crosslineaperture of 1800 m, as used in the previous tests, we now use a crossline apertureof 750 m. The RTM images in figure 10 show that removing the outer cableshas a significant effect on the adaptive double-focusing result. The quality ofthe image in figure 10b has deteriorated and some internal multiples were notpredicted and subtracted. Although the missing outer cables have a large effecton the result of adaptive double-focusing, the image in figure 10b is still ofacceptable quality. This becomes especially clear when comparing it to thestandard redatumed Green’s function in figure 6b, which is constructed from adense and wide azimuth grid of sources and receivers at the acquisition surface.Compared to this image, figure 10b still shows a significant reduction in internalmultiple energy. This is promising for the field data application, since we cannotcompensate for missing outer cables during preprocessing. We remark thatthe effect of removing the outer cables is expected to become more severe ingeological settings with strongly dipping reflectors in the crossline direction.In those cases, the missing outer cables can be a limiting factor that hindersthe application of the adaptive double-focusing method. This observation issupported by reports on the performance of similar multiple prediction andremoval methods (Wang and Hung, 2014; Moore and Dragoset, 2008).9 .6 The combination of all effects
Although the results of the synthetic tests in the previous sections are encour-aging, they do not provide an indication on the feasibility of the application ofthe adaptive double-focusing method to our field dataset, where all acquisitionrestrictions are imposed simultaneously. The negative offsets, the near offsets,the far offsets, some sail lines and the outer cables are all missing. Therefore,we model 32 lines of 3D synthetic data based on the acquisition geometry of ournarrow azimuth streamer data (figure 2b). Next, we reconstruct the negativeoffsets (by applying source-receiver reciprocity) and the near offsets (by inter-polation) and perform interpolation for the sail line spacing (from 450 m to 150m). Figure 11 shows the RTM images of the reflection response “ R , the stan-dard redatumed Green’s function G − ∗ “ f +0 with primaries and internal multiplesand the redatumed Green’s function “ G − + after internal multiple prediction andsubtraction, zoomed in at the target area. We observe an unexpected differ-ence in illumination when comparing figure 11a and figure 11b, especially onthe right side of the images. Figure 11a is constructed by applying an RTMmethod to the reflection response “ R , which uses a finite-difference method toback-propagate the wavefield from the acquisition surface to the redatuminglevel. In contrast, figure 11b is constructed by first back-propagating the reflec-tion response R using convolutions with the modeled direct downgoing focusingfunction, according to “ f +0 ∗ Ψ R ∗ “ f +0 (see equation 3), to obtain “ G − ∗ “ f +0 , whichis subsequently back-propagated from the redatuming level into the target us-ing the RTM method. In principle, back-propagation using a multidimensionalconvolution is equivalent to back-propagation using an RTM method (Esmersoyand Oristaglio, 1988). However, in practice, these are only equivalent when thesame numerical method is used. We use an Eikonal solver to model the directwave “ f +0 , which is different from the finite-difference method used in the RTMmethod. As a result, there are slight differences in illumination between the twoimages.Next, we evaluate the effectiveness of internal multiple attenuation. A com-parison of figures 11b and 11c shows that the adaptive double-focusing methodsucceeded in predicting and subtracting internal multiples from the standardredatumed Green’s function. Especially inside the yellow ellipses, the internalmultiple energy is significantly reduced, resulting in a better continuity of thereflectors. Again, we observe that the internal multiples seem to mainly have ahigh frequency content. We conclude that the adaptive double-focusing methodappears to be sufficiently robust for the prediction and adaptive subtraction ofinternal multiples from narrow azimuth streamer data in this geological setting. Based on the results of the synthetic tests, we continue with the field dataapplication. We have 24 lines of narrow azimuth streamer data, acquired using6 flat streamers with a cable spacing of 150 m and a sail line spacing of 45010. The length of the cables is 6000 m, covering offsets from 250 to 6250 m.The source and receiver spacing is 50 m in the inline direction. The crosslineaperture is 750 m. Deghosting is performed in the f-p domain (Wang et al., 2013)and a designature filter is obtained from the water bottom reflection. Shotand near offset reconstruction are performed using a partial normal moveout(NMO) correction of traces per common depth point (CDP) (e.g. Dragosetet al., 2010). In addition, the data are projected on a regular grid using a τ − p transform (Wang and Nimsaila, 2014). We also remove surface-relatedmultiples, the evanescent wavefield and noise. After preprocessing, we obtain adataset with a sail line spacing of 150 m. Similar to the synthetic tests, we havechosen the redatuming level just above the base of salt. A smoothed version ofthe velocity model in figure 1 and an Eikonal solver are used to model our directdowngoing focusing function “ f +0 (including geometerical spreading), which wesubsequently convolve with a 30 Hz Ormsby wavelet. Convergence is trackedby computing the L norm of the updates of the downgoing focusing function “ f + i .After convolving the individual updates of the wavefields “ f + i and G − j andmigrating them, we obtain the internal multiple predictions in the image domain.Extensive testing shows us that primary reflections are better preserved whensubtracting the internal multiple predictions in the image domain instead of inthe redatumed domain. Prior to subtraction, the predictions are othogonalized.We use the full curvelet transform (for all scales) and tapered windows of 768ms by 256 traces.Figures 12 and 13 show the result of applying the adaptive double focusingmethod to predict and adaptively subtract internal multiples. First, we observethat there is still a slight illumination difference between the RTM migratedimage of the reflection response in figure 12a and the image of the RTM mi-grated redatumed reflection response in figure 12b (see the yellow circle in thetop right), which has been explained above. However, the difference is not aspronounced as in figure 11. Second, a clear difference is visible between fig-ures 12b and figure 12c, indicating that the adaptive double-focusing methodsucceeded in predicting and adaptively subtracting events which are likely in-ternal multiples. The arrow, the ellipses and the lower yellow circle indicateareas in which events were attenuated, while the yellow boxes indicate whatwe interpret to be an improvement in fault definition. The white circle showsconflicting events that are being resolved, but it also shows events that we be-lieve to be remnant internal multiples. These were most likely generated bythe base of salt, which is not part of our overburden. In order to also removethese events and to see if there is perhaps a hidden structure, the redatuminglevel could be placed just below the base of salt. Similar to what we observedin the synthetic data, the image after internal multiple attenuation seems tohave a lower frequency content. The RTM depth slices in figure 13 demonstratethe attenuation of events in the form of ellipses. This example shows that theadaptive double-focusing method is sufficiently robust for the application to our3D narrow azimuth streamer dataset. Naturally, this result can be improved by11he use of more lines of data. In this paper, we identified the key acquisition parameters that affect the ap-plication of our adaptive Marchenko internal multiple attenuation method tonarrow azimuth streamer data. Tests on 3D synthetic data evaluated the effectof removing sail lines, near offsets, far offsets and the outer cables. The resultsof these tests show that the aperture in the crossline direction and the sail linespacing have the strongest effect on the quality of the result. Typically, thesail line spacing can be interpolated, but the aperture in the crossline directioncan possibly be a limiting factor for our method. Surprisingly, the missing nearoffsets and the far offsets only had a modest effect on the result of our method,possibly due to the very deep target area. In addition, we found that theresponses at the negative offsets and the near offsets could be accurately recon-structed. We remark that these tests are only valid for this particular dataset,but they give an impression of the possibilities and limitations of the adaptiveMarchenko double-focusing method. For an Ocean Bottom Node (OBN) ac-quisition geometry, these tests imply that our method will be most sensitiveto the node separation (especially in the direction of the strongest geologicalvariation).Based on the decimation tests, we defined an interpolation strategy thatwas first tested on a realistic synthetic dataset. We reconstructed the negativeoffsets and the near offsets, and interpolated the sail line spacing from 450 m to150 m. When applying the 3D adaptive double-focusing method, an importantaspect that was not visible in earlier 2D applications became visible, therebyshowing that the extension of a method from 2D to 3D is not always trivial.In 3D, when using an Eikonal solver for the modeling of the direct wave and afinite-difference-based RTM method, a slight difference in illumination betweenthe RTM image of the reflection response and the RTM image of the redatumedresponse occurs. Nevertheless, the double-focusing method predicted and adap-tively subtracted internal multiples, thereby improving the image of the targetarea.Next, we applied the adaptive double-focusing method to 24 lines of narrowazimuth streamer data. We reconstructed the negative and the near offsets andinterpolated the sail line spacing. We interpret that internal multiples werepredicted and adaptively subtracted, which resulted in an improved geologicalinterpretation of the target area. Therefore, we conclude that 3D Marchenkointernal multiple attenuation using an adaptive double-focusing method is suf-ficiently robust for the application to narrow azimuth streamer data in a deepmarine setting, provided that there is sufficient aperture in the crossline direc-tion and that the sail lines are interpolated.Note that redatuming is optional for Marchenko methods. The adaptivedouble-focusing method used in this paper includes source-receiver redatuming,which is particularly useful when the aim is to attenuate internal multiples in the12arget area and at the same time reducing the data volume for a next processingstep (for example, a target-oriented full waveform inversion). In contrast, whenthe aim is only to attenuate internal multiples as part of a larger workflow,the adaptive double-focusing method might not be the Marchenko method ofchoice. A direct quality check of the input common source gathers and theredatumed common source gathers is not possible, which is a disadvantage ina general processing workflow. In addition, a quality check on the resultingimages is only possible when the same numerical method is used to obtain thedirect wave for the Marchenko method as for the migration of the original data.Therefore, for the purpose of internal multiple elimination only, we propose theuse of other Marchenko methods that do not include redatuming and thus allowfor an easier quality check. An example is the adaptive overburden eliminationmethod (van der Neut and Wapenaar, 2016), as shown in papers by Pereiraet al. (2018), Krueger et al. (2018) and Pereira et al. (2019). A modified versionof the adaptive double-focusing method is presented by Staring et al. (2019).Other alternatives are the Marchenko multiple elimination scheme (Zhang andStaring, 2018) and the primary-only method proposed by Meles et al. (2016).The Marchenko method used in this paper is acoustic. The synthetic dataare acoustic, but naturally the field data are elastic. A suggestion for furtherresearch is to evaluate the effect of the presence of mode conversions on theacoustic Marchenko method in this geological setting, as was done by Reinickeet al. (2019) for the offshore Middle East. By applying an acoustic Marchenkomethod to elastic synthetic data, and comparing it to the result of applying anacoustic Marchenko method to acoustic synthetic data, Reinicke et al. (2019)evaluated whether the acoustic approximation is valid for structural imaging inthe region. They concluded that the acoustic approximation may be sufficientwhen used for structural imaging in 1.5D geological settings. acknowledgements
This research was performed in the framework of the project ’Marchenko imag-ing and monitoring of geophysical reflection data’, which is part of the DutchOpen Technology Programme with project number 13939 and is financially sup-ported by NWO Domain Applied and Engineering Sciences. The research of K.Wapenaar has received funding from the European Research Council (ERC)under the European Union’s Horizon 2020 research and innovation programme(grant agreement no. 742703). We thank Roberto Pereira, Peter Mesdag andAdel Khalil from CGG for the close collaboration and for providing the fielddata. Also, we would like to thank two anonymous reviewers, Lele Zhang, Gio-vanni Meles and Jan Thorbecke for valuable discussions.
Data Availability Statement
Data sharing is not applicable to this article.13 eferences
Broggini, F., Snieder, R. and Wapenaar, K. (2012) Focusing the wavefield insidean unknown 1D medium: Beyond seismic interferometry.
Geophysics , ,A25–A28.Cypriano, L., Marpeau, F., Brasil, R., Welter, G., Prigent, H., Douma, H.,Velasques, M., Boechat, J., de Carvalho, P., Guerra, C. et al. (2015) Theimpact of interbed multiple attenuation on the imaging of pre-salt targets inthe Santos basin off-shore Brazil. In . European Association of Geoscientists & Engineers.Dragoset, B., Verschuur, E., Moore, I. and Bisley, R. (2010) A perspective on3D surface-related multiple elimination. Geophysics , , 75A245–75A261.Esmersoy, C. and Oristaglio, M. (1988) Reverse-time wave-field extrapolation,imaging, and inversion. Geophysics , , 920–931.Foster, D. J. and Mosher, C. C. (1992) Suppression of multiple reflections usingthe Radon transform. Geophysics , , 386–395.Griffiths, M., Hembd, J. and Prigent, H. (2011) Applications of interbed multipleattenuation. The Leading Edge , , 906–912.Hampson, D. (1986) Inverse velocity stacking for multiple elimination. In SEGTechnical Program Expanded Abstracts 1986 , 422–424. Society of ExplorationGeophysicists.Herrmann, F. J., Wang, D. and Verschuur, D. J. (2008) Adaptive curvelet-domain primary-multiple separation.
Geophysics , , A17–A21.Jakubowicz, H. (1998) Wave equation prediction and removal of interbed multi-ples. In SEG Technical Program Expanded Abstracts 1998 , 1527–1530. Societyof Exploration Geophysicists.Krueger, J., Donno, D., Pereira, R., Mondini, D., Souza, A., Espinoza, J. andKhalil, A. (2018) Internal multiple attenuation for four presalt fields in theSantos Basin, Brazil. In
SEG Technical Program Expanded Abstracts 2018 ,4523–4527. Society of Exploration Geophysicists.Meles, G. A., Wapenaar, K. and Curtis, A. (2016) Reconstructing the primaryreflections in seismic data by Marchenko redatuming and convolutional inter-ferometry.
Geophysics , , Q15–Q26.Moore, I. and Dragoset, B. (2008) General surface multiple prediction: A flexible3D SRME algorithm. First Break , .van der Neut, J., Brackenhoff, J., Staring, M., Zhang, L., de Ridder, S., Slob, E.and Wapenaar, K. (2018) Single-and double-sided Marchenko imaging condi-tions in acoustic media. IEEE Transactions on Computational Imaging , ,160–171. 14an der Neut, J., Vasconcelos, I. and Wapenaar, K. (2015a) On Green’s func-tion retrieval by iterative substitution of the coupled Marchenko equations. Geophysical Journal International , , 792–813.van der Neut, J. and Wapenaar, K. (2016) Adaptive overburden eliminationwith the multidimensional Marchenko equation. Geophysics , , T265–T284.van der Neut, J., Wapenaar, K., Thorbecke, J., Slob, E. and Vasconcelos, I.(2015b) An illustration of adaptive Marchenko imaging. The Leading Edge , , 818–822.Pereira, R., Mondini, D. and Donno, D. (2018) Efficient 3D Internal MultipleAttenuation in the Santos Basin. In . European Association of Geoscientists & Engineers.Pereira, R., Ramzy, M., Griscenco, P., Huard, B., Huang, H., Cypriano, L. andKhalil, A. (2019) Internal multiple attenuation for OBN data with overbur-den/target separation. In SEG Technical Program Expanded Abstracts 2019 ,4520–4524. Society of Exploration Geophysicists.Reinicke, C., Dukalski, M. and Wapenaar, K. (2019) Do we need elastic in-ternal de-multiple offshore Middle East, or will acoustic Marchenko suffice?In
SEG/KOC Workshop: Seismic Multiples - The Challenges and the WayForward, Kuwait . Society of Exploration Geophysicists.Staring, M., Dukalski, M., Belonosov, M., Baardman, R., Yoo, J., Hegge, R.,van Borselen, R. and Wapenaar, K. (2019) Adaptive Marchenko multipleremoval on Arabian Gulf OBC data. In
SEG/KOC Workshop: Seismic Mul-tiples - The Challenges and the Way Forward, Kuwait . Society of ExplorationGeophysicists.Staring, M., van der Neut, J. and Wapenaar, K. (2018a) Marchenko redatumingby adaptive double-focusing on 2D and 3D field data of the Santos basin.In
SEG Technical Program Expanded Abstracts 2018 , 5449–5453. Society ofExploration Geophysicists.Staring, M., Pereira, R., Douma, H., van der Neut, J. and Wapenaar, K. (2018b)Source-receiver Marchenko redatuming on field data using an adaptive double-focusing method.
Geophysics , , S579–S590.Verschuur, D. J. (2013) Seismic multiple removal techniques: past, present andfuture . EAGE publications.Wang, M. and Hung, B. (2014) 3D Inverse Scattering Series Method for Inter-nal Multiple Attenuation. In .European Association of Geoscientists & Engineers.Wang, P. and Nimsaila, K. (2014) Fast progressive sparse Tau-P transform forregularization of spatially aliased seismic data. In
SEG Technical ProgramExpanded Abstracts 2014 , 3589–3593. Society of Exploration Geophysicists.15ang, P., Ray, S., Peng, C., Li, Y. and Poole, G. (2013) Premigration deghost-ing for marine streamer data using a bootstrap approach in Tau-P domain.In
SEG Technical Program Expanded Abstracts 2013 , 4221–4225. Society ofExploration Geophysicists.Wapenaar, K., van der Neut, J. and Slob, E. (2016) Unified double- and single-sided homogeneous Green’s function representations. In
Proc. R. Soc. A , vol.472, 20160162. The Royal Society.Wapenaar, K., Thorbecke, J., Van Der Neut, J., Broggini, F., Slob, E. andSnieder, R. (2014) Marchenko imaging.
Geophysics , , WA39–WA57.Ware, J. A. and Aki, K. (1969) Continuous and discrete inverse-scattering prob-lems in a stratified elastic medium. I. Plane waves at normal incidence. Thejournal of the Acoustical Society of America , , 911–921.Weglein, A. B., Gasparotto, F. A., Carvalho, P. M. and Stolt, R. H. (1997) Aninverse-scattering series method for attenuating multiples in seismic reflectiondata. Geophysics , , 1975–1989.Wu, X. and Hung, B. (2015) High-fidelity adaptive curvelet domain primary-multiple separation. First Break , , 53–59.Zhang, L. and Staring, M. (2018) Marchenko scheme based internal multiplereflection elimination in acoustic wavefield. Journal of Applied Geophysics , , 429–433.Zhou, B. and Greenhalgh, S. A. (1994) Wave-equation extrapolation-based mul-tiple attenuation: 2-D filtering in the fk domain. Geophysics , , 1377–1391.16 ist of Figures “ f +0 , which we obtain using asmooth velocity model and an Eikonal solver. . . . . . . . . . . 214 Cartoons illustrating the redatumed sources at x (cid:48) F and the re-datumed receivers at x F resulting from the adaptive double-focusing method. Internal multiples a) generated by the overbur-den have been removed by double-focusing, while b) later arrivinginternal multiples generated by interactions between the targetarea and the overburden and c) later arriving internal multiplesgenerated by the target area remain. . . . . . . . . . . . . . . . 225 Redatumed common source gathers showing the result of theadaptive double-focusing method applied to 3D synthetic datamodeled on the dense acquisition grid in figure 2a. The gath-ers show: a) the redatumed Green’s function G − ∗ “ f +0 includingprimaries and internal multiples and b) the redatumed Green’sfunction “ G − + after prediction and adaptive subtraction of inter-nal multiples. The yellow ellipses, stripes and arrows indicateareas in which internal multiple attenuation is most visible. . . . 236 RTM images showing the result of the adaptive double-focusingmethod applied to 3D synthetic data modeled on the dense ac-quisition grid in figure 2a. The images show: a) the migratedreflection response “ R , b) the redatumed and migrated Green’sfunction G − ∗ “ f +0 including primaries and internal multiples, andc) the redatumed and migrated Green’s function “ G − + after pre-diction and adaptive subtraction of internal multiples. Below theRTM images are depth slices at 5900 m. The yellow ellipses in-dicate areas in which internal multiple attenuation is most visible. 247 RTM images obtained from the Green’s function “ G − + after inter-nal multiple prediction and subtraction for a reflection responsemodeled with: a) 75 m sail line spacing, b) 150 m sail line spac-ing, c) 300 m sail line spacing and d) 450 m sail line spacing.Below the RTM images are depth slices at 5900 m. The resultsusing 75 m and 150 m sail line spacing are similar, but the imagestarts to deteriorate when moving to 300 m sail line spacing. . . 2517 RTM images obtained from the Green’s function “ G − + after inter-nal multiple prediction and subtraction for a reflection responsemodeled with: a) near offsets (0-250 m) and b) without near off-sets. Below the RTM images are depth slices at 5900 m. Only aslight difference is visible at the arrows, possibly due to the depthof the target zone. . . . . . . . . . . . . . . . . . . . . . . . . . . 269 RTM images obtained from the Green’s function “ G − + after inter-nal multiple prediction and subtraction for a reflection responsemodeled with: a) far offsets (6250-10000 m) and b) without faroffsets. Below the RTM images are depth slices at 5900 m. Onlyminor differences are visible at the arrows, mostly in terms ofamplitude. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2710 RTM images obtained from the Green’s function “ G − + after inter-nal multiple prediction and subtraction for a reflection responsemodeled: a) with outer cables (1.8 km crossline aperture) andb) without outer cables (0.75 km crossline aperture). Below theRTM images are depth slices at 5900 m. The image considerablydeteriorates when removing the outer cables, as indicated by thecircles and arrows. . . . . . . . . . . . . . . . . . . . . . . . . . . 2811 RTM images of the result after applying the adaptive double-focusing method to 3D synthetics with a field data acquistiongeometry. The images show: a) the migrated reflection response “ R , b) the redatumed and migrated Green’s function G − ∗ “ f +0 including primaries and internal multiples, and c) the redatumedand migrated Green’s function “ G − + after prediction and adaptivesubtraction of internal multiples. Below the RTM images aredepth slices at 5900 m. . . . . . . . . . . . . . . . . . . . . . . . 2912 RTM images of the result after applying the adaptive double-focusing method to 3D NAZ streamer data of the Santos Basin,Brazil. The images show: a) the migrated reflection response “ R , b) the redatumed and migrated Green’s function G − ∗ “ f +0 in-cluding primaries and internal multiples, and c) the redatumedand migrated upgoing Green’s function “ G − + after prediction andadaptive subtraction of internal multiples. Below the RTM im-ages are depth slices at 5900 m. Internal multiples were predictedand subtracted, resulting in an improved image of the target area. 3013 Depth slices corresponding to the RTM images in figure 12 at5900 m. The slices show: a) the migrated reflection response “ R , b) the redatumed and migrated Green’s function G − ∗ “ f +0 including primaries and internal multiples, and c) the redatumedand migrated upgoing Green’s function “ G − + after prediction andadaptive subtraction of internal multiples. . . . . . . . . . . . . . 3118igure 1: A 2D slice of the velocity model of the Santos Basin, Brazil.19igure 2: Cartoons showing a) the starting acquisition geometry and b) thefinal acquisition geometry for the synthetic decimation tests in this paper. Thefinal acquisition geometry is based on our narrow azimuth streamer data. Thestars represent sources and the triangles represent receivers.20igure 3: Cartoon illustrating the direct wave “ f +0 , which we obtain using asmooth velocity model and an Eikonal solver.21igure 4: Cartoons illustrating the redatumed sources at x (cid:48) F and the redatumedreceivers at x F resulting from the adaptive double-focusing method. Internalmultiples a) generated by the overburden have been removed by double-focusing,while b) later arriving internal multiples generated by interactions between thetarget area and the overburden and c) later arriving internal multiples generatedby the target area remain. 22igure 5: Redatumed common source gathers showing the result of the adap-tive double-focusing method applied to 3D synthetic data modeled on the denseacquisition grid in figure 2a. The gathers show: a) the redatumed Green’s func-tion G − ∗ “ f +0 including primaries and internal multiples and b) the redatumedGreen’s function “ G − + after prediction and adaptive subtraction of internal mul-tiples. The yellow ellipses, stripes and arrows indicate areas in which internalmultiple attenuation is most visible. 23igure 6: RTM images showing the result of the adaptive double-focusingmethod applied to 3D synthetic data modeled on the dense acquisition gridin figure 2a. The images show: a) the migrated reflection response “ R , b) theredatumed and migrated Green’s function G − ∗ “ f +0 including primaries and in-ternal multiples, and c) the redatumed and migrated Green’s function “ G − + after prediction and adaptive subtraction of internal multiples. Below the RTMimages are depth slices at 5900 m. The yellow ellipses indicate areas in whichinternal multiple attenuation is most visible.24igure 7: RTM images obtained from the Green’s function “ G − + after internalmultiple prediction and subtraction for a reflection response modeled with: a)75 m sail line spacing, b) 150 m sail line spacing, c) 300 m sail line spacing andd) 450 m sail line spacing. Below the RTM images are depth slices at 5900 m.The results using 75 m and 150 m sail line spacing are similar, but the imagestarts to deteriorate when moving to 300 m sail line spacing.25igure 8: RTM images obtained from the Green’s function “ G − + after internalmultiple prediction and subtraction for a reflection response modeled with: a)near offsets (0-250 m) and b) without near offsets. Below the RTM images aredepth slices at 5900 m. Only a slight difference is visible at the arrows, possiblydue to the depth of the target zone. 26igure 9: RTM images obtained from the Green’s function “ G − + after internalmultiple prediction and subtraction for a reflection response modeled with: a)far offsets (6250-10000 m) and b) without far offsets. Below the RTM imagesare depth slices at 5900 m. Only minor differences are visible at the arrows,mostly in terms of amplitude. 27igure 10: RTM images obtained from the Green’s function “ G − + after internalmultiple prediction and subtraction for a reflection response modeled: a) withouter cables (1.8 km crossline aperture) and b) without outer cables (0.75 kmcrossline aperture). Below the RTM images are depth slices at 5900 m. Theimage considerably deteriorates when removing the outer cables, as indicatedby the circles and arrows. 28igure 11: RTM images of the result after applying the adaptive double-focusingmethod to 3D synthetics with a field data acquistion geometry. The imagesshow: a) the migrated reflection response “ R , b) the redatumed and migratedGreen’s function G − ∗ “ f +0 including primaries and internal multiples, and c) theredatumed and migrated Green’s function “ G − + after prediction and adaptivesubtraction of internal multiples. Below the RTM images are depth slices at5900 m. 29igure 12: RTM images of the result after applying the adaptive double-focusingmethod to 3D NAZ streamer data of the Santos Basin, Brazil. The imagesshow: a) the migrated reflection response “ R , b) the redatumed and migratedGreen’s function G − ∗ “ f +0 including primaries and internal multiples, and c) theredatumed and migrated upgoing Green’s function “ G − + after prediction andadaptive subtraction of internal multiples. Below the RTM images are depthslices at 5900 m. Internal multiples were predicted and subtracted, resulting inan improved image of the target area. 30igure 13: Depth slices corresponding to the RTM images in figure 12 at 5900m. The slices show: a) the migrated reflection response “ R , b) the redatumedand migrated Green’s function G − ∗ “ f +0 including primaries and internal multi-ples, and c) the redatumed and migrated upgoing Green’s function “ G − ++