Index of /rozprawy2/11122

Pełen tekst

(1)University of Science and Technology AGH Faculty of Geology, Geophysics and Environmental Protection Department of Fossil Fuels. Doctoral thesis. Spectral decomposition of a seismic signal: thin bed thickness estimation and analysis of attenuating zones Anna Kwietniak ´ Supervisor: prof. dr hab. in˙z. Ryszard Slusarczyk, co-supervisor: dr in˙z. Tomasz Ma´ ckowski. Kraków 2016.

(2) Moim Nauczycielom To my Mentors 1.

(3) Contents Introduction. 4. Notation. 8. 1 Wavelet stability. 10. 1.1. Seismic wavelet extraction techniques — statistical wavelet extraction 11. 1.2. Seismic wavelet extraction techniques — deterministic wavelet extraction, synthetic seismogram for seismic-to-well tie. 1.3. . . . . . . . . . 13. Helpful tools for determining lateral wavelet stability . . . . . . . . . 14. 2 Spectral decomposition for thickness estimation. 17. 2.1. Seismic thickness estimation - state-of-the-art . . . . . . . . . . . . . 17. 2.2. Thin Bed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18. 3 Interference. 34. 3.1. Wolf Ramp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36. 3.2. Various aspects of transition zones. 4 Geological setting. . . . . . . . . . . . . . . . . . . . 42 48. 4.1. Localization of the survey area . . . . . . . . . . . . . . . . . . . . . . 48. 4.2. Geological setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49. 5 Methodology. 60. 5.1. Decomposition based on Fast Fourier transform . . . . . . . . . . . . 60. 5.2. Decomposition based on Continuous Wavelet Transform . . . . . . . . 63. 2.

(4) CONTENTS 5.3. 3. Decomposition based on Complete Ensemble Empirical Mode Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66. 6 Verification of the chosen methodology on model data. 68. 6.1. Wavelet stability verification . . . . . . . . . . . . . . . . . . . . . . . 68. 6.2. Synthetic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71. 6.3. Layer-cake model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72. 6.4. Normalization and wavelet overprint removal techniques. 6.5. Calibration of parameters used for displaying results of spectral de-. . . . . . . . 87. composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.6. Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91. 7 Results. 93. 7.1. Verification of the existence of a Wolf ramp . . . . . . . . . . . . . . . 93. 7.2. Thickness relation of the tuned Silurian and Ordovician deposits . . . 102. 8 Discussion and final conclusions. 120. 8.1. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120. 8.2. Final conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133. 9 Definitions. 136. Acknowledgements. 144. Bibliography. 145. List of Figures. 150. List of Attachments. 153.

(5) Introduction Motivation: All methods that attempt to answer a question concerning how frequency content varies with time can be understood as spectral decomposition. The concept has been studied by many scientists and become the matter of the overriding importance for the last 25 years. Modern algorithms of spectral analysis go far beyond the classical approach of Fourier transform incorporating advanced mathematical tools and various parameters. All of them, however, are based on the fact that at any given time of a signal it consists of many overlapping series whose superposition results in the final function. In this dissertation the author will comment on thin bed thickness estimation by the use of spectral decomposition algorithms and verify if the method can be applicable in the specific geological setting. Also, selective attenuating characteristics of a transition zone will be presented, and via incorporating spectral decomposition the structural interpretation of the zone will be performed. Spectral methods were have been developed and applied in various interpretational works. The subject has been also widely studied in Poland. The author would like to recommend excellent works of Polish scientists, who not only incorporated the concept in software and algorithms development but also made an attempt to translate various English terms to Polish, which is of the highest value. In these work the concept of frequency analysis is presented clearly and the author would suggest that Polish readers, who are not familiar with nomenclature, should get familiar with these papers. Frequency analysis has been the subject of many scholar discussions with the authors’s PhD supervisors and other scientists. These discussion inspired the author 4.

(6) 5 greatly and helped to better understand the frequency analysis. Dissertation layout The dissertation consists of nine chapters in which the author describes the concept of spectral decomposition application is seismic exploration. In Chapter 1 Wavelet stability , the ways wavelet extraction is presented, and its stability determination in the seismic data. Chapter 2 Spectral decomposition for thickness estimation explains the methodology of thin bed thickness estimation on the basis of spectral analysis. Chapter 3 Interference presents the concept of a transition zone (Wolf ramp) as a source of selective attenuation that manifests itself in frequency decrease. In Chapter 4 Geological setting a brief description of geology is given, including the characteristics of specific lithofacies of individual Paleozoic deposits. Chapter 5 Methodology covers the synthetic characteristics of the methodology and gives a brief description of the spectral decomposition algorithms that were implemented. It is followed by an extensive Chapter 6 Verification of the chosen methodology on model data, in which the author describes all the preparatory steps that were performed on the model data in order to determine the optimal set of parameters for frequency analysis. Chapter 7 Results, as synthetic as possible, presents the results and findings of the analysis. The dissertation is summarised in the closing Chapter 8 Discussion and final conclusions, where findings are listed. In this chapter there is a brief discussion where the author asks open questions and makes remarks about the possible sources of the interpretation of anomaly. The thesis is closed by Chapter 9 Definitions where formulas discussed in the main body are evaluated. The list of symbols used throughout the dissertation is included in Section Notation p. 8, at the beginning of the dissertation. This is done on purpose so as not to distract the reader with unnecessary formulas while reading the dissertation. Each chapter begins with short introductory notes that is followed by main body. Figures are included along with the text, bigger figures are added at the end of each chapter for clarity. List of Figures can be found at the end of the dissertation, p. 150. Four maps showing the result of analysis are attached to the dissertation. There are listed at the end of the Chapter 7 Results and at p. 153..

(7) 6. The work was part of the project Blue Gas Polish Shale Gas, a joint undertaking of National Centre for Research and Development and Industrial Development Agency (no.GASLUPSEJSM 17.17.140.86680). Author would like to thank Polish Oil and Gas Company for providing the data and would like to express much appreciation to Halliburton, Tesseral and CGG for their seismic software that was made available to the Department of Fossil Fuels and Department of Geophysics, AGH-UST for scientific and didactic purposes..

(8) 7 Main thesis to be verified: 1. Spectral decomposition is an efficient tool for indicating thin beds and estimating their thickness 2. Correlated results of spectral decomposition and seismic attributes as well as AVO analysis indicate the attenuating zones and gassaturated zones 3. Spectral decomposition enables more precise stratigraphic and structural interpretation of seismic data.

(9) Notation s(t) — seismic trace (Equation 1.1), r(t) — earth reflectivity series (Equation 1.1), w(t) — short pulse wavelet (Equation 1.1), Pf — period of spectral notches (Equation 2.1), T — temporal bed thickness (Equation 2.1), g(t, f ) — seismic response of a thin bed (Equation 2.2), r1 , r2 — reflection coefficients (Equation 2.2), f — frequency (Equation 2.2), g(f ) — complex reflectivity spectrum (Equation 2.2), re , r0 — even and odd component of the reflection coefficients (Equation 2.3, 2.4), fipeak — peak instantaneous frequency (Equation 2.6), Ah — reflection amplitude to a unit-amplitude normal incidence wave (Equation 3.1), r — acoustic impedance ratio (Equation 3.1), h — bed thickness (depth unit) (Equation 3.1), I1 , I2 — acoustic impedance (Equation 3.2), η — the ratio of a bed thickness (h) and the wavelength λ (Equation 3.3), 8.

(10) 9 Rw (f ) — frequency dependent normal-incidence reflection of a transition zone — a Wolf ramp (Equation 3.5), σ, γ — functions of frequency (Equation 3.5), k — velocity ratio (Equation 3.5), Ro — reflection coefficient (Equation 3.8), Rj − 1 — reflection coefficient for N constant velocity layer where N changes from n = 0 to N − 1 (Equation 3.11), Rj — reflection coefficient of the above (i.e. previous) transition layer (Equation 3.11), ω — angular frequency expressed in cycles per octave (Equation 3.11), zj — distance from the start of the zero coordination point, understood as the j t h layer thickness (Equation 3.11), vj − 1 — j t h velocity value (Equation 3.11) V , V0 , V1 — velocities (Equation 3.11) SD — spectral decomposition FFT — Fast Fourier Transform CWT — Continuous Wavelet Transform CEEMD — Complete Ensemble Empirical Mode Decomposition RAP — Relative Amplitude Preservation (processing and display designed to preserve relative amplitudes of seismic reflections).

