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 E L A S T I C M U L T I P L E S F R O M M U L T I C O M P O N E N T 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), The 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 Apri l 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 Earth and Ocean Sciences The University of British Columbia 129-2219 Main Mall 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. 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 ( A V O ) . 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. 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 The Born series 10 2.3 Two dimensional scattering from one dimensional models 14 2.3.1 Scattering by a half space 16 2.3.2 Scattering by a layer over a half space 24 2.4 Conclusions 31 3 Free Surface Multiple Removal for Land Data 34 3.1 Introduction 34 3.2 Forward and inverse elastic wave scattering 34 3.2.1 Homogeneous isotropic reference medium and (P,S) potentials . . . 41 3.3 Free surface (P,S) Green's function 47 in 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 4 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 Removing the reference wavefield 95 4.5 Synthetic data examples 96 4.5.1 Model 5 96 4.5.2 Model 6 102 5 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 Synthetic data example 120 5.4 Discussion and conclusions 121 6 Conclusions 129 Bibliography 134 A Calculated series for 1-D Born scattering 140 iv B Elastic reflection and transmission coefficients 144 C (P,S) decomposition and up-down wave separation 148 C . l Land 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 Data 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 16 2.3 Comparison of the true reflection coefficient to the Born approximate reflec-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 . . . 25 2.8 Order of scattering terms vs. event type 28 2.9 Born series first primary for a plane incident wave 29 2.10 Born series second primary 30 2.11 Born series first internal multiple 31 3.1 Scattering in displacements (ux, uz) 44 3.2 Scattering of ((f>, tp) potentials 45 3.3 Reflected and converted waves for a P wave incident on a free surface . . . 49 3.4 Reflected and converted waves for an S wave incident on a free surface . . 50 3.5 Flowchart of steps in elastic free surface multiple removal 56 3.6 Subevent construction of free surface multiples 58 3.7 Land data Model 3 64 vn 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 Model 3 (P,S) data after multiple removal 68 3.12 Model 3 (P, S) data with velocity error 69 3.13 Model 3 (P, S) data after multiple removal with velocity error 70 3.14 Model 3 (P, S) data with missing components 71 3.15 Model 3 (P,S) data after multiple removal with missing components . . . . 72 3.16 Land data Model 4 73 3.17 Model 4 displacement data, shot gather 446 74 3.18 Model 4 shot gather 446 after (P,S) decomposition 75 3.19 Model 4 shot gather 446 (P, S) data after 2 passes of demultiple 76 3.20 Model 3 displacement data after multiple removal 78 3.21 Model 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 Model 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 Model 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 121 V l l l 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 Bp1], and 124 5.8 Elastic internal multiple suppression for B^l and B ^ 125 ix 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 Dr. Tad Ulrych for supporting me during my years at U B C . This support was financial in 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 Dr. Arthur 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. Dr. Weglein also taught me much outside of science and I am thankful for the many lessons I gleaned from him. Third , I want to thank my committee members Dr. Bruce Buffett and Dr. 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 Graham 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, David Butler, Anton 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 my family who, despite my frequent doubts, always believed that I could do it. The 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. Her support enabled me to devote myself to this work. xi 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. The 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. The 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. INTRODUCTION 2 P r l 0.8 1000 P r 2 F S M 1 I M 1 F S M 2 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. The 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. INTRODUCTION 3 change in 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. The order of the multiple refers to how many times it interacts with the free surface. A n internal multiple has all of its downward reflections below the free surface. The 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 in 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 in 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. As a result, both primaries and multiples are migrated indiscriminately and this leads to confusing or distorted sub-surface 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'Brien 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 mul-tiples in the data. They 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 in shallower portions CHAPTER 1. INTRODUCTION 4 of the earth. Since the velocity generally increases with depth, the multiples have a lower moveout velocity than primaries at the same arrival time. Many methods seek to exploit this differential moveout property and can generally be classified as filtering operations. The second property of multiples that forms the basis of multiple suppression is pe-riodicity. In any seismic experiment, many orders of multiples are generated, whereas primaries only arrive once from each reflector. As 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, Yilmaz (1987)) as it is called in 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 one-dimensional, the velocity is known, or that sufficient differential moveout exists between primaries and multiples. When these assumptions are valid, conventional methods have been and will continue to be effective. When 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. The 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. The 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 Van den Berg (1990) and also by Dragoset and M a c K a y (1993). The method of Verschuur et al. is a surface replacement method based on a Huygen's model of seismic wave propagation. Fokkema and Van 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 in-dependent methods in the forward models which are used to describe seismic data. The latter methods describe seismic waves as propagating with different velocities through the earth and reflecting and transmitting from interfaces where the earth properties change. The 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 in 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 inter-face. Although this method has its strengths; it requires the velocity profile in 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. The 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. INTRODUCTION 6 inversion. It turns out that the full series diverges for all but a small range of earth mod-els (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. The method of Carvalho et al. (1992) assumes that the sources and receivers are located in the water and that the multiples are generated by an acoustic free surface above the sources and receivers. When 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). The 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. By 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 in 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. INTRODUCTION 7 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 in 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. The 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 in 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. In an elastic earth, there are two modes of seismic wave propagation, P-waves and S-waves. These waves respond to the properties of the subsurface in different ways. By 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. The work in 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. The 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. The 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 medium and a perturbation about that reference medium. In order to keep the perturbation (the un-known) 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. The 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 SERIES 9 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 Born 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. The different terms in the perturbation series correspond to the number of scattering inter-actions 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. The 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 about the earth model. 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. The 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 SERIES 11 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. The second section of this chapter gives a brief introduction to the forward scattering Born series. In Sections 2 .3 .1 and 2 .3 .2 , I develop analytic Born series solutions for plane waves incident on simple one dimensional earth models. 2.2 The Born series The 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. The 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 s;f) represents the pressure at point r = (%,z) a n d time t due to a source at point r3 and time t — 0. The waves are emitted by a line source which effectively makes this a two dimensional problem. Throughout this thesis, I will be freely forward and inverse Fourier transforming over time and space variables. To 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 CHAPTER 2. FORWARD SCATTERING SERIES 12 k with subscripts indicating the conjugate spatial variable (i.e. ks —> x3). Time will be denoted by t and w is the angular temporal frequency. Hence, Fourier transforming P(x, z\x3, za; t) over time and horizontal shot and receiver locations gives P(k, z\ks, z3; UJ). 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 = i [ l - a ( r ) ] (2.2) c-2(r) 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 / , u>2\ UJ2 V 2 + -T) P(T\TS; UJ) = S(T - rs) + ^a ( r ) P(r\r3; UJ) (2.4) To find a solution to (2.4), consider the causal free space Green's function Po(r|rs; UJ) which satisfies the equation (V 2 + A;02)Po(r|r,;a;) = ^ ( r - r , ) . (2.5) where k0 = UJ/C0. Although k0 is not actually a spatial transform variable, UJ/CQ has the proper dimensions of a wavenumber. In fact, k0 can be thought of as the spatial transform variable in the direction of wave propagation in the reference medium. Using Po(r|rs; UJ) as a reference wavefield, an integral equation corresponding to (2.4) and its physical boundary conditions is (Weglein, 1985) /oo P0(r|r'; UJ) k20a(r') P(r'\rs; UJ) dr'. (2.6) oo This is called the Lippmann-Schwinger equation and plays a pivotal role in scattering theory. It represents the wavefield in an inhomogeneous medium as the sum of the wave-field in a reference medium and an integral that represents the scattered field due to a CHAPTER 2. FORWARD SCATTERING SERIES 13 perturbation. Iterating the Lippmann-Schwinger equation back into itself generates the Born series /CO P o ( r | r » k20a(r') P0(r'|rs; LJ) dr' OO + r P0(r\r';u>)k20a(r') ( f°° P0(r'\r"; UJ) k20a(r") P0(r»\z„ za;«,) dr") dr'+ .. =P0 + P i + P2 +.... (2.7) Throughout this chapter, I will refer to the terms in this series either according to their order in a, or their place in the sequence. For example, the zeroth order term, P 0 , is the first term in the series. Similarly, the first order term P i , is the second term in the series and so on for the other terms. To facilitate a physical interpretation of the Born series, consider the scattering per-turbation to be composed of a region of closely spaced point scatterers separated by the reference medium. The 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 rs, directly to the measurement point at r. The second term in the Born series contains k^a^r') sandwiched between two Green's functions. The Green's function on the right represents a wave which propagates from the source at r, to a point scatterer at r'. The wave undergoes an interaction with the scatterer at this point and then propagates back to r via the Green's function on the left. The 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 refer-ence 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 inter-actions 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 SERIES 14 velocities in different materials. 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 Born series after the second term leads to the Born approximation. In this approximation, the data (the measured scattered field) is assumed to be linear in the model (the scatterer). When the model is small, the higher order terms in the Born series become less important and the Born 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. When 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 Born 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 scat-tering 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 dimen-sional, I do believe there is value in studying what Arthur 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. The reference medium can then be chosen as a homogeneous whole space composed of the material in the region outside of the scatterers. The 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 CHAPTER 2. FORWARD SCATTERING SERIES 15 is written using the Weyl integral 1 f°° Qik'iX9~Xs) Q^'\Z9~Zs\ G0(xg,zg\xs,zs;uj) = — / — — dks (2.8) ^ J-oo twos where u03 = (UJ2 jc\ — k2)1^2, and c 0 is the wave velocity in the reference medium. Through-out this thesis, v will be used to represent the vertical wavenumber for P wave propagation. The subscripts for v will indicate the velocity c and wavenumber k which are related to v. Rewriting (2.8) as Go(xg,Zg\x3,zs]uj) = — / — o(xg,zg\ks,zs;uj)dks (2.9) ^ J-oo two* where Mxgiz8\k.tzt;u>) = e ^ + ^ — D , (2.10) it is evident that the 2-D Green's function is a superposition of weighted plane wave solu-tions. 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(xg,Zg\x3,zs;uj) = Go(xg,zg\x3,z3;uj) /CO CO , \ G0(xg, zg\x', z'\UJ) V(X', z') G0(X', z'\xa, Z3]OJ) dx'dz' •OO J — OO /CO f OO / G0(xg,zg\x',z';u;)V(x',z'). • CO J — OO / GQ(x', z'\x", z"\ UJ) V(X", z") G0(X", z"\xs, Z3; UJ) dx"dz" + ... CO J — oo (2.11) where V(x,z) = k^a(x,z). For each term in the series, Fourier transform with respect to CHAPTER 2. FORWARD SCATTERING SERIES 16 x3 and multiply by 2IVQS. This gives P(xg,zg\ka,z3;uj) = 0(xg,zg\kt,zs]u) /OO f* OO / G0(xg, zg\x', z'\ UJ) V(x', Z') 4>0(X', z'\k3, z3; UJ) dx'dz' oo J—oo poo poo + I I G0(xg,Zg\xl,z';uj)V(x',z')-G0{x', z'\x\ z"; UJ) V(X", z") o{xg,zg\k„z.;u>) = (2.13) where k3 is the horizontal wavenumber and uo3 is the vertical wavenumber. Different values of k3 then correspond to incident wavefields with different incident angles. To reconstruct the series for a line source, simply multiply each term by 2iv0s and inverse Fourier transform with respect to x3. 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 in the following examples since the source is always located above the scatterers. 2.3.1 Scattering by a half space The simplest case to consider is scattering from the half space as shown in Figure 2.2. I will refer to this as Model 1. The interface between the two media is at z\ and, for simplicity, the source is located at zs = 0. The scattering perturbation for this model is V(x\ z') = V(z') = k2oaiH(z' - Zl) (2.14) where aj = 1 — c\jc\, k0 = UJ/CQ, CI is the wave velocity in the half space, Co is the reference 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. CHAPTER 2. FORWARD SCATTERING SERIES 17 incident wave x z ^ measurement v point (xg,zg) Zl C l Figure 2.2: Model 1: acoustic plane waves incident upon planar interface. The first term in the Born series is just the incident wavefield P0(xg, zg\k; UJ) = e ^ s + ^ s ) . (2.15) For brevity, I have dropped the zs = 0 and will use just k for ks. Keep in mind the difference between k and k0 and that v% = k^ — k2. From (2.12), the second term is ^ /»oo />oo /.oo eikg(xg-x')eivo9\zg-z'\ P1(Xg,Zg\k;UJ) = — / / / — dkg 2TT J_oa 2tu0g k20 ax H(z' - Zl) ei(kx'+'**r> dx'dz' (2.16) 1 f ° ° /"oo POO 7,2 = — / / / e i ( f c- f e 9)x' eikgx9 ei^g\zg-z'\ eiv0z> _o_j_ dz'dx'dkg. 2TT J_00 J_x JZ1 2iu0g Using the Fourier transform /oo eKk-kgydx> = 27r5(kg - k), (2.17) •oo (2.16) becomes / £ ( & s - j f e y ^ e * * . ! * . - * ' ! e - ° 2 ' A; and fog —• giving eikxgeiv0\z>-zg\ e i u 0 z ' d z > ^ CHAPTER 2. FORWARD SCATTERING SERIES 18 When the receiver is on the source side of the interface, zg < z\, and 1 J - " (2-20) = e , f c e'lfe£ aj e_iM,*« / e**22' dz'. 2^o J Z 1 This integral can be solved by considering a small amount of dissipation in the wave prop-agation. 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 0 ( l + ie) where e > 0 for w > 0 (Aki and Richards, 1980). W i t h this modification, (2.20) can be solved directly giving Piix^zg < Zl\k-uj) = ~k2oai e**9eiM^-ze). (2.21) 4z/0 This can also be obtained by a contour integration in the complex z' plane. Adding Po and P i gives the Born approximation PBorniXg, ZG < Zl\k] u) = j ^ ' + w ) + - L * £ 0 l e ^ x J f > o ( 2 Z a - z 3 ) _ ( 2 . 22 ) In this expression, the first term is a wave which propagates outward from the source directly to xg, zg. The second term describes a wave which is reflected from the boundary and propagates with the reference medium velocity. The 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)2 = u^ + k2 where OJ/CQ can be thought of as the wavenumber in the direction of propagation and u0 and k are its vertical and horizontal components. These components can be expressed as k = UJ/C0 sin(#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(#;nc) = sin(0 r e/z). 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 _ [Pil _ k2ai _ 1 ( C L - C Q ) ( C I + C p ) H B o r n - | P q | _ ^ - 4 c o s ^ ) 2 ^ . [2.26) CHAPTER 2. FORWARD SCATTERING SERIES 19 0.5 r co/ci Figure 2.3: Comparison of the true reflection coefficient to the Born approximate reflection coefficient for normal incidence. The velocities of the first and second media are c 0 and C i respectively. For normal incidence, set 9 = 0 and = i fa-»)(•» + »). ( 2 . 2 4 ) 4 c{ This result was also obtained by Weglein and Gray (1983) in a study of the sensitivity of Born 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>2jc\ — & 2 ) 1 / 2 (DeSanto, 1992). For normal incidence, this is (ci — c 0 ) /(ci +c 0 ) which is clearly not in agreement with (2.24). Figure 2.3 shows, however, that the normal incident Born approximate reflection coefficient does approach the proper value for small values of a0 i.e., when CQ/C\ ~ 1. This is the weak scattering convergence criterion for the forward Born series (see, for example, Prosser (1975)). For transmission data, zg > z\ and the integral in (2.19) must be split into two portions: CHAPTER 2. FORWARD SCATTERING SERIES 20 zg < z' and zg > z'. This gives fzg J.2 r O O jL2 P1(xg,zg>z1\k-uj) = e^ei^'^dz' + e ^ e ^ ' e ~ ^ g ^±dz> JZl ™o JZg ,2 25^ = J _ e ^ e - o , 9 f c 2 a i ( 1 + 2iVo(Zl - Zg)). Combining this with the first term gives the Born approximate transmitted field PBorn(xg,Zg > Zl\k]U) = -e^+^ Ik? Avl ^a1(l + 2iu0(z1 -zg)) - 1 ( 2 . 2 6 ) This wave propagates in the proper sense (increasing x and z directions); however it prop-agates with the reference medium velocity, has the wrong amplitude, and in addition, does not have the proper refraction angle. In the following, I show that by adding more and more terms to the forward Born series, the solutions for the reflected and transmitted waves converge to the familiar expressions. The third term in the series is /oo poo / G0(xg, zg\x', Z';UJ) kl a(x', z')Px(x\ z' > zi|fc;u>) dx'dz •CO w —CO Y />co yco yco eikg(xg-x')eiu0g\zg-z''\ = — / / / k% a 1 P i ( a ; / , z' > Zi\k; UJ) dx'dz'dkg 2TT 7_00 J_X JZ1 2iu0g (2 .27 ) Expressing Px(x',z' > z^fyuj) as elkx'P^x', z' > zx\k\uj) , this gives y>oo £oo poo eikg(xg-x')eiv0g\zg-z'\ -> pOO pOO pOO ikg(Xg-X') lV0g\Zg-Z\ ^ P2(xg,zg\k;uj) = ± / / j-A kla^k*Px(x',z'>zx\k-uj)\dx'dz'dK ^ J-oo J-oo JZl ZtvOg /CO pOO T/i. ^ / e****S(kg - ky^-t'^^-Piix^z' > Zl\k;uj)dz'dkg •co JZ1 2lU0g = / eikxgeiv0\z3-z'\^prx^z> > z^.^dz' JZ1 2™o I 2 /»CO = « o ^ i e i t e , / ^\*,-''\p^x'iZ' > Zl\k;uj)dz' 2wo JZ1 (2 .28 ) which can be split into expressions for reflection or transmission data. CHAPTER 2. FORWARD SCATTERING SERIES 21 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: kn 7 2 /»0O Pn(xg,zg < zx\k;w) = ^ i e * " / e^''^Pn^(x\ z' > Zl\k;uj)dz> ^lv0 J zx (2.29) for reflection and Pn(xg,zg > z.lk-u) = ^±eik*° f" C ^<«•-• '>P n _ 1 (x^* , > Zl\k;u;)dz' 2*IVQ JZl fiOO / e i M o ( z ' - Z f l ) P „ - i > zt\k;uj)dz' J Zn k2 n, K0 a l ik e (2.30) 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(xg, zg > Zl\k; UJ) = P i + P 2 + Ps + • • • eikxgeiu0(2z1-zg) laik20 t la\k* | 5 a\kl | 7 qjfcg 4 vl + x 8 v% 64 vi + 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. The reflection coefficient can be obtained using U vl > 2v2 %{vl> l§Kvl> 128 1 vl ' 25Q{vl)^---' (2.32) which leads to the expected result 1 - 4 / 1 -aikl Vn 1 = (2.33) CHAPTER 2. FORWARD SCATTERING SERIES 22 Figure 2.4: For a plane wave incident on a half space with incident angle 0;n c and unit amplitude, the sum of the scattering terms results in a reflected plane wave with reflection angle #;nc and angle dependent reflection coefficient R(9inc). where, as before, v0 = ^u2jc\ - k2 and vx = voy/l - UJ2ax jc20vl = AA> 2 / C I - k2. Figure 2.4 shows the scattered wavefield that results from summing the Born series. For the transmitted wave, the series yields P(xg, zg > Zl\k;uj) = 4>0 + Pi + P2 + Pz + 1 + 1 atk2 , 1 a\kA0 , 5 offcg , 7 a\k* 4 v$> iv0(zg - Z i ) + 8 v$ 1 a-ikl + 64 ven + 128 1 a\k\ 5 a\k% 0 32 " o ( * s - * i ) o 1 a2k? 3 n3k6 I n a. 1 9 (X^/uQ 8 vl 32 vl 128 zv08 + . . . 7 afog 64 /yn8 + + . . . (2.34) CHAPTER 2. FORWARD SCATTERING SERIES 23 which can be expressed as P(xg,zg > zx\k-u) = eikx»eiv^(l + R) 1 + - V0)(zg - ZX) - ~{V\ ~ V0)2(Zg - ^if = (1 + R)eikXg e'VlZfl e 2'( i /i- i /o)(zfl-zij The transmission coefficient, T, is given by + (2.35) 1 -2i/i z'l + v-i (2.36) as expected. In addition, this wave propagates with the proper medium velocity (c\) and has the correct refraction angle. The latter follows from writing k as k = u>/cinc sin(0jn c) — uj/crefr sin(# r e/ r). From this the familiar Snell's law s i n ( ^ n c ) / c j n c = sm(Orefr)/crefr 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). The ratio test can be used to show that this is satisfied for |aifc§/i/Q| < 1. Using u0 = k0cos(8) where 8 is the angle between the direction of propagation and the normal to the interface (see page 18), convergence is ensured for sin(0) < c 0 / c i < ( 1 + cos 2(^)) 1/ 2 . Figure 2.6 shows the region of convergence for a given incident angle. The curve 1 — (crj/ci ) 2 as a function of (cj/co) must remain between the lines dotted lines given by ± cos2(#). With increasing angle, this region 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 \ /2ci . By 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; n c and unit amp-litude, the sum of the scattering terms results in a transmitted plane wave with angle 6refr and angle dependent transmission coefficient T(6inc). The 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 Born series alter the amplitude and adjust the propagation velocity of the scattered wavefield. The 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. Model 2 consists of a single layer with velocity c i on a half-space with velocity c2. As before, the reference 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) + (a2 — a\) H(z — z2) where a\ = 1 — c\jc\, a 2 = 1 — C Q / C | , and Z\ and z2 are the locations of the first and second interfaces respectively. As before, the source is located at zs = 0. 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 CHAPTER 2. FORWARD SCATTERING SERIES 25 Figure 2.6: The Born series for scattering from a half space converges when the curve 1 _ ( c 0 / c i ) 2 is in the region between ± c o s 2 ( 0 ) where 8 is the plane wave incidence angle. 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 zg < z1 is ikxa iun(2z1-zg)kla^ , iu0(2z2-zg) kl(a2 - ai) (2.37) 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 CHAPTER 2. FORWARD SCATTERING SERIES 26 incident wave \ measurement . v point (xg,zg) I C l 2 2 c 2 Figure 2.7: Model 2: acoustic plane waves incident upon a layer over a half space. Meas-urements of incident and scattered waves are made at xg,zg. through the first interface and a wave reflected from the second interface P1(xg, zi P M^ e***» dz> + H M°I E*>(2*'-*,) D Z > JZ1 2luo JZg 2iu0 f°° L% n I 4zv0 4i/Q J (2.38) For zg > z2, the solutions are transmitted waves through both interfaces. P1(xg,z2) = eikx> H dz' + f ' M%*>*« o V + f~ ^jM*'1-*.) dz> (2.39) To first order, the Born 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 SERIES 27 Starting with equations (2.37)-(2.39), the following recursive relations can be derived to calculate successively higher order terms in the Born series: r M i^e-o(,'-,9)pn_i(zi < Z2\k-u)dz' Pn(xg,zg < zi\k;u>) = elkx-+ r ^eiMz'-Z3)Pn-i(z2 < z'\k-w)dz' Jz7 2™o coo r 2 (2.40) + f22 tf^jMz'-^p^zi < z' < z2\k)u>)dz' Jza 2 i v ° r !&ljM*'-*.)pn_1(Z2 < z'\k;u)dz' (2.41) Pn(xg,z2 < zg\k;uj) = e f2 ^e^'-^Pn-iizi <*< z2\k;uj)dz' JZl 2luo + r ^e**<* ' - " '>P„- i (*2 < z'\k;Lo)dz' (2.42) Jz2 2zz/0 POO 7,2 + / ^ - ( ^ ^ ( z ^ z ' l ^ u ; ) ^ ' B y computing a number of terms in 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 < Z l | f c ; U > ) = P0(Xg,Zg < Zl\k\ Uj) + PPTl(Xg,Zg < Z^Uj) + P^ (xg , Zg < z\ | fc J LO) + Pndtl(xg, zg < zx\k\Lo) + P^ixg, Zg < z^k]u) + 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 SERIES 28 Since each term describes a different order of scattering interactions, there is a relation-ship between the terms in the Born 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 sec-ond interface and one from the first. Since this type of path contains three changes in propagation direction, the first order contribution to the internal multiple is contained in the third order term of the Born 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 in 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. The first primary reflection is 7 /Afc.V , 21 (h\a^ + 2 k ( T ) + i i a ^ 1 + " ) ' ( 2 ' 4 4 ) Similar to the previous example, using (2.32) leads to the solution Ppil(xg,zg < Zl\k]uj) = Rieik*geiM2*i-*9) (2.45) where R i = [2 - 2(1 - fc02aiK2)1/2 ~ fcpaiK2] = "o ~ "i ( 2 4 6 ) is the reflection coefficient corresponding to the first interface as shown in Figure 2.9. The second primary reflection is represented by a more complicated expression involving terms in ko, VQ, Z2, Z\, a-i and a2. When certain series expansions are used, these terms reduce to Ppi2(xg, zg < Zl\k-uj) = TlT[R2eikxoeil/^-z^eiv^-z^ (2.47) CHAPTER 2. FORWARD SCATTERING SERIES 29 event direct wave 1st primary 2nd primary 1st order multiple 2nd order multiple number of scattering interactions 1 2 3 4 5 6 1 2 3 4 5 6 3 4 5 6 5 6 Figure 2.8: Events in the reflected wavefield are shown on the vertical axis vs. the different Born scattering terms which contribute to the modelling of that event. where, as before, u\ is the vertical wave number in the second medium (see Appendix A ) . The 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. The reflection coefficient R2 relates to the second interface and can be written as 2 - fcgaa/i/0a - klaxlvl - 2(1 - ^ ^ ( 1 - kjai/v*)1'* - u2 K i = n—r~2——T~2 — —;— • The second primary is shown in Figure 2.10. The corresponding portion of the Born series for the second primary can be obtained from (2.47) by first expanding exp[ifi2(;z2 — zi)\ as a power series in iv-\2(z2 — zi), then substitute the expansion for (1 — A^ax/Vo) 1 / 2 into v\, 2"i, and T[. Next, the expansions for (1 - klax/viyl2, (1 - A ; 2 a 2 / ^ o ) 1 / 2 a n d 1/(kla2/^l - tfax/v2) are substituted into (2.48) which is in turn substituted into (2.47). The convergence of the second primary depends entirely on the convergence of the series expansions for (1 — k\a\lv\y^ and R2. The convergence of (1 — k\axfv^)1//2 was discussed CHAPTER 2. FORWARD SCATTERING SERIES 30 C l ; z2 c 2 Figure 2.9: Born series first primary for a plane incident wave. in the previous example and is ensured for sin(#) < c 0 / c i < (1 + c o s 2 ^ ) ) 1 / 2 . For P t 2 , there are two choices for expanding the denominator in (2.48), namely: Kailvl-^la\lvl tfa2/vo(l - a1/a2) (2.49) and , 2 , / \ • (2-50) K^/^oi1 - « 2 / a i ) These converge for |a2| > | \ 2! n! 7 = — Ri(l — i2 2 )J? 2 e i f c X s e i l / 0 ( 2 ( 2 z 2 - z i ) _ z s ) e i s / o 4 ( z 2 - 2 i ) ( 7 i - i ) = - i ? i ( l - / ^ J ^ e ^ e ^ C ^ ' i - ^ J e ^ 4 ^ - * ! ) (2.51) where 7 i = (1 — k^ai/vl)1/2. From Figure 2.11, it is clear that this expression properly describes the multiple with respect to both the amplitude and the phase. For the second multiple, the terms reduce to PTnlt2(x; k) = R\(l - RDRle^e^^i-^)^!^-^) (2.52) As before, the corresponding portion of the Born series is obtained by expanding out exp[i/Vi6(z2 — zi)], Ri, and R2. A l l of the multiples, including the higher orders not shown here, contain terms in Ri, R2, and exp(ifo(l — k^ai/vl)1!2). Since the primaries are also dependent on the CHAPTER 2. FORWARD SCATTERING SERIES 32 Z2 c 2 Figure 2.11: Born 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. Ultimately, given that sin(#) < C o / c i < (1 + cos 2(^)) 1 / ' 2 and sin(0) < c 0 / c 2 < (1 + c o s 2 ^ ) ) 1 / 2 , the Born series solution for this problem converges. For each velocity and a given angle, a curve similar to Figure 2.6 can be drawn. The curve of 1 — ( c 0 / c 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 The work in 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. The latter represents the way seismolo-gists traditionally think of events as being generated from wave propagation and reflections from interfaces. The 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-CHAPTER 2. FORWARD SCATTERING SERIES 33 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 Born series describes the data, the higher order terms in the series must describe certain data characteristics beyond the Born ap-proximation (see, for example, Stolt and Jacobs (1981)). The calculations in 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 c0 is the reference velocity and cn is the velocity of the n'th layer, the Born series for a plane wave with incident angle 8 converges if c n sin(#) < ci < c n ( l + cos2(#))2. This implies that models with very low velocity layers will lead to divergent Born series which in turn has implications for inverse scattering. In the inverse series originated by Jost and Kohn (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 Born 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. The 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. The 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 The mathematics of elastic wave propagation involve various matrix or tensor quantities. When 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. The notation convention I use is as follows: • Lower case bold (normal and Greek) letters represent vectors with individual compo-nents in lower case italics. I will alternatively use indicial notation or explicitly write column vectors enclosed in brackets. For example, the displacement vector can be r i T written as u, Ui, or ux uy uz • Upper case bold refers to matrix quantities. Individual elements are denoted using in-dices, e.g., Fxx 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(rp|rs; u;). In this thesis, operators are integral operators which involve their associated func-tions. 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. The meaning of the quantities in the above examples will be made clearer as the theory develops. The first step is to derive the scattering equations appropriate for elastic wave prop-agation. 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 - (cijkiukj),j = fi i = 1,2,3 (3.1) where Ui is the displacement of the continuum in the i ' th 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. The directions 1,2,3 are x,y,z respectively. The solution for the displacements at a geophone located at rg due to a body force distribution /m(r,u>) is (see, for example, Morse and Feshbach (1953)) Ui(rg;uj)= / gim(Tg\x'])(Fmj(r'))dr'. (3.7) oo One point of clarification is required. What is the difference between <7V and ^(r)V(r)? The first is Q operating on V which involves an integral over space of the matrix product of their functions Q(r) and V(r). The 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). When 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(rg|r'; UJ) are 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. As 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. The wavefield in the background or reference medium satisfies the operator equation C0g0 = X. (3.8) where Go is the reference Green's operator and C0 is the wave operator for the reference medium. The 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 C0Q=X + Vg. (3.9) Next, I employ Go to invert C0 and arrive at the elastic Lippmann-Schwinger equation G = G0 + G0VG (3.10) or written in terms of functions /oo poo / Goij(rg\r'; u/)Vjfc(r'|r"; UJ) £ f c m ( r " | r 3 ; UJ) CLT'CLT" CO J — oo (3.11) where ( / i m ( r j r , ; u;), Goim(*g\rs',<*>)) and Vjk(r'\r";UJ) 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. The form of V and its im-plications is discussed in more detail in Appendix D and has also been discussed by Stolt and Weglein (1985). Although it is possible that a perturbation with more than four vari-ables 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 wave-field in the actual medium to the wavefield in a reference medium plus a perturbation wavefield which depends on Go, V, and G itself. To isolate G from the other quantities, the Lippmann-Schwinger equation can be iterated to obtain the Born series G = 90 + Q0VQo + 9oV9oVGo + . . • (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> = g0VQQ + Q(S?QoVQo +•• • (3.13) The 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. Of course, using this series to do forward modelling would be terribly inefficient. However, as was also shown in Chapter 2, certain insights can be gained by studying the forward modelling process. What is ultimately of interest to geophysicists is not the forward problem but the inverse problem (although the two are intrinsically related). The 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. The inverse scattering method I employ was originated by Jost and Kohn (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 in more than one dimension (Stolt and Jacobs, 1981; Weglein, 1985). A formulation of the Jost and Kohn 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 U + . . . (3.14) where the n'th order term is n'th order in T>. A simple example illustrates why this approach is reasonable. Suppose A is related to B by A = 75/(1 - B) = B + B2 + B3 + ... CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 40 in analogy with the forward scattering series. The inverse relationship for B given A is the series B - A/(l + A) = A - A2 + As + . . . . 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> = G0V{1)G0 (3.i5) 0 = g0V{2)Go + GoV(1)GoV(1)9o (3.16) 0 = GoVi3)Go + . 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. The 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 at depth. After interacting with \ 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 41 reflection variations. Since sandwiching ^ by Q0 implies upward continuing the sources and receivers to the surface, inverting these operators can be thought of as downward continuation. This can be thought of as migration using the reference medium velocity. Thus, ^ can be thought of as a reflectivity image in scaled time which is obtained by " (2) * (3) migrating the data using the reference medium velocity. The role of the terms V , 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. The 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 in 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) ~ (3) contained in the terms V , 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) poten-tials The question arises as to what to choose for the particular reference medium. The 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 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 42 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), w here g0 = n ^ G o i i (3.18) r - l _ o — l/7o 0 G 0 = GQP 0 0 l//*o 0 Gos (3.19) n dx dz -dz dx n 1 = v dx dz -dz (3.20) 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 - 2 is an integral operator (i.e. if V 2 J 4 = B, then A = V~2B where A is the solution to Poisson's equation). The Green's functions Gop and Gos satisfy the scalar wave equations (v"2 + ^0 G0p(xg,zg\x3,zs]Uj) = S(xg - xs)8(zg - z.), ( V 2 + ^ 2 ) Gos{xg, zg\xs, z3; UJ) = 8(xg - xs)S(zg - z.) (3.21) where a0 = (lo/p)1^2 is the reference P wave velocity and /30 = (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 o " 1 G o contains the independent eigenvalues of the operator Go and II contains the corresponding eigenvectors. The eigenvalues of a system represent its modes, which, in 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 am 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)) u = V

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

dz dx (3.23) Now consider an experiment in the reference medium where a body force f gives rise to P and S waves. The displacements for this experiment are u = Q0f = n ^ r j ^ G o I I f . This can also be written as - 1 u = pu>' - IF G 0 n f (3.24) using (3.21). Comparing this with (3.23), two things can be concluded. First, the Helmholtz potentials are related to G 0IIf by / o u r ^GoIIf . (3.25) Second, II acts to decompose f into its P and S components [, tp), G 0 : represents propagation of (, rp) from source to measurement point, II - 1 : transforms (, VO to (ux,uz) and the outgoing waves from (ux,uz) to (d>,ip). This process has been compressed so that V^ 1 ) simply scatters P and S waves. The result is the simplified scattering diagram in Figure 3.2. Each element of V^ 1 ) has a specific job to do in the scattering process. To understand what these elements do, consider the displacement field due to body force f. Performing the receiver decomposition, (3.5) becomes = r 0 n u jo(dxux + dzuz) f*o(dxuz - dzux) G 0 I I f + G o V W G 0 n f . = G f (3.30) CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 46 Now suppose that f is such that it gives rise to P waves only i.e., I l f = 8 10 . Then GOP 0 + GOP 0 0 Gos T / ( l ) Vpp (1) V. SP V ( i ) PS (1) ss GOP 0 so that (3.31) 4> — GOP + GOP Vpp GQP , (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. The 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 in 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 Born approximation. Applying the source and receiver decomposition to the entire inverse series gives V = V ( 1 ) + V < 2 > + V < 3 ) + . . . (3.34) from which it follows that D = G o V ^ G o 0 = G 0 V ( 2 ) G 0 + G o V ( 1 ) G 0 V ( 1 ) G o 0 = G 0 V ( 3 ) G 0 + G 0 V ( 2 ) G 0 V ( 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 in terms of (P, S) potentials. The reader may question the value of this last result over equations (3.15)-(3.17) con-sidering the effort expended. The point to writing the inverse scattering series in (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 Land 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 displace-ment 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). The solution for Gop (see, for example, DeSanto (1992)) is written as the Weyl integral 1 t°° 1 GA0P{xg)zg\x3,z3]u:) = —J _ c * ( - . - . ) c H « . - . l d f c (3.36) where v = sign(a>)(u;2/a:2 — k2)1/2 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. The sign(u;) ensures that v has the same sign as u> so that the wave propagates forward in time and out from the source. This condition is also known as the Sommerfeld radiation condition. The 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 e%k(xg-x')ew\z'-zg\ 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 1 f°° Gop(xg,zg\xs,z3]uj) = — / elk»x»G^p(kg,zg\xa,za]Uf)dkg (3.37) CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 49 where G k f o , = J - e - * . - . e - . l ' . - * . l (3.38) 2ivg and vg = v(kg). Alternatively, GQP can be represented by an inverse Fourier transform over the shot domain 1 f°° G"op(xg,zg\xsizs;v) = — J e~tk3XsGop(xg,zg\ks,zs;uj)dks (3.39) where G$P{xg,zg\k.,z.;u>) = J - e ^ e - ' " " - ^ ! (3.40) and us = v(ks). This is a powerful formulation which can be used to build Green's functions in more complicated reference media. To 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. The z axis increases downward. Suppose there is a source of plane waves which propagate upward towards the free surface with wavenumber k3. The plane wave which propagates directly from the source and receiver is eiksxgeii/a(zs~zg) (3-41) When this wave strikes the free surface, it creates a reflected P wave of the form ppeikgXgeiva{zs+zg) (3-42) where PP is the free surface P-P reflection coefficient (see Appendix B for the free surface reflection coefficients). The 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;eikaxgei(vszs+r)sZg) (3.43) where rjs = sign(u>)(u;2//32 — A: 2 ) 1 / 2 is the S wave vertical wavenumber. The inclusion of rjs means that the reflected converted wave propagates with the S wave velocity after reflection. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 50 Free Surface (z=0) PS exp(i(ksx, + vazs + V»zg)) 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/2ivs and inverse Fourier transformed with respect to ks. Since each reflected plane wave solution satisfies the free surface boundary condition, so does the weighted superposition of these waves. Hence, scaling the reflected waves by l/2ivs and inverse Fourier transforming with respect to ks builds the reflected wavefields due to a fine source in the presence of a free surface. These solutions can be written in vector format as / - i d ^OP + ^OPP /orfs ^OSP 0 (3.44) where 1 poo ph G%PP(xg,zB\*.,z.;u) = ^ J — e * ' < " ' - - ' > e * ' < * ' + " « > < f c . (3.45) and 1 poo p b Gf0sSP(xg,zg\xa,z3^) = —J — c * . ^ - . ) ^ . . . ^ . . . ) ^ . (3.46) These solutions are also valid when the receivers are below the source location. The su-perscript 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 in an elastic half space. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 51 Free Surface (z=0) exp(i(ksxs + r]s(zs - zg)) ^ incident S wave / \ \ ^ SPexp(i(fc„z3 + r)szs + v,zg)) converted P wave / depth (z) 1 / SS exp(i(ktx, + -q,{z, + zg)) 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 k3 incident on the free surface as in Figure 3.4. This wave is written as &ikaxg eirj3(zs-zg) (g ^ ) The interaction with the free surface creates a reflected S wave Sgeik,xgeir,,(zs+zg) (3_4g) and a reflected converted P wave §peikaxgei(r)az,+vszg) (3.49) where SS and SP are the S-S and S-P elastic reflection coefficients. Scaling the incident wave by l/2ir]s gives Gg0(xg, zg\ka, ZS;UJ). This is then Fourier transformed over ka to obtain the directly propagating S wave due to a line source measured at the geophone i.e. ^ 0 5 ( ^ 9 ' zg\xs, Z*\UJ). Explicitly 1 t°° 1 G*s{xg,zg\x.,z.;u) = — / —eik>(*°-*^-\'>-''\dk. (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{*SP(xg,zg\xSlzs;uj) = - / _ e * . ( - . - - 0 c * > . « . + * . . ) ^ (3.51) l-K J_00 2lVs CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 52 and GQSSS(X9I zg\x ss 2irjs (3.52) Written in matrix form the incident and reflected waves are 0 + /~tfs ^OPS /~ PSeiT>z° SSeir)Z<> 2iv 2ir\ \dk. (3.55) This is the P and S wave Green's function for a purely P or S wave source in the presence of a free surface. The free surface Green's functions can also written in terms of operators as G o = 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. The key to the identification of CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 53 this subseries lies in the special form of the (P, S) free surface Green's function obtained in the previous section. Note that the operator Go has two portions. The first part, G Q , 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. The second portion, G 0 S , exists due to the presence of a free surface in 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 Go did not contain G 0 S , the forward series would create data absent of any free surface effects. This is in 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 in 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'<2> + V ' W + ...) + V 2 ) + V 3 ) + . . . (3.57) V V ' V ' where terms in brackets are the free surface multiple removal subseries V ' and the remain-ing terms perform the other tasks of inversion (converting from time to depth, mapping reflectivity to mechanical earth properties, and removing internal multiples). The quantity V ' represents a series of terms each of which removes a certain order of free surface mul-tiples when added to V W . As in the previous section, V^ 1 ) is the portion of V linear in 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. The CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 54 terms which perform multiple removal are V'( 2 ) = - V ( 1 ) G 0 S V W V ' ( 3 ) = - V ' ^ G ^ v W - V ^ G 0 S V ' ( 2 ) - V ( 1 ) G * V ( 1 ) G £ V ( 1 ) = V ^ ^ G ^ V ^ ^ G ^ V ^ ) (3.58) V'( n ) = — V ' ^ n - 1 ^ G o V ^ 1 ^ . The recursive nature of these terms is due to the special structure of G Q (see Appendix D )• The next step is to reconstruct data using the new V ' series by sandwiching it between two directly propagating Green's functions. The result is a new series of data where each term has an identifiable role, D ' = GjJV'Go 1 = G 2 ( V W + V'<2) + V'( 3 ) + . . . ) Gtf (3.59) = IV1) + D ' ( 2 ) + D '< 3 > + . . . . The first of these terms is D ( i ) = G J J V W G * = G * ( G * + GlT'MGi + GlT'Gt (3.60) = R g D R s where R g performs receiver (geophone) deghosting and R s performs source deghosting. (See Appendix C for a full description of these operators.) Note that the deghosting operators are applied to the original (P, S) decomposed data. The 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 and in general = - ^ ( G ^ G ^ G ^ D W = - D W H D W . ( G o ) _ 1 D ^ ^ ( G o ) _ 1 (3.61) (3.62) Thus, each term in the multiple removal series is obtained from the original multiple con-taminated data L r 1 ) and an operator H which contains elastic free surface reflection co-efficients and appropriate obliquity factors (see Appendix D) . The 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 ) ' th 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. As a result, it is not necessary in practice to calculate all of the terms in the series; there is little use in predicting and removing multiples which arrive after the end of the recording time. Explicitly, the multiples are calculated in the (kg,ka,u>) domain according to the formula HPP(k,uj) HPS(k,u>) HSp(k,uj) HSs{k,uj) D{£,{k,ks,u) D^s{k,ks,u>) D{pp(k,ks,uj) D%>(k,ka,w) dk. (3.63) Note that the operator H becomes a simple matrix multiply by the function H . See Appen-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>(kg, ks;w) Transform to (P, S) and deghost: D( 1)(^,fc s;o;) Calculate n'th order multiples: D(n)(*S)fc,;w) False Sum terms: D'(kg, ks;u>) 3-D inverse Fourier transform over shot, geophone and time: T>'(xg, xs;t) Output multiple free 4 component (P, S) data: B'(xg, xs;t) 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 dimen-sional earth models. In such a case, the data are invariant with respect to the source and receiver midpoint. The spatial transform variable for the midpoint is kg — ks. Since the data are invariant in midpoint, the only nonzero data values are where kg — k3 = 0 or where CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 57 kg = ks. As a result, the data can be written as T)'(n\k; UJ) 8(k — kg). Substituting this into (3.62) gives D'Wfojw) = -B'^\kg]u,)H(kg-,u,)r)V(kg-,u;) (3.64) which is a simple multiplication in (kg, UJ) space. Since the earth doesn't vary with midpoint, one shot is as good as the next; hence only one shot is needed for multiple removal. It is reasonable to ask why G Q V ' G Q creates data when G Q V G Q does not in the forward series. This can be understood by looking at the multiple removal series term by term. The first term is the data with the source and receiver signatures removed. The remaining terms act to remove the free surface; hence they are the multiples with the sign flipped. Since the multiples and the primaries exist in data space, the whole series can be thought of as a sum of terms in the original data which produces new data containing just primaries and internal multiples. This explains why V is more like V « than V . The 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'(2) is a convolution of D^1) with D^1). A good way to 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 in 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 in 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. When the data are convolved with themselves, these subevents are strung CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 58 together with the free surface reflection coefficient to create the multiple PPPP. A minus sign in (3.61) ensures that when this estimate is added to the data, it cancels the actual multiple. The 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. The sources and receivers are represented by * and V respectively. The need for multicomponent source and receiver data becomes clearer under the subevent construction principle. Consider a first order peg-leg multiple which has an S phase (PPSP) (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 PPPP multiple as in Figure 3.6d). As before, this multiple can be split into two PP 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. Data 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. The 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 re-moval The 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. As outlined in Section 3.2.1, the decomposition is done by calculating the divergence and curl of the dis-placement 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. The 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 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 60 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. What follows is my adaptation of this approach to the elastic scattering problem. I start with the equation for elastic displacement data in (3.15), 2> = g0V(1)Go • (3-65) The reference Green's operator can be expressed in the form Go = — ^ n r G 0 n = — - ^ I I T ( G Q + G 0 s)n. (3.66) p = — ^ ( G o 1 + G*yVM{G* + G0 s)n . (3.68) puj2 Factoring out G Q leads to X> = - ± n T { l + G f t G S n G ^ G & l + ( G ^ G ^ I I pu2 (3.69) = N ^ G ^ V W G g M ; 1 where N " 1 = J^TiT[1 + Go(Go) - 1) is the receiver decomposition operator and M (1 + ( G O ) - 1 G Q ) I I is the source decomposition operator. Recall that the quantity G Q V ^ G Q is just i M 1 ) , the deghosted data in (P,S) coordinates. Therefore, the operation D = NgT>Mg (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 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 61 In (kg,ks,u>) space, the decomposition operators are simply matrix functions which are multiplied by the data i.e., D(kg,ka,uj) = TSlg(kg,u>)'D(kg,k„u>)Mg(ktt«>). (3.71) This is a pleasing result since it yields a simple transformation which is applied to displace-ment data in the (kg,ks,w) domain. The required parameters are the source and receivers depths (if nonzero) and the P and S wave velocities of the near surface. When the sources d receivers are located on the surface, the source and receiver functions in (kg,ks,oj) an space are Ng{kg,Uj) -i/32 '9 n 2 - f c 2 Kg Kg 2u„ and M , ( f c „ w ) = — v -2r,g kg n 2 - f c 2 2us (3.72) 2r)s (3.73) 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. The body force which describes these new point sources can be written as ^ (r^)^ = S(r3)^F,A~ Written in (x, z,u>) space, this becomes S(xs) 0 T J x z 0 0 8(zs)_ J~ZX 0 A2(UJ) (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 A2(u>) are the two source wavelets. The data for 8T are given by i> = (g - g0)sr = g0vg0^ + dovgovgoF +. (3.75) 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 = g0vg0 + g0vg0vg0 + = T>. (3.76) 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 in (X,Z,UJ) space is A " 1 (a;) 0 0 A?(u) J~mx-Fzz J~mz-Fz: Tzz ~Tz: (3.77) Unfortunately, the source wavelets are practically never known. They must be either esti-mated 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. The result is "DA'1 = "D where A is the 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. When the corrected field data are used in the (P, S) decomposition, the result is D J 4 _ 1 : the (P, S) CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 63 decomposed data scaled by the inverse of the wavelet. Next, I use DA 1 = 2? to rewrite the multiple removal equations in terms of the (P, S) decomposed field data. 6 ' = 6 ( , ) + 6 ' < " + fi(,, + . . . (3.78) where LS ( 1 ) = GdV^GdA f 0 0 . (3.79) = N g X > M s A = N g X > M s , r V ( 2 ) = - D ( 1 ) H ^ D ( 1 ) , (3-80) and in general 6 ' < " , = - £ , ( - I ) i U - ' £ . ( 1 ) . (3.8i) The 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. The 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. The 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 64 3.7 Synthetic data examples In this section, I demonstrate the removal of elastic free surface multiples using synthetic data examples. Data 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. The data processing was done using the Seismic Unix (SU) system from The Colorado School of Mines. The decomposition and multiple removal codes were written as modules to the S U system. 3.7.1 Model 3 The first case I consider is Model 3 shown in 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). The receiver spacing is 10 meters (m) and the shot geometry is split spread (shot in the middle) with a maximum offset of 640 meters. The trace length is 256 samples at 2 milliseconds (ms) to give a recording time of .512 seconds (s). A minimum phase source wavelet was used as input to the finite difference code. The sources are horizontal and vertical point sources on the surface and contain the same wavelet. The 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. The 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. At 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 65 multicomponent source * multicomponent receivers V V V Elastic Free Surface Elastic solid a0 = 2000 m/s = 1200 m/s po = 1.0 g/cc z = 8 0 m Elastic solid 21 = 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 in the direction of wave propagation. As a result, there is significant cross-talk between different components at mid 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. The resulting data panels in 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)'(2\ These in turn are used as input to calculate the next term D'(3). Multiples are removed order by order by summing the terms calculated from the multiple CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 66 0.4 • Vzz Vxz Vzx Vxx Figure 3.8: Model 3 multicomponent displacement data with the reference wavefield (groundroll) removed. The 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. This history is read from left to right e.g. PPPS begins as a P wave, reflects twice as 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 data. The first panel is the PP portion of the input data. The second panel is the estimate of the first order free surface multiples that begin and end as P. When 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. When 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'( 4\ D'( 5\ etc. It is also evident that the primaries remain intact. This is important CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 67 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. The 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. Like the PP example, the first and second order multiples have been removed, the higher order multiples have been altered, and the primaries are left untouched. Suppression rates for this example are around -20 dB ( 90 %). Effect of near surface velocity errors As I stated previously, the properties of the near surface material (a 0 , /?o, Po) are required as input to the multiple removal algorithm. Unfortunately, the near surface can vary consid-CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 68 Figure 3.10: Model 3 PP multiple removal. The PP portion of the input data is Dpp. First and second order elastic free surface multiples are removed from D p p by D p p and D P P . After summing the first three panels, the resulting data D P P 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. As 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. The 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. The corrupted (P, S) data were used along with the wrong velocities to calculate the CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 69 Figure 3.11: Model 3 (P,S) 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. The 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. When the incident angle is small, this becomes — 2iui / ocA(ui)~l. Now if the velocity ea is used, small angle Hpp becomes —2iuj/a(A(uj)~1/e) = — 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 When 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. The 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 m i s s i n g data c o m p o n e n t s It is rare that full four component data are acquired in the field due to the added expense of orthogonal sources. Can 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. The 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. The 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 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA Figure 3.14: Model 3 (P,S) data with missing Vxx. Separation of P and 5 events degraded. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 73 Figure 3.15: Model 3 (P, S) data after two passes of demultiple with missing Vxx. Multiple suppression is still effective although certain events cannot be estimated and removed. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 74 3.7.2 Model 4 The second case I consider is Model 4 shown in Figure 3.16. Model 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. The 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. The data were modelled using a sample interval of 2 ms and a recording time of .512 seconds. The data were arranged into a cube of 64 shots by 64 receivers by 512 time samples. Much of this data cube is zero padding which is useful when working with Fourier transforms. multicomponent source multicomponent receivers V V V Elastic Free Surface Elastic solid a0 = 2000 m/s /?o = 1200 m/s po = 1.0 g/cc z = 0 z = 146 m at shot 446 Elastic solid ai = 4000 m/s ft = 2500 m/s px = 2.4 g/cc Figure 3.16: Land 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 MULTIPLE REMOVAL FOR LAND DATA 75 Figure 3.17: Model 4 displacement data, shot gather 446. between components. 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. Much 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. The data after two passes of free surface demultiple are shown in Figure 3.19. In this figure, all second order multiples have been removed or severely suppressed. 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 78 3.8 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 sat-isfying 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 apply-ing the inverse (P, S) decomposition giving the original primaries back in (ux,uz) 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* 1 ) = NT> M . Substi-tuting this into equations (3.61) - (3.62) gives f>'(2) = NX> ( 1 ) MHN£ ( 1 ) M , D ' ( 3 ) = N 2 ? ( 2 ) M H N X > ( 1 ) M > (3.82) f)'(") = N 2 ? ( n X ) M H N £ ( 1 ) M . Next, I apply the inverse (P, S) decomposition to these equations to get T>'(2) = £ ( 1 ) M H 1 N P P ( 1 ) = X> ( 1 ) K2> ( 1 ) , (3.83) a n c T>' = + X>' ( 2 ) + T>'{3) + ... + X>' ( n ) . (3.84) where K = M H N . The result is a multiple removal algorithm that uses displacement data as input and gives displacement data as output. The previously explicit (P, S) decompo-sition is implicit in the operator K . The four component displacement data after 2 passes of displacement multiple removal are shown in Figure 3.20. CHAPTER 3. FREE SURFACE MULTIPLE REMOVAL FOR LAND DATA 79 Figure 3.20: 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 equit-able 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). The 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 acous-tic. 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. The 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 4 - C 2 - C 1 - C elastic 1 - C acoustic Figure 3 . 2 1 : Model 3 comparison of elastic vs. acoustic 2 pass demultiple in displacements. The left data panel shows the full four component elastic demultipled data in displacements. The middle left panel uses only vertical source two component elastic data. The 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. B y reciprocity, T>zx can be obtained from Vxz; hence only T>xx are zero. The 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. By 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. The 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). The method shown here removes all orders of elastic free surface multiples without requiring information about the earth which generates the reflections. The formulation iteratively generates a series of terms using the multiple contaminated data as input. Each term in this series removes a certain order of free surface multiples. Specifically, the n'th term removes the (n — l ) ' th 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 bot-tom. The 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 in crustal studies to derive P and S wave velocity models from the analysis of refracted arrivals (Mje-dle, 1992; Mjedle et al., 1995). These velocity models are used to infer rock properties through a/f3 ratios. The course sampling of these experiments has made them unsuitable for more detailed exploration purposes due to the lower spatial (and temporal) resolu-tions. Recently, acquisition systems have been designed specifically for exploration with smaller receiver spacings, broader bandwidth and more robust geophones. Many 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. When the shot lines are completed, the recording cable is moved to a new location and shooting begins again. The second reason for performing O B C measurements is to obtain more complete information about the wavefields propagating in the earth. By measuring both P and S wavefields, a more complete picture of the earth can be inferred. Third , O B C data has the potential for improved reservoir characterization. B y having permanent receivers in 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 mul-tiples. The 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 re-ceivers. 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 zs. I choose an explosive CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 85 source because although multicomponent sources for marine data are possible in principle, they are a long way from being a practical reality. Multicomponent receivers are located at a depth zg which can be anywhere at or below the water bottom. When 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. When this wave hits the bottom, it transmits a P wave into the solid ppeikxggii/i(z-u^-z,,)eiu2{zg-zwb) (^ ^ where PP is the P-P water to elastic solid transmission coefficient, and V\ and v2 are the vertical P wavenumbers in the water and elastic solid respectively. The incident plane wave also transmits a converted S wave p£eikxgeiui{zmb-zs)eir)2(zg-zwb) ^ 2) where PS is the P-S water to elastic solid transmission coefficient and r)2 is the vertical S wavenumber in the elastic solid. The 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 accentv denotes 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 eiv\(zwb-zs)eiv2(zg-zwb) ^ _ pp£ivi2zwb _|_ p p2 eivi4zwb _j_ ^ p pgikxg g i i / i ( z (4.3) (1 + ppeivi2z->»>>) The source ghost creates an additional wavefield which is described by an 'image' plane wave source. This gives ppeikxg^eiui(zwb-zs) _ eiui(zwb+z,)^eiu2(zg-zwb) (1 + PPe^2^") ' ^4'4^ The corresponding transmitted S wave will have PS instead of PP. CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 86 As discussed in Section 3.3, the Weyl integral can be used to construct 2-D wavefields from plane wave solutions. Dividing the plane wave solutions by 2ifi and inverse Fourier transforming over k gives the P and S Green's functions G0p(xg,zg\xs,zs < zwb-u>) = r - ^ f m ( ^ ' ) eik(x9-x,)ei^wbeiM*g-^)dk J-oo ^ (1 + PPe^2^) G0S(xg,zg\xs,z3 < zwb;u) = r -PS^ny^) ^ . . - . . l ^ M e M . , - ^ , J -oo + PPe******) These two expressions represent the P and S reference wavefields measured in the solid when the source is in the water. Next, I develop the Green's functions for sources in 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 in the water column. Summing the P waves for these events and using the Weyl integral gives e * ( » . - » . ) ( e ^ l * . - * . l + ppeiM^-2Zwb) oo ^V2 V f / \ \ PPPP ci\v2(z = d 0 V Go + GoV GoV Go + • • • • (4.13) This equation assumes or implies that the sources are multicomponent. I retain this imprac-tical assumption for the purposes of theoretical development. The strategy is to formulate multiple removal for full four component data and then use the two component data in the algorithm. Two 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 in the previous chapter. Transforming equation (4.13) to (P, S) gives D = G 0 V G 0 + G 0 V G 0 V G 0 + . . . . (4.14) CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 89 The inverse series can be formulated as in infinite series of terms in the data V = + + V<2> + . . . (4.15) where D = G 0 V ( 1 ) G 0 (4.16) 0 = G 0 V ( 2 ) G 0 + G o V ^ G o V ^ G o (4.17) 0 = G o V ^ G o + G 0 V ( 2 ) G 0 V W G O + G o V ^ G o V ^ G o + 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. Each of these terms have a specific role to play in the modelling of data in the forward series. The term G w c creates events due to the presence of the water column. The remaining term, G Q , performs the tasks of data construction exclusive of creating water column multiples. Now, since G Q c creates water column multiples in the forward series, it must be responsible for the removal of these multiples in the inverse series. Based on this idea, I select terms in the inverse series which sandwich G Q c to form the following multiple removal subseries V ' = V W + V ( 2 » ' + V ( 3 » ' + . . . (4.19) where V ( 1 ) = G Q ^ G Q 1 y(3Y = - V ^ ' G ^ V M - \rWG™VW - V W G ^ V ^ G ^ V M (4.20) = V ^ G ^ V ^ G ^ V W CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 90 A recursive relation can be derived from these equations which relates the n'th term in the subseries to the previous term V W = - V ^ ' G ^ V ^ (4.21) where the first term, 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 in Section 3.4, G Q V ' G O creates data since v is more like V W than V . The new series is a sum of data terms D = GDV'GD0 = LV 1 ) + D'( 2) + D ' W + . . . . (4.22) As in the land case, a recursive relation can be used to obtain all of the terms in this series. The first term L V x \ are the original data decomposed into up and downgoing P and S waves via f>W = ftx>]vl (4.23) (see Section 4.4 for a more complete description of this process). The n'th term in the multiple removal series is given by D ' W = - D < 1 > ( G O L ) - 1 G 0 R E ( G 0 1 ) - 1 D ( n - 1 > ' = -LV"-1)'Jr>w. or more explicitly 'DfrVfaM D'£rl\KKu) \DPP(k;k„u) D$l(k,k.,u>) DW(k,k.,«>) DW(k,k.,u>) (4.24) D'^XK^) D'&\kg,k.,«) D'£>{kaX,») D'£\kgX,u) - i dk. (4.25) JPP(k,uj) JPS(k,u>) JSP(k,uj) JSs(k,uj) See Appendix D for the derivation of this equation and the elements of J . When the data contain a single source wavelet A(u>), the inverse wavelet B[u>) = A(ix>)~1 must be incorporated into the multiple removal terms. The 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 DATA 91 As in the free surface case, when the earth model is one dimensional, the above integral becomes a multiplication in the (kg,oj) domain and only a single shot gather is required to calculate the multiple removal terms. Each term in the multiple eHmination subseries has an identifiable role to play in the removal of water column multiples. The 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 in 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. The 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. When D ' ( 2 ) is added to D ^ 1 ) , it cancels all first order water column multiples in the data. For water bottom data, first order refers to the number of times a wave interacts with the water column as a whole. The n'th term in 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. The reference medium parameters must be either estimated from the data or obtained a priori. A n obvious problem lies with the fact that the data have only two components. The solution to this is to use the available data components and set the unknown components to zero. The results of this approach will be demonstrated in the synthetic examples. CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 92 4.4 (P,S) decomposition and up-down going wavefield separation The first two steps in 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. As in 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. At 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). Both source and receiver water bottom decomposition are schem-atically shown in 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. As I will show in 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. The first step in formulating the decomposition is to write the reference Green's oper-ators in the form g0 = nTG0n = nT(GD + &™)n. ( 4 . 2 7 ) The data can then expressed as x> = QoVWQo = nT{6d + Gw c)nv( 1>nT(G^ + G^)II (4.28) = nT(GD + Gjnv^CGo 1 + G ^ n CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 93 original data after source decomposition Figure 4.1: Effect of water bottom source and receiver decomposition. The top diagram shows the effect of the partial source decomposition. The middle shows the effect of the receiver decomposition on a scattered P wave. The bottom digram shows the receiver decomposition for a scattered S wave. where V ^ 1 ) is in (P, S) components. Factoring out G Q gives T> = nr(i + cnGIl)-1) [ G J J V ^ 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 = NgT>Ms (4.30) where = I I r(l + ( G ^ ) " 1 ^ ) - 1 and N f l = (1 + G^G*)-1)-1!!. As in the free surface case, the decomposition in (kg,u) is a simple matrix multiplication of the data by the receiver decomposition function kg n2-k2 2rj 2 n2-k2 ^2 K 9 K + _ pi 1-Z2 2p2vi 1+Z2 Plkg 1-Z2 9 i 2p2v\m 1+Z2 (4.31) where Z = exp(iuizwb), the vertical wavenumbers are related to kg, and the subscripts 1 and 2 refer to quantities in the water and solid respectively. The 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 = o r V2 "9 -hi 2u2 2772 + 0 -0 P i 1-Z2 2P2u1 1+Z2 Plkq 1-Z2 2p2v\t)2 1+Z2 N f s + N w c (4.32) The first matrix is exactly the decomposition matrix for an elastic free surface (see Appendix C) . The 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. The 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 -i(32 1 2 " 2u2 v2-2p2ui Plkg 2rj2 k + 9 ' 2p2vir)2 (4.33) Comparable relations have been derived by Donati and Stewart (1996) and also by Amund-sen and Reitan (1995a, 1995b) in a comprehensive analysis of ocean bottom reflection and transmission. CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 95 A problem arises in 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 (4-34) = Ng D(1 + ( G ^ ) - 1 G D ^ 0 ] 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 in (kg,ks,u)) space by the operator (1 + PPell/1'2Zwb) which is a type of Backus filter in f-k space (Backus, 1959). The source deghosting is performed by the operator 1/ sm(vx,zs) where vls = ((uj/ax)2 — k2)1/2. In summary, the steps in the source decomposition are: 1) Fourier transform the data with respect to the shot domain; 2) apply the operator (1 + ppeiv^2z™b) and 3) inverse Fourier transform back to space. (4.35) 4.4.1 Removing the reference wavefield At 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 Go(xg,zg\x3,zs < Zwb]Ul) = J -P Peiu2(z9~z-»>b) -P Seiv2(Za~Zwb^ + PPeiui2z-">) (4.36) {"Its) CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 96 The total wavefield contains the reference wavefield plus the scattered field. When the source decomposition is applied to the total wavefield, the reference wavefield becomes 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. The 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. The codes for the decomposition and multiple removal were written as modules for the SU processing package. 4.5.1 Model 5 The 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. The impedance contrast between the layers is high to generate many orders of multiples. The source is an explosive source located 6 m below the acoustic free surface. The 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. The sample interval is 2 ms and the recording time is .512 s. A n arbitrary minimum phase wavelet is used as a source wavelet. The 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. (4-37) CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 97 Acoustic Free Surface z =0 # z=6 m explosive source ax — 1500 m/s p! = 1000 k g / m 3 z Water — z=8u m multicomponent receivers Elastic solid a 2 = 2500 m/s & = 1200 m/s p2 = 1200 k g / m 3 =180 m Elastic solid a 3 = 1500 m/s 03 = 900 m/s / 9 3 = 1000 k g / m 3 Figure 4.2: Ocean bottom data Model 5. The multiples as a whole are termed water column multiples. Within this class, I dis-tinguish 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 multi-ples. Events which reflect from the top of the water column are termed water top multiples. Referring to Figure 4.3, Ppp and Pps are primaries, Ppppp and Pppps are first order wa-ter bottom multiples, and PppPPpp are PppPPps water top multiples. Events such as PppPP and PppPPPP 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 (2) removes 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. The 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. As discussed before, the source decomposition is only partial since the source is single component. The left panel in Figure 4.4 shows the vertical component of the scattered field. The reference wavefield is removed by applying the source decomposition operator and then muting out the direct wave. The 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. The 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. The first order water bottom multiple is very faint in the S panels; however, the water top multiple is still prevalent. The effect of multiple removal for P data is shown in Figure 4.6. Here, I show the CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 99 Figure 4.4: Decomposition of Water Bottom Model 5. The panel on the left contains the vertical component of the scattered wavefield. The 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. The 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. The effect of multiple removal on the S data is shown in Figure 4.7. The 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. The panel on the left contains the horizontal component of the scattered wavefield. The 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 in estimating the ocean bottom parameters. In general, these parameters must be estimated from the data or obtained from some other source. As in any measurement, errors are unavoidable. Consequently, the question that needs to be addressed is how sensitive is the method to these errors. In Figure 4.8, I show the effect on P data of using water bottom velocities of a2 = 2200 m/s and ft = 1000 m/s where the true velocities are a2 = 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. As in the land case, an altered source wavelet is used to compensate for the errors in material velocity. Although not shown here, the results CHAPTER 4. MULTIPLE REMOVAL FOR OCEAN BOTTOM DATA 101 D ( i ) D '(2) D '(3) Sum Figure 4.6: 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) which removes first order water column multiples, 3) D^3), 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. To do this, I use Model 5 with a dipping lower reflector. As 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. The 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. Both decomposition and multiple removal were performed using velocities of a 2 = 2200 m/s and ft = 1000 m/s. The true velocities are a 2 = 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. The source and receiver spacing is every 12 m. The 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. The receivers then roll along with the shot in a split spread geometry. Since Model 6 is two dimensional, many shots and receivers are required to remove multiples. The data are arranged in a cube 64 shots by 64 receivers by 512 time samples. Much 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. The time sampling interval is 2 ms and the recording time is .512 s. Acoustic Free Surface z=0 ^ z=6 m explosive source Water V V V a i = 1500 m/s px = 1000 k g / m 3 z multicomponent receivers z=80 m Elastic solid a 2 = 2500 m/s 0 2 = 1200 m/s p2 = 1200 k g / m 3 z = 210 m at shot 506 Elastic solid a 3 = 1500 m/s 03 = 900 m/s p3 = 1000 k g / m 3 Figure 4.9: Ocean bottom data Model 6. The dip of the lower reflector is 10 degrees. In Figures 4.10 and 4 .11, I show a portion of the P and S data taken through the various processing steps. The displayed portion of the data is common shot gather 506 from the middle of the survey. The first data panel contains the vertical component of the total displacement field. The 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. To 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. The result is the second data panel in Figure 4.10. The 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 105 The next stage is to apply the receiver (P, S) decomposition operator. Recall that this transforms the data from displacements into upgoing P and S waves at the receiver. The value of the decomposition process is evident in the third panel in Figure 4.10. B y removing the PppPP event, the first order Ppp 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. The 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. The 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 REMOVAL FOR OCEAN BOTTOM DATA 106 T o t a l f ield Scattered field Figure 4.10: 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 VX VX DP + D'P Total field Scattered field 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. The multiples associated with this new 'datum' 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. The 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. The 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. The task of internal multiple removal has been investigated by Araujo 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. The 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) in terms of the elastic scattering equations used in Chapters 3 and 4. The 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. By 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 The 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 in the forward CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 110 series gives a clue as to how they may be removed in the inverse process. 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 in 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 Born 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. As stated, starting with a medium with no free surface and a plane wave source, a Born series for the data can be CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 111 written as D(Z;UJ)= j G0(z\z';uj)V(z')(f>(z';uj)dz' J—oo /oo poo G0(z\z']u)V(z') / G0(z'\z";w)V(z")4iz"-u)dz"dz' + ... (5.1) -oo J — oo (see Chapter 2) where D are the data, G is the reference Green's function, V is the scattering perturbation, is a n incident plane wave, z is the measurement location, and u> is the angular temporal frequency. The reference Green's function for a homogeneous background medium with velocity c 0 is (DeSanto, 1992) G o ( z | * » = 2 ^ e i f c 0 ' Z " 2 ' 1 ( 5- 2) where k0 = UJ/C0. A n inverse series can be written in the form V = Vx + V2 + V3 + . . . (5.3) where Vn is n'th order in the data. Substituting this inverse series into the forward series yields a series of expressions in successive orders of the data D = G0Vx 0 = GoVxGoV^ + G0V24> (5.4) 0 = GoVxGoVxGoVxtp + GoVxGoVxp + GQV2G0Vl(t> + G0V^ The first equation describes the data in terms of V\. Using

/c0). 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. The third term in the series was selected since it involves three orders of interactions with the data. This term is G0V3 = -GoiVxGoV, + V2G0Vx + VxGoVxGoVx)^ = G 0 (Vsi + %2 + %z)4>- (5-6) CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 112 Analysis of these terms by Araujo has shown that V 3 1 and V32 do not lead to a diagram of the type shown in Figure 5.1. The remaining term is G0V33 = -GoViGoViGoVrf = W33. (5.7) As 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 in the integrals. Each of the four terms imply a relationship between the depths of the three scattering interactions with Vi. One of these terms, W333, contains the desired lower-higher-lower relationship and was chosen for further analysis. This term can be written as -| flOO PZ\ POO WZ33(kQ) = —i / dz1ei2h^V1{zl) / dz2e-ak°"V1{z2) / dz3^k^Vx{z3). (5.8) 4 f e 0 J-OO J-OO J Z2 Using step functions to describe the integration limits and evaluating the integrals using contour integration, (5.8) becomes W333{k0) = ^j^i(ko)W1(-k0)W1(k0) + S residues due to poles in V i ( & o ) . (5.9) 4 K g Using the first part of this expression and (5.5), the following expression was selected as a candidate for internal multiple removal D{f3{z; ko) = D(z- k0)D(z; -k0)D(z; k0) (5.10) where D^(z; k0) = -2-^e%h°zW333{ko). When added to the data, D333 successfully removed first order internal multiples (hence the superscript IM); however it also introduced extra events which corrupted the final output. The 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 /o o eikoZ'D(z; z')dz' (5.11) 0 0 and then breaking the limits of the integral into portions leads to a series of expressions which have a geometrical relationship. The portion which contains a lower-higher-lower relationship is /OO PZ\ POO dzieikaZiD(z; zx) / dz2e-ik°z>D(z; z2) / dz3eik°z>D(z; z3). (5.12) CO J —OO v Z2 CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 113 Hence, instead of V\, the data D are taken through the picture in Figure 5.1. When added to the original data, D3 suppresses all first order internal multiples regardless of 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 multi-ples. For example, the term which attenuates second order internal multiples is /CO fiZl / " C O dzie*°ZlD(z; Zl) / dz2e-ikoZ*D(z;z2) / dz3eik°z> D(z; zs) CO V — CO J Z2 dzAe-ik°ZiD(z; Z4) / dz5eik«z>D(z; z5) . (5.13) •CO J Zi 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) + D3(z; k0) + D5(z; k0) + D7(z; kQ) + . . . . (5.14) 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 in (P, S) components (see Chapter 3) D = G o V ^ G o 0 = G o V ^ G o V W C o + G O V ^ G Q 0 = G o V W G o V W G o V W G o + G o V ^ G o V ^ G o + G Q V W G O V ^ G O + G 0 V ( 3 ) G 0 (5.15) CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 114 Using subscript notation, the data can be written as Da = GoikV^Goij (5.16) 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 x2. For an isotropic homogeneous background medium, the Green's function in (P, S) co-ordinates is diagonal and is given by Goi(zg,zg\z.,z.iU) = —] — e * < - - . > c * « l * . - . l < t t (5.17) where 6u = v and rj for i — P and S respectively. The data are then written as /OO -I —^e^te-'^e^i^-'^dkidxidxidzt (5.18) -oo 2i02j K 1 where depends on three spatial variables. Note that 92j is tied to the horizontal wavenumber k2. Fourier transforming (5.18) over the shot and geophone domains and rec-ognizing that the spatial integrals are Fourier transforms, I obtain the compact expression Dij(h, zg\k2, za; U) = -l-e-«i".^1)(fc1, k2, 6U + *ai)JLc-««" . (5.19) Next, I consider the equation which is third order in 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 CO V£\x3, x,, G 0 m ( z 4 , z4\x5, z6; UJ) V$(X6, X6, Z6) e*** e i k ^ d i (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 OO -I /»oo -I to 7 ^ ^ ( - - ^ ) e ^ | Z 2 - ^ l ^ 3 v£\Xi, Z 4 , Z 4 ) 1 Z"00 1 — / — — e i f e 4 ( x 4 - X 5 ) e i e 4 ™ | z 4 - 2 6 | dfc 4 ^V( : c 5,a ;6 ,^)e^ Z 6 e i f e 2 a : 6 df . (5.21) 27T l / _ 0 0 2l0 4 m Because of the absolute value signs in the exponents, this equation must be split into pieces depending on the relative orientation of the depths z2, z 4 , and z&. I choose the region where z 2 > z 4 and z 6 > z 4 which represents a lower-higher-lower relationship between the scattering depths. Equation (5.21) then becomes /CO /> oc i»CO / / e i 6 ^ V ^ \ h , k 3 , z 2 ) 0 0 J — 0 0 J — 00 1 f°° 1 - / —eid^-^dk3H(z2 - zjvWfa, h , z 4 ) to J-oo 1 poo 1 JT" / j r ^ e ^ ^ - ^ ^ t f ^ - z e j v j ^ (5.22) 27T j * _ O Q 2z0 4 m where I have carried out the x integrations and H is the Heaviside or step function. The Fourier transform of the Heaviside function is H(z2 - z4) = — / -dq e - 0 (5.23) 2TT J^ i{q-te) where e is a small positive quantity which approaches zero. That 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) in (5.22) leads to I poo These expression can be evaluated by performing contour integrations in 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 CO 77 V ^ i h , k2,6U + 02j) = — / V^\k„ k3,6U + 03i)-fV^(k3, fc4,-(03l + 64m)) /oo ^r-vS(k4, k2, 6Am + 62j) + £ residues due to poles in V™ . (5.25) •co ^"im Next, I define a modified form of the data Bij{ki,k2, 6U + 82j) = J — e - ^ V ^ i h , k2, 9U + 8 2 i ) e - i 6 ^ (5.26) which is the original data converted to a plane wave source. Following Araujo, I neglect the extra residues in (5.25) and use the new data type to obtain /oo jr —*-BQ(kA,haie^n + e2j) (5.27) oo 2 l04m where I have converted the portion without the residues, l / ( 3 3 3 ) , to _g(333). This equation is the 2-D elastic version of (5.10) and the 2-D equation in 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 Bij(kuk2,z1)ei^'+e^dzl (5.28) •oo which can be split into /OO flZa Bijih, k2, z 1 ) e i ^ i + e ^ d z 1 = / B^h, k2, z 1 ) e i ^ + e ^ d z 1 •oo «/ —oo fiOO + / Bij{k1,k2,z1)ei^'+e^dz1 (5.29) " Za which correspond to z\ < za and z\ > za respectively. Using this splitting technique, I substitute for each B^ in (5.27) and choose the interval z\ > z2 and z3 > z2 (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. When added to B^\ki,k2,6u + 02j), the role of this term is to suppresses all first order internal elastic multiples. The 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. The 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 = z2 = z 3 which distort the primary reflections. In practice, 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 asso-ciated internal multiples. In other words, internal multiples cannot be suppressed between primaries that cannot be resolved. The 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 ^ G o 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 j jr fz\—a JP- / B ^ h ^ z ^ e ^ ^ ^ e ^ ' ' ^ /OO 77 POO oo ^%"hn Jz2+oo B^{k^z2)e-^6^dz2 I B^j(kuz3)ei^+e^dz3. (5.32) OO J Z2-V) 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 in 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. The variable z is termed pseudo depth and is not the actual depth unless the earth comprises only the reference medium. The data in pseudo depth is then input into (5.27) which computes the first order internal multiples. The multiples are calculated for each set of (ki,0u + 8\j) by integrating the data over a range of pseudo depths z. The resulting estimate of the multiples is then transformed back to (ki,uj) and then to (xg,t). In reality, since multiple suppression is a subtractive process, the suppression need not be done in (xg,t); 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. The 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) coordi-nates. This can be done using the method outlined in Chapters 3 and 4. A consideration must be made to the effectiveness of the (P, S) transform when the velocity in the vicinity of the sources and receivers is in error. The resulting artifacts in 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). The 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 in 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, k2, 9u+ &2j)- The 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. 5.3.1 Synthetic data example 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. The data will be four component (P, S) data in (k, z) space for an incident plane wave and a fixed k. Hence, multiple suppression will be implemented using (5.32) for a single k. I constructed the synthetic data analytically in (k,kz) 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 = 2500m/s 02 = 1200 m / s pi = 2000 k g / m 3 200 m a 3 = 1500m/s 03 = 900m/s ps = 1000 k g / m 3 Figure 5.2: Internal elastic multiple Model 7 The 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 re-flections labelled; other events are first and second order multiples. In these labels, the letters following the underscore denote subscripts (e.g., B_PP is the same as Bpp). For this model, converted internal multiples are more prominent on the BSP and Bps traces than the Bpp and Bss 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 BsZpIM. Four traces contribute to this term and are labelled (1) —•> (4) in Figure 5.4. Also shown are two copies of the data, Bg'p. The second of these will have its multiples suppressed (it will be summed with the elements of B^pIM) while the first will remain untouched for comparison. 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 ^ I M are shown in Figures 5.7 and 5.8. 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. This formulation is applicable to land and ocean bottom multiple component data. Common to all the inverse scattering demultiple methods, this method suppresses all orders of elastic internal multiples independent of the earth which causes the reflections. In addition, primaries are not affected. The 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. The 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. The significance of using an elastic background formulation is that the relative ampli-tudes 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 x 10 B_ps 150 300 450 600 750 900 1050 1200 1350 pseudo depth z (m) B_ss 150 300 450 600 750 900 pseudo depth z (m) 1050 1200 1350 Figure 5.3: Model 7 internal multiple elastic data for a fixed wavenumber CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 124 B_sp, multiple suppression terms 0.01 H 0.005 -0.005 -0.01 -0.015 -0.02 -0.025 -0.03 -0.035 -0.04 n 1 1 1 1 1 1 i r a . vi a . o l -* -^ r*-— * i u — Bl Bl B3(l) B3 (2) B3 (3) B3(4) 150 300 450 600 750 900 1050 1200 1350 pseudo depth z (m) Figure 5.4: Elastic internal multiple suppression for B^. The top two traces are the data Bslp. The bottom four are the contributions to Bs3pIM. CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION 125 B_sp, cumulative sum 0.01 0.005 -0.005 -0.01 -0.015 -0.02 -0.025 -0.03 -0.035 -0.04 h n 1 1 r — B1+B3 (1:4) B3(l:4) B3 (2:4) B3 (3:4) B3(4) 150 300 450 600 750 900 1050 1200 1350 pseudo depth z (m) Figure 5.5: Elastic internal multiple suppression for f?£p- 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 SUPPRESSION x 10 Elastic 'demultiple B3 (1:4), B_sp 150 300 450 600 750 900 1050 1200 1350 1500 pseudo depth z (m) 2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 x 10 Acoustic demultiple B3 (1:4), B_sp -i 1 0 150 300 450 600 750 900 1050 1200 1350 1500 pseudo depth z (m) 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 1050 1200 1350 pseudo depth z (m) x 10 B_sp, Elastic internal demultiple Before After 150 300 450 600 750 900 1050 1200 1350 pseudo depth z (m) Figure 5.7: Elastic internal multiple suppression for BpP and BsLp CHAPTER 5. ELASTIC INTERNAL MULTIPLE SUPPRESSION x 10 B_ps, Elastic internal demultiple 1 1 1 o. o. &- a G O c 1 — a. a o-O H G O G O G O ^ ^ l ^ v ^ ^ I I I I Before I x^ -^^ ljl^ ——-^ K-^ —-*w\Jw^~—Ylr—^ ' L I , A f t e r ^ - " - ^ v * . — * 1 \ w y M - * A J W ' ' " — -~~V"4W--*' —• " _ -i .i i l 1 150 300 450 600 750 900 1050 1200 1350 pseudo depth z (m) B_ss, Elastic internal demultiple 150 300 450 600 750 900 1050 1200 1350 pseudo depth z (m) Figure 5.8: 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. As a preliminary, I presented in Chapter 2 a mapping between the terms in the forward scattering series and both primary and multiple reflections for simple seismic problems. The mapping was obtained by calculating the terms in 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 mul-tiple with n changes in propagation direction is described by the n'th and all higher order scattering terms. In addition, I demonstrated that for a specific event, the terms beyond its Born 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 in 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 in. 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 sub-surface 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. The method also requires data that are regularly sampled and contain near offsets. When the source wave-let is unknown, as it usually is, wavelet estimation techniques used in 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. To my knowl-edge, 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. As in the case of land data, primaries are not CHAPTER 6. CONCLUSIONS 131 affected. Synthetic tests for both one and two dimensional earth models demonstrated the effectiveness of this technique. The 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. The 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. The 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. As 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 descrip-tion 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 formula-tion 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: Data 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. The 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. The 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 in with previous works in inverse scattering demultiple. The works of Carvalho et al. (1992) and Araujo et al. (1994) were designed for marine data where the source and receivers are located in the water CHAPTER 6. CONCLUSIONS 133 column. 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. The forward model for gen-erating 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. The inverse series can be sim-ilarly described as propagation in an acoustic reference medium and scattering from the data. The 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 in 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. The 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 McMaken, H . , 1982, Ray methods for waves in elastic solids: Pitman Publishing Inc, London. A k i , K . , and Richards, P. G . , 1980, Quantitative seismology: W . H . Freeman and Co. , U S A . Alverez, G . , and Larner, K . , 1996, Implications of multiple suppression for A V O analysis and CMP-stacked 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: 64th Annual Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 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 BIBLIOGRAPHY 135 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 us-ing pressure and velocity detectors: 59th Annual Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 653-656. Berg, F . J . , Svenning, B . , and Martin, J . , 1996, S U M I C : a new strategic tool for explo-ration 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 Annual 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 Annual Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1481-1484. Carvalho, P. M . , 1992, Free-surface multiple elimination method based on nonlinear inver-sion 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 Annual 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. BIBLIOGRAPHY 136 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 MacKay, S., 1993, Surface multiple attenuation and subsalt imaging: 63rd Annual 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-Hil l , U S A . Fokkema, J . T . , and Van den Berg, P. M . , 1990, Removal of surface-related wave phenom-ena: 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 Kohn, W . , 1952, Construction of a potential from a phase shift: Physical Review, 87, 977-992. Matson, K . H . , 1996, The 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 Nor-way: 57th Meeting of the European Association of Exploration Geophysicists, Extended Abstracts, P038. BIBLIOGRAPHY 137 Mjelde, R., 1992, Shear-waves from 3-component ocean bottom seismographs off Lofoten, Norway indicative of anisotropy of the lower crust: Geophysical Journal International, 110, 283-296. Morse, P. M . , , and Feshbach, H . , 1953, Methods of theoretical physics: McGraw Hill , New York. Moses, H . E . , 1956, Calculation of the scattering potential from reflection coefficients: Physical Review, 102, 559-567. O'Brien, M . J . , and Gray, S. H . , 1996, Can we image beneath salt?: The 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 Annual Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1531-1533. Prosser, R. T . , 1975, The formal solutions of inverse scattering problems II: Journal of Mathematical Physics, 17, 1775-1779. Prosser, R. T . , 1980, The formal solutions of inverse scattering problems III: Journal of Mathematical Physics, 21, 2648-2653. Razavy, M . , 1975, Determination of the wave velocity in 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, The role of the Born approximation in 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. BIBLIOGRAPHY 138 Stolt, R. H . , and Weglein, A . B. , 1985, Migration and inversion of seismic data: Geophys-ics, 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 Annual Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1017-1020. 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, Decompo-sition 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, Amplitude 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, The sensitivity of Born inversion to the choice of reference velocity: A simple example: Geophysics, 48, 36-38. Weglein, A . B . , and Stolt, R. H . , 1995, I. The wave physics of downward continuation, wavelet estimation and volume and surface scattering. II. Approaches to linear and non-linear migration inversion, in Symes, W . W . , E d . , Mathematical Frontiers in Reflection Seismology: S E G / S I A M publication. BIBLIOGRAPHY 139 Weglein, A . B . , Gasparotto, F . A . , Carvalho, P. M . , and Stolt, R. H . , 1997, A n inverse scattering series method for attenuating multiples in seismic reflection data: Geophysics, (to appear). Weglein, A . B . , 1985, The inverse scattering concept and its seismic applications, in Fitch, 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: 65th Annual Meeting of the Society of Exploration Geophysicists, Expanded Abstracts, 1492-1495. Yilmaz, 0. , 1987, Seismic data processing: Society of Exploration Geophysicists, Tulsa. Appendix A Calculated series for 1-D Born scattering This appendix contains the calculated series expressions for the first and second primaries and the first order multiple in 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. The first primary reflection is given by the series lfcgai , 1 f k l a ^ 2 5 (fcpffli This can be expressed as Pprl(xg,zg < zx\k;uj) = Rieikx9eiM^-zg) (A2j where K-i. — — — = ; ( A . d ) KavK vo + vi is the reflection coefficient corresponding to the first interface. 140 APPENDIX A. CALCULATED SERIES FOR 1-D BORN SCATTERING 141 The series for the second primary reflection is 1. 1 P^ix^Zg < ZilJfejw) = e ^ e - o ( 2 2 2 - , f l ) 2 + ( - 3 2 " 6 4 a 2 -a? 4 a, 4 a 9 H -a? 8 16 2 128 2 256 2 21 1024' + ( 1 ai ^ 8 128 2 : « 2 + i 2 128"" ' l b 2 4 ° 2 ) a i + ( ~ 1 6 " 6 4 a 2 256" a 13 . , 1 27 A . . 4 , 15 , 5 a 2 ) « i + ( - 7 ^ 7 o : ) a i 2048 -o?^ i 2 )a{ 2048 a 2 ) a 3 + (-+ , L 1 - 2 ( « 2 0"> y 8 16 2 128 256 64 2048 y 1 v 2048' 2 ) a i + ( 8 " 3 2 a 2 ~ 1 2 8 a 2 " 2 ^ 6 a 2 ) a i , 3 3 a 2 + ( 1 al ~ V32 512 2 iz/ 0(z 2 - ^ l ) + 1 1 o 5 o , 1 1 A 5 „•> , 1 3 „ . a 4 < " l 6 a 2 - 32" 2 ~ 2 5 6 4 ) a ? + ( 16 " 3 2 ° 2 - 256 ^ ^ + (16 " 256 ^ V 256 ; 1 v 192 ; 1 ^o(z2 - z-x)2 + ivl(z2 - z i ) 3 + I . . . , , , l ) a . (+192«0^ + (- 192' ^o(^2 - ^ i )4 + . . . (A-4) where a = k\ajv\. Factoring this equation gives 1 Ppr2{xg,zg < zx\k;uj) = ei^9eiM^2-zg)2 8 a2 + —a2 + , 1 1 .2 1 -3 7 - , / 1 1 + ( - o + T ^ ° 2 + 7 ^ ° 2 + 7 W ^ ) a i + ( " i 6 - 64 8 128 + ( a 2 V 32 64 128"' ' 1024 13 -. 2\ "3 , i 1 128 « 2 a, + 256 i 2 256 2048 21 1024* a2)a1 27 2048 64 2048 (A.5) 2! + ( _ 2 0 4 8 ) a i + • • • j ( 1 + i z yo2(z 2 - z ! ) ( 7 l - 1) zV 0 38(z 2 - z Q 3 ( 7 i - I ) 3 , ( z V 0 2 ( z 2 - z 1 ) ( 7 i - l ) r 3! + ' n! APPENDIX A. CALCULATED SERIES FOR 1-D BORN SCATTERING 142 which further leads to Ppr2(xg,zg < Zi\k;u>) = eikx°eiv°(2z*-z^ R2T1T[(l + iv02(z2 - z0(7i - 1) v2A{z2 - ^ ) 2 ( 7 i - I) 2 i^0S(z2 - z x ) 3 ( 7 l - l ) 3 3! . •. + 2! (iu02(z2 - z x ) ( 7 l - 1))T (A.6) = ei k x 9 eivo(2z-2-zg) R2 J " g > « / l ) 2 ( 2 2 - 2 i ) ( 7 i - l ) — j[£2 J 1 i / e ^ z 9 e » l / o ( 2 2 i _ z 9 ) e ' l / o 2 ( z 2 - z i ) 7 i _ ^ 2 J , i 2 ^ g i f c 2 f l g J i / o ( 2 z 1 - z 9 ) e i i / i 2(22-21) where 71 = (1 — d i ) 1 / 2 = vi/va- The reflection coefficients from the first and second interfaces are ^> * o / 7 i * \ (A.7) P i = 2 - d x - 2^(1 - a i ) _ zy0 - f i a i ^ 0 + ^ 1 and R2 = « 2 - a x - 2^(1 - a i ) ( l - a 2 ) _ 1^ - ^2 (A.8) d 2 — Ol Vi + ^ 2 respectively. The transmission coefficients through the first reflector in the forward and backward directions are 2\ = 1 — Ri and T[ = 1 + Ri . Thus, the sum of terms in (A.4) leads to a second primary which has the correct amplitude and propagation velocity. The series for the first order multiple is P m l t l ( x g , z g < zi\k;uj) = e ^ e - o ( 2 ( 2 z 2 - 2 1 ) - 2 , ) 1 n2 1 a* ~ e l " 2 " ^ 2 2 64 ~~z 512 a\ )d i + 32 d 2 + 128 1024 -a, a + + 1 1 1 - L A - L ^ 2 \ ' - 3 i 1 a-2 H a 9 a i + 64 32 64 2 ' 1 a2 a, + 128 128 2 ' 1 13 512 ' 1 3 2 a 2 + 3 2 a 2 ' a i 1 5 A i . 4 + 1 & 2 (Xi ^ 1 32 64 ' 1 + 1 32 ivo (z2 - zi) + 1 1 . 1 . 1 ^ 2 1 . 3 — a a 2 H a i 16 128 128 2 ' 1 1 - 2 , 1 " 3 \ - 3 1 16 a2 \ a\ vl{z2 - Zif (A.9) APPENDIX A. CALCULATED SERIES FOR 1-D BORN SCATTERING 143 This can be written as Pmla(xg,zg < Zl\k;co) = - Rx(l - P 2 ) ^ e ^ 9 e - o ( 2 ( 2 , 2 - , 1 ) - z 9 ) ( ^ 1 + 4 ( Z 2 _ 2 l ) ( 7 l _ X)ivo _ [4(z2 - z1)(ll - l)]2vl [ ^ o 4 ( g 2 - C l ) ( 7 i - l ) ] w \ 2! n\ J — _ J2X(1 — ^ 2 ^ 2 e i f c x 9 e i j / o ( 2 ( 2 z 2 - z i ) - z s ) e i i / 0 4 ( 2 2 - 2 l ) ( - Y i - l ) = - P l ( l - R l ) R * e * * 9 e i M ^ - * 9 ) e 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. Each coefficient contains an incident phase followed by a scattered phase. Hence, PS represents an incident P wave and a scattered S wave. A grave accentv denotes 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 The derivation of the reflection coefficients for an free surface has been done previously by various authors. To name a few: Ewing, Jardetsky, and Press (1957), A k i and Richards (1980), and Wapenaar et al. (1990). The derivation shown here follows closely the treat-ment of A k i and Richards (1980). The 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. The 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. The first step is to consider a unit amplitude P wave incident on a free surface. The potential for this incident wave is exp(i(kx + v(za — z))). Reflection from the free sur-face creates both a reflected P wave and a converted S wave. This gives the potentials PPexp(i(kx -f- u(z3 + z))) and PS exp(i(kx + vzs + rjz))) where PP and PS are the P-P and P-S reflection coefficients. The displacements for this experiment can be written as dx -dz ( ei(kx+v(zs-z)) + = dz dx \ 0 p p ei(kx+v(zs+z)) p Cjei{kx+vza+r)z)) (B.l) The tractions are related to the displacements by T J- Xz fidz (j>dx ux T ± zz \dx idz (B.2) where \L US the shear modulus, A is the Lame constant, 7 = A + 2fi is the uniaxial strain modulus, and Txz and Tzz 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 ei(kx+vz3) -2pku /J,(T]2 - k2) - / i ( 7 / 2 - k2) -2\Lki) ppei(kx+uzs) p Cjei(kx+vza) + 2fiku n(r)2 - k2) -fi(V2 - k2) 2fikV 0 = 0. (B.3) This equation can be solved for the reflection coefficients PP and PS provided that the phases in 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. The resulting reflection coefficients for both types of incident waves are PP = SP = (B.4) — 4fc ;y | T7" — fc' I PS 4k2vr) — {r}2 - k2)2 4k2wq + (T/ 2 - k2)2 4:krj(i]2 — k2) 4k2vri + (n 2 - k2)2 Aku(r]2 --k2) 4k2ur) + (n 2 - k2)2 Ak2vq — (rj2 - k2)2 4k2uTj + (rj2 - k2)2 SS 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. At a solid-fluid contact, the horizontal component of displacement may be discontinuous (but not necessarily zero) because the fluid may slip past the solid (Aki and Richards, 1980). Barring cavitation, the vertical component is, however, continuous across the interface. In summary, the three boundary conditions are: 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. The first is an incident P wave in the fluid. The 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 and/or 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 2uJ2p1u1(r]l -k2) where PP + # P 2 * i ( - (vl ~ k2)2 + 4:V2u2k2)) U) PP = — v r -^wlfi + f32P2Vi ((vl - k2)2 + 4V2u2k2)) PP ^-2u>2p2v2{j}2 - k2) PP = PS = PS = Akuiu2uj2 pi ^wb -±kv2vifi\p2{r)2 - A;2) wb Ar)2vxk^2p2{j]2 - k2) 'wb _ -4kj]2v2u>lp2 ^wb -a>VizV/322 + ( ' (Vl - k2)2 + ^2v2k2) S S = n (B.5) SP = - " " " r ^ " (B . 6 ) 0 w b = w W # + ftp2vx ((nl - A; 2) 2 + AV2u2k2)) . (B.7) Appendix C (P,S) decomposition and up-down wave separation C l 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 in (3.15) T> = g0vwg0. (ci) The reference Green's operators can be expressed in the form a, = ^ n T G 0 n = ^ n ^ G * + G*)n. (c.2) This follows from the fact that away from the source, V2Gp0 = —u>2/a2Gpo] hence V~2Gp0 = — a2/u>2Gp0. Similarly V _ 2 G s o = —ft2/UJ2GSO- In operator notation then v- 2 n T r- 1 G 0 = —^n r G 0 . (c.3) pu>z Using this form of the Green's operator, the data can then be written as T> = ~nT(GD0 + G0s)(nv(1)^nT)(G^ + c0s)n ^ p " (C.4) = ~nT(GD0 + Gf)v(Gj + G0s)n. 148 APPENDIX C. (P,S) DECOMPOSITION AND UP-DOWN WAVE SEPARATION 149 Applying T> to a delta function S(ra) and substituting equation (3.55) into the receiver side of (C.4) leads to T>(xg, zs\x3, zs;u>) 1 f°° f°° pw2 J_ co o — o o e - ' v g z g 0 2ivg 0 + 2ir/g p p e > " s z s S Pe™9*3 PSeir><>zs SSeirK>zs 2ivQ e,7>9z 2XT)g dkn /CO fiOO / V^\x', z'|r")Go(r"|r'"; U)US(T'" - r.) dx'dz'dr"dr'" . (C.5) - o o J — 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') / ikg irlg e - ' " 9 z 9 0 x')[ 2iug \ -Wg ikg 0 e~iVgzg 2ir)g ikg irfg Z l/g % kg PSeir>°z° SSeir,'z» IVgi 2iv„ £Vg*_ 2irjg dkg. (C.6) 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 zg = 0 gives /CO eikg(Xg-X>) •OO ikg irjg -ivg ikg ikg -ir/g 11/g t'kg PP SP PS SS e ' v 9 z 2ivg 0 0 2ir)g dk0 -i: Jkg(xg-X') T\T —1 9 2ivg 0 0 2ir)g J dkg. (C.7) Now, except for the N " 1 , the quantity on the right hand side of this equation is exactly the directly propagating Green's function G^Xg^Zg = 0|a:',z';u;) (see equation (3.55)). By eliminating N " 1 , TIT(GQ -f G Q ) is effectively replaced with G „ and the receiver ( P , S) decomposition and deghosting is accomplished in one step. To 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. The result is just upgoing P and S waves measured at the receivers. In the (UJ, kg) domain, the receiver decomposition matrix is -i(32 UJ" "9 Kg 2v„ (C.8) 2Vg > V 9 The source side decomposition can be obtained in a similar fashion. The objective is to convert G0H6 to GdS in (C.5). In this case, the operator II can be treated as if it operates on the quantity to the left. This because when G 0 operates on IIf5, the sifting property of the delta function acts so that the derivatives are applied to the source coordinate of Go-Now / Te i v ' (*"- z >) G0(X",Z"\XS,Z3;UJ] 0i k s ( x " - x s ) j 2iu„ 0 0 + 2iva 0 0 2iv3 0 0 2iVs PP SP PS ss 2 k 0 ^ 2ir]s 0 eir)sZ° \dk3U. (C.9) Setting the source depth to zero gives 0 i k s ( x " ' — x 3 ) 2iva 0 0 iks —iv3 + ik3 2ir), 2iv3 0 0 2ir]s PP SP PS SS 2ivs ^ —iks ivs ^ 2ir)s _ —irjs iks I dk, /oo &iks(x"-Xt • 0 0 2ius 0 0 2ir]s M 7 1 dk. /oo e-IK>X°Gd(x",z"\ks,zz = 0;uj)M;1dks. ( C I O ) • 0 0 This equation says that by Fourier transforming over the shot domain, multiplying by the matrix M3, then applying the inverse Fourier transform, this gives the downgoing P and S waves at the source. In the (uj,k3) domain, the source decomposition matrix is M3 = -i/32 UJ' 2rj„ ~~~2^7~ K s (CU) APPENDIX C (P,S) DECOMPOSITION AND UP-DOWN WAVE SEPARATION 151 C.2 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> = -^nT(Gd0 + G W C ) V ( 1 ) ( G * + G o c ) n , (C.12) where ( G ^ + G ^ 0 ) is the Green's operator for the water bottom reference medium. A particular form of the Green's function related to G 0 will be useful both for decomposition and for formulating multiple removal. For sources and receivers in the solid, this form is 1 r Go(xg,zg\x„zt;u) = — J eiv2 l*s-*gl 0 ^tk^X g Xg o 0 2JJ72 + Rwc e ?v lsp •i/2(z« + z9) p w c e«("2*.+T)2*g) ps 2i^2 1 W C e<12 (*. + * « ) R ss 2im dk (C.13) where -iv22zwb pwc rtpp K P S -0 e-iz~wb(v2+n2) PP + Z2(PPPP - PPPP) 0 SP + Z2(SPPP - SPPP) pwc n s p pwc K s s ( ^ + 1 2 ) r , , (C.14) 0 PS + Z2(PSPP - PPPS) 0 SS + Z2{SSPP - SPPS) in which 0 = 1 + PPZ2, Z = elVl2Zwb, v and rj are the P and S vertical wavenumbers related to the receiver wavenumber, and the subscripts 1 and 2 pertain to the water and solid. The utility of (C.13) is that it isolates the portion of G W C that depends on the source and receiver coordinates. Substituting this equation into the receiver side of (C.12) and APPENDIX C. (P,S) DECOMPOSITION AND UP-DOWN WAVE SEPARATION 152 operating on a delta function leads to />oo «co I /too ^ ' • [ ' " ' " a ) = ^ L L v - J ik —irj2 \IV2 If ("CO fOO PC0 e R iv2(z' + zg) PP 2iu2 w c eH"2z' + V2z9) SP 2iu2 ^wc e R ik ir]2 —iv2 ik '(V2z' + "2zg) PS 2ir,2 0 \dk 0 eiii2(z'-zg) 2iV2 /C />co / V ^ 1 V> z > " , z")G 0 (r"|r" ' ; U)US(T'" - p.) dx'dz'dr"dr'" . (C.15) •CO t/ —CO This portion in curly braces can be expressed as + ifc —irj2 iv2 ik ik irj2 —iv2 ik R^cpeiu22zg R™ei(V2+V2)zg JV2(z'-Zg) eik(xg-x') J 2iv2 0 eiv2(z'-zg) 2iv>2 0 0 e'V2(z'~zg) 2ir)2 0 eiV2(z'-zg) 21772 dk dk (C.16) Recognizing that the integral without N " 1 is just the directly propagating (P, S) Green's function, it follows that N g is the matrix which transforms the data into upgoing P and S waves at the receiver i.e. G Q V ^ G O . As in free surface decomposition, the steps are: 1) Fourier transform the data over time and the geophone domain; 2) apply the matrix N g ; and 3) inverse Fourier transform back to space and time. When the receivers are on the ocean bottom, zg = zwb and the receiver decomposition matrix in the (UJ, kg) domain becomes -iff UJ' •02- 1-Z2 r>22~k2g 2U2 2P2v\ 1+Z2 P\kg_ \-Z2 K + (C.17) 2V2 9 ^~ 2p2fir72 1+Z2 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). Appendix D Elastic interface multiple removal series In this appendix, I derive in detail the equations for both land and ocean bottom multiple removal. The 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. The data are modelled by the forward scattering series as D = G 0 V G o + G o V G o V G o + . . . ( D . l ) where D are the data in (P, S) coordinates, Go 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, 1995) V = V( 1 ) + V( 2 ) + V( 2 ) + . . . (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 0 V « G o (D.3) 0 = G o V ^ G o + G O V W G O V W G O (D.4) 153 APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 154 0 = G o V ( 3 ) G o + G o V ^ G o V W G o + G O V W G O V ^ G O + 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 ( 1 ) , and 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 in operator notation (see equation (3.11)) G = Go + GoVG (D.6) or written in terms of functions /CO />oo / ^oij(rg|r'; UJ)VJI(T'\T"; UJ) £ < m ( r " | r s ; UJ) dr'dr" •oo J—oo (D.7) where £/irn(rg\rs]uj) is the Green's function in the actual medium, ^oim(rg|rs;UJ) is the reference medium Green's function, and VJI(T'\T"; UJ) is the function which pertains to the perturbation operator V = C0 — C. It is not obvious why Vj>(r'|r"; a>) should depend on four spatial variables and not two. The 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. To 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 SERIES 155 where p is the density, K is the bulk modulus, and p is the pressure. Written in terms of an integral wave operator, this equation is tp=r (v • PV)V+wf2)s(r _ =rX(R|R,; W M R ' ) * ' ( D - 9 ) where L(r|r' ; UJ) is the wave function for the scalar wave operator L. The wave equation in a constant reference medium can be written as L0p = 0 and can also be expressed as an integral operator equation L0p = [°° (v • - V + ^-u>2) 8(r - r')P(r')dT' = / ° ° i 0 ( r|r ' ; u,)P(r')dT'. (D.10) . / - c o V PO A O / J-oo Subtracting equations (D.10) and (D.9) gives /CO POO ( V • a p ( r ) V + aK(r)u,2) 6(r - r')p(r')^r' = / V(r\v'; u)p(T')dt'. • o o J — o o ( D . l l ) where V(T\T';UJ) = ( V • ap(r)V + aK(r)u>2) 8(r - r'), a p(r) = l / / 9 0 - l//>(r), and aK(r) = 1/i^o — l / ^ ( 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 in the earth does not vary, then ap = 0 and V(T\T';UJ) looks like a^(r)a; 25(r — r'). 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>). The 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). The 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. Conse-quently, 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^ 2 ) + V ^ 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 it is necessary to cast 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. The form of V^ 1 ) and its implications is discussed in more detail by Stolt and Weglein (1985). Since I am not explicitly performing inversion, I can write a general, albeit impractical, form for V ^ . The 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. APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES D . l Free surface multiple removal 157 Determining Vt 1 ) The 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 DpP DPS DSP Dss PO " T ^PP SP UPS SS ( 1 ) f/(l) VK ' V pp v, PS vSP> v: a ) SS PO ^PP Qis SP G" fs PS Gdso + Gk SS (D.12) (D..13) Expanding out the first term, bpp = (GP0 + G%P) Vpp (Gpo + G%p) + G%s Vgp (GP0 + G%P) + (GPO + Gpp) Vps G*SP + GpS Vgg GfgP , it is clear that Dpp contains contributions from all of the elements in 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 ^ 1 ) . As a preview of things to come, the multiple removal procedure does not require any direct information about the earth properties themselves. A l l information about the earth is contained in V^ 1 ) which is derived from the data. In fact, V^ 1 ) can be thought of as representing reflectivity as opposed to mechanical earth parameters. Returning to the details of the calculation for Dpp, I write /CO />CO /"CO />co / / dxdz dx'dz' • •co J — co J — co J —co (Gpo + G%)(kg, zg\x, z; UJ) V$(x, z\x\ z'- UJ) (GdP0 + G%P){x, z\k„ zs; UJ) + G%s(kg, zg\x, z; UJ) V$(X, z\x', z'; UJ) (GPo + G%P)(x, z\ks, za; UJ) + (GP0 + Gpp)(kg, zg\x, z]u)V$(xt z\x', z'-,u) Gsp(x, z\k., z,\u) + G%(kg,Zg\x,z]u,)V&\x,z\x',z'-,uj)G%(x,z\ks,zs;uj)\ (D.14) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 158 which becomes after substituting equations (3.55) r flOO pOO poo poo Dpp(kg, zg\k3, ZS;UJ) = / / / / dxdzdx'dz'-v — co z> + e_il/sZs) 2*% ( P P 2iV, ,iuazs _ i _ — i i / s 2 s ^ 2ZZA, 2ii/, 2ZZA, + SP- e - i ( * - B - ^ V i i ) ( x , z | ^ z > ) e I ( F C ' , B , + ' ' ' Z ' ) P 5 - . 227/g 2zz/ s (D.15) Recognizing that the integrals in (D.15) represent Fourier transforms, (D.15) becomes Dpp(kg,zg\ks,zs;uj) ( P P e l l / 3 Z 3 I -11/,zs 2iua Vpp(k9, ~Vg\kt,U,\u) (PPeivsZ' + e~ivsZs) ) 2iv, 2iu, + SP:——V^P\kg,-r,g\ks,iy3-,u)) 2irjg + -V£s(ko>-vg\k»>y''>u)PS-2J^ I V , Z , + S P ^ — vss(h, -vg\K,v.;") ps^— (D.16) (D.17) A similar calculation shows that each term in D depends on the quantities vpp(ka> ~vg\k^ u»\u>) Vsp\k9i-v9\ks, vps(k9 i~Vg\ks,fs;^) Vss\kg ,-r)g\ks,rjs;u>) with different coefficients that depend on kg, zg, k3, zs and UJ. This can summarized by-writing T)(kg, zg\k3, ZS;UJ) as the product of three matrices D(kg,zg\k„z.;u) = qg(kg,zg;uj)Y^\kg,k3;uj)Q3(k3,z3;uj), (D.18) where D(kg,Zg\k3,Z3]Uj) DPP(kg, zg\k3, z3; UJ) DPS(kg, zg\k3, z3; UJ) DSp(kg, zg\k3, z3; UJ) DSs(kg, zg\k3, z3; UJ) (D.19) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 159 Qg(kg,zg;uj) (PPe"/9z9+e-"'9z<>) 2iv„ SP(-PS e">9*i> 2iu„ 2ir)g (SSeiv9zs+e~i7>9z9) 2irie (D.20) V ( 1 ) (* , , * , ;< - ) = Vpp (h,-Vg\k*, v,\UJ), Vp1] (kg, -ug\k„ 17,; UJ) V S P {ka,-Vg\ksiVs]u) V^J (kg, -r]g\kA,T]s;uj) (D.21) and Q,(k„za]w) (pp e»'.*. + e-«'«*.) 2iv3 SPe-PS'-2ivs 2irjs 2iris (D.22) Note the symmetry between the subscripts and the Fourier transform variables in 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 kg, ks and 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^ 1 ) to be determined from measurements away from the scattering region. The task remains to find the elements of from (D.18). Since Qg and Qg are both 2 by 2 matrices, their inverses are particularly simple. Pre- and postmultiplying (D.18) with these inverse matrices yields V^(kg,ks;u>) =Q-1(kg,zg;uj) D(kg,zg\ks,zs;u;)Q;1(k3,zs-u>) (D.23) where and det Q s (SSe'v9z9+e-'r>9z9 2ir)g —PS?iVgZg Q - i d e t Q s 2iu„ ( S 5 e " " x ' + e ~ ~ " » x < ) 2irj~s -PS*--SP e'"9z9 2ir)g (PPeiv9z9+e-iv9z9) 2iv„ -SP e'Vsza 2ir}B (PPeiv>>z»+e-i',>z*) (D.24) (D.25) 2iu, 2ius The invertibiHty of Q g and Q3 depends on the properties of determinants of these matrices. If det Q g and det Q s do not have any zeros over a suitable range of k, then APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 160 Q g a n d Q 3 are invertible over that range. By plotting the determinants of these matrices as a function of k, I found that they are nonzero for \k\ < uj/a which is the range for nonevanescent P waves. As 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 in 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, my 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: GD and G0S. The latter is responsible for the creation of free surface 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 Gf0s and write them as a multiple removal subseries. Hence, v'(2) = - V W G ^ V W V ' ( 3 ) = -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$(kg,-ug\kt,u,\w Vs{p\kg,-r]g\ks,us;uj ) VP{s\kg,-vg\ka,r}s;uj) ) Vs(s\kg,-rig\ks,r)a;uj) » 7 , / . . . / \VPP (k9> G%P(X,Z\X',Z';UJ) G] G\ PS SS (a;, z\x', z'\ui (x, z\x', z'\ UJ (D.27) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 161 where I have Fourier transformed over the shot and geophone domains on each side of the equation. This is done so that V'( 2 ) can be added to V^ 1 ) in (D.21). The first component for V'( 2 ) is 2 /*°° f* oo poo fioo poo Vpp(kg,-Ug\ks,utiu) = — / / / / / dkdxdzdx'dz' v — C O v — O O v — OO v — OO v — OO i(k(x-x')+v(z+z')) VpP{1\kg, -ug\x, z-UJ) PP VpPW(x', z'\ka, va; UJ) + VpS^\kg,-ug\x,z;u>)PS IV ei(k{x-x')+vz' +t)z) VpP^\x',z'\ks,us;uj) IV „i(k(x—x')+uz+T)z') + Vpp^(kg, -ug\x, z-UJ) SP : VSP^KX', z'\kt, va; UJ) ITj x /)+rj(^-|-z /)) + VPS{1\kg,-vg\x, z; UJ) SS- : VSP(1\x', z'\ka, VS;UJ). (D.28) IT} Recognizing that the spatial integrals are Fourier transforms, this equation reduces to Vpp{kg,-vg\kayva\uj) = -:— / dkVPP{1\kg> -vg\k,u;uj) -r—VPp^\k, -v\ks,v3;uj) PS + VPS{1)(kg,-vg\k, v, w) — VPPW(k,-v\ka,va; UJ) XV SP + Vpp^Xkg,-vg\k, v; UJ) — VSP{1)(k, -V\ks, va; UJ) 17] SS + Vps{1)(kg,-vg\k,mUJ) — VSPW{k, -r,\k„ va]uj). (D.29) Similar results can be found for VSp, VpS\ and VSKS} leading to the matrix equation /oo V^\kg, k;u>)G{*(k;uj)V^\k, ka;uj)dk (D.30) •oo where V'(2\kg, ka; UJ) is as shown in (D.27), Vppih,-Vg\ks,Vs\u) V$(kg,-vg\ka,r)s;uj) I V£p\kg, -vg\ka,va;uj) Vss'{kg, -r)g\ka,va;uj) (D.31) and PP/iv SP/irj PS/iv SS/ir) (D.32) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 162 It is possible to determine the n'th order multiple removal term from using the recursive relation /oo V ' t " " 1 ^ , k;u,)G0s(fc;u;)V«(fc, ks;oj) dk. (D.33) •oo 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. The reconstruction is done by sandwiching V ' between two directly propagating Green's functions. This gives the data series D ' = GdVGd = G f t V * 1 ) + V'( 2 ) + V'( 3 ) + . . . )Gd (D.34) = £ > « + D'( 2) + D'( 3) + . . . . The equation for the first term is n(i) 'PP ^PS u SP SS rid {rWrid rid \/-V-)rid ^PO VPP ^PO ^PO VPS ^SO rid T/(l)Ad rid T / M A l ^50 VSP ^PO ^SO VSS ^S0 r ( l ) A d (D.35) Fourier transforming over the shot and geophone domains, DPP becomes /O O fiOO flOO fi o o / / / dxdx'dzdz' • O O J—OO J—OO J—OO ei{-kgX + Ug{z-Zg)) — V p p ( a ; , z\x', z ;u>) 2iu, 9 ,i{kax'+us{Z'-z,)) 2%v, (D.36) 2iv„ -V^p\kg,-ug\ka,us;oj)-2iv„ This is possible because the scattering points are always lower than the sources and receiv-ers. A similar calculation for the other elements in D ' yields the matrix expression D'W(kg,kt;u) = Kg(kg,zg;u,)V(1\k]u>)Rt(k.,z.]u>) (D.37) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 163 where V^(k;u) is given by (D.21), Rg(kg, ka;uj) = e—"gzg 2iva 0 e-ir>a'9 2ir]g (D.38) and similarly for R s (Ai s , z3\ UJ). The matrices Rg and Ra can be said to project onto data space. From here, I substitute (D.23) into (D.37) to obtain D ' ( 1 ) ( A ; g , ka; UJ) = Rg(kg, zg;UJ) Qg1(kg, zg;UJ) T)(kg, zg\ks, za;uj) Qj1(k3) za\ UJ) Ra(ka,z3; UJ) (D.39) which effectively removes the source and receiver ghosts from the original data. The source and receiver deghosting matrices have the form [ R . Q J 1 ] (kg,zg;uj) = det O, e~xv9 -(SSeir)sz» + e~irisze) —SP -PS e-ir>9z9 and [Q:lR3] (k3,za;uj) detQa ^(SSeir>°z' + e~iruz°) e — ' v 3 z s 1. iUgVg (PPeil/gZg + e~iVaZg\ iVgT)g \ /_ (D.40) -SP^--PS (PPeiu*Zs + e-iu°z*) (DAI) For the second term in the multiple removal series, f>(2)' = G ^ V ( 2 ) ' G d . In Fourier transform space this becomes T>'{2)(kg, zg\ka, z3; UJ) = Rg(kg, zg; uj)V'W(kg, k3;uj) Ra(ka, ZS;UJ) /oo Rg(kg, zg; UJ) V^(kg, k;UJ) G0s(fc;UJ) V^(k, ka;uj) Ra(ka, z3- UJ) dk • CO /O O z0\Qg\kg, zg; UJ) D(kg, zg\k, za\UJ) Q ; 1 ^ , za; UJ) Ra(k, z3] UJ) O O • R;\k, za; UJ) G0s(k; u))R-\k, zg; UJ) • Rg(k, zg-UJ) Qg^k, zg- UJ)D(fc,ZG|k3,za;UJ) Q 7 1 (A;,, za;u>)R,(As,, za;UJ)dk (D.42) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 164 where I have used (D.23) to substitute for and (D.30) for V'( 2 ) . Expressing the matrix R 7 1 G ^ R ; 1 as H and using (D.39), I obtain /CO B^(kg, zg\k, za;UJ) H(k; UJ) D™{k, zg\ka, za; UJ) dk (D.43) -co where (D.44) H(A3;w) = -7T ivPP iuSP irjPS ir]SS The term D ' ( 2 ) removes all first order elastic free surface multiples from the deghosted data by using the deghosted data to model the events which reflect from the free surface. For the n'th term in the series, I obtain D^n\kg, zg\k„, za;uj) = Rg(kg, z g ; c ^ V ' ^ 1 ^ , ka;UJ)Ra(ka, za;UJ) /CO Rg(kg, zg-UJ) V^~1\kg, k; UJ) G{0s(k; UJ) V^(k, ka; UJ) Ra(ka, za; UJ) dk CO /o o R(kg, Zg' UJ) Q-^kg, Zg, Uj) D'^^kg, Zg\k, Z, J Uj) Qj^k, Z B\ Uj) Ra(k, Z a\ Uj) OO • R~x(k, za-UJ) G%(k-uj)R^(k, zg-UJ) • Rg(k, zg; UJ) Q~1(k, zg; UJ) D ( fc , zg\ka, za; UJ) Q7 1(fc 4, za; UJ) Ra(ka, za; UJ) dk /o o D ' ^ " 1 ^ , zg\k, za-UJ) H(fc; UJ) D « ( f c , zg\ka, ZS;UJ) dk . (D.45) CO In general, the term D ' ( N ) removes all events which undergo n interactions with the elastic free surface. This includes both reflected (P-P) and converted (P-S) waves. As 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 in Section 3.6. It is interesting to see what happens when the reference medium becomes water. In that case, the data elements that contain shear components are zero and PP reduces to APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 165 -1. Equation (D.43) then becomes I poo D'iP\kg,zg\k.,z.-,u,) = — / JDS(fc„ lzB|fc,zJ;W)2.-|/Ci<"+")^(* ! lza|fc.>z.;u;)(ifc 27T J-oo (D.46) which is the equation derived by Carvalho (1992) for marine data. APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 166 D.2 Ocean bottom multiple removal The steps in deriving the multiple removal equations for ocean bottom data are essentially the same as for land data. The 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 /CO poo POO POO I I I dxdz dx" dz' •oo J— oo J— oo J— oo Gi{xg, zg\xt z; UJ) V(1)(x, z\x', z'; UJ) P(x",z"\xa,za;uj) 0 0 ips(x",z"\xa,zs;uj) (D.47 ) where the scattering points are located below the source and receivers. The source side Green's functions describe the P and S energy transmitted throught the water bottom by the P source in 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.47) is therefore an idealized data which will serve as the input into the multiple removal algorithm. The 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.13) and Fourier transforming with respect to both shot and geophone domains, (D.47) becomes D^\kg,zg\ks,zs;uj) = e-<"2g*9 2iu2g 0 0 ~'V2gz9 2iV2g J Vpp -"29\ka, v2a; UJ), Vjiy (kg,-v2g\ka, r/2a; UJ) VSP (kg, -V29\ka, v2a; UJ) V^J (kg, ~V2g\ka, rj2a; UJ) PPe —lV2aZwb APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 1 6 7 or expressed as a matrix equation B(kg,zg\k3,z3;uj) = Rg(kg,zg;oj )V^(kg,k3;uj)R3(k3,z3;uj). Obtaining "V^1) is then a simple operation in (kg,k3,uj) space V ( 1 ) ( f c 3 , k.; UJ) = R^ikg, zg; UJ) B(kg, zg\ks, zs; UJ) R^iK, zs; UJ) ( D . 4 9 ) (D.50) Determining the multiple removal subseries The same logic which led to the free surface multiple removal scheme can be used to obtain a multiple removal series for ocean bottom data. The 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 ( 1 ) G o C V ( 1 ) V ' ( 3 ) = - V ' ^ G J ^ V W - V ( 1 ) G £ C V ' ( 2 ) - V W G ^ V W G J ^ V W ( D . 5 1 ) For this first equation, I write V ' ( 2 ) in the same coordinates as V " W in ( D . 4 8 ) and using G J * from ( C . 1 3 ) , I obtain to J_< dk Vpp\kg,-V29\ks, v23, UJ) VP{s\kg, -u2g\k3,T]2S; UJ) Vsp\kg, -V2g\ks, v2a\UJ) Vs{s\kg, -i]2g\ks,r]2a; UJ) VpP (kg, ~V2g\k, v2, UJ), V^J (kg,-u2g\k, 7/2; UJ) VSP (kg, -V29\k, v2; UJ) VQ (kg, -r]2g\k, r)2; UJ) \ R%P/2iu2 RFs/2iri2 V $ (k, -v2\ks, u2s; UJ), Vp^J (k, -u2\k3,rj2s;uj) RTP/2iu2 RTs/^V2 [V£}(k,-Ti2\k.tv2t;u) V^]) (k, —rj2\k3,rj2s;uj) where I have used the spatial integrals as Fourier transforms. This equation can be written in the compact form /CO dkV^(kg, k-uj)G™(k-a,)V«(fc, k3;uj) ( D . 5 3 ) •CO (D.52) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 168 where 2TT RSP/2iv2 Rss/2iV2 (D.54) A similar process can be performed for all the terms in the subseries leading to the recursive relation /CO dkV'W(kg, k; uj)G^(k; uj)V^(k, ks; UJ) . (D.55) • o o Data reconstruction and multiple removal The next step is to reconstruct data using V and directly propagating reference Green's functions. This gives the data series D ' = G * V ' G * (D.56) = G ^ ( V W + V ' ( 2 ) + V ' ( 3 ) + . . . ) G ^ = D ( I ) + f ) ' ( 2 ) + D ' ( 3 ) + . . . The first term is just the 4 component (P, S) data in (D.49). The second term is D'^\kg, zg\k„ zs-UJ) = Rg(kg, zg; UJ) V'^(kg, ks;u>) R3(k3, z3; UJ) (D.57) which after substituting (D.53) and (D.50) becomes £>'(2\kg,zg\ks,zs;uj) = /o o Rg(kg, zg; UJ) V « ( f c g , k; UJ) G^c(k; UJ) V « ( f c , k3; UJ) R3{k3, z3; UJ) dk • C O /OO D(1)(fe„ zg\k, z3; UJ) R ; \ k , z3; UJ) G^c(k; UJ) R~\k, z3; UJ) D™(k, z3\k3, z3; UJ) dk oo /o o B^(kg, zg\k, z3; UJ) J(k; UJ) T>W(k, z3\k3,z3- UJ) dk (D.58) - o o where 3(k\u) 7T Rpp I PPeiu2(Zg+Zw''} Rpcs IP Pei(r,2Z!>+V2Zwb} Rsp I PSei(V2Z<>+r*iZ™b) R-ss/PSeirl^Zs+z^ (D.59) APPENDIX D. ELASTIC INTERFACE MULTIPLE REMOVAL SERIES 169 When the receiver is located on the water bottom, I use C.14 and the reflection coefficients in Appendix B to obtain J(k;u)) = —!_e^(*'-*»») (PP + Z2SS)/(r,2 - k2) 2p2kvlV2(32 (1 + Z2) p2v^\(k2 - V2) (1 + Z2) (PP + Z2SS)/ku2 (D.60) where 0 = 1 + PPZ2 and Z = eivi2z™\ Proceeding with the data reconstruction for all the terms in (D.56), I obtain the recursive relation /CO D ' * " " 1 ^ , zg\k, z3; UJ) J(fc; UJ) D « ( f c , z3\k3, zs; u>) dk. (DM) •CO