Attenuation of receiver ghosts in variable-depth streamer high-resolution seismic reflection data

Receiver ghosts attenuate marine seismic reflection data at harmonic frequencies that depend on the propagation angle and the streamer depth below the sea surface. The resulting loss of bandwidth is one of the major factors hampering seismic resolution. In near-surface and legacy multi-channel data, receivers depth is often unknown and may vary significantly both along the streamer length and during the survey, making frequency–slowness deghosting techniques unsuitable. In this work, we present a method for the attenuation of receiver ghost reflections in data with an arbitrary streamer depth profile varying during the survey. For each trace, a different deghosting operator is estimated and applied at different two-way-time windows, in order to account for depth-dependent changes in reflection angles. The ghost null frequencies are picked on the time-varying power spectrum via an automatic algorithm, guided by a user-dependent a priori function, and optimised to respect the harmonics’ periodicity. The power of the inverse filter is adjusted by adaptively damping abnormal amplitudes in the deghosted spectra. The algorithm is applied to high-resolution (GI-gun, 20–400 Hz) and ultra-high-resolution (Sparker, 0.2–3.0 kHz) multi-channel datasets, yielding an excellent bandwidth recovery and gain in resolution in the final stacks. Limited computing time and straightforward application make the method widely applicable and cost-effective.