(11) Chapter 1 Wavelet stability The success of applying SD for thickness estimation lays in the following statement: the difference in frequency response between a long-window and short-window amplitude spectrum is significant. The transform from a long window approximates the spectrum of the wavelet, but the transform from a short window compares a wavelet overprint and a local interference pattern representing the acoustic properties and thickness of the geologic layers spanned by the window [Partyka et al. 1999]. Regardless of the chosen method and the goal of interpretation the rate of wavelet impact must be determined for applying SD. Moreover, the seismic wavelet must remain constant throughout the whole seismic volume. Hence, the information of the wavelet itself is of the overriding importance. For minimising its impact common spectral balancing techniques are involved. These techniques rely on spark invariant stationary statistics and originally are based on equalizing each frequency slice after SD according to its average amplitude [Partyka et al. 1999]. Nowadays many spectral balancing techniques are applied with good results. The applications of different algorithms of SD (FFT, CWT, CEEMD) were performed in different programmes so a detailed analysis of wavelet behaviour was done on all the data sets. Structural interpretation requires a high resolution wavelet or more stable phase characteristic [Wardell 1983] and is necessary for the application of subsequent inversion processing for a lithologic interpretation. Before spectral analyses on a given seismic data set is performed, it is of the paramount importance to study seismic wavelet stability and consistency. Frequency 10.

(12) 11 analysis is only possible under condition that the seismic signal is consistent throughout the seismic survey. In the following Chapter a few topics on seismic wavelet are covered: 1. Seismic wavelet extraction techniques. 2. Creating synthetic seismograms for the seismic-to-well tie. 3. Helpful tools for determining lateral wavelet stability.. 1.1. Seismic wavelet extraction techniques — statistical wavelet extraction. Generally, there are two main ways of wavelet extraction. The first is based on purely statistical methods of computing wavelet based on seismic data. The second method, which is considered to be deterministic uses well logs: density and sonic logs by creating and matching synthetic seismograms to seismic data in the close vicinity of the well. Let us comment on them briefly. Statistical approach to wavelet extraction is based mainly on the autocorrelation function. Recently, however, many new algorithms that use more advanced mathematical tools have been proposed, almost all of them rely to some extend on the correlation function or its modification. Classically, the steps of statistical wavelet extraction are as follows: 1. Seismic traces are chosen from a given data set and tapered (i.e. when the z interval defined is part of z range signal outside the zone of interest is tapered, also, traces are usually chosen regularly). 2. The autocorrelation of seismic trace is computed, using the length of the desired wavelet. 3. The frequency spectrum of the autocorrelation function is computed. 4. The square root of the modulus of the frequency spectrum is taken (the zero Hertz component is muted to zero)..

(13) 12 5. The inverse FFT is computed. 6. The zero phase wavelet is the real part of the inverse FFT. In most geophysical software there is an option to extract a wavelet that is nonzero phase wavelet, but essentially zero-phase wavelets are desired: they have high resolution power (resolvability). Nowadays more advanced techniques of statistical wavelet extraction have been introduced and applied with high fidelity. Among many three statistical waveletextraction methods have been widely used [Edgar & van der Baan 2011]: 1. a kurtosis-based approach that produces a stationary, constant phase wavelet, 2. modified mutual-information-rate criterion method that generates stationary wavelet with frequency dependent phase, 3. a nonstationary extension of kurtosis-based approach that gives time-varying constant phase wavelet. Wavelets that are obtained by the use of the above algorithms show a good agreement with a deterministic wavelet derived from seismic to well tie. Regardless of the method chosen and the mathematical advancement of the algorithm applied, a statistical approach has some benefits. Firstly, statistical wavelet can be computed without any well-log data (i.e. for old wells which often have scares or unreliable logs that potentially lead to false deterministic wavelet). Different wells can and sometimes do predict different phase corrections at the nearby location [Edgar & van der Baan 2011], hence they point to suggest a non-consistent wavelet for a given dataset. In view of the above, statistical wavelet extraction enables wavelet estimation away from the well, and this statement alone is of the overriding importance..

(14) 13. 1.2. Seismic wavelet extraction techniques — deterministic wavelet extraction, synthetic seismogram for seismic-to-well tie. Well-log based synthetic seismograms are produced traditionally by time-domain convolution of a seismic wavelet with a well-log-based impedance curve. For such an algorithm two assumptions must be associated: (1) the convolution model is valid, (2) the wavelet is invariant in time. (1) the convolution model The convolution model of earth reflectivity can be described in the equation of seismic trace:. s(t) = r(t) ∗ w(t).. (1.1). where: s(t) — seismic trace r(t) — earth reflectivity function w(t) — short pulse wavelet The model says that the seismic trace is a result of convolution of the earth reflectivity function with a short pulse wavelet. Of course, the equation (1.1) is a simplified expression. In reality the convolution model is far more complicated and more factors should be included in the convolutional model equation, such as: attenuation, layering, multiples, various gradients. Nevertheless, these various effects ideally are reduced at the processing steps and we assume that convolutional model is valid for seismic-to-well tie and at this step no additional corrections are done. But it must always be remembered that such effects may and do affect the final product because it is almost impossible to get rid of them completely. The validation of convolutional model is commonly acceptable. (2) Unfortunately, the second assumption that the wavelet is invariant in time is not fulfilled in the presence of attenuation. For the modelling, usually the attenuation effects are ignored. However, it must be underlined that the process of.

(15) 14 attenuation is complex and its various effects are present in real data. Usually the synthetic seismogram is produced in two steps: first, a reflectivity function is built based on sonic log and density log from the well, and it is converted from depth domain (well) to time domain (seismic) using smoothed velocity function based on sonic log. Second, the time-domain reflectivity function is convolved with the chosen wavelet (synthetic or statistically extracted) using frequencies that are comparable with that of seismic data. In order to match the synthetic seismogram with the seismic data a cross-correlation function and cross-correlation coefficient (CCC) is computed for a given interval. Based on the cross-correlation coefficient, the above procedure is repeated with the use of altered wavelet parameters, until the satisfactory high cross correlation coefficient is obtained. To increase the cross-correlation coefficient the wavelet frequency or/and phase have to be altered in order to smooth the velocity function. It is important though not to arbitrarily stretch or squeeze the synthetic seismogram to find a very high cross correlation coefficient. Typically, in practice CCC of 0.7 is considered very good. The wavelet that is computed on the basis of deconvolution is called deterministic, because well logs determine the wavelet prosperities. Conventional deconvolution attempts to derive this simple wavelet by an analysis of impedance series in comparison with the seismic trace (or for the chosen interval with the respectable portion of trace). It should be mentioned that if the algorithm is performed on a short part of trace this procedure does not always yield a reliable wavelet. i.e. the wavelet is not stable or shows inconsistency. As a result, the reflection character of seismic data cannot be assigned to lithological changes with high fidelity. For wavelet extraction it is advisable to use windows not shorter than 1 second (for a typical seismic survey of the time range about 2 − 4 or 5 seconds).. 1.3. Helpful tools for determining lateral wavelet stability. Wavelet stability has to be verified before any spectral analysis is applied. Various procedures are used to determine how wavelet changes throughout seismic data.

