Science, Faculty of
Earth, Ocean and Atmospheric Sciences, Department of
DSpace
UBCV
Johnson, James
Master of Science - MSc
University of British Columbia
The primary focus of most reflection seismic surveys is to help locate hydro-carbon recourses. Due to an ever increasing scarcity of these re- courses, we must increase the size and quality of our seismic surveys. How- ever, processing such large seismic data volumes to accurately recover earth properties is a painstaking and computationally intensive process. Due to the way reflection seismic surveys are conducted there are often holes in the collected data, where traces are not recorded. This can be due to physical or cost constraints. For some of the initial stages of process- ing these missing traces are of little consequence. However processes like multiple prediction and removal, interferometric ground roll prediction, and migration require densely sampled data on a regular grid. Thus the need to interpolate undersampled data cannot be ignored. Using the fact that reflection seismic data sets obey a reciprocal relation- ship in source and receiver locations, combined with recent advances in the field of compressed sensing, we show that properly regularized the wavefield reconstruction problem can be solved with a high degree of accuracy. We exploit the compressible nature of seismic data in the curvelet domain to solve regularized l1 recovery problems that seek to match the measured data and enforce the above mentioned reciprocity. Using our method we were able to achieve results with a 20.45 dB sig- nal to noise ratio when reconstructing a marine data set that had 50% of its traces decimated. This is a 13.44 dB improvement over using the same method run without taking reciprocity into account.
Seismic wavefield reconstruction using reciprocity by James Johnson B.Sc.Engineering Physics, Queen’s University Kingston Ontario Canada. B.SC. (Gen) Geology, Queen’s University Kingston Ontario Canada. A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in The Faculty of Graduate Studies (Geophysics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) March, 2013 c© James Johnson 2013 Abstract The primary focus of most reflection seismic surveys is to help locate hydro-carbon recourses. Due to an ever increasing scarcity of these re- courses, we must increase the size and quality of our seismic surveys. How- ever, processing such large seismic data volumes to accurately recover earth properties is a painstaking and computationally intensive process. Due to the way reflection seismic surveys are conducted there are often holes in the collected data, where traces are not recorded. This can be due to physical or cost constraints. For some of the initial stages of process- ing these missing traces are of little consequence. However processes like multiple prediction and removal, interferometric ground roll prediction, and migration require densely sampled data on a regular grid. Thus the need to interpolate undersampled data cannot be ignored. Using the fact that reflection seismic data sets obey a reciprocal relation- ship in source and receiver locations, combined with recent advances in the field of compressed sensing, we show that properly regularized the wavefield reconstruction problem can be solved with a high degree of accuracy. We exploit the compressible nature of seismic data in the curvelet domain to solve regularized `1 recovery problems that seek to match the measured data and enforce the above mentioned reciprocity. Using our method we were able to achieve results with a 20.45 dB sig- nal to noise ratio when reconstructing a marine data set that had 50% of its traces decimated. This is a 13.44 dB improvement over using the same method run without taking reciprocity into account. ii Preface This thesis was prepared with Madagascar, a reproducible research soft- ware package available at www.ahay.org. A large amount of code was developed in the Seismic Laboratory for Imaging and Modeling (SLIM). The numerical algorithms and applications are mainly written in Matlab. Early experiments were conducted using SLIMpy (slim.eos.ubc.ca/SLIMpy) a Python interface that exploits func- tionalities of seismic data processing packages, such as MADAGASCAR, through operator overloading. The main work was done using the SPOT and SPGL1 matlab packages developed by Dr. M.P. Friedlander and E.W. Van den Berg at UBC. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . x Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1 Principles of Compressive sensing . . . . . . . . . . . . . . . 6 2.1.1 Choice of the sparsifying transform . . . . . . . . . . 8 2.2 Source-receiver reciprocity . . . . . . . . . . . . . . . . . . . 11 2.2.1 Using reciprocity . . . . . . . . . . . . . . . . . . . . 13 3 Wavefield reconstruction . . . . . . . . . . . . . . . . . . . . . 17 3.1 Definition of the forward problem . . . . . . . . . . . . . . . 18 3.2 Recovery with perfect reciprocity . . . . . . . . . . . . . . . 21 3.2.1 Recovery from undersampled data . . . . . . . . . . . 21 3.2.2 Recovery with symmetric regularization . . . . . . . . 25 3.2.3 Projection to a symmetric subspace . . . . . . . . . . 29 3.3 Recovery with deviations from reciprocity . . . . . . . . . . . 33 3.3.1 Recovery from undersampled data . . . . . . . . . . . 37 3.3.2 Recovery with symmetric regularization . . . . . . . . 37 3.3.3 Projection to a symmetric subspace . . . . . . . . . . 37 3.4 Summary of results . . . . . . . . . . . . . . . . . . . . . . . 47 iv Table of Contents 4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.1 Using reciprocity during wavefield reconstruction . . . . . . . 49 4.2 Open and future research . . . . . . . . . . . . . . . . . . . . 50 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Appendices A Control parameter selection . . . . . . . . . . . . . . . . . . . 55 A.1 Marine data regular sampling . . . . . . . . . . . . . . . . . 56 A.2 Marine data jittered sampling . . . . . . . . . . . . . . . . . 57 A.3 Noisy marine data regular sampling . . . . . . . . . . . . . . 58 A.4 Noisy marine data jittered sampling . . . . . . . . . . . . . . 59 v List of Tables 3.1 Regularly undersamped marine data wavefield reconstruction results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2 Jitter sampled marine data wavefield reconstruction results . 48 3.3 Regularly undersampled noisy marine data wavefield recon- struction results . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.4 Jitter sampled noisy marine data wavefield reconstruction re- sults . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 A.1 Time slice from marine data, regular sampling . . . . . . . . 56 A.2 Time slice from marine data, jittered sampling . . . . . . . . 57 A.3 Time slice from noisy data, regular sampling . . . . . . . . . 58 A.4 Time slice from noisy data, jittered sampling . . . . . . . . . 59 vi List of Figures 1.1 Example showing a) the full data b) the same data under- sampled by a factor of 12 . . . . . . . . . . . . . . . . . . . . . 4 1.2 Reconstruction of the decimated data in 1.1 b, here a) is the data reconstructed without accounting for reciprocity, and b) is the data reconstructed using reciprocity . . . . . . . . . . . 5 2.1 Particulars of curvelet transform a) partitioning of Fourier space b) left curvelets in physical domain, right curvelets in Fourier domain, adapted from Hennenfent, 2008 . . . . . . . . 9 2.2 Alignment of curvelet coefficients with seismic wavefields, adapted from Hennenfent, 2008 . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Decay of seismic wavefield transform coefficient amplitudes, pink = curvelets, red = wavelets, blue = Fourier, adapted from Hennenfent, 2008 . . . . . . . . . . . . . . . . . . . . . . 10 2.4 Three pairs of reciprocal traces that have been overlaid onto each other, top near offsets, middle mid offsets, bottom far offsets adapted from Claerbout 2007 . . . . . . . . . . . . . . 12 2.5 Demonstration of data symmetry due to reciprocity on syn- thetic data. a) A time slice taken from a synthetic data set, b) its asymmetric decomposition, all zeros c) its symmetric decomposition, identical to a). . . . . . . . . . . . . . . . . . . 15 3.1 Different sampling schemes applied to noise free synthetic data. a) the total data for reference, b) regular sampling, c) jittered sampling. . . . . . . . . . . . . . . . . . . . . . . . 20 3.2 F-K spectra for shot gathers taken from a) the total data for reference, b) regular sampling, c) jittered sampling. . . . . . 20 3.3 Recovery from two fold regular undersampling, a) without using reciprocity, SNR 7.01 dB, b) difference between the re- covery and true data. . . . . . . . . . . . . . . . . . . . . . . . 22 vii List of Figures 3.4 Recovery from two fold jittered undersampling, a) without using reciprocity, SNR 9.24 dB, b) difference between the re- covery and true data. . . . . . . . . . . . . . . . . . . . . . . . 23 3.5 F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling, b) its difference from the true data, c) interpolated data from jittered undersampling, d) its difference from the true data. . . . . . . . . . . . . . . . . . . 24 3.6 Recovery from two fold regular undersampling, a) using sym- metric regulariztion, SNR 19.57 dB, b) difference between the recovery and true data. . . . . . . . . . . . . . . . . . . . . . 26 3.7 Recovery from two fold jittered undersampling, a) using sym- metric regularization, SNR 20.09 dB, b) difference between the recovery and true data. . . . . . . . . . . . . . . . . . . . 27 3.8 F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling, using symmetric regularization, b) its difference from the true data, c) interpolated data from jittered undersampling, using symmetric regularization, d) its difference from the true data. . . . . . . . . . . . . . . . . . . 28 3.9 Recovery from two fold regular undersampling, a) with the solution restricted to be symmetric, SNR 20.45 dB, b) differ- ence between the recovery and true data. . . . . . . . . . . . 30 3.10 Recovery from two fold jittered undersampling, a) with the solution restricted to be symmetric, SNR 20.86dB, b) differ- ence between the recovery and true data. . . . . . . . . . . . 31 3.11 F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling, using symmetric projection, b) its difference from the true data, c) interpolated data from jit- tered undersampling, using symmetric projection, d) its dif- ference from the true data. . . . . . . . . . . . . . . . . . . . 32 3.12 Marine data a) as originally acquired and processed b) after being degraded to a SNR of 11.81 c) Skew component of b). d) F-k spectra of b)’s central shot gather. . . . . . . . . . . . 34 3.13 Comparison of a trace from a) the original data, and the same trace after being degraded, and b) the degraded trace and its reciprocal counter part. . . . . . . . . . . . . . . . . . . . . . 35 3.14 Degraded data from previous figure after a) two fold regular undersampling and c) jittered undersampling. Plots b) and d) show the corresponding F-k spectra. . . . . . . . . . . . . . 36 viii List of Figures 3.15 Recovery from two fold regular undersampling of degraded data, a) without using reciprocity, SNR 5.92 dB, b) difference between the recovery and true data. . . . . . . . . . . . . . . 38 3.16 Recovery from two fold jittered undersampling of degraded, a) without using reciprocity, SNR 7.68 dB, b) difference between the recovery and true data. . . . . . . . . . . . . . . . . . . . 39 3.17 F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling of degraded data, b) its differ- ence from the true data, c) interpolated data from jittered undersampling of the same degraded data, d) its difference from the true data. . . . . . . . . . . . . . . . . . . . . . . . . 40 3.18 Recovery from two fold regular undersampling of degraded data, a) with symmetric regularization, SNR 12.13 dB, b) difference between the recovery and true data. . . . . . . . . . 41 3.19 Recovery from two fold jittered undersampling of degraded, a) with symmetric regularization, SNR 11.78 dB, b) difference between the recovery and true data. . . . . . . . . . . . . . . 42 3.20 F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling of degraded data using symmet- ric regularization, b) its difference from the true data, c) in- terpolated data from jittered undersampling of the same de- graded data using symmetric regularization, d) its difference from the true data. . . . . . . . . . . . . . . . . . . . . . . . . 43 3.21 Recovery from two fold regular undersampling of degraded data, a) using symmetric projection, SNR 12.03 dB, b) dif- ference between the recovery and true data. . . . . . . . . . . 44 3.22 Recovery from two fold jittered undersampling of degraded, a) using symmetric projection, SNR 11.77 dB, b) difference between the recovery and true data. . . . . . . . . . . . . . . 45 3.23 F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling of degraded data, using symmet- ric projection, b) its difference from the true data, c) interpo- lated data from jittered undersampling of the same degraded data, using symmetric projection, d) its difference from the true data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 ix Acknowledgments I would like to thank Dr. Felix Herrmann, my supervisor, for providing direction and insight throughout my program of studies. I would also like thank the other members of my committee Dr. Micheal Friedlander, and Dr. Micheal Bostock. All three provided a wealth of knowledge through out my course work and studies. I would also like to thank Dr. Eldad Haber for being one of the examiners for my defence. The entire SLIM group, especially Cody Brown, Sean Ross Ross, and Henryk Modzelewski provided their technical expertise, and support, with out which I could not have completed this research. Dr. Gilles Hennenfent also deserves acknowledgement for providing the foundation on which this work is built. Finally I would like to thank my friends and family for their tireless support and encouragement, which enabled me to complete this project. Thank you all. x Usus est magister optimus xi Chapter 1 Introduction Hydrocarbon resources are becoming increasingly scarce. The scarcity of these valuable resources has led to the drilling of reservoirs at greater depths than previously explored. Extracting oil and gas from reservoirs at such great depths requires high resolution images of the subsurface. However imaging the subsurface at these depths requires large volumes of high qual- ity seismic data. Unfortunately it is possible to spend millions of dollars on a seismic survey, and several months on processing the data without obtain- ing an accurate model of the subsurface, either through a lack of resolution due to not acquiring enough data, or the presence of noise. Many seismic data processing steps, such as migration, require densely sampled data on a regular grid. Some processes can require data that is sampled above the Shannon-Nyquist sampling rate. Conducting reflection seismic surveys is an expansive and labour intensive process. Billions of dol- lars are spent every year on seismic exploration projects. However even with such large expenditures, and intense effort, it is often not feasible to collect all of the information necessary to accurately image the subsurface. Due to monetary and geographic constraints surveys are often conducted with fewer sources and or receivers than required, and the placement of sources and receivers can be imprecise. This creates a situation where it is necessary to reconstruct the measured wavefield on a regular grid, interpolating the missing traces. Figure 1.1 shows a example of the type of under sampled shot records typically acquired in the field and a fully sampled wavefield required by other data processing techniques. The distribution of missing shots —i.e. the size and location of gaps in the shot record greatly affects the accuracy of the reconstruction processes. It has been shown that careful survey design can give a more accurate recon- struction from the same amount of input data (Hennenfent and Herrmann, 2008). Currently there exist many methods of reconstructing undersampled 1 Chapter 1. Introduction seismic wavefields (Abma and Kabir, 2006; Naghizadeh and Sacchi, 2009; Sacchi et al., 1998; Schonewille, 2000; Zwartjes, 2005; Trickett et al., 2010; Cary, 2011; Hollander et al., 2012; Stanton and Sacchi, 2011; Naghizadeh, 2012). Among these there are several methods which use the compressible nature of seismic wavefields in a transform domain and the emerging field of compressed sensing to accurately reconstruct the missing traces. We will examine one such method, and seek to improve upon it (Herrmann and Hen- nenfent, 2007; Hennenfent et al., 2010; Liu et al., 2012). The propagation of seismic wavefields in the earth’s interior is governed by the elastic wave equation. This equation obeys a reciprocal relationship between source and receiver locations (Knopoff and Gangi, 1959). Source re- ceiver reciprocal relationships in wave equations have previously been used in geophysics, in processing electromagnetic data, surface consistent deconvo- lution, migrating ocean bottom node data, and interferometry (Chan et al., 2005; van Vossen et al., 2006b; Alerini et al., 2009; Halliday and Curtis, 2010). There exists an opportunity to use source receiver reciprocity to aid in solving the wavefield reconstruction problem. We investigate two differ- ent methods of doing this. First, we regulate the inversion problem with a penalty term that enforces reciprocity in the processed data. Second, we restrict our solution to obey reciprocity explicitly. Both methods could be adapted to work with many of the interpolation algorithms cited above, however we adapt the previously mentioned sparsity promoting methods of processing seismic data to enforce reciprocity (Hennenfent and Herrmann, 2008). In the following chapters we show how these additions lead to a higher quality output of the inversion algorithm. In our sparsity promoting interpolation scheme we use a combination of a 2D curvelet transform along source/receiver location directions and wavelets along the time axis to sparsely represent seismic data (Candes et al., 2005). With the seismic wavefield transformed to the curvelet domain we attempt to solve an `1 constrained regularization problem. A recent focus on such sparse signal recovery problems, has led to the development of many al- gorithms to solve them. The algorithm we choose to employ is known as SPG`1 (Berg and Friedlander, 2007). Previous work has shown that seismic data processing problems are ami- able to being solved by this approach and that they can lead to accurate inversion results (Hennenfent, 2008; Herrmann et al., 2007; Lin and Her- rmann, 2009). These methods are an application of the emerging field of 2 Chapter 1. Introduction compressed sensing (Candes, 2006). Chapter 2 gives a more thorough explanation of the theory behind our method. In chapter 3 we examine the problem of wavefield reconstruction and show how our contributions affect the quality of the reconstructed wave- field. Figure 1.2 gives an example of this. Figure 1.2 a), is the data recon- structed from two fold regular undersampling with out taking into account reciprocity, figure 1.2 b) ias the same decimated wavefield reconstructed us- ing the same algorithm modified to take reciprocity into account. As the figure shows, using reciprocity leads to a more accurate reconstruction of the undersampled data. 3 Chapter 1. Introduction (a) (b) Figure 1.1: Example showing a) the full data b) the same data undersampled by a factor of 12 4 Chapter 1. Introduction (a) (b) Figure 1.2: Reconstruction of the decimated data in 1.1 b, here a) is the data reconstructed without accounting for reciprocity, and b) is the data reconstructed using reciprocity 5 Chapter 2 Theory We seek to improve the interpolation of undersampled seismic data onto a regular grid. This problem requires the inversion of large systems of equa- tions containing millions of variables, even for relatively small surveys. The accurate solution of such large and possibly under determined problems requires they be regularized. By this we mean the solution is forced to con- tain a priori knowledge. We seek to do this in two ways. First we use the emerging theory of compressed sensing to cast the problem as an sparsity promoting optimization problem (CS- Candes, 2006). Secondly we apply our knowledge of physics to force the recovered wavefields to obey the physical principle of reciprocity (Knopoff and Gangi, 1959). 2.1 Principles of Compressive sensing In noisy `1 optimization problems we seek to invert a operator A that acts on an unknown model vector x, to give the measured data b. In compressed sensing problems the operator A maps the larger model vector x to the much smaller data space b. We formulate the solution of this problem as x̃ = arg min x ‖x‖1 subject to ‖b−Ax‖ ≤ σ. (2.1) The parameter σ is related to the energy, of the noise present in b —i.e. ‖n‖ if b = b0 + n. The `1 norm in which we seek to minimize x is defined as ‖x‖1 = ∑n 1 |xi|, where n is the dimension of x. Compressed sensing assumes x to be sparse or compressible—i.e. few non zeros or coefficients with rapidly decaying sorted magnitudes. If this is true then we could attempt to solve the problem by minimizing the number of nonzero entries in x, however doing so is a NP hard problem (Candes, 2006). The theory of compressed sensing states that for certain classes of operator A, there is a very high probability that minimizing the `1 norm of x is a convex relaxation of this otherwise intractable problem (Candes, 6 Chapter 2. Theory 2006). This formulation is known as basis pursuit denoise BPDN. To solve problem 2.1 we need three ingredients. The first is a sparsity promoting inversion program. The output of the inversion program should match the given data within noise, when forward modelled. This refers to the inequality constraint in 2.1 (‖b−Ax‖ ≤ σ). Sparsity promotion relates to the minimization of x’s `1 norm (arg minx ‖x‖1). In recent years many different algorithms have been developed to do just this. Throughout this paper we use one such algorithm, SPG`1. This algorithm works by solving a series of related lasso problems x̃ = arg min x ‖b−Ax‖ subject to ‖x‖1 ≤ τ. (2.2) The solutions of these lasso problems have a minimal data misfit, Φ(τ) = ‖b−Ax‖, for a given `1-norm on x given by τ . After each lasso problem is solved the algorithm then uses a Newton method to try and find a solution to the equation Φ(τ) = σ (Berg and Friedlander, 2007). Ideally we would let the algorithm run to convergence, in practise this is too computationally expensive for problems of the size dealt with in seismic exploration. Instead we run the algorithm for a limited number of iterations. For noisy data in the case were we cannot run BPDN to convergence, the algorithm is typi- cally run with σ equal to zero this formulation known as basis pursuit BP. The other two ingredients of CS are a sparse representation (x) of the wavefield to be recovered (b), and an operator A that obeys the principles of CS. To ensure that we have these properties we factor the operator A into three parts A = RMSH. The operator SH is the adjoint of a sparsifying transform that ensures the desired compressibility of x. Our choice of S will be covered in the next section. The final ingredients of our CS framework are R and M. The choice of R and M is how we adapt the interpolation problem to the CS framework. In the context of wavefield reconstruction R is a subsampling matrix that represents data acquisition which leads to aliased measurements that need to be interpolated. The operator M is a modelling operator the choice of which can vary. In many cases M is an identity operator, but could equally be a migration or multiple predication operator (Lin and Herrmann, 2012). 7 Chapter 2. Theory 2.1.1 Choice of the sparsifying transform For an `1 recovery program to work we must be able to represent the wave- field that we are trying to recover in a sparse or compressible fashion. This means that we need to be able transform the wavefield to a domain where it has few non zero coefficients, or the coefficients have rapidly decaying ampli- tudes. Seismic wavefields have been shown to be highly compressible in the curvelet domain, —i.e. a complex wavefield can be accurately represented with a relatively small number of curvelets (Candes et al., 2005; Hennenfent and Herrmann, 2008). The curvelet transform decomposes data into multidimensional wave- forms that are formed by a second dyadic partitioning of Fourier space (Can- des et al., 2005). Figure 2.1 a) shows the partitioning of Fourier space and includes examples of a few different curvelet coefficients in both the physical and Fourier domains (Herrmann and Hennenfent, 2007). While individual curvelets are strictly localized in Fourier space, they rapidly decay in the physical domain. This is a function of the partitioning done in Fourier space, which prevents individual curvelets from having all of the frequencies neces- sary to be strictly localized in the physical domain. Figure 2.1 a) shows the 2nd dyadic partitioning of Fourier space into different scales. The coarsest scale is comprised of the lowest frequencies, and the finest scale composed of the highest. Each scale is further divided into angular wedges. At every other scale the number of angular wedges doubles, while the curvelets main- tain a parabolic scaling relationship length = width2 (Candes et al., 2005). Because curvelets are highly anisotropic, and oscillating in one direction and smooth in the other, the only significant coefficients will be those that are aligned with the wavefield. Figure 2.2 illustrates this concept (Herrmann and Hennenfent, 2007). The above properties make curvelets efficient at representing complex wavefields such as those encountered in seismic exploration. Figure 2.3 shows the rapid decay in the amplitude of curvelet coefficients for a ma- rine seismic line, compared to other possible transforms (Hennenfent, 2008; Herrmann and Hennenfent, 2007). Others have shown that curvelets behave as localized Eigen functions for high frequency approximations of the wave equation (Hennenfent et al., 2010), which helps explain how they parsimo- niously represent seismic data. As described above 2D curvelets have 8 times more coefficients than the 8 Chapter 2. Theory k1 k2 angular wedge 2j 2j/2 k1 k2 angular wedge 2j 2j/2 k1 k2 angular wedge 2j 2j/2 (a) -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 (b) Figure 2.1: Particulars of curvelet transform a) partitioning of Fourier space b) left curvelets in physical domain, right curvelets in Fourier domain, adapted from Hennenfent, 2008 data they represent, and 3D curvelets have approximately 24 times the num- ber of coefficients (Herrmann et al., 2009). Seismic data volumes are already quite large, often several Giga-bytes, and the increased redundancy makes using 3D curvelets computationally expensive. To get around this hurdle we use 2D curvelets (C) along the shot receiver location axis. We “Kronecker” the 2D curvelet transform with a non-redundant wavelet transform (W) in the time direction. We define the sparse transform as S ≡ C⊗W. (2.3) 9 Chapter 2. Theory 0 0.5 1.0 1.5 2.0 T im e ( s ) -2000 0 2000 Offset (m) 0 0.5 1.0 1.5 2.0 T im e ( s ) -2000 0 2000 Offset (m) Significant curvelet coefficient Curvelet coefficient~0 Figure 2.2: Alignment of curvelet coefficients with seismic wavefields, adapted from Hennenfent, 2008 Figure 2.3: Decay of seismic wavefield transform coefficient amplitudes, pink = curvelets, red = wavelets, blue = Fourier, adapted from Hennenfent, 2008 10 Chapter 2. Theory This choice in sparsifying transform avoids the 24 times redundancy of 3D curvelets and can still represent seismic data with a very low num- ber of significant coefficients. The combination of curvelets and wavelets also leads to embarrassingly parallel implementations where individual 2D curvelet time slices can be transformed and processed separately on different processing units. 2.2 Source-receiver reciprocity To aid our solution of the wavefield reconstruction problem, we seek to use our knowledge of how waves propagate through the earth. Namely we seek to exploit the reciprocal nature of the Green’s function when solving the minimization problem described above. The wave equation obeys a reciprocal relationship with respect to source and receiver locations. In mathematical terms we define the individual components of the Green’s function as gij(xs, xr, t). The source location is xs and xr is the receiver location. The indexes i and j correspond to the vector wavefield components at the source and receiver, and may be any of the principle coordinates. The source oscillates along the ith co- ordinate, the receiver records vibrations along the jth coordinate, and t signifies the time dependence of the Green’s function. Reciprocity makes gij(a, b, t) = gji(b, a, t). If we convolve (denoted by the symbol ∗) the source signature (qi(xs, t)) and receiver (rj(xr, t)) terms with the Greens function, then the measured trace pij(xs, xr, t) = qi(xs, t) ∗ gij(xs, xr, t) ∗ rj(xr, t) will also obey a reciprocal relationship. The reciprocal counterpart of pij(a, b, t) is pji(b, a, t) = qj(b, t) ∗ gji(b, a, t) ∗ ri(a, t). Reciprocity holds as long as qi(a, t) = qj(b, t) and rj(b, t) = ri(a, t). This means that for reciprocity to hold, the disturbance measured or created at a given location must be in the same direction. That is at location a we are always concerned with the perturbations along the ith coordinate, and at location b we are per- turbing/measuring along the jth coordinate. The above relationship is true for elastic media with arbitrary boundary conditions, inhomogeneities, and anisotropy (Knopoff and Gangi, 1959). For the rest of this paper we will assume single component data—i.e., we are only dealing with the vertical coordinate. This simplifies the notation, however similar results can be de- rived for multi-component data. 11 Chapter 2. Theory It may be difficult to expect reciprocity to hold exactly. It is an ideal- istic assumption for real seismic surveys. However it has been shown that reciprocity holds for reasonable departures from ideal acquisition geometries (Fenati and Rocca, 1984). Figure 2.4 shows 3 sets of reciprocal traces taken from a 2D survey. When the survey was conducted no specific effort was made to obey reciprocity but the three sets of traces show a remarkable amount of overlap (Claerbout, 2007). In practise differences of a few meters in lateral offset between source and receiver locations along a 2D line will not violate reciprocity, however differences on the order of station spacing (20-30m) will (Fenati and Rocca, 1984). Figure 2.4: Three pairs of reciprocal traces that have been overlaid onto each other, top near offsets, middle mid offsets, bottom far offsets adapted from Claerbout 2007 12 Chapter 2. Theory 2.2.1 Using reciprocity Maxwell’s equations also have a reciprocal relationship which has been ex- ploited in resistivity surveys (Chan et al., 2005). Currently source receiver reciprocity of seismic data is only used in a brute force fashion. In industry, if a missing trace’s reciprocal counterpart is available, it will typically be “borrowed” and used to fill in for the missing trace. This method of “bor- rowing” traces can be seen in processes like common depth point binning of seismic data. In many cases, reciprocity is used to generate the negative offsets necessary to make streamer data resemble data acquired with a split spread acquisition. Van Vossen proposed a method of exploiting reciprocity in surface con- sistent deconvolution. His method involves transforming data to the log- Fourier domain to separate the source receiver and Green’s function terms (van Vossen et al., 2006a). For this purpose Van Vossen creates an operator that strictly enforces ideal reciprocity —i.e, it only requires Green’s func- tions for half the shot locations. This operator is used to re-combine, by addition, the estimated source signatures, Green’s functions, and receiver responses before they are transformed back to the physical domain. Using the log-Fourier transform to separate terms can lead to stability problems as any noise in the data, or departures from true reciprocity, are exponentially amplified when the data is transformed, despite its utility in separating the source receiver and Green’s function terms. This method accounts for errors in field data caused by varying near surface conditions that lead to variation in source signature or receiver response along a seismic line. More recently Alerini has used reciprocity in the processing of ocean bottom node data (Alerini et al., 2009). Ocean bottom nodes are relatively expensive to install which leads to data which is often severely under sam- pled, 400 m spacing as opposed to a typical spacing of 25m. To get around this problem Alerini interpolates and re-datums the source locations to coin- cide with the receivers. Doing this allows Alerini to use the densely sampled common receiver gathers in place of the equivalent sparsely sampled com- mon shot gathers when performing stereo tomographic velocity estimation (Alerini et al., 2009). We propose two methods of exploiting reciprocity during the interpola- tion of under sampled data. The first method is a Tykhonov regularization of a method already proven to work on real seismic data. The second is to 13 Chapter 2. Theory directly restrict the solution of the inversion problem to obey reciprocity. Both methods rely on the decomposition of matrices into their symmetric and skew components. To explain these methods we will start with a single time slice taken from a 2D seismic line. We define a single time slice as Pt, a matrix organized with shot locations along the columns and receiver location for the rows. Let pt be the same data arranged as single vector where each column of the matrix has been appended to the end of the previous —i.e. pts = vec(Pt). We define an operator T2d that is a matrix transpose on each time slice Pt, T2dpt ≡ vec ( vec−1 ( pt )T ) , where the T transposes the shot and receiver axis. To extend our transpose operator to the whole data volume p we Kronecker it with an identity operator of dimensions equal to the number of time samples T = T2d⊗Int. Using T, data volumes can be decomposed into symmetric and skew components. This is done by adding or subtracting the transposed data from the original data and dividing by a factor of two. The symmetric component is psym = p + Tp 2 , (2.4) and the asymmetric component is pskew = p−Tp 2 . (2.5) Because of reciprocity, identical traces should be recorded when shot and receivers locations are interchanged. This makes shot records symmetrical about the shot reciever axis. Therefore psym = p and pskew = 0. This holds true for every time slice in a shot record. Figure 2.5 shows a time slice taken from a synthetic data set, and both its asymmetric and skew decompositions. The first way in which we will try to use reciprocity is to use the knowledge that a shot record’s asymmetric component should be zero — i.e. pskew = I−T 2 p = 0, to regularize the inversion process. Going back to the `1 recovery framework of equation 2.1, x̃ = arg min x ‖x‖1 subject to ‖b−Ax‖ ≤ σ, (2.6) we build in the asymmetry term by augmenting the operator A = RMSH. The forward problem remains RMSHx = p, with RM an as yet unspec- ified restriction and modelling operator. The regularization term is incor- porated into the inversion problem by vertically concatenating √ α I−T2 = 0 14 Chapter 2. Theory (a) (b) (c) Figure 2.5: Demonstration of data symmetry due to reciprocity on synthetic data. a) A time slice taken from a synthetic data set, b) its asymmetric decomposition, all zeros c) its symmetric decomposition, identical to a). to the forward problem, —i.e. A ≡ [ RM√ α ( I−T 2 )]SH, (2.7) The parameter α controls how much weight should be given to the asym- metric decomposition term. For a well-designed survey, a higher value can be used. The data vector b must also be augmented to account for the regularization term. We do this by appending zeros of length p to the end of the measured data, —i.e. b ≡ [ p 0 ] (2.8) The augmented system of equations makes 2.1 becomes x̃ = arg min x ‖x‖1 subject to ∥∥∥∥[p0 ] − [ RM√ α ( I−T 2 )]SHx∥∥∥∥ ≤ σ. (2.9) This method of regularization has the advantage of strongly informing the solver to make updates that obey reciprocity. This is because updates that do not obey reciprocity will not have an asymmetric decomposition equal to zero and will be penalized. We will show how this can lead to more accu- rate solutions of the inversion problem in half the iterations. We limit the algorithm to half the iterations of the unregularized solution to account for the increased workload in solving the augmented system of equations. The second way in which we seek to exploit reciprocity is to directly re- strict the solution to be symmetric as a part of the inversion problem, This is done by introducing a symmetric decomposition term (I + T) into the operator A such that A = RM (I + T) SH, and the solution in the physi- cal domain, p̃, becomes p̃ = (I + T) SHx̃. This method can be viewed as 15 Chapter 2. Theory projection of the data into a symmetric subspace. Due to the need for ex- tra matrix additions and transposes imposed by this method it will also be limited to half the iterations of the non-restricted or regularized solutions of equation 2.1. Restricting the solution to be symmetric like this makes equation 2.1 look like x̃ = arg min x ‖x‖1 subject to ∥∥p−RM (I + T) SHx∥∥ ≤ σ. (2.10) In the following chapter we will examine the performance of an algorithm that does not account for reciprocity versus the regularized and restricted methods of enforcing reciprocity. We solve the problem for multiple sampling schemes, R, and for data that does and does not strictly obey reciprocity. 16 Chapter 3 Wavefield reconstruction Due to a variety of factors it is rarely feasible to collect the required vol- ume of seismic data in the field, so data is often measured on a coarse grid, and the missing traces are reconstructed later as a part of data processing (Hennenfent et al., 2010). Accurate wavefield reconstruction is important for other steps in processing seismic data such as migration, interformetric ground roll prediction and multiple estimation. These algorithms tend to be sensitive to errors in the measured data and usually require regular sampling above the Shannon-Nyquist rate. This makes accurate wavefield reconstruc- tion an important and necessary component of seismic data processing. A higher quality reconstruction will lead to a higher quality final image. There are many different approaches to solving the wavefield reconstruc- tion problem (Sacchi et al., 1998; Schonewille, 2000; Zwartjes, 2005). These include F-x prediction filters (Naghizadeh and Sacchi, 2009), or projections onto convex sets (POCS Abma and Kabir, 2006), methods based on rank re- duction and the singular value decomposition (Trickett et al., 2010), and or- thogonal matching pursuit (Hollander et al., 2012). One such method of solv- ing the wavefield reconstruction problem is curvelet reconstruction through sparse inversion (Herrmann and Hennenfent, 2007,CRSI). This method uses the compressible nature of seismic wavefields in the curvelet domain and advances in the theory of compressed sensing to solve the underdetermined wavefield reconstruction inversion problem (Hennenfent, 2008). In this chapter we show how the wavefield reconstruction problem fits into the CS framework, and improve on previous results by using reciprocity induced symmetry of the data as a part of the inversion process. We will show that our contribution improves the quality of the wavefield estimate. 17 Chapter 3. Wavefield reconstruction 3.1 Definition of the forward problem In the previous chapter, we developed the CS framework of solving recovery problems. We write the problem in its standard form Ax = b (3.1) where b is the measured data, x is the vector which we seek to recover by inverting the operator A. As before, we factor the operator A into multiple parts, A = RMSH. The sparsifying transform S, which was previously de- fined, and RM is a sampling/modelling operator. In this section we further explain the structure of RM. To help explain the high dimensional structure of seismic wavefields, shot records from a 2D seismic line are arranged in a data cube where one axis is time, and the other two axises are the shot and receiver locations. We define the fully sampled data set arranged in such a fashion as P0, and its vectorized form as p0 = vec(P0). In the context of wavefield reconstruction we model the measured data p by applying the restriction or sub sampling operator, RM, to p0—i.e., p = RMp0. Our measurements also may be contaminated with noise, n, so that p = RMp0 + n. (3.2) For our problem we define M as the identity operator, M = IN×N , where N is the total number of measured data points, with N = Nt×Nr×Ns and Nt is the number of time samples, Nr is the number of receiver locations, and Ns is the number of shot locations. Our choice of M is a simplification and could easily be replaced by a more complex operator (Lin and Herrmann, 2009). The operator R “sub samples” or restricts the full data set to only those locations where traces have been recorded. To do this we define R as R = Rs ⊗Rr ⊗Rt. The operators Rt, Rr, and Rs, represent a restriction along the time, receiver and shot axis respectively. Modern geophones sam- ple wavefields above the Nyquist rate so we will let Rt = INt . Additionally we will assume full coverage of the receivers and define Rr = INr . This means that we are dealing with sub-sampling along the shot location axis —i.e., RM = Rs ⊗ INr ⊗ INt . (3.3) Here, Rs is a Ns′×Ns matrix, with Ns′ Ns. The operator Rs is composed of ones and zeros. In the case of full coverage of shots, this operator would 18 Chapter 3. Wavefield reconstruction be the identity and Ns′ = Ns. For each shot missing from the survey the corresponding row in Rs is removed. An example of Rs for every other shot missing would be Rs =  1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 . . .  . (3.4) Figure 3.1, shows two different realizations of R applied to a marine data set. For the purposes of display the subsampled data is displayed on the same Cartesian grid as the fully sampled data. Figure 3.1 a) is the total data set for reference b) shows the same dataset regularly under sampled by a factor of 2 —i.e., we only keep every second shot. Figure 3.1 c) contains the same number of measured data as the regular sampling scheme but instead of sampling on a regular grid, we use a jittered sampling scheme (Hennen- fent and Herrmann, 2008). Jittered means we used a sampling scheme that randomly perturbs the shot points around a regular grid. Results show that jittered sampling leads to a higher quality recovery than regular sampling or truly random sampling. Jittered sampling can outperform regular sampling by breaking the coherent aliases that would be introduced by a truly regular under sampling scheme, while limiting the maximal distance between adja- cent shots (Hennenfent and Herrmann, 2008). Figure 3.2 shows F-k spectra for shot gathers taken from all three data sets. In the fully sampled data shown in figure 3.1 a) there are 128 source locations, 128 receiver locations spaced 25 meters apart. The data is sampled at 0.04ms, and there are 256 samples on the time axises. With RM defined we now have all the components necessary to use our sparse recovery program x̃ = arg min x ‖x‖1 subject to ‖b−Ax‖ ≤ σ (3.5) to reconstruct seismic wavefields. The vector x̃ is a set of sparsity domain coefficients that when transformed to the physical domain p̃ = SHx̃ are fully sampled, and have a minimal misfit to the measured data p. In the following sections we will examine how the choice of RM determines the quality of the recovered data p̃, and what affect augmenting the above equations to account for reciprocity induced symmetry has on the recovery. 19 Chapter 3. Wavefield reconstruction (a) (b) (c) Figure 3.1: Different sampling schemes applied to noise free synthetic data. a) the total data for reference, b) regular sampling, c) jittered sampling. (a) (b) (c) Figure 3.2: F-K spectra for shot gathers taken from a) the total data for reference, b) regular sampling, c) jittered sampling. 20 Chapter 3. Wavefield reconstruction 3.2 Recovery with perfect reciprocity To understand the effect of different sampling schemes on the reconstructed wavefield, we solve the recovery problem for both sampling schemes shown in figure 3.1. To make it a fair comparison, both sampling schemes provide the same amount of data to the optimization program. The data set we use was acquired with a towed streamer array in the Gulf of Suez. The data has already had several pre-processing steps applied to it and reciprocity was used to turn the one sided data into the full data cube seen with a split spread acquisition (van Groenestijn and Verschuur, 2009). 3.2.1 Recovery from undersampled data We apply the unmodified algorithm described above to the Gulf of Suez data as a baseline for our experiments. In this case we are solving the problem defined as x̃ = arg min x ‖x‖1 subject to ∥∥p−RMSHx∥∥ ≤ σ. (3.6) Figure 3.1 shows a) the full wavefield, and b) the same wavefield when it has been undersampled with one out of every two traces missing in a regular pat- tern. Figure 3.3 a) shows the results of CRSI, after 500 hundred iterations in SGP`1, b) the difference between a) and the true data. Figure 3.5 a) and b) show F-k spectra for shots gathers taken from the recovered data, and the F-k spectra of difference between the recovered and the true data. The reconstructed wavefield has a signal to noise ratio, SNR = 20 log10 ( p0 p0−p̃ ) , of 7.01 db. Figure 3.1 c) shows the data two fold jitter undersampled. Figure 3.4 show the same plots as figure 3.3 when the wavefield is reconstructed from said jittered undesampling instead of regular spaced samples. Figure 3.5 c) and d) show the F-k spectra of shot gathers taken from the data in 3.4. Using this sampling scheme increases the recovered SNR to 9.24dB for the unmodified CRSI recovery with 500 iterations. 21 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.3: Recovery from two fold regular undersampling, a) without using reciprocity, SNR 7.01 dB, b) difference between the recovery and true data. 22 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.4: Recovery from two fold jittered undersampling, a) without using reciprocity, SNR 9.24 dB, b) difference between the recovery and true data. 23 Chapter 3. Wavefield reconstruction (a) (b) (c) (d) Figure 3.5: F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling, b) its difference from the true data, c) inter- polated data from jittered undersampling, d) its difference from the true data. 24 Chapter 3. Wavefield reconstruction 3.2.2 Recovery with symmetric regularization In the following section we show the results using the recovery program when it has been modified to include a Tykhonov regularization scheme that penalizes any asymmetry in the solution, i.e.– x̃ = arg min x ‖x‖1 subject to ∥∥∥∥[p0 ] − [ RM√ α ( I−T 2 )]SHx∥∥∥∥ ≤ σ (3.7) . Our regularization term introduces a new control parameter α which must be optimized. To do this we apply the algorithm repeatedly to indi- vidual time slices taken from the data, over a wide range of control param- eter values. From these experiments we select the value of α that produces the highest SNR. We use individual time slices to reduce the computational cost. Chapter 5 contains the tables of results generated from individual times slices. Due to the need for additional matrix-vec operations we limit SPG`1 to 250 iterations as opposed to the 500 used in the non-regularized results. Figure 3.6 shows the reconstructed wavefield using this approach, and its difference from the true data. Figure 3.8 a) and b) show the F-k plots of common shot gathers taken from both. The recovered wavefield has a SNR of 19.57 dB. This is a 12.56 dB improvement over the non-regularized results. Figure 3.7 show the results of the regularized recovery, and its differnce plots when using jittered sampling. The jittered results have a SNR of 20.09 dB, which is 10.85 more than the non-regularized results. 25 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.6: Recovery from two fold regular undersampling, a) using sym- metric regulariztion, SNR 19.57 dB, b) difference between the recovery and true data. 26 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.7: Recovery from two fold jittered undersampling, a) using sym- metric regularization, SNR 20.09 dB, b) difference between the recovery and true data. 27 Chapter 3. Wavefield reconstruction (a) (b) (c) (d) Figure 3.8: F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling, using symmetric regularization, b) its difference from the true data, c) interpolated data from jittered undersampling, using symmetric regularization, d) its difference from the true data. 28 Chapter 3. Wavefield reconstruction 3.2.3 Projection to a symmetric subspace As a final approach we cause the solution to be symmetric by including a symmetric projection in our recovery algorithm given below x̃ = arg min x ‖x‖1 subject to ∥∥p−RM (I + T) SHx∥∥ ≤ σ, (3.8) and the recovered data becomes p̃ = (I + T) SHx̃. Once again we limit SPG`1 to 250 iterations. Figure 3.9 shows the reconstructed wavefield when the solution is restricted to be symmetric in the physical domain and the difference between the original data and the interpolated data. The re- constructed wavefield using this method has a SNR of 20.45 dB, a 0.88 improvement over the regularized results and 13.44 dB better than the orig- inal algorithm given twice the number of iterations. Figure 3.10 shows the results of the restricted algorithm when applied to jitter undersampled data. The recovered data has a SNR 20.86 dB, 0.77 dB higher then the regularized results, and 11.62 dB the non-regularized results. The combination of jit- tered sampling and restriction to a symmetric subspace leads to the highest quality reconstructed wavefield. Figure 3.11 shows the F-K spectra of all the data shown in figures 3.9 and 3.10. Tables 3.1-3.2 summarize all of the results shown in this section. 29 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.9: Recovery from two fold regular undersampling, a) with the so- lution restricted to be symmetric, SNR 20.45 dB, b) difference between the recovery and true data. 30 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.10: Recovery from two fold jittered undersampling, a) with the solution restricted to be symmetric, SNR 20.86dB, b) difference between the recovery and true data. 31 Chapter 3. Wavefield reconstruction (a) (b) (c) (d) Figure 3.11: F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling, using symmetric projection, b) its difference from the true data, c) interpolated data from jittered undersampling, using symmetric projection, d) its difference from the true data. 32 Chapter 3. Wavefield reconstruction 3.3 Recovery with deviations from reciprocity The results in the previous section show that our contribution improves the SNR of the recovered wavefield when the input data is perfectly reciprocal. In practise, field data acquired with split spreads will never be acquired on a perfectly regular grid, and for land data differences in near surface conditions can lead to differences in the source/receiver signature in the data along a line (van Vossen et al., 2006a). For ocean bottom node data the sources and receivers are at different datums. To this end, we test our algorithm by artificially degrading the data set used in the previous section. First both sources and receivers were convolved with a damped harmonic oscillator to model differences in the ”coupling” of the sources and receivers. The damped harmonic coupling terms frequency response can be written as C(f) = − [ 1 + ( ffc )ηc ] [ 1− ( ffs )2 + i( f fc )ηc ] , (3.9) where C(f) is the coupling term, ηc is the damping factor, and fc is the resonant frequency. The values of fc, and ηc were chosen from a random distribution with mean and standard deviation given below: fc (Hz) ηc µ 150 1 σ 35 0.2 A different value of fc, and ηc where chosen for each source and receiver location, and the resulting damped harmonic term convolved with each trace corresponding to that source or receiver. Gaussian white random noise with zero mean and a standard deviation equal to 0.05 the data’s standard devi- ation, was also generated, this noise had a 100 Hz low pass filter applied to it before it was added to the degraded data. This results in a SNR of 11.81 for the noisy degraded data. Figure 3.12 shows the original data, the data after it has been noised, the degraded data’s skew decomposition and its F-k spectra. For a visual comparison of the differnces caused by degrading the data, Figure ?? a) shows the same trace from the original data overlane with the degraded version. Figure 3.13 b) shows a comparison of the same noisy trace and its reciprocal counterpart. Figure 3.14 shows the degraded data after being two fold undersampled using both a regular and jittered sampling scheme, and the corresponding F-k spectrum’s. 33 Chapter 3. Wavefield reconstruction (a) (b) (c) (d) Figure 3.12: Marine data a) as originally acquired and processed b) after being degraded to a SNR of 11.81 c) Skew component of b). d) F-k spectra of b)’s central shot gather. 34 Chapter 3. Wavefield reconstruction 0 0.2 0.4 0.6 0.8 1 1.2 1.4 −400 −300 −200 −100 0 100 200 300 400 time (s) (a) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 −400 −300 −200 −100 0 100 200 300 400 time (s) (b) Figure 3.13: Comparison of a trace from a) the original data, and the same trace after being degraded, and b) the degraded trace and its reciprocal counter part. 35 Chapter 3. Wavefield reconstruction (a) (b) (c) (d) Figure 3.14: Degraded data from previous figure after a) two fold regular undersampling and c) jittered undersampling. Plots b) and d) show the corresponding F-k spectra. 36 Chapter 3. Wavefield reconstruction 3.3.1 Recovery from undersampled data As in the previous section we start by solving the unregularized recovery problem for 500 iterations in SPG`1 to set a baseline for our recovery. Fig- ure 3.15 a) shows the results of the unregularized inversion, and b) the difference between the recovered and true data for regularly undersampled and degraded data. The reconstructed wavefield has a SNR of 5.92 dB. We also tested the unregularized recovery on jittered sampled data. The results are shown in Figure 3.16. The interpolated data has a SNR of 7.68 dB. The regularly sampled data has a SNR that is 1.09 dB lower then the recovery when we used noise free data, while the recovery from jittered sampling is 1.56 dB lower then the noise free data. Figure 3.23 shows F-k spectra plots of shot gathers taken from the data in figures 3.21 and 3.22. 3.3.2 Recovery with symmetric regularization We applied the algorithm given in section 3.2.2 to the degraded data. Figure 3.18 a) shows the data that was recovered after 250 iterations in SPG`1, figure 3.18 b) shows the difference between the recovered data in a) and the noise free data . The interpolated data has a SNR of 12.13 dB. Our algorithm has simultaneously interpolated the missing traces and removed some of the noise that was introduced. However the recovery is still 7.44 dB lower then the noise-free ideal recovery. Using our regularization term in the recovery leads to results that have a 6.21 dB higher SNR than not using it. The regularized results for jittered sampling are shown in figure 3.19 a), and the difference with the true data is shown in figure 3.19 b). The recovered SNR is 11.78 dB, almost back to the noise level, and 4.1 dB higher than the non-regularized results. Figure 3.17 shows F-k spectra plots of shot gathers taken from the data in figures 3.18 and 3.19. 3.3.3 Projection to a symmetric subspace As a final test we used the algorithm described in 3.2.3, on the degraded noisy data. Figure 3.21 a) shows the regularly undersampled data that has been interpolated with a symmetric restriction and b) its difference to the true data. The recovered data has a SNR of 12.03 dB. This is 0.1 dB lower then the regularized results. Figure 3.22 a) shows the restricted interpolation of noisy jitter sampled data and b) its difference from the true data. The recovered data has a SNR of 11.77 dB, 0.01 dB lower than the regularized results. Figure 3.23 shows F-k spectra plots of shot gathers taken from the data in figures 3.21 and 3.22. 37 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.15: Recovery from two fold regular undersampling of degraded data, a) without using reciprocity, SNR 5.92 dB, b) difference between the recovery and true data. 38 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.16: Recovery from two fold jittered undersampling of degraded, a) without using reciprocity, SNR 7.68 dB, b) difference between the recovery and true data. 39 Chapter 3. Wavefield reconstruction (a) (b) (c) (d) Figure 3.17: F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling of degraded data, b) its difference from the true data, c) interpolated data from jittered undersampling of the same degraded data, d) its difference from the true data. 40 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.18: Recovery from two fold regular undersampling of degraded data, a) with symmetric regularization, SNR 12.13 dB, b) difference between the recovery and true data. 41 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.19: Recovery from two fold jittered undersampling of degraded, a) with symmetric regularization, SNR 11.78 dB, b) difference between the recovery and true data. 42 Chapter 3. Wavefield reconstruction (a) (b) (c) (d) Figure 3.20: F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling of degraded data using symmetric regulariza- tion, b) its difference from the true data, c) interpolated data from jittered undersampling of the same degraded data using symmetric regularization, d) its difference from the true data. 43 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.21: Recovery from two fold regular undersampling of degraded data, a) using symmetric projection, SNR 12.03 dB, b) difference between the recovery and true data. 44 Chapter 3. Wavefield reconstruction (a) (b) Figure 3.22: Recovery from two fold jittered undersampling of degraded, a) using symmetric projection, SNR 11.77 dB, b) difference between the recovery and true data. 45 Chapter 3. Wavefield reconstruction (a) (b) (c) (d) Figure 3.23: F-k spectra for shot gathers taken from a) interpolated data, from regular undersampling of degraded data, using symmetric projection, b) its difference from the true data, c) interpolated data from jittered un- dersampling of the same degraded data, using symmetric projection, d) its difference from the true data. 46 Chapter 3. Wavefield reconstruction 3.4 Summary of results Tables 3.1-4 summarize the results of our experiments. The first data col- umn shows the number of iterations that SPG`1 was allowed to run for, the second the signal to noise ratio of the recovery, and the third column shows the 2-norm of the skew component of the recovered results, as a fraction of the true data’s 2-norm. For data that perfectly obeys reciprocity, restricting the data so that it must be perfectly symmetrical, as required by reciprocity, leads to the high- est quality results. When the quality of the input data, was degraded to match what might be acquired with a real split spread acquisition the reg- ularization method leads to the highest quality results. By penalizing any asymmetry, rather then strictly requiring symmetry of the output data this algorithm is able to strike a balance between honouring the input data and enforcing reciprocity. We believe this makes the regularization method more robust than the restricted algorithm with respect to interpolating noisy data. However further tests on real data acquired in the field with split spreads would be required to verify this fact. By using the reciprocity induced structure of the data as a priori infor- mation during the reconstruction of missing traces, either through a penalty term, or restricting the data to be truly symmetric we have increased the signal to noise ratio of the interpolated data by as much as 13.44 dB over existing methods. 47 Chapter 3. Wavefield reconstruction Table 3.1: Regularly undersamped marine data wavefield reconstruction results Reciprocity iterations SNR (dB) ‖ p̃−T p̃ 2 ‖2 ‖p0‖2 None 500 7.01 0.2313 Restriction 250 20.45 0 Regularization 250 19.57 0.0098 Table 3.2: Jitter sampled marine data wavefield reconstruction results Reciprocity iterations SNR (dB) ‖ p̃−T p̃ 2 ‖2 ‖p0‖2 None 500 9.24 0.1866 Restriction 250 20.86 0 Regularization 250 20.09 0.0088 Table 3.3: Regularly undersampled noisy marine data wavefield reconstruc- tion results Reciprocity iterations SNR (dB) ‖ p̃−T p̃ 2 ‖2 ‖p0‖2 None 500 5.92 02848 Restriction 250 12.03 0 Regularization 250 12.13 0.0166 Table 3.4: Jitter sampled noisy marine data wavefield reconstruction results Reciprocity iterations SNR (dB) ‖ p̃−T p̃ 2 ‖2 ‖p0‖2 None 500 7.68 0.2374 Restriction 250 11.77 0 Regularization 250 11.78 0.0190 48 Chapter 4 Conclusions The main purpose of this thesis is to improve upon an existing algorithm for interpolating seismic data. Recognizing that reciprocity makes pre-stack 2D seismic data symmetric about the shot and receiver location axis, we de- veloped two methods of exploiting this reciprocity induced symmetry. The first method we use is a Tykhonov regularization term that penalizes any asymmetry in the data. The other method is to strictly enforce the reci- procity induced symmetry of the data. We applied these methods to two different data sets, using two different sampling schemes, and found that both methods lead to a higher SNR of the recovered wavefield. Our work is based on recent advances in using curvelet frames and the theory of compressed sensing to process seismic data (Can- des et al., 2005; Candes, 2006; Herrmann and Hennenfent, 2007; Hennenfent, 2008). 4.1 Using reciprocity during wavefield reconstruction Previous work on the reconstruction of missing seismic traces has been im- proved upon by exploiting the physics of wave propagation through the earth. Namely reciprocity induced symmetry of the data has been used to increase the quality of the inversion results. Two methods of doing so were investigated. The first is a Tykhonov regularization of an existing method, the second is to strictly enforce reciprocity induced symmetry of the data as a part of the same method. For two-fold under sampled data that is truly symmetric both methods increase the SNR of the reconstructed wavefield by more then 10 dB, in half the iterations of the solver. In these tests not accounting for reciprocity resulted in a SNR of 7.01 dB, for regular under sampling and 9.24 for jit- 49 Chapter 4. Conclusions tered spacing of source locations. The regularization method results in a SNR of 19.57 and 20.09 dB respectively for the different sampling schemes. The restriction method led to results of 20.45 and 20.86 dB for each of the sampling methods. The same data set was then degraded by convolving randomized source and receiver terms that vary along the line and adding filtered gaussian white noise. This new dataset has a SNR of 11.63 compared to the true dataset. For this dataset not using reciprocity during the inversion resulted in a recovery of 5.92 dB for regularly under sampled data and 7.68 dB for jitter sampled data. The regularization method resulted in a SNR of 12.13 dB for the regularly sampled data, and 11.78 dB for the jitter sampling. The restriction method has a SNR of 12.03 dB for regularly missing shots, and 11.77 for jitter sampled data. Tables 3.2-5 in the previous chapter sum- marize these results. While both methods improve the quality of the results, the regularization method appears to be more robust with respect to any variations from true reciprocity in the data. We believe that the higher SNR of the recovered wavefield for both regularization and restriction, of the regularly sampled noisy data set compared to the jitter sampled data set is due to the inter- action of the symmetry terms with the noise and the slightly larger gaps in the data caused by the jitter sampling scheme. 4.2 Open and future research We have shown how our contributions lead to a higher quality inversion re- sults when reconstructing seismic wavefields however there are still further areas of enquiry which could be investigated. A study of the control pa- rameter sensitivity in the regularization method would be beneficial. Also our choice of sparsifying transform could be further examined. A combined 2D curvelet Kroneckered with a wavelet transform, is less memory inten- sive than 3D curvelets and easier to implement. However the use of 3D curvelet’s transform has been shown to provide a uplift in SNR of the recov- ered wavefield (Hennenfent, 2008). This is because 3D curvelets can make better use of the higher dimensional structure of seismic wavefields. Our method could also be adapted to other interpolation algorithms. Further- more the work of van Vossen to account for variations in source signature and receiver characteristics could be included in the interpolation scheme, 50 Chapter 4. Conclusions along with enforcement of reciprocity (van Vossen et al., 2006a). This ap- proach could even be combined with recent advances in primary multiple separation which solve directly for the Green’s function, namely estimation of primaries via sparse inversion(van Groenestijn and Verschuur, 2009; Lin and Herrmann, 2009). This would result in a combined designature, inter- polation, and multiple removal algorithm. 51 Bibliography Abma, R., and N. Kabir, 2006, 3d interpolation of irregular data with a pocs algorithm: Geophysics, 71, E91–E97. Alerini, M., B. Traub, C. Ravaut, and E. Duveneck, 2009, Prestack depth imaging of ocean-bottom node data: Geophysics, 74. Berg, E. v., and M. P. Friedlander, 2007, SPGL1: A solver for large-scale sparse reconstruction. (http://www.cs.ubc.ca/labs/scl/spgl1). Candes, E. J., 2006, Compressive sampling. Candes, E. J., D. L. Donoho, and L. Ying, 2005, Fast discrete curvelet transforms: Multiscale Modeling and Simulation, 5, 861–899. Cary, P. W., 2011, Aliasing and 5d interpolation with the mwni algorithm: Presented at the SEG San Antonio 2011 Annual Meeting, SEG. Chan, J., D. W. Oldenburg, and E. Haber, 2005, Reciprocity in electro- magnetics: Application to marine magnetometric resistivity: Physics of the earth and Planetary Interiors, 45–61. Claerbout, J., 2007, Basic earth imaging. Fenati, D., and F. Rocca, 1984, Seismic reciprocity field tests from the italian peninsula: Geophysics, 49, 1690–1700. Halliday, D., and A. Curtis, 2010, An interferometric theory of source- receiver scattering and imaging: Geophysics, 75. Hennenfent, G., 2008, Sampling and reconstruction of seismic wavefields in the curvelet domain: PhD thesis, The University of British Columbia, Vancouver, BC Canada. 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 re- construction via jittered undersampling: Geophysics, 73, no. 3. 52 Bibliography Herrmann, F., D. Wang, G. Hennenfent, and P. Moghaddam, 2007, Curvelet-based seismic data processing: a multiscale and nonlinear ap- proach: Geophysics, 73, A1–A5. Herrmann, F. J., and G. Hennenfent, 2007, Non-parametric seismic data recovery with curvelet frames: Technical report, UBC Earth & Ocean Sci- ences Department. (TR-2007-3). Herrmann, F. J., T. Lin, and Y. A. Erlangga, 2009, Compressive simul- taneous full-waveform simulation: SEG Technical Program Expanded Ab- stracts, SEG, SEG, 2577. Hollander, Y., D. Kosloff, Z. Koren, and A. Bartana, 2012, Seismic data interpolation by orthogonal matching pursuit: Presented at the 74th EAGE Conference and Exhibition, EAGE. Knopoff, L., and A. F. Gangi, 1959, Seismic reciprocity: Geophysics, 24, 681–691. Lin, T., and F. J. Herrmann, 2009, Unified compressive sensing framework for simultaneous acquisition with primary estimation: SEG Technical Pro- gram Expanded Abstracts, SEG, 3113–3117. Lin, T. T., and F. J. Herrmann, 2012, Robust estimation of primaries by sparse inversion via one norm minimization: Presented at the submitted. Liu, Y., C. Liu, D. Wang, X. Feng, and Q. Lu, 2012, Iterative seismic data interpolation beyond aliasing using seislet transform: Presented at the 74th EAGE Conference and Exhibition, EAGE. Naghizadeh, M., 2012, Seismic data interpolation and denoising in the frequency-wavenumber domain: Geophysics, 77. Naghizadeh, M., and M. Sacchi, 2009, F-x adaptive seismic-trace interpo- lation: Geophysics, 74, V9–V16. Sacchi, M., T. Ulrych, and C. Walker, 1998, Interpolation and extrapola- tion using a high-resolution discrete fourier transform: IEEE Transactions on Signal Processing, 46, 31–38. Schonewille, M., 2000, Fourier reconstruction of irregularly sampled seismic data: PhD thesis, Delft University of Technology, Delft, The Netherlands. Stanton, A., and M. Sacchi, 2011, Multicomponent seismic data recon- struction using the quaternion fourier transform and pocs: Presented at the SEG San Antonio 2011 Annual Meeting, SEG. 53 Bibliography Trickett, S., S. L. Burroughs, A. Milton, L. Walton, and R. Dack, 2010, Rank-reduction-based trace interpolation: 81st Annual International Meet- ing, SEG, 1989–1992. van Groenestijn, G. J. A., and D. J. Verschuur, 2009, Estimating primaries by sparse inversion and application to near-offset data reconstruction: Geo- physics, 74, A23–A28. van Vossen, R., A. Curtis, A. Laake, and J. Trampert, 2006a, Surface- consistent deconvolution using reciprocity and waveform inversion: Geo- physics, 71, V19–V30. van Vossen, R., A. Curtis, and J. Trampert, 2006b, Surface-consistent am- plitude corrections for single of multicomponent sources and receivers us- ing reciprocity and waveform inversion: Geophysics Journal International, 165, 311–322. Zwartjes, P., 2005, Fourier reconstruction with sparse inversion: PhD the- sis, Delft University of Technology, Delft, The Netherlands. 54 Appendix A Control parameter selection α = 0 means no regularization of the inversion problem —i.e. not taking into account the symmetry that reciprocity imposes on the data. SNR = 20 ∗ log10 ( ‖P0‖ ‖P0−P̃‖ ) . All problems were run with 5000 maximum iterations and σ = ‖p0‖ ∗ 10−4 55 Appendix A. Control parameter selection A.1 Marine data regular sampling Table A.1: Time slice from marine data, regular sampling α iterations ∥∥P̃−Po∥∥ ∥∥ P̃−TP̃2 ∥∥ SNR (db) 0 5,000 1,805 928.9 5.35 0.10 5,000 442.1 1.04 17.57 0.12 5,000 433.8 0.76 17.75 0.15 5,000 366.3 0.45 19.21 0.18 5,000 369.8 0.12 19.13 0.22 5,000 350.3 0.08 19.60 0.26 5,000 335.9 0.02 19.97 0.32 5,000 331.6 0.00 20.08 0.38 5,000 328.9 0.00 20.15 0.46 5,000 329.8 0.00 20.13 0.56 4,327 329.3 0.00 20.14 0.68 3,285 327.8 0.00 20.18 0.83 3,199 328.1 0.00 20.17 1.00 3,307 327.1 0.00 20.20 1.21 2,973 326.1 0.00 20.22 1.47 5,000 327.0 0.27 20.20 1.78 4,228 328.3 0.00 20.16 2.15 5,000 330.0 0.06 20.11 2.61 5,000 326.3 0.08 20.22 3.16 5,000 328.1 0.00 20.17 3.83 4,954 329.4 0.00 20.14 4.64 5,000 329.9 0.00 20.12 5.62 5,000 328.5 0.03 20.16 6.81 5,000 334.5 0.02 20.00 8.25 5,000 340.4 0.03 19.85 10.00 5,000 351.9 0.00 19.56 56 Appendix A. Control parameter selection A.2 Marine data jittered sampling Table A.2: Time slice from marine data, jittered sampling α iterations ∥∥P̃−Po∥∥ ∥∥ P̃−TP̃2 ∥∥ SNR (db) 0.00 1,792.00 1,604.00 855.86 6.39 0.10 5,000.00 446.45 1.81 17.50 0.12 5,000.00 437.58 0.18 17.67 0.15 5,000.00 385.11 0.10 18.78 0.18 5,000.00 371.42 0.15 19.09 0.22 5,000.00 359.05 0.06 19.39 0.26 5,000.00 357.05 0.01 19.44 0.32 5,000.00 358.38 0.01 19.40 0.38 5,000.00 353.04 0.00 19.53 0.46 4,987.00 351.36 0.00 19.58 0.56 4,883.00 350.50 0.00 19.60 0.68 3,811.00 350.30 0.00 19.60 0.83 3,463.00 348.30 0.00 19.65 1.00 3,522.00 347.45 0.00 19.67 1.21 2,922.00 350.80 0.00 19.59 1.47 5,000.00 350.27 0.10 19.60 1.78 4,499.00 352.47 0.00 19.55 2.15 5,000.00 353.00 0.01 19.54 2.61 5,000.00 355.29 0.02 19.48 3.16 5,000.00 357.84 0.00 19.42 3.83 5,000.00 358.07 0.00 19.41 4.64 5,000.00 360.94 0.00 19.34 5.62 5,000.00 360.71 0.04 19.35 6.81 5,000.00 371.70 0.01 19.09 8.25 5,000.00 373.95 0.03 19.03 10.00 5,000.00 383.17 0.01 18.82 57 Appendix A. Control parameter selection A.3 Noisy marine data regular sampling Table A.3: Time slice from noisy data, regular sampling α iterations ∥∥P̃−Po∥∥ ∥∥ P̃−TP̃2 ∥∥ SNR (db) 0.00 2,499.00 2,050.51 1,126.68 4.38 0.10 2,188.00 993.97 117.96 10.54 0.12 1,869.00 953.93 117.42 10.90 0.15 1,518.00 933.01 116.63 11.09 0.18 1,423.00 915.63 115.49 11.26 0.22 816.00 923.33 113.86 11.18 0.26 855.00 913.10 111.54 11.28 0.32 492.00 913.88 108.31 11.27 0.38 389.00 912.48 103.89 11.29 0.46 394.00 910.43 98.02 11.31 0.56 273.00 909.58 90.52 11.31 0.68 249.00 908.89 81.37 11.32 0.83 251.00 908.03 70.86 11.33 1.00 233.00 906.77 59.57 11.34 1.21 275.00 906.75 48.28 11.34 1.47 289.00 904.87 37.77 11.36 1.78 426.00 906.38 28.62 11.34 2.15 423.00 904.94 21.12 11.36 2.61 414.00 903.48 15.25 11.37 3.16 437.00 902.20 10.83 11.38 3.83 558.00 904.55 7.60 11.36 4.64 753.00 901.73 5.28 11.39 5.62 1,947.00 900.66 3.65 11.40 6.81 1,867.00 901.85 2.51 11.39 8.25 2,417.00 903.37 1.72 11.37 10.00 2,040.00 908.70 1.18 11.32 58 Appendix A. Control parameter selection A.4 Noisy marine data jittered sampling Table A.4: Time slice from noisy data, jittered sampling α iterations ∥∥P̃−Po∥∥ ∥∥ P̃−TP̃2 ∥∥ SNR (db) 0.00 2,284.00 1,843.62 1,062.45 5.18 0.10 2,152.00 1,008.27 141.39 10.42 0.12 1,804.00 975.45 140.74 10.71 0.15 1,527.00 961.16 139.79 10.83 0.18 1,230.00 957.99 138.42 10.86 0.22 926.00 957.58 136.47 10.87 0.26 877.00 955.36 133.69 10.89 0.32 485.00 957.91 129.82 10.86 0.38 462.00 951.65 124.52 10.92 0.46 325.00 950.86 117.49 10.93 0.56 345.00 950.48 108.49 10.93 0.68 296.00 948.57 97.53 10.95 0.83 311.00 946.75 84.94 10.97 1.00 227.00 946.28 71.40 10.97 1.21 244.00 945.67 57.87 10.98 1.47 253.00 943.42 45.27 11.00 1.78 249.00 943.98 34.31 10.99 2.15 449.00 944.87 25.31 10.98 2.61 412.00 944.81 18.28 10.98 3.16 499.00 944.19 12.98 10.99 3.83 663.00 943.64 9.11 10.99 4.64 746.00 942.66 6.33 11.00 5.62 1,945.00 945.82 4.38 10.97 6.81 1,719.00 947.00 3.01 10.96 8.25 2,218.00 947.90 2.07 10.96 10.00 2,361.00 949.05 1.41 10.95 59
Thesis/Dissertation
eng
Geophysics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Seismic wavefield reconstruction using reciprocity
Text