UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Absorption compensation in seismic data processing Zhang, Changjun 2001

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2001-0141.pdf [ 7.52MB ]
Metadata
JSON: 831-1.0089920.json
JSON-LD: 831-1.0089920-ld.json
RDF/XML (Pretty): 831-1.0089920-rdf.xml
RDF/JSON: 831-1.0089920-rdf.json
Turtle: 831-1.0089920-turtle.txt
N-Triples: 831-1.0089920-rdf-ntriples.txt
Original Record: 831-1.0089920-source.json
Full Text
831-1.0089920-fulltext.txt
Citation
831-1.0089920.ris

Full Text

A B S O R P T I O N C O M P E N S A T I O N I N S E I S M I C D A T A P R O C E S S I N G By Changjun Zhang M.E., University of Petroleum in Beijing A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in THE FACULTY OF GRADUATE STUDIES T H E F A C U L T Y O F S C I E N C E D E P A R T M E N T O F E A R T H A N D O C E A N S C I E N C E S We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A November 2000 © Changjun Zhang, 2000 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Earth and Ocean Sciences The University of British Columbia 129-2219 Main MaU Vancouver, Canada V 6 T 1Z4 Date: Abstract As seismic waves propagate through the earth, the anelasticity of the medium will cause energy dissipation and waveform distortion. This phenomenon is refereed to as absorption attenuation. High resolution images require that the effects of absorption are understood and estimated. Estimated attenuation parameters can be used to improve the interpret-ation of seismograms. For these reasons, absorption is an important topic in seismic data processing. The absorptive property of a medium is described by a quality factor which describes the energy decay and determines a velocity dispersion relationship. The quality factor and the velocity govern the propagation of seismic energy in the earth. Velocities determine the arrival times of reflections and quality factors determine their frequency contents. Both of them can be estimated from common mid-point ( C M P ) gathers. Quality factors can be estimated from prestack common midpoint gathers, or a post-stack single trace. B y assuming that the amplitude spectrum of the seismic source signature may be modeled by that of a Ricker wavelet, an analytical relation between a quality factor and the seismic data peak frequency variation with time has been de-rived. This relation plays a central role in quality factor estimation problems. Tests on both synthetic and real data show that the analytical relationship can be used to extract correct quality factors whether the data are one-dimensional or two-dimensional. To remove the effect of absorption, conventionally, inverse Q filtering is used. The stability of inverse Q filtering is always affected by high frequency noise whether it is implemented in the time or the frequency domain. Formulating the absorption compen-sation as an inverse problem, and solving the inverse problem iteratively, the resolution u of seismic data is upgraded step by step. The inverse scheme helps to overcome the instability problem which is a natural drawback of common inverse Q filtering. Absorption compensation can also be considered as a part of the process of seismic data migration. Standard migration is the adjoint of forward modeling. High fidelity images can be obtained by using a least squares migration (LSM) scheme. Theoretical analysis shows that least squares migration is a promising way to implement migration and absorption compensation for a single processing step. m Table of Contents Abstract ii List of Figures vii Acknowledgment x 1 Introduction 1 2 Visco-Elastic Media Models, Mechanisms of Absorption 6 2.1 Physical Hypothesis 6 2.2 Mathematical Description 8 2.2.1 The Voigt Model 9 2.2.2 The Generalized Linear Solid Model 11 2.2.3 Describing Absorption by Quality Factors 12 2.2.4 The Effect of Absorption on Green's Function 16 3 Estimation of Quality Factors from Poststack Seismic Data 19 3.1 Short Window Time Variant Fourier Analysis 21 3.1.1 Spectral Ratio Method 21 3.1.2 Time-Frequency Analysis 23 3.1.3 Windowed Fourier Transform and Wavelet Transform 23 3.1.4 Quality Factor Estimation from WTVS 25 3.2 Quality Factor Estimation by Nonlinear Inversion 32 3.2.1 Quality Factor Estimation from the Viewpoint of the Wave Equation 34 iv 3.2.2 Solutions of the Nonlinear Problem 37 4 Absorption Compensation Theory and Practice 40 4.1 Inverse Q Filtering Revisited 41 4.2 Time-Variant Deconvolution as an Inverse Problem 47 4.2.1 Noise Problem in Deconvolution 47 4.2.2 Basic Theory of Statistical Inversion 48 4.2.3 Absorption Compensation and Inversion . . ; 51 4.2.4 Cauchy-Gauss Model 51 4.2.5 Implementation 53 4.3 Examples 55 5 Estimation of Quality Factors from C M P Gathers 62 5.1 Basic Assumptions 62 5.2 Absorption and Peak Frequency Variation 63 5.2.1 One layer case 65 5.2.2 Multi-layer case 68 5.3 Examples 70 5.4 Observations and Discussion 74 6 Migration and De-Absorption 77 6.1 Migration 78 6.2 Migration, the Adjoint Operation of Modeling 79 6.2.1 Basic Concepts of the Adjoint Operation 79 6.2.2 Gazdag modeling and migration 81 6.3 Least Squares Migration (LSM) 83 7 Conclusions and Discussions 88 v References 91 Appendix 99 A Absorption Described by a Dispersion Relation 99 v i List of Figures 1.1 (a) A seismic wavelet traveling through elastic media, its shape remains unchanged; (b) A seismic wavelet traveling through absorptive media, its high frequency components decrease faster than its low frequency compo-nents, causing the distortion of the wavelet shape 4 1.2 Comparison of the results of depth migration: (a) Elastic media; (b) Ab-sorptive media 5 2.1 Description of visco-elastic model by 3 parameters: velocity v, density p and quality factor Q 14 2.2 Evolution of a unit impulse with time in an absorptive medium 18 3.1 (a) A synthetic trace without absorption; (b) The windowed time variant Fourier analysis. The ridge of the spectrum is almost along a constant frequency line 28 3.2 (a) A synthetic trace with absorption; (b) The time variant Fourier ana-lysis. The decay of peak frequencies with time can be observed from the ridge of the spectrum 29 3.3 (a) A real seismic trace; (b) The WTVS; (c) A wavelet decomposition. . 30 3.4 (a) A WTVS of a synthetic signal, controlling points are marked with circles; (b) The spectrum of the wavelet at time i = 20 ms; (c) The spectrum of the wavelet at time t = 180 ms; (d) The logarithm of the spectral ratio 31 vii 3.5 a) A WTVS of a real seismic trace as in Figure 3.3, the picked points are marked by circles and connected by straight lines. Quality factors are calculated from the slope of each straight line segment, (b) The estimated Q's, the Q value of the last layer is infinite 33 3.6 The variation of the JC\ error vs. Q 39 4.1 De-absorption test on a synthetic noise-free data set: (a) A synthetic sec-tion with absorption (constant Q); (b) The section after absorption com-pensation. For noise-free data, both amplitude attenuation and phase distortion can be corrected completely 43 4.2 (a) A synthetic seismic trace with absorption; (b) After pure phase cor-rection, the minimum phase responses become zero phase 45 4.3 The amplitude of a inverse Q filter and its Gaussian approximation (the curve with circles) 46 4.4 Wavelet removal from a synthetic signal using Cauchy-Gauss model, (e) The original input trace; after 1 iteration the sharper responses is in (a); (b), (c) and (d) are inverse results after iteration 2, 3 and 4 respectively, they are almost the same 56 4.5 (a) A synthetic signal with 20 percent noise; (b) The result of absorption compensation using an inverse method, the rightmost trace is plotted as a reference 57 4.6 Convergence process of the inverse scheme; (e) is a trace from 4.5a. (a), (b), (c) and (d) are inverse results after iteration 1, 2, 3 and 4 respectively. 58 4.7 (a) Synthetic signals with 25 percent noise; (b) The result of absorption compensation using an inverse method 59 4.8 (a) A real seismic section; (b) The output of absorption compensation . . 61 vm 5.1 (a) An event in a CMP gather is generated by a reflector in an absorptive medium, where amplitude of each trace has been normalized, (b) Ampli-tude spectra of original source signature, trace 1, trace 11, and trace 21 which are indicated by [0], [1], [11] and [21] respectively 66 5.2 (a) A synthetic CMP gather of two reflections with absorption, the amp-litude has been normalized, (b) Real Q's in layer 1 and layer 2, it is unknown in layer 3. (c) computed Q values 71 5.3 Variation of peak frequencies of two reflections. Solid line is that of the first layer. Dashed line is that of the second layer 72 5.4 Experiment to show the effect of inverse Q filtering, (a) A synthetic CMP gather with 3 reflections; (b) Extracted Q values, only two layers have been separated; (c) The result of absorption compensation 73 5.5 (a) A common shot gather; (b) The common shot gather after inverse Q filtering. The Q curve has not been shown here 75 ix Acknowledgment I would like to express my deepest thanks to my supervisor, Dr. Tad Ulrych, for his moral and financial support throughout my stay at UBC. His invaluable academic advice contributes to every aspect of this thesis and will contribute to my career. I must also thank Sam Kaplan, Daniel Trad and James Irving. They kindly proofread my manuscripts and helped me correct many errors in both English and science. The courses I have taken from Drs. M. Bostok, B. BufFett and D. Oldenburg at UBC are very useful for my academic pursuit. I should say thanks to them. In particular, I am grateful to Drs. M. Bostok and B. BufFett for their participation as committee members and examiners and their constructive remarks. My stay at UBC and the research that led to this thesis was supported by a UBC international scholarship and grants from the Companies that support the Consortium for the Development of Specialized Seismic Techniques. Chapter 1 Introduction In seismic exploration, we commonly model the Earth as an ideal elastic medium, that is, each layer of the Earth can be fully described by two parameters: velocity and density (if only p-wave is considered), or even simply by just one parameter (in the convolution model): reflectivity. If the earth fits these models, then the frequency contents of the seis-mic wavefield would remain constant with depth; only variation of wavefield magnitude would be observed. In practice, we observe that both the frequency content and amplitude of seismic records vary with depth. The Earth is not an ideal elastic medium, and as a result, seismic waves traveling through the Earth experience attenuation and dispersion; high frequency components become weaker as travel-time increases. This phenomenon is caused by absorption. The attenuation phenomenon has long been studied, a few references being, Ricker (1953), Futterman (1962), White (1983), etc. Absorption is caused by anelasticity of rocks and may be described by quality factors (Q's) of layers. Absorption attenuates the amplitude and distorts the frequency components of prop-agating seismic wavelets. As is illustrated in Figure 1.1, the shape of a wavelet traveling through absorptive media is changed, it becomes broader; in contrast, a wavelet, trav-eling through elastic media, maintains its shape. Absorption decreases the resolution of seismic data. In addition, absorption distorts migrated earth structures. An example is shown in Figure 1.2. The images in Figure 1.2(a) and Figure 1.2(b) are two depth migrated sections. The original velocity models for the two sections are the same. In 1 Chapter 1. Introduction 2 doing forward modeling, two data sets (common shot gathers) have been created. One is without absorption, another one is with absorption attenuation. The image in Figure 1.2(a) is the result of migration from data without absorption and the image in Figure 1.2(b) is the result of migration from data with absorption. We observe a resolution decrease and mis-positioning of the structures in Figure 1.2(b), because of absorption in the corresponding synthetic data. In recent years, the inclusion of the second-order effects, such as absorption and anisotropy, into seismic processing schemes has become more important than ever. Finite-difference modeling in a visco-acoustic medium has been developed (Carcione, 1993), 3-D prestack migration in anisotropic media has been performed (Dong and McMechan, 1993), the inclusion of absorption effects into the prestack migration scheme has been demonstrated by Mittet et al. (1995) in a synthetic data set, and a technique of post-stack migration with attenuation compensation has been given by Song (2000) which has been applied to real data processing using multi-scale resolution propagation focusing. Understanding, estimating and compensating for the attenuation of seismic waves is important in the quest for improving the resolution of seismic images, better interpreting the effects of AVO (amplitude varies with offset) and inverting seismic data for material properties. Absorption of earth materials has long been an important research topic, from ex-plaining its mechanism, estimating absorptive properties, removing its effect, and uti-lizing absorptive properties to locate and describe hydrocarbon reservoirs. As already mentioned, there are many references dealing with these topics. In seismic data processing, a prerequisite for compensating for the effect of attenuation is the knowledge of quality factors ((J's) of the subsurface media. According to Dasgupta and Clark (1998), methods of estimating Q from surface seismic data have not been well developed thus far. Some research has been published concerning the estimation of Q Chapter 1. Introduction 3 from vertical seismic profiling (VSP) or from cross well data. Tonn (1991) investigated ten different approaches of estimating Q from bore-hole seismic data, almost all of them using the amplitude information of received signals. The problem with the use of amplitude information is that, due to the ubiquitous presence of noise and other factors, such as geometrical spreading and scattering effects, such information is not always robust and therefore, not always useful. Quan and Harris (1997) presented a method for estimating seismic attenuation based on the frequency shift inherent in VSP data. Since the centroid of the signal's spectrum experiences a shift to lower frequencies during propagation, they were able to develop a relationship between Q and the centroid of the amplitude spectrum which was represented by a Gaussian, boxcar or triangular shape. To remove the effect of absorption, we can use time variant deconvolution (Griffths, 1977), or inverse Q filtering. Inverse Q filtering is the more popular of these two ap-proaches because in the related algorithms, the dispersion relation is easier to include than in the time variant deconvolution implementation. The dispersion relation describes the relation of velocity with frequency. It is controlled by the absorptive property and velocity of a medium. Inverse Q filtering can be implemented in different ways, e.g. in the frequency domain (Hargreave et al., 1989 ), in the time domain (Varela et al., 1993), or it can be partially implemented by way of a phase correction (Robinson, 1979). In this thesis, firstly, as the basis of this research, absorptive mechanisms are in-troduced which can be applied in different situations. A definition of quality factor Q is presented. Using this Q, a linear dispersion model is adopted because it is appropriate not only for handling the problem of estimating quality factors from both prestack and poststack seismic data, but also for absorption compensation. Secondly two methods are introduced to estimate quality factors from poststack seis-mic data, one using a windowed time-variant Fourier analysis, another using an iterative optimization technique. Chapter 1. Introduction 4 Figure 1.1: (a) A seismic wavelet traveling through elastic media, its shape remains unchanged; (b) A seismic wavelet traveling through absorptive media, its high frequency components decrease faster than its low frequency components, causing the distortion of the wavelet shape. Thirdly, when implementing absorption compensation, I formulate the time-variant deconvolution problem as an inverse problem using statistical inverse theory. Using a sparseness constraint, a stable solution for the inverse problem can be found. As an important step in extending poststack absorption compensation technique to prestack seismic data processing, Chapter 5 concentrates on estimating quality factors from common mid-point ( C M P ) gathers or common shot gathers ( C S G ) . Also, a relation has been developed to connect the quality factors with peak frequency variations. As a future research direction, in Chapter 6, a prestack migration with absorption compensation scheme is briefly introduced which is based on least squares migration theory. This topic is not explored in great detail in the thesis. Chapter 1. Introduction 5 Figure 1.2: Comparison of the results of depth migration: (a) Elastic media; (b) Absorp-tive media. Chapter 2 Visco-Elastic Media Models, Mechanisms of Absorption The propagation of seismic waves in real media is in many respects different from seismic wave propagation in an ideal solid. The anelasticity of the medium will cause dissipation of seismic energy, thus decreasing the amplitude and modifying the frequency content of propagating waves. This attenuation and dispersion of seismic waves is closely related to rock properties (White, 1983). Therefore, these effects are important in exploration geophysics since they may allow us to extract more detailed information about the sub-surface and to construct images with better resolution. In this chapter, two points will be stressed: what causes the absorption phenomenon and how absorption is described in terms of wave theory. 2 .1 Physical Hypothesis A great deal of attention has been given to the study of the attenuation of seismic waves (Sato, 1998). Loss parameters for rocks have been measured by many techniques over a wide range of frequency and environmental conditions. Physical mechanisms for energy loss have been proposed and described mathematically in varying degrees of detail. There are two mechanisms that are usually considered to cause seismic attenuation, scattering and intrinsic mechanisms. Scattering redistributes wave energy within the medium but does not remove energy from the overall wavefield. Conversely, intrinsic at-tenuation refers to various mechanisms that convert vibrational energy into heat through 6 Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 7 friction, viscosity, and a thermal relaxation process. Here intrinsic attenuation or ab-sorption is our research subject. Absorption is a comprehensive effect caused by several factors. Among them, two factors are considered the most important. One is the internal friction resulting from the relative slide of rock matrices along the contact of solid matrices in rocks during wave propagation. Another is the viscosity between the solid matrix and the pore fluid, part of the traveling energy from seismic source is transformed into relative movement and dissipates inside and among the rock grains. Although absorption mechanisms have long been studied (Ricker, 1977, White, 1983, Mavko et al. 1979; Toksoz and Johnston, 1979 ), there is no general consensus as to dominant loss mechanisms . In different environments, a different mechanism may play a dominant role. Many scientists (Ricker, 1977, Aki and Richards, 1980) considered thermo-elasticity as the most viable model to explain intrinsic attenuation at lithospheric temperatures be-cause the required scales for rock grains and cracks along with the amount of attenuation caused by thermo-elasticity are in closest agreement with observations. Rocks have microscopic cracks and pores which may contain fluids. These cracks can have profound influence on the propagation velocity of P and S-waves. Biot (1956a, b) analyzed wave propagation in isotropic porous solids where the coupling of motion between the fluid and the solid matrix was considered. He arrived at expressions for attenuation due to follow of fluids within non-connecting pores initiated by elastic waves. White (1983, P 131) studied Biot's models and concluded that the attenuation predicted by Biot's model is extremely small for frequencies less than 100Hz. Many questions remain unanswered. A major complication is the wide range of properties possessed by earth materials. Conclusions drawn from data of oil-reservoir rocks may not be applicable to materials in the upper mantle. Hence, it is difficult to make a choice among the proposed mechanisms even if this were desirable. In fact, no Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 8 one mechanism can be expected to describe losses in all rocks under all environmental conditions. In spite of these questions and uncertainties, several characteristics of absorption in rocks seem well established. Absorption results in both attenuation and dispersion (Liu et al. 1976). For various rock types and for a wide range of environmental conditions, data indicate constant Q, or attenuation approximately proportional to frequency, for both shear and compressional waves. In many of these experiments, small changes in elastic constants or velocities are observed that are consistent with the measured attenuation and the demands of causality (White, 1983). These characteristics are the basis for a mathematical description of absorption in wave equations which will be introduced in the following section. I have only stated several possible explanations of absorption phenomena. The actual physical mechanism still requires further research. 2.2 Mathematical Description Absorption has been considered an important parameter to measure and characterize in sedimentary rocks for the purpose of petroleum exploration and this has led to efforts to develop models explaining the observed attenuation in sedimentary rocks (Mavko et al. 1979; Toksoz and Johnston, 1979). Many proposed intrinsic attenuation models are relaxation mechanisms having characteristic relaxation times that depend on the physical dimensions of studied objects. In dynamic mechanics, the behavior of wave traveling in a medium is derived from Hooke's law which relates stress and strain. According to Du (1996), generally speaking, there are two kinds of attenuation models built from the viewpoint of dynamic mechanics. One is the Voigt model and the other is the generalized linear solid model. Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 9 2 . 2 . 1 The Voigt Model The Voigt model is a widely investigated method of introducing losses that has the advantage of yielding a linear wave equation which can be solved for arbitrary time dependence. The assumption is made that stresses are directly proportional to strain rates, as well as to the components of strain themselves. This assumption was introduced independently by Stokes, Kelvin and Voigt, and its implications have been investigated by White(1983) and Ricker (1977). This kind of medium is commonly called a Voigt solid and the resulting equation is called Stokes equation (Ricker, 1977). In elastic mechanics, using Einstein notation, Hooke's law is expressed as (Tij = \65ij + 2/zejj where cr^ - is the stress tensor, is the strain tensor and 6ij is the Kronecker delta function. A and u are the two Lame's constants. 6 is the volume strain. If there is no body force, this constitute relation with Newton's Law _ d2Uj lead to the equations of motion in term of particle displacement d2u (\ + ^V9 + aV2u = p—, (2.1) where u is a vector of displacement, | ^ is its second derivative with respect to time and p is density. To a plane wave, if we assume the displacement is parallel to the Z-axis (ux and uy are both zero) and that uz is independent of x and y (In following chapters of this thesis, u is referred as uz when 1-d wave propagation is considered ). Equation (2.1) reduces to Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 10 In absorptive media, the constitutive relation must incorporate the fact that stress is not only proportional to strain, it is also proportional to the time derivative of strain. The stress-strain behavior in absorptive media is described by a modified Hooke's law which includes strain-rate terms: <ri5 = (A + ^Q-^Sij + 2(u + H'-Q^eij, where A' and a' are perturbations of Lame's constants, A and /x, due to absorption. Equation (2.2) in an absorptive medium becomes If we replace Lame's constants by Young's modulus, M = A + 2u, M' = A' + 2u' and transform equation (2.3) into the frequency domain. We have (M + iu;M')^- = -p^Uz, (2.4) where Uz = U(z,u>) designates the space-dependent wavefield at any angular frequency ui. A wavefield satisfying equation (2.4) has the following form U(z,w) = U0e±Gz, anc G pui2 1/2 (2.5) ( M + iuM') This complex propagation constant G may be expressed in terms of attenuation ap and phase velocity vp as G = ap + iuj/vp (2.6) Comparing equation (2.5) and equation (2.6) and defining u>o = M/M', it can be shown that when u>2 < u^, attenuation increases as the square of the frequency and velocity is Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 11 approximately constant (Ricker, 1977). They are expressed as following - r w ° i ^ 2 - r w° I^ 2 ap~[2(M/p)V2l^-[2vpLl and M 1 / 2 P However, many measurements have indicated a first power dependence with frequency. In addition, as pointed out by White (1983, P87), the Voigt model also violates causality. 2.2.2 The Generalized Linear Solid Mode l The constitutive relation between stress a and strain e for viscoelasticity could be for-mulated simply in the frequency domain as (Carcione, 1995) er(u>) = M(u>)e(u>). Here M(u>) is the complex, frequency-dependent, viscoelastic modulus. A generalized linear solid body can be described by (Malvern, 1969) where Mu is the un-relaxed bulk modulus, aj is a constant coefficient, u)j is the jth relax-ation frequency and AMj is the difference between the un-relaxed bulk modulus and the contribution to the bulk modulus due to the j th relaxation mechanism. These parameters can be adjusted to fit any realistic linear absorption law. In synthetic computation, using finite difference method, viscoelastic modeling of this model can be implemented at an expense which is only a little more than that of acoustic modeling (Emmerich and Korn, 1987, Ribodetti et al, 1995, Carcione et al, 1988, Carcione 1995). The two models described in this and the previous section are constructed by assuming that anelasticity of rocks will result in a constitutive relation that is different from that of Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 12 elastic models. Absorption may also be described by directly assuming velocity dispersion relations. 2.2.3 Describing Absorption by Quality Factors Velocities of rocks are determined by rock properties, and variations of medium elastic constants will affect the velocity. When including absorption into the wave equation, it is very straightforward to use complex wavenumbers or complex velocities. The absorptive property determines a real velocity dispersion relation or a complex dispersion relation depending on different applications (Mittet et al., 1995; Claerbout, 1987). The quality factor Q, sometimes also called the internal friction or dissipation factor, is an intrinsic property of rock. Formally, it is defined as a dimensionless measure of the anelasticity which is given by where E is the energy stored at maximum strain in a volume, and — AE is the energy loss in each cycle because of anelasticity. Equation (2.7) means that Q~X is the portion of energy lost during each cycle or wavelength. For a medium with a linear stress-strain relation, the amplitude A of a wave at a frequency is proportional to E1?2 (Aki and Richards, 1980), therefore I = <2-8> Q -KAQ where A0 is the amplitude at the start of a cycle and AA represents the amplitude decay in a cycle. It is observed that amplitude varies with distance because of absorption. We can rewrite A as a function of distance A(z), and AA = (dA(z)/dz)\, (2.9) Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 13 where A is the wavelength given in terms of of u> and phase velocity v(w) by A = 2TTV(UJ)/UJ. (2.10) From equation (2.8), (2.9) and (2.10), it can be shown that A(z) = A0exp (2.11) 2v(u)Q\ Under the constraint of causality, and frequency independent Q, we can obtain the dispersion relation (Aki and Richards, 1980), Hfcal = i + 1 M S L ) . (2.12) From observations of exponentially decaying values of A(x), we use equation (2.11) to define the value of spatial Q. Spatial quality factors are commonly used to study the absorption characteristics of media. From equation (2.11) and (2.12), we can observe that wave propagation in an ab-sorptive medium is completely specified by two parameters, Q and a phase velocity v0 at an arbitrary reference frequency LOQ. To describe absorption, these two equations are necessary. From the above discussion, we can also observe that when describing wave propaga-tion in absorptive media, quality factors and velocities have analogous uses. Using this concept, a layer can be described by three parameters: velocity, density and a quality factor as shown in Figure 2.1. In the following chapters, we will observe as well that quality factors can be estimated in ways very similar to velocity estimation. Another method to formulate the velocity dispersion relation was proposed by Kajar-tanson (1979) and Claerbout (1987). A basic model arises when the dispersion relation is defined by the equation -iui u>0 ('—iu>\1-e , „ -~N Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 14 «i , pi, Qi . . . -* V2, P2, 0,2 . . . •* « 3 , P 3 , ^ 3 . . . ~k V4, PA, QI Figure 2.1: Description of visco-elastic model by 3 parameters: velocity v, density p and quality factor Q. Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 15 Equation(2.13) models the so-called causal, constant Q attenuation also, and it can be derived — ~ tan(7T£j. The estimation of Q is equivalent to the estimation of e. This dispersion relation was used for absorption compensation by Kajartanson (1979). Equation (2.13) models a complex velocity dispersion relation. For e = 0, equation (2.13) gives a real constant velocity which is the case of elastic media. In this thesis I adopt the dispersion relation of equation (2.12). As pointed out by Aki and Richards (1980), equation (2.12) is an important result since it appears to be a good approximation for a variety of attenuation laws in which Q is effectively constant over the seismic frequency range. Besides a reference velocity v0, this visco-elastic model assumes attenuation and dispersion are determined by only one more parameter Q. Theoretically, this linear model (equation (2.11) and (2.12)) is superior to the Voigt model which has been contradicted through practical observations. Implementation is also easier than the generalized linear solid model which requires the adjustment of two parameters when being applied to different media. The linear model of equation (2.11) and (2.12) gives us insight into the development a method to estimate quality factors from CMP records. Based on this model, a method has been developed in Chapter 5 to estimate interval Q's based on the assumption that the amplitude spectrum of a source signature has an analytical expression with a dominant frequency (e.g. that of Ricker wavelet). In the following, we will examine the shape of the impulse response of the dispersion model of equation (2.12). Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 16 2.2.4 The Effect of Absorption on Green's Function As mentioned before, the basic assumptions used to describe absorption behavior are constant Q, linearity and causality. Considering a unit impulse as input at location z = 0, each Fourier component of the impulse, C + oo r+oo / S[t - z/v(uj)}eiajtdt = exp[iu>z/v(w)], J oo will now be attenuated by a factor of exp[—a(u>)z] with a(<*>) = 2V^Q • The propagating pulse shape G(z,t), thus, has a Fourier transform of elfcz, and k is complex which is given by k = - ^ - r + ia(u>). (2.14) V{OJ) Then if G(z,t) - 0; for t < z/v(oo), it can be shown that (Aki and Richards, 1980) U W +H[a(u>)]. (2.15) v(u>) v(oo) This equation is equivalent to the Kramers-Kronig relation in electro-magnetic theory (Futtermann, 1962). Here, v(oo) is the limit of v(UJ) as u> —* oo and H[a(oo)] is the Hilbert transform of the attenuation factor. From a mathematical point of view, there is no Hilbert transform pair for which this relation can be satisfied with constant Q. We must tolerate a frequency dependent Q which is effectively constant over the seismic frequency range. This idea has been developed by Azimi et al. in 1968 ( Aki and Richards, 1980), who derived a dispersion relation of equation (2.12) which satisfies the requirement of causality. Let us now examine the unit impulse response in absorptive media in the one dimen-sional case. Seismic waves traveling through elastic media satisfy the elastic wave equa-tion (2.1). If the medium is visco-elastic, the waves still satisfy equation (2.1), but with Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 17 frequency dependent velocity as described by the dispersion relation ( equation(2.12)) V2U -UJ v(ujy (2.16) In the one dimensional case using the phase shift concept, the wavefield may be extrap-olated from one depth z to another depth z + Az using U(z + AZ,UJ) = U(z,uj)e ikAz (2.17) where k is complex wavenumber as derived from equation (2.14). It has the form (2.18) V(UJ)K 2Q' Substituting equation (2.18) into equation (2.17), and rearranging the terms, we get an approximate equation (see Appendix A) . UJ U(z + AZ,UJ) = U(z,uj)A(Az,uj)ex.j>(i—Az), (2.19) where v0 is a reference velocity, and A(AZ,UJ) = A1(UJ)A2{UJ) (2.20) is an absorption related term in which Ai(u>) = exp( —) U ; VK 2Qv0) and A2{w) = exp UJAZ UJ -z^r—ln( —) Setting we obtain •KQV0 UJ0 J U'(Z,UJ) = U(Z,UJ)A1(UJ)A2(UJ) . UJ U(z + AZ,UJ) = U'(z,uj)exp(i—Az). (2.21) (2.22) (2.23) (2.24) Chapter 2. Visco-Elastic Media Models, Mechanisms of Absorption 18 t=0.0 (s) t=0.8 (s) t=1.6 (s) t=2.4.0 (s) t=3.2 (s) Figure 2.2: Evolution of a unit impulse with time in an absorptive medium Absorption compensation is performed by obtaining U'(z,u>) which is extrapolated from the observed signal U(z,u/) using the operator A(AZ,UJ). As expressed in equation (2.23) ( also in equation (A.8)) absorption compensation includes the recovery of both amplitude and phase. Using this equation, we can examine the impulse response of the absorptive media with a Q = 30. A unit impulse originates at time t = 0 second will become wider as time increases. This is illustrated in Figure 2.2. From this figure we can notice that the responses are causal. The amplitude of the impulse responses at different times has been normalized for convenience of illustration. Chapter 3 Estimation of Quality Factors from Poststack Seismic Data Anelastic absorption attenuates seismic amplitude and distorts the seismic waveform. It is important, therefore, that its effects must be ehminated. Attenuation has also been recognized as a significant seismic parameter which when quantified may be used in improving the interpretation of seismograms. For these reasons, reliable Q-estimates are needed. Because of favorable source-receiver geometry, VSPs are ideal experiments for estim-ating Q (Stainby and Worthington, 1985; Dasgupta and Clark, 1998), and often Q has been computed successfully from such observations. In VSP data, it is easy to compare direct arrivals recorded at two different depths or different reflections recorded at a single depth. From the comparison, quality factors of substructures can be estimated. Several methods have been developed (Tonn, 1989) to estimate quality factors. These methods include the analytical, the rise-time, the spectral ratio and the wavelet fitting methods (Liang 1988; Cheng 1992). The majority of these approaches are used in VSP data processing (Harris et al, 1997). In seismic exploration, as surface seismic data are processed, VSP observations are not always available. Therefore, quality factor estimation techniques are needed for surface seismic data processing. A method for quality factor estimation from prestack common midpoint (CMP) gathers or common shot gathers (CSG) is discussed in Chapter 5. In this chapter, two methods to estimate quality factors from stacked seismic traces are discussed. These are firstly, estimating Q from a windowed time-variant spectrum 19 Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 20 (WTVS) and secondly estimating Q using a nonlinear inverse method. It is not easy to estimate quality factors from a single stacked trace, since the problem of estimating one dimensional parameters v(t) or Q(t) from one dimensional observations u(z = 0,t) does not yield unique solutions. In addition, in seismic data processing, stack processing after normal moveout (NMO) destroys the original absorption character by summing together signals which experience different trajectories. In short, a stacked trace has a distorted attenuation signature and therefore quality factor estimation from prestack data is needed. Nonetheless, to do de-absorption processing on poststack data is often necessary; hence, it is important to estimate quality factors from vertical incid-ence traces. In seismic data processing, absorption compensation is sometimes used to extrapolate the frequency band of poststack seismic data to increase resolution. The aforementioned methods (e.g. the spectral ratio method and the wavelet model-ing method) have also been used for quality factor estimation from surface seismic data (Liang, 1987; Dasgupta and Clark, 1998). Both methods are affected by factors other than absorption (Tonn, 1990). We will explore these factors further in the next section. In this chapter, as already mentioned, two methods will be discussed. In the first method WTVS, the frequency component of a signal at different time locations is ana-lyzed using a short windowed Fourier transform. From the signature of the variation in frequency, quality factors are determined. The second method is nonlinear inversion. It is analogous to using inversion to find velocities from poststack data. Assuming that there are well log data available, using an absorption-free synthetic trace generated from that well log data as a reference, a nonlinear inversion algorithm estimating Q's is presented. A synthetic data example shows that this well-constrained method is powerful for estimating quality factors from a single trace. The restriction is that the well log data are not always available. For the purpose of developing a method that can extract quality factors from poststack seismic Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 21 sections, the discussion below assumes that stacked seismic traces are vertical incidence traces. 3 . 1 S h o r t W i n d o w T i m e V a r i a n t F o u r i e r A n a l y s i s In geophysical exploration, observed signals are often time variant. Waves that reflect from the shallow subsurface are rich in high frequency components, whereas waves that return from the deep subsurface are dominated by low frequency components. In order to analyze such time variant signals, the windowed Fourier transform is commonly used. In the following, a spectral ratio method is presented to illustrate the importance of the time-frequency analysis in quality factor estimation, and then the wavelet transform (WT) and the windowed Fourier transform will be compared. 3 . 1 . 1 S p e c t r a l R a t i o M e t h o d Because absorption increases with frequency, analyzing the variation of the frequency spectrum over time can reveal the absorptive property of the media through which the seismic waves propagate. Assuming that at time t0, the amplitude spectrum of the seismic wavelet is B(f,to), at time t + to, the amplitude spectrum becomes B(f, to + t) and, based on equation (A.11), the relation between B(f,t0) and B(f,to + t) is B(f,t0 + t) = B(f,t0)e-=£, (3.1) where Q is the quality factor of the medium through which the wave propagates. Dividing both sides of equation (3.1) by B(f,to)) gives B(f,t) B(f) = e Ell Q (3.2) Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 22 Taking logarithms of both sides of equation (3.2), we have ln B(f) Setting Ar = ln [ 7 ^ ] , equation (3.2) becomes Q A r = { ~ ^ ) f = pf- ( 3 - 3 ) Plotting the ratio as a function of frequency yields a linear trend whose slope, p, is a function of Q, and thus Q can be determined as Q = =£. (3.4) P The spectral ratio method is based on this principle of using the amplitude spectrum information of seismic traces. Although this method has been used in VSP data pro-cessing since it is possible to pick first arrivals from different traces of VSP data (Tonn, 1990), it is not well established as to how to use it to extract quality factors from surface vertical incidence signals. The difficulty is three fold: the first is in separating reflection arrivals at different time locations, the second is in performing time-frequency analysis (frequency analysis at different times along a trace) and the third aspect restricting the application of this method in practice is that, besides absorption, there are other factors affecting the amplitudes of received signals. These factors include geometrical spreading, reflection and transmission effects, and AGC (automatic gain control) that is applied to balance the trace amplitudes. Even if the amplitude is precisely measured, because of the restricted time-frequency analysis resolution (Mallat, 1998), it is still difficult to cal-culate Q. In the following 1-D case, since we assume vertical ray-paths, the geometrical spreading effect is neglected. Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 23 3.1.2 Time-Frequency Analysis The spectral ratio method requires calculation of the amplitude spectrum at different times. First, we discuss how the accuracy of the calculated spectrum is restricted by the time-frequency analysis. A consequence of Fourier transforms (FT) being built from elwt is that scaling a function to be narrower in one domain scales it to be wider in the other domain. Scaling u> implies inverse scaling of t to keep the product u)t constant. This relation is clearly shown by the scaling property of the Fourier transform, fiat) - |V(-)- (3-5) \a\ a Broadening a signal in the time domain implies sharpen its spectrum in the frequency domain. The localization of a signal in both the time and frequency domains cannot be achieved using the FT. 3.1.3 Windowed Fourier Transform and Wavelet Transform The windowed Fourier transform was introduced in 1946 by Gabor to measure localized frequency components of sound. A real and symmetric window g(t) = g(—t) is translated by u and modulated at frequency u>: guA*) = ^g{t - u). The resulting windowed Fourier transform of a signal /(£) is /+oo f(t)g(t - u)e-^dt. (3.6) -oo This transform is also called the short time Fourier transform because the multiplication by g{t — u) localizes the Fourier integral in the neighborhood of t = u. In the windowed Fourier transform, the window function remains unchanged as it moves along time axis. Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 24 The resolution in time and frequency of the windowed Fourier transform depends on the spread of the window in time and frequency. The time-frequency localization of g can be modified by a scaling factor s. Let ga(t) = ^ 5 ( 7 ) represent the dilation. The choice of a particular scale s depends on the desired resolution trade-off between time and frequency. To analyze signal structures of different size at different localizations, it is necessary to use time-frequency atoms with different time supports. The wavelet transform de-composes signals by means of dilated and translated wavelets. A wavelet is a function Tp 6 L2(R) with zero mean /oo ij>(t)dt = 0. - 0 0 It is normalized (||V>|| = 1), and centered in the neighborhood of t = 0. A family of time-frequency atoms is obtained by scaling ip by s and translating by u. The wavelet transform has been developed for this purpose (Mallat, 1997). The wavelet transform of a signal f(t) at time u and scale s is /OO + 7, f(t)i>( )dt. (3.7) - 0 0 s Wf(u, s) is also called the wavelet coefficient. A set of wavelet bases which are wavelet functions at different scales is equivalent to a set of band pass filters. Orthogonal and complete wavelet bases can be constructed to perform multi-resolution analysis of signals (Mallat, 1997). The projections of a signal onto wavelet bases are equi-valent to outputs from different band pass filters. The wavelet transform has found its application in many areas, for example, image processing (Mallat, 1997), data compres-sion (Mallat, 1997), solution of nonlinear partial differential equation (Beylkin, 1998) and wavefield reconstruction (Song, 2000). Using wavelet analysis to estimate quality factor from ground penetrating radar data has been recently investigated by Irving (2000). In my experience of time variant spectral analysis of seismic signals, Windowed time-variant Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 25 spectrum (WTVS) is superior to wavelet transform (WT) in mapping the tendency of frequency variation with time. 3.1.4 Quality Factor Estimation from W T V S According to the theoretical analysis presented above, the absorption character of seismic wave propagation can be estimated by examining the variation of the spectral content of a seismic trace with time. The output of the wavelet transform is either wavelet coefficients or signal components at different scales. Neither yield frequency spectrum information explicitly. In contrast, windowed Fourier analysis produces the spectrum associated with each time window (Schoepp and Margrave, 1998). In absorptive media, absorption incurred dispersion causes the width of a seismic wavelet to increase with time. To allow the window to enclose the entire seismic wavelet at each time location, a time variant window function g{t) is used. The width of g(t) increases linearly with time. The g(t) functions used in the following numerical tests are boxcar functions. The WTVS computation initially applies a narrow window to the input trace and calculates the ordinary Fourier spectrum of the windowed data. The window is then translated and widened successively with time along the trace and the Fourier transform is calculated for each new position of the window. A WTVS spectrum can be displayed as a grey level or contour plot as in Figure 3.1 with frequency as the horizontal coordinate and time as the vertical coordinate. Each row of the plot corresponds to the spectrum of the input data at a specific time. The relative detail in time and frequency in the WTVS is related to the window size and frequency intervals. Figure 3.1a shows a synthetic seismic trace without absorption, and Figure 3.1b shows a WTVS spectrum of the signal. The bandwidth remains almost constant and the peak Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 26 amplitude is constant with frequency. In contrast, for a synthetic signal with absorption, Figure 3.2a where a constant Q = 30 is used, the bandwidth generally decays with time as shown in Figure 3.2b, and is therefore, non-stationary. A real data example is shown in Figure 3.3, where Figure 3.3a is the real signal from a seismic survey. The windowed time variant spectrum is shown in Figure 3.3b, and Figure 3.3c illustrates a wavelet decomposition. We observe that the WTVS gives a clear indication of the trend of the spectral variation. Although wavelet components also illustrate the frequency evolution with different scales, the trend of frequency variation is less clear in comparison to WTVS. Picking the peak amplitude ridge from WTVS and fitting the ridge with a piece-wise straight line in the coordinates of / versus t , a quality factor can be calculated from each line segment using the formula: where / p i and fp2 denote the peak frequencies which is measured at two different times, £i and t2 on a line segment respectively. Details of the derivation of this formula can be found in Chapter 5. The spectral ratio method can also be used to estimate quality factors from WTVS. First, controlling points in the WTVS are picked. Second, comparing the spectra at the controlling time point locations and plotting the logarithm of the spectral ratio versus frequency yields a linear trend whose slope is a function of the quality factor as shown in equation (3.4). To illustrate the foregoing theory of Q estimation, both synthetic and real data have been considered. Figure 3.4a is a WTVS of a synthetic signal with Q = 30. Controlling points are picked and marked with circles. These points lies along a straight line because Q is constant. Figure 3.4b is the spectrum of the wavelet at time t = 20 ms, and Figure Q = * ( * 2 - * l ) / p 2 / p ] 2(/pi — fp2) (3.8) Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 27 3.4c is the spectrum of the wavelet at time t = 180 ms. These two locations are picked as controlling points. Figure 3.4d, is the logarithm of the spectral ratio of the two spectra. A linear trend can be observed. Either using equation (3.8) or using equation (3.4), a quality factor can calculated which is equal to 32.6. It is very close to the real value quality 30. As an example of quality factor estimation from real data, Figure 3.5a is the same WTVS as that in Figure 3.4b. Five controlling points are picked and marked with circles. The geological structure is divided into 4 layers. The calculated quality factors are plotted in Figure 3.5c. Layer 4 has a very large quality factor (approaching infinity), this means that this layer is elastic (It can be explained that as depth increases, rocks are more compacted and elastic because of increased pressure). Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 28 Figure 3.1: (a) A synthetic trace without absorption; (b) The windowed time variant Fourier analysis. The ridge of the spectrum is almost along a constant frequency line. Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 29 Figure 3.2: (a) A synthetic trace with absorption; (b) The time variant Fourier analysis. The decay of peak frequencies with time can be observed from the ridge of the spectrum. Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 30 (SU) 9WI| Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 31 Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 32 3.2 Quality Factor Estimation by Nonlinear Inversion In the previous section, we discussed using the windowed time variant Fourier analysis to estimate quality factors from a vertical incidence seismic trace. It is observed that the method is effective in finding the variation in the frequency content of a trace. Based on the one-dimensional wave equation, quality factor estimation can be formulated as a nonlinear inverse problem. Assuming a reference trace exists which is absorption free, absorptive parameters contained in a trace can be calculated. In seismic data processing, we often formulate the stacked seismic trace by means of the convolution model x(t) = b(t) * r(t), (3.9) where x(t), b(t) and r(t) are the seismic trace, the seismic wavelet and the reflectivity of the Earth, respectively. To obtain r(i) or b(t) from equation (3.9) is a deconvolution problem. A deconvolution problem can be solved using inverse methods (Claerbout, 1976), since it can be written in the form of a matrix multiplication as equation (3.10) or equation (3.11), x = B r (3.10) or x = R b (3.11) where x, b and r are vectors corresponding to x(t), b(t) and r(t) respectively. B is the wavelet matrix. R is the reflectivity matrix. If the medium is absorptive, a time variant attenuation operator should be included in the convolution model. Equation (3.9) becomes x(r,t) = a(r,t)*b(t)*r(t) (3.12) Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 3.3 0 10 20 30 40 50 60 1 1.2 1.4 1.6 1.8 2 Frequency (Hz) Figure 3.5: a) A WTVS of a real seismic trace as in Figure 3.3, the picked points are marked by circles and connected by straight lines. Quality factors are calculated from the slope of each straight line segment, (b) The estimated Q's, the Q value of the last layer is infinite. Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 34 where a(r, t) is an attenuation operator at a time r. Unfortunately, because a(r, t) is time variant and does not have a simple analytical form, equation (3.12) cannot be written in a matrix form with a(r, t) as a vector at the rightmost, so that, even if the original wavelet and earth reflectivity is known, equation (3.12) can not be solved for attenuation characters using inversion algorithms. In the next section, a nonlinear optimization scheme for Q estimation is devised based on the one-dimensional wave equation rather than the convolutional model of equation (3.12) . 3 . 2 . 1 Q u a l i t y F a c t o r E s t i m a t i o n f r o m t h e V i e w p o i n t o f t h e W a v e E q u a t i o n We assume that a non-absorptive reference trace is available which corresponds to a real seismic trace (at the end of this section, I will explain how to build the reference trace). Based on the one-dimensional wave equation, quality factor estimation can be formulated as a nonlinear inverse problem. The basic equation that we will use is equation (A.7) which is the one-dimensional absorptive wavefield continuation equation in the frequency domain. If depth z is sampled at interval Az, we can denote the depth z, the quality factor and the reference velocity at this depth by z = nAz, Q(Z) = Qn and v0(z) = v0(n), respectively. We also define the wavefield in the frequency domain by U(z,u) = Un{u>) Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 35 and U(z + Az,u) = Un+l{uj). From equation (A.7), we have Un+1(u>) = Un(u>)A(Az,uj)exp I i UJ v0(n) where A(Az,u>) = exp (--J^L-J) e X p tuAz , . UJ . — z — - — — l n ( — ) (3.13) (3.14) 2Qnv0(n)/ is the absorption related operator. Symbols used in equations (3.13) and (3.14) are the same as in Appendix A. Equation (3.13) can also be written in a non-recursive form, UJ Un(uj) = U0(uj)A(z,uj)exp(i zn), (3.15) with vave being the average velocity at frequency u>0 in a section of thickness z, and ujjAz A t \ TT / ujAz A(z,u>) = 11 exp , ^ | exp 3=0 UJ ln ( - ) *QjVo(J) "0 J An(uj). (3.16) Phase shift migration is done in two steps (Gazdag, 1978). Firstly, the surface wave-field UQ(UJ) is downward extrapolated to a depth zn using equation (3.15) to get Un(u>). Secondly, all frequency components of the wavefield Un(uj) are summed together to form un(t = 0), Un(t = 0) = J U0(uj)dw. (3.17) Inserting equation (3.15) into (3.17), we obtain un{t = 0) = I U0{uj)An(uj)ex.p(i-^—zn)duj. (3.18) For constant reference velocities, it can be derived that the sampling rate of a trace At = AZ/VQ, vave = VQ = VQ(J) and tn = zn/vo- Discretizing frequency u> = mAw, we Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 36 have the discrete version of equation (3.18) Mt = 0) - E C / n (^)A n (o; m )exp( i—z n )e { ^Aw. ^ave (3.19) Considering a numerical implementation, if we assume velocity information is available, unknown parameters in term An(uim) are Q's. If we put these <5's into a vector Q n = ( Q l , Q 2 , ' " , Q n ) and define An(n,m) = An(Qn,ujm) = An(u}m)Auj, We obtain ujt = 0) = ^cf n(^m)A n(Q n ,u>m)exp(iwm i n). Furthermore, writing equation (3.22), in matrix form, we have (3.20) (3.21) (3.22) ' ^(1,1)^!! •• A1(l,M)WMi ^ ' u0M ^ it:i > = 0 ) \ A2{2,1)W12 A2(2,2)W22 • •• A2(2,M)WM2 U0(u2) u2(t = 0 ) v AN(N,1)W1N AN(N,2)W2N • •• AN(N,M)WMN , v ttjv(t = ° ) y (3.23) where Wmn = e""mt™. The more explicit form of v 4 „ ( Q „ , u ; m ) is exp An(Qn,um) = Aw JJ exp i——MT—) TTQJVQ LUQ Equation (3.23) may be written in the following concise form, A U R = u. In equation (3.25), A is the coefficient matrix and is a function of Q's. (3.24) (3.25) Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 37 is a vector of a reference signal without absorption in the frequency domain. When velocity information is available, it can be extrapolated to difference depths u = [Ul(t = 0), u2(t = 0), • • •, uN(t = 0)]', whitch is the recorded seismic trace in the time domain. Equation (3.25) is a nonlinear equation, the parameters we want to invert for are in the coefficient matrix. A similar equation was used by Harris et al. (1997) who assumed a quality factor as a function of frequency, and different layers having the same quality factor function. In this way, the problem was formulated as a nonlinear equation to solve for Q(u>) using the SVD (singular value decomposition) method. In fact, the assumption of Harris et al. (1997) conflicts with many experimental results. Evidence indicates that quality factors in different media are different and, in the range of seismic frequencies, the quality factor of a medium remains constant (Turner and Siggins (1994)). Therefore, it is reasonable to assume that quality factors are a function of time or depth. 3.2.2 Solutions of the Nonlinear Problem We can use two methods to solve the nonlinear equation (3.25) for quality factors. One is the nonlinear optimization method. For the quality factors in a layered Earth, we expect to obtain a flat model in the inversion. Nonlinear optimization algorithms can be used to minimize the objective function </>(Q) = | A U r - u | + F Q (3.26) where F is a flatness operator. This nonlinear problem can be solved by many methods, e.g. simulated annealing or genetic algorithms. Another method is simple in its implementation. It is a scanning method, that is, we test a series of quality factors with fixed interval and compute the value of the objective Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 38 function, choosing the quality factor corresponding to the minimum misfit in the solution of the nonlinear problem. Similar methods are used in velocity and noise analysis in seismic data processing. This approach is also suitable for the quality factor estimation problem of equation (3.26). To test this nonlinear optimization scheme, an experiment has been performed on a synthetic trace. The trace used in this test is in Figure 3.2a and its corresponding reference trace is in Figure 3.1a. As shown before, these two traces have the same wavelet and the same reflectivity. The scanning process begins by setting Q = 10 and choosing an increment of 10. The absolute value error is calculated at each Q. The misfit function reaches its smallest value at exactly Q — 30. This value is equal to the actual quality factor used in the modeling. The curve for this scanning process is shown in Figure 3.6. A shortcoming of the scanning method is that it could involve extensive computation. In practice, if well data are available, the reflectivity function could be extracted from well data, and the wavelet could be extracted from close well seismic traces. With the reflectivity and the wavelet, a non-absorptive reference trace can be formed. This nonlinear inverse method could be useful in correlating well data with close well seismic traces. Using this optimization scheme, detailed absorption characters could be extracted from seismic data. Extrapolating the quality factors from close-well traces along horizons of a seismic section, a quality factor section could be obtained as is done in well-constrained acoustic impedance inversion (Martines et al. 1992). This could be very useful for reservoir description. Chapter 3. Estimation of Quality Factors from Poststack Seismic Data 39 Q values F i g u r e 3 .6 : T h e v a r i a t i o n o f t h e L\ e r r o r v s . Q Chapter 4 Absorption Compensation Theory and Practice Inverse Q filtering is used to remove the effect of absorption. The stability of inverse Q filtering will be affected by high frequency noise whether it is implemented in the time or the frequency domain. By formulating the absorption compensation as an inverse problem, a stable solution can be found. In seismic exploration, stacked seismic sections are interpreted in order to find earth structures that could be potential hydrocarbon reservoirs. The detailed interpretation of the earth structure or the sedimentation history relies on high resolution seismic data. Invariably, however, seismic data only have restricted resolution because of band-limited seismic wavelets and absorption. A number of techniques have been used to broaden the frequency band of seismic data. Besides deconvolution, the broad-band constraint inversion of seismic profiles is another widely used method which integrates seismic data with sonic-log data to extrapolate the frequency band to the high and low frequency ends (King, 1988; Martines et al., 1992; Du, 1996 ). Absorption compensation, simply referred to as de-absorption (McGinn and Dui-jndam, 1999), is specifically used to remove the effect of absorption attenuation and distortion of wavelet spectra. It is an important aspect of high resolution processing. De-absorption has received great attention among geophysicists and inverse Q filtering has been gradually incorporated into the seismic data processing flow. The main obstacle that needs to be overcome in the application of inverse Q filtering is its unstable nature. As pointed out by Claerbout (1985), in making earth images, earth dissipation might 40 Chapter 4. Absorption Compensation Theory and Practice 4 1 be compensated for by amplifying high-frequency energy during downward continuation. This might be done, just like migration, by shifting a wavefield by elkzZ, except that kz = yj(uj2/v2 — hi) would be replaced by something like ikz = (—iwY~€ (see equation (2.3)). In practice, however, no one would really do this, since it would severely amp-lify high frequency noise. In fact, many geophysicists have focused their research on developing algorithms that can avoid the amplification of random noise; among them are Robinson (1979), Hale (1987), Hargreaves (1990), Verala et al. (1993) and Song et al. (2000). 4.1 Inverse Q Filtering Revisited There are many methods available to perform absorption compensation, the most com-mon of these being time-variant deconvolution (Margrave, 1998 or multi-window de-convolution (Yilmaz, 1987)). The deconvolution problem involves time-variant seismic wavelets that are governed by the absorption properties of the earth. Different inverse filters are devised in different windows to take into account the non-stationarity of seis-mic wavelets. Another method to compensate for absorption is time-variant spectral whitening (Yilmaz, 1988). This process proceeds in three steps. First, a series of narrow band-pass filters are applied to a trace. Commonly the high frequency components have a faster decay rate than the low frequency components. Second, contents of the different frequency bands are examined, and a series of gain functions are computed to describe the decay rates for each frequency band. This is done by computing the envelope of the band-pass filtered traces. Third, the inverses of these gain functions are applied to each frequency band and the results are summed. The summed trace is the one with absorp-tion compensation. Both time-variant deconvolution and time-variant spectral whitening are done in the time domain, and are widely used in practice. Limitations of these two Chapter 4. Absorption Compensation Theory and Practice 42 methods are obvious. According to equation (A.7), absorption compensation problems can also be handled in the frequency domain. If our observed data are ideally noise-free, the de-absorption 1978). Figure 4.1 shows an example of ideal absorption compensation. The left panel in Figure 4.1 is a synthetic seismic section; absorption makes reflections become broader as time increases. As a result, the resolution of this section becomes worse from top to bottom. The right panel in Figure 4.1 shows the section after de-absorption. Reflections become sharp impulses again and the resolution becomes uniform. Because there is no random noise added to the synthetic data, absorption has been completely compensated for. In practice, since noise is always associated with real signals, this direct de-absorption scheme cannot be used. If an inverse Q filter does not distinguish between reflection signal and noise, it will alter the amplitude of both (Turner, 1994). This can be observed by examining the attenuation operator in equation (A.7). process can be implemented easily, in the same way as phase shift migration (Gazdag, A(Az,u>) = A1(u)A2{u)) (4.1) The absorption operator has two parts: the amplitude attenuation part Ai(u>) and the phase distortion part ^(u;) with the following expressions, respectively: 111 A 7. ^ H = e X p ( - 2 W and A2{yS) = exp — i •KQVQ U>0 The inverse of an attenuation operator A-\LAZ,U) = A-\U>)A-\U>), Chapter 4. Absorption Compensation Theory and Practice 43 Figure 4.1: De-absorption test on a synthetic noise-free data set: (a) A synthetic sec-tion with absorption (constant Q); (b) The section after absorption compensation. For noise-free data, both amplitude attenuation and phase distortion can be corrected com-pletely. Chapter 4. Absorption Compensation Theory and Practice 44 where anc = e x P i——— ln( —) TTQVQ UJQ (4.3) respectively. If Q is known, the phase distortion might be removed separately (Robinson, 1979; Hargreave, 1992) because dispersion correction can be done completely without any in-stability problem. The pure phase correction cannot expand the frequency band of a signal, but it does retrieve a symmetric unit impulse response. An example of phase correction is shown in Figure 4.2. A synthetic signal with con-stant absorption coefficient is shown in Figure 4.2a. The responses are minimum phase. Pure dispersion correction is shown in Figure 4.2b. All of the responses become sym-metric and sharper. Pure phase correction, however, is not enough to compensate for absorption attenuation. To remove the absorptive effect, amplitude compensation is ne-cessary. As mentioned before, if the direct inverse operator is used to compensate for attenuation, high frequency noise will be boosted. To ensure that noise is not unne-cessarily increased, it is important to modify the inverse operator. One way to avoid the instability is to use the band-limited version of the inverse operator i.e., to replace A^iyj) by Wo^A^1^). Here Wo(ui) represents a low-pass filter. An example of this method is to taper the amplitude response in equation (4.2) by multiplying the inverse filter by a time-dependent Gaussian function A1-1(t>a;)-*i4r1HG(t>W)> where G(t,uj) is a time-dependent truncated Gaussian filter (Hargreave, 1993). Figure 4.3 shows the Gaussian approximation to an exponential amplitude compensation operator. Chapter 4. Absorption Compensation Theory and Practice 45 a) 500 550 600 -1 1 1 1 1 1 1 -1 V. , J I . . . J 50 100 150 200 250 300 350 400 450 500 550 600 T i m e (ms) Figure 4.2: (a) A synthetic seismic trace with absorption; (b) After pure phase correction, the minimum phase responses become zero phase. Chapter 4. Absorption Compensation Theory and Practice 46 Frequency (Hz) Figure 4.3: The amplitude of a inverse Q filter and its Gaussian approximation (the curve with circles). To avoid noise amplification, the original monotonically rising function is replaced by a function which decays at high frequencies, and the choice of G(t,co) will be critical to the final results. Undoubtedly this method will compensate for attenuation of the low-frequency parts. However, it can not compensate for the loss of the high frequency energy. Another way of avoiding excessive noise boosting is to project the inverse operator onto a finite orthogonal wavelet basis by eliminating the high frequency component of the inverse operator (Song et al., 2000) where WJ(UJ) represents the wavelet basis at scale j and J is a finite integer. For complete decomposition, J approaches oo. The stability in these two schemes is achieved by Chapter 4. Absorption Compensation Theory and Practice 47 attenuating the inverse filter above a certain frequency, or by just discarding the frequency component above some scale. Using inverse theory, stability in the absorption compensation problem can also be achieved. The details of this will be discussed in following sections of this chapter. 4.2 Time-Variant Deconvolution as an Inverse Problem 4.2.1 Noise Problem in Deconvolution The instability problem is not unique to the absorption compensation problem. It is omnipresent in the implementation of deconvolution. In the following, we discuss this problem and optimum Wiener filter design. In seismic data processing, the recorded data in the one-dimensional case can be modeled by x(t) = w(t)*r(t) + n(t). (4.4) In this equation, x(t), w(t) and r(t) represent a recorded seismic trace, a seismic wavelet and a function of the earth's reflectivity. Assuming that w(i) is known, deconvolution aims at obtaining the reflectivity from the observed seismic trace. Commonly, however, the calculated reflectivity will be affected by the random noise term, whether or not w{t) is band limited. Assuming that f(t) is the inverse filter of w(t) (i. e. f(t) * w(t) = 6(t)), deconvolution results in r'(t) = r(t) + f(t)*n(t). (4.5) Whenever n(t) is not close to zero, the noise effect on the output reflectivity will be very severe because f(t) often has a long duration. To handle the problem of noise in seismic data processing, the inverse filter is often designed using a least squares scheme rather than directly from the original seismic wavelet. The designed filters are called optimum Wiener filters (Robinson and Treitel, Chapter 4. Absorption Compensation Theory and Practice 48 1980). Wiener filtering involves designing the filter f(t) so that the least squares error between the actual and desired output is minimum. Error L is defined as L = T/(d(t)-y(t))2, (4.6) t where d{t) represents the desired output of the filter. The actual output y(t) is the convolution of the filter f(t) with the input seismic trace x(t): y(t) = f(t) * x(t). (4.7) When equation (4.7) is substituted into equation (4.6) , we obtain L = E {^t ~ E /r*t-r) • (4-8) The goal is to compute the filter coefficients (/0, •-,/n-i) so that the error is mini-mized. The filter thus decreases the effect of noise. Optimum Wiener filters are designed to bring the predicted output as close as possible to the desired output using the £ 2 norm. As we will see in the following discussion, in seismic data inversion schemes, model para-meters are adjusted until the difference between the theoretical and observed values is a minimum. The philosophy in deconvolution and inversion is the same. This understand-ing encourages us to find a way to solve the time-variant deconvolution problem using inverse theory. 4.2.2 Basic Theory of Statistical Inversion Actual measurements of geophysical quantities are never exact. It is therefore essential for any practical scheme that solves inverse problems to take this into account. As discussed in the last section, direct deconvolution has inevitable instability. Formulating deconvolution or inverse filtering as an inverse problem has the advantage of allowing us to control what kind of solutions we wish to obtain. Another advantage of this formulation Chapter 4. Absorption Compensation Theory and Practice 49 is that it only uses the forward operator rather than the inverse operator. Many methods have been developed for constructing solutions to linear inverse problems in a way that takes into account the non-uniqueness and instability associated with the incompleteness and inaccuracy of observed data. To solve an inverse problem, two things are now required: a misfit measure, that is, a way of quantifying the disagreement between the observations and their counterparts; and a tolerance, which is a level of misfit that is considered acceptable so that when a model achieves that value, it may be said that the model adequately honors the demands of the measurements. In short, to solve an inverse problem, first an objective function should be constructed, and then a method should be used to minimize or maximize that objective function. This approach is commonly iterative, and the iterations will stop when the solution is acceptable. Inverse problems can be approached from the point of view of probability theory (Ulrych, 1998; Tarantola, 1987). It is common to consider measured data d as uncertain. That is, although there exist true data, we do not know them. The measured data can then be considered as random variables whose mean is that of the true value. Tradi-tionally, the assumption is made that the observed data are random with a Gaussian distribution. This leads to the well known %2 test for goodness of fit and to £2-norm solution. The solution to the inverse problem is the model vector m. It is not unique, but with some constraints we can obtain a reasonable estimate of m. Each of the model elements can be viewed as a random variable, and thus each estimate of the model is just a realization of a random process. Given an estimate of the probability density function (pdf) of the model m, we can, by calculating expected values, obtain an estimate of the model itself. From the above discussion, we can use the following equation to describe the inverse problem, g(m) = d = s + n, (4.9) Chapter 4. Absorption Compensation Theory and Practice 50 where observed data d, noise n and model parameters m can all be considered as realizations of random processes; ^(m) could represent a linear or non-linear function. Before we discuss the de-absorption problem using a statistical inverse method, let us first take a look at some related basic concepts in probability theory. The first is the probability density function (pdf) p(x) that plays a central role. For a random Gaussian process, P(x) = ^ e - ^ ^ 2 . (4.10) In this expression, a is the standard deviation, and a is the mean or mathematical expectation of random process x. Besides the Gaussian distribution, another distribution function that will be used in the de-absorption schemes is the Cauchy distribution: The second basic concept to review is conditional probability. If there are two related random processes, then the conditional probability p(d|m) refers to the probability of a realization d given that m occurs. The third concept is Bayes' Theorem. It states that rfm|d) = (4.12) Statistical inversion theory is based on Bayes' Theorem. Because in an inverse problem, we always have observed data d, p(m|d) is a posterior probability. A common approach in Bayesian inversion aims at maximizing p(m|d). After data are observed, p(d|m) is named the hkelihood and is measure of the possibility of the m that create d. An approach to estimate m from d, in the situation when no additional information about m is available, is to maximize the likelihood p(d|m), giving rise to the MAP (maximum a posteriori) solution. Chapter 4. Absorption Compensation Theory and Practice 51 By assigning a prior probability distribution p(m), the conditional distribution func-tion p(d|m) and maximizing p(m|d), an objective function for the inverse problem can be constructed based on Bayesian theory. 4.2.3 Absorption Compensation and Inversion Throughout this section, we will solve the de-absorption problem using statistical inver-sion theory which is based on Bayes' theorem. We will consider data contaminated by noise that is normally distributed as iV(0,(T2), where the subscript n represents noise. The conditional distribution of the data is given by (from equation (4.10)) / 1 x ( M - l ) / 2 j>(d|m,<rn) = ( ^ - ) e-(i/2. |) | |d-Gm||| > ( 4 1 3 ) where M is the length of the data vector. For a linear system d = G m , where G is called the coefficient or kernel matrix. Let p(m,<7m) indicate a prior distrib-ution for m conditional on a parameter crm. From Bayes' theorem (equation (4.12)), we have: ( \A \ p(djm,o-n)p(m,o-m) p(m|d,<7m,<rn) = — . (4.14) p(d) Here the model parameter m is the reflectivity function, so that we can use sparseness constraint in inversion scheme. 4.2.4 Cauchy-Gauss Mode l The following assumes that a seismic wavelet is available; it is either obtained from check-shot survey, or extracted from a seismic trace. It is believed that the model parameter (reflectivity m) is sparse. We have found that a sparse distribution of m provides not Chapter 4. Absorption Compensation Theory and Practice 52 only a manner to regularize the inversion algorithm, but also to increase the resolution of the inversion results (Sacchi and Ulrych, 1996). The reflectivity function m is a vector, and each element m*. of m is an independent random process. Assuming that all m*. have the same standard deviation <rm, the joint probability density function of all m-k p ( m , c r m ) = Y[p{mk,am). k If p(mk,crm) satisfies a Cauchy distribution (4.11), M-l / 2 \ _ 1 P(m,am)= 1 + r - H • (4-15) k=0 \ / < 7 m / Inserting this Cauchy prior (equation (4.15)) and the data likelihood (equation (4.13)) into equation (4.14) and taking logarithm of the both sides of the equation, we have ln[p(m|d, am,an)} = -c(m) - - ^ ( d - Gm) T(d - Gm), (4.16) where c(m) is a constraint imposed by Cauchy distribution: which is a measure of the sparseness of the model parameter. Furthermore, denoting < c^s(m) = — ln[p(m|d, <rm, crn)\, we can observe that maximiz-ing p(m|d, crm, crn) is equivalent to minimizing 4>cg(m). Therefore, the cost function for the Cauchy-Gauss model is set to be & g(m) = c(m) + r ^ ( d - Gm) T(d - Gm), (4.17) where m is the reflectivity, d is recorded seismic signal in the time domain, G which is composed of time-variant wavelets (bT[t — r)). Both t and r are in the range of a trace Chapter 4. Absorption Compensation Theory and Practice 53 length and b0(t - 0) (4.18) bM{t - M) These time-variant wavelets are calculated from equation (A.7), assuming that the ori-ginal source signature and Q values are known. The time-variant wavelets in the time domain have extremely complex analytical form (Verala et al., 1993). Thus numerical values are computed from their corresponding expressions in the frequency domain. In this way, the inverse Q filtering is formulated as an inverse problem. The model resulting in minimum cost function is the reflectivity function we want to find. The objective of this inverse problem is the same as time-variant deconvolution. The solution of this inverse scheme is stable. This Cauchy-Gauss model has also been used in acoustic impedance inversion (Sacchi and Ulrych, 1997), signal interpolation and interpolation, (Sacchi, Ulrych and Walker, 1998) and high resolution Radon transform (Trad and Ulrych, 1999). 4.2.5 Implementation In this section, we will discuss the implementation of the inverse problem using Cauchy-Gaussian model. The objective function will be minimized is (j)cg(m) in equation (4.17). Taking the derivative of <^ C3(m) with respect to m (Sacchi et al., 1998), we have: d<j>cg(m) _ dc(m) + 1 r _ 3m dm 2<rl (4.19) anc dc(m) dm Chapter 4. Absorption Compensation Theory and Practice 54 where S is a M x M diagonal matrix with elements Skk = l + ^fN A: = 0,1, • • • , M — 1 and depends on m. Equating the derivative in equation (4.19) to zero yields m = (AS - 1 + G T G ) _ 1 G T d , (4.20) where A = This equation is non-linear and must be solved by means of an iterative procedure. To construct an iterative algorithm, equation (4.20) is written as : m = SG r (AI i V + G S G r ) - 1 d = S G T p , (4.21) where the auxiliary vector p is a vector which is obtained from the solution of the system (AIJV + G S G T ) p = d. (4.22) The iterative algorithm starts with setting the observed seismic trace as the initial model parameter m^0). This initial solution is also used to generate the matrix In each iteration, we first compute p ( i t e r - l ) = [ A j ^ + G S ^ - ^ G 7 ] " V 4 6 ' - 1 and then update the model parameter m ( t t e r ) _ g ( t ter - l)Qr ( i t e r - l ) _ The iteration stops when the convergence of solution is reached. To solve the inverse problem iteratively requires a lot of computations, so it is neces-sary to explore the structure of the matrix to find the fastest possible algorithm. In the equation (Al^ y + G S G r ) b = d, the coefficient matrix is Hermitian Toeplitz. Thus a fast solver like the Levinson recursion can be used in the inversion. Chapter 4. Absorption Compensation Theory and Practice 55 4.3 Examples In the following, we show some examples of using this inverse method to perform absorp-tion compensation. The scheme is tested on both synthetic and real data. The first example is shown in Figure 4.4 which is used to examine whether the Cauchy-Gauss model can be used to extract absorption responses by assuming a known wavelet. In the inverse implementation, the kernel matrix (equation (4.18)) is constituted by one time-invariant wavelet. In this figure, the lowest panel (e) is the input synthetic signal without noise. The wavelet is known, and the absorption is unknown. Using the sharpness constraint, the inverse result is almost the same as the original absorption responses after only 2 iterations. When noise is added to the synthetic trace, the extracted absorption responses are much different, the test result is not demonstrated here. The second example shows that when both the wavelet and quality factor are known, the inverse scheme can yield a sharp reflectivity function in different noise environments. Figure 4.4a contains five traces. They are the same except that different noise is added to each trace; the noise level is 20 percent. Figure 4.4b is the inverse result. The rightmost trace is the original reflectivity series plotted as a reference. It can be observed that the reflectivity function is recovered very well with respect to the correct location and value. This shows that Cauchy-Gauss model can help to produce a sharp reflectivity and, using this inverse method, de-absorption may be implemented. Figure 4.6 shows the convergence of the inversion operation. Similar to Figure 4.4, panel (e) is an input trace from Figure 4.5a. Panels (a), (b), (c) and d) are inverse results after 1, 2, 3 and 4 iterations , respectively. The inverse result becomes sharper as the number of iterations increases. As the noise level reaches 25 percent, the inverse result is still acceptable which is shown in Figure 4.7. The locations of extracted reflectivity coefficients are correct; only Figure 4.4: Wavelet removal from a synthetic signal using Cauchy-Gauss model, (e) The original input trace; after 1 iteration the sharper responses is in (a); (b), (c) and (d) are inverse results after iteration 2, 3 and 4 respectively, they are almost the same. Chapter 4. Absorption Compensation Theory and Practice 57 Figure 4.5: (a) A synthetic signal with 20 percent noise; (b) The result of absorption compensation using an inverse method, the rightmost trace is plotted as a reference. Chapter 4. Absorption Compensation Theory and Practice 58 Figure 4.6: Convergence process of the inverse scheme; (e) is a trace from 4.5a. (a), (b), (c) and (d) are inverse results after iteration 1, 2, 3 and 4 respectively. Chapter 4. Absorption Compensation Theory and Practice 59 a) b) 0 2 4 6 0 2 4 6 8 T r a c e Index T r a c e Index Figure 4.7: (a) Synthetic signals with 25 percent noise; (b) The result of absorption compensation using an inverse method. some values are slightly different from the reference. The third example is a real data example. Figure 4.8a, is a stacked seismic section. The structures are not complex, consisting of some horizontal reflections. De-absorption processing can be used to improve the resolution of this section. First, quality factors are estimated from a trace in the section using WTVS method (Chapter 2). Second, using these extracted Q's to perform phase correction, a zero phase wavelet is extracted from the auto-correlation of a trace (only the shallow part of the trace is used). As we have both wavelet and attenuation characters, we can compute the kernel matrix G. The inverse computation begins by setting a value for the ratio of , the inverse results is Chapter 4. Absorption Compensation Theory and Practice 60 updated iteratively. Following the inversion computation, the wavelet is convolved with the extracted reflectivity functions. The final results are shown in Figure 4.8b. A number of reflections, which are originally undifferentiable in Figure 4.8a, can be observed around time t = 920 ms after de-absorption, the lateral continuity can be tracked from trace to trace. The improvement in resolution is clear, a result which could be helpful in the interpretation of this seismic section. Chapter 4. Absorption Compensation Theory and Practice 61 Figure 4.8: (a) A real seismic section; (b) The output of absorption compensation Chapter 5 Estimation of Quality Factors from CMP Gathers It is important to understand the attenuation characteristics of the subsurface because they can help to improve the resolution of seismic data, decrease the possibility of in-correct AVO analysis, and give direct information of lithology. Estimates of the quality factor, Q, are commonly obtained from vertical seismic data or stacked surface seismic data discussed in Chapter 3 and references therein. In this chapter a method is de-scribed that allows Q to be estimated directly from CMP gathers. Attenuation of the wavefield is dependent on three parameters: frequency, traveltime in the medium and medium Q. Assuming that the amplitude spectrum of the seismic source signature may be modeled by a Ricker wavelet spectrum, we derive an analytical relation between Q and seismic data peak frequency variation both along offset and vertical time direction. Q is estimated from CMP gathers in a layer stripping mode. 5.1 Basic Assumptions Despite an early awareness of the effect of intrinsic damping on seismic data, studies of the estimation of Q from multi-fold seismic reflection data remain an important area of research. Conventional methods (e.g. the spectral ratio method which was discussed in Chapter 3) uses amplitude information to estimate quality factor. The problem is that, due to widespread noise and other factors, such as geometrical spreading and scattering effect, such information is not always robust. 62 Chapter 5. Estimation of Quality Factors from CMP Gathers 63 In the following, I will introduce an analytical formula to estimate Q from prestack data, specifically, CMP gathers. The derivation of the formula is based on the two basic assumptions: 1. A CMP section represents multiple observations of an underground structure. It provides information in the time domain and in the offset domain, allowing for the extraction of information concerning structure, lithology and material properties such as velocity and Q. 2. Reflection arrival times are determined by interval velocities and the geometrical structure of the substructure. Attenuation of the received signals is only determined by the interval Q's and traveltimes in each layer. If the amplitude spectrum of a seismic wavelet is assumed to be Ricker-like, interval Q's can be computed solely from the the variation of the peak spectrum frequency as a function of time. We have developed an equation that relates attenuation to spectral peak frequency variation. Using this relationship, we determine interval Q's by means of a layer stripping approach. Tests show that this method allows for the determination of interval Q's with reasonable accuracy. 5.2 Absorption and Peak Frequency Variation In seismic data processing, recorded seismic data are commonly viewed to be the con-volution of a seismic source signature with a reflectivity series. The seismic source is generally unknown, although in some cases, it can either be measured or assumed to be of minimum phase. The effects of anelasticity may be incorporated into the model by means of convolution with an earth filter. This filter is causal, minimum phase, and depends on Q(t). In a routine seismic data processing procedure, inverse Q filtering is often used to remove this attenuation effect (Varela et al., 1993). In this chapter, instead Chapter 5. Estimation of Quality Factors from CMP Gathers 64 of studying details of the attenuation filter, we concentrate on the effect of Q on the wavelet peak frequencies, and construct a relationship between Q and peak frequency translations. First, we consider the amplitude spectrum of the source wavelet to be well represented by that of a Ricker wavelet (Ricker, 1977). The frequency spectrum of the Ricker wavelet is expressed in the frequency domain by equation (5.1): *(/)=4^e_/2//™ (s-1) where fm is called dominant frequeny. In the following, the frequency associated with the maximum of the amplitude spectrum is termed " peak frequency" for convenience. The evolution of the amplitude spectrum through time is now modeled as that of a Ricker wavelet traveling in a viscoelastic medium. Other factors such as geometrical spreading are not considered becausethey do not change the shape of the amplitude spectrum. After traveling a distance corresponding to a time £, the amplitude spectrum becomes mt)=(^)(i)e-"l&mt) (5-2) where H(f, t) is the absorption filter. In the frequency domain, this filter can be expressed by the equation H(f) = eXp(-f a(f,l)dl), (5.3) J ray where the integral is evaluated along the ray-path where Q(l) and c(/) are the quality factor and velocity, respectively, defined at each point of the ray path, and a(f, I) = Q(l)c(l)-Here, as is common practice, we have assumed that Q(l) is independent of frequency. This assumption is justified by many experiments (Ricker, 1977, Kjartasson, 1979, White, 1983). Chapter 5. Estimation of Quality Factors from CMP Gathers 65 5.2.1 One layer case Considering the propagation of a wave traveling in a half space with quality factor Q for t seconds, we determine the spectrum of the received signal as (from equation (A. 11)) B(f,t) = B{f)e-s? (5.4) We can see from this expression that, as time increases, the absorption increases with frequency and results in the translation of the peak frequency. This phenomenon is clearly illustrated in Figure 5.1. Due to absorption, the width of the source wavelet increases and consequently the amplitude spectrum becomes narrower. If the travel time is known, we can obtain Q from the spectral variations. A common method used to accomplish this task is the spectral ratio method which has been discussed in Chapter 3. In the spectral ratio method, the determination of Q is complicated by overlapping wavelets that leads to amplitude spectra that do not reflect the wavelet spectrum. The amplitude of seismic waves is affected by many factors, such as geological struc-ture, geometrical spreading, and AGC (automatic gain control), which is applied to balance the trace amplitudes. However, it is almost exclusively the quality factor that affects the shape of the wavelet amplitude spectrum. Based on this idea, we have devel-oped here a method to estimate Q from the spectral variation of reflections in a CMP gather. Because it contains both time and offset information, a CMP gather is very well suited to determine Q. Including all Q unrelated functions into an amplitude term, we write the spectrum as B(f1t) = A{t)B(f)e-Sir (5.5) where A(t) is an amplitude factor independent of frequency and absorption. The peak frequency, fp, can be determined by equating the derivative of the spectrum to zero ?*m = A ( 0 »fLQ. -¥ + 4, ) * ( / ) . - ? ( - £ ) = 0 (5.6) Chapter 5. Estimation of Quality Factors from CMP Gathers 66 Figure 5.1: (a) An event in a CMP gather is generated by a reflector in an absorptive medium, where amplitude of each trace has been normalized, (b) Amplitude spectra of original source signature, trace 1, trace 11, and trace 21 which are indicated by [0], [1], [11] and [21] respectively. Chapter 5. Estimation of Quality Factors from CMP Gathers 67 From equation (5.1), determining dB(f) (2f e~f2/fm _f_ -2f 9f v^v/Ar ' v^\fij~ \pm/ and inserting this expression into equation (5.5) we obtain the peak frequency (5.7) fp ~ frt \ (£)+ ii TVt (5.8) The relationship between Q and the shift of peak frequency now follows Q = •Ktfpfl (5.9) 2 ( / » - PP) This equation shows that, if the original wavelet spectrum is known, Q may be computed from the CMP gather using only one offset. In practice, of course, we do not know the initial wavelet dominant frequency fm. Let us assume, however, that the amplitude spectrum of the initial source wavelet may be approximated by that of a Ricker wavelet. Designating the peak frequencies at two different times, t\ and t2, by fpl and fp2 respectively Q = 2 ( / £ - / P 2 i ) 2 ( / £ - / £ ) from which one can determine that (5.10) frr fplfp2(t2fpi — t i f p 2 ) (5.11) \ | t2fp2 — ii/pi Equation (5.11) tells us that, if the peak frequencies at two different reflection times can be estimated, the original dominant frequency of the source wavelet can be computed. Equations (5.10) and (5.11) allow us to obtain an average Q thereby removing the effects of surface fluctuations and random noise and consequently, improve the accuracy of Q. They can also be used in estimating quality factors from a vertical-incident trace as in Chapter 3. Chapter 5. Estimation of Quality Factors from CMP Gathers 68 Although we have assumed a Ricker amplitude spectrum in our computations, due to the similarity with observed amplitude spectra, many other models are possible. One example is a Gaussian spectrum i?(/) = e x p ( - ( / ^ / m ) 2 ) (5.12) Following the procedure outlined above, the relationship between Q and the variation of peak frequencies is Q = T^T (5.13) Jm Jp In practice, the amplitude spectrum of a seismic wavelet is not exactly like that of a Ricker wavelet or a Gaussian function. In many situations, however, it can still be approximated by a Ricker spectrum. 5.2.2 Multi-layer case We consider, first, the case of two layers, with quality factors Q\ and Q2 respectively. The travel times of a wave traveling through layer 1 and layer 2 are ti and t2, respectively. Applying the absorption model of equation (5.6) W l *7'2 B(f,t) = A(t)B{f)eSTe-*r (5.14) where t — ti -\- t2. Using the peak frequency fp associated with B(f,t), Q2 can be estimated with knowledge of Q\ and the initial peak frequency fm of the source wavelet nt2Q 1^2 = -a where Qi - 7 r i i VI - Vl a = fpfm We now employ the concept of an average Q. In other words, we let (5.16) e oi e <?2 = e Q (5-17) Chapter 5. Estimation of Quality Factors from CMP Gathers 69 which allows an expression to be derived that relates the interval Q2 to the average Q t2QiQ Qi (5.18) (ti + t 2 )Qi -hQ Since the dominant frequency fm of the initial wavelet and Q\ have already been deter-mined from upper layer arrivals and t\ and t2 can be estimated, in theory, Q2 can be computed using equation (5.15). Given the complexity of the subsurface, in order to use all offset information avail-able in a CMP, some approximations have to be introduced. Considering a multi-layer medium, we write the absorption equation (equation (5.5)) as B(f,t) = A(t)B(f)exp ' N E ,i=l Qi (5.19) where Qi and Ai; are the quality factor and the traveltime in layer i, respectively. Follow-ing the approach used in velocity estimation, we assume straight ray paths and compute the total travel time of a reflection at a particular offset as N ijV = E Ati-i=l We can now write tN to(N) [*o(0 - <o(i-i)] (5.20) (5.21) where io(AT) is the zero-offset travel time of reflection N, and io(i) is the zero-offset travel time of the upper reflections. Defining the Q of layer N to be QN, we can split the attenuation operator in equation (5.21) into two factors as follows B{f,t) = A(t)B(f)exT> N-l E ,t=i -TTfAti Qi exp - 7 T / A i N QN (5.22) which allows us to obtain the following equation for QN 7T A f j v QN — a-B (5.23) Chapter 5. Estimation of Quality Factors from CMP Gathers 70 where a is the same as in equation (5.15), and Proceeding by means of layer-striping, Q values can be calculated layer-by-layer. Because a straight ray path approximation is used, the computed Q is not the actual interval quality factor. Analogous to RMS velocities, such Q values are called RMS Q values. In a manner analogous to that used to determine interval velocities, we can determine interval Q's from the RMS values using the Dix formula. ^interval ~ . _ . K0^) where Q3{ is the average of the calculated Q values at different offsets for layer i, and tSji is the travel time in layer i when ray path is assumed to be straight across interfaces. 5 .3 Examples To estimate Q from CMP gathers, we assume that the arrival times for the main reflec-tion events are known. Fourier transforms are computed in the window containing the reflection at each offset, then each amplitude spectrum is fitted with a Ricker spectrum, and the peak frequency of the spectrum is estimated. Using equations (5.11), (5.10), and (5.23), we can now extract the Q's layer by layer from peak frequency variation. The results of a very simple test are shown in figure 5.2, where a synthetic CMP gather with two events is illustrated (Figure 5.2a). The original Ricker wavelet has a dominant frequency of 60 Hz. Absorption is modeled by using low Q values to emphasize the effect: the Q values in the two layers are 10 and 20. The actual and estimated Q curves are shown in Figures 5.2b and 2c respectively. The corresponding dominant frequency has been estimated as 60.67 Hz. The agreement between the corresponding values shows that the method works well for ideal synthetic data. The variation of peak Chapter 5. Estimation of Quality Factors from CMP Gathers 71 T r a c e Index Input Q va lue Est imated Q va lue Figure 5.2: (a) A synthetic CMP gather of two reflections with absorption, the amplitude has been normalized, (b) Real Q's in layer 1 and layer 2, it is unknown in layer 3. (c) computed Q values. Chapter 5. Estimation of Quality Factors from CMP Gathers 72 2 0 1 Q_ 5 0 5 1 0 1 5 2 0 2 5 Offset Index Figure 5.3: Variation of peak frequencies of two reflections. Solid line is that of the first layer. Dashed line is that of the second layer. frequencies with offset is shown in Figure 5.3. The remarkable decrease is, of course, due to absorption. This example demonstrates very clearly why absorption must be compensated for if resolution of prestack data is an issue. Following the determination of Q we are faced with the task of inverse Q filtering. However, details of the inverse Q filters are not discussed in this chapter (see Chapter 4). An example showing the effect of inverse Q filtering on a prestack CMP gather using obtained Q values is illustrated in Figure 5.4. A synthetic CMP gather composed of three reflections is modeled in Figure 5.4a, where signal amplitudes remain unbalanced. The first two reflections merge at far offset and we consider them as one to estimate the peak frequency variation. The extracted quality factors are shown in figure 5.4b at relative magnitude. Using these Q values, inverse Q filtering is performed. The Q compensated result is shown in figure 5.4c, demonstrating that the reflections associated with the two Chapter 5. Estimation of Quality Factors from CMP Gathers 73 Figure 5.4: Experiment to show the effect of inverse Q filtering, (a) A synthetic CMP gather with 3 reflections; (b) Extracted Q values, only two layers have been separated; (c) The result of absorption compensation. upper layers are separated. These two experiments are employed to simply show the validity of our method for quality factor extraction and the resolution improvement of prestack seismic records using accurate Q values. Although our method is illustrated using CMP geometry, it can also work on common shot gathers (CSG), if the geometry difference can be neglected by considering the substructure as horizontally layered medium. As discussed in many publications (e.g. Futterman, 1962), wave propagation is a causal process, and leads to dispersion in an attenuating medium. In practice, inverse Q filtering should account for both absorption compensation and dispersion compensation at the same time. The next example is to test our extracted Q's in real data processing. Chapter 5. Estimation of Quality Factors from CMP Gathers 74 A CSG gather is shown in figure 5.5a which was obtained from shallow seismic survey, absorption phenomena can be easily observed from the broadened wave shape along time axis. The inverse Q filtering result is shown in Figure 5. 5b, and has been obtained using the Q compensation filter based on the method proposed by Verala et al (1993), because it is fast. Before attenuation compensation was performed, the data were processed by low-pass frequency filtering because inverse Q filter is commonly unstable for high frequency noise. 5.4 Observations and Discussion As illustrated in Figure 5.1, the absorption effects as a function of offset may be mistak-enly interpreted as AVO phenomena. Without Q compensation, the reflection amplitudes may appear to decrease with offset, as illustrated in Figure 5.4. It is clear that compens-ating for absorption in prestack data can be vital in AVO analysis. Figure 5.4 clearly shows that signals attenuate as a function of both traveltime and offset. Conventional deconvolution techniques used to improve resolution must rely on an adaptive formulation of wavelet estimation. The pitfalls of wavelet estimation, even if the wavelet is not time varying, are well known. Absorption compensation, based on a model of the underlying physics of the absorptive process, allows the computation of a section with increased resolution and a stationary source wavelet that can be subsequently used, if desired, for additional resolution by means of deconvolution. Prestack inverse Q filtering is well suited to compensate for absorption. Performed on poststack data, inverse Q filtering is similar to wavelet shaping, because stacking distorts the frequency information of the original seismic data. It is advisable to extract quality factors from prestack seismic data that have not been harshly processed by frequency filtering. Chapter 5. Estimation of Quality Factors from CMP Gathers Figure 5.5: (a) A common shot gather; (b) The common shot gather after inverse filtering. The Q curve has not been shown here. Chapter 5. Estimation of Quality Factors from CMP Gathers 76 Prestack processing provides another dimension, offset, that increases the informa-tion concerning signal and noise separation. Inverse Q compensation, performed in the prestack domain, can provide a balanced signal spectrum near to far offset and over time and will improve the vertical and lateral resolution of the final image. Underground lithology is characterized by its velocity, density and quality factor. The last two do not effect the arrival times of reflections, while the quality factor affects amplitude and the signal's frequency constant. Extracting velocity information from CMP gathers is a common practice. Here we devised similarly a formula to estimate the absorption character from the variation of spectrum with time and offset. A relation that allows the computation of quality factors from observed peak frequency variations was derived in this chapter, based on the reasonable assumption of a Ricker-like amplitude spectrum of the source signature. Chapter 6 Migration and De-Absorption As discussed in Chapter 4 and Chapter 5, absorption compensation and migration are done separately in seismic data processing. Absorption compensation is either performed before migration or after migration in a trace by trace mode. Actual ray-paths, and the nature of lateral continuity of the earth structures are neglected. In addition, because a stack distorts the original absorption characteristics, poststack de-absorption processing can only be considered as a wavelet shaping process to make seismic wavelets uniform along a trace. Considering that seismic waves traveling in the earth satisfy the viscoelastic wave equation, it is necessary to incorporate de-absorption and migration together. Ideally, this process could be done in prestack data; i.e., de-absorption and prestack migration are expected to be implemented in a single process step. The advantage of this process is that it aims at preserving both the amplitude and waveform of seismic signals during imaging. We can call this kind process high fidelity migration. In this chapter, I will review the recent developments of imaging techniques which take into account absorption and discuss the possibility of using a least squares migration (LSM) method to implement high fidelity migration. Prestack seismic data can be sorted into common offset sections. In this chapter, only the zero offset situation will be considered because of its mathematical simplicity. Once we have understood the principle for the zero-offset case, we can extend the technique to non-offset sections without much difficulty. 77 Chapter 6. Migration and De-Absorption 78 6.1 Migration Migration is a procedure in seismic data processing which aims at imaging the subsur-faces. In recent years, because of the developments in computer capacity, 3-D exploration and the requirement of finding small geological bodies, prestack depth and time migration is widely used for imaging more complex geological bodies. The research into lithology, lithofaces, crack density, and fluid content of potential reservoirs in hydrocarbon explo-ration relies heavily on the images obtained from migration. Both the horizontal and the vertical resolution of present migration techniques is expected to be improved. The essen-tials of a high fidelity migration technique are to preserve both amplitude and waveform of seismic signals during imaging (Ma, 1997). Migration methods include mainly ray-tracing-based KirchhofF migration methods, finite-difference methods, phase-shift migration methods and the method developed by Stolt (1978). The key points of these methods are their computation speed, stability and the ability to accommodate spatial velocity fluctuations. Few of these migration approaches take into account the intrinsic absorption in the wave propagation process. However, in recent years, there are increased interests of including de-absorption in seis-mic data migration and inversion. Mittet et al. (1995) proposed a prestack depth mi-gration method for accommodating absorption and dispersion effects in prestack depth migration schemes. Their method was presented in the frequency-wavenumber domain using a linear viscoacoustic medium. An extrapolation operator that compensates for absorption and dispersion was designed using an optimization algorithm which was then employed to devise a finite-difference migration scheme. The images obtained are su-perior to those obtained using an extrapolator without compensation for absorption. However, this method is unstable in the presence of high frequency noise. Causse et al. (1999) presented a viscoacoustic full-waveform inversion scheme which compensates Chapter 6. Migration and De-Absorption 79 for the effects of phase distortion but not for the effects of amplitude distortion using viscoacoustic gradient-based inversion algorithms. Numerical tests were performed on noise free synthetic data. Practical applications of viscoelastic wave inversion have also been reported. Blanck and Symes (1995) applied a 2-D viscoacoustic finite difference linear inversion to a real data set, using iterative linearized inversion to update model parameters. Charara et al. (1996) applied viscoelastic inversion to an offset VSP data also using an iterative scheme. These works show the increased interest and significance of viscoacoustic migration and inversion. The key point of wave theory based absorption compensation method is how to devise an algorithm which is stable to high frequency random noise. This problem has also been addressed in Chapter 4. One dimensional absorption compensation can be solved using inverse theory and migration can be implemented using a least squares scheme. It is a promising direction to formulate the problem of migration with absorption compensation as a model fitting problem by least squares. 6.2 Migration, the Adjoint Operation of Modeling A standard seismic migration operator can be regarded as the adjoint of a seismic forward modeling operator (Claerbout, 1992). It is important to realize that the adjoint is not the exact inverse of forward modeling. 6.2.1 Basic Concepts of the Adjoint Operation A great many of the calculations in science and engineering are really matrix multi-plication in disguise. The adjoint operator back-projects information from data to the underlying model. In matrix manipulation, the adjoint becomes matrix transpose for Chapter 6. Migration and De-Absorption 80 real matrices and the Hermitian conjugate for complex matrices. Geophysical modeling calculations generally use linear operators that predict data from models. Our usual task is to find the inverse of the calculations; i.e., to find models from data. Logically, the adjoint is the first step and a part of all subsequent steps in this inversion process (Tarantola, 1987). However, in practice, the adjoint sometimes does a better job than the inverse. This is because the adjoint operator tolerates imperfections in the data and does not demand that the data provide full information (Claerbout, 1992). Consider a matrix A to be a forward operator which projects a model vector m to a data vector d. Conceptually, the adjoint operation is to obtain m a = L r d , where the definition of matrix transpose is simply aj- = a .^ In practice, however, we often encounter matrices far too large to fit in the memory of any computer. Sometimes, it is also not obvious how to formulate the process at hand as a matrix multiplication. In such cases, the forward operator A and its adjoint A r are commonly represented by different algorithms. The dot-product test enables us to ascertain whether the algorithm for the adjoint of an operator is precisely consistent with the operator itself. To perform the dot-product test, we first load the vectors m and d with random numbers and then, using the forward program for A , we compute the vector d = Am, (6.1) and using the adjoint algorithm for A T , we compute rh = A r d . (6.2) Since d T A m = (A T d) T m, (6.3) the two scalars d r d and m T m are also equal, if the computation of AT is indeed adjoint to the computation of A. Chapter 6. Migration and De-Absorption 81 6.2.2 Gazdag modeling and migration In wave propagation, there are several pairs of adjoint operations (Claerbout, 1992): Kirchhoff migration is the adjoint of diffraction modeling, reverse time migration is the adjoint of the time-extrapolation modeling and, Gazdag phase shift migration is the adjoint of phase shift modeling. This section will demonstrate that Gazdag phase shift modeling and migration are a pair of conjugate operations. We will see that when using Gazdag migration to image absorption attenuated data directly, artifacts will be produced. Gazdag phase shift modeling and migration (Gazdag, 1978) can be expressed by the differential equation, dU dz In the forward modeling application, U(z) = U(kx,z,uj) is the up-going wave, — +iaU = S(z). (6.4) a = ^ / v * - k l , (6.5) S(z) = S(kx, z) is the exploding-reflector source. We take the medium to be layered with constant layer velocity, so that a is constant in a layer, and we put the sources at the layer boundaries. Thus within a layer we have dU/dz + iaU = 0 which has the solution U(zk) = U(zk + Az)eiaAz (6.6) For convenience, let Zk = e t a A z , so the delay of upward propagation is expressed by U(zk) = ZkU(zk + Az) or Uk = ZkU{zk+1). Besides crossing layers, we must cross layer boundaries where the (reflection) sources add to form the up-going wave. Thus we have f / f c _i =ZkUh + Sk.1. (6.7) Recursive use of equation (6.7) across a medium of three layers is expressed in matrix Chapter 6. Migration and De-Absorption 82 form as 1 - Z 0 0 0 1 -Zx 0 0 1 0 0 0 0 0 -z2 1 U0 So ' 5 3 _ (6.8) A recursive solution begins at the bottom with U3 = S3 and propagates upward . The conjugate of the delay operator Z is the time advance operator Z, so that, the conjugate operation to that in equation (6.8) is given by 1 0 -Zo 1 0 -Zi 0 0 0 0 1 0 So Si S2 S3 U0 U3 (6.9) 0 0 - Z 2 1 where S(z) (summed over frequency) is the derived image. Equation (6.8) and (6.9) can be written in matrix form A U = S and A r S = U, and AT is the conjugate transpose of A. Written in explicit form as forward and inverse problems respectively, U = A-'S = A S and ( A T ) _ 1 U = ( A ^ f U = A r U . (6.10) (6.11) The conjugacy can be easily observed from equation (6.10) and (6.11). Let us now examine whether the adjoint operator is the inverse operator. We form A A = 1 Zo ZQZ\ ZQZ\Z2 Zo 2 2 Z\ 2 Z2 Z2 ZQZ\ 1Z\ 3 3^2 ZQZ\Z2 1Z\Z2 2>Z2 4 (6.12) Chapter 6. Migration and De-Absorption 83 which it is not an identity and, therefore S fi S. (6.13) This means that the transpose of matrix A is not its inverse. Standard migration using the adjoint of the forward modeling operator will produce artifacts in migrated images. For absorptive medium, the extrapolator Zk has a form of (equation (A.8)) Zjk — e = exp i, or 2 , / w \ .1 1 - In I — ) +i — irQ \U>QJ Q kl Az. (6.14) Using the 15 degree approximation (Claerbout, 1987), i.e. ^Jk2 — k2 = k — | | and omitting the ^  related terms, we have a =i — TVQ \UJ0J 1 + i 2Q }" TTQ ^ n (wo)] (6.15) For simplicity, denote £ = 1 — ^-ln and separate a into real and imaginary parts (6.16) C > v0kl . a = —— + i VQ IQU) (ui v0kl _2Qv0 2QCu> Therefore, it is seen that the conjugate of Zk is approximately Zk =• exp \ A z ( ^ - V - ^ exp -Az UJ v0kl \2Qv~0 + AQ(UJ/ (6.17) which is a stable operator. This property is very important in constructing a stable de-absorption scheme. In contrast, the inverse, Z^1, is unstable. 6 . 3 Least Squares Migration (LSM) Accurate images are the desired end product of seismic migration. In many situations, the adjoint operation of the forward modeling operator is not a good enough approximation to handle data with absorption attenuation, or with finite aperture. In some cases, the Chapter 6. Migration and De-Absorption 84 exact inverse can be approximated by a generalized inverse which is formulated as a least-squares migration operator, rather than migrating data with the adjoint of the forward modeling operator. In this section, I will introduce the principle of LSM and probe the possibility of using LSM to migrate seismic data and compensate for absorption at the same time. Bamberger et al. (1982), Lailly (1983) and Tarantola (1987) provided the basic form-alism for least squares inversion. These works opened the way to new interpretations of wave equation migration, in particular in regard to iterative migration algorithms, which could be understood as minimization algorithms (Lebras and Clayton, 1988). Recently, least squares migration has received wide interest. Schuster (1993) devel-oped a technique of least-squares cross well migration. Chavent, G. and Plessix, R.(1999) presented a technique of true amplitude least squares depth migration. Nemeth et al. (1999) used the LSM method to migrate incomplete reflection data. These authors demonstrate that LSM can produce more accurate images than standard migration. The standard migration algorithm seeks to image the Earth's structures by applying the adjoint of the forward operator to the reflection data, i.e. the migrated section is given by the adjoint of the forward modeling operator. This, however, is an unneces-sarily restricted definition that can promote artifacts in the migrated section. A more general definition is to recast the standard migration problem as an optimization problem (Schuster, 1993). For an exploding reflector model, the forward modeling operation can be mathematically represented by d = Gm (6.18) where d and m are two vectors representing the predicted seismic data and reflectivity model respectively, and G represents the forward modeling operator. Chapter 6. Migration and De-Absorption 85 Standard migration evaluates m = G r d 0 (6.19) as the migrated image, where d 0 are the observed reflection data. Least squares migration seeks to find m that minimizes a cost function (j>(m) =|| Gm — d 0 ||2 + constraints. This provide a general framework to develop a variety of new migration scheme which might reduce artifacts inherent in standard migration sections. An example is the constrained Conjugate Gradient migration, the benefit of which might be a more accurately migrated image at the expense of increased CPU time. To minimize artifacts, we shall solve for m by minimizing the following objective function: 0(m)= || G m - d n ||2+ || Cm ||2 (6.20) where || • || 2 represents the £ 2 norm, || Cm) || is a regularization function, and C represents a constraint matrix. Expression (6.20) assumes that the C constrained model and the data residual probability distributions are Gaussian. Regularization reduces the uncertainty in the model space by imposing constraints on the model, and hence, assigns more information to the problem. The model that minimizes the cost function is represented by the least squares solution m = [G T G + C r C ] _ 1 G r d 0 . (6.21) If C T C is an identity matrix with large diagonal values then m in equation (6.20) becomes m ss G r d 0 . (6.22) This model is commonly interpreted as the standard migration image. The inverse of the matrix in equation (6.21) may be either too expensive to compute directly or impossible to compute because an operator may not have an explicit matrix Chapter 6. Migration and De-Absorption 86 form. In such cases, we can resort to an iterative technique such as a conjugate gradient method. Minimizing the cost function <^ >(m) in equation (6.20) is equivalent to fitting the following two simultaneous equations G m PS d 0 C m PS 0 Combining these two equations together using block matrices, we have G d 0 m = C 0 (6.23) or, in a concise form where A = G C , x = m and y Ax = y, d 0 0 (6.24) A computation template for the method of conjugate gradients (Claerbout (1992)) for the problem Ax = y is: Set x = A T y , r = Ax — y, 3 = 0, and s = 0 (a symbol in bold face is a vector ), for each iteration, repeat the following: { Ax = A T r if not the first iteration 3 = — ^ A*^ 7 =11 Ax || 2 s = Ax + 3s Ar = As a = - 7 / || Ar || 2 x = x + as r = r + a A r } Chapter 6. Migration and De-Absorption 87 The initial guess for the conjugate gradient method is the conventional migration, which is the image after applying the conjugate operator of the forward modeling. Since this initial guess is very close to the solution except for some artifacts and absorption, it is expected that a fast convergence will occur. Since in the conjugate gradient solution, both the forward and the adjoint operations are stable. Using a LSM scheme, we can perform migration and absorption compensation at the same time to obtain high resolution image. The soundness of this migration with de-absorption scheme needs to be verified by numerical tests. C h a p t e r 7 C o n c l u s i o n s a n d D i s c u s s i o n s The important role that absorption plays in seismic wave propagation has long been recognized by scientists. Absorption is a comprehensive effect caused by several factors. Internal friction along the contact between solid rock matrices and viscosity between the rock matrix and pore fluid are two such factors. There is no general consensus as to which of these factors is dominant. For various rock types and for a wide range of environmental conditions, data indicate that attenuation is approximately proportional to frequency. There are also several mathematical models to describe absorption, such as the Voigt model and Generalized linear solid model. A straightforward way to quantify the absorption phenomenon is to use the quality factor (Q) of a medium. Under the requirement of constant Q, linearity and causality, a velocity dispersion relation was de-picted in equation (2.12). The characteristics of absorption related amplitude attenuation and phase distortion are fully described by Q and the Q related velocity (see equation (A.7) and (A.8)). In this thesis, I emphasize the analogy between quality factors and velocities. They have analogous uses in wave propagation theory and they are estimated from data in similar ways. I have introduced a formula (equation (5.10)) to calculate quality factors using the peak frequency variations of a wave with time (for the one dimensional case, the peak frequency variations were measured along a trace, for the two dimensional case, the vari-ations were measured along a reflection ). The assumption in the theoretical derivation is that the source wavelet has an amplitude spectrum similar to that of an Ricker wavelet 88 Chapter 7. Conclusions and Discussions 89 (equation (5.1)). Although actual amplitude spectra may not, of course, conform to this assumption, in many situations, an amplitude spectrum can still be well approximated by a Ricker spectrum. The implication is that this assumption is not rigorously demanded. For vertically incident data (such as that found in poststack seismic section) ana-lysis, a windowed time-variant Fourier amplitude spectrum was used. The windowed time-variant spectrum (WTVS) analysis is superior to conventional Fourier and wavelet analysis in disclosing the peak frequency variations with time. In WTVS, controlling points can be picked in a way analogous to velocity analysis. The controlling points divide a trace in several layers and their coordinates are used to calculate the quality factors of each layer (Figure 3.5). In many respects, the analytical method is superior to the spectral ratio method or other conventional methods. The analytical method only uses the frequency variation information. Besides absorption, it is not affected by other amplitude factors such as reflection, transmission and geometrical spreading. It is robust to random noise. For the two dimensional case, it uses multi-fold information. For the one dimensional case, it benefits from the frequency trend analysis of WTVS. For estimating quality factors from one dimensional data, a nonlinear optimization scheme is derived (equation (3.25)). An experiment using synthetic data shows promising results. This approach does, however, require a non-absorption trace as a reference, which is often not available in practice and, therefore has restricted applications. If well data are available, it could be useful in correlating well data with close well seismic traces. Using this optimization scheme, detailed absorption characteristics could be extracted from seismic data. Extrapolating the quality factors from close-well traces along horizons of a seismic section, a quality factor section could be obtained as is done in well-constrained acoustic impedance inversion (Martines et al. 1992). This could be very useful for reservoir description. Chapter 7. Conclusions and Discussions 90 Absorption compensation or de-absorption is an important step in seismic data pro-cessing, it helps to improve the resolution of seismic data, and to balance the frequency contents of a signal with time. Using statistical inversion theory, I formulated the inverse Q filtering as an inversion (optimization) problem. The prior information imposed on the model parameters is the sparseness of reflectivity which is represented by a Cauchy model. The inverse problem is solved iteratively. Experiments performed on both synthetic and real data show that the inverse scheme can produce ideal results for absorption compens-ation. The advantage of the scheme is that it overcomes the obstacle of instability which is a natural drawback of common inverse Q filtering. To preserve both the amplitude and waveform of seismic signals during imaging, it is necessary to incorporate de-absorption and migration together (high fidelity imaging). Standard migration is the adjoint operation of forward modeling. When applied in data processing, it will produce artifacts. This thesis concludes with a description of least squares migration which is implemented by recasting the standard migration problem as an optimization problem and solving the optimization problem iteratively. Least-squares migration can produce more accurate images than standard migration. It would be useful to formulate the problem of migration with absorption compensation as a model fitting problem using a least squares migration. References Aki K. and Richards P. G, 1980, Quantitative Seismology-Theory and Methods, W.H. Freeman and Company, pl67-185. Bamberger, A., Chavent, G., Hemon, C. and Lailly, P., 1982, Inversion of normal incid-ence seismograms : Geophysics, 47, 757-770. Beylkin G. and Keiser J. M., A new class of time discretization schemes for the solution of Nonlinear PDEs: Documenta Mathematica, Extra Volume ICM 1998, III, 481-490. Berkhout A. J., 1985, Seismic Migration, Imaging of acoustic energy by wave field extra-polation, A: theoretical aspects, Elsevier Science Publishers. Biot, M. A., 1956a, Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low-frequency range:J. Acous. Soc. Am., Vol. 28, 168-178. Biot, M. A., 1956b, Theory of propagation of elastic waves in a fluid-saturated porous solid. II. higher-frequency range: J. Acoust. Soc. Am., 28 179-191. Blanch, J. O. and Symes, W. W., 1995, Efficient iterative viscoacoustic linearized inver-sion: Annual Meeting Abstracts, Society Of Exploration Geophysicists, 627-630. Bleistein, N., 1984, Mathematical methods for wave phenomena: Academic Press. Carcione, J. M., 1995, Constitutive model and wave equations for linear, viscoelastic, anisotropic media: Geophysics, 60, 537-548. Carcione, J. M., 1993, Seismic modeling in viscoelastic media: Geophysics, 58, 110-120. 91 References 92 Carcione, J. M., Kosloff, D. D. and Kosloff, R., 1988, Viscoacoustic wave propagation simulation in the Earth: Geophysics, 53, 769-777. Causse E. , Rune Mittet and Ursin B., 1999,Preconditioning of full-wave Inversion in visco-acoustic media: Geophysics, 64, no.l, 130-145. Charara, M., Barnes, C. and Tarantola, A., 1996, The state of affairs in inversion of seismic data: An OVSP example: Annual Meeting Abstracts, Society Of Exploration Geophysicists, 1999-2002. Chavent, G. and Plessix R., 1999, An optimal true amplitude least squares prestack depth migration operator: Geophysics, 64, 508-515. Cheng, Q. 1992, Mathematical Foundation of Digital Signal Processing, P718-728, Pet-roleum Industry Publishing Company in Chinese. Claerbout, J., 1976, Fundamentals of geophysical data processing : with applications to petroleum prospecting, New York McGraw-Hill, cl976. Claerbout, J., 1985, Imaging the Earth's Interior, Blackwell Scientific Publications, Bo-ston, p221-235 Claerbout, J., (1992) Earth Soundings Analysis Processing versus Inversion, Blackwell Scientific Publications, Boston, 257-313. Coruh, C , Saatcilar, R., Dimirbag, E. and Costain, J. K., 1990, Estimation of Q-factor for spectral whitening: Annual Meeting Abstracts, Society Of Exploration Geophysi-cists, 1082-1085. Dasgupta, R. and Clark, R. A., 1998, Estimation of Q from surface seismic reflections data; Geophysics, 63, 2120-2128. References 93 Dong, Z. and McMechan, G. A., 1993, 3-D prestack migration in anisotropic media: Geophysics, 5879-90. Du, S., 1996, Joint inversion and its application in seismic sequence analysis: Annual Meeting Abstracts, Society Of Exploration Geophysicists, 1991-1994. Du, S., 1996, Seismic dynamics: University of Petroleum Publishing Company, in Chi-nese. Emerich, H. and Korn, M., 1987, Incorporation of attenuation into time-domain compu-tations of seismic wave fields: Geophysics, 52, no. 09, 1252-1264. Futterman, W. I., 1962, Dispersive Body Waves: Journal of Geophysical Research, 6 7 , 5279-5291. Gabor D., Theory of Communication. J.IEE,93:429-457, 1946. Gazdag, J. , 1978, Wave equation migration with the phase-shift method: Geophysics, 43, 1342-1351. Griffiths, L. J., Smolka, F. R. and Trembly, L. D., 1977, Adaptive deconvolution - A new technique for processing time-varying seismic data: Geophysics, 42, no. 04, 742-759. Hale, D., 1991, Stable explicit depth extrapolation of seismic wavefields: Geophysics, 5 6 , 1770-1777. Hargreaves, N. D., 1992, Similarity and the inverse Q filter: Some simple algorithms for inverse Q filtering (short note): Geophysics,57, 944-947. Hargreves, N. D., and Calvert, A. J., 1991, Inverse Q filtering by Fourier transform, Geophysics, 5 6 , 519-527. References 94 Hargreaves N.D. and Calvert A.J., 1993, Inverse Q filtering by Fourier Transform: Geo-physics, 56, p519-527. Harris, P. E., Kerner, C. and White, R. E. 1997, Multichannel estimation of frequency-dependent Q from VSP data: Geophysical Prospecting, 45, 87-109. Irving, J. D., 2000, Estimation and correction of wavelet dispersion in ground penetrating radar data; Master thesis, the University of British Columbia. Jannsen, D., Voss, J. and Theilen, F., 1985, Comparison of methods to determine Q in shallow marine sediments from vertical reflection seismograms: Geophys. Prosp., 33,479-497. Johnston, D. H., Toksoz, M. N. and Timur, A., 1981, Attenuation of seismic waves in dry and saturated rocks: Geophysics, 44, 691-711. Kjartansson, E., 1979, Constant Q-wave Propagation and Attenuation: Journal of Geo-physical Research, 84, 4737-4748. Kerner, C. and Harris, P. E., 1994, Scattering attenuation in sediments modeled by ARMA processes-validation of simple Q models: Geophysics, 59, 1813-1826. King, M. S., Zimmerman, R. W. and Corwin, R. F., 1988, Seismic and electrical proper-ties of unconsolidated permafrost: Geophys. Prosp., 36 349-364. LeBras, R. and Clayton, R. W., 1988, An iterative inversion of back-scattered acoustic waves: Geophysics, 53, 501-508. Liang, G., et al, 1990, Wave Propagation in Viscoelastic Media and Inverse Q filtering: 60th Ann. Internat. Mtg., Soc. Expl. Geopys., expanded abstracts. P 1079-1081. References 95 Lines, L. R. and Treitel, S., 1984, Tutorial - A review of least-squares inversion and its application to geophysical problems : Geophys. Prosp., 32, 159-186. Liu, H., Anderson D. L. and kanamori H., 1976, Velocity Dispersion due to anelasticity; implications for seismology and mantle composition : Geophys. J. R. astr. Soc, 47, 42-58. Ma, Z., 1997, Computational geophysics (in Chinese): Tongji University Publishing Com-pany Pl-9. Mallat, S., 1998, a wavelet tour of signal processing, Academic Press, P67-85. Malvern L. E., 1969, Introduction to the mechanics of a continuous Media: Prentice-Hall Ltd. Toronto, New Jersey, P313-319. Margrave, G. F. and Ferguson, R. J., 1999, Wavefield extrapolation by non-stationary phase shift: Geophysics, 64, 1067-1078. Margrave, G. F., 1998, Theory of non-stationary linear filtering in the Fourier domain with application to time-variant filtering: Geophysics, 63, 244-259. Martinez, R. D., Cornish, B. E., Rebec, A. J. and Curtis, M. P., 1992, Complex reservoir characterization by multiparameter constrained inversion, in Sheriff, R. E., Ed., Reservoir geophysics: Society Of Exploration Geophysicists, 224-234. Mavko, G. M. and Nur, A., 1979, Wave attenuation in partially saturated rocks: Geo-physics, 44, 161-178. McGinn, A. and Duijndam, B., 1999, Land seismic data quality improvements: The Leading Edge, 17, 1570, 1572, 1574, 1576-1577. References 96 Mittet, R., Sollie, R., and Hokstad, K. 1995, Prestack depth migration with compensation for absorption and dispersion, Geophysics, 60, p 1485-1494. Nemeth, T., Wu, C. and Schuster G. T., 1999, Least-squares migration of incomplete reflection data: Geophysics, 64, 208-221. Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, 1997, Numerical recipes in C: Cambrige University Press Quan, Y. and Harris, J. M., 1997, Seismic attenuation tomography using the frequency shift method: Geophysics, 62, 895-905. Ricker, N., 1977, Transient waves in visco-elastic media: Elsevier Scientific Pub. Co. . Ribodetti, A., Virieux, J. and Durand, S., 1995, Asymptotic theory for viscoacoustic seismic imaging: Annual Meeting Abstracts, Society Of Exploration Geophysicists, 631-634. Robinson, E. A. and Treitel, S., 1980, Geophysical signal analysis: Prentice-Hall, Inc. Robinson, E. A., 1985, Seismic time-invariant convolutional model: Geophysics, 50, 2742-2751. Robinson, J. C , 1979, A technique for the continuous representation of dispersion in seismic data : Geophysics, 44, 1345-1351. Sacchi, M. D., Ulrych, T. J. and Walker C. J., 1998, Interpolation and extrapolation using a high resolution discrete Fourier transform: IEEE Transactions on Signal Processing, 46 31-38. References 97 Sacchi, M. D., and Ulrych, T. J., 1996, Bayesian inversion of acoustic impedance: 1996 annual report of Consortium for the development of specialized seismic techniques (CDSST), UBC, Canada, 152-166. Sato, H. and Fehler, M. C , 1998, seismic Wave propagation and Scattering in the het-erogeneous Earth: Springer-Verlag New York, Inc. P109-145. Schoepp, A. R. and Margrave, G. F., 1998, Improving seismic resolution with non-stationary deconvolution: Annual Meeting Abstracts, Society Of Exploration Geo-physicists, 1096-1099. Schuster, G. T., 1993, Least-squares cross well migration: 63th Ann. Internat. Mtg., Soc. Expl. Geopys., expanded abstracts. P 110-113. Spencer, T. W., Sonnad, J. R. and Butler, T. M., 1982, Seismic-Q - Stratigraphy or dissipation: Geophysics, 4 7 , 16-24 Song, S., Zhang, R. and Ulrych, T. J., 2000, High resolution wavefield reconstruction in seismic migration for viscoelastic media. Journal of Seismic Exploration, 9: 1-18. Stolt, R. H., 1978, Migration by Fourier transform : Geophysics, 4 3 , 23-48. Tarantora, A., 1987, Inverse problem theory. Methods for data fitting and model para-meter estimation: Elsevier Science Publ. Co. Trad, D. S., and Ulrych, T. J., 1999, Radon Transform: 1999 annual report of Consortium for the development of specialized seismic techniques (CDSST), UBC, Canada, 105-147. Turner, G. and Siggins, A. F., 1994, Constant Q attenuation of subsurface radar pulses: Geophysics, 59, 1192-1200. References 98 Tonn, R., 1991, The determination of the seismic Quality Factor Q from VSP data: A comparison of Different Computational method: Geophysical Prospecting 39, 1-27. Ulrych, T. J., 1998, lecture notes of a Bayes tour of inversion at UBC. Ulrych, T. J., Sacchi, M. D. and Graul, J. M., 1999, Signal and noise separation: Art and science: Geophysics, 6 4 1648-1672. Varela, C. L., Rosa A., and Ulrych, T. J., 1993, Modeling of attenuation and dispersion. Geophysics, 58, 1167-1173. White, J. E., 1983, Underground Sound application of seismic waves, Elsevier Science Publishing Company Inc.. Yilmaz, 0., 1988, Seismic data processing : Society Of Exploration Geophysicists: Soci-ety of Exploration Geophysicist Publications. Appendix A Absorption Described by a Dispersion Relation The basic assumptions used to describe absorption behavior are constant quality factor Q , linearity and causality. Under this preposition, different frequency components of a wave travel at different velocities. This is described by a dispersion relation (Aki and Richards, 1980) V(UJ) 1 v0 1 + ~ 7 i l n V / ' (A.l) where Vo is the velocity at a reference frequency UJ0. Commonly u>0 is the highest fre-quency of a wave. Seismic waves traveling in absorptive media will also experience energy attenuation. The wavenumber including an attenuation term is complex 1 \ k UJ 1 + V(UJ)\ 2QJ Inserting equation (A.2) into equation (A.l), we obtain k UJ vo i - -U(^) irQ \UJ0J '2Q (A.2) (A.3) Expanding the right side of equation (A.3), and omitting terms of order equal to or greater than ^ term results in k UJ v0 UJ ( UJ \ , UJ ——In — TVQVQ \UJ0J 2QV0 (A.4) In the two dimensional situation, it is sometimes necessary to calculate k2 or Vk. These quantities are given by k2 = 2 , .1 1 —In I — ) +i~ -KQ \UJ0J Q 1 , / w \ 1 2-KQ \UJOJ AQ (A.5) (A.6) 99 Appendix A. Absorption Described by a Dispersion Relation 100 For one-dimensional wave propagation, if the wavefield U(z,u>) at a depth z is extra-polated downward to a depth z + Az, the wavefield U(z + Az, u>) obtained from equation (A.4) is U(z + Az,u>) = U(z,u>)A(Az,u))exv(i—Az). (A.7) In equation (A.7) A(z,u>) is an absorptive related term . u>Az . , UJ A(Az,w) = exp(-^^)exp 2Qv0' —i „ M — ) ivQvo U)0 = A1(AZ,UJ)A2(AZ,UJ), ( A . 8 ) where u>Az A 1 (A . , c ) = e x p ( - ^ - ) is an attenuation term and A2(Az,u>) — exp . u>Az u> ln(—) (A.9) (A.10) -KQV0 W0 ' is a dispersion term which causes a phase delay. Equation (A.7) is used in Chapter 5 to derive a method of estimating Q's. If a wavelet with an amplitude spectrum B(UJ) at time t = 0 has traveled in an absorptive medium for t seconds, equation (A.7) shows that the amplitude spectrum at time t becomes U)t S(W>0 = B ( a ; ) e x p ( - — ) or B(f,t) = B(f)exp(-^) (A.ll) 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0089920/manifest

Comment

Related Items