(16) 15 volume. The general assumption is to study wavelet prosperities: phase spectrum, amplitude spectrum with respect to time. Wavelet extraction is a repeatable process and requires time and effort to analyse the output critically. In this step of interpretation seismic attributes can be applied, especially instantaneous attributes. Definitions of instantaneous attributes are introduced thoroughly in Definitions. The main idea behind applying instantaneous attributes is that a seismic signal is considered to be non-stationary, i.e. its characteristics vary with time. Several comments are to be mentioned for a non-stationary process with respect to frequency. Based on the classical definition of frequency (frequency is inverse to the period) two major obstacles arise. (1) The definition is oblivious to the fundamental wave conservation law, which requires the wave-number and frequency to be differentiable. In this case it cannot be, because it is constant over the whole wavelength range. (2) The definition fails to work in non-stationary especially in nonlinear processes with intrawave frequency modulations. In a complicated vibration, there may exists multi-extrema in the registered sequence. This observations enabled the concept of instantaneous frequency: the need for frequency as a function of time and the fact that frequency should be the function of time and having an instantaneous value can be justified from both mathematical and physical grounds [Huang et al. 2009]. The most useful application of instantaneous amplitude and phase is based on following property: Fourier phase equals the instantaneous phase at the envelope peak for an isolated band-limited wavelet. This characteristic works well as a tool for detecting lateral wavelet stability [Perz et al. 2004]. It also was stated [Barnes 1991] that for any constant-phase wavelet the instantaneous frequency evaluated at the envelope peak is equal to the average Fourier spectral frequency weighted by its amplitude spectrum. To apply these methodology in practice, first the stable and well-defined seismic reflector has to be found and picked. Then the instantaneous amplitude and the phase are computed and the time of the maximal instantaneous amplitude has to be defined. The next step is to evaluate the instantaneous phase for the maximal instantaneous amplitude, which is good approximation of the Fourier phase (phase spectra). Although the method is straightforward it still has some.

(17) 16 disadvantages: (1) The need for a stable, regional marker, which does not always occur, (2) it assumes that the wavelet is a constant-phase, which is not necessarily true, i.e. in reality the wavelet suffers from frequency-dependent phase and amplitude variations caused by physical processes [Perz et al. 2004]. The method can reveal valuable information about the data quality..

(18) Chapter 2 Spectral decomposition for thickness estimation Spectral decomposition provides means of utilizing seismic data and the discrete Fourier transform (DFT) for imaging and mapping temporal bed thickness and geologic discontinuity over large 3D seismic surveys [Partyka et al. 1999]. The approach suggested in Partyka’s classical paper [Partyka et al. 1999] on the use of SD for bed thickness estimation is described in the following chapters. The methodology herein described was applied to the synthetic 2D model in order to test whether it is applicable in a simple geological setting. These results are discussed in this Chapter. A few more methods of incorporating frequency changes into thin bed thickness estimation are also introduced, as well as different models of thin beds itself.. 2.1. Seismic thickness estimation - state-of-the-art. Seismic layers that are consider to be seismically thin would probably be of a considerable size when examined on the Earth surface. Thicknesses that we deal with during thin bed thickness analysis are in order of tens of meters. The resolvability of the seismic wavefiled depends on acquisition geometry and its parameters as well as on processing procedures, their characteristics and sequence [J¸edrzejowska — Tyczkowska 2012]. The smallest thickness that is observed on any seismic section 17.

(19) 18 is wholly dependant on such seismic wavelet proprieties as wavelength, frequency and signal energy, and because of that the thickness of seismic section cannot be considered without careful examination of wavelet characteristics.. 2.2. Thin Bed. M.B.Widess is considered a founder of thin bed analysis. He himself, as well as many more scientists concerned with answering a question: How thin is a thin bed?. This topic is widely discussed in numerous papers and the scope of this chapter is to introduce the main definitions and discuss them briefly. Generally thin bed models are built from three layers, with different density and velocity values in each of the layer. Due to this reflection coefficients that exist at the top and bottom of the middle layer can vary significantly. Such a model is known as single layer model. Nevertheless, any reflection coefficient pair can be decomposed into even and odd components, with the even components having equal magnitude and sign, and the odd components having equal magnitude and opposite sign [Puryear & Castagna 2008]. The Widess model assumes that reflection coefficient pairs are perfectly odd (see Figure 2.1, A). The second classic model occurs when reflection coefficients are even. A. Example lithology. Acoustic impedance. RC series. Shale Sandstone Shale Increasing. B. Example lithology. Acoustic impedance. Negative. Positive. RC series. Shale Sandstone Limestone Increasing. Negative. Positive. Figure 2.1: Thin bed cases. (A) Odd reflection coefficients, (B) Even reflections coefficients, [Liner 2012].. There are many reasons why thin bed thickness is such an important subject. The motivation can be purely scientific but also it may be of utilitarian value..

(20) 19 Many hydrocarbon reservoirs across the globe are embaded in thinly layered, often intermittent horizons packets of permeable and impermeable deposits. This problem is for instance commonly encountered in the Miocene strata sequence within ˙ lawi´ the Carpathian Fordeep [Zu nski 2010] where hydrocarbon saturations occur in thinely-layered shale-sandstone packets. Hence bed thickness estimation, as well as its detactability has become a subject of the paramount importance. There exist many approaches for quantitative an qualitative predictions of seismically thin beds. The definition of a thin bed has been evolving for many years with advances of seismic acquisition and processing. Form the point of view of information some of the methods will be presented here briefly as well as the conceptual model, whereas the subject will be continued in the next Chapter Interference. Compared with seismically thick beds, thin beds are substantially more difficult to interpret [Zeng & Marfurt 2015]. Thin bed may be understood as a layer of different properties than its upper and lower boundaries, and if thicker, it would produce a visible seismic response at these boundaries but in a situation when it is thin such a pattern is not visible in a seismic response. Thin beds produce reflection at its top and bottom but these two reflections interfere with each other. Being dependent on the exact differences between the pairs of volumes: upper — thin bed and thin bed — lower the interference may be constructive or destructive for the resulting reflection, but only the overlapping result of these two is visible. To address the problem the limit of seismic resolution was studied by many researches. In the case of 2D seismic data two approaches, both referring to the wavelength (λ), defining the limit of seismic resolution are widely known. The classical Widess criterion is that the limit is defined by the minimum distance at which a composite waveform stabilizes as the derivative of the waveform from an individual reflection [Widess 1973]. The Widess proposed λ/8 as the resolution limit. This definition was stated with the use of symmetrical wavelet, hence it is mainly of theoretical value due to the difficulties in determining wavelet stability (symmetry), i.e. the exact shape at the boundaries. A more applicable resolution limit is known as the Rayleigh criterion of peak-to-through separation of λ/4 for a thin bed with the same polarity at the top and bottom of a seismic layer in.

(21) 20 question [Kallweit & Wood 1982]. It also corresponds to the tunning point for which the amplitude reaches its maximum for an opposite-polarity model. Interpretation of more than one thin bed below this limit is barely possible. For 3D data the resolution limit can be significantly increased — 3D data has better seismic resolution (both time resolution and horizontal resolution, and the resulting spatial resolution — the ratio of the previous two). For the time being, the Sheriff limit, which is commonly accepted among interpreters equal to λ/25. This number most probably is going to be smaller with the advances in data acquisition and processing [Zeng & Marfurt 2015] that aim in the widening of the spectral characteristic of a seismic signal [J¸edrzejowska — Tyczkowska 2012].. Spectral decomposition for thin bed thickness estimation Seismic decomposition as a tool for thickness estimation goes beyond the classical tunning thickness estimation proposed by Widess [Widess 1973]. These traditional techniques rely on careful processing and are very sensitive to phase changes, hence require zero-phase signal. SD has two main advantages over this method: It provides a phase-independent result and requires only one guide pick (horizon, event) within the seismic zone of interest [Partyka 2001]. SD used for thickness estimation is based on short window discrete Fourier Transform — SWDFT. The most important parameter for STDFT computation is the length of the window in which the transform is computed. The difference in frequency response between long and short window is significant: transform performed in a long window approximates the spectrum of the wavelet, whereas the transform from a short window comprises a wavelet overprint and a local interference pattern representing the acoustic properties and thickness of the geologic layers [Partyka et al. 1999]. These properties are shown in Figure 5.1. Long window include many geologic features that are statistically random. Hence the interference patterns of individual thin beds give a spectrum which is white [Partyka et al. 1999]. The convolution of a source wavelet with a random geological sequence creates an amplitude spectrum that resembles the wavelet. On the contrary, the spectra obtained from short-windows are sensitive to the local changes of geology and preferably reveal.

