Seismic Absorption Estimation and Compensation by Changjun Zhang M.Sc., The University of British Columbia, 2001 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Geophysics) The University of British Columbia (Vancouver) November, 2008 c Changjun Zhang 2008 Abstract As seismic waves travel through the earth, the visco-elasticity of the earth’s medium will cause energy dissipation and waveform distortion. This phenomenon is referred to as seismic absorption or attenuation. In hydrocarbon reservoir description, seismic absorption can be used as an important attribute to interpret for fluid units. In seismic data processing, the information about seismic absorption can be used to enhance seismic data resolution by means of absorption compensation. The absorptive property of a medium can be described by a quality factor Q, which determines the energy decay and a velocity dispersion relationship. The quality factor and the velocity govern the propagation of seismic energy in the earth. In this thesis, four new ideas have been developed to deal with the estimation and application of seismic absorption. These ideas are detailed in four chapters respectively. By assuming that the amplitude spectrum of a seismic wavelet may be modeled by that of a Ricker wavelet, an analytical relation has been derived to estimate a quality factor from the seismic data peak frequency variation with time. This relation plays a central role in quality factor estimation problems. Quality factors with coarse resolution can be estimated from prestack common midpoint (CMP) gathers, or a post-stack single trace based on event picking. To estimate interval Q for reservoir description, a method called reflectivityguided seismic attenuation analysis is proposed. This method first estimates peak frequencies at a common midpoint location, then correlates the peak frequency with sparsely-distributed reflectivities, and finally calculates Q values from the peak frequencies at the reflectivity locations using an analytical expression. The peak frequency is estimated from the prestack CMP gather using peak frequency variation with offset (PFVO) analysis which is similar to amplitude variation with offset (AVO) analysis in implementation. The estimated Q section has the same layer boundaries of the acoustic impedance or other layer properties. Therefore, the seismic attenuation property obtained with the guide of reflectivity is easy to interpret for the purpose of reservoir description. To overcome the instability problem of conventional inverse Q filterii Abstract ing, Q compensation is formulated as a least-squares (LS) inverse problem based on statistical theory. The matrix of forward modeling is composed of time-variant wavelets. The LS de-absorption is solved by an iterative non-parametric approach. To compensate for absorption in migrated seismic sections, a refocusing technique is developed using non-stationary multi-dimensional deconvolution. A numerical method is introduced to calculate the blurring function in layered media, and a least squares inverse scheme is used to remove the blurring effect in order to refocus the migrated image. This refocusing process can be used as an alternative to regular migration with absorption compensation. Theoretical derivations and numerical tests showing the validity of the new methods are detailed in the chapters. iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Abstract List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgments Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . xi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . 1.1 Thesis Objectives . . . . . . . . 1.2 Seismic Absorption Phenomena 1.3 Seismic Absorption Research . . 1.4 Basic Assumptions . . . . . . . . 1.4.1 Effective Absorption . . . 1.4.2 Frequency Independent Q 1.5 Thesis Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Mechanism and Description of Seismic Absorption . 2.1 Physical Hypothesis . . . . . . . . . . . . . . . . . . . 2.1.1 Scattering Attenuation and Intrinsic Attenuation 2.1.2 Causes of Intrinsic Attenuation . . . . . . . . . 2.2 Mathematical Description . . . . . . . . . . . . . . . . . 2.2.1 The Voigt Model . . . . . . . . . . . . . . . . . 2.2.2 The Generalized Linear Solid Model . . . . . . . 2.2.3 Describing Absorption by Quality Factors . . . . 2.2.4 Choice and Discussion . . . . . . . . . . . . . . 2.2.5 Absorption and Seismic Wave Extrapolation . . 2.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 1 1 3 3 4 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 6 7 7 8 9 11 12 14 16 19 iv Table of Contents 3 Quality Factor Estimation from Seismic Peak Frequency Variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Absorption and Peak Frequency Variation . . . . . . . . . . 3.2.1 Spectral Ratio Method . . . . . . . . . . . . . . . . . 3.2.2 Q and Peak Frequency Variation . . . . . . . . . . . . 3.3 Q Estimation for Seismic Absorption Compensation and Migration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Estimation of Q from CMP Gathers . . . . . . . . . . 3.3.2 Q Estimation from a Poststack Trace . . . . . . . . . 3.4 Summary and Discussions . . . . . . . . . . . . . . . . . . . 4 Seismic Absorption Analysis for Reservoir Description . . 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Review of the Methods for Seismic Attenuation Analysis . . 4.2.1 Iso-frequency Analysis . . . . . . . . . . . . . . . . . 4.2.2 Direct Q Estimation . . . . . . . . . . . . . . . . . . . 4.2.3 Attenuation Tomography . . . . . . . . . . . . . . . . 4.3 Reflectivity-guided Q Analysis, RGQA . . . . . . . . . . . . 4.3.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Peak Frequency Variation with Offset (PFVO) Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.3 Implementation Steps . . . . . . . . . . . . . . . . . . 4.4 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . 4.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . 5 Seismic Absorption Compensation: a Least Squares Inverse Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Issues of Absorption Compensation . . . . . . . . . . . . . . 5.2.1 Method Choice . . . . . . . . . . . . . . . . . . . . . 5.2.2 Inaccurate Q Estimation . . . . . . . . . . . . . . . . 5.2.3 Random Noise . . . . . . . . . . . . . . . . . . . . . . 5.3 Time-variant Deconvolution as an Inverse Problem . . . . . . 5.3.1 Noise and Deconvolution . . . . . . . . . . . . . . . . 5.3.2 Bayes Probabilistic Inference . . . . . . . . . . . . . . 5.3.3 Absorption Compensation and Inversion . . . . . . . 5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 21 22 23 25 27 28 33 35 38 38 39 39 40 41 42 42 44 46 47 51 52 52 54 54 54 55 58 58 59 61 68 v Table of Contents 6 Refocusing Migrated Images in Absorptive Media . . . . . 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Blurring Function caused by Seismic Absorption . . . . . . . 6.2.1 Analytical Analysis . . . . . . . . . . . . . . . . . . . 6.2.2 Numerical Computation . . . . . . . . . . . . . . . . 6.3 Removing Absorption Blurring Function using a Least Squares Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Data Experiments . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Conclusion and Discussion . . . . . . . . . . . . . . . . . . . 69 69 71 72 72 . . . . . . . . . . . . . . . . . . . 86 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 7 Discussion and Conclusions 75 81 85 vi List of Tables 3.1 A pseudo code of prestack Q estimation. . . . . . . . . . . . . 31 6.1 A pseudo code of the refocusing algorithm. . . . . . . . . . . 79 vii List of Figures 2.1 2.2 2.3 2.4 2.5 2.6 3.1 3.2 3.3 3.4 Attenuation is determined by composite rock physical properties, such as lithology, porosity, saturation, pore fluid, etc.. Velocity dispersion comparison with Q = 30. The red curve is for the linear Q model and the green one is for Kjartansson’s model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Amplitude attenuation comparison with Q = 30. The red curve is for the linear Q model and the green one is for Kjartansson’s model. . . . . . . . . . . . . . . . . . . . . . . . . . Seismic wave behavior in absorptive media defined by v, ρ and Q. (a) Without absorption. (b) With absorption. . . . . Evolution of a unit impulse with time in a 1-D absorptive medium (Q = 30). . . . . . . . . . . . . . . . . . . . . . . . . A Green’s function in an absorptive 2-D medium. . . . . . . . (a) A reflection in a CMP gather. The medium is absorptive with Q = 10. Traces have been normalized based on their maximum amplitudes. (b) Amplitude spectra of the original source signature (black), trace 1 (red), trace 11 (green), and trace 21 (gray). . . . . . . . . . . . . . . . . . . . . . . . . . . (a) A synthetic CMP gather of two reflections with absorption and 10% random noise. Traces have been normalized based on their maximum amplitude. (b) Values of input quality factors (Q1 = 10 and Q2 = 20 ). (c) Computed quality factors (Q1 = 10.04 and Q2 = 20.12). . . . . . . . . . . . . . . Variation of peak frequencies of the two reflections in Figure 3.2. Peak frequencies of the first event and the second event are shown in black and red respectively. . . . . . . . . . Q estimation. (a) is a real seismic trace. (b) is its windowed time-variant spectrum. The picked peak frequency points are marked by circles and connected by straight lines. (c) is the inverse of estimated quality factors. . . . . . . . . . . . . . . . 9 15 16 17 18 19 24 32 33 36 viii List of Figures 4.1 4.2 4.3 4.4 4.5 4.6 4.7 5.1 5.2 5.3 5.4 5.5 5.6 Frequency decay caused by absorption. (a) is an attenuation model, (b) is the synthetic trace and (c) is the spectrum of frequency decomposition. . . . . . . . . . . . . . . . . . . . . Frequency decay caused by thin-bed tuning and absorption. (a) shows a reflectivity function with closely-spaced coefficients. (b) is the corresponding absorption property. (c) is the synthetic trace and (d) is its spectrum of time-frequency decomposition. . . . . . . . . . . . . . . . . . . . . . . . . . . Flow diagram of reflectivity-guided Q analysis. . . . . . . . . Peak frequencies at layer boundaries. (a) is the peak frequencies from the trace in Figure 4.1 (b). (b) is the extracted inverse Q curve. . . . . . . . . . . . . . . . . . . . . . . . . . . Peak frequency variation with offset. (a) is a CMP gather built from the model shown in Figure 4.2. (b) is the corresponding peak frequencies at zero offset. (c) is the relative gradient of peak frequencies along offset, and (d) is the estimated Q curve. . . . . . . . . . . . . . . . . . . . . . . . . . . A real data example. (a) is a CMP gather. (b) is the corresponding non-stationary distributed peak frequencies. (c) is the peak frequency curve at zero offset. (d) is the relative gradient curve of peak frequencies. (e) is the reflectivity function and (f) is the estimated Q curve. . . . . . . . . . . . Real data experiment. (a) is the peak frequency section and (b) is the attenuation section (Q−1 ). . . . . . . . . . . . . . . Experiments on amplitude compensation using different Q values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A synthetic trace with three spikes, Q = 30, and its phase correction using different Q values. . . . . . . . . . . . . . . . The amplitude of an ideal inverse Q filter (black) and its approximation (red). The two curves are overlapping before the fork. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Convergence process of the inverse scheme. (e) is a synthetic trace. (a-d) are the inverse results after iteration 1, 2, 3 and 4 respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) is a synthetic trace with 20% noise. (b) is the result of absorption compensation. The rightmost trace is the original reflectivity series plotted as a reference. . . . . . . . . . . . . Q compensation for a real seismic section. (a) shows a real seismic section. (b) is the result of absorption compensation. 40 41 43 44 48 49 50 56 57 58 65 66 67 ix List of Figures 6.1 6.2 6.3 6.4 6.5 6.6 Diagram of the energy distribution of a migration blurring function. Symbol “V” represents a receiver location. . . . . . Examples of the blurring functions in an absorptive medium. (a) at time t = 400 ms. (b) at time t = 700 ms. . . . . . . . . Amplitude of a kernel matrix G. . . . . . . . . . . . . . . . . A synthetic example. (a) shows an absorption-blurred sinefunction shaped structure. (b) is the refocusing result. . . . . Input information. (a) is the velocity section and (b) is the inverse Q section. . . . . . . . . . . . . . . . . . . . . . . . . . Migration results. (a) Time-variant spectral whitening plus phase-shift migration. (b) Phase-shift migration plus refocusing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 76 80 82 83 84 x Acknowledgments I would like to express my deepest thanks to my supervisor Dr. Tadeusz Ulrych for his moral support and scientific enthusiasm on my research. His invaluable academic advice contributes to every aspect of this thesis. My sincere thanks go to Dr. Michael Bostock for his willingness to cosupervise my Ph.D. program. The course of “Theory of Seismology” I took from him is very useful in my academic pursuit. My thanks also go to Dr. Andrew Calvert, Dr. Elizabeth Hearn and Dr. Matthew Yedlin for their participation as committee members. Discussions with and questions from my supervisors and other committee members helped to shape my research on seismic absorption. My special thanks go to Dr. Bill Nickerson. He kindly proofread my manuscripts and helped me correct errors in both English and rock physics. Part of the research included in this thesis was supported by a UBC international scholarship and grants from the companies that supported the Consortium for the Development of Specialized Seismic Techniques, CDSST, which was led by Dr. Tadeusz Ulrych. xi Dedication To the memory of my father Yaoqiang Zhang. xii Chapter 1 Introduction 1.1 Thesis Objectives The research being presented in this thesis deals with signal absorption phenomena observed in seismic exploration. The topics addressed are the estimation of the absorption properties and the compensation of the attenuated signals. 1.2 Seismic Absorption Phenomena Conventionally, in seismic exploration, the earth is modeled as an ideal elastic medium, and seismic wave propagation is explained by means of the elastic (or acoustic) wave equation. In practice, the propagation of seismic waves in the earth is in many respects different from seismic wave propagation in an ideal solid. For example, the earth material is anisotropic, heterogeneous, porous, etc.. The traditional elastic wave equation is not accurate enough to describe the wave behavior for this more complicated medium. Generally, the visco-elasticity of the earth materials causes seismic energy dissipation, and thus decreases the amplitude and modifies the frequency content of propagating waves. This phenomenon of wave energy dissipation is called seismic absorption or seismic attenuation, although attenuation is just one aspect of the wave behavior in visco-elastic media. From the studies of absorptive phenomena observed in seismic data, we may be able to construct images with better resolution in seismic data processing and extract more detailed information about the rock materials in seismic data inversion. 1.3 Seismic Absorption Research Absorption researchers have made contributions in several related areas: from explaining the absorption mechanism, estimating absorptive properties 1 Chapter 1. Introduction and modeling its effect numerically and physically, to removing its effect, and utilizing absorptive properties to locate and describe hydrocarbon reservoirs. In recent years, due to the requirement of obtaining seismic images with higher resolution and higher accuracy, as well as of extracting more information for reservoir description in seismic exploration, seismic absorption has become a field of extensive investigation, similar to other fields that take into account secondary effects in seismic wave propagation, such as anisotropy, angle-dependent reflectivity and wave-mode conversion. Mathematical models, resulting from modification to the elastic wave equation, like the Voigt model, the standard linear solid model, the Futterman model and others have been used for different applications for a long time (Ricker, 1953; Biot, 1956a,b; Futterman, 1962; White, 1975; Toksoz et al., 1979a,b; Aki and Richards, 1980; Johnson, 2001). More complicated seismic attenuation models have also been suggested, such as those by Dvorkin et al. (1995), Pride et al. (2004), and Carcione and Picotti (2006). These models relate seismic attenuation to the porous properties of rock materials, such as porosity, saturation, fluid content, permeability, etc., and try to fit field data and laboratory observations more accurately than allowed by means of traditional methods. An overview of seismic attenuation models can be found in Toverud and Ursin (2005). Numerical simulation of seismic waves in absorptive media is generally implemented through finite-difference modeling. This can either be implemented in the frequency domain or time domain (Day and Minster, 1984; Carcione, 1993; Blanch et al., 1995; Bohlen, 2002). Physical models have also been built to measure seismic absorption in more realistic environments (Toksoz et al., 1979a; Winkler et al., 1979; Sams et al., 1997; Batzle et al., 2005). Both mathematical modeling and physical modeling are important tools to confirm the accuracy of theories and discriminate mechanisms regarding seismic absorption. The seismic absorption property of a medium is usually described by a quality factor designated by Q. Many algorithms to estimate Q values have been published. Among them, the most popular is the spectral ratio method (Spencer et al., 1982; Tonn, 1991). Methods suggested by Quan and Harris (1997), Zhang and Ulrych (2002), Taner and Treitel (2003), Singleton et al. (2006), Rickett (2006) and Guerra and Leaney (2006) are claimed to be more robust than the traditional spectral ratio methods, and/or disclose more detailed attenuation information of rock properties. The objectives of Q estimation can be grouped into two categories, one is for Q compensation, and another one is for reservoir description. For Q compensation, generally, we only require absorptive information for major layers. For reservoir de2 Chapter 1. Introduction scription, we require Q information of higher resolution for the study of fine layers. These two different requirements are similar to the resolution requirements for velocities in seismic data processing and inversion. To remove the effect of absorption, traditionally we can use time-variant spectral whitening or time variant deconvolution (Robinson, 1979; Hargreaves and Calvert, 1991; Bickel, 1993; Varela et al., 1993; Yilmaz, 2000). Computation techniques trying to stabilize inverse Q filtering have been introduced by Margrave et al. (2003), Wang (2003), and Zhang and Ulrych (2007). The inclusion of absorption effects into the migration scheme has also been suggested by Berkhout (1981), Dai and West (1994), Yu et al. (2002), Cui and He (2004) and Wang (2008) for stacked seismic sections and Mittet et al. (1995) for prestack seismic data. Migration plus Q compensation can be considered as multi-dimensional seismic absorption compensation. Many researchers have devoted their efforts to relating seismic attenuation to reservoir description, due to the fact that for porous rocks saturated with fluid, generally, strong absorption can be observed. Parra and Hacket (2002), Castagna et al. (2003), Taner and Treitel (2003), H¨ ubert et al. (2005), Chapman et al. (2005) and Dvorkin and Mavko (2006) use different physical hypotheses to relate attenuation properties to reservoir characterization. There is more and more interest in seismic exploration in using attenuation properties as an attribute in reservoir description, because of reported successful cases. In a medium, shear waves and compressional waves travel at different speeds. They could experience different absorption on their wave-paths depending on the fluid content in rock fractures. Winkler et al. (1979), Klimentos (1995), Deffenbaugh et al. (2000), Bale and Stewart (2002), and Mavko et al. (2005b) have investigated the characteristics of shear wave attenuation in seismic data. 1.4 Basic Assumptions 1.4.1 Effective Absorption Seismic attenuation is usually categorized into scattering attenuation and intrinsic attenuation. Scattering attenuation is generally caused by threedimensional heterogeneities in the subsurface that distribute wave energy in arbitrary directions. Often quoted research about scattering attenuation includes O’Doherty and Anstey (1971), Richards and Menke (1983) and Sato and Fehler (1998). 3 Chapter 1. Introduction Intrinsic attenuation is caused by internal friction among grains in the rock matrix, and the relative movement between the solid rock matrix and the pore fluid, when a seismic wave travels through loosely consolidated or porous media. Internal friction and fluid flow are considered directly related to rock porosity, permeability and the fluid content in the void space. More detailed explanation of the physics of intrinsic attenuation can be found in Biot (1956a), Biot (1956b),White (1975), Ricker (1977), Dvorkin et al. (1995), Johnson (2001) and Pride (2003), etc.. In this thesis, if not stated otherwise, the absorption dealt with is effective absorption, which is the composite effect of scattering and intrinsic attenuation. This is based on the premise that scattering attenuation and intrinsic attenuation have similar features in observed seismic or well log data, and also that, at present, there is no reliable way to separate them. 1.4.2 Frequency Independent Q When using Q to describe seismic attenuation, it is often assumed that Q does not change with frequency in the seismic frequency range (10 − −200Hz). However, there are many mathematical models which describe the change of Q with frequency (Ricker, 1953; Pride et al., 2004; Carcione and Picotti, 2006). Frequency dependent Q models may be theoretically more attractive in visco-elastic wave equation derivation, and/or better fitting laboratory testing data. The usability of these more complicated seismic attenuation models is hindered by their dependence on some unmeasurable parameters and the necessity to use frequency dependent Q values in practice (Dvorkin and Mavko, 2006). When Q is considered as a function of frequency, it would be convenient if the Q function can be possibly separated into a frequency dependent part and a frequency independent part, so that we can use the frequency independent part to describe attention property. In this thesis, we use the simplest convention: absorption property Q does not depend on frequency, while attenuation changes with frequency. The assumption is valid over the signal bandwidth in seismic exploration and conforms to many theoretical models and laboratory observations. 1.5 Thesis Structure This chapter introduces the background of the research on seismic absorption. 4 Chapter 1. Introduction Chapter 2 reviews physical mechanisms and mathematical descriptions of seismic absorption. A definition of quality factor Q and the corresponding dispersion relation between attenuation and velocity are analyzed. The linear dispersion model is adopted in this thesis because it is appropriate not only for handling the problem of estimating quality factors from both prestack and poststack seismic data, but also for absorption compensation. Chapter 3 develops a method to estimate quality factors from wavelet peak frequency variation. The analytical formula derived can be used to estimate Q from a CMP gather, and also from a stacked seismic trace. The estimated Q provides the information required for one-dimensional and multi-dimensional Q compensation. Chapter 4 suggests a scheme called reflectivity-guided Q analysis for reservoir description. It first estimates peak frequencies at a CMP location, then correlates the peak frequency with sparsely-distributed reflectivities, and, finally calculates a Q curve from the peak frequencies at the locations of impedance contrasts using an analytical equation. Chapter 5 formulates seismic absorption compensation as a least-squares (LS) inverse problem based on statistical theory. The purpose of doing so is to overcome the instability problem of regular inverse Q filtering. The kernel matrix of the inversion is composed of time-variant wavelets. The LS de-absorption can be solved iteratively with fast convergence. The results of de-absorption are related to the accuracy of the estimated Q values and also of the seismic wavelets. Chapter 6 investigates the blurring effect in migrated images when using a regular migration algorithm to migrate seismic data with absorption. The deblurring problem is considered as a multi-dimensional time-variant deconvolution and is solved using a LS inverse scheme in time-wavenumber domain to refocus migrated images. The refocusing processing is an alternative to traditional post-stack migration with absorption compensation. Chapter 7 draws conclusions from the research, proposes future research topics, and discusses some other issues in seismic absorption which are not detailed in the previous chapters. 5 Chapter 2 Mechanism and Description of Seismic Absorption Seismic waves traveling in the earth experience amplitude attenuation and velocity dispersion due to the absorptive property of the rocks. Conventionally, in seismic exploration, the earth is modeled as an ideal elastic medium, that is, each layer of the earth can be fully described by two parameters: velocity and density (v, ρ). In practice, the propagation of seismic waves in the earth is in many respects different from seismic wave propagation in an ideal solid. When a seismic wave experiences absorption during its propagating process, the wave behavior can be described by velocity, density and a quality factor (v, ρ, Q), if the absorptive property of a medium is fully determined by a quality factor (Q). There are different physical hypotheses to explain signal attenuation observed in seismic data and there are also numerous mathematical formulas available to describe the attenuation phenomenon. This chapter investigates the causes of the absorption phenomenon and how absorption is described in terms of wave theory by reviewing the appropriate works in the published literature. 2.1 Physical Hypothesis A great deal of attention has been given to the study of the attenuation of seismic waves (Ricker, 1977; White, 1983; Mavko and Nur, 1979; Toksoz et al., 1979b; Kneib and Shapiro, 1995; Carcione, 1995; Dvorkin et al., 1995; Johnson, 2001; Pride et al., 2004). Loss parameters for rocks have been measured by many techniques over a wide range of frequency and environmental conditions. Physical mechanisms for energy loss have been proposed and described mathematically in varying degrees of detail. Absorption is a composite effect caused by the earth media which have a wide range of properties, such as lithology, inhomogeneity, porosity, permeability, saturation, fluid contents, etc.. There is no general consensus as to 6 Chapter 2. Mechanism and Description of Seismic Absorption what the dominant loss mechanism is. In different environments, a different mechanism may play a dominant role. 2.1.1 Scattering Attenuation and Intrinsic Attenuation There are two mechanisms that are usually considered to cause seismic attenuation, scattering and intrinsic attenuation. Scattering redistributes wave energy within the medium but does not remove energy from the overall wavefield. Conversely, intrinsic attenuation refers to various mechanisms that convert vibrational energy into heat through friction, fluid flow, and thermal relaxation processes. The effect of scattering and intrinsic attenuation on seismic data are the same, in that they both attenuate seismic amplitude and distort the seismic waveform. Scattering attenuation is generally caused by relatively simple factors, such as three-dimensional random heterogeneity in rocks, which distributes wave energy in arbitrary directions, and interfaces which reflect and redirect wave propagation. (O’Doherty and Anstey, 1971; Richards and Menke, 1983; Sato and Fehler, 1998). Intrinsic attenuation is caused by more complicated rock properties. Observed seismic attenuation is generally the integrated effects of both scattering and intrinsic attenuation. For porous rocks saturated with fluid, intrinsic absorption generally plays the dominant role (Sams et al., 1997). 2.1.2 Causes of Intrinsic Attenuation As has been stated, seismic absorption is a comprehensive effect caused by many factors. For intrinsic attenuation, two factors are generally considered the most important. One is the internal friction resulting from the relative slide along the contacts of grains in the rock matrix during wave propagation. Another is the fluid flow inside rock pore spaces. Part of the traveling energy from a seismic source is transformed into relative movement and dissipates inside and among rock grains due to friction and fluid flow. Some scientists (Ricker, 1977; Aki and Richards, 1980) considered thermoelasticity as the most viable model to explain intrinsic attenuation at lithospheric temperatures, because the required scales for rock grains and cracks along with the amount of attenuation caused by thermo elasticity are in closest agreement with observations. Rocks have microscopic cracks and pores which may contain fluids. These cracks can have a profound influence on the propagation velocities of P and S-waves. Biot (1956a) and Biot (1956b) analyzed wave propagation in isotropic porous solids where the coupling of motion between the fluid and 7 Chapter 2. Mechanism and Description of Seismic Absorption the solid matrix was considered. He arrived at expressions for attenuation due to fluid flows within non-connecting pores initiated by elastic waves. White (1983) (p.131) studied Biot’s model and concluded that the attenuation predicted by Biot’s model is extremely small for frequencies less than 100Hz. Pride et al. (2004) introduced a model which explains seismic attenuation in porous rocks as resulting from wave-induced fluid flow. Their model also confirms that “squirt-flow” in Biot’s model is incapable of explaining the measured level of attenuation, and describes that “mesoscopic” scale (larger than the grain sizes but smaller than wavelengths) factors are the main contributors of observed seismic attenuation. These factors are lithological variation and fluid-saturated pore spaces in a medium. Research carried out by Carcione and Picotti (2006) shows that the most effective loss mechanisms result from porosity variation and partial saturation in rocks. Rock grain and frame-moduli variations are the next cause of attenuation. Many questions remain unanswered. A major complication is the wide range of properties possessed by earth materials. Conclusions drawn from data of oil-reservoir rocks may not be applicable to materials in the upper mantle. Hence, it is difficult to make a choice among the proposed mechanisms even if this were desirable. In fact, no one mechanism can be expected to describe losses in all rocks under all environmental conditions. Physical mechanism studies help us to understand how rock absorption properties are related to rock physical properties by means of suitable rock physics modelling. Although several such models have been introduced in the cited references, there is no consensus on which model is best. Figure 2.1 illustrates a generic idea that seismic attenuation is the composite effect of rock physical parameters, and that the attenuation property can be denoted by Q. A rock body is very complicated. To describe it needs numerous parameters, such as lithology, porosity, saturation, pore fluid, permeability, etc.. However, the wave behavior inside this rock body can be simply described by density, velocity, and Q. Rock physical models try to link the wave behavior defined by ρ, v, and Q to rock physical properties. Seismic inversion tries to deduce rock properties from seismic attributes. Q is an important attribute for rock property inversion. 2.2 Mathematical Description In dynamic mechanics, the behavior of a wave traveling in a medium is derived by Hooke’s law which relates stress and strain. Generally speaking, there are two kinds of attenuation models built from the point of view of 8 Chapter 2. Mechanism and Description of Seismic Absorption Figure 2.1: Attenuation is determined by composite rock physical properties, such as lithology, porosity, saturation, pore fluid, etc.. dynamic mechanics. One is the Voigt model and the other is the generalized linear solid model (Du, 1996). Mathematical models, derived directly from observed amplitude decay and corresponding attenuation-dispersion relations, describe seismic absorption, generally, by means of a Q parameter (Aki and Richards, 1980; Kjartansson, 1979). 2.2.1 The Voigt Model The Voigt model, also known as the Kelvin-Voigt model, is a widely investigated method of introducing losses that has the advantage of yielding a linear wave equation which can be solved for arbitrary time dependence. The assumption is made that stresses are directly proportional to strain rates, as well as to the components of strain themselves. This assumption was introduced independently by Stokes, Kelvin and Voigt (Ricker, 1977), and its implications have been investigated by Ricker (1977) and White (1983). This kind of medium is commonly called a Voigt solid and the resulting equation is called Stokes equation (Ricker, 1977). In elastic mechanics, using the Einstein notation, Hooke’s law is expressed as σij = λθδij + 2µǫij , (2.1) 9 Chapter 2. Mechanism and Description of Seismic Absorption where σij is the stress tensor, ǫij is the strain tensor, and δij is the Kronecker delta function. λ and µ are the two Lame’s parameters. θ is the volume strain or dilatation. If there is no body force, using Newton’s Law σij,j = ρ ∂ 2 ui . ∂t2 The constitutive relation leads to the equation of motion in term of particle displacement ∂2u (2.2) (λ + µ)(∇ · θ) + µ∇2 u = ρ 2 , ∂t 2 where u is a vector of displacement, ∂∂t2u is its second derivative with respect to time and ρ is density. For a plane wave, if we assume the displacement is parallel to the depth z-axis (ux and uy are both zero) and that uz is independent of x and y, equation (2.2) reduces to (λ + 2µ) ∂ 2 uz ∂ 2 uz = ρ , ∂z 2 ∂t2 (2.3) for one-dimensional wave propagation. In absorptive media, the constitutive relation must incorporate the fact that stress is not only proportional to strain, it is also proportional to the time derivative of strain. The stress-strain behavior in absorptive media is described by a modified Hooke’s law which includes strain-rate terms: σij = λ + λ′ ∂ ∂t θδij + 2 µ + µ′ ∂ ∂t ǫij , where λ′ and µ′ are perturbations of Lame’s parameters, λ and µ, due to absorption. Equation (2.3) in an absorptive medium becomes (λ + 2µ) 3 ∂ 2 uz ∂ 2 uz ′ ′ ∂ uz + (λ + 2µ ) = ρ . ∂z 2 ∂t∂z 2 ∂t2 (2.4) This equation is a Stokes equation without body force. If we replace Lame’s parameters by Young’s modulus, M = λ + 2µ and M ′ = λ′ + 2µ′ , and transform equation(2.4) into the frequency domain, we have (M + iωM ′ ) ∂ 2 Uz = −ρω 2 Uz , ∂z 2 (2.5) where Uz = U (z, ω) designates the space-dependent wavefield at an angular frequency ω. A wavefield satisfying equation (2.5) has the following form U (z, ω) = U0 e±G(ω)z , 10 Chapter 2. Mechanism and Description of Seismic Absorption and G(ω) = − ρω 2 M + iωM ′ 1/2 . (2.6) This complex propagation parameter G(ω) may be expressed in terms of attenuation ap and phase velocity vp as G(ω) = ap + iω/vp . (2.7) By comparing equation (2.6) and equation (2.7) and defining ω0 = M/M ′ , it can be shown that when ω 2 < ω02 , attenuation increases as the square of the frequency, and velocity is approximately constant (Ricker, 1977). These observations are expressed as follows ap = ω2 , 2 M/ρ ω0 1 and vp = M/ρ. The expression for ap means that attenuation increases with the square of frequency. However, many measurements have indicated a first power dependence with frequency. In addition, as pointed out by White (1983)(p.87), the Voigt model also violates causality. Regarding causality, however, there are different observations from that of White (1983). Duren and Heestand (1995) and Buckingham (2005) have derived independently causal solutions to the Stokes equation (equation (2.4)). 2.2.2 The Generalized Linear Solid Model The constitutive relation between stress σ and strain ǫ for visco-elasticity can be formulated simply in the frequency domain as (Malvern, 1969; Carcione et al., 1988) σ(ω) = M (ω)ǫ(ω). Here M (ω) is the complex, frequency dependent, visco-elastic modulus. A generalized linear solid body can be described by N M (z, ω) = Mu (z, ω) + j=1 aj ∆Mj (z) . iω/ωj (z) − 1 (2.8) For one-dimensional P wave, in this equation, Mu (z, ω) is the un-relaxed Young’s modulus at a depth z, N is the number of absorptive mechanisms. 11 Chapter 2. Mechanism and Description of Seismic Absorption aj is a constant coefficient, ωj is the jth relaxation frequency and ∆Mj (z) is the difference between the un-relaxed Young’s modulus and the contribution to the Young’s modulus due to the jth relaxation mechanism. These parameters can be adjusted to fit any realistic linear absorption law, according to Carcione et al. (1988). The two models described in this section and the previous section are constructed by assuming that anelasticity of rocks will result in a constitutive relation that is different from that of elastic models. Absorption may also be described by directly assuming velocity dispersion relations. 2.2.3 Describing Absorption by Quality Factors Wave speeds in rocks are determined by rock properties. When including absorption into the wave equation, it is very straightforward to use complex wave-numbers or complex velocities. The absorptive property determines a real velocity dispersion relation or a complex dispersion relation depending on different applications (Mittet et al., 1995; Claerbout, 1985). Definition of Q The absorptive property is often represented by the quality factor Q, which is an intrinsic property of rocks. Formally, Q is defined as a dimensionless measure of the anelasticity which is given by −∆E 1 = , Q 2πE (2.9) where E is the energy stored at the maximum strain in a volume, and −∆E is the energy loss in each cycle because of anelasticity. Equation (2.9) implies that Q−1 is the portion of energy lost during each cycle or wavelength. For a medium with a linear stress-strain relation, √ the amplitude A of a wave at a particular frequency is proportional to E (Aki and Richards, 1980), therefore −∆A 1 = , (2.10) Q πA0 where A0 is the amplitude at the start of a cycle and ∆A represents the amplitude decay in a cycle. It is observed that amplitude varies with distance because of absorption. We can rewrite A as a function of distance A(z) , and dA(z) , (2.11) ∆A = λ dz 12 Chapter 2. Mechanism and Description of Seismic Absorption where λ is the wavelength given in terms of frequency ω and phase velocity v(ω) by λ = 2πv(ω)/ω. (2.12) From equations (2.10), (2.11) and (2.12), it can be shown that A(z) = A0 exp − ωz , 2v(ω)Q (2.13) which describes amplitude attenuation of a frequency component in an absorptive medium. Attenuation Causing Velocity Dispersion Absorption results in both attenuation and dispersion (Liu et al., 1976). As mentioned before, the basic assumptions used to describe absorption behavior are constant Q, linearity and causality. Considering a unit impulse as input at depth z = 0, each frequency component of this impulse, ∞ −∞ δ[t − z/v(ω)]e−iωt dt = e−iωz/v(ω) , ω 2v(ω)Q . The proptransform eikz , and k is will now be attenuated by a factor e−α(ω)z with α(ω) = agating pulse shape G(z, t), thus, has a Fourier complex and is given by k= ω + iα(ω). v(ω) If the impulse response is causal, G(z, t) = 0, for t < z/v(∞), it can be shown that (Aki and Richards, 1980) ω ω = + H[α(ω)]. v(ω) v(∞) (2.14) Here, v(∞) is the limit of v(ω) as ω → ∞ and H[α(ω)] is the Hilbert transform of the attenuation factor. This equation is equivalent to the KramersKronig relation in electro-magnetic theory (Futterman, 1962). From a mathematical point of view, there is no Hilbert transform pair for which this relation can be satisfied with constant Q. We must tolerate a frequency dependent Q which is effectively constant over the seismic frequency range. 13 Chapter 2. Mechanism and Description of Seismic Absorption Under the constraint of causality, and the assumption of frequency independent Q, a velocity dispersion relation, due to attenuation, can be derived, 1 ω0 v(ω0 ) , =1+ ln v(ω) πQ ω (2.15) From equations (2.13) and (2.15), we can observe that wave propagation in an absorptive medium is completely specified by two parameters, Q and a phase velocity v(ω0 ) at an arbitrary reference frequency ω0 . ω0 often takes the value of the Nyquist frequency in digital signal processing for convenience. Kjartansson’s Absorption Model Another method to formulate the velocity dispersion relation was proposed by Kjartansson (1979). The absorption model is defined by means of a complex velocity dispersion relation −iω v(ω) = v(ω0 ) ω0 γ , (2.16) where γ is a parameter related to absorption. Equation (2.16) also represents the causal, constant Q model, and it can be derived 1 ≈ tan(πγ). Q The estimation of Q is equivalent to the estimation of γ. For γ = 0, equation (2.16) gives a real constant velocity which means no absorption. 2.2.4 Choice and Discussion Besides the foregoing mentioned mathematical models, there are also other models that describe seismic attenuation with differing degrees of precision and utility. A list of these methods can be found in the paper by Toverud and Ursin (2005). In this thesis, the dispersion relation of equation (2.15) is adopted. As pointed out by Aki and Richards (1980), equation (2.15) is an important result since it appears to be a good approximation for a variety of attenuation laws in which Q is effectively constant over the seismic frequency range. Besides a reference velocity v(ω0 ), this visco-elastic model assumes attenuation and dispersion are determined by only one more parameter Q. Theoretically, this linear absorption model (equations (2.13) and (2.15), also 14 Velocity (m/s) Chapter 2. Mechanism and Description of Seismic Absorption 2000 1800 0 20 40 60 80 Frequency (Hz) 100 120 Figure 2.2: Velocity dispersion comparison with Q = 30. The red curve is for the linear Q model and the green one is for Kjartansson’s model. called the Kolsky-Futterman model (Wang and Guo, 2004), is superior to the Voigt model which has been contradicted by practical observations. Application is easier than the generalized linear solid model which requires the adjustment of two parameters when being applied to different media. The linear Q model and Kjartansson’s model define different complex wavenumbers. Their real parts are very close. However, the imaginary parts are different. Figure 2.2 and Figure 2.3 compare the velocity dispersions and amplitude attentions of theses two models. Velocity dispersion and amplitude attenuation are both defined by Q. The red curves result from the linear Q model and the green ones result from Kjartansson’s model. For the same Q value, Kjartansson’s model has stronger attenuation than the linear Q model. In seismic data processing, ideally, Q estimation and Q compensation should be based on the same absorption model. In this research, we only consider the comprehensive effect of attenuation, so that, estimate and use only the effective Q values. This is based on the premise that scattering attenuation and intrinsic attenuation have the same effect on observed seismic data. When using Q to describe seismic attenuation, it is assumed that Q does not change with frequency. This assumption is valid for the seismic data acquired for oil and gas exploration which is in the frequency range of, approximately from 10 to 200 Hz (Sams et al., 1997). Figure 2.4 illustrates the wave behavior in elastic media and in absorp- 15 Amplitude Chapter 2. Mechanism and Description of Seismic Absorption 0 20 40 60 80 Frequency (Hz) 100 120 Figure 2.3: Amplitude attenuation comparison with Q = 30. The red curve is for the linear Q model and the green one is for Kjartansson’s model. tive media. An elastic medium is represented by ρ and v. A visco-elastic medium is represented by ρ, v and Q. In order to extract Q information, we need to examine the spectrum variations in observed seismic data. Without absorption, the wavelet tends to keep its shape during its propagation. With absorption, the wavelet gets distorted during its propagation due to the attenuation of high frequencies. 2.2.5 Absorption and Seismic Wave Extrapolation Absorption causes energy loss and frequency dispersion. From equations (2.13) and (2.15), it can be concluded that for a seismic wave traveling in a visco-elastic medium, the wavenumber including an attenuation term is complex: 1 ω0 1 ω 1+i ln 1+ . (2.17) kc (ω) = v(ω0 ) πQ ω 2Q The real part in this equation is responsible for the wavefront propagation and the phase distortion, the imaginary part is responsible for amplitude attenuation. For one-dimensional wave propagation, if the wavefield U (z, ω) at a depth z propagates downward to a depth z + ∆z, the wavefield U (z + ∆z, ω) at this depth in absorptive media is U (z + ∆z, ω) = U (z, ω)A(∆z, ω) exp −i ω ∆z , v(ω0 ) (2.18) 16 Chapter 2. Mechanism and Description of Seismic Absorption Figure 2.4: Seismic wave behavior in absorptive media defined by v, ρ and Q. (a) Without absorption. (b) With absorption. In equation (2.18), A(∆z, ω) is an absorption related term A(∆z, ω) = A1 (∆z, ω)A2 (∆z, ω), with A1 (∆z, ω) = exp − ω∆z 2Qv(ω0 ) (2.19) (2.20) defining an amplitude attenuation term, and A2 (∆z, ω) = exp −i ω∆z ln πQv(ω0 ) ω ω0 (2.21) expressing a velocity dispersion term which causes a phase delay. For time domain wavefield modeling, equation (2.18) can be written as U (t + ∆t, ω) = U (t, ω)A(∆t, ω) exp(iω∆t), (2.22) which means that the frequency content of a wavefield is time dependant. Using equation (2.22), we can examine the impulse response in an absorptive medium. Figure 2.5 shows five snapshots of an impulse response in a 1-D absorptive medium at five time locations which are 400, 800, 1200, 1600 and 2000 milliseconds respectively. The impulse initiates at time zero and the medium has a quality factor of 30. The shape of the impulse becomes increasingly wider as traveltime increases and its amplitude decreases continuously due to absorption. 17 Chapter 2. Mechanism and Description of Seismic Absorption 5 Trace Index 4 3 2 1 0 400 800 1200 1600 Time (MS) 2000 2400 Figure 2.5: Evolution of a unit impulse with time in a 1-D absorptive medium (Q = 30). 18 Chapter 2. Mechanism and Description of Seismic Absorption 0 -1.0 -0.5 Distance (km) 0 0.5 1.0 Time (sec) 0.5 1.0 1.5 Figure 2.6: A Green’s function in an absorptive 2-D medium. For 2-D modeling, the synthetic records can be obtained by two step calculation. The first step does ray-tracing to calculate ray-paths in the media and the second step extrapolates a wavelet along the ray-paths. Figure 2.6 shows the response of a diffractor in a 2-D absorptive medium with a quality factor which is also equal to 30. Because the medium is uniform, the modeling does not involve ray-tracing, but a zero-phase wavelet is convolved with the impulse response in the records. Without absorption, wavelets on all traces should be the same, because no geometrical spreading, angle-dependent reflectivity or other factors are involved in the calculation. With absorption, the signal at a far offset is weaker and broader than that at a near offset. 2.3 Summary Seismic attenuation is a complicated process, which can be caused by heterogeneity, interfaces, lithology, saturation, fluid content and other properties of the earth media. Mathematically, the absorption equation can be derived using dynamic mechanics or using an attenuation-dispersion relation. In 19 Chapter 2. Mechanism and Description of Seismic Absorption this thesis, the linear Q model (equation (2.13) and (2.15)) is adopted as a reasonable assumption over the exploration seismic frequency bandwidth. 20 Chapter 3 Quality Factor Estimation from Seismic Peak Frequency Variation The fundamental idea to estimate Q from seismic data is, at first to examine the variation of a seismic wavelet with time, and then calculate a Q value from the variation of the wavelet shape or spectrum. By assuming that the amplitude spectrum of a seismic wavelet is like that of a Ricker wavelet, this chapter derives an analytical formula which allows a Q value to be calculated from the spectral peak frequency variation. Theoretically, this analytical approach can be used to estimate Q-factors from seismic data of different acquisition geometries, where peak frequency variations can be reliably examined. Applications of this approach to estimate Q-factors from common midpoint (CMP) gathers and also from stacked seismic traces will be discussed. 3.1 Introduction Seismic waves traveling through the earth experience absorption, i.e., attenuation and dispersion, because of the anelasticity and heterogeneity of the medium (Ricker, 1953; Futterman, 1962; White, 1983; Kneib and Shapiro, 1995). Understanding, estimating, and compensating for absorption of seismic waves is important in the quest to improve the resolution of seismic images. This in turn will allow us to better understand the effects of AVO and, consequently, to realize the quest of inversion for material properties. To compensate for absorption, we require an estimate of the Q-factor. Methods for estimating Q values from surface seismic data are not well developed (Bachrach et al., 2006; Lancaster and Tanis, 2004; Dasgupta and Clark, 1998). However, some research has been published concerning the estimation of Q-factor from vertical seismic profile (VSP) and cross-well data (Spencer et al., 1982; Tonn, 1991). Almost all of these methods use 21 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation the amplitude information of received signals, which information, however, is often inaccurate because of noise, geometrical spreading, scattering, and other effects. Quan and Harris (1997) present a method for estimating seismic absorption based on the frequency shift observed in vertical seismic profiling (VSP) data. Since the centroid frequency of the signal’s spectrum experiences a shift to lower frequencies during propagation, they develop a relation between Q and the centroid of an amplitude spectrum. The amplitude spectrum is represented by a Gaussian, boxcar, or triangular function. In this chapter, an analytical equation will be developed to relate absorption to spectral peak frequency variation. This analytical approach can be used to estimate the Q-factors from prestack common mid-point (CMP) gathers. This approach can also be used to estimate the Q-factors from a single trace, where the spectral peak frequency variation of a seismic wavelet can be picked from a windowed time-variant spectrum (WTVS) or a continuous wavelet transform (CWT) spectrum. 3.2 Absorption and Peak Frequency Variation In seismic data processing, a recorded trace is commonly modeled as the convolution of a seismic source signature with a reflectivity series. The seismic source signature is generally unknown. The effects of anelasticity can be incorporated into the model by convolving with an earth filter. This filter is causal, minimum phase, and depends on the Q-factor (Aki and Richards, 1980). In order to build a relation between Q and peak frequency translation, we begin by assuming that the amplitude spectrum of the source wavelet can be well represented by that of a Ricker wavelet. 2 2 2 f2 B(f ) = √ 2 e−f /fm , π fm (3.1) where fm is the dominant frequency or the main frequency. For convenience, the frequency of maximum amplitude is referred to as the peak frequency, and denoted by fp . For a wavelet at its initial state, the peak frequency is the dominant frequency. One point needs to be emphasized here: the derivation only requires an assumption about the shape of the amplitude spectrum, without any assumption of the phase of a seismic wavelet. The Ricker amplitude spectrum is close to that observed in real data. Approximately, a spectrum which 22 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation rises to the peak amplitude quickly and decays gradually can be fitted by a Ricker spectrum. The evolution of the amplitude spectrum with time is now modeled as a Ricker wavelet traveling in a visco-elastic medium (geometrical spreading and other factors are not considered). After traveling for a time t, the amplitude spectrum is B(f, t) = B(f )H(f, t). (3.2) In equation (3.2), H(f, t) is the absorption filter (Varela et al., 1993), whose frequency response is H(f ) = exp − a(f, l)dl . ray In this expression, the integral is evaluated along the ray-path l, and a(f, l) = πf . Q(l)v(l) Here, Q(l) and v(l) are the Q-factor and velocity, respectively, defined along each point of the ray-path. The value Q(l) is assumed to be independent of frequency (Ricker, 1953; Kjartansson, 1979; White, 1983). By considering the propagation of a wavefield in a half-space with a Q-factor for t seconds, the amplitude spectrum of the received signal is t − πf Q B(f, t) = B(f )e . (3.3) This equation can also be derived from equation (2.22), by only considering the amplitude of the wavefield. As time increases, absorption increases with frequency and results in the peak frequency translating towards lower frequency. This phenomenon is illustrated in Figure 3.1. Because of absorption, the time width of the source wavelet increases and the amplitude spectrum narrows as the wave travels. If the arrival times are known, the Q-factor can be calculated from the spectral variation based on equation (3.3). 3.2.1 Spectral Ratio Method A common method to estimate Q is the spectral ratio method (Spencer et al., 1982; Tonn, 1991). By comparing the amplitude spectra at two arrival times 23 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation Trace index 10 15 5 20 25 Time (s) 0.5 1.0 1.5 Amplitude (a) 0.005 0 0 10 20 30 40 50 Frequency (Hz) 60 70 80 (b) Figure 3.1: (a) A reflection in a CMP gather. The medium is absorptive with Q = 10. Traces have been normalized based on their maximum amplitudes. (b) Amplitude spectra of the original source signature (black), trace 1 (red), trace 11 (green), and trace 21 (gray). 24 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation t1 and t2 respectively, we can take the ratio of the two amplitude spectra, − πf t2 B(f )e Q B(f, t2 ) = . πf t1 B(f, t1 ) B(f )e− Q (3.4) Taking logarithms of both sides, equation (3.4) becomes ln B(f, t2 ) π(t2 − t1 ) f. =− B(f, t1 ) Q (3.5) 2) Denoting Ar = ln B(f,t B(f,t1 ) , the logarithm of the spectral ratio, and plotting Ar as a function of frequency, f , yields a linear trend whose slope p is a function of Q. Q then can be easily expressed in terms of observed slope p as −π(t2 − t1 ) . (3.6) Q= p The spectral ratio method is simple in principle, but in practice determining the Q-factor is complicated by overlapping wavelets which lead to amplitude spectra that do not reflect the wavelet spectrum. Also, the linear trend of Ar to f is not obvious, and linear fitting is required. More seriously, the spectral ratio method will fail if there are other factors affecting the amplitude spectrum besides attenuation. These factors will be detailed in the following discussions. 3.2.2 Q and Peak Frequency Variation As stated above, the magnitude of the amplitude spectrum of seismic waves is affected by many factors besides absorption. These factors include geometrical spreading, reflection and transmission effects, and the methods of amplitude recovery, which may be applied to balance the trace amplitude in seismic data processing. It is almost exclusively the Q-factor, however, that affects the shape of the wavelet spectrum. Based on this idea, a method is developed to estimate the Q-factor from the spectral peak frequency variation of seismic reflections. Including all Q-factor unrelated functions into an amplitude term, the amplitude spectrum can be written as t − πf Q B(f, t) = M (t)B(f )e , (3.7) where M (t) is an amplitude factor independent of frequency and absorption. The peak frequency fp can be determined by equating the derivative of the 25 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation spectrum, with respect to frequency, to zero: πf t ∂B(f, t) ∂B(f ) − πfQt + M (t)B(f )e− Q = M (t) e ∂f ∂f − πt Q = 0. (3.8) Recalling the expression for the Ricker wavelet, from equation (3.1), we obtain ∂B(f ) 2 =√ ∂f π 2f 2 fm − e f2 2 fm 2 + 2 f f2 2 fm − e f2 2 fm −2f 2 fm . (3.9) Finally, by inserting this expression into equation (3.8), the peak frequency at time t is obtained as 2 2 πt 1 πt 2 . − + (3.10) fp = fm 4Q fm 4Q The relation between Q and the shift of peak frequency is, Q= 2 πtfp fm . 2 − f 2) 2(fm p (3.11) This expression shows that if the dominant frequency fm is known, the Qfactor can be computed from the peak frequency at only one time location. In practice, of course, the initial fm is not known. However, it can be estimated if the amplitude spectrum of the initial source wavelet can be approximated by that of a Ricker wavelet. Designating the peak frequencies at times t1 and t2 by fp1 and fp2 , respectively, we have Q= 2 2 πt2 fp2 fm πt1 fp1 fm = . 2 2 −f ) 2 − f2 ) 2(fm 2(fm p1 p2 Thus, the dominant frequency of the source wavelet can be derived from the peak frequencies of a reflection at two different time locations: fm = fp1 fp2 (t2 fp1 − t1 fp2 ) . t2 fp2 − t1 fp1 (3.12) Combining equation (3.11) and (3.12), one obtains 2 π(t2 − t1 )fp2 fp1 Q= 2 − f2 ) . 2(fp1 p2 (3.13) 26 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation Equation (3.13) allows us to determine the absorption property of a layer if an event in the layer is observed at two different time locations. Although a Ricker amplitude spectrum is assumed in the derivation, other wavelet spectral models are also possible, depending on their similarity to the observed amplitude spectra. One example is the Gaussian spectrum of variance σ 2 : −(f − fm )2 . (3.14) B(f ) = exp σ2 Following the above procedure, the relationship between the Q-factor and the variation of peak frequencies is Q= πtσ 2 , fm − fp and the dominant frequency is fm = fp1 t2 − fp2 t1 . t2 − t1 Similarly, observations at two different time locations allow us to determine a unique Q-factor πσ 2 (t2 − t1 ) . (3.15) Q= fp1 − fp2 By assuming the amplitude spectrum of a seismic wavelet has an analytical form which is unimodal, other expressions like equation (3.13) can be derived. For real seismic data, the amplitude spectrum of a seismic wavelet is not exactly that of a Ricker or of a Gaussian wavelet. In many situations, however, it can still be approximated by a Ricker spectrum (Ricker, 1953). This approach may be used to estimate quality factors from both prestack and poststack seismic data. 3.3 Q Estimation for Seismic Absorption Compensation and Migration For absorption compensation and migration in seismic data processing, Q values are only expected for major layers defined by velocities. Therefore, methods to estimate Q values here are based on seismic reflection event picking. 27 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation 3.3.1 Estimation of Q from CMP Gathers A CMP gather is very suitable for Q estimation, for two reasons: First, a CMP gather represents multiple observations of an underground structure. It provides information in the time and offset domains, allowing for the extraction of information concerning structure, lithology, and material properties such as velocity and Q-factor. Secondly, reflection arrival times are determined by interval velocities and the geometrical structure of the subsurface. Absorption of the received signals is only determined by the interval Q-factors and traveltimes in each layer. If the amplitude spectrum of a seismic wavelet is assumed to be Ricker-like, interval Q-factors can be computed solely from the variation of the peak frequency of a spectrum as a function of time. Equation (3.13) allows us to obtain an average Q-factor by using the peak frequency variation along all offsets, thereby allowing us to remove the effects of surface fluctuations and random noise, consequently improving the accuracy of the estimated Q values. Two-layer Case Consider first the case of two layers with quality factors Q1 and Q2 and traveltimes t1 and t2 in each layer, respectively, using equation (3.3), the amplitude spectrum after t = t1 + t2 is − B(f, t) = A(t)B(f )e πf t1 Q1 − e πf t2 Q2 . (3.16) From the peak frequency fp associated with B(f, t) in equation (3.16), Q2 can be estimated by following the same procedure of derivation from equation (3.7) to (3.11), if Q1 and the dominant frequency fm of the source wavelet are known. The expression of Q2 is Q2 = where πt2 Q1 , αQ1 − πt1 (3.17) 2 − 2f 2 2fm p . 2 fp fm α= (3.18) Now, using the concept of an equivalent Q, in other words, letting − e πf t1 Q1 − e πf t2 Q2 = e− πf t Q , (3.19) 28 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation this allows an expression to be derived which relates the interval Q2 to the equivalent Q-factor by Q2 = t 2 Q1 Q . (t1 + t2 )Q1 − t1 Q (3.20) The equivalent Q can be calculated from equation (3.19), using the peak frequency at time t. Since the dominant frequency fm of the initial wavelet and Q1 have already been determined from upper-layer arrivals, if traveltimes t1 and t2 in the two layers are known, Q2 can be computed from equation (3.20). Multi-layer Case Given the complexity of the subsurface, some approximations must be made to use all offset information. For multi-layer media, we may write the amplitude attenuation equation (3.3) as N B(f, t) = A(t)B(f ) exp i=1 − πf ∆ti Qi , (3.21) where Qi and ∆ti are the quality factor and the traveltime in layer i , respectively. Following the approach used in velocity estimation, we assume straight raypaths and compute the total traveltime of a reflection at a particular offset as N ∆ti = tN . i=1 Now define t0 (N ) as the zero-offset traveltime of reflection N , QN as the Qfactor value for layer N , and t0 (i) as a zero-offset traveltime of a reflection above layer N . Using the proportional property of similar triangles, the traveltime in layer i is ∆ti = tN [t0 (i) − t0 (i − 1)] . t0 (N ) Furthermore, the amplitude attenuation operator in equation (3.21) can be split into two factors as follows: B(f, t) = A(t)B(f ) exp N −1 i=1 −πf ∆ti Qi exp −πf ∆tN Qn . (3.22) 29 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation This allows us to obtain the following equation for QN : QN = π∆tN . α−β (3.23) where α is the same as that in equation (3.18) and β= N −1 i=1 π∆ti . Q The Q-factor values can now be calculated layer by layer by means of layer stripping. RMS Q and Interval Q Since a straight raypath approximation is used, the computed QN is not the actual interval Q-factor. Analogous to root mean square (RMS) velocities in seismic velocity analysis, such Q values can be referred to as RMS Q values. In a manner similar to interval velocity analysis, we can determine the apparent interval Q-factor of the ith layer, Qint i , from the RMS values using a relation similar to the Dix formula (Yilmaz, 2000): Qint i = Q2i t0(i) − Q2i−1 t0(i−1) . t0(i) − t0(i−1) (3.24) When using this equation to estimate Q-factors from CMP gathers, it is better to notice that Qi is the average of the calculated Q values at different offsets for layer i, and t0 (i) is the zero-offset arrival time of reflection i when the raypath is assumed to be straight across interfaces. Moreover, since values of Qint are derived from RMS Q values, they are apparent, rather than real, interval Q-factors. Numerical Experiments To estimate Q from CMP gathers, it is assumed that the arrival times for the main reflection events are known. Fourier transforms are computed in the window containing the reflection at each offset, each amplitude spectrum is fitted with a Ricker spectrum and the peak frequency of the spectrum is estimated. Using equation (3.23), we can now calculate Q factors layer by layer from the peak frequency variations. The procedure of using the suggested approach to estimate Q-factors from CMP gathers can be outlined by the the pseudo code in table (3.3.1). 30 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation /* Analyze absorption at selected CMP locations */ { select_a_CMP (); do while (event index < number of events) { pick_an_event (); do while (offset index <number of offsets) { do_windowed_fft (); spectrum_fitting (); check_peak_frequency (); } calculate_Q (); } form_attenuation_field (); } Table 3.1: A pseudo code of prestack Q estimation. Figure 3.2 shows a simple test on a synthetic CMP gather with two events. The original Ricker wavelet has a dominant frequency of 60 Hz. Absorption is modeled by using low Q values to emphasize the absorption effect: the values in the two layers are 10 and 20, respectively. The actual and estimated inverse Q curves are shown in Figures 3.2b and 3.2c, respectively. The two curves are almost the same, because the estimated Q values are close to the actual values. The corresponding dominant frequency is estimated as 60.67 Hz. The agreement between the corresponding values shows that the method works well for ideal synthetic data. The variation of peak frequencies with offset is shown in Figure 3.3. Peak frequency variation with offset of the first event is shown in black and the variation of the second event is shown in red. The remarkable decrease is, of course, due to absorption. This example demonstrates clearly why absorption must be compensated for if resolution of prestack data is an issue. In theory, an event that can be examined at two offset locations allows us to determine the Q value of one layer. By using the peak frequency variation along all offsets, an average Q can be obtained, thereby the effects of surface fluctuations and random noise can be reduced, thus improving the accuracy of Q estimation. 31 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation 1 1/Q 0 0 1 2 2 (a) 1/Q 0.2 0 Time (s) Offset (m) 1000 Time (s) Time (s) 0 0 0 0.2 1 2 (b) (c) Figure 3.2: (a) A synthetic CMP gather of two reflections with absorption and 10% random noise. Traces have been normalized based on their maximum amplitude. (b) Values of input quality factors (Q1 = 10 and Q2 = 20 ). (c) Computed quality factors (Q1 = 10.04 and Q2 = 20.12). 32 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation Peak Frequency (Hz) 20 10 0 0 500 Offset (m) 1000 Figure 3.3: Variation of peak frequencies of the two reflections in Figure 3.2. Peak frequencies of the first event and the second event are shown in black and red respectively. 3.3.2 Q Estimation from a Poststack Trace In seismic data processing, stacking after normal moveout (NMO) destroys the original absorption characteristics by summing together signals which experience different trajectories. In short, prestack CMP gathers before deconvolution are more suitable for attenuation analysis than stacked traces. Nonetheless, to perform de-absorption processing on poststack data is often necessary; hence, it is important to estimate Q-factors from vertically incident traces or stacked traces, if corresponding CMP gathers are not available. The forgoing derived analytical approach, which computes Q values from wavelet peak frequency variation, can be applied to different seismic acquisition geometries. To detect peak frequency variation from a single trace, we can use windowed time-variant spectral (WTVS) analysis or the continuous wavelet transform (CWT) analysis. Windowed Time-variant Spectral Analysis The windowed Fourier transform was introduced by Gabor (1946) to measure localized frequency components of sound. A real and symmetric window 33 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation g(t) = g(−t) is translated by u and modulated with a frequency ω: gu,ω = eiωt g(t − u). The resulting windowed Fourier transform of a signal f (t) is Sf (u, ω) = ∞ −∞ f (t)g(t − u)e−iωt dt. (3.25) This transform is also called the short time Fourier transform because multiplication by g(t − u) localizes the Fourier integral in the neighborhood of t = u. In the windowed Fourier transform, the window function remains unchanged as it moves along the time axis. The resolution in time and frequency of the windowed Fourier transform depends on the spread of the window in time. To examine the variation of a time-variant spectrum, the support of the window function needs to increase with time in order to accommodate the time-variant seismic wavelet. The windowed time-variant Fourier transform is a useful tool for seismic time-frequency analysis. Continuous Wavelet Transform The time-frequency localization of g in equation (3.25) can be modified by a scaling factor s. Let gs (t) = √1s g( st ) represent the dilation of g(t). The choice of a particular scale s depends on the desired resolution trade-off between time and frequency. To analyze signal structures of different size at different localizations, it is necessary to use time-frequency atoms with different time supports. The wavelet transform decomposes signals by means of dilated and translated wavelets. Supposing a mother wavelet is a function ψ(t) ∈ L2 (R) (Mallat, 1999), the wavelet transform of a signal f (t) at time u and scale s is W f (t, s) = ∞ ∞ f (t)ψ t−u s dt. (3.26) A set of wavelet bases which are wavelet functions at different scales is equivalent to a set of band pass filters. Orthogonal and complete wavelet bases can be constructed to perform multi-resolution analysis of signals (Mallat, 1999). It is important to note that the orthogonality and completeness of wavelet bases are not required for the application to find out the peak frequencies of signals. 34 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation Both windowed time-variant spectral analysis and the continuous wavelet transform can be used for time-frequency analysis (Zhang, 2000). The following presents the application of WTVS analysis for Q estimation from a single trace. Quality Factor Estimation from WTVS Analysis To analyze the absorption characteristic of a trace in a poststack seismic section, we use WTVS analysis because WTVS can reveal the variation of frequency content with time. If a trace is assumed to be at zero offset, from the peak frequency variation detected in a WTVS, Q values can be calculated using equation (3.13). For a seismic section, Q can be estimated in many CMP locations. By means of horizontal interpolation, a Q profile can be constructed. The Q profile is like a velocity profile, it can be used in inverse Q filtering and also in Q-migration. Figure 3.4(a) shows a real seismic trace. Figure 3.4(b) shows a windowed time-variant spectrum of this trace. This WTVS gives a clear indication of the trend of the spectral variation. Picking the peak amplitude ridge from WTVS and fitting the ridge with a piece-wise straight line in the coordinates of f versus t, one can calculate a quality factor from each line segment using equation (3.13) for the main events, as shown in Figure 3.4(c). The WTVS computation initially applies a narrow window to the input trace and calculates the conventional Fourier spectrum of the windowed data. The window is then translated and widened successively with time along the trace and the Fourier transform is calculated for each new position of the window. 3.4 Summary and Discussions This chapter has presented an analytical formula which calculates a Q value from the spectral peak frequency variation of a seismic wavelet. The idea is illustrated using the spectrum of a Ricker wavelet. The approach can be applied to prestack CMP gathers or to postack traces for Q-factor estimation. This approach is superior to the traditional spectral ratio method because it only uses the information of frequency variation caused by absorption, instead of using amplitude information which is affected by many factors. Seismic absorption properties are ideally expected to be estimated from prestack CMP gathers, and absorption compensation is ideally implemented on prestack seismic data. The reason is that stacking distorts the frequency information of the original seismic data. 35 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation (a) (b) (c) Figure 3.4: Q estimation. (a) is a real seismic trace. (b) is its windowed time-variant spectrum. The picked peak frequency points are marked by circles and connected by straight lines. (c) is the inverse of estimated quality factors. 36 Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation Underground lithology is characterized by its velocity, density, and Qfactor. While the Q-factor does not affect the arrival times of reflections, it does affect amplitude and the signal’s frequency content. Extracting velocity information from CMP gathers is a common practice. we have devised an analogous approach to estimate the absorption character from the variation of spectrum with both time and offset. The equations that compute Qfactors from observed peak frequency variation are derived based on the reasonable assumption of a Ricker-like amplitude spectrum. Absorption effects as a function of offset may be mistakenly interpreted as AVO phenomena. Without prestack Q-factor compensation, the reflection amplitudes may appear to decrease with offset. Compensating for absorption in prestack data can improve the resolution of seismic data and is important for AVO analysis. Figure 3.2 clearly shows that signals get attenuated with traveltime and offset. Absorption compensation, based on the underlying physics of the absorbing process, lets us compute a section with increased resolution. Subsequently, a stationary source wavelet can be used to improve the resolution by means of deconvolution. It is superior to conventional deconvolution techniques which rely on an adaptive formulation of wavelet estimation. Only when Q values are available, can prestack inverse Q-filtering be applied to compensate for attenuation on prestack seismic data. It is advisable to extract Q-factors from prestack seismic data that have not been processed by frequency filtering methods (e.g. deconvolution) that violate the assumptions of the method. Prestack processing provides another dimension “offset” that increases the information available for signal and noise separation. Q compensation performed in the prestack domain can provide a balanced signal spectrum over offset and over time as well, improving the vertical and lateral resolution of the final image. Although seismic absorption estimation is ideally implemented on a CMP gather, the introduced method applies to seismic data with different acquisition geometries, wherever, the peak frequency of a seismic wavelet can reliably estimated. For poststack data, windowed time-variant spectral analysis or the continuous wavelet transform can be used for time-frequency analysis. In absorptive media, dispersion causes the width of a seismic wavelet to increase with time. To allow the window to enclose an entire seismic wavelet at each time location, WTVS analysis uses a time-variant window function, i.e., the width of the window function increases with time. 37 Chapter 4 Seismic Absorption Analysis for Reservoir Description In reservoir description, seismic attenuation property, Q, is ideally estimated at all prospective layers defined by velocities or acoustic impedance. This chapter introduces a technique called reflectivity-guided seismic attenuation or Q analysis (RGQA) which is based on the idea that Q values can be calculated from the peak frequency variation of a seismic wavelet. We first estimate peak frequencies at a CMP location, then correlate the peak frequency with sparsely-distributed reflectivities, and finally calculate a Q curve from the peak frequencies at the layer interfaces using an analytical equation. The peak frequency is estimated from the prestack CMP gather using peak frequency variation with offset (PFVO) analysis which is similar to amplitude variation with offset (AVO) analysis to reduce the stacking effect on the waveform. PFVO analysis also helps to exclude the effect of event tuning on frequencies. The estimated Q section has the same layer boundaries as the acoustic impedance or other layer properties and can be easily interpreted for the purpose of reservoir description. 4.1 Introduction Seismic attenuation property, usually represented by the quality factor Q, is an important parameter for characterizing rock types. Velocity and attenuation together define the propagating behavior of seismic waves through visco-elastic media. One can expect that attenuation properties estimated in layers defined by velocities or acoustic impedance would be a useful seismic attribute for reservoir description. This chapter first reviews techniques used for seismic attenuation analysis in reservoir description, and then introduces a technique to estimate Q values based on the variation of zero-offset seismic peak frequencies at a CMP location. In order to avoid the distortion in seismic waveforms caused by stacking and to exclude frequency changes caused by event tuning, the 38 Chapter 4. Seismic Absorption Analysis for Reservoir Description peak frequency and its derivative along offset are estimated from a CMP gather using peak frequency variation with offset (PFVO) analysis. Reflectivities are used to define the locations where peak frequencies are required and the layers where Q values will be computed. The technique is called reflectivity-guided seismic Q analysis. It operates on prestack CMP gathers, while in the meantime, it uses the reflectivity (or the interface of acoustic impedance) information, which is usually obtained from stacked seismic section. The estimated Q section can be used as an attribute which helps to delineate prospective hydrocarbon reservoirs. 4.2 Review of the Methods for Seismic Attenuation Analysis There are mainly three categories of techniques to estimate seismic attenuation properties for hydrocarbon reservoir description. They are iso-frequency analysis, direct Q estimation and attenuation or, Q tomography. 4.2.1 Iso-frequency Analysis This technique is regularly implemented using spectral decomposition to decompose a seismic trace into a frequency gather. By examining iso-frequency sections, one can determine frequency dependent energy variations, among which some could be related to hydro-carbon accumulation (Castagna et al., 2003). It is difficult to set a suitable frequency interval for iso-frequency slices so that the low-frequency shadows may be properly identified. In theory, there are numerous mechanisms that might be responsible for the low frequency shadows beneath hydrocarbon reservoirs, as discussed by Ebrom (2004). It is difficult to determine which mechanisms are, in fact, those that introduce the first-order effects. Based on these observations, it may perhaps be concluded that iso-frequency analysis is more a qualitative than quantitative approach. In order to illustrate the above discussion, Figure 4.1 shows a frequency shadow caused by attenuation, whereas Figure 4.2 shows low-frequency shadows caused by attenuation and thin-layer tuning. In Figure 4.1, the only absorptive layer is between 0.4s and 0.6s. The high frequencies are attenuated after the seismic wavelet traveling through the absorptive layer. The frequency decay is irreversible, because there is no physical mechanism to boost the frequency except the interference among waves. In Figure 4.2, there are two closely spaced (0.008s) reflection coefficients at around 0.6s. 39 Chapter 4. Seismic Absorption Analysis for Reservoir Description 0 0 1/Q 0.03 0 0 0 Frequency (Hz) 30 60 90 0.5 Time (s) Time (s) Time (s) 0.03 0.5 0.5 0.02 0.01 0 1.0 1.0 1.0 (a) (b) (c) Figure 4.1: Frequency decay caused by absorption. (a) is an attenuation model, (b) is the synthetic trace and (c) is the spectrum of frequency decomposition. The wavelets from the two layers merge and form a composite wavelet of lower frequency, as shown in Figure 4.2(b). In the time-frequency spectrum of Figure 4.2(c), the composite wavelet appears as a low frequency shadow. This frequency translation is, however, only apparent in that it does not belong to the propagating wavelet. These two figures clearly demonstrate that the low-frequency shadow itself is not accurate enough to define an absorptive zone. 4.2.2 Direct Q Estimation Conceptually, attenuation can be fully defined by the quality factor, Q. Q can be extracted from seismic traces or from sonic logs, and represents a straightforward way to estimate seismic attenuation for reservoir description. A number of techniques have been introduced to estimate Q in reservoir applications. Dvorkin and Mavko (2006) calculated maximum Q values based on high and low frequency P -wave bulk moduli. The extracted Q curve from sonic logs confirms that the gas hydrate area has strong attenuation. A case study presented by H¨ ubert et al. (2005) explains how seismic attenuation and relative acoustic impedance help to identify the gas reservoir and exclude the AVO bright spot caused by Cretaceous shales. Parra et al. (2006) extract intrinsic attenuation from the head P -wave of 40 Chapter 4. Seismic Absorption Analysis for Reservoir Description 0 0 0 0 1/Q 0.03 0 0 0 Frequency (Hz) 30 60 90 0.5 Time (s) 0.5 Time (s) Time (s) Time (s) 0.04 0.5 0.5 0.02 0 1.0 1.0 (a) 1.0 1.0 (b) (c) (d) Figure 4.2: Frequency decay caused by thin-bed tuning and absorption. (a) shows a reflectivity function with closely-spaced coefficients. (b) is the corresponding absorption property. (c) is the synthetic trace and (d) is its spectrum of time-frequency decomposition. a full-waveform sonic log using a spectral ratio approach. The obtained Q log, together with other well logs, helps to discriminate between anomalies which are associated with lithology and those associated with oil and gas saturation. Singleton et al. (2006) introduce a Q estimation method using the Gabor-Morlet spectral ratio of an absorption compensated trace to the original seismic trace. The absorption compensation is achieved by spatial spectral balancing. 4.2.3 Attenuation Tomography This approach attempts to extract attenuation structures from seismic data by minimizing the misfit between model derived data and the observations. Quan and Harris (1997), as well as Plessix (2006), combine traveltime tomography and attenuation tomography together to invert velocity and attenuation from cross-well transmitted data. Seismic attenuation is modeled by means of the variation of the centroid frequency of a seismic wavelet. Pratt et al. (2003) use waveform tomography which is based on a visco-acoustic wave equation. This method helps to find a layered hydrate distribution with attenuation anomalies. Hicks and Pratt (2001), Watanabe et al. (2004) and Gao et al. (2005) also adopt waveform tomography to invert for attenuation. 41 Chapter 4. Seismic Absorption Analysis for Reservoir Description The three classes of techniques for attenuation analysis listed above, have all experienced a degree of success in application. However, the problem of how to obtain an interval intrinsic Q profile from seismic reflection data, or how to obtain a Q that can be easily correlated with acoustic impedance and/or other elastic parameters, remains a challenge in reservoir characterization. 4.3 Reflectivity-guided Q Analysis, RGQA 4.3.1 Theory The following four points outline the theory behind the technique that we have named reflectivity-guided Q analysis, RGQA in abbreviated form: 1. We assume a frequency-independent Q which is the only parameter required to fully describe the attenuation properties of a layer. This Q is related to lithology, porosity, pore fluid, saturation etc.. Together with other available geological and geophysical information, Q provides one more attribute for reservoir characterization. Like velocity, Q is an interval property, which can be derived from the reflections occurring at layer boundaries. 2. A Q curve can be estimated from the frequency variation of a zerooffset trace, which describes a vertically traveling seismic wavelet. Its time-variant spectrum describes the attenuation property of the medium. Frequency variations can be simply described by the peak frequency shift. Assuming the peak frequencies of a wavelet at time t1 and t2 are fp1 and fp2 respectively, if the amplitude spectrum of the seismic wavelet is like that of a Ricker wavelet (the second derivative of a Gaussian function), the inverse Q value can be calculated using (from equation (3.13)) 2 − f2 ) 2(fp1 1 p2 (4.1) = 2 . Q π(t2 − t1 )fp2 fp1 3. The zero-offset trace is not naturally available. Seismic stacking sums up traces in a CMP gather to reduce noise and improve overall data quality. Furthermore, stacking distorts the frequency variation pattern caused by seismic attenuation. The seismic wavelet at a CMP location in a stacked seismic section is different from what is observed on the trace of vertical-incidence. 42 Chapter 4. Seismic Absorption Analysis for Reservoir Description Figure 4.3: Flow diagram of reflectivity-guided Q analysis. 4. Required peak frequencies of seismic wavelets at zero-offset can be estimated by means of modeling the frequency variation along offset in a CMP gather. We call this technique PFVO analysis. Figure 4.3 shows the flow diagram of the RGQA algorithm. Attenuation analysis is implemented on a prestack CMP gather and requires sparsely distributed reflectivities for interface selections. Layer interface information can be obtained from acoustic impedance inversion or sparse reflectivity inversion (Oldenburg et al., 1983; Ulrych and Sacchi, 2005). Nowadays, piece-wise constant acoustic impedance information is generally available for reservoir description (Latimer et al., 2000). RGQA employs poststack acoustic information to help the extraction of attenuation information from the prestack data. Seismic data contain information concerning the interfaces between rock layers along the path of propagation. In general, the information concerning the layers themselves is inferred from that describing the interfaces. Acoustic impedance and AVO inversion are two instances where this strategy is used to infer layer properties. Since layer boundaries defined by reflection coefficients are also representative of the boundaries of the attenuating layers, RGQA derives rock attenuation properties pertinent to the layers themselves. In order to illustrate the general idea behind RGQA, we use the model in Figure 4.1 once again. Peak frequencies at four reflectivity locations 43 Chapter 4. Seismic Absorption Analysis for Reservoir Description 0.5 0 Time (s) Time (s) Frequency (Hz) 20 30 40 50 0 1.0 0 1/Q 0.03 0.5 1.0 (a) (b) Figure 4.4: Peak frequencies at layer boundaries. (a) is the peak frequencies from the trace in Figure 4.1 (b). (b) is the extracted inverse Q curve. are picked from the spectrum in Figure 4.1(c) and these frequencies are transformed to 1/Q values using equation 4.1. Figure 4.4 shows the peak frequencies at these locations, as well as the estimated 1/Q curve. For this simple case, RGQA works perfectly. This method, when applied to the event tuning model shown in Figure 4.2(b), will fail however. The frequency decay or rise caused by event tuning can be excluded by examining the peak frequency variation along offset (PFVO) on CMP gathers. 4.3.2 Peak Frequency Variation with Offset (PFVO) Analysis Considering that NMO correction followed by stacking could distort the shape of the wavelet and thus obscure the sought after absorption effects, we suggest a method to estimate peak frequencies at zero offset directly from prestack CMP gathers based on the frequency variation along offset. This multichannel peak frequency analysis strategy makes the estimation algorithm more robust to random noise and crossing events in a CMP gather. PFVO analysis follows the same implementation procedure as AVO (amplitude variation with offset) (Aki and Richards, 1980). It fits the peak frequency along offset linearly. The intercept, Pf , can be transformed into Q values, and the gradient, Gf , serves the fundamentally important, and thus far unresolved, task of separating Q into real attenuation and event tuning 44 Chapter 4. Seismic Absorption Analysis for Reservoir Description effects. PFVO operates on a CMP gather without NMO correction, because NMO stretch at far offset distorts the frequency spectrum of reflections. There are two main issues in Q estimation using frequency variation with offset. The first is how to determine the peak frequency at each layer boundary. The second is how to determine Q from the peak frequency variation. Regarding the first issue, one general and reliable way is to use spectral decomposition, i.e. to perform time-frequency analysis on each trace, and use the instantaneous amplitude spectrum to determine a peak frequency. Regarding the second issue, we fit the peak frequencies of an event at different times by using a formula like equation (4.1) depending on the spectral shape of a seismic wavelet. For a Ricker-like amplitude spectrum, which is the assumed spectral shape in our work, the peak frequency variation with time follows the following equation (also equation (3.12)) 2 2 1 π∆t πt 2 − + . fp = fm 4Q fm 4Q where fm is the peak frequency of a wavelet at its initial state (time = 0) and fp is the peak frequency after travel time t. It can be observed from this equation that if the peak frequency at time t − ∆t is fp0 after a time interval ∆t, the peak frequency is 2 fp0 π∆t fp π∆t fp = fp0 +1− 0 . (4.2) 4Q 4Q f π∆t In equation (4.2), the term p04Q is almost always less than 1 for small ∆t and weak attenuation (large Q). In a further step, we expand the square root via a Taylor series and keeping only the first term we obtain fp = fp0 1 − fp0 π∆t fp20 π 2 ∆t2 + . 4Q 32Q2 In general, peak frequency fp varies with the difference of arrival times nonlinearly. For weak attention, the Q2 term can be omitted from the above expression, yielding the following approximation: fp = Pf − Gf ∆t, (4.3) πf where Pf = fp0 , Gf = 4Qp0 . For an event in a CMP gather, ∆t is the normal moveout for a particular offset. For an event starting from two-way 45 Chapter 4. Seismic Absorption Analysis for Reservoir Description traveltime t0 at zero-offset, ∆t is a function of t0 , RMS velocity and offset, i.e. ∆t = ∆t(t0 , vrms , offset). Equation (4.3) shows that for a seismic wavelet with a Ricker-like amplitude spectrum, the peak frequency decay of a reflection with the difference of arrival times in absorptive media can be simulated by a linear relation. The intercept, Pf , of the straight line is the peak frequency at zero offset, and the gradient, Gf , is related to seismic attenuation. A first glance at equation (4.3) may give the impression that Q could be computed directly from Gf . Actually, Gf can only provide root-mean square attenuation information, (Zhang and Ulrych, 2002), rather than the required interval Q’s. Instead of providing attenuation information directly, Gf indicates whether there is peak frequency decay along offset. Therefore, we use Gf as an indicator which tells us whether the frequency decay at zero offset is caused by attenuation or by some extrinsic mechanism, e.g. thin-bed scattering. In fact, in our method, Q is calculated from the Pf ’s only. 4.3.3 Implementation Steps The technique is implemented in six steps as outlined below: 1. Build a CMP super-gather to suppress the effect of random noise. A trace in a super-gather is a stacked trace of several neighboring traces (Ostrander, 1984). The CMP gather used for Q-factor estimation needs to be free of NMO stretch, so that if normal moveout is applied when forming a CMP super-gather, reverse NMO needs to be applied afterwards. 2. Do time-frequency analysis on each trace using a continuous wavelet transform or short window Fourier transform. More sophisticated spectral decomposition like exponential pursuit decomposition (Castagna and Sun, 2006) may help to improve the resolution of the time frequency spectrum. Resolution and computation are two factors for the choice. 3. Compute the instantaneous amplitude envelope using Hilbert transform (Taner et al., 1979) to get a smooth amplitude spectrum. 4. Examine for the highest amplitude among frequencies at interface locations. The peak frequency is the frequency corresponding to the maximum amplitude. 46 Chapter 4. Seismic Absorption Analysis for Reservoir Description 5. Fit peak frequencies with normal moveout linearly using l1 norm to get Pf and Gf , and edit Pf by removing the values corresponding to negative or very small Gf ’s. 6. Calculate Q from the effective Pf values. Theoretically, peak frequencies at two offset locations can fully define the trend of peak frequency variation of a reflection. By using multi-channel information along offset, the algorithm is robust to random noise and the effects of crossing events. Forming super-gathers, spectral decomposition and l1 linear fitting are the main parts of the whole PFVO computation. Compared to traditional AVO analysis, PFVO analysis does fewer linear fitting operations, because linear fitting is only performed at layer interface locations, and the extra computation is efficient reverse NMO and the spectral decomposition. The wavelet transform can be regarded as a series of band-pass filters. Because the filters can be precomputed before the convolutions, the computation load of spectral decomposition is limited. PFVO analysis on a CMP gather is computationally comparable to traditional AVO analysis. The reflectivityguided Q analysis is practical even for large 3D prestack data volumes. 4.4 Numerical Experiments First, let’s look back at the model with event tuning and attenuation in Figure 4.2 to see whether the method of the reflectivity-guided attention analysis can help to extract actual attenuation information. Figure 4.5(a) shows a synthetic CMP gather using the reflectivity model (Figure 4.2(a)) and the attenuation model (Figure 4.2(b)). The data at the near offsets on the CMP gather are purposely zeroed to simulate real seismic data acquisition, where zero-offset data are not recorded. Pf and Gf are extracted from the CMP gather. Pf and Gf /Pf are displayed in Figure 4.5(b) and Figure 4.5(c) respectively. The reason for displaying Gf /Pf rather than Gf is that the value of Gf /Pf falls in a relatively narrow range. Gf /Pf is called the relative gradient. The low frequency caused by tuning events at 0.6s is excluded from the Q computation because it does not decay with offset (Gf /Pf ≈ 0). The estimated Q curve is shown in Figure 4.5(d) which matches the original Q model (Figure 4.2(b)). The technique is also tested on a real 2D dataset. Besides the prestack CMP gathers, an acoustic impedance or reflectivity section is required to provide the layer-boundary information. Figure 4.6 shows the test on a 47 Chapter 4. Seismic Absorption Analysis for Reservoir Description 1.0 0.5 1.0 (a) 0 Relative Gradient 0 2 4 0.5 1.0 (b) 0 Time (s) 0.5 Frequency (Hz) 20 30 40 50 0 Time (s) Offset (m) 500 1500 Time (s) Time (s) 0 0 1/Q 0.03 0.5 1.0 (c) (d) Figure 4.5: Peak frequency variation with offset. (a) is a CMP gather built from the model shown in Figure 4.2. (b) is the corresponding peak frequencies at zero offset. (c) is the relative gradient of peak frequencies along offset, and (d) is the estimated Q curve. CMP gather, the super-gather having been formed in a manner similar to that used in AVO analysis (Ostrander, 1984), except with reverse NMO applied to remove NMO stretch. Figure 4.6(b) shows the peak frequency gather. Figure 4.6(c) shows the peak frequencies at zero offset. Figure 4.6(d) illustrates the relative gradients of peak frequency along offset where positive values indicate frequency decay along offset. Figure 4.6(e) is the sparse reflectivity obtained using the Bayesian inversion method described in Ulrych and Sacchi (2005). Figure 4.6(f) is the estimated Q curve which is supposed to be piece-wise constant. In practice, the linear fitting step in PFVO analysis needs only to be done at interface locations. The continuous curves in Figure 4.6(c) and Figure 4.6(d) are computed at all samples. The reason for doing so is only for the purpose of display, and this reason also explains why Figure 4.7(a) shows a peak frequency section. For this 2D dataset, the calculated peak frequencies at all CMP locations are shown in Figure 4.7(a) and the estimated attenuation section is shown in Figure 4.7(b). The attenuation section displays the values of Q−1 . We use the inverse of Q since it lies in the range of (0, 1) and is ideal for the purpose of display. The blue background indicates small Q−1 values which implies smaller absorption than the red areas. 48 Chapter 4. Seismic Absorption Analysis for Reservoir Description Offset (m) 500 1000 0.5 Time (s) 1.0 1.5 2.0 2.5 (a) 10 Pf 20 (b) Gf/Pf -4 -2 0 2 Reflectivity -0.2 0 0.2 0.5 Time (s) 1.0 1.5 2.0 (c) (d) (e) (f) Figure 4.6: A real data example. (a) is a CMP gather. (b) is the corresponding non-stationary distributed peak frequencies. (c) is the peak frequency curve at zero offset. (d) is the relative gradient curve of peak frequencies. (e) is the reflectivity function and (f) is the estimated Q curve. 49 Chapter 4. Seismic Absorption Analysis for Reservoir Description (a) (b) Figure 4.7: Real data experiment. (a) is the peak frequency section and (b) is the attenuation section (Q−1 ). 50 Chapter 4. Seismic Absorption Analysis for Reservoir Description 4.5 Discussion and Conclusion In this chapter, we have introduced a technique of reflectivity-guided Q analysis. This approach uses poststack sparse reflectivity (the interfaces defined by acoustic impedance) to help pick the layer boundaries on prestack CMP gathers, where peak frequencies are selected to compute Q values. This technique basically provides a methodology that allows Q values to be estimated at all prospective layers. The obtained Q section may be computed with the same resolution as other attribute sections. Therefore, the absorption property can be easily incorporated into acoustic impedance interpretation. We do not suggest here the use of attenuation as a direct hydro-carbon indicator but rather as an attribute which could be combined with other parameters to give a better understanding of rock properties. We have also introduced an AVO-like, peak frequency estimation scheme called PFVO. By examining the peak frequency shift along offset, PFVO attempts to estimate a reliable peak frequency at zero offset, and at the same time to separately identify the frequency change caused by thin-bed scattering and that due to absorption. PFVO analysis operates on CMP gathers without NMO to avoid the effect of NMO stretch. If a non-stretch NMO technique (Perroud and Tygel, 2004; Masoomzadeh et al., 2004) were available, PFVO could be done after NMO correction. AVO analysis can be performed on an image gather from prestack time migration. Whether PFVO analysis will benefit from being applied to a common image gather is not clear at present. How prestack migration affects the frequency content also needs to be investigated. The accuracy of our method depends on the accuracy of time-frequency analysis and on the validity of the equation used to calculate Q. Further, we have assumed that frequency translation is caused by absorption and event tuning only. There may certainly also exist other factors which affect the frequency content of seismic data, in which case, such need be identified and quantified. The success of the attenuation analysis is heavily dependent on the accuracy of the interface information, generally represented by the reflectivity or the gradient of acoustic impedance. Without accurate subsurface boundary information, the attenuation property obtained from seismic data will be difficult to interpret. 51 Chapter 5 Seismic Absorption Compensation: a Least Squares Inverse Scheme This chapter addresses the problem of instability that plagues conventional inverse Q filtering. The de-absorption problem is formulated as an inverse problem in terms of least squares and regularization is imposed by means of Bayes’ Theorem. The solution is iterative and non-parametric, and it returns a reflectivity function that has been constrained to be sparse. The inverse scheme is tested on both synthetic and real data and the results obtained demonstrate the viability of the approach. 5.1 Introduction The aim of seismic processing is to obtain a high resolution image of the subsurface. Many steps are involved in the attempt to achieve this objective. One such step consists, or most certainly should consist, of the compensation for the ubiquitous absorption of seismic energy in the earth, a step referred to as de-absorption (McGinn and Duijndam, 1998). To a first approximation, the earth is a linear attenuator, the effect of which is to attenuate the higher frequencies of the seismic signal. The mechanisms and the relevant mathematical details have been described in Chapter 2 and in the literature mentioned there such as Futterman (1962), Kjartansson (1979), Aki and Richards (1980), Wang and Guo (2004) and Brillouin (1960). The most serious setback to the application of de-absorption is the inherent instability of deconvolution filters which, unless very carefully designed, unavoidably amplify high frequency noise. Many authors have considered this problem (Robinson, 1979; Hale, 1991; Hargreaves and Calvert, 1991; Bickel, 1993; Varela et al., 1993; Wang, 2002; Margrave et al., 2003; Wang, 2003, 2006). Although, as a result of the effort of the aforementioned authors, a variety of algorithms are now available to the industry, it is probably 52 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme safe to say that all are very sensitive to the presence of additive noise. The approach introduced here to this noise associated problem is by means of the application of least squares (LS) inversion. This approach attempts to obtain a solution to the inverse problem of de-absorption by minimizing the misfit between observed data and theoretically modeled data. Thus, the instability of applying an inverse operator to the random noise in the data, is avoided. The LS objective function may be constructed in various ways, but the manner which we prefer is based on a probabilistic rationale (Tarantola, 1987; Ulrych and Sacchi, 2005), where the measured and modeled data (which become the solution of the inverse problem) are considered as realizations of random processes. There are two main advantages to the LS approach proposed here as compared to previously published methods. One advantage is that, if the problem is formulated as under-determined, the type of its solution can be controlled as desired. The second advantage is that the LS approach uses only the forward operator, and does not require the unstable inverse operator. In other words, the formulation of de-absorption as an inverse problem, allows a stable solution to be achieved. As pointed out by Claerbout (1985), earth absorption might be compensated for by amplifying high-frequency energy during downward continuation, shifting a wavefield using an absorption operator. This very interesting topic of migration with absorption compensation will be investigated in Chapter 6. This chapter is arranged as follows: First, we discuss very briefly, critical issues in de-absorption, which include Q information, random noise and methods for absorption compensation. Secondly, we describe a deabsorption approach using a LS optimization scheme based on probabilistic inverse theory. Finally, we discuss the results of both synthetic and real data examples using the LS algorithm. Absorption in the earth causes both the attenuation and dispersion of the reflections recorded on the surface. Both phenomena lead to a loss of resolution with resulting loss of information concerning targets of potential interest. This is particularly severe for small , but hopefully profitable, targets at depth. In seismology, absorption is characterized by the quality factor, Q, and many different models have been developed which describe the associated attenuation-dispersion relationship. De-absorption, the process by which the undesirable effects of absorption are attenuated, is an inherently unstable process when applied in the conventional manner. The reason for this is that the model of absorption in the earth, customarily adopted, is one where the seismic trace is convolved with a time varying 53 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme filter. Conventional practice is to deconvolve the trace by means of the inverse of this filter. As is well known, deconvolution is very sensitive to high frequencies (Varela et al., 1993) and so, therefore, is de-absorption. The LS approach to de-absorption that is presented here differs quite radically from the deconvolution techniques in customary use. The problem is cast in terms of probabilistic inference where a priori information is used to shape the solution to our view of the earth. This view is that reflectivities are sparse. 5.2 5.2.1 Issues of Absorption Compensation Method Choice There are many methods available to implement absorption compensation, the most common ones are time-variant spectral whitening, time-variant deconvolution and inverse Q filtering (Yilmaz, 2000). Spectral whitening tries to compensate for absorption caused amplitude decay by applying a gain function at different frequency bands. Time-variant deconvolution implements absorption compensation by means of using a time-variant wavelet in a moving time window. Q information is not used explicitly, because Q information is hard to obtain. In order to compensate for attenuation more accurately and deterministically, Inverse Q filtering or other Q dependent compensation techniques are superior. 5.2.2 Inaccurate Q Estimation Because the estimated values of Q-factors are always different from real values, Q could be under or over estimated. Let us define the actual quality factor as Qo , its estimation as Qe and the difference as ∆Q, where Qe = Qo + ∆Q. (5.1) The inaccuracy ∆Q = 0 in the estimation will affect the amplitude and phase of the compensated signal in a manner outlined below. Effects on the Amplitude If Qe is used to compensate for the amplitude attenuation resulting from Qo (equation (2.20)), after the correction, the amplitude is Ac (ω) = exp −∆Q ω∆z . 2v(ω) (Qo + ∆Q)Qo 54 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme If absorption is under-estimated or ∆Q > 0, Q compensation does not recover the high frequency energy completely. If absorption is over-estimated or ∆Q < 0, Q compensation over-amplifies the high frequency energy. Figure 5.1 shows a test of using inaccurate Q’s in amplitude compensation on a zero-phased attenuation response. Figure 5.1(a) is the input zero-phased absorption response. Original reflection coefficients are three unit spikes at time t = 200 ms, t = 800 ms and t = 1600 ms respectively. Figure 5.1(b) is the result of amplitude compensation using actual Q. The compensation recovers the original spikes. Figure 5.1(c) and Figure 5.1(d) show the amplitude compensated responses using over-estimated and underestimated Q values. Larger Q values compensate part of the lost energy, and thus, the response gets sharper but is not as sharp as a spike. Smaller Q values over-compensate the high frequencies, and thus, the compensated result looks like strong high-frequency noise. Effect on the Phase If Qe is used to compensate for the phase distortion resulting from Qo (equation (2.21)), after the correction, the phase is φc (ω) = −∆Q ω∆z ln πv(ω) (Qo + ∆Q)Qo ω ω0 . If the actual value and estimated value are very close, i.e. |∆Q| ≈ 0, after Q compensation, φc (ω) ≈ 0, which means that the phase has no distortion. Pure phase compensation will make the original minimum-phase absorption response become zero-phase. This is an ideal case. If ∆Q > 0, i.e., Q is over-estimated or absorption is under-estimated, φc (ω) < 0. This means the arrival time of this frequency is delayed. The higher the frequency the more the wave is delayed. If ∆Q < 0, i.e., Q is under-estimated or absorption is over-estimated, φc (ω) > 0. This means the arrival time is advanced. The higher the frequency, the sooner the wave component arrives. Figure 5.2 shows phase corrected responses using under-estimated and over-estimated Q values respectively. The under-corrected signal is minimumphase, and the over-corrected signal is maximum phase. 5.2.3 Random Noise As already mentioned before, inverse Q filtering is inherently unstable since the inverse operator will boost high frequency noise. To ensure that noise is not unnecessarily amplified, it is important to design the inverse operator 55 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme 0.3 0.2 0.1 0 0 500 1000 1500 Time (ms) 2000 (a) Zero-phased impulse response with absorption. 0.8 0.6 0.4 0.2 0 0 500 1000 1500 Time (ms) 2000 (b) Amplitude compensation using actual Q = 30. 0.3 0.2 0.1 0 0 500 1000 1500 Time (ms) 2000 (c) Amplitude compensation using big Q = 40. 14 12 10 8 6 4 2 0 0 500 1000 1500 Time (ms) 2000 (d) Amplitude correction using small Q = 25. Figure 5.1: Experiments on amplitude compensation using different Q values. 56 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme 0.3 0.2 0.1 0 0 500 1000 1500 Time (ms) 2000 (a) Impulse response with absorption. 0.3 0.2 0.1 0 0 500 1000 1500 Time (ms) 2000 (b) Phase correction using Q = 30. 0.3 0.2 0.1 0 0 500 1000 1500 Time (ms) 2000 (c) Phase correction using Q = 20. 0.3 0.2 0.1 0 0 500 1000 1500 Time (ms) 2000 (d) Phase correction with Q = 40. Figure 5.2: A synthetic trace with three spikes, Q = 30, and its phase correction using different Q values. 57 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme 60 Amplitude 50 40 30 20 10 0 10 20 30 40 50 Frequency (Hz) 60 70 80 Figure 5.3: The amplitude of an ideal inverse Q filter (black) and its approximation (red). The two curves are overlapping before the fork. appropriately at high frequencies. One way to avoid the instability is to use a band-limited version of the inverse operator, i.e., to replace the amplitude compensation operator A−1 (ω) by W (ω)A−1 (ω), where W (ω) represents a low-pass filter. This idea is illustrated in Figure 5.3, where an exponential amplitude compensation operator is approximated by a Gaussian function after a certain frequency limit. Of course, as can be well imagined, such a low-pass operation causes an undesirable loss of resolution. The least squares method being introduced is a way to achieve complete compensation while still keeping the stability of the operation of inverse Q filtering. 5.3 5.3.1 Time-variant Deconvolution as an Inverse Problem Noise and Deconvolution The instability issue is not unique to the absorption compensation problem. It is, in fact, ubiquitous in the application of deconvolution in general. In seismic data processing, the recorded data in the one-dimensional case can be modeled as x(t) = w(t) ∗ r(t) + n(t), (5.2) 58 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme where x(t) , w(t), r(t), and n(t) represent the recorded trace, a seismic wavelet, the earth’s reflectivity and random noise, respectively. Seismic data deconvolution aims at obtaining the reflectivity from the observed trace. The estimated reflectivity will be adversely affected by the random noise term, whether or not w(t) is band limited. Defining f (t) as the inverse of w(t), assumed known (i.e. f (t) ∗ w(t) = δ(t) ), the deconvolution result is rˆ(t) = r(t) + f (t) ∗ n(t). Depending on the signal to noise ratio, the effect of the term f (t) ∗ n(t) on the output reflectivity can be very severe and must be compensated for. The required regularization always leads to a loss of resolution. To better handle the problem of additive noise, the inverse filter is often designed using a least squares (LS) approach rather than by directly inverting with the known or estimated seismic wavelet. The resulting filters are called optimum Wiener filters and are well elaborated by Robinson and Treitel (2002). We take a probabilistic approach to the LS solution here, an approach which is believed has a flexibility which leads to improved results. Specifically, we formulate de-absorption as an inverse problem and use a Cauchy-Gauss objective function (Sacchi et al., 1998) to arrive, iteratively, at the desired solution. 5.3.2 Bayes Probabilistic Inference Probabilistic inference is a very powerful approach to the ubiquitous inverse problem (Tarantola, 1987; Ulrych and Sacchi, 2005). It is common to consider the measured data d (where vectors are indicated by bold symbols) as uncertain. That is, true data do exist but are unknown. Measured data can then be considered as random variables whose expectation, in the ensemble sense, is the true value. Traditionally, the assumption is made that the observed data are random variables with a Gaussian distribution. This leads to the well known χ2 test for goodness of fit and to a l2 norm solution. Let us represent the solution to the inverse problem by the model vector m. m is not unique, in the sense that an infinity of such models may be found that fit the data. The reason for this, of course, is that the data, even if noise free, can never describe the continuous model completely. Our problem is under-determined. The problem may indeed be posed as an over-determined one, if we so choose. Wiener inversion is one such possibility. However, we believe, having been well schooled by the late Edwin Jaynes, that an inverse problem posed in such manner, is incorrectly posed (Jaynes, 2004). A discretized inverse problem can never be unique. The over-determined solution 59 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme returns the “smallest model”, i.e., that model whose energy is most distributed in model space. When estimating deconvolution models describing reflectivities of the earth, such models do not appear to be reasonable. An infinity of models is not a useful set. We must impose a priori knowledge, if such exists, feeling or at least hoping, associated with our search. The manner of so doing lies at the heart of Bayes’ theorem, and at the heart of the regularization thus imposed. This is hardly the place for a discussion of the underlying logic, but we feel that a few remarks are in order (for a fairly full discussion vis a` vis Bayes, inversion and a prior information, please see Ulrych and Sacchi (2005)). The central point in the application of Bayesian logic is in how one handles prior information. We begin with Bayes’ Theorem, which states that p(m|d) = p(d|m)p(m) , p(d) (5.3) where p(d|m) is the conditional probability density function(pdf) of d, the data, given that m has occurred. This pdf is called the likelihood and is the function that is maximized in l2 norm problems when in Gaussian form. The Gaussian pdf for a random variable x, that we will have occasion to use, is i2 h 1 − x−µ 2 2σ , (5.4) e p(x) = √ 2πσ where σ 2 is the variance and µ the mean of the pdf, respectively. Returning to equation (5.3), p(d) is the pdf of the observed data and often is not a factor in the inversion. In as much as d is considered to be a random vector, m can be treated in an equivalent manner. p(m|d) is the pdf that we wish to obtain. It is the probability of obtaining m given the data. Maximizing this pdf will obtain an estimated model that is referred to as the maximum a posteriori solution. If p(m) = 1, i.e., a uniform pdf, the formulation of the problem is exactly least squares. Now, however, comes a point of much discord, p(m). In the opinion of a group of statisticians known as “frequentists”, prior probabilities cannot be assigned, they must be measured. Edwin Jaynes, like Laplace before him, fought passionately against this philosophy, and most certainly emerged the victor. So, with victory on our side we press on. p(m) is the pdf that we wish to associate with our expected model. It is our choice based on our belief. It could certainly be inappropriate, but who is to say that the l2 model, the one most distributed in time and/or space, is any more so? What pdf to choose when dealing with de-absorption? We follow the logic of Sacchi et al. (1998) and use the Cauchy distribution, given by 60 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme p(x) = 1+ −1 x2 2σ 2 , (5.5) where σ 2 governs the spread of the distribution (the Cauchy pdf being one that does not possess a theoretical variance). The reason behind this choice is, very basically, that this pdf exhibits long tails and, consequently, is a candidate for the pdf that characterizes sparse reflectivities. Sparseness is of much relevance these days (Hermann, 2005; Cand`es et al., 2005). It is not a concept applicable in all cases, of course. One would not model a Beethoven concerto in this manner, but layers in the earth are, one hopes, sparsely distributed. 5.3.3 Absorption Compensation and Inversion We now formulate our problem in light of the above discussion. Assume the observed data to be contaminated by noise that is normally distributed as N (0, σn2 ), where n represents the noise vector. The Gaussian assumption stems from the principle of maximum entropy (Ulrych and Sacchi, 2005) for a discussion and references) which, informally speaking, states that if nothing is known, use the simplest hypothesis. In this case, the Gaussian pdf follows from the Central Limit Theorem. The conditional distribution of the data is given by (equation (5.3)) p(m|d, σn ) = 1 2πσn2 (M −1)/2 2 2 e−(1/2σn )||d−Gm)||2 , (5.6) where M is the length of the data vector and, for our linear system n = d − Gm, with G the coefficient or kernel matrix. Let p(m|σm ) indicate a prior distribution for m conditional on a parameter σm . From Bayes’ theorem, equation (5.3), we have p(m|d, σm , σn ) = p(d|m, σn )p(m|σm ) .. p(d) (5.7) The model vector, m, is the reflectivity function, and we use a sparseness constraint in the inversion scheme. 61 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme Cauchy-Gauss Model In the following, we will assume that a seismic wavelet is available. It is either obtained from check-shot survey, or from a seismic trace in some manner. We postulate, as reasoned above, that the model, i.e., the reflec2 (a standard tivity m, with elements mk of m i.i.d. with equal variance σm hypothesis). The joint pdf of the mk is p(m|σm ) = k p(mk |σm ), where p(mk |σm ) satisfies a Cauchy distribution given by equation (5.5). Substitution yields p(m|σm ) = 1+ k m2k 2 2σm . By inserting this Cauchy prior and the data likelihood (equation (5.4)) into equation (5.7) and taking logarithms of both sides, we obtain 1 ln[p(m|d, σm , σn )] = −c(m) − 2 (d − Gm)T (d − Gm), (5.8) 2σn where c(m) is a constraint imposed by the Cauchy distribution c(m) = M −1 k=0 ln 1 + m2k 2 2σm . Furthermore, denoting φcg (m) = − ln[p(m|d, σm , σn )], we obtain the cost function for the Cauchy-Gauss model as 1 φcg (m) = c(m) + 2 (d − Gm)T (d − Gm), (5.9) 2σn where G is composed of time variant wavelets, bτ (t − τ ). Both t and τ are in the range of the length of a trace and b0 (t − 0) b1 (t − 1) G= . .. . bM (t − M ) In this manner, one dimensional absorption compensation is formulated as an inverse problem. The model, which is related to the minimum of the cost function, is the sparse reflectivity function we desire. The CauchyGauss model has also been used in acoustic impedance inversion, signal interpolation and extrapolation (Sacchi et al., 1998). 62 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme Implementation The solution of the inverse problem, formulated above, follows the algorithm outlined by Sacchi et al. (1998). The objective function that will be minimized is φcg (m) in equation (5.9). Taking the derivative of φcg (m) with respect to m, we obtain ∂c(m) 1 ∂φcg (m) = + 2 GT (d − Gm) ∂m ∂m 2σn and (5.10) 1 ∂c(m) = 2 S−1 m, ∂m σm m2 where S is a M × M diagonal matrix with elements skk = 1 + σmk , (k = 0, 1, . . . , M − 1) and is m dependent. Equating the derivative in equation (5.10) to zero, yields m = λS−1 + GT G GT d, (5.11) 2σ where λ = mn22 . This equation is nonlinear and must be solved iteratively. To construct an iterative algorithm, equation (5.11) is written as m = SGT λIN + GSGT −1 d = SGT p. (5.12) p is the auxiliary vector which is obtained from the solution of the system λIN + GSGT p = d. (5.13) The iterative computation is initiated by setting the observed seismic data as the initial model, m0 , which is also used to generate the matrix S0 . In each iteration, k, we first compute p(k−1) = λIN + GS(k−1) GT −1 m, and then update the model as m(k) = S(k−1) GT p(k−1) . The algorithm, although iterative, is computationally efficient since the coefficient matrix in equation (5.13) is Hermitian-Toeplitz and is inverted using the Levinson recursion. The processing flow for absorption compensation on a stacked seismic section is implemented in the following five steps: 63 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme 1. Extract Q values by means of the WTVS approach. 2. Apply the required dispersive phase correction. 3. Use a given wavelet or extract a minimum phase wavelet from the early part of the trace. 4. Construct the kernel matrix G. 5. Solve the inverse problem to obtain the reflectivity series. Q and wavelet estimation need not necessarily be performed at every CMP location. Indeed, Q values need only be estimated for selected CMP’s and interpolated to form a Q profile. Computations required, other than those associated with the iterative algorithm, are trivial if both Q and the seismic wavelet do not change dramatically along the CMP direction, and the resulting LS scheme for seismic absorption compensation is computationally practical. Numerical Experiments This section illustrates results obtained using the LS approach on both synthetic and real data. The synthetic test is based on the assumption of known wavelet and Q. Figure 5.4 shows the convergence of the inversion. Panel (e) is the input trace with 20 percent random noise. Panels (a), (b), (c) and (d) are the de-absorption results after 1, 2, 3 and 4 iterations, respectively. After 2 or 3 iterations, the inverted reflectivities are very close to those in the actual model. Figure 5.5 imitates a section of seismic data. The signals in the five traces shown in Figure 5.5(a) are the same, but contain different random noise. The noise level is 20 percent. Figure 5.5(b) is the result obtained by LS timevariant deconvolution. The rightmost trace is the original reflectivity series plotted as a reference. It can be observed that the reflectivity function is recovered very well both as far as locations and amplitudes are concerned. This result shows that the Cauchy-Gauss objective function is a very viable a priori model for obtaining accurate reflectivity inversion. Figure 5.6(a) illustrates a stacked seismic section. The structures are not complex and consist of several flat layers. Processing of the section is initiated by first estimating the Q curve from a trace in the section using the WTVS method. Second, using these extracted Q’s, the dispersive phase correction is computed and is followed by the estimation of the minimum 64 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme Figure 5.4: Convergence process of the inverse scheme. (e) is a synthetic trace. (a-d) are the inverse results after iteration 1, 2, 3 and 4 respectively. 65 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme Figure 5.5: (a) is a synthetic trace with 20% noise. (b) is the result of absorption compensation. The rightmost trace is the original reflectivity series plotted as a reference. 66 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme (a) (b) Figure 5.6: Q compensation for a real seismic section. (a) shows a real seismic section. (b) is the result of absorption compensation. phase wavelet from the early times of a trace. The next step involves computations yielding the kernel matrix G. The inversion begins by setting a value for the parameter λ in equation (5.13) and the result is updated iteratively. Following inversion, the wavelet is convolved with the extracted reflectivities. The final results are shown in Figure 5.6(b). A number of reflections, which are not separated in the section shown in Figure 5.6(a), can be observed around time t = 920 ms after de-absorption and the lateral continuity can be tracked from trace to trace. The improvement in resolution is clear. The de-absorption result would certainly facilitate the interpretation of this seismic section. 67 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme 5.4 Conclusions The LS approach to de-absorption that is presented here differs quite radically from the deconvolution techniques in customary use. The main difference is that the inverse filter is designed using a Bayesian inference approach and is robust with respect to additive noise. The technique described here has very general application. Specifically, since robust Q compensation provides more accurate information concerning both the amplitude and location of the earth’s reflectivity, hydrocarbon reservoir characterization is one obvious target for the introduced method. In particular, the resulting improved resolution may have important consequences in 4D processing where the objective is to delineate changes in reservoir properties. Results that we have obtained on synthetic and real data confirm, we believe, both the sparse assumption that we have made, and the viability of the Cauchy-Gauss prior model that is used to impose it. By taking the advantage of the Hermitian-Toeplitz property of the kernel matrix that describes the forward problem, stable inverse solutions are achieved economically. The synthetic and real data results demonstrate the viability of the approach and the pivotal role of de-absorption in improving the resolution of the processed sections. The LS Q compensation is implemented in a trace by trace manner and does not take into account lateral continuity or the actual ray-paths that are involved. Multi-channel Q compensation and poststack time migration with Q compensation are interesting and important topics which will be discussed in Chapter 6. 68 Chapter 6 Refocusing Migrated Images in Absorptive Media This chapter investigates the blurring effect on migrated images when using a regular migration algorithm to migrate the seismic data with absorption. First, a numerical method is introduced to calculate the blurring functions in layered media, and then, a least squares inverse scheme is used to remove the absorption blurring effect in order to refocus the migrated image. The stability and the convergence property of the refocusing algorithm will be discussed. Experiments on synthetic and real data will be given to show the validity of the introduced technique to compensate for absorption after migration. 6.1 Introduction Conventional Q-compensation is implemented in a trace by trace mode using, typically, time-variant spectral whitening, time-variant deconvolution (Yilmaz, 2000), least squares inversion (Zhang and Ulrych, 2007), onedimensional wavefield extrapolation (Hargreaves and Calvert, 1991) and many others. In reality, the spatial distribution of a Q field is similar to that of a velocity field, i.e., Q is spatially variant. The effect of attenuation on received seismic signals is related to wave-paths, so that absorption compensation is ideally implemented based on the wave equation in the process of seismic imaging. Migration with absorption compensation is referred to as Q-migration. Both prestack Q-migration and poststack Q-migration have been investigated by some researchers. Mittet et al. (1995) proposed a prestack depth migration method for accommodating absorption and dispersion effects in a prestack depth migration scheme. Their method was presented in the frequency-wavenumber domain using a standard linear solid model. An extrapolation operator that compensates for absorption and dispersion was designed using an optimization algorithm, which was then used to devise 69 Chapter 6. Refocusing Migrated Images in Absorptive Media a finite-difference migration scheme. Yu et al. (2002) applied poststack Q-migration to enhance the seismic frequencies on the top of an anticline structure to locate a gas cap, and Wang (2008) presented his tests of a poststack Q-migration technique on a synthetic dataset of the Marmousi model (Versteeg, 1994). The main difference between Q-migration and regular migration is that the amplitude of the wavefield continuation operator used in Q-migration is greater than 1, while the amplitude of the continuation operator used in regular migration is equal to or less than 1. In practice, a Q-migration algorithm must deal with the instability problem caused by the amplification of high frequency random noise in its implementation. Conventionally, if a Q-migration operation is based on wavefield extrapolation, a noise suppression step is required in the wavefield continuation process in order to stabilize the computation. Another option is to use least squares migration (LSM) (Nemeth et al., 1999; Zhang, 2000). A LSM scheme for wavefield Q-compensation may be useful in stabilizing the solution, but to obtain a solution iteratively is computationally expensive. Traditional trace by trace inverse Q filtering neglects the 3D features of wavefield propagation, and regular Q-migration suffers from the amplification of high frequency random noise. This chapter aims to find another approach to do wavefield Q-compensation, which can be considered as a multi-dimensional Q compensation technique. This multi-dimensional Qcompensation problem will be handled as a problem of multi-dimensional deconvolution like other blurring problems in seismic migration. The blurring problem in seismic imaging includes coarse spatial sampling, inaccurate migration operators, and limited acquisition aperture, all of which lead to a blurred migration image. The blurring function is also termed as a point scatter function (PSF) (Jansson, 1997), an illumination function (Xie et al., 2006), a migration Green’s function (Schuster and Hu, 2000), Hessian or a migration resolution function (Rickett, 2003; Guitton, 2004; Toxopeus et al., 2004), etc.. The blurring function has been well investigated in recent years in order to obtain more accurate images of the earth’s structures. While many efforts are focused on trying to design accurate migration operators for the migration algorithm, some issues which affect the quality of seismic imaging can be handled after migration. Schuster and Hu (2000), Hu et al. (2001) and Yu et al. (2006) developed a thread of techniques called migration deconvolution to handle the problem of acquisition footprint in prestack and poststack migrated images. The key in migration deconvolution techniques is how to calculate the blurring function. The next 70 Chapter 6. Refocusing Migrated Images in Absorptive Media important step is how to remove it. This chapter borrows the idea of migration deconvolution to handle the blurring problem caused by using a regular migration operator to migrate attenuated seismic data. In the following, we will discuss the attenuation effect on the migration blurring function. Removing the blurring influence from seismic images is a problem of a 2D or 3D non-stationary deconvolution which can be solved using a least squares inverse scheme. The stability and efficiency of the least squares solver will be discussed. 6.2 Blurring Function caused by Seismic Absorption According to Claerbout (1992), operations of seismic migration and modeling are conjugate pairs. Assuming data acquisition is a modeling process L, the migration process is equivalent to applying the conjugate LT to the recorded data Lm0 , i.e. m = LT Lm0 , (6.1) where m0 is the true image or the real earth structure, and m is the migrated image, or the blurred earth structure, because LT is not the exact inverse of L. For seismic waves traveling in absorptive media, mathematically, a migration blurring function can be defined as Γ(r|r0 ) = LT Lq δ(r0 ), (6.2) where δ(r0 ) stands for a delta function at location r0 and Lq stands for an operator of wave propagation in an absorptive medium. Equation (6.2) implies that when regular migration is applied to the attenuated seismic data, the blurring effect caused by absorption is mixed with other blurring effects in regular migration as those mentioned before. A direct solution to obtain the actual image from equation (6.1) is m0 = LT L −1 m. However, the inverse of LT L is generally difficult to calculate directly. One way to obtain an approximation of the inverse is to use a conjugate gradient algorithm to solve a least squares problem (Nemeth et al., 1999; Sacchi et al., 2006). The same statement applies to the operator LT Lq in equation (6.2). 71 Chapter 6. Refocusing Migrated Images in Absorptive Media 6.2.1 Analytical Analysis Assuming that there is a monochromatic source of angular frequency ω starting from a point scatterer located at r0 in a homogeneous medium with velocity v and quality factor Q, the wavefield at a receiver location rg is the solution of the Helmholtz Green’s function (Schuster and Hu, 2000), G(rg |r0 , ω) = ω |rg −r0 | −j c(ω) e |rg − r0 | . ω In this expression, c(ω) is a complex wavenumber. Absorption caused wavefield attenuation and dispersion are represented through the complex velocity c(ω). Its expression will be given in equation (6.7). The complex conjugate of the Green’s function is the operator used in migration which is equivalent to the wavefield of an impulse starting from a receiver rg down to an image point r, j ω |r −r| e v0 g . G(r|rg , ω) = G (rg |r, ω) = |rg − r| ∗ The migration blurring function in an absorptive medium, if the migration is implemented without absorption compensation, is ω |rg −r0 | j vω |rg −r| −j c(ω) e Γ(r|r0 , ω) = |rg − r0 | rg e 0 drg . |rg − r| (6.3) The blurring function in the time domain is the integral of all frequency components ω −j c(ω) |rg −r0 | j vω |rg −r| e Γ(r|r0 , t) = rg |rg − r0 | e 0 drg dω. |rg − r| (6.4) The closed form of this integral is difficult to derive. However, Schuster and Hu (2000) derived closed forms of both 2D and 3D blurring functions, for acoustic media. For an acquisition with a large aperture, it is assumed that the receivers are symmetrically distributed above a point scatterer. The main energy of a migration blurring function is localized in the shadow area as that shown in Figure 6.1 for homogeneous media. 6.2.2 Numerical Computation Algorithm Choice If an analytical function were derivable, at most, it could be applied to homogeneous media. For layered Q and/or v media, it is difficult to calculate 72 Chapter 6. Refocusing Migrated Images in Absorptive Media Figure 6.1: Diagram of the energy distribution of a migration blurring function. Symbol “V” represents a receiver location. a migration function analytically, but a migration blurring function can be calculated numerically. The calculation is implemented as follows: first compute synthetic data with attenuation Lq δ(r) and then migrate the synthetic to obtain the migration blurring function LT Lq δ(r). The synthetic data can be computed using phase-shift modeling, or using ray-tracing plus wavelet continuation along the ray-paths. Phase-shift modeling is based on the following wavefield extrapolation operator √ 2 2 (6.5) H(ω, z) = e−jz ka −kx , and assumes an impulse begins at a point scatterer. In equation (6.5), ka is a complex wavenumber (refer to equations (2.13 and 2.15)) , and ka = or ka = ω c(ω) , 1 ω ln 1− v0 πQ ω ω0 1+j 1 , Q (6.6) with 1 1 1 = ln 1− c(ω) v0 πQ ω ω0 1+j 1 . Q (6.7) Phase-shift modeling is easy to compute but the synthetic always contains Fourier wrap-around noise. Ray-tracing synthetic computation is expensive due to the cost of ray-bending ray-tracing (Cerveny, 2005). Once a 73 Chapter 6. Refocusing Migrated Images in Absorptive Media ray-path is defined in the media, the waveform at an arrival time can be calculated through extrapolating a wavelet along the ray-path. If the purpose of Q-compensation is purely to enhance the frequency without considering velocity dispersion, a wavelet can be added to a record at an arrival time based on the following expression (equation (3.3)) − B(ω, t2 ) = B(ω, t1 )e ω(t2 −t1 ) 4πQ , where B(ω, t1 ) and B(ω, t2 ) are the wavelet spectra at two different time locations t1 and t2 . Q is the quality factor in the layer between t1 and t2 . Ray-tracing modeling and corresponding Kirchhoff migration (Berkhout, 1981) can apply to different velocity and Q models. Phase-shift modeling and phase-shift migration (Gazdag, 1978) only apply to vertically variant media. The migration algorithm is the adjoint of the modeling algorithm. Therefore, when computing a migration blurring function, if phase-shift migration is used for the imaging, the modeling algorithm is better to use the phase-shift method. If Kirchhoff migration is used for the imaging, correspondingly, the modeling algorithm is better to use ray-tracing. Fast Migration Blurring Function Updating When the phase-shift method is used to calculate the migration blurring function, there exists a fast algorithm for extrapolating a blurring function from one depth step to another depth step or from one time step to another time step. For vertically layered media, vi is the reference velocity of layer i, and ci (ω) is complex velocity of layer i including the effect of seismic absorption. If there is a unit impulse at a depth z0 = n∆z, and its response on the surface is n j∆z e U (kx , z = 0, ω|z0 ) = r ω2 −kx2 ci (ω)2 , i=0 and the response of a unit impulse right below it at depth z0 + ∆z is j∆z U (kx , z = 0, ω|z0 + ∆z) = e r ω2 −kx2 cn+1 (ω)2 U (kx , z = 0, ω|z0 ). To migrate the impulse response U (kx , z = 0, ω|z0 ) from the surface to an image point z = m∆z using downward wavefield continuation, the blurring function is m −j∆z e Γ(kx , z, ω|z0 ) = U (kx , z = 0, ω|z0 ) r ω2 −kx2 v2 i . (6.8) i=0 74 Chapter 6. Refocusing Migrated Images in Absorptive Media To migrate the impulse response U (kx , z = 0, ω|z0 + ∆z) from the surface to the same image point z = m∆z, the blurring function is m −j∆z e Γ(kx , z, ω|z0 + ∆z) = U (kx , z = 0, ω|z0 + ∆z) r ω2 −kx2 v2 i . (6.9) i=0 Examining equations (6.8) and (6.9), we can observe that, in layered media, the blurring function can be derived from one depth to another depth without going through the full process of forward modeling and migration, i.e. r j∆z Γ(kx , z, ω|z0 + ∆z) = e ω2 −kx2 cn+1 (ω)2 Γ(kx , z, ω|z0 ). (6.10) This observation helps to save the calculations of the blurring functions. Figure 6.2 shows two blurring functions in an absorptive medium at two different time locations. The main energy of a blurring function lies in an area, whose shape is like that sketched in Figure 6.1. As time increases, the energy distribution of the migration blurring function extends to a larger area. For layered media, if velocity v and attenuation property Q do not change laterally, the migration blurring function needs only to be computed at one reference location. To allow for lateral v and Q variation, practically we can divide the image into sub-regions or patches (Hu et al., 2001; Yu et al., 2006), and compute a migration blurring function for each sub-region, assuming that in each patch, the lateral variations of velocity and Q are very gentle. 6.3 Removing Absorption Blurring Function using a Least Squares Method As discussed in the foregoing sections, the migration blurring function Γ = LT Lq is not an identity matrix. In our case, the migrated image is a blurred image due to using an inaccurate migration operator, m = LT Lq m0 . (6.11) In order to compensate for absorption after migration to get a solution closer to the right image m0 , the blurring function needs to be removed or at least to be shrunk. Generally, the inverse of LT Lq cannot be solved directly. There are examples of using the diagonal to approximate the effects of the matrix on the image amplitude (Sacchi and Wang, 2007; Hu et al., 2001). In fact, 75 Chapter 6. Refocusing Migrated Images in Absorptive Media Time (s) 1.0 0 1.2 Midpoint (km) 1.4 1.6 1.8 0.5 1.0 (a) Time (s) 1.0 0 1.2 Midpoint (km) 1.4 1.6 1.8 0.5 1.0 (b) Figure 6.2: Examples of the blurring functions in an absorptive medium. (a) at time t = 400 ms. (b) at time t = 700 ms. 76 Chapter 6. Refocusing Migrated Images in Absorptive Media the migrated image can be considered as the real image convolved with the migration blurring function. The blurring function is multi-dimensional (2D or 3D) and non-stationary in time, i.e. equation (6.11) can be written in a form of convolution m(x, t) = γxref ,tj (x, t) ∗md m0 (x, t) + n(x, t). (6.12) In equation (6.12), symbol ∗md represents md dimensional convolution, md is either equal to 2 or 3. x represents 1 or 2 lateral spatial coordinates, and t could be z for depth imaging. γxref ,tj (x, t) represents the migration blurring function at a reference location (xref , tj ). n(x, t) is random noise. The inverse problem could be unstable because of two reasons: the inverse of γ is unstable or even does not exist, and the inverse filtering could raise the energy of random noise. Solving equation (6.12) for the real image m0 (x, t) is a problem of multidimensional non-stationary deconvolution. Time-variant deconvolution can be formulated as a least squares inverse problem (Zhang and Ulrych, 2007). Similarly, the image refocusing process can be solved using the least squares inversion in the wavenumber domain. To expedite the computations, we assume that γ(r|r0 ) is shift invariant in the horizontal coordinates. After Fourier transform on the spatial axes, equation (6.12) becomes ˆ 0 (k, t) + n ˆ (k, t). m( ˆ k, t) = γˆxref ,tj (k, t) ∗ m (6.13) Equation (6.13) changes the multi-dimensional convolution problem to a one dimensional time-variant convolution. Accordingly, the multi-dimensional deconvolution can be solved as a one-dimensional time-variant deconvolution at each wavenumber. Similar to that discussed in Zhang and Ulrych (2007), the time-variant deconvolution can be formulated as a least square problem using Bayes’ theory. The objective function to be minimized is, φ(m ˆ 0 ) = c(m ˆ 0) + 1 ˆ ˆ − Gm (d − Gm ˆ 0 )T (d ˆ 0 ). 2σn2 (6.14) In this equation, the first term c(m ˆ 0 ) on the right-hand side is a model regularizer, and the second term is the misfit which could be interpreted as ˆ represent model and data Gaussian noise in the observed data. m ˆ 0 and d in the wavenumber domain respectively. The objective function can be minimized through solving the following equation, −1 ˆ GT d, (6.15) m ˆ 0 = GT G + λW−1 77 Chapter 6. Refocusing Migrated Images in Absorptive Media where G is the kernel matrix of the inverse problem, W is a regularizer for the image at a wavenumber and λ is a hyper-parameter. The kernel matrix G is formed from a serials of time-variant complex γˆ functions, and γˆ0 (t − 0∆t) γˆ1 (t − 1∆t) .. . (6.16) G= γˆi (t − i∆t) . .. . γˆN (t − N ∆t) In equation (6.16), N is the number of trace samples. γˆi (t − i∆t) represents a migration blurring function in the wavenumber domain at time i∆t with a time shift of i∆t. The weights at wavenumber kj are determined by the weights and the obtained image at lower wavenumber kj−1 Wii [kj ] = αs(mi [kj−1 ]) + (1 − α)Wii [kj−1 ], (6.17) where α ∈ [0, 1] is a “viscosity” parameter (Hugonnet et al., 2001) which means the similarity of energy distribution from low wavenumber to high wavenumber. s(mi ) is a function of model parameters, s(mi ) = (|mi | + ǫ)/(max |mi | + ǫ). i This model constraint implies that the inversion favors strong energy. If the weights are equal, the solution is standard least squares (Robinson and Treitel, 2002). Although the theoretical soundness of this constraint still needs to be further investigated, it works for the two data examples which will be discussed in the following sections. A pseudo code in C language of the refocusing algorithm for 2D migrated time section is listed in table (6.3). The main computation of the refocusing process has two parts: one is the calculation of the blurring functions, another is to solve for m0 using equation (6.15). To extend the algorithm from 2D to 3D is very straightforward. The code needs to be modified to handle one extra spatial dimension y. For phase-shift migration, the blurring function is easy to compute using equation (6.10). For Kirchhoff migration, the computation of the migration blurring function relies on the speed of ray tracing. The algorithm is implemented wavenumber by wavenumber which is easy to be parallelized for distributed computations. In each wavenumber the main computation is matrix multiplication of size (N × N ) which is substantial but not excessive. Therefore, the implementation is practical. 78 Chapter 6. Refocusing Migrated Images in Absorptive Media /* part 1 calculate the blurring function */ { do while (time<maximum_time); { calculate_synthetic (); do_migration (); do_fft_x2kx (a_blurring_function); keep_blurring_function (); time += time_step; } } /* part 2 deconvolution using least squares inversion */ { do_fft_x2kx (a_section); for (ikx=0, ikx<maximum_kx; ikx++) { form_kernel_matrix (); solve_for_m0 (); } do_fft_kx2x (kx_section); } Table 6.1: A pseudo code of the refocusing algorithm. 79 Chapter 6. Refocusing Migrated Images in Absorptive Media Time Index 0 0 Time Index 200 400 200 400 Figure 6.3: Amplitude of a kernel matrix G. The algorithm is stable because GT G is a diagonally dominant Hermitian matrix. The diagonal elements are close to its eigenvalues because the energy of a blurring function is mainly concentrated in a time window centered at the point scatterer. The magnitudes of the elements decay along the diagonal due to amplitude attenuation, but they are greater than zero, which means the kernel matrix is not singular. Figure 6.3 shows an example of the element amplitude of this matrix. The ripples at the top right corner and the bottom left corner are Fourier wrap-around noise in phase shift modeling and migration. They need to be muted out in real data processing, because they do not exist in observed real data. G can also be considered as a band-shaped matrix because of the localization of blurring functions. By considering this characteristic, the computation required by matrix manipulation could be greatly reduced. There are also efforts made to approximate G by a diagonal matrix in order to save computation with the trade-off with accuracy (Sacchi and Wang, 2007). Here, the solution is obtained by solving equation (6.15) using the full matrix of G. Theoretically, the refocusing technique used here is only valid for horizontally layered media. Consequently, it is better to be applied to refocusing a migrated time image without strong lateral v and/or Q variation. Otherwise, the migrated image needs to be processed in patches based on v and/or 80 Chapter 6. Refocusing Migrated Images in Absorptive Media Q variation. These patches are, regularly, overlapping rectangles. In each rectangular area, the refocusing operation is applied independently. Images in overlapping areas are the mathematical average of those involved patches. The final result will depend on the accuracy of the computed blurring functions, the inversion algorithm, the regularization parameter and also patches being divided. 6.4 Data Experiments The refocusing algorithm is first tested on a synthetic model where its migration blurring functions are known. Figure 6.4(a) shows a sine-function shaped structure obtained in an absorptive medium. Figure 6.4(b) shows the refocusing result. Because the refocusing process compensates for the attenuation effect, the resolution of the obtained image is enhanced. The algorithm is also tested on a 2-D field data-set. Its migration blurring functions are computed based on the velocity and Q fields. The velocity field is shown in Figure 6.5(a) and its inverse Q field is in Figure 6.5(b). The inverse Q field is obtained using the procedure described in table 3.3.1, which shows that the major absorption is before 1.4 s with less attenuation afterwards. This absorption profile is generally in line with geologic structures at different environments, because rocks at deep places are relatively harder and better compacted. The refocusing process does not require the Q information, but Q values are required when doing the forward modeling to calculate the absorption blurring functions. Both velocity and Q fields are laterally smoothed and the values at the center CMP (400) location are used to calculate the blurring functions at each time sample. Figure 6.6(a) shows a local area from the image obtained by time-variant spectral whitening (TVSW) plus phase-shift migration. Figure 6.6(b) shows the corresponding area from the image obtained using phase-shift migration plus a migration refocusing process. The section obtained by the refocusing process shows that the algorithm shrinks migration noise and increases the seismic resolution. The overall quality of images has been enhanced. For example, those small vertical faults at time 2.25 s are much clearer in Figure 6.6(b) than that in Figure 6.6(a). The comparison shows the difference between 1D absorption compensation and 2D absorption compensation. The resolution improvement of Figure 6.6(b) compared to Figure 6.6(a) shows the validity of the refocusing algorithm in compensating for the absorption after migration. If the image were processed in patches, the quality of the refocused image would be, hopefully, further improved. 81 Chapter 6. Refocusing Migrated Images in Absorptive Media 0 0 0.2 Midpoint (km) 0.4 0.6 Time (s) 0.5 1.0 1.5 2.0 (a) 0 0 0.2 Midpoint (km) 0.4 0.6 Time (s) 0.5 1.0 1.5 2.0 (b) Figure 6.4: A synthetic example. (a) shows an absorption-blurred sinefunction shaped structure. (b) is the refocusing result. 82 Chapter 6. Refocusing Migrated Images in Absorptive Media (a) (b) Figure 6.5: Input information. (a) is the velocity section and (b) is the inverse Q section. 83 Chapter 6. Refocusing Migrated Images in Absorptive Media 300 350 400 CMP 450 500 550 500 550 Time (s) 1.5 2.0 2.5 (a) 300 350 400 CMP 450 Time (s) 1.5 2.0 2.5 (b) Figure 6.6: Migration results. (a) Time-variant spectral whitening plus phase-shift migration. (b) Phase-shift migration plus refocusing. 84 Chapter 6. Refocusing Migrated Images in Absorptive Media 6.5 Conclusion and Discussion To obtain the real seismic image, the Hessian matrix (LT L) needs to be close to an identity matrix, or the migration blurring function should be a spatial delta function in 2D or 3D domain. Schuster and Hu (2000), Hu et al. (2001) and Yu et al. (2006) investigated the acquisition footprint, or nonuniform illumination resulting from limited recording aperture by means of migration Green’s functions. Here, we have put aside the acquisition geometry blurring function and only discuss the blurring effect caused by an inaccurate migration operator when migrating data without considering the absorption. The numerical implementation of refocusing a migrated image in an absorptive medium includes the following three steps: 1. Make synthetic data with attenuation. 2. Migrate the synthetic data without attenuation to get blurring functions. 3. Remove absorption blurring functions from the image using least squares inversion. Compared to the Q-migration algorithm (Mittet et al., 1995; Yu et al., 2002), the refocusing process is able to address the problem of instability on high frequency noise. Compared to the least squares migration algorithm (Nemeth et al., 1999; Sacchi et al., 2006), the refocusing algorithm is cost effective in implementation. When computation is not a big issue, the refocusing process can be applied to process common offset sections. In doing so, this method can be extended to prestack seismic processing. If the blurring function includes the effect of a seismic wavelet, the refocusing process is equivalent to poststack multi-dimensional deconvolution. In practice, the migration blurring function is related to many factors, such as acquisition geometry, seismic wavelet, velocity structure and absorption properties. This chapter blends the absorption effect with other blurring effects into a comprehensive migration blurring function, and multichannel Q-compensation is achieved in the process to refocus the migrated image by time-variant deconvolution in the wavenumber domain. 85 Chapter 7 Discussion and Conclusions The results of this research provide new techniques for handling signal absorption in seismic data imaging and inversion. The following paragraphs delineate the specific contributions contained in this thesis. The introduced analytical formula (equation(3.13)) relates Q-factor to spectral peak frequency variation. The assumption in the theoretical derivation is that the source wavelet has an amplitude spectrum similar to that of a Ricker wavelet (equation(3.1)). Although an actual amplitude spectrum may not, of course, exactly conform to this assumption, in many cases, the amplitude spectrum can still be well approximated by a Ricker spectrum. This implies that this assumption is not rigorously required. This approach can be used to estimate Q from CMP gathers, from a post-stack trace, and from seismic data acquired from other geometries like VSP (vertical seismic profiling) as well, as long as the peak frequency variation can be reliably examined. To estimate Q-factor from a prestack CMP gather, variation of a wavelet spectrum of an event is analyzed along offset. To obtain spectral peak frequency variation from a poststack trace, we can use windowed time-variant spectral analysis or a continuous wavelet transform analysis. In a WTVS, controlling points can be picked in a way similar to velocity analysis (Figure 3.4). These controlling points divide a trace into layers and their coordinates are used to calculate the quality factors of each layer. In two respects, the analytical method is superior to the spectral ratio method. First, the analytical method only uses the information of frequency variation. Besides absorption, it is not affected by other amplitude factors such as reflection, transmission and geometrical spreading. Second, the analytical method is robust to random noise. For the two-dimensional case, the robustness is achieved by means of using multi-fold information at different offsets. For the one-dimensional case, the robustness is achieved by means of using the frequency trend analysis. Event-based Q estimation can provide the absorption information for Q compensation. For reservoir description, seismic attenuation property, Qfactor, is ideally estimated at all prospective layers defined by velocities or 86 Chapter 7. Discussion and Conclusions acoustic impedance. Techniques to analyze seismic attenuation properties for hydrocarbon reservoir description generally include iso-frequency analysis, direct Q estimation and attenuation or Q tomography. The introduced reflectivity-guided seismic attenuation analysis is potentially another useful technique. It operates on prestack CMP gathers, and uses the location information of structure interfaces to select layer boundaries. The interface information could be obtained from poststack acoustic impedance inversion or well logs. This technique basically provides a methodology to allow Q values to be estimated at all prospective layers. The obtained Q section can have the same resolution as other attribute sections. To estimate spectral peak frequency variation at zero offset, an AVO-like, peak frequency estimation scheme called PFVO which operates on CMP gathers is introduced. The advantage of PFVO over post-stack frequency analysis is that it avoids the distortion caused by seismic stacking on frequency characteristics. PFVO tries to estimate a reliable peak frequency at zero offset, and at the same time to identify the frequency change caused by real attenuation and that caused by event tuning. Absorption compensation or de-absorption is an important step in seismic data processing; it helps to improve the resolution of seismic data, and to balance the frequency contents. In order to stabilize traditional inverse Q filtering, the inverse Q filtering is formulated as a least squares inverse problem based on Bayes’ theory. The prior information imposed on the model parameters is the sparseness of reflectivity, which is represented by a Cauchy model, and random noise in observed data is assumed to be a Gaussian. The inverse problem is solved iteratively. Experiments performed on both synthetic and real data show that the inverse scheme can produce ideal results for absorption compensation. To preserve both the amplitude and waveform of seismic signals during imaging, it is necessary to incorporate de-absorption and migration together. This may be implemented through standard migration with random noise attenuation at each continuation step, or through least squares migration which is implemented by casting the standard migration problem as an optimization problem and solving the optimization problem iteratively. Based on the concept of migration deconvolution, a post-migration refocusing scheme is introduced which could be considered as an alternative to conventional migration with Q compensation. The refocusing scheme can be used to remove the absorption effect on seismic imaging, at the same time to remove the other migration artifacts. Compared to regular migration with Q compensation, the refocusing method handles the instability problem caused by the amplification of the 87 Chapter 7. Discussion and Conclusions high frequency noise by means of inversion. Compared to the least squares migration, the refocusing method is cheaper in implementation. The refocusing method can be applied to common offset sections by using different operators at different offsets. In doing so, this method can be extended to prestack seismic processing. If the blurring function includes the effect of seismic wavelet, the refocusing process is equivalent to multidimensional deconvolution. The work presented in this thesis cover just two areas of seismic absorption research: seismic absorption estimation and compensation. Q-factor estimation from seismic data is the first step to use attenuation properties for seismic data imaging and inversion. To confirm the accuracy of the estimated Q values, numerical modeling, such as reflectivity modeling (Kennett, 1983) or waveform modeling (Carcione, 1995), may be required to check whether synthetic records match observations. For seismic modeling, a mathematical model is expected to describe wave behavior in absorptive media. Generally, a Q model (see equations (2.13) and (2.15)) is enough to describe the absorption phenomenon in visco-elastic media. To use an estimated Q section as an attribute to invert for reservoir properties not only requires accurate Q information, but also depends on reliable rock physics models like those suggested by Biot (1956a), Biot (1956b), Pride et al. (2004) and Dvorkin and Mavko (2006). A rock physics model to relate Q with permeability, porosity, saturation and fluid properties etc. is ideally expected. Rock physical properties cannot be inverted using Q-factor information only. However, integrated inversion of using absorption property, AVO, velocity and well log information is possibly one practical way for obtaining useful parameters for hydrocarbon reservoir description. All aforementioned seismic information could be of P-waves, S-waves or both. Recent developments in multi-component seismic wave acquisition make it possible to obtain S-wave velocities, S-wave AVO effects, and Swave absorption properties as well. There are reports showing applications of using the relation between P-wave and S-wave Q factors as an indicator to identify brine saturated and gas filled layers in seismic sections (Mavko et al., 2005a; Bale and Stewart, 2002). Integrated inversion of using both P-wave and S-wave absorption properties, velocities, AVO effects and other geophysical information for rock properties could be an interesting topic for future research. Accurate Q estimation from CMP gathers provides critical information for prestack Q compensation. With quantitative attenuation information, it may be possible to separate AVO effects from seismic absorption in CMP super-gathers. That is to say, AVO information extracted from seismic data 88 Chapter 7. Discussion and Conclusions would be more reliable if absorption effects on AVO had been taken into account in the analysis accurately. The attenuation effect on AVO needs to be further investigated. Absorption compensation can be performed using the introduced refocusing method after post-stack migration. It can also be implemented on prestack common offset sections. More research is expected to investigate the difference between poststack and prestack refocusing operations when underground structures have strong lateral velocity and Q variations. 89 Bibliography Aki, K. and P. G. Richards, 1980, Quantitative seismology-theory and methods: W. H. Freeman and Company. Bachrach, R., R. Ferber, A. Salama, N. Banik, R. Utech, D. Keller, and M. Bengtson, 2006, Two approaches for Q estimation and its application to seismic inversion: 76th SEG, expanded abstract, 2971–2974, Soc. of Expl. Geophys. Bale, R. and R. Stewart, 2002, The impact of attenuation on the resolution of multi-component seismic data: 72nd SEG, expanded abstract, 1034– 1037, Soc. of Expl. Geophys. Batzle, M., R. Hofmann, M. Prasad, G. Kumar, L. Duranti, and D. Han, 2005, Seismic attenuation: observations and mechanisms: 75th SEG, expanded abstract, 1565–1568, Soc. of Expl. Geophys. Berkhout, A. J., 1981, Wave field extrapolation techniques in seismic migration, a tutorial: Geophysics, 46, 1638–1656. Bickel, S. H., 1993, Similarity and the inverse Q filter: The Pareto-Levy stretch: Geophysics, 58, 1629–1633. Biot, M. A., 1956a, Theory of propagation of elastic waves in a fluidsaturated porous solid. I. low-frequency range: J. Acous. Soc. Am., 28, 168–178. ——–, 1956b, Theory of propagation of elastic waves in a fluid-saturated porous solid. II. high-frequency range: J. Acous. Soc. Am., 28, 179–191. Blanch, J., J. Robertsson, and W. Symes, 1995, Modeling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique: Geophysics, 60, 176–184. Bohlen, T., 2002, Parallel 3-D viscoelastic finite-difference seismic modeling: Computers and Geosciences, 28, 887–899. Brillouin, L., 1960, Wave propagation and group velocity: Academic Press. 90 Bibliography Buckingham, M. J., 2005, Causality, Stokes’s wave equation, and acoustic pulse propagation in a viscous media: Physical Review, E72, 026610–1– 026610–8. Cand`es, E., J. Romberg, and T. Tao, 2005, Stable signal recovery from incomplete and inaccurate measurements: Comm. Pure Appl. Math., 59, 1207–1223. Carcione, J., D. Kosloff, and R. Kosloff, 1988, Viscoacoustic wave propagation simulation in the earth: Geophysics, 53, 769–777. Carcione, J. M., 1993, Seismic modeling in viscoelastic media: Geophysics, 58, 110–120. ——–, 1995, Constitutive model and wave equations for linear, viscoelastic, anisotropic media: Geophysics, 60, 537–548. Carcione, J. M. and S. Picotti, 2006, P-wave seismic attenuation by slowwave diffusion: Effets of inhomogeneous rock properties: Geophysics, 71, O1–O8. Castagna, J. P. and S. Sun, 2006, Comparison of spectral decomposition methods: First Break, 24, 75–79. Castagna, J. P., S. Sun, and R. W. Siegfried, 2003, Instantaneous spectral analysis: Detection of low-frequency shadows associated with hydrocarbons: The Leading Edge, 22, 120–127. Cerveny, V., 2005, Seismic ray theory: Cambridge Univ. Press. Chapman, M., E. Liu, and X. Li, 2005, The influence of abnormally high reservoir attenuation on the AVO signature: The Leading Edge, 11, 1120– 1125. Claerbout, J., 1985, Imaging the earth’s interior: Blackwell Scientific Publications, Boston. ——–, 1992, Earth soundings analysis: Processing versus inversion: Blackwell Scientific Publications, Boston. Cui, J. and J. He, 2004, 2D seismic migration with compensation: a preliminary study: Journal of Geophysical and Engineering, 1, 263–267. Dai, N. and G. F. West, 1994, Inverse Q migration, in 64th SEG, Expanded Abstracts, 1418–1421, Soc. of Expl. Geophys. Dasgupta, R. and R. A. Clark, 1998, Estimation of Q from surface seismic reflection data: Geophysics, 63, 2120–2128. 91 Bibliography Day, S. M. and J. B. Minster, 1984, Numerical simulation of attenuated wavefields using a Pade approximate method: Geophys. J. R. Astr. Soc., 78, 105–118. Deffenbaugh, M., A. Shatilo, B. Schneider, and M. Zhang, 2000, Resolution of converted waves in attenuation media: 70th SEG, expanded abstract, 1189–1192, Soc. of Expl. Geophys. Du, S., 1996, Seismic dynamics (in Chinese): University of Petroleum Publishing Co., Dongying, China. Duren, R. E. and R. L. Heestand, 1995, General solutions and causality for a Voigt medium: geophysics, 60, 252–261. Dvorkin, J. and G. Mavko, 2006, Modeling attenuation in reservoir and non-reservoir rock: The Leading Edge, 25, 194–197. Dvorkin, J., G. Mavko, and A. Nur, 1995, Squirt flow in fully saturated rocks: Geophysics, 60, 97–107. Ebrom, D., 2004, The low-frequency gas shadow on seismic sections: The leading edge, 23, 772. Futterman, W. I., 1962, Dispersive body waves: J. Geophys. Res., 67, 5279–5291. Gabor, D., 1946, Theory of communication: Journal of IEEE, 93, 429–457. Gao, F., G. Fradelizio, A. Levander, G. Pratt, and C. Zelt, 2005, Seismic velocity, Q, geological structure and lithology estimation at a ground water contamination site: 65th SEG, Expanded Abstracts, 1561–1564, Soc. of Expl. Geophys. Gazdag, J., 1978, Wave equation migration with the phase-shift method: Geophysics, 43, 1342–1351. Guerra, R. and S. Leaney, 2006, Q(z) model building using walkaway VSP data: Geophysics, 71, V127–v132. Guitton, A., 2004, Amplitude and kinematic corrections of migrated images for non-unitary imaging operators: Geophysics, 69, 1017–1024. Hale, D., 1991, Stable explicit depth extrapolation of seismic wavefields: Geophysics, 56, 1770–1777. Hargreaves, N. D. and A. J. Calvert, 1991, Inverse Q-filtering by Fourier transform: Geophysics, 56, 519–527. Hermann, F. J., 2005, Robust curvelet-domain data continuation with sparseness constraints: 67th EAGE, Expanded Abstracts, 4248–4251, Eur. Assn. of Geoscientists and Engineers. 92 Bibliography Hicks, G. J. and R. G. Pratt, 2001, Reflection waveform inversion using local descent methods: Estimating attenuation and velocity over a gas-sand deposit: Geophysics, 66, 598–612. Hu, J., G. T. Schuster, and P. Valasek, 2001, Poststack migration deconvolution: Geophysics, 66, 939–952. H¨ ubert, L., U. Strecker, J. Dvorkin, and K. A. Festervoll, 2005, Seismic attenuation and hybrid attributes to reduce exploration risk-North Sea case study: 75th SEG, Expanded Abstracts, 436–439. Hugonnet, P., P. Herrmann, and C. Tibeiro, 2001, High resolution Radon: a review: 63th EAGE, Expanded Abstracts, 3044–3047, Eur. Assn. of Geoscientists and Engineers. Jansson, P., 1997, Deconvolution of images and spetra, 2nd ed.: Academic Press. Jaynes, E. T., 2004, Probability theory: The logic of science: Cambridge University Press. http://bayes.wustl.edu/etj/prob.html. Johnson, D. L., 2001, Theory of frequency dependent acoustics in patchysaturated porous media: J. Acoust. Soc. Am., 110, 682–694. Kennett, B., 1983, Seismic wave propagation in stratified media: Cambridge Univ. Press. Kjartansson, E., 1979, Constant Q wave propagation and attenuation: J. Geophys. Res., 84, 4737–4748. Klimentos, T., 1995, Attenuation of P- and S-waves as a method of distinguishing gas and condensate from oil and water: Geophysics, 60, 555–556. Kneib, G. and S. Shapiro, 1995, Viscoacoustic wave propagation in 2-D random media and separation of absorption and scattering attenuation: Geophysics, 60, 459–467. Lancaster, S. J. and M. C. Tanis, 2004, High density 3D prestack Q estimation: 66th EAGE, Expanded Abstracts, P0004, Eur. Assn. of Geoscientists and Engineers. Latimer, R. B., R. Davidson, and P. van Riel, 2000, An interpreter’s guide to understanding and working with seismic-derived acoustic impedance data: The Leading Edge, 19, 242–256. Liu, H., D. L. Anderson, and H. Kanamori, 1976, Velocity dispersion due to anelasticity; implications for seismology and mantle composition: Geophys. J. R. astr. Soc., 47, 42–58. 93 Bibliography Mallat, S., 1999, A wavelet tour of signal processing: Academic Press. Malvern, L., 1969, Introduction to the mechanics of a continuous media: Prentice-Hall Ltd. Toronto, New Jersey. Margrave, G. F., D. C. Henley, M. P. Lamoureux, V. Iliescu, and J. P. Grossman, 2003, Gabor deconvolution revisited: 73rd SEG, Expanded Abstracts, 714–717, Soc. of Expl. Geophys. Masoomzadeh, H., P. Barton, and S. C. Singh, 2004, Non-stretch stacking in the tau-p domain: Exploiting long-offset arrivals for sub-basalt imaging: 74th SEG, Expanded Abstracts, 2132–2155, Soc. of Expl. Geophys. Mavko, G., J. Dvokin, and J. Wales, 2005a, A rock physics and attenuation analysis of a well from the Gulf of Mexico: 75th SEG, expanded abstract, 1585–1588, Soc. of Expl. Geophys. ——–, 2005b, A theoretical estimate of S-wave attenuation in sediment: 75th SEG, expanded abstract, 1469–1472, Soc. of Expl. Geophys. Mavko, G. and A. Nur, 1979, Wave attenuation in partially saturated rocks: Geophysics, 44, 161–178. McGinn, A. and B. Duijndam, 1998, Land seismic data quality improvements: The Leading Edge, 17, 1570–1577. Mittet, R., R. Sollie, and K. Hokstad, 1995, Prestack depth migration with compensation for absorption and dispersion: Geophysics, 60, 1485–1494. Nemeth, T., C. Wu, and G. T. Schuster, 1999, Least-squares migration of incomplete reflection data: Geophysics, 64, 208–221. O’Doherty, R. F. and N. A. Anstey, 1971, Reflections on amplitudes: Geophysical Prospecting, 19, 430–458. Oldenburg, D. W., T. Scheuer, and S. Levy, 1983, Recovery of the acoustic impedance from reflection seismograms: Geophysics, 48, 1318–1337. Ostrander, W. J., 1984, Plane-wave reflection coefficients for gas sands at nonnormal angles-of-incidence: Geophysics, 49, 1637–1648. Parra, J. O. and C. Hacket, 2002, Wave attenuation attributes as flow unit indicators: The Leading Edge, 21, 564–572. Parra, J. O., C. Hacket, P.-C. Xu, and H. Collier, 2006, Attenuation analysis of acoustic waveforms in a borehole intercepted by a sand-shale sequence reservoir: The Leading Edge, 21, 186–193. Perroud, H. and M. Tygel, 2004, Nostretch NMO: Geophysics, 69, 599–607. 94 Bibliography Plessix, R.-E., 2006, Estimation of velocity and attenuation coefficient maps from crosswell seismic data: Geophysics, 71, S235–S240. Pratt, R. G., K. Bauer, and M. Weber, 2003, Crosshole waveform tomography velocity and attenuation images of arctic gas hydrates: 73rd SEG, Expanded Abstracts, 2255–2258, Soc. of Expl. Geophys. Pride, S. R., 2003, Permeability dependence of seismic amplitudes: The Leading Edge, 22, 518–525. Pride, S. R., J. G. Berryman, and J. M. Harris, 2004, Seismic attenuation due to wave induced flow: Journal of Geophysical Research, 109, B01201. Quan, Y. and J. M. Harris, 1997, Seismic attenuation tomography using the frequency shift method: Geophysics, 62, 895–905. Richards, P. G. and W. Menke, 1983, The apparent attenuation of a scattering medium: Bulletin of the Seismological Society of America, 73, 1005– 1021. Ricker, N., 1953, The form and laws of propagation of seismic wavelets: Geophysics, 18, 10–40. ——–, 1977, Transient waves in visco-elastic media: Elsevier Scientific Pub. Co. Rickett, J., 2003, Illumination-based normalization for wave-equation depth migration: Geophysics, 68, 1371–1379. ——–, 2006, Integrated estimation of interval-attenuation profiles: Geophysics, 71, A19–A23. Robinson, E. A. and S. Treitel, 2002, Geophysical signal analysis: Society of Exploration Geophysicists Publications. Robinson, J., 1979, A technique for the continuous representation of dispersion in seismic data: Geophysics, 44, 1345–1351. Sacchi, M., T. J. Ulrych, and C. Walker, 1998, Interpolation and extrapolation using a high resolution discrete Fourier transform: IEEE Transactions on Signal Processing, 46, 31–38. Sacchi, M. and J. Wang, 2007, Estimation of the diagonal of the migration blurring kernel through a stochastic approximation: 77th SEG, Expanded Abstracts, 2438–2441, Soc. of Expl. Geophys. Sacchi, M., J. Wang, and H. Kuehl, 2006, Regularized migration/inversion: New generation of seismic imaging algorithms: Recorder, CSEG, 55–59. 95 Bibliography Sams, M. S., J. P. Neep, M. H. Worthington, and M. S. King, 1997, The measurement of velocity dispersion and frequency-dependent intrinsic attenuation in sedimentary rocks: Geophysics, 62, 1456–1464. Sato, H. and M. C. Fehler, 1998, Seismic wave propagation and scattering in the heterogeneous earth: Springer-Verlag New York , Inc. Schuster, G. T. and J. Hu, 2000, Migration green’s function: Continuous recording geometry: Geophysics, 65, 167–175. Singleton, S., M. T. Taner, and S. Treitel, 2006, Q estimation using GaborMorlet joint time-frequency analysis techniques: 76th SEG, expanded abstracts, 1610–1614. Spencer, T. W., J. R. Sonnad, and T. M. Butler, 1982, Seismic Qstratigraphy or dissipation: Geophysics, 47, 16–24. Taner, M. T., A. F. Koehler, and R. E. Sheriff, 1979, Complex seismic trace analysis: Geophysics, 44, 1041–1063. Taner, M. T. and S. Treitel, 2003, A robust method for Q estimation: 73rd SEG, Expanded Abstracts, 710–713, Soc. of Expl. Geophys. Tarantola, A., 1987, Inverse problem theory. methods for data fitting and model parameter estimation: Elsevier Science Publ. Co. Toksoz, M. N., D. H. Johnston, and A. T. Timur, 1979a, Attenuation seismic waves in dry and saturated rocks-I laboratory measurements: Geophysics, 44, 681–690. ——–, 1979b, Attenuation seismic waves in dry and saturated rocks-II mechanisms: Geophysics, 44, 691–711. Tonn, R., 1991, The determination of the seismic quality factor Q from VSP data: A comparison of different computational method: Geophysical Prospecting, 39, 1–27. Toverud, T. and B. Ursin, 2005, Comparison of seismic attenuation models using zero-offset vertical seismic profiling VSP data: Geophysics, 70, F17– F25. Toxopeus, G., S. Petersen, J. Thorbecke, and K. Wappenaar, 2004, Simulating high resolution inversion: decomposing the resolution function: 74th SEG, Expanded Abstracts, 3044–3047, Soc. of Expl. Geophys. Ulrych, T. J. and M. D. Sacchi, 2005, Information-based inversion and processing with applications: Elsevier Publishing Company. Varela, C. L., A. L. Rossa, and T. Ulrych, 1993, Modeling of attenuation and dispersion: Geophysics, 58, 1167–11733. 96 Bibliography Versteeg, R., 1994, The marmousi experience: Velocity model determination on a synthetic complex dataset: The Leading Edge, 13, 927–936. Wang, Y., 2002, A stable and efficient approach to inverse Q filtering: Geophysics, 67, 657–663. ——–, 2003, Quantifying the effectiveness of stabilized inverse Q filtering: Geophysics, 68, 337–345. ——–, 2006, Inverse Q-filter for seismic resolution enhancement: Geophysics, 71, V51–V60. ——–, 2008, Inverse Q filtered migration: Geophysics, 74, P1–P8. Wang, Y. and J. Guo, 2004, Modified Kolsky model for seismic attenuation and dispersion: Journal of Geophysics and Engineering, 1, 187–196. Watanabe, T., K. T. Nihei, S. Nakagawa, and L. R. Myer, 2004, Viscoacoustic waveform inversion of transmission data for velocity and attenuation: J. Acoust. Soc. Am., 115, 3059–3067. White, J. E., 1975, Computed seismic speeds and attenuation in rocks with partial gas saturation: Geophysics, 40, 224–232. ——–, 1983, Underground sound application of seismic waves: Elsevier Science Publishing Company Inc. Winkler, K., A. Nur, and M. Gladwin, 1979, Friction and seismic attenuation in rocks: Nature, 277:528. Xie, X., S. Jin, and R. Wu, 2006, Wave-equation-based illumination analysis: Geophysics, 71, S169–S177. Yilmaz, O., 2000, Seismic data analysis: Society of Exploration Geophysicist Publications. Yu, J., J. Hu, G. T. Schuster, and R. Estill, 2006, Prestack migration deconvolution: Geophysics, 71, S53–S62. Yu, Y., R. Lu, and M. Deal, 2002, Compensation for the effects of shallow gas attenuation with viscoacoustic wave-equation migration: 72nd SEG, Expanded Abstracts, 2062–2065, Soc. of Expl. Geophys. Zhang, C., 2000, Absorption compensation in seismic data processing: University of British Columbia, Master thesis. Zhang, C. and T. J. Ulrych, 2002, Estimation of quality factors from CMP records: Geophysics, 67, 1542–1547. ——–, 2007, Seismic absorption compensation: a least squares inverse scheme: Geophysics, 72, R109–R114. 97
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Seismic absorption estimation and compensation
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Seismic absorption estimation and compensation Zhang, Changjun 2008
pdf
Page Metadata
Item Metadata
Title | Seismic absorption estimation and compensation |
Creator |
Zhang, Changjun |
Publisher | University of British Columbia |
Date Issued | 2008 |
Description | As seismic waves travel through the earth, the visco-elasticity of the earth's medium will cause energy dissipation and waveform distortion. This phenomenon is referred to as seismic absorption or attenuation. The absorptive property of a medium can be described by a quality factor Q, which determines the energy decay and a velocity dispersion relationship. Four new ideas have been developed in this thesis to deal with the estimation and application of seismic absorption. By assuming that the amplitude spectrum of a seismic wavelet may be modeled by that of a Ricker wavelet, an analytical relation has been derived to estimate a quality factor from the seismic data peak frequency variation with time. This relation plays a central role in quality factor estimation problems. To estimate interval Q for reservoir description, a method called reflectivity guided seismic attenuation analysis is proposed. This method first estimates peak frequencies at a common midpoint location, then correlates the peak frequency with sparsely-distributed reflectivities, and finally calculates Q values from the peak frequencies at the reflectivity locations. The peak frequency is estimated from the prestack CMP gather using peak frequency variation with offset analysis which is similar to amplitude variation with offset analysis in implementation. The estimated Q section has the same layer boundaries of the acoustic impedance or other layer properties. Therefore, the seismic attenuation property obtained with the guide of reflectivity is easy to interpret for the purpose of reservoir description. To overcome the instability problem of conventional inverse Q filtering, Q compensation is formulated as a least-squares (LS) inverse problem based on statistical theory. The matrix of forward modeling is composed of time-variant wavelets. The LS de-absorption is solved by an iterative non-parametric approach. To compensate for absorption in migrated seismic sections, a refocusing technique is developed using non-stationary multi-dimensional deconvolution. A numerical method is introduced to calculate the blurring function in layered media, and a least squares inverse scheme is used to remove the blurring effect in order to refocus the migrated image. This refocusing process can be used as an alternative to regular migration with absorption compensation. |
Extent | 1419131 bytes |
Subject |
Seismic absorption compensation Quality factor Q estimation Reservoir description Migrated image |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-12-01 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0052450 |
URI | http://hdl.handle.net/2429/2820 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2009-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2009_spring_zhang_changjun.pdf [ 1.35MB ]
- Metadata
- JSON: 24-1.0052450.json
- JSON-LD: 24-1.0052450-ld.json
- RDF/XML (Pretty): 24-1.0052450-rdf.xml
- RDF/JSON: 24-1.0052450-rdf.json
- Turtle: 24-1.0052450-turtle.txt
- N-Triples: 24-1.0052450-rdf-ntriples.txt
- Original Record: 24-1.0052450-source.json
- Full Text
- 24-1.0052450-fulltext.txt
- Citation
- 24-1.0052450.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.24.1-0052450/manifest