A N I N V E R S E S C A T T E R I N G SERIES M E T H O D F O R A T T E N U A T I N G ELASTIC MULTIPLES F R O M MULTICOMPONENT L A N D A N D O C E A N B O T T O M SEISMIC D A T A by Kenneth Howell Matson B.Sc. (Geophysics), T h e University of Alberta, 1989 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF EARTH AND OCEAN SCIENCES We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA April © 1997 Kenneth Howell Matson, 1997 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of E a r t h and Ocean Sciences T h e University of British Columbia 129-2219 M a i n M a l l Vancouver, Canada V 6 T 1Z4 Date: Abstract A method exists for marine seismic data which removes all orders of free surface multiples and suppresses all orders of internal multiples while leaving primaries intact. This method is based on the inverse scattering series and makes no assumptions about the subsurface earth model. T h e marine algorithm assumes that the sources and receivers are located in the water column. In the context of land and ocean bottom data, the sources and receivers are located on or in an elastic medium. This opens up the possibility of recording multicomponent seismic data. Because both compressional ( P ) and shear (S) primaries are recorded in multicomponent data, it has the potential for providing a more complete picture of the subsurface. Coupled with the benefits of the P and S primaries are a complex set of elastic free surface and internal multiples. In this thesis, I develop an inverse scattering series method to attenuate these elastic multiples from multicomponent land and ocean bottom data. For land data, this method removes elastic free surface multiples. For ocean bottom data, multiples associated with the top and bottom of the water column are removed. Internal multiples are strongly attenuated for both data types. In common with the marine formulation, this method makes no assumptions about the earth below the sources and receivers, and does not affect primaries. T h e latter property is important for amplitude variation with offset analysis ( A V O ) . T h e theory for multiple attentuation requires four component (two source, two receiver) data, a known near surface or water bottom, near offsets, and a known source wavelet. Tests on synthetic data indicate that this method is still effective using data with less than four components and is robust with respect to errors in estimating the near surface or ocean bottom properties. n Contents Abstract ii Table of Contents iii List of Tables vi List of Figures vii Acknowledgement x 1 Introduction 1 2 Forward scattering series 8 2.1 Introduction 8 2.2 T h e B o r n series 10 2.3 T w o dimensional scattering from one dimensional models 14 2.3.1 16 2.3.2 2.4 3 Scattering by a half space Scattering by a layer over a half space 24 Conclusions 31 Free Surface Multiple Removal for Land Data 34 3.1 Introduction 34 3.2 Forward and inverse elastic wave scattering 34 3.2.1 3.3 Homogeneous isotropic reference medium and (P,S) potentials Free surface (P,S) Green's function . . . 41 47 in 4 3.4 Elastic free surface multiple removal subseries 52 3.5 Free surface (P,S) decomposition and ghost removal 59 3.6 Incorporating wavelets and source orientation 61 3.7 Synthetic data examples 64 3.7.1 Model 3 64 3.7.2 Model 4 74 3.8 Multiple removal for displacements 78 3.9 Discussion and conclusions 81 Multiple Removal for Ocean Bottom Data 83 4.1 Introduction 83 4.2 (P,S) Green's function for ocean bottom data 84 4.3 Water column multiple removal subseries 88 4.4 (P,S) decomposition and up-down going wavefield separation 92 4.4.1 95 4.5 5 Synthetic data examples 96 4.5.1 Model 5 96 4.5.2 Model 6 102 Elastic internal multiple suppression 108 5.1 Introduction 108 5.2 Review of 1-D internal multiple suppression 109 5.3 Theory for elastic internal multiple suppression 113 5.3.1 120 5.4 6 Removing the reference wavefield Synthetic data example Discussion and conclusions 121 Conclusions 129 Bibliography 134 A Calculated series for 1-D Born scattering 140 iv B C Elastic reflection and transmission coefficients 144 (P,S) decomposition and up-down wave separation 148 C.l L a n d data 148 C. 2 Receiver decomposition for ocean bottom data 151 D Elastic interface multiple removal series 153 D. l Free surface multiple removal 157 D.2 Ocean bottom multiple removal 166 v List of Tables 6.1 D a t a type vs. required multiple removal informat vi List of Figures 1.1 Primary and multiple reflections 2 1.2 Synthetic shot record containing primaries and multiples 2 2.1 Scattering view of reflection from an interface 9 2.2 Model 1: acoustic plane waves incident upon planar interface 2.3 Comparison of the true reflection coefficient to the B o r n approximate reflec- 16 tion coefficient 19 2.4 Plane wave reflective scattering from an interface 21 2.5 Plane wave transmissive scattering from an interface 23 2.6 Convergence criteria for simple scattering problem 24 2.7 Model 2: acoustic plane waves incident upon a layer over a half space 2.8 Order of scattering terms vs. event type 28 2.9 B o r n series first primary for a plane incident wave 29 2.10 B o r n series second primary 30 2.11 Born series first internal multiple 31 3.1 Scattering in displacements (u , u ) 3.2 Scattering of ((f>, tp) potentials 3.3 Reflected and converted waves for a P wave incident on a free surface 3.4 Reflected and converted waves for an S wave incident on a free surface 3.5 Flowchart of steps in elastic free surface multiple removal 56 3.6 Subevent construction of free surface multiples 58 3.7 L a n d data Model 3 64 x . . . 25 44 z 45 vn . . . . . 49 50 3.8 Model 3 multicomponent displacement data 65 3.9 Model 3 (P, S) decomposed data 66 3.10 Model 3 P P data after multiple removal 67 3.11 M o d e l 3 (P,S) 68 3.12 M o d e l 3 (P, S) data with velocity error 69 3.13 M o d e l 3 (P, S) data after multiple removal with velocity error 70 3.14 M o d e l 3 (P, S) data with missing components 71 3.15 M o d e l 3 (P,S) 72 3.16 L a n d data Model 4 73 3.17 Model 4 displacement data, shot gather 446 74 3.18 M o d e l 4 shot gather 446 after (P,S) 75 3.19 Model 4 shot gather 446 (P, S) data after 2 passes of demultiple 76 3.20 M o d e l 3 displacement data after multiple removal 78 3.21 M o d e l 3 elastic vs. acoustic multiple removal comparison 79 4.1 Source and receiver water bottom decomposition 92 4.2 Ocean bottom data Model 5 96 4.3 Model 5 water column events 97 4.4 M o d e l 5 P wave decomposition 98 4.5 Model 5 S wave decomposition 99 4.6 Model 5 P data after water column multiple removal 100 4.7 Model 5 S data after water column multiple removal 101 4.8 Effect of velocity errors on the removal of water column multiples 102 4.9 Ocean bottom data Model 6 103 4.10 M o d e l 6 P data with (P, S) decomposition and water column demultiple . 104 4.11 Model 6 S data with (P, S) decomposition and water column demultiple . 105 5.1 Relationship between successive interactions for internal multiple suppression 108 5.2 Internal elastic multiple Model 7 118 5.3 Model 7 internal multiple elastic data 120 5.4 Elastic internal multiple suppression for B sl data after multiple removal data after multiple removal with missing components . . . . decomposition { Vlll 121 5.5 Elastic internal multiple suppression for B$ 122 5.6 Elastic vs. acoustic internal demultiple term for Bg'p 123 5.7 Elastic internal multiple suppression for Bp ], and 124 5.8 Elastic internal multiple suppression for B^l 125 1 ix and B ^ Acknowledgement No endeavour is solely an individual effort and I have many people to thank for their contributions, scientific and otherwise. I would like to first thank my supervisor D r . T a d U l r y c h for supporting me during m y years at U B C . This support was financial i n part, but more importantly, he provided the moral and psychological support I needed at various points in my graduate career. Second, I would like to thank my co-supervisor and mentor D r . A r t h u r Weglein. Most of the material in this thesis is the result of my collaboration with him and would not have been possible without his guidance. D r . Weglein also taught me much outside of science and I am thankful for the many lessons I gleaned from h i m . T h i r d , I want to thank my committee members D r . Bruce Buffett and D r . Michael Bostock. I would like to gratefully acknowledge the contributions of the scientists at Schlumberger Cambridge Research whom I had the pleasure to work with. Particularly, I would like to thank Richard Coates for his finite difference modelling code; Luc Ikelle and G r a h a m Roberts for their valuable suggestions; James Martin for his support of the water bottom research; and Phil Christie for enabling the research to progress. Thanks goes to Dennis Corrigan for his invaluable advice on practical issues and also for helping me to get my code off the ground. I would also like to thank Fernanda Araujo for her notes and advice and a special thanks to Mauricio Sacchi who always had time to brainstorm. I am grateful to my many friends and colleagues at U B C who made life livable and stimulating. Particularly thanks to: Xingong L i , Shawn Marshall, Gwenn Flowers, D a v i d Butler, A n t o n Kepic, Julian Douglass for teaching me how to throw a frisbee, and Colin Farquharson for showing me that I couldn't kick a soccer ball. Thanks to m y family who, despite m y frequent doubts, always believed that I could do it. T h e last word goes to my wife Shannon who endured the periods when I was absent and was patient with my frustration when I wasn't. myself to this work. xi Her support enabled me to devote Chapter 1 Introduction Seismic data processing provides many challenges for the exploration geophysicist. Many tried and true methods are based on an intuitive understanding of wave phenomena. Others are steeped in theoretical formalism which makes them unreachable to many practitioners and as a result, are met with suspicion and scepticism. It is my hope that this thesis will not represent another in the line of the latter methods. T h e main focus of this thesis is the removal of elastic multiples from multicomponent reflection seismic data. Multiples can be intuitively understood yet difficult to remove. This makes it a tantalizing problem. T h e first step in solving any problem is to define what the problem actually is. Consider a simple seismic experiment with sources and receivers on the surface and two reflectors at depth, as shown in Figure 1 . 1 . Measurements at the receiver would record the arrival of the two primary reflections whose raypaths are shown in Figure 1 . 1 . Primaries are events which undergo a single reflection in the subsurface and are the bread and butter of the seismic method. From the arrival times of these events, a reflection profile can be obtained across the earth giving a rough sketch of the reflector structure. A clearer image of the subsurface can be subsequently obtained in a process called migration. More quantitative and therefore ambitious methods of analysis can also be performed on the primaries. Examples are amplitude variation with offset ( A V O ) analysis and full or partial inversion. 1 CHAPTER 1. 2 INTRODUCTION Prl Pr2 FSM1 IM1 FSM2 0.8 1000 offset (m) Figure 1.2: Synthetic shot record for the model in Figure 1.1. A n automatic gain control has been applied to the data to show the events at later times. T h e problem is that in this simple example and the real world, multiple reflections are observed as well as primaries. Multiples are events which experience more than one CHAPTER 1. 3 INTRODUCTION change i n vertical propagation direction as they propagate from the source to the receiver. Figure 1.1 shows the raypaths of the first order free surface and internal multiple for the two interface model. A free surface multiple is denned as an event which has at least one reflection at the earth's surface. interacts with the free surface. T h e order of the multiple refers to how many times it A n internal multiple has all of its downward reflections below the free surface. T h e order of an internal multiple refers to the number of downward reflections it experiences. A distinction between the two types of multiples is made because the way i n which they are attacked is different (Weglein et al., 1997). The basic problem with multiples is that they can interfere with, or even masquerade as, primaries in the seismic data. Figure 1.2 shows a synthetic shot gather for the model i n Figure 1.1. Without knowing the model, it is difficult to tell which events are primaries and which are multiples. Not only do humans have problems making this distinction, seismic processing algorithms have an even tougher time of it. Most, if not all, migration schemes assume that the data contain only primaries and no multiples. A s a result, both primaries and multiples are migrated indiscriminately and this leads to confusing or distorted subsurface images. Since a migrated image is usually the end product used for interpretation, a multiple free migrated image is particularly important. In a recent article, O ' B r i e n and Gray (1996) demonstrated that even when the velocities of the earth are known exactly, the end product of depth migration is significantly degraded due to the presence of multiples i n the data. T h e y also showed that even a modest suppression of multiples gives a much better depth migration. More sophisticated analysis methods such as A V O are also adversely affected by multiples (see, for example, Alvarez and Larner (1996) and Foster and Mosher (1992)). It has long been recognized that multiples are undesirable, and numerous methods have been devised to remove or suppress them. Weglein (1995) gives a good overview of the present and future state of multiple technology. Multiples have two major properties which distinguish them from primaries and most conventional methods exploit one or both of these properties. First, multiples tend to have more moveout than primaries (i.e., multiples have more curvature as a function of the source-receiver offset). This is because multiples which arrive at the same time as primaries have more of their travel path i n shallower portions CHAPTER 1. 4 INTRODUCTION of the earth. Since the velocity generally increases with depth, the multiples have a lower moveout velocity than primaries at the same arrival time. M a n y methods seek to exploit this differential moveout property and can generally be classified as filtering operations. T h e second property of multiples that forms the basis of multiple suppression is periodicity. In any seismic experiment, many orders of multiples are generated, whereas primaries only arrive once from each reflector. A s a result, multiples tend to be periodic where primaries are not. This has been exploited by various authors to predict and subtract the multiples from the seismic record. A example of such a method is predictive or gapped deconvolution (see, for example, Y i l m a z (1987)) as it is called i n the exploration industry. Since multiples are periodic only for zero or near offsets, this method becomes less effective with increasing offset. Most, if not all, conventional methods typically assume that either the earth is onedimensional, the velocity is known, or that sufficient differential moveout exists between primaries and multiples. W h e n these assumptions are valid, conventional methods have been and will continue to be effective. W h e n these assumptions are violated, conventional methods become less effective. Examples of these situations are dipping reflectors and decrease in velocity with depth. Unfortunately, these situations are not uncommon. T h e limitation of conventional methods of multiple suppression has motivated the search for methods which are independent of the above assumptions. A method developed by Carvalho, Weglein, and Stolt (1992) uses inverse scattering series to remove free surface multiples from two dimensional marine seismic data without making any assumptions about the nature of the subsurface. This method removes all orders of free surface multiples, including those with converted phases, and leaves the primaries untouched. This latter property is particularly important for post multiple suppression analysis such as A V O . Because multiples are removed independent of the subsurface, this method is effective in situations where conventional methods encounter problems or fail entirely. T h e price to be paid is that it requires data which contain near offsets and a known source wavelet. Techniques have been subsequently developed to address these requirements and the marine procedure is currently demonstrating its full potential on real data. Solutions to the marine free surface multiple problem have also been formulated by CHAPTER 1. INTRODUCTION 5 Verschuur et al. (1988), Fokkema and V a n den Berg (1990) and also by Dragoset and M a c K a y (1993). T h e method of Verschuur et al. is a surface replacement method based on a Huygen's model of seismic wave propagation. Fokkema and V a n den Berg and Dragoset and M a c K a y used Green's theorem and reciprocity to formulate free surface multiple removal. Like the inverse scattering method, these methods do not require an earth model to remove free surface multiples. Similar constraints exist for the wavelet and near offsets. Inverse scattering free surface multiple removal differs from the other subsurface i n dependent methods i n the forward models which are used to describe seismic data. T h e latter methods describe seismic waves as propagating with different velocities through the earth and reflecting and transmitting from interfaces where the earth properties change. T h e former models seismic waves using a point scatterer model with reference medium propagation between scattering points. In Chapter 2, I develop a mapping which connects the point scatterer model and the primaries and multiples of reflection seismic data. While all of the above methods for marine data lead to similar formulations for free surface multiple removal, they differ greatly i n the how they approach internal multiples. Under the Verschuur scheme, internal multiple removal is performed by stripping away the layers of the earth from the top down and removing the multiples with each internal interface. Although this method has its strengths; it requires the velocity profile i n the earth. Using a subset of the full inverse scattering series, Araujo et al. (1994) have developed a method which leaves primaries intact and suppresses all orders of internal multiples without requiring a model of the earth. This is where inverse scattering methods stand alone. T h e reason inverse scattering multiple suppression displays this model independence stems from the fact that the inverse scattering series is the only known multidimensional inversion procedure which does not require any a priori information about the earth. In the process of inverting seismic data into an earth model, this series must perform four basic tasks: 1) remove free surface multiples, 2) remove internal multiples, 3) convert from time to depth, and 4) map reflectivity to earth properties. Since the full inversion requires no model information, it follows that each of these tasks must display this same model independence. It is reasonable to ask why the inverse scattering series is not being used to perform full CHAPTER 1. 6 INTRODUCTION inversion. It turns out that the full series diverges for all but a small range of earth models (Carvalho et al., 1992). Fortunately, convergent subseries have been identified which perform free surface multiple removal and internal multiple suppression for marine data. These are the works of Carvalho et al. (1992) and Araujo et al. (1994). A natural extension of marine free surface multiple removal is the removal of free surface multiples from land data. T h e method of Carvalho et al. (1992) assumes that the sources and receivers are located i n the water and that the multiples are generated by an acoustic free surface above the sources and receivers. W h e n seismic data are recorded on land, the sources and receivers are on or in an elastic medium and free surface multiples are generated by an elastic interface. Verschuur (1991) and Wapenaar (1990) have adapted their versions of marine demultiple to remove elastic free surface multiples. In Chapter 3 of this thesis, I extend the inverse scattering concept to remove free surface elastic multiples from multicomponent land data. Similar to the marine case, this method leaves primaries intact and removes all orders of elastic free surface multiples regardless of the subsurface. There are a couple of reasons for using inverse scattering series for elastic free surface demultiple. One is that inverse scattering multiple removal is fully wave theoretic and contains all the required obliquity factors (the works of Vershuur et al. neglect these factors for both land and marine data). T h e other, and perhaps more important reason, is that an inverse scattering method can be developed for land data which suppresses internal multiples independent of the subsurface. New methods of data acquisition create new demands on seismic processing methods. The use of ocean bottom mounted receivers has been gaining interest as of late because of the potential for increased acquisition flexibility, enhanced seismic wavefield information, and the possibility of acquiring time lapse (4-D) seismic surveys. This latter capability is of particular interest because of the potential for improved hydrocarbon recovery through better reservoir characterization. B y seismically monitoring a known oil or gas field while it is being produced, information can be inferred about how the reservoir is being depleted which can i n turn be used to optimize production strategy. Given the high degree of confidence required in development and production projects, the seismic data must be free of multiples to be effective. CHAPTER 1. 7 INTRODUCTION In Chapter 4 of this thesis, I adapt the elastic free surface method to remove multiples from multicomponent seismic data recorded on the ocean bottom. This represents the first multiple removal algorithm designed specifically for water bottom data which removes all orders of water column multiples independent of the subsurface. After the free surface or ocean column multiples have been removed, the stage is set for attenuating internal multiples from multicomponent land and ocean bottom data. In Chapter 5, I use an elastic background medium to extend the method of Araujo et al. (1994) to attenuate internal multiples from multicomponent data. A n elastic background medium is required since the sources and receivers for multicomponent data are located i n an elastic medium. While the internal multiple algorithm for marine data was originally formulated using an acoustic background medium, it has the serendipitous ability to also suppress internal multiples which contain S phases (Coates and Weglein, 1996). Although this suppression is less dramatic than for non-converted phases, it is nevertheless effective. T h e ability to remove converted internal multiples is important in salt structures where significant mode conversion can occur (Ogilvie and Purnell, 1996). In these situations, the method in Chapter 5 creates the potential for enhancing the suppression of converted internal multiples by using multicomponent ocean bottom data and an elastic background internal demultiple scheme. Furthermore, the work i n this thesis provides an avenue to potentially enhance the suppression of converted internal multiples by using standard marine data and an acoustic/elastic background medium. Although it is not necessary that land and ocean bottom data be acquired with more than one component, there are advantages to doing so. two modes of seismic wave propagation, P-waves In an elastic earth, there are and S-waves. These waves respond to the properties of the subsurface in different ways. B y exciting and recording both of these wave types independently, much more subsurface information can, in principle, be obtained from the various P and S primaries. While this is attractive, these P and S primaries are also accompanied by a complex set of P and S free surface and internal multiples. T h e work i n this thesis is designed to attenuate these elastic multiples so that the most reliable information can be obtained from the P and S primaries. Chapter 2 Forward scattering series 2.1 Introduction Forward scattering is concerned with the process of forward modelling. In the context of the seismic problem, it describes how seismic waves propagate through a given earth model. T h e way in which seismic events are modelled in scattering theory is fundamentally different from conventional wave theory where seismic waves propagate through the earth with different wave velocities and are reflected and transmitted at material boundaries. T h e unique formulation of the forward scattering problem leads to some interesting properties in the inverse scattering problem. Specifically, it provides a framework in which internal multiple suppression is possible without requiring a model of the events which generate the multiples. This is not obvious using more conventional models of wave propagation and inversion. In scattering theory, the earth is viewed as a background or reference m e d i u m and a perturbation about that reference medium. In order to keep the perturbation (the unknown) on one side of the source and receivers, the reference medium is chosen so that it agrees with the actual medium at and above the sources and receivers. For example, in the marine case, the background medium is chosen as a free surface bounding water which extends down to infinity. T h e rock layers (including the water bottom) are then simply a perturbation about an infinitely deep ocean. Using this formulation, the wavefield in the 8 CHAPTER 2. FORWARD SCATTERING 9 SERIES water plus rock is the wavefield in water plus a perturbation wavefield which is nonlinearly related to the rock properties. This nonlinear relationship is called the B o r n series and can be interpreted as an infinite series of propagations in the reference medium separated by different orders of scattering interactions with a point scatterer earth model. For the marine case, it is as if the solid earth were made up of point inclusions surrounded by water. In the seismic experiment, waves emanate from a source, propagate through the water, are multiply scattered from the solid bits, then propagate back through the water to a receiver. T h e different terms i n the perturbation series correspond to the number of scattering interactions a wave experiences. Figure 2.1 shows a simple reflection and transmission problem from the point of view of a boundary value problem and scattering theory. T h e goal of inverse scattering is to obtain information about the earth model from the measured seismic data. Like forward scattering, inverse scattering is a nonlinear process which takes the form of an infinite series. Different techniques have been developed for one- and multidimensional inverse scattering which have their relative merits when applied to acoustic waves (Weglein, 1985). These nonlinear methods are generally sensitive to data errors (Snieder, 1990) and are possibly subject to more fundamental convergence problems (Prosser, 1980). Recently the application of inverse scattering series has led to powerful new methods of multiple suppression. Using a subseries of the full inverse series, Carvalho, Weglein and Stolt (1992) developed a method to remove free surface multiples from marine seismic data. In a related work, Araujo et al. (1994) identified an inverse scattering subseries which suppresses internal multiples. These methods suppress multiples independent of the earth model below the sources and receivers. This is a nice property since the efficacy of many conventional multiple suppression methods depends on certain assumptions the earth model. about For example, they require that multiples and primaries have different velocities or that the earth is one dimensional. Inverse scattering demultiple methods as a whole depend on an understanding of how multiples are modelled in forward scattering. T h e purpose of this chapter is to clarify this understanding for internal multiples by analyzing simple forward scattering problems where the solutions are well known. This will establish a mapping between a point scatterer CHAPTER 2. FORWARD SCATTERING SERIES 10 Figure 2.1: a) Reflection and transmission from an interface as a boundary value problem, and b) the same problem from the scattering theory point of view. Some of the contributions to the first, second, and third order scattering terms are shown. CHAPTER 2. FORWARD SCATTERING 11 SERIES model and the primaries and internal multiples of reflection seismic data. In Chapter 3 , I will discuss this relationship for primaries and free surface elastic multiples. T h e second section of this chapter gives a brief introduction to the forward scattering B o r n series. In Sections 2.3.1 and 2.3.2, I develop analytic B o r n series solutions for plane waves incident on simple one dimensional earth models. 2.2 The Born series T h e goal of the forward acoustic problem is to find the wavefield which emanates from a localized source and propagates through a prescribed model. Often the model has a velocity distribution which is constant over intervals and discontinuous at the interval boundaries. T h e earth, for example, is made up of rock layers each having a roughly constant velocity. For such a layered model, the forward problem is traditionally solved as a boundary value problem where constant velocity wave solutions are found which are valid for each layer. To construct the solution which is valid for the entire model, these separate solutions are matched at the layer boundaries according to a set of boundary conditions. For simple models, the total wavefield contains various types of reflected and transmitted waves. I will show how scattering theory can be used to generate the same results. Following closely the work of Weglein (1985), I begin with the two dimensional constant density acoustic wave equation where the scalar P(r|r ;f) represents the pressure at point r = (%, ) z s a n d time t due to a source at point r and time t — 0. T h e waves are emitted by a line source which effectively 3 makes this a two dimensional problem. Throughout this thesis, I will be freely forward and inverse Fourier transforming over time and space variables. T o avoid introducing different symbols to denote Fourier transformed quantities, I will let the form of the variable dependencies denote that a transform has been performed. Spatial variables will be denoted by x, y, z or the vector r with appropriate subscripts g or s to denote geophone and source locations respectively. Spatial transform variables (wavenumbers) will be represented by 2. FORWARD CHAPTER SCATTERING 12 SERIES k with subscripts indicating the conjugate spatial variable (i.e. be denoted by t and w is the angular temporal frequency. k s T i m e will —> x ). 3 Hence, Fourier transforming P(x, z\x , z ; t) over time and horizontal shot and receiver locations gives P(k, z\k , z ; UJ). 3 a s 3 In (2.1), the spatially varying velocity c(r) can be characterized by a constant reference velocity CQ and a perturbation a(r) so that 1 c- (r) 2 =i[l-a(r)] (2.2) or «(r) = 1 - (2.3) A common and convenient choice for reflection data is to choose the reference medium as the actual medium at and above the source and receiver. This places the perturbation on one side of the source and receiver and leads to significant mathematical simplifications. Fourier transforming (2.1) with respect to time and substituting (2.2) gives / ,2 V + u> \ 2 -) UJ 2 P(T\T ; UJ) = S(T - rs) + ^ a ( r ) P(r\r3; UJ) (2.4) S T T o find a solution to (2.4), consider the causal free space Green's function Po(r|r ; UJ) which s satisfies the equation (V + A; )Po(r|r,;a;) = ^ ( r - r , ) . 2 where k 0 = UJ/C . 0 0 (2.5) 2 Although k is not actually a spatial transform variable, UJ/CQ has the 0 proper dimensions of a wavenumber. In fact, k can be thought of as the spatial transform 0 variable i n the direction of wave propagation in the reference medium. Using Po(r|r ; UJ) as s a reference wavefield, an integral equation corresponding to (2.4) and its physical boundary conditions is (Weglein, 1985) oo / P (r|r'; UJ) k a(r') P(r'\r ; UJ) dr'. 0 2 0 s (2.6) oo This is called the Lippmann-Schwinger equation and plays a pivotal role i n scattering theory. It represents the wavefield in an inhomogeneous medium as the sum of the wavefield in a reference medium and an integral that represents the scattered field due to a 2. FORWARD CHAPTER perturbation. SCATTERING 13 SERIES Iterating the Lippmann-Schwinger equation back into itself generates the B o r n series CO P o ( r | r » k a(r') P (r'|r ; LJ) dr' 2 0 0 / + r P (r\r';u>)k a(r') 2 0 0 =P + P i + P 0 2 s OO ( f°° P (r'\r"; UJ) k a(r") P (r»\z„ z ;«,) dr") dr'+ .. 0 2 0 0 a +.... (2.7) Throughout this chapter, I will refer to the terms i n this series either according to their order in a, or their place in the sequence. For example, the zeroth order term, P , is the 0 first term in the series. Similarly, the first order term P i , is the second term i n the series and so on for the other terms. To facilitate a physical interpretation of the Born series, consider the scattering perturbation to be composed of a region of closely spaced point scatterers separated by the reference medium. T h e first term in the Born series is simply the reference medium Green's function. It represents a wave which propagates in the reference medium from the source at r , directly to the measurement point at r. s T h e second term in the Born series contains k^a^r') sandwiched between two Green's functions. T h e Green's function on the right represents a wave which propagates from the source at r, to a point scatterer at r'. T h e wave undergoes an interaction with the scatterer at this point and then propagates back to r via the Green's function on the left. T h e first order term is then the integral over all possible single scattering interactions. Similarly, the third term represents a sum over waves which propagate in the reference medium and undergo two scattering interactions. Following this interpretation of the scattering process, each term in the Born series involves a series of propagations and interactions with points within the scattering region. This is the mathematical foundation for the scattering picture shown in Figure 2.1. It is interesting that all of the Green's functions in the series represent waves that propagate from one point to the next in the reference medium. This is true even within the scattering region which would seem to contradict the fact that waves travel with different CHAPTER 2. FORWARD SCATTERING velocities in different materials. SERIES 14 In fact, all of the scattered solutions together superpose to create a wave which travels at a velocity different from that of the reference medium. Truncating the B o r n series after the second term leads to the B o r n approximation. In this approximation, the data (the measured scattered field) is assumed to be linear in the model (the scatterer). W h e n the model is small, the higher order terms in the B o r n series become less important and the B o r n approximation is valid. This is often referred to as the weak scattering approximation. For the point scatterer model, this corresponds to having the points differ from the reference medium by only a small amount. W h e n the scattering is not weak, the higher order terms play an important role in describing the total wavefield. In the following sections, I will demonstrate how each of the terms in the B o r n series contributes to the modelling of data. 2.3 Two dimensional scattering from one dimensional models In this section, I develop solutions for plane wave scattering from one dimensional scattering models. I consider a two dimensional problem where the earth model varies in the z direction only (see, for example, Figure 2.2). Some take the view that one dimensional problems have little practical value and while I agree that the earth is not truly one dimensional, I do believe there is value in studying what A r t h u r Weglein calls 'baby problems'. These simple problems allow you to build an intuitive understanding of the physics which is, in my opinion, often the key to creating something genuinely innovative. To give the following examples the desired simplicity, I neglect the existence of a free surface. T h e reference medium can then be chosen as a homogeneous whole space composed of the material in the region outside of the scatterers. T h e source will be located in the reference medium and the receivers will be located both outside and inside the scattering region. This describes reflection and transmission data respectively. For two dimensional wave propagation, consider a line source of wave energy. The reference wavefield is then the two dimensional Green's function that satisfies (2.5). This 2. FORWARD CHAPTER SCATTERING 15 SERIES is written using the Weyl integral 1 where u g g s s ik = — / G (x ,z \x ,z ;uj) 0 f°° Q 'i 9~ ) Q^'\ 9~ \ ^ J-oo — X Xs Z Zs — (2.8) dk twos s = (UJ jc\ — k ) ^ , and c is the wave velocity in the reference medium. Through2 03 2 1 2 0 out this thesis, v will be used to represent the vertical wavenumber for P wave propagation. T h e subscripts for v will indicate the velocity c and wavenumber k which are related to v. Rewriting (2.8) as = — / — <f>o(x ,z \k ,z ;uj)dk ^ J-oo two* Go(x ,Zg\x ,z ]uj) g 3 s g g s s (2.9) s where Mx z \k. z ;u>) gi 8 t = e ^ + ^ — D , t (2.10) it is evident that the 2-D Green's function is a superposition of weighted plane wave solutions. This motivates the use of a plane wave source as the incident wave. Solutions for a line source can later be constructed using the plane wave solutions. More formally, the 2-D Born series is P(x ,Zg\x ,z ;uj) g 3 s = Go(x ,z \x ,z ;uj) g CO / •OO CO / \ J— where V(x,z) , G (x , z \x', z'\UJ) V(X', z') G (X', 0 g g 0 3 a G (x ,z \x',z';u;)V(x',z'). 0 g g J — OO / G (x', z'\x", z"\ UJ) V(X", z") G (X", J — oo = k^a(x,z). z'\x , Z ]OJ) dx'dz' f OO Q CO 3 3 OO / • CO g CO 0 z"\x , Z ; UJ) dx"dz" + ... s 3 (2.11) For each term in the series, Fourier transform with respect to 2. FORWARD CHAPTER x 3 SCATTERING 16 SERIES and multiply by 2IVQ . This gives S P(x ,z \k ,z ;uj) g g a 3 <f> (xg,z \k ,z ]u) = 0 OO / g f* OO / J—oo poo poo + I G (x , z \x', z'\ UJ) V(x', Z') 4> (X', z'\k , z ; UJ) dx'dz' 0 oo s t g 0 g 3 3 G (x ,Zg\x ,z';uj)V(x',z')- I 0 l g G {x', z'\x\ z"; UJ) V(X", z") <p {x', z'\k , z ; UJ) dx"dz" + ... 0 = Po + 0 3 3 Pi + P2 + • • • • (2.12) This equation shows that each term in the series is evaluated for an incident plane wave <f>o{x ,z \k„z.;u>) = g (2.13) g where k is the horizontal wavenumber and uo is the vertical wavenumber. Different values 3 3 of k then correspond to incident wavefields with different incident angles. 3 T o reconstruct the series for a line source, simply multiply each term by Fourier transform with respect to x . 3 2iv 0s and inverse For the reference wavefield, care must be taken to preserve the absolute value signs in (2.11) so that the reference wavefield is symmetric about the source. This is not a concern for the scattered field i n the following examples since the source is always located above the scatterers. 2.3.1 Scattering by a half space T h e simplest case to consider is scattering from the half space as shown in Figure 2.2. I will refer to this as Model 1. T h e interface between the two media is at z\ and, for simplicity, the source is located at z = 0. T h e scattering perturbation for this model is s V(x\ z') = V(z') = k H(z' 2 oai where aj = 1 — c\jc\, - ) Zl (2.14) k = UJ/CQ, CI is the wave velocity in the half space, Co is the reference 0 velocity, and H(z' — z\) is the Heaviside function. Clearly, the perturbation is zero where the actual medium is the reference medium and nonzero elsewhere. 2. FORWARD CHAPTER SCATTERING 17 SERIES incident wave x ^ z measurement point (x ,z ) v g g Zl Cl Figure 2.2: Model 1: acoustic plane waves incident upon planar interface. The first term in the B o r n series is just the incident wavefield z \k; UJ) = e ^ s + ^ s ) . P (x , 0 g For brevity, I have dropped the z = 0 and will use just k for k . Keep in m i n d the difference s between k and k (2.15) g s and that v% = k^ — k . 2 0 F r o m (2.12), the second term is ^ P (Xg,Zg\k;U ) 1 J = /»oo />oo /.oo — / / 2TT e x ) f ° ° i /"oo = — / / 2TT J_ J_ Using the Fourier transform 00 dkg 0g kx + r (2.16) dx'dz' 7,2 POO / J x g 2tu e ( ' '** > Zl 9 — oa 2 1 g J_ k a H(z' 0 ik (xg-x') ivo \z -z'\ e / e i(fc - )x' ik x i^ \z -z'\ iv z> fe9 e g 9e g g e _o_j_ 2iu 0 Z1 dz'dx'dkg. 0g oo Kk-k y > e / g dx = - k), 5(k 27r g (2.17) •oo (2.16) becomes / £(& - jfey^e**.!*.-*'! oo Jz! s B y the sifting property of 5(k g — k), k g —> A; and fog —• e ge 0 g e <fc'<flb 2 fl 2zz/ (2.18) 0g giving ikx iv \z>-z \ i u e -° ' 0 z ' d z > ^ 18 CHAPTER 2. FORWARD SCATTERING SERIES W h e n the receiver is on the source side of the interface, z < z\, and g 1 (2-20) J - " e 'lfe£ aj e *« / e** ' dz'. = _iM, ,fce 22 2^o J Z 1 This integral can be solved by considering a small amount of dissipation i n the wave propagation. It is physically reasonable to expect that all real wave systems will display some amount of dissipation such that waves infinitely far from their source will have negligible amplitude. This is satisfied by adding a small imaginary part to the wave slowness such that 1/cn = l / c ( l + ie) where e > 0 for w > 0 ( A k i and Richards, 1980). 0 W i t h this modification, (2.20) can be solved directly giving < \k-uj) = ~k Piix^zg **9 iM^-z ). 2 Zl 4z/ oai e e (2.21) e 0 This can also be obtained by a contour integration i n the complex z' plane. A d d i n g Po and P i gives the Born approximation PBorniXg, Z < \k] u) = j ^ ' + w ) + - L * £ G Zl 0 l e ^x >o(2 -z ) _ In this expression, the first term is a wave which propagates J f Z a ( .22) 3 2 outward from the source directly to x , z . T h e second term describes a wave which is reflected from the boundary g g and propagates with the reference medium velocity. T h e preservation of the horizontal wavenumber before and after reflection means that the angle of reflection is equal to the angle of incidence. To see this, consider the dispersion relation (UJ/CQ) = u^ + k where OJ/CQ 2 2 can be thought of as the wavenumber in the direction of propagation and u and k are its 0 vertical and horizontal components. These components can be expressed as k = UJ/C sin(# ) 0 0 and vo = to/CQCOS(8Q) where 9 is the angle between the propagation direction and the normal to the interface. Since the horizontal wavenumber and velocities are the same before and after reflection, sin(#; ) = sin(0 /z). nc re Consequently the incident and reflected angles are the same. As in traditional wave analysis, the ratio of the amplitudes of the incident and reflected waves is the reflection coefficient. For the Born approximation this is p H B o r n _ [Pil _ k2ai _ _ ^ | P q | 4 c o s 1 ^ ) ( 2 C L - CQ)(CI + ^ Cp) . [2.26) CHAPTER 2. FORWARD SCATTERING SERIES 19 0.5 r co/ci Figure 2.3: Comparison of the true reflection coefficient to the B o r n approximate reflection coefficient for normal incidence. T h e velocities of the first and second media are c and C i 0 respectively. For normal incidence, set 9 = 0 and = ifa-»)(•»+ »). 4 ( 2 . 2 4 ) c{ This result was also obtained by Weglein and Gray (1983) i n a study of the sensitivity of B o r n inversion to the choice of reference velocity. In the corresponding boundary value problem for Model 1, separate solutions on either side of the boundary are matched by equating the pressure and the normal component of velocity across the boundary. This gives the angle dependent reflection coefficient R = — fi/vo + Ui where v\ = (u> jc\ — & ) / 2 2 1 2 (DeSanto, 1992). For normal incidence, this is (ci — c ) / ( c i + c ) which is clearly not in agreement with (2.24). Figure 2.3 shows, however, 0 0 that the normal incident B o r n approximate reflection coefficient does approach the proper value for small values of a i.e., when CQ/C\ ~ 1. This is the weak scattering 0 convergence criterion for the forward Born series (see, for example, Prosser (1975)). For transmission data, z g > z\ and the integral in (2.19) must be split into two portions: CHAPTER 2. FORWARD SCATTERING SERIES z < z' and z g > z'. This gives g fz J.2 g P (x ,z >z \k-uj) 1 g g 20 = 1 e^e ^'^dz' J J _ = + ™o Zl ^ e - o , e 9 f c 2 a jL2 rOO i i ( e J Zg 2i ( Vo 1 + Zl - ^ e ^ ' e ~ ^ ^± > g dz , 2 )). ^ 25 Zg Combining this with the first term gives the Born approximate transmitted field > \k ) = P orn(x ,Z B g Zl g Ik? -e^ ^ ^a (l + 2iu (z -z )) - 1 + ]U 1 Avl 0 1 (2.26) g This wave propagates in the proper sense (increasing x and z directions); however it propagates with the reference medium velocity, has the wrong amplitude, and i n addition, does not have the proper refraction angle. In the following, I show that by adding more and more terms to the forward B o r n series, the solutions for the reflected and transmitted waves converge to the familiar expressions. T h e third term in the series is oo poo / G (x , z \x', Z';UJ) kl a(x', z')Px(x\ z' > 0 g zi|fc;u>) g / •CO w —CO Y dx'dz />co yco yco ikg(x -x') iu \z -z''\ e = —/ 2TT 7_ / e 0g g / J_ 00 g k% a J X 1 Pi(a; , / z' > Zi\k; UJ) dx'dz'dk g 2iu Z1 0g (2.27) Expressing P (x',z' > z^fyuj) as e 'P^x', z' > z \k\uj) lkx x y>oo -> P (x ,z \k;uj) = ± 2 g pOO g ^ £oo poo pOO / x ikg(x -x') iv \z -z'\ ikg(Xg-X') lV g\Zg-Z\ ^ j-A kla^ *P (x',z'>z \k-uj)\dx'dz'dK e pOO / g e 0g g 0 k J-oo J-oo J Og / ****S(k - ky^-t'^^-Piix^z' e g •co J / ikxg iv \z -z'\^pr ^ > e e 0 Z1 I Zl g 3 ™o x > z z^.^dz' 2 J = > \k;uj)dz'dk 0g 2lU Z1 = x T/i. ^ pOO / x Ztv Zl CO , this gives 2 «o^i ite, 2wo e /»CO / ^\*,-''\p^ ' ' x iZ > Zl\k;uj)dz' J Z1 which can be split into expressions for reflection or transmission data. (2.28) 2. FORWARD CHAPTER SCATTERING 21 SERIES B y calculating the higher order terms, I found that each term leads to an expression of the form Pn = elkx'Pn where the horizontal wave number is preserved. From this a recursive relation can be written to generate all the terms in the series: 7 2 n /»0O kn < z \k;w) = ^ i e * " P (x ,z g x g ^ 0 lv / e^''^P ^(x\ (2.29) z' > \k;uj)dz> n Zl Jz x for reflection and P (x ,z n g > z.lk-u) = ^±e *° ik g 2*IVQ k0 n,l ik 2 K a e f" ^ < « • - • ' > P _ ( x ^ * > \k;u;)dz' C n , 1 Zl J (2.30) Zl fiOO / e J Zn '- i M o ( z P „ - i > z \k;uj)dz' Z f l ) t for transmission. Notice that each term in the Born series depends on the transmitted solution for the previous term. It is also interesting that the calculation for two dimensional wave scattering from a half space can be reduced to a series of integrals in depth thus mimicking the one dimensional problem (Matson, 1996). Calculating and summing terms gives the reflected wavefield P(x , z > \k; g g Zl UJ) = P i + P 2 + P s + • • • ikxg iu (2z -z ) e e 0 1 g lk ai 2 0 4 vl t la\k* + x 8 | v% 5 a\kl 64 | + vi 7 qjfcg 128 vl +• (2.31) which is a wave reflected from the boundary at z\ and, as for the first term , has a reflection angle equal to the incident angle. T h e reflection coefficient can be obtained using U vl > 2v 2 % vl> { l§ vl> 128 vl ' K 1 25Q vl ^---' { ) (2.32) which leads to the expected result 1 - 4 / 1 - aikl Vn 1 = (2.33) CHAPTER 2. FORWARD SCATTERING 22 SERIES Figure 2.4: For a plane wave incident on a half space with incident angle 0; nc and unit amplitude, the sum of the scattering terms results in a reflected plane wave with reflection angle #; nc and angle dependent reflection coefficient v = ^u jc\ - k where, as before, 2 0 2 R(9i ). nc v = v /l - UJ a jc vl = and x 2 oy x 2 0 AA> / I 2 C k. 2 Figure 2.4 shows the scattered wavefield that results from summing the B o r n series. For the transmitted wave, the series yields P(x , z > \k;uj) = 4> + Pi + P + Pz + g g 1 + iv (z 0 g Zl 1 4 0 ak t 2 2 , 1 a\kA0 , 5 offcg + + 64 v v$> 8 v$ e - Zi) 1 a-ikl 1 n a\k\ 5 1 a k? 3 2 o 8 vl 7 a\k* + ... a\k% 7 I n 32 nk 3 6 a. vl 8 n 9 1 afog 64 /y 32 0 " o ( * s - * i ) , + 128 (X^/uQ 128 zv 8 0 (2.34) + + ... 23 CHAPTER 2. FORWARD SCATTERING SERIES which can be expressed as P(x ,z g g > z \k-u) = e »e ^(l ikx x 1 + - V )(z 0 - g + R) iv Z) X - ~{V\ ~ V ) (Zg - 2 0 ^if (2.35) + = (1 + R)e ikXg e' VlZfl e 2 '( i- o)( fl- ij i/ i/ z z T h e transmission coefficient, T, is given by 2i/i 1 z'l (2.36) v-i + as expected. In addition, this wave propagates with the proper medium velocity (c\) and has the correct refraction angle. T h e latter follows from writing uj/c f re r sin(# / ). F r o m this the familiar Snell's law s i n ( ^ ) / c j re r nc k nc k as = = u>/ci nc sm(O f )/c f re r sin(0j ) — re r nc can be derived. Figure 2.5 shows the resultant transmitted wavefield. It is amazing that not only do the various scattered waves propagating in the reference medium superpose to create a wave which travels at the proper medium velocity, but that these waves also superpose to cancel the reference wavefield . It is natural to ask under what conditions these series converge. It turns out that this is entirely dependent on the convergent properties of the series in (2.32). T h e ratio test can be used to show that this is satisfied for |aifc§/i/Q| < 1. Using u = k cos(8) where 8 is 0 0 the angle between the direction of propagation and the normal to the interface (see page 18), convergence is ensured for sin(0) < c / c i < ( 1 + cos (^)) / . Figure 2.6 shows the region of 0 2 1 2 convergence for a given incident angle. T h e curve 1 — (crj/ci ) 2 as a function of (cj/co) must remain between the lines dotted lines given by ± cos (#). W i t h increasing angle, this region 2 narrows. Hence the series will diverge for large angles and high velocity contrasts. This is more restrictive than the normal incidence case (8 = 0) where the reference velocity only need be less than \ / 2 c i . B y incorporating nonzero offset, the reference velocity has a lower bound given by sm(#)ci. CHAPTER 2. FORWARD SCATTERING SERIES 24 V receiver z, Figure 2.5: For a plane wave incident on a half space with incident angle 0; nc and unit amp- litude, the sum of the scattering terms results in a transmitted plane wave with angle and angle dependent transmission coefficient T(6i ). nc 6 f re r T h e transmitted wave propagates with the actual medium velocity C\ and with transmission angle given by Snell's law. This simple example illustrates how the higher order terms in the B o r n series alter the amplitude and adjust the propagation velocity of the scattered wavefield. T h e next example demonstrates that the higher order terms also describe interbed multiples. 2.3.2 Scattering by a layer over a half space In this section, I consider scattering from Model 2 as shown in Figure 2.7. M o d e l 2 consists of a single layer with velocity c i on a half-space with velocity c . A s before, the reference 2 velocity is CQ. This model will generate multiple reflections within the layer in addition to the simple reflected and transmitted waves. For this model, the perturbation part of the velocity is written as a(z) = a\H(z — zi) + (a 2 a 2 — a\) H(z — z ) where a\ = 1 — c\jc\, 2 = 1 — C Q / C | , and Z\ and z are the locations of the first and second interfaces respectively. 2 As before, the source is located at z = 0. s In practice, measurements are often made for reflection data which includes both the direct and backscattered wavefields. In the following example, both the transmitted and 25 CHAPTER 2. FORWARD SCATTERING SERIES Figure 2.6: T h e Born series for scattering from a half space converges when the curve 1 _ (c /ci) 0 2 is i n the region between ± c o s ( 0 ) where 8 is the plane wave incidence angle. 2 For wide incident angle and large velocity contrasts, the series diverges. reflected wavefields are calculated for each term; full series solutions will be presented for just the reflection data. Using (2.12) and (2.17), the first order term in the Born series for the region z g ikx a iu (2z -zg) l ^ n 1 k a , iu (2z -z ) kl(a - ai) 0 2 g < z is 1 (2.37) 2 which describes waves reflected from the first and second boundaries with reflection angle equal to the incident angle. For the region between the interfaces, the first order term gives a wave transmitted 26 CHAPTER 2. FORWARD SCATTERING SERIES incident wave \ v measurement point (xg,z ) . I g Cl 22 c 2 Figure 2.7: Model 2: acoustic plane waves incident upon a layer over a half space. Measurements of incident and scattered waves are made at x ,z . g g through the first interface and a wave reflected from the second interface P (x , zi<z < 1 g g P M^***» > + z \k; w) = e *> ik 2 e J Z1 f°° L% I dz o J 2lu Zg H M°I *>(2*'-*,) E D Z > 2iu 0 n 4zv 4i/ 0 J Q (2.38) For z > z , the solutions are transmitted waves through both interfaces. g 2 P (x ,z <z \k;u>) 1 g 2 g = e > ikx H dz' + f ' M%*>*« o V + f~ ^jM*' -*.) 1 > dz (2.39) To first order, the B o r n approximate solutions for this model predict the expected types of waves. Beyond that, however, they display errors in amplitude, velocity, and propagation angle similar to the previous example. Multiple reflections are also conspicuously absent. CHAPTER 2. FORWARD SCATTERING 27 SERIES Starting with equations (2.37)-(2.39), the following recursive relations can be derived to calculate successively higher order terms in the B o r n series: < zi\k;u>) P (xg,z n g = e lkx M^ie-o(,'-, )p _ r 9 + r r2 iMz 2 7 Z3) f tf^jMz'-^p^zi 22 + Jz 2 i v a n f J < z \k;uj) = e 2 2 g + Z2 k r Jz 2 POO Zl < z' < z \k)u>)dz' (2.41) 2 ° r !&ljM*'-*.)p _ ( g \ -u)dz' < < z'\k-w)dz' ^e '- Pn-i(z2 ™o Jz n > i(zi <z (2.40) coo P (x ,z n 1 Z2 < z'\k;u)dz' ^e^'-^Pn-iizi o <*< 2lu z \k;uj)dz' 2 ^ e * * < * ' - " ' > P „ - i ( * 2 < z'\k;Lo)dz' 2zz/ (2.42) 0 7,2 ^ - ( ^ ^ ( z ^ z ' l ^ u ; ) ^ ' + / B y computing a number of terms i n this series, I found that they could be grouped according to what type of event they contribute to. For reflection data, this is written as P(Xg,Zg + < Zl|fc;U>) P (x , ndtl g = z < g P (Xg,Z 0 g < z \k\Lo) + x Zl\k\ Uj) P^ixg, P (Xg,Zg PTl + Z g < Z^Uj) < z^k]u) + + P^ (xg , Zg < z\ |fcJ LO) higher order multiples . (2.43) where the superscripts p r l , pr2, m l t l , and mlt2 represent the first and second primary reflections, and the first and second order multiple reflections respectively. Events and their series constituents have the same travel path except the latter propagate with only the reference velocity. This phase signature makes it possible to identify which terms contribute to which events. CHAPTER 2. FORWARD SCATTERING 28 SERIES Since each term describes a different order of scattering interactions, there is a relationship between the terms in the B o r n series and interlayer multiple reflections. For reflection data, multiples do not appear until a certain odd number of scattering interactions are summed. For example, a first order multiple involves three reflections: two from the second interface and one from the first. Since this type of path contains three changes i n propagation direction, the first order contribution to the internal multiple is contained in the third order term of the B o r n series. Similarly, a second order multiple requires terms with at least five scattering interactions before it begins to appear. It should be empha- sized that it is not only the terms in the series involving odd numbers of interactions which create multiples. For a given multiple, its first term and all of the following terms contain information about that multiple. For example, the first order multiple requires the third, fourth, fifth, and all higher order terms for its construction. So not only do the higher order terms contribute to the primary reflections and transmissions as i n the previous example, they progressively describe higher and higher order multiple reflections. Figure 2.8 charts the different scattering terms and the events which they contribute to. T h e first primary reflection is + 7 /Afc.V ,+ 21 (h\a^ 1 + ) 2k(T) iia^ " ' ' ( 2 4 4 ) Similar to the previous example, using (2.32) leads to the solution P (x ,z pil g < \k]uj) = g Zl ik* iM2*i-* ) Rie ge (2.45) 9 where R i = [2 - 2(1 - fc aiK ) 2 2 ~ fcpaiK ] 1/2 2 0 "o ~ " i = ( 2 4 6 ) is the reflection coefficient corresponding to the first interface as shown i n Figure 2.9. T h e second primary reflection is represented by a more complicated expression involving terms i n ko, VQ, Z , Z\, a-i and a . 2 2 W h e n certain series expansions are used, these terms reduce to P (x , pi2 g z < \k-uj) = T[R e o ^- ^e ^- ^ g Zl Tl 2 ikx il/ e z iv z (2.47) CHAPTER 2. FORWARD SCATTERING 29 SERIES number of scattering interactions event direct wave 1st primary 1 2 3 4 5 6 2nd primary 1 2 3 4 5 6 3 4 5 6 5 6 1st order multiple 2nd order multiple Figure 2.8: Events in the reflected wavefield are shown on the vertical axis vs. the different B o r n scattering terms which contribute to the modelling of that event. where, as before, u\ is the vertical wave number in the second medium (see A p p e n d i x A ) . T h e transmission coefficients T\ and T[ are given by 1 + Ri and 1 — Ri respectively and correspond to waves transmitted through the first interface in the + and — directions. T h e reflection coefficient R 2 2 - fcga /i/ a K i a 0 relates to the second interface and can be written as - kjai/v*) '* 1 - 2(1 - ^ ^ ( 1 - kla lvl x n—r~2——T~2 = — —;— - u 2 • T h e second primary is shown in Figure 2.10. T h e corresponding portion of the B o r n series for the second primary can be obtained from (2.47) by first expanding exp[ifi2(;z2 — i)\ z as a power series in iv-\2(z 2 — zi), then substitute the expansion for (1 — A ^ a x / V o ) / into v\, 2"i, and T[. Next, the expansions for 1 (1 - klax/viyl , 2 (1 - A; a /^o) 2 2 1/2 a n d 2 /( l 2/^l 1 k a - tfax/v ) 2 are substituted into (2.48) which is in turn substituted into (2.47). T h e convergence of the second primary depends entirely on the convergence of the series expansions for (1 — k\a\lv\y^ and R . 2 T h e convergence of (1 — k\axfv^) 1//2 was discussed 30 CHAPTER 2. FORWARD SCATTERING SERIES Cl ; c z 2 2 Figure 2.9: Born series first primary for a plane incident wave. in the previous example and is ensured for sin(#) < c / c i < (1 + c o s ^ ) ) / . For P t , there 2 0 1 2 2 are two choices for expanding the denominator i n (2.48), namely: K il l-^l \l l a v a - a /a ) tf 2/vo(l v a 1 (2.49) 2 and , , 2 / K^/^oi - \ • (2-50) «2/ai) 1 These converge for |a | > |<xij and \a \ < |oi| respectively. Interestingly, (2.47) yields the 2 2 same solution using either expansion. Consequently, the convergence of R 2 on the convergence of ( l — ^ o i / ^ o ) a 1 / / 2 a n d ( l ^ o 2 / ' o ) ' - T h e two primary reflection series — a Z/ are then convergent for sin(0) < c / c i < ( l + c o s ( 0 ) ) 0 depends entirely 2 1/ 1/2 2 and sin(6) < c /c < 0 2 (l+cos (0)) 2 1/2 where B is the incident angle. Like the second primary, the first order multiple P " ' 1 terms i n k, x , x\, a\ and a 2 2 1 1 is also a messy looking sum of (see Appendix A ) . B y employing the same series expansions CHAPTER 2. FORWARD SCATTERING 31 SERIES C2 Figure 2.10: B o r n series second primary (plane incident wave). This wave has the correct amplitude and propagates with the proper angle and velocity through the layer. as for the reflected waves, this portion of the series becomes P^ {x ,z a g < Zi\k;u) g = - Ri(l - R\)R\j^'e^W *-^- j(\ 2 _ [4(z - z Q ( 2 + 4(z - z 2 - 7 l n! = — Ri(l — i 2 ) J ? e = - i?i(l - 2 ifcXs 7 l e i l / 0 ( 2 ( 2 z 2 - z i) _ z s) e i s / l)iu - [^o4(a2-»i)(7i-l)] 2! 2 Z i ) ( 0 w > \ 7 o4(z -2i)(7i-i) 2 /^J^e^e^C^'i-^Je^ ^-*!) 4 (2.51) where 7 i = (1 — k^ai/vl) / . From Figure 2.11, it is clear that this expression properly 1 2 describes the multiple with respect to both the amplitude and the phase. For the second multiple, the terms reduce to P (x; k) = R\(l - RDRle^e^^i-^)^!^-^) Tnlt2 (2.52) A s before, the corresponding portion of the B o r n series is obtained by expanding out exp[i/Vi6(z — zi)], Ri, and R . 2 2 A l l of the multiples, including the higher orders not shown here, contain terms i n Ri, R , and exp(ifo(l — k^ai/vl) ! ). 2 1 2 Since the primaries are also dependent on the CHAPTER 2. FORWARD SCATTERING 32 SERIES Z2 c Figure 2.11: 2 B o r n series first internal multiple (plane incident wave). Similar to the second primary, this wave has the correct amplitude and propagates with the proper angle and velocity through the layer. convergence of these same series expansions, the multiples and primaries have the same region of convergence. sin(0) < c /c 0 < 2 Ultimately, given that sin(#) < Co/ci < (1 + cos (^)) ' 2 1/ 2 and (1 + c o s ^ ) ) / , the B o r n series solution for this problem converges. 2 1 2 For each velocity and a given angle, a curve similar to Figure 2.6 can be drawn. T h e curve of 1 — ( c / c ) 0 n 2 must then lie in the region of convergence. Hence for sufficiently low angle, the velocity contrast between the reference medium and the inner layers determines the convergence of the series. 2.4 Conclusions T h e work i n this chapter has established a mapping between a volume point scatterer model of seismic reflection data (and a series representation of that model) and well known solutions for primary and multiply reflected events. T h e latter represents the way seismologists traditionally think of events as being generated from wave propagation and reflections from interfaces. T h e former is the scattering theory point of view which forms the basis of recent work in multiple attenuation. Inverse scattering is a step by step process for in- 33 CHAPTER 2. FORWARD SCATTERING SERIES verting seismic data into earth parameters and the removal of multiples can be thought of as an intermediate step in that process. This work enhances our understanding of inverse scattering multiple attenuation by understanding how earth properties are mapped into primary and multiple reflections in the forward problem. Specifically, I have explored the role of the terms in the forward Born series by analytically demonstrating that the series properly describes all the reflected and transmitted waves for a simple one dimensional layered scattering problem. It is known that since the convergent forward B o r n series describes the data, the higher order terms in the series must describe certain data characteristics beyond the B o r n approximation (see, for example, Stolt and Jacobs (1981)). T h e calculations i n this chapter have not only confirmed that the role of the higher order terms is to alter the amplitude of reflected and transmitted waves, adjust the propagation velocity of transmitted waves, and describe interlayer multiples, but have also shown the details of how this occurs. For re- flection data, primary reflections depend on all of the terms in the series excluding the first term; whereas multiple reflections depend on all of the terms past a certain order. More specifically, a multiple which contains n changes in propagation direction (i.e., reflections) requires the T i ' t h and all higher order terms in the Born series for a full construction. For a simple layered model where c 0 is the reference velocity and c n is the velocity of the n ' t h layer, the Born series for a plane wave with incident angle 8 converges if c sin(#) < ci < c ( l + cos (#)) . n n 2 2 This implies that models with very low velocity layers will lead to divergent B o r n series which in turn has implications for inverse scattering. In the inverse series originated by Jost and K o h n (1951) and later developed by Moses (1956), the convergence of the forward series is necessary for the convergence of the inverse series (Prosser, 1980). This means that earth models which lead to divergent B o r n series have no hope of being recovered using this inverse series. Nevertheless, even when the forward series is divergent, certain convergent inverse subseries have been identified for attenuating multiples. T h e application of these series for multicomponent seismic data is the subject of the coming chapters. Chapter 3 Free Surface Multiple Removal for Land Data 3.1 Introduction This chapter deals with the application of inverse scattering methods to remove multiples associated with an elastic free surface. Although data contain both free surface and internal multiples, I delay the treatment of internal multiples until Chapter 5. T h e work in this chapter is motivated by a method which uses inverse scattering series to remove acoustic free surface multiples from marine seismic data (Carvalho et al., 1992; Carvalho, 1992). This previous work assumes that multiples are generated by an acoustic free surface situated above sources and receivers located in the water column. In land data, sources and receivers are located on or in an elastic medium and the unwanted multiples are generated by an elastic free surface. M y approach is to adapt the marine free surface demultiple method to remove elastic free surface multiples from multicomponent land data. 3.2 Forward and inverse elastic wave scattering T h e mathematics of elastic wave propagation involve various matrix or tensor quantities. W h e n written explicitly, these quantities often create unnecessary clutter which obscure 34 CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 35 the physics contained in the equations. To avoid this, I employ linear operator notation to represent the different elastic quantities. Wherever warranted, I will explicitly show salient portions of the equations but I defer the full detail to the appendices. T h e notation convention I use is as follows: • Lower case bold (normal and Greek) letters represent vectors with individual components in lower case italics. I will alternatively use indicial notation or explicitly write column vectors enclosed in brackets. r written as u, Ui, or u i x u y For example, the displacement vector can be T u z • Upper case bold refers to matrix quantities. Individual elements are denoted using indices, e.g., F xx is an element of F. Matrix or tensor quantities will also be represented using indicial notation, e.g., Fjj is equivalent to F. • Upper case calligraphic letters refer to displacement quantities, e.g., Q. • Upper case normal letters refer to quantities decomposed into compressional ( P ) and shear (S) wave components, e.g., G . These components are essentially the divergence and curl of the displacement field and are discussed in more detail in Section 3.5. • Upper case quantities written with a hat are operators, e.g., Q, G , and Gpp. Quan- tities without a hat are functions which are related to the operator having the same symbol, e.g., Q is matrix operator which is related to the matrix function £7(r |r ; u;). p s In this thesis, operators are integral operators which involve their associated functions. Explicit examples will be discussed below. • Upper case bold Greek letters (II) also represent matrix operators; however, there is no distinction between vertical and horizontal displacements and (P, S) components. T h e meaning of the quantities in the above examples will be made clearer as the theory develops. T h e first step is to derive the scattering equations appropriate for elastic wave propagation. I begin with the elastic wave equation for displacements in an inhomogeneous, CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 36 anisotropic and unbounded medium (see, for example, A k i and Richards (1980)) pui - (cijkiu j),j k = fi i = 1,2,3 (3.1) where Ui is the displacement of the continuum in the i ' t h direction, p is the material density, Cijki is the fourth order elastic stiffness tensor, and fi is a body force per unit volume which represents the source of wave energy. T h e directions 1,2,3 are x,y,z The solution for the displacements at a geophone located at r g respectively. due to a body force distribution / (r,u>) is (see, for example, Morse and Feshbach (1953)) m Ui(r ;uj)= / g v g (T \x']<jj) f (T',io)dT' im (3.2) m g — OO for a fixed temporal frequency UJ which implies a Fourier transform with respect to time. T h e quantity Qi is the Green's displacement tensor (or matrix) which satisfies the equation m -pW Gim{Tg\T'; U}) - {Ci lGkm,l{ 9\ 'i 2 where £ tTO T jk T = im6(r S g ~ T') (3.3) is the Kronecker delta. T h e Green's tensor represents the z'th component of the displacement due to a delta function source in the m'th direction. Equation (3.2) then gives the displacement due to an arbitrary distributed body force as the spatial and temporal convolution of the impulse response (Green's function) with the body force. A more compact way of writing equations (3.1)-(3.3) is to use operator notation. T h e elastic wave equation (3.1) i n terms of operators is Cu = f (3.4) where C is the elastic wave matrix operator which operates on the displacement vector u. For the reader unfamiliar with operators, Stolt (1981) gives a good introduction to linear operators and their properties. T h e main thing to remember is that operators operate on the quantity to their right. In addition, the operators i n this thesis will be described in terms of integrals over functions. This is discussed i n more detail i n A p p e n d i x D . The solution for displacements in (3.4) is u = 0f (3.5) CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 37 which is just (3.2) in operator notation. T h e Green's operator Q satisfies the relation Cg=± where I (3.6) is the identity matrix operator; hence the Green's operator is the inverse of the wave operator. Comparing (3.5) with (3.2) reveals that Q implies an integral over all space which becomes explicit when the operators are written in terms of the Green's function. Since displacement is a vector, the operators are matrices. In terms of the Green's function, the Green's operator can be thought of as an integral matrix operator which operates on some vector (or even tensor) function of space F j(r') m oo (3.7) g (r \T';u>)(F (r'))dr'. im / mj g oo One point of clarification is required. W h a t is the difference between <7V and ^(r)V(r)? T h e first is Q operating on V which involves an integral over space of the matrix product of their functions Q(r) and V(r). T h e second is simply the matrix product of two functions. In further developments, it is convenient to work with Q rather than u because solutions for any u can be constructed using (3.5). W h e n working with actual displacement data, the effect of distributed sources must be incorporated using f. Strictly speaking, operators are incomplete by themselves and need to operate on something to become a function which has physical meaning. Hence, actual measured data must be described in term of Green's functions. However, since the Green's operator Q and the Green's function ^j(r |r'; UJ) are g intimately related, (the Green's function is obtained by applying Q to a delta function, Q = QS) this justifies being somewhat loose in referring to operator quantities as data. The elastic Green's function represents the wavefield in a heterogeneous and possibly anisotropic medium; hence evaluating it can be difficult. A s discussed in Chapter 2, the scattering theory approach is to recast the wave problem into a simpler form by describing the earth as a background medium with a perturbation. T h e wavefield in the background or reference medium satisfies the operator equation Cg = X. (3.8) where Go is the reference Green's operator and C is the wave operator for the reference 0 0 0 medium. T h e perturbation describes how the actual earth deviates from the reference CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 38 medium. Mathematically, it takes the form V ~ Co — C which is the difference between the wave operators in the two media. Physicists often refer to V as the scattering potential. I will use the term perturbation to avoid confusion with the Helmholtz potentials I will be using later. Substituting V into (3.6) and rearranging terms leads to C Q=X 0 Next, I employ Go to invert C + Vg. (3.9) and arrive at the elastic Lippmann-Schwinger equation 0 G = G + G VG 0 (3.10) 0 or written in terms of functions oo / poo / Goij(r \r'; u/)Vjfc(r'|r"; UJ) £ f c ( r " | r ; UJ) J — oo m g CO 3 CLT'CLT" (3.11) where ( / i ( r j r , ; u;), Goim(*g\ s',<*>)) and Vjk(r'\r";UJ) m r are the functions related to the op- erators G, Go- and V- In (3.11), the perturbation depends on four spatial variables to accomodate a wide variety of earth models. This is done because, for example, in addition to earth properties changing in space, the additional degrees of freedom may be required to describe angle dependent scattering and anisotropic effects. T h e form of V and its implications is discussed i n more detail i n A p p e n d i x D and has also been discussed by Stolt and Weglein (1985). Although it is possible that a perturbation with more than four variables may be needed to describe very complex earth behaviour, I will assume that four are enough. This is still an open issue and warrants further investigation. The Lippmann-Schwinger equation is a fundamental equation which relates the wavefield i n the actual medium to the wavefield in a reference medium plus a perturbation wavefield which depends on Go, V , and G itself. T o isolate G from the other quantities, the Lippmann-Schwinger equation can be iterated to obtain the B o r n series G = 9 0 + Q VQo + 9oV9oVGo 0 + .. • (3.12) In this form, the perturbation wavefield is nonlinearly related to the reference wavefield and the perturbation. Generally, the perturbation represents the unknown quantity or the CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 39 earth model. For this reason, the data are defined as the portion of Q which depends on V . Thus T> = g VQ 0 + Q Q S?QoVQo + • • ( • (3.13) T h e quantity X> is, of course, a perturbation wavefield operator and yields displacement data in space and time (or frequency) when *D operates on an impulsive body force matrix. For simplicity, I will use the term data interchangeably to refer to either the operator or the function. To avoid confusion, I will refer to (3.13) as the forward scattering series or just forward series and (3.12) as the Born series. As discussed in Chapter 2, the forward scattering series describes the forward modelling of 2? using a reference wavefield and a perturbation V . O f course, using this series to do forward modelling would be terribly inefficient. However, as was also shown i n Chapter 2, certain insights can be gained by studying the forward modelling process. W h a t is ultimately of interest to geophysicists is not the forward problem but the inverse problem (although the two are intrinsically related). T h e goal of the inverse problem is to obtain information about the earth given the data. In scattering theory, this is called inverse scattering and many methods have been devised. T h e inverse scattering method I employ was originated by Jost and K o h n (1952) and later developed by Moses (1956) and also Razavy (1975). Although other inverse scattering methods exist that are applicable to 1-D acoustic wave problems (e.g., the Gelfand-Levitan method (Ware and A k i , 1968)), this is the only known method which can accomodate earth models that vary i n more than one dimension (Stolt and Jacobs, 1981; Weglein, 1985). A formulation of the Jost and K o h n method using an elastic background medium has been previously presented by Weglein and Stolt (1995) and I will follow closely their development. To formulate the inverse problem, an inverse series is postulated of the form V = V U + V U + V where the n ' t h order term is n ' t h order in T>. U + ... (3.14) A simple example illustrates why this approach is reasonable. Suppose A is related to B by A = 75/(1 - B) = B + B 2 + B 3 + ... CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA in analogy with the forward scattering series. the series B - A/(l + A) = A - A + A 2 s 40 T h e inverse relationship for B given A is + . . . . Consequently, since the forward process describes the modelling of data as an infinite series, it is reasonable that the inversion of data can be also expressed as an infinite series. In the first case, the data are given by an infinite series in the model. In second case, the model is given by an infinite series in the data. This simple example also clarifies why the data are the difference in the actual and reference wavefields. If, in the above example, A = 1 / ( 1 — 5 ) = 1 - f 5 + 5 2 + B 3 + ..., then there is no simple inverse series expansion for B in terms of A. Substituting (3.14) into the forward series gives a set of equations in various orders of T> T> = G V G 0 0 = g V Go {2) 0 0 = GoV Go i3) {1) (3.i5) 0 + GoV GoV 9o (1) (3.16) (1) + <?oV <?oV <? + <7oV <?oV <?o + GoV GoV QoV Go (2) (1) 0 (1) (2) (1) (1) {1) (3.17) These equations form the basis of the inversion process. T h e y suggest that in principle, it is possible to construct the various terms in V starting with "D and the reference wavefield Go- T h e first step is to obtain ^ from (3.15). This is then used i n (3.16) to obtain \ Continuing this process, all of the terms in V can be found iteratively starting with T>. To understand this inversion, it is important to make the distinction between V and The former is the difference in the wave operators between the actual and reference media and represents the difference in mechanical earth properties between these two media. T h e latter is the first term in the inverse series for V and is derived linearly from the data. To further interpret \ consider (3.15). Here, the data are modelled by waves ~ (!) propagating in the reference medium from the source to V with at depth. After interacting \ the waves then return via the reference medium to the receiver. Consequently, contains all of the characteristics of data: multiples and primaries with angle dependent CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA reflection variations. Since sandwiching and 41 ^ by Q implies upward continuing the sources 0 receivers to the surface, inverting these operators can be thought of as downward continuation. Thus, This can be thought of as migration using the reference medium velocity. ^ can be thought of as a reflectivity image in scaled time which is obtained by " (2) migrating the data using the reference medium velocity. T h e role of the terms V * (3) , V , etc., is to change this reflectivity image to V : the relative change in earth properties as a function of depth. This is conceptually how the inverse series performs inversion. A t this point, it is useful to take a step back and look at the larger picture. T h e job of the forward series is to model seismic data for a specified earth model. In this process, the forward series must perform four main tasks: 1) convert from depth to time; 2) transform earth properties into reflectivity (or some quantity related to it); 3) create free surface multiples; and 4) create internal multiples. It follows that the job of the inverse series is to undo the tasks performed i n the forward series. Specifically, the inverse series must 1) transform from time to depth; 2) convert reflectivity to earth properties; 3) remove free surface multiples; and 4) remove internal multiples. Moreover, these inverse tasks are * ( ) 2 contained i n the terms V ~ ( ) 3 , V , etc. Unlike the tasks in the forward series, it has been shown that tasks in the inverse series are separable to a greater or lesser extent. Carvalho et al. (1992) showed that the removal of marine free surface multiples is completely separable from the other tasks of inversion. For internal multiples, Araujo (1994) found a separable subseries that performs internal multiple suppression as opposed to outright removal. It may be that the removal of internal multiples can be formulated as a completely separable task. However, it is a difficult problem and a path has not yet been found. M y aim is to find a separable subseries which removes elastic free surface multiples. 3.2.1 Homogeneous isotropic reference medium and (P,S) potentials The question arises as to what to choose for the particular reference medium. T h e scattering formulation essentially takes a difficult problem and recasts it into a simpler problem with a perturbation. In this way, much of the complexity of the original problem is transferred CHAPTER 42 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA to the perturbation; therefore, the simpler the reference medium, the better. W i t h this in mind, I consider for now that the reference medium is an isotropic homogeneous whole space. I also consider only two dimensional wavefields. In two dimensions, the reference Green's displacement operator can be decomposed as (Weglein and Stolt, 1995), g (3.18) = n ^ G o i i 0 w here r - l 0 l/7o _ o — 0 0 = l//*o d d x n G n z -d d z 1 x =v GQP 0 0 Gos (3.19) -d d z x (3.20) d z Here, 70 is the uniaxial strain modulus, //o is the shear modulus, Gop and Gos are scalar P and S wave Green's operators, and V A = V~ B 2 - 2 is an integral operator (i.e. if V J 4 = B, then 2 where A is the solution to Poisson's equation). T h e Green's functions Gop and Gos satisfy the scalar wave equations (v" + ^0 2 ( where V a = (lo/p) ^ 2 1 0 2 G p(x ,z \x ,z ]Uj) 0 g g 3 = S(x - x )8(z g s s - z.), g + ^ 2 ) os{x , z \x , z ; UJ) = 8(x - x )S(z G g is the reference g P s 3 g s wave velocity and /3 = 0 g (3.21) - z.) (f-o/p) ^ 1 2 is the reference S wave velocity. The form of the Green's operator in (3.18) is reminiscent of the diagonalized form of a matrix. In fact, the diagonal matrix r " G o contains the independent eigenvalues of the operator Go and II o 1 contains the corresponding eigenvectors. T h e eigenvalues of a system represent its modes, which, i n this case, are decoupled P and SV waves. In three dimensions, an additional SH mode exists and the eigenvalues are slightly more complex. Since I a m considering two dimensional wave propagation, only the P and SV (or simply S) modes appear. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 43 Equation (3.18) can also be understood using the well known Helmholtz representation for displacements in two dimensional elastic media (see, for example, A k i and Richards (1980)) (3.22) u = V<p + V x (0,-0,0) where u is the particle displacement, cj> 1S the P wave potential, and ip is the S wave potential. B y taking the divergence and curl of (3.22), independent expressions for <p and ip can be found in terms of u. Equation (3.22) can also be written as d x u = d z -d 4> z (3.23) d x Now consider an experiment in the reference medium where a body force f gives rise to P and S waves. T h e displacements for this experiment are u = Q f = n ^ r j ^ G o I I f . 0 This can also be written as u = -1 -IF G 0 (3.24) n f pu>' using (3.21). Comparing this with (3.23), two things can be concluded. First, the Helmholtz potentials are related to G IIf by 0 ^GoIIf. (3.25) /our Second, II acts to decompose f into its P and S components [<pf,ipf]T giving u The lesson is that -1 pur II transforms (ux,uz) G (pf 0P (3.26) GosTpf displacements to (P, S) potentials while II -1 transforms them back. From this, the displacement u can be described as resulting from: 1) (P, S) transforming the f(x, z) source components; 2) propagation of these components via the decoupled (P, S) modes; and 3) transforming back to displacements at the receivers. This model of (P, S) wave propagation will be useful in the coming forward and inverse scattering series. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA I comment that translating V _ 2 G to — l/(puj )G 2 0 44 using (3.21) is valid only when 0 the measurement point is not coincident with the source point. This is not worrisome i n practice because when a receiver is placed on or next to a source, the physics of what is being measured is not modelled by the wave equation anymore. Moving back to the elastic scattering equations, substituting (3.18) into (3.15) gives for the displacements g = n ^ G c i i + n- To ' G o n v ^ n ^ G o i i . x (3.27) Using the lesson of the (P, S) decomposition, the individual operators have the following roles to play in the elastic scattering process II: transforms (u ,u ) x z to (<^>, tp), G : represents propagation of (</>, rp) from source to measurement point, 0 I I : transforms (<p,ip) to -1 V ^ : scatters displacements (u ,u ), x z (u ,u ). x z T h e scattering in (3.27) can then be represented by the diagram i n Figure 3.1. This may seem a lengthy interpretation of the scattering process; however, it leads to a simplified picture using only P and S waves. Premultiplying (3.27) by the operators To II and postmultiplying by the operator L T - 1 transforms the displacement Green's operator matrix to G = Tonga = 1 GPP Gps Gsp Gss (3.28) where the rows (subscript 1) represent the type of measurement and the columns (subscript 2) are the source type. Hence the elements of G are the P and S measurements of the displacement field due to purely P and S sources. After (P, S) decomposition, (3.28) is then G = Go + Go (nv^n-'ro ) G x Go + G V « G 0 0 0 (3.29) CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA where "V"W is ^ in (P, S) coordinates. Previously, the role of 45 ^ was to scatter the vertical and horizontal components of displacement where the sandwiching II operators converted the incoming waves from (</>, to (d>,ip). VO to (u ,u ) x z and the outgoing waves from This process has been compressed so that V ^ ) simply scatters P and 1 (u ,u ) x S z waves. The result is the simplified scattering diagram in Figure 3.2. Each element of V ^ ) has a specific job to do in the scattering process. To understand 1 what these elements do, consider the displacement field due to body force f. Performing the receiver decomposition, (3.5) becomes = r n 0 jo(d u + d u ) x u x z f*o(d u - d u ) x z z G IIf +GoVWG nf. 0 z 0 x = Gf (3.30) CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA Now suppose that f is such that it gives rise to P waves only i.e., I l f = 8 1 0 GOP 0 + GOP 0 T/(l) Gos (1) Vpp (i) GOP (1) 0 PS V. SP Vss 4> — GOP + GOP Vpp GQP , 0 46 . Then (3.31) so that (3.32) i> = Gos Vfp GOP • (3-33) This says that the measured P waves propagate as P from the source, are scattered into P by Vpp and then propagate as P to the measurement point. Also included is the P wave which propagates directly from the source. T h e measured S waves propagate as P from the source, are scattered and converted into S by Vgp and then propagate as S to the measurement point. Hence the role of the diagonal elements i n is to perform like CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 47 scattering ( P —*• P and S —> S) and the role of the off-diagonal elements is to perform converted scattering ( P —> S and S —> P ) . This interpretation has been made previously by Clayton (1980) in inverting the elastic B o r n approximation. A p p l y i n g the source and receiver decomposition to the entire inverse series gives V = V ( 1 ) + V< > + V < ) + . . . 2 3 (3.34) from which it follows that D = GoV^Go 0 = G 0 V ( 2 ) G 0 + GoV 0 = 0 V ( 3 ) G 0 + G G + 0 V ( ( 1 2 ) ) G V 0 G 0 V Go ( 1 ) ( 1 ) G 0 + G Q V ^ G Q V ^ G O (3.35) G O V W G O V W G O V ^ G O where D are data with sources and receivers i n terms of (P, S) potentials. T h e reader may question the value of this last result over equations (3.15)-(3.17) considering the effort expended. T h e point to writing the inverse scattering series i n (P, S) potentials is that it allows a more intuitive picture of wave scattering using decoupled P and S waves. In addition, it allows a ready analogy with the marine free surface multiple removal method. In the next section, I consider an isotropic homogeneous half space reference medium and derive the appropriate Green's function. This will prepare the way for obtaining a multiple removal subseries. 3.3 Free surface (P,S) Green's function L a n d data are generally acquired with sources and receivers on the earth's surface and a target located at depth. For this geometry, an elastic half space is a better choice for the reference medium than a whole space because the shallow earth and its surface are well approximated by a half space bounded by a free surface. This is also mathematically better CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 48 since the scatterers can be placed on one side of the measurement surface thus eliminating the need to keep track of scattering from many directions. Recall that the scattering problem transfers the complexity of the inhomogeneous earth problem onto a perturbation about a simple reference medium. Although a half space is not the simplest of reference media, there are significant gains in making the reference medium slightly more complex. In this section I derive the Green's function appropriate for a homogeneous isotropic half space which will serve as the reference wavefield for land seismic data. A similar development can be found in Achenbach et al. (1982). This derivation sets the stage for developing a scheme to remove elastic free surface multiples. In the previous section, I discussed the forward and inverse scattering series for displacement data. I showed that these series can be decomposed into P and S waves given an isotropic homogeneous background medium. Recall that the decomposed reference Green's function for an isotropic homogeneous reference medium contains both P and S wave Green's functions GQP and Gos respectively. These scalar quantities satisfy the scalar wave equations in (3.21). T h e solution for Gop (see, for example, DeSanto (1992)) is written as the W e y l integral 1 G {x z \x ,z u:) = A g) 0P where v= sign(a>)(u; /a: 2 2 3 g —k)/ 2 1 2 3] t°° 1 —J _ *(-.-.) H«.-.ldfc c c (3.36) is the vertical wavenumber (some refer to this quantity as the obliquity factor since it describes an angle dependent amplitude) and a is the P wave reference velocity. T h e sign(u;) ensures that v has the same sign as u> so that the wave propagates forward i n time and out from the source. This condition is also known as the Sommerfeld radiation condition. T h e superscript d indicates that the solution represents propagation directly from the source to receiver without any interactions. W i t h the source and receiver fixed, (3.36) can be considered to be the superposition of plane wave solutions ( g- ') \ '- g\ %k e x x w e z z scaled by l/2iv for a fixed wavenumber k. It can also be thought of as a spatial inverse Fourier transform over the geophone domain f°° Gop(xg,zg\xs,z3]uj) = — / 1 e » »G^p(k ,z \x ,z ]Uf)dk lk x g g a a g (3.37) CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 49 where Gkfo, = J- -*.-. -.l'.-*.l e (3.38) e 2iv g and v = v(kg). Alternatively, g GQ can be represented by an inverse Fourier transform over P the shot domain f°° 1 G"op( g, g\ si s;v) x z x = — J z e~ Gop(x ,z \k ,z ;uj)dk tk3Xs g g s s (3.39) s where G$ {x ,z \k.,z.;u>) = P and u = v(k ). s s g J-e^e-'""-^! g (3.40) This is a powerful formulation which can be used to build Green's functions in more complicated reference media. T o construct the Green's function appropriate for an elastic half space, consider Figure 3.3 which depicts an elastic medium with P and S wave velocities a and j3 and bounded at z = 0 by a free surface. T h e z axis increases downward. Suppose there is a source of plane waves which propagate upward towards the free surface with wavenumber k . T h e 3 plane wave which propagates directly from the source and receiver is ik xg ii/ (z ~z ) e s e a s (3-41) g W h e n this wave strikes the free surface, it creates a reflected P wave of the form pp ik Xg iv {z +z ) e g e a s g (3-42) where PP is the free surface P-P reflection coefficient (see Appendix B for the free surface reflection coefficients). T h e change in vertical propagation can be interpreted as a P wave image source. In addition to a reflected P wave, the incident P wave gives rise to a converted S wave p Z; ik xg i(v z +r) Zg) e where rj s = sign(u>)(u; //3 — A : ) / 2 2 2 1 2 a e s s s is the S wave vertical wavenumber. (3.43) T h e inclusion of rj means that the reflected converted wave propagates with the S wave velocity after s reflection. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 50 Free Surface (z=0) PS exp(i(k x, + v z + s a V» g)) z s converted S wave Figure 3.3: Reflected and converted waves for a P wave incident on a free surface. To build the incident wave due to a line source, the incident plane wave is scaled by l/2iv s and inverse Fourier transformed with respect to k . Since each reflected plane wave s solution satisfies the free surface boundary condition, so does the weighted superposition of these waves. Hence, scaling the reflected waves by l/2iv s and inverse Fourier transforming with respect to k builds the reflected wavefields due to a fine source in the presence of a s free surface. These solutions can be written in vector format as /-id ^OP + 0 ^OPP (3.44) /orfs ^OSP where poo 1 G%P (x ,z \*.,z.;u) P g B = ^ J ph — *'<"'--'> *'<*' "«><fc. e e (3.45) + and 1 G (x ,z \x ,z ^) f s 0 SP g g a 3 = poo p —J b —c*.^-.)^...^...)^ . (3.46) These solutions are also valid when the receivers are below the source location. T h e superscript fs denote waves which interact with the free surface. Together, these solutions represent the P and S wave Green's functions due to a pure P wave source i n an elastic half space. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 51 Free Surface (z=0) exp(i(k x + r] (z - z )) s s s s ^ g incident S wave / \ \ ^ SPexp(i(fc„z 3 / depth (z) / 1 + r) z + v,z )) s s g converted P wave SS exp(i(k x, + -q,{z, + z )) t g reflected S wave Figure 3.4: Reflected and converted waves for an 5" wave incident on a free surface. To build the free surface Green's functions for incident S waves, consider a plane S wave travelling with the velocity (3 and wavenumber k incident on the free surface as i n 3 Figure 3.4. This wave is written as ik x irj (z -z ) & a g e 3 s (g ^ ) g T h e interaction with the free surface creates a reflected S wave Sg ik,x ir,,(z +z ) (3_ g) §p ik xg i(r) z,+v z ) (3.49) e ge s g 4 and a reflected converted P wave e a e a s g where SS and SP are the S-S and S-P elastic reflection coefficients. Scaling the incident wave by l/2ir] s gives Gg (x , 0 g z \k , Z ;UJ). g a This is then Fourier transformed over k S a to obtain the directly propagating S wave due to a line source measured at the geophone i.e. ^ 0 5 ( ^ 9 ' g\ s, z x *\UJ). Explicitly Z t°° 1 1 G* {x ,z \x.,z.;u) s g g = —/ —e >(*°-*^-\'>-''\dk. ik (3.50) Scaling and transforming the reflected waves builds the P and S wave Green's function for an S wave line source due to the free surface. These are G * (x ,z \x z ;uj) { SP g g Sl s = - / _ *.(-.--0 *>.«.+*..)^ e l-K J_ 00 c 2lV s (3.51) CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 52 and GQ SS( 9I S ss g\ X z x (3.52) 2irj s Written i n matrix form the incident and reflected waves are 0 /~tfs ^OPS + Gs (3.53) /~<fs "-"OSS d Finally, the two P and S incident wave experiments can be written together as columns in a matrix Gopp GQSP G 0 GQPS d OP — Goss 0 G _ d •~ffs + U 0 P P ^OSP os ^OPS ^oss (3. Go — G g + G Q or more explicitly 1 G (x ,z \x ,z ;uj) 0 g g s s /°°° = — J-oo e » zk{x if\x,-Zg\ e 0 0 iv\*3-*g\ 2ir) e + pp ivz e g SPe <> iuz PSe > ° SSe <> iT z ir)Z 2iv \dk. (3.55) 2ir\ This is the P and S wave Green's function for a purely P or S wave source i n the presence of a free surface. T h e free surface Green's functions can also written i n terms of operators as Go = G 0 + G o • (3.56) In the next section, this will be used to identify the free surface demultiple subseries. 3.4 Elastic free surface multiple removal subseries In this section, I obtain a subseries which performs free surface elastic multiple removal by selecting a portion of the full inverse scattering series. T h e key to the identification of CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 53 this subseries lies i n the special form of the (P, S) free surface Green's function obtained in the previous section. Note that the operator G o has two portions. T h e first part, GQ, represents propagation directly from the source to receiver in the reference medium and is the same (P, S) Green's function as for a homogeneous isotropic background. T h e second portion, G 0 S , exists due to the presence of a free surface i n the reference medium. In the forward scattering series, G Q is responsible for creating events associated with a free surface i.e. ghosts and free surface multiples. If G o did not contain G 0 S , the forward series would create data absent of any free surface effects. This is i n direct analogy with the marine case where the reference Green's function has two portions: one which represents direct propagation and another which represents reflections from an acoustic free surface. Understanding how free surface events are generated in the forward series gives a clue as to how these events are removed in the inverse series. Since G 0 S is solely responsible for the creation of free surface events, it must be responsible for the removal of these events i n the inverse process (Weglein et al., 1997). Therefore, terms which contain G„ are separated from the full inverse series as a candidate subseries whose task is to remove free surface events. Conceptually, the inverse series is partitioned into V = V( V « + V'< > + V ' W + ...) 2 V V' ' +V 2 ) + V 3 ) + ... (3.57) where terms i n brackets are the free surface multiple removal subseries V ' and the remaining terms perform the other tasks of inversion (converting from time to depth, mapping reflectivity to mechanical earth properties, and removing internal multiples). T h e quantity V ' represents a series of terms each of which removes a certain order of free surface multiples when added to V W . A s i n the previous section, V ^ ) is the portion of V linear in 1 the model; it contains all primaries, internal and free surface multiples. It can be related to reflectivity since it displays all the characteristics of data yet lies in model space. T h e CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 54 terms which perform multiple removal are V'( ) 2 V = -V ( 1 ) G 0 S VW '( ) = - V ' ^ G ^ v W - V ^ G 3 = S 0 V' ( 2 ) - V ( 1 ) G * V ( 1 n n - 1 G £ V ( 1 ) (3.58) V ^ ^ G ^ V ^ ^ G ^ V ^ ) V'( ) = — V ' ^ ) ^Go V^ ^. 1 The recursive nature of these terms is due to the special structure of G Q (see A p p e n d i x D )• The next step is to reconstruct data using the new V ' series by sandwiching it between two directly propagating Green's functions. T h e result is a new series of data where each term has an identifiable role, D ' = GjJV'Go 1 = G 2 ( V W + V'< ) + V ' ( ) + . . . ) Gtf 2 = IV ) 1 3 + D ' ( ) + D'< > 2 3 (3.59) +.... The first of these terms is D (i) = = G J J V W G * G * ( G * + = Rg where R g D R GlT'MGi + GlT'Gt (3.60) s performs receiver (geophone) deghosting and R s (See Appendix C for a full description of these operators.) performs source deghosting. Note that the deghosting operators are applied to the original (P, S) decomposed data. T h e result is a new data set D W with the source and receiver ghosts removed. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 55 Using (3.60), the second term in the demultiple series is D'(2) -G*VMG*VWG = (Go) = _ 1 D^^(Go) _ 1 (3.61) - ^ ( G ^ G ^ G ^ D W = - D W H D W . and i n general (3.62) T h u s , each term i n the multiple removal series is obtained from the original multiple contaminated data L r ) and an operator H which contains elastic free surface reflection co1 efficients and appropriate obliquity factors (see Appendix D ) . T h e operator H is directly analogous to the acoustic free surface reflection coefficient and the obliquity factor term 2iv in the marine algorithm. Each term in the multiple removal series removes a certain order of elastic free surface multiples. In general, the re'th term removes all (re — l ) ' t h order free surface multiples and alters higher order multiples in preparation for their removal by higher order demultiple terms. Since the data are limited in time and space, there are not an infinite number of multiples on each record. A s a result, it is not necessary in practice to calculate all of the terms in the series; there is little use i n predicting and removing multiples which arrive after the end of the recording time. Explicitly, the multiples are calculated in the (k ,k ,u>) g a domain according to the formula H (k,uj) H (k,u>) D £,{k,k ,u) H p(k,uj) H s{k,uj) D p (k,k ,uj) D%>(k,k ,w) PP S PS S { s { p D^ {k,k ,u>) s s s dk. (3.63) a Note that the operator H becomes a simple matrix multiply by the function H . See A p p e n dix D for the details of how this occurs. A flowchart in Figure 3.5 summarizes the processing stream for elastic free surface multiple removal. This includes (P, S) decomposition and deghosting which, as discussed in Section 3.5, are performed in one step. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 56 Start Input 4 component displacement data with source deconvolved and groundroll removed: T>(£ , £«; 3 3-D Fourier transform over shot, geophone and time: T>(k , k ;w) g s Transform to (P, S) and deghost: D( )(^,fc ;o;) 1 s Calculate n'th order multiples: D( )(* fc,;w) n S) False Sum terms: D'(k , k ;u>) g s 3-D inverse Fourier transform over shot, geophone and time: T>'(x , x ;t) g Output multiple free 4 component (P, S) data: B'(x , g s x ;t) s End Figure 3.5: Flowchart of steps in elastic free surface multiple removal. A significant decrease in computation time can be realized using data from one dimensional earth models. In such a case, the data are invariant with respect to the source and receiver midpoint. T h e spatial transform variable for the midpoint is k g — k. data are invariant in midpoint, the only nonzero data values are where k — k s g 3 Since the = 0 or where CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 57 kg = k . A s a result, the data can be written as T)'( \k; UJ) 8(k — k ). Substituting this into n s g (3.62) gives D'Wfojw) = -B'^\k u,)H(kg-,u,)r)V(kg-,u;) (3.64) g] which is a simple multiplication in (k , UJ) space. Since the earth doesn't vary with midpoint, g one shot is as good as the next; hence only one shot is needed for multiple removal. It is reasonable to ask why series. G Q V ' G Q creates data when G Q V G Q does not in the forward This can be understood by looking at the multiple removal series term by term. T h e first term is the data with the source and receiver signatures removed. T h e remaining terms act to remove the free surface; hence they are the multiples with the sign flipped. Since the multiples and the primaries exist i n data space, the whole series can be thought of as a sum of terms i n the original data which produces new data containing just primaries and internal multiples. This explains why V is more like V « than V . T h e ra'th term in the multiple removal subseries represents a modelling of events which have undergone n — 1 interactions with the elastic free surface. This term can be thought of as an n — l ' t h order temporal and spatial convolution of the multiple contaminated data with themselves. For example, D'( ) is a convolution of D^ ) with D^ ). A good way to 2 1 1 understand how this process models first order free surface multiples is to use the concept of subevents. Basically, a given multiple can be decomposed into a series of subevents. A l l of these subevents are recorded as separate and distinct events i n the original data and can be primaries, lower order multiples, or even internal multiples. In the multiple removal algorithm, a multiple is estimated by reconstructing its propagation path and amplitude using its subevents. This tool has been used by Weglein et al. (1997) for marine free surface and internal multiple removal and also by Verschuur (1991) and Riley and Claerbout (1976). To illustrate the principle of subevents, I consider a simple elastic model consisting of a free surface and two reflectors. For the sake of simplicity, surface sources and receivers emit and measure P and/or S waves only. Suppose that there is a first order PPPP peg- leg multiple which can be decomposed into two PP subevents as i n Figure 3.6a). These subevents exist in the data as PP primaries with a source and a receiver at the end points of each subevent. W h e n the data are convolved with themselves, these subevents are strung CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 58 PPPP. A minus together with the free surface reflection coefficient to create the multiple sign in (3.61) ensures that when this estimate is added to the data, it cancels the actual multiple. T h e principle of subevent construction requires that each subevent is recorded in the data - no subevent, no estimate. This becomes important when near offsets are missing in the data. In general, the subevents for a multiple at a given offset are recorded at a nearer offset. Hence, if there is a gap in the near offsets, subevents which would otherwise be recorded in the gap are missing and multiples which require these subevents cannot be constructed and removed. Therefore, near traces must be either recorded or interpolated. Figure 3.6: Subevent construction of free surface multiples. T h e sources and receivers are represented by * and V respectively. The need for multicomponent source and receiver data becomes clearer under the subevent construction principle. phase (PPSP) Consider a first order peg-leg multiple which has an S (see Figure 3.6b)). This is partitioned into a PP and an SP subevent CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 59 which must be recorded in the data for effective multiple modelling and removal. Thus it becomes necessary to measure just the S waves at the receivers. B y extension, S waves sources must be used in order to remove multiples with SS subevents (see Figure 3.6c)). Another aspect of the multiple removal algorithm can be understood using subevents. Equation (3.63) states that in order to estimate a single multiple at a single shot and receiver wavenumber, the entire range of shot and receiver wavenumbers must be used. Consider now a single dipping layer and a free surface and a first order as in Figure 3.6d). A s before, this multiple can be split into two PP PPPP multiple subevents. In order to construct the estimate of the multiple, subevents with a specific source and receiver location must be used. Since the dip and location are unknown, the whole range of shots and receivers are used. Those subevents which do not construct the multiple contribute zero. D a t a are, of course, not measured in terms of purely P and S waves and must be transformed into the P and S domain prior to multiple removal. T h e next section shows how (P, S) decomposition may be done in tandem with deghosting in an elegant and efficient way. 3.5 Free surface (P,S) decomposition and ghost removal T h e first two steps in preparation for multiple removal are: 1) perform a decomposition from (x, z) to (P, S) coordinates and 2) remove the source and receiver ghosts. A s outlined in Section 3.2.1, the decomposition is done by calculating the divergence and curl of the displacement field. This becomes a simple process in the frequency domain where derivatives translate to multiplication by an appropriate wavenumber. Differentiating the wavefield in the horizontal direction is natural since the wavefield is sampled in that direction. T h e problem lies in performing the vertical derivatives. First, the velocity of the near surface must be known. Second, the sign of the derivative will change with up or down going waves. This suggests that it is necessary to first separate the two wavefields by deghosting. This is CHAPTER 60 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA problematic because the deghosting is tied to the wave type (P or S). Perhaps then (P, S) decomposition should be done first but this depends on deghosting ... Fortunately, a paper by Wapenaar et al. (1990) shows how to avoid this sort of chicken and egg relationship by performing the (P, S) decomposition and deghosting in one step. W h a t follows is m y adaptation of this approach to the elastic scattering problem. I start with the equation for elastic displacement data in (3.15), 2> = g V Go • (1) 0 (3-65) The reference Green's operator can be expressed in the form = —^n G n = Go r 0 p<jj* G )n. + —-^II (GQ T pw' (3.66) s 0 This follows from the fact that away from the source, V GOP V~ G p In matrix notation then 2 0 = — a /uj G p2 2 Similarly V~ G s 2 0 v - 2 n T r 0 = —(3 /v G s- - i £ o 2 = 2 z l pu n 0 r 2 G 0 2 = —UJ j O?GQP\ hence 2 . (3.67) Using this form of the Green's function, the data can be written as i> = — ^ ( G o + G*yVM{G* 1 puj 2 + G )n . s 0 (3.68) Factoring out G Q leads to X> = { l + G f t G S n G ^ G & l + ( G ^ G ^ I I pu = N^G^VWGgM; - ± n T 2 (3.69) 1 where N " 1 = ^Ti [1 J (1 + ( G O ) - 1 T + Go(Go) ) is the receiver decomposition operator and M -1 G Q ) I I is the source decomposition operator. Recall that the quantity G Q V ^ G Q is just i M ) , the deghosted data in (P,S) coordinates. Therefore, the operation 1 D = N T>M g g (3.70) transforms the original displacement data into data with upgoing P and S waves at the receiver and downgoing P and S waves at the source. These data are used as input into the multiple removal algorithm. CHAPTER In 3. FREE (k ,k ,u>) g s SURFACE MULTIPLE REMOVAL FOR LAND 61 DATA space, the decomposition operators are simply matrix functions which are multiplied by the data i.e., D(k ,k ,uj) g = a (3.71) TSl (k ,u>)'D(k ,k„u>)M (k «>). g g g g tt This is a pleasing result since it yields a simple transformation which is applied to displacement data in the (k ,k ,w) g domain. T h e required parameters are the source and receivers s depths (if nonzero) and the P and S wave velocities of the near surface. W h e n the sources an d receivers are located on the surface, the source and receiver functions i n (k ,k ,oj) g s space are 2 Ng{k ,Uj) g 2u„ (3.72) 2r) (3.73) K -i/3 '9 g g 2r, n -fc K 2 2 g and M,(fc„w) = — kg - v s n -fc 2 2 2u s See Appendix C for the derivation of the decomposition matrices. Synthetic examples of the decomposition results are shown in Section 3.7. 3.6 Incorporating wavelets and source orientation In practice, multicomponent sources may not be aligned along the x or z axis or may not even be orthogonal to each other. Furthermore, each source has a source time function which must be incorporated into the multiple removal algorithm for real data. This source time function is commonly called the wavelet. T h e body force which describes these new point sources can be written as ^(r^)^ = S(r )^F,A~ 3 S(x ) s 0 T 0 8(z )_ s J J~ZX Written i n (x, z,u>) space, this becomes 0 x z 0 A (UJ) 2 (3.74) CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 62 where the columns of T indicate the vector orientation of each source, and Ai(u>) and A (u>) are the two source wavelets. T h e data for 8T are given by 2 i> = (g - g )sr 0 = g vg ^ 0 (3.75) + dovgovgoF 0 +. where 2> are the data for arbitrarily oriented point sources. Note that the delta function has transformed the operators to functions. Postmultiplying this equation by the inverse matrix T gives 1 = g vg 0 0 + g vg vg 0 0 0 + (3.76) = T>. This is the same equation which led to both (P, S) decomposition and free surface multiple removal in terms of "D. It follows that the acquired field data "D must be transformed to "D prior to decomposition and multiple removal. This transformation is not difficult provided the wavelets and the relative amplitude and orientation of the sources are known. Using (3.74), the transform i n (X,Z,UJ) space is A " (a;) 1 0 0 A?(u) J~mx-Fzz J~mz-Fz: Tzz (3.77) ~Tz: Unfortunately, the source wavelets are practically never known. T h e y must be either estimated or recorded; these approaches have problems of their own. One possible solution to the wavelet problem is to assume that the wavelets are unknown but the same. In fact, the wavelets will generally not be same; but their differences can be reduced by applying spiking deconvolution (see, for example, Yilmaz (1987)) to both data sets as a preprocessing step. Another lever may be applied if the sources and receivers are vertical and horizontal. In this case, reciprocity and a shaping filter can be used to match the xz and zx components. Assuming that the source wavelets are the same or have been equalized, the sources may be oriented in the desired x and z directions. T h e result is "DA' = "D where A is the 1 remaining wavelet and T) is the field data corrected for source orientation. This gives the relationship between the ideal data and the field data with a source wavelet. corrected field data are used i n the (P, S) decomposition, the result is D J 4 _ 1 W h e n the : the (P, S) CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA decomposed data scaled by the inverse of the wavelet. Next, I use DA 1 63 = 2? to rewrite the multiple removal equations in terms of the (P, S) decomposed field data. 6' = 6 ( , ) + 6' " + < fi (,, + ... (3.78) where LS = G V^G A d (1) d f . = N X>M A = N X>M , 0 0 g rV s = -D ( 2 ) ( 1 ) g H^D (3.79) s ( 1 ) , (3-80) and in general 6' " < The , = -£, ( I ) iU-'£. ( 1 ) . (3.8i) result is a set of equations which are identical to the equations for ideal data except that each iteration requires the removal of the source wavelet. T h e advantage to this is that it is possible to first calculate the multiple removal terms without knowing the wavelet, and then apply an inverse wavelet to eliminate the multiples. This is in direct analogy with the marine case. Another aspect that can be borrowed from the marine case is to use the multiple removal method to estimate the multiples. T h e idea is that the data contain less energy without the multiples than with them. Using a minimum energy criteria, the wavelet is then found which best removes the multiples from the data (Verschuur, 1991). This is the true source wavelet. Various realizations of this concept have been demonstrated for marine seismic data (Carvalho et al., 1994; Ikelle et al., 1995) and are also applicable to the land case. In practice, the wavelet derived from the demultiple algorithm will not be the actual source wavelet. This is caused by various things: source and receiver array effects; source and receiver variations from location to location; and the data are 3-D where the algorithm assumes they are 2-D. Therefore, it is just as well that the source wavelet is not known exactly since it would likely not give the best multiple suppression. To first order, the source wavelet should be similar to the optimum multiple removal wavelet. A useful approach is to use an a priori estimate of the source wavelet as a starting point for finding the optimum multiple removal wavelet. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 3.7 64 Synthetic data examples In this section, I demonstrate the removal of elastic free surface multiples using synthetic data examples. D a t a were generated using a two-dimensional elastic finite difference code. Simple models are used since my purpose is to demonstrate the principles contained in the theory. T h e data processing was done using the Seismic U n i x (SU) system from T h e Colorado School of Mines. T h e decomposition and multiple removal codes were written as modules to the S U system. 3.7.1 Model 3 T h e first case I consider is Model 3 shown i n Figure 3.7. Model 3 has multicomponent point sources and receivers affixed to an elastic free surface. Since there is a single elastic reflector at 80 m depth, the data contain only a single primary; any other events are necessarily free surface multiples. A high contrast in elastic properties across the reflector ensures many orders of free surface multiples. Since this model is one dimensional, only a single four component shot gather is needed to calculate and remove the multiples (see Section 3.4). T h e receiver spacing is 10 meters (m) and the shot geometry is split spread (shot in the middle) with a maximum offset of 640 meters. T h e trace length is 256 samples at 2 milliseconds (ms) to give a recording time of .512 seconds (s). wavelet was used as input to the finite difference code. A m i n i m u m phase source T h e sources are horizontal and vertical point sources on the surface and contain the same wavelet. T h e measurements are vertical and horizontal velocities on the surface. Instead of treating the data as velocity measurements, I consider them to be displacements convolved with a source wavelet which is the time derivative of the original wavelet. T h e two representations are equivalent. Figure 3.8 shows the 4 component displacement shots arrayed in panels. Since the theory for multiple removal requires that the scattered wavefield is used as input for the multiple removal scheme, the reference wavefield i.e., direct waves and the groundroll, has been removed in these panels. A t first glance, it is striking that for such a simple model, the data are complex by comparison. I have labelled the four primaries (PP, SP, PS, SS) and some of the more important multiples. Read from left to right, the label tells the CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA multicomponent multicomponent source receivers * 65 V V V Elastic Free Surface Elastic solid a = 2000 m / s = 1200 m / s po = 1.0 g/cc 0 z = 80m Elastic solid 2 1 = f«« ft = 2500 m / s />! = 2.4 g/cc Figure 3.7: Land data Model 3. history of an event as it propagates through the model. Since the sources and receivers are aligned along the x and z axis, the measurements are not always i n the direction of wave propagation. A s a result, there is significant cross-talk between different components at m i d to far offsets. A n f-k filter has been applied to these data with a cutoff velocity of cto to exclude evanescent P wave energy. Since the evanescent/nonevanescent ranges of P and S waves overlap, the cutoff velocity is at the expense of some nonevanescent S energy. Fortunately, this effect is not severe. The first step is to decompose the data into downgoing P and S waves at the source and upgoing P and S waves at the receiver. T h e resulting data panels i n Figure 3.9 show that the cross-talk has been significantly reduced by rotating the different primaries into their respective (P, S) components. Next, the decomposed data, are input into an algorithm which calculates the next multiple removal data term T)'( \ These i n turn are 2 used as input to calculate the next term D'( ). 3 Multiples are removed order by order by summing the terms calculated from the multiple 66 CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 0.4 • V V zz Figure 3.8: xz V V zx xx Model 3 multicomponent displacement data with the reference wavefield (groundroll) removed. T h e first and second subscripts in T> indicate the receiver and source orientation and the maximum offset is 640 m. The primaries and some of the prominent multiples are labelled according to their P and 5 history from source to receiver. history is read from left to right e.g. PPPS begins as a P wave, reflects twice as This P, is converted to S, and is then measured on the surface. contaminated data. In Figure 3.10, I demonstrate the removal of multiples from PP The first panel is the PP portion of the input data. data. T h e second panel is the estimate of the first order free surface multiples that begin and end as P. W h e n the second panel is added to the first, all first order free surface multiples are removed, higher order multiples are altered, and the primaries are left untouched. Similarly, the data in the third panel removes all second order free surface multiples. W h e n all three data are summed, the result is D ' in the fourth panel. These data are free of first and second order multiples and also demonstrate how the higher order multiples are altered in preparation for their removal by D'( \ D'( \ etc. It is also evident that the primaries remain intact. This is important 4 5 FREE SURFACE MULTIPLE 3. Figure 3.9: Model 3 data decomposed into downgoing P and S waves at the source and upgoing P and S waves at the receiver. REMOVAL FOR LAND DATA 67 CHAPTER T h e first and second subscripts in indicate the receiver and source type. for processes such as amplitude variation with offset ( A V O ) analysis which depends on primaries which are free of residual multiples or noise introduced in the multiple removal. O f course, the demultiple algorithm removes multiples for PS, SP, SS data as well. A l l four components after two passes of multiple removal are shown in Figure 3.11. PP Like the example, the first and second order multiples have been removed, the higher order multiples have been altered, and the primaries are left untouched. this example are around -20 d B ( Suppression rates for 90 %). Effect of near surface velocity errors A s I stated previously, the properties of the near surface material ( a , /?o, Po) are required as 0 input to the multiple removal algorithm. Unfortunately, the near surface can vary consid- CHAPTER Figure 3.10: 3. FREE SURFACE M o d e l 3 PP MULTIPLE multiple removal. REMOVAL T h e PP FOR LAND portion of the input data is First and second order elastic free surface multiples are removed from D D P P 68 DATA . After summing the first three panels, the resulting data D P P p p by D p Dp . p p and contain multiples higher than second order which are removed by higher order multiple removal terms. erably which means that errors in estimating its properties are inevitable. To understand the sensitivity of elastic demultiple to these errors, I reran the previous examples using the wrong near surface velocities. Instead of the actual P and S velocities of 2000 m / s and 1200 m / s , I used 1600 m / s and 900 m / s . Velocity errors have a noticable effect on the (P, S) decomposition. A s shown in Figure 3.12, significant crosstalk remains between components. This is to be expected since the transformation depends on the vertical and horizontal derivatives of the displacement field. In Fourier transform space, the horizontal derivative is simply a multiplication by ik, the horizontal wavenumber. T h e vertical derivative, however, is a multiplication by iv (or ir/) which depends on the a (or f3). Therein lies the source of error in (P, S) decomposition. T h e corrupted (P, S) data were used along with the wrong velocities to calculate the CHAPTER 3. Figure 3.11: FREE SURFACE Model 3 (P,S) MULTIPLE REMOVAL FOR LAND 69 DATA decomposed data after two passes of multiple removal. Re- maining multiples are removed by higher order demultiple terms. multiples. Despite the obvious error in decomposition, the predicted multiples are correct in time but wrong in amplitude. Nevertheless, the multiples can be largely removed by using a scaled version of the source wavelet. T h e results are shown in Figure 3.13. Compensating for velocity errors by altering the wavelet can be understood by looking at the components of H . For example, the first component is PP2ivA(u>)~ . 1 incident angle is small, this becomes — 2iui / ocA(ui)~ . l angle H pp becomes —2iuj/a(A(uj)~ /e) 1 = W h e n the Now if the velocity ea is used, small — liuj/aA^uj) . -1 Therefore, when the angle is small enough, a scaled wavelet compensates for errors in near surface velocity thus giving effective multiple suppression. A slight complication arises for components of H which depend on rj instead of v. These components are sensitive to errors in (3 and therefore require a differently scaled wavelet. This is because the percentage errors in a and 0 are unlikely to be the same. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 70 W h e n the incident angle is large, the wavelet compensation becomes less effective and the multiples at far offsets are suppressed instead of removed. In the end though, the results in Figure 3.13 are very good. Figure 3.12: Model 3 (P, S) data with velocity error. T h e true P and S wave velocities of 2000 m / s and 1200 m / s are misestimated as 1600 m / s and 900 m / s . Separation of P and S events is degraded. Effect of missing data components It is rare that full four component data are acquired in the field due to the added expense of orthogonal sources. C a n anything be done with data having less than four components? To answer this question, I consider data with multicomponent receivers but only a vertical component source. Using reciprocity, these data lack only the xx component. Encouraged by the robustness of the demultiple algorithm given errors in reference velocities, I set the missing data component to zero and ran the incomplete data through the algorithm. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 71 Figure 3.13: Model 3 (P, S) data after two passes of multiple removal with velocity error. T h e true P and S wave velocities of 2000 m / s and 1200 m / s are misestimate d as 1600 m / s and 900 m / s . Compensation for velocity error is done using an altered source wavelet. T h e resulting multiple suppression is most effective at near offsets. Not unexpectedly, the (P, S) decomposition suffers as a result of having incomplete data. Significant cross talk is present across the different components in Figure 3.14. Despite the degraded ( P , S) data, the multiple removal algorithm is still effective. Figure 3.14 shows the four data components after first and second order multiple removal. While the overall results are very good, some of the multiples have not been removed. Recall that the method constructs an estimate of a multiple using subevents recorded in the data. Since some of the subevents are missing, multiples which contain them cannot be removed. CHAPTER Figure 3.14: degraded. 3. FREE SURFACE Model 3 (P,S) MULTIPLE REMOVAL data with missing V . xx FOR LAND Separation of P DATA and 5 events CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA Figure 3.15: Model 3 (P, S) data after two passes of demultiple with missing V . xx 73 Multiple suppression is still effective although certain events cannot be estimated and removed. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 3.7.2 74 Model 4 T h e second case I consider is Model 4 shown in Figure 3.16. M o d e l 4 is same as Model 3 except with the reflector dipping at 10 degrees. Since this model is two dimensional, many shots and receivers are needed to calculate and remove the multiples. T h e sources and receivers are located every 12 m and are referenced according to absolute coordinates; hence the ends of the model survey are at 62 m and 818 m with stations every 12 m . A receiver spread of 32 geophones was laid out split spread giving a maximum offset of 180 m . Since the spread is not allowed to go beyond the edges of the survey, shots at the beginning and end of the survey roll into the spread. T h e data were modelled using a sample interval of 2 ms and a recording time of .512 seconds. T h e data were arranged into a cube of 64 shots by 64 receivers by 512 time samples. M u c h of this data cube is zero padding which is useful when working with Fourier transforms. multicomponent multicomponent source receivers V V V Elastic Free Surface Elastic solid a 0 z = 0 = 2000 m / s /?o = 1200 m / s po = 1.0 g/cc z = 146 m at shot 446 Elastic solid ai = 4000 m / s ft = 2500 m / s px = 2.4 g/cc Figure 3.16: L a n d data Model 4. In Figure 3.17, I show shot 446 of the modelled data in (x,z) coordinates with the groundroll removed. Considerable crosstalk of the different P and S arrivals can be seen CHAPTER 3. FREE SURFACE Figure 3.17: between components. MULTIPLE REMOVAL FOR LAND 75 DATA Model 4 displacement data, shot gather 446. Due to the nature of the dipping reflector, the multiples are non periodic and increase their delay time as the reflector gets deeper. Figure 3.18 shows this same shot gather after source and receiver (P, S) decomposition and up-down going wave separation. M u c h of the crosstalk has been eliminated in these data particularly for the primaries. These data are then used as input into the multiple removal algorithm. T h e data after two passes of free surface demultiple are shown in Figure 3.19. figure, all second order multiples have been removed or severely suppressed. In this Residual events are higher order multiples or small processing artifacts at the edges of the data. These edge effects are introduced in the decomposition process as a result of finite offset. Since conventional methods become less effective for dipping earth layers, this example demonstrates the power of the inverse scattering method. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 3.8 78 Multiple removal for displacements In Section 3.7, I showed that multiple suppression is still effective with errors in the near surface velocity despite the degraded results of the (P, S) decomposition. While it is satisfying that the multiples are removed, the resultant primaries are distorted due to the erroneous decomposition. A solution to this problem is to remove the distortion by applying the inverse (P, S) decomposition giving the original primaries back in (u ,u ) x z space. Instead of doing this inverse transform numerically, I reformulate the multiple suppresion algorithm so that the forward and inverse decomposition is done implicitly. This avoids numerical errors associated with applying the forward and inverse decompositions. Referring to Section 3.5, the (P, S) decomposition is given by D* ) = NT> 1 M . Substi- tuting this into equations (3.61) - (3.62) gives f>'() = N X > M H N £ M , 2 (1) D'(3) = N2? ( 2 ) f)'(") = N2? ( n (1) MHNX> X ) ( 1 ) MHN£ M > ( 1 ) (3.82) M . Next, I apply the inverse (P, S) decomposition to these equations to get T>' (2) = £ ( 1 ) MH1NPP = X> ( 1 ) (1) K2> (1) , (3.83) anc T>' = where K = M H N . + X>' (2) + T>' {3) + ... + X>' (n) . (3.84) T h e result is a multiple removal algorithm that uses displacement data as input and gives displacement data as output. T h e previously explicit (P, S) decomposition is implicit i n the operator K . T h e four component displacement data after 2 passes of displacement multiple removal are shown in Figure 3.20. CHAPTER Figure 3.20: 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 79 Model 3 displacement data after two passes of multiple removal. Instead of performing the forward and inverse decomposition explicitly, it is done implicitly with the input and output data in displacement. One of the advantages of using the displacement formulation is that it allows an equitable comparison between elastic free surface demultiple and acoustic free surface multiple removal on single component data. In other words, one can compare apples to apples. In Figure 3.21, I compare three cases of elastic demultiple with acoustic demultiple run on single component data (vertical source and receiver). T h e first panel shows the result of two pass elastic demultiple using full four component (4-C) data. In comparison with the single component acoustic result in the fourth panel, the acoustic method performs admirably considering the data are recorded on an elastic medium and the free surface is not acoustic. Using subevents, the performance of the acoustic method is understandable since the vertical component data does in fact contain all the needed subevents for multiple removal. T h e acoustic method fails to remove events which are converted at the free surface. This makes sense because the acoustic method assumes that an upgoing wave reflects from the CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 80 0.1 0.2 • 0.3 • 0.4 0.5 1 1 - C elastic 2-C 4-C 1 - C acoustic Figure 3 . 2 1 : Model 3 comparison of elastic vs. acoustic 2 pass demultiple in displacements. T h e left data panel shows the full four component elastic demultipled data in displacements. T h e middle left panel uses only vertical source two component elastic data. T h e middle right panel shows single component data treated as one component of an elastic wavefield. In the far right panel, single component data are treated as an acoustic wavefield. free surface with an amplitude of - 1 and that the angle of reflection is equal to the angle of incidence. Since this is certainly not the case for converted waves, the acoustic method fails to remove them. Nevertheless, acoustic free surface reflection is approximately correct for events without conversion and near vertical incidence. Thus, the acoustic method gives good but not perfect results. In the second panel of Figure 3 . 2 1 , I show the results of elastic demultiple using two component ( 2 - C ) receiver data with a single vertical source component. T> zx can be obtained from V ; xz hence only T> xx B y reciprocity, are zero. T h e 2 - C results are very similar to the 4 - C results with small differences at far offsets. As a final test, I consider only single component data. B y treating these data as one CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 81 component of an elastic wavefield, I obtained the demultipled data shown in the third panel. These compare favourably with the acoustic result. T h e main conclusion from these tests is that the elastic formulation gives comparable (perhaps slightly better) results for multiple removal using single component data and superior results using multicomponent data . 3.9 Discussion and conclusions In this chapter, I have developed a method to remove free surface elastic multiples from multicomponent land data using inverse scattering series. This work is an adaptation of a method developed for marine seismic data (Carvalho et al., 1992). T h e method shown here removes all orders of elastic free surface multiples without requiring information about the earth which generates the reflections. T h e formulation iteratively generates a series of terms using the multiple contaminated data as input. E a c h term in this series removes a certain order of free surface multiples. Specifically, the n ' t h term removes the (n — l ) ' t h multiple. Certain data requirements must be met for this process. Specifically the data must: contain 4 components (2 source and 2 receiver), be 2 dimensional (fine source), be regularly sampled in both shot and geophone, contain all near offsets, have surface waves (groundroll) removed, display no variation in shot and receiver amplitude, and contain a known source wavelet. Information about the near surface is also needed. Other than the requirement for multicomponent data, these are the same requirements for removal of marine free surface multiples. Significant effort has been expended meeting the data requirements for multiple removal in the marine case using a combination of acquisition and processing technology. These efforts have contributed to the successful application of this method to real data and underline the importance of being able to remove multiples without needing to know the subsurface. Since the data requirements for land are very similar to the marine case, it is expected that similar efforts are needed to make this technology viable. Fortunately, techniques such as wavelet estimation are portable to the land case. In this chapter, I have addressed the requirements for multicomponent data and a known near surface. Synthetic data tests CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 82 indicated that the method is robust given data with less than four components and errors in near surface information. It is interesting that the data can be placed in a suitable form for multiple removal by applying a series of pre- and postprocessing steps many of which are stochastic in nature. Thus, even though the theory behind multiple suppression is entirely deterministic, the success of the method is dependent on both deterministic and stochastic methodologies. Chapter 4 Multiple Removal for Ocean Bottom Data 4.1 Introduction In this chapter, I explore the application of the inverse scattering series for the removal of interface related multiples from multicomponent seismic data recorded on the ocean bottom. T h e idea of measuring seismic arrivals on the ocean bottom is not new. Numerous experiments have been made to measure multicomponent wave displacements by affixing instruments to the ocean bottom. This technique has been used extensively i n crustal studies to derive P and S wave velocity models from the analysis of refracted arrivals (Mjedle, 1992; Mjedle et al., 1995). These velocity models are used to infer rock properties through a/f3 ratios. T h e course sampling of these experiments has made them unsuitable for more detailed exploration purposes due to the lower spatial (and temporal) resolutions. Recently, acquisition systems have been designed specifically for exploration with smaller receiver spacings, broader bandwidth and more robust geophones. M a n y systems incorporate both vertical velocity and pressure measurements into a cable mounted to the sea bottom (Barr and Sanders, 1989). Later designs measure both vertical and horizon- tal velocity components (Berg et al., 1996). One of the more challenging aspects to ocean bottom recording is proper placement and coupling of the receivers to the ocean bottom. 83 CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 84 Gradually, better instrument engineering has made ocean bottom technology viable though still expensive. The reasons for using ocean bottom seismometers are basically threefold. First, existing production facilities often make it difficult to operate conventional streamers towed behind a single vessel. W i t h ocean bottom cable ( O B C ) data, the cable is placed on the bottom using a recording vessel. A shooting vessel then sails the shot lines which can be much closer to existing infrastructure than with a conventional streamer. W h e n the shot lines are completed, the recording cable is moved to a new location and shooting begins again. T h e second reason for performing O B C measurements is to obtain more complete information about the wavefields propagating in the earth. B y measuring both P and S wavefields, a more complete picture of the earth can be inferred. T h i r d , O B C data has the potential for improved reservoir characterization. B y having permanent receivers i n place over a producing field, different surveys can be acquired over time as the field is being depleted. This time lapse information can be used to help characterize the movement of fluids through the reservoir which in turn can be used to develop an optimal production strategy. In Chapter 3, an inverse subseries was used for the removal of elastic free surface multiples. T h e terms in this subseries were selected based on the reference Green's function which is responsible for the creation of multiples in the forward scattering process. To apply this same principle to ocean bottom data, a few modifications are needed. First, the reference medium should agree with the actual medium at and above the sources and receivers. Since the receivers are assumed to be affixed to the bottom, the reference medium is modelled as a water layer over an elastic half space which has the same properties as the actual material at the ocean bottom. I assume that the ocean bottom is an elastic medium, although, in many cases, the shear modulus of a muddy bottom may range from low to negligible. In the next section, I construct the Green's functions for this reference medium. 4.2 (P?S) Green's function for ocean bottom data Consider a water layer bounded by a free surface on top and an elastic halfspace on the bottom. A n explosive source is located in the water at some depth z . s I choose an explosive CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 85 source because although multicomponent sources for marine data are possible i n principle, they are a long way from being a practical reality. Multicomponent receivers are located at a depth z g which can be anywhere at or below the water bottom. W h e n the receivers are on the water bottom, I assume that they are coupled to the elastic medium. Using this model, consider a downward going plane P wave incident on the water bottom. W h e n this wave hits the bottom, it transmits a P wave into the solid pp ikx gii/i(z-u^-z,,) iu {z -z ) e g e 2 g (^ wb where PP is the P-P water to elastic solid transmission coefficient, and V\ and v 2 ^ are the vertical P wavenumbers in the water and elastic solid respectively. T h e incident plane wave also transmits a converted S wave p£ ikx iui{z -z ) ir) (zg-z ) e ge mb s e 2 ^ wb 2) where PS is the P-S water to elastic solid transmission coefficient and r) is the vertical S 2 wavenumber i n the elastic solid. T h e reflection and transmission coefficients refer to the water-solid interface and follow the nomenclature of A k i and Richards (1980). Specifically, the incident phase is on the left, the scattered phase is on the right, a grave accent denotes v a downgoing wave, and an acute accent 'denotes an upgoing wave. Hence, PP and PP are the P-P solid to water and water to solid transmission coefficients respectively. Including the effects of the free surface, the water column will create a series of multiples which will each transmit into the elastic solid. For P waves, the series is pp ikxg & iv\(z -z ) iv2(z -z ) e wb p pgikxg s e g ^ wb _ pp ivi2z £ _|_ p p ivi4z 2 wb e wb _j_ ^ gii/i(z (4.3) (1 + ppe i ->»>>) iv 2z T h e source ghost creates an additional wavefield which is described by an 'image' plane wave source. This gives pp ikxg^ iui(z -z ) e e wb s _ iui(z +z,)^ iu2(z -z ) e (1 + PPe^ ^") 2 wb e g wb ' T h e corresponding transmitted S wave will have PS instead of PP. ^'^ 4 4 CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 86 A s discussed i n Section 3.3, the Weyl integral can be used to construct 2-D wavefields 2ifi from plane wave solutions. Dividing the plane wave solutions by and inverse Fourier transforming over k gives the P and S Green's functions G p(x ,z \x ,z 0 g g s s G (x ,z \x ,z 0S g g s 3 r PPe^ ^) < z ;u) = r - ^ y^) < z -u>) = - ^ J-oo ^(1 + f wb m e 9 e wbe g dk 2 PS wb J-oo ik(x -x,) i^ iM* -^) ( ^ ' ) n ^ . . - . . l ^ M e M . , - ^ , + PPe******) These two expressions represent the P and S reference wavefields measured i n the solid when the source is i n the water. Next, I develop the Green's functions for sources i n the solid. Consider a line source of P waves in the solid and a measurement point also in the solid. I do not consider measurement points in the liquid since both the actual measurements and the scattering perturbation exist in the solid. For this experiment, three basic waves will be measured: one which propagates directly from the source to receiver, a reflection from the bottom of the water column, and a series of downgoing events which represents the reverberations i n the water column. Summing the P waves for these events and using the Weyl integral gives e oo ^2 *(».-».)( ^l*.-*.l e f / il/l2z where G d PP Zwb i\v2(z +z -2z )+vi 2z \\ c <l i wb L wb ' k jfc (4.6) /roc i — ^OPP e \ \ PPPP (1 + P Pe ™ ) /-^d pp iM^-2 ) + V V ^OPP T represents propagation directly from the source to receiver and G^pp represents the events associated with the water column: the reflections from the bottom and top of the water column as well as the reverberations within the water transmitted back down into the solid. For a line source of P waves i n the solid, the measured S waves are -oo oo 1 / / v v ik{z - )tp§ irn(z -z ) iv {z,-z ) PPPS "2 (1 + PPe ^ ^) / e oo GO S P 'WC g e X3 iv 3 wb e 2 wh V 2l 2z _ i[V2(z +Zs-2z )+v 2z , ] J } e g mb 1 71 b S d ; (4-7) CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM 87 DATA where again the superscript wc denotes contributions from water column events. For an S wave source in the solid, similar equations exist CO -I , / SPPP i\v2{z<,+zs-1z , )+v 2z } c v h 1 wh (1 -f pp i"l2z-wb) (4.8) dk e OPS and CO / I , *(..-.)( *al*.-«.l + 4$ M».+..-2«.0 PPS +z - 2 z ) 2 z ] ^ -co SV2 e c Zl c V } s wb (4.9) wb (1 -I- PPe'"i »») 2z Gd i nwc OSS T o s s • u Combining the two experiments into columns of a matrix gives the Green's function G 0 = OPP I ^OPP r i /Owe GnOSS .<?.<? + T ^^O S S ^OSP GQP P 0 OPS d /Owe 0 Go GO P P wc + 55 (4.10) /mwc ^ O P S /Owe /Owe = G + Gl d 0 ^OSP ^ O S S or i n terms of operators 0 — /-<d I ^ I W 0 (4.11) wc • W i t h this form, a portion of the Green's operator can be identified which is responsible for the creation of events related to the water column - namely G Q . A s in free surface c multiple removal, this will be used to identify an inverse subseries which removes water column multiples. The Green's operator matrix can also be written i n terms of displacements using g = n-To G n 0 1 0 (4.12) CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 88 however, since the propagation and scattering of P and S waves is more intuitive, I will be using the scattering equations written in (P, S) potentials. T h e water bottom reflection coefficients are derived in Appendix B . Although many authors have previously derived elastic reflection coefficients for various situations ( A k i and Richards, 1980; Ewing et al., 1957; deHoop and van der Hijden, 1985), I have chosen to rederive them to avoid both errors i n transcribing notation and errors i n existing publications. 4.3 Water column multiple removal subseries T h e scattering equations for ocean bottom seismic data are essentially the same as for land data. T h e data are described by a forward scattering series in terms of the reference wavefield and a perturbation written in (P, S) components. displacements (u ,u ) x z T h e measured data are generally and, as i n the land case, the multiple removal can be formulated for displacements. Since the data already contain a natural P source, the data are transformed into up and down going P and S waves prior to multiple removal. Section 4.4 shows that this transformation also eliminates certain types of water column reverberations. Following the arguments presented in Chapter 3, the modelling of ocean bottom data can be represented by T> = d V Go + GoV GoV Go + • • • • 0 (4.13) This equation assumes or implies that the sources are multicomponent. I retain this impractical assumption for the purposes of theoretical development. T h e strategy is to formulate multiple removal for full four component data and then use the two component data in the algorithm. T w o component data must be preprocessed to fit the theoretical input data as close as possible. This approach proved effective for land data with missing components i n the previous chapter. Transforming equation (4.13) to (P, S) gives D = G V G + G V G V G 0 0 0 0 0 + ... . (4.14) CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 89 The inverse series can be formulated as i n infinite series of terms i n the data V = + + V< > + . . . (4.15) 2 where D = G V 0 0 = G V 0 0 = GoV^Go + G 0 V ( 2 ) G 0 ( 2 ) V W G O G 0 ( 1 ) G (4.16) 0 + GoV^GoV^Go + GoV^GoV^Go+ (4.17) G O V ^ G O V W G O V W G O (4.18) In the previous section, I showed that the reference Green's operator can be split into two terms. E a c h of these terms have a specific role to play i n the modelling of data in the forward series. T h e term G w c creates events due to the presence of the water column. T h e remaining term, G Q , performs the tasks of data construction exclusive of creating water column multiples. Now, since G Q creates water column multiples in the forward series, it c must be responsible for the removal of these multiples i n the inverse series. Based on this idea, I select terms in the inverse series which sandwich G Q to form the following multiple c removal subseries V ' = V W + V( »' + V( »' + . . . 2 3 (4.19) where V ( 1 ) = GQ ^ G Q 1 y( Y = - V ^ ' G ^ V M - \rWG™VW 3 = V ^ G ^ V ^ G ^ V W - V W G ^ V ^ G ^ V M (4.20) CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM 90 DATA A recursive relation can be derived from these equations which relates the n ' t h term in the subseries to the previous term V W where the first term, (4.21) = - V ^ ' G ^ V ^ acts as a kernel to obtain all the higher order terms. This new inverse subseries can be transformed into data by sandwiching V between two directly propagating Green's functions. In effect, this is a projection of V ' back onto data space. As discussed i n Section 3.4, G Q V ' G O creates data since v is more like V W than V . T h e new series is a sum of data terms D = G V'G D = L V ) + D'( ) + D ' W + . . . . D 1 0 (4.22) 2 A s i n the land case, a recursive relation can be used to obtain all of the terms in this series. T h e first term L V \ are the original data decomposed into up and downgoing P and S x waves v i a ftx>]vl f>W = (see Section 4.4 for a more complete description of this process). (4.23) T h e n ' t h term i n the multiple removal series is given by D'W = -D< >(G 1 O L )- G 1 0 R E (G 0 1 )- D( - >' 1 n 1 (4.24) = -LV"- )'Jr>w. 1 or more explicitly D'^XK^) D'£>{k X,») a D'&\k ,k.,«) g D'£\k X,u) -i J (k,uj) J (k,u>) g PP J (k,uj) SP PS J s(k,uj) S D'£r \KKu) 'DfrVfaM l \DP (k;k„u) P DW(k,k.,«>) D$l(k,k.,u>) dk. DW(k,k.,u>) See Appendix D for the derivation of this equation and the elements of J . (4.25) data contain a single source wavelet A(u>), the inverse wavelet B[u>) = A(ix>)~ W h e n the 1 must be incorporated into the multiple removal terms. T h e multiple free data are then given by D ' = D « + D ' ^ B H + D'^\B(UJ)) 2 + .... (4.26) CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM 91 DATA A s i n the free surface case, when the earth model is one dimensional, the above integral becomes a multiplication i n the (k ,oj) g domain and only a single shot gather is required to calculate the multiple removal terms. E a c h term i n the multiple eHmination subseries has an identifiable role to play in the removal of water column multiples. T h e first term, are the (P, S) decomposed data with the water column effects at the source and receiver removed. This is analogous to removing the source and receiver ghosts from conventional marine data except that the ghosts i n this case are a train of elastic reflected and converted waves. In other words, the objective is to remove the impulse response of the water layer from the source and receiver. T h e second term, T)'( \ 2 represents a modelling of events which propagate up from the scattering region, interact once with the water column, and then propagate back down towards the scatterers. Similar to deghosting, this procedure is more involved than the conventional marine or land case. Instead of a single reflected or converted event coming from the free surface, the water column emits a train of events which are its impulse response. In this way, a train of water column multiples are modelled. W h e n D ' ( ) is added 2 to D ^ ) , it cancels all first order water column multiples in the data. 1 For water bottom data, first order refers to the number of times a wave interacts with the water column as a whole. T h e n ' t h term i n the multiple removal series represents a modelling of events which undergo (n — 1) interactions with the water column. B y summing all these terms order by order, successive orders of water column multiples are removed. Equations (4.23) and (4.26) show that the only pieces of information needed to remove water column multiples are the data, the reference medium parameters, and the source wavelet. This means that the elastic properties of the ocean bottom must be known as well as the thickness and velocity of the water layer. In addition, the water layer is assumed to be flat. T h e reference medium parameters must be either estimated from the data or obtained a priori. components. A n obvious problem lies with the fact that the data have only two T h e solution to this is to use the available data components and set the unknown components to zero. synthetic examples. T h e results of this approach will be demonstrated i n the CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 4.4 92 (P,S) decomposition and up-down going wavefield separation T h e first two steps i n preparation for removing the unwanted multiples are to perform a decomposition from (x, z) to (P, S) coordinates and then remove the source and/or receiver 'ghosts'. This last process removes the effect of downgoing waves at the receiver and the upgoing waves at the source. A s i n the free surface case, the two tasks can be performed at the same time. Unlike the free surface case, the ocean bottom experiment contains only one source component; hence a full source separation into downgoing P and S waves is not possible. Instead, the source can be decomposed into downgoing P waves at the ocean bottom where they are translated into both P and S waves propagating down into the solid. A t the receiver side, (x,z) measurements of up and downgoing waves can be completely transformed into upgoing P and S waves because there are two independently measured components. This process has also be described for multicomponent velocity measurements by Orsen et al. (1996). B o t h source and receiver water bottom decomposition are schematically shown i n Figure 4 . 1 . Theoretically, this procedure requires a knowledge of the P and S wave properties of the water and ocean bottom as well as the source and receiver depths. A s I will show i n section 4 . 5 , the decomposition is relatively insensitive to errors in the elastic properties of the reference medium. Conversely, it is particularly sensitive to the receiver depths since the reverberations within the water column are estimated and removed. T h e first step in formulating the decomposition is to write the reference Green's operators in the form =nGn g T 0 = n (G T 0 D + &™)n. (4.27) T h e data can then expressed as x> = QoVWQo = n (G T D = n {6 + G )nv( >n (G^ T d wc 1 + Gjnv^CGo + G ^ n 1 T + G^)II (4.28) CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA original data Figure 4.1: 93 after source decomposition Effect of water bottom source and receiver decomposition. T h e top diagram shows the effect of the partial source decomposition. T h e middle shows the effect of the receiver decomposition on a scattered P wave. T h e bottom digram shows the receiver decomposition for a scattered S wave. where V ^ ) is in (P, S) components. Factoring out G Q gives 1 T> = n(i + cnGIl)- ) r 1 [GJJV^G;}] (i + (GjJ-^nn (4.29) where the quantity in square brackets represents D : the data decomposed into downgoing CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 94 P and S waves at the source and upgoing P and S waves at the receiver. B y inverting the operators sandwiching D , this gives the desired transformation D = where = II (l + ( G ^ ) " ^ ) r 1 1 and N (4.30) N T>M g f l s = (1 + G^G*)- )- !!. 1 As in the free surface case, the decomposition in (k ,u) 1 is a simple matrix multiplication g of the data by the receiver decomposition function n -k ^2 2 where Z = exp(iuiz b), w _ 9 2 (4.31) 2 2 2 2rj K 2 n -k 2 2 pi 1-Z 2p vi 1+Z Plkg 1-Z K9 +i 2p v\m 1+Z kg 2 2 2 the vertical wavenumbers are related to k, g and the subscripts 1 and 2 refer to quantities i n the water and solid respectively. T h e full derivation of (4.31) is given in Appendix C . A n interesting result can be obtained by splitting (4.31) into the two matrices N = 9 2u "9 or V2 + 2 -hi 2772 0 Pi - 0 P2 2 1 2 q 2 1-Z 1+Z 1-Z 1+Z 2 2 u Plk 2p v\t) 2 N f s + N w c (4.32) 2 T h e first matrix is exactly the decomposition matrix for an elastic free surface (see Appendix C ) . T h e second matrix exists as a result of the water layer. Note that it uses the vertical measurement component only. This is not surprising since the water layer should not affect the horizontal component of the displacements in the solid. T h e decomposition matrix for a water/solid interface with no acoustic free surface can be found by setting Z = 0. This effectively makes the free surface go to infinity giving 12" -i(3 2 2u 2 v2 2rj 2 k + 9 ' 2p ui Plkg 2 (4.33) 2p vir) 2 2 Comparable relations have been derived by Donati and Stewart (1996) and also by A m u n d sen and Reitan (1995a, 1995b) in a comprehensive analysis of ocean bottom reflection and transmission. CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM 95 DATA A problem arises i n the source decomposition since the data are due to only a P source in the water column. Using equations (3.5) and (4.30), the displacement data can be written as the vector d = I ^ D I ^ f = N g (4-34) D(1+ ( G ^ ) - G D ^ 0] . 1 r where f is the body force for an explosion in water. It is clear that the source side operator for this equation cannot be inverted since it is a vector. It is possible, however, to eliminate some downgoing energy from the source side by removing the source ghosts and the multiple reflections within the water column. This latter process can be thought of as deconvolving the P wave impulse response of the water layer from the source. Mathematically, this is done by multiplying the data i n (k ,k ,u)) g space by the operator (1 + PPe ' ) which ll/1 2Zwb s is a type of Backus filter i n f-k space (Backus, 1959). T h e source deghosting is performed by the operator 1/ sm(v ,z ) x s where v ls = ((uj/ax) 2 —k)/. 2 1 2 In summary, the steps i n the source decomposition are: 1) Fourier transform the data with respect to the shot domain; 2) apply the operator (1 + pp ^ ™ ) iv 2z e b (4.35) and 3) inverse Fourier transform back to space. 4.4.1 Removing the reference wavefield A t this point, it is convenient to discuss the removal of the reference wavefield which is a necessary first step before decomposition. Referring to equation (4.5), the reference wavefield is -P Go(x ,z \x ,z g g 3 s < Z ]Ul) wb = J Pe ( 9~ -»>b) iu2 z -P Se ( ~ iv2 Za {"Its) z ^ Zwb + PPe -">) iui2z (4.36) CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA T h e total wavefield contains the reference wavefield plus the scattered field. 96 W h e n the source decomposition is applied to the total wavefield, the reference wavefield becomes (4-37) Other than refracted waves at farther offsets, this wave is generally the first to arrive at the ocean bottom receivers. Consequently, the remaining portion of the reference wavefield can be removed by muting. T h e result is the scattered field with the source decomposition already applied. A synthetic example of this is shown in Section 4.5.2. 4.5 Synthetic data examples In this section, I test the removal of water column multiples using synthetic finite difference data. Like the free surface tests, the data were generated using an elastic finite difference code and simple models. T h e codes for the decomposition and multiple removal were written as modules for the S U processing package. 4.5.1 Model 5 T h e first model, shown in Figure 4.2, consists of a water layer in contact with an elastic layer on top of an elastic half space. T h e impedance contrast between the layers is high to generate many orders of multiples. T h e source is an explosive source located 6 m below the acoustic free surface. T h e receivers are multicomponent velocity geophones affixed to the ocean bottom and are spaced 10 m apart. Each shot is located in the center of a spread of 129 ocean bottom receivers; hence the maximum offset is 640 m . T h e sample interval is 2 ms and the recording time is .512 s. A n arbitrary minimum phase wavelet is used as a source wavelet. T h e data for this model are more complex than for the land examples. To keep track of the different events, Figure 4.3 illustrates the major primaries and multiples that I will be concerned with. In this diagram, the source reverberations have already been removed. Lower case letters indicate propagation through the solid and upper case through the fluid. CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA Acoustic Free Surface # z=6 z m explosive source a x — 1500 97 =0 m/s p! = 1000 k g / m 3 z Water — multicomponent receivers z=8u m a = 2500 m / s & = 1200 m / s p = 1200 k g / m 2 Elastic solid 2 3 =180 m a = 1500 m / s 03 = 900 m/s 3 Elastic solid / 9 = 1000 3 Figure 4.2: kg/m 3 Ocean bottom data Model 5. T h e multiples as a whole are termed water column multiples. W i t h i n this class, I distinguish between multiples that reflect from the bottom and the top of the water column. Events which reflect from the bottom of the water column are termed water bottom multiples. Events which reflect from the top of the water column are termed water top multiples. Referring to Figure 4.3, Ppp ter bottom multiples, and PppPP and PppPPPP and Pps PppPPpp are primaries, are PppPPps Ppppp and Pppps are water top multiples. first order wa- Events such as are downgoing waves at the receiver and are removed in the (P, S) decomposition. As mentioned in Section 4.3, each iteration of water bottom demultiple removes a different order of water column multiples. For example, the second term D ( ) removes 2 events which have one interaction with the water column. This includes the reflection from the water bottom and a train of events which reflect inside the water column. These latter events reflect many times from the acoustic free surface but each transmits only once from the water to the solid. Hence, the order of water column multiples is the number of times CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 98 an event interacts with the water column between interactions with the scatterers at depth. Figure 4.3: Schematic diagram of water column events with source reverberation removed. Lower case letters denote propagation through the solid and upper case through the water. T h e first step is to decompose the data into upgoing P and S waves at the receivers and downgoing P and S waves at the source. A s discussed before, the source decomposition is only partial since the source is single component. T h e left panel in Figure 4.4 shows the vertical component of the scattered field. T h e reference wavefield is removed by applying the source decomposition operator and then muting out the direct wave. T h e panel in the right shows the P portion of the data after decomposition. It is clear that downgoing waves have been severely suppressed from the data. These waves are not removed because stabilizing the decomposition operators introduces small imperfections. Figure 4.5 shows the horizontal component of the scattered field and the corresponding S portion of the decomposed field. T h e events are labelled according to their zero offset arrival times. Notice that the Ppp event has been fully separated from the S component. Since all of the S events originate from a converted P wave at the ocean bottom, their amplitudes are zero at zero offset. T h e first order water bottom multiple is very faint in the S panels; however, the water top multiple is still prevalent. T h e effect of multiple removal for P data is shown i n Figure 4.6. Here, I show the MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 99 CHAPTER 4. Figure 4.4: Decomposition of Water Bottom Model 5. T h e panel on the left contains the vertical component of the scattered wavefield. T h e panel on the right is the P portion of the data after (P, S) decomposition and up-down going wave separation. multiple contaminated data and the terms which remove first and second order water column multiples. Note that the second order multiples are weak in amplitude as compared to the first order multiples. T h e results of multiple removal show that the first and second order multiples are well suppressed. Certain events are not removed due to missing data components i.e., missing subevents. T h e effect of multiple removal on the S data is shown in Figure 4.7. T h e removal of S events is less dramatic than for P . This is expected since the data have missing components. This is similar to the land example in Figure 3.15 where the missing horizontal source component gave good results for P and modest results for S. Nevertheless, when compared against the original horizontal displacement data, the data after multiple suppression are definitely better. CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 100 Figure 4.5: Decomposition of water bottom Model 5. T h e panel on the left contains the horizontal component of the scattered wavefield. T h e panel on the right is the 5" portion of the data after (P, S) decomposition and up-down going wave separation. Effect of errors in water bottom parameters Here I discuss the sensitivity of (P, S) decomposition and multiple removal to errors i n estimating t h e o c e a n bottom parameters. In general, these parameters must be estimated from the data or obtained from some other source. unavoidable. bottom velocities of a 2 2 errors are Consequently, the question that needs to be addressed is how sensitive is the method to these errors. a A s i n any measurement, In Figure 4.8, I show the effect on P data of using water = 2200 m / s and ft = 1000 m / s where the true velocities are = 2500 m / s and ft = 1200 m / s . In the decomposition process, the downgoing events are attenuated but less so than before. Despite the degraded decomposition, the results after multiple suppression are good. A s in the land case, an altered source wavelet is used to compensate for the errors i n material velocity. Although not shown here, the results CHAPTER 4. MULTIPLE REMOVAL (i) D D Figure 4.6: FOR OCEAN '(2) D BOTTOM '(3) 101 DATA Sum Removal of first and second order water column multiples for Model 5. From left to right the data panels contain: 1) P portion of the decomposed data, 2) D^ 2) removes first order water column multiples, 3) D^ , 3) which and 4) sum of panels 1 through 3. for the S component are similar though less dramatic. Tests using various cases indicate that altering the wavelet compensates for either under or overestimating the water bottom velocities. For the reasons discussed in Section 3.7, this compensation is most effective for near offsets. In addition to velocity errors, the method is also sensitive to errors in the water bottom density. One area where the method is very sensitive is in the depth of the water bottom. Since both the decomposition and demultiple processes use the water depth to effectively predict and then subtract certain events, errors in water depth have deleterious effects. Hence, the water depth must be known to a high degree and errors due to topography will be important. A possible extension to the method is to use a reference medium which incorporates variable water depth. CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 102 Figure 4.7: Removal of first and second order water column multiples for Model 5 S data. From left to right the data panels contain: 1) horizontal component of the displacement, 2) S portion of the decomposed data, and 3) S data after removing first and second order water column multiples. In the land case, I showed that errors in decomposition can be avoided by formulating multiple removal in displacements. In principle, a similar procedure could be performed for water bottom demultiple. However, since the decomposition attenuates certain types of downgoing events even when velocity errors occur, it is preferable to have decomposed dat cl CIS Ci final output. 4.5.2 Model 6 In this section, I demonstrate the effectiveness of ocean bottom demultiple using data from a two dimensional earth model. T o do this, I use Model 5 with a dipping lower reflector. A s before, the water layer is flat, the sources are at 5 m depth in the water, and the receivers CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 103 Figure 4.8: Effect of velocity errors on the removal of water column multiples. T h e vertical component of the displacements are shown on the left, the P decomposed data are shown in the middle panel, and the data after two passes of demultiple are shown on the right. B o t h decomposition and multiple removal were performed using velocities of a m / s and ft = 1000 m / s . T h e true velocities are a 2 2 = 2200 = 2500 m / s and ft = 1200 m / s . A n altered wavelet was used to compensate for the velocity errors. are affixed to the ocean bottom. T h e source and receiver spacing is every 12 m . T h e source and receivers locations are referenced by horizontal surface coordinates where one side of the survey is at 62 m and the other at 818 m . Each shot is measured by a spread of 129 receivers which are laid out split spread (shot in the middle) where possible. Since receivers cannot go beyond the edge of the survey, shots at the beginning and end of the survey roll into the receiver spread. T h e receivers then roll along with the shot i n a split spread geometry. Since Model 6 is two dimensional, many shots and receivers are required to remove multiples. T h e data are arranged in a cube 64 shots by 64 receivers by 512 time samples. M u c h of the data cube is zero padding, which is useful when working with Fourier CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 104 transformed data. T h e time sampling interval is 2 ms and the recording time is .512 s. Acoustic Free Surface ^ z=6 z=0 m explosive source a i = 1500 m / s p = 1000 k g / m x 3 z Water V V V multicomponent receivers Elastic solid z=80 a 0 p 2 2 2 = 2500 m / s = 1200 m / s = 1200 k g / m m 3 z = 210 m at shot 506 a = 1500 m / s 03 = 900 m / s 3 Elastic solid p 3 = 1000 kg/m 3 Figure 4.9: Ocean bottom data Model 6. T h e dip of the lower reflector is 10 degrees. In Figures 4.10 and 4 . 1 1 , I show a portion of the P and S data taken through the various processing steps. T h e displayed portion of the data is common shot gather 506 from the middle of the survey. T h e first data panel contains the vertical component of the total displacement field. T h e second panel is the same data with the reference wavefield removed. This is done by applying the source decomposition operator discussed in section 4.4. T o review briefly, the source decomposition involves standard marine source deghosting and then the application of a type of Backus operator which removes the reverberations in the water column. This transforms the reference wavefield into a single downgoing wave recorded on the water bottom. Since this is the first wave to arrive at the receivers, it can be removed using a mute. T h e result is the second data panel in Figure 4.10. T h e source decomposition must be applied to the total field prior to muting so that all of the events in the reference wavefield are removed. For example, on the far left side of the first panel, the event at about .19 s is removed in the source decomposition. CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA The next stage is to apply the receiver (P, S) decomposition operator. 105 Recall that this transforms the data from displacements into upgoing P and S waves at the receiver. T h e value of the decomposition process is evident in the third panel i n Figure 4.10. the PppPP event, the first order Ppp B y removing water bottom multiple can be seen more clearly. In addition, the first order water top multiple PPppPPpp is more continuous. Hence, the decomposition removes some of the unwanted multiples. T h e decomposed S data in Figure 4.11 also show the certain events are attenuated by the decomposition. Note also the Ppp primary has been removed from the S panel. The decomposed data are then used as input into the multiple removal algorithm. After one pass of multiple removal, the data for shot 506 are shown in the fourth panel of both figures. Clearly, the first order water top and bottom multiples have been suppressed or effectively removed from the P data. T h e removal of events from the S data is also apparent although less dramatic than for the P data. This is a consequence of the lack of proper S source data. Nevertheless, even with missing data components, the multiples are significantly suppressed from both components and, best of all, no information about the earth below the water bottom was required. This example demonstrates that the method is effective when the earth is not one-dimensional which is the pitfall of many conventional methods. CHAPTER 4. MULTIPLE T o t a l field Figure 4.10: REMOVAL Scattered FOR OCEAN BOTTOM DATA 106 field Removal of first and second order water column multiples for Model 6 P data, shot 506. From left to right the data panels contain: 1) vertical component of total displacement field, 2) vertical component of scattered displacement field, 3) P portion of (.P, S) decomposed data, 4) P data after one pass of water column demultiple. CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 107 o.i 0.4 0.5 V X Total field V X Scattered field DP + D'P Figure 4.11: Removal of first and second order water column multiples for Model 6 S data, shot 506. From left to right the data panels contain: 1) horizontal component of total displacement field, 2) horizontal component of scattered displacement field, 3) S portion of (P, S) decomposed data, 4) S data after one pass of water column demultiple. Chapter 5 Elastic internal multiple suppression 5.1 Introduction Free surface multiples are generally more prevalent than internal multiples because the free surface is a nearly perfect reflector whereas reflectors below the free surface are weaker. Nevertheless, internal or peg-leg multiples pose a significant hazard for interpreters as they attempt to squeeze more information from seismic data. Conventional methods suffer from the same limitations in attenuating internal multiples as they do in attenuating free surface multiples. Their effectiveness can be further compromised because there is often very little differential moveout between primaries and short period internal multiples. Verschuur et al. (1991) have adapted their surface multiple removal scheme to remove internal multiples. Their approach is to first remove all the free surface multiples, then downward continue the data to the next reflector. T h e multiples associated with this new ' d a t u m ' are then removed in analogy with free surface multiple removal. This layer stripping approach requires the velocities of the earth and the strength of the reflectors as a function of depth. While this method has its merits, it is desirable to remove multiples without needing to know the earth. T h e inverse scattering method provides just such an avenue. Recall from Chapter 3 that the inverse scattering series takes measured surface data as input and gives the earth model (the perturbation) as output. It performs all the necessary tasks of inversion using just the information derived from the data. T h e task of 108 CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 109 free surface multiple removal can be formulated as a completely separable subseries of the full inverse series. T h e task of internal multiple removal has been investigated by A r a u j o et al. (1994) who used an acoustic background reference medium appropriate for marine data and found a separable subseries which performs internal multiple suppression as opposed to removal. Fortunately, this suppression is highly effective (up to 90 %) for realistic earth models. This subseries is always convergent and suppresses all orders of internal multiples independent of the subsurface reflectors which create the multiples. This is a remarkable algorithm and the results usually need to be seen to be believed. Almost as surprising is its ability to suppress converted internal multiples. Using synthetic models, Coates et al. (1996) showed that internal multiples with converted phases could be suppressed using the marine data formulation. T h e attenuation of these converted multiples is not as effective as for nonconverted waves; suppression rates of about 20% were observed. In this chapter, I formulate the equations for internal multiple suppression using an elastic background medium which is appropriate for multicomponent data. M y approach will be to cast the development of Araujo (1994) i n terms of the elastic scattering equations used in Chapters 3 and 4. T h e logic will be the same. Throughout this chapter, I may occasionally refer to the marine method as the acoustic algorithm. It should be made clear that this method is not a formulation for acoustic media only. B y acoustic I mean that it uses an acoustic background medium. I will begin with a brief overview of one dimensional marine internal demultiple and then develop the equations for two dimensional elastic wavefields using an elastic reference medium. 5.2 Review of 1-D internal multiple suppression T h e starting point for internal multiples is to consider data which have no free surface effects i.e., free surface multiples have been removed. In this case, the scattering equations can be written using a perturbation and a homogeneous reference medium. Recall that the subseries which performs marine and land free surface multiple removal was selected using a portion of the reference wavefield which is responsible for the creation of free surface events in the forward series. Thus, understanding how events are generated i n the forward CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION series gives a clue as to how they may be removed in the inverse process. 110 For internal multiple suppression, no such mechanism exists in the reference Green's function; hence a more subtle selection criterion is needed. It has been postulated (Weglein et al., 1997) that a particular type of wave scattering is responsible for the creation of internal multiples in the forward scattering series. This type of wave scattering, shown in Figure 5.1, has the same number of interactions as the actual multiple has reflections i n the earth and shares the same geometric relationship between reflections (lower-higher-lower). It was thought that if scattering of this type creates internal multiples in the forward series, perhaps it removes them in the inverse series. This was the motivation for selecting certain terms in the inverse series. This has a connection to the work in Chapter 2 where I showed that the B o r n approximation to the first order internal multiple was contained in the third order scattering term. It is interesting that the internal multiple algorithm predates this analysis of the forward series. Figure 5.1: Relationship between successive interactions for internal multiple suppression. Here I present an overview of the 1-D internal algorithm which provides the motivation for later deriving the expressions for elastic internal demultiple. A s stated, starting with a medium with no free surface and a plane wave source, a B o r n series for the data can be CHAPTER 5. ELASTIC INTERNAL MULTIPLE 111 SUPPRESSION written as D(Z;UJ)= G (z\z';uj)V(z')(f>(z';uj)dz' j 0 J—oo oo / poo G (z\z']u)V(z') / G (z'\z";w)V(z")4iz"-u)dz"dz' + ... J — oo 0 -oo 0 (5.1) (see Chapter 2) where D are the data, G is the reference Green's function, V is the scattering perturbation, <j> is a n incident plane wave, z is the measurement location, and u> is the angular temporal frequency. T h e reference Green's function for a homogeneous background medium with velocity c is (DeSanto, 0 1992) Go(z|*» = where k = 0 2 ^ ' " ' e i f c 0 Z 2 (-) 1 5 2 UJ/C . 0 A n inverse series can be written in the form V = Vx where V n + V + V 2 3 + ... (5.3) is n ' t h order in the data. Substituting this inverse series into the forward series yields a series of expressions in successive orders of the data D = G V <f> 0 x 0 = GoVxGoV^ + G V 4> 0 2 (5.4) 0 = GoVxGoVxGoVxtp + GoVxGoVxp + G V G V t> + G V^ Q 2 0 l( 0 T h e first equation describes the data in terms of V\. Using <p = exp(ifc z') and (5.2) yields 0 D(zwhere Wi(k ) 0 ko) = -i-e^y1(2&0)= = Vi(2k ) and D(z;k ) 0 ^e^'Wxiko) = D(z;u>/c ). 0 0 (5.5) Similar to the forward series, the terms in the inverse series can be interpreted as a series of reference medium propagations and 'interactions' with V\ which is related to the data. T h e third term in the series was selected since it involves three orders of interactions with the data. This term is G V <f> 0 3 = -GoiVxGoV, + V G Vx + VxGoVxGoVx)^ 2 0 = G ( V s i + %2 0 + %z)4>- (5-6) CHAPTER 5. ELASTIC INTERNAL MULTIPLE 112 SUPPRESSION Analysis of these terms by Araujo has shown that V 3 1 and V 32 do not lead to a diagram of the type shown in Figure 5.1. T h e remaining term is G V <f> = -GoViGoViGoVrf = W . 0 33 (5.7) 33 A s a consequence of the absolute value signs in the exponentials of Go, this equation is split up into a series of four terms which have different limits i n the integrals. E a c h of the four terms imply a relationship between the depths of the three scattering interactions with W , contains Vi. One of these terms, the desired lower-higher-lower relationship and was 333 chosen for further analysis. This term can be written as -| flOO W (k ) = —i / Z33 4 f e 0 PZ\ dz e ^V {z ) Q 1 i2h 1 / l J-OO POO dz e- °"V {z ) ak 2 1 / 2 J-OO dz ^ ^V {z ). k 3 x (5.8) 3 J Z 2 Using step functions to describe the integration limits and evaluating the integrals using contour integration, (5.8) becomes W {k ) = ^j^i( o)W (-k )W (k ) 333 k 0 4 1 0 1 + 0 S residues due to poles in V i ( & o ) . (5.9) Kg Using the first part of this expression and (5.5), the following expression was selected as a candidate for internal multiple removal D{f {z; ko) = D(z- k )D(z; -k )D(z; k ) 3 where 0 D^(z; k ) = - -^ ° W {ko). 0 e%h 2 z 0 W h e n added to the data, 333 (5.10) 0 D successfully removed 333 first order internal multiples (hence the superscript IM); however it also introduced extra events which corrupted the final output. T h e continued search for a multiple removal scheme that did not produce unacceptable artifacts led to further analysis of D^. Substituting into (5.10) the expression oo e 'D(z; z')dz' (5.11) ikoZ / 00 and then breaking the limits of the integral into portions leads to a series of expressions which have a geometrical relationship. T h e portion which contains a lower-higher-lower relationship is OO / PZ\ dz D(z; ikaZi ie CO z) x / J —OO POO dz e- ° >D(z; 2 ik z z) / 2 v Z2 dz e ° >D(z; 3 ik z z ). 3 (5.12) CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 113 Hence, instead of V\, the data D are taken through the picture in Figure 5.1. added to the original data, D When suppresses all first order internal multiples regardless of 3 the geometry of the reflectors in the earth. A l l the necessary information is obtained from the data which is used as input to the algorithm. One of the strengths of this method is that it suppresses an entire class of multiples which have a lower-higher-lower geometry. Moreover, since this method is subtractive, primaries are not affected. Further terms were subsequently found which suppress higher orders of internal multiples. For example, the term which attenuates second order internal multiples is CO / fiZl dz *° D(z; ) Zl ie CO Zl / V /"CO dz e- *D(z;z ) / dz e ° > D(z; J Z2 ) / dz e « >D(z; J Zi ikoZ 2 ik 2 — CO dz e- ° D(z; ik z 3 Zi ik A Z4 •CO z) s z 5 z) . 5 (5.13) B y summing the multiple suppresion terms, a series can be written which attenuates all orders of internal multiples D\z- ko) = D(z- ko) + D (z; k ) + D (z; k ) + D (z; k ) + . . . . 3 0 5 0 7 (5.14) Q It turns out that more than the first three terms are rarely needed to adequately suppress all of the internal multiples since higher order multiples tend to be weak in amplitude. 5.3 Theory for elastic internal multiple suppression In this section, I follow the development of Araujo (1994) for two-dimensional earth models to derive an internal multiple suppression method using an elastic background medium. I start with the equations for forward and inverse elastic scattering i n (P, S) components (see Chapter 3) D = GoV^Go 0 = GoV^GoVWCo + G O V ^ G Q 0 = GoVWGoVWGoVWGo + G o V ^ G o V ^ G o + G Q V W G O V ^ G O + G V( 0 3 ) G 0 (5.15) CHAPTER 5. ELASTIC INTERNAL MULTIPLE 114 SUPPRESSION Using subscript notation, the data can be written as (5.16) Da = GoikV^Goij where the subscripts take on values of P or S to indicate P or S wave propagation and scattering. I employ subscript notation because it is easier to keep track of P and S modes of propagation. Numbered subscripts will be used to distinguish between different variables e.g., xi and x. 2 For an isotropic homogeneous background medium, the Green's function i n (P, S) coordinates is diagonal and is given by Goi(z ,z \z.,z.i ) g g — *<--.> *«l*.-.l<tt = —] U e (5.17) c where 6u = v and rj for i — P and S respectively. T h e data are then written as OO -I —^e^te-'^e^i^-'^dkidxidxidzt / -oo 2i0 (5.18) K 2j where depends on three spatial variables. 1 Note that 9 j is tied to the horizontal 2 wavenumber k . Fourier transforming (5.18) over the shot and geophone domains and rec2 ognizing that the spatial integrals are Fourier transforms, I obtain the compact expression Dij(h, z \k , z ; ) = g 2 a U -l- -«i".^ )(fc , k , 6 1 e 1 2 + U *)JL-««" ai c . (5.19) Next, I consider the equation which is third order i n the data and select just the term = GQIV^GomVff. Sandwiching this function between two reference Green's functions and Fourier transforming over shot and geophone yields oo / V£\x , 3 CO x,, G 0 m ( z , z \x , 4 4 5 z ; UJ) V$(X , 6 6 X , Z ) e*** e 6 6 i k ^di (5.20) where the integration symbol with a tilde and f denote integration over all the x and z CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION . 115 spatial variables. Substituting (5.17) into (5.20) gives CO / to Z" 1 — / 27T / _ l OO /»oo -I 0 0 -I 7 ^ ^ ( - - ^ ) 1 00 — — e 2l0 ife4(x4 - X5) e e ^ | ie4 ™ Z 2 |z4 - - ^ l ^ 26| v£\ , 3 Xi Z 4 , Z 4 ) dfc ^V( 5,a;6,^)e^ e :c 4 Z6 ife2a:6 df . (5.21) 4m Because of the absolute value signs in the exponents, this equation must be split into pieces depending on the relative orientation of the depths z , z , and z&. I choose the region 4 2 where z > z and z > z which represents a lower-higher-lower relationship between the 2 4 6 4 scattering depths. Equation (5.21) then becomes CO / 00 1 i»CO i 6 ^V^\h,k ,z ) 1 f°° - /> oc / / e J — 0 0 J — 00 / —e ^-^dk H(z J-oo id to 3 2 3 2 - zjvWfa, h,z ) 4 1 poo 1 JT" / j r ^ e ^ ^ - ^ ^ t f ^ - z e j v j ^ 27T j * _ 2z0 OQ (5.22) 4m where I have carried out the x integrations and H is the Heaviside or step function. T h e Fourier transform of the Heaviside function is - z)= — / H(z 4 2 2TT J^ -dq e- 0 (5.23) i{q-te) where e is a small positive quantity which approaches zero. T h a t this expression gives the Heaviside function can be easily seen using a contour integration where a single pole exists in the upper half q plane at q = ie. Using (5.23) i n (5.22) leads to I poo These expression can be evaluated by performing contour integrations i n the complex q and r planes using a infinite half circle contour which runs through the real axis. For CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 116 each integration, there is a single pole at ie and also poles in the function which are unknown in general. Using the residue theorem, (5.24) gives -i V^ih, k ,6 2 +0 ) = — U / 2j oo ^r- S(k , v / •co 77 CO V^\k„ k ,6 3 k, 6 4 2 + 0 i)-fV^(k , U 3 + fc ,-(0 3 3l 4 6 )) 4m + 6 ) + £ residues due to poles i n V™ . Am 2j (5.25) ^"im Next, I define a modified form of the data Bij{ki,k , 6 2 U + 8 ) = J — e - ^ V ^ i h , k, 9 2j 2 + 8 U which is the original data converted to a plane wave source. )e- 2 i i 6 ^ (5.26) Following Araujo, I neglect the extra residues in (5.25) and use the new data type to obtain jr oo A / +e ) —*-BQ(k ,h e^ ai 2j n (5.27) oo 2 l 0 4 m where I have converted the portion without the residues, l / ( 333 ) , to _g( ). This equation 333 is the 2-D elastic version of (5.10) and the 2-D equation i n Araujo (1994). It should have the same properties as the 1-D and 2-D versions: effectively removes multiples but generates extra events which distort the final data. Following the path taken in the marine formulation, the new data are written as the Fourier transform oo B (k k ,z )e ^' ^dz ij / u 2 i 1 +e (5.28) l •oo which can be split into OO flZ a Bijih, / k, z )e ^ 2 1 i i + e ^dz •oo 1 = / B^h, k, z )e ^ 2 i 1 + e ^dz 1 «/ —oo fiOO + / B {k ,k ,z )e ^' ^dz ij 1 2 1 i +e (5.29) 1 " Za which correspond to z\ < z a substitute for each B^ and z\ > z respectively. a Using this splitting technique, I in (5.27) and choose the interval z\ > z 2 and z 3 > z 2 (i.e., a CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 117 lower-higher-lower relationship). This gives This expression is the elastic version of (5.12) and the 2-D equation derived by Araujo. W h e n added to B^\ki,k ,6u 2 + 0 j), 2 the role of this term is to suppresses all first order internal elastic multiples. T h e only real difference between this equation and the marine formulation is the matrix nature of the data and the wavenumber 77. In fact, the expression for marine data is contained in this equation when all the subscripts are P. T h e remaining terms describe the mixing of different data types to estimate different types of multiples. Note that I have introduced the small parameter a into the integration limits. This ex- cludes contributions from Zi = z = z which distort the primary reflections. In practice, 2 3 a must be greater than the length of the wavelet in the data to avoid primaries interacting with themselves. Thus, a limit is placed on the spacing between primaries and their associated internal multiples. In other words, internal multiples cannot be suppressed between primaries that cannot be resolved. T h e thin bed response is an example of this situation. The same process used to obtain (5.30) was employed to obtain an expression which suppresses second order internal multiples. Sparing the reader the detail, I start with the fifth order term G O V W G O V W G Q V W G O V ^ G o V^Go and invoking a lower-higher-lower- CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 118 higher-lower relationship between the depths of successive interactions, I obtain oo / jjr JP- OO / 77 fz\—a / B ^ h ^ z ^ e ^ ^ ^ e ^ ' ' ^ POO oo ^ "hn Jz +<r % 2 oo ^ 6?6p p y_oo r -B^^Kz^e^+^dz,. (•OO (5.31) 1-D earth models A t first sight, (5.30) may seem like a horrible mess. A few simplifications will make this expression much more readable. T h e first is to consider data from one-dimensional earth models. In this situation, one shot gather is the same as all the others; hence the data vary only i n offset. Mathematically, this means that the only nonzero data values are for Aii = k ; thus 2 the data can be written as B^\ki, k , 0u + 9 j) 2 2 — B^\ki,0u + &ij)5(ki — k ). 2 Substituting this expression into (5.30) eliminates the integrals over horizontal wavenumber. I will also consider data that has been redatumed so that z = g z 3 = 0. T h e resulting expression is z\—a 2 / />oo B^{k^z )e-^ ^dz 6 2 OO I B^ (k z )e ^ ^dz . j u 3 i +e (5.32) 3 J Z -V<T 2 Comparing this to (5.12), it can be seen that this equation describes the one-dimensional internal multiple algorithm for a fixed k\. T h e type of data this algorithm requires is specific but understandable. First, start with data which has been modified from a line source to a plane wave source. This is a simple multiplication by 2iflij i n (fci,u>) space {k\ is the Fourier transform wavenumber over geophone ). From (ki,uj) coordinates, the data are then transformed to (An,fli; + 9\j) coordinates. This is essentially an elastic version of the transform used i n Stolt migration CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 119 (see, for example, Claerbout (1985)). From here, the data are inverse Fourier transformed to (&i,z) which is, in fact, data migrated using a constant reference velocity. T h e variable z is termed pseudo depth and is not the actual depth unless the earth comprises only the reference medium. T h e data in pseudo depth is then input into (5.27) which computes the first order internal multiples. T h e multiples are calculated for each set of (ki,0u + 8\j) by integrating the data over a range of pseudo depths z. T h e resulting estimate of the multiples is then transformed back to (ki,uj) and then to (x ,t). g In reality, since multiple suppression is a subtractive process, the suppression need not be done in (x ,t); g any suitable domain will do. Elastic internal and free surface multiple attenuation share similar data requirements: line sources, four components, full and regularly sampled data, and known source wavelets. T h e only real difference is that the internal demultiple algorithm assumes that the free surface multiples are removed. Hence, internal multiple attenuation follows free surface multiple removal. Another requirement is that the data be transformed to (P, S) coordinates. This can be done using the method outlined i n Chapters 3 and 4. A consideration must be made to the effectiveness of the (P, S) transform when the velocity i n the vicinity of the sources and receivers is i n error. T h e resulting artifacts i n the (P, S) data will have an effect on the estimate of the internal multiples. It is not obvious that equation (5.32) is independent of the reference velocities. To see this, it is helpful to consider only the marine formulation ( P subscripts only). T h e reference medium velocity appears to be imbedded in the transform to vertical wavenumber 8. However, if B\^IM is likewise Fourier transformed to (ki,z), the algorithm describes a series of 1-D dimensional internal multiple calculations for each k\ with the data i n pseudo depth. Since the 1-D algorithm doesn't depend on an earth model, it will generate a proper estimate of the multiples regardless of the velocity used to transform (or migrate) the data to pseudo depth. In other words, using different reference velocities will stretch the pseudo depth axis of the input data. However, the 1-D algorithm will properly estimate the internal multiples regardless of the relationship between the pseudo depth and travel time. For 2-D earth models, equation (5.30) computes internal multiples for a fixed (ki, k , 9u+ 2 &2j)- T h e internal workings of the algorithm are a little more complicated since a series of CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 120 integrations over horizontal wavenumber are required. Synthetic data example 5.3.1 To demonstrate the attenuation of internal elastic multiples, I consider the one-dimensional model shown in Figure 5.2. This model (Model 7) contains two interfaces, no free surface and multicomponent sources and receivers. T h e data will be four component (P, S) in (k, z) space for an incident plane wave and a fixed k. data Hence, multiple suppression will be implemented using (5.32) for a single k. I constructed the synthetic data analytically in (k,k ) z coordinates using the same plane wave construction I discussed in Section 4 . 2 . A n inverse Fourier transform then places the data in (k, z) space where z is pseudo depth. multicomponent multicomponent sources receivers * V a i = 1500 m/s 0x = 900 m/s I 100 m Px = 1000 k g / m 3 a 2 0 2 = 2500m/s = 1200 pi = 2000 k g / m a 0 3 3 200 m m/s 3 = 1500m/s = 900m/s ps = 1000 k g / m 3 Figure 5.2: Internal elastic multiple M o d e l 7 T h e data for Model 7 were calculated for k = —.0074 1 / m and include the first and second order internal multiples. Figure 5.3 shows the data with the various primary reflections labelled; other events are first and second order multiples. In these labels, the letters following the underscore denote subscripts (e.g., B _ P P is the same as Bpp). model, converted internal multiples are more prominent on the BSP and Bps the Bpp and Bss For this traces than traces. Since the elastic formulation is likely to have the most effect on the cross components (i.e., PS, SP), I will look closely at Bsp- CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 121 Using (5.32), I calculated the first term in the multiple attenuation series for Four traces contribute to this term and are labelled (1) —•> (4) in Figure 5.4. are two copies of the data, Bg'p. Z s IM Also shown T h e second of these will have its multiples suppressed (it will be summed with the elements of B^p ) IM for comparison. B p . while the first will remain untouched Predictably, most of the energy comes from the traces which have Bpp and Bg'g sandwiched between other data components. This is simply because these traces contain more energy. In Figure 5.5, I show the same plot with a cumulative sum running from traces 6 through 2. Trace labels such as B3(2:4) denote a sum from B3(2) to B2(4). This figure shows how the four multiple suppression traces interact with each other. In the final sum (trace 2) the internal elastic multiples have been suppressed at rates ranging from -10 to -20 db (about 60 to 90 %). These rates are comparable to the performance of the marine algorithm. In Figure 5.6, I compare the elastic demultiple term to an acoustic result which uses just B^p data. Although the acoustic algorithm does very well in predicting the time of many of the multiples, the amplitudes are very low. Scaling prior to suppression may alleviate these differences but the relative scaling of the multiples is inaccurate. As a final display, all four data components after adding B ^ 5.7 and 5.8. are shown i n Figures I M In all of these data, multiples of all phases are consistently suppressed and primaries are left untouched. 5.4 Discussion and conclusions Elastic internal multiple suppression is an adaptation of the marine method developed by Araujo et al. (1994) except using an elastic background medium. is applicable to land and ocean bottom multiple component data. This formulation C o m m o n to all the inverse scattering demultiple methods, this method suppresses all orders of elastic internal multiples independent of the earth which causes the reflections. are not affected. In addition, primaries T h e data requirements are the same as for free surface elastic multiple removal. Using a synthetic data example, I have demonstrated that the elastic formulation effectively suppresses both reflected and converted internal multiples. T h e suppression can likely be improved by scaling the multiple suppression terms according to a energy CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 122 minimization criteria after multiple suppression. This may be necessary in practice since the wavelet is never known exactly. T h e significance of using an elastic background formulation is that the relative amplitudes of the reflected and converted multiples are more correct than if the marine method is used. Since, in the marine method, these relative amplitudes contain more errors, the scaling procedure will not improve the suppression of both converted and reflected waves and a tradeoff is therefore necessary. This suggests that when there is a strong mixture of both P and S internal multiples, an effective multiple suppression technique would be to record multicomponent elastic measurements either on land or the ocean bottom, remove the free surface or water column multiples, then apply elastic internal multiple suppression. CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION B_ps x 10 150 300 450 600 750 900 pseudo depth z (m) 1050 1200 1350 1050 1200 1350 B_ss 150 300 450 600 750 900 pseudo depth z (m) Figure 5.3: Model 7 internal multiple elastic data for a fixed wavenumber CHAPTER 5. ELASTIC INTERNAL MULTIPLE 124 SUPPRESSION B_sp, multiple suppression terms 0.01 H n 1 1 1 a. a. 1 1 1 i r vi ol 0.005 Bl -0.005 Bl -0.01 B3(l) -0.015 -0.02 -0.025 B3 (2) B3 (3) -0.03 * -0.035 -^r*-— * i u —B3(4) -0.04 150 300 450 600 750 900 pseudo depth z (m) 1050 Figure 5.4: Elastic internal multiple suppression for B^. B p. l s T h e bottom four are the contributions to B p . 3 s IM 1200 1350 T h e top two traces are the data CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 125 B_sp, cumulative sum n 0.01 1 1 r 0.005 -0.005 B1+B3 (1:4) -0.01 B3(l:4) -0.015 -0.02 B3 (2:4) — -0.025 B3 (3:4) -0.03 -0.035 B3(4) -0.04 h 150 300 450 600 750 900 pseudo depth z (m) 1050 1200 Figure 5.5: Elastic internal multiple suppression for f?£p- 1350 A cumulative sum has been applied from traces 6 through 2 in Figure 5.5. Trace 1 is the original data for comparison. CHAPTER 5. ELASTIC INTERNAL MULTIPLE Elastic 'demultiple B3 (1:4), B_sp x 10 150 300 450 600 750 900 pseudo depth z (m) 1050 1200 1350 1500 1200 1350 1500 Acoustic demultiple B3 (1:4), B_sp x 10 2.5 SUPPRESSION -i 1 2 1.5 1 0.5 0 -0.5 -1 -1.5 0 150 300 450 600 750 900 pseudo depth z (m) 1050 Figure 5.6: Elastic (top) vs. acoustic (bottom) internal demultiple term for B CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION B_pp, Elastic internal demultiple 150 300 450 600 750 900 pseudo depth z (m) 1050 1200 1350 B_sp, Elastic internal demultiple x 10 Before After 150 300 450 600 750 900 pseudo depth z (m) 1050 1200 Figure 5.7: Elastic internal multiple suppression for Bp 1350 P and B p L s CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION B_ps, Elastic internal demultiple x 10 1 1 1 — 1 o. o. a c &- GO a. a OH GO GO GO o- I I I Before I I x^^-^ljl^——-^K-^—-*w\Jw^~—Ylr—^ ^^l^v^^ L I , ' A f t e r ^-"-^v*.—* 1 \ y -*AJW''"— -~~V"4W--*' M w —• " _ i 150 300 450 .i 600 750 900 pseudo depth z (m) i l 1050 1200 1350 1200 1350 1 B_ss, Elastic internal demultiple 150 Figure 5.8: 300 450 600 750 900 pseudo depth z (m) 1050 Elastic internal multiple suppression for Bpg and Bg'g. Chapter 6 Conclusions In this thesis, I have investigated various aspects of the forward and inverse scattering series - particularly the application of the inverse series for the attenuation of elastic multiples. A s a preliminary, I presented i n Chapter 2 a mapping between the terms i n the forward scattering series and both primary and multiple reflections for simple seismic problems. T h e mapping was obtained by calculating the terms i n the forward scattering series and comparing them with well known solutions for simple reflection problems. This analysis revealed that primaries are described by all the terms in the forward series, whereas a multiple with n changes in propagation direction is described by the n ' t h and all higher order scattering terms. In addition, I demonstrated that for a specific event, the terms beyond its B o r n approximation change reflection and transmission amplitude and alter propagation velocity. B y quantitatively identifying those terms which build internal multiples in the forward series, this work aids in the understanding of internal multiple suppression using inverse scattering series. A common consideration when talking about any series is its region of convergence. To this end, I have also identified the convergence criteria for simple scattering solutions. For 1-D models, the region of convergence is dependent on the contrast between the actual and the reference medium and on the incident angle; the region of convergence shrinks with increasing angle. This has implications for performing full inversion using the inverse scattering series since the convergence of this series is dependent on the convergence of the 129 CHAPTER 6. CONCLUSIONS 130 forward series. Despite this fundamental limitation, the series which perform free surface multiple removal and internal multiple attenuation, including the work in this thesis, are convergent. This convergence can be said to be dependent on the data being small in. some way. This can be addressed by considering that the input data are scaled by the source wavelet; hence, their amplitudes are relative to the reflection coefficients i n the earth. Because these reflection coefficients are always less than one, the data are also, in a sense, less than one. This is another way of saying that the earth doesn't give back more energy than was initially put i n . In practice, the input data can be scaled to avoid instability and, since the data are limited in time and space, the entire series doesn't need to be computed anyway. In Chapter 3, I developed a 2-D method which removes all orders of free surface elastic multiples from multicomponent land data without requiring information about the subsurface reflectors. Moreover, primaries are not affected. This method is an adaptation of the inverse scattering method developed for marine seismic data (Carvalho et al., 1992). The effectiveness of this algorithm was demonstrated using synthetic seismic data for both one and two-dimensional earth models. Theoretically, the method requires the use of four component data, a known near surface and a known source wavelet. Synthetic data tests indicated that the method is effective using data with less than four components and is robust with respect to errors in estimating the near surface velocities. T h e method also requires data that are regularly sampled and contain near offsets. W h e n the source wavelet is unknown, as it usually is, wavelet estimation techniques used i n the marine case are portable to the land case. Based on synthetic data tests, one of the main conclusions is that when data are recorded on or in an elastic medium, elastic free surface multiple removal using multicomponent data gives better results than the acoustic formulation using single component data. In Chapter 4, I adapted the elastic free surface method to develop a procedure which removes multiples from multicomponent data recorded on the ocean bottom. T o my knowledge, this is the first method developed specifically for ocean bottom data that removes all multiples associated with the top and bottom of the water column without requiring a model of the reflectors below the ocean. A s in the case of land data, primaries are not CHAPTER 6. 131 CONCLUSIONS affected. Synthetic tests for both one and two dimensional earth models demonstrated the effectiveness of this technique. T h e method makes demands on the data type similar to the land case: regular sampling, known wavelet, and near offsets. In addition, the water top and bottom are assumed to be flat and the properties of the water column and ocean bottom material are required. It is assumed that the source is a pressure source in the water column and that multicomponent measurements are made of the displacement (or velocity) at the ocean bottom. Hence, the data need not be four component. T h e tradeoff for using two component data is that certain types of multiples cannot be removed. Tests indicated that the method is robust given errors in estimating the material at the ocean bottom but is sensitive to errors in water velocity and depth. A n information criterion is common to marine, land, and ocean bottom free surface (or ocean column) multiple removal. Each method requires information about the reference medium appropriate to the recording geometry. Table 6 summarizes the information re- quirements for different data types. T h e underlying principle is that while these methods make no assumption about the earth at depth, they rely on being able to correctly describe interactions with an interface (or interfaces) which creates the multiples. For example, since the algorithm for marine data uses a liquid half space reference medium, this method requires the water velocity and the source and receiver depths. A s stated before, all of the free surface multiple removal methods are sensitive to some parameters but robust with respect to others. Common to all these methods are the recorded data as input and knowledge of the source wavelet(s). In Chapter 5, I extended the theory for inverse scattering internal multiple suppression (Araujo et al., 1994) to an elastic background medium. Unlike the free surface and ocean bottom methods, the internal multiple suppression scheme does not depend on the description of an interface which causes the multiples. A l l reflector information is obtained from the input seismic data. In common with the method for marine data, the elastic formulation assumes that the data do not contain free surface or water column multiples. Hence, the work in Chapters 3 and 4 sets the stage for implementing internal elastic multiple suppression. While no subsurface information is required, internal elastic multiple suppression makes CHAPTER 6. CONCLUSIONS 132 data type reference medium required parameters marine water half space source and receiver depths; water velocity land elastic half space source and receiver depths; P and S wave near surface velocities ocean bottom water layer over source, receiver, and water depths; elastic half space water velocity and density; P and S wave velocities and density of water bottom Table 6.1: D a t a type and the corresponding reference medium determines the needed information for free surface multiple removal demands on the data that are similar to those for free surface multiple removal. T h e data must contain near offsets, four components, be regularly sampled, and contain known source wavelets. A n important property of inverse scattering multiple suppression for land and marine data is the ability to suppress multiples without affecting the primaries. I have mentioned this numerous times but the significance of this property cannot be overstated. T h e problem with many multiple suppression techniques, particularly filtering methods, is that while they do a good job suppressing the multiples, they can also adversely affect the primaries. In areas with severe multiple problems, for example, areas with a hard water bottom, filtering techniques are like applying chemotherapy to the data: hopefully the disease dies before the patient does. For this reason, prediction and subtraction methods are desirable. The point to all this is that leaving the primaries intact minimizes a source of error for true amplitude processes such as A V O and inversion which are applied to the data after multiple suppression. Finally, let me reiterate where the work in this thesis fits i n with previous works i n inverse scattering demultiple. T h e works of Carvalho et al. (1992) and Araujo et al. (1994) were designed for marine data where the source and receivers are located i n the water CHAPTER 6. column. CONCLUSIONS 133 These methods use an acoustic background medium; however, it should not be construed that they represent an acoustic approximation to an elastic earth. For example, the marine free surface multiple removal properly removes free surface multiples which contain shear propagation in the elastic subsurface. There is no extension to elastic theory which would improve upon this. For internal multiples, the situation is more complicated. T h e forward model for generating marine data can be described as propagation in an acoustic reference medium and scattering from points within a perturbation model. This forward model must somehow generate both P and S primaries and internal multiples. T h e inverse series can be similarly described as propagation in an acoustic reference medium and scattering from the data. T h e internal multiple attenuation method is based on a subseries of this inverse series which, in principle, outputs an elastic earth model. Hence, the mechanics of elastic internal multiple removal are contained somewhere within the full inverse series. It happens that the subseries of Araujo et al. attenuates both P and S internal multiples, but is more effective for P. It may be possible to find terms in the inverse series for marine data which attenuate both types of multiples equally. This remains an open question. In contrast with the above works, the methods in this thesis are concerned with data that are recorded on or i n an elastic medium. In that case, the appropriate reference medium is elastic. Moreover, I have investigated the removal of multiples from multicomponent data (although the data need not be multicomponent for these procedures to be effective). Multicomponent data have the benefit of providing access to both P and S primaries which contain more information about the subsurface than just P primaries. Coupled with this extended set of primaries are a complicated set of P and S multiples. T h e work in this thesis addresses the removal of these multiples so that the most reliable information can be obtained from multicomponent data. Bibliography Achenbach, J . D . , Gautesen, A . K . , and M c M a k e n , H . , 1982, Ray methods for waves in elastic solids: P i t m a n Publishing Inc, London. A k i , K . , and Richards, P. G . , 1980, Quantitative seismology: W . H . Freeman and C o . , USA. Alverez, G . , and Larner, K . , 1996, Implications of multiple suppression for A V O analysis and C M P - s t a c k e d data: 66th Annual Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1518-1521. Amundsen, L . , and Reitan, A . , 1995a, Decomposition of multicomponent sea-floor data into upgoing and downgoing P- and S-waves: Geophysics, 60, 563-572. 1995b, Extraction of P- and S-waves from the vertical component of the particle velocity at the sea floor: Geophysics, 60, 231-240. Araujo, F . V . , Weglein, A . B . , Carvalho, P. M . , and Stolt, R. H . , 1994, Inverse scattering series for multiple attenuation: an example with surface and internal multiples: A n n u a l Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 64th 1039- 1041. Araujo, F . V . , 1994, Linear and non-linear methods derived from scattering theory: backscattering tomography and internal multiple attenuation: P h . D . thesis, Universidade Federal da Bahia, Brazil, (in Portuguese). 134 135 BIBLIOGRAPHY Backus, M . M . , 1959, Water reverberations - their nature and elimination: Geophysics, 24, 233-261. Barr, F . J . , and Sanders, J . I., 1989, Attenuation of water-column reverberations using pressure and velocity detectors: 59th A n n u a l Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 653-656. Berg, F . J . , Svenning, B . , and M a r t i n , J . , 1996, S U M I C : a new strategic tool for exploration and reservoir mapping: 56th Meeting of the European Association of Exploration Geophysicists, Extended Abstracts, G055. Carvalho, P. M . , Weglein, A . B . , and Stolt, R. H . , 1992, Non-linear inverse scattering for multiple attenuation: Applications to real data, Part I: 62nd A n n u a l Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1093-1095. Carvalho, P. M . , Weglein, A . B . , and Stolt, R. H . , 1994, Wavelet estimation for surface multiple attenuation using a simulated annealing algorithm: 64th A n n u a l Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1481-1484. Carvalho, P. M . , 1992, Free-surface multiple elimination method based on nonlinear inversion of seismic data: P h . D . thesis, Universidade Federal da Bahia, Brazil, (in Portuguese). Claerbout, J . F . , 1985, Imaging the earth's interior: Blackwell scientific publications, Oxford. Clayton, R. W . , 1980, A n inversion method for elastic wavefields: Stanford Exploration Project Report, 24, 81-92. Coates, R. T . , and Weglein, A . B . , 1996, Internal multiple attenuation using inverse scattering: Results from prestack 1 and 2-D acoustic and elastic synthetics: 66th A n n u a l Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1522-1525. deHoop, A . T . , and van der Hijden, J . H . M . T . , 1985, Seismic waves generated by an impulsive point source in a solid/fluid configuration with a plane boundary: Geophysics, 5 0 , 1083-1090. 136 BIBLIOGRAPHY DeSanto, J . A . , 1992, Scalar wave theory: Springer Verlag, Berlin. Donati, M . S., and Stewart, R. R., 1996, P- and S-wave separation at a liquid-solid interface: Journal of Seismic Exploration, 5, 113-127. Dragoset, W . H . , and M a c K a y , S., 1993, Surface multiple attenuation and subsalt imaging: 63rd A n n u a l Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1099-1102. Ewing, W . M . , Jardetsky, W . S., and Press, F . , 1957, Elastic waves in layered media: McGraw-Hill, U S A . Fokkema, J . T . , and V a n den Berg, P. M . , 1990, Removal of surface-related wave phenomena: the marine case: 60th Annual Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1689. Foster, D . J . , and Mosher, C . C . , 1992, Suppression of multiple reflections using the Radon transform: Geophysics, 57, 386-395. Ikelle, L . T . , Roberts, G . , and Weglein, A . B . , 1995, Source signature estimation based on the removal of first order multiples: 65th Annual Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1478-1481. Jost, R . , and K o h n , W . , 1952, Construction of a potential from a phase shift: Physical Review, 87, 977-992. Matson, K . H . , 1996, T h e relationship between scattering theory and the primaries and multiple of reflection seismic data: Journal of Seismic Exploration, 5, 63-78. Mjedle, R . , Berg, E . W . , Strom, A . , Simamura, H . , Kanazawa, T . , and Fjellanger, J . P., 1995, Three-component ocean bottom seismographs used in prospecting off northern Norway: 57th Meeting of the European Association of Exploration Geophysicists, Extended Abstracts, P038. 137 BIBLIOGRAPHY Mjelde, R . , 1992, Shear-waves from 3-component ocean bottom seismographs off Lofoten, Norway indicative of anisotropy of the lower crust: 110, Geophysical Journal International, 283-296. Morse, P. M . , , and Feshbach, H . , 1953, Methods of theoretical physics: M c G r a w Hill, New York. Moses, H . E . , 1956, Calculation of the scattering potential from reflection coefficients: Physical Review, 102, 559-567. O ' B r i e n , M . J . , and Gray, S. H . , 1996, C a n we image beneath salt?: T h e Leading Edge, 15, no. 1, 17-22. Ogilvie, J . S., and Purnell, G . W . , 1996, Effects of salt-related mode conversions on subsalt prospecting: Geophysics, 61, 331-348. Orsen, A . , Amundsen, L . , Reitan, A . , and Helgesen, H . K . , 1996, Removal of water-layer multiples from multicomponent sea-bottom data: 66th A n n u a l Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1531-1533. Prosser, R. T . , 1975, T h e formal solutions of inverse scattering problems II: Journal of Mathematical Physics, 17, 1775-1779. Prosser, R. T . , 1980, T h e formal solutions of inverse scattering problems III: Journal of Mathematical Physics, 21, 2648-2653. Razavy, M . , 1975, Determination of the wave velocity i n an inhomogeneous medium from the reflection coefficient: Journal Acoustic Society of America, 58, 956-963. Riley, D . C , and Claerbout, J . F . , 1976, 2-D multiple reflections: Geophysics, 41, 592-620. Snieder, R., 1990, T h e role of the B o r n approximation i n nonlinear inversion: Inverse Problems, 6, 246-266. Stolt, R. H . , and Jacobs, B . , 1981, A n approach to the inverse seismic problem: Stanford Exploration Project Report, 25, 121-132. 138 BIBLIOGRAPHY Stolt, R. H . , and Weglein, A . B . , 1985, Migration and inversion of seismic data: Geophysics, 50, 2458-2472. Stolt, R. H . , 1981, Some basic stuff about linear operators and vector spaces: Stanford Exploration Project Report, 25, 227-249. Verschuur, D . J . , Herrmann, P., Kinneging, N . A . , Wapenaar, C . P. A . , and Berkhout, A . J . , 1988, Elimination of surface-related multiple reflected and converted waves: 58th A n n u a l Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 10171020. Verschuur, D . J . , Berkhout, A . J . , and Wapenaar, C . P. A . , 1992, Adaptive surface-related multiple elimination: Geophysics, 57, 1166-1177. Verschuur, D . J . , 1991, Surface-related multiple elimination: an inversion approach: P h . D . thesis, Delft University of Technology. Wapenaar, C . P. A . , Herrmann, P., Verschuur, D . J . , and Berkhout, A . J . , 1990, Decomposition of multicomponent seismic data into primary P- and S-wave responses: Geophysical Prospecting, 38, 633-661. Wapenaar, C . P. A . , Verschuur, D . J . , and Herrmann, P., 1992, A m p l i t u d e preprocessing of single and multicomponent seismic data: Geophysics, 57, 1178-1188. Ware, J . A . , and A k i , K . , 1968, Continuous and discrete inverse-scattering problems in a stratified elastic medium. I. Plane waves at normal incidence: Journal Acoustic Society of America, 45, 911-921. Weglein, A . B . , and Gray, S. H . , 1983, T h e sensitivity of B o r n inversion to the choice of reference velocity: A simple example: Geophysics, 48, 36-38. Weglein, A . B . , and Stolt, R. H . , 1995, I. T h e wave physics of downward continuation, wavelet estimation and volume and surface scattering. II. Approaches to linear and nonlinear migration inversion, in Symes, W . W . , E d . , Mathematical Frontiers in Reflection Seismology: S E G / S I A M publication. 139 BIBLIOGRAPHY Weglein, A . B . , Gasparotto, F . A . , Carvalho, P. M . , and Stolt, R. H . , 1997, A n inverse scattering series method for attenuating multiples i n seismic reflection data: Geophysics, (to appear). Weglein, A . B . , 1985, T h e inverse scattering concept and its seismic applications, in F i t c h , A . A . , E d . , Developments in Geophysical Exploration methods: Elsevier L t d . , 6, 111-138. Weglein, A . B . , 1995, Multiple attenuation: recent advances and the road ahead: A n n u a l Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 65th 1492- 1495. Yilmaz, 0., 1987, Seismic data processing: Society of Exploration Geophysicists, Tulsa. Appendix A Calculated series for 1-D B o r n scattering This appendix contains the calculated series expressions for the first and second primaries and the first order multiple i n Model 2. Using equations (2.40) - (2.42), I calculated the terms in the series up to eight orders and grouped the terms according to the type of event that they represent. Here, I present only the first five orders of terms. T h e first primary reflection is given by the series lfcgai , 1 fkla^ 5 2 (fcpffli This can be expressed as P (x ,z prl g g < z \k;uj) x = R ikx iM^-z ) ie 9e ( j g A2 where K-i. — — — = K vK a is the reflection coefficient corresponding to the first interface. 140 ; vo + vi (A.d) APPENDIX A. CALCULATED SERIES FOR 1-D BORN 141 SCATTERING T h e series for the second primary reflection is < ZilJfejw) P^ix^Zg +( 1 ^ 8 + 128 + , L y 8 ( V + V where a = g ) 1. -a? 8 2 16 + ( a, 4 a H 9 128 2 2 1 256 2 ) a 2 a a i + 4 y 5 1 8"32 ( a a v a 2 ~128 a 2 "2^6 0 5 o , 1 o 2 4 ) a ? + ( 1 5 A 16 " 3 2 ° 2 „•> - 256 ^ ^ + , 1 ( I . . . , , , l ( + 1 9 2 « 0 ^ + (- 192' ivl(z - z i ) + 1 ) 3 2 ; 16 " 256 a . 4 1 i^9 iM^ -z ) e e 2 8 a2 g 2 + —a + 2 128 a, + 7 1 2048 256 i + • • • j ( 1+ zV 8(z - z Q ( i - I ) 3 2 3! 7 , 3 + izy o2(z 2 - z!)( - 1) 7 l 2! (zV 2(z -z )( i-l)r 0 ' 2 n! 1 7 a )a 2 1 (A.5) V a 1024* 256 i2 1 i 21 2 3 ^ ^o(^2 - ^ i ) + . . . 2 0 a 4 Factoring this equation gives a ) 3 „ . 2 1 2048 (A-4) 2 ^o(z - z-x) + 1 x ( _ a 2 ) a i iz/ (z - ^ l ) , 1 1 .2 -3 , / «2 + ( - o8 + T ^ ° 2 + 7128"' ^ ° 2 +' 71024 W ^ ) + ( " i 6 - 64 128 27 13 -. \ "3 , i 1 + ( a 64 2048 2048 32 64 + 1024' 2 al ~ 2 21 -a? 256 i2 -o?^)a{ 2048 ~16"64 256" 27 . . , 15 , 2 ) « i + ( 7 ^ 7 o : ) i 2048 2048' A 3 128 2 1 4 i 2 2 ) a i - 32" ~ 2 5 6 a 2 f l + < z \k;uj) = g -, 2 2 a 2 k\ajv\. P {x ,z pr2 2 2 ; 192 v -o(2 :«2 0"> 32 1 512 256 e 128"" ' l b 2 4 ° 13 . , 1 a ) a + (64 2048 16 1 <"l6 ^ 1-2 3 +( ai a 2 «2 , 3 e 2 -32"64 ( = APPENDIX A. CALCULATED SERIES FOR 1-D BORN 142 SCATTERING which further leads to P (x ,z g g v A{z 2 e °e °( *- ^R2T1T[(l Zi\k;u>) = < pr2 ikx - ^ ) ( i - I) 2 2 7 iv 2z i^ S(z 2 0 - z ) ( 2 (iu 2(z . •. + i = k 9 x e - z )( x 2 R o( -2- g) 2z iv e - l) 7 l 2 0 z0(7i - 1) 3 3! 2! 0 3 x iv 2(z - + z (A.6) T — J"g>« l)2(22-2i)(7i-l) 2 z - 1)) 7 l / j[£ J 2 1 i / e ^z 9 e » l / o( 2 2 i _ z 9) ' e l / o2(z2-zi)7i ^ J 2 ^ g i f c 2 f l g J i o ( 2 z - z ) i i / i 2(22-21) _ , 2 / i where 71 = (1 — d i ) / 1 1 9 = 2 e T h e reflection coefficients from the first and second vi/va- interfaces are 2 - d - 2 ^ ( 1 - a i ) _ zy - f i ^> Pi * o / 7 i * \ x = (A.7) 0 ai ^0 + ^1 and respectively. - «2 R2 = a 2^(1 - ai)(l - a ) _ - x 2 d 2 — Ol ^1 - ^2 (A.8) Vi + ^ 2 T h e transmission coefficients through the first reflector i n the forward and backward directions are 2\ = 1 — Ri and T[ = 1 + Ri . Thus, the sum of terms i n (A.4) leads to a second primary which has the correct amplitude and propagation velocity. T h e series for the first order multiple is P m l t l (x ,z g < zi\k;uj) = g e^e-o(2(2z -2 )-2,) 2 n a* ~ e l " 2 " 64 ^ ~~z 1 2 1 2 1 + 1 64 1 - L A a-2 H 32 1- 64 L a\ ) d i + 32 1 + 1 ^ 32 1 a 2 32 + 5 64 A a 2 i . 4 & 2 (Xi ' 1 1 ' 32 d + 2 9 a ' 2 a 1 i+ —a a i ivo (z 2 16 - zi) + -a, 1024 128 ^ 2 \ ' - 3 i 1 . + + 512 2 1 a a, + 2 128 128 1 . a H 1 2 ' a , 512 ' 1 i 128 ' 1 "3 \-3 1 2 -2 13 ^2 1 . 3 2 128 1 1 a 1 16 a \ a\vl{z - Zif 2 2 1 32 (A.9) APPENDIX A. CALCULATED SERIES FOR 1-D BORN SCATTERING 143 This can be written as P (x ,z mla g g < \k;co) = - R (l - P ) ^ ^ Zl 2 x e _ [4(z - z )( 2 1 9 e -o(2(2, -, )-z )(^ 2 1 9 - l)] vl + 2! J2 (1 — ^2^2 ifcx = Pl(l - - ( g 2 Z 2 - _ C l 2 l ) ( 7 l e R l ) R * 9 e e X) _ ivo )(7i-l)] \ n\ — _ X 4 [^o4( 2 ll 1 w J ij/o(2(2z2-zi)-z ) ii/ 4(22-2l)(-Yi-l) s * * 9 e i M ^ - * 9 ) e e 0 i ^ 2 - Z l ) (A.10) giving a first order multiple with the proper amplitude and propagation velocity. Appendix B Elastic reflection and transmission coefficients In this appendix, I derive the elastic reflection and transmission coefficients for both an elastic free surface and an elastic-acoustic boundary. These are needed for the reference medium Green's functions for land and ocean bottom data. I employ the nomenclature of A k i and Richards (1980) for scattering (reflection and transmission) coefficients from a plane interface. E a c h coefficient contains an incident phase followed by a scattered phase. Hence, PS represents an incident P wave and a scattered S wave. A grave accent denotes v a downgoing wave, and an acute accent ' denotes an upgoing wave. Thus, PP and PP are P-P reflection and transmission coefficients respectively. Free surface reflection coefficients T h e derivation of the reflection coefficients for an free surface has been done previously by various authors. T o name a few: Ewing, Jardetsky, and Press (1957), A k i and Richards (1980), and Wapenaar et al. (1990). T h e derivation shown here follows closely the treatment of A k i and Richards (1980). T h e reflection and transmission of waves from an elastic free surface can be solved as a boundary value problem where solutions for incident and reflected waves are matched according to the free surface boundary conditions. T h e boundary conditions for any elastic 144 APPENDIX B. ELASTIC REFLECTION AND TRANSMISSION COEFFICIENTS 145 interface are continuity of traction across the interface. In the case of a free surface, the material on one side of the interface is a vacuum. Consequently, the tractions go to zero on a free surface. Although the atmosphere is certainly not a vacuum, the relative properties between the air and the solid earth are large enough so that the free surface approximation is valid. T h e first step is to consider a unit amplitude P wave incident on a free surface. T h e exp(i(kx + v(z — z))). potential for this incident wave is Reflection from the free sur- a face creates both a reflected P wave and a converted S wave. PPexp(i(kx -f- u(z + z))) 3 PS exp(i(kx + vz + rjz))) and s This gives the potentials where PP and PS are the P-P and P-S reflection coefficients. T h e displacements for this experiment can be written as d = x d z ( -d z i(kx+v(z -z)) e s \ d x + 0 p p p i(kx+v(z +z)) e s (B.l) Cj i{kx+vz +r)z)) e a T h e tractions are related to the displacements by T (j>d fid T zz u x z z J- X \d x x (B.2) idz ± where \L US the shear modulus, A is the Lame constant, 7 = A + 2fi is the uniaxial strain modulus, and T xz and T zz are the shear and normal tractions respectively across a horizontal plane. B y combining these equations, the sum of the tractions at the free surface is -2pku -/i(7/ 2 /J,(T] 2 - k) 2 -2\Lki) k) 2 pp i(kx+uz ) e p s Cj i(kx+vz ) e a + 2fiku -fi( 2 V n(r) - k ) 2 - k) 2 2fik 2 V i(kx+vz ) e 3 0 = 0. (B.3) This equation can be solved for the reflection coefficients PP and PS provided that the phases i n the exponentials are equal. This latter condition means that k is the same before and after reflection, which is simply an expression of Snell's law. A similar calculation can be performed for an incident S wave of unit amplitude giving a reflected S wave and a APPENDIX B. ELASTIC REFLECTION AND TRANSMISSION COEFFICIENTS 146 converted P wave. T h e resulting reflection coefficients for both types of incident waves are 4k vr) — {r}- k ) 4k wq + (T/ - k ) 4:krj(i] —k ) SP = 4k vri + (n - k ) — 4 f c ; y | T7" —-k —Aku(r] -fc') I PS 4k ur) + (n - k ) Ak vq — (rj - k ) SS 4k uTj + (rj - k ) 2 PP = 2 2 2 2 2 2 2 2 2 2 2 2 2 (B.4) 2 2 2 2 2 2 2 2 2 2 2 2 2 2 Water bottom reflection and transmission coefficients To derive the water bottom reflection and transmission coefficients, an appropriate model to consider is reflection and transmission from an elastic-acoustic boundary. Continuity of traction across the interface is again one of the boundary conditions. Since the shear traction is zero in a fluid, the shear component of the traction is zero at the interface. Another boundary condition is continuity of displacement. For two media in welded contact, both the vertical and horizontal displacements must be continuous. A t a solid-fluid contact, the horizontal component of displacement may be discontinuous (but not necessarily zero) because the fluid may slip past the solid ( A k i and Richards, 1980). Barring cavitation, the vertical component is, however, continuous across the interface. boundary conditions are: In summary, the three continuity of normal traction and normal displacement across the interface and zero shear traction at the interface. There are three types of incident waves to consider. T h e first is an incident P wave in the fluid. T h e second and third are incident P and S waves in the solid. A l l of these cases create reflected and/or transmitted P and S waves in the solid and reflected a n d / o r transmitted P waves in the fluid. Using the same techniques as for the free surface, a set of equations can be obtained for the three boundary conditions. Solving these equations APPENDIX B. ELASTIC REFLECTION AND TRANSMISSION COEFFICIENTS 147 yields the following reflection and transmission coefficients 2u p u (r]l 2 J PP -k ) 2 1 + #P2*i ( - U) = — PP v 1 -^wlfi r + PP = SP 2 V2 ((vl - k ) + f3 2Vi P 2 2 2 4 u k )) V2 2 (B.5) 2 ^2u> p v {j} - k ) 2 PS 2 2 2 PP PS (vl ~ k ) + 4: u k )) = = 2 2 2 2 Akuiu uj pi 2 2 ^wb -±kv vifi\p {r) 2 2 2 - A; ) 2 wb Ar) v k^ p {j] - k ) 2 2 x 2 2 2 (B. ) = - " " " r ^ " 'wb 6 _ -4kj] v u>lp 2 2 2 ^wb -a>VizV/3 + ( ' (Vl - k ) + 2 2 S S 2 = ^vk) 2 2 2 2 n where 0 w b = wW # + ftp v ((nl - A; ) + A u k )) . 2 x 2 2 V2 2 2 (B.7) Appendix C (P,S) decomposition and up-down wave separation Cl Land data Here I the obtain the source and receiver (P, S) decomposition matrices for land data. This derivation has been done previously by Wapenaar et al. (1990). I begin with elastic displacement data as i n (3.15) T> = g vwg . 0 The (ci) 0 reference Green's operators can be expressed in the form a, = ^ n G n = ^ n ^ G * + G*)n. T This follows from the fact that away from the source, V Gp 2 V~ Gp 2 0 = — a /u> Gp . 2 2 Similarly V 0 (c.2) 0 _ 2 G s o = —ft /UJ GSO2 0 = —u> /a Gpo] T 1 r 0 hence 2 In operator notation then 2 v- n r- G = — ^ n G . 2 2 (c.3) 0 pu> z Using this form of the Green's operator, the data can then be written as T> = ~n (G T D 0 + G )(nv( )^n )(G^ + c )n 0 s 1 ^ p = ~n (G T D 0 T 0 " (C.4) + Gf)v< >(Gj + G )n. l 148 s 0 s APPENDIX C. (P,S) DECOMPOSITION AND UP-DOWN WAVE SEPARATION 149 A p p l y i n g T> to a delta function S(r ) and substituting equation (3.55) into the receiver side a of (C.4) leads to T>(x , z \x , z ;u>) g s 3 s f°° f°° 1 pw J_ o —oo 2 co e -' g g v z 0 2iv g 0 CO / -oo + 2ir/ V^\x', e > " s z S Pe™ * 9 3 s 2iv n e >9 PSe ><> s SSe K> s z ir z'|r")Go(r"|r'"; )US(T'" U dk Q ,7 ir g fiOO / J— p p z z 2XT)g dx'dz'dr"dr'" . - r.) (C.5) co Here, the vertical wavenumbers refer to propagation in the reference medium. Focussing on the quantity in curly braces, applying I I T to the receiver coordinates gives different results depending on which way the waves propagate. Explicitly, this gives / n T Q Jkg(Xg-x') x')[ °=£ ikg \ -Wg lg ir e -'"9 9 g irfg Z l/g % kg ~iVg g z e 0 ikg ikg 0 z 2iu 2ir) g IVgi 2iv„ PSe >° ° ir z dkg. £Vg*_ SSe ' » ir, z (C.6) 2irj g Next, the receiver depths are set to zero to simulate surface mounted sources. This reflects a realistic practical consideration however is not necessary; the decomposition can also be done with buried receivers (and sources). Setting z = 0 gives g ik CO / g ikg irjg ikg(Xg-X>) PP SP e' 9 t'kg PS SS 0 e -ivg •OO 11/g ikg -i: Now, -ir/g Jkg(x -X') g T\T —1 9 v 0 z 2iv dk g 0 2ir) g 0 2iv dkg. g 0 (C.7) 2ir) J g except for the N " , the quantity on the right hand side of this equation is exactly 1 the directly propagating Green's function G^Xg^Zg eliminating N " 1 , TI (GQ T = 0|a:',z';u;) (see equation (3.55)). B y -f G Q ) is effectively replaced with G „ and the receiver ( P , S) decomposition and deghosting is accomplished i n one step. T o reiterate, the steps are: 1) Fourier transform the data over time and the geophone domain; 2) apply N G ; and 3) APPENDIX C. (P,S) DECOMPOSITION AND UP-DOWN WAVE SEPARATION 150 inverse Fourier transform back to space and time. T h e result is just upgoing P and S waves measured at the receivers. In the (UJ, k ) domain, the receiver decomposition matrix is g -i(3 2 2v„ "9 UJ" g (C.8) 2 K > Vg V 9 T h e source side decomposition can be obtained in a similar fashion. T h e objective is to convert G H6 to G S in (C.5). In this case, the operator II can be treated as if it operates d 0 on the quantity to the left. This because when G operates on IIf5, the sifting property of 0 the delta function acts so that the derivatives are applied to the source coordinate of GoNow / G (X",Z"\X ,Z ;UJ] 0 S + i k s ( x " - x 0 s Te '(*"- >) i v 3 0 0 2iv a 0 z 2iu„ j ) PP SP 2iv 0 0 2i 3 0 PS Vs 0 2 k ss ^ 0 2ir] s \dk U. (C.9) 3 e ° ir)sZ Setting the source depth to zero gives 0 i k s ( x " ' — x 3 0 2iv ) a 0 ik s —iv 3 ik 3 2ir), 2iv 0 0 2ir] 3 + s PP SP 2iv PS SS ^ 2iu s ik (x"-Xt & / s M7 s 0 •00 1 iv s 2ir) _ 0 oo —ik ^ s s I dk, —irjs ik s dk. 2ir] s oo e- > °G (x",z"\k ,z IK / X d s z = 0;uj)M; dk . 1 s (CIO) •00 This equation says that by Fourier transforming over the shot domain, multiplying by the matrix M , then applying the inverse Fourier transform, this gives the downgoing P and 3 S waves at the source. In the (uj,k ) 3 M = 3 domain, the source decomposition matrix is -i/3 2 UJ' 2rj„ ~~~2^7~ K (CU) s APPENDIX C.2 C (P,S) DECOMPOSITION AND UP-DOWN WAVE SEPARATION 151 Receiver decomposition for ocean bottom data For ocean bottom data, the equations are similar to land data. Starting with the equation for displacement data T> = -^n (G T + G d 0 W ) V C ( 1 ) ( G * Go )n, + (C.12) c where ( G ^ + G ^ ) is the Green's operator for the water bottom reference medium. 0 particular form of the Green's function related to G 0 A will be useful both for decomposition and for formulating multiple removal. For sources and receivers in the solid, this form is r 1 Go(x ,z \x„z ;u) g g = — J iv l*s-*gl t e ^tk^X g Xg o 0 2 0 R wc + p 2JJ72 ? w v c e e •i/ (z« + z ) 2 9 2i^ ps «("2*.+T)2*g) R ss 1 W C sp l 2 e <12 (*. + * « ) dk (C.13) 2i m where -iv22z wb PP + Z (PPPP pwc - 2 rtpp 0 -iz~wb(v2+n2) PPPP) e K - P S 2 (^+12) pwc 0 in which 0 = , , PS + Z (PSPP SS + Z {SSPP pwc K r 2 sp n SP + Z (SPPP 0 ss 2 0 1 + PPZ , 2 Z = e - - SPPP) - PPPS) (C.14) SPPS) , v and rj are the P and S vertical wavenumbers lVl2Zwb related to the receiver wavenumber, and the subscripts 1 and 2 pertain to the water and solid. T h e utility of (C.13) is that it isolates the portion of G and receiver coordinates. W C that depends on the source Substituting this equation into the receiver side of (C.12) and APPENDIX C. (P,S) DECOMPOSITION AND UP-DOWN 152 WAVE SEPARATION operating on a delta function leads to ^ ' • [ ' " ' " a ) = ik />oo «co I /too ik ir]2 ^ L L v - J —iv ik 0 2 \IV w c V^ •CO t/ —CO 2 g V> 1 ^wc e z e z 0 z z 2iV2 z 2ir,2 \dk R z > " , z " ) G ( r " | r " ' ; )US(T'" 0 z '(V2 ' + "2 g) PS 2 2 ("CC O O fOO/>co / iv (z' + z ) e 2iu H"2 ' + V2 9) R SP 2iu If 2 / PC PP —irj 2 0 iii2( '- g) e - p.) U dx'dz'dr"dr'" . This portion i n curly braces can be expressed as ik irj 2 —iv ik 2 + ifc iv 2 ^c iu 2z —irj R 2 pe 2 JV2(z'-Zg) R™ i(V2+V2)z g e 0 g 2iv2 e 0 ik 0 g 2iv>2 J g iV2( '- g) 0 Recognizing that the integral without N " function, it follows that N g 1 z 2ir)2 iv2(z'-z ) e ik(x -x') e dk 'V2( '~ g) z e z z dk (C.16) 21772 is just the directly propagating (P, S) Green's is the matrix which transforms the data into upgoing P and S waves at the receiver i.e. G Q V ^ G O . A s i n free surface decomposition, the steps are: 1) Fourier transform the data over time and the geophone domain; 2) apply the matrix N ; and 3) inverse Fourier transform back to space and time. W h e n the receivers are on g the ocean bottom, z = z b and the receiver decomposition matrix in the (UJ, k ) domain g w g becomes •022U2 -iff UJ' r> 2~k g 2V2 2 2 K9 ^~ + 1-Z 1+Z 2 2 v\ P2 2 P\kg_ \-Z 2p fir72 1+Z 2 2 (C.17) 2 It is remarkable that such a simple expression comes out considering the complexity of the water bottom reflection coefficients. Comparable expressions have also been obtained by Orsen et al. (1996). (C Appendix D Elastic interface multiple removal series In this appendix, I derive i n detail the equations for both land and ocean bottom multiple removal. T h e steps in the two derivations are basically the same and differ only in the Green's function which creates the unwanted multiples. I start with the elastic scattering equations in (P, S) coordinates. T h e data are modelled by the forward scattering series as D = G V G o + GoVGoVGo + ... (D.l) 0 where D are the data in (P, S) coordinates, G o is the reference (P, S) Greens' function for an elastic half space, and V is the perturbation or the model. A n inverse series for V can be formulated as (Weglein and Stolt, V = V( ) + V( ) + V( ) + . . . 1 2 2 1995) (D.2) where the rc'th term is n'th order in the data. Substituting (D.2) into ( D . l ) and equating terms that are equal order in the data gives D = G V«Go 0 (D.3) 0 = GoV^Go + G O V W G O V W G O (D.4) 153 APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL 0 = GoV( )Go + G o V ^ G o V W G o + G O V W G O V ^ G O 3 + SERIES 154 G O V W G O V W G O V W G , , (D.5) This is the inverse series which will be used to obtain an inverse subseries which removes multiples from both land and ocean bottom data. The form of V, V , and (1) Before I proceed to solve for terms in the inverse series, certain aspects of the forward and inverse scattering equations need to be explained or clarified - specifically the form of the perturbation V . Consider the elastic displacement Lippmann-Schwinger equation i n operator notation (see equation (3.11)) G = Go + GoVG (D.6) or written i n terms of functions />oo / ^oij(rg|r'; UJ)VJI(T'\T"; UJ) £ < ( r " | r ; UJ) / •oo J—oo CO m s dr'dr" (D.7) where £/ (r \r ]uj) is the Green's function i n the actual medium, ^oim(r |r ;UJ) is the irn g s g s reference medium Green's function, and VJI(T'\T"; UJ) is the function which pertains to the perturbation operator V = C — C. It is not obvious why Vj>(r'|r"; a>) should depend on four 0 spatial variables and not two. T h e reason lies in the fact that V is a differential operator which means that it has a non-local form i.e. it depends on more than just position. T o see why a differential operator is non-local, it is easiest to consider the variable density acoustic wave equation APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL 155 SERIES where p is the density, K is the bulk modulus, and p is the pressure. Written in terms of r( • an integral wave operator, this equation is tp= v PV) wf ) 2 V+ s(r r = _ X(R|R,; W M R ')*' where L(r|r'; UJ) is the wave function for the scalar wave operator L. T h e wave equation in a constant reference medium can be written as L p 0 = 0 and can also be expressed as an integral operator equation Lp 0 = • - V + ^-u> ) [°° (v ./-co 8(r - r') ( ')dT' 2 PO V A O P r = / ° ° i ( r | r ' ; u,) (r')dT'. 0 P (D.10) J-oo / Subtracting equations (D.10) and (D.9) gives CO ( V • a ( r ) V + a (r)u, ) p / K 2 6(r - r')p(r')^r' = POO / J •oo V(r\v'; u)p(T')dt'. — oo (D.ll) where V(T\T';UJ) = ( V • ap(r)V + a (r)u>2) 8(r - r'), a (r) K 1/i^o — l / ^ ( ) r p = l//9 0 - l//>(r), and aK(r) = Now, since the derivative of a delta function is not, strictly speaking, a local function, a differential operator is non-local; hence V depends on more than just position. This has also been discussed by Stolt (1981). If the mass density i n the earth does not vary, then a p = 0 and V(T\T';UJ) looks like a^(r)a; 5(r — r'). 2 Since this function has values at only r = r', it is purely a function of position and V(r|r'; oo) can be written in the simpler form V(T;U>). T h e lesson is that when wave operators contain both functions of position and spatial derivatives of functions, this implies a more complex earth model and a more general form of the perturbation is required. These arguments can be extended to the elastic case where the wave operator contains spatial derivatives. This justifies writing the elastic perturbation function as V(r'|r"; UJ). T h e same can be said for its (P, S) counterpart V(r'|r";u;). While it is true that V £ ( r — r') is a non-local function, it is very nearly local. Consequently, it may be possible to use a local perturbation operator in practice. However, by writing V(r'|r";w) in as general a form as possible, the perturbation function can account for very complex earth behaviour not described by a simple differential wave equation. Recall that V is simply the difference in wave operators between an actual medium and a APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 156 reference medium. Since the form of C need not be actually specified, it can represent very complex wave behaviour. In the following development of multiple attenuation, I also express "V^(x, z\x', z'; u>) with five independent variables: four spatial and one temporal. This follows from the form of the inverse series which is written as V = + V^ ) + V^ ^ + . . . . 2 3 Since the most general form of V is V ( a ; , z\x', z'; u>), it is reasonable that its constituents should have the same dependencies. In practice, the additional degrees of spatial freedom may be needed to describe angle dependent reflectivity and anisotropic effects as a function of position. In general, since the data depend on only three independent variables, only a portion of can be 'illuminated'. To perform an actual inversion of the data into necessary to cast it is in a form which depends on three independent variables V ^ ^ a ; , x', z). This form corresponds to uncollapsed migration. Since migration can be thought of as downward continuation of the data followed by an imaging condition, uncollapsed migration involves downward continuation only. T h e form of V ^ ) and its implications is discussed in 1 more detail by Stolt and Weglein (1985). Since I am not explicitly performing inversion, I can write a general, albeit impractical, form for V ^ . T h e saving grace is that V W and the multiple removal terms which are derived from it are transformed back into data space, thus avoiding the issue of how to access a function which depends on five degrees of freedom. I should point out that in the following formulation of elastic free surface multiple suppression, a non- local V W is not a necessity - a local operator would work just fine and the formulation would be the same. However, by using a general perturbation operator, the multiple suppression algorithm is valid for a wide class of earth models. 157 APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES D.l Free surface multiple removal Determining V t ) 1 T h e first step constitutes a linear elastic inversion where V ^ is solved for in terms of the data D . To begin, I rewrite (D.3) using the elastic free surface Green's operator Dp DPS D D P SP PO "T (1) V Vpp ' K PS ^PP U SP fs G" PS ^PP PO v > v: SS SP ss f/(l) V v, PS a) SS G Qis SP d so + GSS k (D.12) Expanding out the first term, bpp = (G P0 + G% ) + G% Vgp (G + Gpp) Vps + GpS Vgg G gP , P + (GPO it is clear that Dpp + G% ) Vpp (G po p G*SP s P0 + G% ) P (D..13) f contains contributions from all of the elements i n V ^ . This is due to the conversion between P and S energy at the free surface. Equation (D.13) also shows that the scattering history can be determined from the subscripts by reading them from right to left. Since all of the elements in D depend on all of the elements in V ^ \ all elements in D will be required in order to find just one element in V ^ ) . As a preview of things to 1 come, the multiple removal procedure does not require any direct information about the earth properties themselves. is derived from the data. A l l information about the earth is contained in V ^ ) which 1 In fact, V ^ ) can be thought of as representing reflectivity as 1 opposed to mechanical earth parameters. Returning to the details of the calculation for Dpp, CO />CO / / •co /"CO / />co I write dxdz dx'dz' • J — co J — co J —co (Gpo + G%)(k , z \x, z; UJ) V$(x, z\x\ z'- UJ) (G g d g + G% (k , z \x, z; UJ) V$(X, z\x', z'; UJ) (G s + (G P0 g g Po + Gpp)(k , g z \x, z]u)V$(x g t P0 + G% ){x, z\k„ z ; UJ) P + G% )(x, z\k , z ; UJ) P s ] a z\x', z'-,u) Gsp(x, z\k., z,\u) + G%(k ,Zg\x,z u,)V&\x,z\x',z'-,uj)G%(x,z\k ,z ;uj)\ g s s s (D.14) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 158 which becomes after substituting equations (3.55) flOO r Dpp(k , z \k , Z ;UJ) = / g g 3 pOO poo poo / S / / dxdzdx'dz'- v — co </ - oo J — oo J—oo (PPe > > + iv z ) _il/sZs e 2iV, (PP _i_ ,iu z a s —ii/ 2 ^ s s 2ZZA, 2*% 2ZZA, 2ii/, e- (*- -^Vii (x,z|^z>)e ( + SP- i 227/ B ) I F C ' , B , + ''' ' P5Z (D.15) . ) 2zz/ g s Recognizing that the integrals in (D.15) represent Fourier transforms, (D.15) becomes Dpp(kg,z \k ,z ;uj) g s s pp( 9, ~Vg\k ,U,\u) V 2iu a k (PPe ' + ~ ivsZ + SP ——V^ \kg,-r,g\k ,iy -,u ) : 2irjg P s + 3 (PPe ll/ Z 3 t e 3 I -11/,z ) s 2iv, ) ivsZs 2iu, ) -V£s( o>- g\ »>y''> ) -2J^ k v k u PS IV, + S P -v \K,v.;") ^ — ss(h, v s^— p g Z, (D.16) A similar calculation shows that each term in D depends on the quantities pp( a> v k ps( 9 v k »\ ) u ~ g\ ^ v k sp\ 9i-v \ s, u> V k k 9 (D.17) V \k ,-r)g\k ,rj ;u>) i~Vg\ks,f ;^) ss s g s s with different coefficients that depend on k , z , k , z and UJ. This can summarized byg g 3 s writing T)(kg, z \k , Z ;UJ) as the product of three matrices g S 3 D(k ,z \k„z.;u) g g = q (k ,z ;uj)Y^\k ,k ;uj)Q (k ,z ;uj), g g g g 3 3 3 (D.18) 3 where D(kg,Zg\k ,Z ]Uj) 3 3 D (kg, z \k , z ; UJ) D (kg, z \k , z ; UJ) PP g 3 3 PS g 3 3 D p(k , z \k , z ; UJ) D s( g, z \k , z ; UJ) S g g 3 3 S k g 3 3 (D.19) 159 APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES (PPe" 9 9+e-"'9 <>) 2iv„ / Q (k ,z ;uj) g V (1) g (*,,*,;<-) g V z PS Vpp = z SP - 2ir) ( iv 2iu„ k 1 z g g 17,; UJ) (D.21) V^J (k , -r] \k ,T] ;uj) k a i7 Vp ] (k , -u \k„ S P { ,-Vg\ siVs]u) k z e v,\UJ), (h,-Vg\ *, (D.20) g (SSe 9 s+e~ >9 9) 2iri e">9*i> g g A s and (pp »'.*. -«'«*.) e Q,(k„z ]w) SP - 2irj +e e 2iv (D.22) s 3 a PS'-2iv 2iri s s Note the symmetry between the subscripts and the Fourier transform variables i n V ^ ) . 1 This is related to the types of scattering performed by each of the elements of V ^ . Equation (D.18) shows that the data, which have three independent variables k , k and g s UJ, are related to the product of three matrices each of which also depends on these same three independent variables. It is this property that allows V ^ ) to be determined from 1 measurements away from the scattering region. T h e task remains to find the elements of from (D.18). Since Q and Q are both 2 g g by 2 matrices, their inverses are particularly simple. Pre- and postmultiplying (D.18) with these inverse matrices yields V^(k ,k ;u>) g =Q- (k ,z ;uj) 1 s g D(k ,z \k ,z ;u;)Q; (k ,z -u>) g g g 1 s s 3 s (D.23) where (SSe' 9 9+ -' >9 9 v z r e 9 -SP e'"9 2ir) z z 2ir) g det Q —PS s g (D.24) B (D.25) (PPe 9 9+ - 9 9) 2iv„ iv ?iVgZg 2iu„ z e iv z and (S5e"" '+e~~"» ) x Q detQ T h e invertibiHty of Q matrices. If det Q g 2irj~s - i g -PS*- 2iu, s and Q and det Q 3 s x< 'Vs a 2ir} (PPe >> »+e- ' > *) -SP iv z z e i , z 2ius depends on the properties of determinants of these do not have any zeros over a suitable range of k, then APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES Q g andQ 160 are invertible over that range. B y plotting the determinants of these matrices 3 as a function of k, I found that they are nonzero for \k\ < uj/a which is the range for nonevanescent P waves. A s discussed by Claerbout (1985), this is a reasonable restriction on k. In addition to invertibility constraints, the inverse matrices must be stabilized to avoid zeros i n the denominators of these functions. Determining the multiple removal subseries Carvalho, Weglein and Stolt (1992) derived an inverse subseries which performs surface multiple removal for marine data. M y purpose is to formulate an analogous subseries for elastic data. Similar to the acoustic case, m y aim is to derive a formula for the entire subseries in terms of the data. In the elastic inverse series, the reference Green's function contains two portions: G and G . T h e latter is responsible for the creation of free surface D S 0 multiples in the forward series. It follows that G^ must therefore be responsible for the removal of these events in the inverse series. I select out those terms in the inverse series which involve only G v'( ) = V'( ) = f s 0 2 3 and write them as a multiple removal subseries. Hence, - V W G ^ V W -V'^GfivV - V « G 0 S V ' ( 2 ) - V W G J V W G J V W (D.26) In functional form, this first equation can be written as Vp$(k ,-u \k ,u,\w ) V s\k ,-v \k ,r} ;uj) V p\k ,-r] \k ,u ;uj ) V s\k ,-ri \k ,r) ;uj) g { s g g { P t g s ( s » G% (X,Z\X',Z';UJ) P g s g 7 , / . . . / g a g s \ PP V (a;, G] PS z\x', z'\ui (x, G\ SS z\x', z'\ UJ s a ( 9> k (D.27) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL 161 SERIES where I have Fourier transformed over the shot and geophone domains on each side of the equation. This is done so that V ' ( ) can be added to V ^ ) in (D.21). T h e first component 2 1 for V ' ( ) is 2 2 Vpp(k ,-Ug\k ,u iu) g s = — / t /*°° / f* oo poo / / fioo poo / dkdxdzdx'dz' v — C O v — O O v — OO v — OO v — OO i(k(x-x')+v(z+z')) Vp \k , {1 P z-UJ) PP -u \x, g g + Vp ^\kg,-u \x,z;u>)PS S Vp W(x', z'\k , v ; UJ) P IV i(k{x-x')+vz' +t)z) a Vp ^\x',z'\k ,u ;uj) e g a IV P s s „i(k(x—x')+uz+T)z') + Vpp^(k , -u \x, z-UJ) SP g VSP^KX', z'\kt, va; UJ) : g ITj x )+rj(^-|-z )) / + V \kg,-v \x, {1 PS z; UJ) SS- g / V \x', z'\k , V ;UJ). (1 : SP IT} a (D.28) S Recognizing that the spatial integrals are Fourier transforms, this equation reduces to pp{k ,-v \k v \uj) V g g ay = -:— / a + V {1 PP -v \k,u;uj) -r—V p^\k, g> g P PS (kg,-v \k, {1) PS dkV \k v, w) — V W(k,-v\k ,v ; g PP a g 3 (k, - \k , v ; UJ) {1) SP 17] s UJ) a XV SP + Vpp^Xkg,-v \k, v; UJ) — V -v\k ,v ;uj) V s a SS + Vps (k ,-v \k, UJ) {1) Similar results can be found for S / V'( \k , k ; 2 g a g V p, Vp \ S oo where g and V SP a] (D.29) leading to the matrix equation K } S S V^\k , k;u>)G *(k;uj)V^\k, k ;uj)dk { g •oo V W{k, -r,\k„ v uj). — m a (D.30) UJ) is as shown in (D.27), Vppih,-Vg\k ,Vs\u) V$(k ,-v \k ,r) ;uj) V£p\k , -v \k ,v ;uj) V '{kg, -r) \k ,v ;uj) s g g a a g g ss a g s a I (D.31) a and PP/iv SP/irj PS/iv SS/ir) (D.32) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES It is possible to determine the n'th order multiple removal term from 162 using the recursive relation oo V ' t " " / 1 k ;oj) dk. ^ , k;u,)G0s(fc;u;)V«(fc, s •oo (D.33) This is a consequence of the fact that the exponentials in G Q do not contain absolute value signs. Data reconstruction and multiple removal Now that the multiple removal subseries has been found, the next step is to reconstruct data using this subseries. T h e reconstruction is done by sandwiching V ' between two directly propagating Green's functions. This gives the data series D' = G VG d d G f t V * ) + V ' ( ) + V ' ( ) + . . . )G 1 = 2 3 (D.34) d = £ > « + D'( ) + D'( ) + . . . . 2 3 The equation for the first term is n(i) 'PP u ^PS SP rid {rWrid r(l)Ad rid \/-V-)rid ^PO PP ^PO ^PO V ^50 SP ^PO ^SO V V rid T/(l)Ad SS V PS ^SO SS ^S0 Fourier transforming over the shot and geophone domains, D becomes PP OO / fiOO / flOO fi / / (D.35) rid T / M A l oo dxdx'dzdz' •OO J—OO J—OO J—OO ,i{k x'+u { '-z,)) i{-kgX + Ug{z-Zg)) a e 2iu,9 2iv„ — V p p ( a ; , z\x', z ;u>) -V^p\k ,-u \k ,u ;oj)g g a s s Z 2%v, (D.36) 2iv„ This is possible because the scattering points are always lower than the sources and receivers. A similar calculation for the other elements in D ' yields the matrix expression D'W(k ,k ;u) g t = K (k ,z ;u,)V( \k u>)R (k.,z. u>) g g g 1 ] t ] (D.37) 163 APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES V^(k;u) where is given by (D.21), —"g g Rg(k , k ;uj) = g 0 z e 2iv a a (D.38) - >a'9 2ir] ir e g and similarly for R ( A i , z \ UJ). T h e matrices R s s 3 and R g can be said to project a onto data space. From here, I substitute (D.23) into (D.37) to obtain D' (1) (A; , g k ; UJ) = R (k , a g g z ;UJ) g Q (k , 1 g g T)(k , z \k , z ;uj) Qj (k z ;UJ) g g g s 1 a z \ UJ) R (k ,z ; 3) a a a UJ) 3 (D.39) which effectively removes the source and receiver ghosts from the original data. T h e source and receiver deghosting matrices have the form e~ 9 -(SSe s » xv [R.QJ ] 1 (k ,z ;uj) = g g ir) det O, + e~ s e) z iri —SP iUgVg z (PPe - >9 9 -PS ir e il/gZg z iVgT)g + e~ \ iVaZg \ /_ (D.40) and e [Q: R ] l (k ,z ;uj) 3 3 a z ^(SSe >° ' —'v s 3 1. detQ ir -SP^- + e~ °) z iruz (PPe * + e- ° *) iu -PS a Zs iu (DAI) For the second term in the multiple removal series, f>( )' = G ^ V ( ) ' G . 2 2 In Fourier d transform space this becomes T>' (k , z \k , z ; UJ) = R (k , z ; uj)V'W(k , k ;uj) R (k , {2) g g a oo 3 g g g g g 3 a k;UJ) G (fc;UJ) s 0 a Z ;UJ) S V^(k, k ;uj) R (k , z - UJ) dk a a a 3 • CO O O z \Qg\k , 0 / g R (k , z ; UJ) V^(k , g / g g z ; UJ) D(k , z \k, z \UJ) Q ; ^ , z ; UJ) R (k, z ] UJ) g g g 1 a a a 3 O O • R;\k, z ; UJ) G (k; u))R-\k, z ; UJ) s a 0 g • R (k, z -UJ) Qg^k, z - UJ)D(fc,Z |k ,z ;UJ) g g g G 3 a Q 7 (A;,, z ;u>)R,(As,, za;UJ)dk 1 a (D.42) z APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES where I have used (D.23) to substitute for R7 1 G ^ R ; 164 and (D.30) for V ' ( ) . Expressing the matrix 2 as H and using (D.39), I obtain 1 CO B^(k , z \k, z ;UJ) H(k; UJ) D™{k, g / -co g a z \k , z ; UJ) dk g a (D.43) a where ivPP H(A3;w) = - iuSP (D.44) irjPS ir]SS 7T T h e term D ' ( ) removes all first order elastic free surface multiples from the deghosted data 2 by using the deghosted data to model the events which reflect from the free surface. For the n ' t h term in the series, I obtain D^ \k , n z \k„, z ;uj) = R (k , g g a g z g ; c ^ V ' ^ g 1 ^ , k ;UJ)R (k , a a z ;UJ) a a CO R (k , g / g z -UJ) V^~ \k , 1 g g k; UJ) G (k; UJ) V^(k, k ; UJ) R (k , { s 0 a a z ; UJ) dk a a CO oo / R(kg, Zg' UJ) Q-^kg, Zg, Uj) D'^^kg, Z \k, Z, J Uj) Qj^k, Z \ Uj) R (k, g B a Z \ Uj) a OO • R~ (k, z -UJ) G%(k-uj)R^(k, x • R (k, g z ; UJ) Q~ (k, 1 g z -UJ) a g z ; UJ) D ( f c , z \k , z ; UJ) Q7 (fc , z ; UJ) R (k , g g a a 1 4 a a a z ; UJ) dk a oo / D ' ^ " 1 ^ , z \k, z -UJ) H(fc; UJ) D « ( f c , z \k , Z ;UJ) dk . g a g a S (D.45) CO In general, the term D ' ( ) removes all events which undergo n interactions with the elastic N free surface. This includes both reflected (P-P) and converted (P-S) waves. A s in marine demultiple, the strength of this method is that the only needed information is the near surface elastic parameters, the source and receiver depths, and the seismic data with the source wavelet deconvolved. If the source wavelet is not removed from the data, it must be removed from each term prior to multiple removal. This is discussed i n Section 3.6. It is interesting to see what happens when the reference medium becomes water. that case, the data elements that contain shear components are zero and PP In reduces to APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 165 -1. Equation (D.43) then becomes I D'i \k ,z \k.,z.-,u,) P g g poo = — / JDS(fc„ z |fc,z ; )2.-|/ <" ")^(* z |fc. z.;u;)(ifc l B J W i + C 27T J-oo which is the equation derived by Carvalho (1992) for marine data. !l a > (D.46) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES D.2 166 Ocean bottom multiple removal T h e steps i n deriving the multiple removal equations for ocean bottom data are essentially the same as for land data. T h e difference lies in using the Green's function for a water layer on an elastic half space. Determining I start with an equation for four component (P, S) data poo CO POO I / Gi{x , z \x g g •oo z; UJ) V (x, I dxdz dx" dz' J— oo J— oo J— oo z\x', z'; UJ) (1) t POO I 0 <f> (x",z"\x ,z ;uj) P a a 0 (D.47) ip (x",z"\x ,z ;uj) s where the scattering points are located below the source and receivers. a s T h e source side Green's functions describe the P and S energy transmitted throught the water bottom by the P source i n the water column. Physically, the P and S waves are transmitted into the solid at the same time; different explosive sources will give rise to the same relative mix of P and S waves. Equation ( D . 4 7 ) is therefore an idealized data which will serve as the input into the multiple removal algorithm. T h e actual two component data will be processed to approximate the ideal data using the source and receiver ( P , S) decomposition. Using G Q from ( C . 1 3 ) and Fourier transforming with respect to both shot and geophone domains, ( D . 4 7 ) becomes -<"2g*9 0 e D^\k ,z \k ,z ;uj) g g s s = 2iu 2g ~'V2g 9 z 0 2iV2g J -"2 \k , v ; UJ), Vjiy (k ,-v \k , VSP (k , -V2 \k , v ; UJ) V^J (kg, ~V2 \k , Vpp 9 g a 9 a 2a 2a g PPe 2g g —lV2aZwb r/ ; UJ) a 2a a rj ; UJ) 2a 167 APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES or expressed as a matrix equation B(k ,z \k ,z ;uj) g g 3 = R (k ,z ;oj 3 g g )V^(k ,k ;uj)R (k ,z ;uj). g g 3 3 3 (D.49) 3 Obtaining "V^ ) is then a simple operation in (kg,k3,uj) space 1 V ( 1 ) (fc , 3 k.; UJ) = R^ikg, z ; B(k , z \k , z ; UJ) g g g s UJ) s R^iK, z ; s UJ) (D.50) Determining the multiple removal subseries T h e same logic which led to the free surface multiple removal scheme can be used to obtain a multiple removal series for ocean bottom data. T h e portion of the reference Green's function which builds multiples in the forward series is used to select a multiple removal series from the full inverse series. Hence, V' ( 2 ) = V ' ( 3 ) = - V ( 1 ) G o V C ( 1 ) - V ' ^ G J ^ V W - V ( 1 ) G £ V ' C ( 2 ) - V W G ^ V W G J ^ V W (D.51) For this first equation, I write V ' ( ) in the same coordinates as V " W i n ( D . 4 8 ) and using 2 GJ* from ( C . 1 3 ) , I obtain pp\kg,-V2 \k , v , UJ) V s\k , -u V 9 s { 23 P g \k ,T] ; 2g 3 UJ) 2S sp\kg, -V2g\k , v \UJ) V s\k , -i] \k ,r] ; UJ) V s to J_< { 2a dk s VpP (kg, RF /2iri2 RT /2iu RTs/^V2 P 2 2 2g s ~V2g\k, 2a v , UJ), V^J (k ,-u \k, 2 g 2g 7/ ; UJ) 2 SP (k , -V2 \k, v ; UJ) VQ (k , -r] \k, r) ; UJ) \ V R% /2iu P g g 9 V $ s 2 g (k, -v \k , u ; 2 s 2s [V£}(k,-Ti2\k. v ;u) t 2t UJ), 2g 2 Vp^J (k, -u \k ,rj ;uj) 2 3 2s (D.52) V^]) (k, —rj \k ,rj ;uj) 2 3 2s where I have used the spatial integrals as Fourier transforms. This equation can be written in the compact form CO dkV^(k , k-uj)G™(k-a,)V«(fc, g / •CO k ;uj) 3 (D.53) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 168 where R /2iv 2TT SP (D.54) R /2i 2 ss V2 A similar process can be performed for all the terms in the subseries leading to the recursive relation CO / dkV'W(k , k; uj)G^(k; g k ; UJ) . uj)V^(k, (D.55) s •oo Data reconstruction and multiple removal T h e next step is to reconstruct data using V and directly propagating reference Green's functions. This gives the data series D' = G * V ' G * = G ^ ( V W + V'( ) + V'( ) + = D(I) + 2 f)'(2) + D 3 ' ( 3 ) (D.56) . . . ) G ^ . . . + T h e first term is just the 4 component (P, S) data i n (D.49). T h e second term is D'^\k , z -UJ) = R (k , z \k„ g g s g z ; UJ) V'^(k , g g g z ; UJ) k ;u>) R (k , s 3 3 (D.57) 3 which after substituting (D.53) and (D.50) becomes £>'( \k ,z \k ,z ;uj) = 2 g g s s oo / z ; UJ) V « ( f c , k; UJ) G^ (k; R (k , g g UJ) V « ( f c , k ; UJ) R {k , c g g 3 3 z ; UJ) dk 3 3 •CO OO D (fe„ z \k, z ; UJ) R ; \ k , z ; UJ) G^ (k; oo UJ) R~\k, c (1) / g 3 3 z ; UJ) D™(k, z \k , 3 3 z ; UJ) dk 3 3 oo / B^(k , g z \k, z ; UJ) J(k; UJ) T>W(k, z \k ,z g 3 3 3 3 UJ) dk (D.58) -oo where Rpp I PPe ( iu2 3(k\u) 7T ''} Zg+Zw Rsp I PSe ( <> * ™ i V2Z +r iZ b) Rp IP Pe ( i c s R-ss/PSe ^ ^ irl Zs+z > } (D.59) r,2Z! +V2Zwb APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 169 W h e n the receiver is located on the water bottom, I use C.14 and the reflection coefficients in A p p e n d i x B to obtain J(k;u)) = —!_ ^(*'-*»») e (PP + Z SS)/(r, 2 PPZ 2 and Z= - k) 2 p v^\(k - ) (1 + Z ) 2 2 where 0 = 1 + 2 2 2 V 2p kv (3 (1 + Z ) 2 lV2 (PP + 2 2 Z SS)/ku 2 (D.60) 2 e ™\ ivi2z Proceeding with the data reconstruction for all the terms i n (D.56), I obtain the recursive relation CO D'*"" ^, 1 / •CO z \k, z ; UJ) J(fc; g 3 UJ) D « ( f c , z \k , z ; u>) dk. (DM) 3 3 s
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- An inverse scattering series method for attenuating...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
An inverse scattering series method for attenuating elastic multiples from multicomponent land and ocean… Matson, Kenneth Howell 1997
pdf
Page Metadata
Item Metadata
Title | An inverse scattering series method for attenuating elastic multiples from multicomponent land and ocean bottom seismic data |
Creator |
Matson, Kenneth Howell |
Date Issued | 1997 |
Description | A method exists for marine seismic data which removes all orders of free surface multiples and suppresses all orders of internal multiples while leaving primaries intact. This method is based on the inverse scattering series and makes no assumptions about the subsurface earth model. The marine algorithm assumes that the sources and receivers are located in the water column. In the context of land and ocean bottom data, the sources and receivers are located on or in an elastic medium. This opens up the possibility of recording multicomponent seismic data. Because both compressional (P) and shear (S) primaries are recorded in multicomponent data, it has the potential for providing a more complete picture of the subsurface. Coupled with the benefits of the P and S primaries are a complex set of elastic free surface and internal multiples. In this thesis, I develop an inverse scattering series method to attenuate these elastic multiples from multicomponent land and ocean bottom data. For land data, this method removes elastic free surface multiples. For ocean bottom data, multiples associated with the top and bottom of the water column are removed. Internal multiples are strongly attenuated for both data types. In common with the marine formulation, this method makes no assumptions about the earth below the sources and receivers, and does not affect primaries. The latter property is important for amplitude variation with offset analysis (AVO). The theory for multiple attentuation requires four component (two source, two receiver) data, a known near surface or water bottom, near offsets, and a known source wavelet. Tests on synthetic data indicate that this method is still effective using data with less than four components and is robust with respect to errors in estimating the near surface or ocean bottom properties. |
Extent | 19696269 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-03-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0087987 |
URI | http://hdl.handle.net/2429/6665 |
Degree |
Doctor of Philosophy - PhD |
Program |
Environmental Sciences |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1997-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1997-196178.pdf [ 18.78MB ]
- Metadata
- JSON: 831-1.0087987.json
- JSON-LD: 831-1.0087987-ld.json
- RDF/XML (Pretty): 831-1.0087987-rdf.xml
- RDF/JSON: 831-1.0087987-rdf.json
- Turtle: 831-1.0087987-turtle.txt
- N-Triples: 831-1.0087987-rdf-ntriples.txt
- Original Record: 831-1.0087987-source.json
- Full Text
- 831-1.0087987-fulltext.txt
- Citation
- 831-1.0087987.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.831.1-0087987/manifest