(22) 21 the acoustic properties and thickness of the geological body. This phenomenon gave reasons for applying a short time window for computing DFT. The window used in the process of decomposition has a fixed length; it starts from the first sample of a seismic trace; DFT is computed in the window and then, with the defined step the window moves down (forward in time) and another computation is performed. The DFT is computed until the window reaches the last segment of the trace, which is equal to its length. The whole process results in the volume, that shows local frequency spectra, i.e. spectra for a given time. Apart from its interpretative value this, algorithm enables, at least theoretically, to estimate a bed thickness by inverting the frequency spectra using the high frequency components. Qualitative analysis Thin bed thickness estimation is addressed by two approaches that are based on SWDFT [Partyka et al. 1999]. The first uses spectral-decomposition-derived dominant frequency and dominant amplitude mapping, while the second incorporates spectral-decomposition-derived thickness via discrete Fourier components. Dominant frequency and dominant amplitude The method requires two attributes to be isolated after the decomposition process: maximum spectral amplitude (1) and frequency at which the first frequency peak occurs (1st frequency mode (2)). The maximum spectral amplitude has to be understood as the first amplitude peak that occurs in amplitude frequency spectra. This analysis corresponds to the classic peak-to-through separation but requires only one seismic event to be mapped, not two horizons as in the classical analysis. By looking closely to the dominant amplitude and 1st frequency mode changes, thickness estimation can be performed. However, 1st frequency mode can be used to quantify thickness that is greater than tunning thickness, and, dominant amplitude can be used to quantify thickness that is less than the tuning thickness. Discrete Fourier components Analysis of discrete Fourier components enables robust estimation of the bed thickness. The factor that should be closely mapped by this approach is the period of spectral notches that occurs for different thicknesses. The choice of the use transform can not be accidental — the characteristic notches in the frequency spectra occur.

(23) 22 ˙ lawi´ only for Fourier transform of short window (SWDFT) [Zu nski 2012a]. It should be mentioned here, in order not to confuse terms period and frequency that period of spectral notches is expressed in frequency units (Hz). Roughly speaking, the longer the period the thinner the seismic bed — spacing between the peaks and notches is precisely the inverse of the layer thickness. The above-mentioned attributes have to be carefully separated after the decomposition process. To verify this methodology the model of a thin bed (single layer model) was created and by the concept of convolution its seismic response to the Ricker wavelet was obtained. It should be mentioned that the original algorithm is incorporated in SpectralDecon plugin that is available in OpenWorks platform (Halliburton). Only recent spectral decomposition has become more popular and can be found in numerous commercial software. The analysis presented herein was mainly done in OpendTect. It is important to identify correctly dominant amplitude and 1st frequency mode. Once they are properly defined and their values are known, the thin bed estimation can be continued. For a real data, a step of the overriding importance is connected with wavelet overprint removing. As was described in the previous Chapter Wavelet Stability, the type and characteristics of the seismic wavelet can greatly influence the seismic data. We have to be certain that the parameters and their changes are related directly to geology. Approaches of removing the wavelet overprint that are widely used are called spectral balancing techniques. For the use of the a/m algorithms the seismic wavelet has to be assumed constant at all locations. In order to verify the methodology described above the synthetic model of a thin wedge was created. Thickness of the wedge increased linearly from zero to fifty milliseconds and model was of Type I reflectivity. In Figure 2.2 the concepts of thin bed thickness estimation are compared and summarised. Figure in A and B come from the orginal paper of Partyka [Partyka et al. 1999], whereas Figures C and D show the result of the experiment. Figure 2.2 presents the dominant amplitude and 1st frequency mode after SWDFT for a thin bed model as a function of temporal bed thickness. With the increasing temporal thickness dominant amplitude increases up to the tunning thickness, then stabilises and afterwards increases slowly. The 1st.

(24) 23 frequency mode decreases slowly from rather high frequencies (in order of @ 35 Hz) and when it reaches tunning thickness it decreases rapidly to the frequencies in order of about 10 Hz. Izo-frequency panels available after spectral decomposition process were carefully studied to map the dominat amplitude and 1st frequency mode changes as a function of the thin bed. Figure C shows 1st frequency mode as a function of a temporal bed thickness. It can clearly be shown that 1st frequency mode decreases with the increasing temporal thickness. Result presented in Figures C is consistent which proves the method to be valid and applicable for a simple 2D model. With the use of discrete Fourier Components technique thickness estimation was possible. As is visible in Figure D the method correctly estimated the temporal. 1st frequency mode [Hz]. Dominant amplitude. thickness of the model.. (B) Spectral-decomposition-based approach uses amplitude variability at discrete frequencies to determine thickness [Partyka 2001].. (C) 1st frequency mode as a function of bed thickness.. (D) Thickness estimation basing on discrete Fourier components technique.. 1st frequency mode [Hz]. (A) Dominant amplitude and 1st frequency mode as a function of bed thickness [Partyka 2001].. Figure 2.2: Concept on thin bed thickness estimation via the use of algorithm presented by Partyka [Partyka et al. 1999], (A) Frequency and amplitude mapping, (B) Discrete Fourier Components, (C) 1st frequency mode as a function of thin bed thickness, (D) Thickness estimation basing on discrete Fourier components approach..

(25) 24 Quantitative analysis Quantitative analysis, i.e. the actual value of the temporal bed thickness is also possible by means of decomposition algorithm. Information that is available after spectral decomposition can be directly used for thin bed thickness quantitative estimations. Two equations, both introduced by Partyka in his first work on the application of spectral decomposition [Partyka et al. 1999], state that the following relations are valid:. Pf =. 1 , T. (2.1). where: Pf — the period of the notches in the amplitude spectrum with respect to frequency [Hz], T — thin bed thickness [s] (one-way time) [Partyka et al. 2000]. The equation (2.1) can be used where the period of the notches with respect to frequency can be defined. Izo-frequency panels that are available in OpendTect may provide such information with respect to frequency. In some cases, however, the exact period is difficult to detect, and in such cases we can estimate the period by the occurrence of the 1st frequency mode. In most cases the period of the notches is simply twice the 1st frequency mode. This kind of approach works well for the model data, but not necessarily for the real data. Usually, the period of the notches is visible when the amplitude contrast of a thin bed is significant, i.e. the differences in acoustic impedance is distinguishable. When a weak seismic response is produced by the thin bed such notches may not be present. The applicability of the algorithm was tasted on the same wedge model presented earlier. The results agreed well with the real bed thickness what is presented in Figure 2.2, D. In Figure 2.3 the influence of the window length is visible — too short windows fail to estimate the temporal thickness properly. For the case window length of 60 ms gave the best-fitting results. The 60 ms window encloses the whole wedge model..

(26) 25. Figure 2.3: Thicknesses of the wedge estimated by spectral decomposition with the use of FFT method. Window length influences the results.. Different thin bed models The algorithm described above assumes that the reflection in the top and bed of the thin layer has the same strength and opposite polarity, with the positive reflection in its top and negative in its bed. It is just one possibility in a continuum of the possible reflection configuration pairs [Puryear & Castagna 2008]. Nevertheless, such a concept is widely incorporated for the model studies on thin bed thickness estimation. The very useful characteristics of any given reflection pair is that regardless of the value and sign of the reflection coefficient at the top and the bottom they can always be understood as a superposition of even and odd components (see Figure 2.4) [Puryear & Castagna 2008]. In this case even components have equal magnitudes and signs and odd components have equal magnitudes and opposite signs. The method of splitting the original reflection pair into even and odd components is not only helpful but also elegant from the mathematical point of view. If we consider the seismic response from a thin bed and start the analysing it at the centre, we will result in the equation [Puryear & Castagna 2008] expressed as:. g(t, f ) = r1 exp[−i2π(t − where: r1 — top reflection coefficient, r2 — base reflection coefficient,. T T )] + r2 exp[−i2π(t + )]; 2 2. (2.2).

