Seismic Absorption Estimation andCompensationbyChangjun ZhangM.Sc., The University of British Columbia, 2001A THESIS SUBMITTED IN PARTIAL FULFILMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinThe Faculty of Graduate Studies(Geophysics)The University of British Columbia(Vancouver)November, 2008c Changjun Zhang 2008AbstractAs seismic waves travel through the earth, the visco-elasticity of the earth?smedium 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 asan important attribute to interpret for fluid units. In seismic data process-ing, the information about seismic absorption can be used to enhance seismicdata resolution by means of absorption compensation. The absorptive prop-erty of a medium can be described by a quality factor Q, which determinesthe energy decay and a velocity dispersion relationship. The quality factorand 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 infour chapters respectively. By assuming that the amplitude spectrum of aseismic wavelet may be modeled by that of a Ricker wavelet, an analyti-cal relation has been derived to estimate a quality factor from the seismicdata peak frequency variation with time. This relation plays a central rolein quality factor estimation problems. Quality factors with coarse resolu-tion can be estimated from prestack common midpoint (CMP) gathers, ora 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 estimatespeak frequencies at a common midpoint location, then correlates the peakfrequency with sparsely-distributed reflectivities, and finally calculates Qvalues from the peak frequencies at the reflectivity locations using an ana-lytical expression. The peak frequency is estimated from the prestack CMPgather using peak frequency variation with offset (PFVO) analysis whichis similar to amplitude variation with offset (AVO) analysis in implementa-tion. The estimated Q section has the same layer boundaries of the acousticimpedance or other layer properties. Therefore, the seismic attenuationproperty obtained with the guide of reflectivity is easy to interpret for thepurpose of reservoir description.To overcome the instability problem of conventional inverse Q filter-iiAbstracting, Q compensation is formulated as a least-squares (LS) inverse problembased on statistical theory. The matrix of forward modeling is composedof time-variant wavelets. The LS de-absorption is solved by an iterativenon-parametric approach.To compensate for absorption in migrated seismic sections, a refocusingtechnique is developed using non-stationary multi-dimensional deconvolu-tion. A numerical method is introduced to calculate the blurring functionin layered media, and a least squares inverse scheme is used to remove theblurring 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 validityof the new methods are detailed in the chapters.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . xiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Seismic Absorption Phenomena . . . . . . . . . . . . . . . . 11.3 Seismic Absorption Research . . . . . . . . . . . . . . . . . . 11.4 Basic Assumptions . . . . . . . . . . . . . . . . . . . . . . . . 31.4.1 Effective Absorption . . . . . . . . . . . . . . . . . . . 31.4.2 Frequency Independent Q . . . . . . . . . . . . . . . 41.5 Thesis Structure . . . . . . . . . . . . . . . . . . . . . . . . . 42 Mechanism and Description of Seismic Absorption . . . . 62.1 Physical Hypothesis . . . . . . . . . . . . . . . . . . . . . . 62.1.1 Scattering Attenuation and Intrinsic Attenuation . . 72.1.2 Causes of Intrinsic Attenuation . . . . . . . . . . . . 72.2 Mathematical Description . . . . . . . . . . . . . . . . . . . . 82.2.1 The Voigt Model . . . . . . . . . . . . . . . . . . . . 92.2.2 The Generalized Linear Solid Model . . . . . . . . . . 112.2.3 Describing Absorption by Quality Factors . . . . . . . 122.2.4 Choice and Discussion . . . . . . . . . . . . . . . . . 142.2.5 Absorption and Seismic Wave Extrapolation . . . . . 162.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19ivTable of Contents3 Quality Factor Estimation from Seismic Peak FrequencyVariation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.2 Absorption and Peak Frequency Variation . . . . . . . . . . 223.2.1 Spectral Ratio Method . . . . . . . . . . . . . . . . . 233.2.2 Q and Peak Frequency Variation . . . . . . . . . . . . 253.3 Q Estimation for Seismic Absorption Compensation and Mi-gration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.3.1 Estimation of Q from CMP Gathers . . . . . . . . . . 283.3.2 Q Estimation from a Poststack Trace . . . . . . . . . 333.4 Summary and Discussions . . . . . . . . . . . . . . . . . . . 354 Seismic Absorption Analysis for Reservoir Description . . 384.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.2 Review of the Methods for Seismic Attenuation Analysis . . 394.2.1 Iso-frequency Analysis . . . . . . . . . . . . . . . . . 394.2.2 Direct Q Estimation . . . . . . . . . . . . . . . . . . . 404.2.3 Attenuation Tomography . . . . . . . . . . . . . . . . 414.3 Reflectivity-guided Q Analysis, RGQA . . . . . . . . . . . . 424.3.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 424.3.2 Peak Frequency Variation with Offset (PFVO) Anal-ysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.3.3 Implementation Steps . . . . . . . . . . . . . . . . . . 464.4 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . 474.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . 515 Seismic Absorption Compensation: a Least Squares InverseScheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 525.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 525.2 Issues of Absorption Compensation . . . . . . . . . . . . . . 545.2.1 Method Choice . . . . . . . . . . . . . . . . . . . . . 545.2.2 Inaccurate Q Estimation . . . . . . . . . . . . . . . . 545.2.3 Random Noise . . . . . . . . . . . . . . . . . . . . . . 555.3 Time-variant Deconvolution as an Inverse Problem . . . . . . 585.3.1 Noise and Deconvolution . . . . . . . . . . . . . . . . 585.3.2 Bayes Probabilistic Inference . . . . . . . . . . . . . . 595.3.3 Absorption Compensation and Inversion . . . . . . . 615.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68vTable of Contents6 Refocusing Migrated Images in Absorptive Media . . . . . 696.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 696.2 Blurring Function caused by Seismic Absorption . . . . . . . 716.2.1 Analytical Analysis . . . . . . . . . . . . . . . . . . . 726.2.2 Numerical Computation . . . . . . . . . . . . . . . . 726.3 Removing Absorption Blurring Function using a Least SquaresMethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 756.4 Data Experiments . . . . . . . . . . . . . . . . . . . . . . . . 816.5 Conclusion and Discussion . . . . . . . . . . . . . . . . . . . 857 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . 86Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90viList of Tables3.1 A pseudo code of prestack Q estimation. . . . . . . . . . . . . 316.1 A pseudo code of the refocusing algorithm. . . . . . . . . . . 79viiList of Figures2.1 Attenuation is determined by composite rock physical prop-erties, such as lithology, porosity, saturation, pore fluid, etc.. 92.2 Velocity dispersion comparison with Q = 30. The red curve isfor the linear Q model and the green one is for Kjartansson?smodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.3 Amplitude attenuation comparison with Q = 30. The redcurve is for the linear Q model and the green one is for Kjar-tansson?s model. . . . . . . . . . . . . . . . . . . . . . . . . . 162.4 Seismic wave behavior in absorptive media defined by v, rhoand Q. (a) Without absorption. (b) With absorption. . . . . 172.5 Evolution of a unit impulse with time in a 1-D absorptivemedium (Q = 30). . . . . . . . . . . . . . . . . . . . . . . . . 182.6 A Green?s function in an absorptive 2-D medium. . . . . . . . 193.1 (a) A reflection in a CMP gather. The medium is absorptivewith Q = 10. Traces have been normalized based on theirmaximum amplitudes. (b) Amplitude spectra of the originalsource signature (black), trace 1 (red), trace 11 (green), andtrace 21 (gray). . . . . . . . . . . . . . . . . . . . . . . . . . . 243.2 (a) A synthetic CMP gather of two reflections with absorptionand 10% random noise. Traces have been normalized basedon their maximum amplitude. (b) Values of input qualityfactors (Q1 = 10 and Q2 = 20 ). (c) Computed qualityfactors (Q1 = 10.04 and Q2 = 20.12). . . . . . . . . . . . . . . 323.3 Variation of peak frequencies of the two reflections in Fig-ure 3.2. Peak frequencies of the first event and the secondevent are shown in black and red respectively. . . . . . . . . . 333.4 Q estimation. (a) is a real seismic trace. (b) is its windowedtime-variant spectrum. The picked peak frequency points aremarked by circles and connected by straight lines. (c) is theinverse of estimated quality factors. . . . . . . . . . . . . . . . 36viiiList of Figures4.1 Frequency decay caused by absorption. (a) is an attenuationmodel, (b) is the synthetic trace and (c) is the spectrum offrequency decomposition. . . . . . . . . . . . . . . . . . . . . 404.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) isthe synthetic trace and (d) is its spectrum of time-frequencydecomposition. . . . . . . . . . . . . . . . . . . . . . . . . . . 414.3 Flow diagram of reflectivity-guided Q analysis. . . . . . . . . 434.4 Peak frequencies at layer boundaries. (a) is the peak frequen-cies from the trace in Figure 4.1 (b). (b) is the extractedinverse Q curve. . . . . . . . . . . . . . . . . . . . . . . . . . . 444.5 Peak frequency variation with offset. (a) is a CMP gatherbuilt from the model shown in Figure 4.2. (b) is the corre-sponding peak frequencies at zero offset. (c) is the relativegradient of peak frequencies along offset, and (d) is the esti-mated Q curve. . . . . . . . . . . . . . . . . . . . . . . . . . . 484.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 reflectivityfunction and (f) is the estimated Q curve. . . . . . . . . . . . 494.7 Real data experiment. (a) is the peak frequency section and(b) is the attenuation section (Q-1). . . . . . . . . . . . . . . 505.1 Experiments on amplitude compensation using different Qvalues. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 565.2 A synthetic trace with three spikes, Q = 30, and its phasecorrection using different Q values. . . . . . . . . . . . . . . . 575.3 The amplitude of an ideal inverse Q filter (black) and itsapproximation (red). The two curves are overlapping beforethe fork. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 585.4 Convergence process of the inverse scheme. (e) is a synthetictrace. (a-d) are the inverse results after iteration 1, 2, 3 and4 respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 655.5 (a) is a synthetic trace with 20% noise. (b) is the result ofabsorption compensation. The rightmost trace is the originalreflectivity series plotted as a reference. . . . . . . . . . . . . 665.6 Q compensation for a real seismic section. (a) shows a realseismic section. (b) is the result of absorption compensation. 67ixList of Figures6.1 Diagram of the energy distribution of a migration blurringfunction. Symbol ?V? represents a receiver location. . . . . . 736.2 Examples of the blurring functions in an absorptive medium.(a) at time t = 400 ms. (b) at time t = 700 ms. . . . . . . . . 766.3 Amplitude of a kernel matrix G. . . . . . . . . . . . . . . . . 806.4 A synthetic example. (a) shows an absorption-blurred sine-function shaped structure. (b) is the refocusing result. . . . . 826.5 Input information. (a) is the velocity section and (b) is theinverse Q section. . . . . . . . . . . . . . . . . . . . . . . . . . 836.6 Migration results. (a) Time-variant spectral whitening plusphase-shift migration. (b) Phase-shift migration plus refocus-ing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84xAcknowledgmentsI would like to express my deepest thanks to my supervisor Dr. TadeuszUlrych for his moral support and scientific enthusiasm on my research. Hisinvaluable 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 tookfrom him is very useful in my academic pursuit. My thanks also go to Dr.Andrew Calvert, Dr. Elizabeth Hearn and Dr. Matthew Yedlin for theirparticipation as committee members. Discussions with and questions frommy supervisors and other committee members helped to shape my researchon seismic absorption.My special thanks go to Dr. Bill Nickerson. He kindly proofread mymanuscripts and helped me correct errors in both English and rock physics.Part of the research included in this thesis was supported by a UBCinternational scholarship and grants from the companies that supported theConsortium for the Development of Specialized Seismic Techniques, CDSST,which was led by Dr. Tadeusz Ulrych.xiDedicationTo the memory of my fatherYaoqiang Zhang.xiiChapter 1Introduction1.1 Thesis ObjectivesThe 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 attenuatedsignals.1.2 Seismic Absorption PhenomenaConventionally, in seismic exploration, the earth is modeled as an ideal elas-tic medium, and seismic wave propagation is explained by means of theelastic (or acoustic) wave equation. In practice, the propagation of seismicwaves 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 notaccurate enough to describe the wave behavior for this more complicatedmedium. Generally, the visco-elasticity of the earth materials causes seis-mic energy dissipation, and thus decreases the amplitude and modifies thefrequency content of propagating waves. This phenomenon of wave energydissipation is called seismic absorption or seismic attenuation, although at-tenuation is just one aspect of the wave behavior in visco-elastic media. Fromthe studies of absorptive phenomena observed in seismic data, we may beable to construct images with better resolution in seismic data processingand extract more detailed information about the rock materials in seismicdata inversion.1.3 Seismic Absorption ResearchAbsorption researchers have made contributions in several related areas:from explaining the absorption mechanism, estimating absorptive properties1Chapter 1. Introductionand modeling its effect numerically and physically, to removing its effect, andutilizing absorptive properties to locate and describe hydrocarbon reservoirs.In recent years, due to the requirement of obtaining seismic images withhigher resolution and higher accuracy, as well as of extracting more informa-tion for reservoir description in seismic exploration, seismic absorption hasbecome a field of extensive investigation, similar to other fields that take intoaccount secondary effects in seismic wave propagation, such as anisotropy,angle-dependent reflectivity and wave-mode conversion.Mathematical models, resulting from modification to the elastic waveequation, like the Voigt model, the standard linear solid model, the Fut-terman model and others have been used for different applications for along 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 byDvorkin et al. (1995), Pride et al. (2004), and Carcione and Picotti (2006).These models relate seismic attenuation to the porous properties of rockmaterials, such as porosity, saturation, fluid content, permeability, etc., andtry to fit field data and laboratory observations more accurately than al-lowed by means of traditional methods. An overview of seismic attenuationmodels can be found in Toverud and Ursin (2005).Numerical simulation of seismic waves in absorptive media is generallyimplemented 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 havealso 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 importanttools to confirm the accuracy of theories and discriminate mechanisms re-garding seismic absorption.The seismic absorption property of a medium is usually described bya quality factor designated by Q. Many algorithms to estimate Q valueshave been published. Among them, the most popular is the spectral ratiomethod (Spencer et al., 1982; Tonn, 1991). Methods suggested by Quan andHarris (1997), Zhang and Ulrych (2002), Taner and Treitel (2003), Singletonet al. (2006), Rickett (2006) and Guerra and Leaney (2006) are claimed tobe more robust than the traditional spectral ratio methods, and/or disclosemore detailed attenuation information of rock properties. The objectives ofQ 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-2Chapter 1. Introductionscription, we require Q information of higher resolution for the study offine layers. These two different requirements are similar to the resolutionrequirements for velocities in seismic data processing and inversion.To remove the effect of absorption, traditionally we can use time-variantspectral 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 alsobeen suggested by Berkhout (1981), Dai and West (1994), Yu et al. (2002),Cui and He (2004) and Wang (2008) for stacked seismic sections and Mittetet al. (1995) for prestack seismic data. Migration plus Q compensation canbe 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 andHacket (2002), Castagna et al. (2003), Taner and Treitel (2003), H?bertet al. (2005), Chapman et al. (2005) and Dvorkin and Mavko (2006) usedifferent physical hypotheses to relate attenuation properties to reservoircharacterization. There is more and more interest in seismic exploration inusing attenuation properties as an attribute in reservoir description, becauseof reported successful cases.In a medium, shear waves and compressional waves travel at differentspeeds. 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 Mavkoet al. (2005b) have investigated the characteristics of shear wave attenuationin seismic data.1.4 Basic Assumptions1.4.1 Effective AbsorptionSeismic attenuation is usually categorized into scattering attenuation andintrinsic attenuation. Scattering attenuation is generally caused by three-dimensional heterogeneities in the subsurface that distribute wave energyin arbitrary directions. Often quoted research about scattering attenuationincludes O?Doherty and Anstey (1971), Richards and Menke (1983) and Satoand Fehler (1998).3Chapter 1. IntroductionIntrinsic attenuation is caused by internal friction among grains in therock matrix, and the relative movement between the solid rock matrix andthe pore fluid, when a seismic wave travels through loosely consolidated orporous media. Internal friction and fluid flow are considered directly relatedto rock porosity, permeability and the fluid content in the void space. Moredetailed explanation of the physics of intrinsic attenuation can be foundin 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 intrinsicattenuation. This is based on the premise that scattering attenuation andintrinsic attenuation have similar features in observed seismic or well logdata, and also that, at present, there is no reliable way to separate them.1.4.2 Frequency Independent QWhen using Q to describe seismic attenuation, it is often assumed thatQ does not change with frequency in the seismic frequency range (10 --200Hz). However, there are many mathematical models which describe thechange of Q with frequency (Ricker, 1953; Pride et al., 2004; Carcione andPicotti, 2006). Frequency dependent Q models may be theoretically moreattractive in visco-elastic wave equation derivation, and/or better fittinglaboratory testing data. The usability of these more complicated seismicattenuation models is hindered by their dependence on some unmeasurableparameters and the necessity to use frequency dependent Q values in practice(Dvorkin and Mavko, 2006).When Q is considered as a function of frequency, it would be convenientif the Q function can be possibly separated into a frequency dependentpart and a frequency independent part, so that we can use the frequencyindependent part to describe attention property. In this thesis, we use thesimplest convention: absorption property Q does not depend on frequency,while attenuation changes with frequency. The assumption is valid over thesignal bandwidth in seismic exploration and conforms to many theoreticalmodels and laboratory observations.1.5 Thesis StructureThis chapter introduces the background of the research on seismic absorp-tion.4Chapter 1. IntroductionChapter 2 reviews physical mechanisms and mathematical descriptionsof seismic absorption. A definition of quality factor Q and the correspond-ing dispersion relation between attenuation and velocity are analyzed. Thelinear dispersion model is adopted in this thesis because it is appropriatenot only for handling the problem of estimating quality factors from bothprestack and poststack seismic data, but also for absorption compensation.Chapter 3 develops a method to estimate quality factors from waveletpeak frequency variation. The analytical formula derived can be used toestimate Q from a CMP gather, and also from a stacked seismic trace.The estimated Q provides the information required for one-dimensional andmulti-dimensional Q compensation.Chapter 4 suggests a scheme called reflectivity-guided Q analysis forreservoir 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 locationsof 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 sois to overcome the instability problem of regular inverse Q filtering. Thekernel matrix of the inversion is composed of time-variant wavelets. The LSde-absorption can be solved iteratively with fast convergence. The resultsof de-absorption are related to the accuracy of the estimated Q values andalso of the seismic wavelets.Chapter 6 investigates the blurring effect in migrated images when usinga regular migration algorithm to migrate seismic data with absorption. Thedeblurring problem is considered as a multi-dimensional time-variant decon-volution and is solved using a LS inverse scheme in time-wavenumber domainto refocus migrated images. The refocusing processing is an alternative totraditional post-stack migration with absorption compensation.Chapter 7 draws conclusions from the research, proposes future researchtopics, and discusses some other issues in seismic absorption which are notdetailed in the previous chapters.5Chapter 2Mechanism and Descriptionof Seismic AbsorptionSeismic waves traveling in the earth experience amplitude attenuation andvelocity 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, rho). In practice, the propagation of seismic waves inthe earth is in many respects different from seismic wave propagation in anideal solid. When a seismic wave experiences absorption during its propa-gating process, the wave behavior can be described by velocity, density anda quality factor (v, rho, Q), if the absorptive property of a medium is fullydetermined by a quality factor (Q).There are different physical hypotheses to explain signal attenuationobserved in seismic data and there are also numerous mathematical formulasavailable to describe the attenuation phenomenon. This chapter investigatesthe causes of the absorption phenomenon and how absorption is describedin terms of wave theory by reviewing the appropriate works in the publishedliterature.2.1 Physical HypothesisA great deal of attention has been given to the study of the attenuationof seismic waves (Ricker, 1977; White, 1983; Mavko and Nur, 1979; Toksozet al., 1979b; Kneib and Shapiro, 1995; Carcione, 1995; Dvorkin et al., 1995;Johnson, 2001; Pride et al., 2004). Loss parameters for rocks have beenmeasured by many techniques over a wide range of frequency and environ-mental conditions. Physical mechanisms for energy loss have been proposedand described mathematically in varying degrees of detail.Absorption is a composite effect caused by the earth media which have awide range of properties, such as lithology, inhomogeneity, porosity, perme-ability, saturation, fluid contents, etc.. There is no general consensus as to6Chapter 2. Mechanism and Description of Seismic Absorptionwhat the dominant loss mechanism is. In different environments, a differentmechanism may play a dominant role.2.1.1 Scattering Attenuation and Intrinsic AttenuationThere are two mechanisms that are usually considered to cause seismic atten-uation, scattering and intrinsic attenuation. Scattering redistributes waveenergy within the medium but does not remove energy from the overallwavefield. Conversely, intrinsic attenuation refers to various mechanismsthat convert vibrational energy into heat through friction, fluid flow, andthermal relaxation processes. The effect of scattering and intrinsic atten-uation on seismic data are the same, in that they both attenuate seismicamplitude and distort the seismic waveform. Scattering attenuation is gen-erally caused by relatively simple factors, such as three-dimensional randomheterogeneity in rocks, which distributes wave energy in arbitrary directions,and interfaces which reflect and redirect wave propagation. (O?Doherty andAnstey, 1971; Richards and Menke, 1983; Sato and Fehler, 1998). Intrin-sic attenuation is caused by more complicated rock properties. Observedseismic attenuation is generally the integrated effects of both scattering andintrinsic 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 AttenuationAs has been stated, seismic absorption is a comprehensive effect caused bymany factors. For intrinsic attenuation, two factors are generally consideredthe most important. One is the internal friction resulting from the relativeslide 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 energyfrom a seismic source is transformed into relative movement and dissipatesinside 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 cracksalong with the amount of attenuation caused by thermo elasticity are inclosest agreement with observations.Rocks have microscopic cracks and pores which may contain fluids. Thesecracks can have a profound influence on the propagation velocities of Pand S-waves. Biot (1956a) and Biot (1956b) analyzed wave propagation inisotropic porous solids where the coupling of motion between the fluid and7Chapter 2. Mechanism and Description of Seismic Absorptionthe solid matrix was considered. He arrived at expressions for attenuationdue 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 than100Hz. Pride et al. (2004) introduced a model which explains seismic at-tenuation in porous rocks as resulting from wave-induced fluid flow. Theirmodel 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 arethe main contributors of observed seismic attenuation. These factors arelithological variation and fluid-saturated pore spaces in a medium. Researchcarried out by Carcione and Picotti (2006) shows that the most effective lossmechanisms 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 widerange of properties possessed by earth materials. Conclusions drawn fromdata of oil-reservoir rocks may not be applicable to materials in the uppermantle. 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 expectedto describe losses in all rocks under all environmental conditions.Physical mechanism studies help us to understand how rock absorptionproperties are related to rock physical properties by means of suitable rockphysics modelling. Although several such models have been introduced inthe cited references, there is no consensus on which model is best. Figure 2.1illustrates a generic idea that seismic attenuation is the composite effect ofrock physical parameters, and that the attenuation property can be denotedby Q. A rock body is very complicated. To describe it needs numerousparameters, 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 thewave behavior defined by rho, v, and Q to rock physical properties. Seismicinversion tries to deduce rock properties from seismic attributes. Q is animportant attribute for rock property inversion.2.2 Mathematical DescriptionIn 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 of8Chapter 2. Mechanism and Description of Seismic AbsorptionFigure 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 directlyfrom observed amplitude decay and corresponding attenuation-dispersionrelations, describe seismic absorption, generally, by means of a Q parameter(Aki and Richards, 1980; Kjartansson, 1979).2.2.1 The Voigt ModelThe Voigt model, also known as the Kelvin-Voigt model, is a widely inves-tigated method of introducing losses that has the advantage of yielding alinear wave equation which can be solved for arbitrary time dependence. Theassumption is made that stresses are directly proportional to strain rates,as well as to the components of strain themselves. This assumption wasintroduced independently by Stokes, Kelvin and Voigt (Ricker, 1977), andits 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 theEinstein notation, Hooke?s law is expressed assigmaij = lambdathetadeltaij + 2?epsilon1ij, (2.1)9Chapter 2. Mechanism and Description of Seismic Absorptionwhere sigmaij is the stress tensor, epsilon1ij is the strain tensor, and deltaij is the Kroneckerdelta function. lambda and ? are the two Lame?s parameters. theta is the volumestrain or dilatation. If there is no body force, using Newton?s Lawsigmaij,j = rhopartialdiff2uipartialdifft2 .The constitutive relation leads to the equation of motion in term of particledisplacement(lambda + ?)(nabla? theta) + ?nabla2vector = rhopartialdiff2vectorupartialdifft2 , (2.2)where vector is a vector of displacement, partialdiff2vectorupartialdifft2 is its second derivative with respectto time and rho is density. For a plane wave, if we assume the displacementis parallel to the depth z-axis (ux and uy are both zero) and that uz isindependent of x and y, equation (2.2) reduces to(lambda + 2?)partialdiff2uzpartialdiffz2 = rhopartialdiff2uzpartialdifft2 , (2.3)for one-dimensional wave propagation.In absorptive media, the constitutive relation must incorporate the factthat stress is not only proportional to strain, it is also proportional to thetime derivative of strain. The stress-strain behavior in absorptive media isdescribed by a modified Hooke?s law which includes strain-rate terms:sigmaij =parenleftbigglambda + lambdaprime partialdiffpartialdifftparenrightbiggthetadeltaij + 2parenleftbigg? + ?prime partialdiffpartialdifftparenrightbiggepsilon1ij,where lambdaprime and ?prime are perturbations of Lame?s parameters, lambda and ?, due toabsorption. Equation (2.3) in an absorptive medium becomes(lambda + 2?) partialdiff2uzpartialdiffz2 + (lambdaprime + 2?prime) partialdiff3uzpartialdifftpartialdiffz2 = rhopartialdiff2uzpartialdifft2 . (2.4)This equation is a Stokes equation without body force. If we replace Lame?sparameters by Young?s modulus, M = lambda + 2? and Mprime = lambdaprime + 2?prime, andtransform equation(2.4) into the frequency domain, we have(M + iomegaMprime)partialdiff2Uzpartialdiffz2 = -rhoomega2Uz, (2.5)where Uz = U(z,omega) designates the space-dependent wavefield at an angularfrequency omega. A wavefield satisfying equation (2.5) has the following formU(z,omega) = U0e?G(omega)z,10Chapter 2. Mechanism and Description of Seismic AbsorptionandG(omega) =bracketleftbigg- rhoomega2M + iomegaMprimebracketrightbigg1/2. (2.6)This complex propagation parameter G(omega) may be expressed in terms ofattenuation ap and phase velocity vp asG(omega) = ap + iomega/vp. (2.7)By comparing equation (2.6) and equation (2.7) and defining omega0 = M/Mprime,it can be shown that when omega2 < omega20, attenuation increases as the square ofthe frequency, and velocity is approximately constant (Ricker, 1977). Theseobservations are expressed as followsap = 12radicalbigM/rho omega2omega0 ,andvp =radicalbigM/rho.The expression for ap means that attenuation increases with the square offrequency. 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, thereare different observations from that of White (1983). Duren and Heestand(1995) and Buckingham (2005) have derived independently causal solutionsto the Stokes equation (equation (2.4)).2.2.2 The Generalized Linear Solid ModelThe constitutive relation between stress sigma and strain epsilon1 for visco-elasticity canbe formulated simply in the frequency domain as (Malvern, 1969; Carcioneet al., 1988)sigma(omega) = M(omega)epsilon1(omega).Here M(omega) is the complex, frequency dependent, visco-elastic modulus. Ageneralized linear solid body can be described byM(z,omega) = Mu(z,omega) +Nsummationdisplayj=1ajdeltaMj(z)iomega/omegaj(z) -1. (2.8)For one-dimensional P wave, in this equation, Mu(z,omega) is the un-relaxedYoung?s modulus at a depth z, N is the number of absorptive mechanisms.11Chapter 2. Mechanism and Description of Seismic Absorptionaj is a constant coefficient, omegaj is the jth relaxation frequency and deltaMj(z) isthe difference between the un-relaxed Young?s modulus and the contributionto the Young?s modulus due to the jth relaxation mechanism. These pa-rameters can be adjusted to fit any realistic linear absorption law, accordingto Carcione et al. (1988).The two models described in this section and the previous section areconstructed by assuming that anelasticity of rocks will result in a constitu-tive relation that is different from that of elastic models. Absorption mayalso be described by directly assuming velocity dispersion relations.2.2.3 Describing Absorption by Quality FactorsWave speeds in rocks are determined by rock properties. When includingabsorption into the wave equation, it is very straightforward to use complexwave-numbers or complex velocities. The absorptive property determines areal velocity dispersion relation or a complex dispersion relation dependingon different applications (Mittet et al., 1995; Claerbout, 1985).Definition of QThe absorptive property is often represented by the quality factor Q, whichis an intrinsic property of rocks. Formally, Q is defined as a dimensionlessmeasure of the anelasticity which is given by1Q =-deltaE2piE , (2.9)where E is the energy stored at the maximum strain in a volume, and -deltaEis the energy loss in each cycle because of anelasticity. Equation (2.9) impliesthat Q-1 is the portion of energy lost during each cycle or wavelength. Fora medium with a linear stress-strain relation, the amplitude A of a waveat a particular frequency is proportional to radicalE (Aki and Richards, 1980),therefore 1Q =-deltaApiA0 , (2.10)where A0 is the amplitude at the start of a cycle and deltaA represents theamplitude decay in a cycle. It is observed that amplitude varies with distancebecause of absorption. We can rewrite A as a function of distance A(z) ,anddeltaA = lambdadA(z)dz , (2.11)12Chapter 2. Mechanism and Description of Seismic Absorptionwhere lambda is the wavelength given in terms of frequency omega and phase velocityv(omega) bylambda = 2piv(omega)/omega. (2.12)From equations (2.10), (2.11) and (2.12), it can be shown thatA(z) = A0 expbracketleftbigg- omegaz2v(omega)Qbracketrightbigg, (2.13)which describes amplitude attenuation of a frequency component in an ab-sorptive medium.Attenuation Causing Velocity DispersionAbsorption results in both attenuation and dispersion (Liu et al., 1976).As mentioned before, the basic assumptions used to describe absorptionbehavior are constant Q, linearity and causality. Considering a unit impulseas input at depth z = 0, each frequency component of this impulse,integraldisplay infinity-infinitydelta[t -z/v(omega)]e-iomegatdt = e-iomegaz/v(omega),will now be attenuated by a factor e-alpha(omega)z with alpha(omega) = omega2v(omega)Q. The prop-agating pulse shape G(z,t), thus, has a Fourier transform eikz, and k iscomplex and is given byk = omegav(omega) + ialpha(omega).If the impulse response is causal,G(z,t) = 0, for t < z/v(infinity),it can be shown that (Aki and Richards, 1980)omegav(omega) =omegav(infinity) + H[alpha(omega)]. (2.14)Here, v(infinity) is the limit of v(omega) as omega arrowrightinfinity and H[alpha(omega)] 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 thisrelation can be satisfied with constant Q. We must tolerate a frequencydependent Q which is effectively constant over the seismic frequency range.13Chapter 2. Mechanism and Description of Seismic AbsorptionUnder the constraint of causality, and the assumption of frequency inde-pendent Q, a velocity dispersion relation, due to attenuation, can be derived,v(omega0)v(omega) = 1 +1piQ lnparenleftBigomega0omegaparenrightBig, (2.15)From equations (2.13) and (2.15), we can observe that wave propagationin an absorptive medium is completely specified by two parameters, Q anda phase velocity v(omega0) at an arbitrary reference frequency omega0. omega0 oftentakes the value of the Nyquist frequency in digital signal processing forconvenience.Kjartansson?s Absorption ModelAnother method to formulate the velocity dispersion relation was proposedby Kjartansson (1979). The absorption model is defined by means of acomplex velocity dispersion relationv(omega)v(omega0) =bracketleftbigg-iomegaomega0bracketrightbigggamma, (2.16)where gamma is a parameter related to absorption. Equation (2.16) also representsthe causal, constant Q model, and it can be derived1Q approxequal tan(pigamma).The estimation of Q is equivalent to the estimation of gamma. For gamma = 0, equation(2.16) gives a real constant velocity which means no absorption.2.2.4 Choice and DiscussionBesides the foregoing mentioned mathematical models, there are also othermodels that describe seismic attenuation with differing degrees of precisionand utility. A list of these methods can be found in the paper by Toverudand Ursin (2005).In this thesis, the dispersion relation of equation (2.15) is adopted. Aspointed out by Aki and Richards (1980), equation (2.15) is an importantresult since it appears to be a good approximation for a variety of atten-uation laws in which Q is effectively constant over the seismic frequencyrange. Besides a reference velocity v(omega0), this visco-elastic model assumesattenuation and dispersion are determined by only one more parameter Q.Theoretically, this linear absorption model (equations (2.13) and (2.15), also14Chapter 2. Mechanism and Description of Seismic Absorption0 20 40 60 80 100 120Frequency (Hz)18002000Velocity (m/s)Figure 2.2: Velocity dispersion comparison with Q = 30. The red curve isfor the linear Q model and the green one is for Kjartansson?s model.called the Kolsky-Futterman model (Wang and Guo, 2004), is superior tothe Voigt model which has been contradicted by practical observations. Ap-plication is easier than the generalized linear solid model which requires theadjustment of two parameters when being applied to different media.The linear Q model and Kjartansson?s model define different complexwavenumbers. Their real parts are very close. However, the imaginary partsare different. Figure 2.2 and Figure 2.3 compare the velocity dispersionsand amplitude attentions of theses two models. Velocity dispersion andamplitude attenuation are both defined by Q. The red curves result fromthe linear Q model and the green ones result from Kjartansson?s model. Forthe same Q value, Kjartansson?s model has stronger attenuation than thelinear Q model. In seismic data processing, ideally, Q estimation and Qcompensation 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 onthe premise that scattering attenuation and intrinsic attenuation have thesame effect on observed seismic data. When using Q to describe seismicattenuation, it is assumed that Q does not change with frequency. Thisassumption is valid for the seismic data acquired for oil and gas explorationwhich is in the frequency range of, approximately from 10 to 200 Hz (Samset al., 1997).Figure 2.4 illustrates the wave behavior in elastic media and in absorp-15Chapter 2. Mechanism and Description of Seismic Absorption0 20 40 60 80 100 120Frequency (Hz)AmplitudeFigure 2.3: Amplitude attenuation comparison with Q = 30. The red curveis for the linear Q model and the green one is for Kjartansson?s model.tive media. An elastic medium is represented by rho and v. A visco-elasticmedium is represented by rho, v and Q. In order to extract Q information, weneed 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 tothe attenuation of high frequencies.2.2.5 Absorption and Seismic Wave ExtrapolationAbsorption causes energy loss and frequency dispersion. From equations(2.13) and (2.15), it can be concluded that for a seismic wave traveling ina visco-elastic medium, the wavenumber including an attenuation term iscomplex:kc(omega) = omegav(omega0)bracketleftbigg1 + 1piQ lnparenleftBigomega0omegaparenrightBigbracketrightbigg parenleftbigg1 + i 12Qparenrightbigg. (2.17)The real part in this equation is responsible for the wavefront propagationand the phase distortion, the imaginary part is responsible for amplitudeattenuation.For one-dimensional wave propagation, if the wavefield U(z,omega) at a depthz propagates downward to a depth z + deltaz, the wavefield U(z + deltaz,omega) atthis depth in absorptive media isU(z + deltaz,omega) = U(z,omega)A(deltaz,omega)expparenleftbigg-i omegav(omega0)deltazparenrightbigg, (2.18)16Chapter 2. Mechanism and Description of Seismic AbsorptionFigure 2.4: Seismic wave behavior in absorptive media defined by v, rho andQ. (a) Without absorption. (b) With absorption.In equation (2.18), A(deltaz,omega) is an absorption related termA(deltaz,omega) = A1(deltaz,omega)A2(deltaz,omega), (2.19)withA1(deltaz,omega) = expparenleftbigg- omegadeltaz2Qv(omega0)parenrightbigg(2.20)defining an amplitude attenuation term, andA2(deltaz,omega) = expbracketleftbigg-i omegadeltazpiQv(omega0)lnparenleftbigg omegaomega0parenrightbiggbracketrightbigg(2.21)expressing a velocity dispersion term which causes a phase delay. For timedomain wavefield modeling, equation (2.18) can be written asU(t + deltat,omega) = U(t,omega)A(deltat,omega)exp(iomegadeltat), (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 ina 1-D absorptive medium at five time locations which are 400, 800, 1200,1600 and 2000 milliseconds respectively. The impulse initiates at time zeroand the medium has a quality factor of 30. The shape of the impulse be-comes increasingly wider as traveltime increases and its amplitude decreasescontinuously due to absorption.17Chapter 2. Mechanism and Description of Seismic Absorption0 400 800 1200 1600 2000 2400Time (MS)12345Trace IndexFigure 2.5: Evolution of a unit impulse with time in a 1-D absorptivemedium (Q = 30).18Chapter 2. Mechanism and Description of Seismic Absorption00.51.01.5Time (sec)-1.0 -0.5 0 0.5 1.0Distance (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 stepcalculation. The first step does ray-tracing to calculate ray-paths in themedia 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 witha 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 waveletis 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 thatat a near offset.2.3 SummarySeismic attenuation is a complicated process, which can be caused by hetero-geneity, interfaces, lithology, saturation, fluid content and other propertiesof the earth media. Mathematically, the absorption equation can be derivedusing dynamic mechanics or using an attenuation-dispersion relation. In19Chapter 2. Mechanism and Description of Seismic Absorptionthis thesis, the linear Q model (equation (2.13) and (2.15)) is adopted as areasonable assumption over the exploration seismic frequency bandwidth.20Chapter 3Quality Factor Estimationfrom Seismic Peak FrequencyVariationThe fundamental idea to estimate Q from seismic data is, at first to examinethe variation of a seismic wavelet with time, and then calculate a Q valuefrom the variation of the wavelet shape or spectrum. By assuming that theamplitude spectrum of a seismic wavelet is like that of a Ricker wavelet,this chapter derives an analytical formula which allows a Q value to becalculated from the spectral peak frequency variation. Theoretically, thisanalytical approach can be used to estimate Q-factors from seismic dataof different acquisition geometries, where peak frequency variations can bereliably examined. Applications of this approach to estimate Q-factors fromcommon midpoint (CMP) gathers and also from stacked seismic traces willbe discussed.3.1 IntroductionSeismic waves traveling through the earth experience absorption, i.e., atten-uation and dispersion, because of the anelasticity and heterogeneity of themedium (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 seismicimages. This in turn will allow us to better understand the effects of AVOand, 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 welldeveloped (Bachrach et al., 2006; Lancaster and Tanis, 2004; Dasgupta andClark, 1998). However, some research has been published concerning theestimation of Q-factor from vertical seismic profile (VSP) and cross-welldata (Spencer et al., 1982; Tonn, 1991). Almost all of these methods use21Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variationthe amplitude information of received signals, which information, however,is often inaccurate because of noise, geometrical spreading, scattering, andother 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 relationbetween 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 canbe 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 asingle trace, where the spectral peak frequency variation of a seismic waveletcan be picked from a windowed time-variant spectrum (WTVS) or a con-tinuous wavelet transform (CWT) spectrum.3.2 Absorption and Peak Frequency VariationIn seismic data processing, a recorded trace is commonly modeled as theconvolution of a seismic source signature with a reflectivity series. Theseismic source signature is generally unknown. The effects of anelasticitycan be incorporated into the model by convolving with an earth filter. Thisfilter is causal, minimum phase, and depends on the Q-factor (Aki andRichards, 1980).In order to build a relation between Q and peak frequency translation,we begin by assuming that the amplitude spectrum of the source waveletcan be well represented by that of a Ricker wavelet.B(f) = 2pi f2f2m e-f2/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 isthe dominant frequency.One point needs to be emphasized here: the derivation only requires anassumption about the shape of the amplitude spectrum, without any as-sumption of the phase of a seismic wavelet. The Ricker amplitude spectrumis close to that observed in real data. Approximately, a spectrum which22Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variationrises to the peak amplitude quickly and decays gradually can be fitted by aRicker spectrum.The evolution of the amplitude spectrum with time is now modeled asa Ricker wavelet traveling in a visco-elastic medium (geometrical spreadingand other factors are not considered). After traveling for a time t, theamplitude spectrum isB(f,t) = B(f)H(f,t). (3.2)In equation (3.2), H(f,t) is the absorption filter (Varela et al., 1993), whosefrequency response isH(f) = expparenleftbigg-integraldisplayraya(f,l)dlparenrightbigg.In this expression, the integral is evaluated along the ray-path l, anda(f,l) = pifQ(l)v(l).Here, Q(l) and v(l) are the Q-factor and velocity, respectively, defined alongeach point of the ray-path. The value Q(l) is assumed to be independent offrequency (Ricker, 1953; Kjartansson, 1979; White, 1983).By considering the propagation of a wavefield in a half-space with aQ-factor for t seconds, the amplitude spectrum of the received signal isB(f,t) = B(f)e-piftQ . (3.3)This equation can also be derived from equation (2.22), by only consideringthe amplitude of the wavefield.As time increases, absorption increases with frequency and results inthe peak frequency translating towards lower frequency. This phenomenonis illustrated in Figure 3.1. Because of absorption, the time width of thesource wavelet increases and the amplitude spectrum narrows as the wavetravels. If the arrival times are known, the Q-factor can be calculated fromthe spectral variation based on equation (3.3).3.2.1 Spectral Ratio MethodA common method to estimate Q is the spectral ratio method (Spencer et al.,1982; Tonn, 1991). By comparing the amplitude spectra at two arrival times23Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation0.51.01.5Time (s)5 10 15 20 25Trace index(a)0 10 20 30 40 50 60 70 80Frequency (Hz)00.005Amplitude(b)Figure 3.1: (a) A reflection in a CMP gather. The medium is absorptivewith Q = 10. Traces have been normalized based on their maximum ampli-tudes. (b) Amplitude spectra of the original source signature (black), trace1 (red), trace 11 (green), and trace 21 (gray).24Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variationt1 and t2 respectively, we can take the ratio of the two amplitude spectra,B(f,t2)B(f,t1) =B(f)e-pift2QB(f)e-pift1Q. (3.4)Taking logarithms of both sides, equation (3.4) becomeslnbracketleftbiggB(f,t2)B(f,t1)bracketrightbigg= -pi(t2 -t1)Q f. (3.5)Denoting Ar = lnbracketleftBigB(f,t2)B(f,t1)bracketrightBig, the logarithm of the spectral ratio, and plottingAr as a function of frequency, f, yields a linear trend whose slope p is afunction of Q. Q then can be easily expressed in terms of observed slope pasQ = -pi(t2 -t1)p . (3.6)The spectral ratio method is simple in principle, but in practice determiningthe Q-factor is complicated by overlapping wavelets which lead to amplitudespectra that do not reflect the wavelet spectrum. Also, the linear trendof 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 theamplitude spectrum besides attenuation. These factors will be detailed inthe following discussions.3.2.2 Q and Peak Frequency VariationAs stated above, the magnitude of the amplitude spectrum of seismic wavesis affected by many factors besides absorption. These factors include geo-metrical spreading, reflection and transmission effects, and the methods ofamplitude recovery, which may be applied to balance the trace amplitudein seismic data processing.It is almost exclusively the Q-factor, however, that affects the shapeof the wavelet spectrum. Based on this idea, a method is developed toestimate the Q-factor from the spectral peak frequency variation of seismicreflections. Including all Q-factor unrelated functions into an amplitudeterm, the amplitude spectrum can be written asB(f,t) = M(t)B(f)e-piftQ , (3.7)where M(t) is an amplitude factor independent of frequency and absorption.The peak frequency fp can be determined by equating the derivative of the25Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variationspectrum, with respect to frequency, to zero:partialdiffB(f,t)partialdifff = M(t)partialdiffB(f)partialdifff e-piftQ + M(t)B(f)e-piftQparenleftbigg-pitQparenrightbigg= 0. (3.8)Recalling the expression for the Ricker wavelet, from equation (3.1), weobtainpartialdiffB(f)partialdifff =2radicalpiparenleftbigg2ff2mparenrightbigge-f2f2m + 2f2parenleftbiggf2f2mparenrightbigge-f2f2mparenleftbigg-2ff2mparenrightbigg. (3.9)Finally, by inserting this expression into equation (3.8), the peak frequencyat time t is obtained asfp = f2mbracketlefttpbracketleftbtradicalBiggparenleftbiggpit4Qparenrightbigg2+parenleftbigg 1fmparenrightbigg2- pit4Qbracketrighttpbracketrightbt. (3.10)The relation between Q and the shift of peak frequency is,Q = pitfpf2m2(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 beestimated if the amplitude spectrum of the initial source wavelet can beapproximated by that of a Ricker wavelet. Designating the peak frequenciesat times t1 and t2 by fp1 and fp2, respectively, we haveQ = pit1fp1f2m2(f2m -f2p1) =pit2fp2f2m2(f2m -f2p2).Thus, the dominant frequency of the source wavelet can be derived from thepeak frequencies of a reflection at two different time locations:fm =radicalBiggfp1fp2(t2fp1 -t1fp2)t2fp2 -t1fp1 . (3.12)Combining equation (3.11) and (3.12), one obtainsQ = pi(t2 -t1)fp2f2p12(f2p1 -f2p2) . (3.13)26Chapter 3. Quality Factor Estimation from Seismic Peak Frequency VariationEquation (3.13) allows us to determine the absorption property of a layer ifan 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 similarityto the observed amplitude spectra. One example is the Gaussian spectrumof variance sigma2:B(f) = expbracketleftbigg-(f -fm)2sigma2bracketrightbigg. (3.14)Following the above procedure, the relationship between the Q-factor andthe variation of peak frequencies isQ = pitsigma2fm -fp ,and the dominant frequency isfm = fp1t2 -fp2t1t2 -t1.Similarly, observations at two different time locations allow us to determinea unique Q-factorQ = pisigma2(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 bederived.For real seismic data, the amplitude spectrum of a seismic wavelet isnot 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 prestackand poststack seismic data.3.3 Q Estimation for Seismic AbsorptionCompensation and MigrationFor absorption compensation and migration in seismic data processing, Qvalues are only expected for major layers defined by velocities. Therefore,methods to estimate Q values here are based on seismic reflection eventpicking.27Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation3.3.1 Estimation of Q from CMP GathersA 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 forthe extraction of information concerning structure, lithology, and materialproperties such as velocity and Q-factor.Secondly, reflection arrival times are determined by interval velocitiesand the geometrical structure of the subsurface. Absorption of the receivedsignals is only determined by the interval Q-factors and traveltimes in eachlayer. If the amplitude spectrum of a seismic wavelet is assumed to beRicker-like, interval Q-factors can be computed solely from the variation ofthe peak frequency of a spectrum as a function of time.Equation (3.13) allows us to obtain an average Q-factor by using thepeak frequency variation along all offsets, thereby allowing us to remove theeffects of surface fluctuations and random noise, consequently improving theaccuracy of the estimated Q values.Two-layer CaseConsider first the case of two layers with quality factors Q1 and Q2 andtraveltimes t1 and t2 in each layer, respectively, using equation (3.3), theamplitude spectrum after t = t1 + t2 isB(f,t) = A(t)B(f)e-pift1Q1 e-pift2Q2 . (3.16)From the peak frequency fp associated with B(f,t) in equation (3.16), Q2can 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 sourcewavelet are known. The expression of Q2 isQ2 = pit2Q1alphaQ1 -pit1, (3.17)wherealpha = 2f2m -2f2pfpf2m . (3.18)Now, using the concept of an equivalent Q, in other words, lettinge-pift1Q1 e-pift2Q2 = e-piftQ , (3.19)28Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variationthis allows an expression to be derived which relates the interval Q2 to theequivalent Q-factor byQ2 = 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 andQ1 have already been determined from upper-layer arrivals, if traveltimest1 and t2 in the two layers are known, Q2 can be computed from equation(3.20).Multi-layer CaseGiven the complexity of the subsurface, some approximations must be madeto use all offset information. For multi-layer media, we may write the am-plitude attenuation equation (3.3) asB(f,t) = A(t)B(f)expparenleftBigg Nsummationdisplayi=1-pifdeltatiQiparenrightBigg, (3.21)where Qi and deltati are the quality factor and the traveltime in layer i , re-spectively. Following the approach used in velocity estimation, we assumestraight raypaths and compute the total traveltime of a reflection at a par-ticular offset asNsummationdisplayi=1deltati = 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 reflectionabove layer N. Using the proportional property of similar triangles, thetraveltime in layer i isdeltati = tNt0(N)[t0(i) -t0(i -1)] .Furthermore, the amplitude attenuation operator in equation (3.21) can besplit into two factors as follows:B(f,t) = A(t)B(f)expparenleftBiggN-1summationdisplayi=1-pifdeltatiQiparenrightBiggexpparenleftbigg-pifdeltatNQnparenrightbigg. (3.22)29Chapter 3. Quality Factor Estimation from Seismic Peak Frequency VariationThis allows us to obtain the following equation for QN :QN = pideltatNalpha -beta . (3.23)where alpha is the same as that in equation (3.18) andbeta =N-1summationdisplayi=1pideltatiQ .The Q-factor values can now be calculated layer by layer by means of layerstripping.RMS Q and Interval QSince a straight raypath approximation is used, the computed QN is not theactual interval Q-factor. Analogous to root mean square (RMS) velocitiesin seismic velocity analysis, such Q values can be referred to as RMS Qvalues. In a manner similar to interval velocity analysis, we can determinethe apparent interval Q-factor of the ith layer, Qinti , from the RMS valuesusing a relation similar to the Dix formula (Yilmaz, 2000):Qinti =radicalBiggQ2i 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 isbetter to notice that Qi is the average of the calculated Q values at differentoffsets for layer i, and t0(i) is the zero-offset arrival time of reflection i whenthe raypath is assumed to be straight across interfaces. Moreover, sincevalues of Qint are derived from RMS Q values, they are apparent, ratherthan real, interval Q-factors.Numerical ExperimentsTo estimate Q from CMP gathers, it is assumed that the arrival times for themain reflection events are known. Fourier transforms are computed in thewindow containing the reflection at each offset, each amplitude spectrumis fitted with a Ricker spectrum and the peak frequency of the spectrumis estimated. Using equation (3.23), we can now calculate Q factors layerby layer from the peak frequency variations. The procedure of using thesuggested approach to estimate Q-factors from CMP gathers can be outlinedby the the pseudo code in table (3.3.1).30Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation/* Analyze absorption at selected CMP locations */{select_a_CMP ();do while (event index < number of events){pick_an_event ();do while (offset index <number of offsets){do_windowed_fft ();spectrum_fitting ();check_peak_frequency ();}calculate_Q ();}form_attenuation_field ();}Table 3.1: A pseudo code of prestack Q estimation.Figure 3.2 shows a simple test on a synthetic CMP gather with twoevents. The original Ricker wavelet has a dominant frequency of 60 Hz.Absorption is modeled by using low Q values to emphasize the absorptioneffect: the values in the two layers are 10 and 20, respectively. The actualand estimated inverse Q curves are shown in Figures 3.2b and 3.2c, respec-tively. The two curves are almost the same, because the estimated Q valuesare close to the actual values. The corresponding dominant frequency isestimated as 60.67 Hz. The agreement between the corresponding valuesshows that the method works well for ideal synthetic data. The variation ofpeak frequencies with offset is shown in Figure 3.3. Peak frequency varia-tion with offset of the first event is shown in black and the variation of thesecond event is shown in red. The remarkable decrease is, of course, dueto absorption. This example demonstrates clearly why absorption must becompensated for if resolution of prestack data is an issue.In theory, an event that can be examined at two offset locations allows usto determine the Q value of one layer. By using the peak frequency variationalong all offsets, an average Q can be obtained, thereby the effects of surfacefluctuations and random noise can be reduced, thus improving the accuracyof Q estimation.31Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation012Time (s)0 1000Offset (m)(a)012Time (s)0 0.21/Q(b)012Time (s)0 0.21/Q(c)Figure 3.2: (a) A synthetic CMP gather of two reflections with absorptionand 10% random noise. Traces have been normalized based on their maxi-mum amplitude. (b) Values of input quality factors (Q1 = 10 and Q2 = 20). (c) Computed quality factors (Q1 = 10.04 and Q2 = 20.12).32Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation0 500 1000Offset (m)01020Peak Frequency (Hz)Figure 3.3: Variation of peak frequencies of the two reflections in Figure 3.2.Peak frequencies of the first event and the second event are shown in blackand red respectively.3.3.2 Q Estimation from a Poststack TraceIn seismic data processing, stacking after normal moveout (NMO) destroysthe original absorption characteristics by summing together signals whichexperience different trajectories. In short, prestack CMP gathers before de-convolution are more suitable for attenuation analysis than stacked traces.Nonetheless, to perform de-absorption processing on poststack data is of-ten necessary; hence, it is important to estimate Q-factors from verticallyincident traces or stacked traces, if corresponding CMP gathers are not avail-able. The forgoing derived analytical approach, which computes Q valuesfrom wavelet peak frequency variation, can be applied to different seismicacquisition geometries. To detect peak frequency variation from a singletrace, we can use windowed time-variant spectral (WTVS) analysis or thecontinuous wavelet transform (CWT) analysis.Windowed Time-variant Spectral AnalysisThe windowed Fourier transform was introduced by Gabor (1946) to mea-sure localized frequency components of sound. A real and symmetric window33Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variationg(t) = g(-t) is translated by u and modulated with a frequency omega:gu,omega = eiomegatg(t -u).The resulting windowed Fourier transform of a signal f(t) isSf(u,omega) =integraldisplay infinity-infinityf(t)g(t -u)e-iomegatdt. (3.25)This transform is also called the short time Fourier transform because mul-tiplication by g(t - u) localizes the Fourier integral in the neighborhood oft = u. In the windowed Fourier transform, the window function remainsunchanged as it moves along the time axis. The resolution in time andfrequency of the windowed Fourier transform depends on the spread of thewindow in time.To examine the variation of a time-variant spectrum, the support of thewindow function needs to increase with time in order to accommodate thetime-variant seismic wavelet. The windowed time-variant Fourier transformis a useful tool for seismic time-frequency analysis.Continuous Wavelet TransformThe time-frequency localization of g in equation (3.25) can be modified bya scaling factor s. Let gs(t) = 1radicalsg(ts) represent the dilation of g(t). Thechoice of a particular scale s depends on the desired resolution trade-offbetween time and frequency.To analyze signal structures of different size at different localizations, itis necessary to use time-frequency atoms with different time supports. Thewavelet transform decomposes signals by means of dilated and translatedwavelets. Supposing a mother wavelet is a function psi(t) element L2(R) (Mallat,1999), the wavelet transform of a signal f(t) at time u and scale s isWf(t,s) =integraldisplay infinityinfinityf(t)psiparenleftbiggt -usparenrightbiggdt. (3.26)A set of wavelet bases which are wavelet functions at different scales is equiv-alent to a set of band pass filters. Orthogonal and complete wavelet basescan be constructed to perform multi-resolution analysis of signals (Mallat,1999). It is important to note that the orthogonality and completeness ofwavelet bases are not required for the application to find out the peak fre-quencies of signals.34Chapter 3. Quality Factor Estimation from Seismic Peak Frequency VariationBoth windowed time-variant spectral analysis and the continuous wavelettransform can be used for time-frequency analysis (Zhang, 2000). The fol-lowing presents the application of WTVS analysis for Q estimation from asingle trace.Quality Factor Estimation from WTVS AnalysisTo analyze the absorption characteristic of a trace in a poststack seismicsection, we use WTVS analysis because WTVS can reveal the variation offrequency content with time. If a trace is assumed to be at zero offset,from the peak frequency variation detected in a WTVS, Q values can becalculated using equation (3.13). For a seismic section, Q can be estimatedin many CMP locations. By means of horizontal interpolation, a Q profilecan be constructed. The Q profile is like a velocity profile, it can be used ininverse Q filtering and also in Q-migration.Figure 3.4(a) shows a real seismic trace. Figure 3.4(b) shows a windowedtime-variant spectrum of this trace. This WTVS gives a clear indication ofthe trend of the spectral variation. Picking the peak amplitude ridge fromWTVS and fitting the ridge with a piece-wise straight line in the coordinatesof f versus t, one can calculate a quality factor from each line segment usingequation (3.13) for the main events, as shown in Figure 3.4(c).The WTVS computation initially applies a narrow window to the inputtrace and calculates the conventional Fourier spectrum of the windoweddata. The window is then translated and widened successively with timealong the trace and the Fourier transform is calculated for each new positionof the window.3.4 Summary and DiscussionsThis chapter has presented an analytical formula which calculates a Q valuefrom the spectral peak frequency variation of a seismic wavelet. The idea isillustrated using the spectrum of a Ricker wavelet. The approach can be ap-plied to prestack CMP gathers or to postack traces for Q-factor estimation.This approach is superior to the traditional spectral ratio method becauseit only uses the information of frequency variation caused by absorption,instead of using amplitude information which is affected by many factors.Seismic absorption properties are ideally expected to be estimated fromprestack CMP gathers, and absorption compensation is ideally implementedon prestack seismic data. The reason is that stacking distorts the frequencyinformation of the original seismic data.35Chapter 3. Quality Factor Estimation from Seismic Peak Frequency Variation(a) (b) (c)Figure 3.4: Q estimation. (a) is a real seismic trace. (b) is its windowedtime-variant spectrum. The picked peak frequency points are marked bycircles and connected by straight lines. (c) is the inverse of estimated qualityfactors.36Chapter 3. Quality Factor Estimation from Seismic Peak Frequency VariationUnderground lithology is characterized by its velocity, density, and Q-factor. While the Q-factor does not affect the arrival times of reflections, itdoes affect amplitude and the signal?s frequency content. Extracting velocityinformation from CMP gathers is a common practice. we have devised ananalogous approach to estimate the absorption character from the variationof spectrum with both time and offset. The equations that compute Q-factors from observed peak frequency variation are derived based on thereasonable assumption of a Ricker-like amplitude spectrum.Absorption effects as a function of offset may be mistakenly interpretedas AVO phenomena. Without prestack Q-factor compensation, the reflec-tion amplitudes may appear to decrease with offset. Compensating for ab-sorption in prestack data can improve the resolution of seismic data and isimportant for AVO analysis. Figure 3.2 clearly shows that signals get atten-uated with traveltime and offset. Absorption compensation, based on theunderlying physics of the absorbing process, lets us compute a section withincreased resolution. Subsequently, a stationary source wavelet can be usedto improve the resolution by means of deconvolution. It is superior to con-ventional deconvolution techniques which rely on an adaptive formulation ofwavelet estimation. Only when Q values are available, can prestack inverseQ-filtering be applied to compensate for attenuation on prestack seismicdata.It is advisable to extract Q-factors from prestack seismic data that havenot been processed by frequency filtering methods (e.g. deconvolution) thatviolate the assumptions of the method. Prestack processing provides an-other dimension ?offset? that increases the information available for signaland noise separation. Q compensation performed in the prestack domaincan provide a balanced signal spectrum over offset and over time as well,improving the vertical and lateral resolution of the final image.Although seismic absorption estimation is ideally implemented on a CMPgather, the introduced method applies to seismic data with different acqui-sition geometries, wherever, the peak frequency of a seismic wavelet canreliably estimated.For poststack data, windowed time-variant spectral analysis or the con-tinuous wavelet transform can be used for time-frequency analysis. In ab-sorptive media, dispersion causes the width of a seismic wavelet to increasewith time. To allow the window to enclose an entire seismic wavelet at eachtime location, WTVS analysis uses a time-variant window function, i.e., thewidth of the window function increases with time.37Chapter 4Seismic Absorption Analysisfor Reservoir DescriptionIn reservoir description, seismic attenuation property, Q, is ideally estimatedat all prospective layers defined by velocities or acoustic impedance. Thischapter introduces a technique called reflectivity-guided seismic attenua-tion or Q analysis (RGQA) which is based on the idea that Q values canbe calculated from the peak frequency variation of a seismic wavelet. Wefirst estimate peak frequencies at a CMP location, then correlate the peakfrequency with sparsely-distributed reflectivities, and finally calculate a Qcurve from the peak frequencies at the layer interfaces using an analyticalequation. The peak frequency is estimated from the prestack CMP gatherusing peak frequency variation with offset (PFVO) analysis which is similarto amplitude variation with offset (AVO) analysis to reduce the stackingeffect on the waveform. PFVO analysis also helps to exclude the effect ofevent tuning on frequencies. The estimated Q section has the same layerboundaries as the acoustic impedance or other layer properties and can beeasily interpreted for the purpose of reservoir description.4.1 IntroductionSeismic attenuation property, usually represented by the quality factor Q, isan important parameter for characterizing rock types. Velocity and atten-uation together define the propagating behavior of seismic waves throughvisco-elastic media. One can expect that attenuation properties estimated inlayers defined by velocities or acoustic impedance would be a useful seismicattribute for reservoir description.This chapter first reviews techniques used for seismic attenuation anal-ysis in reservoir description, and then introduces a technique to estimateQ values based on the variation of zero-offset seismic peak frequencies at aCMP location. In order to avoid the distortion in seismic waveforms causedby stacking and to exclude frequency changes caused by event tuning, the38Chapter 4. Seismic Absorption Analysis for Reservoir Descriptionpeak frequency and its derivative along offset are estimated from a CMPgather using peak frequency variation with offset (PFVO) analysis. Reflec-tivities are used to define the locations where peak frequencies are requiredand the layers where Q values will be computed. The technique is calledreflectivity-guided seismic Q analysis. It operates on prestack CMP gathers,while in the meantime, it uses the reflectivity (or the interface of acousticimpedance) information, which is usually obtained from stacked seismic sec-tion. The estimated Q section can be used as an attribute which helps todelineate prospective hydrocarbon reservoirs.4.2 Review of the Methods for SeismicAttenuation AnalysisThere are mainly three categories of techniques to estimate seismic attenua-tion properties for hydrocarbon reservoir description. They are iso-frequencyanalysis, direct Q estimation and attenuation or, Q tomography.4.2.1 Iso-frequency AnalysisThis technique is regularly implemented using spectral decomposition to de-compose a seismic trace into a frequency gather. By examining iso-frequencysections, one can determine frequency dependent energy variations, amongwhich some could be related to hydro-carbon accumulation (Castagna et al.,2003). It is difficult to set a suitable frequency interval for iso-frequencyslices so that the low-frequency shadows may be properly identified. In the-ory, there are numerous mechanisms that might be responsible for the lowfrequency shadows beneath hydrocarbon reservoirs, as discussed by Ebrom(2004). It is difficult to determine which mechanisms are, in fact, thosethat introduce the first-order effects. Based on these observations, it mayperhaps be concluded that iso-frequency analysis is more a qualitative thanquantitative approach.In order to illustrate the above discussion, Figure 4.1 shows a frequencyshadow caused by attenuation, whereas Figure 4.2 shows low-frequencyshadows caused by attenuation and thin-layer tuning. In Figure 4.1, theonly absorptive layer is between 0.4s and 0.6s. The high frequencies areattenuated after the seismic wavelet traveling through the absorptive layer.The frequency decay is irreversible, because there is no physical mechanismto boost the frequency except the interference among waves. In Figure 4.2,there are two closely spaced (0.008s) reflection coefficients at around 0.6s.39Chapter 4. Seismic Absorption Analysis for Reservoir Description00.51.0Time (s)0 0.031/Q(a)00.51.0Time (s)0(b)00.51.0Time (s)30 60 90Frequency (Hz)00.010.020.03(c)Figure 4.1: Frequency decay caused by absorption. (a) is an attenuationmodel, (b) is the synthetic trace and (c) is the spectrum of frequency de-composition.The wavelets from the two layers merge and form a composite wavelet oflower frequency, as shown in Figure 4.2(b). In the time-frequency spectrumof Figure 4.2(c), the composite wavelet appears as a low frequency shadow.This frequency translation is, however, only apparent in that it does notbelong to the propagating wavelet. These two figures clearly demonstratethat the low-frequency shadow itself is not accurate enough to define anabsorptive zone.4.2.2 Direct Q EstimationConceptually, attenuation can be fully defined by the quality factor, Q.Q can be extracted from seismic traces or from sonic logs, and represents astraightforward way to estimate seismic attenuation for reservoir description.A number of techniques have been introduced to estimate Q in reser-voir applications. Dvorkin and Mavko (2006) calculated maximum Q valuesbased on high and low frequency P-wave bulk moduli. The extracted Qcurve from sonic logs confirms that the gas hydrate area has strong at-tenuation. A case study presented by H?bert et al. (2005) explains howseismic attenuation and relative acoustic impedance help to identify the gasreservoir and exclude the AVO bright spot caused by Cretaceous shales.Parra et al. (2006) extract intrinsic attenuation from the head P-wave of40Chapter 4. Seismic Absorption Analysis for Reservoir Description00.51.0Time (s)0(a)00.51.0Time (s)0 0.031/Q(b)00.51.0Time (s)0(c)00.51.0Time (s)30 60 90Frequency (Hz)00.020.04(d)Figure 4.2: Frequency decay caused by thin-bed tuning and absorption.(a) shows a reflectivity function with closely-spaced coefficients. (b) is thecorresponding absorption property. (c) is the synthetic trace and (d) is itsspectrum of time-frequency decomposition.a full-waveform sonic log using a spectral ratio approach. The obtained Qlog, together with other well logs, helps to discriminate between anomalieswhich are associated with lithology and those associated with oil and gassaturation. Singleton et al. (2006) introduce a Q estimation method usingthe Gabor-Morlet spectral ratio of an absorption compensated trace to theoriginal seismic trace. The absorption compensation is achieved by spatialspectral balancing.4.2.3 Attenuation TomographyThis approach attempts to extract attenuation structures from seismic databy minimizing the misfit between model derived data and the observations.Quan and Harris (1997), as well as Plessix (2006), combine traveltime to-mography and attenuation tomography together to invert velocity and atten-uation from cross-well transmitted data. Seismic attenuation is modeled bymeans of the variation of the centroid frequency of a seismic wavelet. Prattet al. (2003) use waveform tomography which is based on a visco-acousticwave equation. This method helps to find a layered hydrate distribution withattenuation anomalies. Hicks and Pratt (2001), Watanabe et al. (2004) andGao et al. (2005) also adopt waveform tomography to invert for attenuation.41Chapter 4. Seismic Absorption Analysis for Reservoir DescriptionThe three classes of techniques for attenuation analysis listed above, haveall experienced a degree of success in application. However, the problem ofhow to obtain an interval intrinsic Q profile from seismic reflection data, orhow to obtain a Q that can be easily correlated with acoustic impedanceand/or other elastic parameters, remains a challenge in reservoir character-ization.4.3 Reflectivity-guided Q Analysis, RGQA4.3.1 TheoryThe following four points outline the theory behind the technique that wehave named reflectivity-guided Q analysis, RGQA in abbreviated form:1. We assume a frequency-independent Q which is the only parameterrequired to fully describe the attenuation properties of a layer. ThisQ is related to lithology, porosity, pore fluid, saturation etc.. Togetherwith other available geological and geophysical information, Q providesone more attribute for reservoir characterization. Like velocity, Q is aninterval property, which can be derived from the reflections occurringat layer boundaries.2. A Q curve can be estimated from the frequency variation of a zero-offset trace, which describes a vertically traveling seismic wavelet.Its time-variant spectrum describes the attenuation property of themedium. Frequency variations can be simply described by the peakfrequency shift.Assuming the peak frequencies of a wavelet at time t1 and t2 are fp1and fp2 respectively, if the amplitude spectrum of the seismic waveletis like that of a Ricker wavelet (the second derivative of a Gaussianfunction), the inverse Q value can be calculated using (from equation(3.13))1Q =2(f2p1 -f2p2)pi(t2 -t1)fp2f2p1 . (4.1)3. The zero-offset trace is not naturally available. Seismic stacking sumsup traces in a CMP gather to reduce noise and improve overall dataquality. Furthermore, stacking distorts the frequency variation patterncaused by seismic attenuation. The seismic wavelet at a CMP locationin a stacked seismic section is different from what is observed on thetrace of vertical-incidence.42Chapter 4. Seismic Absorption Analysis for Reservoir DescriptionFigure 4.3: Flow diagram of reflectivity-guided Q analysis.4. Required peak frequencies of seismic wavelets at zero-offset can beestimated by means of modeling the frequency variation along offsetin a CMP gather. We call this technique PFVO analysis.Figure 4.3 shows the flow diagram of the RGQA algorithm. Attenuationanalysis is implemented on a prestack CMP gather and requires sparselydistributed reflectivities for interface selections.Layer interface information can be obtained from acoustic impedanceinversion or sparse reflectivity inversion (Oldenburg et al., 1983; Ulrych andSacchi, 2005). Nowadays, piece-wise constant acoustic impedance informa-tion is generally available for reservoir description (Latimer et al., 2000).RGQA employs poststack acoustic information to help the extraction ofattenuation information from the prestack data.Seismic data contain information concerning the interfaces between rocklayers along the path of propagation. In general, the information concerningthe layers themselves is inferred from that describing the interfaces. Acous-tic impedance and AVO inversion are two instances where this strategyis used to infer layer properties. Since layer boundaries defined by reflec-tion coefficients are also representative of the boundaries of the attenuatinglayers, RGQA derives rock attenuation properties pertinent to the layersthemselves.In order to illustrate the general idea behind RGQA, we use the modelin Figure 4.1 once again. Peak frequencies at four reflectivity locations43Chapter 4. Seismic Absorption Analysis for Reservoir Description00.51.0Time (s)20 30 40 50Frequency (Hz)(a)00.51.0Time (s)0 0.031/Q(b)Figure 4.4: Peak frequencies at layer boundaries. (a) is the peak frequenciesfrom the trace in Figure 4.1 (b). (b) is the extracted inverse Q curve.are picked from the spectrum in Figure 4.1(c) and these frequencies aretransformed to 1/Q values using equation 4.1. Figure 4.4 shows the peakfrequencies at these locations, as well as the estimated 1/Q curve. For thissimple case, RGQA works perfectly. This method, when applied to theevent tuning model shown in Figure 4.2(b), will fail however. The frequencydecay or rise caused by event tuning can be excluded by examining the peakfrequency variation along offset (PFVO) on CMP gathers.4.3.2 Peak Frequency Variation with Offset (PFVO)AnalysisConsidering that NMO correction followed by stacking could distort theshape of the wavelet and thus obscure the sought after absorption effects,we suggest a method to estimate peak frequencies at zero offset directlyfrom prestack CMP gathers based on the frequency variation along offset.This multichannel peak frequency analysis strategy makes the estimationalgorithm more robust to random noise and crossing events in a CMP gather.PFVO analysis follows the same implementation procedure as AVO (am-plitude variation with offset) (Aki and Richards, 1980). It fits the peak fre-quency along offset linearly. The intercept, Pf, can be transformed into Qvalues, and the gradient, Gf, serves the fundamentally important, and thusfar unresolved, task of separating Q into real attenuation and event tuning44Chapter 4. Seismic Absorption Analysis for Reservoir Descriptioneffects. PFVO operates on a CMP gather without NMO correction, becauseNMO stretch at far offset distorts the frequency spectrum of reflections.There are two main issues in Q estimation using frequency variationwith offset. The first is how to determine the peak frequency at each layerboundary. The second is how to determine Q from the peak frequencyvariation. Regarding the first issue, one general and reliable way is to usespectral decomposition, i.e. to perform time-frequency analysis on eachtrace, and use the instantaneous amplitude spectrum to determine a peakfrequency. Regarding the second issue, we fit the peak frequencies of anevent at different times by using a formula like equation (4.1) depending onthe spectral shape of a seismic wavelet.For a Ricker-like amplitude spectrum, which is the assumed spectralshape in our work, the peak frequency variation with time follows the fol-lowing equation (also equation (3.12))fp = f2mbracketlefttpbracketleftbtradicalBiggparenleftbiggpit4Qparenrightbigg2+parenleftbigg 1fmparenrightbigg2- pideltat4Qbracketrighttpbracketrightbt.where fm is the peak frequency of a wavelet at its initial state (time = 0)and fp is the peak frequency after travel time t. It can be observed fromthis equation that if the peak frequency at time t - deltat is fp0 after a timeinterval deltat, the peak frequency isfp = fp0bracketlefttpbracketleftbtradicalBiggparenleftbiggfp0pideltat4Qparenrightbigg2+ 1 - fp0pideltat4Qbracketrighttpbracketrightbt. (4.2)In equation (4.2), the term fp0pideltat4Q is almost always less than 1 for small deltatand weak attenuation (large Q). In a further step, we expand the squareroot via a Taylor series and keeping only the first term we obtainfp = fp0bracketleftBigg1 - fp0pideltat4Q + f2p0pi2deltat232Q2bracketrightBigg.In general, peak frequency fp varies with the difference of arrival times non-linearly. For weak attention, the Q2 term can be omitted from the aboveexpression, yielding the following approximation:fp = Pf -Gfdeltat, (4.3)where Pf = fp0, Gf = pifp04Q . For an event in a CMP gather, deltat is thenormal moveout for a particular offset. For an event starting from two-way45Chapter 4. Seismic Absorption Analysis for Reservoir Descriptiontraveltime t0 at zero-offset, deltat is a function of t0, RMS velocity and offset,i.e.deltat = deltat(t0,vrms,offset).Equation (4.3) shows that for a seismic wavelet with a Ricker-like amplitudespectrum, the peak frequency decay of a reflection with the difference ofarrival times in absorptive media can be simulated by a linear relation. Theintercept, Pf, of the straight line is the peak frequency at zero offset, andthe gradient, Gf, is related to seismic attenuation.A first glance at equation (4.3) may give the impression that Q couldbe computed directly from Gf. Actually, Gf can only provide root-meansquare attenuation information, (Zhang and Ulrych, 2002), rather than therequired interval Q?s. Instead of providing attenuation information directly,Gf indicates whether there is peak frequency decay along offset. Therefore,we use Gf as an indicator which tells us whether the frequency decay atzero offset is caused by attenuation or by some extrinsic mechanism, e.g.thin-bed scattering. In fact, in our method, Q is calculated from the Pf?sonly.4.3.3 Implementation StepsThe technique is implemented in six steps as outlined below:1. Build a CMP super-gather to suppress the effect of random noise.A trace in a super-gather is a stacked trace of several neighboringtraces (Ostrander, 1984). The CMP gather used for Q-factor estima-tion needs to be free of NMO stretch, so that if normal moveout isapplied when forming a CMP super-gather, reverse NMO needs to beapplied afterwards.2. Do time-frequency analysis on each trace using a continuous wavelettransform or short window Fourier transform. More sophisticatedspectral decomposition like exponential pursuit decomposition (Castagnaand Sun, 2006) may help to improve the resolution of the time fre-quency spectrum. Resolution and computation are two factors for thechoice.3. Compute the instantaneous amplitude envelope using Hilbert trans-form (Taner et al., 1979) to get a smooth amplitude spectrum.4. Examine for the highest amplitude among frequencies at interface lo-cations. The peak frequency is the frequency corresponding to themaximum amplitude.46Chapter 4. Seismic Absorption Analysis for Reservoir Description5. Fit peak frequencies with normal moveout linearly using l1 norm toget Pf and Gf, and edit Pf by removing the values corresponding tonegative or very small Gf?s.6. Calculate Q from the effective Pf values.Theoretically, peak frequencies at two offset locations can fully definethe trend of peak frequency variation of a reflection. By using multi-channelinformation along offset, the algorithm is robust to random noise and theeffects of crossing events.Forming super-gathers, spectral decomposition and l1 linear fitting arethe main parts of the whole PFVO computation. Compared to traditionalAVO analysis, PFVO analysis does fewer linear fitting operations, becauselinear fitting is only performed at layer interface locations, and the extracomputation is efficient reverse NMO and the spectral decomposition. Thewavelet transform can be regarded as a series of band-pass filters. Becausethe filters can be precomputed before the convolutions, the computation loadof spectral decomposition is limited. PFVO analysis on a CMP gather iscomputationally comparable to traditional AVO analysis. The reflectivity-guided Q analysis is practical even for large 3D prestack data volumes.4.4 Numerical ExperimentsFirst, let?s look back at the model with event tuning and attenuation inFigure 4.2 to see whether the method of the reflectivity-guided attentionanalysis can help to extract actual attenuation information.Figure 4.5(a) shows a synthetic CMP gather using the reflectivity model(Figure 4.2(a)) and the attenuation model (Figure 4.2(b)). The data at thenear offsets on the CMP gather are purposely zeroed to simulate real seismicdata acquisition, where zero-offset data are not recorded. Pf and Gf are ex-tracted from the CMP gather. Pf and Gf/Pf are displayed in Figure 4.5(b)and Figure 4.5(c) respectively. The reason for displaying Gf/Pf rather thanGf is that the value of Gf/Pf falls in a relatively narrow range. Gf/Pf iscalled the relative gradient. The low frequency caused by tuning events at0.6s is excluded from the Q computation because it does not decay withoffset (Gf /Pf approxequal 0). The estimated Q curve is shown in Figure 4.5(d) whichmatches the original Q model (Figure 4.2(b)).The technique is also tested on a real 2D dataset. Besides the prestackCMP gathers, an acoustic impedance or reflectivity section is required toprovide the layer-boundary information. Figure 4.6 shows the test on a47Chapter 4. Seismic Absorption Analysis for Reservoir Description00.51.0Time (s)500 1500Offset (m)(a)00.51.0Time (s)20 30 40 50Frequency (Hz)(b)00.51.0Time (s)0 2 4Relative Gradient(c)00.51.0Time (s)0 0.031/Q(d)Figure 4.5: Peak frequency variation with offset. (a) is a CMP gatherbuilt from the model shown in Figure 4.2. (b) is the corresponding peakfrequencies at zero offset. (c) is the relative gradient of peak frequenciesalong offset, and (d) is the estimated Q curve.CMP gather, the super-gather having been formed in a manner similar tothat used in AVO analysis (Ostrander, 1984), except with reverse NMOapplied to remove NMO stretch. Figure 4.6(b) shows the peak frequencygather. Figure 4.6(c) shows the peak frequencies at zero offset. Figure 4.6(d)illustrates the relative gradients of peak frequency along offset where posi-tive values indicate frequency decay along offset. Figure 4.6(e) is the sparsereflectivity obtained using the Bayesian inversion method described in Ul-rych and Sacchi (2005). Figure 4.6(f) is the estimated Q curve which issupposed to be piece-wise constant. In practice, the linear fitting step inPFVO analysis needs only to be done at interface locations. The continuouscurves in Figure 4.6(c) and Figure 4.6(d) are computed at all samples. Thereason for doing so is only for the purpose of display, and this reason alsoexplains why Figure 4.7(a) shows a peak frequency section.For this 2D dataset, the calculated peak frequencies at all CMP locationsare shown in Figure 4.7(a) and the estimated attenuation section is shownin Figure 4.7(b). The attenuation section displays the values of Q-1. Weuse the inverse of Q since it lies in the range of (0,1) and is ideal for thepurpose of display. The blue background indicates small Q-1 values whichimplies smaller absorption than the red areas.48Chapter 4. Seismic Absorption Analysis for Reservoir Description0.51.01.52.02.5Time (s)500 1000Offset (m)(a) (b)0.51.01.52.0Time (s)10 20Pf(c)-4-2 0 2Gf/Pf(d)-0.20 0.2Reflectivity(e) (f)Figure 4.6: A real data example. (a) is a CMP gather. (b) is the correspond-ing non-stationary distributed peak frequencies. (c) is the peak frequencycurve at zero offset. (d) is the relative gradient curve of peak frequencies.(e) is the reflectivity function and (f) is the estimated Q curve.49Chapter 4. Seismic Absorption Analysis for Reservoir Description(a)(b)Figure 4.7: Real data experiment. (a) is the peak frequency section and (b)is the attenuation section (Q-1).50Chapter 4. Seismic Absorption Analysis for Reservoir Description4.5 Discussion and ConclusionIn this chapter, we have introduced a technique of reflectivity-guided Qanalysis. This approach uses poststack sparse reflectivity (the interfacesdefined by acoustic impedance) to help pick the layer boundaries on prestackCMP gathers, where peak frequencies are selected to compute Q values.This technique basically provides a methodology that allows Q values to beestimated at all prospective layers.The obtained Q section may be computed with the same resolution asother attribute sections. Therefore, the absorption property can be easilyincorporated into acoustic impedance interpretation. We do not suggesthere the use of attenuation as a direct hydro-carbon indicator but ratheras an attribute which could be combined with other parameters to give abetter understanding of rock properties.We have also introduced an AVO-like, peak frequency estimation schemecalled PFVO. By examining the peak frequency shift along offset, PFVOattempts to estimate a reliable peak frequency at zero offset, and at thesame time to separately identify the frequency change caused by thin-bedscattering and that due to absorption. PFVO analysis operates on CMPgathers without NMO to avoid the effect of NMO stretch. If a non-stretchNMO technique (Perroud and Tygel, 2004; Masoomzadeh et al., 2004) wereavailable, PFVO could be done after NMO correction. AVO analysis canbe performed on an image gather from prestack time migration. WhetherPFVO analysis will benefit from being applied to a common image gather isnot clear at present. How prestack migration affects the frequency contentalso needs to be investigated.The accuracy of our method depends on the accuracy of time-frequencyanalysis and on the validity of the equation used to calculate Q. Further,we have assumed that frequency translation is caused by absorption andevent tuning only. There may certainly also exist other factors which affectthe frequency content of seismic data, in which case, such need be identifiedand quantified. The success of the attenuation analysis is heavily dependenton the accuracy of the interface information, generally represented by thereflectivity or the gradient of acoustic impedance. Without accurate subsur-face boundary information, the attenuation property obtained from seismicdata will be difficult to interpret.51Chapter 5Seismic AbsorptionCompensation: a LeastSquares Inverse SchemeThis chapter addresses the problem of instability that plagues conventionalinverse Q filtering. The de-absorption problem is formulated as an inverseproblem in terms of least squares and regularization is imposed by meansof Bayes? Theorem. The solution is iterative and non-parametric, and itreturns a reflectivity function that has been constrained to be sparse. Theinverse scheme is tested on both synthetic and real data and the resultsobtained demonstrate the viability of the approach.5.1 IntroductionThe aim of seismic processing is to obtain a high resolution image of the sub-surface. Many steps are involved in the attempt to achieve this objective.One such step consists, or most certainly should consist, of the compensationfor the ubiquitous absorption of seismic energy in the earth, a step referredto as de-absorption (McGinn and Duijndam, 1998). To a first approxima-tion, the earth is a linear attenuator, the effect of which is to attenuate thehigher frequencies of the seismic signal. The mechanisms and the relevantmathematical details have been described in Chapter 2 and in the literaturementioned there such as Futterman (1962), Kjartansson (1979), Aki andRichards (1980), Wang and Guo (2004) and Brillouin (1960).The most serious setback to the application of de-absorption is the inher-ent instability of deconvolution filters which, unless very carefully designed,unavoidably amplify high frequency noise. Many authors have consideredthis problem (Robinson, 1979; Hale, 1991; Hargreaves and Calvert, 1991;Bickel, 1993; Varela et al., 1993; Wang, 2002; Margrave et al., 2003; Wang,2003, 2006). Although, as a result of the effort of the aforementioned au-thors, a variety of algorithms are now available to the industry, it is probably52Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Schemesafe to say that all are very sensitive to the presence of additive noise. Theapproach introduced here to this noise associated problem is by means ofthe application of least squares (LS) inversion. This approach attempts toobtain a solution to the inverse problem of de-absorption by minimizing themisfit between observed data and theoretically modeled data. Thus, theinstability of applying an inverse operator to the random noise in the data,is avoided.The LS objective function may be constructed in various ways, but themanner which we prefer is based on a probabilistic rationale (Tarantola,1987; Ulrych and Sacchi, 2005), where the measured and modeled data(which become the solution of the inverse problem) are considered as re-alizations of random processes.There are two main advantages to the LS approach proposed here ascompared to previously published methods. One advantage is that, if theproblem is formulated as under-determined, the type of its solution can becontrolled as desired. The second advantage is that the LS approach usesonly the forward operator, and does not require the unstable inverse opera-tor. In other words, the formulation of de-absorption as an inverse problem,allows a stable solution to be achieved. As pointed out by Claerbout (1985),earth absorption might be compensated for by amplifying high-frequency en-ergy during downward continuation, shifting a wavefield using an absorptionoperator. This very interesting topic of migration with absorption compen-sation will be investigated in Chapter 6.This chapter is arranged as follows: First, we discuss very briefly, crit-ical issues in de-absorption, which include Q information, random noiseand methods for absorption compensation. Secondly, we describe a de-absorption approach using a LS optimization scheme based on probabilisticinverse theory. Finally, we discuss the results of both synthetic and realdata examples using the LS algorithm.Absorption in the earth causes both the attenuation and dispersion ofthe reflections recorded on the surface. Both phenomena lead to a loss ofresolution with resulting loss of information concerning targets of poten-tial interest. This is particularly severe for small , but hopefully profitable,targets at depth. In seismology, absorption is characterized by the qualityfactor, Q, and many different models have been developed which describethe associated attenuation-dispersion relationship. De-absorption, the pro-cess by which the undesirable effects of absorption are attenuated, is aninherently unstable process when applied in the conventional manner. Thereason for this is that the model of absorption in the earth, customarilyadopted, is one where the seismic trace is convolved with a time varying53Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Schemefilter. Conventional practice is to deconvolve the trace by means of the in-verse of this filter. As is well known, deconvolution is very sensitive to highfrequencies (Varela et al., 1993) and so, therefore, is de-absorption.The LS approach to de-absorption that is presented here differs quiteradically from the deconvolution techniques in customary use. The problemis cast in terms of probabilistic inference where a priori information is usedto shape the solution to our view of the earth. This view is that reflectivitiesare sparse.5.2 Issues of Absorption Compensation5.2.1 Method ChoiceThere are many methods available to implement absorption compensation,the most common ones are time-variant spectral whitening, time-variant de-convolution and inverse Q filtering (Yilmaz, 2000). Spectral whitening triesto compensate for absorption caused amplitude decay by applying a gainfunction at different frequency bands. Time-variant deconvolution imple-ments absorption compensation by means of using a time-variant waveletin a moving time window. Q information is not used explicitly, because Qinformation is hard to obtain. In order to compensate for attenuation moreaccurately and deterministically, Inverse Q filtering or other Q dependentcompensation techniques are superior.5.2.2 Inaccurate Q EstimationBecause the estimated values of Q-factors are always different from realvalues, Q could be under or over estimated. Let us define the actual qualityfactor as Qo, its estimation as Qe and the difference as deltaQ, whereQe = Qo + deltaQ. (5.1)The inaccuracy deltaQ negationslash= 0 in the estimation will affect the amplitude andphase of the compensated signal in a manner outlined below.Effects on the AmplitudeIf Qe is used to compensate for the amplitude attenuation resulting from Qo(equation (2.20)), after the correction, the amplitude isAc(omega) = expbracketleftbigg omegadeltaz2v(omega)-deltaQ(Qo + deltaQ)Qobracketrightbigg.54Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse SchemeIf absorption is under-estimated or deltaQ > 0, Q compensation does notrecover the high frequency energy completely. If absorption is over-estimatedor deltaQ < 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 inputzero-phased absorption response. Original reflection coefficients are threeunit 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. Thecompensation 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, andthus, 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 resultlooks like strong high-frequency noise.Effect on the PhaseIf Qe is used to compensate for the phase distortion resulting from Qo (equa-tion (2.21)), after the correction, the phase isphic(omega) = omegadeltazpiv(omega) -deltaQ(Qo + deltaQ)Qolnparenleftbigg omegaomega0parenrightbigg.If the actual value and estimated value are very close, i.e. |deltaQ| approxequal 0, afterQ compensation, phic(omega) approxequal 0, which means that the phase has no distortion.Pure phase compensation will make the original minimum-phase absorptionresponse become zero-phase. This is an ideal case. If deltaQ > 0, i.e., Q isover-estimated or absorption is under-estimated, phic(omega) < 0. This means thearrival time of this frequency is delayed. The higher the frequency the morethe wave is delayed. If deltaQ < 0, i.e., Q is under-estimated or absorption isover-estimated, phic(omega) > 0. This means the arrival time is advanced. Thehigher the frequency, the sooner the wave component arrives.Figure 5.2 shows phase corrected responses using under-estimated andover-estimated Q values respectively. The under-corrected signal is minimum-phase, and the over-corrected signal is maximum phase.5.2.3 Random NoiseAs already mentioned before, inverse Q filtering is inherently unstable sincethe inverse operator will boost high frequency noise. To ensure that noiseis not unnecessarily amplified, it is important to design the inverse operator55Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme0 500 1000 1500 2000Time (ms)00.10.20.3(a) Zero-phased impulse response with absorption.0 500 1000 1500 2000Time (ms)00.20.40.60.8(b) Amplitude compensation using actual Q = 30.0 500 1000 1500 2000Time (ms)00.10.20.3(c) Amplitude compensation using big Q = 40.0 500 1000 1500 2000Time (ms)02468101214(d) Amplitude correction using small Q = 25.Figure 5.1: Experiments on amplitude compensation using different Q val-ues.56Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme0 500 1000 1500 2000Time (ms)00.10.20.3(a) Impulse response with absorption.0 500 1000 1500 2000Time (ms)00.10.20.3(b) Phase correction using Q = 30.0 500 1000 1500 2000Time (ms)00.10.20.3(c) Phase correction using Q = 20.0 500 1000 1500 2000Time (ms)00.10.20.3(d) Phase correction with Q = 40.Figure 5.2: A synthetic trace with three spikes, Q = 30, and its phasecorrection using different Q values.57Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme0 10 20 30 40 50 60 70 80Frequency (Hz)102030405060AmplitudeFigure 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 usea band-limited version of the inverse operator, i.e., to replace the amplitudecompensation operator A-1(omega) by W(omega)A-1(omega), where W(omega) represents alow-pass filter. This idea is illustrated in Figure 5.3, where an exponentialamplitude compensation operator is approximated by a Gaussian functionafter a certain frequency limit. Of course, as can be well imagined, such alow-pass operation causes an undesirable loss of resolution. The least squaresmethod being introduced is a way to achieve complete compensation whilestill keeping the stability of the operation of inverse Q filtering.5.3 Time-variant Deconvolution as an InverseProblem5.3.1 Noise and DeconvolutionThe instability issue is not unique to the absorption compensation problem.It is, in fact, ubiquitous in the application of deconvolution in general. Inseismic data processing, the recorded data in the one-dimensional case canbe modeled asx(t) = w(t) asteriskmath r(t) + n(t), (5.2)58Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Schemewhere x(t) , w(t), r(t), and n(t) represent the recorded trace, a seismicwavelet, the earth?s reflectivity and random noise, respectively. Seismic datadeconvolution aims at obtaining the reflectivity from the observed trace. Theestimated 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) asteriskmath w(t) = delta(t) ), the deconvolution result is?(t) = r(t) + f(t) asteriskmath n(t).Depending on the signal to noise ratio, the effect of the term f(t) asteriskmath 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 directlyinverting with the known or estimated seismic wavelet. The resulting filtersare called optimum Wiener filters and are well elaborated by Robinson andTreitel (2002). We take a probabilistic approach to the LS solution here, anapproach which is believed has a flexibility which leads to improved results.Specifically, we formulate de-absorption as an inverse problem and use aCauchy-Gauss objective function (Sacchi et al., 1998) to arrive, iteratively,at the desired solution.5.3.2 Bayes Probabilistic InferenceProbabilistic inference is a very powerful approach to the ubiquitous inverseproblem (Tarantola, 1987; Ulrych and Sacchi, 2005). It is common to con-sider the measured data d (where vectors are indicated by bold symbols) asuncertain. That is, true data do exist but are unknown. Measured data canthen be considered as random variables whose expectation, in the ensemblesense, is the true value. Traditionally, the assumption is made that the ob-served data are random variables with a Gaussian distribution. This leadsto the well known chi2 test for goodness of fit and to a l2 norm solution.Let us represent the solution to the inverse problem by the model vectorm. m is not unique, in the sense that an infinity of such models may be foundthat fit the data. The reason for this, of course, is that the data, even if noisefree, can never describe the continuous model completely. Our problem isunder-determined. The problem may indeed be posed as an over-determinedone, if we so choose. Wiener inversion is one such possibility. However, webelieve, having been well schooled by the late Edwin Jaynes, that an inverseproblem posed in such manner, is incorrectly posed (Jaynes, 2004). A dis-cretized inverse problem can never be unique. The over-determined solution59Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Schemereturns the ?smallest model?, i.e., that model whose energy is most dis-tributed in model space. When estimating deconvolution models describingreflectivities of the earth, such models do not appear to be reasonable.An infinity of models is not a useful set. We must impose a prioriknowledge, if such exists, feeling or at least hoping, associated with oursearch. The manner of so doing lies at the heart of Bayes? theorem, and atthe heart of the regularization thus imposed. This is hardly the place for adiscussion of the underlying logic, but we feel that a few remarks are in order(for a fairly full discussion vis a`is Bayes, inversion and a prior information,please see Ulrych and Sacchi (2005)). The central point in the applicationof Bayesian logic is in how one handles prior information. We begin withBayes? Theorem, which states thatp(m|d) = p(d|m)p(m)p(d) , (5.3)where p(d|m) is the conditional probability density function(pdf) of d, thedata, given that m has occurred. This pdf is called the likelihood and isthe 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 touse, isp(x) = 12pisigmae-h x-?2sigma2i2, (5.4)where sigma2 is the variance and ? the mean of the pdf, respectively.Returning to equation (5.3), p(d) is the pdf of the observed data andoften is not a factor in the inversion. In as much as d is considered to bea random vector, m can be treated in an equivalent manner. p(m|d) is thepdf that we wish to obtain. It is the probability of obtaining m given thedata. Maximizing this pdf will obtain an estimated model that is referred toas the maximum a posteriori solution. If p(m) = 1, i.e., a uniform pdf, theformulation of the problem is exactly least squares. Now, however, comesa point of much discord, p(m). In the opinion of a group of statisticiansknown as ?frequentists?, prior probabilities cannot be assigned, they mustbe measured. Edwin Jaynes, like Laplace before him, fought passionatelyagainst this philosophy, and most certainly emerged the victor. So, withvictory on our side we press on. p(m) is the pdf that we wish to associatewith our expected model. It is our choice based on our belief. It couldcertainly be inappropriate, but who is to say that the l2 model, the onemost distributed in time and/or space, is any more so? What pdf to choosewhen dealing with de-absorption? We follow the logic of Sacchi et al. (1998)and use the Cauchy distribution, given by60Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Schemep(x) =parenleftbigg1 + x22sigma2parenrightbigg-1, (5.5)where sigma2 governs the spread of the distribution (the Cauchy pdf being onethat does not possess a theoretical variance).The reason behind this choice is, very basically, that this pdf exhibitslong tails and, consequently, is a candidate for the pdf that characterizessparse reflectivities. Sparseness is of much relevance these days (Hermann,2005; Cand`s et al., 2005). It is not a concept applicable in all cases, ofcourse. One would not model a Beethoven concerto in this manner, butlayers in the earth are, one hopes, sparsely distributed.5.3.3 Absorption Compensation and InversionWe now formulate our problem in light of the above discussion. Assume theobserved data to be contaminated by noise that is normally distributed asN(0,sigma2n), where n represents the noise vector. The Gaussian assumptionstems from the principle of maximum entropy (Ulrych and Sacchi, 2005)for a discussion and references) which, informally speaking, states that ifnothing is known, use the simplest hypothesis. In this case, the Gaussianpdf follows from the Central Limit Theorem.The conditional distribution of the data is given by (equation (5.3))p(m|d,sigman) =parenleftbigg 12pisigma2nparenrightbigg(M-1)/2e-(1/2sigma2n)||d-Gm)||22, (5.6)where M is the length of the data vector and, for our linear systemn = d -Gm,with G the coefficient or kernel matrix.Let p(m|sigmam) indicate a prior distribution for m conditional on a param-eter sigmam. From Bayes? theorem, equation (5.3), we havep(m|d,sigmam,sigman) = p(d|m,sigman)p(m|sigmam)p(d) .. (5.7)The model vector, m, is the reflectivity function, and we use a sparsenessconstraint in the inversion scheme.61Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse SchemeCauchy-Gauss ModelIn the following, we will assume that a seismic wavelet is available. It iseither obtained from check-shot survey, or from a seismic trace in somemanner. 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 sigma2m (a standardhypothesis). The joint pdf of the mk isp(m|sigmam) =productdisplaykp(mk|sigmam),where p(mk|sigmam) satisfies a Cauchy distribution given by equation (5.5).Substitution yieldsp(m|sigmam) =productdisplaykparenleftbigg1 + m2k2sigma2mparenrightbigg.By inserting this Cauchy prior and the data likelihood (equation (5.4)) intoequation (5.7) and taking logarithms of both sides, we obtainln[p(m|d,sigmam,sigman)] = -c(m) - 12sigma2n(d -Gm)T (d -Gm), (5.8)where c(m) is a constraint imposed by the Cauchy distributionc(m) =M-1summationdisplayk=0lnparenleftbigg1 + m2k2sigma2mparenrightbigg.Furthermore, denoting phicg(m) = -ln[p(m|d,sigmam,sigman)], we obtain the costfunction for the Cauchy-Gauss model asphicg(m) = c(m) + 12sigma2n(d -Gm)T (d -Gm), (5.9)where G is composed of time variant wavelets, btau(t -tau). Both t and tau arein the range of the length of a trace andG =bracketlefttpbracketleftexbracketleftexbracketleftexbracketleftbtb0(t -0)b1(t -1)...bM(t -M)bracketrighttpbracketrightexbracketrightexbracketrightexbracketrightbt.In this manner, one dimensional absorption compensation is formulatedas an inverse problem. The model, which is related to the minimum of thecost function, is the sparse reflectivity function we desire. The Cauchy-Gauss model has also been used in acoustic impedance inversion, signalinterpolation and extrapolation (Sacchi et al., 1998).62Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse SchemeImplementationThe solution of the inverse problem, formulated above, follows the algorithmoutlined by Sacchi et al. (1998). The objective function that will be min-imized is phicg(m) in equation (5.9). Taking the derivative of phicg(m) withrespect to m, we obtainpartialdiffphicg(m)partialdiffm =partialdiffc(m)partialdiffm +12sigma2n GT (d -Gm) (5.10)and partialdiffc(m)partialdiffm =1sigma2m S-1m,where S is a M ? M diagonal matrix with elements skk = 1 + m2ksigmam , (k =0,1,... ,M - 1) and is m dependent. Equating the derivative in equation(5.10) to zero, yieldsm = parenleftbiglambdaS-1 + GT Gparenrightbig GT d, (5.11)where lambda = 2sigman2m2 . This equation is nonlinear and must be solved iteratively.To construct an iterative algorithm, equation (5.11) is written asm = SGT parenleftbiglambdaIN + GSGT parenrightbig-1 d = SGT p. (5.12)p is the auxiliary vector which is obtained from the solution of the systemparenleftbiglambdaIN + GSGTparenrightbig p = d. (5.13)The iterative computation is initiated by setting the observed seismic dataas the initial model, m0, which is also used to generate the matrix S0. Ineach iteration, k, we first computep(k-1) =bracketleftBiglambdaIN + GS(k-1)GTbracketrightBig-1m,and then update the model asm(k) = S(k-1)GT p(k-1).The algorithm, although iterative, is computationally efficient since thecoefficient matrix in equation (5.13) is Hermitian-Toeplitz and is invertedusing the Levinson recursion.The processing flow for absorption compensation on a stacked seismicsection is implemented in the following five steps:63Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme1. 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 theearly 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 everyCMP location. Indeed, Q values need only be estimated for selected CMP?sand interpolated to form a Q profile.Computations required, other than those associated with the iterativealgorithm, are trivial if both Q and the seismic wavelet do not change dra-matically along the CMP direction, and the resulting LS scheme for seismicabsorption compensation is computationally practical.Numerical ExperimentsThis section illustrates results obtained using the LS approach on both syn-thetic and real data. The synthetic test is based on the assumption of knownwavelet 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 inthe actual model.Figure 5.5 imitates a section of seismic data. The signals in the five tracesshown in Figure 5.5(a) are the same, but contain different random noise. Thenoise level is 20 percent. Figure 5.5(b) is the result obtained by LS time-variant deconvolution. The rightmost trace is the original reflectivity seriesplotted as a reference. It can be observed that the reflectivity function isrecovered very well both as far as locations and amplitudes are concerned.This result shows that the Cauchy-Gauss objective function is a very viablea priori model for obtaining accurate reflectivity inversion.Figure 5.6(a) illustrates a stacked seismic section. The structures arenot complex and consist of several flat layers. Processing of the section isinitiated by first estimating the Q curve from a trace in the section usingthe WTVS method. Second, using these extracted Q?s, the dispersive phasecorrection is computed and is followed by the estimation of the minimum64Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse SchemeFigure 5.4: Convergence process of the inverse scheme. (e) is a synthetictrace. (a-d) are the inverse results after iteration 1, 2, 3 and 4 respectively.65Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse SchemeFigure 5.5: (a) is a synthetic trace with 20% noise. (b) is the result ofabsorption compensation. The rightmost trace is the original reflectivityseries plotted as a reference.66Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme(a) (b)Figure 5.6: Q compensation for a real seismic section. (a) shows a realseismic 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 avalue for the parameter lambda in equation (5.13) and the result is updated it-eratively. Following inversion, the wavelet is convolved with the extractedreflectivities. 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), canbe 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 isclear. The de-absorption result would certainly facilitate the interpretationof this seismic section.67Chapter 5. Seismic Absorption Compensation: a Least Squares Inverse Scheme5.4 ConclusionsThe 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 approachand is robust with respect to additive noise. The technique described herehas very general application. Specifically, since robust Q compensation pro-vides more accurate information concerning both the amplitude and locationof the earth?s reflectivity, hydrocarbon reservoir characterization is one obvi-ous target for the introduced method. In particular, the resulting improvedresolution may have important consequences in 4D processing where theobjective is to delineate changes in reservoir properties.Results that we have obtained on synthetic and real data confirm, webelieve, 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 takingthe advantage of the Hermitian-Toeplitz property of the kernel matrix thatdescribes the forward problem, stable inverse solutions are achieved econom-ically. The synthetic and real data results demonstrate the viability of theapproach and the pivotal role of de-absorption in improving the resolutionof the processed sections.The LS Q compensation is implemented in a trace by trace manner anddoes not take into account lateral continuity or the actual ray-paths thatare involved. Multi-channel Q compensation and poststack time migrationwith Q compensation are interesting and important topics which will bediscussed in Chapter 6.68Chapter 6Refocusing Migrated Imagesin Absorptive MediaThis chapter investigates the blurring effect on migrated images when usinga regular migration algorithm to migrate the seismic data with absorption.First, a numerical method is introduced to calculate the blurring functionsin layered media, and then, a least squares inverse scheme is used to removethe absorption blurring effect in order to refocus the migrated image. Thestability and the convergence property of the refocusing algorithm will bediscussed. Experiments on synthetic and real data will be given to showthe validity of the introduced technique to compensate for absorption aftermigration.6.1 IntroductionConventional 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) andmany others. In reality, the spatial distribution of a Q field is similar tothat of a velocity field, i.e., Q is spatially variant. The effect of attenuationon received seismic signals is related to wave-paths, so that absorption com-pensation is ideally implemented based on the wave equation in the processof 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 depthmigration method for accommodating absorption and dispersion effects ina prestack depth migration scheme. Their method was presented in thefrequency-wavenumber domain using a standard linear solid model. An ex-trapolation operator that compensates for absorption and dispersion wasdesigned using an optimization algorithm, which was then used to devise69Chapter 6. Refocusing Migrated Images in Absorptive Mediaa finite-difference migration scheme. Yu et al. (2002) applied poststackQ-migration to enhance the seismic frequencies on the top of an anticlinestructure 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 thatthe amplitude of the wavefield continuation operator used in Q-migrationis greater than 1, while the amplitude of the continuation operator used inregular migration is equal to or less than 1. In practice, a Q-migration al-gorithm must deal with the instability problem caused by the amplificationof high frequency random noise in its implementation. Conventionally, ifa Q-migration operation is based on wavefield extrapolation, a noise sup-pression step is required in the wavefield continuation process in order tostabilize the computation. Another option is to use least squares migration(LSM) (Nemeth et al., 1999; Zhang, 2000). A LSM scheme for wavefieldQ-compensation may be useful in stabilizing the solution, but to obtain asolution iteratively is computationally expensive.Traditional trace by trace inverse Q filtering neglects the 3D featuresof wavefield propagation, and regular Q-migration suffers from the amplifi-cation of high frequency random noise. This chapter aims to find anotherapproach to do wavefield Q-compensation, which can be considered as amulti-dimensional Q compensation technique. This multi-dimensional Q-compensation problem will be handled as a problem of multi-dimensionaldeconvolution 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, allof which lead to a blurred migration image. The blurring function is alsotermed as a point scatter function (PSF) (Jansson, 1997), an illuminationfunction (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 theearth?s structures.While many efforts are focused on trying to design accurate migrationoperators 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 acquisitionfootprint in prestack and poststack migrated images. The key in migrationdeconvolution techniques is how to calculate the blurring function. The next70Chapter 6. Refocusing Migrated Images in Absorptive Mediaimportant step is how to remove it. This chapter borrows the idea of migra-tion deconvolution to handle the blurring problem caused by using a regularmigration operator to migrate attenuated seismic data.In the following, we will discuss the attenuation effect on the migrationblurring function. Removing the blurring influence from seismic images isa problem of a 2D or 3D non-stationary deconvolution which can be solvedusing a least squares inverse scheme. The stability and efficiency of the leastsquares solver will be discussed.6.2 Blurring Function caused by SeismicAbsorptionAccording to Claerbout (1992), operations of seismic migration and mod-eling are conjugate pairs. Assuming data acquisition is a modeling processL, the migration process is equivalent to applying the conjugate LT to therecorded data Lm0, i.e.m = LT Lm0, (6.1)where m0 is the true image or the real earth structure, and m is the migratedimage, or the blurred earth structure, because LT is not the exact inverseof L.For seismic waves traveling in absorptive media, mathematically, a mi-gration blurring function can be defined asG(r|r0) = LT Lqdelta(r0), (6.2)where delta(r0) stands for a delta function at location r0 and Lq stands foran operator of wave propagation in an absorptive medium. Equation (6.2)implies that when regular migration is applied to the attenuated seismicdata, the blurring effect caused by absorption is mixed with other blurringeffects in regular migration as those mentioned before.A direct solution to obtain the actual image from equation (6.1) ism0 = bracketleftbigLT Lbracketrightbig-1 m.However, the inverse of LT L is generally difficult to calculate directly. Oneway to obtain an approximation of the inverse is to use a conjugate gradientalgorithm to solve a least squares problem (Nemeth et al., 1999; Sacchi et al.,2006). The same statement applies to the operator LT Lq in equation (6.2).71Chapter 6. Refocusing Migrated Images in Absorptive Media6.2.1 Analytical AnalysisAssuming that there is a monochromatic source of angular frequency omegastarting from a point scatterer located at r0 in a homogeneous medium withvelocity v and quality factor Q, the wavefield at a receiver location rg is thesolution of the Helmholtz Green?s function (Schuster and Hu, 2000),G(rg|r0,omega) = e-j omegac(omega) |rg-r0||rg -r0| .In this expression, omegac(omega) is a complex wavenumber. Absorption caused wave-field attenuation and dispersion are represented through the complex veloc-ity c(omega). Its expression will be given in equation (6.7).The complex conjugate of the Green?s function is the operator used inmigration which is equivalent to the wavefield of an impulse starting from areceiver rg down to an image point r,G(r|rg,omega) = Gasteriskmath(rg|r,omega) = ej omegav0 |rg-r||rg -r| .The migration blurring function in an absorptive medium, if the migrationis implemented without absorption compensation, isG(r|r0,omega) =integraldisplayrge-j omegac(omega) |rg-r0||rg -r0|ej omegav0 |rg-r||rg -r| drg. (6.3)The blurring function in the time domain is the integral of all frequencycomponentsG(r|r0,t) =integraldisplay bracketleftBiggintegraldisplayrge-j omegac(omega) |rg-r0||rg -r0|ej omegav0 |rg-r||rg -r| drgbracketrightBiggdomega. (6.4)The closed form of this integral is difficult to derive. However, Schuster andHu (2000) derived closed forms of both 2D and 3D blurring functions, foracoustic media. For an acquisition with a large aperture, it is assumed thatthe receivers are symmetrically distributed above a point scatterer. Themain energy of a migration blurring function is localized in the shadow areaas that shown in Figure 6.1 for homogeneous media.6.2.2 Numerical ComputationAlgorithm ChoiceIf 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 calculate72Chapter 6. Refocusing Migrated Images in Absorptive MediaFigure 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 becalculated numerically. The calculation is implemented as follows: first com-pute synthetic data with attenuation Lqdelta(r) and then migrate the syntheticto obtain the migration blurring function LT Lqdelta(r).The synthetic data can be computed using phase-shift modeling, or us-ing ray-tracing plus wavelet continuation along the ray-paths. Phase-shiftmodeling is based on the following wavefield extrapolation operatorH(omega,z) = e-jzradicalk2a-k2x, (6.5)and assumes an impulse begins at a point scatterer. In equation (6.5), ka isa complex wavenumber (refer to equations (2.13 and 2.15)) , andka = omegav0bracketleftbigg1 - 1piQ lnparenleftbigg omegaomega0parenrightbiggbracketrightbigg bracketleftbigg1 + j 1Qbracketrightbigg, (6.6)or ka = omegac(omega) , with1c(omega) =1v0bracketleftbigg1 - 1piQ lnparenleftbigg omegaomega0parenrightbiggbracketrightbigg bracketleftbigg1 + j 1Qbracketrightbigg. (6.7)Phase-shift modeling is easy to compute but the synthetic always containsFourier wrap-around noise. Ray-tracing synthetic computation is expen-sive due to the cost of ray-bending ray-tracing (Cerveny, 2005). Once a73Chapter 6. Refocusing Migrated Images in Absorptive Mediaray-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 purposeof Q-compensation is purely to enhance the frequency without consideringvelocity dispersion, a wavelet can be added to a record at an arrival timebased on the following expression (equation (3.3))B(omega,t2) = B(omega,t1)e-omega(t2-t1)4piQ ,where B(omega,t1) and B(omega,t2) are the wavelet spectra at two different timelocations 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 modelingand phase-shift migration (Gazdag, 1978) only apply to vertically variantmedia. 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 thephase-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 UpdatingWhen the phase-shift method is used to calculate the migration blurringfunction, there exists a fast algorithm for extrapolating a blurring functionfrom one depth step to another depth step or from one time step to anothertime step. For vertically layered media, vi is the reference velocity of layeri, and ci(omega) is complex velocity of layer i including the effect of seismicabsorption. If there is a unit impulse at a depth z0 = ndeltaz, and its responseon the surface isU(kx,z = 0,omega|z0) =nproductdisplayi=0ejdeltazromega2ci(omega)2 -k2x,and the response of a unit impulse right below it at depth z0 + deltaz isU(kx,z = 0,omega|z0 + deltaz) = ejdeltazromega2cn+1(omega)2 -k2xU(kx,z = 0,omega|z0).To migrate the impulse response U(kx,z = 0,omega|z0) from the surface to animage point z = mdeltaz using downward wavefield continuation, the blurringfunction isG(kx,z,omega|z0) = U(kx,z = 0,omega|z0)mproductdisplayi=0e-jdeltazromega2v2i -k2x. (6.8)74Chapter 6. Refocusing Migrated Images in Absorptive MediaTo migrate the impulse response U(kx,z = 0,omega|z0 + deltaz) from the surfaceto the same image point z = mdeltaz, the blurring function isG(kx,z,omega|z0 + deltaz) = U(kx,z = 0,omega|z0 + deltaz)mproductdisplayi=0e-jdeltazromega2v2i -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 depthwithout going through the full process of forward modeling and migration,i.e.G(kx,z,omega|z0 + deltaz) = ejdeltazromega2cn+1(omega)2 -k2xG(kx,z,omega|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 twodifferent time locations. The main energy of a blurring function lies in anarea, whose shape is like that sketched in Figure 6.1. As time increases, theenergy distribution of the migration blurring function extends to a largerarea.For layered media, if velocity v and attenuation property Q do not changelaterally, the migration blurring function needs only to be computed at onereference location. To allow for lateral v and Q variation, practically we candivide 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, assumingthat in each patch, the lateral variations of velocity and Q are very gentle.6.3 Removing Absorption Blurring Functionusing a Least Squares MethodAs discussed in the foregoing sections, the migration blurring function G =LT Lq is not an identity matrix. In our case, the migrated image is a blurredimage due to using an inaccurate migration operator,m = LT Lqm0. (6.11)In order to compensate for absorption after migration to get a solution closerto the right image m0, the blurring function needs to be removed or at leastto be shrunk.Generally, the inverse of LT Lq cannot be solved directly. There areexamples of using the diagonal to approximate the effects of the matrix onthe image amplitude (Sacchi and Wang, 2007; Hu et al., 2001). In fact,75Chapter 6. Refocusing Migrated Images in Absorptive Media00.51.0Time (s)1.0 1.2 1.4 1.6 1.8Midpoint (km)(a)00.51.0Time (s)1.0 1.2 1.4 1.6 1.8Midpoint (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.76Chapter 6. Refocusing Migrated Images in Absorptive Mediathe migrated image can be considered as the real image convolved with themigration blurring function. The blurring function is multi-dimensional (2Dor 3D) and non-stationary in time, i.e. equation (6.11) can be written in aform of convolutionm(vector,t) = gammavectorxref,tj(vector,t) asteriskmathmd m0(vector,t) + n(vector,t). (6.12)In equation (6.12), symbol asteriskmathmd represents md dimensional convolution, mdis either equal to 2 or 3. vector represents 1 or 2 lateral spatial coordinates,and t could be z for depth imaging. gammavectorxref,tj(vector,t) represents the migrationblurring function at a reference location (vectorref,tj). n(vector,t) is random noise.The inverse problem could be unstable because of two reasons: the inverseof gamma is unstable or even does not exist, and the inverse filtering could raisethe energy of random noise.Solving equation (6.12) for the real image m0(vector,t) is a problem of multi-dimensional non-stationary deconvolution. Time-variant deconvolution canbe formulated as a least squares inverse problem (Zhang and Ulrych, 2007).Similarly, the image refocusing process can be solved using the least squaresinversion in the wavenumber domain. To expedite the computations, weassume that gamma(r|r0) is shift invariant in the horizontal coordinates. AfterFourier transform on the spatial axes, equation (6.12) becomes? (vector,t) = ?vectorxref,tj(vector,t) asteriskmath ? 0(vector,t) + ?(vector,t). (6.13)Equation (6.13) changes the multi-dimensional convolution problem to a onedimensional time-variant convolution. Accordingly, the multi-dimensionaldeconvolution can be solved as a one-dimensional time-variant deconvolutionat each wavenumber.Similar to that discussed in Zhang and Ulrych (2007), the time-variantdeconvolution can be formulated as a least square problem using Bayes?theory. The objective function to be minimized is,phi( ?0) = c( ?0) + 12sigma2n(? -G? 0)T (? -G ?0). (6.14)In this equation, the first term c( ?0) on the right-hand side is a modelregularizer, and the second term is the misfit which could be interpreted asGaussian noise in the observed data. ?0 and ? represent model and datain the wavenumber domain respectively.The objective function can be minimized through solving the followingequation,?0 = parenleftbigGT G + lambdaW-1parenrightbig-1 GT ?, (6.15)77Chapter 6. Refocusing Migrated Images in Absorptive Mediawhere G is the kernel matrix of the inverse problem, W is a regularizer forthe image at a wavenumber and lambda is a hyper-parameter. The kernel matrixG is formed from a serials of time-variant complex ? functions, andG =bracketlefttpbracketleftexbracketleftexbracketleftexbracketleftexbracketleftexbracketleftexbracketleftexbracketleftexbracketleftbt?0(t -0deltat)?1(t -1deltat)...?i(t -ideltat)...?N(t -Ndeltat)bracketrighttpbracketrightexbracketrightexbracketrightexbracketrightexbracketrightexbracketrightexbracketrightexbracketrightexbracketrightbt. (6.16)In equation (6.16), N is the number of trace samples. ?i(t-ideltat) representsa migration blurring function in the wavenumber domain at time ideltat witha time shift of ideltat. The weights at wavenumber kj are determined by theweights and the obtained image at lower wavenumber kj-1Wii[kj] = alphas(mi[kj-1]) + (1 -alpha)Wii[kj-1], (6.17)where alpha element [0,1] is a ?viscosity? parameter (Hugonnet et al., 2001) whichmeans the similarity of energy distribution from low wavenumber to highwavenumber. s(mi) is a function of model parameters,s(mi) = (|mi| + epsilon1)/(maxi |mi| + epsilon1).This model constraint implies that the inversion favors strong energy. Ifthe weights are equal, the solution is standard least squares (Robinson andTreitel, 2002). Although the theoretical soundness of this constraint stillneeds to be further investigated, it works for the two data examples whichwill be discussed in the following sections.A pseudo code in C language of the refocusing algorithm for 2D migratedtime section is listed in table (6.3). The main computation of the refocusingprocess has two parts: one is the calculation of the blurring functions, an-other is to solve for m0 using equation (6.15). To extend the algorithm from2D to 3D is very straightforward. The code needs to be modified to handleone extra spatial dimension y. For phase-shift migration, the blurring func-tion is easy to compute using equation (6.10). For Kirchhoff migration, thecomputation of the migration blurring function relies on the speed of raytracing. The algorithm is implemented wavenumber by wavenumber whichis easy to be parallelized for distributed computations. In each wavenum-ber the main computation is matrix multiplication of size (N ?N) which issubstantial but not excessive. Therefore, the implementation is practical.78Chapter 6. Refocusing Migrated Images in Absorptive Media/* part 1 calculate the blurring function */{do while (time<maximum_time);{calculate_synthetic ();do_migration ();do_fft_x2kx (a_blurring_function);keep_blurring_function ();time += time_step;}}/* part 2 deconvolution using least squares inversion */{do_fft_x2kx (a_section);for (ikx=0, ikx<maximum_kx; ikx++){form_kernel_matrix ();solve_for_m0 ();}do_fft_kx2x (kx_section);}Table 6.1: A pseudo code of the refocusing algorithm.79Chapter 6. Refocusing Migrated Images in Absorptive Media0200400Time Index 0 200 400Time IndexFigure 6.3: Amplitude of a kernel matrix G.The algorithm is stable because GT G is a diagonally dominant Hermi-tian matrix. The diagonal elements are close to its eigenvalues because theenergy of a blurring function is mainly concentrated in a time window cen-tered at the point scatterer. The magnitudes of the elements decay alongthe diagonal due to amplitude attenuation, but they are greater than zero,which means the kernel matrix is not singular. Figure 6.3 shows an exam-ple of the element amplitude of this matrix. The ripples at the top rightcorner and the bottom left corner are Fourier wrap-around noise in phaseshift modeling and migration. They need to be muted out in real data pro-cessing, because they do not exist in observed real data. G can also beconsidered as a band-shaped matrix because of the localization of blurringfunctions. By considering this characteristic, the computation required bymatrix manipulation could be greatly reduced. There are also efforts madeto approximate G by a diagonal matrix in order to save computation withthe trade-off with accuracy (Sacchi and Wang, 2007). Here, the solution isobtained by solving equation (6.15) using the full matrix of G.Theoretically, the refocusing technique used here is only valid for hori-zontally layered media. Consequently, it is better to be applied to refocusinga migrated time image without strong lateral v and/or Q variation. Other-wise, the migrated image needs to be processed in patches based on v and/or80Chapter 6. Refocusing Migrated Images in Absorptive MediaQ variation. These patches are, regularly, overlapping rectangles. In eachrectangular area, the refocusing operation is applied independently. Imagesin overlapping areas are the mathematical average of those involved patches.The final result will depend on the accuracy of the computed blurring func-tions, the inversion algorithm, the regularization parameter and also patchesbeing divided.6.4 Data ExperimentsThe refocusing algorithm is first tested on a synthetic model where its mi-gration blurring functions are known. Figure 6.4(a) shows a sine-functionshaped structure obtained in an absorptive medium. Figure 6.4(b) showsthe refocusing result. Because the refocusing process compensates for theattenuation effect, the resolution of the obtained image is enhanced.The algorithm is also tested on a 2-D field data-set. Its migration blur-ring functions are computed based on the velocity and Q fields. The velocityfield is shown in Figure 6.5(a) and its inverse Q field is in Figure 6.5(b). Theinverse Q field is obtained using the procedure described in table 3.3.1, whichshows that the major absorption is before 1.4 s with less attenuation after-wards. This absorption profile is generally in line with geologic structuresat different environments, because rocks at deep places are relatively harderand better compacted. The refocusing process does not require the Q in-formation, but Q values are required when doing the forward modeling tocalculate the absorption blurring functions. Both velocity and Q fields arelaterally smoothed and the values at the center CMP (400) location are usedto calculate the blurring functions at each time sample. Figure 6.6(a) showsa local area from the image obtained by time-variant spectral whitening(TVSW) plus phase-shift migration. Figure 6.6(b) shows the correspondingarea from the image obtained using phase-shift migration plus a migrationrefocusing process. The section obtained by the refocusing process showsthat the algorithm shrinks migration noise and increases the seismic resolu-tion. The overall quality of images has been enhanced. For example, thosesmall vertical faults at time 2.25 s are much clearer in Figure 6.6(b) thanthat in Figure 6.6(a). The comparison shows the difference between 1Dabsorption compensation and 2D absorption compensation. The resolutionimprovement of Figure 6.6(b) compared to Figure 6.6(a) shows the validityof the refocusing algorithm in compensating for the absorption after migra-tion. If the image were processed in patches, the quality of the refocusedimage would be, hopefully, further improved.81Chapter 6. Refocusing Migrated Images in Absorptive Media00.51.01.52.0Time (s)0 0.2 0.4 0.6Midpoint (km)(a)00.51.01.52.0Time (s)0 0.2 0.4 0.6Midpoint (km)(b)Figure 6.4: A synthetic example. (a) shows an absorption-blurred sine-function shaped structure. (b) is the refocusing result.82Chapter 6. Refocusing Migrated Images in Absorptive Media(a)(b)Figure 6.5: Input information. (a) is the velocity section and (b) is theinverse Q section.83Chapter 6. Refocusing Migrated Images in Absorptive Media1.52.02.5Time (s)300 350 400 450 500 550CMP(a)1.52.02.5Time (s)300 350 400 450 500 550CMP(b)Figure 6.6: Migration results. (a) Time-variant spectral whitening plusphase-shift migration. (b) Phase-shift migration plus refocusing.84Chapter 6. Refocusing Migrated Images in Absorptive Media6.5 Conclusion and DiscussionTo obtain the real seismic image, the Hessian matrix (LT L) needs to beclose to an identity matrix, or the migration blurring function should be aspatial delta function in 2D or 3D domain. Schuster and Hu (2000), Hu et al.(2001) and Yu et al. (2006) investigated the acquisition footprint, or non-uniform illumination resulting from limited recording aperture by meansof migration Green?s functions. Here, we have put aside the acquisitiongeometry blurring function and only discuss the blurring effect caused by aninaccurate migration operator when migrating data without considering theabsorption. The numerical implementation of refocusing a migrated imagein an absorptive medium includes the following three steps:1. Make synthetic data with attenuation.2. Migrate the synthetic data without attenuation to get blurring func-tions.3. Remove absorption blurring functions from the image using least squaresinversion.Compared to the Q-migration algorithm (Mittet et al., 1995; Yu et al.,2002), the refocusing process is able to address the problem of instability onhigh frequency noise. Compared to the least squares migration algorithm(Nemeth et al., 1999; Sacchi et al., 2006), the refocusing algorithm is costeffective in implementation.When computation is not a big issue, the refocusing process can be ap-plied to process common offset sections. In doing so, this method can beextended to prestack seismic processing. If the blurring function includesthe effect of a seismic wavelet, the refocusing process is equivalent to post-stack multi-dimensional deconvolution. In practice, the migration blurringfunction is related to many factors, such as acquisition geometry, seismicwavelet, velocity structure and absorption properties. This chapter blendsthe absorption effect with other blurring effects into a comprehensive mi-gration blurring function, and multichannel Q-compensation is achieved inthe process to refocus the migrated image by time-variant deconvolution inthe wavenumber domain.85Chapter 7Discussion and ConclusionsThe results of this research provide new techniques for handling signal ab-sorption in seismic data imaging and inversion. The following paragraphsdelineate the specific contributions contained in this thesis.The introduced analytical formula (equation(3.13)) relates Q-factor tospectral peak frequency variation. The assumption in the theoretical deriva-tion is that the source wavelet has an amplitude spectrum similar to that ofa Ricker wavelet (equation(3.1)). Although an actual amplitude spectrummay not, of course, exactly conform to this assumption, in many cases, theamplitude spectrum can still be well approximated by a Ricker spectrum.This implies that this assumption is not rigorously required. This approachcan be used to estimate Q from CMP gathers, from a post-stack trace, andfrom seismic data acquired from other geometries like VSP (vertical seismicprofiling) as well, as long as the peak frequency variation can be reliablyexamined.To estimate Q-factor from a prestack CMP gather, variation of a waveletspectrum of an event is analyzed along offset. To obtain spectral peak fre-quency variation from a poststack trace, we can use windowed time-variantspectral analysis or a continuous wavelet transform analysis. In a WTVS,controlling points can be picked in a way similar to velocity analysis (Fig-ure 3.4). These controlling points divide a trace into layers and their coor-dinates are used to calculate the quality factors of each layer.In two respects, the analytical method is superior to the spectral ratiomethod. First, the analytical method only uses the information of frequencyvariation. Besides absorption, it is not affected by other amplitude factorssuch as reflection, transmission and geometrical spreading. Second, the ana-lytical method is robust to random noise. For the two-dimensional case, therobustness is achieved by means of using multi-fold information at differentoffsets. For the one-dimensional case, the robustness is achieved by meansof using the frequency trend analysis.Event-based Q estimation can provide the absorption information for Qcompensation. For reservoir description, seismic attenuation property, Q-factor, is ideally estimated at all prospective layers defined by velocities or86Chapter 7. Discussion and Conclusionsacoustic impedance. Techniques to analyze seismic attenuation propertiesfor hydrocarbon reservoir description generally include iso-frequency analy-sis, direct Q estimation and attenuation or Q tomography.The introduced reflectivity-guided seismic attenuation analysis is poten-tially another useful technique. It operates on prestack CMP gathers, anduses the location information of structure interfaces to select layer bound-aries. The interface information could be obtained from poststack acous-tic impedance inversion or well logs. This technique basically provides amethodology to allow Q values to be estimated at all prospective layers.The obtained Q section can have the same resolution as other attributesections. To estimate spectral peak frequency variation at zero offset, anAVO-like, peak frequency estimation scheme called PFVO which operateson CMP gathers is introduced. The advantage of PFVO over post-stack fre-quency analysis is that it avoids the distortion caused by seismic stacking onfrequency characteristics. PFVO tries to estimate a reliable peak frequencyat zero offset, and at the same time to identify the frequency change causedby real attenuation and that caused by event tuning.Absorption compensation or de-absorption is an important step in seis-mic data processing; it helps to improve the resolution of seismic data, andto balance the frequency contents. In order to stabilize traditional inverseQ filtering, the inverse Q filtering is formulated as a least squares inverseproblem based on Bayes? theory. The prior information imposed on themodel parameters is the sparseness of reflectivity, which is represented by aCauchy model, and random noise in observed data is assumed to be a Gaus-sian. The inverse problem is solved iteratively. Experiments performed onboth synthetic and real data show that the inverse scheme can produce idealresults for absorption compensation.To preserve both the amplitude and waveform of seismic signals dur-ing imaging, it is necessary to incorporate de-absorption and migration to-gether. This may be implemented through standard migration with randomnoise attenuation at each continuation step, or through least squares migra-tion which is implemented by casting the standard migration problem asan optimization problem and solving the optimization problem iteratively.Based on the concept of migration deconvolution, a post-migration refocus-ing scheme is introduced which could be considered as an alternative toconventional migration with Q compensation. The refocusing scheme canbe used to remove the absorption effect on seismic imaging, at the sametime to remove the other migration artifacts.Compared to regular migration with Q compensation, the refocusingmethod handles the instability problem caused by the amplification of the87Chapter 7. Discussion and Conclusionshigh frequency noise by means of inversion. Compared to the least squaresmigration, the refocusing method is cheaper in implementation. The refo-cusing method can be applied to common offset sections by using differentoperators at different offsets. In doing so, this method can be extendedto prestack seismic processing. If the blurring function includes the effectof seismic wavelet, the refocusing process is equivalent to multidimensionaldeconvolution.The work presented in this thesis cover just two areas of seismic absorp-tion research: seismic absorption estimation and compensation. Q-factorestimation from seismic data is the first step to use attenuation proper-ties for seismic data imaging and inversion. To confirm the accuracy ofthe estimated Q values, numerical modeling, such as reflectivity modeling(Kennett, 1983) or waveform modeling (Carcione, 1995), may be required tocheck whether synthetic records match observations. For seismic modeling,a mathematical model is expected to describe wave behavior in absorptivemedia. Generally, a Q model (see equations (2.13) and (2.15)) is enough todescribe the absorption phenomenon in visco-elastic media.To use an estimated Q section as an attribute to invert for reservoir prop-erties not only requires accurate Q information, but also depends on reliablerock physics models like those suggested by Biot (1956a), Biot (1956b), Prideet al. (2004) and Dvorkin and Mavko (2006). A rock physics model to relateQ with permeability, porosity, saturation and fluid properties etc. is ide-ally expected. Rock physical properties cannot be inverted using Q-factorinformation only. However, integrated inversion of using absorption prop-erty, AVO, velocity and well log information is possibly one practical wayfor obtaining useful parameters for hydrocarbon reservoir description.All aforementioned seismic information could be of P-waves, S-waves orboth. Recent developments in multi-component seismic wave acquisitionmake it possible to obtain S-wave velocities, S-wave AVO effects, and S-wave absorption properties as well. There are reports showing applicationsof using the relation between P-wave and S-wave Q factors as an indicatorto identify brine saturated and gas filled layers in seismic sections (Mavkoet al., 2005a; Bale and Stewart, 2002). Integrated inversion of using bothP-wave and S-wave absorption properties, velocities, AVO effects and othergeophysical information for rock properties could be an interesting topic forfuture research.Accurate Q estimation from CMP gathers provides critical informationfor prestack Q compensation. With quantitative attenuation information,it may be possible to separate AVO effects from seismic absorption in CMPsuper-gathers. That is to say, AVO information extracted from seismic data88Chapter 7. Discussion and Conclusionswould be more reliable if absorption effects on AVO had been taken intoaccount in the analysis accurately. The attenuation effect on AVO needs tobe further investigated.Absorption compensation can be performed using the introduced refo-cusing method after post-stack migration. It can also be implemented onprestack common offset sections. More research is expected to investigatethe difference between poststack and prestack refocusing operations whenunderground structures have strong lateral velocity and Q variations.89BibliographyAki, K. and P. G. Richards, 1980, Quantitative seismology-theory andmethods: W. H. Freeman and Company.Bachrach, R., R. Ferber, A. Salama, N. Banik, R. Utech, D. Keller, and M.Bengtson, 2006, Two approaches for Q estimation and its application toseismic inversion: 76th SEG, expanded abstract, 2971?2974, Soc. of Expl.Geophys.Bale, R. and R. Stewart, 2002, The impact of attenuation on the resolutionof multi-component seismic data: 72nd SEG, expanded abstract, 1034?1037, Soc. of Expl. Geophys.Batzle, M., R. Hofmann, M. Prasad, G. Kumar, L. Duranti, and D. Han,2005, Seismic attenuation: observations and mechanisms: 75th SEG, ex-panded abstract, 1565?1568, Soc. of Expl. Geophys.Berkhout, A. J., 1981, Wave field extrapolation techniques in seismic mi-gration, a tutorial: Geophysics, 46, 1638?1656.Bickel, S. H., 1993, Similarity and the inverse Q filter: The Pareto-Levystretch: Geophysics, 58, 1629?1633.Biot, M. A., 1956a, Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low-frequency range: J. Acous. Soc. Am., 28,168?178.???, 1956b, Theory of propagation of elastic waves in a fluid-saturatedporous solid. II. high-frequency range: J. Acous. Soc. Am., 28, 179?191.Blanch, J., J. Robertsson, and W. Symes, 1995, Modeling of a constantQ: Methodology and algorithm for an efficient and optimally inexpensiveviscoelastic technique: Geophysics, 60, 176?184.Bohlen, T., 2002, Parallel 3-D viscoelastic finite-difference seismic model-ing: Computers and Geosciences, 28, 887?899.Brillouin, L., 1960, Wave propagation and group velocity: Academic Press.90BibliographyBuckingham, M. J., 2005, Causality, Stokes?s wave equation, and acousticpulse propagation in a viscous media: Physical Review, E72, 026610?1?026610?8.Cand`s, E., J. Romberg, and T. Tao, 2005, Stable signal recovery fromincomplete and inaccurate measurements: Comm. Pure Appl. Math., 59,1207?1223.Carcione, J., D. Kosloff, and R. Kosloff, 1988, Viscoacoustic wave propa-gation simulation in the earth: Geophysics, 53, 769?777.Carcione, J. M., 1993, Seismic modeling in viscoelastic media: Geophysics,58, 110?120.???, 1995, Constitutive model and wave equations for linear, viscoelastic,anisotropic media: Geophysics, 60, 537?548.Carcione, J. M. and S. Picotti, 2006, P-wave seismic attenuation by slow-wave diffusion: Effets of inhomogeneous rock properties: Geophysics, 71,O1?O8.Castagna, J. P. and S. Sun, 2006, Comparison of spectral decompositionmethods: First Break, 24, 75?79.Castagna, J. P., S. Sun, and R. W. Siegfried, 2003, Instantaneous spectralanalysis: Detection of low-frequency shadows associated with hydrocar-bons: The Leading Edge, 22, 120?127.Cerveny, V., 2005, Seismic ray theory: Cambridge Univ. Press.Chapman, M., E. Liu, and X. Li, 2005, The influence of abnormally highreservoir attenuation on the AVO signature: The Leading Edge, 11, 1120?1125.Claerbout, J., 1985, Imaging the earth?s interior: Blackwell Scientific Pub-lications, Boston.???, 1992, Earth soundings analysis: Processing versus inversion: Black-well Scientific Publications, Boston.Cui, J. and J. He, 2004, 2D seismic migration with compensation: a pre-liminary study: Journal of Geophysical and Engineering, 1, 263?267.Dai, N. and G. F. West, 1994, Inverse Q migration, in 64th SEG, ExpandedAbstracts, 1418?1421, Soc. of Expl. Geophys.Dasgupta, R. and R. A. Clark, 1998, Estimation of Q from surface seismicreflection data: Geophysics, 63, 2120?2128.91BibliographyDay, S. M. and J. B. Minster, 1984, Numerical simulation of attenuatedwavefields using a Pade approximate method: Geophys. J. R. Astr. Soc.,78, 105?118.Deffenbaugh, M., A. Shatilo, B. Schneider, and M. Zhang, 2000, Resolutionof converted waves in attenuation media: 70th SEG, expanded abstract,1189?1192, Soc. of Expl. Geophys.Du, S., 1996, Seismic dynamics (in Chinese): University of Petroleum Pub-lishing Co., Dongying, China.Duren, R. E. and R. L. Heestand, 1995, General solutions and causality fora Voigt medium: geophysics, 60, 252?261.Dvorkin, J. and G. Mavko, 2006, Modeling attenuation in reservoir andnon-reservoir rock: The Leading Edge, 25, 194?197.Dvorkin, J., G. Mavko, and A. Nur, 1995, Squirt flow in fully saturatedrocks: Geophysics, 60, 97?107.Ebrom, D., 2004, The low-frequency gas shadow on seismic sections: Theleading edge, 23, 772.Futterman, W. I., 1962, Dispersive body waves: J. Geophys. Res., 67,5279?5291.Gabor, D., 1946, Theory of communication: Journal of IEEE, 93, 429?457.Gao, F., G. Fradelizio, A. Levander, G. Pratt, and C. Zelt, 2005, Seismicvelocity, Q, geological structure and lithology estimation at a ground watercontamination site: 65th SEG, Expanded Abstracts, 1561?1564, Soc. ofExpl. Geophys.Gazdag, J., 1978, Wave equation migration with the phase-shift method:Geophysics, 43, 1342?1351.Guerra, R. and S. Leaney, 2006, Q(z) model building using walkaway VSPdata: Geophysics, 71, V127?v132.Guitton, A., 2004, Amplitude and kinematic corrections of migrated imagesfor non-unitary imaging operators: Geophysics, 69, 1017?1024.Hale, D., 1991, Stable explicit depth extrapolation of seismic wavefields:Geophysics, 56, 1770?1777.Hargreaves, N. D. and A. J. Calvert, 1991, Inverse Q-filtering by Fouriertransform: Geophysics, 56, 519?527.Hermann, F. J., 2005, Robust curvelet-domain data continuation withsparseness constraints: 67th EAGE, Expanded Abstracts, 4248?4251, Eur.Assn. of Geoscientists and Engineers.92BibliographyHicks, G. J. and R. G. Pratt, 2001, Reflection waveform inversion usinglocal descent methods: Estimating attenuation and velocity over a gas-sanddeposit: Geophysics, 66, 598?612.Hu, J., G. T. Schuster, and P. Valasek, 2001, Poststack migration decon-volution: Geophysics, 66, 939?952.H?bert, L., U. Strecker, J. Dvorkin, and K. A. Festervoll, 2005, Seismicattenuation and hybrid attributes to reduce exploration risk-North Seacase study: 75th SEG, Expanded Abstracts, 436?439.Hugonnet, P., P. Herrmann, and C. Tibeiro, 2001, High resolution Radon:a review: 63th EAGE, Expanded Abstracts, 3044?3047, Eur. Assn. of Geo-scientists and Engineers.Jansson, P., 1997, Deconvolution of images and spetra, 2nd ed.: AcademicPress.Jaynes, E. T., 2004, Probability theory: The logic of science: CambridgeUniversity Press. http://bayes.wustl.edu/etj/prob.html.Johnson, D. L., 2001, Theory of frequency dependent acoustics in patchy-saturated porous media: J. Acoust. Soc. Am., 110, 682?694.Kennett, B., 1983, Seismic wave propagation in stratified media: Cam-bridge Univ. Press.Kjartansson, E., 1979, Constant Q wave propagation and attenuation: J.Geophys. Res., 84, 4737?4748.Klimentos, T., 1995, Attenuation of P- and S-waves as a method of distin-guishing gas and condensate from oil and water: Geophysics, 60, 555?556.Kneib, G. and S. Shapiro, 1995, Viscoacoustic wave propagation in 2-Drandom media and separation of absorption and scattering attenuation:Geophysics, 60, 459?467.Lancaster, S. J. and M. C. Tanis, 2004, High density 3D prestack Q estima-tion: 66th EAGE, Expanded Abstracts, P0004, Eur. Assn. of Geoscientistsand Engineers.Latimer, R. B., R. Davidson, and P. van Riel, 2000, An interpreter?s guideto understanding and working with seismic-derived acoustic impedancedata: The Leading Edge, 19, 242?256.Liu, H., D. L. Anderson, and H. Kanamori, 1976, Velocity dispersion due toanelasticity; implications for seismology and mantle composition: Geophys.J. R. astr. Soc., 47, 42?58.93BibliographyMallat, S., 1999, A wavelet tour of signal processing: Academic Press.Malvern, L., 1969, Introduction to the mechanics of a continuous media:Prentice-Hall Ltd. Toronto, New Jersey.Margrave, G. F., D. C. Henley, M. P. Lamoureux, V. Iliescu, and J. P.Grossman, 2003, Gabor deconvolution revisited: 73rd SEG, Expanded Ab-stracts, 714?717, Soc. of Expl. Geophys.Masoomzadeh, H., P. Barton, and S. C. Singh, 2004, Non-stretch stackingin the tau-p domain: Exploiting long-offset arrivals for sub-basalt imaging:74th SEG, Expanded Abstracts, 2132?2155, Soc. of Expl. Geophys.Mavko, G., J. Dvokin, and J. Wales, 2005a, A rock physics and attenuationanalysis of a well from the Gulf of Mexico: 75th SEG, expanded abstract,1585?1588, Soc. of Expl. Geophys.???, 2005b, A theoretical estimate of S-wave attenuation in sediment:75th SEG, expanded abstract, 1469?1472, Soc. of Expl. Geophys.Mavko, G. and A. Nur, 1979, Wave attenuation in partially saturated rocks:Geophysics, 44, 161?178.McGinn, A. and B. Duijndam, 1998, Land seismic data quality improve-ments: The Leading Edge, 17, 1570?1577.Mittet, R., R. Sollie, and K. Hokstad, 1995, Prestack depth migration withcompensation for absorption and dispersion: Geophysics, 60, 1485?1494.Nemeth, T., C. Wu, and G. T. Schuster, 1999, Least-squares migration ofincomplete reflection data: Geophysics, 64, 208?221.O?Doherty, R. F. and N. A. Anstey, 1971, Reflections on amplitudes: Geo-physical Prospecting, 19, 430?458.Oldenburg, D. W., T. Scheuer, and S. Levy, 1983, Recovery of the acousticimpedance from reflection seismograms: Geophysics, 48, 1318?1337.Ostrander, W. J., 1984, Plane-wave reflection coefficients for gas sands atnonnormal angles-of-incidence: Geophysics, 49, 1637?1648.Parra, J. O. and C. Hacket, 2002, Wave attenuation attributes as flow unitindicators: The Leading Edge, 21, 564?572.Parra, J. O., C. Hacket, P.-C. Xu, and H. Collier, 2006, Attenuation analy-sis of acoustic waveforms in a borehole intercepted by a sand-shale sequencereservoir: The Leading Edge, 21, 186?193.Perroud, H. and M. Tygel, 2004, Nostretch NMO: Geophysics, 69, 599?607.94BibliographyPlessix, R.-E., 2006, Estimation of velocity and attenuation coefficientmaps from crosswell seismic data: Geophysics, 71, S235?S240.Pratt, R. G., K. Bauer, and M. Weber, 2003, Crosshole waveform tomog-raphy velocity and attenuation images of arctic gas hydrates: 73rd SEG,Expanded Abstracts, 2255?2258, Soc. of Expl. Geophys.Pride, S. R., 2003, Permeability dependence of seismic amplitudes: TheLeading Edge, 22, 518?525.Pride, S. R., J. G. Berryman, and J. M. Harris, 2004, Seismic attenuationdue to wave induced flow: Journal of Geophysical Research, 109, B01201.Quan, Y. and J. M. Harris, 1997, Seismic attenuation tomography usingthe frequency shift method: Geophysics, 62, 895?905.Richards, P. G. and W. Menke, 1983, The apparent attenuation of a scat-tering medium: Bulletin of the Seismological Society of America, 73, 1005?1021.Ricker, N., 1953, The form and laws of propagation of seismic wavelets:Geophysics, 18, 10?40.???, 1977, Transient waves in visco-elastic media: Elsevier Scientific Pub.Co.Rickett, J., 2003, Illumination-based normalization for wave-equationdepth migration: Geophysics, 68, 1371?1379.???, 2006, Integrated estimation of interval-attenuation profiles: Geo-physics, 71, A19?A23.Robinson, E. A. and S. Treitel, 2002, Geophysical signal analysis: Societyof Exploration Geophysicists Publications.Robinson, J., 1979, A technique for the continuous representation of dis-persion in seismic data: Geophysics, 44, 1345?1351.Sacchi, M., T. J. Ulrych, and C. Walker, 1998, Interpolation and extrapola-tion using a high resolution discrete Fourier transform: IEEE Transactionson Signal Processing, 46, 31?38.Sacchi, M. and J. Wang, 2007, Estimation of the diagonal of the migrationblurring kernel through a stochastic approximation: 77th SEG, ExpandedAbstracts, 2438?2441, Soc. of Expl. Geophys.Sacchi, M., J. Wang, and H. Kuehl, 2006, Regularized migration/inversion:New generation of seismic imaging algorithms: Recorder, CSEG, 55?59.95BibliographySams, M. S., J. P. Neep, M. H. Worthington, and M. S. King, 1997, Themeasurement of velocity dispersion and frequency-dependent intrinsic at-tenuation in sedimentary rocks: Geophysics, 62, 1456?1464.Sato, H. and M. C. Fehler, 1998, Seismic wave propagation and scatteringin the heterogeneous earth: Springer-Verlag New York , Inc.Schuster, G. T. and J. Hu, 2000, Migration green?s function: Continuousrecording geometry: Geophysics, 65, 167?175.Singleton, S., M. T. Taner, and S. Treitel, 2006, Q estimation using Gabor-Morlet joint time-frequency analysis techniques: 76th SEG, expanded ab-stracts, 1610?1614.Spencer, T. W., J. R. Sonnad, and T. M. Butler, 1982, Seismic Q-stratigraphy or dissipation: Geophysics, 47, 16?24.Taner, M. T., A. F. Koehler, and R. E. Sheriff, 1979, Complex seismic traceanalysis: Geophysics, 44, 1041?1063.Taner, M. T. and S. Treitel, 2003, A robust method for Q estimation: 73rdSEG, Expanded Abstracts, 710?713, Soc. of Expl. Geophys.Tarantola, A., 1987, Inverse problem theory. methods for data fitting andmodel parameter estimation: Elsevier Science Publ. Co.Toksoz, M. N., D. H. Johnston, and A. T. Timur, 1979a, Attenuationseismic waves in dry and saturated rocks-I laboratory measurements: Geo-physics, 44, 681?690.???, 1979b, Attenuation seismic waves in dry and saturated rocks-IImechanisms: Geophysics, 44, 691?711.Tonn, R., 1991, The determination of the seismic quality factor Q fromVSP data: A comparison of different computational method: GeophysicalProspecting, 39, 1?27.Toverud, T. and B. Ursin, 2005, Comparison of seismic attenuation modelsusing zero-offset vertical seismic profiling VSP data: Geophysics, 70, F17?F25.Toxopeus, G., S. Petersen, J. Thorbecke, and K. Wappenaar, 2004, Simu-lating high resolution inversion: decomposing the resolution function: 74thSEG, Expanded Abstracts, 3044?3047, Soc. of Expl. Geophys.Ulrych, T. J. and M. D. Sacchi, 2005, Information-based inversion andprocessing with applications: Elsevier Publishing Company.Varela, C. L., A. L. Rossa, and T. Ulrych, 1993, Modeling of attenuationand dispersion: Geophysics, 58, 1167?11733.96BibliographyVersteeg, R., 1994, The marmousi experience: Velocity model determina-tion on a synthetic complex dataset: The Leading Edge, 13, 927?936.Wang, Y., 2002, A stable and efficient approach to inverse Q filtering:Geophysics, 67, 657?663.???, 2003, Quantifying the effectiveness of stabilized inverse Q filtering:Geophysics, 68, 337?345.???, 2006, Inverse Q-filter for seismic resolution enhancement: Geo-physics, 71, V51?V60.???, 2008, Inverse Q filtered migration: Geophysics, 74, P1?P8.Wang, Y. and J. Guo, 2004, Modified Kolsky model for seismic attenuationand dispersion: Journal of Geophysics and Engineering, 1, 187?196.Watanabe, T., K. T. Nihei, S. Nakagawa, and L. R. Myer, 2004, Viscoacous-tic waveform inversion of transmission data for velocity and attenuation:J. Acoust. Soc. Am., 115, 3059?3067.White, J. E., 1975, Computed seismic speeds and attenuation in rocks withpartial gas saturation: Geophysics, 40, 224?232.???, 1983, Underground sound application of seismic waves: ElsevierScience Publishing Company Inc.Winkler, K., A. Nur, and M. Gladwin, 1979, Friction and seismic attenua-tion in rocks: Nature, 277:528.Xie, X., S. Jin, and R. Wu, 2006, Wave-equation-based illumination anal-ysis: Geophysics, 71, S169?S177.Yilmaz, O., 2000, Seismic data analysis: Society of Exploration Geophysi-cist Publications.Yu, J., J. Hu, G. T. Schuster, and R. Estill, 2006, Prestack migrationdeconvolution: Geophysics, 71, S53?S62.Yu, Y., R. Lu, and M. Deal, 2002, Compensation for the effects of shallowgas attenuation with viscoacoustic wave-equation migration: 72nd SEG,Expanded Abstracts, 2062?2065, Soc. of Expl. Geophys.Zhang, C., 2000, Absorption compensation in seismic data processing: Uni-versity of British Columbia, Master thesis.Zhang, C. and T. J. Ulrych, 2002, Estimation of quality factors from CMPrecords: Geophysics, 67, 1542?1547.???, 2007, Seismic absorption compensation: a least squares inversescheme: Geophysics, 72, R109?R114.97
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Seismic absorption estimation and compensation
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Seismic absorption estimation and compensation Zhang, Changjun 2008-12-01
pdf
Page Metadata
Item Metadata
Title | Seismic absorption estimation and compensation |
Creator |
Zhang, Changjun |
Publisher | University of British Columbia |
Date Issued | 2008 |
Description | As seismic waves travel through the earth, the visco-elasticity of the earth's medium will cause energy dissipation and waveform distortion. This phenomenon is referred to as seismic absorption or attenuation. The absorptive property of a medium can be described by a quality factor Q, which determines the energy decay and a velocity dispersion relationship. Four new ideas have been developed in this thesis to deal with the estimation and application of seismic absorption. By assuming that the amplitude spectrum of a seismic wavelet may be modeled by that of a Ricker wavelet, an analytical relation has been derived to estimate a quality factor from the seismic data peak frequency variation with time. This relation plays a central role in quality factor estimation problems. To estimate interval Q for reservoir description, a method called reflectivity guided seismic attenuation analysis is proposed. This method first estimates peak frequencies at a common midpoint location, then correlates the peak frequency with sparsely-distributed reflectivities, and finally calculates Q values from the peak frequencies at the reflectivity locations. The peak frequency is estimated from the prestack CMP gather using peak frequency variation with offset analysis which is similar to amplitude variation with offset analysis in implementation. The estimated Q section has the same layer boundaries of the acoustic impedance or other layer properties. Therefore, the seismic attenuation property obtained with the guide of reflectivity is easy to interpret for the purpose of reservoir description. To overcome the instability problem of conventional inverse Q filtering, Q compensation is formulated as a least-squares (LS) inverse problem based on statistical theory. The matrix of forward modeling is composed of time-variant wavelets. The LS de-absorption is solved by an iterative non-parametric approach. To compensate for absorption in migrated seismic sections, a refocusing technique is developed using non-stationary multi-dimensional deconvolution. A numerical method is introduced to calculate the blurring function in layered media, and a least squares inverse scheme is used to remove the blurring effect in order to refocus the migrated image. This refocusing process can be used as an alternative to regular migration with absorption compensation. |
Extent | 1419131 bytes |
Subject |
Seismic absorption compensation Quality factor Q estimation Reservoir description Migrated image |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2008-12-01 |
Provider | Vancouver : University of British Columbia Library |
DOI | 10.14288/1.0052450 |
URI | http://hdl.handle.net/2429/2820 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2009-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2009_spring_zhang_changjun.pdf [ 1.35MB ]
- Metadata
- JSON: 24-1.0052450.json
- JSON-LD: 24-1.0052450-ld.json
- RDF/XML (Pretty): 24-1.0052450-rdf.xml
- RDF/JSON: 24-1.0052450-rdf.json
- Turtle: 24-1.0052450-turtle.txt
- N-Triples: 24-1.0052450-rdf-ntriples.txt
- Original Record: 24-1.0052450-source.json
- Full Text
- 24-1.0052450-fulltext.txt
- Citation
- 24-1.0052450.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0052450/manifest