Seismic groundroll prediction by interferometry and separation in curvelet domain by Jiupeng Yan B.Sc, Peking University, Beijing, China A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Geophysics) THE UNIVERSITY OF BRITISH COLUMBIA May, 2011 c Jiupeng Yan, 2011 Abstract Groundroll is a type of surface wave that propagates along the Earth’s surface. Groundroll usually has low frequency, low velocity and high amplitude. Due to its high amplitude, groundroll almost always dominates reflected body waves in land seismic data and covers important reflection information. Therefore, removing groundroll noise is a very important step before seismic imaging. The most common methods used in industry to remove groundroll are the Fourier domain filtering methods based on the different characteristics of groundroll and reflections, i.e. the low frequency and low velocity properties of groundroll. However, groundroll and reflection usually have large overlap in both physical and frequency domain. Also groundroll is spatially aliased at normal receiver intervals causing additional processing difficulties. Therefore, a good separation of groundroll by Fourier domain filtering method is challenging. In this thesis, we propose a data-driven workflow to remove groundroll. Our workflow is motivated both by SRME (Surface Related Multiple Elimination) method and a recently proposed interferometry method for the prediction of groundroll. It consists of a prediction step based on interferometry and a robust separation step that involves curvelet domain matched filtering and sparsity promotion. Tests of our workflow on synthetic data show clear removal of large amplitude groundroll and preservation of seismic reflection events. Test of our separation step on real data shows improvement over conventional Fourier domain filtering methods. ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi . . . . . . . . . . . . . . . . . . . . . . . . . . . ix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Abstract List of Tables Acknowledgments Dedication 1 Introduction . . . . . . . . . . . . . 1.1 Surface waves - groundroll . . . 1.2 Thesis overview . . . . . . . . . 1.3 Theoretical background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 4 6 2 Groundroll prediction . . . . . . . . . . . . . . . . . . . . . 2.1 Seismic interferometry and stationary phase analysis . . . 2.2 Groundroll prediction by cross-correlation interferometry 2.3 Data matrix formulation of interferometry . . . . . . . . . . . . . . . . 10 10 19 23 3 Groundroll separation by curvelet-domain matched filtering 3.1 Least-square adaptive subtraction strategies . . . . . . . . . 3.2 Curvelet-domain matched filtering . . . . . . . . . . . . . . . 3.3 Synthetic examples and discussions . . . . . . . . . . . . . . 28 29 31 33 4 Bayesian separation for groundroll . . . . . . . . . . . . . . . 4.1 Bayesian wavefield separation algorithm . . . . . . . . . . . . 4.2 Application to synthetic data . . . . . . . . . . . . . . . . . . 38 38 42 iii Table of Contents 5 Conclusion and future research . . . . . . . . . . . . . . . . 5.1 A data-driven workflow for groundroll removal . . . . . . . . 5.2 Open and future research . . . . . . . . . . . . . . . . . . . . 45 45 46 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Appendix A A real example for curvelet matching and Bayesian separation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 iv List of Tables 3.1 SNR of the reflections produced by subtracting curvelet matched groundroll for different parameter selections. . . . . . . . . . 34 4.1 The Bayesian-separation algorithm. Table from Yarham and Herrmann (2008) . . . . . . . . . . . . . . . . . . . . . . . . . SNR of the Bayesian separated reflections for different parameter selections. . . . . . . . . . . . . . . . . . . . . . . . . 4.2 42 44 v List of Figures 1.1 1.2 1.3 1.4 1.5 1.6 2.1 2.2 2.3 Particle motion of Rayleigh waves. . . . . . . . . . . . . . . . A shot gather of land seismic data. . . . . . . . . . . . . . . . A comparison of our workflow for interferometry groundroll removal and conventional SRME workflow. . . . . . . . . . . Six real curvelets examples of different scales and angles in the physical (left) and frequency (right) domain. Curvelets are strictly localized in the Fourier domain and have rapid decay in the physical domain. Image from Herrmann and Hennenfent (2008). . . . . . . . . . . . . . . . . . . . . . . . Partitioning for 2D discrete curvelet transform. The 2D Fourier plane is partitioned into second dyadic coronae and the coronae are sub-partitioned into angular wedges. Image from Hennenfent et al. (2010). . . . . . . . . . . . . . . . . . . . . . . Example of correlation of curvelets with curved wavefronts. Image from Herrmann and Hennenfent (2008). . . . . . . . . (a) A Ricker wavelet with 20 Hz central frequency at 0.5s (b) (a) shifted for 1.5s (c) correlation of (a) and (b). . . . . . . . Simple acoustic model for testing seismic interferometry. Letters A and B indicate the two receiver positions. Black dots indicate the source locations numbered from 1 to 200. Acoustic wave velocity of the top layer is 2.0 km/sec and the bottom layer is 4.0 km/sec. . . . . . . . . . . . . . . . . . . . . . . . . (a) Responses at receiver A. (b) Responses at receiver B from different serial sources at the circle boundary. . . . . . . . . . 3 3 5 8 9 9 12 13 15 vi List of Figures 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 2.13 2.14 3.1 3.2 Ray path analysis for acoustic interferometry example in Figure 2.2. Sources locate around green circle are stationary sources for the direct wave. Correlation of traces at A and B from sources around green circle eliminates the common ray path and reveals the difference, i.e. raypath for the direct wave. Sources locate around purple circle are stationary sources for the reflected wave. Correlation of traces at A and B from sources around purple circle eliminates the common ray path and reveals the difference, i.e. raypath for the reflected wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Correlated traces for different serial sources (b) Accumulative sum of correlated traces from different portions of all the sources. . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Actual response at A from a source at B by finite difference (b) Response at A from B constructed by interferometry . . P-wave velocity model for generating (a) reflections and (b) groundroll. . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) a synthetic shot gather at 1.0 km (b) true reflections contained in (a) . . . . . . . . . . . . . . . . . . . . . . . . . . . . Correlated traces at A and B with auto gain control from difference source locations at the surface. . . . . . . . . . . . (a) Example of 2D line data cube in time domain (b) Example of 2D line data cube in frequency domain. . . . . . . . . . . . Frequency slice of data volume. . . . . . . . . . . . . . . . . . Prediction of groundroll for a monochromatic matrix. . . . . Example of R matrix. . . . . . . . . . . . . . . . . . . . . . . (a) Interferometry prediction of the groundroll for the shot gather at 1.0 km. (b) True groundroll for the shot gather at 1.0km. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Global Fourier matched groundroll for shot gather at 1.0 km (b) reflection produced by subtracting (a) from the total data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Curvelet matched groundroll with γ = 0.02 for shot gather at 1.0 km. (b) Reflection produced by subtracting (a) from original data. (c) difference between (b) and true reflection . 16 17 18 21 22 23 24 25 25 26 27 34 35 vii List of Figures 3.3 3.4 4.1 (a) Comparison of reflection produced by global Fourier matching with the true reflection at offset 0.4km (b) Comparison of reflection produced by curvelet-domain matching with the true reflection at offset 0.4km. Shot location is at 1.0km. Red lines represent the true reflection and blue lines represent reflection results produced by subtracting global Fourier matched groundroll and curvelet matched groundroll. . . . . . (a) Mosaic and (b) 2D plots of the curvelet domain filter for shot gather at 1.0 km and γ = 0.02. In (a), the low frequency components (i.e. coarse scales) of the curvelet domain filter are in the center and high frequency components (i.e. fine scales) are at the perimeter. The angle increases clockwise and the number of angles double every other scale. The color represents the amplitudes of the curvelet coefficients. . . . . . (a) Groundroll produced by Bayesian separation (b) Reflection produced by Bayesian separation for λ1 = 2.0, λ2 = 8.0, and η = 2.0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1 (a) A real shot gather. (b) groundroll prediction produced by Fourier domain threshold and filtering. . . . . . . . . . . . . . A.2 (a) Reflection produced by direct subtraction of the Fourier domain groundroll prediction. (b) Separated reflection through curvelet domain matched filtering and Bayesian separation. . 36 37 43 51 52 viii Acknowledgments First of all, I would like to thank my supervisor Felix Herrmann. Without his support, encouragement and creative ideas, none of this could be possible. I would like to thank all the team members of SLIM group for their help and encouragement. I am very fortunate to be part of this group and to witness the amazing development that it has gone through. I am very fortunate to have Michael Bostock and Doug Oldenburg to be part of my supervisory committee. I have learned a great deal from their insightful courses at UBC. The knowledge and inspiration they provided were invaluable. I would like to thank Christian Schoof for reading through my thesis and being part of my defense committee. Finally, I would like to thank my family for their support and encouragement. Thank you all. ix To my parents x Chapter 1 Introduction Seismic surveying is the most common method for hydrocarbon exploration both on land and in the ocean. The basic principle for a survey is to send seismic waves into the Earth using sources such as a seismic vibrator on land or an air-gun at sea. The artificial seismic waves sent out during a survey are reflected by layer structures in the Earth and are recorded using hydrophones in water or geophones on land. The recorded data is processed to produce images of the subsurface structures, which will be interpreted by geoscientists to determine the potential structural traps that may contain hydrocarbons. During acquisition of seismic data, both signal and unwanted energy are recorded. Such undesired energy, more commonly known as noise, contaminates seismic data, leading to lower quality images and misinterpretations. Thus, enhancing the reflection signal by suppressing noise is one of the main goals in seismic data processing. Seismic noise is generally classified into two categories: incoherent and coherent noise (Yilmaz, 2001). Incoherent noise is uncorrelated between adjacent traces and can often be attenuated by stacking traces and by filtering during processing. Coherent noise is unwanted seismic energy that is consistent trace by trace, which makes it more difficult to remove. The most common types of coherent noise include: multiples, air waves, and surface waves (Yilmaz, 2001). This thesis deals with separating seismic signal from coherent groundroll noise. It focuses on removing surface waves in land seismic data. 1.1 Surface waves - groundroll Surface waves are waves propagating along the surface of the Earth. The amplitudes of surface waves decay exponentially with depth but slowly along the surface. Therefore surface waves have large amplitudes. The most typical surface waves are Love waves and Rayleigh waves. Love waves are named after A.E.H. Love, a mathematician who first predicted the existence of these waves in 1911. Love waves are horizontally polarized shear waves 1 Chapter 1. Introduction (SH waves) trapped in a low-velocity elastic layer over an elastic half-space. The particle motion of a Love wave is parallel to the surface and perpendicular to the direction of propagation. Love waves are generated by the constructive interference of reverberations. In common seismic exploration, Love waves are inconsequential because the seismic sources used in exploration do not generate large amount of Love waves. Common geophones that respond to vertical ground motion cannot detect Love waves, which have horizontal ground motion (Telford et al., 1990). In seismic exploration, the more important surface waves are Rayleigh waves because they commonly exist in land seismic data. Rayleigh waves were first predicted by British physicist John Willisam Strutt (Lord Rayleigh) in 1885. Rayleigh waves have both vertical and horizontal particle motions in the plane normal to the surface and parallel to the direction of propagation. The vertical and horizontal motions have a definite phase relation and thus produce overall elliptical motions in the vertical plane parallel to the direction of propagation. Due to the elliptical particle motion, Rayleigh waves ’roll’ along the ground, which is why they are commonly known as groundroll. Figure 1.1 shows a simple schematic diagram of the particle motion of groundroll. Rayleigh waves are generated by the constructive interference of reflected and transmitted body waves (both P and SV waves) at the free surface. Their particle motions decay exponentially with depth but slowly along the surface. Therefore Rayleigh waves live for a long time in seismic record. The depth of significant particle motion is at the same scale of the wavelength. Rayleigh wave velocity depends on the near surface elastic parameters and is less than the S wave velocity. Because the depth of significant motion of a Rayleigh wave depends on the wavelength, for media that have densities or velocities changing with depth, the velocity of Rayleigh waves changes with wavelength (or frequency) and Rayleigh waves become dispersive. For theoretical analysis and more detailed discussion on characteristics of Rayleigh waves, we refer the reader to Aki and Richards (2002). To summarize, Rayleigh waves travel along the Earth surface, are low frequency, low group velocity, dispersive and have elliptical particle motions that decay exponentially with depth. Groundroll is the most common coherent noise source in land seismic data. As a result of its large amplitude, groundroll almost always dominates seismic reflection events that may exist in the recorded data. Figure 1.2 shows a land shot gather contaminated by groundroll. The high amplitude groundroll corrupts the data, causing the primary reflections to be practically invisible. This is one of the major reasons for low signal to noise ratio (SNR) of land seismic data. 2 Chapter 1. Introduction Propagation direction of Rayleigh wave Figure 1.1: Particle motion of Rayleigh waves. Figure 1.2: A shot gather of land seismic data. 3 Chapter 1. Introduction The most common methods to remove groundroll are based on different characteristics of groundroll and reflected body waves, which we refer to as reflections from this point on. The low group velocity and low frequency characteristics of groundroll are most utilized to separate them from higher frequency reflections. The most popular methods used for groundroll removal include band-pass frequency filtering (Yilmaz, 2001), f − k domain filtering, multi-channel filtering (Galbraith and Wiggins, 1968), and other Fourier based filtering with different variants including wide-band velocity filtering (Embree et al., 1963) and fan filtering (Fail and Grau, 1963; Treitel et al., 1967). The general idea of these methods is to transfer the data into the Fourier domain, remove characteristic frequencies of groundroll and then transfer the data back to the physical domain. However, a good separation by the Fourier domain filtering methods is challenging because groundroll and reflection usually have significant overlaps both in the physical domain and in the Fourier domain. Also groundroll is spatially aliased at typical receiver intervals causing additional processing problems. Several papers have shown that these methods have a tradeoff of distorting reflections (McMechan and Sun, 1991; Liu, 1999). To address these issues, other methods beyond the Fourier domain have been developed to remove groundroll and these include Wiener-Levinson algorithm (Karslı and Bayrak, 2004), Radon transform (Russel et al., 1990; Trad et al., 2003), and one dimensional or two dimensional wavelet transform (Deighan and Watts, 1997; Zhang and Ulrych, 2003). Unfortunately, challenges still remain, particularly when dealing with aliased groundroll with poor frequency domain separation. Moreover, because the velocity and frequency of groundroll depends on near-surface velocity structure and density structure, most of these filtering methods require user interpretation and careful selection of parameters for different datasets. Therefore, a data-driven method that avoids selection of sensitive parameters and user interpretation for different datasets is valuable. 1.2 Thesis overview The main theme of this thesis is a data-driven workflow to remove groundroll. The workflow is motivated both by Surface Related Multiple Elimination method (SRME- Verschuur et al., 1992) and recently developed seismic interferometry method for the prediction of groundroll (Curtis et al., 2006; Halliday et al., 2007; Dong et al., 2006; Vasconcelos et al., 2008). 4 Chapter 1. Introduction Instead of using filtering methods which are based on different characteristics of groundroll and reflection, our workflow is based on a data-driven prediction of groundroll by seismic interferometry, coupled with a robust subtraction step that involves curvelet-domain matched filtering (Herrmann et al., 2008) and sparsity promotion (Yarham and Herrmann, 2008; Wang et al., 2008). During the prediction, we exclude reflected body waves by muting the corresponding stationary sources in seismic interferometry (Vasconcelos et al., 2008). Since the data itself has been used as a prediction operator, our method does not require prior information. In the separation step, the matching process forces the predicted groundroll to fit with the true groundroll in the input data. The underlying assumption for matching is that data after subtraction has minimum energy. Finally, we utilize a Bayesian algorithm developed by Yarham and Herrmann (2008) to replace direct subtraction and thus enhance the overall separation of groundroll and reflection. Figure 1.3 shows a comparison of our workflow for interferometry groundroll removal and conventional SRME workflow. Our Workflow Conventional SRME Workflow input prediction (multidimensional correlation) input prediction (multidimensional convolution) conservative Fourier matching conservative Fourier matching curvelet-domain matching Windowed Fourier matching Bayesian separation Direct subtraction Figure 1.3: A comparison of our workflow for interferometry groundroll removal and conventional SRME workflow. Fourier domain matching is most utilized for adaptive subtraction. Despite the major advances by Fourier domain filtering methods, it is also known that prediction errors that vary as a function of offset, time and 5 Chapter 1. Introduction dip pose serious challenges. To achieve better results than conventional Fourier matching techniques, we use the curvelet transform developed by Cand`es and Donoho (2004). The curvelet transform is a data-independent, multiscale and multidirectional transform. Curvelets are localized in the frequency domain and have rapid decay in the physical domain. In the physical domain, curvelets look like localized plane waves. The curvelet transform has been utilized to approximate piece-wise smooth images with curved edges, e.g., complex wavefronts, with few significant curvelet coefficients (Cand´es et al., 2000; Herrmann and Hennenfent, 2008). In this thesis, we utilize curvelets to build a multidirectional and multiscale filter to match predictions to true groundroll. We also utilize the sparsity provided by curvelets to improve groundroll separation. In the last section of this chapter, we discuss the properties of curvelets in more detail as a theoretical background of the thesis. In Chapter 2, we give an introduction on seismic interferometry. Next we discuss how to predict surface waves using seismic interferometry and the problems associated with interferometry prediction of surface waves. In Chapter 3 and 4, we discuss the curvelet domain matching technique and Bayesian separation algorithm as mentioned earlier. Finally, in Chapter 5, we summarize our workflow. We give conclusions and discuss possibilities for future research. 1.3 Theoretical background As mentioned earlier, we use curvelets in the separation step of our method for matched filtering and sparsity promotion. Before we present our workflow, we first give a general introduction to the curvelet transform and its advantages for representing seismic data over other transforms. Real seismic data contains wavefronts coming from different directions with various frequencies and may have conflicting dips (Hennenfent and Herrmann, 2008). Finding a transform that can detect and optimally represent (i.e. using only a few significant coefficients) these complex wavefronts is highly sought after by geophysicists. In seismology, one of the mostly utilized transforms is the Fourier transform, which decomposes the signal into monochromatic plane waves with constant frequency and dip. The Fourier transform is good at representing periodical signals that are sums of simple sinusoidal waves, but it is not efficient at representing signals with discontinuities, i.e. it requires many plane waves to represent a curved wavefront in seismic data. One solution to overcome the difficulty of representing signals 6 Chapter 1. Introduction with discontinuities is to use the wavelet transform. The wavelet transform decomposes the signal into a set of basis functions called wavelets. The basic wavelet in the wavelet transform is a normalized function with zero average, which is centered and well-concentrated in the neighborhood of time zero. A set of time-frequency atoms is obtained by dilation and translation of the basic wavelet function. The wavelet transform is better at representing nonperiodic signals that have discontinuities compared to the Fourier transform. However, the basis functions of 2D/3D wavelet transform are isotropic thus in high dimensions the wavelet transform is non-directional. 2D/3D wavelet transform decompose the signal into ’fat’ isotropic wave packets. Since seismic data contains curved wavefronts, many wavelet coefficients are required to properly represent seismic data. To resolve the issues of representing wave-front like features, curvelets have been utilized by Cand´es et al. (2000). Examples of curvelets in the physical domain and in the frequency domain are given in Figure 1.4. In the frequency domain, curvelets are strictly localized. In the physical domain, curvelets have fast decay and are smooth in one direction and oscillatory in the perpendicular direction. Figure 1.5 shows the partitioning for 2D discrete curvelet transform. The 2D Fourier plane is partitioned into different scales and the different scales are sub-partitioned into angular wedges. The number of angular wedges double at every second scale (Cand´es et al., 2000). Curvelets have the parabolic scaling: width ∝ length2 , which makes them act as localized eigenfunctions. Such construction gives curvelets following properties: multiscale, multidirectional and anisotropic. With these properties, the curvelet transform becomes an optimal choice for representing seismic data (Herrmann and Hennenfent, 2008). Curvelets can handle the multiscale and multidirectional properties of seismic data. Curvelets can detect curved wavefronts due to their special shape, as shown in Figure 1.6. Moreover, because of its specific partitioning, the curvelet transform is able to differentiate between different signal components on the basis of location, angle and frequency, making it suitable for separation of seismic signals. Curvelets also provide the sparsest representation of complex seismic data (Herrmann and Hennenfent, 2008; Cand`es et al., 2006). We refer the reader to Herrmann and Hennenfent (2008) for a comparison of the sparsity of seismic data obtained by using the Fourier transform, the wavelet transform and the curvelet transform. To summarize, we utilize the curvelet transform in our workflow to improve the separation of groundroll for the two major reasons: • The curvelet transform is a localized, multiscale and multidirectional 7 Chapter 1. Introduction transform, which decomposes signals into different location, angle and frequency. Since groundroll and reflection have different frequencies and velocities, they are more separable in the curvelet domain than the Fourier domain or the wavelet domain; • The sparsity provided by curvelets increases the probability of separation (sparse signals have more chances to be separated) and allows us to utilize effective signal separation algorithms through l1 minimization (Yarham and Herrmann, 2008; Wang et al., 2008). -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 Figure 1.4: Six real curvelets examples of different scales and angles in the physical (left) and frequency (right) domain. Curvelets are strictly localized in the Fourier domain and have rapid decay in the physical domain. Image from Herrmann and Hennenfent (2008). 8 Chapter 1. Introduction wavenumber FINEST SCALE FINEST SCALE frequency FINEST SCALE 2j 2j/2 FINEST SCALE Offset (m) -2000 0 Figure 1.5: Partitioning for 2D discrete curvelet transform. The 2D Fourier Complete curvelet tilling of the f-k domain 0 One scale One angular wedge plane is partitioned into second dyadic coronae and the coronae are subpartitioned into angular wedges. Image from Hennenfent et al. (2010). 0 -2000 Offset (m) 0 2000 0.5 Time (s) 0.5 1.0 1.5 Significant curvelet coefficient Curvelet coefficient~0 2.0 Time (s) Figure 1.6: Example of correlation of curvelets with curved wavefronts. Image from Herrmann and Hennenfent (2008). 1.0 9 2000 Chapter 2 Groundroll prediction In the previous chapter, the concept of groundroll was introduced. Groundroll is a type of coherent noise generated by Rayleigh waves, which have low velocity, low frequency and high amplitudes. Groundroll almost always dominates reflection energy in land seismic data and degrades overall data quality due to its high amplitude. The most common methods to remove groundroll are the Fourier domain filtering methods, which utilize the low velocity and low frequency characteristics of groundroll. However, groundroll and reflection usually have significant overlap both in the physical domain and the frequency domain. Also, groundroll is spatially aliased at normal receiver intervals. Thus, it is challenging to achieve a good separation for groundroll and reflection by Fourier domain filtering methods alone. Moreover, variations of near-surface velocity structures cause groundroll to change in magnitude, velocity and frequency, requiring user interpretation for different datasets. Our method to remove groundroll is based on a datadriven prediction of groundroll by seismic interferometry coupled with a separation step that involves curvelet-domain matching and sparsity promotion. Since our method does not require prior information about near-surface velocity structure, it avoids user interpretation and sensitive parameter selection. In this chapter, we will discuss groundroll prediction method by seismic interferometry, which has recently been introduced by Curtis et al. (2006); Halliday et al. (2007); Dong et al. (2006); Vasconcelos et al. (2008). We first introduce the concept of seismic interferometry and then explain how to predict surface waves using interferometry. 2.1 Seismic interferometry and stationary phase analysis In exploration geophysics, seismic interferometry is beneficial since it can extract inter-receiver responses from seismic data generated by active exploration sources or passive random noise sources. The term interferometry was first used in radio astronomy, which means extracting information by cor10 Chapter 2. Groundroll prediction relating radio signals from distant objects. Claerbout (1968) proposed the idea of daylight imaging, a method that mimics an active source/receiver response from correlating passive random noise traces recorded at two different locations. Schuster et al. (2004) applied the correlation method to active sources and proposed interferometry imaging. More recently, Wapenaar and Fokkema (2006) developed a theory based on seismic representation theorem and formally solidified the relationship in seismic interferometry for both acoustic and elastic media. This representation theory forms the theoretical basis for seismic interferometry using passive or active sources. Following Wapenaar and Fokkema (2006), we use monochromatic Green’s function to express the cross-correlation interferometry in acoustic media: 2 {G(ra , rb , ω)} ≈ 2 ρc G(ra , s, ω)∗ G(rb , s, ω)d2 s. (2.1) ∂V In this expression, all the quantities are in the Fourier domain. The vectors ra and rb represent the locations of two receivers A and B inside the domain V . The boundary of the domain V is denoted by ∂V . The acoustic wave speed of the medium is denoted by c and the density is denoted by ρ. The symbol expresses the operation of taking the real part of a complex value and ∗ denotes complex conjugate. For an impulsive source of volume injection rate at rb , Green’s function G(ra , rb , ω) represents observation of the pressure response at A for a single frequency ω. The correlation of two signals in the time domain is equivalent to the element wise product of the complex conjugate of the first signal and the second signal in the frequency domain (i.e. f g = F −1 (F(f ))∗ (F(g)), where f g denotes the correlation of two signals f and g). Since Equation 2.1 is in the frequency domain, the product G(ra , s, ω)∗ G(rb , s, ω) in Equation 2.1 represents the correlation of the two Green’s functions in the time domain. For two continuous functions f and g, their cross-correlation in time domain is defined as: (f g)(t) = ∞ f ∗ (τ )g(t + τ )dτ (2.2) −∞ where * denotes complex conjugate. According to Equation 2.2, crosscorrelations measures the similarity of two waveforms as a function of a time-lag applied to the second waveform (i.e. it measures the time alignment of two signals). Suppose there are two signals a and b, where b is a with a time shift δ. In that case their correlation a b has maximum value at time δ. Figure 2.1 shows a simple example of correlation for two signals. 11 Chapter 2. Groundroll prediction Figure 2.1(a) is a plot of a 20Hz central frequency Ricker wavelet at 0.5s. Figure 2.1(b) has the same wavelet as Figure 2.1(a) but it has been shifted from 0.5s to 2.0s. The correlation of the these two signals (Figure 2.1(c)) shows a wavelet at 1.5s. In the seismic context, two traces recorded at different receivers from the same source give us two signals. The difference between the two signals that share the common source depends on the ray paths from the common source to the two receivers. Cross-correlation of the two traces eliminates the ray path that seismic waves have in common and reveals the difference. (a) (b) (c) Figure 2.1: (a) A Ricker wavelet with 20 Hz central frequency at 0.5s (b) (a) shifted for 1.5s (c) correlation of (a) and (b). Equation 2.1 allows us to approximate the response at A from a source at B by integrating the correlated responses at A and B over sources at 12 Chapter 2. Groundroll prediction ∂V . To derive Equation 2.1, several assumptions and approximations are made including the condition that the medium is non-attenuating inside domain V and homogeneous at and outside the boundary ∂V . For more detailed discussion on approximations and assumptions needed to derive Equation 2.1, we refer the reader to Wapenaar and Fokkema (2006). 100 125 75 A B 150 50 175 25 01 Figure 2.2: Simple acoustic model for testing seismic interferometry. Letters A and B indicate the two receiver positions. Black dots indicate the source locations numbered from 1 to 200. Acoustic wave velocity of the top layer is 2.0 km/sec and the bottom layer is 4.0 km/sec. To explain how interferometry works, we consider an acoustic example with two horizontal layers. The top layer has a velocity of 2.0 km/sec and the bottom layer has a velocity of 4.0 km/sec. Figure 2.2 includes the model, where the positions of two receivers are indicated by letter A and B. The black dots indicate the 200 source positions at a circle boundary of 1km radius. Source 1 is located in the middle at the bottom of the circle and the remaining sources distribute along the circle in an anti-clockwise fashion. 13 Chapter 2. Groundroll prediction We simulate the scalar wavefield by solving the acoustic wave-equation: ∆U − 1 ∂2U = f (t). c(x)2 ∂t2 (2.3) In this expression, ∆ is the Laplacian operator, f (t) is the source wavelet, c(x) is the acoustic wave velocity which is a variable of spatial postion, and U is a scalar wavefield. We use a fourth-order finite difference scheme for the Laplacian operator in space and second order scheme for time stepping. The time step is set to 0.25 ms and the spatial grid is set to 5 m for both spatial directions. We use a 20 Hz central frequency Ricker wavelet as the source wavelet. Figure 2.3 shows the responses at A and B from the 200 different sources at the circle boundary. Since receivers A and B are enclosed by the sources located at the circle boundary, we can use Equation 2.1 to reconstruct the response at A from a source at B. Figure 2.5(a) shows correlation of traces at A and B from different sources at the boundary without stack. Figure 2.5(b) shows the stack of the correlated traces in Figure 2.5(a) from different portions of the 200 sources. For example, trace 100 in Figure 2.5(b) is the sum of correlated traces from source 1 to source 100 in Figure 2.5(a). Figure 2.6(b) is the stack of correlated traces from all the 200 sources, i.e. the last trace in Figure 2.5(b). Figure 2.6(a) is the response at A from a Ricker wavelet source at B produced by finite difference modeling. For this simple two layer acoustic model, there are two main events: the direct waves (indicated by green color) and reflected waves (indicated by purple color) , which are clearly visbile in Figure 2.6(a). The direct arrival occurs around 0.5s and the reflected wave occurs around 1.2s. By comparing the stacked trace in Figure 2.6(b) with the actual response in Figure 2.6(a), we find that they are approximately equal to each other as predicted by Equation 2.1. However, the response by interferometry (Figure 2.6(b)) is not the same as the original response by finite-difference modeling (Figure 2.6(a)) for several reasons. First, Equation 2.1 applies to Green’s function and not to data that carries the imprint of a source wavelet (i.e. the 20Hz central frequency Ricker wavelet). Second, Equation 2.1 is based on several assumptions and approximations including the condition that the medium is homogeneous at the boundary of all the sources. However, in our simulation, the circle boundary of sources is heterogeneous due to the interface of the two layers. The sources located at the interface that indicated by the small red circle in Figure 2.4 cause the non-physical arrival before the direct wave in Figure 2.6(b). In addition, there are other approximations that we discussed earlier in this chapter. 14 Chapter 2. Groundroll prediction (a) (b) Figure 2.3: (a) Responses at receiver A. (b) Responses at receiver B from different serial sources at the circle boundary. 15 Chapter 2. Groundroll prediction 100 A B 150 50 10 Figure 2.4: Ray path analysis for acoustic interferometry example in Figure 2.2. Sources locate around green circle are stationary sources for the direct wave. Correlation of traces at A and B from sources around green circle eliminates the common ray path and reveals the difference, i.e. raypath for the direct wave. Sources locate around purple circle are stationary sources for the reflected wave. Correlation of traces at A and B from sources around purple circle eliminates the common ray path and reveals the difference, i.e. raypath for the reflected wave. 16 Chapter 2. Groundroll prediction (a) (b) Figure 2.5: (a) Correlated traces for different serial sources (b) Accumulative sum of correlated traces from different portions of all the sources. 17 Chapter 2. Groundroll prediction (a) (b) Figure 2.6: (a) Actual response at A from a source at B by finite difference (b) Response at A from B constructed by interferometry Figure 2.5(a) and Figure 2.5(b) represents the contributions of different sources in the response by interferometry. Figure 2.5(a) shows the contribution of a single source while Figure 2.5(b) shows the cumulative contribution from source 1 to the corresponding source number at the horizontal axis. Both Figure 2.5(a) and Figure 2.5(b) show that different sources contribute differently in the response by interferometry. We refer to sources that contribute constructively to the interferometry response as stationary sources. In our simple experiment, there are two major kinds of stationary sources contributing respectively to the direct and reflected waves. Other sources that sum up destructively in the interferometry response are nonstationary sources. As mentioned earlier, the interface on the boundary of the sources breaks the homogeneous boundary condition which is needed to derive Equation 2.1 and causes the nonphysical arrival before direct wave in Figure 2.6(b)). The red circle in Figure 2.4 indicates the sources at the interface while the red arrows in Figure 2.5(a) and Figure 2.5(b) show their effect on the interferometry response. Sources located at approximately No.70 (indicated by green circle in Figure 2.4 and green arrows in both Figure 2.5(a) and Figure 2.5(b)) are the major contributors to the direct wave that arrives at 0.5s roughly. In Figure 2.4, sources around the small green circle are stationary sources for 18 Chapter 2. Groundroll prediction the direct wave. Correlation of traces at A and B from these sources eliminates the ray path they have in common (dashed green line) and reveals the difference (solid green line), i.e. the direct wave traveling from B to A. The green arrows in Figure 2.5(a) and Figure 2.5(b) indicate that these sources contribute constructively for the direct wave, as predicted by the ray path analysis. Similarly, sources around the small purple circle are stationary sources for the reflected wave. Correlation of traces at A and B from these sources eliminates the ray path they have in common (dashed purple line) and reveals the difference (solid purple line), i.e. reflected wave traveling from B to A. The purple arrows in Figure 2.5(a) and Figure 2.5(b) indicate that sources around the purple circle contribute constructively for the reflected wave. 2.2 Groundroll prediction by cross-correlation interferometry In the previous section, we discussed seismic interferometry for acoustic media and performed the stationary phase analysis for a simple two layer example (Figure 2.2). The stationary phase analysis plays a key role in the interferometry prediction for groundroll. By performing a similar stationary phase analysis for realistic land seismic data, we show that all the sources located at surface act as stationary sources for surface waves. We start our analysis by introducing the interferometry for elastic media. In reality land seismic data is the response of Earth’s elastic material from seismic sources such as dynamite or a seismic vibrator. The observation of geophones is not the pressure of a simple scalar wavefield but particle velocity in different directions. Seismic interferometry theory for elastic media is more complex than acoustic media. Our objective is using seismic interferometry to predict groundroll and it only exists in elastic media. To predict groundroll using seismic interferometry, we must consider seismic interferometry theory for elastic media. Wapenaar and Fokkema (2006) show a similar representation of seismic interferometry for elastic media using Green’s functions: 2 {Gv,f p,q (xa , xb )} ≈ On left side, Gv,f p,q (ra , rb ) 2 ρcK ∗ v,Φ 2 Gv,Φ p,K (xa , x)Gq,K (xb , x) d x (2.4) ∂V is the Green’s function observed at ra for an im19 Chapter 2. Groundroll prediction pulsive source at rb in the frequency domain. The superscript v indicates that the measured quantity is particle velocity and f indicates the source is a point source of force. Subscript p indicates the direction of particle velocity and q indicates direction of point force source. On the right side of Equation 2.4, Gv,Φ p,K (ra , s) also denotes the Green’s function. The difference between the Green’s function on left side and right side is that subscript Φ denotes sources of potential and K denotes the different potential source types, i.e. P-wave source (K = 0) and S-wave sources with different polarization directions (K = 1, 2, 3). cK is the velocity of different waves and ρ is the density of the media. The repeated subscript K denotes Einstein summation from 0 to 3, i.e. summation over different source types. From a number of numerical comparisons, Blonk et al. (1995) show that the vertical components from vertical point force dominate land seismic data. In reality, usually only vertical components and vertical sources are available. Therefore, we use only vertical particle velocity from a vertical point force source as an approximation for seismic interferometry in elastic media instead of summation over different types of sources. To clearly explain how to perform interferometry for real seismic data and predict groundroll, we built an elastic model to generate synthetic data. We generate the groundroll and reflection events separately and then combine them into a synthetic dataset. Such generation of synthetic data allows more accurate modeling of groundroll at a finer grid and gives us a ground truth for the reflection. Figure 2.7(a) shows the P-wave velocity model for generation of reflection events. S-wave velocity and density are not shown in Figure 2.7(a) but have similar structures. Instead of horizontal layers, this model has more complex structures and produces more complicated reflections. Groundroll is generated from a model with linear increase of P and S wave velocities from the surface to 50m depth over a half space, as shown in Figure 2.7(b). A line of receivers is placed horizontally on the free surface with 8m sampling intervals. A single vertical source is put at the location of a receiver then moved along the receiver line in sequential order to produce shot gathers. This setting simulates realistic exploration surveys. We only consider vertical sources that are restricted to the free surface. We generate the synthetic data by elastic finite difference modeling. The modeling code we use has fourth order accuracy in space and second order accuracy in time. We use a 20Hz central frequency Ricker wavelet as the source wavelet. We set the spatial grid to 2m in both directions (i.e. x and z) and time step to 0.2ms for generating reflections and we set the spatial grid to be 0.5m and time step to 0.05ms for generating groundroll. Figure 2.8(a) shows a combined synthetic shot gather with a source at 1.0km 20 Chapter 2. Groundroll prediction A B (a) (b) Figure 2.7: P-wave velocity model for generating (a) reflections and (b) groundroll. 21 Chapter 2. Groundroll prediction at the surface. Figure 2.8(b) shows the true reflections for the shot gather in Figure 2.8(a). (a) (b) Figure 2.8: (a) a synthetic shot gather at 1.0 km (b) true reflections contained in (a) Letter A and B in Figure 2.7(a) indicate two receiver positions. According to Equation 2.4, we can reconstruct the response at B from A if we have all types of sources (i.e. all the sources including K = 0, 1, 2, 3 in Equation 2.4) located at a boundary that encloses A and B. However, we only have the vertical force sources and they are located at the free surface, so we will only have an approximation to the exact Green’s function. Similar to previous analysis for the acoustic example, we perform the stationary phase analysis for sources at the surface. For a source located at the surface to the left of A, the correlation of the responses at A and B eliminates the common ray path (dash line) and shows residual from A to B (solid line). Since groundroll travels along the surface from A to B, sources to the left of A are stationary sources for groundroll observed at B. For a source close to A, the correlation of responses at A and B is close to correlation of the source wavelet and the response at B. Therefore sources close to A are stationary sources for all the events in the response at B Figure 2.9 shows the 22 Chapter 2. Groundroll prediction correlated traces at A and B from different source at the free surface. All the sources located to the left of A contribute constructively to groundroll as indicated by the purple arrows. Sources close to A contribute to all the events as indicated by the green arrows. If we then remove the sources close to A by muting and sum only part of the sources ( i.e. not summing the sources close to A as in standard interferometry), we get a response dominated by groundroll. We can use this result as a prediction for groundroll. The advantage of this interferometry prediction is that it does not require any prior information about near surface velocity and it is fully data-driven. Figure 2.9: Correlated traces at A and B with auto gain control from difference source locations at the surface. 2.3 Data matrix formulation of interferometry As discussed in the previous section, interferometry with partial summation provides a data-driven prediction for groundroll. To simplify our notation, we rearrange the interferometry operation into the form of matrix multiplication, using Berkhout’s monochromatic matrix notation. For a 2D seismic line with fixed receivers, we arrange the shot gathers into a 3D time-receiver-shot data volume as done within SRME (Verschuur 23 Chapter 2. Groundroll prediction (a) (b) Figure 2.10: (a) Example of 2D line data cube in time domain (b) Example of 2D line data cube in frequency domain. 24 Chapter 2. Groundroll prediction Primary operator [Berkhout & Verschuur ‘96] Shots Receivers Shots ∆P Receivers Frequency Frequency slice from data cube Figure 2.11: Frequency slice of data volume. Source Sg˘ T Receiver X ˆ ∗ )T R(P Source = Receiver Receiver Source ˆ RP Figure 2.12: Prediction of groundroll for a monochromatic matrix. 25 Chapter 2. Groundroll prediction 0 0 1 1 1 1 0 0 0 1 1 1 1 0 0 0 1 1 1 1 0 0 0 1 1 1 1 0 0 0 1 1 1 1 0 0 Figure 2.13: Example of R matrix. and Berkhout, 1997). Figure 2.10(a) shows an example of the 2D seismic line data cube. The three axises are time, receiver position and shot position. The data volume is then Fourier transformed along the first axis (i.e. time) into the frequency domain, as shown in Figure 2.10(b). Then we take a single frequency slice out of the frequency data volume, as shown in Figure 2.11, which forms the monochromatic data matrix P. The rows of P correspond to different shots and the columns of P correspond to different receivers. P∗ denotes taking the complex conjugate for each element of P and PT denotes the transpose of the matrix P. We can also arrange the groundroll prediction of a 2D seismic line into the same data volume structure and take the equivalent monochromatic data matrix out, denoted by S. For each monochromatic matrix of original data and groundroll prediction, we have the following relationship: S = R(P∗ )T RP (2.5) where R is a matrix with 0s in the diagonals and 1s elsewhere e.g. Figure 2.13. R removes the stationary sources for all events, leaving only stationary sources for groundroll. We choose the width of the muting window to be the same scale with the wavelength. Since S and P are in the frequency domain, the matrix multiplication between R(P∗ )T and RP is mathematically equivalent to the correlation and stack operation in interferometry. We apply the interferometry prediction for the synthetic data generated from Figure 2.7(a). Figure 2.14(a) shows the prediction of groundroll. In 26 Chapter 2. Groundroll prediction Figure 2.8(a) groundroll is dominant in the original data, but there are still visible reflections in the image. In the correlation prediction shown in Figure 2.14(a), there are barely any reflections. Comparing the prediction with true groundroll, we find that the prediction has the same dip with true groundroll. We also note that the amplitude of the prediction is not at the same scale of the original data. The phase and position are also not exactly the same. This is mainly due to the imprint of the source wavelet and vertical component approximation we use. Therefore, we have to adaptively match the prediction instead of directly subtracting it from the original data. We will discuss the separation step in the next two chapters. (a) (b) Figure 2.14: (a) Interferometry prediction of the groundroll for the shot gather at 1.0 km. (b) True groundroll for the shot gather at 1.0km. 27 Chapter 3 Groundroll separation by curvelet-domain matched filtering In the previous chapter, we put emphasis on how to predict groundroll. The next step is to separate the prediction from the total data that contains true groundroll and reflection. However, as discussed in the previous chapter, the prediction is not accurate and therefore we cannot directly subtract the prediction from the total data. There are several reasons why the prediction is inaccurate. First, the interferometry prediction is based on the correlation of two waveforms (instead of Green’s function) and consequentially correlates the imprint of the source signature. Second, several approximations and assumptions are used in the cross-correlation interferometry. For example, sources are assumed to be uniformly distributed at a boundary that encloses the receivers. The medium is assumed to be nonattenuating in the domain and homogeneous at and outside the boundary of sources. For elastic seismic interferometry, we need displacement responses at all directions from different type of sources (K = 0, 1, 2, 3 in Equation 2.5). Since we only use the vertical displacements from vertical sources and lack other components, the cross-correlation interferometry we have is only an approximation to the full elastic case. We only have sources located at the surface, which is also an approximation. One of the results of using only sources at the surface is that these surfaces sources are stationary sources for groundroll, therefore the reconstructed shot gathers are dominated by groundroll i.e. it provides a prediction for groundroll, although it is not accurate. Due to the approximations and inaccuracy, we must match the prediction to the true groundroll in the original data before subtraction, i.e. the prediction needs to be adaptively subtracted from the original data. In this chapter, we will discuss the adaptive subtraction methods and focus on the curvelet domain matched filtering (Herrmann et al., 2008), which has 28 Chapter 3. Groundroll separation by curvelet-domain matched filtering advantages over conventional Fourier domain matched filtering. 3.1 Least-square adaptive subtraction strategies As mentioned earlier, the interferometry prediction of groundroll is usually imperfect and inaccurate. Therefore, we rely on adaptive subtraction to fit the prediction to the true groundroll contained in the original data. The common adaptive subtraction method is Fourier domain matched filtering. The underlying assumptions for the Fourier domain matched filtering is that the errors in the prediction are stationary and signal after subtraction is orthogonal to the noise and therefore has minimum energy. The Fourier domain filtering operation can be viewed as a scaling of the Fourier coefficients of the prediction: mtrue ≈ Fmpredicted with F = F H diag ˆf F (3.1) where mpredicted is the inaccurate prediction that needs to be matched. mtrue is the true groundroll, which is usually unknown. F is the Fourier domain filtering operation that matches the prediction to the true groundroll. F is the Fourier transform which decomposes the prediction into the frequency domain. Vector ˆf is the Fourier domain filter that corrects the coefficients of the prediction. F H is the inverse Fourier transform that reconstructs the filtered prediction in the physical domain. The filter ˆf is estimated by solving a least-square minimization problem: ˆf = arg min 1 d − F H diag (ˆ g) Fm ˆ 2 g 2 2 ˆ 22 . + λ LF g (3.2) Since we do not know mtrue , total original data d that contains mtrue is used as a substitute for the reference of matching. To avoid over-fitting the total data (d also contains reflection events), smoothness in the Fourier domain i.e. fast decay in the physical domain is imposed as a regularization. The ˆ between adjacent frequenoperator LF computes the difference of the filter g 2 ˆ 2 as a penalty term in the minimization objective cies. By adding λ LF g function, the minimization process in Equation 3.2 seeks for a solution that not only minimizes the misfit of matching, but also minimizes the difference of the filter between adjacent frequencies, i.e. promotes smoothness of the filter in the Fourier domain. The regularization parameter λ controls ˆ in the frequency domain. For a the amount of smoothness of the filter g larger λ, there is more emphasis on the penalty of the difference between 29 Chapter 3. Groundroll separation by curvelet-domain matched filtering ˆ , thus the solution is more smooth in the adjacent frequencies of the filter g Fourier domain. Because smoothness in the frequency domain is equivalent ˆ 22 can be to fast decay in the time domain, the regularization term λ LF g incorporated implicitly by adding a short-time window to constrain the filter ˆ. g Verschuur and Berkhout (1997) argue that the least square subtraction is best applied in two stages: firstly, global subtraction per shot record to remove the overall source signature; secondly, local subtraction for each shot in overlapping space-time windows. The reason behind this approach is that in SRME method, the errors in the prediction for multiples come from the convolution of the source signature and non-ideal local mismatch that depends on offset, time and dip. Similarly, the errors in the interferometry prediction for groundroll come from the correlation of the source signature and local imperfections which depend on offset, time and dip. Therefore, the adaptive subtraction for groundroll should also be carried in two stages: global and local Fourier domain matched filtering. The global Fourier filters handle the spectrum mismatch and global kinematic errors which mainly occur due to the correlation of the source signature. In actuality, because of 3D geology of the Earth, directivity effects in sources and receivers, and dispersive effects of groundroll, it is better to estimate the global filters per offset instead of per shot gather. The windowed Fourier domain matched filtering mitigates the local imperfections (as a function of offset, time and dip) in the groundroll prediction. However, the windowed Fourier matched filtering runs the risk of over fitting the data (i.e. it removes part of the reflection together with groundroll) because it has limited control over window-to-window variations. It is also known that local imperfections that depend on the location and dip in the prediction is particularly challenging for the Fourier domain matched filtering. When applied to multiple primary separation problem, the least square Fourier matching technique leads to residual multiple energy, highfrequency clutter and deterioration to the primaries (Herrmann et al., 2007). With the similarities of multiple primary separation and groundroll reflection separation, we expect analogous problems will arise when applying the least square Fourier matching technique to groundroll matching. Due to the large amplitudes and dispersion of groundroll, these problems will be even more severe for groundroll separation. A small mismatch in the groundroll prediction may cause coverage of important reflection information. To avoid the problems in windowed Fourier domain matched filtering, Herrmann et al. (2008) has proposed an alternative — the curvelet domain matched filtering and applied it successfully to multiple separation. Since the groundroll separation problem has similarities to the multiple separation 30 Chapter 3. Groundroll separation by curvelet-domain matched filtering problem (i.e. offset, time and dip dependent errors in the prediction), we apply curvelet-domain matched filtering to match the groundroll prediction. Our matching approach consists of two stages: first global Fourier domain matched filtering and second curvelet-domain matched filtering. The imprint of the source signature i.e. kinematic shifts and amplitude spectral mistakes are removed by a conservative global Fourier domain matched filtering procedure per offset. Remaining dip and offset dependent errors are addressed by the curvelet domain matched filtering. 3.2 Curvelet-domain matched filtering After we remove the correlation of source signature and directivity by global Fourier domain matching, the prediction of groundroll can be viewed as a non-stationary scaling of the true groundroll. Following Herrmann et al. (2008), this non-stationary scaling can be mathematically modeled as a zero order pseudo-differential operator. This operator links the prediction to the true groundroll in the original data: mtrue = Bm0 with B ≈ CT diag (b) C, {b}µ∈M > 0 (3.3) where B is a full positive definite matrix, approximating the zero order pseudo-differential operator. Through the global Fourier matched filtering, we gain m0 = F H diag ˆf Fmpredicted — predicted groundroll without source signature and directivity. The true groundroll is denoted by mtrue . The action of B is approximated by a curvelet domain filtering operation , which has a similar form compared to the Fourier domain matched filtering. The matrix C represents the curvelet transform, vector b is the curvelet domain scaling vector and M is the index set of curvelet coefficients. The curvelet transform C decomposes m0 into the curvelet domain. After filtered by b, the inverse curvelet transform CT reconstructs the curvelet matched groundroll in the physical domain. In this approximate model, the groundroll prediction is connected to the true groundroll through a scaling of the curvelet coefficients. The scaling of the curvelet coefficients fixes the amplitude errors in the groundroll prediction which vary smoothly as a function of location, scale and dip. Following Herrmann et al. (2008), we estimate the scaling vector b in Equation 3.3 through the least-square optimization that minimizes the l2 misfit between the true groundroll mtrue and curvelet-domain matched groundroll 31 Chapter 3. Groundroll separation by curvelet-domain matched filtering CT diag Cm0 b. Since we do not know the true groundroll mtrue , as a reference we use the total data d which includes the true groundroll. According to Herrmann (2008), there are several issues that make the estimation of the scaling vector difficult: (i) the forward model is undetermined due to the redundancy of the curvelet transform (i.e. CCT is not identity); (ii) it has the risk of over-fitting the total data which contains reflections; (iii) it is requisite the scaling vector is positive. As with the Fourier domain matching, which is regularized by imposing the Fourier domain smoothness, the curvelet domain matching imposes smoothness in the curvelet domain to address issue (i) and (ii). To impose curvelet domain smoothness (i.e. to control curveletto-curvelet variations) for the scaling vector b in Equation 3.3, we add the penalty term γ LC b 22 in Equation 3.4. The curvelet-domain sharpening operator LC consists three terms, namely DT1 ,DT2 , and DTθ , which compute the differences for the curvelet-domain scaling vector in space (x1 , x2 ) and angle (θ) directions for each scale: b = arg min b 1 d − CT diag Cm0 b 2 2 2 + γ LC b 22 , LC = [DT1 DT2 DTθ ]. (3.4) (3.5) The parameter γ in Equation 3.4 controls the amount of smoothness of the scaling vector. For a larger γ, there is more emphasis on penalizing the variations between different curvelets, and the scaling vector becomes more smooth in the curvelet domain. To address issue (iii), we substitute the scaling vector b with exponentiation ez , where the exponentiation is taken element-wise for vector z. Finding the solution for Equation 3.4 is equivalent to solving the following augmented system: d CT diag{Cm0 } = b with b = ez > 0 0 γL (3.6) by the minimizing the functional: 1 ˜ d − Hγ ez 22 2 CT diag{Cm0 } H= . γL Jγ (z) = ˜= d where d 0 (3.7) This nonlinear system can be solved using the limited-memory BFGS solver (Nocedal and Wright, 1999) with the gradient: 32 Chapter 3. Groundroll separation by curvelet-domain matched filtering ˜ gradJ (z) = diag{ez }[HTγ (Hγ )ez − d)]. 3.3 (3.8) Synthetic examples and discussions We test our matching approach on the synthetic data which has been produced in the previous chapter. (See Figure 2.7 in Chapter 2). Figure 3.1 shows the result of global Fourier matched groundroll and the reflection produced by subtracting global Fourier matched groundroll. Figure 3.1(a) demonstrates that the global Fourier matching fixes the global spectral error by matching the groundroll prediction to the right scale and phase. Figure 3.1(b) shows that subtracting this global Fourier matched groundroll from the original data reveals some of the reflection events but there is still a significant portion of residual groundroll remaining. The SNR of the reflection produced by subtracting global Fourier matched groundroll is 8.66 dB. Figure 3.2 shows the result of curvelet matched groundroll and reflection produced by subtracting the curvelet matched groundroll from the original data. We test a few parameters for γ from 0.001 to 1.0 and find γ = 0.02 to be an optimal choice. Table 3.1 shows the SNR of the reflections produced by subtracting curvelet matched groundroll for different parameter selections. From the table we can see results of γ = 0.02 produced the highest SNR, which has 3.33 dB improvement compared with global Fourier matched result. Figure 3.2(a) is the curvelet-domain matched groundroll with input being the global Fourier matched groundroll. Figure 3.2(b) illustrates that subtracting the curvelet-domain matched groundroll from the original data removes much of the groundroll residual found in Figure 3.1(b). Figure 3.2(c) shows the difference between true reflection and reflection produced by subtracting the curvelet matched groundroll. From Figure 3.2(c) we can see there are still some groundroll residuals left in reflection produced by subtracting the curvelet matched groundroll. There are also some leakages to reflections at places where groundroll and reflections overlap with each other, which is moderate. Figure 3.3 shows the comparison of our separation result and true reflection at offset 0.4km. Figure 3.3(a) is the comparison of true reflection and reflection from subtracting global Fourier matched groundroll. Figure 3.3(b) is the comparison of true reflection and reflection from subtracting the curvelet-domain matched groundroll. Figure 3.3 demonstrates that our 33 Chapter 3. Groundroll separation by curvelet-domain matched filtering two stage matching procedure removes most of groundroll while conserves reflections. Figure 3.4 shows the curvelet domain matched filter for γ = 0.02. In Figure 3.4(a), the low frequency components (i.e. coarse scales) of the curvelet domain filter are in the center and high frequency components (i.e. fine scales) are at the perimeter. The angle increases clockwise and the number of angles double every other scale. The mosaic plot shows that the curvelet domain matched filter is smooth along two space directions and the angle direction. (a) (b) Figure 3.1: (a) Global Fourier matched groundroll for shot gather at 1.0 km (b) reflection produced by subtracting (a) from the total data SNR (dB) γ = 1.0 8.99 γ = 0.1 10.92 γ = 0.02 11.99 γ = 0.01 11.23 γ = 0.001 6.25 Table 3.1: SNR of the reflections produced by subtracting curvelet matched groundroll for different parameter selections. Ideally, the result after curvelet-domain matching should be the appro34 Chapter 3. Groundroll separation by curvelet-domain matched filtering (a) (b) (c) Figure 3.2: (a) Curvelet matched groundroll with γ = 0.02 for shot gather at 1.0 km. (b) Reflection produced by subtracting (a) from original data. (c) difference between (b) and true reflection 35 Chapter 3. Groundroll separation by curvelet-domain matched filtering (a) (b) Figure 3.3: (a) Comparison of reflection produced by global Fourier matching with the true reflection at offset 0.4km (b) Comparison of reflection produced by curvelet-domain matching with the true reflection at offset 0.4km. Shot location is at 1.0km. Red lines represent the true reflection and blue lines represent reflection results produced by subtracting global Fourier matched groundroll and curvelet matched groundroll. 36 Chapter 3. Groundroll separation by curvelet-domain matched filtering 18 16 14 12 10 8 6 4 2 0 0 2 4 6 8 10 12 14 5 x 10 (a) (b) Figure 3.4: (a) Mosaic and (b) 2D plots of the curvelet domain filter for shot gather at 1.0 km and γ = 0.02. In (a), the low frequency components (i.e. coarse scales) of the curvelet domain filter are in the center and high frequency components (i.e. fine scales) are at the perimeter. The angle increases clockwise and the number of angles double every other scale. The color represents the amplitudes of the curvelet coefficients. priate estimate for groundroll. However, other kinematic and phase errors that could not be removed by global Fourier domain matching may affect our curvelet-domain filtering, making the direct subtraction less effective. As shown in Figure 3.2(b) and Figure 3.2(c), there ares still residual errors left in the reflection produced by subtracting curvelet domain matched groundroll from the original data. To handle these errors, in the next chapter, we introduce a Bayesian method that separates the curvelet matched groundroll and reflection robustly by promoting sparsity in curvelet domain. 37 Chapter 4 Bayesian separation for groundroll In the previous chapter, we discussed the curvelet domain matched filtering that fits the prediction to the true groundroll in the original data. Ideally, the scaled groundroll, yielded by the solution of the curvelet-domain matched filtering, can be subtracted directly from the total data. Unfortunately, noises in seismic data and imperfections in the prediction may affect our curvelet domain filtering, making direct subtraction less effective. As shown in the previous chapter, there are still residual errors left in the reflection produced by subtracting the curvelet domain matched groundroll from the original data. Herrmann et al. (2008) has shown that a similar problem arises in multiple separation. Due to the large amplitudes of groundroll, the imperfection in the matched groundroll could be a more severe problem, i.e. a relatively small residual error in the matched groundroll can cover important reflection information. Recently a Bayesian separation algorithm has been proposed to address this problem and has been successfully applied to multiple separation (Saab et al., 2007; Wang et al., 2008). Yarham and Herrmann (2008) applied this Bayesian method to improve the groundroll separation when a prediction for groundroll from conventional methods is provided. In this chapter, we apply the Bayesian algorithm to robustly separate the curvelet matched groundroll instead of direct subtraction. 4.1 Bayesian wavefield separation algorithm The Bayesian algorithm (Wang et al., 2008) is proposed to separate two coherent wavefield components, when a prediction is given for one of the components. This thesis addresses how to separate groundroll and reflection. The prediction of groundroll which contains moderate residual errors (not the interferometry prediction in Chapter 2 ) is given by the two stage 38 Chapter 4. Bayesian separation for groundroll matching process in Chapter 3 (first conservative global Fourier matching and then curvelet matching) and we refer to it as curvelet matched groundroll. We form our forward model of signal separation problem as: b = s1 + s2 + n (4.1) where b is the total data. Vectors s1 and s2 are reflection and groundroll respectively. We assume the noise n in the total data is Gaussian noise with a standard deviation σ. We denote the curvelet matched groundroll by b2 . Since there are moderate noise residuals left in the curvelet matched groundroll, we assume: b2 = s2 + n2 (4.2) b1 = s1 + n1 (4.3) and we obtain where b1 = b − b2 and n1 = n − n2 . Since we do not have accurate knowledge about the noises, we simply assume the noise terms ni to be independent white Gaussian noises with possibly different standard deviation σi . Such assumption of Gaussian noises is consistent with previous work about matched filtering (Herrmann et al., 2008; Verschuur et al., 1992). We decompose the two unknown signal components into the curvelet domain i.e. xi = Csi , i = 1, 2, where C is the curvelet transform. And we get the following equations: b1 = CT x1 + n1 b2 = CT x2 + n2 (4.4) where CT represents the inverse curvelet transform. The unknown curvelet coefficients for reflection (i.e. x1 ) and groundroll (i.e. x2 ) are linked to the direct subtraction b1 and the curvelet matched groundroll b2 . Using the system equations above, we are able to utilize the sparsity in the curvelet domain to further improve groundroll separation. Since curvelet matched groundroll and reflection prediction (i.e. b1 and b2 ) are known, we aim to find the curvelet coefficients x1 and x2 that maximize the conditional probability P (x1 , x2 |b1 , b2 ). Using Bayes’ rule, we choose to maximize: 39 Chapter 4. Bayesian separation for groundroll P (x1 ,x2 )P (b1 |x1 ,x2 )P (b2 |b1 ,x1 ,x2 ) P (b1 ,b2 ) P (x1 ,x2 )P (n)P (n2 ) . P (b1 ,b2 ) P (x1 , x2 |b1 , b2 ) = = (4.5) In order to get x1 and x2 , we make additional assumption that x1 and x2 have weighted Laplacian prior distributions. Such an assumption acts as a sparsity-promoting prior and is consistent with the high compression rate of seismic data in the curvelet domain (Wang et al., 2008). Following previous work by Saab et al. (2007); Wang et al. (2008); Yarham and Herrmann ˆ 1 and x ˆ 2 by solving the optimization problem: (2008), we find estimations x arg max P (x1 , x2 |b1 , b2 ) = arg max P (x1 , x2 )P (n)P (n2 ) x1 ,x2 x1 ,x2 = arg max exp x1 ,x2 − λ1 ||x1 ||1,w1 − λ2 ||x2 ||1,w2 − ||CT x2 − b2 ||22 σ22 − ||CT (x1 + x2 ) − (b1 + b2 )||22 σ2 = arg max − λ1 ||x1 ||1,w1 + λ2 ||x2 ||1,w2 x1 ,x2 + ||CT x2 − b2 ||22 σ22 + ||CT (x1 + x2 ) − (b1 + b2 )||22 . σ2 (4.6) where ||xi ||1,wi = µ∈M |wi,µ xi,µ | is the weighted l1 norm of the curvelet coefficients xi and M is the index set for the curvelet coefficients. According to Saab et al. (2007); Wang et al. (2008), the maximization of the conditional probability can be cast into: arg max P (x1 , x2 |b1 , b2 ) = arg min f (x1 , x2 ) x1 ,x2 x1 ,x2 = arg min λ1 ||x1 ||1,w1 + λ2 ||x2 ||1,w2 x1 ,x2 +||CT x2 − b2 ||22 +η||CT (x1 + x2 ) − (b1 + b2 )||22 (4.7) 40 Chapter 4. Bayesian separation for groundroll Following previous work by Herrmann et al. (2007); Wang et al. (2008) and their empirical findings, we choose the weights using curvelet matched groundroll and reflection produced by direct subtraction: w1 = |Cb2 |, w2 = |Cb1 | (4.8) This choice of weights ensures they are strictly positive. Furthermore, it ˆ1) makes it less likely that the curvelet coefficient vector for reflection (i.e. x ˆ 2 ). Therefore, will overlap with the curvelet coefficient for groundroll (i.e. x the two signal components: reflection and groundroll will be driven apart by the weights. Equation 4.7 minimizes the sum of four components: the weighed l1 norm of curvelet coefficients for groundroll and reflection, the l2 misfit of the curvelet matched groundroll, and l2 misfit of the total data. The two parameters λ1 , λ2 in Equation 4.7 controls the relative sparsity of groundroll and reflection. If we expect sparser reflection compared to groundroll, we set λ1 > λ2 , putting more weight on the penalty of the l1 norm of reflections. The parameter η controls the confidence level for the curvelet matched groundroll. Increasing η will put more emphasis on the penalty of misfit of total data and less emphasis on the misfit of curvelet matched groundroll, i.e. a larger η corresponds to a worse prediction of groundroll. Equation 4.7 leads to solutions that are not only sparse, but also fit the curvelet matched groundroll and total data. The minimization of f (x1 , x2 ) can be solved by soft iterative thresholding. From the initial starting conditions x01 , x02 , the separation algorithm computes the nth iteration as: xn+1 = T λ1 w1 [Cb2 − CCT xn2 + Cb1 − CCT xn1 + xn1 ] 1 2η xn+1 2 = T λ2 w2 [Cb2 − CCT xn2 + xn2 + 2(1+η) η (Cb1 − CCT xn1 )] (4.9) 1+η where Tu is the element-wise soft thresholding operator on a vector v: Tuµ := sgn(vµ ) max(0, |vµ | −| uµ |), µ ∈ M. (4.10) Here uµ is the threshold value and M is the index set for v. With the use of strictly positive weights w1 and w2 (Equation 4.8), the algorithm above converges to the global minimum (Daubechies et al., 2004; Wang et al., 2008). Table 4.1 demonstrates the algorithm in detail. 41 Chapter 4. Bayesian separation for groundroll Initialize: x01 = 0,x02 = 0, m = 0, Choose: mmax , η, λ1 , λ2 Define Threshold: w1 = (λ1 |Cb2 |)/(2η), w2 = (λ2 |Cb1 |)/(2(η + 1)) while m < mmax do m=m+1 T m m ˜ 1 = Cb2 −CCT xm x 2 +Cb1 −CC x1 +x1 ; {Coefficient update} η T m m ˜ 2 = Cb2 − CC x2 + x2 + η+1 (Cb1 − CCT xm x 1 ); {Coefficient update} xm+1 = sign(˜ x1 ) · max(0, |˜ x1 | −| w1 |); {Soft threshold} 1 xm+1 = sign(˜ x ) · max(0, |˜ x2 | −| w2 |); {Soft threshold} 2 2 end while Table 4.1: The Bayesian-separation algorithm. Table from Yarham and Herrmann (2008) 4.2 Application to synthetic data We apply the Bayesian separation method to the synthetic data we used in the previous chapters. We test several pairs of the parameters and find λ∗1 = 2.0, λ∗2 = 8.0, and η ∗ = 2.0 to be a optimal choice. The optimal choice of parameters leads to a SNR of 13.23 dB. The bigger value for λ2 is consistent with sparser groundroll than reflections. Table 4.2 shows SNR of Bayesian separated reflection for different parameter choice. From the SNR numbers, we can see the algorithm is quite stable for different parameter choice. Fig 4.1 shows the groundroll and reflection produced by the Bayesian separation algorithm for optimal parameter choice λ∗1 = 2.0, λ∗2 = 8.0, and η ∗ = 2.0. Compared with reflection produced by subtracting curvelet matched groundroll ( i.e. Fig 3.2(b) ), our Bayesian algorithm removed partial amount of groundroll residual that could not be resolved by the curvelet domain matched filtering. The SNR also shows a 1.24 dB improvement compared with the curvelet matched reflection for γ = 0.02 which has a SNR of 11.99 dB. 42 Chapter 4. Bayesian separation for groundroll (a) (b) Figure 4.1: (a) Groundroll produced by Bayesian separation (b) Reflection produced by Bayesian separation for λ1 = 2.0, λ2 = 8.0, and η = 2.0 43 Chapter 4. Bayesian separation for groundroll SNR (dB) η∗ 1 ∗ 2 ·η ∗ 2·η λ∗1 , λ∗2 13.23 13.14 13.15 2 · λ∗1 , λ∗2 13.10 12.62 13.06 λ∗1 , 2 · λ∗2 13.20 13.14 13.21 Table 4.2: SNR of the Bayesian separated reflections for different parameter selections. 44 Chapter 5 Conclusion and future research 5.1 A data-driven workflow for groundroll removal Groundroll is one of most common coherent noises in land seismic data. Groundroll and reflection seismic events usually have significant overlap in both time and frequency domain. How to remove groundroll while preserving reflection remains an active research topic in exploration geophysics. Our goal for this thesis is to develop a data-driven workflow to remove groundroll noise in land seismic data. In this thesis, we have proposed a workflow that is motivated both by by surface related multiple elimination (SRME- Verschuur et al., 1992) and recently developed interferometry method for groundroll prediction (Curtis et al., 2006; Halliday et al., 2007; Dong et al., 2006; Vasconcelos et al., 2008). Our workflow to remove groundroll consists of two steps: prediction and adaptive subtraction. For prediction, we exclude stationary sources for reflections and use only stationary sources for groundroll in the inteferometry. We show that such a cross-correlation based interferometery can serve as a data-driven but imperfect prediction for groundroll, which need to be adaptively matched to the true groundroll in the original data. We show that imperfections in the inteferometry prediction for groundroll mainly come from two parts: (1) the stationary differences which mainly occur due to the correlation of source signature during interferometry, and (2) non-stationary difference which vary smoothly as a function of dip and offset along the wavefronts. For adaptive subtraction, we remove stationary differences, i.e., global kinematical errors and amplitude spectra mismatch, through a global Fourier matching process. Then we address the remain45 Chapter 5. Conclusion and future research ing dip and offset dependent errors by curvelet domain matched filtering. We show that position, scale and dip dependent amplitude errors in the prediction can be successfully corrected through a diagonal scaling of curvelet coefficients. Overfitting to reflections is avoided by imposing smooth constraints among neighboring curvelet coefficients of the diagonal scaler. Finally, we apply a Bayesian separation algorithm to further improve the separation result over the output of typical direct subtraction. We show that the sparsity of seismic data provided by the curvelet transform increases the separation of coherent wavefield and allows us to utilize robust signal separation through l1 mininization. 5.2 Open and future research Recently proposed interfeometry method to predict groundroll adds a new family of groundroll removal methods based on the predictability of groundroll, different from the traditional method based on the low velocity and low frequency properties of groundroll. In this thesis, we demonstrate that for a 2D seismic line, interferometry using inline sources can effectively predict inline groundroll that has linear moveout. For 3D seismic survey, groundroll has hyperbolic moveout when there is non-zero distance from the source to the receiver line. Crossline groundroll with hyperbolic moveout is particularly challenging for traditional f-k filtering method due to its significant overlap with reflections. One possible future research direction is to expand the interferometry prediction to 3D seismic survey, using out of line sources to predict crossline groundroll and then adaptively subtract it. 46 Bibliography Aki, K., and P. G. Richards, 2002, Quantitative seismology. University Science Books. Blonk, B., G. C. Herman, and G. G. Drijkoningen, 1995, An elastodynamic inverse scattering method for removing scattered surface waves from field data: Geophysics, 60, 1897–1905. Cand´es, E., D. Donoho, et al., 2000, Curvelets: A surprisingly effective nonadaptive representation for objects with edges. Cand`es, E., J. Romberg, and T. Tao, 2006, Stable signal recovery from incomplete and inaccurate measurements: Communications on Pure and Applied Mathematics, 59, 1207. Cand`es, E. J., and D. L. Donoho, 2004, New tight frames of curvelets and optimal representations of objects with C 2 singularities: Communications on Pure and Applied Mathematics, 57, no. 2, 219 – 266. Claerbout, J. F., 1968, Synthesis of a layered medium from its acoustic transmission response: Geophysics, 33, 264–269. Curtis, A., P. Gerstoft, H. Sato, R. Snieder, and K. Wapenaar, 2006, Seismic interferometry—turning noise into signal: The Leading Edge, 25, 1082–1092. Daubechies, I., M. Defrise, and C. De Mol, 2004, An iterative thresholding algorithm for linear inverse problems with a sparsity constraints: 57, 1413– 1457. Deighan, A., and D. Watts, 1997, Ground-roll suppression using the wavelet transform: Geophysics, 62, 1896–1903. Dong, S., R. He, and G. T. Schuster, 2006, Interferometric predcition and least squares subtraction of surface waves: SEG Technical Program Expanded Abstracts, 2783–2786, SEG. Embree, P., J. P. Burg, and M. M. Backus, 1963, Wide band velocity filtering-the pie-slice process: Geophysics, 28, 948–974. 47 Bibliography Fail, J. P., and G. Grau, 1963, Les filters on eventail: Geophys. Prospect., 11, 131–163. Galbraith, J. N., and R. A. Wiggins, 1968, Characteristics of optimum multichannel stacking filters: Geophysics, 33, 36–48. Halliday, D. F., A. Curtis, J. O. A. Robertsson, and D.-J. van Manen, 2007, Interferometric surface-wave isolation and removal: Geophysics, 72, A69–A73. Hennenfent, G., L. Fenelon, and F. J. Herrmann, 2010, Nonequispaced curvelet transform for seismic data reconstruction: a sparsity-promoting approach: Technical Report TR-2010-2, UBC-Earth and Ocean Sciences Department. Hennenfent, G., and F. J. Herrmann, 2008, Simply denoise: wavefield reconstruction via jittered undersampling: Geophysics, 73, no. 3. Herrmann, F. J., 2008, Curvelet-domain matched filtering: SEG Technical Program Expanded Abstracts, 3643–3647, SEG. Herrmann, F. J., U. Boeniger, and D. J. Verschuur, 2007, Non-linear primary-multiple separation with directional curvelet frames: Geophysical Journal International, 170, no. 2, 781–799. Herrmann, F. J., and G. Hennenfent, 2008, Non-parametric seismic data recovery with curvelet frames: Geophysical Journal International, 173, 233–248. Herrmann, F. J., D. Wang, and D. J. Verschuur, 2008, Adaptive curveletdomain primary-multiple separation: Geophysics, 73, no. 3, A17–A21. Karslı, H., and Y. Bayrak, 2004, Using the Wiener–Levinson algorithm to suppress ground-roll: Journal of Applied Geophysics, 55, 187–197. Liu, X., 1999, Ground roll spression using the Karhuren-Loeve transform: Geophysics, 64, 564–566. McMechan, G. A., and R. Sun, 1991, Depth filtering of first breaks and ground roll: Geophysics, 56, 390–396. Nocedal, J., and S. J. Wright, 1999, Numerical optimization: Springer. Russel, B., D. Hampson, and J. Chun, 1990, Noise elimination and the radon transform; part 1: The Leading Edge, 9, no. 10, 18–23. Saab, R., D. Wang, O. Yilmaz, and F. Herrmann, 2007, Curvelet-based primary-multiple separation from a bayesian perspective: Presented at the SEG International Exposition and 77th Annual Meeting. 48 Bibliography Schuster, G. T., J. Yu, J. Sheng, and J. Rickett, 2004, Interferometric/daylight seismic imaging: Geophysical Journal International, 157, 838–852. Telford, W. M., L. P. Geldart, and R. E. Sheriff, 1990, Applied geophysics. Cambridge University Press. Trad, D., T. Ulrych, and M. Sacchi, 2003, Latest views of the sparse radon transform: Geophysics, 68, 386–399. Treitel, S., J. L. Shanks, and C. W. Fraster, 1967, Some aspects of fan filtering: Geophysics, 32, 789–800. Vasconcelos, I., J. Gaiser, A. Calvert, and C. Calder´on-Mac´ıas, 2008, Retrieval and suppression of surface waves using interferometry by correlation and by deconvolution: SEG Technical Program Expanded Abstracts, 2566– 2570, SEG. Verschuur, D. J., and A. J. Berkhout, 1997, Estimation of multiple scattering by iterative inversion, part II: practical aspects and examples: Geophysics, 62, 1596–1611. Verschuur, D. J., A. J. Berkhout, and C. P. A. Wapenaar, 1992, Adaptive surface-related multiple elimination: Geophysics, 57, 1166–1177. Wang, D., R. Saab, O. Yilmaz, and F. J. Herrmann, 2008, Bayesian wavefield separation by transform-domain sparsity promotion: Geophysics, 73, no. 5. Wapenaar, K., and J. Fokkema, 2006, Green’s function representations for seismic interferometry: Geophysics, 71, SI33–SI46. Yarham, C., and F. J. Herrmann, 2008, Bayesian ground-roll seperation by curvelet-domain sparsity promotion: SEG Technical Program Expanded Abstracts, 2576–2580, SEG. Yilmaz, O., 2001, Seismic data analysis. Curves and Surfaces. Society of Exploration Geophysicists. Zhang, R., and T. J. Ulrych, 2003, Physical wavelet frame denoising: Geophysics, 68, 225–231. 49 Appendix A A real example for curvelet matching and Bayesian separation In reality, depends on the survey geometry and data quality, the interferometry prediction may not be available or at a satisfactory level. In this appendix, we show that our curvelet domain matching and Bayesian separation can still serve as a powerful tool for separating groundroll when other prediction for groundroll is available. Figure A.1(a) shows a real shot gather from a 2D seismic line. The receivers are placed with a distance of 25 meters and the sampling interval in time is 4ms. Figure A.1(b) shows the groundroll prediction produced by Fourier domain filtering and threshold. Figure A.2(a) shows the reflections produced by directly subtract the Fourier domain groundroll prediction. Figure A.2(b) shows the reflection produced by curvelet domain matching and Bayesian separation. Notice the improvement around near offset. 50 Appendix A. A real example for curvelet matching and Bayesian separation (a) (b) Figure A.1: (a) A real shot gather. (b) groundroll prediction produced by Fourier domain threshold and filtering. 51 Appendix A. A real example for curvelet matching and Bayesian separation (a) (b) Figure A.2: (a) Reflection produced by direct subtraction of the Fourier domain groundroll prediction. (b) Separated reflection through curvelet domain matched filtering and Bayesian separation. 52
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Seismic groundroll prediction by interferometry and...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Seismic groundroll prediction by interferometry and separation in curvelet domain Yan, Jiupeng 2011
pdf
Page Metadata
Item Metadata
Title | Seismic groundroll prediction by interferometry and separation in curvelet domain |
Creator |
Yan, Jiupeng |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | Groundroll is a type of surface wave that propagates along the Earth’s surface. Groundroll usually has low frequency, low velocity and high amplitude. Due to its high amplitude, groundroll almost always dominates reflected body waves in land seismic data and covers important reflection information. Therefore, removing groundroll noise is a very important step before seismic imaging. The most common methods used in industry to remove groundroll are the Fourier domain filtering methods based on the different characteristics of groundroll and reflections, i.e. the low frequency and low velocity properties of groundroll. However, groundroll and reflection usually have large overlap in both physical and frequency domain. Also groundroll is spatially aliased at normal receiver intervals causing additional processing difficulties. Therefore, a good separation of groundroll by Fourier domain filtering method is challenging. In this thesis, we propose a data-driven workflow to remove groundroll. Our workflow is motivated both by SRME (Surface Related Multiple Elimination) method and a recently proposed interferometry method for the prediction of groundroll. It consists of a prediction step based on interferometry and a robust separation step that involves curvelet domain matched filtering and sparsity promotion. Tests of our workflow on synthetic data show clear removal of large amplitude groundroll and preservation of seismic reflection events. Test of our separation step on real data shows improvement over conventional Fourier domain filtering methods. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-05-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0053425 |
URI | http://hdl.handle.net/2429/34685 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2011-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2011_fall_yan_jiupeng.pdf [ 2.41MB ]
- Metadata
- JSON: 24-1.0053425.json
- JSON-LD: 24-1.0053425-ld.json
- RDF/XML (Pretty): 24-1.0053425-rdf.xml
- RDF/JSON: 24-1.0053425-rdf.json
- Turtle: 24-1.0053425-turtle.txt
- N-Triples: 24-1.0053425-rdf-ntriples.txt
- Original Record: 24-1.0053425-source.json
- Full Text
- 24-1.0053425-fulltext.txt
- Citation
- 24-1.0053425.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0053425/manifest