UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Aspects of spatial wavelets and their application to modelling seismic reflection data 1986

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
UBC_1987_A6_7 N38.pdf [ 6.44MB ]
UBC_1987_A6_7 N38.pdf
Metadata
JSON: 1.0052976.json
JSON-LD: 1.0052976+ld.json
RDF/XML (Pretty): 1.0052976.xml
RDF/JSON: 1.0052976+rdf.json
Turtle: 1.0052976+rdf-turtle.txt
N-Triples: 1.0052976+rdf-ntriples.txt
Citation
1.0052976.ris

Full Text

A S P E C T S OF SPATIAL W A V E L E T S A N D T H E I R A P P L I C A T I O N T O M O D E L L I N G SEISMIC R E F L E C T I O N D A T A by ATUL NAUTIYAL B.A.Sc, The University of Toronto, 1984 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Geophysics and Astronomy We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 1986 ©Atul Nautiyal, 1986 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 G e o p h y s i c s and Astronomy The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date December 17, 1986 A B S T R A C T The propagation of seismic waves may be described in the space-frequency domain by the Rayleigh-Sommerfeld convolution integral. The kernel of this integral is called a spatial wavelet and it embodies the physics and geometry of the propagation problem. The concepts of spatial convolution and spatial wavelet are simple and are similar to other topics studied by geophysicists. With a view to understanding these concepts, some aspects of spatial wavelets and their application to two-dimensional, zero-offset, acoustic seismic modelling were investigated. In studying the spatial wavelet, two topics in particular were examined: spatial aliasing and wavelet truncation. Spatial aliasing arises from the need to compute a discrete wavelet for implementation on a computer. This problem was solved by using an analytic expression for the spatial wavelet in the Fourier (wavenumber) domain. In the wavenumber domain the wavelet was windowed by a fourth order Butterworth operator, which removed aliasing. This technique is simple and flexible in its use. The second problem of wavelet truncation is due to the necessity of having a wavelet of finite length. A length limiting scheme based upon on the energy content of a wavelet was developed. It was argued that if that if a large portion of the wavelet energy was contained in a finite number of samples, then truncation at that sample would incur a minimal loss of information. Numerical experiments showed this to be true. The smallest length wavelet was found to depend on temporal frequency, medium velocity and extrapolation increment. The combined effects of these two solutions to the practical problem of computing a spatial wavelet resulted in two drawbacks. First, the wavelets provide ii modelling capabilities up to structural dips of 30 degrees. Second, there is a potential for instability due to recursive application of the wavelet. However, neither of these difficulties hampered the modelling of fairly complex structures. The spatial wavelet concept was applied to seismic modelling for media of vary- ing complexity. Homogeneous velocity models were used to demonstrate diffraction evolution, dip limitations and imaging of curved structures. The quality of modelling was evaluated by migrating the modelled data to recover the time-image model of the reflection structure. Migrations of dipping and synform structures indicated that the modelled results were of a high calibre. Horizontally stratified velocity models were also examined for dipping and synform structures. Modelling these reflection structures showed that the introduction of a depth variable velocity profile has a tremendous influ- ence on the synthetic seismic section. Again, migration proved that the quality of the data was excellent. Finally, the spatial wavelet algorithm was extended to the case of laterally varying velocity structures. The effects of space variant spatial convolution in the presence of a smoothed velocity field were examined. Smoothed velocity fields were computed by a simple weighted averaging procedure. The weighting function used was a decaying exponential whose decay rate determined the amount of smoothing. Seis- mograms computed for this case showed that the algorithm gave smoother and more continuous reflection signatures when the velocity field has been smoothed so that the largest lateral velocity gradient corresponded to the lower end of the temporal frequency band of the spatial wavelets. In this respect, the results are similar to those of geometric ray theory. Also, the traveltimes of these models compared favourably with those of ray tracings iii T A B L E OF C O N T E N T S ABSTRACT ii LIST OF FIGURES vi ACKNOWLEDGMENTS viii DEDICATION ix 1. INTRODUCTION 1 1.1 Seismic reflection modelling . . 1 1.2 Representation of seismic reflection data 2 1.2.1 The common midpoint stack 3 1.2.2 The exploding reflector hypothesis 7 1.3 A review of modelling methods 10 1.3.1 Analogue methods 10 1.3.2 Analytic methods 11 1.3.3 Numerical methods 12 1.3.4 Transform methods 13 1.3.5 Hybrid methods 14 1.4 The spatial convolution approach 14 1.4.1 Convolution with spatial wavelets 15 1.4.2 Recursive versus nonrecursive extrapolation 19 1.5 Scope of this study 21 2. ASPECTS OF SPATIAL WAVELETS 23 2.1 Theoretical development 23 2.1.1 The Kirchhoff-Helmholtz integral 23 2.1.2 The Rayleigh-Sommerfeld integral 27 2.1.3 Seismic approximations and limitations of theory 33 2.2 Spatial aliasing 33 2.2.1 Analysis by asymptotic approximations 34 2.2.2 Analysis by Fourier transform 47 2.2.3 Frequency spectra 57 2.3 Wavelet truncation 65 2.3.1 Wavelet energy level 66 2.3.2 Minimum aperture wavelet 67 2.3.3 Seismograms, spectra and stabilty 72 3. MODELLING SEISMIC REFLECTION DATA 80 3.1 Homogeneous acoustic media 80 3.1.1 A truncated model 80 i v 3.1.2 Dipping models 84 3.1.3 Antiform and synform models 93 3.2 Horizontally stratified acoustic media 100 3.2.1 A dipping model 100 3.2.2 A synform model 103 3.3 Laterally varying acoustic media 106 3.3.1 Accomodating lateral velocity changes 107 3.3.2 A simple experimental model 109 4. CONCLUSION 119 4.1 Summary 119 4.2 Future investigations 124 REFERENCES 126 APPENDIX Butterworth windows and windowing strategy 131 v LIST OF F IGURES 1-1. Field geometries for seismic data acquisiton 4 1-2. Common midpoint gather geometry 6 1-3. Exploding reflector model for zero-offset data 8 1-4. Physical interpretaion of spatial convolution 16 1-5. Forward extrapolation in the space-frequency domain 18 1- 6. Recursive and nonrecursive extrapolation 20 2- 1. Geometry for derivation of Kirchhoff-Helmholtz integral 25 2-2. Geometry for derivation of Rayleigh-Sommerfeld integral . . . . . . . 28 2-3. Seismograms satisfying aliasing conditions 40 2-4. Seismograms violating frequency limit 42 2-5. Seismograms violating step size limit 45 2-6. Seismogram violating frequency and step size limits 46 2-7. Complex spectrum of spatial wavelets at 30 Hz 49 2-8. Complex spectrum of spatial wavelets at 50 Hz 51 2-9. Seismogram violating frequency limit with bandlimited wavlets 54 2-10. Seismogram violating step size with bandlimited wavelet 55 2-11. Seismogram violating frequency and step size limits with bandlimited wavelets 56 2-12. Frequency spectra for seismograms in 2-3 59 2-13. Spectral comparison of space and wavenumber domain wavelets . . . . 61 2-14. Frequency spectrum for seismogram in 2-9 62 2-15. Frequency spectrum for seismogram in 2-10 63 2-16. Frequency spectrum for seismogram in 2-11 64 2-17. Wavelet energy variation with extrapolation step size 69 vi 2-18. Wavelet energy variation with frequency 70 2-19. Wavelet energy variation with velocity 71 2-20. Seismograms computed using band and aperture limited wavelets . . . 73 2-21. Frequency spectra for seismograms in 2-20 77 2- 22. Wavenumber spectra for seismograms in 2-20 78 3- 1. A truncated model 81 3-2. Model with a 15 degrees dip . 85 3-3. Model with a 30 degrees dip 87 3-4. Model with a 45 degrees dip 89 3-5. Model with a 60 degrees dip 91 3-6. Antiform model 94 3-7. Synform model with d < r 96 3-8. Synform model with d — r 97 3-9. Synform model with d > r 98 3-10. Model with a 30 degrees dip in horizontally layered media 101 3-11. Synform model with d > r in horizontally layered media 104 3-12. Waves generated at a velocity discontinuity 110 3-13. Seismogram for a lateral velocity discontinuity 112 3-14. Seismograms for smooth lateral velocity transitions 115 A - l . Dependence of Butterworth characteristics on the order . . . . . . . 133 A-2. Dependence of Butterworth characteristics on recursion 134 A-3. Seismograms for a varying cutoff, fourth order windowing strategy . . . 137 A-4. Seismograms for a fixed cutoff, fourth order windowing strategy . . . . 138 A-5. Seismograms for a fixed cutoff, box-car windowing strategy 139 vii A C K N O W L E D G E M E N T S There are many people who made the completion of this study possible and I would like to express my gratitude. I thank my parents, Malti and Jagdish, and my little sister, Nandita, for their never-ending supply of love and their well-timed shipments of German coffee cakes. It is a pleasure to acknowledge my advisor, Doug Oldenburg, for his unwavering enthusiasm for this work and for his many helpful suggestions. Equally important, I thank Doug for his friendship. Tad Ulrych provided me with inspiration when I was in search of a thesis topic. Indeed, this study originated from a paper I wrote for a course with him. Matt Yedlin helped me with numerous aspects of the mathematics and the physics of wave theory. My interest in seismology was spurred by Chris Chapman while I was an undergraduate student at the University of Toronto. I thank him especially for having given me an opportunity to learn something new. Needless to say, my fellow graduate students and I shared in the agony of our research and in the thrill of our eventual successes. You know who you are darlings and you are marvelous. In particular, I thank Tim Scheuer and Ken Whittall. Also, Sonya Dehler kindly allowed me to pirate her T^X files which eased the task of thesis production. To my mates at Barbarian House, Erick Baziw, Erik Blake and Gord Williams, thanks for feeding and clothing me. Words are not enough. This project was supported by Natural Sciences and Engineering Research Council grant A4270 to Doug Oldenburg and by a scholarship to myself from the Canadian Society of Exploration Geophysicists and Pembina Resources Limited. viii D E D I C A T I O N Rachna, j te amo ! ix 1. I N T R O D U C T I O N 1.1 Seismic reflection modelling The interpretational quality of seismic reflection data is largely dependent on the skill and experience of the interpreting geophysicist or geologist. In petroleum exploration during the past, this task might have been considered more straightforward than it is today. This is perhaps because the obvious oil and gas fields were easy to locate, be- ing large in spatial extent as well as being found in relatively uncomplicated geology at shallow depths. As these known reserves dwindle and the need for petroleum re- sources increases continuously, geophysicists and geologists are challenged by having to probe new regions of this planet where the geology of oil and gas bearing rocks may be very complex and the deposits are of small size. With exploration extending to areas of more complicated geology, there arises the difficulty of interpreting correspondingly complicated reflection seismograms. For example, the seismic responses in thrust belt and fault zones are extremely difficult to understand, even for the seasoned interpreter. Therefore, it becomes highly desirable to simulate the processes of seismic wave propa- gation and data collection under the controlled conditions of a laboratory (which may be a computer) and thus gain insight by the investigation of known situations. This line of thought leads to the concept of seismic reflection modelling for the purpose of inferring subsurface geologic structures. The results of modelling seismic wave propagation can be used in many ways by the geophysicist to understand and predict seismic wave phenomena. These uses fall into the main categories of (l) illustrating the consistency of interpretation of real data; 1 (2) providing synthetic data sets for testing processing methods, acquisition parameters and new concepts; and (3) educating geophysicists in wave propagation (Kelly et al., 1982; Berkhout, 1984). One of the most obvious uses of modelling is in the creation of synthetic data sections for a proposed model for comparison with the observed data. It is unreasonable to expect the synthetic data to match the real data precisely, although skilled interpreters can, by the degree of similarity, increase their confidence in the proposed geologic model. The insight given by the modelling procedure often allows the geophysicist to refine the interpretation. This process may be repeated until an acceptable match is obtained between the simulated data from the model and the true data. 1.2 Representation of seismic reflection data The fundamentals of seismic reflection data acquisition and processing have been dis- cussed in numerous texts. Dobrin (1976), Waters (1978), Sheriff and Geldart (1982a), and McQuillin et al. (1984) have dealt in depth with acquisition, while Robinson and Treitel (1980), Kanasewich (1981), Sheriff and Geldart (1982b), Claerbout (1985a) and Hatton et al. (1986) have provided details of processing. Instead of presenting a thor- ough review of the contents of these references I offer a brief synopsis of the essential points that are needed later for studying aspects of synthetic seismograms. Concise, yet complete, descriptions of the salient material may be found in Robinson (1983, p.17-33) and Gazdag and Sguazzero (1984). I summarize these two works in the remainder of this section. 2 1.2.1 The common midpoint stack In land seismic reflection surveys two types of signal emitting sources are commonly employed. First, there is dynamite which produces a short, intense pulse. Second, there are vibratory sources such as Vibroseisf which send a long duration (7 to 35 seconds) sinusoidal signal into the earth. The latter type of source has come into widespread use over the last decade. Velocity sensitive detectors at the earth's surface, called geophones record the signals due to the seismic waves reflected from the subsurface geology. At sea, an airgun source is used to emit an explosive signal into the earth and pressure sensitive hydrophones, placed in a cable, record the information carried by reflected waves. The sampled time series data associated with a single source and receiver location is called a seismogram, a seismic trace or more simply a trace. There are various geometries for shot and receiver which are employed in field work. The basis for all coordinates in field geometry are the locations of the source and receiver. Let s and r represent the horizontal coordinates of the source and receiver respectively. Also define two other descriptive coordinates, m — (s + r ) /2 and h — (r — s)/2, which are known as the midpoint and half-offset coordinates respectively. Note that m and h define a set of axes rotated 45 degrees with respect to the axes r and s. A simple illustration of the different geometries in the field procedure is shown in Figure 1-1. Although the scene depicts a marine application, the principles apply on land as well. Each point is representative of a seismic trace for a standard survey. The ship moves in the x direction along a straight line and tows a source and a cable of hydrophones. After the vesselmoves one-half receiver interval, a shot is fired and the | Trademark of Conoco Inc. 3 F i g u r e 1 - 1 . R e l a t i o n s h i p a m o n g the ho r i zon ta l coord ina tes r , s , m a n d h. A l l axes represent d i s tances o n the p a t h of the sh ip w h i c h moves i n the +x d i r e c t i o n . 1-1 is a c o m m o n hal f -of fset ga the r , 2 - 2 is a c o m m o n m i d p o i n t ga ther , 3 - 3 is a c o m m o n source g a t h e r a n d 4 - 4 is a c o m m o n receiver ga ther . E a c h po in t represents a poss ib le se ismic t r ace (af ter R o b i n s o n , 1983, p.22). 4 hydrophones along the cable record the pressure. The geometry indicates that there are four types of common point gathers of traces. (1) Common half-offset point (CHP) gather: sets of traces with the same half-offset coordinate, h. (2) Common midpoint (CMP) gather: sets of traces with the same midpoint coordinate, rn. (3) Common source point (CSP) gather: sets of traces with the same source coordinate, s. (4) Common receiver point (CRP) gather: sets of traces with the same receiver coordi- nate, r. The current commercial standard for subsequent processing is the CMP gather. There are many processing procedures involved in constructing an image of the subsurface from the acquired seismic data, but I limit the discussion to stacking and normal moveout (NMO) corrections. Stacking consists of the summation of the traces of each CMP gather after correcting them to compensate for the offset between source and receiver. This compensation is known as the NMO correction and the process is illustrated in Figure 1-2. The two-way traveltime of the signals from a horizontal reflector is given by m=ti+%, (1.1) V J . where the velocity t; above the reflector is constant. In the case where there are many horizontal layers of differing velocities, t; in (1.1) is replaced by the root-mean-square velocity. The first arrival at h = 0 has the reflection time, to- As can be seen, (1.1) 5 Figure 1-2. (a) Common midpoint gather geometry indicates that sources are located to the left of the midpoint, m and receivers are to the right, (b) Data display in the form of a space-time (x - i) section shows that the locus of the reflectors from the given horizontal reflecting plane is a hyperbola (after Robinson, 1983, p.26). 6 describes a hyperbola. The difference between the reflection time observed at an offset 2h and that corresponding to zero-offset is the NMO correction and is given by A ^NMO = t{h) - t0. (1.2) The time shift defined by (1.2) is then applied to the signals with the result that they appear as if they were recorded with coincident source and receiver, i.e., a zero-offset gather. The NMO corrected traces are then stacked. The collection of stacked traces along the midpoint axis is called a CMP stacked section. Stacked data have the advan- tages of reduced random noise and attenuation of multiple reflections in comparison to unstacked data. 1.2.2 The exploding reflector hypothesis A CMP stacked section may be considered a zero-offset section in which the path trav- elled by a seismic ray from source to reflector is identical to that from reflector to receiver (Figure l-3a). In this model the assumptions are that all the shots are fired simultaneously and that each receiver records only the signal from its coincident source. However, this view does not correspond to data resulting from a single seismic experi- ment. To hypothesize a physical experiment which would produce the zero-offset data Loewenthal et al. (1976) proposed an exploding reflector model (ERM) in which the en- ergy sources are not at the surface, rather they are located along the reflecting surfaces (Figure l-3b). That is, the reflectors are represented by buried sources, which are set off at the same time t = 0. With such a model, one need only be concerned with upward travelling waves. Since a record section involves two-way time, it must be converted to 7 Figure 1-3. (a) The common midpoint section may be regarded as data obtained with coincident (zero-offset) sources and receivers, (b) The exploding reflector model of the common midpoint data for modelling sections replaces the reflector surface by buried sources (after Gazdag and Sguazzero, 1984). 8 one-way time. In practice, the time scale of the C M P section is left unchanged while the wave propagation velocity is scaled down by a factor of two. This is known as the half-velocity substitution. There are some important limitations to E R M which should be discussed. To begin with, the concept is not applicable to nonzero-offset geometries and extensions for this case have yet to be formulated. Therefore, its use is restricted to the C M P stacked section. In fact, E R M has shortcomings for zero-offset sections which make it quanti- tatively incorrect. In particular, there are three obvious failings in E R M (Claerbout, 1985b, p. 10). These are reiterated below. First, waves travelling through media with strong lateral velocity variations cannot be modelled, although they appear on zero-offset sections. That is, E R M cannot account for ray bending except if the velocity structure is horizontally stratified. Second, E R M fails in its prediction of multiple reflections. For example, a flat sea bottom causing a primary reflection at a two-way traveltime t\ will generate multiples at times of 2 £ 1 , 3 r 1 , 4 r 1 , etc. According to E R M the first multiple goes from sea bottom to surface, then from surface to sea bottom, then from sea bottom to surface once again, for a total traveltime of 3r,. Other multiples follow at times 5ti,7ti, etc. Thus, multiple reflections generated on a zero-offset section are not correctly forcast by E R M . A third drawback of E R M is in its prediction that waves emitted in opposite directions (upgoing and downgoing) have the same polarity. However, the physical concepts of determining reflection coefficients implies that the downgoing waves will see reflection coefficients of opposite polarity to those seen by upgoing waves. 9 1.3 A review of modelling methods There are numerous techniques for modelling the propagation of seismic waves. These methods span the spectrum from analogue modelling in laboratories, to evaluation of analytical solutions of the wave equation, to numerical modelling of the wave equation or its variants using computers. I offer a brief overview of these modelling techniques. 1.3.1 Analogue methods Laboratory seismic models can provide a bridge to close the gap between theoreti- cal studies of wave motion in the simplest of geometries and the more complicated behaviour of waves in realistic situations. Scaled seismic models were the subject of in- tensive research during the 1950s and a large variety of problems were solved by model techniques. Despite the evolution of other approaches, scale models may still present the closest approximations to what may be considered the true solution. There are three categories into which seismic models may be classified depending on the number of dimensions used. A one-dimensional (l-D) model may be built in various forms such as the sound tube used by Woods (1956) or the plastic rod described by Berryman et al. (1958a, 1958b). The l-D model is ideal for the analysis of many reflectors and the problem of associated ghosts and multiples which can be described by the convolutional model of Robinson (1957). This assumes that waves have normal incidence on flat, horizontal reflection surfaces. The two-dimensional (2-D) model was introduced by Oliver et al. (1954) and used later by Angona (1960) and Kosloff and Baysal (1982). In this procedure thin sheets of elastic materials are cut into the form of discs and strips and bonded together to form 10 a layered cross-section. In geologic regions where the stratigraphy does not change ap- preciably in the horizontal direction perpendicular to the seismic line, the 2-D approach may be applied without loss of generality. The third type of seismic model is used for describing three-dimensional (3-D) ge- ometry and has been detailed by Howes et al. (1953), Evans et al. (1954), Hilterman (1970) and French (1974). Construction of such models for all but the most simple of geometries is more difficult than for the 1-D and 2-D cases. Also, the materials that may be used for the construction are limited since few materials are sufficiently homogeneous in their constitution. The extension to 3-D is most useful when the stratigraphy varies significantly in the direction perpendicular to the seismic line. 1.3.2 Analytic methods In exploration seismology, analytic methods have essentially concentrated upon the Kirchhoff-Helmholtz integral solution to the wave equation. A first order approximation to the Kirchhoff-Helmholtz description of wave propagation is geometrical ray tracing. However, ray tracing alone cannot explain the common phenomenon of diffraction. Keller (1962) extended the ray concept with the development of the geometrical theory of diffraction. The application of extended geometric ray theory to seismology was made by Cerveny et al. (1977) who examined wave propagation in inhomogeneous media. With this approach, seismograms can be computed which predict the travel times and approximate the amplitudes of the signals travelling through complex models. The Kirchhoff-Helmholtz approach has enjoyed considerable popularity in industrial applications. Trorey (1970) used this method to study seismic diffractions and Hilter- man (1970) used it to explain results from 3-D scale modelling experiments. Both these 11 works considered coincident source-receiver geometry and it was not until the research of Berry hill (1977) and Trorey (1977) that extensions for arbitrary source-receiver ge- ometries were derived. Hilterman (1975) found that the Kirchhoff-Helmholtz theory was useful in explaining seismic amplitude anomalies such as bright spots. The con- tributions of Hilterman and Larson (1975), Berryhill (1979), Carter and Frazer (1983), and Deregowski and Brown (1983) extended this method to the case of laterally varying velocity. 1.3.3 Numerical methods The need to interpret seismic responses in increasingly complex geologic media necessi- tated the evolution of more complicated methods which could simulate the reflection pro- cess. Numerical modelling of the wave equation and its variants by the finite-difference method (FDM) was developed during the late 1960s and early 1970s. Alterman and Karal (1968) used FDM to study waves in an homogeneous, isotropic, elastic half-space overlain by a layer with similar properties. Diffracted waves due to quarter and three- quarter planes were studied by Alterman and Loewenthal (1970). Later, Boore (1972) investigated applications of FDM to heterogeneous media. The accuracy of FDM was analyzed by Alford et al. (1974) by comparing FDM results with analytic means for simple geometries. These investigators found the results to be similar. Kelly et al. (1976) contributed the first discussion on FDM applications to exploration seismics and interpretation of sections. In a review of FDM, Kelly et al. (1982) demonstrated its use in complicated acoustic and elastic media and through the various stages of acquisition and processing. Recently, Kelamis and Kjartansson (1985) extended FDM to laterally varying acoustic velocity structures using a space-frequency domain implementation. 12 The use of the frequency domain allows the modelling of certain frequency dependent phenomena, such as attenuation, dispersion and spatial bandwidth. Another numerical approach is the finite-element method (FEM). Bamberger et al. (1980) and Marfurt (1984) have analyzed the accuracy of both FDM and FEM and compared them. They found that in homogeneous media, FEM and FDM are comparable when 0.30 is the upper limit on Poisson's ratio. In contrast, where Poisson's ratio is between 0.30 and 0.45 the FEM approach is superior. This is an important factor to consider when looking for petroleum indicators such as the ratio of compressional and shear wave velocities (Vp/Vs) (Tatham and Stoffa, 1976). It has been noted that there is a relationship between Vp/Vs and Poisson's ratio, so that using either FEM or FDM will have consequences on the interpretation. 1.3.4 Transform methods Fourier transform methods of modelling the wave equation have also been used. Gazdag (1981) used the frequency-wavenumber domain to model primary compressional reflec- tions in homogeneous and horizontally stratified media. In the frequency-wavenumber domain propagation of waves may be accomplished by a multiplication of the source wavefield by a pure phase operator. The time-wavenumber domain has been used by Kosloff and Baysal (1982) for application in laterally varying velocity structures. These investigators Fourier transformed the wave equation from the space domain to the wavenumber domain and used the FDM to propagate waves forward in time. This method was extended to the case of elastic wave computations by Kosloff et al. (1984) and was applied to certain reflection interpretation problems by Reshef and Kosloff (1985). 13 1.3.5 Hybrid methods Each group of methods has its advantages and limitations, so that the choice of method used for computing synthetic seismograms is dependent on the specific situation. It is for these reasons that Shtivelman (1984, 1985) has developed a hybrid method of analytical (ray tracing) and numerical (FDM) techniques for acoustic media. 1.4 The spatial convolution approach An interesting approach to the modelling of seismic wave propagation was suggested by Berkhout and van Wulfften Palthe (1979, 1980). They proposed that waves could be propagated through a medium by a recursive series of spatial convolutions or de- convolutions in the space-frequency domain by using the principles of Kirchhoff scalar diffraction theory. The deconvolution procedure represents a fundamental step in the seismic imaging technique known as migration. A comprehensive discussion of deconvo- lution in the light of linear inverse theory has been given by Oldenburg (1981), although this approach has not yet been applied to the migration problem. In contrast, spatial convolution is an expression of the forward problem of seismic modelling. The spatial convolution approach is novel in that it may be able to accomodate complicated velocity structures by a series of recursive extrapolation steps. Berkhout (1981, 1982, 1984) did much to develop the theoretical concepts for both seismic modelling and migration, with an emphasis on the latter, using this approach. In this section I introduce the concepts involved in the spatial convolution method for seismic reflection modelling of zero-offset data. 14 1.4.1 Convolution with spatial wavelets An important result from scalar diffraction theory is that for monochromatic pressure waves, propagation may be cast as a convolutional process over the space domain. For the most general geometries this can be expressed in terms of the monochromatic Kirch- hoff-Helmholtz integral. Consider the simplifying assumption that the geometry of the source and recording surfaces is planar, then one arrives at a variant of the Kirchhoff- Helmholtz integral, known as the Rayleigh-Sommerfeld integral. A futher simplification is to assume two-dimensional geometry so that all source-receiver geometries and ma- terial properties are considered invariant along one of the space coordinates. The result is that the propagation convolution may be expressed by P{x0,zl+1,fn) = / W(x0- x,Az,fn)P(x,zl,fn)dx. (1.3) The quantity P represents the pressure field of the waves, x0 indicates the locations of receivers on the recording plane at a depth level Zi+i, x indicates the location of emitters on the source surface at a depth level 2 t and fn is the nth frequency component. When the wave is composed of many frequencies, then the convolution over space must be done for each frequency. Assume a monochromatic point source of unit amplitude at x = 0 to be the wavefield on depth level Z{, then the response at depth level zt'+i will D e given by W(x, Az, f n ) . Therefore, W(x, Az, fn) is the the spatial impulse response, or the spatial wavelet, for the temporal frequency component fn and extrapolation step size Az. The spatial wavelet is a symmetric or zero-phase function about the point source when the medium is homogeneous. I present a rigorous account of scalar diffraction theory in the next chapter. 15 I 1 1 I I I I I I 1 I b Figure 1-4. Physical interpretation of spatial convolution, (a) The wavelet w gives the response of a linear array of sensors to a point source, (b) The wavelet gives the value of the wavefield at a single sensor due to a linear array of point sources. 16 I interpret the mathematical operation of convolution for wave propagation with a simple physical model. Consider a point source embedded in a homogeneous medium and an array of sensors on a recording surface. When the source produces an acoustic emission, the wavefield radiates outwards and is detected by recorders. The portion of the wavefield received at a particular geophone is determined by the spatial wavelet. Thus, the spatial wavelet can be regarded as a weighting function that governs the distribution of the wavefield, due to a point source, everywhere in space (Figure l-4a). Since the sensors are located somewhere in this space, the spatial wavelet determines the wavefield values at the sensor array. Consequently, the spatial wavelet can also give the total response at a geophone due to a collection of buried point sources (Figure l-4b). Determining the wavefield at each receiver involves sliding the spatial wavelet over the collection of sources and summing up the source contributions according to the weighting function. Thus, convolution computes the wavefield on a recording surface by using linear combinations of the wavefield on the source surface. The situation described above was for a single frequency. In reality, signals are composed of many frequencies and each frequency requires its own spatial wavelet for a single extrapolation step. In the space-frequency domain a one-dimensional convo- lution is involved along the space axis for each frequency component. Inverse Fourier transformation from the frequency domain to the customary time domain for each trace completes the procedure. Thus, the space-time seismic section is obtained. This is illustrated in Figure 1-5. 17 *, K3 i f, f 3 T — X — • J <« f T a f N 1 — b t y G Figure 1-5. Forward extrapolation in the (a) space-frequency domain involves (b) reordering from traces to frequency components and one-dimensional convolution along the x axis for each frequency then (c) a one-dimensional inverse Fourier transformation (/ —* t) to obtain the space-time section (after Berkhout, 1982 p.138). 18 1.4.2 Resursive versus nonrecursive extrapolation Wavefield continuation or extrapolation may be done in two ways: recursive and non- recursive. In recursive extrapolation the output of the previous extrapolation step is used as input for the next extrapolation step. In contrast, for nonrecursive extrapola- tion the input for each extrapolation step is the wavefield on each source plane and the continuation step is always directly to the observation plane (Figure 1-6). When the velocity in the medium is constant, the nonrecursive approach has the ad- vantage over the recursive method in that it can compute the seismogram with greater accuracy. However, as soon as the next level of complexity, the case of horizontally stratified media, is introduced the recursive method immediately gains the upper hand. Recursive computation implies that each extrapolation takes place in a confined portion of a homogeneous medium whose velocity differs from the surrounding strata. Thus, vertical velocity variations may be accounted for with greater ease than with the non- recursive method. In addition, when the more complicated case of laterally varying velocity structure is considered, only recursive application can provide a possible means of reflection modelling It appears that the best approach to extrapolate is by recursive means if compli- cated velocity structure are to be accomodated. However, the points made above apply in the space-frequency domain and implementation in other domains may introduce dif- ferent advantages and disadvantages for these extrapolation schemes. I refrain from a discussion of these, as I am only interested in an inquiry regarding the space-frequency application. 19 Figure 1-6. Recursive (indicated by R) versus nonrecursive (indicated by NR) wavefield continuation. The recursive approach can model complex structure more easily than can the nonrecursive method. 20 1.5 Scope of this study The spatial wavelet approach is appealing to geophysicists as it expresses the complex phenomena of wave propagation as a simple one-dimensional convolution. Therefore, to understand aspects of wave propagation one need only gain some knowledge of the spatial wavelet and the operation of convolution. Berkhout (1981, 1982, 1984) has discussed the spatial wavelet concept in considerable depth, although those works deal primarily with the theoretical aspects of seismic migration. However, to date there have been no explicit studies in the literature which demonstrate the application of the spatial wavelet to migration or modelling. Here, I investigate the spatial wavelet concept from a different perspective and apply the method to the forward problem of seismic reflection modelling. Numerical examples in the form of synthetic reflection seismograms demonstrate the application of this approach. In the second chapter of this thesis I give an in depth discussion of the spatial wavelet and some of its aspects. I begin with the fundamentals of Kirchhoff scalar diffraction theory and show that the spatial wavelet is the kernel of the Rayleigh-Sommerfeld integral. I then address two practical problems in wavelet computation and provide original solutions which overcome these difficulties. The first problem, that of aliasing the spatial wavelet, is pursued using analyses by asymptotic approximations and by Fourier transformation. A new spectral domain method of eliminating spatial aliasing is proposed as a result of these examinations. The second difficulty in the practical computation of synthetic seismograms is the problem of wavelet truncation due to a finite number of observations. I analyse truncation effects based on the consideration of the energy contained in a truncated spatial wavelet. It is demonstrated that when a 21 finite number of samples contains a large portion of the wavelet energy, then errors due to truncation after that number of samples may be suitably minimized. The spatial wavelet method is evaluated in chapter three. Performance is judged by using synthetic, zero-offset seismic sections which are computed for earth models with velocity structures of varying complexity: homogeneous acoustic media, horizon- tally stratified acoustic media and laterally varying acoustic media. To the best of my knowledge these results are the first using the the spatial wavelet approach as I have described it. In the fourth and final chapter I conclude with a summary of the investigation and give some suggestions for further inquiries into the topic of pressure wave propagation using spatial wavelets in the space-frequency domain. 22 2. A S P E C T S OF SPATIAL WAVELETS 2.1 Theoretical development The spatial convolution approach is essentially a Kirchhoff-Helmholtz technique. Tere- fore, I review the classical diffraction theory of Kirchhoff which provides the solution to the scalar wave equation using a point source. This particular aspect of the theory can be made applicable to seismic reflection analysis by modelling reflections as the re- sponse due to an alignment of point diffractors. Although the exploding reflector model is rather simplistic, it allows analyses of seismic reflection data that otherwise might not be possible. Kirchhoff scalar diffraction theory for optics has been discussed in great depth in the scientific literature (Goodman, 1968, p.30-56; Born and Wolf, 1975, p.370-458). Therefore, I review only the salient points in the derivation of the Kirchhoff-Helmholtz integral and its variants in two and three dimensions. 2.1.1 The Kirchhoff-Helmholtz integral I begin the discussion of the three-dimensional case with the second theorem of Green: In the above equation G = G(x, y,z,u) and P — P(x,y, Z,UJ) may be any two harmonic complex valued functions. S is a surface surrounding the volume V, n is an outward (2.1) 23 pointing unit vector with respect, to S , u> is the angular frequency and x,y,z are spatial coordinates. Consider a pressure wave disturbance p(x,y,z,t) with some source distribution and an observation point ro surrounded by an arbitrary, source free, closed surface S (Figure 2-1). Then p(x,y,z,t) satisfies the inhomogeneous wave equation everywhere in space. That is, 1 d2 N V 2 P - - - - 4 7 T Y, Sj{t) 6(x - Xj)6{y - yj)8(z - Zj). (2.2) V 3 = 1 In the above, Sj(t) is the jth source strength function, v is the medium velocity and 6 is the symbol for the Dirac delta function. For the closed surface the homogeneous equation applies ^2 1 d2p , . V * P - - J = 0 . (2.3) Fourier transformation with respect to time to equation (2.3) gives the Helmholtz equa- tion V 2 P + k2P = 0- (2.4) where k is the angular wavenumber. The application of Green's theorem to (2.4) requires that G(x,y, z,u>) also be de- fined. I select the free space Green's function for a monopole source G= , (2.5) 24 Figure 2-1. Geometry for derivation of Kirchhoff-Helmholtz integral. The pressure field at r 0 may be computed by knowing the pressure field on the surface S which surrounds the closed volume V. 25 which is a solution of V 2 G - k2G = -4n 6{x - x0)6{y - y0)6{z - z0). (2.6) Here r = \J(XQ — x)2 + (t/o — y ) 2 + (ZQ — z)2 and i = \f—l. Although other choices for G are possible (Kuhn and Alhilali, 1976), I restrict my study to the use of (2.5). Substitution of (2.4) and (2.6) into (2.1) gives j (PVG- GVP) •ndS = -Air j P(r) 6(x - x0)6{y - y0)6{z - z0) dV, (2.7) where I have chosen to write P(x,y,z,OJ) as P[r). The sifting property of the delta functioni yields / ( PVG — GVP) • n dS = -A%P{r0). (2.8) J s Substituting (2.5) into (2.8) results in If d e~ikr dP(r) e~ikr P ^ = - h n p { r )  + V i d s - ( 2 - 9 ) 47r Js dn r dn r Expression (2.9) is the Kirchhoff-Helmholtz integral and it is a solution to the one- way wave equation governing the radiation of waves away from their source. Through this equation it is possible to compute the pressure field at any interior point r 0 by evaluating the surface integral. However, the solution given by the Kirchhoff-Helmholtz integral is uniquely determined by specifying either P (Dirichlet boundary condition) or j £ (Neumann boundary condition) on the surface 5; but both should not be specified because they will not, in general, be consistent with one another. The inconsistency 26 presented by the Kirchhoff-Helmholtz integral boundary conditions was removed by Sommerfeld who eliminated the requirement of imposing boundary conditions on both P and on 2.1.2 The Rayleigh-Sommerfeld integral The imposition of boundary conditions on the pressure field and its directional derivative may be circumvented by considering an extended Green's function given by G= + G*. (2.10) r where G* has no singularities in and on surface S and is to be determined. The extended Green's function is still a solution to equation (2.6) and had it been used, it would give an alternate Kirehhoff-helmholtz integral: If d e~ikr dPlr) e"lkr Piro) — - jf | P(r) _ (_ + + - JU I— + G*)\dS. (2.11) For the seismic reflection case, ^ is not measured directly and therefore one would like to remove it from this expression. This is accomplished by choosing G* such that —— + G* =0 (2.12) on S. To do so, choose a closed surface S defined by the plane 2 = 0 and a hemisphere in the upper half-space z > 0 (Figure 2-2). 27 Figure 2-2. Geometry for derivation of Rayleigh-Sommerfeld integral. The pressure field at ro may be computed knowing the pressure field on the surface S\ which bounds the closed volume V along 2=0. The pressure field on £>2 does not contribute to the solution when the Sommerfeld radiation condition is invoked. 28 With this choice of S one may rewrite (2.11) as P ( R ° ) = - 7 - / P{r)£(—+G-')dS1-j-f P[r)£-(— + G*)dS2. (2.13) 47r Js an r 4n Js^ on r Application of the Sommerfeld radiation condition allows one to drop the integral over S2 as its contribution is naught. One is then left with dn — — dz and (2.13) becomes P(r°) = " i / PW^C-^+ (2.14) J Si Select G" as -ikr' G* = - — (2.15) such that r' defines a point that is the mirror image through the z 0 plane of the point defined by r 0 , i.e., r' = yj(x0 _ x)2 + (no ~ u)2 + (zo + z)2- Clearly, the extended' Green's function is zero on S for this choice of G*. The derivative of the extended? Green's function is dz[ r r> > r* r ' ' where Az — ZQ — z. Therefore, (2.14) becomes Equation (2.17) is the Rayleigh-Sommerfeld integral. This implies that if P is known everywhere on S\, then the pressure field P(ro) at a point r 0 away from S\ 29 (2.16) (2.17) may be determined. This result defines the equation for computing upward continued wavefields in three dimensions. Until this point the discussion has concerned the theoretical developments for three dimensions (3-D). In reflection seismology it is standard to acquire data along seismic lines rather than seismic grids. For this case it is customary to consider the data to be independent of a lateral coordinate, either x or y, in order to make interpretation easier. This study assumes independence in the y direction to make the transformation to two-dimensional (2-D) geometry. The result is an expression that may be used for seismic lines in the x direction. Consider the surface S\ in Figure 2-2 to be of infinite extent, then equation (2.17) may be rewritten as (2.18) To transform (2.18) to 2-D consider that d e~l dz r — ikr = -Az ikr (2.19) which is just a result of (2.16). Then (2.18) may be written as (2.20) 30 Interchanging the derivative and the integral over y gives pM = hLkpwl,cr*'i*- (221) where the pressure function has explicitly been written as independent of y so that P(r 0) becomes P(XQ) and P(r) becomes P(x). Letting the integral over y range from — oo to oo leads to the Heine integral (Magnus and Oberhettinger, 1949, p.27) / oo — ikr dy= -t7xH^](krx), (2.22) -oo r where rx = \J(xo — x)2 + Az2 and HQ is the Hankel function of the second kind and zeroth order. Equation (2.21) becomes P(*o) = -%- / ^Hf\krx) P(x) dx. (2.23) 2 Jx dz Evaluating the derivative using Kinsler et al. (1982, p.450) gives „, . ikAz f H\2)(krx) . , ^ ( * o ) = — / 1 1 X > P(x)dx, (2.24) ^ Jx rx (2) where H\ is a Hankel function of the second kind and of the first order. The range of integration for x represents the size of the aperture in the diffraction slit problem of classical optics. Equation (2.24) is the 2 - D Rayleigh-Sommerfeld integral and it will be the basis for computations. 31 Equation (2.24) may be cast as a convolution where the kernel of the 2-D Rayleigh- Sommerfeld integral is defined by W(x) ikAz HJ2)(kVx2 + Az2) (2.25) 2 sjx2 + Az2 and depends on the frequency and extrapolation step size when the medium velocity has been specified. Consequently the 2-D Rayleigh-Sommerfeld integral becomes the following convolution integral (cf., Figure 1-4): The spatial wavelet is defined by equation (2.25) and is dependent upon the parmeters Az and k. The convolution in (2.26) is very similar to the more familiar one-dimensional convolutional model of seismograms (Robinson, 1957). However, now the seismic source wavelet has been replaced by a spatial wavelet derived from wave theory. The spatial in the y direction) as perceived on a planar recording surface. Equation (2.25) repre- sents the wavefield value at a recording point (x, Az) due to a point source located at (x, Az) = (0,0). Thus the wavelet is essentially a weighting function which transforms Section 1.4.1). Also, Robinson's (1957) reflectivity series gives way to a series repre- senting the pressure source distribution along a line of the source plane, oriented in the x direction. Most geophysicists are comfortable with one-dimensional convolution, so that this approach forms an excellent viewpoint from which to examine the phenomena of wave propagation. (2.26) wavelet is the impulse response for a 2-D point diffraetor (or equivalently, a line source the pressure from the point source to a pressure distribution elsewhere in space (cf., 32 2.1.3 Seismic approximations and limitations of theory The application of scalar diffraction theory to reflection seismology requires that ap- proximations be made in the seismic reflection procedure (Trorey, 1970). First, a scalar theory demands that the seismic wavefields be either compressional or shear waves, but not both. This is equivalent to assuming that the seismic waves impinge with normal incidence on the reflection surface. For a compressional source normal incidence implies that there is no mode conversion, so that the generation of shear waves will be neg- ligible and the final wavefield consists of compressional waves only. In reality, normal incidence seismic sections may be approximated by the source-receiver geometry of a common midpoint stacked profile. The second approximation is that energy is treated as travelling only in one direction. That is, only the upcoming wavefield is considered, as a consequence of the exploding reflector model (cf., Section 1.2.2). Third, the reflec- tion coefficients of the reflectors are considered to be small so that transmission losses may be ignored. Fourth, the velocity in the medium between the point diffractors and receivers is assumed to be constant. Lastly, I mention that the scalar theory, as an approximate theory, has a shortcoming worth noting (Jackson, 1962, p.283; Goodman, 1968, p.32). The theory works best in the short wavelength or high frequency limit where the dimension of the diffraction opening is much larger than the wavelength, but deteriorates in its predictions otherwise. 2.2 Spatial aliasing I begin this discussion by briefly defining the phenomena of aliasing. When a signal is aliased it has not been sampled sufficiently often and a subsequent analysis of the signal will not be correct. Suppose a sinusoidal waveform, sin(kx), is sampled at a uniform 33 interval Ax. The wavenumbers k (where k = 2n/X and A is wavelength) within the interval 0 < A; < k^ are said to be in the principle interval. The highest wavenumber, fcjv, in the interval is referred to as the Nyquist wavenumber (Kanasewich, 1981, p.110). The relationship between k^ and A x is kN = (2.27) A x This expression exists because it is necessary to have at least two samples per wavelength to detect that harmonic component of a signal. A consequence of discrete sampling and the Nyquist criterion above, is that any k greater than kjy is impossible to distinguish from those in the principle interval. This means that if any signals are present with A less than Ajy (A,v corresponding to k^), their spectral amplitudes will be reflected back, or aliased, into the amplitude spectrum over the principle interval of wavenumbers (Kanasewich, 1981, p.111). The fact that the integral equation (2.24) must be diseretized for practical purposes demands that great care be taken to not alias the spatial wavelet when it is sampled in the x domain. Following the example of Berkhout and van Wulfften Palthe (1979) for the 3-D wavelet, some insight is gained into the 2-D case by studying the asymptotic behaviour of the wavelet in* the near and far fields. 2.2.1 Analysis by asymptotic approximations The near field behaviour of (2.25) can be appreciated by examining the case for krx < 1. This condition is satisfied if either k or the extrapolation step size, Az, is very small or 34 if both k and Az are small. In any case Az must be small. Then the Hankel function may be approximated by (Kinsler et al., 1982, p.449) Hp)(jbrx) = ^ _ ( * ! J ^ + t- 2 ( 2 2 g ) 1 v ; 2 16 ivkrx v 1 As a result the wavelet may be rewritten as where the subscript "nf" denotes near field. Equation (2.29) contains a real part inde- pendent of k and an imaginary component dependent on fc. In the near field, the ratio of the real part to the imaginary part is 4 /(7rfc 2 r 2 ) > 1, so that the dominant portion of the wavelet is given by the real component. As Az goes to zero, the real part becomes a Dirae delta function lim Wn{(x) ~^6{x). (2.30) Az—»0 Since the spectrum of a delta function is unity, if Az is very small there will be large spectral values near fc^ and beyond so that spatial sampling will result in aliasing of the wavelet. This necessitates the determination of a lower bound on Az to prevent aliasing in the x domain sampling. A minumum value on Az may be obtained by analysing the Fourier transform of the asymptotic wavelet (2.29). As shown, the real component of (2.29) becomes dominant and so one need examine its Fourier transform only. Using Gradshteyn and Ryzhik 35 (1965, p.312) the result is oo Az 7r(z2 + Az 2) >xdx — oo (2.31) That is, when the field point is close to the source, the spectrum is a decaying exponential with the decay rate governed by the extrapolation step size and it is independent of k. If Az is large, then the decay is rapid, whereas small Az values correspond to a near white spectrum. Clearly Az must be sufficiently large if the spectral values of the near field wavelet are to be substantially small at wavenumbers equal to the Nyquist wavenumber. Suppose that at kx — k^j the value of (2.31) is e which is a small fraction of the peak amplitude. Then acceptably small aliasing can be guaranteed if Using e = 10 - 2 or -40 dB gives A z m i n > 1.5Ax. This limiting rule may often be too large to model or migrate the geologic fine structure unless the sampling interval is reduced. However, a less stringent rule (i.e., picking e larger) would allow a smaller Az, but this < e. (2.32): Manipulation of (2.32) and (2.27) yields A z m i n > Ax (--Ine). (2.33) would result in an aliased wavelet whose ill effects might not be realized until a large 36 number of extrapolations had been performed. Therefore, in the near field, there exists a tradeoff between aliasing and the ability to accomodate detailed structure. Now I examine far field case where krx > 1 is the defining condition. The far field approximation of the Hankel function is given by (Kinsler et al., 1982, p.450) H > 2 ) { k r x ) = ^ ^ ' ^ ( 2 - 3 4 ) i and the asymptotic wavelet can be expressed as [Y -i(ky/x? + Az*-§ 7T) W f f ( x ) = - t W — • (2.35) V 27F ( x2 + / \ 2 2 ) t The subscript "ffi" in equation (2.35) denotes far field. As with the 3-D case (Berkhout and van Wulfften Palthe, 1979), fast variations in VKff(x) are determined by the expo- nential term. To ensure unaliased sampling of the exponential the Nyquist criterion is employed: k(rx,i+i ~ rx>i) < Tr. (2.36) Setting ( r I ) t + i — r I t ) = Arx and differentiating rx with respect to x gives \/x2 + Ax2 Substitution of (2.37) into (2.36) yields A \/x 2 + Az2 Ax < . (2.38) 2 oo 37 This inequality shows that the selection of Ax depends on the maximum x value used in defining the wavelet. The maximum x value is called the aperture of the wavelet, L. Since the wavelet is symmetric about x=0, I use L to indicate the half-length of the wavelet. Therefore, L can replace x on the right-hand side of (2.38). If L extends to infinity then Ax becomes independent of L and A Ax < - - 2 However, if L is shortened, then Ax can be larger than that given by (2.39). Equiva- lently, truncating the wavelet at a shorter L and keeping Ax constant allows A to be smaller while avoiding aliasing. Although truncation of the wavelet may mean that aliasing effects are reduced, truncation will lead to errors if not done carefully. Expres- sion (2.39) is a restatement of the Nyquist criterion for a harmonic wave with wavelength A. Therefore, the same rule applies to the spatial wavelet with the corresponding fre- quency / = v/X. Since there are many A values, one can select Ax corresponding to the minimum A. This avoids different aperture considerations for different frequencies. The use of inequalities (2.33), (2.38) and (2.39) is illustrated with examples of zero- offset synthetic seismograms due to a point source embedded in a constant velocity medium. The point diffractor is placed at a depth of 500 m at the centre trace position of a seismic line 1600 m long where the spatial sampling (geophone interval) is a uniform interval of 25 m (so the number of traces is 65). The velocity of the medium is 4000 m/s. According to these parameters and the Nyquist criterion, waves with frequencies above 80 Hz will be aliased by the array of geophones. Since the exploding reflector model is employed to generate two-way traveltime seismograms, the frequencies must be restricted to half that mentioned above (i.e., 40 Hz) if aliasing of the spatial wavelet 38 (2.39) is to be prevented. This factor of two arises from the half-velocity substitution of the exploding reflector hypothesis into equation (2.25). Clearly, this is a drawback if the field seismogram to be modelled contains data with frequencies greater than this threshold. Also, the near field asymptotic analysis implies that the extrapolation step size is limited by Az > 1.5Ax = 37.5 m, if spatial aliasing of the wavelet is to be circumvented. This limitation implies that vertical resolution can be no better 37.5 m, and consequently, there may be difficulty in modelling fine earth structure. The lowest frequency used was 10 Hz. Figure 2-3 shows the correct seismograms with the frequency and extrapolation step size limits are satisfied. In Figure 2-3a the seismogram was computed using Az = 500 m so that the diffraction response observed was propagated upwards in a single step. This example represents the exact answer. As can be seen the first arrival is at 0.25 s and the last arrival is at 0.47 s which are the expected two-way traveltimes. It should be noted that the temporal wavelets on each trace are not quite zero-phase as anticipated. In fact, these wavelets have been phase shifted by — n/4 radians due to the formulation of wave propagation in 2-D (c.f., equation (2.35)). This can be removed for interpretive reasons, but I do not do so, as it does not interfere with the observations. The seismogram in Figure 2-3b was computed using Az = 50 m and thus required ten iterative steps for its generation. The seismograms are nearly identical with those in Figure 2-3, except for the extremely weak diffraction tails which appear to be reflected back from the boundary formed by the last traces. This artifact may be viewed more clearly if viewed across the figure at an acute angle. The fact that there are a finite number of traces and a correspondingly finite number of samples in what should be an infinite spatial wavelet, leads to the artifact. For the point diffractor, only the first step in the recursion is exact and all subsequent steps are slightly inaccurate because 39 Figure 2-3. Correct seismograms. (a) Az=500 m and (b) Az=50 m. Parameters used were Ax=25 m, u=4000 m/s and frequency content of 10-40 Hz. Horizontal axis is surface offset coordinate and vertical axis is two-way traveltime in seconds. 40 the wavefield at the edge traces is not properly compensated for when the traces end abruptly. Therefore, the wavefield near the edges acts as a set of diffractors whose diffraction tails are not correctly cancelled by neighbouring diffractors. This artifact can be removed by extending the number of traces by a factor of two (padding each side of the original trace set by an equal number of traces), generating the seismogram and displaying only the centre traces. However, this adds to the expense of computation. The artifact is not a hinderence to the present observations and I refrain from removing it. The effect of extending the frequency range over the 40 Hz limit to 60 Hz is illus- trated in Figure 2-4. These seismograms were computed in one step using Az = 500 m. Although the antialiasing criterion has not been met, a very good seismogram is generated. The result compares favourably with that in Figure 2-3a, except the addi- tional frequencies have added to the sharpness and resolution of the temporal wavelet. The reason for this good result is that a single step extrapolation does not significantly distort the aliased frequency components of the seismogram. In contrast, Figure 2-4b shows the seismogram when Az = 50 m was used for its generation. Note that the char- acteristic hyperbola has been lost and in its place there appears a large checkerboard artifact between 0.30 s and 0.50 s with tails radiating out with increasing traveltime. The cause of this phenomena is the aliasing of the spatial wavelet combined with sev- eral extrapolation steps. The result is an amplification of aliased contributions which results in a very poor seismogram. However, equation (2.38) implies that a decrease in the spatial wavelet length or aperture can help acommodate the frequencies beyond the upper limit. The variation of L with frequencies beyond 40 Hz is displayed in Figure 2-4c and in Figure 2-4d I show the seismogram generated using this scheme. There is a clear improvement over the seismogram in Figure 2-4b, as one is now able to recognize a 41 Figure 2-4. Seismograms with upper frequency limit extended beyond 40 Hz to 60 Hz for (a) Az=500 m and (b) Az=50 m. Other parameters used were Ax=25 m and u=4000 m/s. 42 APERTURE VARIATION cn ZD J— Ql Ul Q_ < CO Ul < in 40 44 48 52 56 FREQUENCY (Hz) 60 0 . 0 1 0.1 1 0 . 2 0 . 3 0 . 4 1 0 . 5 0 . 6 0 . 7 I 0 . 8 0 . 9 ± 1.0 Figure 2-4 cont'd. Attempted reconstruction of (a) from (b) using aperture limitation of wavelets in frequency band 40-60 Hz. (c) Aperture variation as a function of frequency when the frequency corresponding to the Nyquist wavenumber is exceeded. The dashed line shows the number of samples needed in the aperture according to (2.38) and the solid line indicates the actual integer used for computation, (d) Only partial recovery is possible. 43 weak diffraction hyperbola. Nonetheless, the result is still poor in comparison with that in Figure 2-4a. This illustrates that equation (2.38) can be used to solve the aliasing problem partially, but it is still inadequate. The reason for this failure is that in limiting the aperture, the wavelet was truncated abruptly and significant information contained in those portions of the wavelets was rejected. The seismograms in Figure 2-5 demonstrate the effects of violating (2.33) and select- ing Az too small. Figure 2-5a shows the surprising result of using Az = 25 m. Although, Az is significantly less than the minimum value of 37.5 m (required by equation (2.33) if e = 10 - 2), the constructed seismogram is remarkably good. This is due to picking e very small and suggests that a larger e or smaller Az could be used. Picking Az = 10 m corresponds to e = 10 - 0 ' 5 and the result of this case is shown in Figure 2-5b. It is evident that this selection was a poor one and the seismogram shows the checkerboard pattern characteristic of aliasing. Therefore, selection of Az between 10 m and 25 m (closer to Az=10 m) will lead to considerable deterioration of the seismogram. The effects of extended frequency range and selecting Az too small have been ex- amined independently. Now I investigate the case where both these factors are acting simultaneously. To do so, I select the frequency range 10 to 60 Hz and Az = 25 m. Pre- viously Az — 25 m (Figure 2-5a) gave adequate results on its own while the frequency range of 10 to 60 Hz allowed some recovery of the diffraction hyperbola when the aper- ture was limited using Az = 50 m (Figure 2-4c). When the frequency range is extended and Az is made too small, the aperture limitation fails to recover any portion of the diffraction hyperbola (Figure 2-6). I conclude that when each aliasing factor, either extended frequency or selection of Az too small, acts independently, the antialiasing 4 4 Figure 2-5. Seismograms an extrapolation step size which is too small, (a) Selections of Az=25 m and (b) Az—10 m violate the antialiasing criterion, but the former generates an acceptable result. Other parameters used were Ax=25 m, u=4000 m/s and frequency content of 10-40 Hz. 45 Figure 2-6. S e i s m o g r a m s for w h i c h u p p e r f requency l im i t is 60 H z a n d A z = 2 5 m. O t h e r pa rame te rs used were A x = 2 5 m a n d v=4000 m / s . 46 criterion can be relaxed slightly. However, when these aliasing conditions are violated together there can be no hope of generating an adequate seismogram. 2.2.2 Analysis by Fourier transform The investigations of near and far field approximations to the spatial wavelets have provided rules for preventing spatial aliasing without really explaining what causes the problem. Considerable insight may be gained by examining the spectrum of the spatial wavelet in the wavenumber domain. To do so, a Fourier transformation from the space variable x to the correponding wavenumber domain variable, kx, is required. Fourier transformation of the spatial wavelet (2.25) is easily accomplished by W(kx) = - (2.40): The integral is given by Erdelyi (1954, p.56) as The Hankel function of one-half order may be expressed as an exponential (Magnus and Oberhettinger, 1940, p.l8) H[2)(Azv/k2~^kl) = i\P^(k2 - k2x)L* e - A V ^ I . (2.42) 2 V TTAZ 47 Substitution of equations (2.41) and (2.42) into (2.40) immediately yields the analytic expression for the spectrum w{kx) = E - I W F C 2 - F E ? . (2.43) Equation (2.43) is the phase shift operator used by Gazdag (1981) in Fourier domain modelling. The phase term contains the dispersion relation of the wave equation and may be expressed as Thus the spatial wavelet and the phase shift operator form a Fourier transform pair. A similar result was discussed in detail for the 3-D case by Schneider (1978) and Berkhout and van Wulfften Palthe (1979), but they did not explicitly present the 2-D case as I have done. The spectrum of the spatial wavelet, as given by (2.43) indicates that the wavelet is an operator having infinite bandwidth, extending to all wavenumbers. The portion of the spectrum for kx < k is made up of complex sinusoids and represents the wavefield propagating in the vertical direction. In contrast, the region kx > k is indicative of evanescent waves which decay exponentially in the vertical direction and which diffuse in the lateral coordinate. In Figure 2-7, the Fourier transform of the spatial wavelet in the kx domain is shown for various extrapolation step sizes. For this particular case the frequency used was 30 Hz, the velocity was 4000 m/s and the spatial sampling interval was 25 m. The corresponding propagation and Nyquist wavenumbers are 0.094 and 0.125 radians/m, respectively. The real componenent of the spectral wavelet is shown in Figure 2-7a. kz = yjk* - kj. (2.44) 48 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250 I 1 1 1 1 1H 1 , 1 1 a 1 I 1 1 1 1 1 —F 0.000 0.025 0.050 0.075 0.100 0.125 D.150 0.175 0:200 0.225. 0.250^ I 1 1 1 1 1 1 1 1 1 b 1 1 1 1 1 1 1 1- Figure 2-7. (a) Real and (b) imaginary components of wavenumber spectrum of spatial wavelet at various Az values. Parameters used were f=30 Hz and u=4000 m/s. Az values used were 5, 10, 20, 50, 100, 250 and 500 m arranged from top to bottom on each panel. Horizontal axis is wavenumber, kx in radians/m. Positive values are shaded. Computations were performed using equation (2.43). 4 9 Observe that for small Az that the amplitude is near one as expected for a cosine function of small argument. As Az increases, the propagating wave portion exhibits the oscillations of the cosine function. The abrupt change in character occurs at the propagation wavenumber, where there is the transition from propagation to evanescence. At small Az values, a large portion of the spectrum lies in the evanescent zone, but as Az increases the contributions this zone become smaller exponentially. The imaginary component (Figure 2-7b) illustrates the sine function variations in the wavelet. Again there is an abrupt change at the transition from propagation to evanescence. However, here the imaginary evanescent zone is zero. This example shows that for all Az the imaginary component will not be aliased because it has no contri- butions above k^. In contrast, this will only be true for the real component at large Az values where the evanescence is nearly naught at the propagation wavenumber. If Az is made small, then there are large contributions due to evanescence beyond k^. Therefore, aliasing arises due to evanescent waves when Az is smalll Now I consider the characteristics of the spectral wavelet if the operating frequency is raised to 50 Hz which corresponds to a propagation wavenumber of 0.157 radians/m (Figure 2-8). The results shown are similar to those in Figure 2-7, but the transition point from propagation to evanescence has moved beyond k^. In this case, both real and imaginary components have contributions above the Nyquist wavenumber. There- fore the wavelet will always be aliased if the operating frequency selected exceeds that corresponding to kj<j, regardless of the Az size. A solution to the problem posed in the discussion of Figure 2-7 (i.e., picking Az very small) was claimed by Berkhout (1981, 1982) for the 2-D problem. In those works, he suggested using a Taylor series expansion of the wavelet which was truncated after a few 50 Figure 2-8. Wavenumber spectra as in Figure 2-7, but with f=50 Hz. (a) Real and (b) imaginary components. 51 terms. Also, he stated that this approach removed the evanescent waves, thus allowing the computation of the wavelet in the space domain. However, I found Berkhout's (1981, 1982) formulation rather confusing. A solution to the second problem (cf., Figure 2-8) of propagation frequency exceeding the Nyquist was demonstrated by the use of equation (2.38) (aperture limitation) in Figure 2-4c, but the results were far from adequate. To circumvent both of these aliasing problems, I propose using the phase shift op- erator or spectral wavelet (2.43). This approach is conceptually simple when compared with Berkhout's (1981, 1982) contributions. Since an analytic expression for the spa- tial wavelet in the wavenumber domain exists, it is easily computed in that domain. Contributions beyond the Nyquist wavenumber may be removed by the application of a lowpass filter (0-0.125 radians/m). This ensures that when the wavelet is Fourier trans- formed back to the space domain, it will' be unaliased for all Az and all frequencies. The selection of a lowpass filter is straightforward for this problem. The preservation of the amplitude and phase information of the wavelet as much as possible is required and this may be achieved in the spectral domain multiplying a zero-phase box-car win- dow (Oppenhiem and Schafer, 1975, p.239) with the phase shift operator. However, an abrupt truncation of the spectrum will cause the generated spatial wavelet to ring. The ringing, known as Gibbs phenomenon, may be minimized by selecting a lowpass filter which does not have a sharp discontinuity like the box-car window. A better choice is the Butterworth window which is maximally flat with unit amplitude in the passband and maximally flat with zero amplitude in the stopband. The transition between pass- band and stopband is determined by the order of the window. High order will provide a sharper transition than low order, but high order will also induce a more pronounced Gibbs phenomenon. In the examples that follow a Butterworth window of order four was used. The cutoff wavenumber used for the window was 0.8A; when k < and 0.8fcjv for 52 k > kflj- Windows implemented with these parameters will smooth the spectrum near the Nyquist wavenumber and reduce the amount of energy in the evanescent wavefield. There are other choices of windows that could work equally well, but I do not explore their use in this study. The selection of the window is important, as is the strategy with which to implement the window. However, in order to preserve the continuity of this section I do not elaborate on these topics here. Instead, a vital and detailed discussion of Butterworth windows and windowing strategy is given in the Appendix. The success of the new method of wavelet computation is demonstrated in Figures 2-9, 2-10 and 2-11 where the parameters used were identical to those in Figures 2-4b, 2-5b and 2-6. Those latter figures demonstrated the limitations of computing the spatial wavelet in the space domain while attempting simultaneously to eliminate aliasing. A comparison of Figures 2-4b and 2-9 shows that the bandlimited spatial wavelet recovers the hyperbola (which was lost in Figure 2-4b due to an extended frequency range) in a manner superior to that of the aperture limiting scheme of equation (2.38). Thus, the Butterworth bandlimiting approach safely allows use of frequencies beyond those corresponding to the Nyquist wavenumber. However, the hyperbola tails decay faster in Figure 2-9 than those observed using the space domain method (Figure 2- 4c). This happens because the removal of the higher wavenumbers in the bandlimiting procedure. Since the Hankel function high wavenumber (short wavelength) components (c.f. equation 2.25) occur at large offsets, the Butterworth smoothing leaves only very small amounts of low wavenumber energy at those positions. This also results in the removal of the boundary reflections that were observed earlier and implies that traces need not be padded as initially suggested. The result of a violation of step size was shown in Figure 2-5b. This situation is remedied by the bandlimited wavelet as seen in Figure 53 0 . 0 O J 0 . 2 1 0 . 3 I 0 . 4 1 0 . 5 0 . 6 0 1 7 0 . 8 0 . 9 T.O Figure 2-9. Seismogram generated with spectral bahdlimited wavelets. Upper fre- quency limit is 60 Hz and Ar=50 m. Other parameters used were Ax=25 m and v=4000 m/s (cf., Figure 2-4b). 5 4 0 . 0 I 0.1 j l . 0 . 2 jl 0 . 3 Jl 0 . 4 j_ 0 . 5 j _ 0 . 6 JL 0 . 7 JL 0 . 8 I 0 . 9 1 1.0 Figure 2 - J L 0 _ S e i s m o g r a m genera ted w i t h s p e c t r a l b a n d l i m i t e d wavelets us ing A z = 1 0 m . O t h e r pa rame te r s used were A i = 2 5 m a n d t>=4000 m / s (c.f., F i g u r e 2 -5b) . 55 0 . 0 0 .1 t 0 . 2 0 . 3 0 . 4 0 .5 1 0 . 6 JL 0 . 7 j_ 0 . 8 I 0 . 9 1.0 \\\ Figure 2-11. Seismogram generated with spectral bandlimited wavelets. Upper fre- quency limit is 60 Hz and Az=25 m. Other parameters used were Ai=25 m and u=4000 m/s (cf., Figure 2-6). 56 2-10. This result is most encouraging as it indicates that fine geologic structure may be modelled using the spatial wavelet method. Another failing of the spatial sampling scheme was considered in the discussion of Figure 2-6 where both frequency limit and Az antialiasing conditions were not satisfied. As expected, the result was a very poor seismogram. In contrast, the spectral method of computation provides an excellent seismogram with higher frequency content at the desired step size increments (Figure 2-11). These results indicate that bandlimiting the spatial wavelet in the wavenumber domain using a Butterworth window provides a new and improved method of avoiding aliasing. However, as a consequence of bandlimiting with this scheme this is a limitation on the steep dips that can be modelled (Section 2.3 pursues this matter in futher detail). 2.2.3 Frequency spectra The spectral bandlimiting procedure has overcome the difficulties encountered in the space domain computation of synthetic seismograms when the parameters are extended beyond theoretical limits. Although it is now possible to generate good seismograms under conditions where previously impossible, one would still like to investigate how well the seismograms preserve the parameters of the spatial wavelets. In particular, the frequency content of the seismogram is an important element in understanding some effects that may be due to spectral bandlimiting. Rather than calculating the frequency spectra along all traces of a seismic section, I select a few representative traces. Since the seismograms are symmetric about the 33rd trace I use this centre trace along with the 49th and 65th traces which give information off the axis of symmetry. A brief note on aesthetics. The seismograms shown so far were frequency filtered by a cosine window (also called a Hanning window). The cosine tapering was implemented 57 in the manner of a trapezoidal window, but with the cosine determining the attenuation of frequency components rather than a simple ramp. This resulted in the production of clean (as opposed to noisy) seismograms. In the spectral plots that follow the cosine window corner frequencies are given in the figure captions. In Figures 2-12a and 2-12b I show the spectra for the seismograms of Figures 2-3a and 2-3b respectively. Since the seismogram in Figure 2-3a was generated nonrecursively it represents the exact solution and its spectra represent the exact spectra. Amplitudes decrease from inner trace spectra to outer trace spectra due to a l /r decay and the Butterworth windowing. Unlike the exact spectra, Figure 2-12b exhibits a perturbing ripple phenomena on an otherwise exact spectra. This effect is from reflections at the boundaries of the data due to recursive computation (c.f. Section 2.2.1). Since the reflections bounce back onto the inner traces, they modify the traces (outer traces most drastically) and their corresponding spectra. Proof of this explanation can be given with a brief analysis. On trace 33 of Figure 2-3b the initial pulse arrives at 0.250 s and the reflected pulse arrives at 0.838 s. The spacing between signal is 0.588 s and so the spacing between the peaks of the ripples in the spectrum will be the reciprocal of this time separation, or 1.70 Hz. Close inspection of of the spectrum in Figure 2-12b shows that there are about five or six peaks for every 10 Hz which corresponds to the expected peak separation. For the 49th trace the first and boundary reflected arrival times are 0.320 s and 0.650 s respectively. This implies a spectral peak separation of 3.03 Hz. Again, this expectation is borne out by Figure 2-12b where observations yield approximately four peaks per 10 Hz. As mentioned earlier, the boundary reflections may be attenuated by padding traces. This will also improve the spectra, but will also raise the expense of computation. 58 Figure 2-12. Frequency spectra (a) for Figure 2-3a and (b) for Figure 2-3b respectively. Parameters used were Az=500 m (for 2-12a) and Az=50 m (for 2-12b). Also Ax=25 m, v=4000 m/s and the frequency content was 10-40 Hz. Spectra of 33rd, 49th and 65th traces are shown from top to bottom on each panel. Horizontal axis is frequency in Hz. Cosjne window corner frequencies used were 15 and 35 Hz. 59 Next I consider the spectra of seismograms which violate the upper frequency limit of 40 Hz. Figure 2-4a is such a seismogram, computed by the space domain method in a single step for up to 60 Hz frequency content. Even though the frequency limit was violated a good seismogram was observed. Not surprisingly, a correspondingly good set of spectra are shown in Figure 2-13a. I previously discussed this issue and concluded that the aliasing errors were not significant for nonrecursive implementation. Figure 2-13b shows spectra for the nonrecursive implementation of the spectral bandlimiting method. In contrast to the spectra in Figure 2-13a, these spectra show a slight deterioration after 40 Hz. There are two possible reasons which may explain the apparently poor performance of the spectral method. First, there is the Gibbs phenomena due to the shape of the Butterworth window used in the bandlimiting process. Although for a fourth order Butterworth window the Gibbs phenomena is small, it is sufficient to cause significant distortion to wavelet samples far from the zero lag sample. The second perturbation could arise from bandlimiting itself. The true wavelet contains frequencies higher than 40 Hz, but at higher frequencies greater portions of significance lie beyond the Nyquist wavenumber and are being removed. In the place of correct information for the band 40 Hz to 60 Hz, one obtains a distorted representation because significant information has been rejected. Despite this fact, the procedure stabilized computation of the seismograms in a recursive manner (Figure 2-9) where it was previously not possible (Figures 2-4b and 2-4c). Finally, I examine the spectra for the seismograms displayed in Figures 2-9 through 2-11. Figure 2-14 shows the spectra for Figure 2-9 which is well behaved and exhibits a very rapid amplitude decay for the outer traces. In fact, there appears to be almost no energy for the 65th trace spectrum. The same observation applies to Figures 2-15 and 2-16 which are the spectra for the seismograms in Figures 2-10 and 2-11 respectively. 60 0 10 20 30 40 50 60 70 80 90 100 110 120 1 1 1 1 1 1 1 1 1 1 1 r Figure 2-13. Frequency spectra (a) for Figure 2-4a by space domain computation and (b) for spectral bandlimited counterpart. Parameters used were Az=500 m, Ax=25 m, t>=4000 m/s and frequency content of 10-60 Hz. Horizontal axis is frequency in Hz. Cosine window corners frequencies used were 15 and 55 Hz. 61 Figure 2-14. Frequency spectra for seismogram in Figure 2-9. Parameters used were Az=50 m, A x = 2 5 m, t ;=4000 m/s and frequency content of 10-60 Hz. Horizontal axis is frequency in Hz. Cosine window corner frequencies used were 15 and 55 Hz. 62 '- 1 \ 1 K - h h 1 I 1. h 1 h — r— Figure 2-15. Frequency spectra for seismogram in Figure 2-10. Parameters used were Az="10 ni, Ax=25 m, u=4000 m/s and frequency content of 10-40 Hz. Horizontal axis is frequency in Hz. Cosine window corner frequencies used were 15 and 35 Hz. 63 Figure 2-16. Frequency spectra for seismogram in Figure 2-11. Parameters used were Az=25 m, Ax=25 m, v=4000 m/s and frequency content 10-40 Hz. Horizontal axis is frequency in Hz. Cosine window corner frequencies used were 15 and 55 Hz. 64 Regardless of slightly disturbed spectra, spectral bandlimiting provides a method of extending the limits of operation beyond those defined by space domain computation with greater stability for recursive implementation. However, the penalty paid for this increase in flexibility is the sacrifice of reliability in the regions of extension. 2.3 Wavelet truncation The problem of wavelet truncation was briefly alluded to, when I discussed aperture lim- itation as a means of attenuating distortions due to spatial aliasing. However, truncation effects appear in two other important considerations for the computation of synthetic seismograms. First, the limited number of traces restricts the number of samples that can be used to represent the spatial wavelet. Thus, although the wavelet is infinite in extent, it must be truncated for practical computation. All the wavelets in the seis- mograms computed thus far were limited by the number of traces (i.e., the number of samples used for the aperture was 65). However, it was shown that when seismograms were computed using the spectral bandlimiting scheme the wavelet amplitudes at far offsets were very small. This suggests that if wavelets are computed with this method, then an approach to wavelet truncation may be found so that an adequate represen- tation of the wavelets exists in only a few discrete samples. The second consideration arises when attempting to accomodate lateral velocity variations in the geology. Here, zones of lateral variation may be small and if any modelling is to be done, then the spatial wavelet must also be small. Otherwise, use of a larger wavelet will lead to the smearing out of the lateral heterogeneity and a poor representation of the structure will result. In addition to the possibility of accomodating lateral velocity, another advantage gained by wavelet truncation is the reduction in computation cost as the convolution 65 procedure will be less time consuming using a shorter wavelet. Thus, it is imperative to determine a scheme by which the wavelet aperture can be limited to a minimum, yet still contain a maximal amount of information for wavefield computation. To the best of my knowledge neither Berkhout (1981, 1982, 1984) nor Berkhout and van Wulfften Palthe (1979, 1980) discussed the effects of wavelet truncation in their investigations. Therefore, the following analysis, which quantifies the effects of wavelet truncation and which presents an approach to remedy this problem, is a new contribution. 2.3.1 Wavelet energy level In this section, I develop an approach to limiting the wavelet aperture based on the energy content of the computed bandlimited wavelet. The general problem to truncate an infinite length signal after N discrete samples such that the loss of information due to missing samples is negligable. I argue that if an infinite signal contains a large portion of its total energy within the first TV discrete samples, then the representation of the signal by those samples is sufficient. Therefore, truncation of the signal at i V samples is justified and the error induced by doing so is expected to be minimal. In the space domain, the bandlimited wavelet is an infinite length signal and com- putation of its total energy in this domain is not always possible. However, due to the fact that the signal is bandlimited an approach exists by which the energy may be determined in the wavenumber domain, where an analytic expression for the wavelet is known. Parseval's theorem (Kaplan, 1981, p.244) states that (2.45) 66 where f(x) is the signal and F(kx) is the corresponding Fourier transform. The factor of 1 /2TT arises from the definition of the asymmetric Fourier transform which normalizes the inverse transform only. The right hand side of equation (2.45) for the wavelet problem can be computed easily by replacing the infinite limits of integration by limits that are the lower and upper bounds on the wavenumbers, 0 and kN respectively. Thus, the total energy of a bandlimited spatial wavelet becomes a known quantity. The fraction of the total energy in the first TV samples of a signal is the JVth partial energy and for a symmetric signal, like the spatial wavelet, it is defined by {*" \W(x)\2dx , . Ew = 2 T T ^ — 1 — — . (2.46 j0kN \W(kx)\2dkx W(kx) is the Butterworth smoothed phase shift operator and W(x) is the corresponding space domain representation. When E^ reaches some threshold value, say 0.95, then truncation of the signal1 at N samples will contain 95 percent of the total energy. Only 5 percent of the total1 energy is contained in the infinite number of samples outside the TVth sample and it is safe to say that their effect will be negligable over a single extrapolation step. T examine this more rigorously later (Section 2.3.3). 2.3.2 Minimum aperture wavelet The determination of the minimum length or aperture of a spatial wavelet based on energy considerations is highly dependent on several parameters. The extrapolation step size, the frequency and the velocity all play important roles in this quest. Some heuristic arguments, given below, provide insight into how these three parameters may affect the energy content in a spatial wavelet. 67 Earlier in this work I showed that the spatial wavelet approached a delta function as Az became small. The energy in delta-like functions tends to be concentrated in the first sample or first few samples of the discrete representation. Thus, a small number of samples may be adequate in their description of the function. This implies that for small Az the wavelet may be well represented by only a few samples. This is exactly what is needed when modelling complicated cases, that is, for media with laterally varying velocities and fine geologic structure. Therefore, selection of Az to be small will give the the smallest aperture wavelet that maximizes the energy content. The wavelet expression in equation (2.25) indicates that the Hankel function argu- ment is krx or urx/v, using k = OJ/V. It is easy to see that u> and v are contracting and stretching factors respectively, for the Hankel function. If w increases and v is kept constant, the Hankel function contracts. As a consequence, more energy is concentrated in the first few samples of a discrete representation. Similarly, if v decreases while u> is held fixed, the energy is again concentrated in the first few samples. Thus, the minimum aperture of the wavelets can be found when the frequency is at its maximum value and the velocity is at its minimum value. Accordingly, every wavelet for every Az-frequency-velocity combination will have a different minimum aperture. Therefore, it becomes necessary to select a standard minimum aperture from all of these. I use the maximum of the minimum apertures for this purpose. This aperture represents the worst case of all the given possibilites and, therefore, it satisfies the energy threshold criterion for all the frequencies and all the velocities once the minimum extrapolation step size has been specified. Theoretically, this aperture will be a property of the spatial wavelet for the lowest frequency value and highest velocity value of interest. I illustrate the relationships between Az, u and 68 >- o 1 .0 cc UJ ^ 0 . 8 £ 0 . 6 < Q 0 . 4 U l M ^ 0 . 2 or o NORMALIZED PARTIAL ENERGY CURVES Tf/ i i i 1 1 M i i i i i i I I 1 1 i i i i 1 1 1 1 1 0 1 1 0 2 LOG SAMPLE NUMBER 1 0 3 or UJ m < to o (/) UJ CC o o THRESHOLD SAMPLE VARIATION 1 0 0 2 0 0 3 0 0 AZ (m) 4 0 0 5 0 0 Figure 2-17. Wavelet energy variation with extrapolation step size. Fixed parameters are /=30 Hz, v=4000 m/s and Ax=25 m. (a) Normalized partial energy curves are for Az=5, 25, 50, 100, 250 and 500 m are arranged from top to bottom in this panel. The energy for small Az rises to the saturation level faster than for larger Az values, (b) Variation of wavelet aperture satisfying a threshold value of 99.99 percent for various Az values. 69 NORMALIZED PARTIAL ENERGY CURVES 6 8 10 SAMPLE NUMBER 12 14 THRESHOLD SAMPLE VARIATION TO 20 30 40 50 FREQUENCY (Hz) 60 Figure 2-18. Wavelet energy variation with frequency. Fixed parameters are Az=5 m, v=4000 m/s and Ax=25 m. (a) Normalized partial energy curves for /=60, 50, 40, 30, 20 and 10 Hz are arranged from top to bottom on this panel. The energy for high / rises to the saturation level faster than for low / values, (b) Variation of wavelet aperture satisfying a threshold value of 99.99 percent for various f values. 70 N O R M A L I Z E D P A R T I A L E N E R G Y C U R V E S 4 6 8 S A M P L E N U M B E R 10 T H R E S H O L D S A M P L E V A R I A T I O N 2000 4000 V E L O C I T Y ( m / s ) 6000 Figure 2-19. Wavelet energy variation with velocity. Fixed parameters are A2==5 m, /=30 Hz and Az=25 m. (a) Normalized partial energy curves for v=1000, 2000, 3000, 4000, 5000 and 6000 m/s are arranged from top to bottom on this panel. The energy for low v rises to the saturation level faster than for higher v values, (b) Variation of wavelet aperture satisfying a threshold value of 99.99 percent for various v values. 71 v with the minimum length in Figures 2-17, 2-18 and 2-19 respectively. Computations were performed using single precision and integration was done using a trapezoidal rule. These diagrams demonstrate that computations support the theoretical expections. All of these diagrams were computed for a threshold value of 99.99 percent and a Ax=25 m. Although, I do not show what happens for other choices of threshold value and Ax, I discuss these possibilities based on deductions from Figures 2-17 through 2-19. Selection of smaller threshold values give the result that smaller aperture sizes may be adequate. This leads to seismograms poorer than those computed with larger aperture wavelets. Since the 99.99 percent threshold gives suitably small values for the aperture size, I think that it is a good choice and smaller threshold values should be avoided. Smaller Ax values yield more accurate integration results and require that the wavelet aperture be larger than that for choices of larger Ax. Thus, the normalized partial energy curves reach their saturation limits at slightly slower rates. 2.3.3 Seismograms, speetra and stability With a quantitative theory for limiting the wavelet aperture laid out, I demonstrate its efficacy with seismograms for point diffractors as before. Again the point source is located at a depth of 500 m and it is embedded in a homogeneous medium with v=4000 m/s. The trace spacing is Ax=25 m. Figure 2-20 shows seismograms computed using the band and aperture limiting schemes that I have developed. The energy threshold used was 99.99 percent as in the experimental curves and the frequency range was from 10-80 Hz. An extrapolation step size of 50 m (Figure 2-20a) used a wavelet aperture of 42 samples to generate the 72 a 0.0 I 0.1 j - 0.2 j l 0.3 J l 0.4 J L 0.5 J L 0.6 JL 0.7 j l 0.8 JL 0.9 1 1.0 0.0 OJ I 0.2 0.3 0.4 1 0.5 0.6 0.7 t 0.8 t 0.9 1.0 Figure 2-20. Seismograms computed with band and aperture limited wavelets, (a) Az=50 m and (b) Az=25 m. The hyperbolas in both panels have the same size, the one in (b) decays in amplitude faster than the one in (a). Other parameters used were Ax=25 m,;v=4000 m/s and frequency content of 10-80 Hz. Horizontal axis is offset coordinate and vertical axis is two-way traveltime in seconds. 73 0 . 0 O J 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 1.0 0 . 0 O J 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 1.0 Figure 2-20 cont'd. Seismograms computed with band and aperture limited wavelets, (c) Az=10 m and (d) A z = 5 m. The amplitude in (d) decays faster than that in (c). These hyperbolas are of the same approximate size as those in (a) and (b). 74 hyperbola. The emergence angle of waves can be determined using 0 = tan-1\{Ln-l)Ax/Az}, (2.47) where Ln is the number of samples corresponding to an aperture of size L. Theoret- ically, using Az=50 m and Ln = 42 corresponds to the extrapolation being able to accomodate waves travelling at emergence angles of up to approximately 87 degrees. Similarly, for extrapolation step sizes of 25, 10 and 5 m the Ln was 20, 8 and 7 samples respectively. These all correspond to wavelets that describe rays with emergence angles up to 87 degrees. However, on observing the seismograms, I can only measure angles up to approximately 45 degrees on the tapered hyperbolas. It appears that the theoretical expectation is not fullfilled. This happens becouse of the combined effects of Butter- worth windowing (to avoid aliasing), wavelet truncation (to avoid large expeses) and recursive wavefield eontiouation (to allow modelling of fine geologic structure)1. Thus, an interpreter would not be able to discern structures dipping at angles greater than 45 degrees. Even this is optimistic, since most of the energy in the hyperbolas of Figure 2-20 is found near the apex. A more conservative estimate for the steepest dip that can be easily observed is 30 degrees. Another characteristic of these seismograms is the complete lack of hyperbola tail reflections from the boundary defined by the last traces. The reason for this apparent improvement is due to the weak diffraction tails from the aperture limiting procedure. The boundary reflections are there, but well below the recognizable level. Also, the location of the point source at a position far from both boundaries enhances this illusion. 75 As in previous sections on circumventing the aliasing problem, I examine the fre- quency spectra (Figure 2-21). Spectra are displayed for only two of the seismograms in Figure 2-20, namely 2-20a and 2-20c. I use Az=10 m as it is accepted as a lower limit of seismic resolution and Az=50 m just as a comparison. Both spectral diagrams show that the amplitudes decay with distance from the apex of the hyperbola. However, the spectra of the Az=10 m hyperbola (Figure 2-21b) are not as smooth as those of the other hyperbola. This happens because the small errors in the wavelet for the smaller extrapolation step size eventually become significant enough to distort the spectrum. I found that increasing the wavelet aperture reduced the effect of these perturbing ripples. Setting the aperture equal to the number of traces resulted in frequency spectra very much like that in Figure 2-21a. In spite of small spectral distortion the interpreter's eye cannot detect such subtleties on the seismograms and therefore it is not a cause for great alarm. In the case where the number of extrapolations becomes very large, however, the distortion of frequency spectra will inevitably destroy the seismogram. The computation of these seismograms with minimum aperture selection was four to five times faster than their counterparts using apertures as large as the number of traces. Given the savings in cost, I forsake accurate spectra in preference of speed, provided that the spectral distortion does not have a deprecating effect on the seismogram and its interpretational quality. In addition to investigating the frequency spectra I analyse the wavenumber spectra to examine the effects of recursive extrapolation on wavelet aperture limitation. Again, I use the seismograms in Figures 2-20a and 2-20c. These results are displayed in Figure 2-22. An immediate conclusion made from this figure is that the spectra from the Az=10 m extrapolation result are more distorted than those of A2=50 m. This is 76 0 1Q 20 30 40 50 60 7D 80 90 100 110 120 4 1 1 1 1 -+ 1 1 1 1 1 r a 1 (: 1 1 1 1 1 1 1 1 1 h 0. 10; 20' 30: 40s 50; 60 70' 80- 90 100. 110, 120 7 r H h H f 1 \r H 1 1 : Y 1- b 1 1 1 y 1 1 1 1 1 1 1 1—I Figure 2-21. Frequency spectra for seismograms computed using (a) Az=50 m and (b) Az=10 m. Spectra are for the 33rd, 49th and 65th traces of each seismograms and are arranged from top to bottom in that order for each panel. The horizontal axis is frequency in Hz. Cosine window corner frequencies used were 25 and 65 Hz. 77 FREQUENCY=10.7 Hz Q 3 Q_ < .0 0.04 0.08 0.12 WAVENUMBER ( rad/m) FREQUENCY=50.4 Hz Q ZD} CL. < 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0 - —\ f— ~ \ X V \ \ \ \ v - \ v G i i \ v , 0 0.04 0.08 0.12 WAVENUMBER ( rad/m) FREQUENCY=30.0 Hz < .0 0.04 0.08 0.12 WAVENUMBER (rad/m) FREQUENCY=70.7 Hz D_ < 1.2 1.0 0.8 0.6 0.4 0.2 0.0 0 -A — 1— \ \ r \ x \ v \ v -• \ v \ v d \ v 0 0.04 0.08 0.12 WAVENUMBER (rad/m) Figure 2-22. Wavenumber spectra for seismograms computed using Az=50 (dashed line) and 10 m (solid line) at various frequencies, (a) /=10.7 Hz, (b) /=30.0 Hz, (c) /=50.4 Hz and (d) /=70.7 Hz. 78 again due to the errors in wavelet truncation becoming significant after many recursive extrapolations when Az is small. In partcular, the result for the /=10.7 Hz spectra appears worse than the others. This is caused by the wavelet aperture being too small, even though it contained a large portion of the energy when Az=10 m. I think this is related to the frequency distortion problem described in the discussion of Figure 2-21b. The effect can be reduced by increasing the wavelet aperture in a judicious manner. Another interesting observation is the reduced band range of the spectra for the smaller extrapolation step size. This is expected from the Butterworth smoothing win- dow applied in the bandlimiting procedure to attenuate aliasing. Each extrapolation will reduce the spectral bandwidth and since small Az will require more steps to propa- gate a wave the same distance as larger Az, the wavenumber spectra will suffer greater distortion as a consequence. The inevitable distortion in the frequency and wavenumber spectra of the spatial wavelet lead to instabilty in the extrapolation procedure if it is done too often. The transition from stabilty to instabilty is difficult to analyse in advance of seismogram computation. In computing seismograms for this study I have not encountered insta- bility and I do not pursue it here. 79 3. M O D E L L I N G S E I S M I C R E F L E C T I O N D A T A 3.1 Homogeneous acoustic media A homogeneous acoustic medium is one in which only compressional travelling waves are considered and their velocity of propagation is constant throughout the space of interest. Since seismic reflections arise from variations in acoustic impedance, this implies that only the densities of the media involved are changing. Although this is not a very realistic model, it is adequate as a tool for demonstrating the modelling procedure. In the synthetic seismograms that follow I show that the spatial convolution ap- proach using the spatial wavelets may be successfully applied with a view to modelling seismic reflection data in homogeneous acoustic media. All seismograms simulate a field experiment conducted over the frequency range 10-80 Hz for u=4000 m/s, Ax=25 m (over 1600 m) and Az=10 m. Also. I assign a value of one to all reflection coefficients for the purpose of illustration. 3.1.1 A truncated model The first example is that of three truncated, horizontal reflectors, each of which is at a different depth of burial. The earth model is shown in Figure 3-la and the seismic section is displayed in Figure 3-lb. This computation used a model defined by 65 points sampled at a uniform separation of 25 m. The sampling interval for the spatial wavelets was also 25 m. Therefore, in addition to a small Az choice, the frequencies above 40 Hz would cause aliasing of the wavelets if computations were done in the space domain. However, with the spectral method of wavelet computation these hazards are removed. 80 EARTH MODEL x i— Q_ U J Q 400 800 1200 LINE COORDINATE (m) 1600 0 . 0 0.1 0 . 2 1 0 . 3 0 . 4 1 0 . 5 0 . 6 0 . 7 X 0 . 8 0 . 9 1.0 Figure 3-1. A truncated model, (a) Earth model with constant velocity of 4000 m/s. (b) Synthetic reflection seismogram for (a). Vertical axis is two-way traveltime in seconds and horizontal axis is offset coordinate with Ax=25 m. Arrow indicates termination of reflectors. 81 11 m m m m ™ 0 . 0 0.1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 . 1.0 Figure 3-1 cont'd, (c) Expanded view of diffractions from the left half of (b). Arrows from left to right indicate positions at the end of diffractions from the bottom to top reflectors respectively, (d) Stolt migration of modelled data shows that even the poor diffractions are collapsed well. 82 The seismogram correctly reproduces the two-way traveltimes of the flat reflectors and of the diffractions due to the truncation. It is known that the diffraction branch under the reflection is the inverse of the branch that extends beyond the reflector. However, the seismogram does not produce the diffraction branch beneath the reflector. This arises because the diffractions generated by the termination of the reflectors are small and the low amplitude of negative polarity event is hidden by the large amplitude primary reflection. Despite this shortcoming, the spatial wavelet modelling algorithm appears to correctly show that the curvature of the diffractions becomes smaller with depth. The diffractions at greater depth are of larger lateral extent than those at shallower depths (Figure 3-lc). This happens because the diffractions from deeper depths have more extrapolation steps over which to evolve than do diffractions closer to the recording surface. Spatial wavelets are computed to have small apertures. When the wavelet is convolved with a point source a wavefield represented by a diffraction hyperbola (on a seismic section) is generated. Each successive convolution expands the wavefield and the diffraction response. The extent of wavefield and diffraction evolution is limited by the number of extrapolation steps between the source and the recording surface. Therefore, if a source is located at a shallow depth the diffraction development will be stopped before a more complete response can evolve. Consequently, deeper sources generate larger diffraction hyperbolas. The overall result, although not ideal, is certainly of good interpretational quality. An important test of the spatial wavelet algorithm is to migrate the modelled data and recover the original structural geometry. Migration was accomplished using the Stolt (1978) method in the frequency-wavenumber domain and the result is shown in Figure 3-ld. Even though the diffractions are not of full extent in the modelled data, they are sufficient, to produce a very good migrated section. 83 3.1.2 Dipping models Next, I show seismograms for a suite of reflectors dipping at angles of 15, 30, 45 and 60 degrees (Figures 3-2, 3-3, 3-4 and 3-5 respectively). These examples demonstrate the dip limitations of the spatial wavelet algorithm. These computations represented the models by 129 points which corresponded to a spatial sampling interval of 12.5 m. I doubled the number of points so that the density of point reflectors was sufficient to simulate a continuous dipping reflector. The seismograms show every other trace so that the field parameters were mimicked. The spatial wavelets were sampled at the same interval. Therefore, frequencies up to 80 Hz were accomodated without danger of aliasing. So in these examples the use of the spectral method of wavelet computation served only to eliminate aliasing due to the small choice of Az. This note applies to other parts of Section 3.1 with regards to computations. In Figure 3-2 the 15 degree dipping reflector model and seismogram are shown. The modelled data correctly show that at the sharp corner at the top of the ramp the amplitudes are reduced and that the sharp character of the structure becomes smoothed. This is expected due to diffraction at a sharp corner where only scattered waves are registered on the seismogram. In contrast, the corner at the bottom of the ramp in the model produces a high amplitude response on the synthetic seismogram. This happens because the signal is composed of reflections from the dipping and lower flat planes plus diffraction at the corner. A structural dip of a degrees in the geology will be transformed to an apparent dip (/?) in the seismic section according to P = tan~l(sin a). (3-1) 84 EARTH MODEL 400 800 1200 LINE COORDINATE (m) 1600 0.0 0.1 J L 0.2 JL o.3 0.4 X 0.5 0.6 0.7 0.8 i 0.9 1.0 Figure 3-2. Model with 15 degree dip. (a) Earth model with constant velocity of 4000 m/s. (b) Synthetic seismogram for (a). Vertical axis is two-way traveltime in seconds and horizontal axis is offset coordinate with Ax=25 m. Arrows indicate positions of corner points in the model. 85 Consequently, a true dip of 15 degrees becomes an apparent dip of 14.5 degrees in theory. The observed dip on the seismogram is 14.5 degrees by measurement. Therefore, the spatial wavelet method correctly maps a true dip to the apparent dip. A 30 degree dipping structural model and its corresponding seismic section are dis- played in Figure 3-3a and 3-3b respectively. Like the 15 degree model the upper corner has reduced amplitudes and the lower corner has reinforced amplitudes. The theoretical and observed dips compare favourably, being 26.6 and 26.8 degrees respectively. Unlike the shallower dip example, this example shows that there is another reflection response from below the lower flat reflector. This is event is from the dipping segment as can be confirmed by the continuous character across the fiat reflector. In addition to the expected features there appears a linear artifact just below the lower flat reflector at a time of 0.62 seconds. An explanation for this feature is discussed throughly in the Appendix. Similar artifacts persist in further examples, but fortunately, these are small effects and leave the seismograms essentially undisturbed. Migration was accomplished' using the Stolt method in the frequency-wavenumber domain. Since the migration algorithm produces a corrected vertical traveltime-image model (the time-image model for Figure 3-3a is shown in Figure 3-3c), it is computed by simply dividing the distances by the half-velocity. The recovered model (Figure 3-3d) gives a basis for judging the quality in the modelled data. The migrated result shows that the time-image model is reconstructed very well. However, the amplitudes near the downdip corner are rather weak. This may be from the poor representation of the diffracted energy at that corner. Also, the presence of the modelling artifact is noted in the migrated section. 86 EARTH MODEL 0 200 - I T 400 X 1 600 1 Q_ ^^^^ LU 800 ^^^^ Q 1000 - a 1200 i i i i i i i 0 400 800 1200 LINE COORDINATE (m) 1600 lb 0 . 0 0.1. J L 0 .2 -J i_ 0 .3 _i _ 0 . 4 _: _ 0 . 5 _J L 0.6 J. -0 . 7 J - 0 . 8 _J L 0 . 9 i • 1.0 Figure 3-3. Model with 30 degree dip. (a) Earth model with constant velocity of 4000 m/s. (b) Synthetic seismogram for (a). Vertical axis is two-way traveltime in seconds and horizontal axis is offset coordinate with Ax=25 m. Arrows indicate positions of corner points in the model. 87 Figure 3-3 cont'd, (c) Time-image model and (d) Stolt migration result. Arrows indicate positions of corner points in the model. Q_ LU O 0 EARTH MODEL 400 800 1200 LINE COORDINATE (m) 1600 0.0 O J 0.2 0.3 X 0.4 0.5 1 0.6 1 0.7 0.8 J _ 0.9 1.0 Figure 3-4. Model with 45 degree dip. (a) Earth model with constant velocity of 4000 m/s. (b) Synthetic seismogram for (a). Vertical axis is two-way traveltime in seconds and horizontal axis is offset coordinate with Ax=25 m. Arrows indicate positions of corner points in the model. 89 Figure 3-4 cont'd, (c) Time-imaged model and (d) Stolt migration result. Arrows indicate positions of corner points in the model. 90 EARTH MODEL 0 2 0 0 I T 4 0 0 . X 1 6 0 0 Q_ LJJ 8 0 0 Q 1000 1200 0 4 0 0 8 0 0 1200 LINE COORDINATE (m) 1600 0.0. 0.1 1 0.2 I 0.3 0.4 0,5 0.6 1 0.7 0.8 1 0.9 1.0 Figure 3-5. Model with 60 degree dip. (a) Earth model with constant velocity of 4000 m/s. (b) Synthetic seismogram for (a). Vertical axis is two-way traveltime in seconds and horizontal axis is offset coordinate with Ax=25 m. 91 0 4 0 0 8 0 0 1 2 0 0 1 6 0 0 0.2 0,3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Figure 3-5 cont'd, (c) Time-imaged model and (d) Stolt migration result. Arrows indicate positions of corner points in the model. 92 Figures 3-4a and 3-4b show that a structural dip of 45 degrees is poorly modelled. The seismogram barely gives an outline of the dipping segment of the reflector. Studies in the second chapter predicted this breakdown and so it does not come as a surprise. Only weak amplitude, steep dip portions of the diffraction hyperbolas interfere con- structively to image the dipping segment so that a poor model is produced. Migration of the model (Figures 3-4c and 3-4d) show that the reconstructed time-image is still quite good, but the dipping portion of the reflector is of weak amplitude. An extreme case where the inclined portion of the structure is completely lost occurs when the structural dip is 60 degrees (Figures 3-5a and 3-5b). The modelled data would be incorrectly interpreted as two truncated reflectors at opposite ends of the seismic section. Stolt migration of the synthetic section (Figure 3-5d) cannot reconstruct the time-image model in Figure 3-5e. However, there does appear to be some recovery of structure between the two flat segments, although this is at the noise level. 3.1.3 Antiforrn and! synform models The next level of model complexity is that of antiforrn and synform structures. These features are of particular interest in petroleum exploration as they have been associated with potential zones of hydrocarbon accumulation. For example, antiforms may be oil bearing carbonate reefs and synforms may indicate gas charged channel sands. Figure 3-6 displays an antiforrn model and the corresponding synthetic seismogram. The modelled data show that the antiforrn is broadened beyond its true structural dimension. Incorrect interpretation of this feature would lead to an over estimate of the potential hydrocarbon reserves contained within the structure. The sharp corners of the structure give rise to diffractions. These results are expected. 93 E A R T H M O D E L x i— Q_ U J Q 4 0 0 8 0 0 1 2 0 0 L I N E C O O R D I N A T E (m) 1 6 0 0 o.o 0.1 1 0.3 0.4 0.5 0.6 X 0.7 1 0 . 8 0.9 1.0 Figure 3-6. Antiforrn model, (a) Earth model with constant velocity of 4000 m/s. (b) Synthetic seismogram for (a). Vertical axis is two-way traveltime in seconds and horizontal axis is offset coordinate with Ai=25 m. Arrows indicate positions of corner points in the model. 94 I now examine seismic sections for a suite of synform structures which have different depths of burial. When the model (Figure 3-7a) is such that the depth of burial (d) of the synform is less than the radius of curvature (r) of the synform, then the synthetic seismogram (Figure 3-7b) reflects the geometry of the structure. The amplitudes on the steep flanks of the synform are not modelled very well and diffractions from the sharp corners are difficult to observe. The latter deficiency occurs because of the shal- low depth. When the section is viewed from the side at an acute angle the large dip response is weak, but is definitely interpretable as a synform where d < r. In the case d = r (Figure 3 -8a) , reflections from the synform will all arrive at the receiver at the centre of the synform simultaneously. The seismogram will then show a very strong amplitude response at this offset position and the remainder of the synform will not be imaged. Again diffractions are expected from the sharp edges of the structure. These expectations are fulfilled quite well in Figure 3-8b. An interesting effect occurs on the synthetic seismogram for the synform model where d > r (Figure 3-9a) . In Figure 3-9b the synform is shown to be inverted. This is called the bow-tie effect and it is familiar to reflection seismologists. Again I examine the quality of the spatial wavelet approach by migration of the synthetic seismogram. I do this only for the synform with d> r. The time-image earth model is displayed in Figure 3-9c. The spatial wavelet modelled data were migrated (Figure 3-9d) with the Stolt method and an excellent reconstruction of the time-image structural model was obtained. The steeply dipping portions of the synform are not recovered completely and this is likely due to the short diffraction tails observed in Figure 3-9b. 95 EARTH MODEL 400 800 1200 LINE COORDINATE (m) 1600 0.0 4 0 J JL 0.2 1 0.3 0.4 0.5 1 0.6 0.7 0.8 1 0.9 1.0 Figure 3-7. Synform model with d < r. (a) Earth model with constant velocity of 4000 m/s. (b) Synthetic seismogram for (a). Vertical axis is two-way traveltime in seconds and horizontal axis is offset coordinate with Ai=25 m. Arrows indicate positions of corner points in the model. 96 EARTH MODEL Q_ UJ Q 400 800 1200 LINE COORDINATE (m) 1600 0 . 0 0.1 0 . 2 0 . 3 X 0 . 4 j _ 0 . 5 X 0,6 0 . 7 X 0 . 8 1 0 . 9 1.0 Figure 3-8. Synform model with d - r . (a) Earth model with constant velocity of 4000 m/s. (b) Synthetic seismogram for (a). Vertical axis is two-way traveltime in seconds and horizontal axis is offset coordinate with Ax-25 m. Arrows indicate positions of corner points in the model. 97 Q 0 EARTH MODEL 400 800 1200 LINE COORDINATE (m) 1600 0 . 0 0.1 0 . 2 0 . 3 0 . 4 0 . 5 X 0 . 6 0 . 7 0 . 8 0 . 9 1.0 Figure 3-9. Synform model with d > r. (a) Earth model with constant velocity of 4000 m/s. (b) Synthetic seismogram for (a). Vertical axis is two-way traveltime in seconds and horizontal axis is offset coordinate with Ax=25 m. Arrows indicate positions of corner points in the model. 98 0.0 O J 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Figure 3-9 cont'd, (c) Time-image model and (d) Stolt migration result. Arrows indicate positions of corner points in the model. 99 3 . 2 Horizontally stratified acoustic media A more realistic model of the subsurface velocity field is that of a horizontally stratified structure. The velocity may vary with depth, but it is laterally invariant, i.e., v = v(z). This is called the layered earth model and is supported by the empirical work of Faust (1951) which showed that compressional wave velocity varied with the burial depth and geologic age of the reflecting structure. Therefore, the layered earth model is a justifiable approximation to the real earth for the purpose of preliminary analysis at a level of complexity higher than that of homogeneous media. The viability of the spatial wavelet method in horizontally stratified acoustic media is illustrated briefly in this section. Reflector models in this portion are the same as those presented for homogeneous acoustic media and so the effects of vertical velocity inhomogeneity can be studied by comparison. The parameters used were as in the last section: frequency range over 10-80 Hz, Az=10 m and Ax=25 m. Like the dipping and synform examples I define the model with 129 points (Ax=12.5 m) and display seismic sections with every other trace so that there are only 65 traces. 3 . 2 . 1 A dipping model 1 modified the 30 degree dipping model (Figure 3-3a) by introducing a velocity jump at a depth of 730 m. The velocities above and below this depth are 2000 m/s and 4000 m/s respectively (Figure 3-10a). This velocity change occurs at the middle of the dipping segment in the model. The synthetic seismogram (Figure 3-10b) shows that the result is markedly different from that without the velocity stratification. Now a new apparent segment is observed with a high apparent dip. The large dip of this segment is due to the low velocity region in the upper part in the model. This low 100 EARTH MODEL 0 CL UJ o 1000 0 400 800 1200 LINE COORDINATE (m) 1600 Figure 3-10. Model with 30 degree dip. (a) Earth model with velocities of u1=4000 m/s and v2=2000 m/s below and above the 730 m depth, (b) Synthetic seismogram for (a). Vertical axis is two-way traveltime in seconds and horizontal axis is offset coordinate with Ax=25 m. Arrows indicate positions of corner points in the model. 101 Figure 3 - 1 0 cont'd, (c) Time-image model and (d) phase shift migration result. Arrows indicate positions of corner points in the model. 102 velocity causes the point diffractors to have steeply dipping hyperbola flanks which do not cancel out completely. As a result noise is observed below the steep segment. Increasing the number of points defining the model can reduce this effect. Initially, I computed this seismogram using only 65 points to define the model and found that doubling the number of points improved the seismogram quality. However, there was a large increase in the cost of computation. The cost of computing the section shown was approximately 130 seconds of central processing unit time. Therefore, I refrained from increasing the number of points any further. A time-image of the model is shown in Figure 3-10c. The modelled data were migrated using the phase shift method (Gazdag, 1978). Figure 3-10d illustrates that the reconstructed time-image of the model matches Figure 3-10c very well. However, the steep dipping segment is not reconstructed as well as the other segments of the model. As previously mentioned, this is probably due to the poor representation of waves with large emergence angles. 3.2.2 A synform model Earlier, in Figure 3-9b I showed the well known bow-tie effect. Now I consider a model like that in Figure 3-9a, but with a velocity change at a depth of 630 m. Above this depth the velocity is 2000 m/s, and below it, the velocity is 4000 m/s. The velocity discontinuity occurs through the synform structure (Figure 3-lla). In Figure 3-llb the seismic model shows that with the introduction of the velocity stratification the bow- tie effect has been lost. In its place an image that might be interpreted as a synform with depth of burial less than the radius of curvature of the bowl, is observed. Clearly, 103 0 EARTH MODEL 200 - 400 - 800 - 1 0 0 0 L ^ _ J 1 1 L 1 1 1 0 400 800 1200 1600 LINE COORDINATE (m) Figure 3 - 1 1 . Synform model with d > r. (a) Earth model with velocities of u1=4000 m/s and u2=2000 m/s below and above the 630 m depth, (b) Synthetic seismogram for (a). Vertical axis is two-way traveltime in seconds and horizontal axis is offset coordinate with Ax=25 m. Arrows indicate positions of corner points in the model. 104 0 400 800 1200 Y 1600 0.0 J . 0.1 1 0 .2 JL 0 .3 JL O.4 j . 0 .5 JL 0.6 1 0.7 JL 0.8 JL 0 .9 1 1.0 Figure 3-11 cont'd, (c) Time-image model and (d) phase shift migration result. Arrows indicate positions of corner points in the model. 105 vertical variations in the velocity may have important implications on the observed reflection structure. A time-image model of Figure 3-1 la is given in Figure 3-llc. Reconstruction of the time-image model from the modelled data was again accomplished by the phase shift method. The migration result is shown in Figure 3-lld. The reconstruction is excellent, but with the loss of detail in areas of steeply dipping structure. 3.3 Laterally varying acoustic media A more complicated velocity model employed by geophysicists is one in which the veloc- ity may vary vertically and laterally, i.e., v = v(x, z). Computing seismograms in these types of media is a difficult problem and has been a topic of ongoing research over the past decade. Usually mathematical formulations are limited in their ability to describe wave propagation in complex velocity structures. Therefore, approximate approaches become necessary if a solution is to be obtained. In the first chapter, I mentioned various methods by which lateral velocity inhomo- geneities can be modelled in the seismic reflection process. These methods are usually applicable when the velocity inhomogeneity is weak. When large velocity contrasts oc- cur over small distances, all methods produce poor results. In this section I present a method of accomodating lateral velocity variations which uses the spatial wavelet concept and discuss some of the limitations with this approach. 106 3.3.1 Accomodating lateral velocity changes Berkhout and van Wulfften Palthe (1979, 1980) and Berkhout (1981, 1982, 1984) pro- posed an approximate method to model wave propagation in laterally varying media using spatial convolution. They suggested that the convolution be dependent on the lo- cal velocity. That is, the convolution should be made space variant with wavelet changes being governed by a localized velocity function. However, these works described neither how to find a localized velocity function nor how to implement space variant convolution for the purpose of seismic reflection modelling, although they did intimate that such a procedure might only be valid when the velocity contrast is small so that transmission losses and internal reflections can be ignored. A localized velocity function may be calculated in many ways. The particular method used to compute such a function is not crucial provided that the approach smooths the transition between large lateral velocity contrasts. However, some methods may be simpler and/or more flexible than others. I compute localized velocity functions using a simple mean, weighted by a decaying exponential function. The eomputaional expression is f_L v(x)ws(xo — x)dx J-L ws(zo - x)dx M*o) = L r L , 7 7 — (3.2) where vs(x) is the smoothed velocity function, ws(x) is the weighting function and L is the aperture length of the wavelet. The smoothing integral is over the aperture because that is the extent to which the spatial wavelet is attempting to model wave propagation. The reason for selecting the exponential function was the ease of designing a flexible window. For example, a ramp function would be adequate, but the decay rate (slope) would be fixed for a given wavelet aperture. In contrast, an exponential function can 107 have different rates of decay when the wavelet aperture is constant. The particular weighting function that I used is given by ws(x) = e-^rf, (3.3) where n is the decay rate or smoothing factor. Small values of 77 indicate fast decay and minimal smoothing whereas large values express the opposite case. With this weighting function the localized or smoothed velocity function is determined by, Given a smoothed velocity function, space variant convolution becomes simple to implement. The variation of the wavelet is determined by the velocity function along the x direction. • Based on the results of the velocity smoothing, all the necessary wavelets for wavefield continuation can be computed in advance of convolution. This constitutes an overhead cost. The computed wavelets are stored and recovered as needed for the space variant or point-wise convolution. Therefore, the increase in cost of computation depends only on the number of wavelets that are needed for the different velocities. The examples that I compute in this section are small and therefore wavelets for all the necessary velocities can be computed. In larger examples this may not be possible due to computer limitations and expense. For such cases a new strategy is needed to decide on velocity selection. The expense of the wavefield continuation procedure is unaffected by the space variant nature of the the convolution and so a simple, approximate approach to modelling complex structure becomes available. vs{x0]r)) f^Lv{x)e-Ux<>-xVrfdx (3.4) I-L e-\{x»-*)h\2dx 108 3.3.2 A simple experimental model In modelling lateral velocity variations it is extremely important to consider the gradi- ent of the velocity transition. The gradient is simply the rate of change of velocity with distance and it has units of frequency (Hz). If seismograms with lateral velocity transi- tions are to be modelled by spatial wavelets, then the lowest frequency component of the wavelets must be greater than the largest local velocity gradient. Otherwise, wavelets with too low a frequency will not be able to describe the velocity change. In cases where large velocity gradients exist, a solution may not be possible. However, if the velocity transition is smoothed sufficiently to reduce the gradient, then adequate seismograms may be computed. In this section I invesitgate the effects of lateral velocity smoothing. I employ a simple experimental model to study seismic modelling in laterally varying velocity structures. The model consists of a point diffraetor embedded in a velocity structure that varies laterally. There model space is bisected into two velocity zones, v\ and V2 with the source being located near the velocity boundary in either the vv or v2 material. In Figure 3-12 the theoretical expectation of the types of' waves this geometry produces when the source is in the slower medium are shown. The spatial wavelet formulation cannot model the critically refracted wave (3) or the internally reflected wave (2), but it does model the direct wave (l) and the transmitted wave (4). Velocity smoothing was investigated using an earth model like that shown in Figure 3-12. In particular, the source is located 50 m away from the velocity boundary at a depth of 500 m in the slow velocity medium. The lower velocity is ui=4000 m/s and the higher velocity is ̂ 2=5000 m/s. In computing the seismograms I used a frequency range of 10-40 Hz, trace separation of 25 m over an 800 m line and an extrapolation increment of 10 m. 109 Figure 3-12. Waves generated near a lateral velocity discontinuity. When the source (S) is in the slower velocity medium four types of compressional waves are expected: (1) the direct wave, (2) the internally reflected wave, (3) the critically refracted wave and (4) the transmitted wave. If the source were in the faster medium the type (3) wave would not be generated. The spatial wavelet method is able to describe only wave types (1) and (4). 110 In Figures 3-13a and 3-13b the unsmoothed velocity model and computed seismic section are shown, respectively. The velocity transition occurs over 25 m (the spacing of a single trace) and this results in a lateral velocity gradient of Vzt;(x)=40 Hz. Since the frequency range used is 10-40 Hz, the diffraction formulation presented here is not expected to adequately describe the rapid velocity variation. For proper modelling, frequencies higher than 40 Hz would have to be used. Regardless, some interesting effects can be seen and explained even in this case. First, there is the marked change in the shape of the characteristic hyperbola. In a constant velocity medium the hyperbola is symmetric with the location of the source being identified with the position of the hyperbola apex. However, in the presence of a lateral velocity change, the hyperbola becomes distored and it is no longer symmetric. Nor is the position of the earliest traveltime correlated with the lateral location of the point diffractor (Figure 3-13b). Second, there is a change in the amplitude of the seismic signature on either side of the velocity boundary. According to equation (2.25) the waves in the faster medium will have lower amplitudes due to a l / r factor which will effect the waves in the fast medium because of the increased travel path attributed to refraction at the velocity discontinuity (cf., wave type (4) shown in Figure 3-12). Also, the traveltime reconstruction is fairly good (Figure 3-13e). I made the comparison with the results from simple ray tracing of the directly scattered and refracted waves. Finally, an artifact of the processing is the second apparent reflector directly beneath the primary response. This occurs due to the smearing of wavelets as they cross the velocity boundary. Initially I thought that this might be due to internal reflections, but simple ray tracing considerations for this geometry indicated no such arrivals as observed in the seismogram. I concluded that this was a spurious effect of the computation. I l l VELOCITY MODEL 4800 4600 >- 4400 O o _ l 4200 L J > 4000 3800 0 200 400 600 LINE COORDINATE (m) 800 F V S Y Y Y 0.0 0.1 4 0.2 J L 0 .3 1 0.4 0.5 Figure 3-13. Seismogram for a point source located near a lateral velocity discontinu- ity, (a) Lateral velocity model without smoothing. Maximum local velocity gradient of 40 Hz. Boxes indicate actual velocities used in computations, (b) Synthetic seismogram. F, V and S indicate positions of first arrival, velocity boundary and source respectively. 112 0.0 Figure 3-13 cont'd, (c) Comparison of traveltimes. The result of (b) is shown with a traveltime curve obtained by ray tracing (solid line). The spatial wavelet approach does well in traveltime reconstruction, but displays the familiar -TT/4 phase shift. The direct wave (l) and the transmitted wave (4) are marked as indicated. 113 An alternative to increasing the lowest frequency limit exists, where the same fre- quency band may be retained and the velocity transition is smoothed instead. Physi- cally, this approach is appealing since geology rarely changes rapidly as depicted in the previous example, except for faults. It is often that one medium will gradually merge with another over a large distance. To model this effect I smoothed the velocity model using smoothing factors of n = l (Figures 3-14a and 3-14b), r}=2 (Figures 3-14c and 3- 14d) and n=3 (Figures 3-14e and 3-14f). In Figure 3-14a (f?=l) the average value of Vxv(x) is 8 Hz, but the maximum local value is 22.4 Hz (at the centre of the transition). Therefore, the spatial wavelets frequency band 10-20 Hz will not be able to describe the velocity transition sufficiently well. The seismogram (Figure 3-14b) confirms this by showing that there is little change from the case without smoothing. For a larger amount of smoothing (?7=2) the average Vxv(x) value is 3.6 Hz, and the maximum local Vxv(x) is 11.6 Hz (Figure 3-14c). Now most of the wavelets in the frequency band can describe the velocity change sufficiently well as can be seen in Figure 3-14d. The velocity model with further smoothing (?7=3) is shown in Figure 3-14e. The average VIt>(x)=2.6 Hz and the local maximum Vxv(x) is 7.2 Hz.-In this case all the spatial wavelets in the frequency band can contribute to the description of waves travelling across the lateral boundary. The smooth transition is reflected in the high interpretational quality of the seismogram in Figure 3-14f. An expectation from geometrical ray theory indicates that good seismograms can be computed when the largest local velocity gradient is less than the the lowest component of the frequency band (Yedlin, pers. comm.). The examples shown in this section indicate that the limitations of the spatial convolution approach to modelling wave propagation in laterally varying media are similar to those for geometrical ray theory. I found that although the spatial convolution method was able to give what appeared to 114 VELOCITY MODEL 3800 0 200 400 600 LINE COORDINATE (m) 800 Figure 3-14. Seismograms for a point source near a smooth lateral velocity transition, (a) Lateral velocity model with r/=l. Maximum local velocity gradient of 22.4 Hz. Boxes indicate actual velocities used in computations, (b) Synthetic seismogram. F , V and S indicate positions of first arrival, original velocity boundary and source respectively. Outer arrows indicate extent of smoothed velocity transition zone. 115 VELOCITY MODEL 4800 4600 >- 4400 O o _ J 4200 LJ > 4000 3800 0 200 400 600 LINE COORDINATE (m) 800 Figure 3-14 cont'd, (c) Lateral velocity model with r?=2. Maximum local velocity gradient of 11.6 Hz. (d) Synthetic seismogram. 116 VELOCITY MODEL >- O O _ J UJ > 4800 4600 4400 4200 4000 3800 0 200 400 600 LINE COORDINATE (m) 800 Figure 3-14 cont'd, (e) Lateral velocity model with r/=3. Maximum local velocity gradient of 7.2 Hz. (f) Synthetic seismogram. 117 be reasonable results for a velocity boundary without smoothing, the best seismograms were for the cases where there was significant smoothing. The measures of a reasonable seismogram that I used were continuity and smoothness of the seismic signature, which are important in terms of the interpretational quality. Unlike modelling in homogeneous and horizontally stratified velocity structures, it is extremely difficult to evaluate the quality of seismograms for waves in laterally varying media. The primary difficulty at the time of this writing was the lack of an analytic or exact numerical solution by other means with which to compare the spatial convolution method. I terminate my investigation at this stage, but I provide suggestions for further research on this topic in the next chapter. 118 4. C O N C L U S I O N 4.1 Summary The spatial convolution approach of describing wave propagation is a conceptually sim- ple method based on the classical diffraction theory of Kirchhoff. In this thesis I have investigated two main themes concerning this method: the computational aspects of spatial wavelets and their application to modelling seismic reflection data for the 2-D zero-offset case. Either of these topics, when studied in great depth could have made excellent projects in themselves. However, one without the other would have left a con- siderable gap between theory and practice. Therefore, I have endeavoured to present an overall viewpoint combining theoretical discussion, computational detail and appli- cation. The advantages of this approach were that a large amount of material was covered and that an integrated platform from which to launch new investigations was established. Investigations of the first theme produced two new developments. The first result was a fast, flexible and conceptually simple approach to computing spatial wavelets while preventing spatial aliasing. I showed that computing spatial wavelets in the space domain leads to aliasing if the extrapolation step size was chosen too small or if temporal frequencies beyond that corresponding to the Nyquist wavenumber were considered. In earlier work Berkhout (1981) also recognized these difficulties and proposed a method to circumvent the hazards of selecting an extrapolation step that was too small. He formulated an approach using a truncated Taylor series expansion of a wavelet. I found this method difficult to follow and as an alternative I proposed the spectral method 119 of computation. The spectral method relies on the fact that the spatial wavelet and the phase shift operator (Gazdag, 1978) are a Fourier transform pair. Therefore, when potential aliasing conditions arise the wavelet may be bandlimited in the wavenumber domain to prevent aliasing. In my study I found that a fourth order Butterworth window gave a suitable bandlimiting operator. The effect of using a low order window was to smooth the wavenumber representation of the wavelet so that the information near the Nyquist wavenumber is zero. As well as preventing aliasing this minimized truncation errors that may occur during Fourier transformation from the wavenumber domain to the space domain. However, such windowing also rejected information at wavenumbers below the Nyquist limit and when the wavelet is used recursively this effect is compounded. As a consequence the seismograms computed in this manner suffered from dip limitations. Thus, a tradeoff exists between flexibility and ease of computation with restrictions on the imaging of steeply dipping structures. The second development was the determination of the minimum length or aperture of a wavelet. By necessity the wavelet must be truncated at a certain number of samples; for computational purposes. The maximum number of samples in a wavelet is equal to the number of traces that are to be used in a seismic section. Often the number of traces is large and recursive convolution with a long wavelet becomes very expensive. In addition, if laterally varying velocity structures are to be modelled, a short wavelet is necessary. With a view to finding the shortest wavelet, I argued that if the discrete representation of a continuous function contained a large portion of the total energy in a finite number of samples, then truncation of the series at that number of samples would be sufficient to minimize information loss. I found that the minimum aperture wavelet depended on three parameters: the extrapolation step size, the medium velocity and the temporal frequency. The energy variation with each of these parameters as a 120 function of the number of samples was studied. Results showed that the energy rise (to a certain preselected threshold value) with the number of samples was; (l) faster for small extrapolation steps sizes than for large step sizes; (2) faster for low velocities than for high velocities; and (3) faster for high temporal frequencies than for low temporal frequencies. Therefore, the minimum aperture wavelet could be computed using the smallest desired extrapolation step size, the lowest frequency of interest and the highest velocity of interest. For practical purposes a standard minimum aperture was selected so that one length was applied to all wavelets regardless of the frequency or velocity. The maximum of the minimum aperture sizes was used for this purpose. This aperture length was determined by selecting the lowest frequency of interest and the highest velocity of interest so that all wavelets satisfied the energy threshold criterion. Tests with wavelet truncation indicated that this procedure reduced computational expenses by a factor of four to five without severe side-effects. The second theme of this thesis demonstrated the successful application of spatial convolution to seismic reflection modelling in media of varying complexity. To the best of my knowledge these results are new contributions. I examined seismograms for models at three levels of complexity: homogeneous acoustic media, horizontally stratified acoustic media and laterally varying media. An acoustic medium is one that supports the propagation of compressional waves, but not of shear waves. Homogeneous acoustic media represent the lowest level of model complexity. Al- though this is not a very realistic model of the earth, it has been successfully employed in demonstrating the fundemental aspects of seismic modelling. I used homogeneous ve- locity models to explore simple reflector configurations. In particular, reflector models 121 representing truncated beds, structures of varying dips and antiform/synform struc- tures were investigated. The truncated reflectors example showed that the evolution of the wavefield as a function of extrapolation steps affected the diffraction response on the seismograms. Diffractions due to deep sources are more evolved than those from shallow sources. Migration of the modelled data demonstrated that even poorly evolved diffractions collapse to a focus fairly well. The dipping reflector models served to illus- trate the dip imaging limitations of the spatial wavelet method. Examples with dips of 15 and 30 degrees were easily accomodated, but cases with 45 and 60 degree dips gave poor results. These observations were consistent with the developments in the second chapter. Seismic migration of the shallow dip example (30 degrees) gave an excellent traveltime model reconstruction. Even for the steep dip cases the modelled data were migrated well enough to image portions of the 45 degree dipping segment and the 60 degree segment to a much lesser extent. The antiform/synform examples showed the applicability of the recursive spatial wavelet algorithm' in modelling curved reflectors. The antiforrn example produced the exaggerated doming effect clearly. Synelines buried at various depths also gave excellent results. When the depth of burial was less than the radius of curvature of the synform bowl, the entire reflection structure was imaged. In the case where the burial depth equalled the radius of curvature, the bowl was imaged as a focus beneath the centre of the array. The well known bow-tie effect was observed when the depth of burial exceeded the radius of curvature of the bowl. This example was migrated and a very good traveltime reconstruction was obtained. However, the steeply dipping poritons of the bowl were not imaged as well as other parts of the reflector. Horizontally stratified acoustic media represent the next level of model complexity. This layered earth model is more realistic than the homogeneous case and is founded on the empirical work of Faust (1951). I used horizontally stratified velocity models 122 to demonstrate the effect of such velocity distributions on the reflection characteristics of modelled data. This was done by introducing a low velocity layer to homogeneous models. The homogeneous models used for comparison were the 30 degree dipping reflector and the synform with depth of burial greater than the radiaus of curvature of the bowl. The 30 degree dipping model was changed by introducing a low velocity (2000 m/s compared to 4000 m/s) layer half way through the dipping segment. The seismogram showed that this modification resulted in a kink in the dipping portion of the modelled data. Migration of this example by the phase shift method produced an excellent time-image model of the structure. The synform case gave equally dramatic results. Initially, the lack of velocity stratification produced the bow-tie effect, but with the introduction of a low velocity layer half way through the bowl, the bow-tie is annihilated. The modelled data were migrated to give a good reconstructed image. However, as with previous examples, the steeply dipping parts of the structure were not imaged' well'. Finally, an extension of the spatial convolution approach to laterally varying media. This model of the earth is, perhaps, the most realistic type and the most difficult to use. However, it is becoming more commonly employed in the analysis of seismic reflection data, while remaining an active area of research (Claerbout, 1985b). Berkhout and van Wulfftem Palthe (1979, 1980) and Berkhout (1981, 1982, 1984) suggested that lateral velocity variations could' be accomodated by smoothing the lateral velocity field and' that the convolution be made space variant. The wavelet used in the convolution would be governed by the smoothed or local velocity. I used a simple mean, weigthed over the aperture length by a decaying exponential function, to compute smoothed velocity profiles. The earth model used to examine this formulation was equally simple. The velocity model simulated a large vertical discontinuity with differing velocities on either 123 side. A point source was embedded in the slower velocity medium. Seismograms were computed for cases with no velocity smoothing and for various amounts of smoothing. The traveltimes for the seismogram without smoothing compared favourably with those generated by simple ray tracing. However, the character of the seismic response was nei- ther smooth nor continuous across the boundary. Experiments with velocity smoothing clearly demonstrated that when the largest velocity gradient corresponded to the lower end of the temporal frequency band, the seismograms produced smooth and continu- ous reflection signatures. Also, increased smoothing results in lower velocity gradients. These observations are consistent with the results of geometrical ray theory. The work on using spatial convolution in laterally varying media shows that this for- mulation is worth pursuing even though it handles only upgoing wavefield's and produces results for zero-offset sections. With this in mind, I recommend further investigaton. 4.2 Further investigations Clearly, in this first attempt at understanding some basic concepts and demonstrating their use, I have not answered all the questions posed by earlier analyses. Indeed, this project has raised new questions, some of which may warrant further study. I give a brief indication of some possible topics for later investigations. First, stability studies might prove useful to better understand the relationships be- tween the number of extrapolation steps, the type of antialiasing window implemented, wavelet truncation, dip limitation and how these affect the computation of seismograms. The basis for this was laid in the second chapter. I think this is a difficult problem to subject to a quantitative analysis, but qualitative cataloging of the effects (although 1 2 4 tedious) may be helpful in providing a deeper appreciation of the spatial wavelet's properties in a homogeneous medium. Second, an attempt at making the computation more efficient would be fruitful. Presently the cost of computing large models (e.g., Figure 3-10b) is about 130 seconds of central processing unit time on an Amdahl 580/V6 system. I am uncertain as how any improvement can be made, since the main cost is incurred by the recursive convolution. If possible, implementation of the algorithm on an array processor might prove fruitful. Third, I think that extensions of the spatial convolution approach to laterally varying media could be useful. I have shown that the method has some potential, but limitations on computing power restricted my investigation. Even so, it might be useful to compare the spatial wavelet approach with various other algorithms (e.g., finite-difference or Kirchhoff if they are easily available) for simple cases. Also, development of a more sophisticated algorithm for purposes of comparison with scale model studies would be most beneficial. 125 R E F E R E N C E S Alford, R. M. , Kelly, K. R. and Boore, D. M., 1974. Accuracy of finite-difference modeling of the acoustic wave equation, Geophysics, 3 9 , 843-842. Alterman, Z. S. and Karal, F. C , 1968. Propagation of elastic waves in layered media by finite-difference methods, Bull, seism. Soc. Am., 5 8 , 367-398. Alterman, Z. S. and Loewenthal, D., 1970. Seismic waves in a quarter and three- quarter plane, Geophys. J. R. astr. Soc, 2 0 , 101-126. Angona, F. A., 1960. Two-dimensional modeling and its application to seismic problems, Geophysics, 25, 468-482. Bamberger, A., Chavent, G. and Lailly, P., 1980. Etude de schemas numeriques pour les equations de l'elastodynamique hneaire, Institut National de Recherche en Informatique et en Automatique, 0 3 2 - 7 9 . Berkhout, A. J., 1981. Wave field extrapolation techniques in seismic migration: a tutorial, Geophysics, 4 6 , 1638-1656. Berkhout, A. J., 1982. Seismic Migration, Developments in Solid Earth Geophysics, 14A, Elsevier, Amsterdam. Berkhout, A. J.,1984. Seismic Migration, Developments in Solid Earth Geophysics, 1 4 B , Elsevier, Amsterdam. Berkhout, A. J. and van Wulfften Palthe, D. W., 1979. Migration in terms of spatial deconvolution, Geophys. Prosp., 2 7 , 261-291. Berkhout, A. J. and van Wulfften Palthe, D. W., 1980. Migration in the presence of noise, Geophys. Prosp., 2 8 , 372-383. Berryhill, J. R., 1977. Diffraction response for nonzero separation of source and receiver, Gephysics, 4 2 , 1158-1176. Berryhill, J. R., 1979. Wave-equation datuming, Geophysics, 44,1329-1344. Berryman, L. H., Goupillaud, P., L. and Waters, K. H., 1958a. Reflections from multiple transition layer, part I: theoretical results, Geophysics, 23, 223-243. 126 Berryman, L. H., Goupillaud, P., L. and Waters, K. H., 1958b. Reflections from multiple transition layer, part II: experimental results, Geophysics, 2 3 , 244-252. Boore, D. M., 1972. Finite-difference methods for seismic wave propagation in het- erogeneous materials, in Methods in Computational Physics, 1 1 . Academic Press, New York. Born, M. and Wolf, E. , 1975. Principles of Optics, 5th ed., Pergamon Press, Oxford. Carter, J. A. and Frazer, L. N., 1983. A method for modeling reflection data from media with lateral velocity changes, J. geophys. Res., 8 8 , 6469-6476. Cerveny, V., Molotkov, I. A. and Psencik, I., 1977. Ray Method in Seismology, Univerzita Karlova, Praha. Claerbout, J. F., 1985a. Fundementals of Geophysical Data Processing, Blackwell, Oxford. Claerbout, J. F., 1985b. Imaging the Earth's Interior, Blackwell, Oxford. Deregowski, S. M. and Brown, S. M. , 1983. A theory of acoustic diffractors applied to 2-D models, Geophys. Prosp., 3 1 , 293-333. Dobrin, M. B., 1976. Introduction to Geophysical Prospecting, 3rd. ed., McGraw- Hill, New York. Evans, J . F., Hadley, C. F., Eisler, J . D. and Silverman, D., 1954. A three- dimensional wave model with both electrical and visual observation of waves, Geo- physics, 1 9 , 220-236. Faust, L. T., 1951. Seismic velocity as a function of depth and geologic time, Geophysics, 1 6 , 192-206. French, W. S., 1974. Two-dimensional and three-dimensional migration of model- experiment reflection profiles, Geophysics, 3 9 , 265-277. Gazdag, J., 1978. Wave equation migration with the phase shift method, Geo- physics, 4 3 , 1342-1351. Gazdag, J. , 1981. Modeling of the acoustic wave equation by transform methods, Geophysics, 4 6 , 854-859. 127 Gazdag, J. and Sguazzero, P., 1984. Migration of seismic data, Proc. I.E.E.E., 72, 1302-1315. Goodman, J. W., 1968. Introduction to Fourier Optics, McGraw-Hill, New York. Gradshteyn, I. S. and Ryzhik, I. M. , 1965. Table of Integrals, Series and Products, 4th ed., Academic Press, New York. Hatton, L., Worthington, M. H. and Makin, J. , 1986. Seismic Data Processing, Blackwell, Oxford. Hilterman, F. J., 1970. Three-dimensional seismic modeling, Geophysics, 35, 1020- 1037. Hilterman, F. J., 1975. Amplitudes of seismic waves: a quick look, Geophysics, 40, 745-762. Hilterman, F. J. and Larson, D. E. , 1975. Kirchhoff wave theory for multi-velocity models, paper presented at 45th annual Society of Exploration Geophysieists meet- ing, Denver, Colorado. Howes, E. T., Tejada-Flores, L. H. and Randolph, L., 1953. Seismic model study, JL acoust. Soc. Am-., 25, 915-921. Jackson, J. D., 1962. Classical Electrodynamics, John Wiley and Sons, New York. Kanasewieh, El R., 1981. Time Sequence Analysis in Geophysics, 3rd ed., University of Alberta Press, Edmonton. Kaplan, W., 1981. Advanced Mathematics for Engineers, Addison-Wesley, Reading. Kelamis, P. G. and Kjartansson, E. , 1985. Forward modeling in the frequency-space domain,, Geophys. Prosp., 33, 252-262. Keller, J. B., 1962. Geometrical theory of diffraction, J . opt. Soc. Am., 52, 116-130. Kelly, K. R., Ward, R. W., Treitel, S. and Alford, R. M., 1976. Synthetic seis- morams: a finite-difference approach, Geophysics, 41, 2-27. Kelly, K. R., Alford, R. M. and Whitmore, N. D., 1982. Modeling: the forward method, in Concepts and Techniques in Oil and Gas Exploration, Society of Explo- ration Geophysieists. Tulsa. 128 Kinsler, L. E. , Frey, A. R., Coppens, A. B., and Sanders, J. V., 1982. Fundementals of Acoustics, 3rd ed., John Wiley and Sons, New York. Kosloff, D. D. and Baysal, E., 1982. Forward modeling by a Fourier method, Geo- physics, 47, 1402-1412. Kosloff, D., Reshef, M. and Loewenthal, D., 1984. Elastic wave calculations by the Fourier method, Bull, seism. Soc. Am., 74, 875-891. Kuhn, M. J. and Alhilali, K. A., 1976. Weighting factors in the construction and reconstruction of acoustical wavefields, Geophysics, 42, 1183-1198. Loewenthal, D., Lu, L., Roberson, R. and Sherwood, J. , 1976. The wave equation applied to migration, Geophys. Prosp., 24, 380-399. Magnus, W and Oberhettinger, F., 1949. Formulas and Theorems for the Special Functions of Mathematical Physics, Chelsea , New York. Marfurt, K. J., 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations, Geophysics, 49, 533-549. McQuillin, R., Bacon, M. and Barclay, W., 1984. An Introduction to Seismic Inter- pretation, 2nd ed., Gulf, Houston. Oldenburg, D. W., 1981. A comprehensive solution to the linear deconvolution problem, Geophys. J. R. astr. Soc, 65, 331-357. Oliver, J. , Press, F. and Ewing, M., 1954. Two-dimensional seismology, Geophysics, 19, 202-219. Oppenheim, A. V. and Schafer, R. W., 1975. Digital Signal Processing, Prentice- Hall, Englewood Cliffs. Reshef, M. and Kosloff, D., 1985. Applications of elastic forward modeling to seismic interpretation, Geophysics, 50, 1266-1272. Robinson, E. A., 1957. Predictive decomposition of seismic traces, Geophysics, 22, 767-778. Robinson, E. A., 1983. Migration of Geophysical Data, International Human Re- sources Development Corp., Boston. 129 Robinson,E. A. and Treitel, S., 1980. Geophysical Signal Analysis, Prentice-Hall, Englewood Cliffs. Schneider, W. A., 1978. Integral formulation for migration in two and three dimen- sions, Geophysics, 43, 49-76. Sheriff, R. E. and Geldart, L. P,. 1982a. Exploration Seismology, vol. 1, Cambridge University Press, Cambridge. Sheriff, R. E. and Geldart, L. P., 1982b. Exploration Seismology, vol. 2, Cambridge University Press, Cambridge. Shtivelman, V., 1984. A hybrid method for wave field computation, Geophys. Prosp., 32, 236-257. Shtivelman, V., 1985. Two-dimensional acoustic modeling by a hybrid method, Geophysics, 50, 1273-1284. Stolt, R., 1978. Migration by Fourier transform, Geophysics. 43, 23-48. Tatham, R. H. and Stoffa, P. L., 1976. Vp/V8: A potential hydrocarbon indicator, Geophysics, 41, 837-849. Trorey, A. W., 1970. A simple theory for seismic diffractions, Geophysics, 35, 762- 784. Trorey, A. W., 1977. Diffractions for arbitrary receiver-source locations, Geophysics, 42, 1177-1182. Waters, K. H., 1978. Reflection Seismology, John Wiley and Sons, New York. Woods, J. P., 1956. Composition of reflections, Geophysics, 21, 261-176. 130 A P P E N D I X Butterworth windows and windowing strategy The term Butterworth window is somewhat unconventional, while the usage of But- terworth filter is more common. These two are related and perform the same tasks, but their implementations are very different. Butterworth filters perform filtering by a recursive formulation in the space domain. In contrast, Butterworth windows perform filtering by multiplication in the wavenumber domain. A Butterworth window is related to a Butterworth filter in that the Fourier amplitude spectrum of the later gives the former. Butterworth windows may be implemented for lowpass, highpass and bandpass filr tering purposes. In this thesis, only the lowpass operation was required. The general form of the Butterworth window in the wavenumber domain for the lowpass application is given by B(kx\ = 1 = , (A.I) where kx is the horizontal wavenumber, kc is the cutoff wavenumber and N is the order of the window. Butterworth windows are defined by the property that their amplitude is maximally flat in the passband. For an Nth order lowpass window, this means that the first 27V-] derivatives of the amplitude function are zero at kx—0. Another property is that the amplitude is monotonic in the passband and in the stopband (Oppenheim and Schafer, 1975, p.21l). 131 As the order, N, in equation (A.i) increases, the window characteristics become sharper; i.e., the amplitudes remain closer to unity over the passband and become closer to zero more rapidly in the stopband. The amplitude at the cutoff wavenumber, kc, will always be l/\/2 because of the nature of equation (A.i). The dependence of Butterworth window characteristics on the order is indicated in Figure A-1. There are two aspects to be considered when applying the lowpass filter. First, the order of the window must be selected. Second, the mode of implementation for the chosen window must be determined. Figure A-1 provides a guide to choosing window order. The amplitude characteristics show that high order windows are desirable because of their sharpness in the wavenumber domain. However, these windows also produce large oscillations in the space domain. This could produce large errors in the space domain computation procedure. In contrast, the low order windows have large transition regions in the wavenumber domain, but they have small oscillations in the space domain. With this qualitative guide, a judicious choice of filter order must be made. I decided to use JV/'=4, and found that this choice gave fairly good results. Another consideration that must be accounted for is the dependence of the window characteristics on the number of extrapolation steps as the wavefield is upward continued. This is illustrated for a fourth order window in Figure A-2. Every step that the wavefield is extrapolated through implies that the window is multiplied by itself in the wavenumber domain. Consequently, as the number of extrapolation steps becomes large, window amplitudes in the transition region and the stopband decay. If this decay occurs in the region where waves are propagating (i.e., kx < k), then the result is a loss in the ability to image steep dips. Selection of a large order window will reduce this problem, but it may enhance the likelihood of accumulating errors due to wavelet truncation in the space domain convolution operation and possibly provide a worse image. Thus, the 132 W I N D O W V A R I A T I O N WITH O R D E R 0 . 0 0 . 0 4 0 . 0 8 0 . 1 2 W A V E N U M B E R ( r a d / m ) W I N D O W V A R I A T I O N WITH O R D E R UJ Q ZD CL < 6 0 1 2 0 D I S T A N C E ( m ) 1 8 0 Figure A-1. Dependence of lowpass Butterworth amplitude characteristic on the order, N. (a) Butterworth windows with A;c=0.10 rad/m in the wavenumber domain. Increas- ing N produces sharper characteristics. The values of A' used were 2 (smoothest), 4, 6, 8, 10 and 390 (sharpest). As N becomes very large, the window resembles a box-car function, (b) Space domain representation of windows in (a). Oscillations with larger amplitudes .correspond to larger order windows. These windows were computed for a spatial sampling interval of Ax — 25 m as used for the experiments in the body of the thesis. The Nyquist wavenumber is kN =0.125 rad/m. 133 WINDOW VARIATION WITH RECURSION 0.0 0.04 0.08 0.12 WAVENUMBER ( r a d / m ) Figure A-2. Dependence of fourth order Butterworth amplitude characteristic on the number of recursive extrapolation steps. As the number of recursions increases, the amplitudes in the transition zone and the stopband decay. Curves were computed for 1, 50 and 100 extrapolation steps. The cutoff wavenumber initially selected was fcc=0.10 rad/m, but with an increasing number of recursions the effective cutoff is near 0.06 rad/m. These windows were computed for a spatial sampling interval of Ax = 25 m as used for the experiments in the body of the thesis. The Nyquist wavenumber is fcjv =0.125 rad/m. 134 tradeoff between computational stability and dip resolution must be considered when making a best guess at the order to use. I reiterate the tradeoff. Large order windows are desirable because they do not winnow away a large portion of the propagating wavefield, but they may cause instability. Small order windows are desirable because they tend to be computationally stable, but they may cause excessive decay of the propagating wavefield. The implementation of the chosen window constitutes the second stage of windowing strategy. I found that in attempting to preserve the consistency of wavelet energy de- terminations for minimum aperture, that a varying cutoff, fourth order window worked well. Consistency refers to the ability to numerically confirm that the smallest aperture corresponds to the wavelet with the lowest temporal frequency and the highest velocity (cf., Section 2.3). In such an implementation the cutoff wavenumber changes with the propagation wavenumber of the wavelet being computed. In paticular, the rule that I employed was kc — 0.8k if k < and kc — 0.8k otherwise. This effectively removed1 the evanescent portion of the wavefield described by the spatial wavelet. The results of this implementation are shown with examples in Figure A-3. Figure A-3a shows the characteristic hyperbola for a point source at a depth of 500 m. However, there is a small artifact, directly beneath the main hyperbola, which appears as two mini-hyperbolas. The origin of this effect will become apparent with further examples. The modelling used the following parameters: Az=10 m, Az=25 m and a frequency range of 10-40 Hz. The tapered appearance of the hyperbola is immediately seen. This implies dip limi- tations and verification is provided by Figure A-3b where a 60 degrees dipping model (cf., Figure 3-5) is shown. As expected, the steeply dipping segment is not imaged well. In addition, the linear artifact observed below the flat segments is generated by the mini-hyperbola phenomena mentioned previously. The mini-hyperbolas from the point 135 diffractors stack in to form the linear artifact, just as the diffractors stack in to image the flat segment. Figure A-4 shows the results of using the same models as in Figure A-3, but with a different filter implementation. Here, the fourth order window cutoff was fixed at kc = 0.8kN. Figure A-4a shows that the hyperbola is of greater extent than in the previous case, implying that large dips may be accomodated. The mini-hyperbolas in thei case are different from those in Figure A-3a. Here, they are much weaker in amplitude than before and they are not as close to one another. Also, the distance from the main hyperbola apex to the the level of the mini-hyperbola peaks is greater. As indicated, large dips can be observed in the 60 degrees dipping model of Figure A-4b. The mini-hyperbolas again contribute to the presence of a linear artifact beneath the flat reflector segments. In addition there is noise form Fourier transform (frequency domain to time domain) wraparound which appears as linear dipping features of low amplitude. Overall, the result is a good one. Finally, I illustrate the effects of using a 390th order window (effectively a box-car) with a fixed cutoff as for Figure A-4. In Figure A-5a an excellent hyperbola with large flanks is seen. The mini-hyperbola artifact is no longer present. With the observa- tions from the previous examples, the suggests that the artifact arises from window implementations that filter out portions of the non-evanescent (i.e., propagating) wave- field. Also, the steeply dipping segment of the 60 degree model is imaged very well despite the presence of some artifacts and noise. Even after 102 extrapolation steps, the computation was stable. The aperture selection in the last two examples was inconsistent with the theoretical expectations and implies that too large or too small an aperture might be chosen. 136 Figure A -3 . Seismograms using a varying cutoff, fourth order windowing strategy. A point diffractor model, (b) A 60 degrees dipping model. 137 138 139 However, in the these cases 1 found remarkably good results; the expected instability of many recusive steps was not observed. In retrospect, I now suggest that higher order windows with fixed cutoff (preferably near the Nyquist wavenumber) schemes should be explored further. My decision for using the low order, varying cutoff was primarily influenced by the consistent ability to predict the wavelet aperture. Perhaps too much importance may have been attached to this, thus, leading to rather conservative results in the algortihm's capability to accomodate steeply dipping features. 140

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
United States 10 0
China 6 2
Germany 3 1
France 2 0
City Views Downloads
Ashburn 7 0
Unknown 6 1
Beijing 5 0
Alexandria 2 0
Shenzhen 1 2

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}

Share

Share to:

Comment

Related Items