(27) 26 t — time sample, T — temporal layer thickness, f — frequency. By using trigonometric identities the real part and imaginary part of (2.2) for a complex spectrum g(f )can be defined respectively:. Re[g(f )] = (2re ) cos(πf T ),. (2.3). where: re — even part of the reflection coefficient, and. Im[g(f )] = (2ro ) sin(πf T ),. (2.4). where: ro — odd part of the reflection coefficient. Original reflection coefficient pair. Even part. time domain. r1 T. A. Odd part. r2 t. Fourier transform 0.35. 0.30. 0.30. 0.25. 0.25. 0.20. 0.20. 0.15. 0.15. 0.10. 0.10 0.5. 0.5 0.00. 0. 100. 200. 300. Frequency [Hz]. 400. 0.00. 0. 100. 200. 300. Frequency [Hz]. 400. frequency domain. Amplitude. Odd part. Even part. 0.35. B. Figure 2.4: Reflection coefficient pair (r1 is twice as big as r2). Any given pair can be expressed as even and odd part. (A) Time domain — original reflection coefficient, even and odd part of the reflection coefficient, (B) Frequency domain — amplitude versus frequency plots for even and odd components of the reflection coefficient pair (r1 = 0.2 and r2 = 0.1).. For thin bed thickness estimation the odd components are the worst case, because when thickness approaches zero, reflections from top and bottom interfere destruc-.

(28) 27 tively. Seismic resolvability improves a lot when reflections are even, because then they interfere constructively. There are four widely accepted reflectivity types that are enumerated with roman numerals: I, II, III and IV . They were widely studied by Chung [Chung 1994], and the results of these studies can be found in Chapter 2 of his doctoral dissertation. Type I reflectivity has negative reflection coefficient in top, positive in bottom, with equal magnitudes. Type II has equal and same — sign reflection coefficients. Types III and IV are neither odd nor equal: they have different reflection coefficients in top and bottom, i.e. Type III has negative reflection coefficient in the bed top, positive in bottom and the top has smaller magnitude than the one from its bottom, Type IV has both positive reflection coefficients, but the one from top is twice as week as the one from bottom. For modelling purposes Type III and IV reflectivity series were decomposed into the sum of Types I and II and superposition was applied to yield their seismic responses. For these models the amplitude reflection (seismic response) for thicknesses less than 1/4 λ differs from the Widdes model. The analytical solution was derived for sinusoidal and Ricker wavelets. The results of the analytical and numerical solutions are included in two papers, simply repeating the solutions described by Chung in his doctoral dissertation. They can easily be found: describing the amplitude response of thin beds by the sinusoidal approximation and Ricker approximation can be found in Geophysics [Chung & Lawton 1995a], and frequency characteristics of seismic reflections from thin beds [Chung & Lawton 1995b]. The equations derived in these papers are widely accepted and quite accurately agree with the numerical solutions. They were incorporated and used in many other papers, e.g., layer-thickness estimation with the use of spectral inversion [Puryear & Castagna 2008]. This conception leads to two main disadvantages of the model presented by Partyka. One was previously mentioned: the algorithm is dedicated and was tested only for Type I reflectivity model. The other disadvantage is that the notches in the amplitude frequency spectra and the period of their occurrence is very hard to detect in the real data examples. Due to these problems, different methods of detecting.

(29) 28 and estimating thin bed thickness were introduced. These methods respect Chung and Lawton equations [Chung & Lawton 1995b],[Chung & Lawton 1995a], as well as Partyka’s observations on the behaviour of frequency. It was stated afterwards that mainly the frequency response of a thin bed should be used predominantly for thin bed thickness estimation, i.e. the frequency analysis should be in general less sensitive to changes in the reflection coefficient as compared to the amplitude (Widdes) technique [Liu & Marfurt 2006]. It should be mentioned, however, that for modeling purposes it is common to ignore transmission losses and internal multiples. Such effects are negligible provided the acoustic impedance ratio between the thin layer and the surrounding rock lies between the bounds of 0.5 and 2 [Koefoed & de Voogd 1980], which is the range usually encountered in clastic sequences. Thin bed detection — other algorithms Frequency changes have been widely studied for thin bed thickness detection. Different algorithms used for this purpose were listed and briefly commented in the August 2015 issue of the Society of Exploration Geophysicists journal Interpretation. An introduction to the volume Special Section enumerated seven algorithms: instantaneous attributes, spectral decomposition and continuous wavelet transform, frequency spikes and phase residues, singularity/spice, spectral inversion, spectral balancing; the majority of them are associated with frequency changes. Keeping in mind the range of possibilities of thin bed models, as well as frequency characteristics of their seismic response instantaneous attributes, especially instantaneous frequency should be mentioned here. A detailed derivation of its formula as well as the exact definitions are included in chapter Definitions,( no. 4). Instantaneous attributes are classical and widely used tools for a seismic interpretation and have become a part of a standard seismic analysis. From the above listed methods of thin bed identification, the one which is associated with instantaneous attribute analysis will be described here in short. Afterwards a simple method of defining thin bed thickness that incorporates peak instantaneous frequency will be discussed. Thin bed detection and correlation with instantaneous frequency Spikes in instantaneous frequency (sudden increases and decreases in the at-.

(30) 29 tribute)) are huge transient values, and though they may appear unrealistic (they exceed the spectral bandwidth) thus are a direct consequence of the definition of instantaneous frequency. Their magnitude can exceed the Nyquist frequency [Barnes 1991], and also cross the low-frequency limit, even being negative [Zeng et al. 2010]. These abnormal frequency behaviour is linked to the seismic signal interference phenomena, hence these changes in the frequency values can be used as a thin bed detection technique. The existence of spikes was doubted previously, because it was claimed that abnormal frequency increases occur where the signal is poorly defined. Nevertheless the recorded seismic trace is inherently deterministic and is no less defined at envelope minima than at envelope maxima [Barnes 1991]. These spikes may be closely linked to thin beds and have a certain geologicinterpretation meaning. There are some works saying that abnormal instantaneous frequency peaks may be linked to half-wavelength thickness in a wedge model, discontinuities in a gas reservoir succession, lower coherency indicators, flow-unit boundary and stratigraphic terminations in evaluating reservoir compartments [Zeng 2010]. With the use of different thin bed models (Types I — IV ), it can be found that frequency spikes are closely associated with changes in bed thickness. For showing these changes a simple thin bed model will be introduced. In Figures 2.5 and 2.6 the application of instantaneous amplitude and instantaneous frequency is shown for a thin bed model. Thin bed has thinned from right to left from the temporal thickness of 50 ms to 0. The models were created by the convolution of 30Hz-Ricker wavelet with the reflection coefficients of the same magnitudes and odd polarities (Type I) 2.5 and even polarities (Type II) 2.6. For these models instantaneous frequency shows characteristic spikes in its value that are concentrated approximately in the middle of the wedge for Type I reflectivity and in 1/3 for reflectivity Type II. The instantaneous frequency is apparently related to the bed thickness. The spikes that occur after applying instantaneous frequency are classified by Types I and II [Zeng 2010], but in order not to confuse them with the reflectivity types introduced earlier, they will be refer as type A and B respectively (type I denoted by [Zeng 2010] = A, type II denoted by [Zeng 2010] = B )..

(31) 30 Type A reflectivity spikes (Figure 2.5): If we concentrate on the central part of the wedge and investigate the instantaneous frequency changes it can be noted that the instantaneous frequency decreases with increasing thickness and then reverts to the negative values as the thickness approaches the merging point. Then, abruptly its value changes to bigger, normalpositive and afterwards decreases again. However, if the average instantaneous frequency had been taken the negative value would have disappeared, so the spike, treated here as an anomaly would not have occurred. B. 50 ms. A. Amplitude Max. Min. C. D. Frequency Max. Min 0. Figure 2.5: Thin bed model (reflectivity type I) and its instantaneous attributes. (A) Thin bed model convolved with 30 − Hz Ricker wavelet, (B) instantaneous amplitude of the model from (A), (C) instantaneous frequency in wiggle display, (D) instantaneous frequency in colour display. Black arrow indicates instantaneous frequency spikes.. Type B reflectivity spikes (Figure 2.6): Similar behaviour of the instantaneous frequency occur when the central part of the wedge of Type II reflectivity is considered. The difference here is that the spike is visible for thinner part, which is approximately λ/4 — the tunning thickness. Type B spikes correspond to the peak-to-trough separation of a seismically thin bed closer to the tunning thickness..