Introduction
Maximising the frequency bandwidth of recorded seismic reflection data is instrumental to obtaining optimal resolution and penetration depth within the inherent limitations of source and receiver antennas. One of the main factors hampering seismic resolution, therefore compromising the interpretability of the subsurface imaging, is the presence of ghost reflections in the recorded wavefield (e.g., Aytun 1999;Bai et al. 2013;Sun et al. 2015).
A ghost reflection originates at the sea-surface and interferes with the primary up-going wavefield (Fig. 1), producing notches in the recorded frequency spectrum at harmonics whose fundamental is equal to the reciprocal of the relative time delay (Aytun 1999;Yilmaz 2012). The resulting loss of bandwidth depends on the acquisition setting, namely, the deeper the receiver below the sea surface and the narrower the reflection angle, the longer the time delay and the lower the fundamental ghost notch frequency. If the time delay is longer than the wavelet duration, commonly in very-highfrequency seismic data, the receiver ghost may also appear in the time-offset domain as a separate arrival with opposite polarity (e.g., Provenzano et al. 2018).
The cartoon in Fig. 1 shows that up-going plane wavefronts reflected at the seafloor interfere with their ghosts at time delays depending both on the receiver depth and the source-receiver offset. Furthermore, for each source-receiver offset, reflections originating at interfaces deeper than the seafloor interfere with their ghosts at a longer time delay, by virtue of a narrower reflection angle, and hence their bandwidth is attenuated at harmonics of a lower fundamental frequency (Aytun 1999). Though not Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s1100 1-020-09407 -9) contains supplementary material, which is available to authorized users.
1 3 11 Page 2 of 15 shown in Fig. 1, wavefronts propagating upwards from the source can similarly be reflected off the sea-surface, thereby generating source-side ghost reflections (e.g., Kluesner et al. 2018). In particular, this paper focusses on the attenuation of the receiver-side ghosts when the receivers depth is unknown, non-homogeneous across the streamer and nonconstant during acquisition.
Carefully designed acquisition geometries can in principle minimise the loss of bandwidth determined by the presence of ghost reflections (e.g., . For example, slanted streamer geometries increase the post-stack bandwidth by diversifying the ghost notch frequencies over the offset range (Soubaras and Dowle 2010). Simultaneous over-under acquisitions combine two streamers towed at different depths, to fill in the respective ghost notches (Özdemir et al. 2008). Multi-component streamer, on the other hand, combine velocimeters and hydrophones (Carlson 2007) exploiting the different polarities of particle-velocity and pressure reflected at the sea-surface (Yilmaz 2012).
If receiver depth is uniform across the streamer length, single-component hydrophone data can be deghosted by means of an inverse filter in the plane wave domain (Aytun 1999), which can be modified to be applied to data acquired with slanted streamers (e.g., Masoomzadeh et al. 2013;Vrolijk and Blacquiere 2017). More complex streamer geometries, e.g. curved streamers with tail-buoys ( Fig. 1), make the relationship between horizontal slowness and ghost notch frequency channel-dependent. In such cases, ghost attenuation requires more complex frequency-domain strategies (e.g., Masoomzadeh et al. 2015;Zhang and Masoomzadeh 2018), wavefield extrapolation and inversion approaches (Wang et al. 2017;Sun and Verschuur 2018), or deghosting by joint deconvolution of mirror migrated images (Soubaras 2010;Beasley et al. 2013).
High-resolution scientific or engineering-scale seismic acquisitions, are often single component and with little control on the streamer position at depth (Pinson et al. 2008;Provenzano et al. 2017); partially because of budget limitations, but also since depth-controllers may introduce significant levels of noise in short streamer ( ≤ 200m ) highfrequency ( 10 2 −10 3 Hz ) data. This prevents the effective application of both acquisition-based and frequency-slowness methods.
Here we propose a novel time-frequency-offset domain approach for the attenuation of the receiver ghost in data with arbitrarily complex streamer depth profile, in the general case in which the latter is unknown and varies during the survey. The methodology has been developed using multi-channel high-resolution (GI-gun, 20-400 Hz) and ultra-high-resolution (Sparker, 0.2-3.0 kHz) 2D seismic reflection data. We will show that the proposed methodology increases the seismic resolution, improving the quality of velocity analysis and the interpretability of the processed seismic data.

Methodology
In seismic data acquired with a flat streamer towed at depth z r below the sea surface, the time delay between a primary reflection and its ghost can be written as: where is the reflection angle and cos( ) yields the vertical component ( p z ) of the slowness vector . In the frequency domain, the recorded seismogram resulting from the interference between the primary and the receiver ghost, will (1) t = 2z r cos( ), Fig. 1 Multi-channel 2D acquisition. Schematic view of the acquisition geometry for multichannel seismic reflection data with non-homogeneous receiver depth ( z i , vertically exaggerated for clarity). Wave-paths for two primaries (black solid line) and their ghosts (red solid line), and plane wavefronts (dashed lines, same colour code). Power spectra for receivers 1 (grey) and 2 (blue) as they result from the interference between primary and ghost reflection. Note that the period of the ghost harmonics is lower for receiver 1 as a result of the longer time delay between primary and ghost be attenuated at harmonics of f 0 = 1∕ t, according to the equation: where r ss is the sea surface reflection coefficient (see "Appendix 1" for a derivation of Eq. 2). From Eq. 2, an inverse ghost filter can be designed in the frequency-slowness domain, or equivalently frequency-wavenumber, if the receiver depth and sea surface reflection coefficients are known and homogeneous. However, these conditions are often violated both as a consequence of the buoyancy characteristics of the streamer, and of changes in navigation condition e.g. sea currents and swell (Masoomzadeh et al. ( , 2015) during data acquisition. Furthermore, the effective sea surface reflection coefficient may vary significantly due to the wave conditions (Blacquiere and Sertlek 2019).
In the latter case, Eq. 2 is complicated by the variation of z r , and r ss , as a function of the recorded channel. Therefore, a frequency-slowness operator is not sufficient to describe the interaction between primary and ghost reflections, as it would neglect that each plane wave identified by ( , , f ) is in fact recorded by each channel at a different depth. Furthermore, a frequency-offset approach would be inaccurate as it ignores the presence of different slowness components at each trace (e.g., Provenzano et al. 2017). The methodology presented in this paper has been developed to address this issue in 2D multi-channel seismic reflection datasets acquired with an arbitrary streamer depth profile, e.g. curved as a result of the use of a tail-buoy ( Fig. 1), which changes during acquisition, and particularly may vary by a significant fraction of a wavelength.
In order to account for trace-by-trace variations in receiver depth and, for each trace, changes in the recorded propagation angle at different two-way travel times, a timevarying deghosting operator is estimated in the short-time Fourier transform domain (Zhivomirov 2019). For each time window ( t w ) of each trace ( ), the ghost notch frequencies are picked automatically in the short-time power spectrum in the neighbourhood ( ±w ) of an initial estimate [g( )]. The latter can be picked at decimated shot locations on the amplitude spectrum of the data windowed around the seafloor reflection, and then smoothed and interpolated as a function of shot number and absolute offset. In the attempt to reduce the risk of mis-picks, the periodic nature of the ghost reflection is exploited, and a number (n) of harmonics of f 0 are picked at integer multiples of g( ).
The picks are then optimised by line-searching for the optimal fundamental frequency f opt , i.e. the f 0 predicting the higher order harmonics with the lowest mean-square error (4) Note that we make no assumptions about the local receiver depth, nor on the dominant slowness at each time-window, but rather obtain a data-driven estimate of the optimal fundamental ghost notch frequency. This process can be summarised as a time-varying spectral whitening procedure consisting in two main stages: -A first step of ghost notch picking on a coarse grid, to reduce the search range of the optimised automatic picking. -A time-varying trace-by-trace automatic picking that updates the initial estimate of f 0 to account for: (a) variations of the streamer depth during acquisition and (b) changes in the fundamental frequency at different twoway travel time in the data.
Finally, the equation for the inverse filter, for each time window in the recorded trace is: Here, the sea surface is assumed to be a nearly-perfect reflector, with r ss ≃ −1. In order to account for possible overestimation of r ss , the power of the inverse filter is adjusted by multiplying Eq. 5 for a damping function, computed as: where is a user-dependent constant (Aster et al. 2005

Estimation of streamer depth for datuming
The fundamental ghost-notch frequencies picked around the seafloor reflection are also valuable to estimate the receivers depth as a function of shot (s) and channel (c) number, for the purpose of accurate seismogram datuming. Equation 2 can be rewritten for n harmonics as: where V w is the p-wave velocity of sea-water, and can be approximated to the ratio of the seafloor reflection travel time at the 1st and cth channel ( t 1 ∕t c ) (e.g., Pinson et al. 2008;Pinson 2009). Equation 7 can be linearly inverted for the reciprocal of z r (Provenzano et al. 2017), or more conveniently, for a polynomial approximation of the streamer depth profile function of channel number (c) and source number (s): where k j,l are the polynomial coefficients solution of the linear inverse problem. The polynomial approximation has the advantage over a trace-by-trace inversion of reducing the number of unknowns to the polynomial order, therefore making the inverse problem for each shot well posed, even if only one harmonic per trace is used as input data (Menke 1989). Furthermore, by filtering out high-spatial

Application
The method has been developed using a 2D survey, employing a GI-gun source (20-400 Hz, two elements of 160 in. 3 each) and a 120-channel streamer with group spacing equal to 1.56 m and absolute offset ranging from 40 to 230 m. Shot point spacing is equal to approximately 12 m, and common-depth-point (CDP) fold is between 20 and 30 traces in 2 m bins. Ghost notch frequencies are picked every 25 shots and inverted for a 4th order two-dimensional polynomial function of channel and shot number (Eq. 8). In Fig. 2, example common shot gathers sampled at every 50 shots are plotted against the estimated streamer depth. Note that the latter is affected by the use of a tail-buoy (Fig. 1), and that the variation of streamer depth along the survey line is significant with respect to the propagated wavelengths for the frequency band of interest (4-75 m), both within each gather and along the survey line. Before automatic picking is performed in the neighbourhood of the interpolated coarse-picks, data preprocessing is performed in order to remove sources of noise Raw shots: three example GI-gun common shot gathers before deghosting; pre-processing consists in wave-equation datuming, bandpass and dip filtering; frequency-offset power spectra show at least three orders of ghost notches ( f 0 to f 2 ) that might affect the filter parameter estimation: minimum phase bandpass frequency, frequency-wavenumber filtering within a range of plausible apparent wavenumbers, and seafloor multiple attenuation by gapped deconvolution (Yilmaz 2012).
A time-window of 60 ms, sliding with 30 ms overlap and zero-padded to achieve a 1 Hz spectral resolution, is used for the short-time Fourier transform of each trace. The length of the window should be sufficient to include the primary and the ghost for each reflection, as well as short enough to capture the time-variation of the ghost notches. Figure 3 shows for a time-window at different offsets the robustness of the automatic picking procedure against inaccuracies of the guide function. As long as the local minimum of the power spectrum is within the search window (in this case, 40 Hz wide), the automatic picking (Eq. 3) is able to locate the fundamental ghost notch frequency and predict its higher order harmonics. The optimisation stage (Eq. 2) improves such prediction, producing cleaner offset-dependent nullfrequency curves. This is shown with greater detail in Fig. 4, where two power spectra computed on different two-way-time windows of the same trace are compared. The auto-picker repositions the initial guide function frequencies yielding a better prediction of the higher-order notches. It is noteworthy that such update is of different magnitude for the shallow and deep spectra, the latter being affected by notches at lower frequencies, as predicted by Eq. 1 for narrower angle reflections. In Fig. 5, the average spectra for a raw and deghosted shot gather are compared. The broader and flatter frequency band of the deghosted seismogram is apparent; although the two time window are only 100 ms apart, the time-varying approach has a positive impact on bandwidth recovery in the deeper part compared to a simple frequency-offset deghosting.
Finally, before transformation into the time-offset domain, data are filtered in frequency and dip using the same parameters as in the pre-processing stage, in order to attenuate possible deghosting artefacts. The impact of deghosting on the pre-stack dataset can be appreciated by comparing Figs. 6 and 7, both in time-offset and frequency-offset domains. In Fig. 7 the average power spectra overlays show a significant gain in bandwidth between 180 and 400 Hz. The bandwidth recovery, and the correspondent attenuation of ghost reflections in the time-offset domain, improves the quality of semblance velocity analysis on CDP super-gathers (Fig. 8), allowing for a more robust and repeatable stacking velocity estimation.
The combined effect of bandwidth recovery and more robust velocity analysis yields a higher quality stack final image (Fig. 9). Note how the low frequency ringing below the bright spots at CDP intervals 500-600, 1000-1500 and 2000-2300 has been eliminated, therefore improving the interpretability of the seismic image. The close-up in Fig. 10 highlights the increased sharpness of the deghosted stack image in a structurally complex interval, and shows how this corresponds to a broader post-stack bandwidth compared to the raw stack. In particular, the improved resolution allows us to interpret with greater confidence the geometry of the erosional unconformity between 0.3 and 0.4 s (CDPs 2900-3500), as well as providing a clearer image of small scale reflectivity anomalies corresponding gas accumulation pockets within thin stratigraphic horizons (0.46-0.55 s, CDPs 2900-3600).
The method has also been applied to an ultra-high-resolution sparker (0.2-3.0 kHz) survey acquired with a similar streamer configuration in the same area. The source is an Applied Acoustics Squid-2000 Sparker operating at 1750 J with a shot spacing of 4 m. The CDP fold ranges from 45 to 60 traces for a 1.5 m bin size. Pre-processing, in this case, included predictive deconvolution to attenuate the secondary pulse of the sparker signature at ca. 1 ms (e.g., Kluesner et al. 2018). In this dataset, due to the lack of frequencies lower than 0.2 kHz, in most channels the fundamental ghost null frequency was below the noise level, and the auto-picker started at the second order harmonic; however, the broader bandwidth allowed us to introduce up to four higher order harmonics in the optimisation. In Fig. 11, the effect of de-ghosting is apparent in the attenuation of the seafloor ghost reflection, that shadows the true seismic stratigraphy between 0.21 and 0.22 ms twtt. This corresponds to a gain in post-stack bandwidth between 1.0 and 3.0 kHz and a generally flatter frequency response.

Discussion
The deghosting method presented in this paper has been demonstrated to be effective in two high-resolution datasets with unknown receivers depths and complex streamer geometry, in which standard frequency-slowness techniques would not be applicable. The gain in seismic resolution Fig. 7 Deghosted shots: three example GI-gun common shot gathers (same as in Fig. 6) after deghosting; pre-processing consists in waveequation datuming, bandpass and dip filtering; note the bandwidth recovery in the frequency-offset power spectra and in the deghosted average power spectra overlaid to the raw ones ◂ produces a clearer velocity estimation, and a sharper, more interpretable final stack image.
Pre-stack bandwidth recovery is more apparent than post-stack (Figs. 7, 10), because stacking also attenuates ghosts in raw data by averaging traces acquired at different receiver depths (Soubaras and Dowle 2010). However, unlike in conventional slanted streamer acquisitions, in these datasets ghost cancellation by stacking is only partial, since the streamer curvature induced by the tail-buoy reduces the range of ghost notch frequencies across the offset range (Fig. 1).
In the sparker data example, the method effectively attenuates the seafloor ghost reflection (Fig. 11), though the bandwidth recovery is not as apparent as in the GI-gun case (Fig.10). This might be due to an imperfect attenuation of the sparker bubble pulse in the pre-processing stage, as a consequence of the non-repeatability of the primary-bubble delay hampering the effectiveness of predictive deconvolution; the associated frequency-domain notches constitute a source of error in the estimation of the receiver ghost fundamental frequency, and therefore the optimal deghosting operator.
The time-varying approach relies on the simple assumption that each two-way-time window is dominated by arrivals with a specific apparent slowness, and that this variation is captured by the automatic ghost notch picking on the Fig. 8 Impact of deghosting on velocity analysis: semblance velocity analysis of GI-gun data, before deghosting (panel a) shows broader coherency maxima and secondary peaks; after deghosting (panel b) stacking velocity estimation is more precise. Semblance values for the two panels are normalised to one short-time power spectrum. This requires a careful choice of the sliding window size for the short-time Fourier transform. In complex subsurfaces, this assumption is violated by the presence of multi-path arrivals, converted and turning waves; in this scenario, the method could be applied to migrated common reflection point gathers (CRPs) rather than to common shot gathers, in order to approximate a locally one-dimensional wave propagation.
Compared to more complex methodologies, applicable regardless the wavefield complexity, such as deghosting by joint deconvolution of mirror migrated gathers (Soubaras 2010), the proposed approach has the advantage of a low computational cost, since it requires only Fourier-domain operations. Furthermore, the shot-by-shot analysis can easily be parallelised, with dramatic scope for improvement of the computing cost.
The requirements in terms of user-time are limited to the design of a guide-function for the auto-picker. In the examples presented here, a multi-channel picking every 25/30 shots was sufficient. The optimal grid spacing depends on the specific dataset, how variable the navigation conditions were during acquisition and, therefore, how frequently the streamer geometry changes shot-by-shot. In principle, the longer the shot spacing and the rougher the sea conditions, the denser the guide function picking should be. However, as shown in "Application" section, the automatic picker is designed to be robust against inaccurate or noisy guide functions. A careful tuning of the damping factor in Eq. 6 is necessary in order to attenuate artefacts arising from inaccurate values of the sea-surface reflection coefficients ( r ss ). In this work, this is done by quality-checking the results at few selected shot gathers. Iterative estimation of an optimal value for r ss as in (Provenzano et al. 2017) would improve the robustness of the method, while increasing sensibly its computing cost.
Though the methodology has been developed for 2D streamer acquisitions, the impact of lateral streamer movements on the ghost notch frequencies can in principle be captured and corrected for by the trace-by-trace estimated deghosting operators, thereby extending the applicability to datasets with strong streamer feathering. However, it is sensible to assume that vertical movements of streamer sections have a dominant effect on the relative time-delay between primary and ghost reflections.
Under this assumption, as a by-product of the estimation of the optimal deghosting operator, an accurate streamer depth profile can be obtained and used for elevation corrections or wave-equation datuming. It is recommended that only the slowly varying component of the estimated depths are used (e.g. by inverting for a polynomial approximation), so that high-frequency variations of the ghost notch frequencies, due to local changes in sea-surface conditions, are isolated from variations of streamer geometry during acquisition (Masoomzadeh et al. 2015).

Conclusions
In this paper, we have presented a novel deghosting method developed for 2D multichannel data with unknown receiver depths and varying streamer depth profile. The main advantages of this technique are: 1. It requires simple pre-processing, limited to bandpass and dip filtering, 2. It does not make assumptions about the apparent slowness and local receiver depth, 3. It automatically estimates the deghosting filters parameters with little user input, 4. It is applicable to datasets where frequency-slowness/ wavenumber deghosting is unsuitable due to non-constant, non-homogeneous receiver depths, 5. It produces as a by-product an accurate receiver elevation function that can be used for datuming.
The effectiveness of the proposed method has been demonstrated on two high-resolution seismic reflection datasets, in which frequency-wavenumber deghosting would have been inapplicable, showing the gain in pre-stack bandwidth and the improvement in the interpretability of the final stack images. This could be valuable for near-surface marine seismic acquisition with little control on the receiver streamer depth, as well as for the broadband re-processing of legacy seismic reflection data. Multi-channel processing flow for GI-gun and Sparker multichannel reflection data. This is a standard seismic processing flow for 2D multi-channel seismic reflection data, integrated with the timefrequency-offset deghosting presented in this paper