APPLICATION OF WAVELET TRANSFORMS PROCESSING A N D T O SEISMIC INVERSION By Xin-Gong Li B. Sc. (Geophysics) Tong-Ji University, PRC M. Sc. (Geophysics) University of British Columbia A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T O F T H E REQUIREMENTS FOR T H E D E G R E E O F D O C T O R OF PHILOSOPHY in T H E F A C U L T Y O F G R A D U A T E STUDIES EARTH A N D O C E A N SCIENCES We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA 1997 © Xin-Gong Li, 1997 DATA 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. Earth and Ocean sciences The University of British Columbia 129-2219 Main Mall Vancouver, Canada V6T 1Z4 Date: Abstract The goal of this thesis is to use wavelet transforms as tools to analyze simultaneously the time (or location) and frequency (or scale) variant characteristics of seismic data and then to apply these information to two important geophysical problems: noise filtering and traveltime inversion. In the noise filtering problem, the discrete non-orthogonal wavelet transform is used to analyze the time-variant characteristics of signal and noise. This information can be used to filter noise in the transformed domain. This approach is applied to ground roll attenuation for field seismic data according to the localized characteristics. Geophysical inverse problems are, in general, non-linear, underdetermined and illposed. An initial model and some prior information describing the model is needed to find the model perturbation that minimizes an objective function. Traditionally, the perturbation is assumed to be an independent Gaussian or a correlated but stationary process, an assumption which is not physically correct according to the analyses of well logs. Indeed, analysis of logs described in this thesis suggests that log data can be decomposed into (a) a fractal process, usually non-stationary, with wavelet coefficients which follow a power law, and (b) some non-fractal structures including a large scale trend and some spiky details. The inverse problem can be solved using the wavelet transform in two steps. First, I solve an overdetermined problem to invert the large scale trend, then, I solve for a model with fractal constraints. Another important application of the wavelet transform in inverse problems is to reduce the effects of the noise in the input data. Since the noise effect varies with scale, the regularization can be done accordingly. Wavelet transform constrained inversions are tested using 1-D and 2-D synthetic examples. n To the memory of my grandmother m Table of Contents List of Tables viii List of Figures ix Acknowledgment xii 1 INTRODUCTION 1 1.1 A brief history of wavelet transform (WT) 1 1.2 The scope of this study 3 1.3 Some wavelet transform resources 7 1.4 Overview of the thesis 8 2 T H E T H E O R Y OF WAVELET T R A N S F O R M (WT) 10 2.1 Introduction 10 2.2 Non-orthogonal WT - Frame Theory 10 2.2.1 Definitions of frames 11 2.2.2 Gabor Frames 14 2.2.3 Wavelet Frame 20 2.3 Orthogonal wavelet transforms 21 2.4 Summary 26 3 T H E G A B O R T R A N S F O R M IN SEISMIC FILTERING 36 3.1 Introduction 36 3.2 Filtering using the Gabor transform (GT) 37 iv 3.2.1 Algorithm 37 3.2.2 Filtering random noise for signal with localized characteristics 3.2.3 Attenuation of the coherent noise 39 3.2.4 Gabor filtering of real data 40 . . 38 3.3 Relationship between the GT and FT filtering 41 3.4 Gabor filtering for 2-D data 43 3.4.1 Reasons for a 2-D GT 43 3.4.2 Formulation and algorithm 44 3.4.3 Example 45 3.5 Summary 47 4 A T T R I B U T E ANALYSIS USING T H E CONTINUOUS W T 64 4.1 Introduction 64 4.2 The Morlet wavelet transform 65 4.2.1 The wavelet transform in the form of convolution 65 4.2.2 Fourier transforms of the Morlet coefficients and multiscale traces 67 4.3 Understanding wavelet attributes 69 4.3.1 Wavelet attributes of spikes 69 4.3.2 Wavelet attributes of traces 71 4.4 Wavelet attribute analysis for 2-D data 72 4.5 Summary 73 5 ANALYSIS A N D SYNTHESIS OF W E L L LOGS 5.1 5.2 5.3 Introduction 88 . 1/f-type processes 88 89 Diagonalization of the covariance matrix 91 5.3.1 91 Covariance matrix in a transformed space v 5.3.2 5.4 5.5 5.6 5.7 Diagonalization of a F B M covariance matrix using the WT . . . 93 Fractal parameter estimations of logs 94 5.4.1 The wavelet transform approach 94 5.4.2 Examples 95 Synthesis of fractal processes 98 5.5.1 1-D 98 5.5.2 2-D 99 Non-fractal component 101 5.6.1 Definition of the non-fractal component 101 5.6.2 A simple way to evaluate the fractal model 101 Summary 103 6 INVERSION WITH PRIOR INFORMATION 114 6.1 Introduction 114 6.2 Traveltime tomography 115 6.3 The maximum a posteriori probability (MAP) solution 117 6.4 The relationship between regularization and stochastic approaches . . . . 119 6.4.1 On the covariance matrix, damping and smoothing 119 6.4.2 Regularization of a 1/f model using damping and smoothing . . . 121 6.4.3 Conjugate gradient approach with damping and smoothing . . . . 122 6.4.4 The discrete derivative matrices 123 6.5 Summary . 125 7 W A V E L E T T R A N S F O R M INVERSION (WTI) 127 7.1 Introduction 127 7.2 Wavelet transform inversion (WTI) 128 7.2.1 128 Formulation VI 7.2.2 7.3 7.4 On the relationship between WT and Fourier domain constraints 130 Two types of multiscale constrained models 130 7.3.1 The scale-truncated model 130 7.3.2 The fractal model 133 7.3.3 WTI and the damped least squares 134 Example 1: inversion of a blocky model 134 7.4.1 Inversions via damping and smoothing 135 7.4.2 Inversion via constraints in a transformed domain 135 7.5 Example 2: inversion of a F B M model 137 7.6 Example 3: inversion of large problems 138 7.6.1 139 7.7 Model space constrained inversions 7.6.2 WTI solutions 139 7.6.3 Model correlation matrix 140 Summary 141 8 SUMMARY 153 Appendices 165 A The Morlet W T of 6(t; t ) 165 0 vu List of Tables 2.1 Comparison of the three types of Gabor frames 19 5.1 The random processes with different fractal parameters 91 5.2 The fractal parameters estimated from Logl 96 5.3 Estimated fractal parameters from Well 2, 3 and 4 97 5.4 The number of non-fractal outliers with different thresholds 6.1 Estimated damping and smoothing coefficients for a 1/f model when a = 1. 122 7.1 Derivatives used to constrain the 2-D inversion vin 102 140 List of Figures 2.1 The auxiliary function for a dual Gabor frame 28 2.2 The window functions in forward and inverse dual frames 29 2.3 The sampling of self-adjoint frames 30 2.4 Correlated windows in the Gabor transform 30 2.5 The approximate auxiliary function using least squares 31 2.6 Comparison between the Morlet and the Gabor wavelets 32 2.7 The direct implementation and the a trous algorithm for CWT 2.8 The pyramid algorithm 33 2.9 Two Daubechies wavelets 34 33 2.10 2-D orthogonal wavelet transform 35 3.1 Matrix structures of the Gabor transform and Gabor coefficients .... 3.2 Synl: synthetic trace of two Pucker wavelets and random noise 3.3 Gabor transform of synl 3.4 Fourier filtering of synl 51 3.5 Gabor filtering of synl 52 3.6 Syn2: two Packer wavelets and a chirp signal 3.7 Estimated Ricker wavelets using GT and FT approaches . 3.8 Tr70: a real trace and its Gabor and Fourier spectra 3.9 Trace trlO after Gabor and Fourier filtering 49 . . 3.10 Shot2: a real shot record 48 50 53 54 55 56 57 3.11 F-k transform of Shot2 57 ix 3.12 2-D Gabor and f-k bases 58 3.13 Shot2 and its Gabor subbands 58 3.14 The estimated low frequency signal and noise 3.15 The estimated signal after the attenuation of low frequency noise 60 . . . . 61 3.16 The filtering of refracted waves 62 3.17 F-k spectra of estimated reflected waves with and without refraction . . . 63 4.1 Fourier decomposition of a spike 75 4.2 Morlet WT decomposition of a spike 76 4.3 Unwrapped FT and WT phases for a spike in the middle 76 4.4 FT decomposition of two spikes 77 4.5 Morlet WT decomposition of two spikes 78 4.6 Morlet WT attributes and multiscale traces of a Ricker wavelet 79 4.7 Morlet WT attributes and multiscale traces of two Ricker wavelets . . . . 80 4.8 Morlet WT attributes and multiscale traces of Ricker wavelets with noise 81 4.9 Wavelet transform attributes and multiscale traces for a portion of a real trace 82 4.10 The instantaneous amplitudes and phases for a shot2 83 4.11 The amplitudes of the wavelet coefficients at a few scales 84 4.12 The phases of the wavelet coefficients at a few scales 86 5.1 Covariance matrix of a F B M and its wavelet transforms 105 5.2 Log 1 106 5.3 The wavelet coefficients and multiscale decompositions of segment 3 . . . 107 5.4 The multiscale variances for segments 1 and 2 of Log 1 108 5.5(a)The sonic and density logs of Well 2, 3 and 4 109 5.5(b)Three segments of the sonic and density logs of Well 2, 3 and 4 110 x 5.6 Synthesis of 1-D F B M processes Ill 5.7 Synthesis of 2-D F B M images 112 5.8 Decomposition of fraction and non-fractal structures from Log 1 113 6.1 Covariance matrix and weighting matrix for a fractal model 126 7.1 Model 1 and its derivatives, FT and WT 143 7.2 The damped least squares (DLS) solution of smallest model 143 7.3 The flattest and smoothest models 144 7.4 The truncated SVD model 144 7.5 Fourier transform constrained inversion (FTI) 145 7.6 Wavelet transform constrained inversion (WTI) 145 7.7 WTI with location constraints 146 7.8 Model 2: F B M 146 7.9 The trend and F B M components of Model 2 147 7.10 WTI and damped least squares solutions 148 7.11 Wavelet coefficients of the inverted models 148 7.12 Cross borehole geometry 149 7.13 Model 3 149 7.14 Inversion using damping constraints 150 7.15 Inversion using damping and x direction smoothing constraints 150 7.16 Inversion using damping and 2-D smoothing constraints 151 7.17 Inversion using 1-D WT constraints in horizontal direction 151 7.18 Inversion using 2-D WT constraints 151 7.19 Covariance functions of model space and wavelet space constrained inversions 152 XI Acknowledgment I would like to thank my supervisor Dr. Tad Ulrych for his guidance, financial support and for the time and patience he devoted in help the writing of this thesis. My special thanks goes to Mauricio Sacchi for his critical reading of this thesis, for the discussions we shared on various topics about inversion and wavelet transforms and for the laughs he brought to the office which made this study much more enjoyable. My thanks also goes to Drs. Mat Yedlin and Bruce Buffett for their considerable input during the preparation of this thesis. I would also like to thank Dr John Anderson in M E P T E C for the helpful input. I am grateful to Yao-Guo Li, Colin Farguharson, Ken Matson, Denise Long for their proofreading of this thesis and for their tips about English writing. I am also grateful to my friends and colleagues in the department who made life fun and interesting. The research that leads to this thesis was supported in part by an NSERC grant to Dr. Tad Ulrych and by the Consortium for the Development of Specialized Seismic Technologies at UBC. Xll Chapter 1 INTRODUCTION The characteristics of geophysical signals are, in general, variable in both time and frequency. For example, different types of seismic waves appear at different times and have different frequency bands. Examining the significance of the amplitudes of the waves can provide timing or location information, while analyzing the Fourier spectrum can provide frequency information. The time and frequency domain information usually was not analyzed simultaneously in the past because of the lack of mathematical tools to do so. The recent rapidly developed wavelet theory offers such a tool of time and frequency analysis. This thesis studies the application of the wavelet transforms to two common seismic problems: signal processing and inversion. This chapter briefly reviews the development of the wavelet transform and then describes the scope of this thesis. 1.1 A brief history of wavelet transform (WT) The wavelet transform (WT) can be considered as the decomposition and reconstruction of a signal using a localized basis. In terms of orthogonality the WT can be generally divided into orthogonal and non-orthogonal WTs. Orthogonal WT theory, as will be described later, was developed only in the last decade, while the study of non-orthogonal WT dates back to the applications of the short time Fourier transform (STFT) which can be treated as a special WT. The first application of a non-orthogonal WT was the application of the Gabor transform in signal analysis (Gabor, 1946). The Gabor transform is a short time Fourier transform (STFT) with a Gaussian window. Dziewonski 1 Chapter 1. INTRODUCTION 2 et al. (1969) developed a 'multiple filter technique' using a modified Gabor transform with dilated windows to analyze the time-variant features of surface waves. Later, Morlet et al. (1982) used a similar concept to study wave propagation and sample theory in geophysics. The wavelet basis used, named after Morlet, is a fixed complex waveform dilated to different scales and shifted to different locations. The Morlet basis was also used by Goupillaud et al. (1984) to represent seismic traces. The concept initiated by Gabor was quickly combined with theories in engineering, physics, pure mathematics and other fields, and many new theoretical developments have been presented. Except for applications of the non-orthogonal wavelet transforms in signal analysis, more general applications are limited by the fact that information is lost in the process of forward and inverse transformation. This problem is solved with the development of orthogonal WTs. The study of orthogonal WTs by Daubechies (1988) and the development of the pyramid algorithm (Mallat 1989a, b) have increased the efficiency of applications of the WT in image processing, data compression and noise attenuation in the late 80's and early 90's. In recent years, many new orthogonal wavelet transforms have been adopted by different researchers based on different mathematical requirements. The commonly used wavelets, such as the Haar, Daubechies, and Meyer wavelets, are named after their inventors. An excellent review concerning the detailed history and theory of wavelet transforms can be found in the books by Daubechies (1992) and Meyer (1993). A brief description of a fast algorithm is given in the next chapter. In the last decade, the wavelet transform and related localized transforms, have been applied in various branches of science. Much of the research has occurred in the fields of mathematics, and signal and image processing. The rapid progress has benefited greatly from Internet communications between the different fields. Among different applications of orthogonal wavelet transforms, the most important one is probably data compression and one of the successful examples is the compression of finger print images adapted by Chapter 1. INTRODUCTION 3 the Federal Bureau of Investigation of the United States to save computer storage space (See, for example, Wickerhauser, 1991; Strang and Ngyuen, 1996). The orthogonal wavelet transform has shown great potential in compressing seismic data where such data are currently stored on hundreds of magnetic tapes. This application has been studied using wavelet packet decomposition (Luo and Schuster, 1992) and the wavelet transform (Bosman and Reiter, 1993; Donoho et al., 1995). Seismic data compression can also be used to save satellite transmission time in the so-called "almost real-time seismic survey" where the seismic data collected in the field, specially ocean and offshore, is compressed and transmitted to a processing center (Stigant et al., 1995). Because of the increasing need for seismic data compression, the 1996 annual meeting of the Society of Exploration Geophysicists hosted a special workshop dealing with seismic data compression standard and related issues (see 1996 SEG extended abstracts). Another application in seismology which has attracted much attention is the compression of the migration operator to save both computing time and storage space. This application can be achieved by considering the migration as a multiscale process and by ignoring the computations at some of the unimportant scales and locations. Research in this area has been described by Dessing and Wapenaar (1994, 1995) and Wu and McMechan (1995). Wavelet transforms have also shown great potential in other geophysical problems such as filtering and inversion. These aspects will be discussed in the next section and detailed in the thesis. 1.2 The scope of this study- Geophysical data processing and inversion are sometimes implemented by means of the Fourier transform (Stolt, 1978), and other transforms such as the r — p (Noponen and Chapter 1. INTRODUCTION 4 Keeney, 1986) and K-L transforms (Jones and Levy, 1987). The purpose of a particular transform used is dependent on the nature of the problem at hand. For a filtering application, signal and noise usually overlap in the spatial domain and the quest is to transform the data to a new domain where the signal and noise may be optimally separated. In the solution of ill-posed inverse problems, the incorporation of prior model information available in the transformed model and/or data spaces is generally required. The problem with the transforms mentioned above is that the bases are of infinite extent in the continuous case, or periodic in the discrete case. Such transforms are therefore only suitable for data and models with no localized characteristics, which implies stationarity. Unfortunately, since seismic data are, in general, not stationary, many important characteristics related to the localized nature of the data and model cannot be incorporated using these transforms. It is the goal of this thesis to use the localized nature of information to improve the filtering and inversion of seismic data. In particular, the following topics will be studied. Filtering of localized coherent noise Attenuation of various types of noise is an important and ambiguous procedure in geophysical data processing. Usually, as mentioned above, the data are transformed to another space where differences in the signatures of the signal and noise in the transformed space are utilized to attenuate the noise. Since signal, noise, or both are often localized in terms of frequency (scale) and spatial location, the processing results can be improved by the inclusion of the relevant information. The wavelet transform provides a tool for incorporating this localized information in both frequency and location. The problem in WT noise attenuation is to find the right type of wavelet transform to attenuate noise according to the time and location difference between the signal and noise. To attenuate random noise, for example, a denoising technique (Donoho, 1992) can Chapter 1. INTRODUCTION 5 be adapted with the advantage of using fast orthogonal WT algorithms. The subsequent filtered data have been shown to have better resolution for velocity analysis due to the denoising (Olmo et al., 1994). The wavelet transform has been used by a number of investigators for filtering coherent noise in the form of ground roll (Schuster and Sun, 1993; Foster et al., 1994; Miao and Moon, 1994; Miao and Cheadle, 1995). Unfortunately these studies have not fully utilized the location differences between the ground roll and signal. Such differences have been subsequently used in application of the 2-D Gabor transformed domain to attenuate the effects of ground roll (Li and Ulrych, 1996a). Wavelet transform attribute analysis Taner et al. (1979) introduced the use of complex trace attributes to indicate geological changes. A complex trace is computed using the original trace as the real part and the Hilbert transform of the trace as the imaginary part (Claerbout, 1992). The computed attributes are the contribution of all the frequency components at each time (or location). In some situations, for example, in the case that many modes of seismic waves exist in a record, we may wish to analyze the attributes of a certain mode in both time and frequency (or scale). Such time- and scale-variant characteristics can be analyzed using the multiscale attributes defined by complex wavelet coefficients. Application of similar attributes has been used for many years in analyzing the modes of surface waves (Dziewonski et al., 1969). In this thesis, these attributes are used to analyze the distribution of reflections and ground roll (Li and Ulrych, 1996b) in both time and scale. Some preliminary results show that wavelet attributes are sensitive indicators for certain time-frequency characteristics and may be useful in reservoir characterization. Chapter 1. INTRODUCTION 6 Analysis and synthesis of well log data It has long been known that many processes in nature can be modeled as non-stationary fractals (Mandelbrot and Van Ness, 1968). Unfortunately, a non-stationary process cannot be characterized as a 1/f-type process using the Fourier transform as can be done for a stationary fractal process (Flandrin, 1989). Recent studies have shown that the wavelet transform provides a suitable tool for the analysis of the non-stationary nature of fractional Brownian motion, F B M (Flandrin, 1989, 1992; Wornell, 1990; Wornell and Oppenheim, 1992). Applying such concepts in this thesis, well logs are considered as FBMs and their fractal parameters are estimated to represent regional geological characteristics. Well logs are also synthesized as F B M processes from given fractal parameters. Wavelet transform inversion (WTI) For any ill-posed inversion problem, prior information must be incorporated in some manner to constrain the inversion. The orthogonal wavelet transform provides an important tool for this purpose where the prior information may be incorporated in both time and location. This approach was first studied in medical tomography (Bhatia, 1994; Kolaczyk, 1994) and has been extended to seismic traveltime tomography where larger, and seriously ill-posed problems are encountered. In particular, this constrained inversion, or wavelet transform inversion (WTI) in short, has been used to incorporate scale and location information (Li et al., 1994, 1996) and in large 2-D tomographic problems with the help of a conjugate gradient approach (Li and Ulrych, 1995). In the above four applications, the choice of wavelet transform basis is based on the characteristics of the data and the model to be inverted. In contrast to the Fourier transform, which consists only of one basis, the localized transforms can have various bases according to the type of data and the problem being tackled. Although this fact Chapter 1. INTRODUCTION 7 makes things difficult because there is no unique theory to guide the choice of basis, the flexibility in the choice of basis makes the localized transform extremely flexible and powerful. By choosing the right basis, we can achieve a very good signal and noise separation in the transformed domain, or the maximum amount of prior information can be incorporated to constrain the inversion. In this thesis, the Gabor transform is used for filtering of coherent noise, the Morlet wavelet transform is used for localized attribute analysis, and the orthogonal wavelet transform is adopted for the last two applications mentioned above. 1.3 Some wavelet transform resources Since wavelet transform theory and application have been developed very rapidly during the last few years, this thesis does not intend to develop the theoretical aspects in any detail. This section is intended to provide some resources to the theory and some useful Internet addresses where the most recent research results can be obtained. Reviews of theory can be found in the books by Daubechies (1992), Meyer (1993), Chui (1992a, b) and Strang and Nguyen (1996). For applications, there have been two important international conferences on wavelet transforms in 1987 (Combes et al. (eds.), 1990) and 1989 (Meyer (ed.), 1989) in MarsieUe, France. The journal of "IEEE Transactions on Information Theory" has a special issue on the wavelet transform (March, 1992) which provides a wide range of applications. The wavelet transform literature has expanded greatly and it is important to keep up to date on the most recent results of other researchers. The Internet monthly news letter, Wavelet Digest, to be found at http://www.math.sc.edu/~wavelet/ and set up by the Industrial Mathematics Initiative, Department of Mathematics, University of South Carolina, is one of the most helpful Internet resources. Chapter 1. INTRODUCTION 8 There are a few leading groups working on different areas of WT applications. Visiting their home pages and reading their reports have been a great help to this study. The most notable are Mallat's group at New York University working on image processing (Mallat, 1989a, b; Mallat and Hwang, 1992; Mallat and Zhong, 1992); a group at Yale University working on image compression and signal detection using a library of wavelets (Bradley et al., 1993); Donoho and Jonestone and their students working on noise filtering, or denoising (e.g., Donoho and Jonestone, 1992). Software and research reports can be found at the ftp sites of the above three groups. All their ftp and web addresses can be found on the home page of the Wavelet Digest. Algorithms for implementing orthogonal wavelet transforms also can be obtained from various resources. A set of simple Fortran and C routines can be found in Numerical Recipes (Press et al., 1992). A collection of Matlab routines for wavelet and related transforms called WaveLab, which was used for the orthogonal wavelet transform computation in this thesis, was obtained without cost from the Internet site http://playfair.Stanford.edu /~wavelab/. 1.4 Overview of the thesis The thesis is arranged as follows. Chapter 2 provides a brief description of the basic theories of localized transforms. The non-orthogonal wavelet transforms such as the Gabor and Morlet transforms, are introduced in terms of frame theory. For orthogonal wavelet transforms, the fast pyramid algorithm using the quadrature mirror filters is described. In chapter 3, the Gabor transform is used for ground roll attenuation. Also, the relationship between Gabor filtering and conventional Fourier filtering (or f-k filtering in 2-D case) is studied. A field data example is tested and confirms that incorporating Chapter 1. INTRODUCTION 9 location information using Gabor filtering can greatly improve conventional f-k filtering. In order to analyze the time and frequency/scale characteristics of seismic data, multiscale attributes and multiscale traces are defined in Chapter 4 using the concept of the continuous wavelet transform. These multiscale attributes are used to monitor time and frequency features of ground roll and refracted and reflected waves. In Chapter 5, well logs are treated as a combination of a non-stochastic component and a component that can be statistically modeled as non-stationary fractional Brownian motion (FBM). The wavelet transform is used as a tool to analyze and synthesize the F B M component. During analysis, a fractal parameter is computed from a real log to determine the statistics of the appropriate F B M . During synthesis, this model parameter is used to generate a fractal process which is statistically similar to the real log. Chapters 6 and 7 contain perhaps the most significant contribution made in this thesis - the wavelet transform inversion using prior scale and location information. Chapter 6 reviews basic inverse theory as a means of introducing the WTI approach. Chapter 7 applies the WTI approach to incorporate two different types of scale constraints: blocky scale constraints and fractal scale constraints. The use of the conjugate gradient (CG) method for solving large 2-D WTI problems is studied and tested using synthetic examples. In the final chapter, some conclusions of this study are reiterated. Applications which are worth further study are also discussed. Chapter 2 T H E THEORY OF WAVELET T R A N S F O R M (WT) 2.1 Introduction The goal of this chapter is to study the theory of the WT and to prepare the background for the applications of WTs in geophysical data processing and inversion. The two different types of WTs, non-orthogonal and orthogonal WTs, are generally studied from different theoretical points of view. The non-orthogonal WTs such as the Gabor and Morlet transforms are studied in terms of frame theory (Daubechies, 1992). The definition and various forms of frames will be described in the next section. Orthogonal WT basis functions are constructed to satisfy certain mathematical requirements. Instead of pursuing the construction of orthogonal wavelet bases, this chapter will only describe the orthogonal WT in terms of a fast pyramid algorithm (Dautileux et al., 1987; Mallat, 1989a). 2.2 Non-orthogonal W T - Frame Theory The development of frame theory dates back to the study of the Gabor transform (Gabor, 1946), a short time Fourier transform (STFT) with a Gaussian window, in applications to signal communication. In geophysics, the Gabor transform has been used for analyzing the time-variant features of surface waves (Dziewonski et al., 1969) and this type of transform was later used by Morlet et al. (1984) for wave propagation and sampling theory. Since the STFT and the Morlet WT bases are non-orthogonal, the direct inverse 10 Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 11 transform using the conjugate operator, or the self-adjoint operator as in the Fourier transform, cannot recover the signal exactly. The following two approaches can solve this problem with the help of frame theory (Daubechies, 1992). First, the sampling density of the coefficients can be increased to achieve a good approximation using the self-adjoint operator where the quality of the approximation can be evaluated using frame bounds (Daubechies, 1992). Second, a basis which is orthogonal to that used in the forward transform can be constructed. In terms of frame theory, this basis is called a dual frame and can be used to reconstruct the data exactly. The basis function that constructs the dual frame for the Gabor transform is a sinusoid multiplied by an auxiliary function (Bastiaans 1981). Because of its non-physical shape, the auxiliary function has not been found useful in real cases except for image processing (Wang, et. al, 1994). This section starts with the definition of frames, and various forms of frames will be defined. Different methods of computing forward and inverse non-orthogonal wavelet transforms will be discussed using frame theory. 2.2.1 Definitions of frames In order to analyze the time-variant characteristics of a signal, a finite support basis is needed to decompose and reconstruct the signal. Unfortunately, most of the finite support bases such as the STFT and Morlet wavelet bases are not orthogonal. Frame theory (Daubechies, 1992) provides some fundamentals of evaluating a transform in terms of its accuracy and some other features. This subsection presents the definition of a frame and some examples in a convenient matrix form following Daubechies' book (1992). Let {tpj,j = 1, • • •, J} be a family of functions, or a basis (usually non-orthogonal) in Hilbert space 7i. If there exist two constants A and B with 0 < A < B < oo such that for any function / in Ti Chapter 2. THE THEORY OF WAVELET TRANSFORM ^ll/H <El(/»^)| <5||/|| , j 2 2 2 (WT) for all j, where (•, •) denotes inner product, | | / | | = |(/, / ) | , the basis 2 12 (2.1) = 1,..., J} comprises 2 a frame, and A and B are called the frame bounds. The frame is called a tight frame if A = B. For a tight frame, £l(/,Vv)| = ^ll/H j 2 (2.2) 2 In this case, the function can be projected onto a tight frame and reconstructed perfectly from the projected coefficients, i.e., / = ^ E ( / - ^ 3 (2-3) The constant A in this case is called the redundancy rate. If A — 1, and {tpj} is a tight frame, {ipj} is a n orthonormal basis. In the discrete case, the basis {ipj} is a set of vectors. The decomposition and reconstruction of the signal, vector f, can be written in matrix form as f = Ff, (2.4) f = jF f, (2.5) H where f and f are the wavelet coefficients and the reconstructed signal respectively, and the superscript H denotes the complex conjugate transpose. The following simple example can be used to show the form of a frame and the meaning of redundancy. The basis consists of three vectors e\ = (0,1), e-i — ( — | ) and e = ( ^ , | ) and the frame 3 operator can be written as Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 13 0 1 V \/3 2 1 2 \/3 2 1 2 (2.6) It is clear that j^F l? = I, so F is a tight frame with the redundancy constant A = |. H Here we have used three basis vectors to represent a 2-D vector which can, of course, be represented by the orthonormal, non-redundant basis e = (0,1), e = (1,0). It can be x 2 seen by extrapolation that a four vector frame can be chosen with a redundancy of | , and a five vector frame with a redundancy of | , and so on. From a different viewpoint, a 2-D vector can be decomposed into a linear combination of two orthogonal vectors, or two, three and so on non-orthogonal vectors. These non-orthogonal vectors construct nonorthogonal bases. These decompositions are less efficient than the orthogonal basis but may provide useful information. The original vector can be reconstructed using equation (2.5). If A ~ B, the frame is called a 'snug' frame. For a snug frame, we have - F F f + Rf, H A+ B ' < > 27 where the residual Rf is very small. In this case, a signal can only be reconstructed approximately using equations (2.4) and (2.5). Wavelet frames can be treated as adjoint operators. An adjoint operator F of an operator G is defined as (F c,f) = (c,Gf), H (2.8) where F and G are two matrix operators and c is any column vector with the number of elements equal the number of columns of matrices F and G. The operator F is called a Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 14 self-adjoint operator, if F = G. The Gabor and Morlet wavelet frames which will be used in the next two chapters are both self-adjoint operators. In most cases in this thesis, the self-adjoint operators are not orthonormal, and consequently P ff p = i+J* ( I _ R ) > (2 .9) where R is the residual matrix. Many geophysical operators, for example, wave equation modeling and migration, forward and inverse r — p transforms, are also adjoint operators (Claerbout, 1992). The wavelet frames are no different than these operators except that different properties are examined. Adjoint operators are usually easy to construct, but in order to make the residual matrix R less significant and therefore the frames snugger, the redundancy of the frames has to be increased. The dual frame of F = = 1,..., J} is Fd = {ipj, j = 1,..., J} with F F H D = FF H D = I. (2.10) The dual of the Gabor frame will be described in the next section. In many situations, the dual frame is difficult to obtain, or has an undesired structure. So the self-adjoint frame is often used to recover the data. Next, the details of the Gabor and wavelet transform frames will be examined. 2.2.2 Gabor Frames In this subsection, two special forms of the Gabor frames, the dual frame and the selfadjoint frame, are discussed. The functions making up the dual frame are infinite in length which make them unattractive in real situations. To overcome this problem, a least squares approach is developed in this thesis to construct a special type of frame, a Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 15 localized approximate dual frame. Gabor frame and its dual The Gabor basis function can be written as the Fourier basis multiplied by a window function g(t) at time r: ^ r ( t ) = e-^g{t - r). (2.11) The window function for the Gabor transform is 9(t) = V2 L 1 * (2.12) exp ' where the constant L controls the shape of the window function. The Gabor transform maps a 1-D signal to a 2-D space in time and frequency. The sampling of this 2-D domain is arbitrary but will affect the accuracy of the reconstruction. Let D and W be the sampling intervals at time r and frequency u>. The Gabor coefficient at frequency to = mW and time r = nD can be written as c mn = J f(t)g(t - nD)exp(-imWt)dt. (2.13) The above equation can be considered as a decomposition of a function into a number of elementary functions g(t — nD) exp(-imWt). These elementary functions are not orthogonal to each other, so the reconstruction cannot be done using the same basis. Bastiaans (1981) showed that the inverse transform can be found by introducing an auxiliary function 7 with J 7(0 g(t - nD) exp(-imWt) dt = 8 8 . m The signal can then be reconstructed from the coefficients by n (2.14) Chapter 2. THE THEORY OF WAVELET TRANSFORM /(<)= +00 +00 E E (WT) 16 c (t-nD)exp(imWt). (2.15) mnl n= —00 m=—00 The auxiliary function obtained by solving equation (2.14) is * ) = {je> It) t E J (n+l/2)>|t/D| .(-i) B e ^ -7T 'n + 1 2 (2.16) where iiTo = 1.8540746 is a normalization factor (Bastiaans, 1981) and the constant L = D is used for the window function in equation (2.12). Figure 2.1 shows the shape of the auxiliary function. In the discrete case, the number of data is K = MN where M and N are the number of frequencies and locations of the Gabor coefficients. Equations (2.13) and (2.15) can be directly applied to analyze and reconstruct the signal using the explicit forms of g(x) and 7(3;) (Wang et al., 1994) and the computations can be speeded up by using the Zak transform (Daubechies, 1990). A few remarks should be made regarding the Gabor transform with the above dual frames. First, the function g in equation (2.13) and the function 7 in equation (2.15) are interchangeable. In other words, the forward GT can be computed using 7 as the window function and the inverse transform can be evaluated using g as a substitute for the auxiliary function. Second, since the number of frequencies M is much smaller than the number of data, the resolution in frequency is not as good as in the Fourier spectrum. It is this sacrifice of frequency resolution that results in a gain of the resolution in time. Third, since the function 7 is not localized, all the coefficients are used to reconstruct a single sample of the signal. This means that frequency filtering cannot be done for a certain segment without affecting the rest of the signal since the Gabor coefficients of this segment will also be used to construct the remaining part of the signal. The use of Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 17 finite support auxiliary functions will be discussed later in the next section following a description of the matrix forms of the Gabor frames. The matrix forms of Gabor frames In the discrete case, the forward transform, or decomposition, of signal f using frame F can be written as f =Ff (2-17) fmn = J2k=o fk 9k-nD exp(-imWk), where m = [0,M — 1], n = [0, N — 1] with M and N the number of frequencies and locations respectively, and k = [0, K — 1] with K the number of data. The inverse transform, or reconstruction, of the signal can be written as f =F f H (2-18) fk =Y,n=0 £ m = 0 fmn Ik-nD exp(imWk). The structure of F in the forward transform and of F in the inverse transform are controlled by the sampling of the Gabor space, i.e., the number of locations N and of frequencies M, and by the window function 7. For the dual frame, the total number of parameters is MN = K and the window function in the inverse transform is the auxiliary function j(t) in equation (2.16). The form of F = F^ is therefore very different than that of F. For the self-adjoint frame, MN » K and the window function in the inverse transform is 7 = g, and F is simply the complex conjugate of F. Figure 2.2 shows the sampling of the Gabor frame and its dual in the case when K = 256 with N = 16 locations and M — 16 frequencies at each location. Figure 2.3 shows the self-adjoint frame with K = 256, where the frame is oversampled by a factor of 8. Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 18 Approximated auxiliary function using least squares In seismic data processing, it is generally desired to decompose and reconstruct the seismic data locally, so that the localized, or time-variant information can be incorporated into the processing. However, as discussed previously, using a finite support basis to decompose and reconstruct a signal usually requires the redundant sampling of the transformed space. For example, in case of the Gabor transform with 7 = g, we need to oversample the Gabor space, i.e. MN >> K, to achieve a good reconstruction. This subsection tries to find an alternative approach to compute function 7 which can be used to approximately reconstruct the signal when MN is only several times larger than K, i.e., only using a small oversampling rate. The above problem can be reformulated to find an auxiliary function h with finite length L which approximately satisfies equation (2.14): w K-l Yl k 9k- D exp(-imWk) ~ 8 8 fc=o h n m n (2.19) for any location n G [0, N — 1] and frequency m G [0, M — 1]. The above equation can be rewritten as: ~ 68 EkJo h g - cos(mWk) k K 1 Efc=o k nD m x h gk-nDStn(mWk) ~ 0 n ( 2 - 2 0 ) k Now the problem is to find a finite length window function h which satisfies the above equations with minimum least squares error. From Figure 2.4 we see that if the above equations are satisfied in the region between the two vertically dotted lines, they will also be satisfied in the whole region. Assuming that the inverse window series has L w samples, the above equations can be expanded into Chapter 2. THE THEORY OF WAVELET TRANSFORM Lw (WT) —1 hk gu ~ 1 for n = 0; m - 0, fe=o Efc=o k 9k-nD cos(mWk) ~ 0 19 (2.21) _1 h I Ejfeo 1 ^ ^-nfl /or i/ie re.s£ o/ n and m. (2.22) «n(mWfc) ~ 0 In the case of Figure 2.4, the location numbers n = 0 , ± 1 , ± 2 are related to a certain location sample in the Gabor space, so only five windows are considered. Because cos(mWk) = cos[mW(M — k)] and sin(mM^A;) = sin[mW(M — k)), each set above only has M/2 independent equations. There is a total of hL independent equations with L w w unknowns. The problem is overdetermined and the coefficients h can be found by least squares. Figure 2.5 shows the optimal auxiliary window function. Increasing the number of locations increases the number of equations in the system and makes the problem more overdetermined. The shape of the inverse window function becomes closer to Gaussian, and this approximate dual frame closer to the self-adjoint frame. Decreasing the number of locations to the number of data samples results in a well posed system of equations and leads to Bastiaans' equation (2.16). The features of the above three different types of frames for the Gabor transform are summarized in Table 2.1. Dual Self-adjoint Approx. Dual Sample # MN = K MN » K MN > K Basis length infinite finite finite Table 2.1: Comparison of the three types of Gabor frames. Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 2.2.3 20 Wavelet Frame The continuous wavelet transform decomposes a signal f(t) into a linear combination of a wavelet basis which can usually be represented analytically. The basis comprises a fixed waveform that is dilated to different scales a and shifted to different locations b: f(a,b) = 4= / )<&» (2-23) where i\> is the complex conjugate of ij>. The basic wavelet can be of any form, but in order to recover the signal from the coefficients, two admissibility conditions have to be satisfied (Goupillaud et al., 1984; Daubechies, 1992). First, ip(t) has to be square integrable C^= \ W n dco < oo (2.24) where 4> (u) = J iP(t)e- dt FT (2.25) iut is the Fourier transform of ij;. Second, there should be no DC component of ip , i.e., FT / tjj(t)dt = 0. The signal can then be reconstructed from the wavelet coefficients using 0 f(t) = Re ?r f<x poo ^ / 1 / A dn /(a,6)^(—)db-. (2.26) One of the most commonly used wavelets is the Morlet wavelet which will be described in Chapter 4. Figure 2.6 compares the Morlet wavelet basis with the Gabor basis both in the time and the frequency domains. Figure 2.6(a) shows three wavelets with scales a = 0.5,1, 2. For the sampling of scale, two terms, octave and voice are often used. One octave denotes scale from a to 2a, or from 2a to 4a and so on. Voice defines the number of scales between two neighboring octaves. Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 21 The CWT coefficients could be computed by direct implementation of the above equations which can be done in two steps. First, dilate ip to scale moo where a is the 0 scale interval, then correlate the wavelet with the signal. The problem with this approach is that as the scale increases, as shown in Figure 2.6(a), the length of the wavelet increases dramatically and the correlation becomes computationally expensive. For this reason, a fast algorithm called the algorithme a trous or a trous algorithm can be used (Holscheider et al., 1987; Dutilleux, 1987; Shensa, 1992, 1993). The d trous algorithm is designed to avoid convolving the data with long wavelets at larger scales. The left side of Figure 2.7 shows the direct method in which the coefficients at different scales are computed by correlating the wavelet directly with the data. This is very time consuming when the wavelets become longer at large scales. The a trous algorithm avoids the correlation with long wavelets by applying a smoothing filter S and a down sampling operator D. A smoothed data set at scale m is obtained by correlating the smooth filter with the data at scale m — 1. The wavelet coefficients at scale m are then computed by simply correlating the smoothed data with the original wavelet. The down sampling operator is applied to keep less coefficients at larger scales where less spatial resolution is required. 2.3 Orthogonal wavelet transforms A very important discovery in WT studies was that of orthogonal wavelet bases by Daubechies et al., (1986), Daubechies (1988) and Meyer (1990). These localized bases must satisfy certain mathematical conditions and they can be used to transform the signal into the wavelet domain without any loss of information. Since the wavelet bases are orthogonal, the signal can be recovered by means of the inverse WT using the conjugate operator as in the case for the Fourier transform. This section skips the theoretical Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 22 aspects of the orthogonal WT and concentrates on the fast algorithm, known as Mcdlat's pyramid algorithm (Mallat, 1989a, b; Strang, 1989; and Press et al., 1992). In the pyramid algorithm, the WT is performed by a down-sampling from one scale to next larger scale and a convolution of the signal with a pair of quadrature mirror filters (QMFs) at a certain scale. A pair of QMFs consists of a differential and a smoothing filter. These QMFs are convolved with the coefficients at a certain scale to compute the detailed and smoothed coefficients at next larger scale. Figure 2.8 shows the wavelet decomposition and reconstruction of a signal f. The vectors Sj and d j denote the smoothed and the differentiated data, and the subscript indicates the scale. If N = 2 is the number J of elements of the vector f, then d j _ ! and S j _ i represent the smoothed and detailed components of f at the finest scale, respectively. The final decomposition, denoted as f, consists of the following vectors which completely define the signal in the wavelet domain f = tjo> Jo> Jo+i> • • •, dj_i], s where d j is the detailed data at the j d th (2.27) d scale with 2 samples, and S j is the only smooth J 0 coefficient array required at the finest scale jo. In what follows, the oldest and simplest wavelet, the Haar wavelet, is used to demonstrate some properties of the orthogonal wavelet transform. The Haar wavelet is defined by '1 0 < x < 1/2 1>{x)=\-l 1/2<Z<1 -0 (2.28) otherwise. It is possible to compute the Haar wavelet coefficients using equation (2.23) directly, but using the theory of the orthogonal wavelet transform (see Mallat 1989a, b and c; and Daubechies, 1992) a pair of QMFs can be computed and the orthogonal wavelet Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 23 transform can be evaluated by means of the pyramid algorithm. The QMF pair which describes the Haar wavelet consists of a smoothing, or low pass filter [^7j)^7j]> sometimes called the father wavelet, and the differential, or high pass, filter [^j, — ^j], sometimes called the mother wavelet. We now demonstrate, by means of a simple example, the basic operations that are involved in wavelet analysis. Assume that we have a 4 element vector, [2 3 4 5] T which describes some physical property, for example, the velocities in km/sec in a 4 layer earth. Consequently, N = 4 = 2 with J = 2. We assume that the data themselves 2 are at the finest scale, j = 2, and we consider the decomposition at 2 scales, j = 1 and j = 0. The forward WT begins by computing a smoothed and a differentiated, or detailed, version of the original data at scale j = 1- This can be written in matrix form as follows: 1 5 1 I 0 0 2 0 0 1 1 3 0 0 4 -1 5 -1 l 0 -i 0 1 -1 1 9 (2.29) or Li Ha [f] = Sl da where f is the data, L is the matrix composed by the QMF smoothing filter, H is the matrix composed by the QMF differential filter and the subscript denotes the scale. The short bar is used to separate the smoothed and detailed data. To compute the smoothed and detailed components at the next scale, we keep d i and apply the QMFs to S j , where the sampling interval is now twice the original interval. In matrix form Chapter 2. THE THEORY OF WAVELET TRANSFORM 1 72 1 1 1 5 7 1 -1 '72 9 -2 (WT) s 24 0 (2-31) do where s = [7] and d = [—2] are the wavelet coefficients at scale j = 0. s corresponds 0 0 0 to the DC level of the signal in this case. Rather than computing the WT in two steps, as expressed by equations (2.29) and (2.31), a one-step procedure can be implemented as follows: 1 2 1 2 1 2 2 7 1 2 1 2 1 2 3 -2 0 4 I 1 1 0 0 1 v/2 1 V2 5 So (2.32) v/2 1 L v/2 J or f - Wf. (2.33) The transform is orthogonal and complete with W W = I, W T T = W H The data can be recovered by f = W f. (2.34) T Equations (2.33) and (2.34) are used as the notation for the WT in this thesis but they are rarely used in practice because they are expensive to evaluate. The pyramid algorithm, on the other hand, requires about the equivalent computing effort to calculate the WT as does the F F T if short support wavelets are used (Strang, 1989; Press, et al., 1992). Different types of QMFs, as listed in Daubechies' book, can be used just as for the Haar WT. For Daubechies' 4 coefficient wavelet, [ci, c , c , c ] = [1(1 + V3),l(3 + 2 3 \/3), |(3 — Vd), |(1 — v ^)], the matrix form can be written as 7 4 Chapter 2. THE THEORY OF WAVELET TRANSFORM C Ci c 0 2 C 0 c 3 Cl c 2 c 3 c c 0 Cl c 3 —C 3 2 Co C C2 2 25 (WT) c 3 Ci (2.35) Ci —c 0 c 3 -c 2 Ci -c 0 c 3 "C 2 Ci —c 0 c Cl —c 0 3 -c 2 The pyramid approach in Figure 2.8 can be used in the same manner as for the Haar wavelet to compute the wavelet transform coefficients in the form of equation (2.27). Each coefficient computed using the pyramid algorithm described corresponds to a certain wavelet. Although we are not dealing with the wavelet directly, but rather with scales and the QMFs, the wavelet can be obtained by applying an inverse transform to a particular coefficient. Figure 2.9 shows two Daubechies wavelets. So far, only the 1-D transform has been studied. For a 2-D image, we can apply the QMFs to both the rows and columns of the image. A cartoon diagram in Figure 2.10 shows the procedure for a 2-D WT. At each scale, three quarters of the detailed coefficients are kept and one quarter of the smoothed coefficients is used to compute the wavelet coefficients at the next scale. In the diagram, L and H are the QMFs operating c c on the columns, and L and H are those operating on the rows. r r Chapter 2. THE THEORY OF WAVELET TRANSFORM 2.4 (WT) 26 Summary This chapter reviewed some of the theoretical aspects of the wavelet transform. Nonorthogonal WTs such as the Gabor and Morlet transforms were studied using frame theory and orthogonal WTs such as the Haar and Daubechies transforms were described by the pyramid algorithm via the implementation of quadrature mirror niters. The forward and inverse Gabor frames can be constructed in three different pairs, namely, dual frames, self-adjoint frames and least squares approximate dual frames. 1. Dual frames. The forward frame is sampled with MN = K with M, N and K as the number of locations, frequencies and data respectively. The inverse frame is constructed using an infinite length auxiliary function. The advantage of using the dual frame is that with the same number of samples, the data can be reconstructed exactly from their coefficients. However, there are some problems. First, the form of the auxiliary function used to construct the dual frame has infinite length and has no physical meaning. Second, the numbers of frequencies and locations are limited, implying that the frequency and location resolutions are poor. 2. Self-adjoint frames. For self-adjoint frames the inverse transform frame is the (complex conjugate) transpose of the forward frame. The frames are easy to compute, but the signal can only be recovered approximately. In order to achieve a good approximation the transformed space has to be over-sampled, thus increasing computation time and storage space for storing the coefficients. Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 27 3. Approximate dual frames. An finite support inverse window function was computed to construct an approximate frame of the dual frame used for the inverse Gabor transform. The approximate frame is less redundant than the self-adjoint frame. A simple least squares approach was developed to compute the window function for the approximate inverse Gabor transform in the discrete case. The construction of another type of frame, the continuous wavelet transform (CWT) frame was also discussed. The CWT frames can be directly implemented using the correlation of the dilated wavelets with the signal but is very expensive because the correlation of the signal with a long wavelet at large scale for CWT is very time consuming. For this reason, the d trous algorithm was developed by Holscheider et al. (1987) to speed up the transform. Like the Gabor transform, the CWT is also redundantly sampled in order to preserve information. In section 2.3, the pyramid algorithm (Mallat, 1989a) for the computation of the orthogonal WT was described. The transform can be represented by a matrix W with a very important feature: W = W . The applications of WTs in the modeling of well r - 1 logs and in inversion will be explored in Chapters 5 and 7. Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 28 Figure 2.1: The auxiliary function (solid line) in equation (2.16) for a dual Gabor frame. The dotted line is the Gaussian window function of the Gabor frame in the forward transform in equation (2.12). Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 29 (a) / n=16 n=2 /\ =1 -100 i 50 0 -50 i 150 i 200 i 250 i 300 V ii lV !• 100 150 i 200 i 250 300 i 100 i 350 (b) n=2 ILJLJL. n=1 JL JLJL V J II ! If V i V i If (i -100 - 5 0 V r 0 «v ii iV 50 l i L i i 350 Figure 2.2: The window functions in forward and inverse dual frames, (a) The Gaussian windows in the forward frame; (b) the auxiliary window functions in the inverse frame. There are N = 16 locations and M = 16 frequencies at each window. The total number of samples is K = 256, that is, no redundancy is present. The signal exists in between the dotted vertical lines. The samples outside the dotted lines are assumed to be either zero or the periodic extension of the signal. Chapter 2. THE THEORY OF WAVELET TRANSFORM 30 (WT) n=65 ^ n=5 n=1 _| -50 _y L 0 1I I1 I1 I1 50 100 150 200 1L_ 250 11 300 Figure 2.3: The sampling of self-adjoint frames. In this case, the number of locations and frequencies at each window are N = 65, M = 32. The Gabor frame is 8 times oversampled compared with the original data space which has K = 256 samples. Figure 2.4: Correlated windows in the Gabor transform. There are five locations in the Gabor space which affect the inverse of the samples between the dotted lines. A finite length auxiliary window function can be computed to minimize the error in the signal reconstruction caused by the non-orthogonality. See Section 2.2 for details. Chapter 2. THE THEORY OF WAVELET TRANSFORM 31 (WT) 0.4 _Q 21 1 -23.5 1 1 -15.5 -7.5 1 1 0 7.5 1 15.5 I 1 23.5 Figure 2.5: The approximate auxiliary function using least squares for the case in Figure 2.4. Increasing the number of Gabor windows will cause the shape of the auxiliary function to become closer to the Gaussian window as shown by the dotted line. See Section 2.2 for details. Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) (a) (b) Frequency Frequency 32 Figure 2.6: Comparison between the real parts of the complex (a) Morlet and (b) Gabor wavelets. The Morlet wavelets have fixed waveforms and Gabor wavelets have fixed lengths, (c) and (d) show the frequency responses of the wavelets in (a) and (b) respectively. Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 33 SD GM SD Figure 2.7: (a) The direct implementation and (b) the a trous algorithm for the continuous wavelet transform (CWT). In (a), the operators G\, Gi, • • • are the correlations of the signal with wavelets which vary at different scales. The wavelets at large scales are very long, hence the correlations are very expensive. In (b), only three operators: the smoothing operator S, the downsampling operator D and the wavelet operator G are used and they are fixed at all the scales. Since S and G are both short, the algorithm is fast. Figure 2.8: WT decomposition and reconstruction using the pyramid algorithm. See section 2.3 for notation. Chapter 2. THE THEORY OF WAVELET TRANSFORM (WT) 34 Daubechies wavelets Wavelet Coefficients 0.4 0.2 -4 CD TJ _9 4 0 "Q. E -0.2 CO < -8 -10 -0.4 0 0.5 Location -0.6 0.5 Location 0 Figure 2.9: Two Daubechies wavelets, (a) Two wavelet coefficients, and (b) the corresponding wavelets. The total number of data is 1024 = 2 . The finest scale is 9 and has 2 coefficients. The two wavelets are at scale 6 and 4. The location has been scaled from 0 to 1. 10 9 Chapter 2. THE THEORY OF WAVELET TRANSFORM Image 3 H o (WT) 3 H o s © (1) W T C at J - l LrLc I rows lumns Lcl 35 HrLc I (2) He I L r H c I \ HrHc I A I (4) X x HrLc I Y D Y LrHcI D HrHc I Figure 2.10: 2-D orthogonal wavelet transform. At step (1), the QMFs, i.e. the low pass filter L and high pass filter H , are applied to the columns of image I. Then at (2), the QMFs L and H are applied to the rows. The detailed coefficients L H I , H H I and H L I at scale J — 1 are kept and the smoothed coefficients L L I are further decomposed into the detailed coefficients X, Y , D , and a smoothed image at scale J — 2. The final wavelet coefficients are shown as I. c c r P C r r r c c r c Chapter 3 T H E G A B O R T R A N S F O R M IN SEISMIC FILTERING 3.1 Introduction It is well accepted that geophysical data are time-variant, so windowed analysis and processing techniques have been applied to solve various problems such as surface wave velocity analysis (Dziewonski et al., 1969), time-variant deconvolution (Clarke, 1968; Griffiths et al., 1977) and filtering (Scheuer and Oldenburg, 1988; Park and Black, 1995). In these applications, time-variant information is applied in the time domain directly. However, to directly design a filter with time-variant information is difficult especially for 2-D data. This chapter introduces a technique to filter noise in the Gabor transform (GT) space where the time-variant information can be easily incorporated. Since the Gabor transform is a non-orthogonal transform, some information will be lost by using the transform, but this loss can be controlled to a proper level by the use of frame theory. In the next section, an algorithm to compute the Gabor transform is introduced, two synthetic examples are used to show the necessity of GT filtering, and the relationship between GT filtering and Fourier filtering is discussed. When the time-variant information is removed from the GT filter, the proposed approach is the same as applying a Fourier filter with a taper. This relation is studied in section 4 where the 2-D Gabor transform is used to attenuate 2-D coherent noise such as ground roll and head waves. 36 Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 3.2 37 Filtering using the Gabor transform (GT) 3.2.1 Algorithm Traditionally, a seismic trace is assumed to be the convolution of a seismic wavelet with the reflectivity sequence. Assuming that the reflectivity is stationary and that the shape of the seismic wavelet is time invariant, the seismic trace itself is also stationary. In reality, however, neither the reflectivity nor the seismic wavelet are time invariant and, consequently, the Fourier transform cannot be legitimately applied. For this reason, the short time Fourier transform, in this case the GT, is used to analyze the characteristics of signal and noise, and to separate noise from signal. The Gabor transform of a discrete signal f of length K can be computed in two steps as shown in Figure 3.1. First, we compute y = Gf which maps f to a M x N element vector y. Each of the M segments of length N is a truncated signal of f multiplied with a window function g. Second, a fast Fourier transform is computed for each segment of y. The second operation is represented by f = Fy where F is a block diagonal matrix as shown in Figure 3.1. The above steps can written in a matrix form: f = FGf, (3.1) where f is the GT of signal f with M locations and with N frequencies at each location. The vector f can be arranged as a two dimensional M x N matrix and plotted in time and frequency to show the time variant features of the signal. The dimensions of G and F depend on the sampling of location and frequency. Increasing the number of the block diagonal matrices increases the number of locations. We can also increase the size of the block diagonal matrices by padding with zeros or using a longer Gaussian window function. Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 38 According to frame theory described in the previous chapter, the inverse transform can be written as f = G F f = G F*FGf, r H (3.2) r where T denotes transpose and H denotes complex conjugate transpose. Because each small block matrix in F is orthogonal, F is an orthogonal matrix, so F F = I. From H equations (3.1) and (3.2), we obtain f = G Gf. (3.3) T The accuracy to which the signal can be reconstructed depends on how close G G is r to the identity matrix. The range of G G are the bounds of the Gabor frame. More r detailed studies can be found in Daubechies' book (1992). We can design a Gabor filter h(r, OJ), or H in matrix form, to attenuate noise in the Gabor space. The filtered data can be written as f = G F^HFGf. T (3.4) Since all the matrices in above equation are of block form, there is no need to compute the matrices explicitly and the computation can be done block by block to save storage space. 3.2.2 Filtering random noise for signal with localized characteristics Figure 3.2 shows a synthetic trace consisting of two Ricker wavelets with different main frequencies and with additive random noise. The second and third rows of the figure show the Fourier amplitude and phase spectra respectively. Due to the nature of the Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 39 Fourier transform, differences between the two wavelets cannot be seen in the Fourier spectra. Figure 3.3 from the top shows the Gabor amplitude and phase spectra for the data, signal and noise respectively. The difference between the two Ricker wavelets can be clearly seen in the amplitude spectra. Figure 3.4 shows the filtered Fourier spectra and the estimated signals using two different approaches. In the first approach, a bandpass filter is applied in the Fourier space. The filtered spectrum and the estimated signal are shown in the middle row. The second approach considers the spectrum with amplitudes less 2cr as noise and these n coefficients are set to zero. The filtered spectrum and the estimated signal are shown in the bottom row. The results using two Gabor space approaches are shown in Figure 3.5. Compared with the second Ricker wavelet, the first one has a higher frequency and short time duration, thus its spectrum is wider and centered at a higher frequency. In the first approach, two different bandpass filters are applied to different portions of the data. The filtered Gabor coefficients and the estimated signal are shown in the top row. The second approach considers these components with amplitudes less 2<r as noise and these n components are set to zero. The filtered Gabor coefficients and the estimated signal are shown in the bottom row. It is clear that both results have been improved in comparison with the results in the previous figure using conventional Fourier filtering approaches. 3.2.3 Attenuation of the coherent noise In this example, a synthetic trace with two Ricker wavelets is combined with a linear chirp. The chirp is treated as coherent noise. Figure 3.6 shows the data, the Fourier spectrum and the Gabor coefficients. The main frequencies of the two Ricker wavelets are both centered at 35Hz and the frequency of chirp varies from 25 to 0Hz. Because the Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 40 signal and noise are located in time, they are spread in the Fourier space, thus applying a simple bandpass filter cannot achieve a good signal and noise separation. Also, since the basis is of infinite length, although there is no noise in the first portion of the trace, the filtering will cause some damage to the first Ricker wavelet. Figure 3.7 compares the original signal with the estimated signals using Fourier and Gabor filtering approaches. The Gabor technique is clearly superior to the Fourier technique because it attenuates the noise with additional location information about the signal and noise. In this case, a better estimate of the second wavelet is obtained without affecting the first wavelet. 3.2.4 Gabor filtering of real data The Gabor filtering works well for synthetic data where the localized characteristics are very distinct. For real data, localized features may not be as apparent as in the synthetic examples. Figure 3.8 shows a real trace and its Gabor and Fourier spectra. After 0.4 seconds, the trace is dominated by ground roll. Since the ground roll occurs after the first break, it can be attenuated using the Gabor technique without affecting the first break as well as near surface reflections. As will be shown in the next section, without the spatial (time) information about the signal and noise, i.e., applying a frequency bandpass filter in the Gabor spectrum at all locations, the Gabor technique is similar to bandpass Fourier filtering. Incorporating spatial information in the Gabor domain could improve the type of Gabor filtering, but this Gabor filtering is more expensive than Fourier frequency filtering, and improvement has to be evaluated with this consideration in mind. Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 3.3 41 Relationship between the G T and F T filtering As has been shown in the previous section, adding time-variant information improves filtering results. This section shows that when the spatial information is not applied in the Gabor filtering, i.e., h (r,u>) = h (u>), the effect of the Gabor filter is equivalent GT GT to that of a Fourier filter with a Gaussian taper. For simplicity, the continuous case is considered. Assuming f(t) are data, we can define the Fourier transform as f {u>) = f f(t)e-^dt. (3.5) FT Suppose h (uj) is a frequency filter, thus the filtered data are FT = J h^Mf^Me^du,. f (t) FT (3.6) In the continuous case, the forward Gabor transform can be defined as f (r, ) = J f(t)g(t - ry-^dt, GT (3.7) U and the inverse transform is given by f(t) = U ^ ( T , ^ - The Gabor filtered data with a time variant filter h f (t) GT =U (T,UJ) t fe (r,a;)/ (r,a;) (f GT 5 FT are - r)e^drdu:. Gr We now look at the relationship between f {t) (3.8) r)e^drdw. and f (t) GT assuming a Gaussian window function g(t) = je~d) / 2 2 (3.9) when h (r,uj) = h (u>), GT GT with L a constant. We begin by applying the time-variant filter h (r,u>) to a single sinusoid u{t) = e"" . GT Using equation (3.7) the GT of u(t) can be written as ot Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING U (T,U) = GT - « ( ^ o ) r - ( % - . ) ) / 2 (3 2 e e ] In the case of h (r,u>) = h (u>), by applying the GT filter to GT 42 GT U(T,OJ) 1 0 ) and then applying the inverse GT using equation (3.8), we obtain the filtered sinusoid, u (t) = e*'-* j h {u)e-( ("-"°» ? dw. GT GT 2L 2 (3.11) 2 For any function f(t) = J f(u> )e dw , the filtered signal is lWot 0 = J e^J J^{t) 0 h (u,)e-( ("-"°» ? c[^du, . GT 2L 2 (3.12) 2 0 The above equation can be rewritten as f<ZT(t) = J J h {u; )e-W»-»^ l fa du, GT (3.13) 2 2 0 0 by interchanging UJ and UJ . 0 We conclude that applying a time invariant filter Ii (uS) in the Gabor domain correGT sponds to using a filter h (u;) = J h (uJo)e-W"-"°^ du, FT GT (3.14) 2 0 in the Fourier domain which is the GT filter convolved with a window function. In the discrete case, since the sampling of the Fourier and Gabor spaces are usually different, it is not simple to show that a similar relation to equation (3.14) holds. Consequently a real data example is used to show that applying a time-variant filter h GT in the Gabor domain is equivalent to applying h (u>) * e~( ( ')) / , a tapered band pass GT 2L a 2 2 filter, in the Fourier domain. Figure 3.9 shows: (a) the Gabor filtered data with the same bandpass filter applied to all locations; (b) the Fourier filtered data with the same bandpass filter and (c) the Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 43 same filter with a taper. In (a) and (c), the results are relatively similar and support the previous conclusion. Of course, the time variant information such as the start time of the ground roll can be easily incorporated into the Gabor filtering thus improving the result as shown in next section. 3.4 Gabor filtering for 2-D data 3.4.1 Reasons for a 2-D G T In reflection seismology, coherent noise such as ground roll, refracted waves and multiples must be removed in order to obtain a good reflection image. The differences between the reflections and coherent noise are, in general, variable from one location to another. Also, the coherent noise such as ground roll is usually localized in the seismic data. Conventional frequency and f-k filtering cannot incorporate these localized differences. In this section, a 2-D Gabor filtering technique will be used to incorporate the localized information. By applying the GT in both the x and t directions, the 2-D Gabor transform maps data from a 2-D to a 4-D space. In this 4-D space, a 4-D filter can be applied to incorporate the localized differences between the reflections and coherent noise leading to an efficient attenuation of the noise. Figure 3.10 shows a real shot record and Figure 3.11 shows its f-k transform. The difficulties of attenuating ground roll using the conventional f-k method become apparent by examining these figures. First, the signal and noise are not linear events as clearly indicated by the data in the four boxes. Secondly, the signal and the coherent noise are localized. The head waves and near surface multiples are located at smaller times; the ground roll at larger times. Since the 2-D Fourier transform basis functions are of infinite dimension (or in the discrete case, the dimension of the whole image), this location effect will blur the information in the f-k domain, therefore the signal and noise cannot be well Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 44 separated. However, the above problems can be solved by applying the Gabor transform where a localized basis is used. Figure 3.12(a) shows two Gabor basis functions which are compared with a f-k basis function in Figure 3.12(b). 3.4.2 Formulation and algorithm The 2-D Gabor transform maps the 2-D image to a 4-D space and can be written, in continuous form, as S(k,x -uj,t ) = j j s{x,t) e- ^ e{ 0 )2/2 e- -^ e- ikx {t 0 fl2 iwt dxdt (3.15) where s(x,t) and S(k, x$; u>, t ) represent a seismic image and its Gabor transform; x 0 0 and t are the locations in time and offset. The constants L and L control the shape 0 x t of the Gaussian window. By choosing reasonable sample intervals and the constants L x and L , the filtered image can be recovered using the following equation t H(k, x ; u, t ) S{k, x ; w, t ) - ^r e ( 0 0 0 0 e )2/2 ikx e~ ^ e { )212 iwt dk dx dw dt 0 0 (3.16) where A and B are the frame bounds. In the discrete case, the transform can be computed using algorithm similar to the 1-D case. It is obvious that when H — 1, we will simply recover the signal and the accuracy depends on the frame bounds A and B. When H(k,x ;u>,t ) — H(k,u>), i.e., without 0 0 using any location information, the result will be similar to that of f-k filtering. This has been noticed by Womack and Cruz (1994). We can extend relation (3.14) to 2-D where applying H in the Gabor domain is equivalent to applying a Fourier filter which is H(k,u>) convolved with a 2-D window function. We can then predict that with location Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 45 information the GT filtering will improve the result of a simple f-k filtering. This is now demonstrated using the data in Figure 3.10. 3.4.3 Example The most difficult task in filtering ground roll and other coherent noise is that the signal and ground roll overlap at low frequencies. The conventional f-k or frequency filtering cannot overcome this problem because applying a low frequency reject filter will result in the loss of some valuable low frequency information. F-k filtering cannot achieve a good result because of the localized nature of the signal and noise as mentioned previously. In what follows, I will show how the location information concerning the signal and noise can be incorporated in the 2-D GT domain to greatly improve the attenuation of the ground roll and head waves. Step 1: Decomposition into frequency bands using the 1-D Gabor transform. In this case, data S have been decomposed into three frequency bands: low frequency band (0,25] Hz; middle band (25,60] Hz and high frequency band over 60Hz where ( means excluding the number on the right hand side and ] means including the number on the left hand side. The frequency bands are labeled 51,52 and 53 in Figure 3.13. Note that the amplitudes of the low frequency band 51 (in the range of ±5000) are much larger than those of the second and third bands, which are in the ranges of ±1000 and ±500, respectively. This means that although the low frequency signal in SI is not clearly shown, it may have relatively large amplitudes and therefore should not be ignored. Step 2: Recovery of the low frequency information. A 2-D Gabor velocity filter is applied to recover the signal. Although the frequencies of the signal and ground roll overlap, the dips (or the velocities) of the ground roll and signal are quite different. Figure 3.14 shows the recovered low frequency signal by rejecting velocities less than 8000m/.s in the region where the ground roll exists. The Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 46 estimated reflection section with the head waves and surface multiples is shown in Figure 3.15. Step 3: Attenuation of the head waves and their surface multiples. This is done by applying a velocity filter to reject velocities of 2200 to 3200/, m/s at times 0 to lOOOm^ec (samples 1 to 500 in the diagram). The estimated head waves, their surface multiples and the estimated reflected waves are shown in Figure 3.16. For clarity, the notations used in Figures 3.10 to 3.16 are listed below. • S: raw shot; • SI: subband [0, 2S\Hz obtained using 1-D GT; Sis: recovered low frequency signal using 2-D Gabor filtering; 51ra: low frequency noise 51 — Sis; • 52: subband (25,60]#* obtained using 1-D GT; • 53: subband over 60Hz obtained using 1-D GT; • 55: data after the attenuation of low frequency noise: 5 — Sin; • Srfr: estimated refraction and multiples in 55 using 2-D GT; • 555: estimated reflected waves: 55 — Srfr. The f-k transforms of 55 and 555 are shown in Figure 3.17. Comparing with Figure 3.11, the f-k transform of the original image, we find that it is difficult, if not impossible, to directly design any f-k filter that attenuates the coherent noise and achieves the same results as the 2-D Gabor filter. Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 3.5 47 Summary Section 3 showed the correspondence between Gabor and Fourier domain filtering: a bandpass Gabor filter without any spatial information corresponds to a bandpass Fourier filter obtained by convolving the Gabor filter with a window function. When the spatial information about signal and noise are available, they can be used to improve the filtering results using the Gabor approach. Synthetic tests showed that the results using the Gabor filter with spatial information yields a significant improvement to those of the conventional Fourier filtering. The 2-D Gabor filter is used in section 4 to attenuate ground roll, refraction and multiple waves using appropriate location information. The improvement of the GT filtering can be assessed by examining the f-k transform of the estimated signal which would be difficult, if not impossible, to obtain using conventional f-k filtering. Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 48 MN d'ag(g) F T m=l K diag(g) m=2 F T = MN diag(g) m=M F T A f G Figure 3.1: Matrix structures of the Gabor transform and Gabor coefficients. The Gabor transform of a vector f can be computed by multiplying f with a window and resampling matrix G and a blocking Fourier transform matrix F. Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING I 0 49 I 1 Time (sec) 2 0 1 Time (sec) 2 100 50 Freq (Hz) 100 0 50 Freq (Hz) 100 100 50 Freq (Hz) 100 50 Freq (Hz) 100 1 Time (sec) 2 50 Freq (Hz) 50 Freq (Hz) 0 Figure 3.2: Synl: synthetic trace of two Ricker wavelets and random noise. Top row: the synthetic trace (left) which consists of two Ricker wavelets (middle), with different central frequencies, and Gaussian noise (right). They are denoted as data, signal and noise. Middle row: the amplitude spectra of data, signal and noise. The time interval is 2ms so the Nyquist frequency is 250Hz but only up to 100Hz is displayed. The horizontal line on the right indicates the stand deviation of the noise, <r . Bottom row: the phase spectra of data, signal and noise. n Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 50 _80 N cr £40 ll 0 0.5 1 1.5 Time (sec) 0.5 1 1.5 Time (sec) 20 2 0 0.5 1 Time (sec) 1.5 0.5 1 1.5 Time (sec) 0.5 1 Time (sec) 2 80 N £60 cr £40 Ll_ 20 0.5 1 1.5 Time (sec) 1.5 Figure 3.3: Gabor transform of synl. From the top, amplitudes (left) and phases (right) for the data, signal and noise. The frequencies are only displayed up to 100Hz Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING Freq (Hz) Time (sec) Freq (Hz) Time (sec) 51 Figure 3.4: Fourier filtering of synl. Top row: the Fourier spectrum (left) of the original data (right). Middle row: the band passed spectrum and the filtered data. Bottom Row: on the left is the spectrum by setting the amplitude to zero when it is less than 2cr and one the right is the filtered data. n Chapter 3. THE GABOR TRANSFORM 0.5 IN SEISMIC FILTERING 0.5 1 1.5 Time (sec) 52 1 1.5 Time (sec) 2 o 1 t o . - ••• — j L- _ miff E < -1 0.5 1 1.5 Time (sec) -2 0.5 1 1.5 Time (sec) Figure 3.5: Gabor filtering of synl. Top row: the band passed spectrum and the filtered data. Bottom row: the filtered spectrum with threshold 10% of the maximum amplitude and the filtered data. Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 53 Synthetic trace Time (sec) Gabor Coefficients Frequency (Hz) Figure 3.6: From top, Syn2: two Ricker wavelets and a chirp signal; the Gabor spectrum; and the Fourier spectrum. Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 54 Signal I I I I I 0 0.2 0.4 0.6 0.8 1 i 1 : I I 1 1.2 Time (sec) I 1 1 L_ 1.4 1.6 1.8 2 FT filtered with taper 1| 1 I i i i i 0 0.2 0.4 0.6 0.8 1 1 i i 1 1.2 Time (sec) r 1 1 1 1.4 1.6 1.8 L 2 Figure 3.7: From top, The original Ricker wavelets; the estimated Ricker wavelets using Gabor filtering; and using Fourier filtering. Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 55 (a) Raw trace Time (sec) (b) Gabor Coefficients _80£60 £40~ LL 200 0.5 1 Time (sec) 1.5 2 (c) Fourier spectrum 4 i 1 1 1 1 1 1 i r Frequency (Hz) Figure 3.8: From top, tr70, a real trace, and its Gabor and Fourier spectra. The trace is dominated by ground roll from 0.4.sec. 56 Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING (a) G T filter 1 — 0.41 1 1 1 1 1— Time (sec) (c) FT filter with taper 0.41 1 1 | I I 0.2 0.4 0.6 0.8 _ J 1 1 1 1 I I 1 1.2 Time (sec) I I I 1.4 1.6 1.8 L 2 Figure 3.9: Trace tr70 after Gabor and Fourier filtering. ^ (a) Gabor filter = h (u) as box car filter; (b) Fourier filter h same as h (u); (c) Fourier filter with taper. Without the location information the GT filtering is similar to a Fourier filtering with a taper. See section 3.3 for details. I (T,UJ) gt GT FT GT Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 57 0 200 JD Q. 400 E CO w <D E 600 800 1000 20 40 60 Trace number 80 Figure 3.10: Shot2: a real shot record. The four boxes are used to highlight the dip variation of the seismic events. s 200 r 1 1 150 I ^100 0 0.1 0.2 0.3 0.4 k Figure 3.11: F-k transform of Shot2. Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING Trace number 58 Trace number Figure 3.12: (a) 2-D Gabor and (b) f-k bases. Figure 3.13: Shot2, denoted as 5, and its three 1-D Gabor Subbands 51, 52 and 53. 5 = 51 + 52 + 53. The amplitudes of each image have been normalized. The amplitudes of 5 and 51 are in the range of ±5000, much larger than those of 52 in ±1000 and 53 in ±500. Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING S : raw shot 0 20 40 60 Trace number 59 S 1 : Gabor band 0 - 2 5 Hz 80 0 20 40 60 Trace number 80 Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING S1n 0 20 40 60 Trace number 60 S1s 80 0 20 40 60 Trace number Figure 3.14: The estimated low frequency signal Sis and noise Sin for low frequency band SI using a 2-D Gabor transform. 80 Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING 61 0 200 oT 400 Q. E CO a, 600 E 800 1000 0 20 40 60 Trace number Figure 3.15: The estimated signal after subtracting the low frequency noise Sin in Figure 3.14. Chapter 3. THE GABOR TRANSFORM Trace number IN SEISMIC FILTERING 62 Trace number Figure 3.16: Left: estimated refracted waves, Srfr, using 2-D Gabor filtering. Right: estimated reflected waves, SS — Srfr. Chapter 3. THE GABOR TRANSFORM IN SEISMIC FILTERING SS 200 200 63 sss 150 N X O" 100 CD Figure 3.17: F-k spectra of estimated reflected waves with (SS) and without (SSS) refraction. Chapter 4 A T T R I B U T E ANALYSIS USING T H E CONTINUOUS W T 4.1 Introduction Because of the localized nature of earthquake records and seismic wavelets the concept of the wavelet transform has been used to analyze geophysical data for decades (Dziewonski et al., 1969; Morlet et al., 1982; Goupillaud et al., 1984). It has also been widely accepted and studied in other applications such as underwater oceanic acoustic signal analysis (David and Chapron, 1989), speech detection (Lienard and d'Alessandro, 1987), time series analysis (Schiff, 1992), image processing (Mallat, 1989b) and mathematics (Meyer, 1993). In contrast to the Gabor transform where the basis functions are composed of windowed sinusoids with a fixed length, the wavelet basis comprises a single waveform which is dilated at different scales and shifted to different locations. The non-orthogonal wavelet transform studied in this chapter is sometimes used under different names with minor variation in meaning. The continuous wavelet transform (CWT) is commonly used when the analytic form of the wavelet is adopted. Since CWT is often defined using an integral, it is also called the integral wavelet transform, IWT (Chui, 1992a). For computation purposes, the wavelets have to be represented in discrete form and the transform is called the discrete wavelet transform (DWT). An important issue in the case of DWT is how to choose the sampling rates in location and scale. As discussed in Chapter 2, in order to completely transform the information from the spatial domain to the scale-location 64 Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 65 domain, the scale-location space is usually oversampled. If the DWT space is properly sampled, the reconstruction can be quite accurate (Goupillaud et al., 1984). In this chapter, the term CWT will be used in general. If the DWT is used, it is assumed to be sampled properly so it has the same features as the CWT. The goal of this chapter is to understand the features of the CWT, the Morlet WT in particular, and the use of the amplitudes and phases of the complex wavelet coefficients as types of attributes to analyze the time- and scale-variant characteristics of the seismic data. First, the CWT is introduced in terms of the convolution of a signal with the dilated wavelets. In section 3, the Morlet WT is compared with the Fourier transform using a few simple examples. Section 4 analyzes the time- and scale-variant characteristics of different wave types in a shot record using the amplitudes and phases of the Morlet wavelet coefficients. 4.2 4.2.1 The Morlet wavelet transform The wavelet transform in the form of convolution The continuous wavelet transform of a signal f(t) can be considered as a basic wavelet dilated at different scales and then convolved with f(t). First, the basic wavelet can be constructed as with v the main frequency and g(t) a window function. For the Morlet wavelet transform, the window function can be formed by g(t) = ^e-^ (4.2) where L controls the window width. Thus the Morlet wavelet basis is constructed by Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 66 = "TTV-C—), (4-3) where a and b denote scale and location respectively. Denoting Mt) = 4 = l K - ) (4-4) as the dilated wavelet at scale a, the wavelet coefficients at scale a can be written as a function of spatial variable b, fa{b) = f(b)*K(b), (4-5) where * represents convolution and ip denotes complex conjugate of ip. Since f (b) is a a function of both a and b, it can also be written in another form, /«(&) =/(a, 6), {a=.[0,oo); b = (-oo, +oo)}. (4.6) In the continuous case, the signal can be reconstructed exactly from the wavelet coefficients (Goupillaud et al., 1984). This reconstruction can also be written as a convolution, /(*)= ^ReJ °° JZoH^ (t)dadb/a 2 0 (aib) = ±ReSff (b)*MW"/« = RefZ°f (t)da, 2 a (-) 4 7 a where C$ is a constant and the decomposed trace at scale a is /.(*) = - 4 r / a ( * ) * M^. GTC,/, (4-8) In the discrete case, the sampling of a and b cannot be infinitely dense and thus the reconstruction is not exact. The accuracy of the reconstruction depends on the sampling Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 67 rates in scale and location. This accuracy also depends on the sampling methods used in wavelet space (Goupillaud et al., 1984). In this thesis, uniform sampling in both time and frequency are used. The wavelet coefficients at scale a of a discrete signal f are directly computed using equation (4.5) by first computing the discretized wavelet rj> and a then convolving with f. The problem of this direct approach is that the wavelet at a large scale is very long and the convolution becomes very expensive. Since the computing time is not the main concern in this chapter, we will simply use this direct implementation. For large data sets, however, an a trous algorithm (Holscheider et al., 1987; Dutilleux, 1987; Shensa, 1992, 1993) can be used to improve the computational efficiency. The time consuming convolution at large scale can be avoided by the following representation •ft ' )= -fe//W(?)<a , , a 6 vH = (4.9) ^\Jf(at)iP(t- -)dt. b In this case, we can scale f(t) to f(at) at scale a, then convolve f(at) with the same wavelet. 4.2.2 Fourier transforms of the Morlet coefficients and multiscale traces Since at each scale the Morlet wavelet coefficients and multiscale traces are both functions of time, they can be analyzed using the Fourier transform. The Fourier transform of the basic Morlet wavelet in equation (4.3) can be written as V5)H = V ^ ^ - ^ V S and it becomes e^ b (4.10) 68 Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT $™( ) u = ^e- ^- ^/2 L2a2 (4.11) v/a when 6 = 0. The wavelet transform at scale a in equation (4.5) can then be written as /.H^^H^H (4-i2) where^represents two transforms - a wavelet transform then the Fourier transform. The above equation indicates that the Fourier transform of the wavelet coefficients at scale a is the Fourier transform of the signal f (u>) multiplied by the Fourier transform of the FT window function. Similarly, for the multiscale trace in equation (4.8), the Fourier transform of f (t) can a be obtained as = cfcf H\ti H\ = cfcf He- (»- FT FT T 2L2a2 (4-13) 2 v ^74, which is simply f {oj) multiplied by a window function with a center frequency at v/a. FT Equations (4.12) and (4.13) show that the wavelet transform and multiscale trace can be computed using the Fourier transform. However, since the above relations are all derived for the continuous case, the continuous FT and convolution have to be substituted by the discrete FT and circular convolution respectively and these substitutions may cause problems in analyzing the localized characteristics of the seismic data. In the rest of this chapter, the convolution approach will be used to compute the WT coefficients and multiscale traces. Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 4.3 69 Understanding wavelet attributes Like the Fourier transform, the Morlet WT also transforms data to complex coefficients. What makes the WT much more powerful than the FT is that the WT represents the information in both time and frequency or in a localized manner. This section examines the meaning of wavelet attributes - the amplitudes and phases of the wavelet coefficients. 4.3.1 Wavelet attributes of spikes Let us first look at the Fourier transform from the viewpoint of complex signal decomposition. The decomposition of a spike at t can be written as 0 S(t;t ) = J2 ' eiw{t t0) 0 (- ) 4 14 where ; is used to separate the time variable t and parameter to- Figure 4.1 shows the decomposed complex signals and their attributes in time (horizontal) and frequency (vertical). The amplitude is equal to one at all times and frequencies. The phase of u)(t —1 ) for a fixed a; is a straight line in time, and phase lines with different u> intercept 0 at to to indicate the location of the spike. In this case, the phases are wrapped from —7r to 7r, but the phase pattern can still be used to identify the location of the spike. As shown in Appendix A, 8(t;t ) may be approximated using the Morlet basis by 0 5(t;t ) = 0 (4.15) Figure 4.2 displays the same attributes as Figure 4.1 but computed using the WT. The diagram of wavelet amplitudes is also called a scalogram. Note that the phase in the above equation is the same as the phase in the FT decomposition. The scale a is related to frequency u> by u> = ^ where the reference frequency is chosen as v = 60Hz. Since the Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 70 small amplitudes at large \t — to\ are truncated, only part of the phase is kept. As we will see later, it is this truncation that makes the wavelet phase attractive. Both the FT and WT phases contain the location information of the spike. The unwrapped FT phase ui{t — t ) and wavelet phase — to) are plotted in Figure 4.3. The 0 lines in both figures indicate the theoretical phases for frequencies of 20,40,60,80,100 and 120Hz. The horizontal and vertical axes display time and angle respectively. Wavelet phases are the same as the Fourier phases except that they are truncated with different window lengths at different frequencies. All the phase lines intercept at t indicating the 0 location of the input spike. Although the computed phases in Figure 4.1 and 4.2 are wrapped, the phases are consistent at to for all frequencies. Since the phase unwrapping for real data can be very time consuming and often inaccurate, this wrapped phase pattern is used to analyze the location of the signal. For convenience, the state where the phases are consistent for all frequencies is called the in-phase state. As will be demonstrated later, the in-phase state also occurs where a seismic wavelet exists. Therefore, this state can be used as an indicator to locate a seismic reflector. The above decomposition is now applied to a signal with two spikes. Denote the two spikes as the function S(t;to,ti) = S(t;to) + S(t;ti). The FT can be written as S(t; to, h) = ^(e^C*-*") + - e - M ) . e ( . ) 4 16 In Figure 4.4, the signal is decomposed into a combination of sinusoids with different amplitudes and phase delays. Note that in this case, the in-phase state clearly visible in one spike example does not exist. Therefore it is difficult to obtain any information from these diagrams directly. The Morlet WT can be written as Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT sit-M) ^ ( ^ ^ V ^ M = Gxp a + -(^) E 2 E ^(^)). 71 ( 4 . 1 7 ) a, Figure 4.5 shows the decomposition and the attributes using the Morlet WT. The amplitudes show energy distribution of the complex traces in time and frequency. In contrast to the Fourier phases in Figure 4.4, the wavelet phases maintain the pattern exhibited by the single spike, two in-phase patterns indicating the locations of two spikes. Next, the Morlet wavelet transform is applied to some more realistic examples. 4.3.2 Wavelet attributes of traces The Ricker wavelet (Sheriff and Geldart, 1982) shares a very basic feature with all real seismic wavelets, it is band-limited. It has been widely used as synthetic wavelet for seismic modeling. Here, the Ricker wavelet is used to examine and understand the characteristics of WT attributes and trace decomposition. Figure 4.6 shows a Ricker wavelet, its wavelet attributes, multiscale traces, reconstructed signal using the multiscale traces and the error of the reconstruction. As expected, the amplitudes indicate the energy distribution of the Ricker wavelet in time and frequency, while the in-phase state of phases indicates the location of the signal independent of the energy in frequency direction. The real part of the complex multiscale trace shows the waveform of the decomposed Ricker wavelet. Figure 4.7 is equivalent to Figure 4.6 except that two Ricker wavelets with different amplitudes are considered. The amplitude difference can be clearly seen in the wavelet amplitudes. This difference in amplitudes has no effect on the wavelet phases and this implies that identifying the in-phase state may be a good approach to the identification of the location of seismic events with low energy. Figure 4.8 is identical to Figure 4.6 but with the addition of noise with standard Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 72 deviation of 5% of the maximum amplitude of the Ricker wavelet. The amplitudes still show the time and frequency energy distribution. Despite the fact that the phases are very different from those shown in Figure 4.6, the in-phase state at the 50th trace sample is still able to indicate the location of the Ricker wavelet, although the existence of noise has made the in-phase state much narrower in the number of time samples it spans. Figure 4.9 demonstrates results for a portion of a real trace. The envelopes of wavelet coefficients illustrate the variations of reflection strength at different frequency bands. In the diagram of WT amplitudes we can see that the energy of the surface wave starts at approximately the 40th sample and is located at low frequency. In the phase spectrum, some in-phase states can also be recognized at a few locations indicating the locations of a few strong events. 4.4 Wavelet attribute analysis for 2-D data Applying the above approach to a 2-D data set, in a trace by trace manner, we will obtain a 3-D data set in time, frequency and a new variable - offset. These 3-D attributes not only show the location of seismic events but also how they vary with frequency. The next example demonstrates the identification of the reflections, ground roll and head waves in a shot record according to their wavelet attributes in time and location, as well as the frequency bands which they occupy. The data are displayed in Figure 4.10 with attributes as defined by Taner et al. (1979). Since these attributes are the contribution of all frequency components, they cannot reflect variations due to different wave modes in the data. For example, at certain locations where both surface waves and reflections exist, the Taner et al. attributes provide averages for both types of waves collectively for all frequencies. Wavelet attributes, on the other hand, tend to differentiate the wave types into location and frequency. Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 73 Figures 4.11 and 4.12 show the wavelet amplitude and phase attributes at frequencies of 10, 20, 30, 40, 50 and 60Hz. The surface waves are contained mainly in the first two panels, the reflections begin at the third panel, the head waves exist in all panels. The phases show a large amount of events because of phase wrapping. These events are called phase events to differentiate them from the actual seismic events. While it is obvious that the number of phase events does not correspond to the number of seismic events, the existence of phase events does suggest the existence of seismic events. Based on tests in the 1-D case, the actual location of a seismic event should occur where the phase events appear consistently for all frequencies. Here the identification of a particular mode in a certain frequency band, rather than the exact location of that event in time and offset are considered. It is clear that the phases of the reflections appear well defined in the second panel where the seismic reflections are very weak in amplitude and, in fact, can hardly be identified. This leads to the conclusion that the multiscale wavelet amplitudes and phases together provide a powerful tool to study the characteristics of the seismic signal in time and offset, as well as in frequency. 4.5 Summary Morlet WT attributes were used to analyze the variations of the characteristics of seismic data in time and offset, as well as in frequency (scale). Some simple examples were used to understand the meaning of the WT amplitude and phase: the amplitude indicates the energy of the signal in location and frequency and the in-phase state, or a phase consistency over all scales, indicates the location of the signal. Multiscale attribute analysis provides a useful tool for analyzing the frequency bands of both signal and coherent noise such as ground roll and head waves. While the WT amplitudes show the energy distribution over time, offset and frequency, the WT phase Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 74 events show the existence of seismic events independent of their energy. The amplitude and phase attributes together provide a powerful tool for the identification of reflection events at low frequencies where the energy is dominated by ground roll. The time- and frequency-variant characteristics obtained using the multiscale attribute analysis approach can be applied to the design of wavelet transform or Gabor filters. Some of the applications in 1-D wavelet filtering can be found in Miao and Cheadle (1995) and in 2-D Gabor filtering in Li and Ulrych (1996). Multiscale attributes are also expected to have important application in the analysis of features which exhibit time-frequency variations in the post stack domain. Multiscale attributes may also contribute to reservoir delineation, particularly in respect to the fourth dimensional measures which are increasingly being used in reservoir studies. Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 75 Figure 4.1: Fourier decomposition of a spike. Top row: the real and imaginary parts of the ensemble of complex harmonics. Bottom row: the amplitudes and phases of the complex harmonics. The dark amplitude indicates unity at all times and frequencies. For all diagrams, the horizontal axis is in time/location and vertical axis is in frequency/scale. Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 76 Figure 4.2: Morlet WT decomposition of a spike. Top row: the real and imaginary parts of multiscale signals. Bottom row: the amplitudes and phases of the multiscale signals. WT phase FT phase 40 20 g> < 0 f=20Hz -20 0.1 Time (s) 0.2 -40 0 M=120Hz 0.1 Time (s) 0.2 Figure 4.3: Unwrapped FT and WT phases the complex signals in Figures 4.1 and 4.2. Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT \AAAAAAAAAAAAAAAAAAAAAAAAAAAAMAAAAAMAAAAAAAAAAAAAi AAAAAAAAAAAAAAAAAAAAAAMAAAAAAAMAAAAAAAMAAMA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAA 77 mmmmmmmmwmmmmmmmj \AAAMAAMAAAAAAAAAAAAAMAAAAAAAAAAAAAAAAAAAAAA LAAAAAAAAAAAAAAAAAAAAAAAAAA^ \AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA KAAAAAAAAAAAAAAAAAAAAA^J K A A A A A A A A A A A A A A ^ Figure 4.4: FT decomposition of two spikes. Top row: the real and imaginary parts of the ensemble of complex harmonics. Bottom row: the amplitudes and phases of the complex harmonics. Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 78 Figure 4.5: Morlet W T decomposition of two spikes. Top row: the real and imaginary parts of multiscale signals. Bottom row: the amplitudes and phases of the multiscale signals. Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 79 Figure 4.6: Morlet WT attributes and multiscale traces of a Ricker wavelet. Left from top: Ricker wavelet; amplitude (envelope) and phase of wavelet coefficients; Right from top: decomposition of Ricker wavelet at 16 scales; reconstructed signal from 64 scales; the difference between the reconstructed signal and the original signal of 128 samples. Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS 20 40 60 80 100 120 Time (sample) 20 40 60 80 100 120 Time (sample) WT 80 250 Figure 4.7: Morlet WT attributes and multiscale traces of two Ricker wavelets. See Figure 4.6 for explanations. Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 81 Figure 4.8: Morlet WT attributes and multiscale traces of Ricker wavelets with noise. See Figure 4.6 for explanations. Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 82 2000 20 2501 • 20 40 60 80 100 120 Time (sample) • • • • 40 60 80 100 120 Time (sample) Figure 4.9: Wavelet transform attributes and multiscale traces for a portion of a real trace. See Figure 4.6 for explanations. Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 83 Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 84 Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 85 Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 86 Chapter 4. ATTRIBUTE ANALYSIS USING THE CONTINUOUS WT 87 Chapter 5 ANALYSIS A N D SYNTHESIS OF W E L L LOGS 5.1 Introduction Seismic inverse problems are usually underdetermined and some prior knowledge about the geology is required to constrain the inversions (Tarantola, 1987). Since well logs are the most accurate measurement of geophysical properties such as velocity and density, modeling of well logs is an important approach to obtain the prior information. It is well accepted that these geophysical properties shown some self-similarity which can be modeled as a 1/f process (for example, Turcotte, 1992). This self-similarity is the results of erosion, sea floor spreading, sea level changing and other geological phenomenon. Modeling of the well logs relates the above regional geological information with the seismic inverse problems. A 1/f-type process is defined as a random series showing self-similarity. A 1/f process can be statistically stationary or non-stationary. Stationary processes can be modeled using the Fourier power spectrum whereas non-stationary processes can not (Flandrin, 1989). In fact, recent studies show that for a non-stationary 1/f process, localized transforms, such as the wavelet transform, should be used (Flandrin 1992; Wornell, 1990; Wornell and Oppenheim, 1992; Stoksik, 1994). Experimental studies suggest that velocity and density logs can be modeled as non-stationary fractional Brownian motions (FBMs), in which case, the Fourier power spectrum is not suitable to model the process. A correct approach, but an indirect one, as pointed out by Mandelbrot and Van 88 Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 89 Ness (1968), is to infer a F B M from its derivative which should be a stationary process. The wavelet transform provides a direct way to model the FBMs, and therefore, unlike the Fourier transform, can be used as a tool to analyze and simulate the time-variant characteristics of the process. The goal of this chapter is to introduce and discuss some recent research results on modeling and simulation of FBMs and then to apply the results to the study of well logs. Some of the studies related to 1/f processes are reviewed in the following two sections. Section 4 discusses how the orthogonal wavelet transform can be used to estimate the model parameters of a 1/f process for field logs. Synthetic logs are generated in section 5 using the wavelet approach in both 1-D and 2-D cases for given fractal parameters. 5.2 1/f-type processes A stationary 1/f process can be modeled by its Fourier power spectrum P(u) according to P{w) oc l/uj . a (5.1) The scaling parameter a determines the magnitudes of the long and short term correlations. The above definition can only be applied to 0 < a < 1 when the process is stationary, in which case, a time-independent covariance function can be computed from the power spectrum (Pilkington and Todoeschuck, 1992). When a > 1, the process is said to be non-stationary. Specifically, when 1 < a < 3, the process is called fractional Brownian motion (FBM), the covariance function is time dependent and the process has to be defined using a covariance function of the form (for example, Flandrin, 1992) Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS E{b (t)b (r)} H H = ^(\t\ 2H 90 + | r | " - |t - r | * ) , 2 where &#(£) is the F B M process and the constant H = 2 (5.2) is commonly adopted by the signal processing community. Since the meaning of the Fourier power spectrum for non-stationary processes is not clear, the Fourier spectrum cannot be used to describe their characteristics. There are two conventional ways to handle the problem of nonstationarity. One way to overcome the problem is to re-define the power spectrum in equation (5.1) as an average spectrum (Flandrin, 1989) but this average spectrum may be affected by the existence of a trend. Another approach is to define a F B M from its stationary derivative which is a stationary process even when the original process is not. The relationship between the parameter a in a F B M process and the parameter 8 for the derivative is 8 = 2 — a. Assuming P,(u>) is the power spectrum of any signal s(t), either stationary or non-stationary, and Pd(w) is the power spectrum of d{t) = ds(t)/dt, it is clear that P (u) = u) Pd(oj). If the process is assumed to be F B M , in order to 2 s estimated the fractal parameter a, the parameter 8 has to be estimated by fitting 8 to the power spectrum of the derivative, and then to compute a from 8. The drawback of this approach is that we have to assume the process is F B M before analysis. Recent work suggests that the wavelet transform, with localized basis, is a natural choice for the direct analysis of F B M . Before we consider the wavelet transform approach, some summarizing remarks concerning 1/f processes are in order. There are three special cases in denning a random process. The case a = 0.0 results in a white noise process, while a = 1.0 and a = 2.0 give rise to flicker and Brownian noise respectively. A flicker noise process is sometimes called a Joseph process (Todoeschuck and Jensen, 1988). For 0 < a < 1, the process is called Gaussian fractal noise and for 1 < a < 3, as we mentioned earlier, it is F B M . All the processes can be fully characterized Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS by three parameters: the mean 91 the variance o~ and the fractal parameter a. The 2 mean and variance control the amplitude of the process, while the shape of the process is affected by the parameter a. The definitions are summarized in Table 5.1 for 0 < a < 2. White GFN Joseph FBM OBM x(t) a H 0 -\ 0->l - § - > 0 1 0 1 -> 2 0 -» | 2 \ OBM FBM Joseph GFN White Derivative of x(t) P = 2- a H' 2 \ 1 -> 2 0 -> f 1 0 0 l -\ -> 0 0 -| Table 5.1: The random processes with 0 < a < 2. where a = 2iJ + 1 and (3 = 2i7' + 1. GFN: Gaussian fractal noise; OBM: ordinary Brownian motion; F B M : fractional Brownian motion. 5.3 5.3.1 Diagonalization of the covariance matrix Covariance matrix in a transformed space It is well known that the model parameters of a well log are correlated, i.e., any parameter is related to its neighboring parameters. A commonly used approach to analyze this correlation is to transform the model parameters to a new space and hope that, in this space, the relationship between a transformed coefficient and its neighboring ones is simplified. This simplification can be described as the diagonalization of the covariance matrix. For a stationary process, the covariance matrix is Toeplitz and the covariance matrix for the Fourier coefficients is a diagonal matrix. Since, for a F B M process, the relationship between a parameter and its neighboring ones varies with the location of the parameter, the Fourier transform cannot diagonalize the covariance matrix. This section will consider the use of the wavelet transform in the diagonalization of the F B M Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 92 covariance matrix. Given a column vector m which represents a random process, the estimated covariance matrix is C = E{mm } T where E is the expectation operator. If \£ represents an orthogonal transform with in = ^m, since m = S£ m, the covariance matrix for the T transformed model parameters is C =E{mm } T = *£{mm }* r = *C* r (5-3) r This equation shows the relationship between the covariance matrices in the original and the transformed model spaces. So far, in geophysical data analysis, two transforms have been applied to the above equations according to the type of the process involved. First, for a stationary process in the discrete case, C is Toeplitz: each row has the same covariance function but with a cyclic shift. Application of the Fourier transform operator \I> = F T will diagonalize the covariance matrix, i.e., in this case C is diagonal with diagonal elements given by the Fourier transform of the covariance function. Second, for a more general process, the Karhunen-Loeve (K-L) transform can be applied to the above equation. The K-L transform for the covariance matrix C in this case can be defined by choosing \P = U , the eigenvector matrix computed from the spectral decomposition of T C = TJ^AU where A is the eigenvalue matrix. In this case, the covariance matrix C for the transformed coefficients is simply the diagonal matrix A. The drawback to the first approach is that it is only suitable for stationary processes while the drawback to the second approach is that the K-L transform may be time consuming. Recent research suggests that a K-L like transform, such as the wavelet transform, can be used to approximately diagonalize the covariance matrix (Wornell, 1990). The covariance matrix of a F B M in the WT space will now be examined and its Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 93 characteristics will be used to analyze the features of well logs. 5.3.2 Diagonalization of a F B M covariance matrix using the W T The covariance matrix in the wavelet transformed domain can be written as C = WCW , (5.4) r where W denotes the wavelet transform in matrix form. Assuming m is the vector of a F B M process, the WT of m is in = Wm with m^) the coefficient vector at scale j. There are two important characteristics about C (Flandrin, 1989, 1992; Wornell, 1990; Wornell and Oppenheim, 1992). is a random series and its variance follows the power law: vor(in -) = a 2~ , 2 ja J (5.5) with a representing an average standard deviation for the coefficients at all scales. • The off diagonal elements of C have very small values and decay very rapidly. The decay rate depends on the number of vanishing moments of the wavelet, where the decay is rapid for wavelets with a large number of vanishing moments. Figure 5.1 shows the covariance matrices C and C computed using Haar and Daubechies 4 and 10 point wavelets. The Daubechies 10 wavelet with the largest number of vanishing moment among the three wavelets results in the best diagonalization of the covariance matrix. Although a form for an ordinary F B M process using the Haar wavelet is shown by Flandrin (1992), in general, the analytical form of C is difficult to obtain. In this example, the 128 x 128 covariance matrix is first computed, then C is obtained by applying the wavelet transform to both the rows and the columns of C. Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 94 The 1/f behavior has been studied by various authors (e.g., Todoeschuck and Jensen, 1988; Rosa and Ulrych, 1991) and this behavior implies a non-white nature of the earth reflectivity series. The reflectivity series are usually stationary fractals and can be represented by a covariance function or its power spectrum. In some geophysical problems, however, such as traveltime inversion, it is the impedance or velocity which is modeled, quantities which are related to the integral of the reflectivity series and represent nonstationary F B M processes. Conventional Fourier spectrum analysis is, consequently, not applicable. The study of F B M processes can be simplified significantly due to the characteristics of a F B M process in the wavelet domain. A non-stationary process can be represented using a small number of smoothed coefficients to represent the trend and the rest of the detailed coefficients at different scales can be approximately considered as stationary processes. The variances of the wavelet coefficients at different scales can be described using only two parameters W and a. The averages of the detailed coefficients at different scales are considered to be zero because the detailed coefficients can be treated approximately as the discrete derivatives of the data at different scales. The next section will use the features described above to characterize real well logs in the wavelet domain. 5.4 5.4.1 Fractal parameter estimations of logs The wavelet transform approach Suppose a signal f has N = 2 samples. Its wavelet transform, f, consists of a vector J of smoothed coefficients s^°) at the coarsest scale jo, and vectors of detailed coefficients d^) from scale jo to the finest scale J — 1: f = { U) s d( '°) d( ' ) J J 0+1 ••• S ~% J (5.6) Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 95 The vector s^^ contains the trend and will not be studied here. The vectors d ^ can be 0 considered as random series with variances obeying equation (5.5). The parameter a can be calculated in two steps. • Compute the estimated variance of the detailed coefficient at scale j using where Nj = 2 is the number of coefficients at scale j or, in other words, the length J of the array d^. • Compute a by fitting log (o-j) = —j a/2 + constant (i.e., from equation (5.5)). 2 5.4.2 Examples Figure 5.2 shows an impedance log, named 'Logl' for convenience, which contains a total of 3207 samples. The log has been resampled to an interval of 0.3 msec, from 0.82 to 1.78 seconds. Since the characteristics of the first 2/3 portion of the log are very different from the rest, the log is subdivided into three segments as follows: segment 1 includes samples 1 — 2048; segment 2 includes samples 2100 — 3123; and segment 3 includes samples 1100 — 3147 which includes segment 2 and portion of segment 1. The wavelet transforms of Logl are shown in Figure 5.3. Shown on the left are the wavelet coefficients computed using the Haar and Daubechies wavelets for segment 3. Shown on the right are the multi-resolution signals which can be obtained by transforming the wavelet coefficients at each scale back to the time domain. The smoothed signal s^) 4 is shown on the top of the figures on the right and segment 3 is shown at the bottom. The coefficients clearly show the difference between the first and second portions of the segment and accentuate the advantage of using the wavelet transform in allowing the Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 96 localized features to be identified. In contrast, the Fourier transform will not show any of the location information. Figure 5.4 shows the estimated variances of the wavelet coefficients at different scales using the Haar, Daubechies 4 and 10 point, and symmlet 10 point wavelets which are shown in solid, dotted, dashed, and dot-dashed lines respectively. The fitted a's are about 1.3 for segment 1 and about 1.7 for segment 2. The a ~ 1.46 is for segment 3 which partly covers segments 1 and 2. The detailed results are listed in Table 5.2. (samples) Haar Daub4 DaublO SymmletlO Segment 1 (1-2048) 1.18 ± 0.14 1.31 ± 0.07 1.29 ± 0.09 1.33 ± 0.14 Segment 2 (2100-3123) 1.53 ± 0.11 1.80 ± 0.20 1.73 ± 0.19 1.77 ± 0.22 Segment 3 (1000-3047) 1.35 ± 0.12 1.45 ± 0.15 1.55 ± 0.16 1.47 ± 0.22 Table 5.2: The fractal parameters estimated from Logl. The numbers inside the parentheses indicate the ranges of the segments in terms of samples. The second example consists of density and sonic logs from three different wells in the same region. The logs are shown in Figure 5.5(a). The depths of logs start at about 300m to about 1900m with a sampling interval of 0.3048m. The Daubechies 4 and 10 point wavelets are used to compute two estimated a's. The characteristics of the logs vary greatly with depth, and for this reason, the logs are subdivided into three segments with 1024 samples per segment. These segments are shown in Figure 5.5(b). The values of fractal parameter a for most of the segments are between 1 and 2. This suggests that the logs in general show F B M characteristics. The exceptions are the fractal parameters of the second segment for Wells 3 and 4 which are close to, or less than, one and imply that an ordinary Gaussian fractal process applies at that depth interval to both wells. 97 Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS Segment 4, from depth 320 — 1568m with 4096 samples, overlaps all the first three segments. It is obvious that the fractal nature of segment 4 cannot be characterized by a single parameter. Just for the sake of curiosity, fractal parameters for these segments at all wells were computed and, as expected, the estimated a's for the fourth segments are not the averages of the a's of the first three segments. This further confirms that the log should be analyzed by segmentation. Other than the simple segmentation used in this section, more sophisticated methods may be explored in future studies. For example, the variation of the parameter a with depth may be computed over a sliding window. (dO/N) Son 2 Den 2 Son 3 Den 3 Son 4 Den 4 Segment 1 (320/1024) 1.82 ±0.31 1.83 ±0.30 1.10 ±0.21 1.40 ±0.21 1.48 ±0.16 1.50 ±0.14 Segment 2 (800/1024) 1.47 ±0.23 1.85 ±0.11 0.64 ±0.22 0.97 ±0.12 1.17 ±0.13 0.86 ±0.11 Segment 3 (1200/1024) 1.71 ±0.32 1.82 ±0.27 1.34 ±0.22 1.57 ±0.26 2.00 ±0.06 1.77 ±0.16 Segment 4 (320/4096) 1.14 ±0.24 1.28 ±0.21 0.88 ±0.13 1.06 ±0.15 1.08 ±0.16 1.08 ±0.16 Table 5.3: Estimated fractal parameters from Wells 2, 3 and 4. Where dO is the start depth, N is the number of samples. 'Den' and 'Son' represent density and sonic logs respectively. Segments 1, 2 and 3 are shown in Figure 5.5(b). The interpretation of the estimated fractal parameters is an interesting topic but because of the lack of geological information to support any interpretation, this chapter will concentrate on the modeling of the well logs. To confirm the model has fractal properties, some synthetic models will be tested in the next section. Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 5.5 98 Synthesis of fractal processes 5.5.1 1-D Synthesis, as the opposite of analysis, uses the fractal parameter a measured in the analysis to simulate a F B M like process. The synthesis of a 1-D fractal can be done using the following steps. • Generate a series of random vectors with variances following equation (5.5): f = a{2~ s \ joa (jo 2~ n \ joa {ja 2~ ' (j 0+1)a n '' (: 0+1) , ••• 2- (J_1)a n ~ }, (J 1) (5.8) where vector n ^ denotes a random vector of 2 elements and s^°\ is the vector of J the smoothed coefficients which can be chosen randomly or assumed in accordance with prior knowledge. • Compute the inverse WT of f to obtain the fractal time series. Figure 5.6 shows an example of using the above procedure to synthesize some fractals. Figure 5.6(a) shows the coefficients f generated from a Gaussian random series. Figure 5.6(b) and (c) show the synthesized fractals and their multiscale signals. The fractal parameter a can also be considered to vary with depth to construct a fractal with depth variant characteristics. In this example, a = 1.2 is chosen for half of the coefficients on the left, and a = 1.6 for half of the coefficients on the right. As we can see, the right half of the synthetic signals with large a exhibit less small scale structures than the left half. Note that the three fractals generated using different wavelets show some differences at the edges. This edge effect could be reduced by putting some zero coefficients at both edges. Alternatively, wavelets with less edge effects can be used (Cohen, et al., 1993), but this aspect has not been considered here. Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 99 In general, synthetic logs show some similarity to real logs, with the difference being that the real log is more spiky. This spiky feature cannot be modeled as a F B M process and this issue is discussed in section 5.6. 5.5.2 2-D The simulation of 2-D velocity fields is an important topic in some geophysical applications such as deconvolution (Jensen, et al., 1995) and geostatistics (see, for example, Meng, 1995). The methods in the previous subsection for 1-D fractal synthesis can be extended to the 2-D case. As shown in Chapter 2, the WT decomposes a two dimensional image X into a four dimensional image X. The 4-D image consists of two levels. The first level consists of a pyramid structure X {S( ); J0 = ... H^°), D^°\ V( °); H^ \ J 0+1 H^- ), D( - ), V^" )}, 1 J X D^ 0 + 1 \ V( ' ); J 0+1 (5.9) 1 where matrix S( °) contains the smoothed coefficients at the coarsest scale jo and coeffiJ cient matrices Y)^ and \ ^ contain structures in horizontal, diagonal and vertical directions at scale j. The second level is made up of these 2-D coefficient matrices. Each of these matrices can be used to construct an image at a certain scale and with certain directional characteristics. If X is an uncorrelated random image with the structure of equation (5.9) and with a standard deviation W, an image with F B M character can be constructed by scaling the coefficients according to the following power laws: Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS crf )2 af where a^ \ h 2 a (5.10) 2 tr ? { =W 2~i », =a 2-i°"', 2 )2 100 =W 2-i«\ 2 and <T^ are the standard deviations of matrices H , D and V respectively and ah, «d and a are the fractal parameters. In the previous section, we have shown v that for a F B M process the wavelet coefficients can be approximately characterized as uncorrelated random series at different scales with variances following a power law. Equation (5.10) can be considered as the extension of equation (5.5) to 2-D in which case the wavelet coefficients of the F B M image consist of many uncorrelated random images with their variances following power laws for the matrices H's, D's and V s . The above method can be processed in three steps: • Generate a random image X consisting of the wavelet coefficient matrices D( ) and V ^ ) with standard deviation a. J • Apply scaling factors 2~i , 2~^ ah ad and 2~^ '" to the matrices H's, D's and Vs a respectively. The new matrix X follows the power laws in equation (5.10). • Apply the inverse 2-D WT to X . Figure 5.7 shows an example of the synthesis of 2-D fractals. Figure 5.7(a) shows a random Gaussian image which is then scaled according to equation (5.10) and is shown in (b). Figure 5.7(c) and (d) are the inverse wavelet transforms of the image in 5.7(b). Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 5.6 101 Non-fractal component 5.6.1 Definition of the non-fractal component As was discussed before, the statistical features of a fractal can be described using the mean, variance and fractal parameter a. However, some of the features cannot be described statistically by a fractal model. These features include the spiky characteristics and the low frequency, large scale or long time variation representing a non-stationary trend. These features are common in well logs and for these reasons a log should be represented as a combination of a fractal component and a non-fractal component. The non-fractal component can be defined in the wavelet transformed domain as follows. Assume the wavelet transform of a signal f is f as expressed by equation (5.6). The first portion of the non-fractal component in the wavelet domain is defined as the smoothed coefficients s^°) representing the trend. The second portion of the non-fractal component consists of the spiky coefficients that do not follow the power law in equation (5.5). Assuming that <Tj is the standard deviation of the wavelet coefficients at scale j, this portion of the coefficients are those a times larger than the standard deviation at that scale, i.e., \^ \>a<r . ) i (5.11) Clearly, if a is large enough, the non-fractal coefficients consist of only the first part, the trend. In this case, the fractal portion of the signal can be obtained by the inverse wavelet transform of all the detailed coefficients. 5.6.2 A simple way to evaluate the fractal model Two important parameters, the coarsest scale jo and the threshold a, have to be chosen in well log analysis. Scale jo determines the number of scales used to represent the trend Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 102 and the threshold controls the number of outliers that do not statistically fit well into a fractal model. Choosing these two parameters is an arbitrary task and the following example is used to understand the effect of these two parameters. Figure 5.8 shows an example of the separation of the non-fractal signal and the fractal signal using segment 3 of Logl (shown as the dotted lines). Overlapped on the dotted lines are the non-fractal structures calculated with a = 2 and 4. The solid line on the bottom is the difference between the log and the non-fractal data. In this example, the coarsest scale is jo = 4, so there are 16 smoothed coefficients. There are 105 and 19 outliers for a = 2 and 4 respectively. Table 5.4 lists the number of outliers in the cases of a = 1 to 4 for the three segments of Logl. Clearly, compared to the total number of coefficients, 2048, most of the coefficients fit a fractal model. From the viewpoint of information extraction, segment 3 of 2048 samples can be characterized in the wavelet domain using 16 smoothed coefficients representing the trend, 19 outliers (with threshold o = 4) representing significant local structures, and two constants, an average standard deviation a and a fractal parameter a, representing the fractal component. a 1 2 3 4 Segment 1 (1-2048) 261 85 36 15 Segment 2 (2100-3123) 180 63 18 4 Segment 3 (1000-3047) 286 105 47 19 Table 5.4: The numbers of non-fractal outliers with different thresholds. The threshold is defined as constant a multiplying the estimated standard deviation at each scale. Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 5.7 103 Summary The wavelet transform is used as a K-L like transform to approximately diagonalize the covariance matrix of a fractional Brownian motion. The diagonalized covariance matrix relates the wavelet coefficients of a fractal process in two ways: • the detailed coefficients at each scale are Gaussian processes; • the variances of wavelet coefficients at different scales follow a power law. Both properties have been demonstrated by Wornell (1990) and Flandrin (1992) and have been studied by others. The diagonalization is dependent on the particular wavelet basis used and a theoretical formulation of the diagonalized matrices is very complicated. In this chapter a simple numerical approach using matrix operations to diagonalize the covariance matrix of F B M was proposed. Examples using different wavelet bases to compute the diagonalized covariance matrices were shown. The analysis of the fractal characteristics was demonstrated using a few real logs from two different regions. Because of the obvious variation with depth, the logs were subdivided into a few segments. The fractal parameters of these segments were estimated separately by fitting the variances of the detailed coefficients at different scales to power laws. The estimated fractal parameters are in general between 1 and 2 and indicate the F B M nature of the well logs. The synthesis of a well log can be simply accomplished by constructing a series of random vectors with variances following a power law and then by applying an inverse wavelet transform. Since the trend and the spiky features cannot be modeled as fractal, the synthesized fractals are less spiky and with less trend. The fractal analysis and synthesis are based on the assumption that most of the wavelet coefficients at each scale belong to a stochastic process, for example, Gaussian Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 104 in nature. Two portions of the wavelet coefficients cannot be characterized as stochastic processes: the smoothed coefficients and the outliers that are a few times larger than the standard deviation at each scale. These wavelet coefficients comprise the non-fractal component of a well log. Except for the smoothed coefficients and the outliers, the rest of the coefficients can be characterized by a fractal parameter and an average standard deviation. The characterization of well logs is very important in incorporating the correct constraints in inverse problems and these aspects will be studied in the next two chapters. 105 Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS (a) Covariance 20 40 60 80 (b) Haar 100 120 20 (c) Daub4 20 40 60 80 40 60 80 100 120 (d) DaublO 100 120 20 40 60 80 100 120 Figure 5.1: The covariance matrix of a fractional Brownian motion (FBM) and its wavelet transforms, (a) The covariance matrix and its wavelet transforms using (b) Haar, (c) Daubechies 4 point and (d) Daubechies 10 point wavelets. 106 Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS Log1 x 10 500 1000 1500 2000 Sample number 2500 3000 Figure 5.2: Log 1. The horizontal direction is annotated in both time and number of samples. The fractal analysis is done for three segments with the samples from 1 — 2048, 2100 - 3123 and 1000 - 3047 respectively. Segments 1 and 2 have two different types of characters and segment 3 overlaps both segments 1 and 2. Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS Haar: Multi-scale decomp Haar: Multi-scale coeff -4 -5 -6 -7 -8 -9 -10 l. i I • 11 i . 1 . i 'I." \ . -4 -5 -6 -7 -8 -9 -10 •'4' "M'tHM I L . ^ i M H l»H« Daub4: Multi-scale decomp Daub4: Multi-scale coeff -4 -5 -6 -7 -8 -9 -10 T T 107 -4 -5 -6 -7 -8 -9 -10 As* .... •« ' 1 '»• * 1 1 1 .1• ... .. ^J^^~^Afi^^AJ\^^ >. . . J,. . Figure 5.3: The wavelet coefficients and multiscale decompositions of segment 3. Top row: by the Haar wavelet; bottom row: by the Daubechies 4 point wavelet. The multiscale coefficients show clear differences between the left halves and the right halves of the log for both wavelets. Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS 108 Variances Variances 16 14 E o 10 8 4 6 8 10 Scale X 1 Q4 Impedance 1 1.2 Time (sec) x 10 1.4 4 1.5 Impedance 1.6 Time (sec) 1.7 Figure 5.4: The multiscale variances for segments 1 and 2 of Log 1. The estimated variances using the Haar, Daubechies 4, 10, and Symmlet 10 point wavelets are shown in solid, dotted, dashed, and dot-dashed lines respectively. The fractal parameters are fitted using the variances from the sixth scale to the smallest scale. Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS Son2 400 200 200 I J l_ 600 800 1000 1200 1400 1600 1800 2000 400 600 800 1000 1200 1400 1600 1800 2000 L Den2 0 l 200 Son3 400 200 _l 0 2500 2000 1500 I 400 _1 0 2500 2000 1500 109 200 : 4 I I I I L. 400 600 800 1000 1200 1400 1600 1800 2000 400 600 800 1000 1200 1400 1600 1800 2000 400 600 800 1000 1200 1400 1600 1800 2000 400 600 800 1000 1200 Depth (m) 1400 1600 1800 2000 Den3 0 200 Son4 400 200 200 2500 2000 1500 Den4 0 200 Figure 5.5(a): The sonic and density logs of Well 2, 3 and 4. The depths of the logs are in the range of approximately 300m to 1900m and the sampling rate is 0.3048m. Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS Son2 400 200 0 200 400 600 800 1000 1200 1400 1600 1800 2 5 0 0 - Den2 2000 1500 0 200 400 600 800 1000 1200 1400 1600 1800 400 600 800 1000 1200 1400 1600 1800 200 400 600 800 1000 1200 1400 1600 1800 200 400 600 800 1000 1200 1400 1600 1800 400 600 800 1000 Depth (m) 1200 1400 1600 1800 Son3 400 200 0 2500 2000 1500 200 Den3 0 Son4 400 200 h 0 2500 2000 1500 Den4 0 200 Figure 5.5(b): Three segments of the sonic and density logs of Well 2, 3 and 4. estimated fractal parameter a's are listed in Table 5.3. Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS (b) Haar (a) Synthesized coefficients J -3 -4 -5 -6 -7 -8 -9 -10 L J l l L i .'I'.'MI . . I ' I 111 -3 -4 -5 -6 -7 -8 -9 -10 0.5 0.5 (c) Daub4 (d) DaublO -3 -4 -5 -6 -7 -8 -9 10 0.5 0.5 Figure 5.6: Synthesis of 1-D F B M processes, (a) A Gaussian random series scaled according to the power law and and shown in multiscale manner; (b)-(d) Multiscale signals of the coefficients in (a) using three different wavelets. The synthesized fractals are shown at the bottoms of the figures. Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS (c) Daub4 20 40 60 80 112 (d) DaublO 100 120 20 40 60 80 100 120 Figure 5.7: Synthesis of 2-D F B M images from a random Gaussian image, (a) A random image; (b) the scaled image; (c) and (d) synthesized images using the Daubechies 4 and 10 point wavelets. Chapter 5. ANALYSIS AND SYNTHESIS OF WELL LOGS x 10' (a) 2*sigma 1.4 Time (sec) x 10 113 1.5 (b) 4*sigma 1.4 Time (sec) Figure 5.8: Decomposition of fractal and non-fractal structures from Log 1. The non-fractal wavelet coefficients include the smoothed coefficients and the spiky coefficients which are several times, denoted as a, larger than the standard deviation of the wavelet coefficients at each scale. In this case, (a) a = 2; (b) a = 4. The fractal components which are the residuals after the removals of the non-fractal components are shown at the bottoms of the two diagrams. The dotted lines represent the original log. Chapter 6 INVERSION WITH PRIOR INFORMATION 6.1 Introduction Since seismic traveltime tomography is in general an ill-posed and underdetermined problem, some a priori information about the model is required to yield a solution (see, for instance, Jackson, 1972; McMechan, 1982; Oldenburg, 1984; and Lines and Levin, eds., 1988). The inverse problem is to find a model that honours the data to some specified degree, and subject to some prior model information. This goal is achieved by minimizing an objective function that consists of a data misfit term and the model regularization term. The data misfit controls how well the model should fit the data. The model regularization term, or the generalized model norm, can be constructed by a stochastic approach - the Bayesian approach via a model covariance matrix (Menke, 1980; Tarantola, 1987), or by a regularization approach via a combined norm of the model parameters, and/or the first or second derivatives of the model parameters (Scales and Smith, 1994). This chapter discusses a few approaches of incorporating the prior information into the inversion. The limitation of these conventional approaches and the need of a more general approach of incorporating non-stationary prior information will be discussed and a new approach of using the wavelet transform to incorporate the non-stationary constraints will be introduced in the next chapter. This chapter is organized as follows. First, traveltime tomography and the general inverse problem are introduced. Then, the relationship between the regularization and 114 115 Chapter 6. INVERSION WITH PRIOR INFORMATION stochastic approaches is explored by approximating a weighting matrix using a combination of an identity matrix and discrete derivative matrices. With this approximation, a conjugate gradient approach can be used to solve the inversion of a large system of equations. 6.2 Traveltime tomography The seismic traveltime can be related to a slowness field s by / sdl = t(s), Jr(s) (6.1) where T(s) is the raypath connecting source and receiver. Since T is a function of s, the problem is non-linear. In the discrete case, the above equation can be written as Am = t, (6.2) where t is a vector of N observed traveltimes, m is the discrete slowness field with M elements and A is the forward matrix that relates the model parameters to the observations. Let ft m = R A maps an element in fl m M be the model space and fid = R N be the data space. to fi^. The matrix A, in general, depends on m and the problem is non-linear. Some assumptions about the model are usually made to solve the problem (Menke, 1984; Norlet, ed., 1987; Parker, 1996). The non-linear problem is solved by linearization and iteration. For a given initial model W , the matrix A<°> and t<°) are calculated and the model can be improved by S solving (°)(s( )-^ )) = t 1 A 0 o b 3 -t(°), (6.3) Chapter 6. INVERSION WITH PRIOR INFORMATION where t° * is the vector of observed traveltimes. A solution 6 116 is found by minimizing an objective function and is then used as an initial model for the computation of A^ ) for the 1 next iteration. This iterative process is completed when the minimum of the objective function is reached. The problem equation at each iterative step can be rewritten as A W A s W = AtW. (6.4) To simplify the presentation in this thesis, a straight line ray path is assumed so that A is independent of m. The method studied in this thesis can be readily extended to the non-linear case by considering the linear problem as one step in solving the non-linear case. For a matrix A which is M x N, the number of model parameters M is usually assumed to be greater than the number of data N, an assumption which can always be made by refining the model space in the discrete case. Such a problem is said to be underdetermined and is satisfied by an infinite number of solutions, i.e. the problem is non-unique. To find a unique, physically meaningful solution, a priori information is needed to constrain the solution. Sometimes, even if A is of full rank, the condition number which is defined as the ratio of the largest singular value to the smallest one can be very large and the problem is ill-posed. For an ill-posed problem, the noise in the observed data can cause some unrealistic structures in the inverted model. Again, some prior information about the model should be used to regularize the solution. In either case, a particular solution can be obtained by minimizing an objective function which can be either derived using a Bayesian approach or from the viewpoint of regularization. Chapter 6. INVERSION WITH PRIOR 6.3 117 INFORMATION The maximum a posteriori probability (MAP) solution Let us assume that the traveltimes are contaminated with Gaussian noise and the model parameters are also Gaussian random variables. The maximum a posterior probability (MAP) is maximized when the objective function (Tarantola, 1987) ^(m) = ( A r a - t f C - ( A m - t ) + m 1 T n C> is minimized, where C is the estimated data covariance matrix and C n (6.5) m is the a priori model covariance matrix. The inverses of the covariance matrices are called data and model weighting matrices respectively. The minimization of the objective function leads to m = (A C' A T 1 + C- )" A C - \ 1 1 (6.6) T and can also be written as (Tarantola, 1987) m = CA T m (AC A T m + C) n t, (6.7) where the covariance matrix C is applied instead of the weighting matrix C^ to con1 m struct the solution. Ideally, if we construct a covariance matrix according to the available prior information, a solution can be found by solving any of the above systems of equations under the Gaussian assumption. The inversion of C and C can be computed inexpensively under certain conditions. n m Since noise is usually assumed to be independent Gaussian, the covariance matrix is diagonal and C" is computed by simply inverting the diagonal elements. The inverse 1 of C m can also be computed without any difficulty if the model is a stationary Gaussian process, i.e., the model parameters may be correlated but the covariance function is the Chapter 6. INVERSION WITH PRIOR INFORMATION same for all the model parameters. In this 118 Fourier approach can be used to compute the inverse of the covariance matrix (Andrews and Hunt, 1977; Pilkington and Todoeschuck, 1992; Squires and Guillaume, 1992). The difficulty of obtaining the above solutions is the inversion of ^ A C ~ A + C " ) T 1 1 or ^ A C A + C ) in equations (6.6) and (6.7) respectively. For large A , with size T m n over 1000 x 1000, for instance, a direct inverse approach such as Cholesky decomposition (Press et al., 1992) is too slow. In this case, the conjugate gradient approach can be used to solve the system of equations (6.8) given below efficiently by taking the advantage of the fact that A is a sparse matrix (Nolet, 1987; Scales, 1987; Claerbout, 1992; Aldridge and Oldenburg, 1993) 1 C„"*A where C " = C 1 C 2 m 2 TO _ m= and C " = C 1 2 n i cn~n (6.8) 0 C . The equation (6.8) can be easily con2 n firmed to be agreed with the solution in equation (6.6) by multiplying J A C T _ n 2" C ~ 2"] m to both sides of equation (6.8). In the case that the model is independent Gaussan distribution, thus C ^ is diagonal and C 1 2 TO can be computed by taking the square roots of _ the diagonal elements. In more general case where the model is fractal, C 1 2 m is usually not computed directly but approximated using a combination of an identity and a few derivative matrices which will be discussed in the next section. Finally, it is worthwhile emphasizing that the above approach is suitable only for stationary models for which C ^ can be obtained using the Toeplitz structure of C . 1 m The inversion of non-stationary models will be studied in the next chapter using a wavelet related approach. Chapter 6. INVERSION WITH PRIOR INFORMATION 6.4 119 The relationship between regularization and stochastic approaches 6.4.1 On the covariance matrix, damping and smoothing The regularization method uses the general model norm in the following form, = a m m + a 9m 9m + a d m d m r (j) m r 0 2 1 T +••• 2 2 , (6.9) where d denotes partial derivative with respect to space and a , ai, a , • • • are the con0 2 stants that control the relative weights among different terms on the right hand side of the equation. Equation (6.9) can be written in matrix form as <> /m = «oni m + a i m D f D i m + a m D ^ D m + • • • r T T 2 2 , (6.10) where D D , • • • are the discrete matrix differential operators which will be examined at 1 ; 2 the end of this section. In terms of regularization, the first term on the right hand side of equation (6.10) is the damping term and the rest of the terms are related to different degrees of smoothing. If only the first term is used, the inverted model is called the smallest model; if only the second term is used, the model is called the flattest model; if only the third term used, the model is called the smoothest model. By comparing the generalized model norms in equations (6.5) and (6.9), we can represent C " using a series of matrices 1 C - = a I + CiDfDi + a D^D + • • • , 1 0 2 2 (6.11) where I represents an identity matrix. The above equation relates the regularization methods to stochastic inversion. This relation has been shown by Tarantola (1987) and the obvious application is to combine damping and smoothing to invert a correlated model and this may save much time in comparison to computing C ^ directly. The question 1 Chapter 6. INVERSION WITH PRIOR INFORMATION 120 now is how to determine the constants do, ai, and so on. Wang and Lawrence (1994) have applied this relation to traveltime tomography and they have chosen the constants using a trial and error approach. In what follows, a simple approach is introduced to approximate C ^ using a finite number of terms on right hand side of equation (6.11). 1 For a stationary model, the covariance matrix is Toeplitz and thus can be represented by a covariance function vector. The inverse of the covariance matrix, the weighting matrix, can be computed using a Fourier approach as mentioned in the previous section. In this case, the weighting matrix can be represented by a vector q, which will be named as weighting vector or weighting function in short, and equation (6.11) becomes q= + a [---,0,l,0,---] r o aa[-.-,-1,2,-1,-..} T + a [...,-l,4,-6,4,-l,..-] r 2 (6.12) + ••• = Pa, where matrix P relates q to the constants denoted by vector a. i If the elements of q decay rapidly from the center of the vector, a truncated q can be used to approximate the original one. For example, if two elements on each side from the center of q are considered, equation (6.12) can be represented as 9o ?2 1 1 -1 0 -1/2 4/6 0 0 -1/6 a 0 d2 where only half of the matrix is shown, and the constants can then be easily obtained. A trivial example is that when the model is an uncorrelated white process, q = a = 1. The above approach requires the weighting function to decay quickly. Fortunately, Chapter 6. INVERSION WITH PRIOR 121 INFORMATION this is true in many cases, for example in the case of 1/f model which will be discussed next. 6.4.2 Regularization of a 1/f model using damping and smoothing For a stationary 1/f type process, the power spectrum is proportional to l/|u>| , thus a the power spectrum of the weighting function is proportional to |u>| . The weighting a function in the spatial domain is the inverse Fourier transform of \oj\ (Pilkington and a Todoeschuck, 1992) which can be written in a form of q(t) = sin *f T(a + l ) / | * | a/0,2," a + 1 for - (Bracewell, 1978; Rade and Westergren, 1995). For a > 0, q(t) decays much more quickly than the covariance function and it is therefore possible to truncate q(t). In the discrete case, q can be represented using a combination of the column vectors of P in equation (6.12). Therefore C " can be approximated using a combination of an 1 identity and a few derivative matrices. As an example, Figure 6.1 shows the covariance matrix C m when a = 1, its inverse (weighting) matrix Q , and the corresponding covariance and weighting functions. m The inverse covariance matrix can be approximated using the damping and smoothing matrices as in equation (6.11). Table 6.1 shows the constants used for approximating matrix Q m with different orders of smoothing. Since the model norm in equation (6.10) should be minimized, a large constant a implies a greater penalty to the model norm, so 0 that the model will have less structure; large constants a\, ai, • • • imply greater penalties on derivatives and thus the model will be smoother. In the table, as the number of terms increases, a decreases, therefore there will be more structure in the model. Since, in 0 the meantime, the values of ai, a , • • • increase, more smoothing is introduced reducing 2 the structure. This intuitively suggests that using more or less terms should give similar results. Some synthetic tests will be performed in the next chapter to clarify the effect of adding higher order smoothing. Chapter 6. INVERSION WITH PRIOR INFORMATION terms a 2 0.4143 0.2667 3 4 0.2009 0.1638 5 0 ai «3 a 0.5857 1.1762 0.4428 1.7678 1.6262 0.6574 2.3626 3.8566 3.6313 122 2 1.3011 Table 6.1: Estimated damping and smoothing coefficients for a 1/f model when a = 1. 6.4.3 Conjugate gradient approach with damping and smoothing With the constants a computed in the previous subsection, equation (6.8) becomes cn~h C "»A n 0 m= D where C 2 m 2 0 0 is substituted by an identity matrix and a few derivative matrices. The smallest, flattest and smoothest models can be obtained by using only I, D j and D 2 separately. The above equations can also be solved with a certain combination of I, D i , etc. to yield a certain stochastic model. This approach can be applied in both the 1-D and 2-D cases. The forms of 2-D derivative matrices could include 1-D derivatives in both horizontal and vertical directions as well as cross derivatives. Some of the discrete derivatives are listed in the next subsection which will be used in the next chapter. Chapter 6. INVERSION WITH PRIOR 6.4.4 INFORMATION 123 The discrete derivative matrices The derivative matrices are important in solving inverse problems. First, the smoothing model constraints can be directly incorporated using the derivative matrices. Secondly, the covariance matrix can be approximately represented by the combination of the damping and smoothing matrices. Different damping and smoothing constraints will be compared with the wavelet domain constrains in scale and location in the next chapter. Here the forms of the derivative matrices are discussed for both 1-D and 2-D cases. 1-D Suppose N = 8, the first, second and third derivative matrices are 1 D -1 1 1 -1 (6.15) = J 1 2 -1 0 -1 D, and 7x8 0 0 0 0 0 2 -1 0 0 0 0 2 -1 0 0 0 2 -1 0 0 2 -1 0 0 0 -1 0 0 0 -1 0 0 0 0 -1 0 0 0 0 0 -1 2 -1 (6.16) 6x8 Chapter 6. INVERSION WITH PRIOR INFORMATION -1 3 -3 0-1 D = I 3 1 0 0 0 0 1 0 0 0 3 -3 1 0 0 3-3 0 0 -1 0 0 0 -1 0 0 0 3 -3 0 -1 124 (6.17) 1 0 3 -3 1 5x8 Note that the ranks of the ith derivatives D ; and D f D; are M — i and that the above derivative matrices have not included the derivatives at the edges of the model which could be included with certain modifications. 2-D For a 2-D model, if the model is parameterized into M x M elements, the first order derivative in x and z directions can be written as 1 -1 1 Di m = i!B -1 1 (6.18) m -1 21 m2 -1 1 2 -1 (Af-l)MxAfAf MMxl Chapter 6. INVERSION WITH PRIOR 1 0 -1 0 1 0 Di m = z 1 INFORMATION 125 -1 0 ••• -1 ••• 1 -1 0 ••• . (6.19) mi 2 m2 2 (M-l)MxMM J MMxl For a mixed derivative, the derivative can be written as 5 2 m dxjdzj _ (dm\ \9ziJj _ / 3 m \ \8zi)j i (6.20) + = rriij - m j - rriij i+1 +1 + m ji i+1 + where i and j are the indices of an element in x and z directions. 6.5 Summary The main purpose of this chapter was to prepare the theoretical background for inversion of stationary models which will be used in the next chapter for inversion of more general, non-stationary models. The prior model information was incorporated by minimizing a general model norm constructed using a weighting matrix Q = C^ as the inverse of the covariance matrix. 1 m For a stationary model, covariance and weighting matrices can be represented by covariance and weighting functions respectively. Instead of direct inversion of a covariance matrix, the weighting matrix can be constructed by the weighting function which can be computed using a Fourier approach. In cases that the weighting function decays rapidly, Q m can be constructed approx- imated using a combination of an identity and a few derivative matrices which function Chapter 6. INVERSION WITH PRIOR INFORMATION 126 Cm 25 30 Qm 35 40 25 30 35 40 Figure 6.1: Top: covariance matrix C and weighting matrix Q for a fractal model with fractal parameter a = 1. Bottom: the covariance and weighting functions. TO m as damping and different degrees of smoothing. The relative weights on the damping and smoothing can be computed given the number of terms to approximate the weighting function. The covariance matrix of a stationary 1/f type model was used as an example. For a large system of equations, the CG approach can be implemented via damping and smoothing to incorporate the prior model constraints. The approaches discussed in this chapter can only be applied to incorporate stationary information. The next chapter will describe how to incorporate the non-stationary, timevariant information to the inverse problem with the help of the wavelet transform. Chapter 7 WAVELET T R A N S F O R M INVERSION (WTI) 7.1 Introduction The previous chapter reviewed some conventional approaches of incorporating stationary prior model information to constrain the inversion. This was done by minimizing a generalized model norm constructed using a weighting matrix which is the inverse of the covariance matrix. For a stationary model, each model parameter is treated statistically in an identical manner and therefore the covariance and weighting matrices can be represented by the covariance and weighting functions (vectors). But for a non-stationary model, such as the one discussed in Chapter 5, the correlation function changes from one parameter to another and the approaches in the previous chapter cannot be applied to the inversion of this type of model. Recently a new approach that uses the wavelet transform to incorporate non-stationary scale information has been studied in medical image tomography (Bhatia, 1994; Kolaczyk, 1994). This chapter will apply this new approach to incorporate prior scale information in seismic tomography. Two types of scale information will be discussed. The first type of information relates to the size of the objects to be inverted. In order to attenuate the effect of noise in an inverted model, a model with the maximum possible scale is discussed. The second type includes the fractal parameter of a non-stationary fractional Brownian motion process. A large portion of the contents of this chapter have been published recently (Li et al., 1994; Li and Ulrych, 1995; and Li et al., 1996). 127 Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 128 The next section introduces the theory of wavelet transform (constrained) inversion (WTI) as a special way of incorporating prior information in a transformed model space. Section 3 describes two particular types of models constructed using the WTI approach and section 4 discusses three examples of using the WTI. 7.2 Wavelet transform inversion (WTI) 7.2.1 Formulation Assume that slowness m is a F B M process which can be characterized by a covariance matrix C . Since the rows of C m m are not the same, it cannot be represented by a single covariance function as in the case of a stationary process and thus, the approaches discussed in the previous chapter cannot be applied to solve the problem. As discussed in Chapter 5, a F B M covariance matrix can be diagonalized using the wavelet transform and this fact can be used to simplify traveltime tomographic inversion and also to incorporate non-stationary F B M constraints. We first derive the form of the objective function in terms of the wavelet transform coefficients of the model parameters. As shown in Chapter 2, an orthogonal wavelet transform can be represented by a matrix W with W W = I, T W- =W . X T (7.1) According to section 5.3.1 the covariance matrix of the wavelet coefficients is C = WC W , r r o m (7.2) and its inverse therefore can be written as (7.3) Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 129 or, because of the orthogonality of W, c -1 = w c w. r _1 (7.4) Assuming in = Wm, equation (6.2) can be written as A W W m = t or T (7.5) Am = t. The objective function in equation (6.5) can then be rewritten in terms of m as follows, </>(in) = (Am - t ) C ; ( A m - t) + m^C^in. T 1 (7.6) The corresponding MAP solution for the above objective function is (7.7) or in = C A T m (AC A + C ) r m n _ 1 1. (7.8) After obtaining the wavelet coefficients, the model parameters can be computed using the inverse WT. For a large problem, the CG method can be applied to seek the solution of c i 2 m = " cd-h (7.9) 0 As discussed in Chapter 5, the velocity field can be generally modeled as a F B M process. Thus, C is approximately diagonal with the diagonal elements following power m law and C 2 m can be computed by taking square roots of the diagonal elements. As Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 130 another example, if the scale or the size of object to be inverted is approximately known, this information can be incorporated into the inversion by applying a small penalty to the corresponding diagonal elements of the covariance matrix. The above two models will be examined in more detail in the next sections. 7.2.2 On the relationship between W T and Fourier domain constraints The Fourier transform has similar properties to those of the wavelet transform described by equations (7.1), and the equivalent expressions are FF H = I, F _ 1 = F , H (7.10) where F represents the Fourier transform. For a stationary model, the covariance matrix is Toeplitz and the covariance matrix of the Fourier coefficients is a diagonal matrix with the diagonal elements being the Fourier transform of the correlation function (Strang, 1986). This feature has been used in seismic tomography by Pilkington and Todoeschuck (1992) and Squires and Guillaume (1992). However, since the Fourier transform cannot diagonalize the covariance matrix of non-stationary process, it is not suitable as a means of incorporating the non-stationary model constraints. 7.3 7.3.1 Two types of multiscale constrained models The scale-truncated model Choosing the number of model parameters can sometimes be a problem in solving an inverse problem. On the one hand, by choosing a small grid interval, the problem becomes highly ill-posed and unrealistic structures could dominate the solution without the proper constraints. On the other hand, if the grid is not dense enough, we risk the sacrifice of the resolution which the data could provide. Using wavelet transform inversion we should be Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 131 able to find a solution with the largest possible grid interval that satisfies the required data misfit. For convenience, this solution is called the scale-truncated model, or the largest-scale model. The solution can be found by the means of the following steps. First, equation (7.5) is rewritten in a multiscale form, Jo *J0 "•Jo 30 L = t. k-J-1 Jo+l (7.11) dj_i The second step is to find a scale L so that the least squares solution of 'Jo A* ( } A .. = t. 0 0 fits the data to %\ ^ ^> D u t a t s c a ^ L + l the solution of equations e (7.12) Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 132 'Jo *J0 Jo 0 Jo (7.13) 0 0 0 fits data to xi+i < which means that we have overfitted the data by introducing scale L + 1. Equation (7.13) can be rewritten as W t, (7.14) At , (7.15) T . dr+i or A L + 1 d L + 1 = t - A^Si = L where SL is the minimum least squares solution of equation (7.12). Now, d^+i can be obtained by minimizing the following objective function <f> L+1 = (A i + 1 d - At ) C; (A r L + 1 L 1 i + 1 d - At ) + d C^ d , T L + 1 L L+1 L+ L+1 (7.16) subject to misfit % = N. The largest-scale solution can be obtained by the inverse 2 wavelet transform of [s^, d i i , 0, • • •, 0]. In Example 1, a largest-scale solution will be + compared with solutions by other inversion approaches. 133 Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 7.3.2 The fractal model As discussed in Chapter 5, well logs can be decomposed into two components with different features: one component can be characterized by the fractal model and the other comprises the trend and spiky characteristics which cannot be statistically described. In this chapter, to simplify the problem, the spiky features will not be considered and so the smoothed coefficients (defined in Section 2.3) are assumed to represent the trend while the detailed coefficients represent the fractal component. For a non-stationary F B M process, C is dense and non-Toeplitz, thus computing m the inverse matrix, C ^ , is very expensive. Since the wavelet transform can diagonalize 1 the covariance, the inversion of the covariance matrix becomes very simple. Suppose prior F B M model information can be represented using the covariance matrix C . Since m the covariance of the wavelet coefficients is assumed to be diagonal we can write C « diag{<r u , <r u , <r 2 m 2 s jo o 2 jo o+1 u jo+1 , • • •, aj^Uj^}, (7.17) where diag denotes a diagonal matrix, Uj denotes a vector of length 2 with all elements J equal to unit, a is the variance of the smoothed coefficients and a is the variance of 2 2 the detailed coefficients at scale j. Substituting equation (5.5) into the above expression one obtains 2 C « a diag{^2-^u , cr 2~ ^u , 2 - ^ u 2 m a jo a jo , • • •, 2~ ^uj^}, a j o + 1 (7.18) where a = a Vp(a) and Vp as a constant that depends on the chosen wavelet tp (Flandrin, 2 2 1992). The inverse matrix can be computed as C - = ^ diag{^2 "u , 1 a j0 2^u , 2 ^ jo + 1 ) , • • •, 2 " ( - ) _ } , J U j o + 1 1 Uj 1 (7.19) Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 134 which can be used to construct the solution directly. Another possibility is to compute ^ _ i_ Cm and to obtain the solution by the means of the conjugate gradient. 2 7.3.3 W T I and the damped least squares When a = 0 and cr = cr, = tr, all scales are equally weighted, C 2 m = ^fT = /zl. Hence, the solution, equation (7.7), corresponds to the DLS. This can be shown using the relations A = A W , W W = I, W T T =W r - 1 . The solution becomes m = (WA AW + uWW )- WA t ^ T T r V = x J T (7.20) W ( A A + /iI) A t, T _1 T which in the original model space becomes m = ( A A + I)- A t, r 1 Ai :r (7.21) and gives the DLS solution. 7.4 Example 1: inversion of a blocky model A simple vertical seismic profiling example is used to determine how the a priori information can be incorporated into the inversion of the blocky model. In this case, the shot is assumed on the top of the well and the receivers are assumed to be located at the bottoms of modeled layers. The matrix A that relates the model parameters m and the observed traveltimes t, in this case, is lower triangular. Assume M = N = 256, where M is the number of the model parameters and N is the number of the traveltimes. The traveltime is assumed to be contaminated by Gaussian noise with standard deviation o~ . n The existence of this noise has very different effects on the inverted model at different scales. Traditionally, this problem is tackled using eigenspectral (for example, Stork, Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 135 1991) or Fourier domain (Pilkington and Todoeschuck, 1992) constraints. The following example compares the WTI largest-scale model with some conventional models. For all the solutions, the final % misfit is equal to the number of observations. 2 The model and its first and second derivatives are plotted in Figure 7.1(a). The Fourier spectrum and the wavelet transform coefficients are shown in Figures 7.1(c) and (d). In what follows, the model norms will first be constructed using the model parameters and their derivatives. 7.4.1 Inversions via damping and smoothing Figure 7.2 is the damped least squares solution of the smallest model. In this case, if the damping factor is small and the data are overfitted, there exist many unrealistic structures and the model norm is large. If the damping factor is too large and the data are underfitted, less structures are inverted and we lose model resolution. The smallest model is found with \ 2 = N. It is obvious that the small scales are overfitted. To attenuate the small structures, which are assumed to be unrealistic by overfitting the noise, we force the derivatives, the differences between the neighboring cells, to be small by minimizing the model norm of the derivatives. The results are shown in Figure 7.3. The same \ 2 misfit as in Figure 7.2 is achieved. Compared with the smallest model, the flattest and smoothest model have reduced the high frequency structures significantly, but have also reduced the sharpness of the edges of the objects. 7.4.2 Inversion via constraints in a transformed domain Instead of constraining the model parameters directly, we can constrain the frequencies or scales of the model parameters in a transformed model space. Based on the fact that the small structures in the model are affected more seriously by the noise than are the Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 136 * large scales, a big penalty will be applied to small scales, orTiigh frequencies, to yield models with largest possible structures with the same x data misfit. 2 Figure 7.4(a) shows the number of singular values versus the data misfit. Since the singular vectors corresponding to the small singular values will be affected by the noise the high frequency structures constructed by these singular vectors (Stork, 1991) are less reliable. If all the singular values are used, we can find a model that fits the data exactly but with a lot of unrealistic high frequency (small scale) structures. By discarding some small singular values, we can find a model that fits the data with reasonable x error 2 and also with less unrealistic small scale structures. Figure 7.4(b) shows the model with the first 26 singular values with the same data misfit as in the previous models using the model space constrained approaches. Figure 7.5 shows the Fourier space constrained model. There are 26 sinusoids used to construct the model which fits the data to the same % misfit as in Figure 7.4. Because 2 the basis functions are of infinite length, the low frequencies suffer from sidelobe and edge effects. The wavelet space constrained inversion can overcome the above fundamental problems. First in the WT space, the small structures can be reduced by introducing large damping factors to the coefficients at small scales. Since the basis functions are localized, there is no sidelobe effect when representing the local structures. Figure 7.6 shows the wavelet transform constrained model. The scale-truncated model is found in two steps: • Find scale L so that the \ 2 misfits at scales L and L + 1 satisfy XL — N ^ AN< Xl <N. +1 • Solve an overdetermined problem at scale L and keep the model coefficients. Then solve rn( ) by searching for the damping factor p that fits % = N. i+1 2 Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 137 Since each wavelet coefficient is localized, the location constraints can be easily incorporated at the same time as the scale information, while the K-L and Fourier approaches cannot. Figure 7.7 shows the solution with both the prior scale and location information. In this case, a large damping factor is applied to the wavelet coefficients representing the right hand side of the model. To summarize, the solutions in Figures 7.2 and 7.3 can be incorporated with prior spatial but not frequency (scale) information, while the solution in Figures 7.4 and 7.5 can be incorporated with prior frequency but not spatial information. Only using the WTI approach, both the prior scale and spatial information are incorporated to the solution in Figure 7.7. 7.5 Example 2: inversion of a F B M model The study of velocity logs shows that the velocity model consists of two parts: nonfractal and fractal components. The inversions of these two components are done in the following two steps. • Invert the trend at large scale by solving an overdetermined problem. To do this, the smoothed coefficients are inverted by minimizing the data misfit. • Solve for the detailed coefficients by minimizing the model norm with F B M information and fitting % = N. 2 Figure 7.8 shows the synthetic slowness model using the approach introduced in section 5.5. Shown on the left are the wavelet coefficients where the detailed coefficients are generated from a random series scaled by a power law and the 8 smoothed coefficients on top of the diagram are numbers with a fix decrement. Shown on the right are the Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 138 synthetic slownesses and the multiscale components. Figure 7.9 shows the decomposed trend and the F B M components. The solution found using the above steps is shown in Figure 7.10(a). Comparison of the WTI results with the conventional damped least squares (DLS) result in Figure 7.10(b) shows that the WTI solution does not contain unrealistic structures which are associated with small scales and which were not constrained by the DLS inversion. This can be seen more clearly in the scale-location space in Figure 7.11. Comparison with the model coefficients shown in Figure 7.8(a) indicates that, although both methods fail to recover the exact features at the coarse scales, the WTI has constrained the spurious structures at small scales leading to better agreement with the original synthetic model. 7.6 Example 3: inversion of large problems In the 2-D case, the number of model parameters is generally quite large and the problem is usually solved using the conjugate gradient approach. The example in this section compares the WT constrained CG technique described in section 7.2 with the conventional CG approach in section 6.3. The ray geometry of the problem is shown in Figure 7.12. The example simulates a cross borehole survey with surface receivers. There are 16 shots with 32 receivers and 64 x 64 cells. A model with three blocks is shown in Figure 7.13. The models are displayed in both 1-D and 2-D manners. The 1-D diagrams on the left plot all the model parameters as a 1-D vector to show a clear visual comparison between the inverted and the true models. Because we are depicting more than four thousand model parameters, the boxcar functions appear as vertical lines in the diagrams. Inversions with different constraints will be shown in the following diagrams. The model has three velocity anomalies with a background velocity of 1000 m/s which is Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 139 used as the initial velocity for the inversions. All the inversions have the constraints that the cell velocities are known in the two wells and at the surface. This constraint is incorporated by assigning large weights to the appropriate cells. 7.6.1 Model space constrained inversions Figure 7.14 shows the conventional smallest DLS solution. There are many small scale structures in the inverted model. Smoothing constraints, damping the derivatives, can be used to reduce these unrealistic structures. Table 7.1 lists some of the derivatives that can be used as smoothing constraints in the 2-D case. The first column represents the constraints in the a; direction which have been used in Figure 7.15. The third column represents the constraints on the partial derivatives in both x and y directions which have been used in Figure 7.16. By introducing constraints on the derivatives in the horizontal direction the small scale unrealistic structures have been reduced. The better result in Figure 7.16 is achieved by introducing mixed higher order partial derivatives. Again, all the solutions are fitted to % = N. The problem with the model space constrained 2 inversion is that by increasing the degree of the smoothing the edges of the objects become more blurred. This problem is tackled using WTI which can be considered as a multiscale smoothing approach. 7.6.2 WTI solutions Figure 7.17 shows the 1-D WT scale constrained inversion and Figure 7.18 shows the solution with 2-D WT scale constraints. The results described above have been dramatically improved by reducing the effect of noise on the small scale structures and preserving the edges of the objects. A large damping factor has been used at the small scales which are Chapter 7. WAVELET TRANSFORM INVERSION (WTI) X first order second order 140 xy y a_ dx dy- a d dx dy 2 2 2 2 fourth order a 2 dxdy a 4 dx dy 2 2 Table 7.1: Derivatives which could be used as model constraints. known to be affected by the noise more seriously. Both solutions are fitted to the same data misfit as in the damping and smoothing approaches. 7.6.3 Model correlation matrix As was discussed in the previous chapter, a linear combination of the different derivatives can be used to define the correlation function of the model. This correlation function relates one particular model parameter with neighboring ones. By defining the weights at different scales and locations, we also define a covariance matrix for the WTI approach. Figure 7.19 shows examples of correlation functions, the rows of the covariance matrix used in the damping plus smoothing and WTI approaches respectively. The covariance matrices in the damping plus smoothing approaches are obtained by inverting the weighting matrices which are the combination of the damping and smoothing matrices. The covariance matrices in the WTI approach are computed by transforming the wavelet domain covariance matrix back to the original model space. Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 7.7 141 Summary WTI was developed so that different types of scale and location information may be incorporated into a tomographic problem. A largest-scale model was inverted by truncating the small scale wavelet coefficients in model space and fitting the data to % = N. Comparing with the results of damping 2 and smoothing regularizations, the largest-scale model is more reliable and reconstructs edges well. Comparing with the spectral constrained approach, because of the localized basis used in the inversion, the largest-scale model is less affected by sidelobes. The prior information related to the location can be easily incorporated into the largest-scale model to improve the result. Chapter 5 demonstrated that a well log can be decomposed into a combination of a non-stochastic component and a component that can be modeled as fractional Brownian motion. In this chapter, a new procedure was developed to invert these two different components. First, a non-stationary large scale model is inverted by solving an overdetermined problem. Then the detailed coefficients are inverted with wavelet domain constraints. Section 7.2.6 showed that without the scale and location information, the WTI approach becomes a damped least squares approach. With the right scale and location information, however, the WTI approach will, in general, yield a better result than conventional damped least squares. This has been shown in this chapter using three different examples. Both scale and location information can be incorporated into an inversion using WTI. This is superior to both the Fourier domain approach and model space regularization using damping and smoothing constraints. Using model space regularization certain location penalties could be incorporated into the inversion but not the scale or frequency Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 142 information of the model, while the frequency domain approach achieves the opposite in that frequency but not location information can be incorporated. WTI can incorporate both types of information jointly. For a 2-D inversion of large problems, the conjugate gradient approach can be used to incorporate the wavelet constraints more efficiently. A 2-D example was used to test the wavelet domain CG approach and the results were compared to the conventional model space constrained CG approach. Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 143 (b) Derivatives (a) Model 1 6000 D1 lo 5000 E £-4000 o o D2 £3000 100 200 Depth (sample) 100 200 Depth (sample) 2000 (d) WT coefficients (c) Fourier spectrum J I J I I I _g) 4 t o cn 5 6 7 10 Frequency (sample) 10 Figure 7.1: (a) Model 1; (b) its derivatives; (c) its Fourier transform; and (d) its wavelet transform. (b) DLS (a) Misfit 800 1000 2000 3000 Damping factor 4000 100 200 Depth (sample) Figure 7.2: (a) % squared misfit versus damping factor; (b) the smallest model. 2 Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 144 (b) Smoothest (a) Flattest 3000 2000 100 200 Depth (sample) 0 100 200 Depth (sample) Figure 7.3: The (a) flattest and (b) smoothest models. (a) Misfit (b) S V D 320 25 30 35 Number of sigular value 100 200 Depth (sample) Figure 7.4: (a) % squared misfit versus the number of singular values; (b) the model inverted using the first 26 singular values. 2 Chapter 7. WAVELET TRANSFORM INVERSION (WTI) (a) Mistfit 145 (b) FT I 280 260 240 220 20 25 30 35 Number of freq 100 200 Depth (sample) Figure 7.5: Fourier transform constrained inversion (FTI): (a) % squared misfit versus the number of frequencies; (b) the model inverted using the first 26 frequencies. 2 (a) Mistfit (b) WTI Scale Depth (sample) Figure 7.6: Wavelet transform constrained inversion (WTI): (a) % squared misfit versus the number of scales; (b) the model inverted using the first 4 scales. 2 Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 146 WTI with location 6000 v> 5000 E, ^4000 o o £3000 2000 0 100 200 Depth (sample) Figure 7.7: WTI with location constraints. The deeper portion of the model is assumed to have less structure and a large penalty is applied for the corresponding coefficients. Figure 7.8: Model 2: F B M model with a = 1.4. (a) Wavelet coefficients; (b) the multiscale signals and the synthetic model. A linear trend is added to the model. Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 148 Chapter 7. WAVELET TRANSFORM INVERSION (WTI) (a) WTI (b) D L S Figure 7.10: (a) WTI solution; (b) damped least squares solution. The thin lines in both figures show the true model. (b) D L S (a) WTI J I .I I J L 1 I i r- j I ' 1 1 . I, .. •-i"T— 0.5 L i i_ ' | I I I. . . • ... 1 II i 'I r 0.5 Figure 7.11: Wavelet coefficients of the inverted models in the previous figures. Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 149 Distance (m) Figure 7.12: Cross borehole geometry. Model index Distance (cell number) Figure 7.13: Model 3. There are 64 x 64 = 4096 cells and the model is composed of three blocks. On the left, all 4096 model parameters are plotted in 1-D style. The sizes of the two larger blocks are both 16 x 16 cells and the velocity anomalies are —100 m/s and 200 m/s respectively. The size of the smaller block is 8 x 8 cells and the velocity anomaly is —200 m/s. Chapter 7. WAVELET TRANSFORM INVERSION (WTI) model index 150 Distance (cell number) Figure 7.14: Inversion using damping constraints. On the left: the result is plotted in 1-D style and the synthetic model is shown as the dotted line. On the right: the result is plotted in gray scale, from white to dark the amplitude ranges from 850 m/s to 1150 m/s. 1200 model index Distance (cell number) Figure 7.15: Inversion using damping and x direction smoothing constraints. See Figure 7.14 for explanation. Chapter 7. WAVELET TRANSFORM INVERSION (WTI) model index 151 Distance (cell number) Figure 7.16: Inversion using damping and 2-D smoothing constraints See Figure 7.14 for explanation. 1200 model index Distance (cell number) Figure 7.17: Inversion using 1-D WT constraints in horizontal direction. See Figure 7.14 for explanation. model index Distance (cell number) Figure 7.18: Inversion using 2-D WT constraints. See Figure 7.14 for explanation. Chapter 7. WAVELET TRANSFORM INVERSION (WTI) 5 10 15 Distance (cell number) 152 5 10 15 Distance (cell number) Figure 7.19: Covariance functions of model and wavelet space constrained inversions. One row of the covariance matrix (4096 x 4096) is plotted in 2-D. For display purposes, only 16 x 16 out of 64 x 64 cells are plotted. Left: 2-D model space constraints. Right: 2-D WT constraints. Each diagram relates the model parameter at the 4th row and 5th column to neighboring model parameters. Chapter 8 SUMMARY Conventional approaches analyze and process a signal in either location or scale. For example characteristics of a seismic trace can be represented in terms of amplitudes in time or in frequency, but not in both at the same time. The wavelet transform, on the other hand, provides a tool to examine the characteristics in both location and scale at the same time. It turns out that this flexibility of the wavelet transform is of great importance in both the processing and inversion of data in many fields of science. This thesis described the application of the concept of the wavelet transform to two commonly encountered geophysical problems: (1) signal analysis and processing, and (2) velocity inversion. The thesis can be subdivided into two parts according to the types of data dealt with and the wavelets used. The first part consists of Chapters 3 and 4 where the band-limited seismic trace was dealt with directly using the Gabor and Morlet transforms in order to provide location information in addition to frequency/scale information. The second part consists of Chapters 5, 6 and 7 which dealt with the modeling and inversion of seismic velocity fields. Orthogonal wavelet transforms were adopted to map information completely and efficiently from the spatial domain to the scale-location domain. To summarize, this thesis consists of the following contributions. 153 Chapter 8. SUMMARY 154 Coherent noise attenuation with location information In Chapter 3, a localized transform, the Gabor transform, was used to attenuate coherent noise, such as ground roll, using the location information of the signal and noise as prior information. It was shown that a Gabor band pass filter without location information is equivalent to using a tapered Fourier band pass filter. A synthetic example in 1-D demonstrated that the use of location information can dramatically improve results of noise attenuation. A 2-D Gabor transform approach was developed to incorporate the location information of the ground roll and illustrated by application to a shot record. Analyzing the localized seismic attributes In Chapter 4, some wavelet attributes were defined using the Morlet wavelet transform and the meaning of the wavelet attributes and some preliminary applications was explored. Specifically, conventional instantaneous attributes such as those defined by Taner et al. (1979) analyze certain features at a certain time for all frequencies, while wavelet attributes analyze these features at a certain time at various scales (or frequencies). The usefulness of wavelet attributes was shown by analyzing the location and frequency distribution of ground roll, reflection and refraction events. Modeling of well logs Since the Fourier transform cannot be used to analyze variations in location of the characteristics of a signal, it cannot be used to directly model well logs which are highly non-stationary. With the help of wavelet transforms, well logs can be modeled in the wavelet domain directly as a combination of three components: the trend represented by the smoothed coefficients, a non-stationary fractional Brownian motion component represented by most of the detailed coefficients and some very spiky detailed coefficients Chapter 8. SUMMARY 155 which cannot be statistically modeled. Estimates of fractal parameters for some well logs confirmed the F B M nature of well logs and the variation of fractal parameters represents the geological variation associated with the locations of the wells. For a given fractal parameter, both 1-D and 2-D fractals were synthesized from a random series or image. This procedure could be further used for geological simulation. Wavelet transform inversion with scale and location information In Chapter 7 a wavelet transform inversion approach which can incorporate both scale and location information was developed. The WTI approach is superior to both the conventional model space approaches using damping and smoothing constraints, as well as to the Fourier domain approach. Using a conventional damping or smoothing approach, the location information, but not the scale or frequency information, can be applied to the inversion. Frequency, but not location information, can be incorporated into the inversion by the use of a transform approach, for example by means of the Fourier transform. WTI can incorporate both types of information at the same time. Two particular types of models, the largest-scale and F B M models, were discussed in Chapter 7. Three examples were used to test the WTI and they all show a marked improvement over conventional methods. The WTI was also applied to solve a large inverse problem by means of the conjugate gradient method. References Aldridge, D. F., and Oldenburg, D. W., 1993, Two-dimensional tomographic inversion with finite-difference traveltimes: J. of Seis. Explor., 2, 257-274. Andrews, H. C , and Hunt, B. R., 1977, Digital Image Restoration, Prentice-Hall Signal Processing Series. Bastiaans, M.J., 1981, A sampling theorem for the complex spectrogram, and Gabor expansion of a signal in Gaussian elementary signals: Optical Engineering, 20, 4, 594-598. Bhatia, M . , Karl, W. C., and Willsky, A. S., 1996, A wavelet-based method for multi-scale tomographic reconstruction: IEEE Transactions on Medical Imaging, 15, 92-101. Bosman, C., and Reiter, E., 1993, Seismic data compression using wavelet transforms: 63th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1261-1264. Bradley, J., Brislawn, C., and Hopper, T., 1993, The FBI wavelet/scalar quantization standard for gray-scale fingerprint image compression: Tech. Report LA-UR-931659, Los Alamos Nat'l Lab, Los Alamos. Bracewell, R. N., 1986, The Hartley Transform: Oxford Univ. Press. Chui, C. K., 1992a, An Introduction to Wavelets: (volume 1 of a new series) Academic Press, Boston. Chui, C. K., 1992b, Wavelets: Theory and Applications, Academic Press, Boston. Claerbout, J. F., 1992, Earth Sounding Analysis: Blackwell Scientific Publications, Boston. Clarke, G. K. C , 1968, Time-varying deconvolution filters: Geophysics, 33, 936-944. Cohen, A., Daubechies, I., Jawerth B., and Vial, P., 1993, Multiresolution analysis, 156 Chapter 8. SUMMARY 157 wavelets and fast algorithms on an interval: C. R. Acad. Sci. Paris Ser. I Math, I, 417-421. Combes, J.M., Grossman, A., and Tchamitchian, Ph., (editors), 1990, Wavelets: TimeFrequency Methods and Phase Space, Second Edition, Springer-Verlag, New York. Daubechies, I., Grossmann, A., and Y. Meyer, 1986, Painless Nonorthogonal Expansions: J. Math. Phys., 27, 5, 1271-1283. Daubechies, I., 1988, Orthonormal bases of compactly supported wavelets: Communications in Pure and Applied Mathematics, 41, 909-996. Daubechies, I., 1990, The Wavelet Transform, Time-Frequency Localization and Signal Analysis: IEEE Transactions on Information Theory, 36, 5, 961-1005. Daubechies, I., 1992, Ten lectures on wavelets: Society for Industrial and Applied Mathematics, Philadelphia, PA. David, P. M., and Chapron, B., 1989, Underwater acoustic, wavelets and oceanography: in Wavelets and applications, Meyer, Y. Ed., 1989, Springer-Verlag, New York. Dessing, F. J., and Wapenaar, C. P. A., 1994, Wavefield extrapolation using the wavelet transform: 64th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1355-1358. Dessing, F. J., and Wapenaar, C. P. A., 1995, Efficient migration with one-way operators in the wavelet transform domain: 65th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1240-1243. Donoho, D. L., 1992, De-noising by soft-thresholding, TR 409 November 1992, Submitted for publication: Available at playfair.stanford.edu. Donoho, P. L., Ergas, R. A., and Villasenor, J. D., 1995, High-performance seismic trace compression: 65th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 160-163. Dutilleux, P., Aiguier, C. J., and Joliot, R. F., 1987, An implementation of the 'algorithm 158 Chapter 8. SUMMARY a trous' to compute the wavelet transform: in Combes, J. M . , Grossman, A., and Tchamitchian, Ph., Ed.s, 1990, Wavelets: Time-Frequency Methods and Phase Space, Second Edition, Springer-Verlag, New York. Dziewonski, A., Bloch, S., and Landisman, M . , 1969, A technique for the analysis of transient seismic signals: Bull. Seism. Soc. Am., 59, 1, 427-444. Feder, J., 1988, Fractals, Plenum. Flandrin, P., 1989, On the spectrum of Fractional Brownian Motions: IEEE trans, on information theory, 35, 1, 197-199. Flandrin, P., 1992, Wavelet analysis and synthesis of Fractional Brownian Motion: IEEE trans, on information theory, 38, 2, 910-917. Foster, D. J., Mosher, C. C , and Hassanzadadeh, S., 1994, Wavelet transform methods for geophysical applications: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1465-1468. Gabor, D., 1946, Theory of Communication: Journal of the Institute of Electrical Engineers, London, 93, 429-457. Goupillaud, P., Grossmann, P., and Morlet, J., 1984, Cycle-octave and related transforms in seismic signal analysis: Geoexploration, 23, 1984/85, 85-102. Griffiths, L. J., Smolka, F. R., and Trembly, L. D., 1977, Adaptive deconvolution - a new technique for processing time-varying seismic data: Geophysics, 42, 742-759. Grossmann, A., Kronland-Martinet, R., and Morlet, J., 1987, Reading and understanding continuous wavelet transform: in Combes, J. M., Grossman, A., and Tchamitchian, Ph., Ed.s, 1990, Wavelets: Time-Frequency Methods and Phase Space, Second Edition, Springer-Verlag, New York. Holschneider, M . , Kronland-Martinet, R., Morlet, J., and Tchamitchian, Ph., 1987, A real-time algorithm for signal analysis with help of the wavelet transform: in Chapter 8. 159 SUMMARY Combes, J. M . , Grossman, A., arid Tchamitchian, Ph., Ed.s, 1990, Wavelets: TimeFrequency Methods and Phase Space, Second Edition, Springer-Verlag, New York. Jensen, 0. G., Ulrych, T. J., and Gregotski, M. E., 1995, Deconvolution of seismic sections to recover a 2-dimensional fractal reflectivity field: J. of Seismic Exploration, 4, 4560. Jones, I. F., and Levy, S., 1987, Signal-to-noise ratio enhancement in multichannel seismic data via the Karhunen-Loeve transform: Geophys. Prosp., 35, 12-32. Kolaczyk, E. D., 1994, Wavelet methods for the inversion of certain homogeneous linear operators in the presence of noisy data: Ph.D. thesis, Stanford University. Li, X.-G., Sacchi, M . D., and Ulrych, T. J., 1994, Wavelet transform inversion: application to ID tomography: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 972-975. Li, X.-G., and Ulrych, T. J., 1995, Tomography via Wavelet Transform Constraints: 65th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1070-1073. Li, X.-G., Sacchi, M . D., and Ulrych, T. J., 1996, Wavelet transform inversion with prior scale information: Geophysics, 61, 1379-1385. Li, X.-G., and Ulrych, T. J., 1996, Coherent noise filtering using a 2-D Gabor transform: 66th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1180-1183. Lienard, J. S., and d'Alessandro, 1987, Wavelet analysis of speech: in Combes, J. M . , Grossman, A., and Tchamitchian, Ph., Ed.s, 1990, Wavelets: Time-Frequency Methods and Phase Space, Second Edition, Springer-Verlag, New York. Lines, L. R., and Levin, F. K . Ed.s, 1988 Inversion of Geophysical Data: Geophysics reprint series No 9. Luo, Y., and Schuster, G. T., 1992, Wave packet transform and data compression: 62nd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts. Chapter 8. SUMMARY 160 Mallat, S., 1989a, A theory for multiresolution signal decomposition: the wavelet representation, IEEE Trans. Pattern Anal. Machine Intell., 11,7, 674-693. Mallat, S., 1989b, Multifrequency channel decomposition of images and wavelet models: IEEE Transactions on Acoustics, Speech, and Signal processing, 37 No. 12, 2091. Mallat, S., 1989c, Multiresolution approximation and wavelet orthonormal bases of L2: Transactions of the American Mathematical Society.. Mallat, S., and Hwang, W. -L., 1992, Singularity detection and processing with wavelets: IEEE Transactions on Information Theory, 38, 2, Part II, 617. Mandelbrot, B. B., and Van Ness, J. W., 1968, Fractional Brownian motions, fractional noises and applications, SIAM Rev., 10, 4, 422-437. McMechan, G. A., 1982, Seismic tomography in boleholes: Geophys. J. Royal Astro. Soc, 74, 601-612. Meng, H.-Z., 1995, Geostatistical data integration for reservoir management studies: 65th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1554. Menke, W., 1984, Geophysical data analysis: discrete inverse theory: Academic Press Inc., Orlando. Meyer, Y. Ed., 1989, Wavelets and applications: Springer-Verlag, New York. Meyer, Y., 1993, Wavelets: algorithm and applications, translated and revised by R. D. Ryan: SIAM. Miao, X.-G., and Cheadle, S., 1995, Wavelet transform methods for seismic application, 1995 Canadian Soc. Expl. Geophys. National Convention, Expanded Abstracts, 141-142. Miao, X.-G., and Moon, W. M . , 1994, Application of the wavelet transform in seismic data processing: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1461-1464. Morlet, J., and Arens, G., Fourgeau, E., and Giard, D., 1982, Wave Propogation and 161 Chapter 8. SUMMARY Sampling Theory—Part I: Complex Signal and Scattering in Multilayered Media: Geophysics, 47, 203-221. Morlet, J., and Arens, G., Fourgeau, E., and Giard, D., 1982, Wave Propogation and Sampling Theory—Part II: Sampling Theory and Complex Waves: Geophysics, 47, 222-236. Nolet, G., ed., 1987, Seismic tomography with applications in global seimology and exploration geophysics: D. Reidel Publ. Co.. Nolet, G., 1987, Seismic wave propagation and seismic tomography, in Nolet, G., Ed., Seismic tomography with applications in global seimology and exploration geophysics: D. Reidel Publ. Co., 1-23. Noponen, I., and Keeney, J., 1986, Attenuation of waterborne coherent noise by application of hyperbolic velocity filtering during the tau-P transform: Geophysics, 51, 20-33. Oldenburg, D. W., 1984, An introduction to linear inverse theory: IEEE Transactions on Geosecience and Remote Sensing, GE-22, 6, 665-672. Olmo, G., Lo Presti, L., and Spagnolini, TJ., 1994, Wavelet transform (WT) in geophysical signal processing - an application to velocity analysis: 56th Mtg. Eur. Assoc. Expl Geophys., Extended Abstracts, 94, Session H014. Park, C , and Black, R. A., 1995, Simple time-variant, band-pass filtering by operator scaling: Geophysics, 60, 1527-1535. Parker, R. L., 1994, Geophysical inverse theory: Priceton Univ. Press. Pilkington, M., and Todoeschuck, J. P., 1992, Natural smoothness constraints in crosshole seismic tomography: Geophys. Prosp., 40, 227-242. Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P., 1992, Numerical recipes, Second edition: Cambridge Univ. Press. Rade, L., and Westergren, B., 1995, Mathematrics handook: Birkhauser. Chapter 8. SUMMARY 162 Rosa, A. R., and Ulrych T. J., 1991, Processing via spectral modeling: Geophysics, 56, 1244-1251. Scales, J. A., 1987, Tomographic inversion via the conjugate gradient method: Geophysics, 52, 179-185. Scales, J. A., and Smith M . L., 1994, Introductory Geophysical Inverse Theory: Samizdat Press. Scheuer, T. E., and Oldenburg, D. W., 1988, Aspects of time-variant filtering: Geophysics, 53, 1399-1409. Schuster, G. T., and Sun, Y., 1993, Wavelet filtering of tube and surface waves, 63th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts. Sheriff, R. E., and Geldart, L. P., 1982, History, theory, and data acquisition: Exploration seismology, 1, Cambridge Univ. Press. Shensa, M . J., 1992, Wedding the d trous and Mallat algorithms, IEEE Signal Processing Magazine, 2464-2482. Shensa, M . J., 1993, An inverse DWT for nonorthogonal wavelets: NCCOSC Technical Report, TR 1621. Squires, L. J., and Guillaume, C , 1992, A linear filter approach to designing off-diagonal damping matrices for least-squares inversion problems: Geophysics, 57, 948-951. Stigant, J. P., Ergas, R. A., Donoho, P. L., Minchella, A. S., and Galibert, P. Y . , 1995, Field trial of seismic compression for real-time transmission: 65th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 960-962. Stoksik, M . A., Lane, R. G., and Nguyen, 1994, Accurate synthesis of fractional Brownian motion using wavelets: Electronics Letters, 30, 5, 383-384. Stolt, R. H., 1978, Migration by Fourier transform: Geophysics, 43, 1, 23-48. Stork, C , 1991, High resolution SVD of the coupled velocity determination and reflector imaging problem: 61st Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Chapter 8. SUMMARY 163 Abstracts, 981-985. Strang, G., 1986, Introduction to Applied Mathematics: Wellesley-Cambridge Press. Strang, G., 1989, Wavelets and dilation equations: a brief introduction, SIAM Review, 31, 614-627. Strang, G., and Nguyen, T., 1996, Wavelets and Filter Banks: Wellesley-Cambridge Press, Wellesley. Taner, M . T., Koehler, F., and Sheriff, R. E., 1979, Complex seismic trace analysis: Geophysics, 44, 1041-1063. Tarantola, A., 1987, Inverse problem theory: Elsevier, New York. Todoeschuck, J. P., Jensen, O. G., 1988, Joseph geology and seismic deconvolution: Geophyics, 53, 1410-1414. Turcotte, D. L., 1992, Fractal and chaos in geology and geophysics: Cambridge university press. Wang, B., and Lawrence, W., 1994, Stochastic view of damping and smoothing, and applications to seismic tomography: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1343-1346. Wang, L., Chen, C.-T., and Lin, W . - C , 1994, An efficient algorithm to compute the complete set of discrete Gabor coefficients: in IEEE Trans, on Image Processing, 3, 1, 87-92. Wickerhauser, M . V., 1991, Lectures on wavelet packet algorithms: ftp site pascal.math.yale.edu. Womack, J. E., and Cruz, J. R., 1994, Seismic data filtering using a Gabor representation: IEEE Trans, on Geoscience and Remote Sensing, 32, 2, 467-472. Wornell, G. W., 1990, A Karhunen-Loeve-like expansion for 1/f processes via wavelets: IEEE Trans, on Information Theory, 36, 4, 859-861. Wornell, G. W., and Oppenheim, A. V., 1992, Estimation of fractal signals from noisy measurements using wavelets: IEEE Trans. Signal Proc, 40, 3, 611-623. Chapter 8. SUMMARY 164 Wu, Y . , and McMechan, G. A., 1995, Wavefield extrapolation in the wavelet domain with application to poststack migration: Geophys., Expanded Abstracts, 1236-1239. 65th Ann. Internat. Mtg., Soc. Expl. Appendix A The Morlet W T of 6(t;t ) Q From equations (4.1), (4.2) and (4.3), the Morlet WT coefficients are (A.l) S(a,b) = The inverse WT can be written as x * * ; * o = 1 f f 7T- / / 1, U^, J J a ,to — b. , ,t — b. da „ ~ W )-db, a a a* , A . (A.2) which can be further expanded to £(*;*„) = TT I I "T e ^ e _ ( " f c ) e^^dadb. (A.3) (7,/, J J a 3 Integrating over b we obtain (A.4) where the integral can be approximated by a summation in the discrete case. 165
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Application of wavelet transforms to seismic data processing...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Application of wavelet transforms to seismic data processing and inversion Li, Xin-Gong 1997
pdf
Page Metadata
Item Metadata
Title | Application of wavelet transforms to seismic data processing and inversion |
Creator |
Li, Xin-Gong |
Date Issued | 1997 |
Description | The goal of this thesis is to use wavelet transforms as tools to analyze simultaneously the time (or location) and frequency (or scale) variant characteristics of seismic data and then to apply these information to two important geophysical problems: noise filtering and traveltime inversion. In the noise filtering problem, the discrete non-orthogonal wavelet transform is used to analyze the time-variant characteristics of signal and noise. This information can be used to filter noise in the transformed domain. This approach is applied to ground roll attenuation for field seismic data according to the localized characteristics. Geophysical inverse problems are, in general, non-linear, underdetermined and ill-posed. An initial model and some prior information describing the model is needed to find the model perturbation that minimizes an objective function. Traditionally, the perturbation is assumed to be an independent Gaussian or a correlated but stationary process, an assumption which is not physically correct according to the analyses of well logs. Indeed, analysis of logs described in this thesis suggests that log data can be decomposed into (a) a fractal process, usually non-stationary, with wavelet coefficients which follow a power law, and (b) some non-fractal structures including a large scale trend and some spiky details. The inverse problem can be solved using the wavelet transform in two steps. First, I solve an overdetermined problem to invert the large scale trend, then, I solve for a model with fractal constraints. Another important application of the wavelet transform in inverse problems is to reduce the effects of the noise in the input data. Since the noise effect varies with scale, the regularization can be done accordingly. Wavelet transform constrained inversions are tested using 1-D and 2-D synthetic examples. |
Extent | 18242652 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-04-02 |
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.0088134 |
URI | http://hdl.handle.net/2429/6734 |
Degree |
Doctor of Philosophy - PhD |
Program |
Environmental Sciences |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1997-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1997-250946.pdf [ 17.4MB ]
- Metadata
- JSON: 831-1.0088134.json
- JSON-LD: 831-1.0088134-ld.json
- RDF/XML (Pretty): 831-1.0088134-rdf.xml
- RDF/JSON: 831-1.0088134-rdf.json
- Turtle: 831-1.0088134-turtle.txt
- N-Triples: 831-1.0088134-rdf-ntriples.txt
- Original Record: 831-1.0088134-source.json
- Full Text
- 831-1.0088134-fulltext.txt
- Citation
- 831-1.0088134.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0088134/manifest