(32) 31 B. 50 ms. A. Amplitude Min. Max. C. Frequency Max. Min 0. Figure 2.6: Thin bed model (reflectivity type II) and its instantaneous attributes. (A) Thin bed model convolved with 30 − Hz Ricker wavelet, (B) instantaneous amplitude of the model from (A), (C) instantaneous frequency in colour display. Black arrow indicates instantaneous frequency spikes.. The method described above is straightforward and fast. It can be performed at every interpretational platform - instantaneous frequency is commonly available in basic interpretational packages. The exact type and nature of the frequency spike depend on acoustic impedance profile, bed thickness and bed spacing (in the case of a multibeded model). It should be noted here, however, that the abnormal frequency spikes may be related to different features that are observed in seismic data, i.e. discontinuities (such as facies boundaries, flow barriers, etc.), the presence of noise in a seismic section [Zeng 2010]. Due to these difficulties, as much geologic information as possible has to be integrated to define possible sources of abnormal instantaneous frequency spikes. Thin bed detection and correlation with peak instantaneous frequency When we deal with low signal to noise ratio instantaneous attribute may be very difficult to interpret. Therefore, a method of incorporating peak instantaneous frequency was introduced [Liu & Marfurt 2006]. The method was tested on the model data. A trend in changes in peak instantaneous frequency is consistent with the relationships presented earlier, i.e. a thicker thin bed has lower leap instantaneous frequency, while a thinner thin bed has higher peak instantaneous frequency. The relation between peak instantaneous frequency and instantaneous frequency for a model using thin bed (reflectivity Type I) and a Ricker wavelet is defined.

(33) 32 similarly by different authors. Nevertheless, its values are somehow close to one another: The relation between peak frequency and peak instantaneous frequency [Robertson & Nogami 1984]: fipeak ≈ 1.1283fp ;. (2.5). The relation between peak frequency and peak instantaneous frequency [Liu & Marfurt 2006]: fipeak ≈ 1.3293fp ;. (2.6). The difference in the value of peak instantaneous frequency is a result of different assumptions used by the authors. The value in Equation (2.5) was evaluated for single zero-phase Ricker wavelet for a thick bed seismic response. In the contrary, value in Equation (2.6) was derived for a case when time separation between top and bottom reflection coefficients is close to zero — thin bed situation. This behaviour of peak instantaneous frequency gives reason to assume that there exist a relation between thickness and peak instantaneous frequency. In the work of Liu [Liu & Marfurt 2006] the behavior between peak instantaneous frequency and thin bed thickness for a model example was analysed and the authors found that it has an inverse relationship with the wedge thickness. Peak Instantaneous Frequency [Hz] 48 46 44 42. fi peak 1.33fp. Frequency difference. 40 38 36 34. 1/4l. Thickness [m]. Figure 2.7: Peak instantaneous frequency versus wedge thickness for two wedge models (First wedge model with −0.01 and 0.01 RC — solid black line, second wedge model with −0.01 and 0.009 — dashed red line. Wedge thickness 25 meters [Liu & Marfurt 2006])..

(34) 33 For the two models the difference is very small (deviation between the solid and dashed line in Figure 2.7. The phenomena may be used for quantitative thin bed thickness estimation. In the example [Liu & Marfurt 2006] the following algorithm was applied: 1. Instantaneous frequency was computed. 2. In a 9 ms window around the zero-crossing of a selected seismic horizon (Caddo limestone, Fort Worth Basin, USA) the peak instantaneous frequency was computed and displayed as a map view. 3. The thickness of the horizon from 16 well logs within the 3D seismic section was taken. 4. Information from 16 wells was cross-plot with the value of peak instantaneous frequency taken from the exact point, where the wells were situated. 5. The thickness of the bed was computed on the basis of the inverse linear trend. The thicknesses from the wells agree nicely with the thickness taken from the attribute. The method, which is simple and only slightly interferes in the seismic data is considered to be straightforward and easily-applicable. Peak instantaneous frequency may be verification method..

(35) Chapter 3 Interference In the following section the frequency behaviour of a special case of a seismic boundary response is considered — a response from a transition zone. Two models will be described: a model with two neighbouring reflections and a special case of reflection from a gradational transition zone. The ideas presented in this Chapter are taken from work of Liner [Liner 2012] where the author presents the idea of a Wolf ramp — transition zone.. Seismic response from a thin bed The reflection amplitude Ah to a unit-amplitude normal incidence wave on a bed thickness h is: Ah = (1 + r)2 [(2rcot(2πη))2 + (1 + r2 )2 ]−1/2 ,. (3.1). where r — acoustic impedance ratio at the top of the middle layer, η — the ratio of bed thickness to wavelength in the thin bed. These factors are defined respectively: r=. I1 ; I2. 34. (3.2).

(36) 35. η=. h ; λ. (3.3). In equation (3.2) symbols I1 and I2 are acoustic impedance values of the upper and the lower layer, h in (3.3) is thickness of the layer. In the case of long-wavelength, or a very thin bed (when η → 0) the approximate form of the equation (3.1) can be written: π(1 − r)2 η Ah ∼ . = r. (3.4). Interference beyond the thin bed The model of a thin bed presents only one possibility of wave interference that exists in broadband seismic data. One of the most powerful tools to unravel interference patterns is spectral decomposition. Interference phenomenon, as frequencydependent, can be highly sensitive to frequency change. To depict this relationship [Liner 2012] suggested analysing a reflection-coefficient series that shows sonic velocity plotted against reflection time. Generally, such a complex pattern can be created via spectral decomposition of a random seismic trace from a given seismic survey (Figure 3.1). Figure 3.1, as a time-frequency dependent plot shows how the seismic trace appears as a function of bandwidth. This kind of plot might yield information that is not visible on the broadband seismic data. Some features are present across the whole frequency range at almost the same time i.e. they are consistent throughout the whole bandwidth, whereas some peaks may be phase-shifted regarding the specific frequency (event peak time varies with frequency). The main conclusion that can be drawn after analysing of Figure 3.1 is that a geological horizon can be represented by a peak, trough, zero crossing or none of these on various bandwidth versions of data. Secondly, spectral decomposition in general is a tool for unraveling information that is hidden by broadband interference..

(37) 36 A. B 0. 30. C. [Hz] 40. 50. 60. 70. 80. 90 100 110. m/s 0. 5000. Time [s]. 2.95. 3.00. 3.05. 3.10. RC. Bandpass 5 to 110 Hz. Vp. Figure 3.1: (A) Reflection coefficient series, (B) Seismic RC series shown through the application of various bandpass filters, (C) Sonic log velocity (VP ). [Liner 2012].. 3.1. Wolf Ramp. Let us briefly discuss a reflectivity dispersion which refers to frequency dependance of a normal incidence reflection effect. Intrinsic dispersion is present for a single reflection event whereas apparent dispersion is spectral modification associated with interference among several events. A combination of these two types is called mixed dispersion. There are a few causes for reflectivity dispersion, e.g., rough-surface scattering, reflection from an interface porous media, vertical transition zones. For a complete understanding of this phenomenon and development of quantitative seismic interpretation methods, it is necessary to model and identify characteristics that can be used as tools for distinguishing between them. The conclusion suggested by the above description, as well as many interpreting tasks have shown that reflection events that are observed on broadband seismic data are always frequency dependent. In this subsection, a linear velocity transition zone, known as a Wolf ramp, will be introduced. The theoretical discussion results in a ramp-detection method can be applied to migrated seismic data. An exact theoretical review of a Wolf ramp problem and its solution is given in Chapter Definitions. The model of a Wolf ramp is shown in Figure 3.2. Regardless of the values of specific parameters a Wolf ramp is given the following description: it is built of.

(38) 37. Velocity [m/s] V1. h[m] V2. Figure 3.2: Velocity model of a Wolf ramp.. two zones of different apparent velocities separated by a transition zone of a given thickness h, which exhibits linear change of velocity from the upper velocity zone to the bottom zone. The frequency-dependent normal-incidence reflection coefficient for a Wolf ramp is given by equation (3.5):. Rw (f ) =. 1 . 2σ + 2γ coth[γ log(k)]. (3.5). In equation (3.5) k is velocity ratio of the upper and lower part of the transition zone, frequency f is in Herz, transition zone thickness is h and σ and γ depend on frequency by: i2πf h ; (k − 1)V r 1 γ(f ) = + σ2. 4. σ(f ) =. (3.6). (3.7). In equation (3.6) symbol V is a velocity of an upper layer. A classic reflection coefficient can be thought of as a Wolf ramp with zero transition zone thickness so that the upper and lower half-spaces are in direct contact. For this case, a Wolf ramp reflection coefficient is defined: lim Rw (f ) =. (f,h)→0. 1 k−1 = = R0 . coth[1/2 log(k)] k+1. (3.8). In the limit of zero frequency, the Wolf and classic reflection coefficients are again.

(39) 38 equal because the transition zone thickness is vanishingly small compared with the wavelength. The Wolf ramp reflection coefficient is a complex function of frequency and decreases in magnitude as frequency increases (a behaviour that can be detected by spectral decomposition). The expression in equation (3.5) is zero if either σ or coth(γ(logk)) is infinitive, the latter occurring when γ(logk) is a non zero, purely imaginary integer multiple of π. The nth zero-crossing frequency is given by the equation: V (k − 1) fn = 2h. r (. n 2 1 ) . logk 4π 2. (3.9). Let us consider two synthetic examples of a Wolf ramp. The parameters for both models are given in Table 3.1. Table 3.1: Two models of Wolf ramp used for analysis. Parameter Model 1 Model 2 h [m] 10 50 V [m/s] 3500 3500 k [-] 0.8 0.8. A. B. C. Figure 3.3: Frequency behaviour of the exact normal-incidence reflection coefficient Rw (f ) for a Wolf ramp, model 1. (A) Real part of R(w), (B) Imaginary part of R(w), (C) Modulus of R(w).. Figure 3.3 shows the Wolf reflection coefficient by plotting Rw (f ) from 0 to 100 Hz for Model 1. In this case, we can deal with high-velocity cap rock overlying.

(40) 39. A. B. C. Figure 3.4: Frequency behaviour of the exact normal-incidence reflection coefficient Rw (f ), model 2. (A) Real part of R(w), (B) Imaginary part of R(w), (C) Modulus of R(w).. and grading into low-velocity reservoir rock. In Figure 3.3a) a real part of Rw (f ) is presented. It decreases slowly from the value −0.1 and crosses 0 at 80Hz. The imaginary part is zero at zero frequency and slowly grows to the maximum at 55Hz. Superposition of the imaginary part and real part give the total value of Rw (f ) for a given frequency. The amplitude does not exhibit zero. Equation (3.9) indicates a first zero-crossing for this model at 157Hz. The phase is linear with frequency: 180◦ at 0Hz and it decreases moderately to the value of 60◦ . The frequency behaviour for Model 2 is shown in Figure 3.4. For this model all significant change occurs below 20 Hz, above which Rw oscillates around zero (the waves do not see a reflecting interference but smoother, linear medium). Waveform simulation and detection [Liner 2012] proved that the theoretical approach presented above is valid in seismic interpretation and that Wolf ramp can be a possible solution to a specific seismic pattern. The synthetic example of a Wolf ramp is shown in Figure 3.5. The reflection waveform in Figure 3.5 (B) is created by summing the real part of the complex reflection coefficient across frequency. In conclusion, spectral decomposition of the time-domain waveform can image time-frequency behaviour of a Wolf ramp and enables a kind of inversion in which spectral decomposition of seismic trace could be inspected for characteristics of ramp reflectivity. In that way it can indicate some model parameters like ramp thickness and velocity ratio..

