TY - JOUR
AU - Smith, Patrick J
AB - Abstract Since explosive and impulsive seismic sources such as dynamite, air guns, gas guns or even vibroseis can have a big impact on the environment, some companies have decided to record ambient seismic noise and use it to estimate the physical properties of the subsurface. Big challenges arise when the aim is extracting body waves from recorded passive signals, especially in the presence of strong surface waves. In passive seismic signals, such body waves are usually weak in comparison to surface waves that are much more prominent. To understand the characteristics of passive signals and the effect of natural source locations, three simple synthetic models were created. To extract body waves from simulated passive signals we propose and test a Radon-correlation method. This is a time-spatial correlation of amplitudes with a train of time-shifted Dirac delta functions through different hyperbolic paths. It is tested on a two-layer horizontal model, a three-layer model that includes a dipping layer (with and without lateral heterogeneity) and also on synthetic Marmousi model data sets. Synthetic tests show that the introduced method is able to reconstruct reflection events at the correct time-offset positions that are hidden in results obtained by the general cross-correlation method. Also, a depth migrated section shows a good match between imaged horizons and the true model. It is possible to generate off-end virtual gathers by applying the method to a linear array of receivers and to construct a velocity model by semblance velocity analysis of individually extracted gathers. Radon transform, Green's functions, passive seismic signals, cross-correlation, retrieving body waves 1. Introduction In most areas, the origins of passive seismic signals are cultural (e.g. cars and trains) and natural (wind, distant earthquakes and ocean-waves) sources that produce substantial strong propagated surface waves. The lack of deep noise sources in the subsurface hinders the extraction of reflection arrivals from the correlation of passive seismic signals (Forghani & Snieder 2010). To retrieve the complete Green's function from passive signals, a uniform distribution of sources is essential (Wapenaar et al. 2011) otherwise non-physical seismic arrivals are unavoidable (King & Curtis 2012). Cross-correlation (Claerbout 1968; Lobkis & Weaver 2001; Shapiro & Campillo 2004; Bakulin & Calvert 2006; Campillo 2006; Curtis et al. 2006; Godin 2006; Draganov et al. 2007; Boué et al. 2014), time reversal (Van Tiggelen 2003), cross-coherency (Nakata et al. 2011), deconvolution interferometry (Vasconcelos & Snieder 2008), passive inverse filtering (Gallot et al. 2012) and multi-dimensional deconvolution (Wapenaar et al. 2011) have been used to correlate passive seismic signals at spatially distributed receivers to recover coherent seismic waves. These recovered seismic wavefields are highly contaminated by the ground roll (Wapenaar 2003) and need post-correlation processing to highlight the weak reflection events. Also, a non-uniform distribution of sources (as happens in almost all cases) causes the correlation results to be biased by directionality and also causes non-complete reconstruction of the Green's function (Vasconcelos & Snieder 2008; Wapenaar et al. 2011), which creates more difficulties (appearance of non-physical events) in retrieving the reflection arrivals correctly. To avoid bias from the out-of-receiver plane sources it is essential to use beamforming (Behr et al. 2013; Löer et al. 2018) to create a semi-zero azimuth profile. This helps to choose the part of the recorded signal energy that comes from sources that are inline with receivers. As mentioned, the lack of noise sources in the subsurface beneath the reflectors is a substantial reason why only very weak amplitude reflection arrivals are recovered by cross-correlation. 2D FK filtering is a well-known method to attenuate the events inside the reject velocity zone. This filter performs poorly in terms of dealing with spatial aliasing, causes amplitude smearing, can easily miss very weak amplitude reflections and also attenuates not only the ground roll but also the reflection arrivals inside the noise cone, which are not desirable. These pitfalls of FK filtering encouraged us to focus on finding a better method for body wave extraction from passive signals that could be used to image the subsurface. Radon correlation (RC) is a time-spatial correlation that looks for localised maximum distributed energy in a hyperbolic path for body waves with apparent propagation velocities higher than the surface wave velocity. We introduce RC as a method to correlate passive signals to suppress the ground roll, which is the dominant event in the correlated section, and to enhance weak reflection arrivals. It is easy to estimate the true propagation velocity of ground roll from the correlated section. The minimum velocity of the velocity span in RC should be higher than the surface wave velocity. In this way, linear coherent events (ground roll) will be rejected in the forward transform procedure based on the propagation velocities and the amount of move out. The result of RC will be to highlight hyperbolic events with propagation velocities in the chosen velocity range. To mitigate against the creation of artificial events due to amplitude smearing it is vital to stack many processed windows. To test the method, passive signals have been simulated from the propagation of elastic waves through known 2D synthetic velocity models. First, we use a simple two-horizontal-layer model with 301 receivers with a 5 m interstation interval, excited by 80 randomly distributed sources. Second, we use two simple three-layer models including a dipping layer (with and without lateral heterogeneity) with 301 receivers with a 5 m interstation interval, excited by 100 randomly distributed sources and third, we employ the complex Marmousi model (Versteeg 1994) with 381 vertical sensors laid out with a 10 m interstation interval and a 2200 m minimum offset. To simulate passive signals, only sources on the surface of the model are considered, which mimics the dominant source distribution for ambient noise seismic signals. We assumed that no deep sources are available and that the origin of the passive signals could effectively correspond to trains, road traffic etc. A depth migrated section was created by migration of virtual shot gathers created by RC from simulated passive signals due to propagation of elastic waves, excited by 100 arbitrarily distributed sources along surface (minimum offset is 2200 m), through the synthetic Marmousi model (Versteeg 1994). 2. Passive seismic signal simulation The finite difference method is used extensively for the modelling of seismic wave propagation (Aki & Larner 1970; Virieux 1986; Thorbecke & Draganov 2011). It has the advantages of being a simple but powerful technique that is easy to implement in a modelling code. Since passive seismic signals originating from surface sources are dominated mostly by surface waves, to simulate them it is essential to solve the elasto-dynamic equations with appropriate free-surface boundary conditions (e.g. Virieux 1986). To retrieve complete Green's functions from the correlation of passive signals through a 2D profile, the receivers should ideally be surrounded by the seismic sources (mutually uncorrelated) distributed statistically uniformly on a large radius half circle/sphere, and most importantly the noise fields should be diffuse and equipartitioned (Wapenaar et al. 2011). Under these circumstances, deep ‘passive’ sources can lead to the retrieval of valuable reflected body waves in correlated data and relatively weaker surface waves. However, major sources of high frequency passive signals such as railways and road traffic are not uniformly distributed and are biased towards the free surface. To examine the effect of these kinds of sources, we simulate passive signals due to the propagation of elastic waves through different models (which will be described in the following sections) using randomly distributed sources in space and time. The aim here is not to model the specific sources associated with traffic, but more to focus on the effects of bias in their spatial locations (e.g. on the free surface, in the far-field etc.). Vertical component synthetic sensors are used to record the source wavefield, and absorbing boundaries have been used to attenuate the reflections from the model edges (these apply to three of the edges in the 2D model, with a free surface at the top). For stacking purposes, the simulations are repeated many-times, which we called realisations, and each time the source locations are distributed randomly in space and time. The source wavelets are mixed-phase Ricker wavelets with 20 Hz dominant frequencies and vertical forces are used to excite the medium for all simulations. The modelling parameters for all models are summarised in Table 1. Table 1. Summarised modelling parameters for all simulations 1. Vertical sensors are used to record the wavefield. 2. The maximum time delay of shots is randomly chosen between 0 and 100 ms. 3. The source wavelets are phase-shifted Ricker wavelets with a 20-Hz dominant frequency 4. Vertical forces are used to excite the medium 5. Absorbing boundaries have been used to attenuate the reflections from the model edges 6. Free-surface boundary condition is used during the all simulations 7. The sources are highlighted by turquoise asterisks. 8. The receivers are highlighted by green V-shapes 1. Vertical sensors are used to record the wavefield. 2. The maximum time delay of shots is randomly chosen between 0 and 100 ms. 3. The source wavelets are phase-shifted Ricker wavelets with a 20-Hz dominant frequency 4. Vertical forces are used to excite the medium 5. Absorbing boundaries have been used to attenuate the reflections from the model edges 6. Free-surface boundary condition is used during the all simulations 7. The sources are highlighted by turquoise asterisks. 8. The receivers are highlighted by green V-shapes Open in new tab Table 1. Summarised modelling parameters for all simulations 1. Vertical sensors are used to record the wavefield. 2. The maximum time delay of shots is randomly chosen between 0 and 100 ms. 3. The source wavelets are phase-shifted Ricker wavelets with a 20-Hz dominant frequency 4. Vertical forces are used to excite the medium 5. Absorbing boundaries have been used to attenuate the reflections from the model edges 6. Free-surface boundary condition is used during the all simulations 7. The sources are highlighted by turquoise asterisks. 8. The receivers are highlighted by green V-shapes 1. Vertical sensors are used to record the wavefield. 2. The maximum time delay of shots is randomly chosen between 0 and 100 ms. 3. The source wavelets are phase-shifted Ricker wavelets with a 20-Hz dominant frequency 4. Vertical forces are used to excite the medium 5. Absorbing boundaries have been used to attenuate the reflections from the model edges 6. Free-surface boundary condition is used during the all simulations 7. The sources are highlighted by turquoise asterisks. 8. The receivers are highlighted by green V-shapes Open in new tab 3. Radon correlation RC is a time-spatial correlation that looks for localised maximum distributed energy through different hyperbolic paths for body waves with apparent propagation velocities higher than surface wave velocity. 3.1. RC in the frequency domain The forward Radon transform (RT), introduced by Radon (1917), can be expressed as the cross-correlation of data in the time-offset domain with a train of Dirac delta functions, which form the basis functions of the RT (Gholami 2017; Mokhtari et al. 2019): $$\begin{equation}u\!\left( {\tau ,p} \right) = \int\!\!\!\! \int s\left( {t,x} \right)\sigma \left( {t - \sqrt {{\tau ^2} + {{\left( {xp} \right)}^2}} } \right)dtdx,\end{equation}$$(1) where, |$s( {t,\,\,x} )\,\,$|is the original seismogram in the time-offset domain, |$u( {\tau ,\,\,p} )$| is the data in the Radon domain which is called the RT coefficients, |$\sigma \,\,$|is Dirac delta function, t is the two-way time, x is the offset, τ is the intercept two-way time and |$p$| is the slowness. In equation (1) the correlation is applied in a hyperbolic trajectory, and the inverse RT can be written as $$\begin{equation}s\!\left( {t,x} \right) = \int\!\!\!\! \int u\left( {\tau ,p} \right)\sigma \left( {t - \sqrt {{\tau ^2} + {{\left( {xp} \right)}^2}} } \right)d\tau dp.\end{equation}$$(2) Since estimating the RT coefficients in the time domain is computationally expensive, it is much better to apply it in the frequency domain. To be able to perform it correctly in the frequency domain it is essential to apply time-stretching |$( {{t^2} \to t^{\prime}} )$| before domain transformation (Yilmaz 1989, 2001; Gholami & Farshad 2019). After time-stretching, equation (2) can be written as follows: $$\begin{equation}s\!\left( {t^{\prime},x} \right) = \int\!\!\!\!\int u\left( {\tau ^{\prime},p} \right)\delta \left( {t^{\prime} - \tau ^{\prime} - {{\left( {xp} \right)}^2}} \right)d\tau ^{\prime}dp.\end{equation}$$(3) Time-stretching changes the hyperbolic events to semi-parabolic ones (Yilmaz 1989), which helps to estimate the RT coefficients in the frequency domain more efficiently. Applying a Fourier transform on both sides of equation (3) with respect to |$t^{\prime}$| gives: $$\begin{eqnarray}s\!\left( {w,x} \right) &=& \int\!\!\!\! \int\!\!\!\! \int u\left( {\tau ^{\prime},p} \right)\delta \left( {t^{\prime} - \tau ^{\prime} - {{\left( {xp} \right)}^2}} \right)\nonumber\\ &&\times\, ex{p^{\sqrt{ - 1} .wt^{\prime}}}d\tau ^{\prime}dpdt^{\prime},\end{eqnarray}$$(4) where |$w\,\, = \,\,2\pi f$| is the angular frequency of stretched time and |$s( {w,{\rm{\,\,}}x} )$| is the seismic gather in the frequency-offset domain. While noting that $$\begin{equation}\int \delta (t^{\prime} - \tau ^{\prime} - {(xp)^2}) ex{p^{\sqrt { - 1}.wt^{\prime}}} dt^{\prime} = exp{^{\sqrt{-1}_{{.w(\tau}^{\prime}{+{{(xp)}^{2}})}}}},\end{equation}$$(5) it is possible to summarise equation (4) as: $$\begin{equation}s\!\left( {w,x} \right) = \int u\left( {w,p} \right)ex{p^{\sqrt { - 1} *2\pi f{{\left( {xp} \right)}^2}}}dp,\end{equation}$$(6) where |$u( {w,{\rm{\,\,}}p} )$| is the RT coefficients in the frequency-slowness/velocity (⁠|$f - v)$| domain. Equation (6) can be written in matrix form as follows: $$\begin{equation}{S_{Mf}} = \,\,{{\Psi}}{U_{Mf}},\end{equation}$$(7) where, |${S_{Mf}}$| is the matrix form of seismic data in the frequency domain, |${U_{Mf}}$| is the matrix form of the RT coefficients (model space) in the frequency domain and |${{\Psi }}$| is the Kernel matrix. In the Radon domain: $$\begin{equation}{U_{Mf}} = \frac{{{{{\Psi }}^T}.{S_{Mf}}}}{{{{{\Psi }}^T}.{{\Psi }} + \lambda I}}.\end{equation}$$(8) In equation (8) the RT coefficients are obtained by applying a minimum-least-square criterion, where, |${{{\Psi }}^T}$| is the transpose of |${{\Psi }}$|⁠, λ is a pre-whitening percentage and |$I$| is the identity matrix. If in equation (8) instead of |${S_{Mf}}$|⁠, which is the matrix that includes the frequency spectrum of seismic signals, a coherency matrix obtained by correlation of a reference trace with other traces is used, the RC coefficients in the Radon domain can be obtained by $$\begin{equation}FR{C_{Mf}} = \frac{{{{{\Psi }}^T}.{C_{Mf}}}}{{{{{\Psi }}^T}.{{\Psi }} + \lambda I}}\,\,,\end{equation}$$(9) where FRCMf is the forward RC coefficients in the Radon domain, |${C_{{M_f}}} = \mathop \sum \nolimits_{i = 1}^N {C_{{N_i}}}\,\,$|⁠, |${C_{{N_i}}} = \,\,FFT( {{T_s}( {IFFT( {{C_M}} )} )} )\,\,$|where |$FFT$|⁠, |$IFFT$| and |${T_s}$| are the forward Fourier transform, inverse Fourier transform and time-stretching, respectively. Also, |${C_M}\,\,$|is the coherency matrix of the |${i^{th}}$| windowed signals that is estimated by the multi-taper method coherence estimation (MTMC) as proposed by Thomson (1982). MTMC is a periodogram based method that uses a bunch of family group tapers to form independent estimates of the spectral density, before the average of all individual estimates gives the final coherence (Naghadeh & Morley 2016). The inverse RC |$( {IRC} )$| can be written as follows: $$\begin{equation}IRC\,\,\left( {t^{\prime},x} \right) = \int\!\!\!\!\int FR{C_{Mf}}.ex{p^{ - 2\pi f\sqrt { - 1} \left( {t^{\prime} + {{\left( {xp} \right)}^2}} \right)}}dpdf.\end{equation}$$(10) The final step to highlight the hyperbolic events in the section is to conduct time squeezing |$( {t^{\prime} \to {t^2}} )$| in the time domain: $$\begin{equation} IRC\left( {t^{\prime},x} \right)\mathop \Longrightarrow \limits^{time - squeezing} d\left( {t,x} \right).\end{equation}$$(11) 3.2. Discrete RC in the time domain For discrete signals the time-domain RC can be obtained by $$\begin{equation}RC\,\,\left({t,x} \right) = \mathop \sum \limits_p \mathop \sum \limits_x \mathop \sum \limits_{w = 1}^{w = N} \mathop \sum \limits_T d_w^ {={\times}} \left({T,{x_0}}\right).{d_w}\left( {t + T,x} \right),\end{equation}$$(12) where |${{1}/{{{V_{{\rm{max}}}}}}} \le p \le {{1}/{{{V_{{\rm{min}}}}}}},$||${V_{{\rm{min}}}}$| and |${V_{{\rm{max}}}}$| are the minimum and maximum velocities of the velocity range during RC procedure, |$ - nt.dt \le T \le nt.dt$| (⁠|$nt$| is the number of time samples and |$dt$| is the time sample interval), |$t\,\, = \sqrt {{\tau ^2} + {{( {xp} )}^2}} \,\,$|⁠, |$\tau $| is the intercept-time, |${x_0} \le x \le {x_{max}}$|⁠, |${x_0}$| is the distance between the virtual source and the closest receiver (minimum offset) and |${x_{max}}$| is the maximum offset, |${d_w}$| is the windowed data and |$d_w^{={\times}}$| is the complex conjugate of |${d_w}$|⁠. 3.3. The RC procedure as a flowchart In this section, the RC procedure is summarised in a flowchart (figure 1), which shows the main processing steps used to carry out the elimination of the ground roll and highlight the weak reflection arrivals using time-spatial correlation in the Radon domain. Long window lengths of recorded signals are chosen to make N tapered sub-windows using 2N– 1 family group tapers for each individual window. The final coherence spectrum between the reference and desired traces is the average of the estimated coherences for each set of tapered sub-windows. Analysing the coherence section in the |$f - v$| map gives the ability to discriminate between different events based on their different move out and propagation velocity. Choosing a suitable velocity range and amplitude threshold helps to supress and pass some arrivals during the iterative forward RC transform. The domain change from model space to data space (called inverse RC) highlights the most probable reflection arrivals. The beamforming step is added to the flowchart during the processing of real data sets to reduce the effects of bias from sources that are out of plane with the receivers. Beamforming helps to localise the maximum energy sources that can be used to narrow the stationary phase zone, which is useful for retrieving the reflection arrivals. To improve the resolution, it is possible to predefine the dominant frequency and propagation velocity of the ground roll from the correlated section (as surface waves are well recovered in the raw correlograms), as the source characteristics. Figure 1. Open in new tabDownload slide The RC procedure as a flowchart. Choosing a semi-zero azimuth by beamforming is a part of the processing applied to a real data set. Forward and inverse RC are applied in the frequency domain. The amplitude threshold, called scaling, manages the removal of low amplitude focused energies in the |$f - v$| map. Figure 1. Open in new tabDownload slide The RC procedure as a flowchart. Choosing a semi-zero azimuth by beamforming is a part of the processing applied to a real data set. Forward and inverse RC are applied in the frequency domain. The amplitude threshold, called scaling, manages the removal of low amplitude focused energies in the |$f - v$| map. 4. Two-layer synthetic data example To be able to reconstruct reflection events from the passive signals, the reflection energy should be transmitted by a virtual source at the surface to the other receivers in the profile (Snieder 2004, 2007). To evaluate the applicability of RC to our dataset, and to understand the role of multiples in retrieving the reflection arrivals, two different time records of simulated passive signals have been used. The first trial includes 100 s (1-second simulation repeated a hundred times) total time simulated passive signals (figure 2b) resulting from the activity of 80 arbitrarily distributed sources (figure 2a) along the surface between –1000 m and –800 m from the first receiver. The receivers are located in positions between 0 and 1500 m with a 5 m station interval. For each realisation, the locations of the sources are randomly chosen. The result of normalised stacked cross-correlation of 100-realisations is shown in figure 3a. To correlate the first trace with all other traces along the receiver line, we used a correlation window length of 10 s with 90% overlap of two consecutive windows. Only the surface wave (ground roll) and direct P-wave are visible events in the correlated section. Since multiples cannot be seen in a short 1 s record, the cross-correlation (figure 3a) and RC (figure 3b) results do not contain true retrieved reflection events. In the second trial, 500 s of simulated signals (100 realisations) are used to retrieve the Green's function. The same geometry for the sources and receivers has been used with the only difference being the total record time. A 5 s record of each realisation includes PP, PSP and SSP multiples, which gives the virtual source the ability to transmit energy to other receivers. However, the amplitudes of multiples arising from repeated bounces before reaching the virtual source dramatically decrease, but the stacking of numerous analysis windows boosts these amplitude levels. The correlation of long duration signals, where a 15 s window length with 80% overlap between consecutive windows is used to conduct cross-correlation, shows (figure 3c) the importance of multiples (‘scattering’) in aiding retrieval of reflection events from passive signals with a non-uniform distributions of background sources. In figure 3d, the result of RC on simulated signals is shown, where a 50 s window length with 20% overlap between consecutive windows and a velocity range of 3500 to 4500 m s–1 was used. During the correlation, a bunch of family group tapers belonging to discrete prolate spheroidal sequences (Delsarte et al. 1985) are applied to each window of signals to make 10 tapered sub-windows. These sub-windows are then used to estimate the final coherence spectrum between the reference trace and desired one. The final coherence section (from stacking of all of the estimated coherences for all windows) is prepared (equation 8) for analysis in the |$f - v$| map. The maximum amplitude |$({M_A})$| in the |$f - v$| map is defined by $$\begin{equation}{M_A} = {{{\left| {\mathop \sum \limits_{i = 1}^{NT} {A_{\left( {f,i} \right)}}\,.\,{x_i}} \right|}}/{{NT}}},\end{equation}$$(13) where |${A_{( {f,i} )}}$| is the amplitude of |${i^{th}}$| trace at frequency f, |${x_{i\,\,}}{\rm{is\,\,the\,\,offset\,\,of\,\,the\,\,}}{i^{th}}{\rm{\,\,trace}},{\rm{\,\,\,\,\,\,and\,\,}}NT$| is the total number of traces in the section. During forward RC all amplitudes lower than |${M_A}/S,$| where |$S > 1$| (the scale chosen by user), are removed. The result of applying RC to the simulated signals shows that this approach can highlight the reflection events following the elimination of ground roll, as RC looks for the maximum cumulative energy through a hyperbolic path and the move out of ground roll is linear not hyperbolic. During forward RC, all amplitudes along a linear event do not stack together to create large localised cumulative energy. This small cumulative energy in the |$f - v$| map causes the elimination of (surface wave) ground roll (and also the direct body wave) from the correlated section due to scaling. Comparing the result of RC in figure 3d with the active shot record (figure 2d, with a source at 0 m position) demonstrates the ability of RC to remove ground roll and to highlight the correct reflection arrival. Figure 2. Open in new tabDownload slide (a) P-velocity model with sources (light blue asterisks) laid out between the −1000 and −800 m position and receivers (green V-shapes) are located at the 0 to 1500 m position; (b) 100 s total time simulated passive signals; (c) 500 s total time simulated passive signals; (d) active shot record for the model in figure 1 – the shot is located at position (⁠|$x\,\, = {\rm{\,\,}}0{\rm{\,\,m}},{\rm{\,\,}}z\,\, = {\rm{\,\,}}0{\rm{\,\,m}}$|⁠). Several events are recognisable: (i) direct P-arrival, (ii) ground roll, (iii) reflected P-wave from the interface at 800 m depth and (iv) Reflected PS-wave from the interface at 800 m depth. Figure 2. Open in new tabDownload slide (a) P-velocity model with sources (light blue asterisks) laid out between the −1000 and −800 m position and receivers (green V-shapes) are located at the 0 to 1500 m position; (b) 100 s total time simulated passive signals; (c) 500 s total time simulated passive signals; (d) active shot record for the model in figure 1 – the shot is located at position (⁠|$x\,\, = {\rm{\,\,}}0{\rm{\,\,m}},{\rm{\,\,}}z\,\, = {\rm{\,\,}}0{\rm{\,\,m}}$|⁠). Several events are recognisable: (i) direct P-arrival, (ii) ground roll, (iii) reflected P-wave from the interface at 800 m depth and (iv) Reflected PS-wave from the interface at 800 m depth. Figure 3. Open in new tabDownload slide (a) The result of cross-correlation applied to 100 s short record simulations; (b) the RC result applied to 100 s short record simulations. It does not contain true retrieved reflection event and includes the highlighted non-physical event. For short records, due to the lack of multiples, the virtual source is not able to transmit reflection energy to virtual receivers. (c) The result of cross-correlation applied to 500 s long record simulations. Amplitude discontinuities and fluctuations through the hyperbolic path are caused by shadow linear events. (d) The RC result applied to 500 s long record simulations. Amplitude continuity and equal distribution of energy through the hyperbolic path are due to the interpolation power in the Radon domain. Figure 3. Open in new tabDownload slide (a) The result of cross-correlation applied to 100 s short record simulations; (b) the RC result applied to 100 s short record simulations. It does not contain true retrieved reflection event and includes the highlighted non-physical event. For short records, due to the lack of multiples, the virtual source is not able to transmit reflection energy to virtual receivers. (c) The result of cross-correlation applied to 500 s long record simulations. Amplitude discontinuities and fluctuations through the hyperbolic path are caused by shadow linear events. (d) The RC result applied to 500 s long record simulations. Amplitude continuity and equal distribution of energy through the hyperbolic path are due to the interpolation power in the Radon domain. The results in figure 3 emphasise the importance of multiples (‘scattering’) in retrieving the reflection arrivals. The correlation of short duration simulated signals that do not include the multiples does not recover true reflection arrival and contains non-physical reflection events due to the correlation of different arrivals (P, PS, PP, PSP, SSP, S) from the same/different interfaces. In this case, the RC method highlights non-physical events because of the absence of a true reflection arrival. In the second trial in section 6, long records that include multiples allow for the retrieval of the true reflection arrival. Cross-correlation of simulated signals (figure 3c) includes the true reflection and non-physical events. Three steps that are desirable to help RC to retrieve reflection arrivals are: (1) weighting traces by offset during forward RC. This step magnifies the cumulative energy for a hyperbolic event with propagation velocity in the range of the velocity span. (2) Using a high number of iterations to achieve sparsity in the model space during the inversion procedure (here 150 iterations have been used). (3) Finally, the scaling factor. The scaling factor helps to suppress lower cumulative energies than the threshold and so preserves higher local cumulative energies than threshold in the |$f - v$| map. In figure 3c, the retrieved reflection arrival at offsets between 0 and 700 m from the virtual source (0 m position) has strong amplitude because of multiples, from 700 to 1100 m the amplitude of event decreases dramatically indicating a decrease of multiple amplitudes due to many bounces, and after 1100 m there is not enough energy to reach to the far offset receivers. The result of RC (figure 3d) shows a hyperbolic event with constant distributed energy through the full offset. The inverse RC procedure distributes the local cumulative energy in the |$f - v$| map over a hyperbolic event through all the offsets. If lateral velocity changes exist after position 1000 m, due to heterogeneity or anisotropy, the result of RC would be the same as the highlighted event (figure 3d) due to the lack of energy reaching offsets larger than 1000 m. 5. Three-layer model without lateral heterogeneity To further explore and test the RC method, a three-layer model (figure 4a) was chosen where the first layer is a dipping layer with a 3000 m s–1 P-wave velocity. A synthetic-active-shot record with a source located at a 0 m position is shown in Figure 4b. The receiver profile includes 201 receivers with a 5 m interstation interval that are located between 0 and 1000 m in position. To simulate passive seismic signals, 50 sources (turquoise asterisks) are randomly distributed in time (the maximum delay time between the random shots is 100 ms) and space (between the –1000 and −800 m positions). A total simulation time of 500 s is used for 100 realisations to retrieve the reflection arrivals. For each realisation the locations of the sources are randomly chosen. To correlate the signals during the RC procedure, a 25 s long window with 50% overlap between consecutive windows and a 2500 to 5000 m s–1 velocity range is used. Figure 4. Open in new tabDownload slide (a) Three-layer velocity model with P-wave velocities (Vp1 = 3000 m s–1, Vp2 = 4000 m s–1, Vp3 = 5000 m s–1). Source (turquoise asterisks) are randomly distributed between the −1000 and −800 m position and receivers (green V-shapes) are between the 0 and 1500 m position, (b) Synthetic-active-shot record. The source is located at the 0 m position and the source wavelet frequency is 20 Hz. Figure 4. Open in new tabDownload slide (a) Three-layer velocity model with P-wave velocities (Vp1 = 3000 m s–1, Vp2 = 4000 m s–1, Vp3 = 5000 m s–1). Source (turquoise asterisks) are randomly distributed between the −1000 and −800 m position and receivers (green V-shapes) are between the 0 and 1500 m position, (b) Synthetic-active-shot record. The source is located at the 0 m position and the source wavelet frequency is 20 Hz. To show the effect of scaling and iteration number during retrieval of the arrivals, and also to examine the effect of the energy path on the created virtual gathers, 100 gathers were extracted. Each gather includes 201 traces with a virtual source station interval of 5m. Figure 5a shows the virtual gather that was created with a virtual source located at a 0 m position, where the scale and iteration number are 5 and 10, respectively. The extracted gather includes many artificial events that appear earlier than the P-wave reflection arrival from the first interface. Also, a discontinuity appears through the retrieved reflection arrival. Decreasing the scaling factor from 5 to 2 while using the same iteration number (10 iterations) helped to extract the twin continuous reflection arrivals at the corresponding time-space coordinate of the first reflection arrival, with considerably fewer artificial events (figure 5b). Raising the iteration number from 10 to 100 while keeping the same scale (figure 5c and d) improves the highlighting of events (P-wave reflection arrivals from the first and second interfaces) followed by suppressing the artificial events (figure 5d). The created virtual gather with a virtual source at the 500 m position is compared with the active shot record with the same source location of 500 m (figure 6). Due to the longer ray-path, less energy reaches the virtual source at 500 m and this restricts the ability to remove the non-physical events. This is apparent even when using a high (100) iteration number and low scaling factor, with non-physical events appearing at an earlier time than the first-layer P-wave reflection (figure 6b). The reflection arrival from the first interface is detected correctly at the relevant time-space with a low reflection amplitude. The difference between the apparent velocity and the true propagation velocity is the mainly the result of the dipping geometry. The dipping geometry is not a challenge for the retrieval of a true reflection event. A quasi-random wavefield, arriving from various directions at the virtual source locations plays the most important role in correct extraction of the reflection arrival. Figure 5. Open in new tabDownload slide Virtual gather created with a virtual source located at 0 m position by RC due to: (a) scale = 5, iteration number = 10; (b) scale = 2, iteration number = 10; (c) scale = 5, iteration number = 100 and (d) scale = 2, iteration number = 100. Figure 5. Open in new tabDownload slide Virtual gather created with a virtual source located at 0 m position by RC due to: (a) scale = 5, iteration number = 10; (b) scale = 2, iteration number = 10; (c) scale = 5, iteration number = 100 and (d) scale = 2, iteration number = 100. Figure 6. Open in new tabDownload slide (a) Active shot record with a source located at the 500 m position and (b) virtual gather created by RC where the scale = 2, iteration number = 100 and the virtual source located at the 500 m position. Non-physical reflection events were not eliminated by scaling and high iteration number due to the very weak transmitted energy by the virtual source at the 500 m position through the virtual receiver profile. Figure 6. Open in new tabDownload slide (a) Active shot record with a source located at the 500 m position and (b) virtual gather created by RC where the scale = 2, iteration number = 100 and the virtual source located at the 500 m position. Non-physical reflection events were not eliminated by scaling and high iteration number due to the very weak transmitted energy by the virtual source at the 500 m position through the virtual receiver profile. 6. Three-layer model with lateral heterogeneity The existence of lateral heterogeneity forces the linear coherent ground roll event to appear with different linear move outs, and at different offset locations there is also a decrease/increase of acoustic impedance contrast between horizons. In this case, lateral heterogeneity decreases the velocity contrast at the first interface at far offsets which causes low reflection amplitudes, and consequently results in very low multiple (‘scattering’) amplitudes that are harmful to the seismic interferometry. To understand the effect of this heterogeneity on the results given by RC, the lateral heterogeneity was added to the first layer. The velocity of the dipping layer increases latterly from west to east with a 20% gradient. To simulate the passive signals from the model, 50 sources are distributed randomly in space (−1000 to −800 m position) and time (100 ms is the maximum delay between shots), with 301 receivers laid out between positions 0 to 1500 m of the model (figure 7a). A total 500 s simulation time (100 realisations with 5 s recorded signals) was used to create a virtual gather with a virtual source at the 0 m position, which includes the 301 virtual receivers. To perform RC, a 20 s window length with a 50% overlap between two consecutive windows and a velocity range of 2500 to 5000 m s–1 was used. The scale and iteration number are 2 and 100, respectively. The RC result (figure 7c) shows a highlighted hyperbolic event that matches with the first-layer P-wave reflection arrival shown in the active shot record (figure 7b) for offsets between 0 and 1000 m. At far offsets beyond 1000 m (dashed line, figure 7c), where the lateral heterogeneity causes an increase of 500 m s–1 in the P-wave velocity, a maximum of 30 ms deviation in the reflection arrival is visible, occurring at the 1500 m offset. The lateral increase of the velocity in the first layer leads to large reflection angles to the virtual source at the surface and bigger reflection angles for the multiples propagate through the virtual receiver profile. This causes the transmission of reflection energy to offsets beyond 1500 m and an absence of energy between 1000 to 1500 m offset. To ameliorate the retrieved reflection arrival deviation at far offsets, it is essential to extract the virtual gathers for partial offsets and not for full offsets. Figure 7. Open in new tabDownload slide (a) Three-layer velocity model with P-wave velocities (Vp1 = 3000 m s–1, Vp2 = 4000 m s–1, Vp3 = 5000 m s–1). The first layer is a heterogeneous layer with a 20% lateral velocity gradient toward the east. Sources (turquoise asterisks) are randomly distributed between the −1000 and −800 m position and receivers (green V-shapes) are between the 0 and 1500 m position; (b) Synthetic-active-shot record. The source is located at the 0 m position and the source wavelet frequency is 20 Hz. (c) Created virtual shot gather by RC. The virtual source is located at the 0 m position and the gather includes 301 virtual receivers with a 5 m interstation interval. Figure 7. Open in new tabDownload slide (a) Three-layer velocity model with P-wave velocities (Vp1 = 3000 m s–1, Vp2 = 4000 m s–1, Vp3 = 5000 m s–1). The first layer is a heterogeneous layer with a 20% lateral velocity gradient toward the east. Sources (turquoise asterisks) are randomly distributed between the −1000 and −800 m position and receivers (green V-shapes) are between the 0 and 1500 m position; (b) Synthetic-active-shot record. The source is located at the 0 m position and the source wavelet frequency is 20 Hz. (c) Created virtual shot gather by RC. The virtual source is located at the 0 m position and the gather includes 301 virtual receivers with a 5 m interstation interval. 7. Synthetic Marmousi model Imaging the complex Marmousi model using simulated passive signals presents a big challenge for all correlation methods. We used raytracing to determine the minimum recording time needed for multiples from different reflectors to reach the receivers at the surface, and also to identify the area covered by passing energy. In figure 8c and d, sources are located at 0 and 500 m positions, respectively. The red curves show the propagated ray-paths from the source to receivers and rays highlighted in turquoise are the scattered energies from deep scatter points that return energy to the surface at locations of virtual sources with small reflection angles. The amplitudes of these energies are much stronger than multiples, which bounce from different horizons before reaching the virtual sources. Also, the arrival times of these scattered energies from deep scatter points are much earlier than the multiples. On the basis of the raytracing results, it is obvious that the area with the best coverage is between positions 3000 and 4800 m. Figure 8. Open in new tabDownload slide (a) Marmousi P-velocity model. Turquoise asterisks show the location of randomly distributed sources and the green V-shapes are 381 vertical receivers. The minimum offset is 2200 m and the maximum offset is 6500 m; (b) Marmousi S-wave velocity model; (c) the propagated rays (red curves) from the source at the 0 m position to receivers and (d) the paths of rays travelling from the source at the 500 m position to receivers. Figure 8. Open in new tabDownload slide (a) Marmousi P-velocity model. Turquoise asterisks show the location of randomly distributed sources and the green V-shapes are 381 vertical receivers. The minimum offset is 2200 m and the maximum offset is 6500 m; (b) Marmousi S-wave velocity model; (c) the propagated rays (red curves) from the source at the 0 m position to receivers and (d) the paths of rays travelling from the source at the 500 m position to receivers. To estimate the approximate minimum simulation time, multiple raytracing for different reflector depths has been undertaken (figure 9). The arrival times for multiples from a reflector at a 650 m depth to the farthest offset for the PP mode and converted PS mode are about 4.9s and 5.5 s. Also, the arrival time of multiples that bounce between the surface and a 1500 m maximum depth for PP and PS modes are 6.4 and 7.5 s, respectively. To test the ability of RC, and to cover possible additional interlayer trapped energies, we decided to simulate a total 2000 s (200 × 10 s) time record (based on the raytracing) of passive signals. For each 10 s recording, the sources are distributed randomly between the 0 and 500 m position of the model at the surface. The receiver line is located between 2700 and 6500 m and the station interval is 10 m (figure 8a). Here the aim is to extract the virtual gathers and then migrate them in order to image the reflectors. To do so, 47 off-end shot gathers that include 201 channels were created. The sampling frequency is 500 Hz, the shot interval is 50 m and the source wavelet frequency is 20 Hz. Because of the steeply dipping reflectors and the high lateral velocity gradient, a part of the reflection energy could be turning waves that are essential for imaging steeply dipping reflectors. However, these add more difficulties in imaging structures using virtual gathers. In addition, non-physical reflections due to correlation of different reflection phases from different interfaces hinder the imaging and velocity analysis. Figure 10a shows the part of the model where all the virtual sources and receivers are located (between 2700 and 6500 m positions). P-wave, SSP-wave and PSP-wave multiples all play an important role in allowing the virtual sources to transmit the energy to other receivers. To create the gathers by RC, a velocity range of 1300 to 4000 m s–1 was used. The pre-stack reverse time migration method (Yoon & Marfurt 2006) is used to migrate the gathers using a highly smoothed version of the true velocity model (figure 10b). The imaged reflectors (figure 10c) show a good match with the true velocity model (figure 10a), and this demonstrates the ability of the RC method to retrieve the reflection arrivals from passive signals that are highly contaminated by both the non-physical events and surface waves generated by sources at the surface. Figure 9. Open in new tabDownload slide Multiple raytracing for maximum 650 m depth: (a) propagated PP mode multiples and (b) propagated PS mode multiples. Multiple raytracing for a maximum 1500 m depth. (c) Propagated PP mode multiples and (d) propagated PS mode multiples. The maximum arrival time of multiples bouncing between the surface and 1500 m depth is 7.5 s. Figure 9. Open in new tabDownload slide Multiple raytracing for maximum 650 m depth: (a) propagated PP mode multiples and (b) propagated PS mode multiples. Multiple raytracing for a maximum 1500 m depth. (c) Propagated PP mode multiples and (d) propagated PS mode multiples. The maximum arrival time of multiples bouncing between the surface and 1500 m depth is 7.5 s. Figure 10. Open in new tabDownload slide (a) Zoomed-in part of the velocity model (figure 8a) that includes all the virtual sources and receivers; (b) highly smoothed velocity model and (c) depth migrated image after migration of the created virtual gathers. The last virtual source is located at the 4500 m position. The length of each window during the RC procedure is 50 s and the amount of overlap of two consecutive windows is 20%. Figure 10. Open in new tabDownload slide (a) Zoomed-in part of the velocity model (figure 8a) that includes all the virtual sources and receivers; (b) highly smoothed velocity model and (c) depth migrated image after migration of the created virtual gathers. The last virtual source is located at the 4500 m position. The length of each window during the RC procedure is 50 s and the amount of overlap of two consecutive windows is 20%. In the migrated section (figure 10c) some dipping horizons are imaged that match nicely with the original model. Raytracing through the model (figure 8c-d) shows that the most coverage in the area of the receiver profile is between positions 3000 and 4800 m. Also, the pre-stack depth migration result is obtained by stacking all individual migrated virtual gathers that cause low amplitude imaged horizons at the start and the end of the imaged profile due to stacking of only a few gathers. The last virtual source is located at the 4500 m position. The lack of propagated energy beyond 5000 m and the absence of further gathers cause the smiling and discrepancy between the imaged and true horizons. 7.1. Estimating a velocity vector from a virtual gather Figure 11 shows the result of semblance velocity analysis on a virtual common-midpoint (CMP) gather (from the Marmousi model) created using RC, with the CMP located at position 4000 m. In the semblance panel the red and yellow polygons show the P-multiples, reflected S and PS-waves, and non-physical reflections that are recognisable by higher velocity and earlier arrivals than the main reflections, respectively. Also, the black dashed line in the semblance panel shows the picked velocity used to apply the NMO correction. The results indicate that the amplitudes of reflected P-waves, and also PS and S waves (red polygon), are much stronger than the amplitudes of non-physical reflections (yellow polygon). Figure 11. Open in new tabDownload slide Semblance velocity analysis of the CMP gather where the CMP is located at 4000 m. The red polygon shows the region that includes the P-wave multiples and also PS and S reflected waves from deep horizons. The yellow polygon indicates the non-physical reflections with higher propagation velocity than main reflections. In the |${V_{{\rm{stk}}}}$| and |${V_{{\rm{int}}}}$| panels, the true stacking and interval velocities are shown in red and the estimated |${V_{{\rm{stk}}}}$| and |${V_{{\rm{int}}}}$| are shown in blue. Figure 11. Open in new tabDownload slide Semblance velocity analysis of the CMP gather where the CMP is located at 4000 m. The red polygon shows the region that includes the P-wave multiples and also PS and S reflected waves from deep horizons. The yellow polygon indicates the non-physical reflections with higher propagation velocity than main reflections. In the |${V_{{\rm{stk}}}}$| and |${V_{{\rm{int}}}}$| panels, the true stacking and interval velocities are shown in red and the estimated |${V_{{\rm{stk}}}}$| and |${V_{{\rm{int}}}}$| are shown in blue. 7.2. The effect of random noise on the RC results To examine the effect of random noise on the RC results, different levels of Gaussian noise were added to the simulated passive signals from the Marmousi model. The results in figure 12a and b show migrated sections from migration of extracted gathers created using RC from simulated signals with signal-to-noise (S/N) ratios equal to 10 and 1 dB. In figure 12c, the subtraction of migrated sections (figure 12b subtracted from figure 12a) is shown. The results indicate that the method can attenuate not only linear coherent events but also incoherent events (if S/N> = 1) and images the steeply dipping horizons. In this section, the velocity range chosen during the RC procedure is |$1300 \le v \le 4000$| m s–1. Figure 12. Open in new tabDownload slide Migrated sections from migration of created virtual gathers with S/N ratios of: (a) 10 dB; (b) 1 dB and (c) a and b. Figure 12. Open in new tabDownload slide Migrated sections from migration of created virtual gathers with S/N ratios of: (a) 10 dB; (b) 1 dB and (c) a and b. 8. Discussion To simulate passive signals, only sources distributed along the surface and far-field tests were considered. Furthermore, we investigate the special case where the ambient noise sources are clustered and located outside the seismic array (or profile). This is motivated by the situation where cultural noise may represent the main ambient source, such as passing trains, cars and so on (cultural noise sources) in the vicinity of, but outside, the seismic array. The results of correlations made by cross-correlating the reference trace with other traces indicate that the ground roll is the dominant event that masks the reflections. This linear event followed by some scattered surface waves creates a discontinuity in the reflection arrivals, and this can look like the fluctuation of reflection events due to static shifts. In addition, non-physical events due to the correlation of different arrivals (P, PS, PP, PSP, SSP, S) from the same/different interfaces make it difficult to correctly interpret the section for an unknown velocity model. Applying RC to the simulated signals retrieves continuous hyperbolic events that are free of this amplitude fluctuation that is seen in the cross-correlation results. During forward RC amplitudes are weighted by offset, which improves the localisation of the maximum energy correctly in the |$f - v\,\,$|map. The choice of the amplitude threshold plays an important role in eliminating the low cumulative amplitudes during the forward transform in order to avoid artificial events in the RC results. In the RC results, frequency distortion at shallow times occurs due to the stretching effect that requires the application of a Rho filter (Cao et al. 2003) to retrieve high frequencies at shallow time. Also, smearing of the amplitude during reverse RC (because of ground roll, diving waves and non-physical arrivals) is unavoidable. Our RC results from simulated passive signals showed that increasing the number of realisations improves the local focusing of energy and increases the maximum cumulative energy through a hyperbola. Increasing the number of realisations enables the use of large window lengths to improve the performance of the MTMC to estimate the coherence spectrum (here a 50s window length is used to correlate the simulated signals from the Marmousi model). To attenuate ground roll during the procedure, the minimum velocity in the velocity span should be higher than the ground-roll propagation velocity. It is possible to estimate the true velocity of linear coherent events from a correlated section made by general cross-correlation. Additionally, the amount of move out of the linear ground roll event causes low cumulative energy due to stacking of some but not all amplitudes through the event. This small cumulative energy in the |$f - v$| map causes the elimination of linear coherent events in the scaling step. We applied RC to noisy simulated signals with different S/N ratios, and showed that if the amplitude of the coherent signal is not lower than the random noise level, it is possible to highlight the hyperbolic event. Otherwise, at higher noise levels (S/N < 1), artificial events will appear. To be able to retrieve reflections from passive signals it is essential to have reflected energy (emitted by deep scattering points beneath the profile) at the surface at the location of the virtual source reflected back to subsurface media and bouncing from interfaces toward other receivers (Wapenaar et al. 2011). In this situation, the propagated energy suffers amplitude loss due to geometrical spreading, and therefore stacking of a large number of analysis windows from long recorded data sets followed by spherical divergence correction is required. This step helps to increase the S/N ratio of reflection arrivals so as to be easily detectable by RC. During correlation some non-physical reflections could appear at shallow times with high propagation velocity (King & Curtis 2012) due to the correlation of reflections from different interfaces. These spurious hyperbolic events appear in the RC results with lower amplitudes than the main reflections, and it is possible to use subjective RC to reject these events during the forward transform. By subjective RC we mean that the parameters are chosen subjectively by users during the forward RC procedure. Applying a very high iteration number followed by a small scaling factor is most likely to highlight the probable reflection arrivals free of non-physical events, however, this increases the time cost of computation. Trials showed that the scaling factor and iteration number are two important parameters during the retrieval of reflection arrivals. A high iteration number helps to achieve the sparsity during the inversion procedure, but increases the time cost. Using a very small scaling factor sacrifices very low amplitude reflection arrivals from deep horizons, while using a high scaling factor keeps all the amplitudes in the |$f - v$| map that results in the creation of many hyperbolic events highlighted in the time-offset domain during the inverse RC procedure. Those events include reflected PS, S, P, PSP, SS, PP and non-physical reflection events. The semblance velocity analysis of a virtual gather created by RC includes three zones. The first zone includes the reflected S, PS, PP, SP and SS arrivals that either propagate with a lower velocity than the main P-wave reflection arrival or arrive later than the main P-wave reflection arrivals (second zone). The third zone includes the non-physical reflection events due to the correlation of different arrivals (phases) from the same/different interfaces. These events appear at earlier time than the main primary events with higher propagation velocities. The main primary arrivals are located in the second zone and are easily recognisable from the other zones. The best model for retrieving the reflection arrivals is the layer-cake model with an increase of velocity with depth. In this case, the angle of reflection energy reaching the surface is small and enables the virtual source to transmit this energy to other virtual receivers. The reflection arrival retrieved from the three-layer model with heterogeneity showed a deviated reflection curvature at far offsets from the true reflection arrival due to the absence of transmitted energy from virtual source to far offset receivers. The increase of velocity due to lateral heterogeneity causes the decrease of velocity contrast at the first interface (results in low reflection amplitudes and very low multiple amplitudes), and also makes the reflection angles larger before reaching the virtual source locations. These cause the deviated reflection curvature at far offsets. To ameliorate the retrieved reflection arrival deviation at far offsets, it is essential to extract the virtual gathers for partial and not for full offsets. This method can be applied to real recorded passive seismic signals to extract the reflection arrivals. Before highlighting the arrivals from passive signals, some preparation must be undertaken such as beamforming, spectral whitening, geometrical spreading correction, elevation static correction, detrending and mitigation of possible spatial aliasing. 9. Conclusions In this paper, we introduced RC as a time-offset correlation to correlate passive signals and eliminate the ground roll based on the propagation velocity and the amount of move out from the correlated section. We also showed that it can highlight the hyperbolic events produced by reflected arrivals from subsurface interfaces. The results of performing RC on simulated passive signals indicate that the constructed hyperbolic events highlighted by applying RC show a good match with events from a synthetic-active-shot record created using the same velocity model. In addition, the migrated section obtained from the migration of virtual gathers created by RC from the simulated passive signals due to propagation of elastic waves through the Marmousi model shows a good match between the migrated P-wave reflections and the structures in the model. To mitigate against the creation of artificial events due to the correlation of low S/N ratio signals, it is vital to stack a large number of analysis windows using signals with a long recording time. Three steps help RC to highlight reflection arrivals with high certainty which are: (1) weighting amplitudes by offset during forward RC procedure (this is effective in magnifying the reflection arrival cumulative energies in the |$f - v$| map with propagation velocities in the velocity span range), (2) high iteration numbers to achieve sparsity in the model space during the forward RC procedure and (3) scaling helps to keep high cumulative localised energies and reject the low energies in the |$f - v$| map. A high scaling factor keeps all amplitudes in the forward RC procedure, while a very small scaling just keeps the highest energy in the |$f - v$| map, meaning that the low amplitude reflections can be rejected during the forward procedure. Without weighting amplitudes by offsets, different events with different moveouts appear with very close cumulative energies in the |$f - v$| map that make it difficult to segregate them by iteration and scaling. Data availability statement The data that support the findings of this study are available on request from the corresponding author. Acknowledgements We gratefully acknowledge PACIFIC (a European Research Project that has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 776622) for generously supporting this research work. Conflict of interest statement. Authors do not have any conflicts of interest to declare. References Aki K. , Larner K.L., 1970 . Surface motion of a layered medium having an irregular interface due to incident plane SH wave , Journal of Geophysical Research , 75 , 933 – 954 . Google Scholar Crossref Search ADS WorldCat Bakulin A. , Calvert R., 2006 . The virtual source method, theory and case study , Geophysics , 71 , SI139 – SI150 . Google Scholar Crossref Search ADS WorldCat Behr Y. , Townend J., Bowen M., Carter L., Gorman R., Brooks L., Bannister S., 2013 . Source directionality of ambient seismic noise inferred from three-component beamforming , Journal of Geophysical Research , 118 , 240 – 248 . Google Scholar OpenURL Placeholder Text WorldCat Boué P. , Poli P., Campillo M., Roux P., 2014 . Reverberations, coda waves and ambient noise: correlations at the global scale and retrieval of deep phases , Earth and Planetary Science Letters , 391 , 137 – 145 . Google Scholar Crossref Search ADS WorldCat Campillo M. , 2006 . Phase and correlation in random seismic fields and the reconstruction of the green function , Pure and Applied Geophysics , 163 , 475 – 502 . Google Scholar Crossref Search ADS WorldCat Cao Z. , Bancroft J.C., Brown R.J., Xaio C., 2003 . Radon transform and multiple attenuation , CREWES Research Report , Volume 15 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Claerbout J.F. , 1968 . Synthesis of a layered medium from its acoustic transmission response , Geophysics , 33 , 264 – 269 . Google Scholar Crossref Search ADS WorldCat Curtis A. , Gerstoft P., Sato H., Snieder R., Wapenaar K., 2006 . Seismic interferometry – turning noise into signal , The Leading Edge , 25 , 1082 – 1092 . Google Scholar Crossref Search ADS WorldCat Delsarte P.H. , Janssen A.J.E.M., Vries L.B., 1985 . Discrete prolate spheroidal wave functions and interpolation , SIAM Journal on Applied Mathematics , 45 , 641 – 650 . Google Scholar Crossref Search ADS WorldCat Draganov D. , Wapenaar K., Mulder W., Singer J., Verdel A., 2007 . Retrieval of reflections from seismic background-noise measurements , Geophysical Research Letters , 34 , L04305 . Google Scholar Crossref Search ADS WorldCat Forghani F. , Snieder R., 2010 . Underestimation of body waves and feasibility of surface-wave reconstruction by seismic interferometry , The Leading Edge , 29 , 790 – 794 . Google Scholar Crossref Search ADS WorldCat Gallot T. , Catheline S., Roux P., Campillo M., 2012 . A passive inverse filter for Green's function retrieval , The Journal of the Acoustical Society of America , 131 , EL21 – EL27 . Google Scholar Crossref Search ADS PubMed WorldCat Gholami A. , 2017 . Deconvolutive Radon transform , Geophysics , 82 , V117 – V125 . Google Scholar Crossref Search ADS WorldCat Gholami A. , Farshad M., 2019 . Fast hyperbolic Radon transform using chirp-z transform , Digital Signal Processing , 87 , 34 – 42 . Google Scholar Crossref Search ADS WorldCat Godin O.A. , 2006 . Recovering the acoustic Green's function from ambient noise cross correlation in an inhomogeneous medium , Physical Research Letter , 97 , 054301 . Google Scholar Crossref Search ADS WorldCat King S. , Curtis A., 2012 . Suppressing nonphysical reflections in Green's function estimates using source-receiver interferometry , Geophysics , 77 , Q15 – Q25 . Google Scholar Crossref Search ADS WorldCat Lobkis O.I. , Weaver R.L., 2001 . On the emergence of the Green's function in the correlations of a diffuse field , The Journal of the Acoustical Society of America , 110 , 3011 – 3017 . Google Scholar Crossref Search ADS WorldCat Löer K. , Riahi N., Saenger E. H., 2018 . Three-component ambient noise beamforming in the Parkfield area , Geophysical Journal International , 213 , 1478 – 1491 . Google Scholar Crossref Search ADS WorldCat Mokhtari A. , Gholami A., Siahkoohi H.R., 2019 . Fast hyperbolic deconvolutive Radon transform using generalized Fourier slice theorem , Geophysical Prospecting , 67, 408 – 422 . Google Scholar OpenURL Placeholder Text WorldCat Naghadeh D.H. , Morley C.K., 2016 . Propagation source wavelet phase extraction using multi-taper method coherence estimation , Journal of Geophysics and Engineering , 14 , 69 – 76 . Google Scholar Crossref Search ADS WorldCat Nakata N. , Snieder R., Tsuji T., Larner K., Matsuoka T., 2011 . Shear wave imaging from traffic noise using seismic interferometry by cross-coherence , Geophysics , 76 , SA97 – SA106 . Google Scholar Crossref Search ADS WorldCat Radon J. , 1917 . Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten Ber, über Verh. Königlich-Sächsischen Ges , Wiss. Leipzig , 69, 262 – 277 . Google Scholar OpenURL Placeholder Text WorldCat Shapiro N.M. , Campillo M., 2004 . Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise , Geophysical Research Letters , 31 , L07614 . Google Scholar Crossref Search ADS WorldCat Snieder R. , 2004 . Extracting the Green's function from the correlation of coda waves, a derivation based on stationary phase , Physical Review E , 69 , 046610 . Google Scholar Crossref Search ADS WorldCat Snieder R. , 2007 . Extracting the Green's function of attenuating heterogeneous acoustic media from uncorrelated waves , The Journal of the Acoustical Society of America , 121 , 2637 – 2643 . Google Scholar Crossref Search ADS PubMed WorldCat Thomson D.J. , 1982 . Spectrum estimation and harmonic analysis , Proceedings of the IEEE , 70 , 1055 – 1096 . Google Scholar Crossref Search ADS WorldCat Thorbecke J.W. , Draganov D., 2011 . Finite-difference modelling experiments for seismic interferometry , Geophysics , 76 , H1 – H18 . Google Scholar Crossref Search ADS WorldCat Van Tiggelen B.A. , 2003 . Green function retrieval and time reversal in a disordered world , Physical Review Letters , 91 , 243904 . Google Scholar Crossref Search ADS PubMed WorldCat Vasconcelos I. , Snieder R., 2008 . Interferometry by deconvolution, Part 2—theory for elastic waves and application to drill-bit seismic imaging , Geophysics , 73 , S129 – S141 . Google Scholar Crossref Search ADS WorldCat Versteeg R. , 1994 . The Marmousi experience, velocity model determination on a synthetic complex data set , The Leading Edge , 13 , 927 – 936 . Google Scholar Crossref Search ADS WorldCat Virieux J. , 1986 . P-SV wave propagation in heterogeneous media, velocity stress finite-difference method , Geophysics , 51 , 889 – 901 . Google Scholar Crossref Search ADS WorldCat Wapenaar K. , 2003 . Synthesis of an inhomogeneous medium from its acoustic transmission response , Geophysics , 68 , 1756 – 1759 . Google Scholar Crossref Search ADS WorldCat Wapenaar K. , van der Neut J., Ruigrok E., Draganov D., Hunziker J., Slob E., Thorbecke J., Snieder R., 2011 . Seismic interferometry by cross-correlation and by multi-dimensional deconvolution, a systematic comparison , Geophysics Journal International , 185 , 1335 – 1364 . Google Scholar Crossref Search ADS WorldCat Yilmaz Ö. , 1989 . Velocity-stack processing , Geophysical Prospecting , 37 , 357 – 382 . Google Scholar Crossref Search ADS WorldCat Yilmaz O. , 2001 . Seismic data analysis, processing, inversion and interpretation of seismic data , Society of Exploration Geophysicists , Investigations in Geophysics , Vol. 1 . Google Scholar OpenURL Placeholder Text WorldCat Yoon K. , Marfurt K.J., 2006 . Reverse-time migration using the Poynting vector , Exploration Geophysics , 37 , 102 – 107 . Google Scholar Crossref Search ADS WorldCat © The Author(s) 2021. Published by Oxford University Press on behalf of the Sinopec Geophysical Research Institute. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
TI - Retrieving reflection arrivals from passive seismic data using Radon correlation
JF - Journal of Geophysics and Engineering
DO - 10.1093/jge/gxab004
DA - 2021-03-03
UR - https://www.deepdyve.com/lp/oxford-university-press/retrieving-reflection-arrivals-from-passive-seismic-data-using-radon-bIPuo8Ifwc
SP - 1
EP - 15
VL - 18
IS - 2
DP - DeepDyve
ER -