COMPLEX SEISMIC DATA: ASPECTS OF PROCESSING, ANALYSIS, AND INVERSION By Timothy Ellis Scheuer B. Sc. (Geophysics) University of Utah M. Sc. (Geophysics) University of British Columbia A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T ' OF T H E REQUIREMENTS FOR T H E DEGREE OF D O C T O R OF PHILOSOPHY in , * T H E FACULTY OF GRADUATE STUDIES GEOPHYSICS AND ASTRONOMY We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA February 1989 Timothy Ellis Scheuer, 1989 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-eGPhs^fcS k HRjCOHOm The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstract A common theme in this thesis is the use of complex signal attributes to facilitate the processing, analysis, and inversion of seismic data. Complex data are formed from real data by removing the negative frequency portion of the Fourier transform and doubling the positive frequency portion. Removing the dual nature of frequency components in real data allows the development of efficient algorithms for time-variant filtering, local phase velocity estimation, and subsequent velocity-depth inversion of shot gathers using Snell traces. For 1-D seismic data, I develop a computationally efficient time-variant filter with an idealized boxcar-like frequency spectrum. The filter can either be zero phase, or it can effect a time-variant phase rotation of the input data. The instantaneous low-frequency and high-frequency cutoffs are independent continuous time functions that are specified by the user. Frequency-domain rolloff characteristics of the filter can be prescribed, but the achieved spectrum depends on the length of the applied filter and the instantaneous frequency cutoffs used. The primary derivation of this theory depends upon the properties of the complex signal and the complex delta function. This formulation is particularly insightful because of the geometrical interpretation it offers in the frequency domain. Basically, a high-pass filter, can be implemented by shifting the Fourier transform of the complex signal towards the negative frequency band, annihilating that portion of the signal that lies to the left of the origin, and then shifting the truncated spectrum back to the right. This geometrical insight permits inference of the mathematical form of a general time-variant band-pass filter. In addition, I show that the time-variant filter reduces to a Hilbert transform filter when the derivation is constrained to include real ii signal input. Application of the procedure to a spectral function permits frequency-variant windowing of an input time signal. For 2-D arid 3-D seismic data, I propose a new method that uses the concepts of complex trace analysis for the automatic estimation of local phase velocity. A complex seismic record is obtained from a real seismic record by extending complex trace analysis into higher dimensions. Phase velocities are estimated from the complex data by finding trajectories of constant phase. In 2-D, phase velocity calculation reduces to a ratio of instantaneous frequency and wavenumber, and thus provides a measure of the dominant plane-wave component at each point in the seismic record. The algorithm is simple to implement and computational requirements are small; this is partly due to a new method for computing instantaneous frequency and wavenumber which greatly simplifies these calculations for 2-D and 3-D complex records. In addition, this approach has the advantage that no a priori velocity input is needed; however, optimum stability is achieved when a limited range of dipping events is considered. Preconditioning the record with an appropriate velocity filter helps reduce the detrimental effects of crossing events, spatial aliasing, and random noise contamination. Accurate recovery of local phase velocity information about underlying seismic events allows the rapid evaluation of seismic attributes such as rms velocity and maximum depth of ray penetration. I utilize local phase velocity data from a shot gather for the estimation and inversion of Snell traces. The primary Snell trace corresponding to a 1-D velocity model locates all primary reflection energy,corresponding to a fixed emergence angle. Constraints on interval velocity and thickness obtained from several estimated Snell trajectories are inverted using SVD to provide a least squares velocity-depth model. The estimation and inversion is efficiently carried out on an interactive workstation utilizing constraints from a hyperbolic velocity analysis. Finally, Snell trace inversion is extended to an inhomogeneous medium. When dips are small, averaging Snell traces of a common phase iii velocity from forward and reversed shot gathers approximately removes the effects of planar dip. This allows recovery of velocity and depth vertically beneath the midpoint of the source locations used to obtain the reversed information. iv Table of Contents Abstract ii Acknowledgement xviii Preface xix 1 Aspects of time-variant filtering 1 1.1 Introduction . 1 1.2 Highpass filter development 5 1.3 Derivation of general time-variant filters 9 1.4 Analysis of Filter Errors 15 1.5 Relation to the Hilbert transform filter 23 1.6 Frequency-variant windowing 28 1.7 Conclusion 31 2 Local Phase Velocity from Complex Seismic Data 35 2.1 Introduction 35 2.2 Instantaneous frequency 37 2.3 Extension to 2-D •. . . . . . : 41 2.4 General formulation 45 2.5 Effects of curvature and noise . 50 2.5.1 Curvature 50 2 . 5 . 2 Noise ; : 5 3 v 2.6 Examples and related attributes . . . 55 2.6.1 Synthetic example . 55 2.6.2 Estimation of rms velocity 57 2.6.3 Maximum depth bounds 60 2.6.4 Field data example 63 2.7 Image Analysis 66 2.8 Conclusion 67 3 Snell Trace Estimation and Inversion 69 3.1 Introduction . . 69 3.2 Snell trace inversion formulae 73 3.2.1 Discrete formulation 73 3.2.2 Continuous formulation . . 74 3.2.3 Example 76 3.2.4 rms velocity 76 3.3 Least squares algorithm 81 3.4 Inversion of Field Data . 84 3.5 Plane dipping layers 93 3.6 Conclusion 103 4 Discussion 105 4.1 On the time-variant filter 105 4.2 On the calculation of local phase velocity 106 Bibliography 109 A Derivation of the complex bandlimited sine function 114 vi Computation of the digital complex delta function Resolution of Local Phase Velocity Approximate traveltime relations for dipping layers vii List of Figures 1.1 (a) Final stacked reflection section with time-variant predictive deconvolu-tion but without time-variant filtering, (b) Stacked reflection section of (a) after application of the continuously time-variant filter described herein. The band-pass varies linearly from 15 to 80 Hz at time zero to 5 to 25 Hz at 6 s 2 1.2 The time and frequency operations of the time-invariant high-pass filter using complex delta signal input. The cutoff frequency is 50 Hz 10 1.3 The time and frequency operations of the time-invariant low-pass filter using complex delta signal input. The cutoff frequency is 50 Hz 12 1.4 The time and frequency operations of the time-invariant band-pass filter using complex delta signal input. The frequency band-pass is from 20 to 70 Hz , 13 1.5 (a) The instantaneous cutoff frequency functions Q; and £lh and the phase rotation function ip. (b) Trace 1 is the input consisting of seven spikes separated by 250 ms using a sample interval of 4 ms. Traces 2-4 are the high-passed, low-passed, and band-passed output in accordance with the functions given in (a). No phase rotation was applied. Trace 5 is the phase-rotated and band-pass filtered output, (c) a time-variant spectral analysis of the band-pass trace 4 in (b). The desired passband in each case is defined by the outer edges of the darkened portion of the computed spectrum 16 V l l l 1.6 (a) The spectrum of the exact windowed filter of length 200 ms. (b) Real spectra of the approximate filters with cutoff gradients of 25, 50, 100 and 500 Hz/s, while fixing the filter truncation length at 200 ms 19 1.7 (a) Four instantaneous low-pass frequency functions containing frequency gradients of (1) 25, (2) 50, (3) 100, and (4) 500 Hz/s. (b) The output of a time-variant low-pass filter applied to three input spikes using each of the frequency functions in (a). Trace 4b) is the exact output resulting from numerical integration of equation (1.1). . .- 20 1.8 (a) Spectra of the exact windowed filters using truncation lengths of 100, 200, 400, and 800 ms. (b) Real spectra of corresponding approximate filters with a fixed frequency gradient of 100 Hz/s 21 1.9 The output of a time-variant low-pass filter applied to three input spikes using frequency function 3 in Figure B-2a. Four filter half-lengths of (a) 100, (b) 200, (c) 400, and (d) 800 ms were used, (e) is the exact output using a half-length of 800 ms . 22 1.10 The time and frequency operations of the time-invariant, low-pass filter using real input data. The input signal in (a) is a 90 Hz sine function and the cutoff frequency is 50 Hz. The final output signal in (i) is aliased by energy wrapped into its spectrum above 50 Hz. . 27 1.11 Two time-frequency pass fairways that can be utilized in a frequency-variant window or a time-variant filter operation. The time-frequency pass in (a) is more precisely implemented using a frequency-variant win-dow, and (b) is more accurately implemented using a time-variant filter. This follows since the time-cutoff gradient in.(a) and the frequency-cutoff gradient in (b) are respectively small. . . 30 • ix 1.12 Time domain operations of a simple frequency-variant window, a) Input time signal, b) frequency-variant shift, c) annihilation of positive time energy, d) frequency-variant shift, e) annihilation of negative time energy, and f) frequency-variant shift restoration 32 1.13 The time-cutoff window utilized to effect the frequency-variant window shown in Figure 1.12 33 2.1 (a) Real part of a complex 10 to 50 Hz sine function, (b) Instantaneous frequency of the complex sine using three different formulas. The straight line at 30 Hz results from equations 2.5 and 2.7, while the oscillating curve results from a finite difference approximation to equation 2.5 given by Yilmaz. The analytic formula for the complex sine function was used in this example. . . 42 2.2 (a) Input synthetic trace. I have applied a time-variant band-pass filter to a reflectivity series to obtain the result shown, (b) The solid curve is the result of using formula (2.5) and the dashed curve is the result of using my two-point formula (2.7) to calculate instantaneous frequency from the complex version of the trace in (a) 43 2.3 (a) A plane wave of frequency 20 Hz and wavenumber 0.005 cyc/m. (b) A 2-D sine function with frequency passband 0 to 40 Hz and wavenumber passband 0 to 0.01 cyc/m. (c) A seismic event with a frequency band range of 0 to 40 Hz dipping at 4000 m/s. Below each section is the corresponding ideal Fourier a> — k domain amplitude response. The theoretical phase velocity of all three events are identical 46 x 2.4 Local phase velocity images calculated (a) from the accurate data shown in Figure 2.3, and (b) from the same data containing 10% uniformly dis-tributed random noise. Dark areas of the image represent phase velocity values outside the bounds given. Note that the color scale is non-linear at both ends. 51 2.5 Comparison of local phase velocity estimators for a hyperbolic reflection event. In (a), 1) is the forward, 2) is the backward, 3) the forward and backward average, and 4) the 5 point estimator, (b) The hyperbolic re-flection trajectory analyzed in (a) 54 2.6 . (a) A synthetic shot gather generated from the velocity-depth model in Figure 2.7. It is a kinematic representation of the primary (compres-sional) reflective events, (b) A field shot gather from Western Canada, (c) Reflective events extracted from the field gather, after statics corrections, by selective velocity filtering. . 56 2.7 (a) Interval velocity versus depth model used to produce the synthetic shot gather in Figure 2.6a. (b) Corresponding rms velocity,- and (c) av-erage velocity mapped from two-way normal incidence time, to the depth coordinate. Actual model quantities are given in Table 2.1 58 2.8 (a) Raw phase velocity image computed from the synthetic shot gather. (b) Phase velocity of significant events only 59 2.9 (a) Local rms velocity and, (b) maximum depth penetration images derived from the local phase velocity image of the synthetic shot gather 62 2.10 (a) Raw phase velocity image computed from the field shot gather, (b) Phase velocity image obtained by analyzing the field record in separate velocity windows. 64 xi 2.11 (a) rms velocity derived from the local phase velocity image corresponding to the reflective events of Figure 2.6c and (b) Maximum depth bounds extracted from the rms velocity image in (a). The white curve in (a) represents a reflection hyperbola defined by the estimated rms velocity at the cursor location 65 3.1 (a) Constant phase velocity ray path geometry for primary reflections in a horizontally layered earth model; with source location S and receiver positions Xi,x2,x3. (b) Primary reflection traveltime curves observed in an end-on shot gather. The Snell trace locates all energy corresponding to a fixed phase velocity (takeoff/emergence angle) 70 3.2 (a) The positive half of the local phase velocity image shown in Figure 2.8b with several Snell trajectories drawn in. (b) Snell traces determined by skipping over reflection events 1,2,3,6,8, and 9 77 3.3 (a) Comparison of the velocity models obtained from inverting the three Snell traces in Figure 3.2a with the true velocity model shown by the dashed curve 79 3.4 Comparison of the velocity models obtained from inverting three Snell. traces in Figure 3.2b with the true velocity model. Reflectors 1,2,3,6,8, and 9 were bypassed in the Snell trace definition resulting in interval rms velocities between the surface and reflector 4, reflectors 5 and 7, and re-flectors 7 and 10 . 82 3.5 Least squares inversion of all Snell traces in Figure 3.2a (solid curve) com-pared with the true velocity model (dashed, curve) 85 xu 3.6 (a) Shallow water marine shot record from offshore eastern Canada, (b) Shot record after application of velocity filter, (c) phase velocity image corresponding to (b) showing phase velocities between 4000 and 5000 m/s only. A radial trace (with vp = 4500 m/s) corresponding to the water bottom reverberation is shown along with an outer bound on the primary Snell trace 87 3.7 (a) Shot gather from offshore western Canada; data window shown is free of water bottom reverberations, (b) Shot gather after application of -a velocity filter to remove contamination from negative phase velocity noise. 88 3.8 On the left; phase velocity image obtained from the raw shot record. On the right; phase velocity image obtained from the velocity filtered shot record 89 3.9 (a) Graphics screen from the Snell trace estimation and inversion program INSTEIN. On the right is a local phase velocity image corresponding to the shot data in the center. Four Snell trajectories were estimated and inverted to obtain the velocity-depth model shown at the bottom right. Snell trace picks are constrained by the hyperbolic velocity spectrum shown at the far left. Modelled traveltime curves are superimposed on the shot data to determine goodness of fit. (b) Example of phase velocity isolation to aid Snell trace picking 91 3.10 Forward and reverse ray path geometry for a plane dipping layered medium when the emergent angle for both raypaths is identically ip 94 3.11 Velocity models obtained by inverting Snell traces from, (a) the forward Snell traces, (b) the reverse Snell traces given in Table 3.5 above. Also, the velocity model (middle solid curve) obtained from the average Snell trace is compared with the true vertical velocity-depth model (dashed curve). 97 xin 3.12 Velocity models obtained from the data given in Table 3.7 99 3.13 (a) Plane dipping layered test model showing some ray paths arising from a split spread experiment, (b) Resulting traveltime curves from the primary reflections only 101 3.14 (a) Velocity-depth models obtained by inverting forward Snell traces at shot locations of 3, 5, and 7km. (b) Velocity-depth models obtained by inverting reversed Snell traces at the same shot locations, (c) Velocity-depth models obtained by inverting averaged Snell traces from forward and reverse shot data. The dashed curves represent the true vertical velocity-depth models 102 3.15 (a) Velocity-depth models obtained by inverting average Snell information using reversed gathers at (1,5 km), (3,7 km), and (5,9 km) from Figure 3.5, where the first and second distances are the forward and reversed shot locations respectively, (b) Inverted results when the corresponding split-spread Snell traces are included in the Snell trace average. In (a) and (b), the dashed curves represent the true vertical velocity-depth models for the midpoint locations of the shotpoint pairs. 102 B . l The digital complex delta function sampled at 4 ms and truncated at 200 ms . 117 B.2 (a) The spectral response of the digital complex delta function when the exponent of the cosine window is 0, 2, and 8, for a truncation length of 200 ms. (b) The spectral response of the digital complex delta function with truncation lengths of 100, 200, and 800 ms with the window exponent fixed at a = 2. The spectrum is shown around 0 Hz to emphasize the rolloff characteristics; these characteristics also occur at the Nyquist frequency. 119 xiv D. l Geometrical relationships of forward and reverse ray paths used for the development of an approximate average traveltime relation 123 xv List of Tables 2.1 Model parameters used to produce the synthetic shot gather shown in Figure 2.6a. Parameters are listed for each layer with bottom reflecting boundary "R". Corresponding normal-incidence rms and average veloci-ties for each reflector are also listed. 60 2.2 Comparisons of rms velocities (m/s) and depth bounds (m) obtained from the images in Figure 2.9 with the true normal-incidence rms velocities and true depths of the model in Figure 2.7. The measured values repre-sent mean values taken along the traveltime curve of each reflector. "R" indicates reflector number 61 3.1 Snell trace parameters of horizontal offset x (rounded to the nearest meter) and traveltime t (rounded to the nearest milli-second) for three trajectories shown in Figure 3.2a ("R" denotes the reflector number). The inverted ve-locities v (m/s) and thicknesses Az (m) are listed for each layer. Compare these with the true model parameters given in Table 2.1 78 3.2 Snell trace parameters of horizontal offset x (m) and traveltime t (s) for three Snell trajectories shown in Figure 3.2b. Reflectors 1,2,3,6,8, and 9 were not included (skipped) in the Snell trace picks. The inverted veloc-ities v (m/s) are compared with theoretical angle dependent interval rms velocities v^, (m/s) calculated by numerical evaluation of equation 3.15. 81 xv i 3.3 True model parameters on the left. On the right are the velocity (m/s) and depth (m) estimates obtained by least squares inversion of all Snell traces shown in Figure 3.2a 85 3.4 Velocity v and depth z obtained from a least squares inversion of the four Snell traces shown in Figure 3.9a. 90 3.5 Snell trace parameters of offset and traveltime computed by ray tracing forward and reversed rays of the plane dipping model shown in Figure 3.10. The forward and reversed averages are compared to the Snell trace~ parameters of an equivalent horizontal model. The common emergent angle for these results is if; = 20°. The model parameters of vertical depth and dip refer to the model shown in Figure 3.10 97 3.6 Ray traced vertical angles for each reflector shown in Figure 3.10. Averag-ing the cosines of these angles for each boundary results in the.angles 0j, that are compared to the vertical angles obtained from a horizontal model with the vertical depths given by the same model 98 3.7 Interval velocity and thickness estimates from averaged Snell traces when the angles (£) of the dipping layers shown in Figure 3.10 increase. Velocity and depth errors increase with dip 98 3.8 Interval velocity and thickness estimates from averaged Snell traces when the emergence angles (ip) increase. Velocity and depth errors increase with emergence angle 99 X V l l Acknowledgement This thesis is dedicated to my wife Jan. It would have been difficult with out you Jan. I thank the members and faculty of the Department of Geophysics and Astronomy for providing a creative environment to do research. Committee members Garry Clarke, Ron Clowes, and Matt Yedlin deserve special thanks. I appreciate ^ .the assistance of fellow graduate students and research assistants, David Aldridge, John Amor, Andy Calvert, Rob Ellis, David Lumley, and Colin Zelt. Part of this research was inspired by fellow colleagues at Amoco Production Co., Ken Peacock and Richard Crider. I'd like to thank David Morris, Bob Pangle, and Bill Damewood at Amoco for furnishing field seismic data. I am also grateful to Shell Canada Frontier Exploration group for contributing field seismic data. I know that providing proper field data tapes is no easy task. This work was supported by NSERC grant 5-84270. I am grateful for that support. A very special thanks is reserved for my one and only graduate supervisor, Dr. Dou-glas W. Oldenburg. If I have contributed the skeleton, then Doug has added the muscle. xvm Preface Each chapter in this thesis is a self contained study. Indeed, Chapters 1 and 2 are expansions of two separate published papers (Scheuer and Oldenburg, 1988 [35],[36]) and a precis of Chapter 3 will be submitted for publication as well. Since each chapter contains a separate introduction and conclusion, I feel overall introduction and conclusion chapters would be superfluous additions to the thesis. Nevertheless, my research is bound by a common theme; complex data analysis. I therefore provide the following prefatory remarks. Four years experience in the seismic data processing industry convinced me that another look at some of the fundamental processing concepts was warranted. Band-pass filtering is one of the oldest and most useful of signal processing techniques. Nevertheless, I saw a need for new ideas concerning an efficient time-variant filter. Therefore, I returned to UBC in 1985 with an idea for time-variant filtering that stemmed from my exposure to the Hilbert transform filter of Cain 1972. It was well known that the Hilbert transform filter causes aliasing if the input signal contains energy in frequencies above half the Nyquist. This is due to amplitude modulations that cause wrap around effects or aliasing in digital signals. I eventually discovered that by using complex signals in the time-variant convolution integral, aliasing could be completely avoided. This discovery, partly inspired by the work of Pann and Shin 1976, led to the research comprising the body of this thesis. The theory of complex signals is well documented by such authors as Gabor (1946), Oppenheim and Schafer (1975), Bracewell (1978), and Taner et al (1979). It is not my purpose here to reproduce the theoretical development in those publications but to build upon the basic framework provided. Quite simply, complex data are formed from real xix data by removing the negative frequency portion of the Fourier transform and doubling the positive frequency portion. Removing the duality of frequency components from real data allows the development of efficient algorithms for time-variant filtering, local phase velocity estimation from shot gathers, and subsequent velocity-depth inversion using Snell traces. The processing part of the thesis is embodied in Chapter 1; time-variant filtering. Formulation of a time-variant filter using complex signals is particularly insightful because of the geometrical interpretation it offers in the frequency domain. The high-pass filter operations begin by separating the pass band and the reject band of an input complex signal into positive and negative frequency bands. Annihilating the reject band and then restoring the pass band to its original spectral location completes the filter operations. This insight allows inference of a general time-variant band-pass filter. Chapter 2 covers the analysis aspects of the thesis; estimation of seismic dip at-tributes. Extension of complex signal analysis into higher dimensions is straight forward if a frequency domain viewpoint is adopted. Multi-dimensional complex data can be formed by removing the negative frequency portion of the Fourier transform of the input real seismic data. In one dimension, instantaneous frequency can be thought of as the frequency of the best fitting complex harmonic at each point of a complex signal. This provides insight into the generalization of a corresponding higher dimensional attribute; local phase velocity. Local phase velocity estimation can be viewed as finding the best fitting complex harmonic plane wave at each point of a complex seismic record. Elements of Chapter 2 provide basic motivation for the contents of Chapter 3 and, in addition, scratch the surface of possible research into the utilization of multi-dimensional complex data attributes in many fields of science. Chapter 3 covers the estimation and inversion of Snell traces from shot gathers. Once accurate local phase velocity of underlying seismic events are readily available, then xx ideas for the extraction of many other subsurface attributes and parameters come to light. Some ideas are listed in the introduction to Chapter 2 . I chose to concentrate my research efforts on the estimation and inversion of Snell traces from the local phase velocity image of a shot gather. Velocity and thickness information extracted from a Snell trace is very important in the interpretation of seismic data. During the process of restructuring theory into a practical algorithm, it became evi-dent that a local phase velocity image from seismic data provides a massive amount of information to interpret. An expert can extract a proportionate amount of useful seis-mic parameters, but a novice can easily get buried in coloured pixels. The workstation program INSTEIN, presented in Chapter 3, is a result of my attempt to make Snell trace estimation and inversion more efficient for the expert and more accessible to the novice. x x i Chapter 1 Aspects of time-variant filtering 1.1 Introduction The frequency and phase character of seismic reflection data varies with time "mostly due to local variation in the character of the subsurface reflectivity, but also due to the effects of absorption and scattering of acoustic energy. Absorption effects, for instance, lead to a systematic loss of high frequency content and an associated change in phase of the propagating wavelet. Thus, any attempt to recover this attenuated frequency information using adaptive or time-variant deconvolution (Clarke,1969, Griffiths, et.al., 1977), or time-variant spectral whitening (Bickel and Natarajan, 1985), may enhance high frequency noise at larger recording times. This is why it is sometimes necessary, as a final processing stage, to use a time-variant band-pass filter on reflection seismograms to uncover the clearest possible image of stratigraphic boundaries. For example, Figure 1.1a represents a seismic reflection sec-tion in which the processing sequence included time-variant deconvolution. In this case, the enhancement of attenuated high frequencies at larger recording times has decreased reflector coherency because the signal at high frequencies is progressively lost in back-ground noise. Generally, the geophysicist would prefer to make a geological interpretation using the reflection section shown in Figure 1.1b. That section has been processed using a time-variant band-pass filter. 1 Chapter 1. Aspects of time-variant filtering 2 10 20 30 40 10 20 30 40 _ ( b ) Figure 1.1: (a) Final stacked reflection section with time-variant predictive deconvolution but without time-variant filtering, (b) Stacked reflection section of (a) after application of the continuously time-variant filter described herein. The band-pass varies linearly, from 15 to 80 Hz at time zero to 5 to 25 Hz at 6 s. Chapter 1. Aspects of time-variant filtering 3 The seismic section shown in Figure L i b was processed using a continuously time-variant band-pass filter of the type discussed in this thesis. In contrast, the conventional method of time-variant filtering applies a sequence of time-invariant filters to overlapping time gates that are then meshed together using a linear taper in the overlap zones. This approach has two disadvantages. The output signal is distorted in the overlap zones (Pann and Shin, 1976) and the gate-wise stationarity assumption used to justify time-invariant filtering is suboptimal (Clarke, 1969, and Ziolkowski and Fokkema, 1986). This type of filter does not mimic the way nature continuously attenuates higher frequencies as travel time increases. A time-variable band-pass filter in which the high and low frequency cutoffs are allowed to vary continuously with time is thus a desirable alternative if such a filter can meet requirements of accuracy and computational efficiency. Despite the practical implications of a high quality digital time-variant filter, there has been very little written in the geophysical literature on this subject. Nikolic (1975), and. Stein and Bartley (1983) describe recursive methods for implementing a continu-ously time-variant filter. In a recursive or autoregressive moving-average (ARMA) filter, each output value is a weighted sum of input values plus a weighted sum of previously computed output values (Shanks, 1967). The filter weights are a set of ARMA coeffi-cients that depend upon the desired output frequency characteristics. The difficulties with recursive filtering are twofold. If a low order ARMA filter is utilized, its amplitude spectrum may deviate significantly from the desired spectrum, either because of over-shoots within the pass-band or unacceptable rolloff characteristics. This deficit can be improved by designing higher order filters as shown by Stein and Bartley (1983) but the price to be paid is that of additional computation expense. The second drawback for recursive filters is that the input must be recursively filtered in a forward and a reverse time direction to obtain a zero phase response. Shanks (1967) investigated the utility and accuracy of the time-invariant recursive Chapter 1. Aspects of time-variant filtering 4 filter. However, little has been written on error propagation when the ARMA coefficients change with time. The results of Stein and Bartley (1983) indicate that there may be difficulty in obtaining a zero phase response. It is perhaps for these reasons that the continuously time-variant recursive filter has not gained wide acceptance as a replacement for the gate-wise time-invariant filter. Pann and Shin (1976) describe a class of time-variant filters wherein the frequency spectrum of a given filter is shifted, without change in shape, along the frequency axis as a function of time. The design is based on an initial filter with the desired spectral shape, and the time-variant result is computed using an approximate solution to the time-variant convolution integral. This approach is limited in its applications to seismic data processing because the filter bandwidth does not change with time. Nevertheless, the formulation of Pann and Shin suggests a general method by which time-variant filtering can be carried out, and that is the thrust of this work. In this chapter, I develop a computationally efficient time-variant filter with idealized band-pass characteristics; that is, it has a time-variant boxcar-like frequency spectrum. The filter can either be zero phase, or it can effect a time-variant phase rotation to the input data. The instantaneous low and high frequency cutoffs are independent and continuous time functions that can be specified by the user. Rolloff characteristics of the filter can also be prescribed although the achieved spectrum depends on the length of the applied filter and the instantaneous frequency cutoffs used. The primary derivation of the time-variant filter depends upon the spectral properties of the complex signal and also the complex delta function. This formulation is particularly insightful because of the geometrical interpretation it offers in the frequency domain. Basically a high-pass filter can be implemented by shifting the Fourier transform of the complex signal toward the negative frequency band, annihilating that portion of the signal that lies to the left of the origin, and then shifting the truncated spectrum back to Chapter 1. Aspects of time-variant filtering 5 the right. This geometrical insight allows inference of the mathematical form of a general time-variant band-pass filter. In addition, I show that my time-variant filter reduces to the Hilbert transform filter of Cain (1972) and Peacock (1985) when the derivation is constrained to include real signal input. 1.2 Highpass filter development The complex seismic trace (Taner et al., 1979) is defined as s(t) = s(t) — iH[s(t)], where s(t) is the real seismic trace, and H[-} is the Hilbert transform operation (Bracewell, 1978,p. 267). A complex signal of this type contains energy in only positive frequencies and is often called an analytical signal. The general form of the time-variant convolution integral for complex signals is. f(t) = Jj(T,t)s(t-T)dT, (1.1) where f(r,t) is a time-variant complex filter that operates on a complex signal s(t). The computational effort implied by equation (1.1) for a band limiting filter defined in r at each instant of output t is formidable in most cases. However, this integral equation can be greatly simplified if the filter used at each instant has a spectrum that is flat in the pass band and zero in the stop band, i.e., if it is an idealized band-pass filter. I first derive an expression for a time-variant high-pass filter. Geometrical insight gained from this result then allows formulation of expressions for general time-variant band-pass filters. The characteristics of the frequency limiting filter are most easily described in the frequency domain. The Fourier transform of the complex filter is F(u,t) = / f(T,t)e—dr, = .4(a;,i)eM w , t ), (1-2) Chapter 1. Aspects of time-variant filtering 6 where A(u>,t) and tp(u,t) are the time-varying amplitude and phase spectrum respec-tively. The inverse Fourier transform is given by 1 t+°° « f{r,t) = —J_^-F{u,t)e^du;. (1.3) The amplitude spectrum of the time-variant high-pass filter is represented by a Heaviside step function that shifts in time according to a specified instantaneous low-cut frequency. That is, A(u>,t) = u(u - Q(t)), ^ (1.4) where Q.(t) is the instantaneous cutoff frequency function that regulates the time-variable shift. Inverse Fourier transformation of A(w,t) produces a zero phase filter. Using the shift theorem, the resulting high-pass filter is /(r,0 = e i n ( t ) ^(r) /2 , (1.5) where S(T) = S(T) + ijixr shall be referred to as the complex delta function. Equation (1.5) shows that the time-variant high-pass filter consists of the complex delta function modulated by a time-variant phase that depends on the instantaneous cutoff frequency function. This filter is highly localised in r; because of the compact nature of the complex delta function, most of the energy of f(r,t) is concentrated around r = 0. The filter /(r, t) is infinitely long and hence must be truncated in r for digital imple-mentation. A finite impulse response filter /(r, t) can be represented as the product of the desired infinite impulse response filter and a finite duration window w(r) i.e., f(T}t) = f(T,t)-Ar). (1-6) From the convolution theorem, this is equivalent to convolving the Fourier transform of the filter with a shaping operator i.e., F(u>,t)=u(u-Q(t))*U(w), (1.7) Chapter 1. Aspects of time-variant filtering 7 where the shaping operator W(u>) is the Fourier transform of w(r). Equations (1.6) and (1.7) show that a high-pass filter having finite length and predetermined spectral char-acteristics can be designed by suitably choosing w(£) or W(u>). Nevertheless, I continue deriving the time-variant filter utilizing the unwindowed equation (1.5) and defer further discussion of the filter spectral response until after important aspects of the filter which are independent of the truncation window have been elucidated. The compact nature of the filter leads to further simplification. First define 0(0 = ^ , -(1-8) where Q{t) is the instantaneous phase corresponding to the instantaneous cutoff frequency Cl(t). My goal is to write e , n ^ T in terms of {t — r) in order to simplify the time-variant integral equation (1.1). Using a Taylor expansion gives . /i/ \ dB(t) r2d29(t) , . ,(«- T) = , ( 0 - r - £ i + T - s i 2 - . . . . (1.9) Rearranging and using equation (1.8) produces n f , 6{t)-9(t-r) rdSlit) If the term in dQ,(t)/dt and higher orders can be neglected, then . . 6(t)-9(t-r) f 2 ( i ) ~ - ^ ^ '-. (1.11) T The accuracy of this approximation will depend upon the upper limit of the prescribed filter length (i.e., the maximum value of r) and also upon the derivatives of Q,(t). If Q(t) is highly oscillatory or has jump discontinuities, this approximation may be invalid at nonzero values of r. For most seismological problems however, Q(t) is a slowly varying function and hence the approximation will be valid for reasonably large values of r. Substitution of equation (1.11) into (1.5) reduces the high-pass filter to f(r,t) ~ ei[6^-8^S{r)/2. (1.12) Chapter 1. Aspects of time-variant filtering 8 Inaccuracy of approximation (1.11) at larger values of r is offset by the compact nature of the complex delta function. That is, the filter decays rapidly to zero as the approximation is likely to worsen. This property allows the use of most practical instantaneous frequency functions without loss of accuracy. Nevertheless, the validity of equation (1.12) is explored within the context of filter accuracy and the output of a time-variant filter in Section 1.4. Substituting equation (1.12) into equation (1.1) and simplifying gives the following formula for time-variant high-pass filtering, r(t) = eie^[S(t)/2*s(t)e-ie^}. "(1.13) The real part of this complex output is the filtered trace. Equation (1.13) is an efficient simplification of the time-variant convolution integral, since the convolution is now time — invariant. Time variance is transformed to multiplication in time by complex exponential factors dependent on the instantaneous frequency cutoff function. Considerable insight regarding the operation of this high-pass filter can be obtained by applying equation (1.13) to an example in which Q(t) is constant for all time. The time and frequency domain effects of that filter are enlightening. Setting Q(t) — u>c gives 6{t) = u>ct, and thus equation (1.13) becomes f(t) = eiUct[S(t)/2*i(t)e-iu't]. (1.14) The Fourier transform of this equation gives, R(u) = S(ui — u>c) * U(UJ)S(UJ + uic). (1-15) Figure 1.2 shows graphical representations of equations (1.14) and (1.15) when the input signal is a complex delta function. The actual computations for this example shown on the left were done in the time domain using equation (1.14) and a cutoff frequency of 50 Hz or u)c = 27r50. On the right are the resulting frequency domain effects described by Chapter 1. Aspects of time-variant filtering 9 equation (1.15). These figures were produced by taking the Fourier transform of the time domain results. In Figure 1.2a, the real and imaginary parts of the input digital complex delta function are shown on top and bottom respectively. This complex digital signal has a real periodic Fourier transform with energy in the positive frequency spectrum out to the Nyquist frequency of 125 Hz as shown in Figure 1.2b. Application of the negative phase shift in time generates the signals in Figure 1.2c, while Figure 1.2d shows the corresponding frequency domain effect of shifting the entire spectrum to the left by 50 Hz. Because energy has been shifted into the negative frequency spectrum, the resulting signal is no longer defined complex analytic. Figure 1.2e shows the results of convolution with the complex delta function shown in Figure 1.2a. This annihilates the negative frequency spectrum as shown in Figure 1.2f. Finally, application of the positive phase shift in time restores the frequency spectrum to its original location as shown in Figures 1.2g and 1.2h. The resulting complex signal in Figure 1.2g is a high-passed version of the complex delta function. Rolloff characteristics of the spectra shown in Figure 1.2 are due, in part, to the cosine window used to truncate the complex delta function. Appendix B covers this aspect in more detail. 1.3 Derivation of general time-variant filters The high-pass filter development is easily generalized to produce expressions for low-pass and band-pass filters. The difference between high-pass and low-pass filtering is whether the negative or positive frequency spectrum is annihilated after the negative phase shift has been applied. Observe in Figures 1.2c and 1.2d that the negative phase shift serves to separate the low and high frequency bands of the original signal into the negative and positive frequency spectra respectively. The separation is centered on the Chapter 1. Aspects of time-variant filtering 10 s(t) Re 1 Aw. 1_J I I \ J I I L —0.20 0.0 0.20 a SjCtJ^sCtJexpl^-nSOt] J • I I I I I I 1 L -0.20 0.0 0.20 c s2(t)=sf(t)*6(t)/2 i i i i i i i -0.20 0.0 0.20 e r( t)=s2( t Jexp[i2-n50t] U I I I ' ' ' ' ' •0.20 0.0 0.20 time (sec) g 2.0 0.0 -\ \ 1 1 1 1 100 0 b 100 S1(U)=S(U+2TT50) 2.0 0.0 J L J I / 00 0 d 100 S2(CJ)=u(o))S 1 (o) 2.0 0.0 I ' I ' L -100 0 f 100 R(u)=S2(a-2-n50) 2.0 0.0 z^ J I I L 100 0 100 frequency (Hz) h Figure 1.2: The time and frequency operations of the time-invariant high-pass filter using complex delta signal input. The cutoff frequency is 50 Hz. Chapter 1. Aspects of time-variant filtering 11 instantaneous frequency cut-off of 50 Hz which has been shifted to zero frequency. Instead of zeroing the negative frequency spectrum at this point, the positive frequency spectrum may be annihilated by application of the conjugate complex delta function £~(t). Then the subsequent positive phase shift returns the desired low-pass signal. Therefore, by deduction, the formula for time-variant low-pass filtering can be written as p(t) = eie^[S*{t)/2 * s(t)e-ie^}. (1.16) Figure 1.3 shows the time and frequency operations of this filter in time-invariant mode. The above time-variant low and high-pass filters allow complete flexibility in applying zero phase time-variant band-pass filters to complex signals. For example, a time-variant band-pass output may be obtained by computing the difference between the output of two time-variant low/high-pass filters. But, time-variant band-pass filtering may also be accomplished by cascading the low and high-pass filter equations. From the insight gained in the previous derivations, a time-variant band-pass filter that will pass frequencies above Oj(i) and below fl/i(£) can be written explicitly as b{t) = eie'{8/2*eiAe[8'/2*se-ieh}}, (1.17) where, 0,(0= /'flitfK. °h(t)= f u n i m , Jo Jo AO = Oh —and the functional dependence of the right-hand side on t has been dropped for convenience. Figure 1.4 shows the time and frequency operations of the time-invariant mode of this filter. The previous filters were designed to leave the phase of the input signal unaltered. Yet time-variant phase alterations of the propagating seismic wavelet are often evident in final processed sections when compared with synthetic seismograms. Such phase alterations are implied by the physics of wave propagation in a lossy medium and, in addition, Chapter 1. Aspects of time-variant filtering 12 s(t) 1.0 Q h i i i \ i j i L — 0 . 2 0 0.0 0 .20 S(w) 2.0 1.0 0.0 : b \ 1 1 1 1 - 1 0 0 0 100 s 1(t)=s(t)exp[-i27T50t] - C 1.0 0.0 • I I I I l _ J I L -0.20 0.0 0 .20 S 1 ( o j ) = S ( c j + 2 n 5 0 ) 2.0 1.0 0.0 Mill \ L I - 1 0 0 0 100 s 2(t) = S l ( t ) * 6 * ( t ) / 2 1.0 0.0 - e • i i i i i i i - 0 . 2 0 0.0 0 .20 S 2 ( w ) = u ( - a > ) S 1 ( a ) ) Mill 1 1 - 1 0 0 0 100 p(t) = s 2 ( t)exp[i27T50t] 1.0 0.0 g - 0 . 2 0 0.0 0 .20 t ime (sec) P(w) = S 2 ( w - 2 7 T 5 0 ) 2.0 1.0 I II II ' 1 1 1 1 ' - 1 0 0 0 100 f requency (Hz) Figure 1.3: The time and frequency operations of the time-invariant low-pass filter using complex delta signal input. The cutoff frequency is 50 Hz! Chapter 1. Aspects of time-variant filtering 13 s ( t ) 1.0 0.0 -0.20 0.0 0.20 2.0 1.0 0.0 l b \ •100 100 \ s1(t)=s(t)exp[-i2n70t] 1.0 0.0 _ 1 I I I I I L —0.20 0.0 0.20 S , ( W ) = S ( W - I - 2 7 T 7 0 ) 2.0 1.0 0.0 - d : \ 1 1 1 1 -100 100 s2(t)=s1(t)*<5«(t)/2 - e 1.0 0.0 J L • I I I I I I l _ L -0.20 0.0 0.20 S 2 ( C J ) = U ( - C J ) S 1 ( C J ) 2.0 1.0 0.0 ' • ' i i_ -100 100 s 3(t)=s 2(t)exp[i27T50t] 1.0 0.0 -'\/\.— _l I I 1 I I I L -0.20 0.0 0.20 S3(co)=S2(w)*5(w-2n50) 2.0 . 1-0 0.0 : h •100 100 s 4(t)=s 3(t)»<5(t)/2 1.0 0.0 I I I I I L -0.20 0.0 0.20 S 4 (CO)=U(W)S 3 (CJ ) 2.0 1.0 0.0 - j ' 1 1 1 1 -100 100 b(t)=s4(t)exp[i2n20t] - k 1.0 0.0 **^^— — — 3 I I I I I I L _ l -0.20 0.0 0.20 time (sec) B(w)=S4(«)*<5(w - 2TT20) 2.0 1.0 0.0 h _i i i i _ i _ -100 0 100 frequency (Hz) Figure 1.4: The time and frequency operations of the time-invariant band-pass filter using complex delta signal input. The frequency band-pass is from 20 to 70 Hz. Chapter 1. Aspects of time-variant filtering 14 could be caused by incorrect image enhancement procedures applied in the processing stage. Exact characterization of these phase distortions is difficult at best. However, if the resultant phase shift can adequately be approximated as being independent of frequency, then the time-variant filter presented here can be generalized to correct for such shifts. This could prove useful in reconciling processed seismic data with a synthetic trace derived from a well log. Time-variant, but frequency-independent phase rotation may be easily included in the derivation of the time-variant complex filter. For example, consider the frequency spectrum of an all-pass complex filter with time-variable phase rotation F(u>,t) =u{w)ei^t\ (1.18) where. tp(t) is the frequency invariant instantaneous phase spectrum of the filter. The inverse Fourier transform of this expression gives t f(T,t)=eim(T)/2. (1.19) Substituting this into the time-variant convolution integral and simplifying gives q(t) = e^ ( t )i(t). (1.20) The phase rotated output is thus obtained by applying a time-variant phase shift to the input complex signal. Combining time-variant band-pass filtering with time-variant phase rotation leads to the following general relation for most basic time-variant filtering applications • g(t) = ei^+^[6/2*eiAe[S'/2*se-i^]\. (1.21) A computer processing algorithm based on this equation can accommodate all of the foregoing applications of time-variant filtering. Digital implementation of this filter in-volves three complex vector multiplications and three complex convolutions (the third Chapter 1. Aspects of time-variant filtering 15 being s = s •* 6, where s is the real input trace). However, the most efficient way to implement this equation is to perform the convolutions in the frequency domain. Since each convolution involves a form of the complex delta function, only one transfer function needs to be designed for all three frequency domain multiplications. This can reduce the computation time by a factor of three and, thus, the application of this continuously time-variant filter may, in many cases, be more efficient than a gate-wise time-invariant filtering approach (when many gates are required). Figure 1.5 shows examples of the general time-variant filter when the input ^ signal consists of seven unit impulses separated by 250ms and sampled at 4ms (trace 1). Figure 1.5a shows the instantaneous frequency functions used. The band-pass fairway varies linearly from 40 and 100 Hz at t = 0, to 5 and 40 Hz at t = 2s. Time-variant phase rotation varies linearly from 0° at t = .25s, to 180° at t = 1.75s and is shown in Figure 1.5a as the dashed line. The real outputs of a high-pass, a low-pass, a band-pass, and a time-variant phase rotation filter are shown in Figure 1.5b. The performance of the band-pass filter is demonstrated in Figure 1.5c by the am-plitude spectrum of each band-pass response of trace 4 shown in Figure 1.5b. Each spectrum was computed using a 280ms window of data centered on the corresponding response peak. The Fourier magnitude was computed after smoothing the data using a cosine window with an exponent of 1 [see equation (B.5)]. The achieved instantaneous band-pass agrees very well with the desired band-pass shown by the darkened portion of each spectrum. 1.4 Analysis of Filter Errors A primary assumption made in the filter derivation was that higher order terms in equa-tion (1.10) can be neglected. The advantage of that approximation is that it permits Chapter 1. Aspects of time-variant filtering 16 0.0 1.0 time (sec) 2.0 0 . 5 t . O 1 . 5 t i m e ( s e c ) 0.0 0.5 1.5 2.0 1.0 t i m e (sec) Figure 1.5: (a) The instantaneous cutoff frequency functions and Q,H and the phase rotation function I/J'.' (b) Trace 1 is the input consisting of seven spikes separated by 250 ms using a sample interval of 4 ms. Traces 2-4 are the high-passed, low-passed, and band-passed output in accordance with the functions given in (a). No phase rotation was applied. Trace 5 is the phase-rotated and band-pass filtered output, (c) a time-variant spectral analysis of the band-pass trace 4 in (b). The desired passband in each case is defined by the outer edges of the darkened portion of the computed spectrum. Chapter 1. Aspects of time-variant filtering 17 the filter equations to be written as a time-invariant convolution. The disadvantage is that the approximation introduces errors such that the achieved filter is not precisely equal to the desired filter. The discrepancy between the achieved and desired filters is dependent upon the truncation window (as described in Appendix B), the filter length, and the derivatives of the instantaneous frequency cutoff function. The purpose of this section is to assess that discrepancy. It is not feasible to present an error analysis that is completely general because each filter utilizes different functions of Q,h(t) and fij(i). I therefore present an analysis of a time-variant low-pass filter in which the instantaneous frequency cutoff varies linearly with time. Although simple, this analysis provides considerable insight into the close-ness with which more complicated filters achieve their desired goal. It also provides a formalism whereby the user can make a rigorous error analysis and therefore alter the filter parameters to produce optimum results for each particular application. From the insight provided by the high-pass filter development, the low-pass filter can be written /(r ,0 = e i n ^ * ( r ) / 2 . (1.22) Including the truncation window w(r) that provides a desired rolloff for the spectrum I obtain /(r ,0 = e i n ^ ( r ) . w ( r ) / 2 . (1.23) This windowed filter will be used as the starting point in the error analysis; that is, differences in our approximate filter from the exact windowed filter will constitute the error. Both character and performance differences are considered in the error analysis. Using equation (1.10), I expand fi(£)r in terms of the instantaneous phase 9(t) and the higher order derivatives of Q,(t). That is, „ , s „ / x r2dQ(t) r3d2ri(t) m r = « ( 0 - « « - r ) + T ^ - 3 r - S H + - , . Chapter 1. Aspects of time-variant filtering 18 = 6{t)-6(t-T) + e(T,t). (1.24) The approximate time-variant low-pass filter is then, fa(r,t) = ei^-6^S'(T)^(r)/2 = e W ) ^ ^ ) ] ^ ( r ) . w ( r ) / 2 . (1.25) The approximate filter can be used with confidence if the error term e(r, t) is small. However, when e(r, t) becomes large, filter error must be quantified before using fa(T,t) in a time-variant filtering algorithm. In general, the filter error will depend on the magnitude and shape of fi(t) and on r as indicated by the higher order terms in equation (1.24). Oscillatory functions of Q(r) with large gradients increase the error term. And there is also a r variable error to consider which worsens as the truncation length of the filter increases. I next provide an error analysis when.Q(t) is linear. This analysis will suffice for many cases of practical interest in which Q,(t) is specified as a piecewise linear function. Choose Q(t) = 27T7i,.where 7 is the frequency cutoff gradient in Hz/s. The associated instantaneous phase is 9(t) = 7r7£2, and the error term reduces to first order, e(r, t) = TVJT2.. Because the error term does not depend on t, the effect of errors at t = 0 can be investigated without loss of generality. In this case, the exact and approximate filters reduce to /(r,0)=i-(r).w(r)/2 • (1.26) and /.(riOJ = e - ^ ^ ( r ) . i ( r ) / 2 (1.27) respectively. It is thus necessary to investigate the errors due to the exponential term in the approximate filter. I do this by examining the difference between the achieved spectral response given by the Fourier transform of (1.27) and the exact spectral response given by Chapter 1. Aspects of time-variant filtering 19 ( a ) \ \ (b) \ 7=500Hz/s -40 -20 0 20 40 frequency (Hz) Figure 1.6: (a) The spectrum of the exact windowed filter of length 200 ms. (b) Real spectra of the approximate niters with cutoff gradients of 25, 50, 100 and 500 Hz/s, while fixing the filter truncation length at 200 ms. the Fourier transform of (1.26). The error will in general depend on the filter parameters 7, T, and a. In this analysis, I fix a = 2 and observe the effects of varying 7 and T. In addition to the spectral differences, I show the corresponding effects on the output of a time-variant low-pass filter when these parameters are varied. Figure 1.6b shows the real spectrum. of the approximate filter when the frequency cutoff gradient varies from 25 to 500 Hz/s while fixing T=200 ms. The spectrum of the exact windowed filter is shown in Figure 1.6a for comparison. Note that the higher gradients produce more ringing in the filter spectrum. Figure 1.7b shows the output of the approximate time-variant low-pass filter using each of the instantaneous cutoff frequency functions shown in Figure 1.7a, where the in-put consists of three unit spikes at .25, .5, and .755. Trace numbers indicate the frequency function used. For comparison, trace 4b) was calculated using the exact windowed filter Chapter 1. Aspects of time-variant filtering 20 120 100 N (J c <u D CT 4b) i 0.0 0.4 0 .8 t i m e (sec) 0 .0 0 . 2 0 .4 0 .6 t i m e ( s e c ) 0 . 8 1.0 Figure 1.7: (a), Four instantaneous low-pass frequency functions containing frequency gradients of (1) 25, (2) 50, (3) 100, and (4) 500 Hz/s. (b) The output of a time-variant low-pass filter applied to three input spikes using each of the frequency functions in (a). . Trace 4b) is the exact output resulting from numerical integration of equation (1.1). and a numerical integration of equation (1.1). Traces 4a) and 4b) compare, the approxi-mate and exact impulse responses for the extreme 500 Hz/s case. Observing the middle impulse response on both traces, it is evident that the approximate result contains some leakage of undesired frequency information on both sides of the center lobe. That is, the approximate result is lower frequency to the left and higher frequency to the right of the center lobe than is the exact result. Figure 1.8b shows the real spectrum of the approximate filter when the one-sided filter length varies from 100 to 800 ms while fixing 7=100 Hz/s. The spectral response becomes more oscillatory as the filter length increases, whereas, the spectrum of the exact windowed filter becomes sharpened as this length increases as shown in Figure 1.8a. This is due to the r variable error which offsets the presumed advantage of using a longer filter. When large cutoff frequency gradients are desired, a short filter may be advised. Figure 1.9 shows the output of the approximate time-variant low-pass filter using Chapter 1. Aspects of time-variant filtering 21 f requency (Hz) Figure 1.8: (a) Spectra of the exact windowed filters using truncation lengths of 100, 200, 400, and 800 ms. (b) Real spectra of corresponding approximate filters with a fixed frequency gradient of 100 Hz/s. filter half-lengths of a) 100, b) 200, c) 400, and d) 800 ms. The cutoff frequency function used is number 3) in Figure 1.7a (at 100 Hz/s). For comparison, trace e) was calculated using the. exact windowed filter of half-length 800 ms and a numerical integration of equation (1.1). The approximate result in trace d) shows a more oscillatory response than the exact result in trace e), and slight leakage of undesirable frequency information is evident in the approximate result. For linear f2(i), errors in the approximate filter depend on the frequency gradient and on the truncation window. Higher frequency gradients require a shorter filter to control spectral ringing in the filter.response. It• is reasonable to expect that these observations also hold approximately for nonlinear functions of Q(t). However, I suggest the following method of determining optimum filter parameters when using nonlinear instantaneous frequency functions Cli(t) and Q/j(£) to perform'time-variant band-pass filtering. a) Select an initial set of filter parameters (eg., window type and length). Chapter 1. Aspects of time-variant filtering 22 _ j i i i i i 0.0 0.2 0.4 0.6 0.8 1.0 time (sec) Figure 1.9: The output of a time-variant low-pass filter applied to three input spikes using frequency function 3 in Figure B-2a. Four filter half-lengths of (a) 100, (b) 200, (c) 400, and (d) 800 ms were used, (e) is the exact output using a half-length of 800 ms. b) Locate that portion of the band-pass fairway encased by Q,i(t) and ilh(t) where the gradients are most severe. c) Compute the theoretical and approximate impulse responses for a spike situated at that location. d) Compare the time and frequency domain representations of both responses. If they are unacceptably different then alter the filter parameters and return, to c) above. Indeed, for many of the time-variant filtering problems in reflection seismology the frequency cutoff gradients are relatively small (ranging from 0 to 25 Hz/s). Thus, a wider range of filter parameters will provide acceptable filter errors. I have found that using a cosine window.with exponent a = 2 and a length of T=200 ms provides excellent results. For example, the field data example shown in Figure 1.1a is a case where time-variant deconvolution after stack enhanced high frequency noise over seismic events later in time. Chapter 1. Aspects of time-variant filtering 23 The band-pass fairway of an optimum time-variant filter for these data varies linearly from 80 and 15 Hz at time zero to 25 and 5 Hz at 6 s which provides a maximum frequency gradient of about 9 Hz/s. This fairway was determined by analyzing the data in filter panels. In this analysis, filter panels each portray a different narrow-frequency-band of the same section of data. By analyzing the data at a succession of frequency bands, it is possible to roughly determine a time-variant frequency fairway that will optimize the signal to noise ratio of seismic events with time. Figure 1.1b shows the results of applying a time-variant filter with the given frequency cutoffs, and the above recommended filter parameters. The higher frequencies inherent in the shallow events are preserved and the deeper events, though lower in frequency content, are unmasked. 1.5 Relation to the Hilbert transform filter In the above development, the time-variant filter was presented in terms of complex signals. It is possible however, to provide a related derivation in which the input signal is kept real. The reasons for doing this are two-fold. Firstly, the derivation leads directly to a classical formulation of a low-pass filter utilizing the Hilbert transform derived by Cain (1972) and extended to the time-variant case by Peacock (1985). And secondly, the possibility of aliasing arising from redundant information in the negative frequency band occurs when real signals are input to a time-variant Hilbert filter. The formulation presented here will show explicitly how such aliasing arises and under what conditions a real input signal can be used safely. The time-variant convolution integral for mixed signals is Here, the time-variant complex filter / operates on a real input signal s and the real part of the output is desired. (1.28) Chapter 1. Aspects of time-variant filtering 24 Development of only the time-variant low-pass filter satisfies the objectives of this section. The amplitude spectrum of the time-variant low-pass filter can be represented as the difference between two step functions. That is, A{u,t) = 2[U(UJ) - u(u - Q(t))}, (1.29) where Q(t) is the instantaneous low-pass frequency function that regulates the time-variable shift applied to the second step function, and the factor of 2 is necessary for proper scaling of the real output. Inverse Fourier transformation of A(OJ, t) gives the zero phase low-pass filter / ( r ,0 = [ l - e m ^ ( r ) . (1.30) Utilizing relation (1.11) gives f(r,t) = [l - e ^ ) - s ( ' - ) ] ] £ ( r ) . (1.31) Substituting this into equation (1.28) and simplifying gives P(t) = Re{ieie^H{s{t)e-16^}}, (1.32) where H[-] is the Hilbert transform given in Bracewell (1978, p. 267). Further reduction of this result leads to p(t) = cose(t)H[s(t)sin6(t)}-sind(t)H[s(t)cos9(t)]. (1.33) Substituting 6(t) = 2ivfct in equation (1.33) produces the classic low-pass Hilbert trans-form filter originally derived by Cain (1972). Peacock (1985) begins with equation (1.33) and develops criteria for automatic filter parameter selection to help shape the spectral response of the Hilbert transform filter. When this filter is time-invariant, his procedure accurately reproduces the desired spectral response. However, the spectral response of the time-variant Hilbert transform filter also depends on the instantaneous frequency Chapter 1. Aspects of time-variant filtering 25 cutoff function and on the filter length due to approximation (1.11). This leads to errors that should be considered when designing an optimum time-variant filter as shown in Section' 1.4. Digital implementation of equation (1.32) or (1.33) must be done with caution. Be-cause the input signal is real, its Fourier transform has energy in both positive and negative frequency bands. This means that the amplitude modulations induce frequency shifts that will cause aliasing if the input signal contains energy at frequencies above some upper frequency limit. To illustrate this problem, I provide a graphical analysis of the. time and frequency operations of equation (1.32) when the instantaneous low-pass frequency is constant. Equation (1.32) is utilized instead of (1.33) because its Fourier transform is less complicated to interpret graphically. Setting 0(£) = wc, gives 9(t) = u>cty and equation (1.32) becomes p(t) Re{ieiu1'1 II'.s{t)e iuJ'1}. ' (1.34) The Fourier transform of the quantity inside the curly brackets is, P(u>) = S(UJ — wc) * [—sgn(uj)S(u> -f wc)]. (1.35) Figure 1.10 shows the graphical representation of equation (.1.34) on the left, and of equation (1.35) on the right, for a low-pass frequency of 50 Hz. The input function used to exemplify the hazards of applying this filter is a 90 Hz low-pass sine function. This function and its frequency spectrum are shown in Figure 1.10a and 1.10b respectively. Again I use a sample interval of 4 ms leading to a Nyquist frequency of 125 Hz. Figure 1.10c shows the results, of amplitude modulating the input signal with a constant fre-quency of 50 Hz (u>c .= 27r50). The real part of the output is shown above the imaginary part. As shown in.Figure l . lOd, the resulting leftward spectral shift of 50 Hz wraps.part of the negative spectral energy into the positive frequency band. Next, applying the Chapter 1. Aspects of time-variant filtering 26 Hilbert transform combined with multiplication by i, creates the complex signal shown in Figure l.lOe. The effect of this operation in the frequency domain, as shown in Figure l.lOf, is to negate the positive frequency spectrum. The results of amplitude demodula-tion are shown in Figure 1.10g and the corresponding frequency shift of 50 Hz to the right is shown in Figure l.lOh. I then take the real part of the complex time signal in (1.10g) and the even part of the real frequency spectrum in (l.lOh) to obtain the final results shown in Figures l.lOi and l.lOj respectively. Observe that the desired low-pass output has been aliased by energy wrapped into frequencies above, the desired 50 Hz cutoff. The amount of aliasing is dependent on the upper frequency limit of the input data and on the (highest) low-pass frequency of the filter. In order to avoid aliasing effects when using real input signals, the following relation must be satisfied us, + uf < wN, (1.36) where ws is the high frequency limit of the input signal, ujf is the high frequency Hmit of the filter, and up? is the Nyquist frequency. When this relation is satisfied, equation (1.32) or (1.33) can be safely appHed, obviating the need to compute the complex signal. Note however, that if real signals are input with knowledge of only UJN, various procedures must be adopted to ensure that relation (1.36) is satisfied. Procedures such as resampling to increase UJN, or anti-alias filtering to decrease ut, may be necessary when using a filter with a prescribed upper frequency limit. These precautions are, of course, unnecessary if the complex signal and equation (1.16) are utilized in the time-variant low-pass filtering. Equation (1.33) involves four real vector multiplications and two real convolutions (Hilbert transforms). Thus, digital application of equation (1.33) may be slightly more cost effective than equation (1.16) for time-variant low-pass filtering. However, the com-putations are doubled for a band-pass filter (computed as the difference between two Chapter 1. Aspects of time-variant filtering 27 0.4 0.0 s(t) -1~ J I I I I I L -0.20 0.0 0.20 a s^t^s^expl^nSOt] 1.0 0.0 Im ' ' ' ' ' ' ' ' -0.20 0.0 0.20 C s2(t)=iH[Si(t)] 1.0 0.0 -0.20 0.0 0.20 e s3(t)=s2(t)exp[i2n50t] 1.0 0.0 ' ' ' ' I I 1 L —0.20 0.0 0.20 g Re[s3(t)] 0.4 0.0 T ' l ' i I I I L -0.20 0.0 0.20 time (sec) I S(o) 1.0 0.0 - \ \ 1 1 > 1 -100 0 100 b S1(U)=S(U+2T\50) 1.0 0.0 1 L J L — 100 100 S2(u)=-sgn(o)S,(a) 1.0 -I _ J U • i -100 0 100 f S3(a)=S2(a-2i\50) 1.0 0.0 J l_ t r _J L. —100 0 100 h Even[S3(cj)] 1.0 0.0 1 1 u 1 1 100 0 100 frequency (Hz) J Figure 1.10: The time and frequency operations of the time-invariant, low-pass filter using real input data. The input signal in (a) is a 90 Hz sine function and the cutoff frequency is 50 Hz. The final output signal in (i) is aliased by energy wrapped into its spectrum above 50 Hz. Chapter 1. Aspects of time-variant filtering 28 low-pass outputs). Nor is phase rotation easily included in this formulation. There-fore, the use of complex signals and equation (1.21) are recommended as a basis for a production algorithm to perform time-variant digital band-pass filtering. 1.6 Frequency-variant windowing Time-variable frequency cutoffs may also be effected by filtering the spectral function of an input signal. I use the term frequency-variant window to distinguish that operation from the time-variant filter operation described in previous sections. The time-variant filter described in Section 1.3 permits a different frequency band at each instant of time, while the frequency-variant window described in this section permits a different time window for each frequency component of an input signal. Both operations can lead to an output signal that contains time-variable frequency cutoffs. The above time-variant filtering algorithm can be generalised for spectral functions to achieve frequency-variant windowing. The derivation is essentially identical to that presented for complex time signals except that operating in the frequency domain intro-duces conjugations that must be accounted for in the equations. I next briefly discuss the important considerations of a frequency-variant windowing algorithm. Equation (1.1) may be rewritten for complex spectral signals as function S. The time domain characteristics of this operator are obtained by applying the inverse Fourier transform (1.37) where F(K;<JJ) is a frequency-variant time limiting operator that filters a complex spectral = A(r,v)e (1-38) Chapter 1. Aspects of time-variant filtering 29 where A and ip are the frequency-variant amplitude and phase characteristics of the filter respectively. Using the analogy of the time-variant high-frequency-pass filter previously described, the amplitude of a frequency-variable long-time-pass operator is represented by a step function that shifts along the time axis with frequency according to a specified short-cut time. That is, A(T,U)=U(T-T(U>)), (1.39) where T(UJ) is the time cutoff function that regulates the frequency-variable shift. Using the time shifting theorem, the forward Fourier transform of equation (1.39) leads~~to the following long-time-pass operator F ( « > W ) = C - ^ C " ) ^ > (1.40) where S(K) = S(K) — Z2/K is the frequency domain equivalent of the complex delta function. Following the development of the time-variant filter I first define a phase function 0 associated with the frequency dependent time cutoff as duj Carrying out a similar derivation to that done for the time-variant filter leads to a general expression for frequency-variant windowing with phase compensation; G H = e - i ( e ' - p ) [ ^ l * e-^e{tM * S(u)ei@>}], (1.42) where, 0j = 0;(u>) = Jo Ti(f)df, is the short-time-cut phase, 0^ = Oh((jj) = So Th(f)df, is the long-time-cut phase, A 0 = Qh(w) - ©,(», Chapter 1. Aspects of time-variant filtering 30 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 078 1.0 Time (sec) Time (sec) Figure 1.11: Two time-frequency pass fairways that can be utilized in a frequency-variant window or a time-variant filter operation. The time-frequency pass in (a) is more precisely implemented using a frequency-variant window, and (b) is more accurately implemented using a time-variant filter. This follows since the time-cutoff gradient in (a) and the frequency-cutoff gradient in (b) are respectively small. and, P = P(w), is a frequency-variable phase compensation function. Error analysis, identical to that provided in Section 1.4, shows that errors in this spectral filter are a function of. the time-cutoff gradient (in s/Hz); smaller gradients in the given time-pass fairway trajectories result in a more precise window response for a given set of filter parameters. It then follows that errors in the frequency-variant window and the time-variant filter are inversely related. In Section 1.4, I showed that errors in the time-variant filter increase with the frequency cutoff gradient of the filter in Hz/s. For example, both types of operations are capable of filtering an input time signal according to the time-frequency pass regions shown in Figure 1.11. However, the time-frequency pass shown in Figure 1.11a is most accurately carried out using the frequency-variant window formulation, where as the time-frequency pass shown in Figure 1.11b is most Chapter 1. Aspects of time-variant filtering 31 accurately carried out using the time-variant filtering formulation. This follows since the time-cutoff gradients in Figure 1.11a and the frequency-cutoff, gradients in Figure 1.11b are correspondingly small. A simple example of this frequency-variant window is shown in Figure 1.12. The input signal (Figure 1.12a) consists of three sinusoids of 5, 20, and 60 Hz, and the desired time windows for each are (0,1), (.2,.8), and (4,.6)s respectively. The complete time-frequency pass window used for this example is shown in Figure 1.13. Figure 1.12 shows each of the operations of equation (1-42) in the time domain with the final result showing the effect of frequency-variant windowing on the input signal. No phase compensation is included in this example. This example concisely illustrates the key to parameter-variable windowing/filtering described in this chapter; each frequency component is independently shifted, operated on, and then restored. Further significance of a frequency-variant window in signal processing/analysis re-mains to be seen: Oppenheim and Schafer (1975) discuss the importance of traditional window functions in signal processing theory. However, to my knowledge, window func-tions given in the literature all treat frequency components equally. It is reasonable to assume that a frequency-variant window could have applications in many areas of signal processing. I do not investigate those possibilities here. 1.7 Conclusion Time-variant filtering is efficiently and concisely formulated using complex signals. This formulation is particularly insightful because of the geometrical interpretation it offers in the frequency domain. The filter essentially separates the pass-band and the reject band of an input complex signal into positive and negative frequency bands. Application of the (conjugate) complex delta function removes the reject band, and subsequent restoration Chapter 1. Aspects of time-variant filtering 32 a ) IJUM 1 UnllliMliiMlli -i i e ) 0 1 1 - 1 . 0 - 0 . 5 0.0 0.5 1.0 time (sec) Figure 1.12: Time domain operations of a simple frequency-variant window, a) Input time signal, b) frequency-variant shift, c) annihilation of positive time energy, d) fre-quency-variant shift, e) annihilation of negative time energy, and.f) frequency-variant shift restoration. Chapter 1. Aspects of time-variant filtering 33 120 -100 -0 I i i i i i i I 0.0 0.2 0.4. 0.6 0.8 1.0 TIME (sec) Figure 1.13: The time-cutoff window utilized to effect the frequency-variant window shown in Figure 1.12. of the pass-band to its proper location completes the filter operations. The filter can be zero-phase or it can effect a time-variant phase rotation of the input data. Imposing a tapering window on the (conjugate) complex delta function allows control of the spectral characteristics of the filter, although spectral rolloff also depends on each instantaneous frequency function used. Higher frequency gradients require a shorter filter to produce acceptable spectral rolloff. However, with these cautions in mind, this time-variant filter can produce excellent, interpretable results. When the derivation is constrained to include real input signals, I show that the filter is related to the Hilbert transform filter of Cain (1972) and Peacock (1985). As a consequence, I quantified two sources of error, when using the Hilbert transform filter in a time-variant mode, that have not been previously elucidated. First, the spectral rolloff depends on the instantaneous frequency cutoff function used. Second, filter aliasing can Chapter 1. Aspects of time-variant filtering 34 be avoided if the highest signal frequency plus the highest filter frequency is less than the Nyquist frequency. Frequency-variant windowing is introduced as a complementary method of effecting a time-frequency pass function that contains large frequency-cutoff gradients, and hence, small time-cutoff gradients. The inverse relation between the errors in the frequency-variant window arid the time-variant filter, allows a user to improve filter performance by choosing the appropriate -method for each desired time-frequency pass function. Considering 2-D seismic data, a space-variant filter is an obvious extension-to the above derivations. And some limited1 kinds of time-space-variant filter operations: can be achieved by cascading time and space-variant filters. It was my hope that by considering complex signals in higher dimensions, a more versatile multi-dimensional filter could be developed. The pursuit of this problem (which remains unsolved) led me to the formulation of local phase velocity from multi-dimensional complex data described in chapter 2. Chapter 2 Local Phase Velocity from Complex Seismic Data 2.1 Introduction Phase velocity measured at a specific location on a seismic record is defined as the velocity of any given phase such as a trough or peak. It is also thought of as the local slope of a travel time surface produced by a seismic event such as a reflection arrival. In exploration seismology, local phase velocity is an indicator of seismic attributes such as head wave velocity (Grant and West, 1965), rms and interval velocity (Gonzales and Claerbout, 1984), and phase alteration with offset (McMechan, 1983). Local dip estimates can be utilized in the migration process (Chun and Jacewitz, 1981), for inverse ray tracing procedures (May and Covey, 1981) and trace interpolation (Bardan, 1987). Indeed, local phase velocity is one of the most important characteristics of a seismic event. The measurement of local phase velocity from seismic records was historically done by hand. Recently, however, automatic methods based on localized slant stacks have been proposed by McMechan (1983) and Milkereit (1987) using semblance analysis along a sweep of linear trajectories of limited aperture defined by a set of trial velocities. The analysis is carried out at each point on the seismic record and the trajectory with maximum semblance is picked as local phase velocity. This technique is robust in the presence of noise but requires the computational intensity of a stacking velocity algorithm (Neidell and Taner, 1971). I propose a different method for the automatic estimation of local phase velocity 35 Chapter 2. Local Phase Velocity from Complex Seismic Data 36 that uses the concepts of complex trace analysis (Taner et al., 1979). A complex seis-mic record is an extension of complex trace analysis into higher dimensions and phase velocity is estimated by finding trajectories . of constant phase. In two dimensions, the calculation of phase velocity reduces to a ratio of instantaneous frequency and wavenum-ber and thus, provides a measure of the dominant plane wave component at each point of a seismic record. The algorithm is simple to implement and computational require-ments are small partly due to a new method of computing instantaneous frequency and wavenumber which greatly simplifies these calculations for two-dimensional (2-D) and three-dimensional (3-D) complex records. In addition, my approach has the advantage that no a priori velocity input is needed; however, optimum stability is achieved when a limited range of dipping events are considered. Preconditioning the record with an appropriate velocity filter (Treitel et al., 1967) helps reduce the detrimental effects of crossing events, spatial aliasing, and random noise contamination. Computation of local phase velocity from complex seismic data is dependent on fast and accurate estimation of instantaneous frequency and wavenumber. I therefore begin with a basic review of this one-dimensional (1-D) problem and show how it relates to the estimation of phase velocity from multidimensional signals. This heuristic view is centered around simple functions such as the sine function that have well defined in-stantaneous properties. Analyzing the sine function in one and two dimensions provides a strong basis for intuitive understanding of local phase velocity estimation using field seismic records. Chapter 2. Local Phase Velocity from Complex Seismic Data 37 2.2 Instantaneous frequency A seismic trace can be thought of as a superposition of harmonic components. Thus, the simplest seismic signal is represented by a single harmonic of the form s(t) = aQcos{ujQt + (/>}, (2-1) where UJ0 is the radial frequency, <f> is the initial phase, and a0 is the amplitude. The Fourier transform of this signal is a0/2{el*£(a> — a;0) + e~l^8{ij + ^ o)}- At any instant of time, this harmonic is described by the two frequencies ±u;0- Thus, instantaneous frequency of a real harmonic appears to be dual valued. This ambiguity is avoided by using the complex exponential form s{t) = o 0e i^ o t +*> > (2.2) where s(t) = s(t) — iH[s(t)] and H[-] is the Hilbert transform (Bracewell, 1978). The Fourier transform of equation (2.2) is simply a0el^8(uj — u0) and the instantaneous fre-quency is now the single value co0. This is consistent with the traditional definition of instantaneous frequency as the time derivative of the phase {co0t + (j)}. Although UJ0 can be either positive or negative without loss of generality in the equations, I shall consider complex data with only positive instantaneous frequencies since they are traditionally more acceptable (Gabor, 1946; Taner et al., 1979). I do not impose the same restriction on wavenumber content. Indeed, I later rely on both positive and negative wavenumbers to estimate vector valued phase velocity. An important function in signal processing theory is the bandlimited sine function, often used as an idealized seismic wavelet. The analytic expression for a complex ban-dlimited sine function shows that it has a simply defined instantaneous frequency. The expression for a such a function containing spectral energy between ui and o>2 and with Chapter 2. Local Phase Velocity from Complex Seismic Data 38 initial phase (j) is s(t) = a{t)ei{Qt+*}, (2.3) where, a(t) =' 2sin\!^f^\t/7rt, <I> = (u>2+^i)/2, s(0) = e**(u/2 - W O / T T . The instantaneous frequency of the complex sine function is simply ui or the cen-ter frequency of its passband. This is consistent with the work of.Mandel (1974.) and Robertson and Nogami (1984), who show that instantaneous frequency is a local measure of the spectral centroid. In particular, Taner et al. (1979) show that the instantaneous frequency of a Ricker wavelet is a good estimate of the Fourier centroid when the analysis is located around the envelope maximum. Boashash and Whitehouse (1986) show that for an arbitrary signal, instantaneous frequency is viewed as the spectral centroid of the Wigner-Ville distribution. A complex seismic trace y can be computed by convolving the real trace y with a complex sine function of suitable bandrange, y(t) = s(t)*y(t)=A(t)eie(t\ (2.4) where A{t) is the amplitude, and 6{t) is the phase of the complex trace. When the pass-band of s extends between zero and the Nyquist frequency for discrete signals, y becomes the complex trace of Taner et al., (1979), while s becomes the complex delta function described in Chapter 1. Instantaneous frequency of y typically exhibits variations from the center frequency of s due to local structure of seismic events, changes in spectral character of the propagating source wavelet, and noise. Several methods have been proposed for estimating an instantaneous frequency func-tion from the complex signal. Taner et al. (1979) estimate instantaneous frequency by Chaptef 2. Local Phase Velocity from Complex Seismic Data 39 computing the local time derivative of instantaneous phase, , v d6(t) r 1 dy(t), , x ^=-r = I^m^]- (2'5) Boashash and Whitehouse (1986) suggest calculating the spectral centroid of the Wigner-Ville distribution. However, these methods become computationally restrictive in higher dimensions. I therefore present another viewpoint on this subject which leads to simpler calculations of higher dimensional attributes. Theoretical work done on the relationship between the time derivative of instanta-neous phase of the complex signal and the Fourier spectrum and Wigner-Ville distribu-tions leads to a pragmatic definition of instantaneous frequency. Relating local frequency measures to a spectral distribution aids intuitive understanding of complex signal charac-teristics. I further define instantaneous frequency to be the frequency of a single complex harmonic that locally fits the complex signal in an arbitrary time span around the point of interest. This definition leads to an alternate, yet equivalent method of calculating instantaneous frequency in addition to furthering the pragmatic understanding of this attribute. An estimate of instantaneous frequency at a point t0 can be obtained by locally fitting a complex harmonic of the form of (2.2) to the complex signal. Thus, I desire to find an u0 such that y(t) = a0e^<"+*>. (2.6) in a region around t — t0. There are clearly many ways to effect this fit. The simplest, and one that is appbcable if the data are accurate is to require that (2.6) passes exactly through the points y(to) and y(to + St). This yields a pair of equations y(t0) = a 0 e^ + *>. y(t0 + 8t) = o^Wo+rtH*}. Chapter 2. Local Phase Velocity from Complex Seismic Data 40 which may be solved to obtain the desired instantaneous frequency . ' (2.7) where for the complex number z = x + iy, arg(z) = arctan(?//a:) with the branch of the arctan chosen so that arg(z) lies between ±7r. Thus, upper and lower frequency bounds on u>o are ±7r / £ t . Choosing 8t equal to the signal sample interval provides limits on u>0 of the Nyquist frequency. Equation (2.7) avoids the computation of dy/dt necessary in equation (2.5) and, in addition, retains the simplicity of a finite difference approximation to (2.5) given by Yilmaz (1987,p.521). It is important to verify that the two-point formula (2.7), and therefore, my definition of instantaneous frequency, is consistent with previous work. Figure 2.1 shows an example comparing (2.7) with the formulae of Taner and Yilmaz using an input complex sine function with a band range of 10 to 50 Hz. My formula and that of Taner produce the true result of 30 Hz, however, Yilmaz' finite difference approximation shows significant errors that render it unsuitable for my purposes. This example was carried out using an analytic expression for the complex sine function. The comparison using a more typical seismic trace is shown in Figure 2.2. A complex trace was first computed from the real trace shown in Figure 2.2a before computing the instantaneous frequency functions using formulae (2.5) and (2.7). In Figure 2.2b,-the solid curve results from Taner's formula (2.5) and the dashed curve results from application of my formula (2.7). The two curves agree very well for the most part; however, disagreement exists at locations where estimation of instantaneous frequency becomes unstable1. In those areas, I truncated the output from equation (2.5) at the zero and Nyquist frequency boundaries. Truncation of the output from equation (2.7) is 1 Formula (2.5) often produces instantaneous frequencies well outside the Nyquist bounds at locations where the signal envelope is small. Unstable results are also encountered when the instantaneous phase function contains unusually rapid variations. Hence, I have truncated the anomalous results from this formula in Figure 2.2b to facilitate comparison in areas of stable estimation. Chapter 2. Local Phase Velocity from Complex Seismic Data 41 unnecessary because instantaneous frequency is automatically constrained within those boundaries2. The input data are sampled at 4 ms intervals. From these and other examples, I conclude that my definition of instantaneous fre-quency as the frequency of the best fitting local complex harmonic at each point of a complex seismic trace is consistent with previous work done by Taner et al. (1979) on this subject. Reconciliation in those areas where estimation of instantaneous frequency becomes unstable remains a subject of further research. However, philosophical differ-ences can be characterized. The formulation of Taner et al (1979) provides little insight in how to generalize instantaneous attributes to higher dimensions. On the other hand, thinking of instantaneous frequency as representing a local fit to a complex harmonic function is immediately extendable. As I discuss next, in 2-D, one can think of fitting to a complex harmonic plane wave. 2.3 Extension to 2-D In two dimensions, a seismic record is decomposed into plane wave components by the Fourier transform. Thus, the simplest complex seismic record is represented by a single complex plane wave. Such a plane wave of frequency UJ0 and wavenumber k0 is given by s(x,t) = a 0 e i { , u o t - k ° x + * } , (2.8) where a0 is the amplitude and <f> is the initial phase. The Fourier transform of s(x,t) is a0e1'^8(k + k0,uj — u>0). This complex plane wave exhibits moveout or dip that can be characterized in terms of phase velocity, vp = uj0/k0, (2.9) 2 Negative frequency estimates sometimes arise from my formula at locations where frequency calcu-lation becomes unstable. This is undesirable since the input complex signal does not contain negative frequency information. However, I have found that such anomalous values can be satisfactorily un-wrapped by the addition of IT to the calculated radial frequency. -io]-J— 1 — — —' : -1 L j 0.0 0.2 0.4 0.6 0.8 1.0 Time (sec) Figure 2.1: (a).Real part of a complex 10 to 50 Hz sine function, (b) Instantaneous frequency of the complex sine using three different formulas. The straight line at 30 Hz results from equations 2.5 and 2.7, while the oscillating curve results from a finite difference approximation to equation 2.5 given by Yilmaz. The analytic formula for the complex sine function was used in this example. Chapter- 2. Local Phase Velocity from Complex Seismic Data 43 0.8 ( a) 0.4 0.0 - 0 . 4 - 0 . 8 i i i i t i 0.0 0.2 0.4 0.6 0.8 1.0 TIME (sec) Figure 2.2: (a) Input synthetic trace. I have applied a time-variant band-pass filter to a reflectivity series to obtain the result shown, (b) The solid curve is the result of using formula (2.5) and the dashed curve is the result of using my two-point formula (2.7) to calculate instantaneous frequency from the complex version of the trace in (a). Chapter 2. Local Phase Velocity from Complex Seismic Data 44 which is a simple ratio of frequency and wavenumber. Plane waves whose frequency and wavenumber vary proportionately are considered locally similar since they display the same phase velocity. My purpose here is to distinguish seismic events on the basis of phase velocity. The intuitive advantages of the bandbmited sine function can be extended into two dimensions, where it can be thought of as an idealized 2-D seismic wavelet. The 2-D complex sine function is defined by a rectangular pass region in the Fourier k — ui domain and has a simply defined phase velocity. The analytic expression for such a function containing spectral energy in the frequency band between UJ1 and UJ2 a n < i the wavenumber band between ki and k2 is given by (see Appendix A) s(x,t) = a ( a ; , i ) e i { ^ ^ x + * } , (2.10) where, a(x,t) = 4 : s i n \ ^ ^ \ x • sin\^^\t/TT2xt, UJ = (a; 2+^i) /2, k = (k2 + k1)/2, i(0,0) = e**|u;2-Wil • \k2 - k^/n2. Recognizing that instantaneous frequency is the partial time derivative, and instanta-neous wavenumber is the negative partial space derivative of the phase {uit — kx + (/>}, the phase velocity of this function is simply u>/k, where (k,UJ) is the centroid of the rect-angular pass region in its Fourier k — UJ spectrum. Local phase velocity is still envisioned as a measure of the Fourier centroid when more realistic two dimensional seismic wavelets are considered that contain nonuniform k — UJ spectra. As I show later, calculation of phase velocity using estimates of instantaneous frequency and wavenumber is viewed as finding the best fitting complex harmonic plane wave component at each point in the complex seismic record. Chapter 2: Local Phase Velocity from Complex Seismic Data 45 A complex seismic record can be created by convolving the real seismic record with an appropriate two dimensional complex sine function, y(x,t) = s(x,t) * *y(x,t), (2.11) where ** denotes a two dimensional convolution. Estimates of local phase velocity from y may deviate from the phase velocity of s due to local structure of seismic events, space and time variation in the propagating wavelet, and noise. However, Figure 2.3 shows the special case where a harmonic plane wave, a two dimensional sine function,-and a dipping event all have the same phase velocity. The spectral centroid of each event is (.005 cyc/m, 20 cyc/s) producing a common phase velocity of 4000 m/s. Hence, these events are considered locally similar despite differences in global character. The above development provides a basis for intuitive understanding of phase velocity computation using a complex seismic record but does not lend itself to more complicated situations. A more general treatment of the problem in three dimensions is rigorously handled by use of the Fourier transform. In the analysis that follows, I assume that the recording plane is z = 0 and that seismic signals are therefore represented as functions of x, y, and t. 2.4 General formulation The space-time Fourier transform of a 3-D seismic record g(x,y,t) is defined as . OO OO OO (2.12) — OO —OO —OO and the inverse Fourier transform is OO OO OO — OO — OO —OO (2.13) Chapter 2. Local Phase Velocity from Complex Seismic Data 46 0.44 plane wave 2-D sine dipping event 0.006 0.012 wavenumber 0.006 . 0.012 Wavenumber 0.006 0.012 Wavenumber • Figure 2.3: (a) A plane wave of frequency 20 Hz and wavenumber 0.005 cyc/m. (b) A 2-D sine function with frequency passband 0 to 40 Hz and wavenumber passband 0 to 0.01 cyc/m. (c) A seismic event with a frequency band range of 0 to 40 Hz dipping at 4000 m/s. Below each section is the corresponding ideal Fourier ui — k domain amplitude response. The theoretical phase velocity of all three events are identical. Chapter 2. Local Phase Velocity from Complex Seismic Data 47 The three dimensional Fourier transform decomposes a seismic wavefield into its plane wave components. Each of the four independent positive frequency octants of the Fourier transform contains plane wave components of differing emergence angle (dip direction). The four negative frequency octants are conjugate symmetric with respect to the posi-tive frequency octants and thus, provide redundant information about local dip. Con-sequently, reconstructing a space-time function from only the positive frequency compo-nents of an input seismic record provides a complex record suitable for phase velocity estimation. Inverse Fourier transformation over the positive frequency axis produces OO OO OO g{x,y,t) = — j J J G^ky^y^-^'-Wdktdkydw 0 — o o — oo = g(x,y,t)-iHt[g(x,y,t)}. (2.14) where, Ht[-] is the Hilbert transform with respect to t. Furthermore, velocity filtering is often necessary to reduce unwanted noise and to enhance signal over a limited dip range. This may be accomplished while creating a complex record by modifying the integration limits in the wavenumber portions of the inverse Fourier transform. That is, OO w/vy2 w/Vv2 f(x,y,t) = -^J J j G{kx,ky,u)ei^t-^x-k^dkxdkydw = R{x,y,t)+il(x,y,t), (2.15) represents a complex seismic record with a limited range of dipping events, where vxi,vx2, and vyi,vy2 define velocity pass fans in the x and y directions respectively and R and / denote the real and imaginary parts of the complex record f. The above complex representation permits expression of the seismic record in terms of local amplitude and phase as r(x,y,t) = A(x,y,t)eie^t\ (2.16) Chapter 2. Local Phase Velocity from Complex Seismic Data 48 where, A(x,y,t) = y/R*(x,y,t) + P(x,y,t), (2.17) is the amplitude and ='^-'{fel)}' (2-18) is the local phase. Local phase velocity may be derived by first noting that a seismic event lies along a surface of constant phase (neglecting dispersion and other possible phase alterations incurred by the propagating energy). Along such a surface S, the total phase is conserved, that is, f | , = £ + * . V , = 0. (2.19) where, v = vxi+vyj, is a horizontal phase velocity field, and, V# = f j^ + f^'i) i s the spatial phase gradient. The horizontal phase velocity I seek is normal to wavefront propagation along the recording surface. It naturally follows that the local phase velocity pointing in the direction of the spatial gradient is the key result. Therefore, substituting v = aV# into the above equation and solving for a provides the following expression for local phase velocity _ -do/at ve *=mm\- (2-20) When the phase surface is smooth, the above phase velocity is a single valued vector function. However, at points where the phase surface is discontinuous or is crossed by a secondary event surface, ambiguities in phase velocity computation arise. These com-plications, along with the problems of spatial aliasing and random noise contamination, will be considered in Section 2.5. In two dimensions, expression (2.20) reduces to -d6{x,t)/dt V p ( M ) = de(x,t)/dx- ( 2 2 1 ) Chapter 2. Local Phase Velocity from Complex Seismic Data 49 Here, phase velocity is seen to be the ratio of instantaneous frequency and wavenumber with a sign reversal that provides the correct dip direction. Generalizing my previ-ous definition of instantaneous frequency, I now define the instantaneous frequency and wavenumber to those values of a complex harmonic plane wave that locally yield a fit to the complex seismic data. Again, when the data are accurate, a simple formula for estimating the best fitting plane wave near the point (x0, t0) can be generated by requiring f(x0 + Ax,t0) = a o e i ( -oto-fco[xo+Ax]) r(x0,to + At) = a 0 e i { m l t o + A t ] - k o x o l These relations are easily manipulated to yield v(xt) ^arg[r(x,t + At)/r(x,t)} p [ X ' } ~ Atarg[r(x + Ax,t)/r(x,t)Y K ] where Aa; and At are the spatial and temporal sample intervals. Figure 2.4a shows the results of applying this formula to the three records shown in Figure 2.3. Some edge effects are observed in the plane wave and dipping events; this is due to a smoothing window used to taper the trace ends. Those areas will not be included in the following analysis. Figure 2.4a shows an accurate recovery of the true local phase velocity (4000 m/s) for each type of event, however, the numerical estimation is not exact. A standard error for each event is defined as cr = 1 N M ^ E ( ^ - 4 0 0°) 2> (2-23) where V{ is the estimated local phase velocity value, and the sum is taken over data values where the event envelope is greater than 90% of its maximum to analyze the most significant portions of the events. Standard errors calculated from the phase velocity Chapter 2. Local Phase Velocity from Complex Seismic Data 50 images of the plane wave, the sine function, and the dipping event shown in Figure 2.4a are 5, 20, and 2 m/s respectively (rounded to the nearest integer). Numerical evaluation of the localized sine function is evidently more difficult. This example illustrates that estimation of local phase velocity using numerical tech-niques is not exact. Within this context, Appendix C attempts to assess the question of limits on the resolution of local phase velocity given that the data are accurate yet finite in spatial and temporal extent and discretely sampled. Event curvature and noise are other considerations which impede accurate calculation of local phase velocity of seismic events. These considerations are discussed next. 2.5 Effects of curvature and noise 2.5.1 Curvature Equation (2.22) is an explicit three point forward estimator of local phase velocity. When the signal is a 2-D sine function or a linear event, this formula gives reasonably accurate values of the phase velocity as shown by Figure 2.4a. However, if the trajectory of the event is curved, then (2.22) will produce estimates that are too high or too low depending on the curvature of the event in the x — t domain. To minimize this bias, additional neighboring points are included in the formula. I illustrate the curvature problem by considering the following estimators: (1) the 3 point forward estimator equation (2.22), (2) the 3 point backward estimator vp(x,t) Ax arg[s(x,t)/s(x,t — At] (2.24) At arg[s(x,t)/s(x — Ax,t)] (3) the forward and backward average, and Chapter 2. Local Phase Velocity from Complex Seismic Data 51 a . mm 7 - - 1 ( b ) B B S _ •• c Hfifl P . J • . » w" mm mMWw 1 1 • Figure 2.4: Local phase velocity images calculated (a) from the accurate data shown in Figure 2.3, and (b) from the same data containing 10% uniformly distributed random noise. Dark areas of the image represent phase velocity values outside the bounds given. Note that the color scale is non-linear at both ends. Chapter 2. Local Phase Velocity from Complex Seismic Data 52 (4) the 5 point forward/backward estimator equation • • A fflrarf(x't+At) + r*(a-t? i vp(x,t) = -¥ ['[fn{TT (2-25) The above estimators are compared using a 1-D complex sine function with hyperbolic trajectory simulating a reflection event from a horizontal subsurface boundary s(x,t) = V 0 1 >.ei«{t-^tl+x*iv*) ( 2 2 g ) *(t-y/tZ + z>/Vi) where to is the normal incidence travel time, x is the offset, and V is a constant medium velocity. I compare both sides of the equation3 V = ^vp(xtt) (2.27) at locations along the hyperbolic trajectory. The right hand side is numerically calculated using estimators (1) through (4) on the complex function (2.26) above for a given velocity V =2000 m/s. And, I use a hyperbolic sine function with a frequency band of 5 to 40 Hz. The results are illustrated in Figure 2.5. All estimators show greatest error at near offsets where the curvature is greatest, that is near the apex of the hyperbola. Estima-tor (3), the forward and backward average, shows the most stability in those areas of greatest curvature. Nevertheless, I have used equation (2.25) for many of the examples in this thesis. I found this formula to provide a good compromise between accuracy and computational efficiency for use on synthetic and field data, cases. Multi-point estimators can be envisioned that include more neighboring data points, however, this may be defeating the objective of getting a localized, estimate. In examples I have tested that contain localized events (such as an isolated dipping event), phase veloc-ity estimates are severely degraded by considering a 9 point forward/backward estimator 3See Section 2.6.2 on rms velocity for the derivation of this relation Chapter 2. Local Phase Velocity from Complex Seismic Data 53 or average. Even with accurate data, instabilities in the computation of instantaneous parameters (as seen in Figure 2.4a) suggests caution with respect to increasing the order of multi-point estimators. Median type operators may be more appropriate in reducing the effect of such outliers (Bednar, 1983). 2.5.2 Noise The types of noise that affect the calculation of local phase velocity for a single seismic event are (a) crossing events, (b) spatial abasing, and (c) random and systematic noise. Each type of noise distorts the phase surface of the primary event and this in turn leads to fluctuations in the estimated phase velocity. For instance, at locations where events of opposite dip cross, both positive and neg-ative wavenumbers exist simultaneously; the estimated spectral centroid is biased in the crossover region and this leads to anomalous estimates of phase velocity there. In fact, local phase velocity is multi-valued at any crossing point, regardless of the sense of the interfering events. It is well known however, that crossing events can often be effectively separated by appbcation of an appropriate velocity filter (Treitel, et. al., 1968). Precon-ditioning the seismic data by selective velocity filtering can often preserve information in the crossover region. . For the case of spatial aliasing, higher frequency components of the primary event appear to dip in an opposing direction and thereby create the effect of crossing events. Appbcation of an appropriate velocity \filter allows the effects of spatial aliasing to be locally analyzed. However, a steeply dipping seismic event bkely contains a lower band of frequency information that is not spatially abased. Thus, a highcut frequency filter can often remove this problem. The effects of random noise on the estimation of local phase velocity are illustrated in Figure 2.4b. Uniform random noise of maximum ampbtude 10% of the maximum data Chapter 2. Local Phase Velocity from Complex Seismic Data 54 2150 2100 2050 & 3 2000 1950 1900 -\1) CURVATURE EFFECTS -\2) ( a ) I I I i 0 3000 0.8 o CO 1000 2000 OFFSET (meters) 3000 Figure 2.5: Comparison of local phase velocity estimators for a hyperbolic reflection event. In (a), 1) is the forward, 2) is the backward, 3) the forward and backward average, and 4) the 5 point estimator, (b) The hyperbolic reflection trajectory analyzed in (a). Chapter 2. Local Phase Velocity from Complex Seismic Data 55 amplitude were added to the data sections shown in Figure 2.3 before utilizing formula (2.25) to calculate the phase velocity. The standard errors of phase velocity calculation using the noisy plane wave, sine function, and dipping events are 379, 341, and 422 m/s respectively, or roughly, 10% standard errors in each case. I used data above 90% envelope threshold, described previously, to estimate these errors. Two possible methods of minimizing the effects of random noise are velocity filtering and utilizing multi-point estimators. After a velocity filter with a pass range of 3000-5000 m/s was applied to the input noisy data, the standard errors of the calculated phase velocities were reduced to 54, 215, and 36 m/s respectively. By contrast, use of an order 9 multi-point estimator (of the form given in (2.25)) produces standard errors of 18, 2123, and 1615 m/s respectively. This example indicates that the amount of stability achieved by each approach depends on the type of event being analyzed. Pervasive harmonic plane wave type events respond well to multipoint estimators, whereas localized events respond much better to velocity filtering prior to use of the simpler formula (2.25). Static shifts are an example of a systematic noise problem. This type of noise in-troduces a phase velocity error that extends forward in time at a fixed receiver location (static shifts are evident in Figure 2.6b). Static corrections must be made prior to the estimation of local phase velocity if analysis of underlying reflection events are desired. Further examples of phase velocity stabilization in the presence of noise are shown in the next Section. 2.6 Examples and related attributes 2.6.1 Synthetic example Figure 2.6a shows a synthetic shot gather based on the velocity-depth model shown in Figure 2.7 and listed in Table 2.1. Two-way travel times to all reflectors were computed Chapter 2. Local Phase Velocity from Complex Seismic Data 56 28 «3 60 80 100 120 60 SO 100 120 iSliiiS f ; f iM:; fJ • • • H i i l l l l l-5T§ a ) 28 48 60 88 100 120 Figure 2.6: (a) A synthetic shot gather generated from the velocity-depth model in Figure 2.7. It is a kinematic representation of the primary (compressional) reflective events, (b) A field shot gather from Western Canada, (c) Reflective events extracted from the field gather, after statics corrections, by selective velocity filtering. by the ray path integral equations (Grant and West, 1965, p. 134) that exactly describe wave propagation in a horizontally layered earth. For this example, the reflection travel time curves within the receiver spread are approximately hyperbobc. There are 60 traces on either side of the shot point (symmetrical split spread field geometry) with trace spacing of 33 m. I placed a unit reflection coefficient (with true polarity) at each two-way traveltime location and convolved the resulting spike traces with a 0-35 Hz bandpass filter operator to produce the kinematic representation of reflective events shown in Figure 2.6a. Figure 2.8a shows a local phase velocity image obtained from the synthetic shot gather in Figure 2.6a using equation (2.25). In that example I computed the complex record using equation (2.14). As shown, the estimation of instantaneous parameters is noisy in Chapter 2. Local Phase Velocity from Complex Seismic Data 57 areas of low signal energy, hence, I employed amplitude thresholding in Figure 2.8b to display only the significant events. Although the threshold level is arbitrary in any case, this operation greatly enhances the interpretability of the image. It is clearly seen that the resultant phase velocity image is a very accurate representation of local dip displayed on the input section; for example note the straight line drawn using the estimated phase velocity at the tangent point. Once accurate phase velocity is available from a seismic shot gather, many other seismic attributes are immediately available. For example, note that a trajectory of constant phase velocity defines a Snell trace (Claerbout, 1985). I have drawn in a few Snell traces in Figure 3.2. Theoretically, each Snell trace can be inverted for subsurface velocity-depth information (Gonzales and Claerbout, 1984). This subject is covered in detail in Chapter 3. However, prior to formal velocity-depth inversion, some approximate relations may be utilized to provide a rough analysis of subsurface parameters using local phase velocity estimates. I next give a brief discussion of two such calculations; rms velocity and maximum depth estimates (Milkereit,1987). 2.6.2 Estimation of rms velocity Here, I provide an approximate relation between rms velocity and local phase velocity. This connection is made more rigorous in Chapter 3. Figure 2.9a shows the result of converting the phase velocity image in Figure 2.8b into local rms velocity using the hyperbolic travel time approximation t2 = t\ + x2/Vl^ns, where t0 is the normal incidence two-way travel time and x is the source-receiver offset. Differentiation of this relation produces the following dependence on local phase velocity (2.28) Chapter 2. Local Phase Velocity from Complex Seismic Data 58 velocity (km/s) 4 6 Figure 2.7: (a) Interval velocity versus depth model used to produce the synthetic shot gather in Figure 2.6a. (b) Corresponding rms velocity, and (c) average velocity mapped from two-way normal incidence time to the depth coordinate. Actual model quantities are given in Table 2.1. Figure 2.8: (a) Raw phase velocity image computed from the synthetic shot gather, (b) Phase velocity of significant events only. Chapter 2. Local Phase Velocity from Complex Seismic Data 60 Model Parameters R v (m/s) Az (m) z (m) VTms (m/S) Vave (W S) 1 2000 100 100 2000 2000 ' .2 2742 137 237 2400 2371 3 2689 133 370 2500 2476 4 4153 207 577 3000 2896 5 3291 329 906 3100 3028 6 4187 628 1534 3500 3415 7 4399 220 1754 3600 3513 8 6141 307 2061 3900 3753 9 4969 248 2309 4000 3854 10. 4655 465 2774 4100 3968 11 4842 484 3258 4200 4078 Table 2.1: Model parameters used to produce the synthetic shot gather shown in Figure 2.6a. Parameters are listed for each layer with bottom reflecting boundary "R". Corre-sponding normal-incidence rms and average velocities for each reflector are also listed. As shown in Table 2.2, measured rms velocities compare very well with the normal-incidence rms velocities computed from the model in Figure 2.7. The measured values listed are the result of averaging rms velocities along the traveltime curve for each event. But for those events where traveltime curves cross, only the near offset traces were included in the average. Theoretically, such estimates of local rms velocity could be used to perform an offset dependent moveout correction; a necessary operation before common depth point stacking. I do not investigate that possibility here. 2.6.3 Maximum depth bounds Rigorous depth inversion is possible using local phase velocity estimates from a shot gather as shown in Chapter 3. However, a maximum depth bound associated with each reflection event can be quickly calculated from local phase velocity estimates. Assuming horizontal layering and general velocity increase with depth, I begin with the straight Chapter 2. Local Phase Velocity from Complex Seismic Data 61 Measured Normal-Incidence Measured True R Rms Velocity Rms Velocity Depth Bound Depth 1 2006 2000 100 100 2 2427 2400 243 . 237 3 2510 2500 377 370 4 3031 3000 606 577 ' •5 3091 3100 927 906 6 3513 3500 1586 1534 7 3595 3600 1798 1754 8 3935 3900 2167 2060 9 4009 4000 2407 2309 10 4110 4100 2878 2774 11 4203 4200 3363 3258 Table 2.2: Comparisons of rms velocities (m/s) and depth bounds (m) obtained from the images in Figure 2.9 with the true normal-incidence rms velocities and true depths of the model in Figure 2.7. The measured values represent mean values taken along the traveltime curve of each reflector. "R" indicates reflector number. ray path approximation, (x/2)2 + z2 = (t/2Va)2, relating average velocity Va and depth z to offset x and two-way travel time t. Then, rms velocity is used as an upper bound on average velocity (Al-Chalabi, 1974) to arrive at the following relation between maximum depth of reflection and local phase velocity; zmax{x,t)^^^\vp{x,t)\-l. (2.29) Figure 2.9b shows the maximum depth image associated with the synthetic shot gather. This result provides a quick look at approximate reflector depths if one keeps in mind that the correspondence with true reflector depths worsens with increasing depth. Ta-ble 2.2 compares the measured depth bounds with the true depths of each reflector. The measured depth bounds listed represent average values along each traveltime curve determined as given above for rms velocity averages. Figure 2.9: (a) Local rms velocity and, (b) maximum depth penetration images derived from the local phase velocity image of the synthetic shot gather. Chapter 2. Local Phase Velocity from Complex Seismic Data 6 3 2.6.4 Field data example The field shot gather shown in Figure 2.6b contains the same field geometry and trace parameters given above for the synthetic example. The associated local phase velocity image, shown in Figure 2.10a, was computed using formula (2.25). Phase velocities corre-sponding to the first break arrivals are well determined and suitable for useful headwave analysis. And estimated local phase velocities are in general a good representation of lo-cal dip present in the original data. However, the input gather contains ground roll and refraction wave trains that cross reflection events causing distortion of the local phase velocity in crossover zones. Effects of random noise and static shifts are also evident in these data. Although equation (2.25) can provide accurate local dip information of the input data (signal plus noise), phase velocities of underlying seismic events are biased by distortions in the field data. Appbcation of appropriate velocity filters provides the best solution to this problem. The stabilized phase velocity image shown in Figure 2.10b was obtained by analyzing the original data using four separate velocity analysis windows. I separated the positive dipping events from the negative dipping events, and I also separated the reflection events from the ground roll and refraction events, using appropriate velocity windows. Ground roll and headwave arrivals were passed with a velocity filter that ranged from ±300 to ±5000 m/s, while the reflection events were passed using a velocity range from ±5000 m/s on up to ±oo m/s. I then combined the results from each analysis using amplitude thresholding (the event with the highest ampbtude is displayed). Using this approach, I attempt to close in on the local phase velocities of underlying global seismic events by reducing variations in the phase trajectories caused by crossing events and the local effects of noise. Figure 2.11a shows a local rms velocity image created using only the reflection events Figure 2.10: (a) Raw phase velocity image computed from the field shot gather, (b) Phase velocity image obtained by analyzing the field record in separate velocity windows. Chapter 2. Local Phase Velocity from Complex Seismic Data 65 Figure 2.11: (a) rms velocity derived from the local phase velocity image corresponding to the reflective events of Figure 2.6c and (b) Maximum depth bounds extracted from the rms velocity image in (a). The white curve in (a) represents a reflection hyperbola defined by the estimated rms velocity at the cursor location. Chapter 2. Local Phase Velocity from Complex Seismic Data 66 extracted from the raw shot gather in Figure 2.6b. The reflection events shown in Figure 2.6c were extracted from the raw gather, after statics corrections, by velocity filtering (with a pass range of-5000 to 5000 m/s). Variability in rms velocity along the reflection travel time curves could be due to failure of the hyperbolic travel time assumption, but is more likely due to remnant noise effects. Averaging velocities along selected reflection curves provides rms velocity estimates of 2600 m/s at .3 s, 3350 m/s at .9 s, and 3930 m/s at 1.3 s of normal incidence travel time. These results agree, in general, with an independent stacking velocity analysis performed on the raw data, and I therefore remain optimistic about the possibilities of using this approach as an alternate method of automatic rms velocity analysis. Figure 2.11b shows a maximum depth image derived from the local phase velocity image derived from the reflection data in Figure 2.6c. Average depth bounds range from 390 m at .3 s to 1518 m at .9 s to 2560 m at 1.3 s of normal incidence travel time. These values were obtained by taking the mean value along a portion of the reflection curve corresponding to the given normal incidence time. 2.7 Image Analysis The estimation of local phase velocity from seismic data as presented in this chapter is scale independent. That is, the local dip of a seismic event is not dependent on the frequency (u>) and wavenumber (kx) bandwidth used to carry out the analysis. A notable exception occurs if the propogating seismic energy exhibits dispersion. When the velocity of a particular seismic phase varies with frequency, computation of local phase velocity will differ with each frequency component analyzed. Similarly, a 2-D spatial image such as a photograph generally contains variations of Chapter 2. Local Phase Velocity from Complex Seismic Data 67 scale that must be considered in a local dip analysis. A higher wavenumber4 band region of an object in the photograph will Hkely exhibit different local slope characteristics than a lower wavenumber band. Therefore, when utilizing the procedure outbned above for computing local dip in an arbitrary image, one must have in mind what scale or wavenumber bandwidth is of most interest. Simple appbcation of equation (2.22) to a complex version of the original broad band image may not provide sensible answers to specific questions about local dip structure. This is because the calculated slope represents an estimation of the centroid of local wavenumber content, arid therefore, measures an average slope for all scales present. With this concern in mind, local-slope, calculation by use of complex data analysis should provide fertile ground for research in the field of pattern recognition. 2.8 Conclusion I have developed a new method for the automatic estimation of local phase velocity that uses the concepts of complex trace analysis. Defining instantaneous frequency as the frequency of the best fitting complex harmonic at each point of a complex trace provides insight into the generabzation of instantaneous attributes to higher dimensions. In two dimensions, local phase velocity is envisioned as finding the best fitting complex harmonic plane wave at each point of a complex seismic record. The calculation reduces to a ratio of instantaneous frequency and wavenumber with a change of sign. The algorithm is simple to implement and computational requirements are small partly due to a new method of computing instantaneous frequency and wavenumber which greatly simphfies these calculations for 2-D complex seismic records. In addition, my approach has the advantage that no a priori velocity input is needed; however, optimum stabibty is achieved 4 Spatial wavenumber ky replaces frequency as the second spectral dimension in a spatial image Chapter 2. Local Phase Velocity from Complex Seismic Data 68 when a limited range of dipping events are considered. Preconditioning the record with appropriate velocity filters helps reduce the detrimental effects of crossing events and random noise contamination. Accurate local phase velocity estimates of underlying seismic events lead to the cal-culation of many other seismic attributes. Maximum depth bounds and rms velocity are two such attributes easily extracted from shot data. A more complex appbcation of local phase velocity calculation from a shot gather is the subsequent estimation and inversion of Snell traces. A Snell trace is a curve connecting constant phase velocity values on successive primary reflective events. This curve provides information on the traveltime and horizontal travel distance of a ray path within layers bounded by reflections and is directly invertable for interval velocity and thickness of subsurface layers. I provide the details of this subject in Chapter 3. Chapter 3 Snell Trace Estimation and Inversion 3.1 Introduction The efficient recovery of local phase velocity described in Chapter 2 leads naturally to the estimation of Snell trajectories from shot records. A Snell trace1 is defined as a path connecting constant phase velocity in the x — t plane of a shot gather (Claerbout, 1985). For a horizontally layered subsurface, the Snell trace corresponding to primary reflection events (see Figure 3.1) has a significant physical interpretation. Such a trajectory locates all primary reflection energy corresponding to a fixed takeoff/emergence angle. In terms of ray theory, at time t/2, the tip of a ray is at the horizontal position x/2. When v — v(z), Snell's law dictates that the horizontal projection of the ray velocity is the constant 1 vi . Vp=~ = — p sino1 where, vp is termed horizontal phase velocity, 6^ is the initial launch angle and p is the ray parameter or horizontal slowness. Hence, energy leaving the source at a given takeoff angle 6i emerges at the surface with a common phase velocity. The time and space distribution of this energy in the recorded shot gather defines a Snell trace. In this Chapter, I show that Snell trace estimates provide data for the inversion of velocity and depth in a plane layered earth system; interval velocity and layer thickness can be expressed in terms of the time derivative of a Snell trace. 1 I also use the terms "trajectory", "curve" and "path" as post modifiers of Snell; these are meant to be synonymous with "trace". 69 Chapter 3. Snell Trace Estimation and Inversion 70 Figure 3.1: (a) Constant phase velocity ray path geometry for primary reflections in a horizontally layered earth model; with source location S and receiver positions Xi,X2,%3-(b) Primary reflection traveltime curves observed in an end-on shot gather. The Snell trace locates all energy corresponding to a fixed phase velocity (takeoff/emergence angle). Chapter 3. Snell Trace Estimation and Inversion 71 There has been a limited amount of previous work done on estimating Snell trajecto-ries. Gonzales and Claerbout (1984) extracted a Snell trace by wave equation migration of a linear moveout corrected common midpoint (CMP) gather. Their objective was a velocity analysis of the input CMP gather but with no emphasis on the problems of subsequent extraction of velocity-depth information. An advantage of their wave equa-tion based approach is that amplitude information can be accounted for in addition to traveltimes. Low and high amplitude primary reflective arrivals in the CMP gather can be focused onto a single Snell trajectory. Nevertheless, assessment of the final results is made difficult by a tortuous path to the final result. Retarded Snell midpoint coordinate transformation is followed by a 15° finite difference migration, and finally, focused energy maxima must be picked to define the desired Snell trajectory. Diebold and Stoffa (1981) describe a "r-sum" method of inverting reflection traveltime information from CMP gathers. They preferred to analyze CMP data because inherent averaging of takeoff and emergent ray parameters in this coordinate frame smoothed the effects of dip on the resulting velocity-depth inversions when carried out in the r — p domain. However, steeply dipping reflectors give rise to significant dip moveout (DMO) in a CMP gather, and thus, inverting Snell traces obtained from this seismic record results in biased velocity and depth estimates. In such cases, a DMO correction (Hale, 1984) should be applied prior to Snell trace estimation and inversion. I develop an algorithm to estimate and invert Snell traces from the local phase velocity image of a shot gather. My reasons for analyzing the shot gather instead of the CMP gather are twofold. Firstly, the CMP gather does not represent a physical experiment, and thus, hinders the evaluation of simplifying assumptions and coordinate transformations necessary to estimate the velocity function. Secondly, the shot gather often exhibits denser spatial sampling over the same total spread length. This leads to better resolution of local phase velocity and less spatial aliasing at the far offsets of shallower reflection Chapter 3. Snell Trace Estimation and Inversion 72 and refraction events. Estimation and inversion of Snell traces, from a local phase velocity image offers advantages and presents difficulties. An advantage is that there are many Snell traces that can be denned from a single shot gather (e.g., Figure 3.2a). Yet all trajectories arise from a unique velocity-depth distribution. Thus, multiple data sets are available in a shot gather for the inversion of a single v(z). This adds stabibty to the inverse problem. On the other hand, the inverse problem is dictated by difficulties in specifying the data. Snell paths from reverberatory events, ground roll, and noise can lead to first order uncertainties in determining a primary Snell trajectory. Therefore, the inversion algorithm must be flexible. Flexibibty and robustness are achieved through two mechanisms. First, the esti-mation and inversion procedure is interactive; decisions on Snell trace selection need to be immediately tested. For example, erroneous Snell trace estimates can be quickly identified by inversion and then forward modelling to compare calculated and observed traveltime curves. And second, the estimation of Snell traces can be made more robust by incorporating additional information. For example, picking Snell traces only from pre-determined primary reflective events maintains consistency while constraining the depth of reflector boundaries. In addition, constraints from semblance velocity analysis (Neidell and Taner, 1971) are helpful. A hyperbobc velocity spectrum can be used to constrain Snell trace picks where Snell information is most uncertain. Before outbning practical details of the inverse algorithm, I first derive the equations relating Snell trace information to both discrete and continuous velocity structures. Chapter 3. . Snell Trace Estimation and Inversion 73 3.2 Snell trace inversion formulae 3.2.1 Discrete formulation Relating subsurface velocity and depth information to a Snell trace from a shot gather is accompbshed by considering the ray path geometry shown in Figure 3.1. I assume a horizontally layered subsurface with constant velocity in each layer. Then, the total horizontal distance covered by a reflecting ray in the ith layer is given by Axi = ViAtisinOi, (3-1) and the vertical thickness can be expressed as Azi = ^—^cosdi, (3.2) 2 where Vi is the interval velocity, and Ati = ^ — ti_1 is the two-way traveltime in the zth layer. Along any ray path, Snell's law requires a constant horizontal phase velocity vp = Vi/sinOi, (3-3) which provides cosOi = ^-(vi/vpf. (3.4) Substituting relations (3.3) and (3.4) into equations (3.1) and (3.2) gives interval velocity and thickness in terms of the Snell trace parameters Axi, At{, and vp, that is, and, Vi = , (3.5) ( 3 6 ) The total depth to the ith boundary can be obtained by summing interval thicknesses if the complete Snell trajectory is known above this boundary ^ = ±Az: = l±AxJvp^-l. (3.7) Chapter 3. Snell Trace Estimation and Inversion 74 These formulae provide the basis of an algorithm for least squares inversion of interval velocity and depth described in Section 3.3. However, they are cumbersome for further development of important Snell trace concepts; thus, I present the following continuous formulation. 3.2.2 Continuous formulation I consider a continuous analogue of the above layered system to aid theoretical develop-ment of Snell trace concepts. First define a continuous Snell trajectory as a function of phase velocity and time x = s(vp,t). For a given v(z), each phase velocity describes an independent Snell curve in x — t space. Then, for a fixed phase velocity, the time derivative of the Snell trace becomes Ax; dsivjj.t) , ' , 2 ^ = * ( 3- 8 ) In the limit At{ —> 0, equations (3.5) and (3.6) become v(vp,t) = y/v^sJ^T), (3.9) and, dz(vp,t) = ^ f ^ - - 1, (3.10) 2 y s{vp,t) where, v(vp,t) represents instantaneous velocity along ray paths defined by vp. The total ray penetration depth is given by integrating (3.10) in time '(*p»o = £ fK^n.l^-T-idr (3.H) 2 Jo y s(vp,t') And thus, the 1-D subsurface velocity-depth function can be written v(z(vpit)) = ^Ups(vp,t). (3.12) Chapter 3. Snell Trace Estimation and Inversion 75 For each given phase velocity, there exists a corresponding Snell trace that gives rise to a common subsurface velocity-depth distribution down to the maximum depth penetration of the corresponding ray. Depth penetration of the Snell trace terminates when the critical angle is reached for a particular subsurface boundary. In this case, the Snell trace extends along the head wave branch of the traveltime curve for the bottoming boundary. Furthermore, a Snell trajectory terminates when a ray exceeds critical angle propagation, truncating deeper recovery of the velocity function. Equation (3.12) shows that velocity is determined directly by evaluating the time derivative of the Snell trace. Since differentiation is a roughening operation, the resulting velocity inversion will be sensitive to small fluctuations in the estimated Snell path2. Thus, constraining the variability of the estimated Snell trace is desirable. My algorithm constrains the Snell path to be linear between reflective events picked and is therefore vabd for a layered medium with constant velocity layers. Smoothness constraints are possible which allow velocity gradients between picks, but I do not consider them here. The reader is referred to Oldenburg et al (1984) for methods of recovering smooth velocity structures. On the other hand, equation (3.11) shows that depth information is determined by integrating a function which involves the time derivative of the Snell trace. Because integration is a smoothing operation, the resulting depth estimates given by (3.11) may be less sensitive to small fluctuations in the estimated Snell path. However, it is expected that smoothness constraints on the estimated Snell trace would provide a more stable depth recovery as well. 2 The square root operation in equation (3 .12) somewhat dampens the roughening caused by differentiation. Chapter 3. Snell Trace Estimation and Inversion 76 3.2.3 Example Figure 3.2a shows the positive phase velocity image of the synthetic shot gather in Figure 2.6a. I have picked four partial Snell traces spanning the eleven reflective events present in the original data. Table 3.1 shows the results of inverting three of those Snell traces using equations (3.5) and (3.6). Horizontal offsets and time picks along with the corresponding interval velocity and thickness estimates are given in that Table for each Snell trace. The true model parameters are listed in Table 3.3. And Figure 3.3a illustrates the comparison of the inverted velocity-depth models with the true subsurface velocity model. The shallowest reflection event was not picked in the Snell traces corresponding to phase velocities of 18000 and 13000 m/s because resolving power was lost at the near offsets. Thus, the corresponding velocities and thicknesses given for the second layer represent averages over the upper two layers. In addition, Snell trace penetration for the phase velocities of 13000 and 8000 m/s is limited by the spatial aperture of the shot gather and, hence, picks could not be made for the deeper reflectors. This example illustrates that Snell trace inversion produces good velocity and thick-ness estimates when the data are accurate and successive reflections are included in the analysis. Often, however, there are low amphtude reflections in field seismic data that cannot be rebably utibzed to define a Snell trajectory. Utilization of rms velocity con-cepts helps clarify the inversion in those cases when only major reflective events are used to define the Snell trace. This is discussed next. 3.2.4 rms velocity Along a given Snell trace, rms velocity is defined as VL.(vP,t) = ±£v\vp,t')dt\ (3.13) Chapter 3. Snell Trace Estimation and Inversion 77 Figure 3.2: (a) The positive half of the local phase velocity image shown in Figure 2.8b with several Snell trajectories drawn in. (b) Snell traces determined by skipping over reflection events 1,2,3,6,8, and 9. Chapter 3. Snell Trace Estimation and Inversion 78 vp = 18000 m/s VP = 13000 m/s vp = 8000 m/s R X t V Az X t V Az X t V Az 1 - - - - - - - - 60 .100 2194 106 2 88 .202 2792 279 99 .200 2538 249 158 .210 2666 138 3 126 .300 2671 130 161 .304 2792 142 263 . .317 2795 141 4 224 .404 4107 208 298 .408 4130 203 516 .436 4125 210 5 337 .603 3201 312 457 .614 3175 317 804 .652 3272 321 6 640 .912 4201 633 905 .933 4270 644 1567 1.003 4169 625 7 746 1.014 4307 213 1034 1.029 4164 190 1789 1.111 4061 188 8 971 1.126 6040 317 1365 1.148 6013 317 - - - -9 1100 1.230 4719 237 1544 1.258 4610 236 - - - -10 1357 1.434 4762 468 1933 1.477 4807 489 - - -11 1610 1.643 4662 472 - - - - - - - -Table 3.1: Snell trace parameters of horizontal offset x (rounded to the nearest meter) and traveltime t (rounded to the nearest milh-second) for three trajectories shown in Figure 3.2a ("R" denotes the reflector number). The inverted velocities v (m/s) and thicknesses Az (m) are listed for each layer. Compare these with the true model parameters given in Table 2.1. Substituting (3.9) into (3.13) gives a simple formula relating local rms and phase velocity V2n.(vp,t)=vps{vp,t)/t=vp^. (3.14) This formulation represents rms velocity along the refracted ray path. The continuous formulation is helpful in elucidating theoretical aspects of Snell traces. However, it is difficult to identify reflections from a subsurface velocity model that is con-tinuously variable. Only the major velocity and density contrasts will provide reflections that can be utibzed to define a Snell trace. Velocity derived from a Snell trace picked from only the major reflecting boundaries actually represents phase-velocity-dependent rms velocity between those reflective events picked. To prove this, consider two Snell trace picks corresponding to separate reflective events from the bottoms of layers i and where n > 1. Assume that reflection events from small impedance contrasts be-tween those layers provide unreliable Snell trace picks. The rms velocity over the time Chaptef 3. Snell Trace Estimation and Inversion 79 V E L O C I T Y M r o 0 4 ( j j - £ - - t ^ c n c n c r > o c n o c n o c n o c n o o o o o o o o o o o o o o o o o o o Figure 3.3: (a) Comparison of the velocity models obtained from inverting the three Snell traces in Figure 3.2a with the true velocity model shown by the dashed curve. interval t{ —> ti+n for a given phase velocity is defined as VLsiint} = 1 f t,+nv 2(vp,t')dr. W'+n *i J ti The velocity calculated between those picks using equation (3.5) is explicitly (3.15) (3.16) Differencing the rms velocity corresponding to the common phase velocity located on each traveltime curve provides ^ ™ . { * + »}-**V r™.M= v 2(vp,t')dt' (3.17) Then, utilizing equation (3.14) transforms this to vp(xi+n - Xi) = f t,+nv*(vp,t')dt~. (3.18) Chapter 3. Snell Trace Estimation and Inversion 80 Substituting (3.18) into (3.15) produces V L A ™ t } = v p X - ± ^ . (3.19) Correlation of this result with equation (3.16) establishes that the velocity calculated between Snell trace picks is actually a phase-velocity-dependent rms velocity over the interval. Thickness estimates are correspondingly affected by this velocity averaging between Snell trace picks. For example, Table 3.2 shows the results of inverting three Snell traces shown in Figure 3.2b. Reflectors 1, 2, 3, 6, 8, and 9 were skipped to simulate the effects of Snell trace inversion when reflection events are missed, a definite possibility with field data sets. Horizontal offsets and traveltimes are given along with corresponding interval velocity estimates for each Snell trace. The true interval rms velocities are also listed for comparison. Figure 3.4 illustrates the comparison of inverted velocity-depth models with the true velocity-depth model. As indicated, this procedure produces angle dependent interval rms velocity estimates and the inverted depths are somewhat biased from the true depths of the reflectors picked. Nevertheless, "interval rms" velocity estimates are useful subsurface parameters if one keeps in mind their relation to the true interval velocity model-Synthetic examples show that several Snell traces from a single shot gather can be inverted to obtain an estimate of v(z). However, such examples also reveal that it is not always possible to define a complete Snell curve over the desired range of primary subsurface reflectors. At small source-receiver offsets, lack of resolution and high sen-sitivity to picking errors make Snell trace estimation at near surface reflectors difficult for larger phase velocities (smaller incident angles). In addition, post critical reflection and finite spread aperture restrict the depth penetration of measurable Snell trajectories, Chapter 3. Snell Trace Estimation and Inversion 81 V P = 18000m/s vp = 13000m/s vp = 8000m/s R X t V X t V "rms X t V "rmi 4 208 .404 3045 3006 294 .408 3060 3013 516 .436 3075 3039 5 345 .610 3452 3291 461 .616 3237 3291 808 .650 3309 3291 7 746 1.014 4226 4241 1045 1.035 4255 4241 1809 1.114 4151 4242 10 1341 1.439 5023 5149 1922 1.473 5098 5157 - - -11 1614 1.640 4950 4842 - - - - - - -Table 3.2: Snell trace parameters of horizontal offset x (m) and traveltime t (s) for three Snell trajectories shown in Figure 3.2b. Reflectors 1,2,3,6,8, and 9 were not included (skipped) in the Snell trace picks. The inverted velocities v (m/s) are compared with theoretical angle dependent interval rms velocities Vrm, (m/s) calculated by numerical evaluation of equation 3.15. leaving deeper reflectors unconstrained at smaller phase velocities (higher incident an-gles). Therefore, in practice, only portions of Snell traces are rebably obtainable from a field shot gather. However, if overlapping Snell trace portions that sufficiently constrain the stratified medium illuminated by the shot are estimated, then the procedure outbned in the next Section can be utibzed to invert for a least squares solution to v(z). 3.3 Least squares algorithm The objective of this section is to describe a simple algorithm for the least squares in-version of constraints provided by several partial Snell trace estimates. I assume that the subsurface consists of a relatively few significant formations that give rise to distinct primary reflections. The total number of reflections utibzed is initially specified. Then, partial Snell traces are picked from the predetermined reflections by locating a common phase velocity on as many events as possible. Many Snell curves can be estimated (there are an infinite number of them), however, each must be consistent with the same sub-surface velocity-depth distribution. Each Snell trace corresponding to primary reflection Chapter 3. Snell Trace Estimation and Inversion 82 VELOCITY ro C-J -r^ -p- cn cn CD O cn o cn o cn O cn O O o o o o O O O O O o o o o O O O o Figure 3.4: Comparison of the velocity models obtained from inverting three Snell traces in Figure 3.2b with the true velocity model. Reflectors 1,2,3,6,8, and 9 were bypassed in the Snell trace definition resulting in interval rms velocities between the surface and reflector 4, reflectors 5 and 7, and reflectors 7 and 10. Chapter 3. Snell Trace Estimation and Inversion 83 energy provides constraints on local subsurface velocity and thickness. Velocity between Snell trace picks is constrained by equation (3.5) and thickness is constrained by equation (3.6). These constraints provide a system of linear equations Ap = 6, (3.20) where, p = [vi, v2, • • •, vn, Z\,z2r- • •, zn]T is the unknown parameter vector, b is the known data vector, and A is a sparse system matrix. As an example, system entries for a single Snell trace corresponding to a three layered model are simply 1 0 0 0 0 0. 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 - 1 1 0 0 0 0 0 -1 1 A = P = [Vl,V2,V3,Z1}Z2,ZS\ , and, b = y/vyAxi/At! yJvpAx2/At2 ^/vpAx3/At3 lAx^VpAtJAa?i - 1 \Ax2^JvpAt2lAx2 - 1 \Ax3^vpAt3/Ax3 - 1 where the data vector 6 consists entirely of Snell trace parameters extracted from the local phase velocity image of a shot gather. This linear system becomes overdetermined when many partial Snell traces are in-cluded that provide overlapping constraints. Singular Value Decomposition (Lines and Chapter 3. Snell Trace Estimation and Inversion 84 Treitel, 1984) provides a stable least squares inverse to this problem. Furthermore, lin-ear velocity and depth contraints from other analyses could be included in this system (Oldenburg et al, 1984), in addition to a weighting matrix that provides statistical nor-mabzation and parameter emphasis (Jackson 1972, and Oldenburg 1984). But I do not consider the effects of those additional constraints here; my purpose is to outbne the importance of Snell trace constraints. Further discussion on the solution of discrete bn-ear inverse problems is given by Menke (1984) and treatment of the philosophies and practicalities of solving general inverse problems is provided by Tarantola (1987). As an example, I have simultaneously inverted the four Snell traces in Figure 3.2a. Table 3.3 compares the least squares output to the true model parameters while Figure 3.5 illustrates a good match between the constructed and the true velocity model. I note that in general, the depth parameter is more accurately recovered than the interval velocity in synthetic tests. This is presumeably due to the relative stabibty of the respective inverse relations (3.11) and (3.12) to errors in Snell trace estimates as suggested previously. 3.4 Inversion of Field Data Estimation and inversion of Snell traces from field data are more difficult. Multiple valued phase velocities along a single reflection event can arise from crossing events, and random and systematic noise as discussed in Chapter 2. In addition, first order surface multiples give rise to radial SneU traces and pegleg multiple Snell paths branch off bnearly from the primary Snell trace (Gonzales and Claerbout, 1984), thus hindering the accurate estimation of primary Snell traces. For example, Figure 3.6 shows a marine shot gather that exhibits significant water bottom reverberations. These data were obtained offshore eastern Canada in water depths around 150 m. Figure 3.6a shows a raw shot record and Figure 3.6b shows Chapter 3. Snell Trace Estimation and Inversion 85 True Model. Inverted Model R V Az z V z 1 2000 100 100 2194 106 2 2742 137 237 2666 243 3 2689 133 370 2753 381 4 4153 207 577 4121 588 5 3291 329 906 3216 905 6 4187 628 1534 4203 1538 7 4399 220 1754 4323 1743 8 6141 307 2061 5982 2059 9 4969 248 2309 4680 2293 10 4655 465 2774 4724 2763 11 4842 484 3258 4852 3252 Table 3.3: True model parameters on the left. On the right are the velocity (m/s) and depth (m) estimates obtained by least squares inversion of all Snell traces shown in Figure 3.2a. V E L O C I T Y Figure 3.5: Least squares inversion of all Snell traces in Figure 3.2a (solid curve) compared with the true velocity model (dashed curve). Chapter 3. Snell Trace Estimation and Inversion 86 the record after velocity filtering has been applied to better isolate reflective arrivals. Figure 3.6c shows the phase velocity image of 3.6b where only the narrow band of phase velocities from 4000-5000 m/s is shown. Multiple events and noise broaden this phase velocity window, making it more difficult to identify a primary Snell path. I have drawn in the radial Snell trajectory (for vp = 4500 m/s) corresponding to first order water bottom multiples. Nevertheless, for a subsurface model that exhibits general velocity increase with depth, the primary Snell trace can often be distinguished along the upper edge of the phase velocity window. In this case, the time derivative of a primary Snell path increases with time causing the trajectory to curve upward with increasing offset, while the Snell paths from reverberations maintain briear trajectories. For example, a primary Snell trace corresponding to 4000 m/s might be defined by following a trajectory near the upper phase velocity bound shown in Figure 3.6c. A note of caution is necessary if a low velocity zone is encountered. Then, the time gradient of the primary Snell path decreases, possibly causing interference of primary and multiple Snell paths. In those cases, it may be necessary to first remove the contaminating multiple events (e.g., Wiggins, 1988). As an example of inverting field data, I choose a marine shot gather that does not contain significant multiple contamination. Figure 3.7a shows a field shot gather obtained in deep (2.5 km) water offshore Vancouver Island, British Columbia. The time window of data shown does not include the first water bottom multiple. For the purposes of Snell trace estimation and inversion presented at the end of this section, I assume that other reverberatory events are insignificant in this data window as well. The first step in preparing these data for Snell trace estimation is the appbcation of a velocity filter to remove contaminating negative phase velocity noise. Figure 3.7b shows the field gather after appbcation of a velocity filter with pass range of 2000 m/s to infinite phase velocity. Figure 3.8a shows phase velocity output using the raw shot gather, and Chapter 3. Snell Trace Estimation and Inversion 87 Phase V e l o c i t y Figure 3.6: (a) Shallow water marine shot record from offshore eastern Canada, (b) Shot record after application of velocity filter, (c) phase velocity image corresponding to (b) showing phase velocities between 4000 and 5000 m/s only. A radial trace (with vp = 4500 m/s) corresponding to the water bottom reverberation is shown along with an outer bound on the primary Snell trace. Chapter 3. Snell Trace Estimation and Inversion 88 TSN 30 GO 90 120 150 180 210 240 T5N (a) (b) Figure 3.7: (a) Shot gather from offshore western Canada; data window shown is free of water bottom reverberations, (b) Shot gather after appbcation of a velocity filter to remove contamination from negative phase velocity noise. Figure 3.8b shows phase velocity output using the velocity filtered shot gather. The image obtained after velocity filtering shows an overaU signal-to-noise enhancement. Further signal-to-noise enhancement is possible by stacking the phase velocity images of adjacent shot records. Stacking techniques tend to reduce the effects of random noise and the effects of small scale inhomogeneities in subsurface structure. More sophisticated image enhancement techniques, using adjacent shot records, could be utibzed to improve phase velocity estimates. Principal component analysis (Ready and Wintz, 1973) is one such approach. Once an optimum image of local phase velocity is available, Snell trace estimation and inversion is most efficiently carried out on an interactive workstation. Figure 3.9a illustrates a graphics screen from the program INSTEIN (for LNteractive Snell Trace Estimation and INversion) that I have developed on a Micro V A X II workstation. The Chapter 3. Snell Trace Estimation and Inversion 89 Phase V e l o c i t y Figure 3.8: On the left; phase velocity image obtained from the raw shot record. On the right; phase velocity image obtained from the velocity filtered shot record. primary image in that figure is the local phase velocity data obtained by stacking the phase velocity images of four adjacent velocity filtered shot gathers of the Vancouver Island data. One raw shot gather is shown in wiggle trace form in the center of that Figure. On the far left is a hyperbolic velocity spectrum obtained by semblance analysis of the input shot gather. And on the bottom right is an area reserved for the velocity-depth functions obtained by Snell trace inversion. As an example, four partial Snell traces were picked covering seven reflectors and inverted to obtain the velocity-depth model shown. Table 3.4 lists the velocity and depth parameters obtained from the least squares inversion. The important features of the interactive program INSTEIN are listed below. 1. Snell trace picking: (a) picking Snell traces from a local phase velocity image while utilizing constraints Chapter 3. Snell Trace Estimation and Inversion 90 Least Squares Solution R v (m/s) z (m) 1 1486 2529 2 1750 2713 3 2050 3157 4 2166 3348 5 2337 3759 6 2304 3948 7 3000 4583 Table 3.4: Velocity v and depth z obtained from a least squares inversion of the four Snell traces shown in Figure 3.9a. from the original shot data and a hyperbobc velocity spectrum, (b) isolation of individual phase velocity windows to aid picking (see Figure 3.9b). 2. Inversion of individual Snell traces (see Figure 3.9b). 3. Least squares inversion of several Snell traces. 4. Forward modelling of inverted velocity models: (a) forward Snell trace modelling, (b) forward reflection traveltime modelling, (c) forward headwave traveltime modelbng. As indicated in 1(a) above, SneU trace picks are constrained by two additional cursors; one located on the wiggle trace data, and the other located on the hyperbobc velocity spectrum by utihzing V ^ m s = <Jvpx/t, and t0 = \Jt2 — {x/Vrms )2- Each observed phase velocity defines an rms velocity; and I use the hyperbolic approximation of normal in-cidence reflection, time t0 to locate a particular pick in the velocity spectrum3. It is 3 Hyperbolic velocity maxima provide an approximation to VTrrl, as shown by Taner and Koehler (1969). Chapter 3. Snell Trace Estimation and Inversion 91 COLtMPJ BGSUH ADVEL IPLTVEL BGPIK ADPIK INVERT 1 2 3 4 5 6 SNELL V(Z) EXIT ' Figure 3.9: (a) Graphics screen from the Snell trace estimation and inversion program IN-STEIN. On the right is a local phase velocity image corresponding to the shot data in the center. Four Snell trajectories were estimated and inverted to obtain the velocity-depth model shown at the bottom right. Snell trace picks are constrained by the hyperbobc velocity spectrum shown at the far left. Modelled traveltime curves are superimposed on the shot data to determine goodness of fit. (b) Example of phase velocity isolation to aid Snell trace picking. Chapter 3. Snell Trace Estimation and Inversion 92 important to make picks on the same reflection each time, and it is often desirable for each pick to maximize hyperbolic semblance information (particularly when picking at near offsets where Snell information is weakest). Ideally, the inversion could be driven by maximizing semblance along the exact trav-eltime curves given by the forward parametric equations However, the hyperbolic approximation provides adequate constraints with less compu-tational effort. Many algorithms are available for rapidly computing hyperbolic velocity spectra from seismic reflection data (Neidell and Taner, 1971). Nevertheless, I use the parametric equations to determine the goodness of fit between the observed reflection traveltime curves and those provided by the inverted results. Figure 3.9a shows the modelled reflection traveltime curves superimposed on the observed data. There is good agreement between modelled and observed traveltime curves. In addition, head wave traveltime branches can be modelled and tested for goodness of fit with those observed in field data (head waves are not present in this deep water data). A good match in both cases strengthens the acceptability of the final velocity-depth model. Snell trace theory provides a very physical means of approaching 1-D velocity in-version. Utilization of these concepts leads to an efficient program (e.g. INSTEIN) for estimating and inverting Snell traces from the local phase velocity image of a shot gather. It is therefore worthwhile investigating the possible extension of Snell trace concepts to laterally heterogeneous media. I next discuss this extension when the subsurface struc-ture contains planar dips. and, Chapter 3. Snell Trace Estimation and Inversion 93 3.5 Plane dipping layers Snell trace theory loses much of its physical significance when the subsurface velocity model is laterally inhomogeneous. Snell's law requires a variable horizontal phase veloc-ity along the ray path in a plane dipping model, and thus a single Snell trace no longer accurately constrains layer velocity and thickness. Downgoing ray paths of a common emergent phase velocity are not necessarily coincident. However, I show that for small reflector dips, averaging Snell traces of a common phase velocity from forward and re-versed shot gathers can provide approximate velocity and thickness information about the subsurface model vertically beneath the midpoint of the two shots. The accuracy of this approximation depends on the magnitude of subsurface dips and on the emergence angle of the average Snell trace. Errors increase with reflector dip and with emergence angle. Acceptable errors must be defined by the user. I show that a rough quantification of expected errors is possible by simple forward modelling tests, if a rough estimate of subsurface dips is given. The traveltime relation for a 2-D plane dipping model of n layers has been given by Johnson (1976) *(«P) = - ^ + E -{coscti+cospi), (3.21 where a and (3 are ray angles referenced to a vertical axis (see Figure 3.10), i>; is the interval velocity, Azi is the thickness of the ith layer vertically beneath the shotpoint, and vp is the horizontal phase velocity of the emergent energy at the surface. This relation, which is closed form in terms of the linear headwave traveltime branches, has also been verified for reflection traveltimes by Diebold and Stoffa (1981). However, forward modelling for the nonlinear reflection traveltime curves must involve ray tracing techniques (e.g., Hubral and Krey, 1980). Johnson (1976) describes a recursive algorithm that attempts to unravel the velocity Chaptei 3. Snell Trace Estimation and Inversion 94 vi = 2000m/a u2 = 3000m/s v3 = 2000m/« u4 = 3000m/« Figure 3.10: Forward and reverse ray path geometry for a plane dipping layered medium when the emergent angle for both raypaths is identically rjj. and thickness of dipping subsurface structure by using headwave data from forward and reversed shot gathers. Unfortunately, velocity reversals are undetectable using headwave data. It is conceivable that a layer stripping operation could be designed that utibzes local phase velocity data from forward and reversed reflection data. Such an approach could then recover low velocity zones. However, the recursive nature of layer stripping operations leads to cumulative errors. Uncertainty in the recovered model parameters increases as the layer stripping proceeds. As an alternative, I investigate ah approximate solution obtained by simply averaging Snell traces of a common emergent phase velocity from forward and reversed shot gathers. It is my experience that errors are not necessarily cumulative using this approach. Averaging traveltime relations with a common emergent phase velocity for a spbt Chapter 3. Snell Trace Estimation and Inversion 95 spread shot gather (see Figure 3.10) gives tf + tT Xf + xr 1 J ^ Az; — - — = — h - V (coscti + cos(3i + cosvi + costi). (3.22) This result can be simply written as tove= — +2'£^cos8i, (3.23) vP ~[ ^ where, cosdi = (coscti + c o s / ? i + cos£i + cosj] i)/A. (3-24) In Appendix D, I show that equation (3.23) can be approximated by x 1 n Az-* = - + - E —-{cos{9i - u-i) + cos(8i -ei) + cos(6i + pi) + cos{6i + e i)}, (3.25) vp 2. = 1 «i where e and fi are perturbation angles described in Appendix D, and Q{ is the refraction angle in the ith layer corresponding to a horizontal subsurface model with velocity-depth model parameters vertically beneath the shotpoint. When the angles e and /z are small, cos9 — cos9, and thus, equation (3.23) provides a rough approximation to the exact horizontal traveltime relation . , x(vv) Azv „ . t(vp) = + 2T —-cos9i. 3.26 VP i = i vi This suggests that an average Snell trace formed from t a v e and xave, for sequential dipping reflectors, may provide a good approximation to the Snell trace corresponding to an equivalent horizontal model; and that subsequent inversion of an average Snell trace may provide a good estimate of the true vertical velocity model. As an example, Table 3.5 shows ray traced traveltimes and offsets of common-emergent-angle, forward and reversed, reflections from the four boundaries shown in Figure 3.10. An emergent angle of ip = 20° was used giving a phase velocity of 5848 m/s. This table also compares the average Snell trace parameters with Snell trace parameters from an Chapter'3. Snell Trace Estimation and Inversion 96 equivalent horizontal model containing the same vertical depth and velocity parameters. Figure 3.11 compares velocity functions obtained by inverting the Snell traces bsted in Table 3.5. Inversion of the average Snell trace data results in velocity and thickness errors that are less than 1%. For completeness, Table 3.6 bsts the ray angles for each reflector and compares the cosine average angles 9{ to the equivalent horizontal angles 6{. As shown, the correspondence between 6i and 9{ worsens when the angle or /X; is large for a particular layer. In fact, e{ and pi are dependent on refraction effects along both forward and reverse ray paths, and thus, errors depend on the configuration of the subsurface model. The above example illustrates the utibty of averaging forward and reversed Snell traces to reduce the effects of dip on subsequent Snell inversion. However, errors remain that are dependent upon the size of layer dip and on the angle of emergence chosen in addition to the configuration of subsurface dips. For example, Table 3.7 bsts interval velocity and thickness obtained from averaged Snell traces as dip angles are increased for the reflectors shown in Figure 3.10. Reflector 2 was held horizontal while the dip of reflector 1 was increased positively and reflectors 3 and A were increased negatively according to the parameter 8. And Figure 3.12 provides illustration of these errors in the resulting velocity-depth models. For dips of 6° and 10°, the interval velocity and thickness errors are roughly 1% and 5% respectively4. Furthermore, Table 3.8 bsts interval velocity and thickness obtained from averaged Snell traces as the emergence angle is increased. In general, smaller emergence angles produce less error as verified in that Table. A further example emphasizes the utibty of this approximation. Figure 3.13a shows a dipping subsurface model (model dips range from 0 to about 3.5°) with several ray paths eminating from a source at 5 km. Primary reflection traveltime curves, appropriate for 4 I t is interesting to note that recovery of velocity in the upper layer is unaffected by dip yet the thickness estimate is biased. A more accurate depth estimate for the upper layer could be obtained by using 3.21 if gathers are available that provide reciprocal source and receiver locations. Chapter 3. Snell Trace Estimation and Inversion 97 Ray Traced Snell Traces Horizontal Snell Trace Vertical Depth Reflector Dip R tf tT tave X f %T Xave t X z 8 1 .528 .528 .528 274 448 361 .532 364 500 5° 2 .914 .930 .922 843 1120 981 .920 962 1000 0° 3 1.443 1.450 1.447 1517 1174 1345 1.452 1326 1500 -5° 4 1.830 1.837 1.833 2211 1693 1952 1.841 1923 2000 . -5° Table 3.5: Snell trace parameters of offset and traveltime computed by ray tracing forward and reversed rays of the plane dipping model shown in Figure 3.10. The forward and reversed averages are compared to the Snell trace parameters of an equivalent horizontal model. The common emergent angle for these results is ip = 20°. The model parameters of vertical depth and dip refer to the model shown in Figure 3.10. Velocity Figure 3.11: Velocity models obtained by inverting Snell traces from, (a) the forward Snell traces, (b) the reverse SneU traces given in Table 3.5 above. Also, the velocity model (middle sobd curve) obtained from the average Snell trace is compared with the true vertical velocity-depth model (dashed curve). Chapter 3. Snell Trace Estimation and Inversion 98 Vertical Angles Vertical Angles Dipping Model Horizontal Model R i Pi Vi 6 8i Bi 1 1 10.00 20.00 30.00 20.00 21.16 20.00 R i A Vi 6 Bi 2 1 16.20 20.00 24.07 20.00 20.25 20.00 - 2 27.84 27.84 34.34 34.34 31.24 30.87 R i Pi Vi ti Bi Bi 3 1 25.73 20.00 13.83 20.00 20.31 20.00 - 2 45.03 27.84 18.31 34.34 32.70 30.87 3 28.14 18.14 12.09 22.09 20.9.1 20.00 R i «t Pi Vi ti- Bi 4 1 25.73 20.00 13.83 20.00 20.31 20.00 - 2 45.03 27.84 18.31 34.34 32.70 30.87 - 3 28.14 18.14 12.09 22.09 20.91 20.00 - 4 41.12 31.12 21.16 31.16 31.85 30.87 Table 3.6: Ray traced vertical angles for each reflector shown in Figure 3.10. Averaging the cosines of these angles for each boundary results in the angles 0{, that are compared to the vertical angles obtained from a horizontal model with the vertical depths given by the same model. Velocity an d Thickness Variations with Dip 8 = 0° 8 = 2° 8 - 4° 8 = 6° 8° 8 = 10° R v Az v Az V Az v Az V Az V Az 1 2000 500 2000 499 2000 498 2000 495. 2000 490 2000 485 2 3000 500 3005 502 3022 507 3051 516 3093 529 3153 547 .3 2000 500 2002 499 2008 497 2022 495 2055 495 2146 507 4 3000 500 3005 500 3018 501 3041 502 3074 503 3119 506 Table 3.7: Interval velocity and thickness estimates from averaged Snell traces when the angles (8) of the dipping layers shown in Figure 3.10 increase. Velocity and depth errors increase with dip. Chapter 3. Snell Trace Estimation and Inversion 99 o o o Velocity -» K> K> C-J C-4 > cn o cn o en o o o o o o o o o o o o o — 1 1 1 1 1 cn o o o CD o o o o o NO o o o Figure 3.12: Velocity models obtained from the data given in Table 3.7. Velocity and Thickness Variations with Emerg ence Angle ip = 5° 4> = io° Tp = 15° i> = 20° ip = 25° V = 30° R S v Az v Az v Az v Az v Az V Az 1 5° 2000 496 2000 496 2000 496 2000 496 2000 496 2000 496 2 0° 3020 508 3022 508 3027 509 3035 511 3048 514 3073 523 3 -5° 1998 492 2001 492 2005 494 2014 496 2034 502 2102 524 4 --5° 3020 499 3021 500 3024 500 3028 501 3035 503 3047 507 Table 3.8: Interval velocity and thickness estimates from averaged Snell traces when the emergence angles (ip) increase. Velocity and depth errors increase with emergence angle. Chapter 3. Snell Trace Estimation and Inversion 100 a split-spread experiment, are shown in Figure 3.13b. Synthetic shot gathers were then generated from ray traced reflection times at selected source locations for the purpose of Snell trace estimation and inversion. Figures 3.14a and 3.14b show the results of inverting the forward and reversed gathers separately for source locations at 3, 5, and 7 km. Snell trace estimation and least squares inversion was carried out with INSTEIN using local phase velocity images of the synthetic shot gathers. As shown, estimated velocity functions (sobd curves) are in error when compared with the true models (dashed curves), confirming that Snell trace inversion is highly sensitive to dip. Figure 3.14c shows the results of inverting averaged Snell traces formed from each spbt-spread shot gather. The resulting velocity-depth functions are within 1% of the true subsurface functions. This example strengthens the premise that (3.23) represents an adequate approximation to (3.26) for the purposes of Snell trace inversion when reflector dips are small. Empirical studies also indicate that the utibty of averaging forward and reversed Snell traces of a common phase velocity is not restricted to the spbt-spread experiment. Any two reversed shot gathers may be used to produce averaged Snell trace information. The subsequent inversion provides approximate velocity-depth information vertically beneath the midpoint of the two shotpoints. For example, Figure 3.15 shows the results of in-verting averaged SneU traces from reversed gathers at shotpoint pairs from the model in Figure 3.13. The shotpoint pairs are located at (1,5 km), (3,7 km), and (5,9 km), where the first and second distances are the forward and reversed shot locations respectively. Again the velocity and thickness errors are less than 1%. The midpoints of the above pairs correspond to the split-spread shot locations above. Hence, corresponding split-spread Snell traces can be averaged with the above Snell traces to improve the statistical measure. Figure 3.15b illustrates the results of inverting a Snell trace average containing two Snell trace pairs referenced to each midpoint. - In fact, the typical spbt-spread seismic reflection survey provides multi-fold, forward and Chapter 3. Snell Trace Estimation and Inversion 101 DISTANCE -(KM) Q 0 2 4 6 8 10 DISTANCE (KM) Figure 3.13: (a) Plane dipping layered test model showing some ray paths arising from a split spread experiment, (b) Resulting traveltime curves from the primary reflections only. Chapter 3. Snell Trace Estimation and Inversion 102 Figure 3.14: (a) Velocity-depth models obtained by inverting forward Snell traces at shot locations of 3, 5, and 7km. (b) Velocity-depth models obtained by inverting reversed Snell traces at the same shot locations, (c) Velocity-depth models obtained by inverting averaged Snell traces from forward and reverse shot data. The dashed curves represent the true vertical velocity-depth models. Figure 3.15: (a) Velocity-depth models obtained by inverting average Snell information using reversed gathers at (1,5 km), (3,7 km), and (5,9 km) from Figure 3.5, where the first and second distances are the forward and reversed shot locations respectively, (b) Inverted results when the corresponding split-spread Snell traces are included in the Snell trace average. In (a) and (b), the dashed curves represent the true vertical velocity-depth models for the midpoint locations of the shotpoint pairs. Chapter 3. Snell Trace Estimation and Inversion 103 reverse, subsurface coverage. This recording geometry allows definition of several Snell trace pairs that can be referenced to a specific midpoint location. As a result, Snell trace estimation could be designed to facilitate common midpoint. Snell trace stacking before inversion. This is a subject of further investigation. Exact quantification of errors obtained by inverting averaged Snell traces is difficult since errors are model dependent. The amount of tolerable error will depend upon the application and on user specifications. When some idea of subsurface dip is known, a rough quantification of expected, errors can be obtained by simple forward modelling tests such as I have presented. In any case, errors due to dip should ,be kept in mind when interpreting velocity models obtained from averaged Snell traces. Nevertheless, reasonably accurate velocity and depth inversions can be obtained for a plane dipping subsurface when the dips are sufficiently small. When the dips are more severe, I speculate on possible improvements in the shot record inversion by first applying a DMO correction. Documented improvements in CMP velocity analysis after such a correction (Hale, 1984) suggests that improvement may also be achieved in the analysis of Snell traces. I have not investigated this possibility. 3.6 Conclusion -I have developed a procedure to estimate and invert Snell traces from the local phase velocity image of a shot gather. The Snell trace corresponding to primary reflection events has a significant physical interpretation when the subsurface being explored is horizontally layered. Such a trajectory locates all primary reflection energy corresponding to a fixed takeoff/emergence angle. Snell trace estimates provide data for the inversion of velocity and depth in a layered earth system; interval velocity and layer thickness can be expressed in terms of the time derivative of a Snell trace. Chapter 3. Snell Trace Estimation and Inversion 104 Estimation and inversion of Snell traces from a local phase velocity image offers advantages and presents difficulties. An advantage is that there are many Snell traces that can be denned from a single shot gather. Yet each trajectory arises from a unique velocity-depth distribution. Thus, multiple data sets are available in a shot gather for the inversion of a single v(z). This adds stability to the inverse problem. On the other hand, the inverse problem is dictated by difficulties in specifying the data. Snell paths from reverberatory events, ground roll, and noise can lead to first order uncertainties in determining a primary Snell trajectory. Therefore, the estimation is constrained by velocity semblance data and a priori reflector information. Finally, Snell trace inversion can be extended to a planar dipping subsurface when the dips are small. The estimation and inversion is most efficiently carried out on an interactive workstation in real time. Chapter 4 Discussion This chapter is reserved for speculations on possible research topics stemming from the ideas put forth herein. It is not meant to be a rigorous outline complete with references of support. As a result, suggestions and questions provided below should be taken with a grain of salt. 4.1 On the time-variant filter I have done a bmited amount of research into the possibibty of a space-time-variant filter using the principles outbned in Chapter 1. Calculation of a 2-D phase surface corresponding to input frequency and wavenumber cutoff functions proved to be non-trivial. The problem is analogous to constructing a surface from a grid of given surface gradients. One problem is that the solution is nonunique. Another problem is that any proposed solution appears to be computationally intensive. It is not clear whether further research in this area would lead to the production of an efficient filter. On the plus side, this research led to development of the theory for Chapter 2. It would be interesting to further investigate the frequency-variant window and pos-sible appbcations. One possible appbcation is a time-variant notch filter. This could be used to remove a quasi-harmonic noise component that slowly varies in frequency with time. Real time (or close to real time) appbcation of the time-variant filter may be possible. Complex signal computation and convolutions with the complex delta function could be 105 Chapter 4. Discussion 106 hard coded for seismic applications. Desired frequency cutoff functions (translating to phase functions in the complex exponential factors) could then be programmed by the field operator for each specific time-variant filter appbcation. However, it is hard to think of a good reason to. pre-time-variant filter raw seismic data. In general, it is safer to recover raw data in the field and filter in the office. Nevertheless, data that require a time-variable digitization interval would be one motivation for a pre-time-variant filter. There may be a need for a high quabty time-variable anti-abasing filter. Finally, given the speed and storage capacity of today's computers, it would be worth-while to perform benchmark tests on the relative efficiency of my approximate filter and integration of the full blown time-variant convolution integral. 4.2 On the calculation of local phase velocity One likely appbcation of local phase velocity calculation is spatial interpolation. Mea-sured seismic records are sometimes spatially abased (SA). Spatial interpolation reduces SA artifacts produced by processing algorithms such as migration. Local slopes of seis-mic events are necessary parameters for the interpolation. A local phase velocity image calculated from data that have been low-pass frequency filtered could provide that infor-mation. Low-pass filtering reduces resolution of local slope but eliminates phase velocity bias caused by spatial abasing. I feel more research is warranted on the numerical evaluation of instantaneous fre-quency and wavenumber attributes, particularly in those areas where the estimation becomes unstable. If the recovery of local frequency and wavenumber can be viewed as a formal inverse problem, then the methods of inverse theory may be used to characterize and perhaps stabilize those problems. In addition, the formulation of a directional phase velocity estimator, one that can be tuned to a specified range of dips, would be a great Chapter 4. Discussion 107 benefit. This could obviate the application of velocity filters to stabilize the calculations in the presence of noise. Dispersion analysis of surface waves is another possible application of local phase velocity calculation. The calculation of local phase velocity images at a succession of narrow frequency bands could expose a variation of phase velocity with frequency. Of course, subsurface earth parameters could then be inverted from the recovered phase velocity vs frequency curves. Considering the complex data formulation of local phase velocity using the Fourier transform; could a similar local slope calculation algorithm be formulated using other orthogonal transforms such as the Walsh transform (Wood, 1974)? Is the complex for-mulation avoidable using an alternate to the Fourier transform? Local slope calculation may prove useful in artificial vision. Evaluation of the envi-ronment is determined to a large extent by patterns in the variation of light intensity. Edge detection along with local slope estimation could be used as a basis for pattern recognition. What about 3-D target location using a spherical radar source and the principles of local phase velocity estimation? Is real time estimation of local phase velocity possible? These questions are probably worth answering. I have done limited work on using local phase velocity data for calculation of Fresnel zone widths and for constant velocity map migration. The latter application is more of a curiosity than a practical method of migration. Inverse ray tracing techniques may be a more beneficial use of local phase velocity information for the purpose of depth migration. In addition, T-p information can be obtained at each point in a shot record without resorting to a full T-p transform. Thus, many data are available for input to traditional velocity inversion algorithms that utilize T-p data. I feel that interactive data evaluation Chapter 4. Discussion 108 is an essential element in any inversion scheme utilizing local phase velocity data. The large amount of information available is most efficiently handled using the tools of a workstation. Bibliography [1] Al-Chalabi, M. , 1,974, An. analysis of stacking, rms, average, and interval velocities of a horizontally layered ground: Geop. Prosp., 22, 458-475. [2] Bardan, V., 1987, Trace interpolation in seismic data processing: Geoph. Prosp., 35, 345-348. [3] Bednar, J.B., 1983, Appbcations of median filtering to deconvolution, pulse estima-tion, and statistical editing of seismic data: Geophysics, 48, 1598-1610. [4] Bickel, S.H., and Natarajan, R.R., 1985, Plane-wave Q deconvolution: Geophysics, 50, 1426-1439. [5] Boashash, B. and Whitehouse, H.J., 1986, Seismic appbcations of the Wigner-Ville distribution: Proceedings of the 1986 IEEE International Circuits and Systems Sym-posium, San Jose, USA, 1, 34-37. [6] Bracewell, R.N., 1978, The Fourier Transform and its Appbcations: New York, McGraw-Hill Book Co., Inc. [7] Cain, G.D., 1972, Hilbert-transform description of bnear filtering: Electron. Lett., 8, 380-382. [8] Chun, J.H. and Jacewitz, C.A., 1981, Fundamentals of frequency domain migration: Geophysics, 46, 717-733. [9] Claerbout, J.F., 1976, Fundamentals of data processing: New York, McGraw-Hill Book Co. 109 Bibliography 110 [10] Claerbout, J.F., 1985, Imaging the earth's interior: Blackwell Scientific Pub. [11] Clarke, G.K.C., 1969, Time-varying deconvolution filters: Geophysics, 33, 936-944. [12] Diebold, J, and, Stoffa, P, 1981, The traveltime equation, tau-p mapping and inver-sion of common midpoint data: Geophysics, 46, 238-254. [13] Gabor, D., 1946, Theory of communication: J. IEE(London) 93, 429-457. [14] Gonzalez-Serrano, A., and Claerbout, J.F.,1984, Wave equation velocity analysis: Geophysics, 49, 1432-1456. [15] Grant, F.S. and West, G.F., 1965, Interpretation theory in applied geophysics: McGraw-Hill Book Co. [16] Griffiths, L.J., Smolka, F.R., and Trembley, L.D., 1977, Adaptive deconvolution: A new'technique for processing time-varying seismic data: Geophysics, 42, 742-759. [17] Hale, D.j 1984, Dip-moveout by Fourier transform: Geophysics, 49, 741-757.. [18] Hubral, P., and Krey, T., 1980, Interval velocities from seismic reflection time mea-surements: Society of Exploration Geophysics Publications. [19] Jackson, D.D., 1972, Interpretation of inaccurate, insufficient and inconsistent data: Geophysical Journal of the Royal Astronomical Society, 28, 97-109. [20] Johnson, S.H., 1976, Interpretation of split-spread refraction data in terms of plane dipping layers: Geophysics, 41, 418-424. [21] Lines, L.R., and Treitel, S., 1984, Tutorial: A review of least-squares inversion and its application to geophysical problems: Geop. Prosp. 32, 159-186. Bibliography 111 [22] Mandel, L., 1974, Interpretation of Instantaneous frequencies: American J. Physics, 42, 840-846. [23] May, B.T. and Covey, J.D., 1981, An inverse ray method for computing geologic structures from seismic reflections - zero offset case: Geophysics, 46, 268-287. [24] McMechan, G.A., 1983, P-x imaging by locabzed slant stacks of T-x data: Geophys-ical Journal of the RAS, 72, 213-221. [25] Menke, W., 1984, Geophysical data analysis: discrete inverse theory; Academic Press Inc. [26] Milkereit, B., 1987, Decomposition and inversion of seismic data - an instantaneous slowness approach: Geop. Prosp., 35, 875-894. [27] NeideU, N.S., and Taner, M.T., 1971, Semblance and other coherency measures for multichannel data: Geophysics, 36, 482-497. [28] Nikobc, Z.J., 1975, A recursive time-varying bandpass filter: Geophysics, 40, 520-526. [29] Oldenburg, Doug W., 1984, An introduction to bnear inverse theory: IEEE Trans-actions on Geoscience and Remote Sensing, GE-22, 665-674. [30] Oldenburg, D.W., Levy, S., and Stinson, K., 1984, Root-mean-square velocities and recovery of the acoustic impedance: Geophysics, 49, 1653-1663. [31] Oppenheim, A.V., and Schafer, R.W., 1975, Digital Signal Processing: Englewood Cliffs, N.J., Prentice Hall. [32] Pann, K., and Shin, Y., 1976, A class of convolution time-varying filters: Geophysics, 41, 28-43. Bibliography 112 [33] Peacock, K.L., 1985, Kaiser-Bessel weighting of the Hilbert transform high-cut filter: IEEE Trans. Acoust., Speech, Signal Processing, ASSP-33, 329-331. [34] Ready, R.J., and Wintz, P.A., 1973, Information extraction, SNR improvement, and data compression in multispectral imagery: IEEE Transactions on Communications, COM. 21(10), 1123-1130. [35] Robertson, J.D., and Nogami, H.H., 1984, Complex seismic trace analysis of thin beds: Geophysics, 49, 344-352. [36] Scheuer, T.E., and Oldenburg, D.W., 1988, Aspects of time-variant filtering: Geo-physics, 53, 1399-1409. [37] Scheuer, T.E., and Oldenburg, D.W., 1988, Local phase velocity from complex seis-mic data: Geophysics, 53, 1503-1511. [38] Shanks, J.L., 1967, Recursion filters for digital processing: Geophysics, 32, 33-51. [39] Stein, R.A., and Bartley, N.R., 1983, Continuously time-variable recursive digital bandpass filters for seismic signal processing: Geophysics, 48, 702-712. [40] Taner, M.T., and Koehler, F., 1969, Velocity spectra - digital computer derivation and appbcation of velocity functions: Geophysics, 34, 859-881. [41] Taner, M.T., Koehler, F., and Sheriff, R.E., 1979, Complex trace analysis: Geo-physics, 44, 1041-106,, [42] Tarantola, A., 1987, Inverse Problem Theory: Elsevier Science Pub. [43] Treitel, S., Shanks, J.L., and Frasier, C.W., 1967, Some aspects of fan filtering: Geophysics, 32, 789-800. Bibliography 113 [44] Wiggins, J.W., 1988, Attenuation of complex water-bottom multiples by wave-equation-based prediction and subtraction: Geophysics, 53, 1527-1539. [45] Wood, L.C., 1974, Seismic data compression methods: Geophysics, 39, 499-525. [46] Yilmaz, O., 1987, Seismic data processing: SEG publications, Tulsa, Oklahoma. [47] Ziolkowski, A., and Fokkema, J.T., 1986, The progressive attenuation of high-frequency energy in seismic reflection data: Geop. Prosp., 34, 981-1001. Appendix A Derivation of the complex bandlimited sine function The idealized complex bandlimited sine function is represented in the frequency domain as a boxcar of magnitude 2 between the frequencies u>\ and u>2. Explicit inverse Fourier transformation between two positive frequencies where o>2 > ^ 1 , gives 1 rW2 ^ ^ 27T iw 2elwtdu>, -—[cosu2t — cosw-J, + i(sinuj2t — sinuit)], iirt 2sinujt id where u> = (UJ2 +o>1)/2, and u> = (o>2 — u\~)j2. A complex sine function defined between the negative frequencies — Ui and —u)2 is given by s(t) = ——e-"*. (A.l) It then follows that a general expression for a complex sine function with any given band range of w\ to u2 is . 2sin\u>\t •-, . „. A low-pass complex sine function (u;1 = 0) with positive frequency energy is given as s(t) = 2sm^t/2eZUj2t/2 = —{sinuj2t-\-i(l — cosuj2t)}. The 2-D complex sine function, denned by a rectangular region in the Fourier k — u domain, can be derived utilizing equation (A.2). Multiplication of a spatial complex sine 114 Appendix A. Derivation of the complex bandhnaited sine function 115 function by a temporal complex sine function is equivalent to convolving the respective Fourier domain boxcars and, thus, produces the desired rectangular pass region. The analytic expression for such a function containing energy in the frequency band between UJI and u>2 and the wavenumber band between k\ and A;2 is i{x,t) = a(x,t)ei{Qt-'kx\ (A.3) where, a[x,t) — (4/7T xt)sin\ \x • sm\ —-\t, i(0,0) = \u>2 - Wi | • \k2 - kx\l-K2. The sign reversal on the spatial wavenumber results from the definition of the Fourier transform over spatial coordinates (see equation (2.12)). This definition produces a positive phase velocity given positive frequency and wavenumber components. Appendix B Computation of the digital complex delta function The complex delta function may be defined in terms of a generabsed function as the asymptotic bmit of a low-pass complex sine function (see Appendix A), i.e., 6(t) = bm —[sinuit + i(l — cosutfj. (B-l) w •oo The low-pass complex sine function reproduces the spectral characteristics of the com-plex delta function within its pass-band. Therefore, an unabased digital complex delta function with sample interval At is computed by limiting the upper frequency in (B.l) to the Nyquist value of U>N — TV/ At. The resulting expression is 8(mAt) = — (^sinmir + z(l — cosmivj^, m = 0, ±1, ± 2 , . . . , (B-2) which reduces to S(mAt) = < it> m = 0, i2 ... m = ±1 , ±3, ±5 , . .. , 0, otherwise. The imaginary part of this result is an unabased Hilbert transform operator. Although the digital complex delta function is infinite in length, it decays rapidly to zero away from t = 0. It therefore can be truncated for practical implementation. Figure B . l shows a normabzed digital complex delta function with a one-sided length of 200 ms for At = 4 ms. In addition to simple truncation, the complex delta function may also be weighted by a window function w(t). The type of truncation window used, and the truncation length, determine important spectral characteristics of the digital 116 Appendix B. Computation of the digital complex delta function 117 Re Im i i i i i i i i i -0.20 -0.10 0.0 0.10 0.20 time (sec) Figure B . l : The digital complex delta function sampled at 4 ms and truncated at 200 ms. complex delta function. This in turn, helps determine the spectral response of the time-variant filter. Consider the high-pass filter in equation (1.14) when the window function is included r{t) = eiUct[S(t) • w(t)/2 * i( i)e _ i " c t ] . (B.3) Taking the Fourier transform produces R(w) = 8(UJ - ue) * (\U(LU) * W(w)] • S(u + w c)) , (B.4) where W(a>) is the Fourier transform of w(i). The term u(u)) * W(o>) is the smoothed spectrum of the windowed complex delta function. It is a modified step function that exhibits deviations from unity in the positive frequency band and most likely leaks energy into the negative frequency band. Since S(u> + UJC) in equation (B.4) is multiplied by U{OJ) * W(a>), this product exhibits rolloff characteristics which are preserved in the final output of the filter. For a time-invariant filter then, the spectral characteristics of the filter are identical to the spectral characteristics of the windowed digital complex delta Appendix B. Computation of the digital complex delta function 118 function. In the time-variant case however, there is an additional complication that arises due to gradients in the instantaneous cutoff frequency function. That complication is considered in section 1.5. Oppenheim and Schafer (1975, chap. 5) discuss several types of truncation windows and their spectral shaping characteristics. Any one of them could be used to satisfactorily shape the spectrum of the digital complex delta function. I use a simple digital cosine window of the form, •77-777 A / . w(mAi) = c o ^ l ^ ) , (B.5) where T is the one-sided truncation length and a is an exponent that allows a family of windows. . The effects on the spectrum of the digital complex delta function by varying the exponent a and the truncation length T are shown in Figure B.2. Figure B.2a shows ampbtude spectra of 6(mAt) • w(mA£) when a=0,2,and 8, fixing T=200 ms and Figure B.2b shows the effect of varying T from 100 to 200 to 800 ms, fixing a = 2. It is evident that choosing an optimum cosine window involves a trade-off of leakage into the negative frequency spectrum against smooth spectral response and computational efficiency. Increasing the value of a gives a smoother spectral response and decreasing the length of T leads to faster computation speeds, but both lead to more rolloff into the negative frequency band. Nevertheless, the spectral characteristics of the windowed complex delta function play a role in describing the spectral characteristics of the output from the time-variant filter. Appendix B. Computation of the digital complex delta function 1 1 9 — I . 1 I I L_ - 4 0 - 2 0 0 2 0 4 0 f r e q u e n c y (Hz) Figure B.2: (a) The spectral response of the digital complex delta function when the exponent of the cosine window is 0, 2, and 8, for a truncation length of 200 ms. (b) The spectral response of the digital complex delta function with truncation lengths of 100, 200, and 800 ms with the window exponent fixed at a = 2. The spectrum is shown around 0 Hz to emphasize the rolloff characteristics; these characteristics also occur at the Nyquist frequency. Appendix C Resolution of Local Phase Velocity Seismic data are finite in temporal and spatial extent, discretely sampled, and bandlim-ited in frequency content. These conditions must impose limitations on the resolution of phase velocity using equation (2.20). The following analysis attempts to provide maxi-mum error bounds on the computation of phase velocity in terms of the above limitations. If the total recording time is denoted T, and the spatial extents of the seismic record are denoted X and Y, then by the uncertainty principle (Claerbout, 1976), the best resolution obtainable in the wavenumber-frequency domain is given by the increments, dw = 2TT/T, dkx = 2TT/X, dky = 2TT/Y. (C.l) This principle can be extended to provide a maximum uncertainty in the computation of phase velocity using estimates of local frequency and wavenumber as given by equation (2.20). We begin by taking the magnitude of equation (2.20) ' v = = v*7F (C'2) where zv is the instantaneous frequency and A and rj are the instantaneous wavenumber estimates at a particular location in the record. The differential of this result gives K K,3 where K = \ / A 2 + if 120 Appendix C. Resolution of Local Phase Velocity 121 To obtain a maximum bound on this differential we take the absolute value of each term giving \dV\'< K^idw+ —jXdX+71dr]}), (CA) where all variables are taken positive. Using the uncertainty relations (C.l) to define maximum uncertainty in measurements of the instantaneous parameters we have dzu = 2n/T, dX •= 2TT/X, di]=2Tr/Y. (C.5) Substituting these into (C.4) produces the following upper bound on the resolution of phase velocity, 2TT IK VIA! V\i\\\ ,„ x W S « ( T + -T + -T)- ( C 6 ) In two dimensions this expression reduces to This result shows expbcitly how the resolution of phase velocity depends on the temporal and spatial extents of the data and on the sought phase velocity. Not surprisingly, higher phase velocities have higher maximum bounds. Imphcit in this result is the dependence on the spatial sample interval Aa; since it bmits the maximum possible A. And for phase velocities greater than Ax/ At, X is also limited by the upper frequency bmit of the data (where At is the temporal sample interval). In practice, we have found that estimation of local phase velocity from noise free synthetics using equation (2.22) to be well within the bounds imposed by equations (C.6). Appendix D Approximate traveltime relations for dipping layers Experience from ray tracing dipping models such as the one shown in Figure 3.10, suggests that simple, yet approximate, relationships exist between ray angles di, /? ; , 77;, & of the dipping model, and the ray angle from an equivalent horizontal model. The approximate relations are given by d ; ~ 6i + /2; • Vi - &i ~ Pi ii -Oi+ii ' (D.l) where fi and ept are perturbation angles that complete the relations. I show that these relations are more accurate when layer dips are relatively small. As shown in Figure D. l , d{ refers to a refraction angle referenced to the interface normal, while a; refers to the refraction angle referenced to a vertical axis. These angles are simply related by the interface dip 8i as given below; ai — di — 8i rn = rji - f 8i ft • Pi ~ 8, .6 ^ 1 8{. (D.2) 122 Appendix D. Approximate traveltime relations for dipping layers 123 Figure D. l : Geometrical relationships of forward and reverse ray paths used for the development of an approximate average traveltime relation. Appendix D. Approximate traveltime relations. for dipping layers 124 The validity of relations (D.l), for small angles p and e, can be established by in-duction. I proceed by inverse ray tracing downward along the upcoming branches of the forward and reverse ray paths. First, by assuming that k-i = ^ - 1 - ^ - 1 (D.3) &_! = i i l t i , (D.4) hold identically in the i — 1th layer, I show that for rays impinging upon the zth interface, Snell's law implies that A (D.5) Ci ^Oi + ii- (D.6) Choosing forward and reverse ray paths with a common emergence angle satisfies (D.3) and (D.4) in the first layer and hence forms the basis for the following inductive reasoning. Applying Snell's law at the zth interface for the forward ray gives . sinfli'= (vi/vi-i)sinf3i-i (D.7) Substituting (D.3) into this result gives sinfii .= (vi/vi^l)sin(9i^1 - e;_i) = sinOiCosii^i — (vi/vi-i)cos8i-isinii-\, (D.8) In addition, Snell's law provides sin(6i ~. ii) = (yilvi-l)sin(6l-1 ~ Q-J). (D.9) Expanding the sine functions gives sinOiCOsii — cosOisinii = sin9iCosii^{ — (vi/vi_i)cos9i-.isinei_i. (D.10) Appendix D. Approximate traveltime relations for dipping layers 125 When e, and are small (or when ii ~ ei-i), I obtain the following approximation cosdisinii ~ (vi/vi_x)cos9i_isinii^i. (DAI) Substituting this into equation (D.8) and assuming e,_i is small provides sinfli ~ sinOi — cos8isinii. (DA2) Taking the inverse sine leads to fli ~ sin~1[sin9i — cosBisinii]. (D.13) Then, expanding the right hand side of equation (D.13) in a Taylor series produces sin2e"i _ Pi~9i- smii + ^ — ^ + • • • . (D.14) When ii is small, sinii ~ ii, and the higher order terms in the Taylor expansion can be dropped to give the desired result Pi~9i-ii. (D.15) The derivation for & proceeds similarly. Applying SneU's law at the zth interface for the reverse ray gives sin(i = (vi/vi-i)sin£i-i. (D.16) Substituting (D.4) into this equation and expanding gives sin£i = sin9iCosii-i + (yi/vi_i)cos9i-\sinii_i. (DAI) Utibzing relation (D.ll) and assuming i{_i is small leads to sin£i ~ sinBi + cos9isinii. (D.18) Taking the inverse sine of this result gives ~ sin~1[sin9i + cos9isinBi\. (D.19) Appen'dix D. Approximate traveltime relations for dipping layers 126 Expanding the right hand side of (D.19) in a Taylor series produces ? --* SZTh 6' £i~6i + sinii + - ± + •••. (D.20) Assuming small ii and dropping the higher order terms in sinii produces the desired approximation ii ^ Bi + ii. (D.21) A similar approximate relationship holds for a; and iji] but an independent set of perturbation angles \ii is needed. In the layer of reflection, I have shown that 8 ~ 9 - i Hn — "n cn> i n ~&n + in. (D.22) Upon reflection Oin = 0n -On - £n, 7 j n = £„ - f l n + t V (D.23) Therefore, the inductive arguement can be extended back to the surface to yield approx-imate relations a; ~ Oi-fn, (D.24) Vi ^ Oi+Pi, (D.25) A new set of perturbation angles are needed since the ray tracing procedes in an opposite direction back to the surface. Equations (D.2) show how to convert the relations (D.15), (D.21), (D.24), and (D.25) into angles referenced to the vertical axis. The results are Appendix D. Approximate traveltime relations for dipping layers 127 77; ~ 6{ - pi, (3{ ~ 9i - eiy . 6 ^Oi + ei, ' (D.26) where, pi = jii + 6i, and e; = e,- — These results allow equation (3.22) in the text to be approximated by tave ~ ~ + lf] —{cos(6i - p^ + cos(6i - ei) + cos(6i + p{) + cos{6i + e,-)}.' (D.27) VP 2 i=l vi Consider the term r { = cos(6i - pi) + cos(6i - e{) + cos(6i + m) + cos(6i + ei). (D.28) This can be rewritten as Ti = 2cos9i(cospi + cose,). (D.29) Substituting this simpbfication into equation (D.27) gives tave = — + 2 2 — C O S ^ ^ 1 • (D.30) Thus, when ^ and e,- are small equation (D.27) can be approximately written as t a v e ^ ^ + 2 J 2 — ^ei, (D.31) which resembles the traveltime relation for a horizontally layered system. Expbcit in the above approximations is that e;,e; and Pi,Pi must be small. However, the size of these angles are model dependent. This makes a definitive error analysis dif-ficult. Factors such as layer dip, emergence angle, and the configuration of subsurface boundaries, all affect the degree to which equation (D.31). approximates traveltimes in a corresponding horizontal model. Nevertheless, examples shown in Section 3.5 demon-strate that the above relations provide a good foundation for assessing the results of average Snell trace inversion in the presence of planar dips.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Complex seismic data : aspects of processing, analysis,...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Complex seismic data : aspects of processing, analysis, and inversion Scheuer, Timothy Ellis 1989
pdf
Page Metadata
Item Metadata
Title | Complex seismic data : aspects of processing, analysis, and inversion |
Creator |
Scheuer, Timothy Ellis |
Publisher | University of British Columbia |
Date Issued | 1989 |
Description | A common theme in this thesis is the use of complex signal attributes to facilitate the processing, analysis, and inversion of seismic data. Complex data are formed from real data by removing the negative frequency portion of the Fourier transform and doubling the positive frequency portion. Removing the dual nature of frequency components in real data allows the development of efficient algorithms for time-variant filtering, local phase velocity estimation, and subsequent velocity-depth inversion of shot gathers using Snell traces. For 1-D seismic data, I develop a computationally efficient time-variant filter with an idealized boxcar-like frequency spectrum. The filter can either be zero phase, or it can effect a time-variant phase rotation of the input data. The instantaneous low-frequency and high-frequency cutoffs are independent continuous time functions that are specified by the user. Frequency-domain rolloff characteristics of the filter can be prescribed, but the achieved spectrum depends on the length of the applied filter and the instantaneous frequency cutoffs used. The primary derivation of this theory depends upon the properties of the complex signal and the complex delta function. This formulation is particularly insightful because of the geometrical interpretation it offers in the frequency domain. Basically, a high-pass filter, can be implemented by shifting the Fourier transform of the complex signal towards the negative frequency band, annihilating that portion of the signal that lies to the left of the origin, and then shifting the truncated spectrum back to the right. This geometrical insight permits inference of the mathematical form of a general time-variant band-pass filter. In addition, I show that the time-variant filter reduces to a Hilbert transform filter when the derivation is constrained to include real signal input. Application of the procedure to a spectral function permits frequency-variant windowing of an input time signal. For 2-D arid 3-D seismic data, I propose a new method that uses the concepts of complex trace analysis for the automatic estimation of local phase velocity. A complex seismic record is obtained from a real seismic record by extending complex trace analysis into higher dimensions. Phase velocities are estimated from the complex data by finding trajectories of constant phase. In 2-D, phase velocity calculation reduces to a ratio of instantaneous frequency and wavenumber, and thus provides a measure of the dominant plane-wave component at each point in the seismic record. The algorithm is simple to implement and computational requirements are small; this is partly due to a new method for computing instantaneous frequency and wavenumber which greatly simplifies these calculations for 2-D and 3-D complex records. In addition, this approach has the advantage that no a priori velocity input is needed; however, optimum stability is achieved when a limited range of dipping events is considered. Preconditioning the record with an appropriate velocity filter helps reduce the detrimental effects of crossing events, spatial aliasing, and random noise contamination. Accurate recovery of local phase velocity information about underlying seismic events allows the rapid evaluation of seismic attributes such as rms velocity and maximum depth of ray penetration. I utilize local phase velocity data from a shot gather for the estimation and inversion of Snell traces. The primary Snell trace corresponding to a 1-D velocity model locates all primary reflection energy,corresponding to a fixed emergence angle. Constraints on interval velocity and thickness obtained from several estimated Snell trajectories are inverted using SVD to provide a least squares velocity-depth model. The estimation and inversion is efficiently carried out on an interactive workstation utilizing constraints from a hyperbolic velocity analysis. Finally, Snell trace inversion is extended to an inhomogeneous medium. When dips are small, averaging Snell traces of a common phase velocity from forward and reversed shot gathers approximately removes the effects of planar dip. This allows recovery of velocity and depth vertically beneath the midpoint of the source locations used to obtain the reversed information. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-10-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052967 |
URI | http://hdl.handle.net/2429/29283 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1989_A1 S33.pdf [ 12.49MB ]
- Metadata
- JSON: 831-1.0052967.json
- JSON-LD: 831-1.0052967-ld.json
- RDF/XML (Pretty): 831-1.0052967-rdf.xml
- RDF/JSON: 831-1.0052967-rdf.json
- Turtle: 831-1.0052967-turtle.txt
- N-Triples: 831-1.0052967-rdf-ntriples.txt
- Original Record: 831-1.0052967-source.json
- Full Text
- 831-1.0052967-fulltext.txt
- Citation
- 831-1.0052967.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052967/manifest