(41) 40. Time [s]. A. B Amplitude -0.02 0 0.02 -0.10 -0.10 -0.05. -0.05. 0. 0. 0.05. 0.05. Frequency [Hz] 10 20 30 40 50 60 70 80 90. C Velocity [m/s] 3500. 50m. 0.10 Summed. 2800. 0.10 Time-Frequency (real). Model. Figure 3.5: Convolutional 1-t 100 Hz waveform modeling of a Wolf ramp. (A) Timeddomain response by summation over frequency, (B) real part of time-frequency Wolf reflection coefficient formed by convolving the complex Wolf reflection coefficient with a unit spike. (C) Velocity model [Liner 2012].. The field data that test a given methodology for the Wolf ramp comes from the Andaman Basin, offshore Thailand [Liner 2012], example is visible in Figure 3.6. The Wolf ramp is thought to be a seafloor — it can exhibit a gradational boundary below the seafloor from mud to poorly consolidated sediments and finally lithified material. In order to verify the hypothesis of existence of a Wolf ramp the following steps were undertaken: 1. Classical data processing; result: migrated seismic section, 2. Choosing a representative seismic trace from a section that depicts the features connected with the seabed. 3. Spectral decomposition of the chosen seismic trace; result: time-frequency plot. 4. Finding a synthetic solution that would give a result similar as real timefrequency plot To proceed step 4 the velocities of P wave in water and in the seafloor were applied. P wave velocity in water was assumed to be 1500 m/s, whereas in the seafloor 2100 m/s. With those parameters it was possible to test various rampthickness values. The best fit between real and synthetic data was obtained for.

(42) 41 h = 17m. Throughout the analysis other possible explanations were studied. The most interesting comparison was made for the model described above with pure interference of two reflections in the shallow seafloor (no transition zone). The results were similar but they had one different feature: unlike real data, the Wolf ramp model spectral decomposition raveled notches that have quick recovery of energy beyond the notch frequency. That is always true of interference in which frequencies are not lost or attenuated but are phased out locally to form notches. In a Wolf ramp there is a clear decay above the notch and notches are clear and repetitive. Trace. A 1.0. 1800 1850 1900 1950 2000 2050 2100 2150. Time [s]. 1.1 1.2 1.3 1.4 Migration data B. Amplitude. Time [s]. 1.20. -5 0 5. C 1.20. 1.25. 1.25. 1.30. 1.30. 1.35. Frequency [Hz] 10 20 30 40 50 60 70 80 90. 1.35 Trace 1950. TF of trace 1950 (real). Figure 3.6: Offshore migrated field data from Andaman Basin, Thailand. (A) Seismic time section showing the seafoor reflection, (B) Seismic trace, (C) Timefrequency decomposition from 0 to 100 Hz. [Liner 2012].. The theory discussed in this chapter proves that even in the simple case of normal incidence and constant density, frequency dependent reflectivity has long been predicted. Moreover, as was presented briefly in subsection Interference beyond the thin bed, frequency-dependence of amplitude can be attributed to many factors and the main task of the discussion is how to develop tools for quantitative interpretation that unravel competing frequency-dependent phenomena of the local scale..

