UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

An inverse scattering series method for attenuating elastic multiples from multicomponent land and ocean… Matson, Kenneth Howell 1997

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1997-196178.pdf [ 18.78MB ]
Metadata
JSON: 831-1.0087987.json
JSON-LD: 831-1.0087987-ld.json
RDF/XML (Pretty): 831-1.0087987-rdf.xml
RDF/JSON: 831-1.0087987-rdf.json
Turtle: 831-1.0087987-turtle.txt
N-Triples: 831-1.0087987-rdf-ntriples.txt
Original Record: 831-1.0087987-source.json
Full Text
831-1.0087987-fulltext.txt
Citation
831-1.0087987.ris

Full Text

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

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0087987/manifest

Comment

Related Items