@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix dc: . @prefix skos: . vivo:departmentOrSchool "Science, Faculty of"@en, "Earth, Ocean and Atmospheric Sciences, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Zhang, Changjun"@en ; dcterms:issued "2008-12-01T22:00:51Z"@en, "2008"@en ; vivo:relatedDegree "Doctor of Philosophy - PhD"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms: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."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/2820?expand=metadata"@en ; dcterms:extent "1419131 bytes"@en ; dc:format "application/pdf"@en ; skos:note "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 phe- nomenon 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 process- ing, the information about seismic absorption can be used to enhance seismic data resolution by means of absorption compensation. The absorptive prop- erty 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 es- timation 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 analyti- cal 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 resolu- tion 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 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 using an ana- lytical 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 implementa- tion. 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 filter- ii 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 deconvolu- tion. 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 pro- cess can be used as an alternative to regular migration with absorption com- pensation. Theoretical derivations and numerical tests showing the validity of the new methods are detailed in the chapters. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Seismic Absorption Phenomena . . . . . . . . . . . . . . . . 1 1.3 Seismic Absorption Research . . . . . . . . . . . . . . . . . . 1 1.4 Basic Assumptions . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4.1 Effective Absorption . . . . . . . . . . . . . . . . . . . 3 1.4.2 Frequency Independent Q . . . . . . . . . . . . . . . 4 1.5 Thesis Structure . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Mechanism and Description of Seismic Absorption . . . . 6 2.1 Physical Hypothesis . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1 Scattering Attenuation and Intrinsic Attenuation . . 7 2.1.2 Causes of Intrinsic Attenuation . . . . . . . . . . . . 7 2.2 Mathematical Description . . . . . . . . . . . . . . . . . . . . 8 2.2.1 The Voigt Model . . . . . . . . . . . . . . . . . . . . 9 2.2.2 The Generalized Linear Solid Model . . . . . . . . . . 11 2.2.3 Describing Absorption by Quality Factors . . . . . . . 12 2.2.4 Choice and Discussion . . . . . . . . . . . . . . . . . 14 2.2.5 Absorption and Seismic Wave Extrapolation . . . . . 16 2.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 iv Table of Contents 3 Quality Factor Estimation from Seismic Peak Frequency Variation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2 Absorption and Peak Frequency Variation . . . . . . . . . . 22 3.2.1 Spectral Ratio Method . . . . . . . . . . . . . . . . . 23 3.2.2 Q and Peak Frequency Variation . . . . . . . . . . . . 25 3.3 Q Estimation for Seismic Absorption Compensation and Mi- gration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3.1 Estimation of Q from CMP Gathers . . . . . . . . . . 28 3.3.2 Q Estimation from a Poststack Trace . . . . . . . . . 33 3.4 Summary and Discussions . . . . . . . . . . . . . . . . . . . 35 4 Seismic Absorption Analysis for Reservoir Description . . 38 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.2 Review of the Methods for Seismic Attenuation Analysis . . 39 4.2.1 Iso-frequency Analysis . . . . . . . . . . . . . . . . . 39 4.2.2 Direct Q Estimation . . . . . . . . . . . . . . . . . . . 40 4.2.3 Attenuation Tomography . . . . . . . . . . . . . . . . 41 4.3 Reflectivity-guided Q Analysis, RGQA . . . . . . . . . . . . 42 4.3.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.3.2 Peak Frequency Variation with Offset (PFVO) Anal- ysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.3.3 Implementation Steps . . . . . . . . . . . . . . . . . . 46 4.4 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . 47 4.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . 51 5 Seismic Absorption Compensation: a Least Squares Inverse Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.2 Issues of Absorption Compensation . . . . . . . . . . . . . . 54 5.2.1 Method Choice . . . . . . . . . . . . . . . . . . . . . 54 5.2.2 Inaccurate Q Estimation . . . . . . . . . . . . . . . . 54 5.2.3 Random Noise . . . . . . . . . . . . . . . . . . . . . . 55 5.3 Time-variant Deconvolution as an Inverse Problem . . . . . . 58 5.3.1 Noise and Deconvolution . . . . . . . . . . . . . . . . 58 5.3.2 Bayes Probabilistic Inference . . . . . . . . . . . . . . 59 5.3.3 Absorption Compensation and Inversion . . . . . . . 61 5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 v Table of Contents 6 Refocusing Migrated Images in Absorptive Media . . . . . 69 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.2 Blurring Function caused by Seismic Absorption . . . . . . . 71 6.2.1 Analytical Analysis . . . . . . . . . . . . . . . . . . . 72 6.2.2 Numerical Computation . . . . . . . . . . . . . . . . 72 6.3 Removing Absorption Blurring Function using a Least Squares Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 6.4 Data Experiments . . . . . . . . . . . . . . . . . . . . . . . . 81 6.5 Conclusion and Discussion . . . . . . . . . . . . . . . . . . . 85 7 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . 86 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 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 Attenuation is determined by composite rock physical prop- erties, such as lithology, porosity, saturation, pore fluid, etc.. 9 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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Amplitude attenuation comparison with Q = 30. The red curve is for the linear Q model and the green one is for Kjar- tansson’s model. . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4 Seismic wave behavior in absorptive media defined by v, ρ and Q. (a) Without absorption. (b) With absorption. . . . . 17 2.5 Evolution of a unit impulse with time in a 1-D absorptive medium (Q = 30). . . . . . . . . . . . . . . . . . . . . . . . . 18 2.6 A Green’s function in an absorptive 2-D medium. . . . . . . . 19 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 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 3.3 Variation of peak frequencies of the two reflections in Fig- ure 3.2. Peak frequencies of the first event and the second event are shown in black and red respectively. . . . . . . . . . 33 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 viii List of Figures 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. . . . . . . . . . . . . . . . . . . . . 40 4.2 Frequency decay caused by thin-bed tuning and absorption. (a) shows a reflectivity function with closely-spaced coeffi- cients. (b) is the corresponding absorption property. (c) is the synthetic trace and (d) is its spectrum of time-frequency decomposition. . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.3 Flow diagram of reflectivity-guided Q analysis. . . . . . . . . 43 4.4 Peak frequencies at layer boundaries. (a) is the peak frequen- cies from the trace in Figure 4.1 (b). (b) is the extracted inverse Q curve. . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.5 Peak frequency variation with offset. (a) is a CMP gather built from the model shown in Figure 4.2. (b) is the corre- sponding peak frequencies at zero offset. (c) is the relative gradient of peak frequencies along offset, and (d) is the esti- mated Q curve. . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.6 A real data example. (a) is a CMP gather. (b) is the cor- responding non-stationary distributed peak frequencies. (c) is the peak frequency curve at zero offset. (d) is the rela- tive gradient curve of peak frequencies. (e) is the reflectivity function and (f) is the estimated Q curve. . . . . . . . . . . . 49 4.7 Real data experiment. (a) is the peak frequency section and (b) is the attenuation section (Q−1). . . . . . . . . . . . . . . 50 5.1 Experiments on amplitude compensation using different Q values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.2 A synthetic trace with three spikes, Q = 30, and its phase correction using different Q values. . . . . . . . . . . . . . . . 57 5.3 The amplitude of an ideal inverse Q filter (black) and its approximation (red). The two curves are overlapping before the fork. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 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 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 5.6 Q compensation for a real seismic section. (a) shows a real seismic section. (b) is the result of absorption compensation. 67 ix List of Figures 6.1 Diagram of the energy distribution of a migration blurring function. Symbol “V” represents a receiver location. . . . . . 73 6.2 Examples of the blurring functions in an absorptive medium. (a) at time t = 400 ms. (b) at time t = 700 ms. . . . . . . . . 76 6.3 Amplitude of a kernel matrix G. . . . . . . . . . . . . . . . . 80 6.4 A synthetic example. (a) shows an absorption-blurred sine- function shaped structure. (b) is the refocusing result. . . . . 82 6.5 Input information. (a) is the velocity section and (b) is the inverse Q section. . . . . . . . . . . . . . . . . . . . . . . . . . 83 6.6 Migration results. (a) Time-variant spectral whitening plus phase-shift migration. (b) Phase-shift migration plus refocus- ing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 co- supervise 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 phe- nomena observed in seismic exploration. The topics addressed are the esti- mation 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 elas- tic 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 prop- agation 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 seis- mic 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 at- tenuation 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 informa- tion 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 Fut- terman model and others have been used for different applications for a long time (Ricker, 1953; Biot, 1956a,b; Futterman, 1962; White, 1975; Tok- soz et al., 1979a,b; Aki and Richards, 1980; Johnson, 2001). More compli- cated 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 al- lowed 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 imple- mented 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 re- garding 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 de- 2 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; Harg- reaves and Calvert, 1991; Bickel, 1993; Varela et al., 1993; Yilmaz, 2000). Computation techniques trying to stabilize inverse Q filtering have been in- troduced 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 attenu- ation to reservoir description, due to the fact that for porous rocks satu- rated with fluid, generally, strong absorption can be observed. Parra and Hacket (2002), Castagna et al. (2003), Taner and Treitel (2003), Hübert 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 de- pending on the fluid content in rock fractures. Winkler et al. (1979), Klimen- tos (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 three- dimensional 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 ef- fective 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 dependentQ 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 absorp- tion. 4 Chapter 1. Introduction Chapter 2 reviews physical mechanisms and mathematical descriptions of seismic absorption. A definition of quality factor Q and the correspond- ing 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 decon- volution 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. Convention- ally, 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 propa- gating 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 environ- mental 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, perme- ability, 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 atten- uation, 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 atten- uation on seismic data are the same, in that they both attenuate seismic amplitude and distort the seismic waveform. Scattering attenuation is gen- erally 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). Intrin- sic 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 ab- sorption 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 thermo- elasticity as the most viable model to explain intrinsic attenuation at litho- spheric 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 attenu- ation predicted by Biot’s model is extremely small for frequencies less than 100Hz. Pride et al. (2004) introduced a model which explains seismic at- tenuation in porous rocks as resulting from wave-induced fluid flow. Their model also confirms that “squirt-flow” in Biot’s model is incapable of ex- plaining 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 mecha- nisms 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 de- scribed 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 de- rived 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 general- ized 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 inves- tigated 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 equa- tion 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 = ρ ∂2ui ∂t2 . The constitutive relation leads to the equation of motion in term of particle displacement (λ+ µ)(∇ · θ) + µ∇2~u = ρ∂ 2~u ∂t2 , (2.2) where ~u is a vector of displacement, ∂ 2~u ∂t2 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µ) ∂2uz ∂z2 = ρ ∂2uz ∂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µ) ∂2uz ∂z2 + (λ′ + 2µ′) ∂3uz ∂t∂z2 = ρ ∂2uz ∂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 ′) ∂2Uz ∂z2 = −ρω2Uz, (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, ω) = U0e ±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 < ω20 , attenuation increases as the square of the frequency, and velocity is approximately constant (Ricker, 1977). These observations are expressed as follows ap = 1 2 √ M/ρ ω2 ω0 , and vp = √ M/ρ. The expression for ap means that attenuation increases with the square of frequency. However, many measurements have indicated a first power de- pendence 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 M(z, ω) =Mu(z, ω) + N∑ 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 pa- rameters 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 constitu- tive 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 1 Q = −∆E 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 1 Q = −∆A πA0 , (2.10) 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 ∆A = λ dA(z) dz , (2.11) 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 ab- sorptive 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ωtdt = e−iωz/v(ω), will now be attenuated by a factor e−α(ω)z with α(ω) = ω2v(ω)Q . The prop- agating pulse shape G(z, t), thus, has a Fourier transform eikz, and k is complex and is given by k = ω v(ω) + iα(ω). If the impulse response is causal, G(z, t) = 0, for t < z/v(∞), it can be shown that (Aki and Richards, 1980) ω v(ω) = ω v(∞) +H[α(ω)]. (2.14) Here, v(∞) is the limit of v(ω) as ω →∞ and H[α(ω)] is the Hilbert trans- form of the attenuation factor. This equation is equivalent to the Kramers- Kronig relation in electro-magnetic theory (Futterman, 1962). From a math- ematical 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 inde- pendentQ, a velocity dispersion relation, due to attenuation, can be derived, v(ω0) v(ω) = 1 + 1 πQ ln (ω0 ω ) , (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 v(ω) v(ω0) = [−iω ω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 Q ≈ tan(πγ). 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 atten- uation 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 Chapter 2. Mechanism and Description of Seismic Absorption 0 20 40 60 80 100 120 Frequency (Hz) 1800 2000 Ve lo ci ty (m /s) 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. Ap- plication 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 Chapter 2. Mechanism and Description of Seismic Absorption 0 20 40 60 80 100 120 Frequency (Hz) Am pl itu de 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. With- out 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: kc(ω) = ω v(ω0) [ 1 + 1 πQ ln (ω0 ω )]( 1 + i 1 2Q ) . (2.17) 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 ω v(ω0) ∆z ) , (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, ω), (2.19) with A1(∆z, ω) = exp ( − ω∆z 2Qv(ω0) ) (2.20) defining an amplitude attenuation term, and A2(∆z, ω) = exp [ −i ω∆z πQv(ω0) ln ( ω ω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 ab- sorptive 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 be- comes increasingly wider as traveltime increases and its amplitude decreases continuously due to absorption. 17 Chapter 2. Mechanism and Description of Seismic Absorption 0 400 800 1200 1600 2000 2400 Time (MS) 1 2 3 4 5 Tr ac e In de x 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 0.5 1.0 1.5 Ti m e (se c) -1.0 -0.5 0 0.5 1.0 Distance (km) 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. Fig- ure 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 uni- form, 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 hetero- geneity, 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., atten- uation 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 seis- mic 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 ab- sorption based on the frequency shift observed in vertical seismic profiling (VSP) data. Since the centroid frequency of the signal’s spectrum experi- ences a shift to lower frequencies during propagation, they develop a relation between Q and the centroid of an amplitude spectrum. The amplitude spec- trum is represented by a Gaussian, boxcar, or triangular function. In this chapter, an analytical equation will be developed to relate ab- sorption 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 con- tinuous 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. B(f) = 2√ π f2 f2m e−f 2/f2m , (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 as- sumption 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 ( − ∫ ray a(f, l)dl ) . 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 B(f, t) = B(f)e −pift Q . (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 0.5 1.0 1.5 Ti m e (s) 5 10 15 20 25 Trace index (a) 0 10 20 30 40 50 60 70 80 Frequency (Hz) 0 0.005 Am pl itu de (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 ampli- tudes. (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, B(f, t2) B(f, t1) = B(f)e −pift2 Q B(f)e− pift1 Q . (3.4) Taking logarithms of both sides, equation (3.4) becomes ln [ B(f, t2) B(f, t1) ] = −π(t2 − t1) Q f. (3.5) Denoting Ar = ln [ B(f,t2) 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 Q = −π(t2 − t1) p . (3.6) 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 geo- metrical 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 B(f, t) =M(t)B(f)e −pift Q , (3.7) whereM(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: ∂B(f, t) ∂f =M(t) ∂B(f) ∂f e− pift Q +M(t)B(f)e− pift Q ( −πt Q ) = 0. (3.8) Recalling the expression for the Ricker wavelet, from equation (3.1), we obtain ∂B(f) ∂f = 2√ π ( 2f f2m ) e − f2 f2m + 2 f2 ( f2 f2m ) e − f2 f2m (−2f f2m ) . (3.9) Finally, by inserting this expression into equation (3.8), the peak frequency at time t is obtained as fp = f 2 m   √( πt 4Q )2 + ( 1 fm )2 − πt 4Q   . (3.10) The relation between Q and the shift of peak frequency is, Q = πtfpf 2 m 2(f2m − f2p ) . (3.11) This expression shows that if the dominant frequency fm is known, the Q- factor 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 = πt1fp1f 2 m 2(f2m − f2p1) = πt2fp2f 2 m 2(f2m − f2p2) . Thus, the dominant frequency of the source wavelet can be derived from the peak frequencies of a reflection at two different time locations: fm = √ fp1fp2(t2fp1 − t1fp2) t2fp2 − t1fp1 . (3.12) Combining equation (3.11) and (3.12), one obtains Q = π(t2 − t1)fp2f2p1 2(f2p1 − f2p2) . (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: B(f) = exp [−(f − fm)2 σ2 ] . (3.14) 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 = fp1t2 − fp2t1 t2 − t1 . Similarly, observations at two different time locations allow us to determine a unique Q-factor Q = πσ2(t2 − t1) fp1 − fp2 . (3.15) By assuming the amplitude spectrum of a seismic wavelet has an analyti- cal 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 struc- ture. 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 −pift1 Q1 e −pift2 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 equa- tion (3.7) to (3.11), if Q1 and the dominant frequency fm of the source wavelet are known. The expression of Q2 is Q2 = πt2Q1 αQ1 − πt1 , (3.17) where α = 2f2m − 2f2p fpf2m . (3.18) Now, using the concept of an equivalent Q, in other words, letting e −pift1 Q1 e −pift2 Q2 = e− pift 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 = t2Q1Q (t1 + t2)Q1 − t1Q. (3.20) The equivalent Q can be calculated from equation (3.19), using the peak fre- quency 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 am- plitude attenuation equation (3.3) as B(f, t) = A(t)B(f) exp ( N∑ i=1 −πf∆ti Qi ) , (3.21) where Qi and ∆ti are the quality factor and the traveltime in layer i , re- spectively. Following the approach used in velocity estimation, we assume straight raypaths and compute the total traveltime of a reflection at a par- ticular offset as N∑ i=1 ∆ti = tN . Now define t0(N) as the zero-offset traveltime of reflection N , QN as the Q- factor 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(N) [t0(i)− t0(i− 1)] . 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, Qinti , from the RMS values using a relation similar to the Dix formula (Yilmaz, 2000): Qinti = √ Q2i t0(i) −Q2i−1t0(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 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 compen- sation 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 under- estimated 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 val- ues 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 (equa- tion (2.21)), after the correction, the phase is φc(ω) = ω∆z πv(ω) −∆Q (Qo +∆Q)Qo ln ( ω ω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 minimum- phase, 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 500 1000 1500 2000 Time (ms) 0 0.1 0.2 0.3 (a) Zero-phased impulse response with absorption. 0 500 1000 1500 2000 Time (ms) 0 0.2 0.4 0.6 0.8 (b) Amplitude compensation using actual Q = 30. 0 500 1000 1500 2000 Time (ms) 0 0.1 0.2 0.3 (c) Amplitude compensation using big Q = 40. 0 500 1000 1500 2000 Time (ms) 0 2 4 6 8 10 12 14 (d) Amplitude correction using small Q = 25. Figure 5.1: Experiments on amplitude compensation using different Q val- ues. 56 Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme 0 500 1000 1500 2000 Time (ms) 0 0.1 0.2 0.3 (a) Impulse response with absorption. 0 500 1000 1500 2000 Time (ms) 0 0.1 0.2 0.3 (b) Phase correction using Q = 30. 0 500 1000 1500 2000 Time (ms) 0 0.1 0.2 0.3 (c) Phase correction using Q = 20. 0 500 1000 1500 2000 Time (ms) 0 0.1 0.2 0.3 (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 0 10 20 30 40 50 60 70 80 Frequency (Hz) 10 20 30 40 50 60 Am pl itu de Figure 5.3: The amplitude of an ideal inverse Q filter (black) and its ap- proximation (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 Time-variant Deconvolution as an Inverse Problem 5.3.1 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 of- ten 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 con- sider 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 ob- served 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 dis- cretized 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 dis- tributed 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 av̀is 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 p(x) = 1√ 2πσ e − h x−µ 2σ2 i2 , (5.4) 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 + x2 2σ2 )−1 , (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ès 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, σ2n), 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πσ2n )(M−1)/2 e−(1/2σ 2 n)||d−Gm)||22, (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 form conditional on a param- eter σ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 reflec- tivity m, with elements mk of m i.i.d. with equal variance σ 2 m (a standard 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) = ∏ k ( 1 + m2k 2σ2m ) . By inserting this Cauchy prior and the data likelihood (equation (5.4)) into equation (5.7) and taking logarithms of both sides, we obtain ln[p(m|d, σm, σn)] = −c(m)− 1 2σ2n (d−Gm)T (d−Gm), (5.8) where c(m) is a constraint imposed by the Cauchy distribution c(m) = M−1∑ k=0 ln ( 1 + m2k 2σ2m ) . Furthermore, denoting φcg(m) = − ln[p(m|d, σm, σn)], we obtain the cost function for the Cauchy-Gauss model as φcg(m) = c(m) + 1 2σ2n (d−Gm)T (d−Gm), (5.9) where G is composed of time variant wavelets, bτ (t − τ). Both t and τ are in the range of the length of a trace and G =   b0(t− 0) b1(t− 1) ... 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 Cauchy- Gauss 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 min- imized is φcg(m) in equation (5.9). Taking the derivative of φcg(m) with respect to m, we obtain ∂φcg(m) ∂m = ∂c(m) ∂m + 1 2σ2n GT (d−Gm) (5.10) and ∂c(m) ∂m = 1 σ2m S−1m, where S is a M ×M diagonal matrix with elements skk = 1 + m 2 k σm , (k = 0, 1, . . . ,M − 1) and is m dependent. Equating the derivative in equation (5.10) to zero, yields m = ( λS−1 +GTG ) GTd, (5.11) where λ = 2σ n 2 m2 . This equation is nonlinear and must be solved iteratively. To construct an iterative algorithm, equation (5.11) is written as m = SGT ( λIN +GSG T )−1 d = SGTp. (5.12) p is the auxiliary vector which is obtained from the solution of the system( λIN +GSG T ) 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)GTp(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 dra- matically 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 syn- thetic 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 time- variant 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 com- putations yielding the kernel matrix G. The inversion begins by setting a value for the parameter λ in equation (5.13) and the result is updated it- eratively. Following inversion, the wavelet is convolved with the extracted reflectivities. The final results are shown in Figure 5.6(b). A number of re- flections, 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 con- tinuity 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 radi- cally from the deconvolution techniques in customary use. The main differ- ence 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 pro- vides more accurate information concerning both the amplitude and location of the earth’s reflectivity, hydrocarbon reservoir characterization is one obvi- ous 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 viabil- ity 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 econom- ically. 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 us- ing, typically, time-variant spectral whitening, time-variant deconvolution (Yilmaz, 2000), least squares inversion (Zhang and Ulrych, 2007), one- dimensional 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 com- pensation 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 investi- gated 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 ex- trapolation 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 post- stack 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 al- gorithm 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 sup- pression 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 amplifi- cation 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 Q- compensation 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 sam- pling, 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 in- vestigated 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 qual- ity 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 tech- niques 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 migra- tion 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 mod- eling 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 = LTLm0, (6.1) wherem0 is the true image or the real earth structure, andm 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 mi- gration blurring function can be defined as Γ(r|r0) = LTLqδ(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 = [ LTL ]−1 m. However, the inverse of LTL 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 LTLq 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, ω) = e −j ω c(ω) |rg−r0| |rg − r0| . In this expression, ωc(ω) is a complex wavenumber. Absorption caused wave- field attenuation and dispersion are represented through the complex veloc- ity 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, G(r|rg , ω) = G∗(rg|r, ω) = e j ω v0 |rg−r| |rg − r| . The migration blurring function in an absorptive medium, if the migration is implemented without absorption compensation, is Γ(r|r0, ω) = ∫ rg e −j ω c(ω) |rg−r0| |rg − r0| e j ω v0 |rg−r| |rg − r| drg. (6.3) The blurring function in the time domain is the integral of all frequency components Γ(r|r0, t) = ∫ [∫ rg e −j ω c(ω) |rg−r0| |rg − r0| e j ω v0 |rg−r| |rg − r| drg ] dω. (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 ho- mogeneous 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 func- tion. 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 com- pute synthetic data with attenuation Lqδ(r) and then migrate the synthetic to obtain the migration blurring function LTLqδ(r). The synthetic data can be computed using phase-shift modeling, or us- ing ray-tracing plus wavelet continuation along the ray-paths. Phase-shift modeling is based on the following wavefield extrapolation operator H(ω, z) = e−jz √ k2a−k2x , (6.5) 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 = ω v0 [ 1− 1 πQ ln ( ω ω0 )] [ 1 + j 1 Q ] , (6.6) or ka = ω c(ω) , with 1 c(ω) = 1 v0 [ 1− 1 πQ ln ( ω ω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 expen- sive 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 cal- culated 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) 4piQ , 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 mi- gration is used for the imaging, the modeling algorithm is better to use the phase-shift method. If Kirchhoff migration is used for the imaging, corre- spondingly, 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 U(kx, z = 0, ω|z0) = n∏ i=0 e j∆z r ω2 ci(ω) 2−k2x , and the response of a unit impulse right below it at depth z0 +∆z is U(kx, z = 0, ω|z0 +∆z) = e j∆z r ω2 cn+1(ω) 2−k2x 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 Γ(kx, z, ω|z0) = U(kx, z = 0, ω|z0) m∏ i=0 e −j∆z r ω2 v2 i −k2x . (6.8) 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 Γ(kx, z, ω|z0 +∆z) = U(kx, z = 0, ω|z0 +∆z) m∏ i=0 e −j∆z r ω2 v2 i −k2x . (6.9) Examining equations (6.8) and (6.9), we can observe that, in layered me- dia, 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. Γ(kx, z, ω|z0 +∆z) = e j∆z r ω2 cn+1(ω) 2−k2x Γ(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 propertyQ 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 Γ = LTLq is not an identity matrix. In our case, the migrated image is a blurred image due to using an inaccurate migration operator, m = LTLqm0. (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 LTLq 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 0 0.5 1.0 Tim e ( s) 1.0 1.2 1.4 1.6 1.8 Midpoint (km) (a) 0 0.5 1.0 Tim e ( s) 1.0 1.2 1.4 1.6 1.8 Midpoint (km) (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 multi- dimensional 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 m̂(~k, t) = γ̂~xref ,tj ( ~k, t) ∗ m̂0(~k, t) + n̂(~k, t). (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 2σ2n (d̂−Gm̂0)T (d̂−Gm̂0). (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 Gaussian noise in the observed data. m̂0 and d̂ represent model and data in the wavenumber domain respectively. The objective function can be minimized through solving the following equation, m̂0 = ( GTG+ λW−1 )−1 GT d̂, (6.15) 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 G =   γ̂0(t− 0∆t) γ̂1(t− 1∆t) ... γ̂i(t− i∆t) ... γ̂N (t−N∆t)   . (6.16) 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 i |mi|+ ǫ). 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, an- other is to solve form0 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 func- tion 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 wavenum- ber 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