(43) 42. 3.2. Various aspects of transition zones. The methodology presented above is an elegant solution that is valid for a simple one-layer model. The same results are visible in real data examples. Both numerical solution as well as real data examples will be presented afterwards. Nevertheless, a few more aspects of transition layers have to be discussed. Multiple transition layers In practice, it would be difficult to imagine a well-separated, well-defined transition layer. Such a situation, however distinguishable in seismic data, is rather rare. It is especially preferable to have velocity logs from a well to define the separated transition layer with high fidelity. The theory of transition zones enables to approximate the velocity log by a series of straight-lines segments of the linear velocity change. By doing so, the velocity log can be represented by the multiple transition zone layers [Berryman et al. 1958]. In this article Berryman et al. (1958) presented a theoretical solution to such a complex system of multiple transition layers. The result is a reflection coefficient equation that is a function of frequency, expressed in circles per period, and contrast of velocities between separated layers, the value which is constant for a pair of layers. The equation assumes that density is constant throughout the section and the solution respects Wolf’s assumption/condition that all waves are plane waves and are at normal incidence to the layers considered. The reflection coefficient for a separated layer is given by. R( j − 1) =. rj + Rj (−2iω(zj −zj−1 ))/vj −1 e 1 + rj Rj. (3.10). The equation 3.10 needs some comments. The variables are as follows: 1. Rj − 1 - reflection coefficient for N constant velocity layer where N changes from n = 0 to N − 1, 2. Rj - reflection coefficient of the above (i.e. previous) transition layer, 3. ω - angular frequency expressed in cycles per octave,.

(44) 43 4. zj - distance from the start of the zero coordination point, understood as the j t h layer thickness, 5. vj − 1 - j t h velocity value rj depends on velocity contrast:. rj =. vj−1 − vj vj−1 + vj. (3.11). and Rj is a function of frequency and scaling factor that represents velocity change between the layers Careful derivation of equation 3.10 is given in Chapter 9 Definitions. Equation 3.10 can be solved recursively for all defined transition layers. The results of solutions to a single transition layer are consistent with the Wolf solution [Wolf 1937], and for the model values both equations 3.10 and 3.5 yield the same outcome ([Berryman et al. 1958], (p.229, Fig. 5)). The result of computation is the reflection response of a system of transition layers of variable velocity but constant density. Originally [Berryman et al. 1958], the function was named a transfer function, but it must be taken into account that in seismic exploration we do not deal with a transmission experiment but reflection experiment conducted over a layered geological body and hence isolated transmission signals are not recorded [Mateeva 2003]. In the light of the above statement, a transfer function computed in the manner described above cannot be compared with the seismic section in terms of absolute similarity i.e. one to one. Processing procedures that may have impact on transition zone analysis Transition zones, so-called Wolf ramps find application in various aspects of seismic interpretation. such as permafrost analysis [Justice & Zuba 1986] and carbon dioxide saturated transition zones [Gómez & Ravazzoli 2012]. These applications lead to reveal very interesting aspects of transition layers. Nevertheless, not all assumptions took into account the linear velocity changes in the transition zone [Gómez & Ravazzoli 2012]. Even though the concept of transition zones was introduced in 1937, the usual study of seismic reflections has been limited to those.

(45) 44 from sharply defined contrasts in acoustic impedance for a long time. However, there exist a wide variety of examples that differ from the sharp boundary and show the frequency-selective nature of attenuation and phase distortion that may lead to pattern recognition in order to identify such reflections [Justice & Zuba 1986]. The corresponding formulas for the reflection coefficient in time and frequency domain are given by equations 3.12 and 3.13 respectively. S(τ ) =. ln(V1 /V0 ) w(τ ), 2. (3.12). where V1 and V0 are velocities at the bottom and at the top of the transition layer, τ expresses the temporal position of the wavefiled within the transition zone.. R(f ) =. ln(V1 /V0 ) sin πf T , 2 πf T. (3.13). where T is the total temporal thickness of the transition zone, expressed in milliseconds, unlike h that is expressed in meters. The transition zone reflection is prone to modifications of amplitude of no less than a factor of 1/2ln(V − 1/V0 ) if V1 < V0 . In this case, a polarity reversal of amplitude will occur provided that the wavelet is band-limited to the frequency range comparable with the temporal thickness of the transition zone (− T1 ≤ f ≤ T1 ) [Justice & Zuba 1986]. Wave-shape modifications and frequency-selective attenuation will occur if the wavelet has frequency components outside this range. Nevertheless, if the velocity ratio is close to one the reflected wavelet may still not be detectable. Another important observation depicted by the study of permafrost [Justice & Zuba 1986] is that if the time-width of the transition zone is large compared with the time-width of the wavelet, due to the integration effect there will be a reflection at each end of the transition zone, connected by the constant zone where trace amplitude is just the integral of the wavelet. In this context one significant remark arises. The deconvolution process which uses as an operator the very wavelet that is present in seismic data may severely distort the characteristics of the transition zone and make it invisible on seismic section. The results of these studies are expressed by the relationship in Figure 7.3 which.

(46) 45 shows a seismic response obtained by the convolution of wavelets of the specific frequencies with a 50-meters thick transition zone. If sufficiently low frequencies are absent from the wavelet, the transition zone can be seismically invisible, i.e. response obtained for 60 Hz has very low amplitude. Although such wavelets may be useful, they may represent the poorest choice for delineating transition zones. The results of numerical computation that are presented in Chapter 7 Results, proved that the analytical solution is valid and effects discussed above occur. Variations of transition zone reflections with offset Moreover it was proven that there exists the dependence of the transition zone response due to a different angle of incidence. The said modification may manifest in a large increase in reflection parameter or/and wave-shape modification due to the frequency dependance of the reflection coefficient with the offset [Gupta 1966]. These two effects exert a lot of pressure on a stacking process i.e. only meticulous stacking process may sustain the presence of the transition zone (see [Justice & Zuba 1986] Fig. 4 Justie which is reproduced in 3.7).. B CDP. VELOCITY (Pwave), m/s. 2 000. WAVELET. A. CHANNEL 3 500. 4. 8. 12 16. 0.00. 0.10. 0.20. 0.30. 0.50. TIME [s]. TIME [s]. 0.40. 0.60. 0.70. 0.80. 0.90. 1.00. 1.10. Figure 3.7: Amplitude versus offset relationship for a transition zone. (A) Velocity model, (B) Amplitude increase with increasing offset [Justice & Zuba 1986]..

(47) 46 Discussion Generally, every model of a layer that is not homogeneous in terms of velocity or/and density can be referred as a transition layer. Changes in these two parameters can significantly modify the seismic response. The models described above , which consist of one or multiple transition layers are most simplistic, and real data analysis shows that more complicated cases can be encounter. In Figure 3.8, different models of transition zones are presented. The first three cases presented in E — G to the transition zones that correspond to the classical case, they differ only in the transition zone thickness. Others, however, differ not only in the transition zone thickness, but also in the character of the velocity function at the transition zone boundary, i.e. the velocity function has one (at its top) or two (at the top and bottom) discontinuities. The equations presented herein show solutions to the most simplistic case of a transition zone. It should be underlined that in real data more complex cases occur, for instance those presented in Figure 3.8. It is common to assume that the reflection amplitude should look elegant: it should be as symmetrical as possible. Processing procedures often aim to modulate reflection amplitudes in order to make them more interpreter-friendly. This may not always be desirable — nature of many reflections patterns may be different, e.g. not looking as zero-phase reflections, which is perfectly true. By making such a reflection more symmetrical and zero-phase we may loose theirs real character and misinterpret true geological signature of the data. The theory presented in this chapter is only a rough approximation of all phenomena that are encountered and have influence on seismic data. Conclusions To bring this Chapter to a close, I summarise the main points: 1. Reflection from the transition zone — a Wolf ramp is a function dependent on frequency. It suffers wave-shape and amplitude alternation as well as phaseshift. 2. Reflection from the broad (thick) transition zone, the one whose temporal limits are much larger than wavelet time-width will result in two separate.

(48) 47 A. B. C. D. E. F. G. H. I. J. K. L. M. N. Model. Seismic response. Figure 3.8: Seismic response of a layer modified by the seismic impedance changes: (A), (B) thick layer models, (C), (D) thin bed models, (E) — (N) transition zone models, [Kasina 1998].. reflections at the top and bottom (if the wavelet is zero-phase). The opposite relation may cause the relative or complete invisibility of a Wolf ramp on a seismic section. 3. Modification of the amplitude and wave shape is offset-dependent, i.e. the amplitude changes with the changes in angle of incidence. For further offset, an increase in amplitude can be observed. This effect is more likely to occur for broad transition zones. 4. Deconvolution and the process of stacking have a significant influence on the visibility of potential transition zones..