A GENERALIZED KIRCHHOFF-WKBJ DEPTH MIGRATION THEORY FOR MULTI-OFFSET SEISMIC REFLECTION DATA: REFLECTIVITY MODEL CONSTRUCTION BY WAVEFIELD IMAGING AND AMPLITUDE ESTIMATION By David Edward Lumley B. Sc. (Hon. Geophysics) University of British Columbia A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T O F THE REQUIREMENTS FOR T H E DEGREE OF MASTER OF SCIENCE in THE FACULTY O F GRADUATE STUDIES DEPARTMENT O F GEOPHYSICS A N D ASTRONOMY We accept this thesis as conforming to the required standard THE UNIVERSITY O F BRITISH COLUMBIA October, 1989 © David Edward Lumley, 1989 In presenting this degree at the thesis in University of partial fulfilment of of this thesis for department or by his requirements British Columbia, I agree that the freely available for reference and study. I further copying the her representatives. an advanced Library shall make it agree that permission for extensive scholarly purposes may be granted or for It is by the understood that head of copying my or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada Date DE-6 (2/88) O c - T b p g g . \o t figcj Abstract This thesis embodies a mathematical, physical, and quantitative investigation into the imaging and amplitude estimation of subsurface earth reflectivity structure within the framework of pre stack wave-equation depth migration of multi-offset seismic reflection data. Analysis is performed on five prestack depth migration reflectivity "imaging conditions" with respect to image quality and quantitative accuracy of recovered reflectivity amplitudes. A new computationally efficient and stable prestack depth migration imaging method is proposed which is based upon a geometric approximation to the theoretically correct, but unstable, "dynamic" imaging condition. The "geometric" imaging condition has the desirable property of true-amplitude reflectivity recovery in regions of both 1-D and 2-D velocity variation, while fully retaining and optimizing the favorable imaging characteristics of current, nontrue amplitude formulations. The currently predominant "crosscorrelation" and "excitationtime" migration imaging methods are shown to possess significantly less accurate imaging and amplitude-recovery characteristics relative to the proposed geometric migration. The respective signal-to-noise recovery of their imaged amplitudes deteriorates approximately linearly (excitation-time) and quadratically (crosscorrelation) with depth. As a necessary prerequisite to the imaging analysis, a true-amplitude prestack depth migration equation is derived which appears to be new to the literature. This result is obtained in the form of a 2.5-D farfield Kirchhoff integral solution to the acoustic wave equation, after the application of a dynamic imaging condition to the reconstructed upgoing and downgoing ii wavefields. This solution is in harmony with zeroth order asymptotic ray theory (ART) assumptions, and depends upon WKBJ Green's functions which can be numerically evaluated for arbitrary migration models by raytracing methods. A new and "generalized" Kirchhoff prestack depth migration equation is subsequently obtained by the introduction of a weighting function into the true-amplitude migration integral. The weight is a function of both the reconstructed upgoing and downgoing wavefields, and is determined analytically by a mathematical application of each specific reflectivity imaging condition. This generalized equation is significant in that it provides a common mathematical, physical and computational basis for the comparative analytical and quantitative analysis of reflectivity image quality and amplitude recovery among current prestack migration philosophies and variants of those migration themes. In addition, three ancillary research objectives are achieved. The first achievement is the development of a Kirchhoff prestack depth migration computer algorithm to implement the generalized imaging of surface-recorded seismic reflection data. This algorithm can be readily modified to perform seismic wavefield imaging for other recording geometries such as crossborehole or vertical seismic profiling, and may be suitable to non-seismic applications such as the imaging of electromagnetic wavefields and satellite-acquired synthetic aperture radar data. The second result is the development of a fast two-point raytracing computer algorithm which provides accurate computation of a subsurface grid of traveltimes and 1.5-D zeroth order ART amplitudes in a 1-D acoustic medium. This algorithm is useful for subsurface wavefield reconstruction and imaging, and for inversion applications such as geotomography. The third objective is the detailed quantitative examination of migration imaging quality and true relative-amplitude normal-incidence reflectivity recovery from numerically migrated depth images. This is achieved successfully in an extensive 2.5-D synthetic data analysis, iii using a challenging 2-D structural model, synthetic multifold reflection seismogram shot gathers, and the numerical imaging and modelling algorithms developed as part of this thesis research. iv Table of Contents Abstract ii List of Tables xii List of Figures xv Acknowledgment xx 1 Introduction 1 2 The Fundamentals of Kirchhoff Migration Theory 6 2.1 Introduction 6 2.2 An Overview 2.3 The Scalar Wave Equation 2.4 The Green's Function Problem 2.5 Physical Interpretation of the Green's Function Solutions 7 . 8 13 15 3 4 2.6 The Kirchhoff Integral Solution to the Scalar Wave Equation 18 2.7 Defining the ZSR Reflectivity Imaging Condition 25 2.8 The Conventional Kirchhoff Post-Stack Migration Equation 27 2.9 Interpreting the Kirchhoff Migration Equation 30 A True-Amplitude Kirchhoff-WKBJ Migration Theory 33 3.1 Introduction 33 3.2 A True-Amplitude Migration Philosophy 34 3.3 First Principles 36 3.4 The Upgoing Wavefield Solution 42 3.5 The Downgoing Wavefield Solution 49 3.6 Applying the Reflectivity Imaging-Inversion Condition 50 3.7 The True-Amplitude Kirchhoff Migration Solution 51 The Migration Reflectivity Imaging Conditions 55 4.1 Introduction 55 4.2 The Generalized Kirchhoff Migration Solution 56 vi 4.3 The Dynamic Imaging Condition 57 4.4 The Crosscorrelation Imaging Condition 58 4.5 The Excitation-Time Imaging Condition 59 4.6 The Geometric Imaging Condition 60 4.7 The Kinematic Imaging Condition 61 5 Mathematical Analysis of Reflectivity Amplitude Recovery 5.1 Discussion on "Velocity Models" 63 63 5.2 Geometric Migration 66 5.3 Kinematic Migration 71 5.4 Excitation-Time Migration 75 5.5 Crosscorrelation Migration 78 5.6 Summary 81 6 Zeroth Order Asymptotic Ray Theory 83 6.1 Introduction 83 6.2 The Ray Series Expansion 84 vii 7 6.3 Solving the Eikonal Equation 87 6.4 Solving the Transport Equation 89 6.4.1 Transforming to Ray-Centred Coordinates 90 6.4.2 A General Solution to the Transport Equation 92 6.4.3 Evaluating the Source Term 93 6.4.4 Evaluating the Jacobian Term 94 6.4.5 Energy Partitioning 98 6.4.6 Q-Attenuation 100 6.4.7 Transport Equation Amplitude Summary 100 1-D Two-Point Raytracing 103 7.1 Introduction 103 7.2 The Ray Kinematics 104 7.2.1 The Newton-Raphson Method 107 7.2.2 A Discrete Layer Formulation 109 7.2.3 Aspects of Convergence 110 viii 7.3 The Ray Amplitudes 113 7.4 Summary 116 8 A 2-Dimensional Test Model 117 8.1 The "Arc" Test Model 117 8.2 The Synthetic Data 120 8.3 ° The Migration Model 122 9 Raytracing the WKBJ Green's Functions 127 9.1 Introduction 127 9.2 The Raytraced Imaging Times 128 9.3 The Raytraced WKBJ Amplitudes 130 9.4 Computational Requirements 138 10 Quantitative Analysis of Migration Data Examples 140 10.1 Introduction 140 10.2 Single Trace Depth Migration Responses 141 10.2.1 Dynamic Migration 141 ix 10.2.2 Geometric Migration 144 10.2.3 Kinematic Migration 146 10.2.4 Excitation-Time Migration 146 10.2.5 Crosscorrelation Migration 149 10.3 Single Shot Gather Depth Migration 149 10.3.1 Dynamic Migration 151 10.3.2 Geometric Migration 153 10.3.3 Kinematic Migration 156 10.3.4 Excitation-Time Migration 158 10.3.5 Crosscorrelation Migration 158 10.4 Multifold Prestack Depth Migration 161 10.4.1 Dynamic Migration 161 10.4.2 Geometric Migration 170 10.4.3 Kinematic Migration 173 10.4.4 Excitation-Time Migration 177 10.4.5 Crosscorrelation Migration 182 x 10.5 Geometric Migration Stability in the Presence of Noise 186 10.6 Geometric Migration Stability in the Presence of Velocity Model Errors . . . 191 10.7 Computational Requirements 197 10.8 Summary of the Migration Imaging Condition Quantitative Analysis 199 11 Conclusion 203 References 206 xi List of Tables 8.1 The 2-D Arc Test Model Material Property Specification 119 8.2 The 1-D Migration Model Material Property Specification 125 9.1 139 CPU Times for Green's Function Raytracing Evaluation 10.1 The True Absolute and Normalized 1-D Horizontal-Layer Normal-Incidence Acoustic Reflection Coefficients as Calculated from the Arc Test Model. . . 164 10.2 Normalized 1-D Horizontal-Layer Reflection Coefficients from the True Arc Test Model and the Dynamic Migration Depth Section 165 10.3 The True Absolute and Normalized 2-D Dipping Arc Segment Normal-Incidence Acoustic Reflection Coefficients as Calculated from the Arc Test Model. 166 10.4 Normalized 2-D Dipping Arc Segment Reflection Coefficients from the True Arc Test Model and the Dynamic Migration Depth Section 167 10.5 Normalized 1-D Horizontal-Layer Reflection Coefficients from the True Arc Test Model and the Geometric Migration Depth Section xii 172 10.6 Normalized 2-D Dipping Arc Segment Reflection Coefficients from the True Arc Test Model and the Geometric Migration Depth Section 173 10.7 Normalized 1-D Horizontal-Layer Reflection Coefficients from the True Arc Test Model and the Kinematic Migration Depth Section 175 10.8 Normalized 2-D Dipping Arc Segment Reflection Coefficients from the True Arc Test Model and the Kinematic Migration Depth Section 176 10.9 Normalized 1-D Horizontal-Layer Reflection Coefficients from the True Arc Test Model and the Excitation-Time Migration Depth Section 179 lO.lOAbsolute and Normalized Normal-Incidence Geometric Divergence Correction Factors F Calculated from the True Arc Test Model n 179 10.11 Normalized 2-D Dipping Arc Segment Reflection Coefficients from the True Arc Test Model and the Excitation-Time Migration Depth Section 181 10.12Normalized 1-D Horizontal Layer Reflection Coefficients from the True Arc Test Model and the Crosscorrelation Migration Depth Section 184 10.13Normalized 2-D Dipping Arc Segment Reflection Coefficients from the True Arc Test Model and the Crosscorrelation Migration Depth Section 185 10.14Normalized 1-D Horizontal Layer Reflection Coefficients from the True Arc Test Model and the Clean and Noisy Geometric Migration Depth Traces. . . 191 10.15Migration Model Velocities and Depths from the True Arc Test Model and the Velocity Semblance Analysis of a Noisy Shot Gather xiii 193 10.16 Normalized 1-D Horizontal Layer Reflection Coefficients from the True Arc Test Model and the Geometric Migration Depth Traces xiv List of Figures 2.1 The Geometry for the Scalar Wave Equation 11 2.2 The Geometry for the Green's Function Solution 16 2.3 The Raypath Geometry for ZSR Reflections 26 2.4 The Exploding Reflector Model 28 3.1 An Idealized Shot Gather Geometry 37 3.2 An Idealized Reflection Raypath Geometry 39 7.1 A Two-Point Raypath in a 1-D Horizontally Layered Medium 105 7.2 A Graphical Representation of the Newton-Raphson Method for Finding the Ray Parameterp* Ill 8.1 The 2-D "Arc" Test Model 118 8.2 A Representative Shot Gather Raytraced at x = 1200 m along the Arc Test s Model 121 xv 8.3 The Pre-Processed Shot Gather at x - 1200 m along the Arc Test Model.. . 123 s 8.4 The 1-D Migration Model Obtained from the 2-D Arc Test Model 124 9.1 A 2-D Subsurface Map of the Raytraced Green's Function Reflectivity Imaging Times t(x0) 129 9.2 A 2-D Subsurface Map of the Raytraced Green's Function Geometric Divergence Amplitude Factor An (x0) 131 9.3 A 2-D Subsurface Map of the Raytraced Green's Function Energy Partitioning Amplitude Factor Ap(x0) 133 9.4 A 2-D Subsurface Map of the Raytraced Surface-Incident Obliquity Amplitude Factor cos0r(x<>) 134 9.5 A 2-D Subsurface Map of the Raytraced WKBJ Green's Function Amplitude Factor A(x0) 136 9.6 A 2-D Subsurface Map of the Raytraced Obliquity-Modified WKBJ Green's Function Amplitude Factor A(x0) cos 9r(x0) 137 10.1 Single Trace Depth Migration Response to the Dynamic Reflectivity Imaging Condition 142 10.2 Single Trace Depth Migration Response to the Geometric Reflectivity Imaging Condition 145 xvi 10.3 Single Trace Depth Migration Response to the Kinematic Reflectivity Imaging Condition 147 10.4 Single Trace Depth Migration Response to the Excitation-Time Reflectivity Imaging Condition 148 10.5 Single Trace Depth Migration Response to the Crosscorrelation Reflectivity Imaging Condition 150 10.6 A Single Shot Gather Depth Migration using the Dynamic Reflectivity Imaging Condition 152 10.7 Ray Diagram Illumination Pattern of Reflected Arrivals for a Shot Gather at xs = 1200 m 154 10.8 A Single Shot Gather Depth Migration using the Geometric Reflectivity Imaging Condition 155 10.9 A Single Shot Gather Depth Migration using the Kinematic Reflectivity Imaging Condition 157 10.10A Single Shot Gather Depth Migration using the Excitation-Time Reflectivity Imaging Condition 159 10.11 A Single Shot Gather Depth Migration using the Crosscorrelation Reflectivity Imaging Condition 160 10.12Full-Fold Kirchhoff-WKBJ Prestack Depth Migration using the Dynamic Reflectivity Imaging Condition 162 xvii 10.13Migration Depth Trace Stability Analysis of the Dynamic Reflectivity Imag2 ing Condition for Variable Stabilization Parameter Values e 169 10.14Full-Fold Kirchhoff-WKBJ Prestack Depth Migration using the Geometric Reflectivity Imaging Condition 171 10.15Full-Fold Kirchhoff-WKBJ Prestack Depth Migration using the Kinematic Reflectivity Imaging Condition 174 10.16Full-Fold Kirchhoff-WKBJ Prestack Depth Migration using the ExcitationTime Reflectivity Imaging Condition 178 10.17Full-Fold Kirchhoff-WKBJ Prestack Depth Migration using the Crosscorrelation Reflectivity Imaging Condition 10.18The Muted Shot Gather at x = 1200 m with 10% Additive Noise s 183 187 10.19 Geometric Migration of a Noisy Shot Gather 188 10.20Normalized Depth Traces Corresponding to the True Normal-Incidence Acoustic Reflection Coefficients of the Arc Test Model, and the Geometric Shot Gather Migration of Clean and Noisy Synthetic Data 190 lO.HGeometric Migration of a Noisy Shot Gather with an Erroneous Migration Velocity Model 194 10.22Normalized Depth Traces Corresponding to the True Normal-Incidence Acoustic Reflection Coefficients of the Arc Test Model, and the Geometric Shot Gather Migrations 195 xviii 10.23Full-Fold Kirchhoff-WKBJ Prestack Depth Migration Traces for each of the Five Reflectivity Imaging Conditions xix 200 Acknowledgment I would like to take this opportunity to thank those who contributed vital research funding and service support during this thesis preparation. My work was generously supported by a research grant from Esso Resources Canada Limited, and operating grant 5-84270 to Prof. Doug Oldenburg from the National Science and Engineering Research Council of Canada. I am very grateful to the Amoco Canada Petroleum Company Ltd., and the Canadian Society of Exploration Geophysicists for their significant level of scholarship support. In addition, I thank Marvin Bloomquist, Bob Keys, Keh Pann, and Mark Willis of Mobil Research and Development Corp., Dallas, for the opportunity to work in an industrial research environment, and for their constructive criticism and advice during the initial stages of this project. I would like to thank the many who contributed to an environment of intellectual stimulation during my MSc program and thesis research. My acquaintances at the Dept. of Geophysics and Astronomy made UBC a very enjoyable place in which to pursue this degree. A special note of utmost regard goes to Professor Doug Oldenburg for his ongoing inspiration and faith in my ability. Finally, I thank my family, who kindled within me an awe of the physical phenomena that pervade our existence, and the endless curiosity to explore and understand them. And most important of all, I thank my wife Glynis, and my daughter Brianne, for their love, support, and sharing of a very pleasant life indeed. xx Chapter 1 Introduction An increase in the application of geophysical inversion theory to quantitative seismic exploration research problems is currently in evidence. This rise in interest is exemplified by the impedance inversion of surface seismograms (Oldenburg et al., 1986), cross-borehole tomographic inversion methods (Bregman et al., 1989), and the application of generalized linear inversion techniques to amplitude-offset data (Hampson, 1989), for example. Recent advances in the theoretical and implementational aspects of seismic wave-equation inversion (Tarantola, 1988, and Crase et al., 1988) are heightening interest in the future prospects of seismic inversion theory. Quantitative impedance models, produced as direct results of these seismic inversion techniques, can be useful in collaboration with other geophysical data to better determine the structure and distribution of subsurface geophysical properties. For example, such impedance models can be used with direct rock property data, obtained from borehole cores, logs and other geophysical/geologic control, to assist in the prediction of detailed rock and interstitial fluid distributions in the subsurface. Seismic inversion can be useful for the refinement of the spatial delineation of reservoir rock structure, and the physical parameter estimates of the fluid and gaseous constituents in their pore space. Improvement in the quality of these inversions would be of major consequence in the exploration, production, or monitoring of 1 subsurface hydrologic, hydrothermal, or hydrocarbon resources, for example. The research documented in this thesis forms a part of my larger global research goal to develop a multi-dimensional elastic inversion theory and implementation method for seismic reflection data. My intention has been to approach such an inversion by developing an iterative hybrid method of elastic parameter model construction consisting of successive steps of: seismic wavefield reconstruction, reflectivity depth imaging and amplitude-estimation, reflectivity inversion to obtain the elastic parameter model perturbations, forward modelling of synthetic seismograms from the updated elastic model, and optimization of the "fit" between the observed and synthetic data. Some of the fundamental inversion theory related to this problem has been discussed by Tarantola (1986), for example. Although not presented in this thesis, I have done a considerable amount of research work regarding 3-D acoustic seismic inversion theory. I have also implemented a 2-D iterative acoustic seismic inversion using the Born Approximation of inverse scattering theory, with very encouraging results. However, this preliminary research exposes several obstacles which confront the quantitative construction of geophysically reliable impedance models from multifold seismic reflection data. These include the critical task of recovering reliable true-amplitude reflectivity information during the subsurface reflectivity imaging (migration) process. Amplitude accuracy in the reflectivity imaging step is an essential prerequisite to thefinalcomponent of the inversion process, where the preliminary reflectivity and subsequent impedance distribution values are computed. By approaching seismic inversion within a true-amplitude prestack depth migration theory, vital wavefield traveltime, amplitude, and raypath information, destroyed in a post-stack approach, may be used advantageously to yield enhanced subsurface geophysical parameter estimations. 2 The combined requirement of accurate prestack reflectivity imaging and true-amplitude recovery for the generalized wave-equation inversion of seismic reflection data motivates the primary goal of this thesis research: reflectivity model construction. Philosophically, the problem of reflectivity model construction can be decomposed into two related problems: reflectivity imaging, and reflectivity amplitude estimation. Traditionally, seismic wave- equation migration research has focussed on the imaging aspect of reflectivity model construction. Migration has proved to be valuable for qualitative analysis of seismic reflection data, especially in the area of seismic structural interpretation. However, my present research is concerned with the quantitative, rather than qualitative, analysis and inversion of seismic data. Hence, this thesis addresses the more complete reflectivity model construction problem, by developing methods which perform simultaneous reflectivity imaging and amplitude-estimation, within the context of seismic wavefield reconstruction. On this basis,fivedistinct reflectivity imaging conditions are examined in a comprehensive mathematical, physical, and quantitative analysis within a common Kirchhoff-WKBJ spacetime integral solution domain. The fundamentals of Kirchhoff post-stack migration theory arefirstexamined, in order to introduce the mathematical constructs and physical intuition necessary to ascend to the more complete Kirchhoff-WKBJ prestack migration theory. Following this, a true-amplitude Kirchhoff-WKBJ prestack depth migration equation is derived, which appears to be new to the literature. This result is obtained in the form of a 2.5-D farfield Kirchhoff integral solution to the acoustic wave equation, after the application of a "dynamic" imaging condition to the reconstructed upgoing and downgoing wavefields. The solution is dependent upon WKBJ Green's functions, which can be numerically evaluated for heterogeneous migration models by raytracing methods, and is therefore fundamentally related to the eikonal and transport equation solutions of zeroth order asymptotic seismic ray theory (ART). 3 A new and "generalized" Kirchhoff-WKBJ prestack depth migration equation is subsequendy obtained by the introduction of a weighting function into the true-amplitude migration integral. The weight is a function of both the reconstructed upgoing and downgoing wavefields, and is determined analytically by a mathematical application of each specific reflectivity imaging condition. This generalized equation is significant in that it provides a common mathematical, physical and computational basis for the comparative analytical and quantitative analysis of reflectivity image quality and amplitude accuracy among current prestack migration philosophies and variants of those migration themes. The generalized KirchhoffWKBJ equation is then subjected to an extensive mathematical and physical analysis, in order to determine the leading-order behaviour of the amplitude-estimation properties inherent in each prestack depth migration reflectivity imaging condition. The relationship between the WKBJ Green's functions of the generalized migration equation and zeroth order asymptotic ray theory is discussed. The discussion of ART includes the traveltime and WKBJ amplitude solutions to the eikonal and transport equations, and leads to the development of an iterative 1.5-D two-point raytracing computer algorithm. The raytracing algorithm is subsequently implemented to numerically evaluate the WKBJ Green's function components of the generalized Kirchhoff-WKBJ migration equation. Finally, a computer algorithm is implemented to perform the generalized Kirchhoff-WKBJ prestack depth migration imaging and amplitude-estimation of subsurface reflectivity, given a set of observed seismic reflection data. A detailed quantitative examination of migration imaging quality and true relative-amplitude normal-incidence acoustic reflectivity recovery from numerically migrated depth images is conducted. This is achieved successfully in an extensive 2.5-D synthetic data analysis, using a challenging 2-D structural model, synthetic 4 multifold reflection seismogram shot gathers, and the numerical imaging and modelling algorithms developed as part of this thesis research. As afinalremark, the Kirchhoff-WKBJ wavefield imaging and amplitude recovery methods presented in this thesis may also be applicable to other seismic wavefield imaging and inversion problems. For example, I have recently entertained an idea regarding a waveequation "migration" offirstarrival waveforms in seismic data, including the "head-wave" arrivals corresponding to critically refracted, or diving-wave energy. I have implemented this for a simple test case, using a modified version of the Kirchhoff-WKBJ migration algorithm developed in this thesis, with extremely encouraging results. Such a procedure could be useful as a new method for constructing a subsurface structural image or parametric model, using a wave-theoretic inversion of seismicfirstarrival waveforms. This mayfindapplication in the quantitative analysis of seismic refraction survey data, and also for VSP and crossborehole data inversions directed toward the enhanced delineation of subsurface reservior material property distributions. Furthermore, the Kirchhoff-WKBJ research presented in this thesis may lead to the possibility of somewhat non-traditional geophysical applications of migration imaging and inversion theory. This includes the migration of electromagnetic (EM) wavefield data, as discussed by Zhdanov and Frenkel (1983), for example, in which it may be possible to apply the Kirchhoff-WKBJ formalism to the EM wave equation and employ the use of EM diffusion, rather than seismic wave propagation Green's functions. Other possible applications of my thesis research arise from the increasing potential for interaction between the seismic and the remote sensing communities, regarding the application of seismic wave-theoretic imaging and inversion theory to remotely sensed satellite data, such as that of synthetic aperture radar for example, as discussed by Rocca (1987). 5 Chapter 2 The Fundamentals of Kirchhoff Migration Theory 2.1 Introduction The purpose of this chapter is to introduce and examine the fundamental theory leading to the conventional Kirchhoff integral solution of the scalar wave equation, and its application to the reflection seismic migration problem. The problem of reconstructing the subsurface seismic reflection wavefield from the observed boundary data is posed in a rigorous mathematical framework. The auxiliary problem for the Green's function impulse response to the scalar wave operator is presented, and the Green's function solution for a constant-velocity halfspace is examined from a physically intuitive viewpoint. Using the Green's function solution, the Kirchhoff integral solution to the scalar wave equation is obtained by an application of Green's Second Identity. Subsequently, the zero-offset source-receiver (ZSR) reflectivity imaging condition is developed and applied to the Kirchhoff integral solution to obtain a solution to the conventional post-stack seismic migration problem. Finally, the Kirchhoff post-stack migration equations are physically interpreted to enhance the understanding of the seismic wave-equation migration procedure. This chapter serves as an introduction to the more complete Kirchhoff-WKBJ prestack depth migration theory that forms the core component of this thesis research. By understanding 6 the conventional Kirchhoff theory presented in this chapter, the reader can more easily transcend to, and acquire an appreciation of, the Kirchhoff-WKBJ theory derived later. In particular, the conventional Kirchhoff theory assumes piece-wise constant-velocity media, and is traditionally developed from a post-stack ZSR perspective. In contrast, the KirchhoffWKBJ theory derived later, admits heterogeneous media and prestack multi-offset wavefield reconstruction by the use of WKBJ ray-theoretic Green's functions. Of primary importance, the Kirchhoff-WKBJ theory allows for the recovery of true-amplitude reflectivity, and offers a generalized framework for mathematical, physical, and quantitative analysis of various reflectivity image construction and amplitude recovery methods. 2.2 An Overview Migration is a two-part problem requiring: (1) the reconstruction of the subsurface upgoing (reflected) wavefield using the observed surface seismic data as a boundary condition to the wave equation, and (2) the physical definition and mathematical application of a reflectivity imaging condition as a function of the reconstructed wavefields. Kirchhoff migration is so named because the method used for reconstructing the upgoing wavefield throughout the subsurface medium is based upon the Kirchhoff integral solution to the scalar wave equation. Two early papers that relate the Kirchhoff solution to wave equation migration are given by French (1975) and with somewhat more generality, Schneider (1978). An exhaustive treatment of the subject is given by Berkhout (1984). The Kirchhoff integral solution is a boundary value method which allows the wavefield within a given volume V to be obtained by a space-time integration of the boundary values 7 that are observed on the surface that bounds the volume V. In contrast, the f-k mi- gration methods of Stolt (1978) and Gazdag (1978) involve a frequency domain phase-shift downward continuation of the Fourier-transformed surface data to obtain the upgoing wavefield as a function of the subsurface coordinates. The w-x migration method of Claerbout (1971) employsfinitedifference operators to approximate the scalar wave equation, and thus reconstructs the subsurface upgoing wavefield by wave-extrapolation of the observed surface wavefield. The mathematics involved in the Kirchhoff solution are the most difficult of any of the migration solutions. However, the end result is the most intuitive from a physical point of view. Finally, once the upgoing wavefield has been reconstructed, a reflectivity imaging condition is defined and applied independent of the particular migration wavefield reconstruction method used. This chapter examines the Kirchhoff theory that forms the framework for conventional poststack Kirchhoff migration, in contrast with the prestack Kirchhoff-WKBJ migration theory developed in the body of this thesis. This section attempts to clearly derive mathematical results which are presented as starting, or intermediate points in the two classic post-stack Kirchhoff migration papers by French and Schneider. 2.3 The Scalar Wave Equation We begin by assuming that the observed seismic reflection data represent the surface measurement along d'V of an upgoing pressure wavefield L7(x, t) propagating through the Earth's volume V. In general, the surface data are recorded in shot gathers \P(xr, t; xs) for given shot locations xs each containing many receiver locations xr. However, for the purpose of our subsequent post-stack migration analysis, the shot gathers can be regathered into a common 8 midpoint (CMP) domain, where the midpoint coordinate x m is defined as x m = (xs + xr)/2. The traces within each CMP gather can be approximately corrected for the normal move-out (NMO) of reflected events due to increasing source-receiver offset, and then summed, or stacked, together. Each stacked trace resulting from a CMP gather approximates the zero-offset source-receiver (ZSR) trace that would have been recorded at that surface CMP location. The CMP-ZSR approximation is exact for the situation in which the subsurface background velocity structure and reflectivity structure are 1-D. Hence, the "observed" upgoing reflected wavefield that we desire to post-stack migrate is really a pre-processed, CMP-stacked modification of the originally observed shot gathers, and approximates the upgoing wavefield that would have been observed from many hypothetical simultaneous ZSR experiments conducted over the surface dl/ of a 1-D earth in the volume V. We begin the mathematical analysis by assuming the NMO-corrected CMP-stacked data *F approximate a surface-measured ZSR upgoing reflected wavefield U(x, t) which is governed by the homogeneous, constant-velocity scalar wave equation: V U(x, t) - ^d U{% 0 = 0; x e V, 2 tt (2.1) satisfying the spatial boundary conditions: U(x,t) = ¥(x,t) x 2 ; xedV, lim ||L7(x,0|| =0 ; x e dl/, 9 (2.2) (2.3) and satisfying the temporal boundary conditions: U(x, t) = 0 ; t < 0, / \\U(x,t)\\ dt 2 < oo Jo (2.4) (2.5) ; *e(0,«.). The geometry associated with the system of equations (2.1) - (2.5) is given in Figure 2.1. The temporal boundary condition given by (2.4) ensures that the wave potential U must be causal. Boundary condition (2.5) ensures that the wave potential U must be stable (i.e., U(x, t) must be square-integrable over all time). The combined conditions of causality and stability define a physically realizable quantity. This is a most desirable property to require of any possible scalar wave potential solution U(x, t) satisfying (2.1). The spatial boundary condition (2.2) specifies the value of U(x, t) everywhere on the surface ~dV that bounds the solution domain U. The enclosing surface d'V can be decomposed into r two surfaces: a lower hemispherical surface 9l and the upper planar areal surface A upon which it is assumed we record seismic data. Boundary condition (2.3) is an expression of the requirement for conservation of energy, and is intimately related to the Sommerfeld radiation condition (Zauderer, 1983, p.432). In basic terms, as the radius ^,of the hemispherical surface Of increases, the surface area of 9i increases in proportion to the square of radius, causing the energy density flux through a unit area on the expanding surface 9f to 1 decrease in inverse proportion to the square of the radius,'Sz .In the limit as ^approaches infinity the energy densityfluxapproaches zero as stated in (2.3). We may notice that conditions (2.2) and (2.3) are Dirichlet boundary conditions since, in general, any spatial boundary conditions can be classified as Dirichlet, Neumman, or Cauchy: 10 air Figure 2.1: The Geometry for the Scalar Wave Equation. The solution domain for the scalar wave equation is a hemispherical volume V of radius % . The spatial boundary conditions for the scalar wave equation are specified on the surface dV, which bounds IS, and consists of an upper planar surface A and a lower hemispherical surface !H. 11 Dirichlet: Neumman: dn U(x,t) = ¥N (x,t) ; xzdV, Cauchy: a(x)L7(x, t) + $(x)dn U(x, t) = Y (x, t) ; x e 9V. c The operator 3 is the spatial derivative in the direction h, normal to the surface dV, as RA shown in Figure 2.1. For the seismic migration problem at hand we only have data available on the surface Si and not on 0(. For this reason we will assume that !^is sufficiently large that boundary condition (2.3) applies to the data measured on the hemispherical surface Oi, but not on !A, of the total surface measurements on dV. In other words, we assume that the boundary data contribution from Oi has decayed sufficiently, U(x,t) = 0 ; x e Oi, (2.6) such that boundary condition (2.2) reduces to x U(x,t) = ¥(x,t) ; x e A, where ^(x, t) is, in general, the observed upgoing wavefield measured at the surface (2.7) in shot gathers. However, for the post-stack analysis to follow, ^P(x, t) is assumed to be the 12 approximate ZSR section represented by the NMO-corrected CMP-stacked seismic reflection data. 2.4 The Green's Function Problem Given the partial differential equation system represented by (2.1) and its boundary conditions (2.3), (2.4), (2.5), (2.6), and (2.7), we can set up an auxiliary problem to solve for the impulse response to the linear partial differential wave operator L of (2.1). The auxiliary problem is known as the Green's function problem, and its solution is the Green's function impulse response to the scalar wave operator (Morse and Feshbach, 1953, §7.3). To set upthe auxiliary problem we replace the scalar wave potential U(x, t) with a Green's function G(x, t) in (2.1) and drive the equation with a body forcing function that is a unit impulse Dirac delta function: V <G(x, t) - ^3«G(x, t) = -5(x - x')8(£ -1') ; x, x' e V , 2 c (2.8) where (x', t') is the impulse source coordinate and (x, t) is the observation coordinate. We observe the Green's function impulse response wavefield G at the observation coordinate (x, t) due to the impulsive source -8(x)S(£) placed at the source location (x', t'). The Green's function boundary conditions are the homogeneous versions of the spatial boundary conditions (2.3), (2.6), (2.7), for the primary problem (2.1): G(x, t) = 0 ; x e %, 13 (2.9) and 2 (2.10) lim ||G(x,*)|| = 0 ; x e M. The solution to the problem posed by equations (2.8) - (2.10) can be obtained by Fourier transform methods to solve for the Green's function impulse response to the wave operator in a whole-space, followed by an application of the method of images to match the half-space boundary condition (2.9) on the planar data surface !A. (Morse and Feshbach, 1953, §7.3). The result is: G{x,t,x,t)-Gx+ G2 - — — , where the radial impulse source and image distances are given by 2 2 2 n = y/(x- x') + (y -y') + (z- z') , (2.12) and 2 2 2 r2 = \/(x - x') + (y -y') + (z + z') . (2.13) We remark that G\ is the Green's function solution to the scalar wave operator in a wholespace, and that G2 is the negative polarity image source Green's function of G\. The combination of the whole-space Green's function G\ and its negative polarity image G2 creates the half-space Green's function G which satisfies the wave equation (2.8) and the homogeneous boundary conditions (2.9) and (2.10). To see that the Green's function solution does indeed satisfy the wave equation and its boundary conditions, we reference the coordinate origin to the surface 14 In this case, z = 0 on ft, which makes G\ = -G2 on A, which causes the half-space Green's function G to vanish on the boundary surface ft. as required by the boundary condition (2.9). Also as r -> °°, as is the case on the surface Oi, both G\, G2, and hence G —> 0, as required by boundary condition (2.10). Finally, the Green's functions Gi and G2 are obviously solutions to the scalar wave equation since they have the familiar D'Alembertian functional form, f(r±ct), and the 1/r spherical divergence factor. And by the principle of superposition, G is also a solution to the scalar wave equation, since it is a linear combination of the two wave solutions, Gi and G2. 2.5 Physical Interpretation of the Green's Function Solutions The wholespace Green's functions can be physically interpreted. The wholespace Green's functions G\ and G2 match no boundary conditions alone, other than the energy conservation requirement that they vanish as r — » 00. G\ represents the wholespace response to an impulsive source placed at x', G2 is merely a negative polarity image of G\ positioned by reflecting x' about the plane ft where the boundary condition needs to be satisfied. The Green's function geometry is illustrated in Figure 2.2. The impulse source at x' is ignited at time t = t' and its energy takes a finite amount of time x = ri/c to propagate from the source location x' to the observation location x where the impulse response G is measured. This temporal bookkeeping and impulsive nature of the source is imbedded in the delta function component of the Green's function: 8(t-t'-ri/c). The wave energy is expected to propagate as spherical waves in a constant velocity medium, and hence, G\ has an amplitude component proportional to the geometric divergence of spherical waves in a constant velocity medium: 1/471^. In summary, the wholespace Green's functions 15 image observer Figure 2.2: The Geometry for the Green's Function Solution. An impulsive source is ignited at x', which is subsequendy observed at x. The Green's function G\ corresponds to the direct arrival, which travels the distance r\ from x' to x. The Green's function G2 corresponds to the reflected arrival, which travels the greater distance r%, represented by the image-source to observer raypath. 16 represent the measured response of an impulse propagating as a spherical wave with constant velocity c and spherical amplitude wavefront divergence 1/r, traveling from source location x' to an observation point x in. afinitetime X = r/c, having been ignited at an initial time t'. The half-space Green's function can be interpreted in a similar fashion. G is composed of two such wholespace arrivals as described above. Thefirstarrival to be measured at x is due to G\, and is the direct arrival. This arrival has all of the characteristics mentioned above for an impulse propagating from the source location x' direcdy to the observation location x. The second arrival measured at x is due to the image source Green's function G2, and is the reflected arrival. This arrival accounts for the reflection of directly arriving incident energy from the impulse source location x', off of the surface ft, and back down toward the observation location x. Since the reflected travel path distance r 2 is longer than the direct travel path distance r l s the reflected event arrives at x later than the direct arrival by atimedifference of Ax = fo-rO/c). Naturally, the spherical divergence of the reflected arrival is larger than the direct arrival, and is expressed by the G 2 attenuation l/r2 versus 1/n of the larger amplitude direct arrival Gi. The boundary condition (2.9) defines ft as a free surface. The free surface reflection coefficient is (-1.), and is represented by the fact that the reflected arrival G2 is negative polarity compared with the direct arrival G\. Finally, for observation points very near the free surface ft, the direct and reflected arrivals are nearly of equal arrival times and mutually negative amplitudes. On the surfaceft,the two events cancel each other exactly, and no response is measured, as required by the boundary condition (2.9). A very small distance away from the surface ft, the arrivals imperfectly cancel, giving the appearance of a dipole, with the leading positive lobe attributable to the direct arrival. In the farfield situation, at observation points very far from both the source 17 and the free surface, the two arrivals again appear as a dipole. In the nearfield situation, at observation points close to the source but not the free surface, the reflected event has essentially vanished due to spherical divergence, and only the direct arrival has significant energy. This corresponds to the whole-space assumption, in which the boundaries of a region are so far away that they have essentially no effect on the distant interior wave solutions. 2.6 The Kirchhoff Integral Solution to the Scalar Wave Equation In this section we use the Green's function problem (2.8) and its solution (2.11) to derive the Kirchhoff integral solution to the scalar wave equation as posed in the primary problem (2.1). To begin, we form the difference of two product terms, U(x',t')L'G(x',t';x,t) - G(Z',t';x,t)L'U(x',t'), (2.14) where L! is the wave operator of the scalar wave equation (2.1) with respect to the primed impulse source coordinates (x', t'), L '= K ~ ^,t - (2-15) We then integrate (2.14) over the primed solution domain (x', t'), in both space and time: 18 1 = J JJJ \. &'> *') ' ^'' u t'' ' ')" L G( '' x G(x x ' ') rLr(x '' ] d x ' ( 1 1 6 ) Substituting for L'U and £'G from equations (2.1) and (2.8) respectively, relation (2.16) becomes 1 = J JJJ [- ^> ' ' W * ' " W " ') " ] 0 U{ rfx '' '' dt (- ) 2 1? which can be integrated out such that (2.16) reduces to I = -U(x,t), (2.18) which is the upgoing wavefield we seek to reconstruct. Hence, we can legitimately make the equation U(x, t) = - J JJJ^l f) £G(x', t'\ x, t) - G(x', t'\ % t) L'U(x', t') ] dx' dt'. (2.19) Next we apply Green's Second Identity, or as it is often called, Green's theorem (Morse and Feshbach, 1953, p.804), to the right hand side of (2.19) resulting in j JJJ^i u &> *') ' &> M .V *> *) ~ &> ?'> *> t) ' (*'> ?) ] *'> dt' = L G G L u d [ U(x', t') d >G(x', t'\ x, t) - G(x', t'; x, t) d , U(x', t') ] dx' dt', n n 19 (2.20) where d i is the derivative, with respect to the primed spatial source coordinates x', in n the direction of the unit vector h, which is normal to the bounding surface d'V. Using the boundary conditions (2.6) and (2.10), the contribution to the dV surface integrals (2.20) from the hemispherical surface of which x' e 9-(' vanishes identically. Hence, the only remaining contribution to (2.20) is from the planar surface such that, using relation (2.19), U(% t) = - f If [¥(*', t') oVG(x', t'; x, t) - G(x', t'\ % t) d ,¥(x', t')} dx'dt', (2.21) Jt' J JA. n where T(x', t') are the ZSR data, approximated by the NMO-corrected, CMP-stacked seismic reflection data, measured on the surface J3, as required by the original boundary condition (2.7): V(x',t') = tf(x',0|t' * 6 Equation (2.21) is equivalent to the starting point equation (2) of Schneider (1978). We can further reduce (2.21) by applying boundary condition (2.9), which states that the Green's function must vanish on the surface J?L U(x, t) = - f ff ¥(x', *') d G{x', t'; x, t) dx' dt' Jt' J JA n> (2.22) To proceed further, we need to evaluate the normal derivative of the Green's function, d >G(x', t'). This evaluation requires an explicit definition of the Green's function G(x', t') n in the new primed source coordinate system, which we shall now consider. 20 Recall the definition of G(x, t) from (2.11): Sjt-t-ri/c) G{%t\x',t') = 4nr { _ §(t-t'-r /c) 2 4nr 2 Since G depends on x and x' only in the radial distance terms r\, r , given by relations (2.12) 2 and (2.13), G is reciprocal with respect to the spatial coordinates: G(%t;x',t') = G(x',t;%t'), (2.23) However, G is not reciprocal in the temporal coordinates! Let us consider the impulse arrival at the observation coordinates (x, t) due to an impulsive source at located at (x', t'). The arrival at x occurs at atimet, which is equal to the time of source ignition t', plus thefinitetravel time from x' to x, x = rjc. We can describe this arrival as G(x, t = t' + rjc). Hence the delta function numerator of the Green's function must be specified as 5(i -1' - rjc), as is borne out in (2.11). Now let's perform the same physical analysis with respect to the primed source coordinates. The same arrival at x occurs now at a time t' in the new source coordinates, which is equal to the time of the arrival t in the old observational coordinates, minus the traveltime x = r\/c. In otherwords, at the moment of source ignition f, the arrival at the observation coordinate x will occur at a later time t in the future. We can describe this same arrival, but with respect to the source coordinates, as G(x', t' — t — rjc). Hence the delta function in the new source coordinate Green's function must be specified as b(t' — t + rjc). Notice that this change of coordinates is spatially reciprocal, but the temporal coordinates must obey the directional 21 "arrow of time." Given the discussion above, we can now define the coordinate-transformed Green's function G(x', t') as S( Of!, f.%t) = ''7 + r ' /C) - «?-* + '#). (2 . 4) 2 Now that we have defined the coordinate-transformed Green's function G(x', t'), we can return our attention to equation (2.22). If we perform the normal derivative of the Greens function, d iG(x', t'), using the definition (2.24), and evaluate the derivative on the surface n A, as indicated by (2.22), we find that d >G(x', Ok'e* = 23 ,G (x', OIl'eA n B 1 (2.25) where Gi(x', t') is the whole-space Green's function in the primed source coordinates, as given by the first term of (2.24): G,(xU';M) = 8 ( < '-/ + r i f c > . (2.26) Hence, equation (2.22) becomes U{%t) = -2 I ff^(x',t')d >Gi(x',t'-%t)dx'dt'. n (2.27) Since h is the unit, outward directed, normal to the surfaceft,and z' is positive downward 22 into the Earth, oV = -dz'. (2.28) Substitution of relations (2.25), (2.26), and (2.28) into equation (2.27) yields: U(x, t) = ^ J t ) JJWP, t')d , { z 5 ( f "^ + r / c ) } dx'dt', (2.29) where r is the radial distance r\ given by (2.12). This equation (2.29) corresponds to equation (4) of Schneider (1978). Equation (2.29) is the Kirchhoff integral solution for the reconstruction of the upgoing wavefield U{% t) which we seek. The wavefield reconstruction at each point in the subsurface domain is obtained by performing an integration over the boundary wavefield at the surface Si, (z=0), where, again, W(x',t') = l/(x',nii<E* is the recorded, NMO-corrected CMP-stacked surface seismic reflection data, as indicated by the boundary condition (2.7). Equation (2.29) could be used for migration if we were to apply a reflectivity imaging condition to the reconstructed wavefield. However, let us rearrange (2.29) into a more convenient form that avoids the explicit depth derivative and removes the time integration. Examining the depth derivativefirst,wefindthat 23 (2.30) d i = d d <r = -cos0(x')oV, z r z where the radial distance r is given by r\ of (2.12), and 0(x') is the angle between the ZSR raypath and the vertical, measured at the planar surface A Performing the radial derivative of the Green's function, as indicated by (2.29), we obtain b(t'-t + r/c)} bjt'-t + r/c) 2 + r h'jt'-t + rlc) rc ( } Substituting (2.30), (2.31) into (2.29), we obtain 8(t'-t + r/c) b'(t'-t + r/c) rc dx'dt'. (2.32) X After carefully transferring the delta-function derivative 8' to the observed wavefield F, and then performing the indicated t' integration, (2.32) reduces to U(x,t) = ± y^cos0(x') T(x', t - r/c) r 2 dtVjx'.t-r/c) rc dx'dy'. (2.33) Equation (2.33) corresponds to equation (5) of Schneider (1978), to which he subsequently applies the ZSR reflectivity imaging condition and obtains a 3-D Kirchhoff migration equation. 24 2.7 Defining the ZSR Reflectivity Imaging Condition Now that we have a method of reconstructing the upgoing wavefield U(x, t) given by the Kirchhoff integral solution (2.33), or (2.29), we need to apply a reflectivity imaging condition. A ZSR imaging condition, appropriate for the assumed CMP-stacked data we wish to migrate, is given by Claerbout (1971) and is interpreted by the exploding reflector model of Loewenthal et al.(1976). The ZSR reflectivity imaging condition states that the subsurface reflectivity distribution R(x) will be imaged by the reconstructed wavefield U(x, x), when the temporal argument x is evaluated for time x = 0, and at half the true medium velocity: R(x) = U(x,x) and c' = cll. (2.34) Let us briefly derive the ZSR imaging condition. We are considering the migration of ZSR reflection data, approximated by the CMP-stacked data. A ZSR reflection has the important property that the downgoing portion of the reflection path, is exactly coincident with the upgoing portion of the reflection path, as illustrated in Figure 2.3. Similarly, the downgoing traveltime from the source to the reflector is exactly equal to the upgoing traveltime from the reflector to the ZSR receiver location. Since the downgoing and upgoing raypaths and traveltimes are identical, we can replace the true experiment with a hypothetical exploding reflector experiment as follows. First, we 25 air Figure 2.3: The Raypath Geometry for ZSR Reflections. Due the the zero-offset nature of the source-receiver geometry, the upgoing and downgoing portions of reflected raypaths are coincident. Hence, the upgoing and downgoing traveltimes are equal. 26 place a source at the true reflection point, whose strength is equal to the local reflection coefficient. The hypothetical propagating medium is identical to the true medium, except that the velocity is halved. Then we ignite the source on the reflector at the recording time x = 0, which simulates an exploding reflector. The measured response at the surface receiver occurs at a time x = r/c' = 2r/c later, since the velocity has been halved. The geometry of the exploding reflector model is illustrated in Figure 2.4. The hypothetical exploding reflector experiment produces a reflected event measured at the ZSR surface location which has the exact same arrival time compared with the true downgoing/upgoing experiment. Furthermore, the raypath and exploding reflector source strength are the same as in the true experiment. These equalities suggest a way of imaging the subsurface reflectivity by using the exploding reflector model. To image the reflector in the exploding reflector model scenario, we would need to reconstruct the upgoing wavefield L7(x, T) from the hypothetical surface data, which we don't have, and image at time % = 0, just at the moment when the reflection sources were ignited. Since the surface data for both experiments are identical, we can reconstruct L7(x, T) by substituting the true ZSR upgoing/downgoing data *F for the unavailable hypothetical data, and then performing the wavefield reconstruction at half the true medium velocity. Finally, we can image the reflectivity from this reconstructed result at the exploding reflector time x = 0. Hence, we have derived the ZSR reflectivity imaging condition stated by (2.34). 2.8 The Conventional Kirchhoff Post-Stack Migration Equation Now we can apply the ZSR reflectivity imaging condition (2.34) to the Kirchhoff wavefield reconstruction integrals given by (2.33), which results in the full 3-D constant-velocity 27 air Figure 2.4: The Exploding Reflector Model. A hypothetical model in which point-sources are positioned along subsurface reflecting boundaries. The point-source amplitudes are given by the local reflectivity values, and the sources are ignited at time x = 0. The velocity c' of the model is half the true medium velocity. Under these conditions, the hypothetical trace data resulting from the exploding reflector model are equivalent to the true ZSR experimental data. 28 Kirchhoff post-stack migration equation: cos 0(x') x', t = 2r/c) 2df¥(x', t = 2r/c) dx'dy'. (2.35) The two terms in solution (2.35) can be recognized as the nearfield (first) and farfield (second) terms of the scalar wave potential (Berkhout, 1984). The nearfield term dominates in the region where the radial distance is much less than a few wavelengths of seismic energy, and the farfield term dominates the region in which the radial distance is much greater than a few wavelengths. In the typical seismic reflection acquisition case, the farfield term is largely predominant. In that case, we obtain a farfield approximation to the 3-D constant-velocity Kirchhoff post-stack migration solution: (2.36) Equation (2.36) is the common Kirchhoff post-stack migration equation. Note that, strictly speaking, (2.36) requires seismic data recording coverage over the entire planar surface A. In practice, the ^'-integration is ignored reducing (2.36) to a 2-D migration equation; (2.37) only the ^'-integration is computed since typical seismic profiles are recorded along a line, rather than a plane, of receivers. In the prestack Kirchhoff-WKBJ migration theory that follows, the ^'-integration will be analytically performed, resulting in a more accurate 2.5-D migration equation. 29 2.9 Interpreting the Kirchhoff Migration Equation The farfield Kirchhoff migration equation (2.37) is now in a form that can provide some physical intuition into the migration process. The first remark of note is that the CMPstacked seismic reflection trace data ^(x, t) are summed spatially over thetimetrajectories t = 2r/c. In the 2-D (x, t) data plane thesetimetrajectories sweep out hyperbolae. In the 3-D (x, y, t) data volume they sweep out hyperboloid surfaces. The hyperbolae correspond to the potential diffraction-hyperbolae that would arise due to any point-diffractors, or pointscatterers in the subsurface medium. Equation (2.37) suggests that we sum along hyperbolae and place the resulting energy of each hyperbola at its corresponding diffraction-apex location in the reflectivity image space R(x). Successive integrations over many traces and many hyperbolae will tend to reinforce the reflectivity image, by constructive wave interference, in regions that truly contain pointdiffractors, and wash out the reflectivity image, by destructive wave interference, at locations where point-diffractors do not actually exist. Recalling the exploding reflector model of Figure 2.4, we see that an image of the subsurface point-diffracting sources which generate the upgoing wavefield U{x, t) will indeed image subsurface reflecting interfaces. This intuitive approach also suggests why diffracted as well as reflected energy tends to be correcdy imaged by the migration process. Prior to the advent of wave-equation migration, this intuitive approach of summing along diffraction-hyperbolae was used to "migrate" CMP-stacked data under the name of diffraction-stack migration. The diffraction-stack process is based on the concept that, given an appropriate velocity structure, the measured surface wavefield can be mapped back into the depth domain and re-focussed to image the diffracting sources of the reflected wavefield. 30 The main difference between a diffraction-stack operator and the Kirchhoff migration equation (2.37), is that the latter attempts to respect the intrinsic physics of wave propagation which is inherent in the seismic reflection data. Most importantly, (2.37) attempts to migrate the data in a manner which satisfies the amplitude and phase requirements placed upon the operator by the scalar wave equation. The 1/r amplitude factor in (2.37) represents the geometric divergence of a spherical wavefront in a constant velocity medium. The cos 9(x') amplitude term is often called the obliquity factor, and represents the incident-angle dependency of wave energy arriving obliquely at the free surface A and being subsequently projected onto vertical component geophones (in the case of land data), or pressure-sensitive hydrophones which measure the normal derivative of ambient pressure at the recording surface. The factor of (-2), which arises due to the evaluation of the normal derivative of the Greens function, is interpreted as the acoustic free-surface conversion coefficient. Finally, the time derivative applied to the surface seismic data *F corresponds to the well-known 90° phase rotation of the source waveform caused by its propagation from the source region (nearfield) to a distant (farfield) observation point. In particular, note that the full solution (2.35) allows for variable phase along the Kirchhoff wave-equation migration operator's hyperbolic trajectory, ranging from 0° (nearfield) to 90° (farfield) depending on the radial distance and frequency content of the surface observation. As afinalimplementational remark, the 3-D Kirchhoff post-stack migration integrals, (2.35) and (2.36), and the 2-D integral (2.37), can be discretized into summations, to accommodate the discretely sampled seismic traces that exist at discrete locations on the recording surface. The resulting summations do not require regular trace spacing intervals, and are thus wellsuited to the often,irregular recording geometries present in real seismic data surveys. The Kirchhoff migration solutions can also be generalized for variable velocity models c(x) and 31 prestack geometries, as will be developed next in the main body of this thesis. 32 Chapter 3 Derivation of a True-Amplitude Kirchhoff-WKBJ Prestack Depth Migration Equation 3.1 Introduction This chapter forms the core theory component of my thesis research. The central theme is the derivation of a true-amplitude reflectivity migration solution. The main motivation for true-amplitude migration is the requirement for accurate amplitude recovery in the reflectivity imaging step of a more complete formalism for the generalized inversion of seismic reflection data, which is my larger ongoing research interest. Of course, true-amplitude migration could also be valuable for the interpretational analysis of seismic reflection data, since velocity structure and geophysically interesting reflection-amplitude anomalies may be more accurately estimated and hence more readily subject to "direct" inversions. It is my contention that, in general, true-amplitude reflectivity is not recoverable from poststack processed data. Valuable offset and raypath information, destroyed in the stacking process, is vital to the reconstruction of reflectivity amplitudes, and therefore suggests the requirement for a prestack methodology. Furthermore, true-amplitude reflectivity is dependent upon the incorporation of realistic amplitude-attenuation recovery in the Kirchhoff wavefield reconstruction step. This is achieved by extending the Kirchhoff theory with the use 33 of WKBJ ray-theoretic Green's functions in place of the conventional simplistic constantvelocity Green's functions. The WKBJ Green's functions are extremelyflexiblein the presence of 2-D migration velocity model heterogeneity, and can be numerically evaluated for rxaveltime and amplitude response by raytracing orfinitedifference techniques. Finally, the requirement for migrated depth models of reflectivity arises from physical considerations. Time migration images are non-physical from an inverse theory perspective, from which migration can be viewed as an approximate inverse mapping resulting in the construction of a reflectivity model, as a function of spatial coordinates, which satisfies the observed seismic reflection data. In summary, the above considerations regarding the recovery of true-amplitude reflectivity from seismic reflection data lead to the requirement for a prestack, amplitude-preserving depth migration theory. In this chapter, such a theory is invoked by developing a hybrid KirchhoffWKBJ solution to the scalar wave equation and applying physically-defined reflectivity image construction and reflectivity amplitude estimation principles. 3.2 A True-Amplitude Migration Philosophy The derivation of a true-amplitude prestack depth migration equation is based upon the following conventional migration philosophy: Migration = wavefield reconstruction + reflectivity imaging. Virtually all migration methods can be viewed as different approaches to the same two-part problem. Thefirstpart of the migration problem is the reconstruction of the downgoing 34 (source) and upgoing (reflected) wavefields throughout the subsurface given the recorded surface data as a boundary condition to the wave equation. Conventionally, the second part of the migration problem requires the definition and application of a reflectivity imaging condition, in order to effectively image the subsurface reflectivity distribution R(Xo). Since we are interested in true-amplitude reflectivity model construction rather than imaging, we need to physically define and mathematically apply a reflectivity inversion condition as a function of the reconstructed wavefields, to both image and amplitude-estimate the reflectivity distribution R(x ). 0 First, we will examine the true-amplitude recovery of reflectivity fromfirstprinciples. This examination emphasizes migration as an inversion of the observed seismic data in the construction of a model describing the subsurface spatial distribution of the physical quantity: acoustic reflectivity. Migration posed as an inverse problem represents a strong philosophical departure from the conventional, and somewhat less demanding perspective of migration as an image enhancement process. However, it is precisely the inversion perspective which establishes a physical migration philosophy capable of guiding the mathematical derivation to its goal of true-amplitude reflectivity model construction. Next, wavefield reconstruction solutions for the upgoing and downgoing wavefields are obtained by extending the conventional Kirchhoff integral theory with the use of WKBJ raytheoretic Green's functions. Based on the section regardingfirstprinciples of migration from an inversion perspective, a true-amplitude reflectivity imaging condition is physically defined and mathematically applied as a function of the reconstructed wavefields. The result is the desired true-amplitude Kirchhoff-WKBJ prestack depth migration equation. 35 3.3 First Principles Suppose a given seismic experiment produces a shot gather of data, ^(xr,fr,xs), where: *F is the observed pressure wavefield, xs = (x ,y >Zs) is the (fixed) shot coordinate, and s s xr = (x , y , z ) is the (variable) receiver coordinate. r r r For simplicity, we assume that the shot and receivers are coplanar, although this is not required by the method. The planar data surface may be the actual physical recording surface, or as is often the case, may represent some convenient imaginary datum plane. The idealized shot gather geometry is illustrated in Figure 3.1. From the shot to the subsurface the source energy propagates as downgoing waves, D(\o,fr,xs), where Xo = (x , yo, z ) is the subsurface model coordinate. 0 0 Upon reflection, a portion of the downgoing energy is converted into upgoing waves, denoted by U(\ ,fr,xs). When the upgoing waves impinge upon the recording surface, we record the 0 resulting observed seismic data: V(x ,fr,Xs) = L/(Xo,fr,%)\^ . r The upgoing wavefield U, and hence r (3.1) contains energy from primary and multiple re- flections, but is also contaminated by upgoing diffractions, mode conversions and refracted 36 source receivers Figure 3.1: An Idealized Shot Gather Geometry. The flag symbol locates the shot source position xs, and the inverted triangle symbols locate the receiver positions xr. 37 energy. We will assume that the upgoing wavefield U and the observed data *F can be processed to ideally represent energy due to primary reflections alone. The objective of conventional migration is to obtain an image of the subsurface earth reflectivity distribution i?(x ), which upon secondary excitation by the incident source field 0 D(x , t; x ), gives rise to the upgoing reflection wavefield L7(x,fr,x ), subsequently mea0 s 0 s sured at the recording surface as the observed seismic reflection data ^(x,.,fr,x ). s Unconventionally, the objective of this thesis is to obtain both a reflectivity image compo- nent and a true-amplitude reflectivity-estimate component of the true subsurface reflectivit distribution. In effect, we seek an inversion of the surface seismic reflection data in order to construct a physically reliable model of the true subsurface reflectivity. Let us examine the problem fromfirstprinciples. If at any reflecting point in the subsurface there exists a downgoing wave incident upon the reflector, then at an incremental time later, a reflected upgoing wave will be emergent as shown by the simplified ray diagram in Figure 3.2. If the propagation distances involved from the point of incidence to the point of emergence are very small, then the amplitude of the upgoing wave is approximately equal to the amplitude of the downgoing wave, scaled by the local angle-dependent reflectivity: U = R(x ,Q)D. 0 (3.2) This assumes that the geometric divergence, transmission loss, and phase rotation along the incremental raypaths are negligible. For the moment let us omit the explicit angle dependence of the reflectivity and assume its implicit notational dependence: 38 air Figure 3.2: An Idealized Reflection Raypath Geometry. An ideal primary reflection occurs when a downgoing wave D is incident upon a subsurface reflector with local reflectivity magnitude R. A reflected upgoing wave U is generated which, in the local vicinity of the reflection point, can be represented by the downgoing wave scaled by the reflectivity magnitude: U = RD. 39 #(x ,e)->E(xo). 0 In the spectral domain, (3.2) implies that (3.3) R = U/D where U is the temporal Fourier transform of U. Usually seismic reflection data are observed over a range of frequencies, which would generalize (3.3) to a reflectivity function dependent upon position, angle of incidence and frequency: R(x , 0, w). In practice, we will assume 0 that R is independent of frequency but representative of the dominant-frequency behaviour of the bandlimited seismic reflection data. We desire to transform (3.3) back into the time domain. Multiplying numerator and denominator by the complex conjugate of the downgoing wavefield, D*, results in „ UD R = DD* (3.4) v ' We recognize the numerator of (3.4) as the spectral equivalent of a time-domain crosscorrelation. Since we do not desire R to be a function of crosscorrelation lag, we take the zero-lag contribution of the reflectivity. Hence, the numerator of (3.4) transforms to UD* -> U*D 40 (3.5) and 7(x ;x ) 0 s = U*D\p=o (3.6) where the * symbol denotes correlation and the |^=o indicates the zero-lag evaluation. The term I represents the image component of the recovered reflectivity distribution, since the phase of (3.4) is contained in the upgoing-downgoing crosscorrelation. The denominator of (3.4) is shown below to control the relative amplitude within the reflectivity image. The combination of reflectivity structural imaging, and reflectivity amplitude estimation within the structural image, constitutes an inversion for a physical model estimate of the true reflectivity distribution itself. The denominator of (3.4) is the spectral equivalent of a time-domain autocorrelation. The zero-lag contribution is effectively an average measure of the local energy of the downgoing 2 wavefield: | D | . Analogously, this average energy may be thought of as representative of a volume integration of the downgoing energy over a tiny sphere centered about x„. In practice, we calculate this value using dynamic raytracing amplitudes, as discussed later. The denominator of (3.4) transforms to DD* ->D*D (3.7) D*D (3.8) and \D\ = 2 Finally, in the space-time domain (3.4) becomes / ( x ; Eg) 0 R = 41 (3.9) We shall later call (3.9) the dynamic imaging condition, in the sense that it is a true-amplitude equivalent of a conventional migration reflectivity imaging condition. However, we may now remark that (3.9) should be classified as a reflectivity inversion condition, since it constructs both a reflectivity image and a model estimate of reflectivity amplitude. The form of condition (3.9) implies that we reconstruct the upgoingfield,U(x , t; xs), and the downgoingfield, 0 D(x , t; xs), perform their zero-lag crosscorrelation, normalize by the energy of the downgoing 0 field, and hence obtain the reflectivity model estimate R(x ). Note that the implementation of 0 (3.3), (3.4) or (3.9) may require stabilization in regions of weak source illumination energy, as discussed in later sections of this thesis. 3.4 The Upgoing Wavefield Solution Let us examine how to reconstruct the upgoing wavefield L7(x0, t; xs) needed to solve for the reflectivity model estimate R(x ) in Equation (3.9). 0 Claerbout (1971) uses the upgoing wavefield recorded at the surface as a boundary condition to the scalar wave equation. This assumes an acoustic wavefield as opposed to the true elastic nature of the data. Usingfinitedifference operators to approximate the wave equation, Claerbout reverse-time extrapolates the boundary wavefield in the space-time domain. This approach forms the basis offinite-differencewave-equation migration. A major advantage of this approach is that an arbitrary 2-D migration velocity model can be readily specified on thefinite-differencegrid. Another advantage is that full wavefield effects such as postcritical reflections, critical refractions, and multiples may be directly incorporated in the migration without pre-processing the recorded data to remove these effects. Furthermore, thefinite-differencealgorithm can be extended to the elastic wave equation to produce an 42 elastic wave-equation migration of both compressional and shear wave energy. Another common approach to reconstructing the upgoing wavefield is to use Fourier Transform methods. This alternative is given by Stolt (1978) for a constant-velocity migration model, and by Gazdag (1978) for vertical velocity variation. The boundary data are Fouriertransformed into the spectral domain whereupon a phase-shift downward-continuation operator is applied to reconstruct the entire wavefield. The greatest advantage of this spectral method is its rapid computational speed due to the use of the Fast Fourier Transform (FFT), and the fact that wave extrapolation reduces to simple multiplication in the Fourier domain. Alternatively, we have used a Kirchhoff method to reconstruct the upgoing wavefield by applying the Kirchhoff boundary integral solution to the homogeneous scalar wave equation (e.g., Berkhout, 1986). The appropriate wave equation governing scalar wave propagation in a linear, acoustic, homogeneous and isotropic medium, in the absence of body forces, is the homogeneous Scalar Wave Equation: 2 V C/(x0, t; Xs) - \dt U(Xo, t c l t; xs) = 0, (3.10) subject to boundary and initial conditions. The parameter c represents the phase velocity of waves propagating in the medium. Strictly speaking, c is constant in the derivation of (3.10), since the acoustic parameters of the medium that define c are constant. In practice, we allow these parameters to be spatially variable in a "smooth, slow-varying" manner with respect to the dominant seismic wavelengths to be considered. Allowing a heterogeneous acoustic medium effectively introduces extra terms to (3.10) which depend on the gradient of the acoustic material properties defining c. If the acoustic parameterfield,and hence c, 43 is spatially smooth and slow-varying as assumed, then the gradient terms become negligible and (3.10) remains the appropriate equation of motion, to zeroth order. A classical solution to (3.10) in a halfspace is the Kirchhoff Integral Solution: U(x0 , t; xs) = -2 / / / ¥ ( x r , t'; xs) d G (xr, t'\ x0, t)dxr dyr dt', n r derived previously as Equation (2.27), where d is the derivative normal to the recording n surface ft, and G\ is an appropriate wholespace Green's function as yet to be specified. Note that (3.11) assumes that the complete halfspace Green's function G^, of which the wholespace Green's function Gi is one component, vanishes on the observation plane ft. Equation (3.11) is identical to Equation (5) of Carter and Frazer (1984), except that their equation is Fourier-transformed with respect to time. Equation (3.11) is also identical to Schneider's (1978) Equation (4) in which he uses the constant-velocity reciprocal halfspace Green's function: G h ^ M = G 1 + G 2 = ^ ' - 4jcri ^ . * - « ' - ' * > 47tr2 the radial source and image distances are given by 2 + r\ = {Xr-Xof + iyr -y0 ) (zr -z0 f r\ = (xr -x0 f + (yr -y0 f + (zr +z0 f. and ; ( 3 . 1 2 , (3.1 If the velocity is a function of x and z , i.e. c(x , z ), then we can look for a WKBJ Green's 0 a 0 Q function solution by assuming a frequency domain representation of the form: ±l M; G(x , w'; x ) = A(x ; x )S(M;')e ' ' r 0 r 0 t(Xr;Xo) , (3.13) where the G\ subscript has been omitted for subsequent notational simplicity (i.e., G\ —> G). Equation (3.13) describes a wholespace WKBJ Green's function which is suitable for heterogeneous media, in contrast with the piece-wise constant velocity Green's function Gi of (3.12). We will see later that the WKBJ solution corresponds to the geometrical optics, or zeroth order asymptotic ray theory, approximation of wave propagation. The WKBJ approach to the seismic reflection problem is discussed by Clayton and Stolt (1981) in the context of inversion and by Castle (1982) in the context of migration. The WKBJ Green's function of Equation (3.13) has an amplitude term A(xr;Xo) which represents the amplitude attenuation (largely due to geometric divergence) from the subsurface location XQ to any surface receiver location x . The phase term t(x ; x ) represents the traveltime from a given r r 0 subsurface location to the receiver locations . The function S(w') represents the temporal Fourier transform of atime-domainsource function S(t'). The choice of sign in the complex exponential depends upon the time-direction of wave propagation. Since we are considering energy propagating from a particular subsurface location x to all possible surface receiver locations x , the choice of negative phase is correct for 0 r the desired upgoing waves. This is consistent for a depth axis defined as positive downward and a negative complex exponential as the kernel of the forward Fourier Transform. In thetimedomain, (3.13) inverse-transforms to 45 G(x , t'\ x ) = A(x ; x0 )S(t' - x(xr ; x )). r 0 r 0 (3.14) We can time-shift S(t') by a value t corresponding to the "ignition-time" of the Green's function source impulse. This generalizes (3.14) to G(xr, t'\ x0, t) = A(xr; x0 )S(t' - T(xr; x0) -1). (3.15) Using the form of the Green's function (3.15), the normal derivative in (3.11) can be reexpressed as d G = VGh n (3.16) where V is the gradient operator with respect to the subsurface coordinates x0, and h is the unit vector normal to the recording surface. Expanding the spatial gradient of the Green's function, we obtain VG = SVA+AS VT, T (3.17) where Sz is the derivative of <S with respect to x. Since the Green's function is a 3-D pointsource impulse response, its amplitude A decays approximately as 1/r, and hence the gradient 2 of A decays approximately as 1/r. Because of this the second term in (3.17) is dominant for large r. We neglect thefirstterm in (3.17) and anticipate a farfield approximation to the full solution: 46 VG~AS Vx (3.18) VGn~AS Vx-n. (3.19) x and T From (3.13) we see that x(xr;x0) is the phase of the WKBJ Green's function and as such can be shown to satisfy the eikonal equation from asymptotic ray theory (ART) considerations (Cerveny et al., 1977): ^ = s k y ( 3 - 2 0 ) Also, since the loci of constant x(xr; x0) map out a wavefront surface, Vx is the ray normal to that wavefront, and hence V x .=COS0(Xo) c(x0) where 0(xo) is the angle between the ray propagation direction incident at a subsurface location x0 and the local normal to the ray-emergent surface location xr. Note that from (3.14) we can re-express Sx as Sz = St > (3.22) where Sp is the derivative of S with respect to t'. Combining Equations (3.13) through (3.22) into (3.11) results in 47 U(x ,t;%) = 2 III 0 Jt'JjA ^^ACxrjXoJSr^'-x^XoJ-O^^Xs)^^^' c (3.23) ( r) x where 6(x ) and c(x ) have been evaluated at the surface location x as required by the surface 0 0 r integral. We assume that the form of the source time function S(t') can be represented by a Dirac delta function: 8(f')- Then we use the property of the delta function: J8 (x) q>(*)dx = (-1)" J 8(a) <? \x)dx, (n) (3.24) {n to transfer the time derivative from &(t') to W. Also, noting from (3.15) that dfS = -d S (3.25) t and performing the indicated time integration in (3.23) results in r r c o s 0( ) x t/(x , t; x ) = 2 / / — - i - I i A ( x ; x„) 3^(x , t + x(x ; x ); x ) dx dyr . 0 s JJj? c(x ) r r r G s r (3.26) r Equation (3.26) is in the correct form but requires seismic data recorded over the entire (%r, yr) plane rather than the more commonly available in-line profile for a given shot. To overcome this difficulty, we assume that the subsurface geology is constant in the y-direction, and integrate the ^-dependence out of (3.26), e.g., Bleistein (1986). The ^-integration can also be acheived by consideration of its physical interpretation. The integration essentially converts the 3-D point-scatterer response of a reflecting segment at x 48 0 to that of a 2-D line-scatterer response. This results in the amplitude A(xr; %) diverging as l/y/r instead of 1/r as the reflected energy propagates from the reflection location XQ to the receiver location xr. The nearfied-to-farfield spatial phase shift of such a line-source is 45° rather than the 90° shift due to a point-source. This phase-shift property is controlled in the space-time domain by the point-source spatial derivative dtlc integrating to the line-source l2 spatial derivative d] /y/c. The term d] is known as the "half-time-derivative operator" and 12 is defined as the inverse Fourier transform of the spectral operator \/iw. The net result of the ^-integration is the 2.5-D solution for the upgoing wavefield: L7(x0, t; xs) = 2 J cos 6(x)y 1 / 2 a, ¥(xr, t + x(x ; x ); x ) dx r r 0 s r (3.27) in analogy with Equations (13) and (14) of Carter and Frazer (1984). 3.5 The Downgoing Wavefield Solution Now we calculate D(x t;x ) as required by the reflectivity inversion (condition 3.9). If we 0) s assume that the downgoing wavefield is composed of the direct arrival only from the source location and does not contain scattered or multiply reflected energy, then the solution for D simply takes the form of the Green's function response to the medium. Assuming the form of a WKBJ Green's function again: D(x0, w; xs) = A(x0; x )S(w)e- ^ ™ iw s 49 0 (3.28) Notice that the phase shift is required to be negative for downgoing, outward-propagating waves, as opposed to positive phase shifts for reverse-time depropagating waves. In the time domain, (3.28) is D(x , t; Sg) = A(x0; x )S(t - t(x0; xs)). 0 (3.29) s Thus, D(x , t;x ) is the original source wavelet S(t) time-shifted by the traveltime 0 s T(X ;X ) 0 S from xs to x0 and amplitude-attenuated by the 3-D point-source divergence factor A(x0; xs). 3.6 Applying the Reflectivity Imaging-Inversion Condition We have now reconstructed the upgoing wavefield everywhere in the subsurface medium assuming 2-D geology and a 3-D point-source illuminationfield,hence we obtain a 2.5-D migration solution. Now we use the reflectivity imaging condition (3.6), and crosscorrelate the upgoing and downgoingfieldsretaining only the zero-lag reflectivity image contribution Z(x0;xs), (3.30) Substituting U and D from equations (3.27) and (3.29), (3.30) becomes f 7(x0; xs) = / A(x0; xs)S(* - T ) • S 2 cos Qr A(xr;x0) Bl^ixr, c(xr) 50 t + T ; Sg) dx • dt. (3.31) r r Assuming again that the source function S(t) takes the form of a delta function, we perform thetimeintegration in (3.31) resulting in 1/2x 7(Xo;xs) = 2AS f cos9 r J — at J y cy P(xr, t;xs)|^s+Trdx . (3.32) r To complete the reflectivity inversion condition (3.9), we must normalize (3.32) by the energy of the downgoing wavefield which we associate with the square of the Green's function amplitude term in (3.29): 2 2 | D(x ; Xs) | = A (x0; xs) = A 2 (3.33) 0 Substituting (3.32) and (3.33) into (3.9) we obtain the migrated reflectivity estimate for a fixed shot gather: Rbo,Zo\V) = .n?''*! ] ,2 = / ^ 1 I •L'l.Xo) s ) I x 3.7 J y/Cr ^ ( * r , t;%)\ ^ dx . t r r (3.34) As The True-Amplitude Kirchhoff Migration Solution Equation (3.34) gives a true absolute-amplitude estimate of the subsurface reflectivity for a single shot gather at xs. In practice, the absolute source amplitude A is never known, which s reduces (3.34) to a measure of true relative-amplitude reflectivity. The reflectivity estimate is an average with respect to the reflection-incidence angle 9, since, within a shot gather, (3.34) averages the true angle-dependent nature of the reflectivity over many shot-receiver 51 offsets and hence over many reflection raypath angles. Equation (3.34) provides a robust reflectivity estimate in subsurface regions of high signalto-noise ratio. In areas of weak source illumination, little or no reflectivity information is recoverable. This suggests performing shot-gather migrations at different xs locations, and compositing the resulting reflectivity estimates. Mathematically, we achieve this composite by integrating (3.34) over a finite range of recorded shot locations to obtain an estimate of the subsurface reflectivity independent of shot location: R(x z )= m 0 ff d]*V&r,t Z )\t dx dx . t a Wr r a (3.35) Equation (3.35) provides a true-amplitude reflectivity estimate, if subsequently normalized by the shot line-source length and a true relative-amplitude reflectivity estimate otherwise. The composite full-fold multioffset migrated reflectivity estimate of (3.35) is effectively averaged over all individual shotdependent ray-incident angles and shot-recording-geometry subsurface-illumination patterns. It is important to realize that the total seismic experiment may not sufficiently illuminate certain portions of the subsurface resulting in the loss of reflectivity resolution and accuracy in thefinalmigrated depth section. Hence, the illumination dependence and the specular dependence affect one's ability to recover true-amplitude reflectivity. Equation (3.35) gives the solution for the reflectivity in terms of a Kirchhoff prestack depth 52 migration which incorporates both imaging traveltimes and WKBJ amplitudes. The traveltimes and amplitudes can be computed by a raytracing algorithm given an arbitrary migration velocity model. Specifically, then, for the implementation of the migration Equation (3.35), we repeat the notational definitions: xr is the receiver coordinate, xs is the shot coordinate, Xo is the subsurface depth-point coordinate, R(x , z ) is the migrated depth estimate of subsurface reflectivity, 0 0 d] is the half-time-derivative operator, 12 ^P(xr, t, %)\t=% +i s r are the recorded shot-gather traces evaluated at the total shot- subsurface-receiver traveltime, xs = T(X , X s ) is the raytraced traveltime from xs to x0, 0 %r - i(xr, x0) is the raytraced traveltime from x0 to xr, A — A(x0, xs) is the 3-D amphtude attenuation from xs to x0, s A = A(xr, x0) is the 3-D amplitude attenuation from x0 to xr, r c = c(xr) is the velocity evaluated at the receiver location xr, and r 0 = 0(xr) is the ray incident angle at xr with respect to the surface normal for r a ray traveling from x0 to xr. Finally, an implementational remark concerning the Kirchhoff-WKBJ prestack depth migration theory just derived is in order. No special pre-migration data processing is required other than standard trace-editing, bandpassfilteringand the appropriate muting or attenuation of post-critical reflections, critical refractions, converted waves, surface waves and multiples: all of which are not included in the fundamental mathematical assumptions of the method. Any specific sub-portion of the reflectivity model may be estimated as specified by the user. 53 This contrasts with other migration methods which may intrinsically require that the full subsurface model be wave-extrapolated and imaged at every step. Lastly, Kirchhoff migration is amenable to non-uniform trace spacing and irregular trace gathering. This is advantageous for the common situation in which adverse recording geometries are necessary to acquire seismic data, and for the avoidance of large, costly computer sorting of trace data in the processing flow. 54 Chapter 4 The Migration Reflectivity Imaging Conditions 4.1 Introduction The previous chapter developed a detailed derivation of a true-amplitude prestack depth migration equation using Kirchhoff-WKBJ theory. The solution fundamentally involves a reconstruction of the necessary upgoing and downgoing seismic wavefields, followed by the application of a physically-defined reflectivity imaging condition. Since we emphasized a migration inversion philosophy, the true-amplitude reflectivity imaging condition may be more accurately termed a reflectivity inversion condition. In this chapter, we examine several "imaging" conditions, in order to compare the true-amplitude results of this thesis with other current non-true-amplitude migration imaging principles. To compare various migration imaging conditions with the new reflectivity inversion condition, it isfirstnecessary to cast all of the relevant migration formalisms into a common mathematical framework. Once this is accomplished, the reflectivity imaging and inversion conditions can be subjected to comparative mathematical, physical, and quantitative analysis on a common footing. For this reason, a new generalized Kirchhoff prestack depth migration equation is derived containing a weighting kernel function W, which governs the reflectivity imaging-inversion condition specification. 55 The generalized Kirchhoff migration equation is subsequently used to consider five distinct imaging conditions: dynamic, excitation-time, crosscorrelation, and the new geometric and kinematic conditions. A special emphasis is placed on the physical interpretation of the reflectivity imaging-inversion conditions, in terms of their potential for image resolution and amplitude accuracy. 4.2 The Generalized Kirchhoff Migration Solution Following the reconstruction of both the upgoing and downgoing wavefields U and D, the migration philosophy can be completed by the physical definition and mathematical application of a reflectivity imaging condition. By applying each imaging condition as a function of the reconstructed wavefields, and assuming a 2.5-D farfield geometry as before, the following new generalized Kirchhoff prestack depth migration equation is obtained: R{x ,z )= 0 0 I f ^ Cr^ ^ y ^A s J JA V dl V(x ,t,x )\ ^ dx dx . a r s t r r s (4.1) Notice that the generalized solution (4.1) is equivalent to the true-amplitude reflectivity solution (3.35), except that we have now introduced a weighting function, *W, into the integration kernel. The exact form of W is obtained by performing the mathematical application of a given imaging condition to the upgoing and downgoing wavefield solutions (3.27) and (3.29), as was done in the previous section. In general, the weight TV is a function of both the upgoing and downgoing wavefields, and its specification governs the selection of the particular reflectivity imaging condition to be implemented as described below. 56 4.3 The Dynamic Imaging Condition • R (x ) = | [U(x , t) *Z)(xo, t)] I [D(x0, t) * D(x , t)] |, . d • Wd 0 0 0 (=0 = 1. Assuming that the upgoing reflected wavefield is approximated by the downgoing wavefield convolved with the reflectivity distribution, then a rearrangement of terms defines the dynamic imaging condition: R = U/D, (Claerbout, 1971). This imaging condition corresponds to d the physical notion that local reflectivity is estimated by dividing (deconvolving) the local upgoing wavefield by the local downgoing wavefield in the spectral (temporal) domain. The dynamic imaging condition has the correct formulation to recover true-amplitude reflectivity, and hence qualifies as a relectivity inversion condition, but requires stabilization in subsurface regions where the energy of the downgoing wavefield D is small. The dynamic migration can be implemented by setting = 1 in (4.1). The dynamic imaging condition formed the basis for the derivation of (3.35), which was derived in order to recover true-amplitude reflectivity. Accordingly, (4.1) reduces to (3.35) for the special case that Ufa = 1. Equations (3.35) and (4.1) appear to be new forms of the Kirchhoff migration equation and as yet unpublished in the geophysical literature. Keho and Beydoun (1988, Equations (5)-(7)) obtain a frequency-domain result somewhat similar to (3.35), however, their equation does not correctly compensate for the 2.5-D effect of outof-plane reflector excitation, and hence they overestimate the amplitude divergence along the reflected portion of a raypath, and incur a 45° phase error. Theoretically, we expect the dynamic imaging condition to provide the best recovery of true relative-amplitude normal-incidence reflectivity in a depth-migrated section. The dynamic 57 imaging method is derived from strong physical considerations which emphasize an inversion for the absolute-amplitude reflectivity distribution values rather than a simple image of reflector locations. In practice, the dynamic imaging condition may produce unstable reflectivity estimates in regions of weak source-illumination energy, i.e., where A —> 0. In the s frequency domain this problem arises when the spectrum of the downgoing field contains (near) zeroes. The instability can be overcome by adding a "stabilization parameter" to the source amplitude term such that: 2 A s -> A s + e . 2 The parameter e is generally taken to be constant with its magnitude being a small fraction of the average value of the source amplitude A . For an optimal stabilization, the parameter s should be set to the estimated background-noise amplitude-level, which in general is a function of x„. Obviously, the suggested amplitude decomposition of the data into signal and noise components may not be trivial. 4.4 The Crosscorrelation Imaging Condition • .Rc(xo) = U(x ,t)-kD(xo,t)\t>^o. 0 • <W = A . 2 C Claerbout (1971) overcomes the dynamic imaging condition instability resulting from the spectral division U/D, by defining a crosscorrelation imaging condition: R c = U*D, evaluated at zero-lag. Crosscorrelation imaging is equivalent to dynamic imaging except that 58 the required normalization by the energy of the downgoing field is subsequently ignored. Physically, the crosscorrelation method retains the correct dynamic imaging characteristics, since its imaging phases (times) are identical to those of the dynamic image. However, the accuracy of the amplitude information within the image is sacrificed in return for the enhanced image stability. The crosscorrelation imaging condition migration can be implemented by setting the weighting factor equal to the square of the source amplitude term: W = A , in 2 c S (4.1). Claerbout usesfinite-differenceoperators to downward-continue the upgoing wavefield and forward-propagate a source wavefield. In this way neither the source nor receiver amplitudes need be calculated separately since they are intrinsically computed as part of the wavefield extrapolations. However, to implement a Kirchhoff crosscorrelation migration with Equation (4.1), the A amplitudes must be determined by external means such as raytracing. r 4.5 The Excitation-Time Imaging Condition • -Re(Xo) = L7(X0,£)|*=I S • W = A. e 8 If the downgoing wavefield D is assumed to be everywhere a unit amplitude delta function, then the crosscorrelation imaging condition reduces to the excitation-time imaging condition: R e = L7(xo, t — T s ). Physically, this is equivalent to the assumption that the reflectivity is equal to the downward-continued upgoing wavefield evaluated at the source-to-subsurface ray traveltime x(x ; x ). This assumption would be valid for a delta-function source of unit 0 s amplitude which does not attenuate nor disperse along the raypath from the source location to the reflection point. 59 This imaging principle is used in a reverse-time prestack depth-migration finite-difference algorithm (Chang and McMechan, 1986), and in a prestack Kirchhoff migration method (Etgen, 1987). Most post-stack migration methods based on the "exploding reflector" model implicitly obey the excitation-time imaging condition since ZSR reflection is treated as oneway upgoing propagation at half the true medium velocity. The ZSR imaging times are correct but only the A amplitude divergence is compensated along the upgoing raypath; the r A term along the downgoing raypath is effectively ignored. The excitation-time imaging s condition migration can be implemented by choosing the integrand weight of (4.1) equal to the source amplitude component: *M4 = A . s 4.6 The Geometric Imaging Condition U\x , 0 t) * D\x , 0 t)] I \D\X , 0 t) • D\x , 0 t) «'=0 W = A r ly/ArT g s s r A stable and efficient geometric approximation to the dynamic imaging condition is now proposed: R g = ff/D\ The terms if and & are identical to U and D except that the dynamic amplitude factors A and A are replaced with their constant-velocity divergence r s approximations l/r> and l/rs respectively. This avoids singularities due to the downgoing wavefield D while retaining a very useful approximation to the true average wavefield attenuation. Since the amplitudes are approximately correct, and the image phases are identical to those of the dynamic imaging method, the geometric imaging condition presents a very simple and robust method for accurate migration imaging and amplitude recovery. As a result, the geometric imaging condition can be classified as a reflectivity inversion method. The geometric imaging condition migration is implemented by setting Wg = 60 A r l\/A r s s r r in the generalized Kirchhoff equation (4.1). Furthermore, a simple geometric estimation of reasonable minimum and maximum r values along reflection raypaths provides an upper and lower bound on the resulting amplitude ratio r /y/f^. These bounds constrain the amplis tude ratio to any desired degree of stability at the gradual expense of 2-D true-amplitude recovery. This allows for an optimal combination of the competing demands for a high quality amplitude-balanced image and the accurate quantitative recovery of true-amplitude reflectivity distribution values. This method is also relatively efficient since the computationally-intensive dynamic amplitude variations need not be raytraced, but rather reduce to simple algebraic calculations of the radii r and r . Furthermore, for an arbitrary 2-D migration model, the geometric amplitude r s factors need be calculated only once, whereas the other imaging conditions require extensive raytracing of the dynamic amplitude maps A and A for each surface (shot and receiver) r s location. 4.7 The Kinematic Imaging Condition • R (Xo) = • ( X r , * ) | * w k In the extreme case that we constrain the variation of the amplitude ratio r /y/ry such that s it is constant for all subsurface locations, then the geometric imaging condition reduces to the kinematic imaging condition: = ^(x,., t = x + x ). The kinematic condition obtains s r a reflectivity estimate by evaluating a modified version of the surface-recorded wavefield ^(xr, t) at the total shot-subsurface-receiver traveltime: t = T(x0;xs)-K(xr;x0). The reflected data are modified such that 61 *(x,o r = ^^ar*F(xr,o. As the name implies, the kinematic imaging condition uses only the raytraced imaging times xs and xr , and does not explicitly use any amplitude information (other than the recorded r trace amplitudes), to perform the migration. Setting M4 = A /y/A^ s removes all dynamic amplitude dependence from the generalized solution (4.1), and thus implements the kinematic imaging condition migration. Although the kinematic imaging condition does not utilize any wavefield amplitude information, the migrated depth images tend to recover true relative-amplitude reflectivity estimates in regions of 1-D to mildy 2-D velocity variation. This aspect is discussed later in greater detail. As a result, the kinematic imaging condition is proposed as an efficient true-amplitude migration method in regions of 1-D velocity variation. If true-amplitude recovery is not required, then the kinematic migration is recommended as a practical reflectivity imaging tool in terms of computational efficiency and image quality. The dynamic, geometric, and kinematic imaging conditions can be thought of as single members in a series of true-amplitude reflectivity imaging-inversion conditions. The series represents a tradeoff in qualitative reflectivity image resolution, versus quantitative reflectivity amplitude accuracy. The dynamic imaging condition is an end member of the series representing the best in amplitude accuracy, but the worst in image stability. The kinematic imaging condition tends toward the opposite end of the tradeoff curve, in that it possesses extremely good image stability, at the expense of amplitude accuracy. The geometric imaging condition occupies the middle ground, and represents an optimal compromise of both reflectivity image stability and amplitude accuracy. 62 Chapter 5 Mathematical Analysis of Reflectivity Amplitude Recovery In this chapter we seek a physical intuition for the first order reflectivity amplitude recovery behaviour of the various imaging condition migrations. In particular, it is useful to understand the specific approximations and circumstances for which they approximate the dynamic imaging condition amplitude recovery, which we subsequently assume to be correct. First, however, we discuss some aspects of the rather ubiquitous term velocity model. This discussion will be of subsequent help when we examine velocity model considerations for true relative-amplitude reflectivity recovery among the various imaging conditions. 5.1 Discussion on "Velocity Models" In general, we can separate the true subsurface velocity distribution, c , into a reference true velocity component c ,, and a residual velocity component Ac, such that true(*o) C Usually, the reference velocity c velocity field c true = Cre/r(XD) + AC(X 0 ). (5.1) is specified to incorporate the major features of the true while maintaining some desirable properties such as spatial smoothness 63 and small spatial gradients. This situation is applicable to the case of a strongly scattering medium, in which c , is the first order approximation to the true field c , and the residual true field Ac represents a second-order variation such that: true C ~ ref C and \te\<\c \. nf (5.2) In the case of a weakly scattering medium, the residual velocity field is on the order of a perturbation to, and of much smaller magnitude than, the reference field: Ac —> &, such that ^true — '-ref and l & K l c J . (5.3) In a scattering medium, the reflected wavefield is due primarily to scattering of the incident wavefield by the residual velocityfielddistribution 6b or Ac. The reference velocity does not contribute to the reflected energy since its spatially smooth and slow-varying properties rule out the possibility rapid velocity contrast required for reflection of incident energy. However, the reference velocity is important with regard to the propagation of incident and reflected wavefields in the subsurface scattering medium. In weakly scattering media, the 64 propagation of both incident and reflected energy is strongly controlled by the reference velocity c ,. The residual velocity plays no major role in weakly scattered propagation since & <C c . In strongly scattering media, the reference velocity is a significant factor in wave f propagation, however, variations in Ac are also large enough to affect the propagation of the incident and reflected wavefields. In the context of migration, we need to specify a migration velocity model c mig in order to accurately reconstruct the incident and reflected wavefields. Ideally, we should define m mi c ~ truec g practice, we do not know the true velocity model c true a priori, and so usually estimate the reference velocity component c of the true velocityfieldsuch that nf W * o ) - creixc). (5.4) If the true subsurface model can be classified as a weakly scattering medium, satisfying conditions (5.3), then the migration velocity c migration result, provided that c mig mig can be expected to provide an excellent is a reasonable estimate of c as implied by (5.4). In this nf case, the migration wavefield reconstruction integrals can be implemented accurately with c mig since the true wavefield propagation is governed primarily by c , and not 8?. Hence, nf with the correct wave propagation accounted for, the application of the imaging condition as a function of the correctly reconstructed upgoing and downgoing wavefields will produce an excellent estimate of the residual velocity distribution Sc that gives rise to the reflected wavefield. If the true subsurface is strongly scattering, as suggested by conditions (5.2), some difficulties arise. In this case, the true wave propagation is governed by both c tion of the migration wavefield reconstruction integrals with c 65 mig and Ac. Implementa- ~ c will deteriorate the nf migrated result as the approximation c mig ~c true deteriorates. The reconstructed wavefields will not be correctly propagated, and hence will not correctly focus near reflective regions caused by the residual velocity variation Ac. The migrated reflectivity images will be skewed and distorted in proportion to the propagation raypath direction, traveltime, and amplitude errors incurred during wavefield reconstruction. The deleterious situation caused by a strongly scattering subsurface can be improved by trying to estimate and incorporate some component of the strong residual velocity Ac into the migration velocity model c . However, this may be difficult to do in an a priori fashion. mig A residual velocity estimation scheme may be implemented by an iterative migration method, in which the residual velocity Ac and migration velocity model c mig are continually estimated and updated respectively between successive migration iterates. This process may prove to be unstable in that convergence to a migration velocity model that maximizes the coherency of the reflectivity image may not be attainable. Furthermore, such an iterative process may be extremely computationally intensive. 5.2 Geometric Migration Now we turn our attention to the mathematical analysis of the amplitude recovery behaviour of each imaging condition. The new geometric imaging condition is expected to give the closest approximation to the theoretically correct dynamic imaging condition, and hence we examine itfirst.Let us begin with the geometric version of the generalized migration solution (4.1), 66 R (xo) = 2 JJ?^^M{x ,%\t)dxrdxa, g (5.5) r where M(x , xs; t) are the recorded surface data modified by the half-time-derivative operator r and evaluated at the imaging times: I/2 M(x , xs; t) = 3, T(xr, xs, t) | < = T s + T r . (5.6) r The correct dynamic version of (4.1) is given similarly by: R (x ) = 2 JJ^°-?jt^ d M(xr, xs; t) dx dx . 0 r s (5.7) It is evident that the geometric migration given by (5.5) is a good approximation to the correct dynamic migration (5.7) when A r ~ l/r r and A s ~ l/r . s In general, (5.5) approaches (5.7) when the WKBJ Green's function amplitude factor A is approximately equivalent to the constant-velocity geometric-divergence factor: At ~ Mn, 67 (5.8) where A;(X0;XJ) is the WKBJ amplitude response evaluated at any subsurface location x0 having been excited by a unit amplitude impulse located at Xj, and ri is the radial (straightray) distance from xi to x0. Condition (5.8) holds identically for a constant velocity model, and we will now show that (5.8) is a good approximation for models with weak 2-D velocity variation. Let us assume that, in general, most true subsurface velocity distributions can be represented, at least to first order, by a 1-D reference model, as discussed in the previous section. This assumption is often quite valid, since velocity increases with depth into the Earth's subsurface due to the horizontal deposition and vertical compaction of sedimentary layers in the lithification process. The residual velocity component allows the subsurface to possess 2-D variation in a weak or strong scattering sense. The amplitude recovery of the geometric imaging condition is dependent upon the choice of a migration velocity model. As per the discussion in the previous section, and given the generalization above, we henceforth assume a 1-D migration model. The migration amplitude recovery is controlled by the WKBJ amplitude factors A;. The A; amplitude variations, in turn, are largely dominated by the geometric divergence component of the total wavefront amplitude function, as will be demonstrated in later examples. This is because the amplitude attenuation associated with geometric divergence ranges over several more orders of magnitude than the effects of transmission loss or Q attenuation, for the class of subsurface models examined in this text. Hence, to examine the first order behaviour of the geometric migration amplitude recovery, we need to examine the geometric divergence factor in a 1-D medium. O'Brien and Lucas (1971, p.13) give a discrete expression for the geometric amplitude wavefront divergence Arj in the n layer of a 1-D medium composed of N constant-velocity horizontal plane layers: th 68 ( n<N \ /n<N \ XJWcij i^hckcos QJcrcos Q J, 2 (5.9) 2 k where ^ is the length of the ray in the k layer, is the velocity in the k layer, and 0^ is th th the angle between the ray and the vertical in the k layer. If we assume that, to first order, th the downgoing source energy (or the upgoing reflected energy) travels vertically through the subsurface, as is the general case for deep seismic reflection data, then cos 0*. — 1 and (5.9) reduces to: n<N A \z ) D n ~ £Z*c*/c,. (5.10) If we consider the true 1-D reference velocity to be comprised of an infinite number of thin constant-velocity plane horizontal layers, i.e., as N —> °°, then 1^ —> dl, —» c(z) = c , mig and remarking that the constant factor C\ reduces to the velocity c at the source region, then s (5.10) reduces to a continuous relative-amplitude 1-D geometric divergence factor: fl(Zo) i A \z ) D 0 ~ - / c(z)dl. s Jl(z )=0 (5.11) c s The integration (5.11) is performed along the raypath length I from the source depth z to s the subsurface point z . In a constant velocity medium, the raypath is a straight line of total 0 length l = r, hence (5.11) reduces to AQ = 1/r as expected. However, we consider a 1-D subsurface velocity variation. Continuing thefirst-orderanalysis, we assume that c(z) is approximately linear with depth, as is common in the Gulf of Mexico, for example, (Sheriff and Geldart, 1982): 69 (5.12) c(z) ~ Co + az. The velocity at the surface, c0, may assume an average value of about 2 km/s, and the velocity gradient a may be reasonably specified as 0.3/s in offshore sedimentary sections (Sheriff and Geldart, 1982,1, p.87). If the gradients are not excessively strong, then dl ~ dr and dr = 2zdz/Vx + z = IcosQdz -> dz, 2 2 with the previous assumption of predominantly vertically-propagating energy. In this case, (5.11) becomes A-J( ) Zo ~ - f Z°c(z)dz. (5.13) s Jzs C Substituting (5.12), performing the z-integration, and noting that c = c , reduces (5.13) to s 0 (5.14) Given the values for c0 and a above, and considering a one-way divergence distance range: (0 < r < 5) km, it is easy to show that thefirstterm of (5.14) is of leading order: A D ~ 1/r (5.15) Hence we have shown that, in a medium with a reasonable vertical velocity gradient, the 70 geometric divergence AD, and hence the WKBJ amplitude factors AR andAs, is approximately equivalent to the 1/r-type divergence of a constant velocity medium. As this is the condition required by (5.8), we can say, tofirstorder, that the geometric imaging condition amplitude recovery is a good approximation to the true dynamic amplitude recovery in the presence of realistic 2-D earth models. Furthermore, in the case that the true subsurface strongly satisfies the above assumptions of a weakly scattering medium with a mildly variable 2-D reference (migration) velocity, and the reflected data are composed of primarily vertically propagated energy, then the geometric imaging condition is expected to recover a highly accurate estimate of the true amplitude reflectivity distribution, as the migration examples will later demonstrate. 5.3 Kinematic Migration The kinematic imaging condition is a regularized version of the geometric imaging condition. The previous section demonstrates that the geometric migration is, in general, a very good first-order approximation to true dynamic migration reflectivity amplitude recovery. We now consider the conditions necessary for the recovery of true relative-amplitude reflectivity when using the kinematic imaging-condition migration. Let us begin with the kinematic version of (4.1): (5.16) Recall that M(x , xs; t) are the recorded surface data modified by the half-time-derivative r operator and evaluated at the imaging times, as given by relation (5.6). 71 We assume the true subsurface geology to be a mildly variable 2-D weakly scattering medium in the sense of conditions (5.3). That is, the reference (migration) velocity is essentially 1-D, and the residual velocity variation is mildly 2-D and of a relatively small magnitude. In this case, the data M(x , x ; t) do not depend (much) on the shot location x , but rather on the r s s shot-receiver offset within a given shot gather. Hence, the independent variable x becomes s a parameter such that M(x , x ; t) -» M(\ ; x , t). Now M(x ; x , t) can be thought of as a r s r s r s single "generalized" shot gather, attached to no specific surface shot location, since all shot gathers are identical as a function of source location over a 1-D earth. With this assumption, the integral over the shot coordinate in (5.16) reduces to an overall scale factor multiplying the receiver integral. However, since we are ultimately interested in relative rather than absolute amplitudes, the shot-integral scale factor, and the factor of (2.), can be ignored without loss of generality. Equation (5.16) reduces to (5.17) Now consider the dynamic imaging-condition migration equation: (5.18) Re-ordering (5.18) to make the shot integration the inner integral, and neglecting the factor of (2.) for the relative-amplitude formulation, results in 72 The previous assumption of 1-D or mildly 2-D geology allows us to bring the modified data, M, out from underneath the shot integral: (5.20) The inner shot integral is an integration over the shot point-sources, each of which have a 1/r-type pseudo-spherical amplitude divergence term A ( x , x ). The shot-point integration s 0 s sums over a line of closely spaced point-sources, effectively creating a line-source response. Hence, the integration sums the many l/r(x , xs)-type point-source divergence responses to 0 a single l/y/r(x ; xs)-type cylindrical line-source divergence. The line-source integration 0 converts the shot-variable WKBJ point-source amplitude A ( x , x ), to the shot-independent s 0 s WKBJ line-source amplitude factor \ZA (x ; x ). Hence, Equation (5.20) integrates to s 0 s (5.21) Both of the WKBJ amplitude terms A ; behave like 1/r; as previously discussed, especially within the mildly 2-D assumptions of this section. Since we have assumed 1-D to mildly 2-D geology, the generic shot-gather source-to-subsurface distance r of the total reflection s path is approximately equal to the receiver-to-subsurface reflection distance r . This in turn r makes A approximately equal to A , and the WKBJ amplitude ratio in (5.21) reduces to r s unity; A~ 1/rs> , A~ 1/rr> , s r 73 r s r, ^ r which implies ~ 1, resulting in: (5.22) This result is obtained from the true dynamic migration equation (5.18) under the assumptions of a mildly 2-D, weakly scattering subsurface medium. Comparing this result (5.22) with the kinematic migration equation (5.17) obtained for the same assumptions, demonstrates their mutual equivalence, (5.23) under the assumption of 1-D to mildy 2-D velocity structure. Thus, the kinematic imaging condition approximates the dynamic and geometric imaging conditions, under the assumption of a mildly 2-D weakly scattering subsurface, and hence tends to recover true relative-amplitude reflectivity in those circumstances, as the migration examples later demonstrate. This analysis also helps to explain the apparent amplitude-balancing property of kinematic 74 migration, as later demonstrated. If reflected energy originates from steeply-dipping 2-D reflectivity structure, then the source and receiver amplitude factors A and A will be s r significantly different due to distincdy unequal raypath-distances rs and rr . This implies that the WKBJ migration amplitude ratio of the dynamic or geometric imaging condition, y/A^lA , s will have a significant dynamic range from values much less than, to values much greater than, unity. However, in this strongly 2-D situation, the kinematic migration neglects the true WKBJ amplitude-variation, and regularizes all extremal values to unity. This has the effect of balancing the amplitudes in the resulting kinematic image in regions of strong 2-D structure, yet preserving the correct relative-amplitude variation in areas of 1-D or mild 2-D structure. The net result is an imaging method which possesses a desirable compromise of the conflicting requirements for image stability and amplitude recovery. 5.4 Excitation-Time Migration The excitation-time imaging condition has been used for prestack migration in afinitedifference implementation by Chang and McMechan (1986) and in a Kirchhoff formulation by Etgen (1987). Also, the excitation-time imaging condition implicitly forms the basis of many of the current post-stack migration algorithms in use today: the f-k domain phase-shift methods of Stolt (1978) and Gazdag (1978), the x-t domain Kirchhoff diffraction summation methods of French (1975) and Schneider (1978), and the w-x domain finite-difference method developed by Claerbout (1976), for example. In very basic terms, this imaging condition correctly compensates for amplitude attenuation during the reconstruction of the surface-measured upgoing reflected wavefield, but does not compensate for the downgoing source wavefield attenuation, as detailed in the previous chapter. Hence, we anticipate an irrecoverable amplitude error in the imaged reflectivity, which can be described as a function 75 of the uncompensated source wavefield WKBJ amplitude factor A, The excitation-time migration equation is given by (4.1) for *W = A : e -R (x ) = e 2 cos 6r 0 s A M(x , x ; t) dx dx. r r s r (5.24) where M are the modified data given by (5.6). We will assume the subsurface to be a mildly 2-D weak scatterer, as was done in the kinematic section previously. Thus, the migration velocity model is essentially 1-D, and the residual reflectivity structure is mildly 2-D. In this case, the recorded data do not depend much on the shot coordinate, reducing the shot integration to an overall constant scale factor. Since we are practically concerned with relative amplitudes only, we can ignore the integration constant, and the factor of (2.), reducing (5.24) to: (5.25) We desire to compare this result with the correct dynamic or geometric migration result. The dynamic equation was derived under identical velocity assumptions in the previous section, and is given by Equation (5.21) as (5.26) where the A factor has been brought out from underneath the integral for clarity, and does s not depend on the receiver coordinate xr. Comparing (5.25) to (5.26), we see that the amplitude recovery of the excitation-time imaging Condition is in error by a factor of \fA~ : s 76 R (x ) ~ e \/A (x ;x )R (x ). Q s 0 s d (5.27) 0 In previous sections we have shown that the WKBJ amplitude A can be usefully approximated by a 1/r-type attenuation factor: A ~ l/r . s s In this case, (5.27) can be interpreted as R (xo) ~ e r ] _R (x ). v s(,x0; xs) d (5.28) 0 This clearly illustrates the nature of the reflectivity amplitude error inherent in the excitationtime migration. The excitation-time amplitude recovery error increases as the square root of the radial distance (depth and offset) of the reflector increases, with respect to the shot location. This implies that only the shallowmost reflectors positioned vertically beneath a shot location may have meaningfully recovered amplitude estimates. The deeper reflectivity, which is often of greater interest, may be imaged successfully, but the relative amplitudes within the image are definitely unreliable. Relation (5.28) was developed for an essentially 1-D earth, and for the more general 2-D case, the anticipated amplitude recovery error would be naturally much more severe. As an ancillary remark, relation (5.28) suggests methods of partially correcting the amplitude recovery error inherent in the excitation-time imaging migration. A partial 2-D correction could be attempted by applying an individually-designed y/r (x ) weighting matrix to each s 77 0 migrated shot gather before stacking into the final composite estimate of R (*o)- Also, if the e downgoing source energy can be shown to propagate nearly vertically, as is the case for near offset or very deep reflection surveys, then r ~ (z0 -zs ) and a partial 1-D correction can be s implemented. This 1-D correction could also be implemented prior to migration and applied directly to the shot gather traces. This could be achieved by converting the depth-variable gain function to a time-variant trace gain function, using the migration velocity model for the implicit depth-time conversion. Such amplitude correction schemes may be of practical use since most of the post-stack migration algorithms implicitly use the excitation-time imaging condition. The above analysis quantifies the expected amplitude recovery error of these algorithms, and suggests approaches for compensating and minimizing that error. Such corrections, however, are somewhat ad hoc since (5.28) is only afirst-orderapproximation based on strict velocity assumptions. Also, such a correction becomes increasingly inaccurate as the reflectivity structure (residual velocity) becomes increasingly 2-D. Finally, any serious attempt at reflectivity amplitude compensation must take into account the amplitude and waveform variation along each individual source-subsurface-receiver reflection raypath. For these reasons, the dynamic and geometric imaging conditions still maintain superior amplitude recovery compared to the corrected, or uncorrected results of an excitation-time migration. 5.5 Crosscorrelation Migration The crosscorrelation imaging condition forms the basis of the pioneering work of Jon Claerbout (1971) in seismic wavefield migration. In an attempt to stabilize the volatility of the dynamic imaging condition, Claerbout decided to sacrifice amplitude recovery for image 78 quality. Crosscorrelation migration is indeed inherendy stable, but its amplitude recovery is so poor that it risks washing out any significant image information as well. In basic terms, crosscorrelation migration has the same amplitude characteristics as excitation-time migration, except an additional amplitude attenuation is imposed by the crosscorrelation source amplitude term A . For this reason, we expect crosscorrelation migration to have an amplis tude recovery error that can also be expressed as a function of the WKBJ source amplitude A , but even more severe than the excitation-time imaging method. s The crosscorrelation implementation of the generalized Kirchhoff migration equation (4.1) can be obtained by setting W =A .: 2 C -Rc(xo) 2 c / / = °^ A \fA~ M(x , x ;t)dx dx , r s r r s r (5.29) s y/Cr JJA where M are the modified surface seismic data given by (5.6). Re-ordering (5.29) to make the shot integration the inner integral, and neglecting the factor of (2.) for the relative-amplitude formulation, results in R (x ) = J c y/A • J A M(x , x ; t) dx • dx . 0 r s r s s r (5.30) We assume a mildly 2-D weakly scattering subsurface again, which implies an essentially 1-D migration velocity model. In this case, the data M do not depend much on the shot coordinate x , and can be brought out from underneath the shot integral: s C / O S Q — —j==- y/A M(x ;x , r r 79 s r *)• JAsdxs-dx . r (5.31) As in previous sections, the shot integration converts a line of point-sources into an effective line-source response, and so the point-source amplitude divergence 1/r is integrated to 1/^/r. Hence, A —» >/A7, and (5.31) reduces to s R (x ) ~ y/A c 0 s / — ~ y/A M(x -,x ,t)dx . r r s r (5.32) Comparing (5.32) with the dynamic migration equation (5.26) developed for the same velocity assumptions, we obtain a relation quantifying the reflectivity amplitude recovery error inherent in the crosscorrelation imaging method: R (x ) ~ A (x ;x )R (x ). c 0 s 0 s d 0 (5.33) If we again make the correspondence of A ~ l/r , (5.33) can be interpreted as s s Rcfro) i? d (x 0 ). (5.34) Relation (5.34) demonstrates the severe attenuation in reflectivity amplitude recovery expected from the crosscorrelation imaging condition migration. The amplitude deterioration increases as the radial distance (offset and depth) of imaged reflectors increases with respect to the shot location. Note that this amplitude recovery is more severe than the excitation-time migration by a factor of y/A~ . Relation (5.34) suggests amplitude correction schemes similar s to those discussed in the previous excitation-time section, and as such, will not be further elaborated upon here. 80 5.6 Summary This chapter has examined the reflectivity amplitude recovery behaviour of the various migration imaging conditions from simple, but reasonable, mathematical and physical analysis. Ranked from best to worst, in terms of amplitude recovery accuracy, the imaging conditions are: geometric, kinematic, excitation-time, and crosscorrelation. The new geometric migration is shown to be capable of good amplitude recovery in regions of moderate 2-D velocity structure, and excellent amplitude recovery in regions of mild 2D structure. The new kinematic imaging condition is shown to have very good amplitude recovery in essentially 1-D regions of weakly scattering velocity structure. The excitation-time imaging condition, which forms the basis of many current post-stack as well as prestack migration algorithms, is shown to have a shot gather migration amplitude deterioration proportional to the square root of the WKBJ source amplitude term: A / A J . The crosscorrelation imaging condition supports an even less desirable shot gather migration amplitude deterioration, being linearly proportional to the WKBJ source amplitude term: A . S Finally, as an anticipatory remark, the dynamic imaging condition, which is theoretically expected to recover correct reflectivity amplitudes, suffers in practical implementation from its inherent instability, as will be shown in later examples. Regularization of the dynamic imaging condition can produce good image quality, but the volatility incurred by uncertainty in the estimation of the stabilization parameter degrades the accuracy of the amplitude estimates more than the limiting assumptions of the geometric imaging condition. Hence, in terms of practical implementation, and certainly efficiency, the geometric imaging condition possesses superior performance over all the imaging conditions, including the dynamic 81 imaging condition. Chapter 6 Zeroth Order Asymptotic Ray Theory 6.1 Introduction This chapter is an attempt to present an understanding of the ray theory considerations that fundamentally underlay the raytracing implementation of the generalized Kirchhoff-WKBJ migration equation (4.1) derived earlier. In this thesis, the W K B J Green's functions are evaluated by raytracing techniques. Raytracing algorithms are useful for this application since the direct arrival ray family can be easily specified, and the subsequent traveltime and desired amplitude calculations can be performed efficiently. This is in direct contrast to finite difference evaluations, for example, which are more computationally intensive and, as full-wavefield techniques, lack the corresponding degree of control in calculating separate wavefield components. However, zeroth order ART is limited to models in which material properties vary smoothly and slowly between boundaries. Migration models that violate this limitation may require an alternate approach to the Green's function evaluation. Fortunately, most migration models should respect the ART validity conditions, given the velocity model discussion in the previous chapter, and the equivalent assumptions that are fundamentally required to justify the use of the scalar wave equation itself. 83 6.2 The Ray Series Expansion We begin with an exposition of zeroth order asymptotic ray theory, and the ray equations that result from the eikonal and transport equation approximations to the scalar wave equation. Let us start with the scalar wave equation, 2 V C/(x, t) - ^-a«L7(x, t) = 0. c (x) z (6.1) To solve this wave equation we try a series solution for L7(x, t), oo L7(x,r) = ^L7A(x)F^(r-T(x)), (6.2) k=o where the functions Fk(Q satisfy the relation **(0 = k-dO, F (6.3) as given by Equations (2.3)-(2.4) of Cerveny et al.(1977). A trial solution solution of the form (6.2) is known as a ray series expansion. Note that the amplitude coefficients Uk, and the ray travel time % depend on the spatial coordinates x only. We further choose the basis functions F^(C), in accordance with relation (6.3), to represent high-frequency harmonic waves: 84 Fk (0 exp A (-ico) = (6.4) Substituting (6.4) into trial solution (6.2) yields the ray series expansion: t/(x,0 = e x p - ^ ) ) f : ^ . (6.5) Substituting the ray expansion (6.5) into the scalar wave equation (6.1) results in iQ,( ) (ico)exp- ^ £(i^( "< v2 2 V y to)<v u k (ico) + < -^) )2 + V = 0. (6.6) The zeroth order asymptotic ray theory solution is obtained by solving (6.6) for k = 0 as co — > oo. The zeroth order ART solution reduces to two important equations: the Eikonal Equation, VT(X) | 2 _ : C (X) ' 2 (6.7) and the Transport Equation, 2 A(x)V x(x) + 2Vx(x)-VA(x) = 0, where A - U^. 85 (6.8) The solution of the eikonal equation (6.7) is the ray traveltime function x(x). The traveltime function is accurate (identical) for all terms in the ray series expansion since T(X) does not depend on k. Hence, the zeroth order ART solution is perfectly accurate with respect to travel time. The solution of the transport equation (6.8) is the leading order ray amplitude function A(x). This amplitude is strictly valid only for k = 0, but is generally accurate enough for most applications. Notable exceptions include head wave and diffraction amplitude calculations. In effect, A(x) is a generalized geometrical spreading function for heterogeneous media. The zeroth order ART solution essentially corresponds to the geometrical optics solution of WKBJ theory, since, for k = 0 in the ray series expansion (6.5), i<0( (S)) U(x, t) = A(x) exp- ^ , (6.9) where, again, x(x) are the ray traveltimes obtained by solving the eikonal equation (6.7) and A(x) are the leading order ray amplitudes given by solving the transport equation (6.8). Comparing (6.9) with the WKBJ Green's function (3.13) of previous sections, we see that the zeroth order ART traveltimes and amplitudes are in complete agreement with the KirchhoffWKBJ migration theory. Hence, evaluating the WKBJ Green's functions by zeroth order ART texhniques is completely justified. The validity conditions for the ray method outiined above, in terms of slow-varying, smooth material properties, are discussed in detail by Cerveny et al. (1977), Chapter 8. 86 6.3 Solving the Eikonal Equation The eikonal equation (6.7) is afirstorder nonlinear partial differential equation (PDE): I VT(X) I= 2 2 1 c (x) It can be converted to a system offirstorder linear ordinary differential equations (ODE's) using the method of characteristics (Cerveny et al., 1977, p.58). The resulting system of ODE's can be solved numerically with a Runga-Kutta algorithm in discretized velocity models, for example. For a heterogeneous 2-D velocityfield,c = c(x, z), the raytracing system is given by three ODE's (Cerveny et al, Equation (3.16)): dx . — = c sin 9 dx dz — = c COS0 dx n dQ dx — = -c cos9 + c sin9, x z (6.10) where c = d c(x, z), etc., and 9 is the angle the ray makes with respect to the vertical axis x x k. For a heterogeneous 1-D velocityfield,c = c(z), two raytracing ODE's must be solved: dx Tx = c 87 2 p and 2 ^ 2 = ±cy/l -p c , (6.11) where p is the ray parameter defined as sin6(z) p = " ^ p (6-12) Equations (6.11) can be rewritten into a system for x{z) and t(z): dx dx dx dz dx dz cp 2 2 A/1 -p c and dxdx _ dx _ 1 2 2 dx dz dz c A/1 -p c (6.13) Integration of Equations (6.13) results in z x(z)= f pc(z')dz' p K / [2 2 + *(zs) (6.14) ./*• V1-P C (2') and = / 2 2 7* c(*') y/1 -P C (Z>) + x(zs ) (6.15) These are the familiar parametric ray equations for offset and traveltime as a function of depth in a 1-D medium, as per Cerveny et al. (1977, equations (3.20')), or AM and Richards 88 (1980, p.404), for example. They can also be derived independently from geometrical considerations. Equation (6.14) describes the path of a ray through the 1-D velocity medium c(z) as a function of the ray parameter p, which is constant for each specified initial ray take-off angle 0. Equation (6.15) gives the total traveltime x along the raypath described by (6.14). These equations are solved by an iterative Newton-Raphson approach in a 1.5-D raytracing algorithm developed to evaluate the WKBJ Green's function traveltimes x and x for a s r 1-D migration velocity model. The details of the iterative solution and examples of the WKBJ Green's function traveltime maps necessary to implement the generalized Kirchhoff migration equation (4.1) are discussed later. 6.4 Solving the Transport Equation We desire a ray amplitude solution A(x) to the transport equation (6.8): A(x)V x(x) + 2Vx(x)- VA(x) = 0. 2 To facilitate a solution, we will transform from the present Cartesian coordinate system to a new ray-centred coordinate system. 89 6.4.1 Transforming to Ray-Centred Coordinates To assist in determining a zeroth order ART amplitude solution A(x) to the transport equation, we make the following transformation from Cartesian to ray-centred coordinates (Cerveny et al., 1977, §2.4): x = (x,y,z) = x(Yi,y2,x), (6.16) where T is the (time) coordinate within a given ray, and Yi, Y2 are orthogonal to T and lay within the wavefront surface. For afixedray traveltime % , x{y\, Y2, tr) sweeps out a r wavefront surface in the spatial medium. Expressing the Laplacian operator in the new curvilinear ray-centred coordinates we obtain (Babich and Alekseev, 1958): V * = J T S J ] * { \ B J ] ) - ( 6 - 1 7 ) where J is related to the Jacobian of the transformation from Cartesian to ray-centred coordinates (Cerveny et al., 1977, Equation (3.25)): J = c J (6.18) The term \Tij\ is the determinant of the coordinate transformation partial derivative matrix, 90 (6.19) where Xi = (x, y, z) and ji = (yi, Y2, t). The "Jacobian term" J is also related to the crosssectional area da of an elementary ray tube by Equation (3.23) of Cerveny et al. (1977): da = (6.20) \J\dyidy . 2 Re-expressing the o\x term of (6.17), dx dx 9Tx I = 1 dx dx =1 vx r = c, (6.21) the Laplacian becomes z 1 d (J Vx = — — , , . Jc 9T \c J (6.22) Also, the transport equation dot-product can be evaluated in ray-centred coordinates: V A Vx - — — - ^ i ^ L ^ l - ^ (Hl\ - 1 — 2 dxi dxi dx dxi dxi dx \dxi J c dx' 2 (6.23) Combining (6.17)-(6.23), the transport equation in ray-centred coordinates becomes 2d,A c 2 + A_ d Jc T fj_ \c 91 0. (6.24) 6.4.2 A General Solution to the Transport Equation Now that we have transformed the transport equation (6.8) into ray-centered coordinates, a general solution is more easily facilitated. Rearranging (6.24), A 2 (J/c) (6.25) and integrating both sides of (6.25) results in In A = - i In (- } + do. 2 \c (6.26) Exponentiating both sides results in a solution of the transport equation for A, Ao " ]fj' (6.27) where A 0 is the amplitude of the initial source, A D = A(x = xs). Equation (6.27) is related to the geometric divergence of the cross-sectional area of a ray-tube in ray-centred coordinates. In an analogous manner, although somewhat more complex, Babich and Alekseev (1958, Equation (4.5)) solve the elastodynamic wave equation for the elastic counterpart of (6.27). They obtain the result, corroborated by Cerveny et al., 1977, Equation (2.30), as A = -^=, y/jpb" 92 (6.28) where \)/ is a source amplitude function, p(x) is the material density function, and c(x) may 0 take the value of oc(x) or p\x) for compressional or shear amplitudes respectively. 6.4.3 Evaluating the Source Term To determine the source amplitude term \(/ , it is necessary to examine the wave solution in 0 a small spherical homogeneous volume centered about the source location x . Cerveny et al. s (1977, Equations (2.40)-(2.41)), give the results in an elastic medium for a unit-amplitude impulse point-source: = g(®s, §s) A / c p s i n 0 , (6.29) Vo = #(9 ) VeTpI, (6.30) s s s and line-source: S where g is a function describing the directional characteristics of the source, and the s subscript indicates evaluation at the source location x . For the seismic migration application s at hand, we assume isotropic shot gather point-sources and point-receivers. In this case, g(&s, <t>s) = 1, and (6.28) reduces to (6.31) At this point we note that thefirstfactor of the amplitude solution is related to the geometric divergence of the raypaths, and the second factor is an impedance compensation between the 93 source impedance and the ray endpoint impedance. 6.4.4 Evaluating the Jacobian Term The factor J in the amplitude solution for A is related to the Jacobian of the transformation from ray-centered coordinates to Cartesian coordinates. Geometrically, J is related to the change in the ray-tube cross-sectional area resulting from a small change in the ray take-off angle 9S. In general, J can be factored into two components such that (6.32) J = where J\\ represents the in-plane spreading, J|l e (ixk), (6.33) and J±_ represents the out-of-plane spreading, (6.34) The in-plane spreading J\\ is given by Cerveny et al. (1977, Equation (3.60)) as (6.35) 94 where the partial derivative is evaluated for fixed z. For 2-D and 3-D media, the partial derivative in the Jy expression (6.35) is often approximated by ray differencing. For a given ray, a ray-tube is created by raytracing adjacent rays for very small increments dQ of the ray take-off angle 0S. The change in offset or range dx s is then measured across the ray-tube cross-sectional area. The cos 6 term is the projection of the dx quantity onto the ray-tube cross-sectional area. Note that if J\\ —> 0, e.g., near a turning point (cos 6 —> 0), or near a caustic (dx —» 0), then the amplitude solution for A is invalid since A —> °°. Hence, zeroth order ART breaks down, and higher order terms of the ray series expansion (6.5) in the wave equation (6.6) are required to accurately predict amplitudes in these troublesome regions. For the migration problem at hand, this means that upgoing or downgoing segments of reflection raypaths that contain turning points or caustics will not be correctly imaged by the zeroth order ART evaluation of the WKBJ functions. Practically speaking, this migration error corresponds to regions of large gradients in a migration velocity model. In such cases,finitedifference evaluation of the Green's functions may be necessary. Fortunately, a migration velocity model is generally required to be as smooth as possible while still representing the estimated subsurface velocity structure, as per the discussion of a previous chapter, and hence often satisfies the zeroth order ART limitations. For 1-D media, the derivative in (6.35) can be evaluated analytically. From the ray-equation solution to the eikonal equation (6.14), we have 95 where p is the ray parameter defined by (6.12) as P _ sin 0(z) ~ ~cW' Rearranging the required partial derivative (6.36) we can evaluating thefirstterm from the ray equation (6.14) for x(z), c(z')dz' cos 0(z')' (6.37) 3 Zs where cos 0 is defined by the ray parameter as cos 0 2 2 = Vl-sm 0 = V l -P c 2 (6.38) Also, we can make the evaluation cos 0 dp\ Combining (6.36)-(6.39), S we obtain an analytic expression for the JM derivative, 96 (6.39) \dQsJ c z 3 cos J s Zs 6(2') Equation (6.40) agrees with Equation (6) of Newman (1973), except for his factor of (2.) relating to two-way raypaths. Combining (6.40) with (6.35), we obtain the full 1-D analytic expression for the in-plane spreading factor, j j D , _ cos8,cos8(*) c 11 f#l_„ J s . (6 41) 3 Z s cos 9(z') The out-of-plane spreading factor, J±, is given by Equation (5.22) of Cerveny and PSenCfk (1979) as an integral along the length I of the ray from the ray initial point Z; to the ray endpoint If, ^ = sine, f c s lf c{l)dl ( 6 4 2 ) hi For general 2-D or 3-D media, this quantity is again readily available from the raytracing kinematics. The total general Jacobian term is defined by (6.35) and (6.42) as JO. U' J = '*1*2!M) (* ( | ( / c(l)dl) . For the 1-D medium, Jj. is given simply from (6.42) as Cs Jzs COs6(z') 97 (6.43) The in-plane and out-of-plane spreading factors J | | and Jj_ given by (6.41) and (6.44) can be combined to yield an analytic expression for the total spreading factor (6.32) in a 1-D medium, sin 0 cos 0 cos 0(z) S c(z') dz' cos 6(z') S (6.45) As a final remark, the J equations listed above may need to be reparameterized as a function of offset x instead of depth z for near-horizontal ray propagation. This occurs for the situation in which 0(z) -> 0°. 6.4.5 Energy Partitioning So far we have been considering the amplitude of a ray traveling through a smooth and slowly-variable heterogeneous medium. We also need to consider the partition of energy amplitude effects at discontinuous layer boundaries. This includes incorporating the Zoeppritz reflection-transmission amplitude coefficients and matching impedance jumps from the divergence amplitude A in (6.31) at a boundary: (6.46) where the impedance terms £ = pc are defined as: £s at the ray starting point (source), C, at r the ray endpoint (receiver), of emergence at the i th at the point of incidence at the i th layer, and ^ at the point layer. Also, n is the number of boundaries encountered by the ray, and Zi is the Zoeppritz displacement amplitude coefficient at the i th 98 boundary. Unfortunately, the amplitude divergence as written now also exhibits a non-physical amplitude discontinuity across layer boundaries. Cerveny et al. (1977) add a compensation factor so that the resulting amplitude expression becomes A = A D A P , (6.47) where AD is the geometrical divergence amplitude factor given by and Ap is the energy partition amplitude factor given by ^•te)"n(g)% In general, the Zoeppritz coefficients Zi can accommodate reflection, transmission and conversion at an interface between two elastic media. For the migration application at hand, we need only the transmission loss component across each layer. Since we will eventually consider acoustic WKBJ amplitudes, we will use the specular transmission coefficient, T;, which is a plane-wave, incident-angle dependent acoustic transmission coefficient, JI. 1 = 2p£Cf/cosef pi c /cosQ + pjCi/cosGj" +1 i+l i+1 99 6.4.6 Q-Attenuation As an option, we may desire to calculate an approximate Q-attenuation amplitude factor, AQ. Given the equation for Q-attenuation as a function of distance traveled in the propagation direction, i.e., along the ray, from Aki and Richards (1980, p. 169), we may write 71 AQ = Y[ exp(-7t/ Zj/ciQj), i=i 0 (6.51) where f is the dominant frequency of the propagating pulse, U is the length of the ray 0 segment, c; is the average velocity, and Q{- is the average Q value along the raypath within the i th layer. This relation (6.51) is approximate since it ignores the true dispersive nature of Q-attenuation, and assumes that the source wavelet is spectrally narrowband and centered about a dominant frequency f . Fortunately, this approximation is often very good for seismic reflection data. Q 6.4.7 Transport Equation Amplitude Summary An amplitude solution for the transport equation (6.8) has been developed for the general case of multidimensional elastic media: A = A D A 100 P A Q . (6.52) The full solution A has an amplitude component AD representing geometric wavefront divergence due to an isotropic point-source which is defined by (6.48) as where the Jacobian J is defined by (6.43) as ' - pr^) (a (I;H • Also, part of the full solution A corresponds to the partitioning of energy amplitude component Ap arising at layer boundary interfaces, defined by (6.49) as where, in general, the Z; are the Zoeppritz elastic reflection-transmission amplitude coefficients. Furthermore, an approximate form of Q-attenuation is accommodated by the amplitude component AQ, defined by (6.51) as n A Q = II zxv(-rfok/ciQi). i=\ These are the amplitude quantities that may be computed by a raytracing algorithm in order to evaluate the Kirchhoff migration WKBJ Green's function amplitude factors AS and AR , necessary to implement the generalized Kirchhoff-WKBJ migration equation (4.1). Details 101 of the 1.5-D two-point acoustic raytracing algorithm developed for this purpose will be discussed next. 102 Chapter 7 A Fast and Accurate Two-Point Raytracing Algorithm for 1-D Media 7.1 Introduction At this point it is convenient to briefly describe the raytracing algorithm developed to compute the Green's function traveltime and amplitude variations, A(XQ) and x(x0), as required to implement the generalized migration equation (4.1). A detailed treatment of the zeroth order asymptotic ray theory that forms the foundation for this algorithm is discussed in the previous chapter. In this chapter, an iterative Newton-Raphson solution is developed and implemented for the 1-D raytracing equations previously derived. To implement the 1-D iterative ray solutions, a two-point raytracing algorithm is developed to calculate traveltimes and amplitudes from a given surface location to all subsurface points specified on a finite difference grid of the migration model. This algorithm accepts an unlimited number of horizontal layers, each specified by a depth, constant P-wave velocity, density, and Q value. The algorithm computes direct transmission times, geometric divergence, specular transmissivity, Q-attenuation, and surface-incident obliquity from the source location to the gridpoint nodes. These quantities may be calculated over the entire model or in specific subportions of the model specified by a source aperture cone or by zones of pre-critical reflection. The amplitudes obtained are 1.5-dimensional, since the algorithm 103 emulates a 1-D model illuminated by a 3-D isotropic point source. Specifically, the algorithm uses an iterative Newton-Raphson approach to solve the 1-D eikonal and transport equations in a two-point raytracing scheme. Two-point raytracing formulations are based on the premise that two ray endpoints, the initial andfinallocations, are given, and the ray that joins the two points must be found such that the eikonal equation, which embodies Snell's Law and Fermat's Principle, is satisfied in the given velocity model. The raytracing method presented is fast since the two-point iterative ray solution possesses quadratic convergence and is accurate since it avoids interpolation of a ray fan onto a finite difference grid. The main disadvantage of Newton-type algorithms is in the required selection an initial trial solution. However, for the 1-D two-point raytracing problem an initial estimate for the ray parameter can always be found for which convergence is guaranteed. 7.2 The Ray Kinematics We are interested in calculating the traveltime t along a ray from a given source coordinate (x ,z ) to any arbitrary depth point (x*,z*) in a 1-D medium, as shown in Figure 7.1. Since s s the velocity model is 1-D, the traveltime is a function of the source-receiver offset, | x -x* |, s and the depth of the subsurface point relative to the shot depth, | z -z* |. We may set 8 x , z = 0 and use x* and z* as the relative offset and depth without loss of generality. s s A common approach to the problem is one in which the ray properties at any desired endpoint (x*, z*) are interpolated from a dense fan of rays traced from a given initial ray point (x , z ) s s and a specified set of take-off angles 0S, (Aki and Richards, 1980, p.727). This method can be slow and inaccurate since many rays must be traced for each subsurface point, and if the rays are not spacedfinelyenough, interpolation error can be significant. 104 (Xs, Zs) Figure 7.1: A Two-Point Raypath in a 1-D Horizontally Layered Medium. The example 1-D medium consists of N horizontal layers each specified with a constant velocity c;. The object of two-point raytracing is to define a ray which departs from a given source location (Xs> z ) and arrives at a specified subsurface location (x*, z*), such that the appropriate ray equations are satisfied. s 105 Using the Newton-Raphson method, we solve a nonlinear integral equation for the ray parameter p* given the initial andfinalboundary conditions for the ray, (x , z ) and (x*,z*). s s The well-known offset and traveltime parametric ray equations for a 1-D medium are given by (6.14) and (6.15) as p c{z') dz' M ^\-p C {z>) = LL-JI 2 (7.1) 2 s and f dz' z Jzs C(Z') y/\ -P C (Z') 2 2 where the ray parameter p is constant along any given ray and defined by (6.12) as " sin Q(z) " (7 - > 3 Using Equation (7.1), p* is found such that X i P ) k Jl-P* c\z>) 2 ( 7 ' 4 ) Equation (7.4) can be solved iteratively by the Newton-Raphson Method which has quadratic convergence. Once we have solved forp*, we have defined the ray which travels from (x ,z ) to (x*,z*) s and whose take-off angle Gs is given by: 106 s 9* = arcsin(csp*), (7.5) where c is the velocity at the source location. We can also calculate the incident angle at s the ray endpoint by substituting c* = c(x*, z*) for c in (7.5), s 9* = arcsin(c*p*), (7.6) where c* is the velocity at the desired subsurface depth-point. This angle will have an application for determining the ray amplitudes, as discussed later. After p* has been determined from (7.4), the traveltime x(p*) can be computed directly from (7.2). Hence, we desire to solve (7.4) with the Newton-Raphson method for p*, thus determining the kinematic ray traveltime %, which is the desired traveltime along the ray from the source location (x ,z ) to any given subsurface point (x*,z*). s s 7.2.1 The Newton-Raphson Method We now turn our attention to the Newton-Raphson solution of the ray equation (7.4). Let the correct ray parameter solution p* be represented such that, p* = p + &>, 0 (7.7) where p is an initial estimate, or trial solution, that is "close" in some sense to p*. If p 0 is reasonably close to p*, then the correction term Sp is on the order of a perturbation such 107 0 that dp < p*, . (7.8) Po We can expand x(p*) of (7.4) in a Taylor Series Expansion: x(p*) = x(p ) + x'(p )$p + x"(p )bp /2\ + •••. (7.9) 2 0 0 0 If condition (7.8) is valid, then to first order, x(p*) ~ x(p ) +x'(p )8p, 0 (7.10) 0 and the ray parameter correction Sp to the initial estimate p is found to be 0 op [x(p*)-x(p )] — — . X'(p ) 0 (7.11) 0 If x(p) is a linear function, then (7.10) gives the exact solution for p*. However, (7.4) shows the true nonlinear nature of x(p). Hence, we use the linearized version of x(p) given by (7.10) to iteratively determine the successive ray parameter corrections Sp from (7.11), in order to continually update and ameliorate the initial ray parameter estimate p . We can halt 0 the iterative process when p is arbitrarily close to p*, in the sense that Q 2 \\x(p*)-x(p )\\ < e , 0 108 (7.12) 2 where e is a user specified error criterion. Obviously, the convergence of such an iterative solution depends upon the initial ray parameter estimate p . 0 We need to evaluate Equation (7.11) in order to calculate the ray parameter perturbation, 5p. To evaluate (7.11), the quantity x(p*) is given simply as the x* coordinate of the subsurface gridpoint in concern. The quantity x(p ) is given by the evaluation of (7.1) using 0 the subsurface gridpoint depth coordinate z* and the initial ray parameter estimate p . We Q also need to calculate the derivative x'(p), which is defined analytically by (6.37) as, (7.13) where cos 9 is determined from the ray parameter definition (7.3) and is given by (6.38) as 2 cosG = Vl - sin Q = A/1 -p c . 2 2 (7.14) This completes the requirements for iteratively evaluating the recursive solution for p* using equations (7.4), (7.11) and (7.12). A solution thus obtained immediately provides all desired ray kinematic information, including the ray traveltime % which can be obtained by direct evaluation of (7.2) given the solution p*. 7.2.2 A Discrete Layer Formulation It is worthwhile to note that the iterative Newton-Raphson solution of the kinematic raytracing equations (7.1) and (7.2) can be recast into a discrete formulation. This is suitable for 1D model parameterizations which consist of a sequence of N horizontal, constant-velocity 109 layers. In this case, the ray equation (7.1) for offset becomes, 71-1 x(p) = tanGj + (7.15) (z-z )tane„, B i=i the offset derivative equation (7.13) becomes, X'ip) = > s-r- + ^ — ' COS 0; f=l T"7\ > (7-16) T(P) = V — + , ^—' cj, cos 0tc„ COS 0ra f=l (7.17) 3 COS 0ra and the ray traveltime equation (7.2) becomes, where n < N is the number of the layer in which the subsurface gridpoint (x, z) is located, hi is the thickness of the i th top of the n i th th layer, c; is the velocity in the i th layer, z is the depth to the n layer, and 0j is the angle that the ray makes with respect to the vertical in the layer, defined from (7.5) as, 0; = arcsin(pc;). (7.18) 7.2.3 Aspects of Convergence Figure 7.2 represents a graphical illustration of the Newton-Raphson method for finding the ray parameter p*. Examination of Equation (7.13) or (7.16) shows that the offset derivative 110 Graph of X ( p ) t Ray Parameter Figure 7.2: A Graphical Representation of the Newton-Raphson Method for Finding the Ray Parameter p*. Ill x'(p) is always positive since 0e(O,rc/2).Hence, for a 1-D medium, the graph of x(p) is always concave upward and increasing with p values ranging between zero and l/c waje . As a point of interest, this positivity also gaurantees that solution (7.11) for the ray parameter perturbation bp is never singular. At p = 0, the ray travels vertically such that the offset x(p) from (7.1) is zero for all z, as would be expected. As p approaches the ray parameter limit l/c , max the ray approaches a horizontal (turning point) trajectory, such that the offset x(p) becomes unbounded. According to the algorithm, for an initial ray parameter estimate p , the updated estimate p\ a is found at the intersection of the horizontal line x = x*, with the line that passes through the tangent x'(p ). From the graph in Figure 7.2, it can be seen that successive iterations 0 converge to the ray parameter p* when the initial estimate p is selected such that x(p ) > x*. 0 0 It is also apparent that the region of convergence decreases as p* —> Vc . max Physically, this implies that more iterations are required tofinda ray which is near to critical incidence at the desired endpoint (x*,z*). A proof has been demonstrated by Lumley and Keys (1988), showing that the 1-D NewtonRaphson two-point raytracing algorithm sequence of Newton iterates, Pn l = Pn+ + *~ [X X(Pn)] X'(p ) ; 71 = 0,1,2,..., (7.19) n always converges to the unique ray parameter solution p* which satisfies x* = x(p*), provided that the initial ray parameter estimate p is chosen such that, Q 0 < p 0 < 1/Cmax, 112 (7.20) and (7.21) x(p ) > x*. 0 An initial estimate which always satisfies this condition is the starting value, Po = where c max (x Cmax y/(x* Xg) (7.22) 1 -Xsf+H ' is the maximum layer velocity, and H is the thickness of the layer having that velocity. In most applications, the initial ray parameter estimate p must be determined for an entire 0 suite of offsets. In such cases, the initial ray parameter may be estimated more accurately than (7.22) by using an adjacent existing ray parameter solution p*. Of course, the better the estimate p , the faster the convergence to the correct solution p*. Therefore, the most Q efficient procedure for estimating the initial ray parameter values p , is to use (7.22) for the 0 first offset, and then use the resulting ray parameter solution p* as an initial estimate for the next offset. In this manner, one can bootstrap along with increasing offset. As long as the initial ray parameter estimate satisfies the convergence-zone tests (7.20) and (7.21), convergence is gauranteed. 7.3 The Ray Amplitudes Once the ray kinematics have been determined for a given subsurface gridpoint, using the method outlined above, then the appropriate zeroth order asymptotic ray theory amplitudes can also be directly evaluated. Recall the general ray amplitude solution (6.52) for A from 113 the previous chapter, A = A D A P A Q (7.23) . The amplitude component AD represents the geometric wavefront divergence and is defined analytically by (6.48) as, , cos 0, 1=1 x 1 where the Jacobian J is defined for a 1-D medium by (6.45) as Equation (7.24) can be rearranged into a discrete 1.5-D expression for the point-source geometric amplitude wavefront divergence AD in a 1-D medium composed of N constantvelocity horizontal plane layers. The discrete version of (7.24) is obtained by discretizing the Jacobian (7.25), and after some simplification, results in Expression (7.26) agrees exactly with expression (5.9), and also that of O'Brien and Lucas (1971, p.13). 114 The amplitude component represents the energy partitioning at layer boundaries, and is defined by (6.49) as, /„ \ 1/2 re \ 1/2 where, in general, the Z; are the Zoeppritz elastic reflection-transmission amplitude coefficients. For the migration application at hand, the raytracing evaluation of the WKBJ Green's function amplitudes A(x0) for implementation of the generalized Kirchhoff migration equation (4.1), we are concerned only with the transmission amplitude component. Hence, we can substitute the acoustic specular transmission coefficient T; given by (6.50) in place of the Zoeppritz coefficient Z;, TM = 2p^/cose, pi Cj ,/cOS0i + 1 + + p;C;/cOS0;' + 1 Furthermore, an approximate form of Q-attenuation is accommodated by the amplitude component AQ, defined by (6.51) as AQ (p) = Yl exp(-7t/ Aj/QjCjCos0j). o (7.29) In summary, Equations (7.23)-(7.29) represent the amplitude quantities that may be computed by the two-point Newton-Raphson raytracing algorithm once the ray kinematics have been determined by the solution for the ray parameters p*. These zeroth order ART amplitudes are required in order to evaluate the Kirchhoff migration WKBJ Green's function amplitude factors As and Ar , necessary to implement the generalized Kirchhoff-WKBJ migration 115 equation (4.1). 7.4 Summary An algorithm has been developed that calculates raytraced direct transmission times in a 1-D medium using an iterative Newton-Raphson two-point method. Furthermore, once the ray kinematics have been determined, 1.5-D zeroth order ART amplitudes may also be calculated. The method is fast and allows the user to specify arbitrary accuracy to within machine precision. The method avoids the computationally slow "ray fan" method, its interpolation inaccuracies, and "dead zones" near the critical angle where ray-take-off angles become unstable. The algorithm can be easily modified to accommodate various acquisition geometries such as subsurface sources for borehole or VSP work. The traveltimes and amplitudes are especially useful for geophysical imaging and inversion procedures that require the numerical evaluation of Green's functions. This has applications in seismic migration, geotomography, generalized inversion, and non-traditional ray-theoretic problems such as acoustics and synthetic aperture radar imagery. 116 Chapter 8 A 2-Dimensional Test Model A detailed example is examined later which illustrates the performance of the true-amplitude Kirchhoff migration Equation (3.35) and thefiveimaging conditions related by (4.1). The example is based on a 2-D test model which contains six horizontal layers truncated by a dipping arc structure. Multifold synthetic seismic reflection data are generated from the test model to provide 51 shot gathers for migration. A 1-D migration model is derived from the test model, which is raytraced to provide the Green's function imaging times and WKBJ amplitudes necessary for implementation of the generalized migration equation (4.1). Details regarding the test model, migration model, and the synthetic data are discussed below. 8.1 The "Arc" Test Model Figure 8.1 displays the 2-D test model that forms the basis of the synthetic examples to be discussed later. The color legend values represent P-wave velocity in km/s. The model consists of six horizontal layers truncated by an arc structure, and is representative of sedimentary layering adjacent to a continental margin, salt dome flank, or listric fault. Table 8.1 lists the material property specifications of the arc test model. The compressional velocities vary from 1.5 km/s at the surface to 6.0 km/s at depth. Layer 6 at 117 Figure 8.1: The 2-D "Arc" Test Model. This test model is used to generate synthetic seismogram data for the prestack depth migration examples discussed later. The arc test model consists of six horizontal layers truncated by a steeply dipping arc structure. A low-velocity zone is present at 2000 m depth, which may be representative of fluid or gaseous phase reservoir potential. 118 Density (g/cc)Arc Dip Layer .Base De/Jf/z (km) Velocity (km/s) 1 2 3 4 5 6 7 0.600 0.900 1.200 1.500 1.800 2.100 3.000 1.500 2.300 3.500 4.500 5.500 2.500 6.000 1.929 2.147 2.384 2.539 2.670 2.192 2.728 82° 69° 60° 50° 38° 16° - Table 8.1: The 2-D Arc Test Model Material Property Specification. 2 km depth is a 2.5 km/s low-velocity zone, possibly indicative of hydrocarbon potential, for example. Densities vary as a function of velocity according to the Gardner relation, Equation (8.1), which is determined empirically from sedimentary basin rock samples (Gardner et al., 1974). The arc consists of a quarter of a pseudo-circle centered at 1 km on the model surface (zero depth) with a radius of 2.1 km. Thus, the arc is tangent to the base of layer 6 at the (x , z ) 0 Q coordinate (1.0 km, 2.1 km), and outcrops at the model surface at (3.1 km, 0.0 km). The dip of the six linear segments comprising the arc ranges from approximately 15° at depth to 80° near the surface. This model is very challenging in terms of the large velocity contrasts and range of steep dips present. Furthermore, the model contains a low-velocity zone which often poses troublesome imaging problems, and may be interesting from a geophysical explorationist's perspective. Finally, the model may be similar to the subsurface geology encountered in regions of continental margin sedimentary deposition, or to the salt dome geologic structures common in the Gulf of Mexico, for example. 119 8.2 The Synthetic Data Multifold synthetic seismic reflection trace data are generated from the test model to provide 51 shot gathers for migration. A representative shot gather is shown in Figure 8.2. The synthetic data are raytraced using the algorithm of Zelt and Ellis (1988). Only the primary reflection arrivals are raytraced in this example, and only the acoustic amplitudes and phase rotations due to a 3-D isotropic point source are calculated in order to satisfy the basic assumptions of the Kirchhoff migration theory developed in the preceding sections. Subsequent to calculating the raytraced synthetic data impulse responses, the "spike" traces are convolved with a 1-second zerophase 30 Hz peak-frequency Klauder wavelet to generate the synthetic seismograms. The 51 shot gathers of synthetic data represent a shot spacing of 60 m, and a trace (reciever) spacing of 15 m. The receiver geometry is fixed between 0.0 km and 3.0 km along the surface, while the shot locations vary over the same range. Hence, the recording geometry of the shot gathers ranges from a forward end-on spread at the origin, to a split-spread at 1.5 km, to a reverse end-on spread at 3.1 km. This recording geometry provides 200 traces per shot gather, each trace consisting of 3 seconds of digital data sampled at 4 msec intervals. Thefiveprominent hyperbolic events in the shot gather of Figure 8.2 correspond to the five shallowest horizontal reflectors in the test model. Only a small portion of the sixth horizontal reflector is illuminated by this shot gather recording geometry, and is visible in the synthetic data on trace 40 at 1.75 sec. The linear dipping events in the seismic section correspond to the linear segments that comprise the arc structure in the test model. At this point it is worthwhile noting that the synthetic data do not contain diffraction responses that realistically arise due to the truncation corners of the horizontal layering against the arc structure. Their 120 Distance/15 (m) Figure 8.2: A Representative Shot Gather Raytraced at x = 1200 m along the Arc Test Model. The star symbol marks the shot location. The receiver traces are deployed in an asymmetric split-spread, ranging from trace 1 at 0.000 km to trace 200 at 2.985 km. The record length is 3 seconds. s 121 inclusion is not readily available from the raytracing algorithm. This absence of diffractions has subsequent consequences in the imaging of the truncation corners as will be seen later. Figure 8.3 shows the pre-processed synthetic shot gather ready to be migrated. Since the synthetic data are generated with a bandlimited source wavelet and do not require trace editing, the only pre-processing required is the application of a standard mute pattern. The mute pattern consists of a mask with mute boundaries defined by line segments joining the critical offset-time pairs calculated manually from the test model. A mute ramp of 100 msec was further applied. The mute application is necessary to eliminate post-critical energy, which violates the basic assumptions of the migration method. The synthetic data do not include multiples, critical refractions, surface waves or converted waves, all of which would otherwise need to be removed during the pre-processing step. Finally, it is mandatory that the pre-processing be done in such a way as to preserve the true-relative amplitude information in the data. This means that all trace-dependent scaling processes, such as Automatic Gain Control or deconvolution, should be avoided if true-amplitude recovery is desired in the final migrated depth images. 8.3 The Migration Model Figure 8.4 shows the 1-D migration model obtained by subtracting the 2-D arc structure from the true test model of Figure 8.1. The migration model is otherwise identical to the arc test model in terms of velocity, density, and horizontal layer structure specification. Hence, the 1-D migration model is specified completely by thefirstfour columns of Table 8.1, as shown in Table 8.2. A migration model estimate of the true arc test model is required for subsequent raytracing evaluation the WKBJ Green's function amplitudes and imaging times required to 122 Distance/15 (m) Figure 8.3: The Pre-Processed Shot Gather at x = 1200 m along the Arc Test Model. This figure is identical to Figure 8.2, except for the application of a standard mute pattern to suppress post-critical energy. The pre-processing is required prior to migration to eliminate seismic energy which violates the fundamental assumptions of the Kirchhoff-WKBJ migration method. s 123 Figure 8.4: The 1-D Migration Model Obtained from the 2-D Arc Test Model. The 1-D migration model is obtained by subtracting the 2-D arc structure from the test model of Figure 8.1. Hence, the migration model is identical to the 1-D structural component of the arc test model, in terms of layer positioning and material property specification. 124 Density (glee) Layer Ztase De/tf/j (km) Velocity (km/s) 1 2 3 4 5 6 7 1.500 2.300 3.500 4.500 5.500 2.500 6.000 0.600 0.900 1.200 1.500 1.800 2.100 3.000 1.929 2.147 2.384 2.539 2.670 2.192 2.728 Table 8.2: The 1-D Migration Model Material Property Specification. implement the generalized migration equation (4.1). A 1-D model is an acceptable migration model for the true 2-D nature of the synthetic data, since all reflected energy propagates only through the 1-D portion of the 2-D test model. An algorithm was developed in order to parameterize a 2-D model in general, and a 1-D migration model in particular, on a finite difference grid. The algorithm generates a finite difference migration model to be used as input to the Newton-Raphson two-point raytracing algorithm developed in the previous chapter to calculate the WKBJ Green's function amplitudes and traveltimes. In the 1-D mode of model parameterization, the algorithm accommodates an essentially unlimited number of horizontal layers. Velocities are specified as constants within each layer. Densities may be specified similarly by layer, or defaulted to vary as a function of velocity according to the Gardner relation: p(z) = p0{/c(i), (8.1) which is determined empirically for sedimentary basin rocks (Gardner et al., 1974), and where the constant pG is assigned the value 1.74 for velocity units of km/s. An approximate 125 Q-attenuation can be accommodated in the migration model, based on the assumption that the seismic data attenuate exponentially as a function of one dominant frequency w . This 0 assumption is reasonably valid for the common case of narrowband seismic reflection data. The Q-attenuation thus implemented does not compensate for the true dispersive nature of the data, but rather, attempts to incorporate an amplitude factor during the migration procedure that corrects for the amplitude attenuation of reflected energy due to the estimated subsurface Q structure. For model parameterization, Q values can be specified by layer, or defaulted as a function of velocity using the Waters (1981) relation: Q(z) = Q c {z), 2 0 where the constant QQ may be specified as 18.83 for velocity units of km/s. 126 (8.2) Chapter 9 Raytracing the WKBJ Green's Functions 9.1 Introduction The WKBJ Green's functions play an extremely important role in the development and implementation of the generalized migration solution (4.1). As the kernel functions of the migration integrals, the Green's functions are essentially responsible for mapping the recorded seismograms from observational data space back to the original reflectivity model space. Examination of the integral solution (4.1) demonstrates this in that the Green's functions discriminantly select trace data at the imaging times t' = % + xs, implement amplitude r recovery with the dynamic amplitude weights A and A , and successively map each portion r s of the observed seismic data as a contribution to the total estimate of the subsurface reflectivity function R(x ). As such, it is evident that accuracy in the evaluation of the Green's function 0 imaging times T(X0) and WKBJ amplitudes A(x0) is critical to the eventual quality of the recovered reflectivity estimate i?(x0). An iterative two-point raytracing algorithm was developed to compute the Green's function traveltime and amplitude variations, A(x„) and x(x0), as required to implement the generalized migration equation (4.1). A detailed and extensive treatment of the zeroth order asymptotic ray theory and the iterative Newton-Raphson method of solution developed for this 1.5-D case 127 is presented in Chapters 6 and 7 respectively. The two-point raytracing algorithm calculates traveltimes and amplitudes from a given surface location to all subsurface points specified on afinitedifference grid of the migration model. The algorithm model parameterization admits an unlimited number of horizontal layers, each specified by a depth, P-wave velocity, density, and Q value. Direct transmission times, geometric divergence, specular transmissivity, Qattenuation, and surface-incident obliquity, zeroth order ART amplitudes can be computed. These quantities can be calculated throughout the entire model or in a specific portion of the medium determined by a source aperture cone or by zones of pre-critical reflection. The amplitudes obtained are 1.5-dimensional, since the algorithm assumes a 1-D model illuminated by a 3-D isotropic point source. 9.2 The Raytraced Imaging Times Figure 9.1 shows the raytraced transmission times x(x0) necessary for reflectivity imaging, specified on afinitedifference grid of the complete migration model of Figure 8.4. The color legend specifies traveltime in seconds. The source (dark blue) is located at the origin. The color boundaries represent surfaces of constanttimeand hence delineate direct-transmission wavefronts. It is evident that the wavefronts are circular within thefirstconstant-velocity layer. This corresponds to the well-known fact that spherical wavefronts are obtained from an isotropic point-source in a homogeneous, isotropic medium. At far offsets, it is notable that the wavefront curvature approaches the classical plane-wave approximation. Furthermore, it is observed that as the wavefronts propagate vertically with depth, they are noticably stretched (blue), due to the general velocity increase in that direction. Furthermore, we can note that past the critical distance along the top of a layer boundary, 128 Figure 9.1: A 2-D Subsurface Map of the Raytraced Green's Function Reflectivity Imaging Times x(x0). Each color pixel represents the direct one-way traveltime x(x0; xs) from the source position xs (dark blue), located at the origin, to the pixel subsurface coordinate location x0. The color legend specifies traveltime in seconds. 129 at 3.0 km offset and 0.6 km depth for example, the transmission wavefield splits discontinuosly (red-yellow) across the boundary; incident wave energy cannot refract through the layer boundary at post-critical incidence. This phenomenon results in transmission-time discontinuities across post-critically illuminated sections of layer boundaries. These traveltime discontinuities are obviously a potential source for migration imaging instability. Hence, it is advisable to mute the recorded reflection shot-gather data in order to remove any such post-critical arriving energy. A mute pattern based on the zone of post-critical arrivals can be determined by the migration model velocity distribution. Since the present migration model is 1-D, both the source-subsurface and receiver-subsurface traveltimes, xs and xr, may be obtained from the single traveltime map in Figure 9.1. However, this is not the case for a heterogeneous 2-D model, in which individual xs and xr maps must be obtained separately, from a suite of generic x(xQ) maps raytraced for each unique source and/or receiver surface location. 9.3 The Raytraced WKBJ Amplitudes The raytraced geometric divergence Aj)(x ) due to an isotropic point source excitation in the 0 migration model is displayed in Figure 9.2. The plot displays geometric divergence amplitude attenuation with increasing distance from the source, and is plotted using a logarithmic scale. At the origin, the blue zero dB reading indicates the location of a unit-amplitude source. The amplitude attenuation ranges from zero dB at the source to about 40 dB (red) at far offsets. Such a large dynamic range suggests that compensation for geometric divergence will be a most significant factor in the recovery of true relative-amplitude information. In thefirstlayer the spreading is inversely proportional to the wavefront radius of curvature. 130 Figure 9.2: A 2-D Subsurface Map of the Raytraced Green's Function Geometric Divergence Amplitude Factor Arj(x ). Each color pixel represents the geometric divergence A j ; ( x ; x ) from the unit-amplitude source position x (dark blue), located at the origin, to the pixel subsurface coordinate location X Q . The color legend specifies amplitude attenuation on a logarithmic (dB) scale. 0 0 s 1 3 1 s This is the familiar 1/r geometric-spreading term in a constant-velocity medium. For subsequent layers the geometric divergence changes in a complicated fashion. At post-critically illuminated sections along a boundary, the transmission amplitude is zero (dark red) since the wave cannot transmit beyond critical incidence. To obtain information beneath postcritical sections of layer boundaries, the source energy must transmit pre-critically at small offsets, and then leak out horizontally underneath the boundary, as exemplified by the blue in layer 2 at 0.75 km depth. This results in large wavefront radii of curvature tending toward vanishingly small red amplitudes beneath post-critical sections of layer boundaries. Figure 9.3 displays the cumulative raypath specular transmissivity, plotted on a linear scale, and represents the transmission amplitude component of the general Zoeppritz energy partitioning amplitude factor Ap(x0), which arises at layer boundaries. This amplitude factor is calculated assuming an acoustic plane-wave specular-transmission coefficient T; at each subsurface location. Successive coefficients are multiplied together along a given raypath resulting in a cumulative transmissivity estimate from the source location to any subsurface location. Note that the transmissivity is constant along an individual ray within a given layer, and changes across layer boundaries as a function of incident angle and impedance contrast. Within the shallowest layer, the source energy is perfectly transmitted, Ap = 1., which plots as a constant hue of red. Successively less energy is transmitted as a function of raypath depth and offset, as suggested by the dark blue plot regions. As an interesting remark, the light blue transmission energy is relatively large in the low-velocity layer, since there is no post-critical-incidence barrier impeding transmitted energy incident from shallower regions of the model. Figure 9.4 demonstrates the variation of the cos 6r obliquity factor of the generalized migration equation (4.1), as plotted on a linear scale. Note that 0r is the incident angle, or 132 Figure 9.3: A 2-D Subsurface Map of the Raytraced Green's Function Energy Partitioning Amplitude Factor Ap(Xo). Each color pixel displays the cumulative specular transmission amplitude factor Ap(x0;xs) along a raypath from the unit-amplitude source position xs, located at the origin (dark red), to the pixel subsurface coordinate location x0. The color legend specifies cumulative amplitude on a linear scale. 133 Figure 9.4: A 2-D Subsurface Map of the Raytraced Surface-Incident Obliquity Amplitude Factor cos8 (x ). Each color pixel displays the obliquity amplitude factor cos 6 ( x ; x ) , detenruned by the surface-incident angle of the direct transmission ray joining the source position x , located at the origin (dark red), to the pixel subsurface coordinate location Xo. The color legend specifies obliquity amplitude on a linear scale. r 0 r s 134 0 s reciprocally, the take-off angle, with respect to the vertical axis, of a ray at the surface. This term is not of leading order significance in situations of reasonable increasing-withdepth velocity structure, since in such cases it is evident that its amplitude is nearly unity (red) over most of the model. To first order, only arrivals originating at large offsets from within the shallow near-surface layers experience any significant obliquity variation. If amplitude recovery in this zone is not critical, the obliquity term can be safely ignored without deterioration in the amplitude estimation of deeper reflectivity structure. Combining the amplitude variations, Figure 9.5 represents the total WKBJ Green's function amplitude factor A(x0), obtained by multiplying the geometric divergence and specular transmissivity factors Ajr> and Ap at each subsurface location x0. The plot displays amplitude attenuation, calculated on a logarithmic (dB) scale. The dominant contribution of the geometric divergence Ap, Figure 9.2, to the total WKBJ amplitude factor A is immediately obvious. Since the migration model is 1-D, the source and receiver WKBJ amplitude components A andAr are both obtainable, as a function of offset, from the single general WKBJ s amplitude map A(x0) shown in Figure 9.5. However, this is not the case for a heterogeneous 2-D migration model, in which the individual A and A maps must be obtained separately, s r from a suite of generic A(x0) maps raytraced for each unique source and/or receiver surface location. Similarly, Figure 9.6 shows the WKBJ Green's function amplitude factor A(x0) of Figure 9.5 modified by the surface-incidence obliquity factor cos 8r. Figures 9.5 and 9.6 are identical except for the obliquity factor cos 0r which is seen to only significantly affect the first-layer response of the latter, as predicted. In summary, Figure 9.1 displays the Green's function reflectivity imaging times Ts(x0) and tr(x0) required for implementing the generalized migration equation (4.1). Figures 9.5 and 135 Figure 9.5: A 2-D Subsurface Map of the Raytraced W K B J Green's Function Amplitude Factor A ( x ) . Each color pixel displays the total W K B J Green's function amplitude factor Ai(x ; Xj), determined by multiplying the 2-D maps of geometric divergence AD (Figure 9.2), and specular transmissivity Ap (Figure 9.3), at each subsurface location Xo, with respect to the "source" location Xj at the origin. The color legend specifies amplitude attenuation on a logarithmic (dB) scale. 0 0 136 Figure 9.6: A 2-D Subsurface Map of the Raytraced Obliquity-Modified W K B J Green's Function Amplitude Factor A(x ) cos 8 (x ). Each color pixel displays the total W K B J Green's function amplitude factor A;(Xo; Xj) of Figure 9.5, modified by the surface-incident obliquity amplitude factor cos 0 (x ; x ) of Figure 9.4. The color legend specifies amplitude attenuation on a logarithmic (dB) scale. 0 r r o 0 r 137 9.6 display the necessary WKBJ Green's function amplitude terms As(x0) and Ar(x0), and the surface-incidence obliquity amplitude variation cos0r(xo). Since the migration model is 1-D, the source-to-depth-point and receiver-to-depth-point traveltimes and amplitudes are merely a function of source-receiver offset. This greatly reduces the migration raytracing computational effort compared to the case of a 2-D migration model. For an arbitrary 2-D migration model, raytraced traveltime and amplitude maps are required to be computed for each source and receiver location, or equivalently, for every unique source and/or receiver surface position, assuming shot and receiver reciprocity. 9.4 Computational Requirements The Newton-Raphson two-point raytracing algorithm affords both rapid and accurate numerical evaluation of the WKBJ Green's functions over a 2-Dfinitedifference grid of the subsurface migration model coordinate xQ. Several computational examples were performed to construct the CPU time listing of Table 9.1. Each example consisted of a single (N ,N ) X Z gridded map of raytraced values, where N = N = 101. Hence, in effect, a ray is traced x z from a single shot location, N = 1, to each of N • N = 10,201 subsurface "receiver" s x z locations. In general, the computation time required for the raytracing component of the migration is proportional to N • N • N , where N is the number of raytraced maps required m x z m to implement the migration. All raytracing computations were performed on a MicroVax II with 9 Mb RAM. The respective CPU times for each raytraced map item are given in Table 9.1, and do not include I/O processing. These CPU times can be expected to decrease by at least one order of magnitude with the use of an array processor, and at least two orders of magnitude with the use of a supercomputer. 138 CPU Time (s) Raytraced Quantity Transmission Times Geometric Divergence A Energy Partitioning A Q Attenuation Surface Obliquity cos 0r All Amplitudes At All Amplitudes and Times Ai, Xi D P 19.3 16.6 24.2 24.5 15.7 25.2 29.1 Table 9.1: CPU Times for Green's Function Raytracing Evaluation. Since the migration model is 1-D, the Green's functions depend only on the offset and depth coordinates, and not independently on the source and receiver coordinates. For this reason, only one raytraced map, N = 1, is required to implement the generalized Kirchhoffm WKBJ migration equation (4.1). Hence, the CPU times listed in Table 9.1 represent the total computational requirement for the raytracing component of the Kirchhoff prestack depth migration examples discussed later. However, for the general 2-D case, the minimum required raytracing effort will be on the order of the quoted times for a single map, multiplied by the number of surface locations 7VS considered. Since the shot spacing is usually greater than the receiver-group spacing, the number of surface positions is roughly equivalent to the number of independent receiver positions along a survey line. Finally, it may be more efficient, for large survey data sets, to raytrace a sparse set of traveltime and amplitude maps, and then interpolate between these maps to obtain the required intermediate subsurface location Green's function evaluations. 139 Chapter 10 A Quantitative Analysis of Seismic Reflection Data Migration Examples 10.1 Introduction This chapter is concerned with strengthening the physical understanding of the generalized Kirchhoff-WKBJ prestack depth migration equation (4.1), and the quantitative analysis of the the five reflectivity imaging conditions in terms of image resolution and amplitude accuracy. A particular emphasis is placed on the quantitative accuracy of the migration reflectivity amplitudes, since this is the central theme of the thesis. In light of this, numerous carefully constructed synthetic data examples are considered. Some analysis has also been performed with "real" field data, resulting in conclusions similar to that of the synthetic data analysis. However, field data examples will not be presented here since the uncertainty regarding the true nature of the subsurface makes it difficult to realistically quantify those results. At this point, we have discussed the philosophical and theoretical development leading to a true-amplitude migration formalism, and the resulting generalized Kirchhoff-WKBJ migration equation (4.1). A blend of mathematical and physical analysis has been applied to give insight to the first-order behaviour of the reflectivity amplitude recovery inherent to each of thefiveimaging conditions. An examination of zeroth order asymptotic ray theory, and the development and implementation of an iterative 1-D two-point raytracing solution to 140 numerically evaluate the migration Green's functions, exposes a fundamental link between the Kirchhoff-WKBJ migration theory and the solutions to the eikonal and transport equations of seismic ray theory. The generalized Kirchhoff-WKBJ migration can now be implemented since the raytraced Green's function traveltime and amplitude maps have been evaluated. We will presently examine various migration examples resulting from the implementation of Equation (4.1). 10.2 Single Trace Depth Migration Responses To begin the migration analysis, we consider the single trace migration responses of the synthetic arc test model data, Figure 8.3, to the five reflectivity imaging conditions. The single trace migrations are equivalent to computing the integral kernel function of the migration equation (4.1) for a single (xs, xr) trace location, and represent a simple starting point for understanding the complete prescription. In this way, we will gradually build up a physical intuition for the complete Kirchhoff-WKBJ prestack depth migration process,firstby examining the single trace migrations, then by summing reflectivity contributions from the receiver locations xr to obtain a single migrated shot gather at xs, and thenfinallyby summing over all shot gathers xs, as prescribed by (4.1), to obtain the full-fold migrated depth sections. 10.2.1 Dynamic Migration Figure 10.1 shows the single trace migration response to the dynamic reflectivity imaging condition, obtained by migrating a far-offset trace from the muted shot gather of Figure 8.3. The shot location is marked by the star symbol at 1200 m, and the receiver location is marked by the inverted triangle at 0 m. This shot and reciever location corresponds to the left-most 141 Dynamic Distance/15 (m) Figure 10.1: Single Trace Depth Migration Response to the Dynamic Reflectivity Imaging Condition. The single trace data correspond to a shot position at xs = 1200 m, marked by the star symbol, and a receiver position at xr = 0 m, marked by the triangle symbol. trace of the muted shot gather displayed in Figure 8.3. The migrated result is displayed in distance and depth coordinates, as opposed to CMP location and vertical traveltime, and has the same dimensions as the arc test model of Figure 8.1. The migrated depth section displays all possible subsurface reflecting locations that have the potential to reproduce the reflection events observed on the recorded seismic trace for the given shot-receiver geometry and migration velocity model. Note that only four migration ellipses, corresponding to four out of a possible six horizontal reflection trace events, are in evidence. This is because the two shallowest reflection events, at this large offset, have been muted from the post-critical zone of the shot gather data, as shown in Figure 8.3. Also, no portion of the arc is present in the migrated result, since the arc reflections arrive too late to be recorded at this large offset. Coherent energy can be seen to be aligned at the true horizontal reflector depths of 1200 m, 1500 m, 1800 m, and 2100 m. The coherent energy focusses midway between the shot and receiver location, as expected, since this coincides with the reflection depth point for this trace in a 1-D medium. In a constant velocity medium, for a ZSR trace, a migrated event would sweep out a circle in the depth domain. Due to the non-zero offset of this trace, the migrated events sweep out portions of ellipses whose foci coincide with the source and receiver locations. Due to the non-constant migration velocity, the ellipses are stretched downward in depth. At subsurface locations corresponding to super-critical reflection, the migration ellipses are forced to be discontinuous across migration velocity model boundaries, in compliance with the physics of wave propagation through the given migration velocity structure. The migration ellipses demonstrate a significant amplitude asymmetry, as predicted by the migration equation (4.1) in the case of the dynamic imaging condition. The dominant dynamic 143 amplitude factor, \jA~i.lA , which is approximated to first order by the geometric amplitude s factor, r /y/ry, is relatively small near the source location (r < ry), and relatively large near s s the receiver location (r > r ). Relatively moderate migration ellipse amplitudes occur near s r the central midpoint offset region where r ~ r>. s The migration ellipses demonstrate the good steep-dip resolution of the algorithm where they sweep up to near-vertical dip and maintain good signal-to-noise quality. The effect of the cos 9 obliquity factor is evidenced by the attenuation of the migration ellipses near the r model surface. One can also observe that the phase of the single trace migration operator lags the zero-phase of the input trace wavelets by 45°, which arises from the effect of the half-time-derivative operator in (4.1). If we migrate over many trace offsets within a shot gather, simulating a line integral over many point-receiver responses, we incur an opposite 45° phase-shift. The net result is a shot-gather migration which has the desired zero-phase property. 10.2.2 Geometric Migration The single trace migration response to the geometric reflectivity imaging condition is displayed in Figure 10.2. The geometric migration ellipse locations coincide with those of the dynamic imaging condition, as indeed all of the imaging conditions do, since the phase, or equivalently, the imaging times Tj, are identical among the various imaging conditions. The amplitude variation along the migration ellipses, however, differs considerably among the five imaging conditions. The amplitude variation along the geometric migration ellipses is seen to have the same asymmetry, although not as pronounced, as the dynamic migration response of Figure 10.1. 144 Geometric TSN 20 40 60 80 100 1.20 140 160 180200220240260 TSN Distance/15 (m) Figure 10.2: Single Trace Depth Migration Response to the Geometric Reflectivity Imaging Condition. The single trace data correspond to a shot position at x = 1200 m , marked by the star symbol, and a receiver position at x = 0 m, marked by the triangle symbol. s r The geometric amplitude factor r /^/ry has been constrained to vary only within a range s determined by estimating the minimum and maximum one-way reflection raypath distances r as 600 m and 3000 m respectively. This constraint decreases the allowable dynamic range of the geometric migration amplitudes in order to achieve stability in the amplitude recovery. Other than the stabilization effect, the geometric amplitude variation is a good approximation to the correct dynamic migration response. 10.2.3 Kinematic Migration Figure 10.3 displays the single trace migration response to the kinematic imaging condition. Recall that this imaging condition is a highly regularized version of the dynamic and geometric imaging conditions, such that the true migration amplitude factor y/A~^/A is constrained s to be constant, i.e., no source-receiver amplitude variation is represented at all. As a result, the amplitude variation along each migration ellipse is symmetric and fairly constant. Some residual amplitude asymmetry is present near the model surface, due to the cos 0 obliquity r factor. One may observe that relative amplitude variation of the kinematic migration ellipses is similar to that of the dynamic and geometric migration responses, near the true reflector locations at 600 m offset. This reinforces previous analysis which indicates that the kinematic imaging condition recovers accurate amplitudes in 1-D reflectivity structure. 10.2.4 Excitation-Time Migration Figure 10.4 displays the single trace migration response to the excitation-time imaging condition. The excitation-time imaging condition exhibits a similar amplitude asymmetry to that of the dynamic and geometric imaging conditions. However, the kinematic variation 146 Kinematic TSN 20 40 BO 8 0 1 0 0 ) 2 0 140 160 180 2 0 0 2 2 0 2 4 0 2 6 0 TSN Distance/15 (m) Figure 10.3: Single Trace Depth Migration Response to the Kinematic Reflectivity Imaging Condition. The single trace data correspond to a shot position at x = 1200 m, marked by the star symbol, and a receiver position at x = 0 m, marked by the triangle symbol. s r Excitation-Time Distance/15 (m) Figure 10.4: Single Trace Depth Migration Response to the Excitation-Time Reflectivity Imaging Condition. The single trace data correspond to a shot position at x = 1200 m, marked by the star symbol, and a receiver position at x = 0 m, marked by the triangle symbol. s r is erroneously manifested by a \fA factor, as opposed to the correct \fA lA r r s factor. The uncompensated source divergence A remains buried in the excitation-time reflectivity ims age, and appears in Figure 10.4 as a gradual reflectivity amplitude loss with shot-reflector distance in the migration ellipses. As a result, the excitation-time migration ellipses appear to get progressively weaker with depth as compared to the dynamic or geometric migration response. 10.2.5 Crosscorrelation Migration Figure 10.5 displays the single trace migration response to the crosscorrelation imaging condition. The crosscorrelation imaging condition exhibits a dramatically different amplitude behaviour compared with the other imaging conditions. The kinematic amplitude variation is controlled by a A \/A^ factor, which is in error by a factor of Aj compared with the dynamic s or geometric imaging conditions. Hence, the reflectivity amplitude recovery is attenuated by the uncompensated source divergence, such that the amplitude loss is proportional to the square of the shot-reflector distance. This appears in Figure 10.5 as an over-accentuation of near-surface reflectivity amplitudes, and a drastic attenuation of deeper reflectivity amplitudes. In this way, the crosscorrelation imaging condition represents an accelerated amplitude-loss version of the excitation-time migration. 10.3 Single Shot Gather Depth Migration To continue the migration analysis, we follow the migration equation (4.1), which indicates that we must integrate (sum) the single trace migration responses for all receiver locations xr within a shot gather. This provides a single shot gather migration, and results in a reflectivity 149 Crosscorrelation 20 TSN OoO 500a0 4 0 6 0 8 0 1007,20 1 4 0 1 6 0 1 8 0 200 220 2 4 0 2 6 0 I I " i 1 * i I I I I I | L, TSN "TT OaO 500a0 -.T liliB 1000 "5 o lOOOoO aO"^ •ir 1 5 0 0 = 0^- 1500a0 a O 2 0 0 0 „ 0 2 0 0 0 a 0 + 2500.(Hi 2500o0 3 0 0 0 . , 0- 3000u0 Distance/15 (m) Figure 10.5: Single Trace Depth Migration Response to the Crosscorrelation Reflectivity Imaging Condition. The single trace data correspond to a shot position at x = 1200 m, marked by the star symbol, and a receiver position at x = 0 m, marked by the triangle symbol. s r model estimate for that portion of the subsurface which is illuminated by the particular shot gather location and geometry in question. Since the reflectivity estimate consists of single trace migrations from separate, and unique, offsets within a single shot gather, information regarding the variation of i?(x0,8) as a function of reflection angle may be recoverable. True-amplitude shot gather migration may have important applications for amplitude-offset (AVO) analysis and inversion, for example. 10.3.1 Dynamic Migration Figure 10.6 shows a single shot-gather Kirchhoff-WKBJ prestack depth migration obtained by migrating all 200 offset traces in the gather of Figure 8.3, corresponding to a surface shot location at 1200 m, using the dynamic reflectivity imaging condition. A considerable portion of the horizontal and arc reflectivity structure has been imaged well, although some residual migration ellipse and wavefront discontinuity artifacts remain uncancelled at the edges of the subsurface shot coverage. Furthermore, both the horizontal and arc reflectivity reflectivity structure are positioned correctly in the depth image, compared with the true arc test model of Figure 8.1. At this shot location, most portions of the arc are weakly illuminated by reflected energy. Since the arc structure is very distant from the shot location, the source amplitude A attenus ates severely before reflection. For this shot location, A has attenuated to an amplitude level s 2 on the order of or less than, the dynamic stabilization parameter e = 0.005, as can be seen in the WKBJ amplitude ray map of Figure 9.5. In this case, the stabilized dynamic imaging condition amplitude variation y/A^/(A +£ ) 2 s reduces to constant • \fA~i-, which is proportional to the excitation-time amplitude recovery. Hence, the severe source amplitude attenuation and subsequent stabilization results in the arc being only weakly imaged, as opposed to the 151 Dynamic TSN 20 OnO -m I 40 1 GO 80 1 1 1 0 0 : 2 0 140 160 180 200 220 240 260 I * l I J 1 1 1 L 1 TSN ^0 = 0 500a0 5 0 0 » 0 -.1 ~'r lOOOaO lOOOaO- 1 500 • 0 L500n0^r 2000.(HE- :r 2 0 0 0 a 0 2500-O^r -.r llllllllllllWIMIlllllllllllllllI 3000aO- 2500„0 3000o0 Distance/15 (m) Figure 10.6: A Single Shot Gather Depth Migration using the Dynamic Reflectivity Imaging Condition. The single shot gather data correspond to a shot position at x = 1200 m, marked by the star symbol, and 200 receiver positions, ranging from x = 0 m (migrated depth trace 34) to x = 2985 m (migrated depth trace 233), in 15 m increments. s r r anticipated true-amplitude reflectivity recovery expected of the dynamic imaging condition. Later, it will be shown that this overstabilization phenomenon has somewhat compromised the amplitude recovery on the deep horizontal reflectors as well. Figure 10.7 shows a ray diagram indicating the subsurface illumination pattern of primary reflected energy which arrives within the given surface receiver array at this shot location. The close correspondence between the ray diagram illumination of reflectors in Figure 10.7, and the recovered migration reflectivity estimate of Figure 10.6, indicates that the kinematic physics of the Kirchhoff-WKBJ prestack migration equation are correct. The correspondence between the ray diagram illumination and reflectivity estimate resolution is particularly remarkable along the portion of the arc in the low-velocity layer at 2000 m depth, for example. Also, only a small portion of the sixth reflector is imaged since this particular shot gather geometry does not illuminate much of that structure, as shown in the ray diagram of Figure 10.7, and by the notable lack of associated reflected arrivals in the synthetic shot gather data of Figure 8.3. 10.3.2 Geometric Migration Figure 10.8 displays the Kirchhoff-WKB J prestack depth migration of the shot gather in Figure 8.3 using the geometric reflectivity imaging condition. Again, the horizontal and arc reflectivity structure are correctly positioned in depth, as governed by the accuracy of the imaging times T; and the migration velocity model. The geometric reflectivity amplitudes seem to increase with depth, which is the correct trend of the true reflection coefficients derived from the arc test model of Figure 8.1. This is in direct contrast to the general decrease in reflectivity amplitude with depth of the dynamic migration result, as shown in Figure 10.6. Furthermore, the geometric reflectivity estimates along the arc structure are 153 DISTANCE CJ-0.5 0.0 0.5 1-0^ 1-5 (Krl) 2.0 2.5 3.0 3.5 63 Figure 10.7: Ray Diagram Hlumination Pattern of Reflected Arrivals for a Shot Gather at x = 1200 m. Only those raypaths are displayed corresponding to primary reflections which arrive back at the recording surface within the model boundaries. s Geometric TSN OoO - S Cu 20 I 40 60 l I 80 -I IOC !?0 140 160 180 200 220 240 260 I * i 1 500-0 -.'r 500o0 1000-O^r ir [500 = 0-^ ~.r 1500-0 2000-0-$ •r 2000o0 2500o(Hr ir 2500-0 1000-0 Q 3000-0 3000-0- Distance/15 (m) Figure 10.8: A Single Shot Gather Depth Migration using the Geometric Reflectivity Imaging Condition. The single shot gather data correspond to a shot position at x = 1200 m , marked by the star symbol, and 200 receiver positions, ranging from x = 0 m (migrated depth trace 34) to x = 2985 m s r (migrated depth trace 233), in 15 m increments. r relatively larger in amplitude than the dynamic case, which is again of significance compared to the correct model trend. In fact, it will be shown in the full-fold migration section that the geometric imaging condition recovers a superior reflectivity amplitude estimate compared with the stabilized dynamic imaging condition. Evidently, the geometric amplitude factor r /^/ry contains sufficient infors mation regarding the true dynamic amplitude variation T/A^/A , S and is adequately stable, in order to recover superior reflectivity amplitude estimation when compared with the formally correct, but implementationally somewhat unstable, dynamic imaging condition. 10.3.3 Kinematic Migration Figure 10.9 displays the Kirchhoff-WKBJ prestack depth migration of the shot gather in Figure 8.3 using the kinematic reflectivity imaging condition. The kinematic migration reflectivity amplitude estimation displays a highly stable blend of the dynamic and geometric amplitude recovery properties. The 1-D horizontal reflectivity structure appears to have the correct reflectivity amplitude variation, similar to the geometric result of Figure 10.8. This is as expected, since the physical analysis of Chapter 5 suggests that the kinematic imaging condition provides accurate and stable amplitude recovery in regions of 1-D reflectivity structure. The reflectivity amplitudes along the 2-D arc structure, however, have a balanced appearance compared to the geometric or dynamic results. As a result, although the true reflectivity amplitudes decrease along the upper five segments of the arc, as accurately portrayed in the dynamic and geometric results, the kinematic migration amplitudes appear to be more constant along the arc. This is the amplitude balancing effect predicted by the analysis of 156 Kinematic 20 TSN 40 GO 80 IOC ! 20 140 160 180 200 220 240 2 6 0 OaO - 5 0 0 a 0 -_~r TSN 3 - OaO V lOOOaO lOOOaO— 53 a 500a0 l500 Oj- -ir 1500 = 0 2000a(Ht -.r 2 0 0 0 a 0 2500a(Hr ~~ 2 5 0 0 a 0 o Q 3000a0 3000a0" Distance/15 (m) Figure 10.9: A Single Shot Gather Depth Migration using the Kinematic Reflectivity Imaging Condition. The single shot gather data correspond to a shot position at x = 1200 m, marked by the star symbol, and 200 receiver positions, ranging from x = 0 m (migrated depth trace 34) to x = 2985 m (migrated depth trace 233), in 15 m increments. s r r Chapter 5. Although the true reflectivity amplitude variation is proportional to the factor y/A^/A , the kinematic imaging condition forces this ratio to be a constant. Hence, although s the true-amplitude reflectivity recovery factor has a dynamic range of much greater than, to much less than unity, the kinematic factor always takes the same intermediate conservative value. This results in an amplitude-balanced image in regions of 2-D velocity variation, as opposed to the true-amplitude 2-D reflectivity estimation of the dynamic and geometric imaging conditions. This effect becomes more pronounced in the full-fold migration examples. 10.3.4 Excitation-Time Migration Figure 10.10 displays the Kirchhoff-WKBJ prestack depth migration of the shot gather in Figure 8.3 using the excitation-time reflectivity imaging condition. Figure 10.10 clearly displays the loss of reflectivity amplitude recovery inherent in the excitation-time imaging condition. It is evident that the excitation-time imaging condition could be useful as a reflectivity imaging process, since an application of AGC to Figure 10.10 would definitely highlight the correct reflectivity structure. However, the reflectivity amplitude estimation, necessary for the complete reflectivity model construction, is definitely unreliable. 10.3.5 Crosscorrelation Migration Figure 10.11 displays the Kirchhoff-WKBJ prestack depth migration of the shot gather in Figure 8.3 using the crosscorrelation reflectivity imaging condition. Figure 10.11 dramatically displays the loss of reflectivity amplitude recovery inherent in the crosscorrelation imaging condition. Recalling the analysis of Chapter 5, the amplitude error (loss) is proportional to 158 Excitation-Time TSN 0.0 T 20 I 40 I 60 I 80 100 '20 140 160 180 200 220 240 260 I I* J 1 1 1 1 1 1 1 r TSN OoO 500o0 :r 500-0 lOOOoO-^r :r lOOOaO 1500=0 1500a0" -.'r 2000o0 2000., (Hr 2500-CHr 2500o0 3000aO- 3000a0 Distance/15 (m) Figure 10.10: A Single Shot Gather Depth Migration using the Excitation-Time Reflectivity Imaging Condition. The single shot gather data correspond to a shot position at xs = 1200 m, marked by the star symbol, and 200 receiver positions, ranging from xr = 0 m (migrated depth trace 34) to xr = 2985 m (migrated depth trace 233), in 15 m increments. Crosscorrelation TSN OaO - TSN 4 0 6 0 8 0 100 ; 2 0 140 160 180 2 0 0 220 2 4 0 2 6 0 i | | I AI l l l I I I I OaO 1 L 1 : — —; I ——T-—' ;'"; ' ^ — ™ T T ^ — — : • •: r; •;•!. , ••!..! • |- ] ; ; 20 11 l :r 500a0 500a0 -|- 1000a0-=r lOOOaO 1500-0$- ~.~ I 5 0 0 a 0 2000a0 2000a0- -.r 2 5 0 0 a 0 2500a(Ht 3000aO" i . i 3000a0 ; Distance/15 (m) Figure 10.11: A Single Shot Gather Depth Migration using the Crosscorrelation Reflectivity Imaging Condition. The single shot gather data correspond to a shot position at x = 1200 m, marked by the star symbol, and 200 receiver positions, ranging from x = 0 m (migrated depth trace 34) to x = 2985 m (migrated depth trace 233), in 15 m increments. s r r the shot-reflector distance, r. The implied amplitude attenuation with depth is readily observable. Only the two shallowest reflectors are imaged at all, and the arc is notably absent. This result illustrates the unfavourable image resolution and amplitude estimation properties of the crosscorrelation reflectivity imaging condition. 10.4 Multifold Prestack Depth Migration At this point, we have examined the single trace and the single shot gather KirchhoffWKB J prestack depth migrations obtained by partially implementing Equation (4.1) as a function of the five reflectivity imaging conditions. Now we examine the complete full-fold implementation of (4.1) for each imaging condition by integrating (summing) the reflectivity contributions from all available shot gather traces. A total of 51 shot gathers, 200 traces per gather are migrated into each full-fold reflectivity depth image. The shot spacing of 60 m gives approximately 3 km length of surface coverage over the arc test model. Further details regarding the complete synthetic survey data are discussed in Chapter 8. 10.4.1 Dynamic Migration Figure 10.12 displays the full-fold Kirchhoff-WKB J prestack depth migration for the dynamic reflectivity imaging condition, obtained by migrating all 51 synthetic shot gathers individually and stacking the results. The positioning of the 1-D horizontal reflectivity structure and the 2-D arc structure accurately images the reflective boundaries of the true arc test model, Figure 8.1, used to generate the synthetic shot gather data. Some residual migration-ellipse artifacts distort the image, 161 Dynamic TSN 20 40 60 80 100 120 140 160 180 200 220 240 260 TSN Distance/15 (m) Figure 10.12: Full-Fold Kirchhoff-WKBJ Prestack Depth Migration using the Dynamic Reflectivity Imaging Condition. A total of 51 shot gathers are migrated for shot locations ranging from xs = 0 m (migrated depth trace 34) to xs = 3000 m (migrated depth trace 234). The 200-receiver array is fixed independent of shot location, and spans nearly the same range as the shot locations: xr = 0 m to xr = 2985 m. especially near the truncation corners at the abuttment of the horizontal layering against the arc. These artifacts result from the absence of diffraction energy in the synthetic seismograms, and thefinitedomain of the total survey subsurface illumination pattern. However, although some artifact distortion blurs the image, the truncation corners are still clearly located and imaged. Such image clarity near truncation corners may be important, for example, in the structural delineation of hydrocarbon reservoirs and traps, especially in low-velocity zone settings, like that of layer 6 at 2000 m depth. Some spatial aliasing is evident in the image of the two shallowest arc segments. The extreme dip of these two reflectors, approximately 70° and 80° respectively, is too steep to be correctiy imaged with the current shot spacing of 60 m. A simple calculation, using a maximum frequency content of 50 Hz and the near-surface velocity of 1500 m/s, shows that the minimum anti-aliasing spatial sampling increment needs to be about 15 m or less. Although the migration depth model and receiver sample intervals, dx = 15m, dx = 15m, r 0 and dz = 7.5m, meet the anti-aliasing requirement, the horizontal shot spacing, dx = 60, 0 s does not. This implies that individual shot gather migrations will not be spatially aliased, as we have seen, but that steeply dipping structure in the full-fold stacked image may be horizontally distorted. As the individually migrated shot gathers are stacked together, the steeply dipping migration structures are not spacedfinelyenough to achieve constructive interference. Rather, the large horizontal shot gather spacing shifts the near-vertical structural images out of phase, causing the destructive interference observed in Figure 10.12. Further deterioration in the image of the near-surface reflectors arises from the necessary application of a post-critical mute to the shot gathers. Reflected arrivals from shallow structure often fall within the post-critical mute zone in the shot gather data, thus reducing the 163 Layer 1 2 3 4 5 6 7 1.500 2.300 3.500 4.500 5.500 2.500 6.000 P(z) 1.929 2.147 2.384 2.539 2.670 2.192 2.728 Rabs Rnorm 0.261 0.256 0.156 0.125 -0.456 0.498 1.000 0.982 0.597 0.478 -1.748 1.909 - - Table 10.1: The True Absolute and Normalized 1-D Horizontal-Layer Normal-Incidence Acoustic Reflection Coefficients as Calculated from the Arc Test Model. ability to accurately image and amplitude-estimate near-surface, or steeply dipping, reflectivity structure. This type of mute pattern distortion has significantly affected the image and amplitude estimate of the shallowest reflector in Figure 10.12. Table 10.1 lists the true absolute and normalized normal-incidence acoustic reflection coefficients calculated manually from the horizontal reflectors in the arc test model. Table 10.2 compares the true normalized reflection coefficients with the corresponding peak-amplitude reflectivity estimates obtained directly from the dynamic migration depth section. The depths of the migrated reflectivity peaks are accurate to within half of one depth sample increment, or + 3.5 m. The coefficients have been normalized relative to thefirstpeak amplitude at 600 m depth. The dynamic reflection coefficients are obtained by performing a partial stack of 11 depth traces along a small portion of each horizontal reflector, in order to reduce the possibility of spurious amplitudes. The 11-trace partial stack represents a 150 m average along each horizontal reflector, from x r = 685 m to x r = 835 m. The stacked depth trace is then searched for peak values within a specified set of depth ranges. A peak reflection coefficient amplitude is obtained from within a depth window of 200 m centered about each true reflector depth. 164 Layer True 1 2 3 4 5 6 Dynamic 1.00 1.00 ± .31 0.98 1.06 ± .16 0.60 0.56 + .06 0.48 0.39 ± .03 -1.75 -0.87 ± .08 1.91 0.78 ±.11 Table 10.2: Normalized 1-D Horizontal-Layer Reflection Coefficients from the True Arc Test Model and the Dynamic Migration Depth Section. The peak amplitude error is calculated by comparing the two depth-trace samples adjacent to the peak amplitude value sample, and then by halving the minimum peak-neighbour amplitude difference. Hence, the possibility of peak amplitude measurement error due to traveltime, amplitude and sampling inaccuracies, is greater for more discontinuous reflectivity averaging functions like the one at 600 m depth (± 31%), compared to smoother averages like at 1500 m depth (± 8%). Although this error measurement is somewhat unconventional, it suits the application presented here. However, the extremal error values computed by this method are probably somewhat over-pessimistic and over-conservative depending on which extremal bound is considered. Table 10.2 illustrates the reflectivity amplitude recovery behaviour of the dynamic imaging condition. The three shallowest reflection coefficients are all accurately recovered to within 2 the allowed error tolerances. At 1500 m depth, the stabilization parameter e = 0.005 begins to have an effect on the amplitude recovery of the fourth reflection coefficient, which narrowly misses the allowable error tolerance. The stabilization parameter value 0.005 corresponds approximately to the 25 dB contour in the A raytraced amplitude map of Figure 9.5. s The 25 dB contour is encountered near the 1500 m depth, and thus amplitudes below this depth become increasingly erroneous from an application of the stabilized dynamic imaging 165 Arc c(z) 1 2 3 4 5 6 7 1.500 2.300 3.500 4.500 5.500 2.500 6.000 PC*) Rabs Rnorm 1.929 2.147 2.384 2.539 2.670 2.192 2.728 0.700 0.536 0.325 0.178 0.054 0.498 2.680 2.055 1.244 0.681 0.208 1.909 - - Table 10.3: The True Absolute and Normalized 2-D Dipping Arc Segment Normal-Incidence Acoustic Reflection Coefficients as Calculated from the Arc Test Model. condition. The reflectivity amplitude error increase with depth is borne out by the remaining dynamic reflectivity estimates tabulated in Table 10.2. Table 10.3 lists the true absolute and normalized normal-incidence acoustic reflection coefficients calculated manually from the dipping arc segments in the arc test model. Table 10.4 compares the true normalized dipping arc segment reflection coefficients with the corresponding peak-amplitude reflectivity estimates measured direcdy from the dynamic migration depth section. The coefficients have been normalized relative to the same peak amplitude as Table 10.2, corresponding to the shallowest horizontal reflector at 600 m depth. The dynamic migration arc segment reflection coefficients are obtained by direct measurement (pencil and ruler) from the depth section. This method is convenient, since dipping events cannot be easily partially stacked, and the amplitude variation along each arc segment makes it difficult to obtain a reasonable amplitude estimate by computer event detection, as was done for the horizontal reflectors. The arc segment peak amplitude error is assigned an estimate of ± 10%, which accounts for reasonable measurement error and amplitude variability within a given arc segment. In a qualitative sense, the migration-amplitude variation along the arc structure corresponds 166 Arc True Dynamic 1 2 3 4 5 6 2.68 2.06 1.24 0.68 0.21 1.91 1.82 ± .18 1.21 + .12 0.65 ± .07 0.39 ± .04 0.11 + .01 0.78 + .08 Table 10.4: Normalized 2-D Dipping Arc Segment Reflection Coefficients from the True Arc Test Model and the Dynamic Migration Depth Section. well with the true relative-amplitude normal-incidence reflection-coefficient variation, including thefirstlayer where the shot-gather mute pattern has largely dip-filtered this obliquelyreflected energy. However, in a quantitative sense, the dynamic amplitude recovery along the arc is not very encouraging. Most amplitude values are in error on the order of 50%, which is not acceptable from a quantitative inversion point of view. The reason for this poor performance is probably due to the fact that the arc structure is too distant from many of the shot gather locations. For example, consider the shot location at x = 1200 m. At this remote distance, the source amplitude A has decayed significantly, s s as can be seen in the ray map of Figure 9.5, to the point that the stabilization parameter 2 e becomes significant. As previously noted, whenever A —> 0, the stabilized dynamic s 2 migration amplitude factor y/A^/(A + e ) approaches the relative amplitude recovery of the s excitation-time factor y/A^, which causes the dynamic amplitude recovery in regions of stabilization to be attenuated relative to the true reflectivity amplitudes. Of course, some shot locations are close enough to recover true amplitudes along the arc, but the net effect of stacking all the shot gather migrations is an overall underestimation of the true arc reflectivity amplitude variation. This is reinforced in Table 10.4 by the fact that the estimated arc amplitudes become more erroneous with increasing raypath distance from the shot location. 167 Figure 10.13 demonstrates the inherent instability involved in using the dynamic imaging condition. Thisfigureconsists of a dynamic migration partial-stack depth trace corresponding 2 to each of the various values of the stabilization parameter e . Trace 3 shows the true relativeamplitude normal-incidence acoustic reflection coefficients as calculated manually from the test model. All traces in thisfigureare normalized relative to the last peak amplitude at 2100 m depth. The stabilization parameter is necessary to stabilize the dynamic amplitude factor, given by the ratio \/A^/A amplitudes A . 8 2 S in the migration equation (4.1), in regions of vanishingly small source The noise parameter variation represented in Figure 10.13 ranges from -1 2 -6 _1 e = 10 on trace 1, to e = 10 on trace 7, in factors of 10 . The stabilized dynamic amplitude factor T/A^/(A 2 S + e ) has the effect of stabilizing the dynamic amplitude division for data imaged from regions of the test model where the source amplitude A has decayed S below the stabilization parameter level. The ray map in Figure 9.5 implies that a stabilization _1 parameter amplitude level of 10 (10 dB atten) effectively implies a stabilization of the dynamic migration amplitude recovery at a radius of 300 m from the source. Hence, true relative-amplitude reflectivity can only be hoped to be recovered within this radius. A 3 stabilization parameter value of 10" (30 dB atten) offers the potential for relative-amplitude recovery in a rectangular area of offset and depth from the source ranging about 1000 m in width and 3000 m in depth. Based on the results of Figure 10.13, a stabilization parameter was simplistically chosen to 2 be e = 0.005, or about 25 dB in the ray map of Figure 9.5. This value is assumed to approximately recover the true relative-amplitude reflectivity of trace 3, based on a simple comparison of the amplitudes and parameter values of trace 2 and trace 4. This stabilization value was used in the dynamic migration of Figure 10.12. Note that for large values of 168 TSN TSN 0 = 0 - - 500o0 r 1000 ? x - - 5 0 0 . 0 - - 1000=0 \ = 0-- f 1 5 0 0 . 0 — 0 . 0 ^ 1500=0 4 2000=0 2000 = 0-- r f f ii r - - 2500.0- 2500=0 3000=0 3000o0- Figure 10.13: Migration Depth Trace Stability Analysis of the Dynamic Reflectivity Imaging Condition for Variable Stabilization Parameter Values e . Trace 3 corresponds to the true relative-amplitude normal-incidence acoustic reflection coefficients derived manually from the arc test model. The remaining traces correspond to the reflectivity amplitude recovery of the dynamic imaging condition for varying stabilization parameter values, ranging from e = 10" on Trace 1, to e = 10 on Trace 7, in factors of 10" . 2 2 1 2 -6 1 the stabilization parameter, e.g., traces 1 and 2, shallow reflectivity amplitudes tend to be correctly recovered, but appear over-gained with respect to the deeper highly-attenuated (stabilized) amplitudes. For small stabilization parameter values, e.g., traces 6 and 7, the deeper coefficients are over-gained relative to the true shallow amplitudes due to the singular nature of the dynamic amplitude factor A division. This instability is also evident in that s the shallower events are subject to serious waveform degradation, which arises at regions of vanishingly-weak source illumination beneath post-critical layer-boundary sections of the migration velocity model. 10.4.2 Geometric Migration Figure 10.14 displays the full-fold Kirchhoff-WKBJ prestack depth migration for the geometric reflectivity imaging condition. Again, the positioning of the 1-D horizontal reflectivity structure and the 2-D arc structure accurately images the reflective boundaries of the true arc test model, Figure 8.1, used to generate the synthetic shot gather data. Similar residual migration-ellipse artifacts distort the reflectivity image as compared with the dynamic migration result of Figure 10.12. The artifacts appear to be slightly larger in amplitude than in the dynamic result, as do the recovered reflectivity amplitudes on the deeper horizontal reflectors. Table 10.5 lists the true reflection coefficients of the horizontal reflectors in the arc model, and the corresponding peak-amplitude reflectivity estimates obtained directly from the geometric migration depth section. The depths of the migrated reflectivity peaks are accurate to within half of one depth sample increment, or + 3.75 m. The reflectivity amplitudes tabulated in Table 10.5 represent a dramatic improvement in amplitude recovery compared with the dynamic migration results of Table 10.2. Especially significant is the recovery of the large 170 Geometric TSN 20 40 GO 80 100 120 I 40 1 GO 180 200 220 240 2G0 Distance/15 (m) Figure 10.14: Full-Fold Kirchhoff-WKBJ Prestack Depth Migration using the Geometric Reflectivity Imaging Condition. A total of 51 shot gathers are migrated for shot locations ranging from xs = 0 m (migrated depth trace 34) to xs = 3000 m (migrated depth trace 234), in increments of 60 m. The 200-receiver array is fixed independent of shot location, and spans nearly the same range as the shot locations: xr = 0 m to xr = 2985 m, in increments of 15 m. TSN Layer True Geometric 1 2 3 4 5 6 1.00 1.00 ± .32 0.98 1.00 + .14 0.60 0.68 ± .05 0.48 0.62 ± .03 -1.75 -1.74 ± .16 1.91 1.51 ± .02 Table 10.5: Normalized 1-D Horizontal-Layer Reflection Coefficients from the True Arc Test Model and the Geometric Migration Depth Section. low-velocity zone reflection coefficient at the base of layer 5, and (nearly) layer 6. Only reflection coefficients 4 and 6 fall markedly outside the estimated error tolerance. However, these erroneous reflectivity amplitudes still maintain the correct trend of the true reflectivity, and at + 25%, are not at all unpalatable. And, compared with the corresponding dynamic amplitude values from Table 10.2, the geometric reflectivity amplitude errors are eminently acceptable. Table 10.6 lists the true normal-incidence acoustic reflection coefficients calculated manually from the dipping arc segments in the arc test model, and the corresponding peak-amplitude reflectivity estimates measured directly from the geometric migration depth section. Again, the geometric reflectivity amplitude estimates of Table 10.6 represent a dramatic improvement compared to the dynamic arc amplitudes of Table 10.4. The geometric amplitude estimates are very accurate considering the extremely challenging 2-D nature of the test model. Even the amplitude in the shallowest layer is quite respectable, considering that most of the segment image has been obliterated by spatial aliasing and the application of the mute pattern and obliquity cosine. These results reinforce the analysis of Chapter 5, which predicts that the geometric imaging condition can achieve accurate reflectivity amplitude recovery in situations where the reference or migration model has mild 2-D variation, even though the reflectivity, 172 Arc True Geometric 1 2 3 4 5 6 2.68 1.86 ± .19 2.06 2.00 ± .20 1.24 1.03 ± .10 0.68 0.69 ± .07 0.21 0.21 ± .02 1.91 1.51 ± .15 Table 10.6: Normalized 2-D Dipping Arc Segment Reflection Coefficients from the True Arc Test Model and the Geometric Migration Depth Section. or "scattering" structure may be strongly 2-D. In conclusion, the geometric reflectivity amplitude recovery results, for both the 1-D and 2-D structures in this example, prove to be superior to all other imaging conditions examined. Hence, Tables 10.5 and 10.6 verify the previous analysis which predicts that the new geometric reflectivity imaging condition, which arises from this thesis research, represents an optimal compromise between reflectivity image resolution (stability) and reflectivity amplitude estimation accuracy. The combined requirements for accurate imaging and amplitude recovery make the geometric imaging method a prime candidate for optimal reflectivity model construction and impedance inversion. 10.4.3 Kinematic Migration Figure 10.15 displays the full-fold Kirchhoff-WKBJ prestack depth migration for the kinematic reflectivity imaging condition. The horizontal and arc reflectivity structure are accurately imaged. The migration-ellipse artifacts appear to be slightly decreased in amplitude 173 Kinematic TSN 20 40 GO 80 1 00 1 20 t 40 1 60 1 80 200 220 240 260 Distance/15 (m) Figure 10.15: Full-Fold Kirchhoff-WKBJ Prestack Depth Migration using the Kinematic Reflectivity Imaging Condition. A total of 51 shot gathers are migrated for shot locations ranging from \ = 0 m (migrated depth trace 34) to xs • 3000 m (migrated depth trace 234). The 200-receiver array is fixed independent of shot location, and spans nearly the same range as the shot locations: xr = 0 m to xr = 2985 m. s TSN Layer True Kinematic 1 2 3 4 5 6 1.00 1.00 ± .32 0.98 0.97 ± .13 0.60 0.57 + .04 0.48 0.47 ± .02 -1.75 -1.21 ± .11 1.91 0.97 ± .02 Table 10.7: Normalized 1-D Horizontal-Layer Reflection Coefficients from the True Arc Test Model and the Kinematic Migration Depth Section. relative to the geometric migration result of Figure 10.14. Table 10.7 lists the true reflection coefficients of the horizontal reflectors in the arc model, and the corresponding peakamplitude reflectivity estimates obtained directly from the kinematic migration depth section. The depths of the migrated reflectivity peaks are accurate to within half of one depth sample increment, or ± 3.75 m. Table 10.7 demonstrates that the kinematic migration possesses better 1-D reflectivity amplitude recovery properties than the stabilized dynamic imaging condition, but is inferior to the geometric imaging condition. In fact, the kinematic imaging condition slightly outperforms the geometric imaging condition over thefirstfour reflection coefficients, but then seems to deteriorate rapidly for the last two estimates. The kinematic migration probably performs well for the first four depth levels, because the 1-D straight-ray approximation, for which the kinematic migration is suitably accurate, is fairly valid. However, when the velocity structure changes rapidly at the low velocity zone bracketed by the last two reflection coefficients, the straight-ray assumption is severely violated, and the kinematic imaging condition seems to fail accordingly. Table 10.8 lists the true normal-incidence acoustic reflection coefficients calculated manually 175 Arc True Kinematic 1 2.68 0.88 ± .09 2 2.06 1.25 ± .13 1.24 0.73 ± .07 3 4 0.68 0.52 ± .07 5 0.21 0.16 + .02 6 1.91 0.94 + .09 Table 10.8: Normalized 2-D Dipping Arc Segment Reflection Coefficients from the True Arc Test Model and the Kinematic Migration Depth Section. from the dipping arc segments in the arc test model, and the corresponding peak-amplitude reflectivity estimates measured directly from the kinematic migration depth section. The amplitude variation along the arc structure is fairly good in a qualitative sense. In a quantitative sense, the kinematically recovered amplitudes are not as good as the geometric imaging condition, but better than the dynamic imaging, which suffers from its inherent instability. Hence, the kinematic amplitude recovery is quite good for 1-D reflectivity structure, as suggested by the results in Table 10.7, but not reliable for 2-D structure, as evidenced by Table 10.8. This is just as predicted by the analysis of Chapter 5. A close scrutiny of the kinematic results in Table 10.8 suggests that the recovered 2-D reflectivity amplitude variation along the arc appears to be more conservative, or balanced than either of the geometric or dynamic migration results. This balanced amplitude quality 2 can be quantified somewhat, by calculating the variance a of the recovered arc amplitude values separately for each column 3 of Tables 10.4,10.6 and 10.8 respectively. The kinematic 2 arc amplitude variance, o\ = 0.141, is much less than that of the geometric variance, a , = 0.487, the dynamic variance, ad = 0.373, or the variance of the true arc segment reflection 2 coefficients, a = 0.854. These results reinforce previous analysis which predicts that the kinematic imaging condition migration provides a useful reflectivity image, consisting of 176 fairly accurate amplitude recovery along 1-D structure, and a gradual amplitude-balancing trade-off between amplitude accuracy and image resolution along 2-D structure. 10.4.4 Excitation-Time Migration Figure 10.16 displays the full-fold Kirchhoff-WKBJ prestack depth migration for the excitation-time reflectivity imaging condition. The horizontal reflectivity structure is fairly well imaged, although the image tends to fade with depth. The depth attenuation of the imaged amplitudes contributes to a weak, and in some places, barely discernable image of the arc structure. This is the effect of the uncompensated source divergence A that is buried in the s excitation-time image, as predicted by the analysis of Chapter 5. Table 10.9 lists the true reflection coefficients of the horizontal reflectors in the arc model, and the corresponding peak-amplitude reflectivity estimates obtained directly from the excitationtime migration depth section. The depths of the migrated reflectivity peaks are accurate to within half of one depth sample increment, or ± 3.75 m. The fourth column represents the directly measured excitation-time amplitudes of the third column, after application of a geometric divergence correction. The normal-incidence divergence correction is calculated using a simple relation derived from Equation (5.9), or that of O'Brien and Lucas (1971, p.13), n<N A \zn) = J2 * k > k=\ l C /Cl ( i a i ) D which is valid for a 1-D earth consisting of N horizontal constant-velocity layers, and where Ik and Ck are the thickness and velocity respectively of the k 177 th layer. Table 10.10 lists the Excitation-Time TSN 20 40 60 80 100 120 140 160 180 200 220 240 260 TSN Distance/15 (m) Figure 10.16: Full-Fold Kirchhoff-WKBJ Prestack Depth Migration using the Excitation-Time Reflectivity Imaging Condition. A total of 51 shot gathers are migrated for shot locations ranging from xs = 0 m (migrated depth trace 34) to xs = 3000 m (migrated depth trace 234). The 200-receiver array is fixed independent of shot location, and spans nearly the same range as the shot locations: xr - 0 m to xr = 2985 m. Layer True Excitation-TimeDiv. Corrected 1 2 3 4 5 6 1.00 0.98 0.60 0.48 -1.75 1.91 1.00 ± .35 0.62 ± .13 0.23 ± .03 0.14 ± .01 -0.29 ± .03 0.27 ± .05 1.00 + .35 1.09 ± .22 0.68 + .09 0.62 ± .06 -1.84 ± .17 1.95 ± .33 Table 10.9: Normalized 1-D Horizontal-Layer Reflection Coefficients from the True Arc Test Model and the Excitation-Time Migration Depth Section. Div. Factor Absolute Normalized 600. 1.000 F F 1060. 1.767 1760. 2.933 F F 2660. 4.433 3760. 6.267 F 4260. 7.100 Fe 1 2 3 4 5 Table 10.10: Absolute and Normalized Normal-Incidence Geometric Divergence Correction Factors F Calculated from the True Arc Test Model. n normal-incidence geometric divergence correction factor, F =A (z ), calculated manually 1 n D n from the true arc test model at the six depths corresponding to the reflection coefficients of Table 10.9. Table 10.9 demonstrates the potential for significant error in reflectivity amplitude estimation incurred with the uncorrected excitation-time imaging condition. Even thefirstpiece of independent uncorrected amplitude information, at the second reflector, is unreliable, and the uncorrected amplitude accuracy is subject to further deterioration with increasing depth. The analysis of Chapter 5 predicts that excitation-time amplitude estimates may be corrected post-migration to remove the effect of the source amplitude term A . In light of this, column s 179 4 of Table 10.9 lists the amplitudes of column 3, after a true normal-incidence geometric divergence correction, calculated directly from the test model, has been applied. The corrected excitation-time amplitudes are about as accurate as the geometric imaging condition results, except that the deepest reflection coefficient at 2100 m is better estimated in the former case. This is a very encouraging result which implies that accurate amplitude estimates can be recovered, in the vicinity of 1-D reflectivity structure, by using the non-true-amplitude excitation-time imaging condition migration supplemented by a judicious geometric divergence correction. Table 10.11 lists the true normal-incidence acoustic reflection coefficients calculated manually from the dipping arc segments in the arc test model, and the corresponding peak-amplitude reflectivity estimates measured directly from the excitation-time migration depth section. The uncorrected amplitude estimates are dramatically incorrect. The accuracy of the divergence-corrected amplitude estimates along the arc structure is also discouraging, although the deepest three reflection coefficients are approaching somewhat reasonable values, in a qualitative sense. This is probably due to the fact that, as the reflected arrivals return from deeper, less obliquely-dipping structure, the raypath geometry becomes less strongly 2-D. The results of Tables 10.9 and 10.11 suggest that reflected events which have a nearly 1-D raypath geometry, tend to be accurately amplitude-estimated by the excitation-time imaging condition supplemented by standard geometric divergence correction. Evidently, dipping reflectors which exhibit clearly 2-D raypath geometry cannot be accurately amplitude-estimated with the excitation-time imaging condition or a standard geometric divergence correction. Hence, the combination of excitation-time amplitude recovery and geometric divergence correction is quite good for 1-D reflectivity structure, as suggested by the results in Table 10.9, but not reliable for 2-D structure, as evidenced by Table 10.11. 180 Arc 1 2 3 4 5 6 True Excitation-TimeDiv. Corrected 2.68 2.06 1.24 0.68 0.21 1.91 0.67 0.37 0.20 0.12 0.03 0.19 ± .07 + .04 ± .02 ± .01 ± .00 + .02 0.22 0.51 0.46 0.42 0.18 1.26 ± .02 ± .05 + .05 + .04 ± .02 ± .13 Table 10.11: Normalized 2-D Dipping Arc Segment Reflection Coefficients from the True Arc Test Model and the Excitation-Time Migration Depth Section. As an ancillary remark, it can be noted that Figure 10.16, and column 3 of Table 10.9, indicate the approximate signal-to-noise ratio expected from the excitation-time imaging condition. Gaining up the excitation-time migration results with a geometric divergence correction factor will increase both signal and noise levels in the corrected results. The presence of noise in the data may seriously degrade even the corrected excitation-time amplitude estimates in regions of low signal-to-noise ratio. In contrast, the results tabulated for the geometric, kinematic, and dynamic imaging conditions do not require a gain correction, and so their values indicate better signal-to-noise recovery than the corrected excitation-time results. The geometric divergence factors in Table 10.10 indicate the approximate s/n recovery of the corrected excitation-time amplitudes relative to the true-amplitude geometric, kinematic and dynamic formulations. For example, Table 10.10 suggests that the excitation-time s/n recovery of the deepest reflector amplitude at 2100 m depth may be about seven times worse than the associated geometric, kinematic, or dynamic result. In conclusion, nearly all true relative-amplitude information is destroyed by the uncorrected excitation-time imaging condition. However, it may be possible to approximately compensate for the source-divergence term A s by the application of a pre- or post-migration divergence correction, in order to produce a better structural image, and a fairly accurate estimate of 1-D, 181 but not 2-D, reflectivity amplitudes. Also, as most of the true reflective structure appears in the uncorrected image, an application of AGC, for example, may produce a useful map of subsurface reflectors for structural interpretation applications. Finally, the excitation-time imaging condition suffers from inherently poor signal-to-noise recovery, which may further degrade the reflectivity image and corrected amplitude estimates relative to the geometric, kinematic, or dynamic imaging conditions, depending upon the ambient noise level in the recorded seimic data. 10.4.5 Crosscorrelation Migration Figure 10.17 displays the full-fold Kirchhoff-WKB J prestack depth migration for the crosscorrelation reflectivity imaging condition. This dramaticfigureillustrates the poor image resolution and amplitude accuracy inherent in the crosscorrelation imaging condition. Only the shallowest horizontal reflectivity structure is imaged, and even this tends to fade rapidly with depth. The extreme depth attenuation of the imaged reflectivity amplitudes results in the complete loss of the arc structure from the migrated depth section. This deleterious effect is due to the lack of compensation for the square of the source divergence Al, which remains imbedded in the crosscorrelation image, as predicted by the analysis of Chapter 5. Table 10.12 lists the true reflection coefficients of the horizontal reflectors in the arc model, and the corresponding peak-amplitude reflectivity estimates obtained directly from the crosscorrelation migration depth section. The depths of the migrated reflectivity peaks are accurate to within half of one depth sample increment, or ± 3.75 m. The fourth column represents the direcdy measured crosscorrelation amplitudes of the third column, after application of a geometric divergence correction. The correction factor attempts to compensate for the square of the source amplitude term, A , and is implemented by applying the square of the F values 2 S n 182 Crosscorrelation TSN 20 40 60 80 100 120 140 160 180 200 220 240 260 TSN Distance/15 (m) Figure 10.17: Full-Fold Kirchhoff-WKBJ Prestack Depth Migration using the Crosscorrelation Reflectivity Imaging Condition. A total of 51 shot gathers are migrated for shot locations ranging from xs = 0 m (migrated depth trace 34) to xs = 3000 m (migrated depth trace 234). The 200-receiver array is fixed independent of shot location, and spans nearly the same range as the shot locations: xr = 0 m to xr = 2985 m. Layer 1 2 3 4 5 6 True Crosscorrelation Div. Corrected 1.00 0.98 0.60 0.48 -1.75 1.91 1.00 ± .40 0.25 ± .07 0.04 ± .01 0.01 + .00 -0.02 ± .00 0.02 ± .01 1.00 ± .40 0.77 ± .23 0.34 ± .06 0.26 ± .04 -0.67 ± .08 1.11 + .35 Table 10.12: Normalized 1-D Horizontal Layer Reflection Coefficients from the True Arc Test Model and the Crosscorrelation Migration Depth Section. of Table 10.10 to each associated uncorrected crosscorrelation amplitude in column 3. Table 10.12 demonstrates the dramatic error in reflectivity amplitude estimation incurred with the uncorrected crosscorrelation imaging condition. Evidently, none of the uncorrected crosscorrelation reflectivity amplitude estimates are even remotely reliable. Such significant amplitude error represents an accelerated degradation of the reflectivity estimates with depth as compared to the uncorrected excitation-time migration of Table 10.9. The analysis of Chapter 5 predicts that crosscorrelation amplitude estimates may be corrected post-migration to remove the effect of the square of the source amplitude term, A .. In light 2 of this, column 4 of Table 10.9 lists the amplitudes of column 3, after a squared normal-incidence geometric divergence correction has been applied. Qualitatively, the divergencecorrected crosscorrelation reflectivity amplitudes recover the approximate trend of the true test model reflectivity. However, in a quantitative sense, the corrected crosscorrelation amplitudes are not as accurate as the associated geometric, kinematic, dynamic, or excitation-time results. The corrected crosscorrelation amplitude estimates are certainly not accurate enough for subsequent quantitative analysis or impedance inversion. Hence, these results suggest that a corrected crosscorrelation migration may be useful as an qualitative imaging process, but 184 Arc 1 2 3 4 5 6 True Crosscorrelation Div. Corrected 2.68 2.06 1.24 0.68 0.21 1.91 0.77 ± .08 0.03 ± .00 0.01 ± .00 0.01 ± .00 0.00 ± .00 0.01 + .00 0.09 ± .01 0.05 + .00 0.07 + .01 0.08 ± .01 0.04 ± .00 0.52 ± .05 Table 10.13: Normalized 2-D Dipping Arc Segment Reflection Coefficients from the True Arc Test Model and the Crosscorrelation Migration Depth Section. not a quantitative inversion process. Table 10.13 lists the true normal-incidence acoustic reflection coefficients calculated manually from the dipping arc segments in the arc test model, and the corresponding peak-amplitude reflectivity estimates measured directly from the crosscorrelation migration depth section. Both the uncorrected and divergence-corrected amplitude estimates along the arc structure are dramatically incorrect. Evidently, the errors involved in estimating the 1-D crosscorrelation amplitudes are even further magnified in the presence of strong 2-D structure, and cannot be unravelled by the application of a standard divergence correction. Finally, in the presence of noise, the crosscorrelation migration images may tend to suffer from the extremely poor s/n recovery of this method, indicated by Figure 10.17 and Table 10.12. Based on the squared correction factors of Table 10.10, the crosscorrelation s/n recovery of the deepest reflector amplitude at 2100 m depth may be about 50 times worse than the associated geometric, kinematic, or dynamic result. This extremely poor level of s/n recovery may seriously degrade even the imaging aspect of the crosscorrelation migration method, relative to the other imaging conditions. 185 10.5 Geometric Migration Stability in the Presence of Noise The purpose of this section is to examine the stability, or robustness, of the geometric imaging condition migration method in the presence of noise. A single shot gather migration, using the geometric imaging condition, will be tested for robustness of reflectivity imaging and amplitude recovery. In light of previous results, the single shot gather results can be extrapolated to the analogous full-fold situation. A full-fold migration example in the presence of noise is not presented here due to the large additional computational effort required, and the somewhat redundant nature of the shot-gather and multi-fold results. Figure 10.18 displays a "noisy" synthetic shot gather located at x = 1200 m along the arc s test model. The noisy shot gather is obtained by adding evenly-distributed random noise to the "clean" shot gather used in the previous migration examples, shown in Figure 8.3. The standard deviation of the evenly-distributed random noise is set at 10% of the peak absolute amplitude in the gather. Hence, the s/n level ranges from about 10:1 at the earliest reflected event, to less than 1:1 at the later, weaker arrivals. This large level of additive noise represents a challenging test of the stability of the 2-D imaging and amplitude recovery properties of the geometric imaging condition, and is clearly comparable, or greater than, the noise levels observed in typical seismic reflectionfielddata. Figure 10.19 displays the geometric migration of the noisy shot gather in Figure 10.18. This noisy migration is directly comparable with the analogous clean shot gather migration of Figure 10.8. To the experienced eye, it is immediately evident that the migrated result contains a severe amount of noise contamination. The migrated noise is swept out along migration ellipses corresponding to potential, but erroneous, reflector locations. This effect is familiar to data processors and interpreters as "migration smile." The different appearance 186 Distance/15 (m) Figure 10.18: The Muted Shot Gather at x = 1200 m with 10% Additive Noise. This noisy data is obtained by adding uniformly distributed random noise to the clean shot gather data of Figure 8.3. The standard deviation of the noise is specified to be 10% of the peak absolute amplitude value in the clean data. Hence, the s/n level ranges from about 10:1 on the shallowest reflection event, to less than 1:1 on deeper reflection events. s 187 Distance/15 (m) Figure 10.19: Geometric Migration of a Noisy Shot Gather. The single shot gather data correspond to a shot position at x = 1200 m, marked by the star symbol, and 200 receiver positions, ranging from x = 0 m (migrated depth trace 34) to x = 2985 m (migrated depth trace 233), in 15 m increments. s r r of the migrated noise at 2000 m depth is due to the large contrast of the low-velocity zone in the migration model at this depth. This layer has a velocity of 2500 m/s, and is bound by high-velocity layers at 5500 m/s and 6000 m/s above and below respectively. The banded appearance arises from the strong change in wavefront direction, in accordance with Snell's Law, across the low-velocity zone boundaries. Also, the spatial wavelength of waveforms, given by the relation X = c/f, is much smaller in a low-velocity region than in a high-velocity region. Hence, recalling that Figure 10.19 is a depth migration, the spatial frequency of noise in the low-velocity zone is much higher than that of adjacent layers. The image of the 1-D reflectors, and to a somewhat lesser degree, the image of the arc structure, is readily observable. This is extremely encouraging considering that the reflected trace arrivals corresponding to the deep 1-D and 2-D structure in the shot gather of Figure 10.18, are barely discernable. Evidently the intrinsic wavefield divergence correction given by the geometric amplitude factor rs/^/7> is extremely robust in that it amplifies the trace data just enough to recover the correct reflectivity image signal, but not enough to destroy the coherency of the image by overemphasizing the noise. Trace 3 of Figure 10.20 displays a depth trace from the noisy geometric migration, along with a corresponding clean migration depth trace, and the true normal-incidence acoustic reflection coefficients. The traces have been normalized to thefirstpeak amplitude at 600 m on each trace. The geometric migration depth traces are obtained by a partial stack along the 1-D reflectivity structure, in a manner similar to that described previously. Table 10.14 lists the peak reflectivity amplitudes directly measured from each depth trace. The depths of the migrated reflectivity peaks are accurate to within half of one depth sample increment, or ± 3.75 m. In a qualitative sense, the depth traces of Figure 10.20 show the extreme imaging robustness 189 500o0 - - 500.0 1000 = 0-- 1000.0 1500.0-- — 1500.0 2000.0 2000.0 2500.0 -- 2500.0 3000.0 3000.0 Figure 10.20: Normalized Depth Traces Corresponding to the True Normal-Incidence Acoustic Reflection Coefficients of the Arc Test Model, and the Geometric Shot Gather Migration of Clean and Noisy Synthetic Data. Trace 1 corresponds to the true arc test model reflection coefficients, Trace 2 to the geometric migration of tlie clean shot gather, and Trace 3 to the geometric migration of the noisy shot gather. All traces are normalized individually to thefirstreflection event at 600 m depth. Layer True 1 2 3 4 5 6 1.00 0.98 0.60 0.48 -1.75 1.91 Clean Noisy 1.00 + .31 0.97 ± .14 0.68 ± .05 0.62 ± .03 -1.74 ± .16 1.51 ± .02 1.00 ± .35 0.93 ± .10 0.70 ± .06 0.62 ± .04 -1.76 ± .14 1.39 ± .08 Table 10.14: Normalized 1-D Horizontal Layer Reflection Coefficients from the True Arc Test Model and the Clean and Noisy Geometric Migration Depth Traces. of the geometric migration method in the presence of large amounts of random noise. The s/n recovery of the shot gather depth traces, obtained by partial stacking, is indicative of the the full-fold s/n recovery. The large degree of s/n recovery by the geometric migration method is very encouraging indeed. In a quantitative sense, Table 10.14 demonstrates the inherent stability of the geometric imaging condition, which has recovered the true relative-amplitude reflectivity variation of the test model with a high degree of accuracy comparable to that of the clean geometric migration. In conclusion, the geometric migration method is extremely robust in terms of reflectivity imaging and amplitude recovery in the presence of noise. 10.6 Geometric Migration Stability in the Presence of Velocity Model Errors The purpose of this section is to examine the stability of the geometric migration method in the presence of migration velocity model errors. A single shot gather migration, using the geometric imaging condition, will be tested for robustness of reflectivity imaging and amplitude recovery. Again, in light of previous results, the single shot gather results may be extrapolated to the analogous full-fold situation, which will not be presented here. 191 The input data to the analysis of this section is again the noisy shot gather of Figure 10.18. However, instead of using the correct migration velocity model described in Chapter 8, we desire to implement the geometric migration of the noisy shot gather with an erroneous migration model. It is difficult to decide how erroneous the migration model can be, and yet still represent a fair test of stability. It was decided that the most reasonable specification of an erroneous migration velocity model would be obtained from a velocity semblance scan of a noisy shot gather. In the velocitysemblance approach, a shot gather, or CMP gather, is scanned along a suite of hyperbolic NMO trajectories each specified by a depth and stacking velocity. The stacking velocities are then manually "picked" as a function of depth from the peak amplitudes in the velocity-depth semblance map. Assuming that the stacking velocity function is approximately equivalent to the rms velocity function, interval velocity can be computed as a function of depth by an application of the Dix Equation, for example. The velocity semblance approach to determining a migration velocity model is probably the most commonly used method by the practicing seismic data processor, in the face of "real"fielddata and gross uncertainty regarding the true nature of the Earth's subsurface. For this reason, obtaining a migration velocity model by the semblance method provides both a realistic and reasonable test of migration stability in the presence of migration velocity model errors. Table 10.15 lists the results of performing the described velocity semblance analysis on the noisy shot gather of Figure 10.18. The velocities have units of meters/second and the depths are listed in meters. The erroneous migration velocity model is constructed from the Csemb interval velocities in the Table resulting from the semblance analysis. It is evident that significant error has been incurred in the semblance estimates of the correct interval velocities, and to a somewhat lesser degree, the depths to the base of each constant-velocity 192 Layer 1 2 3 4 5 6 7 Ctrue Csemb 1500. 2300. 3500. 4500. 5500. 2500. 6000. 1500. 2336. 3712. 4480. 5338. 2480. AC 0. 36. 212. -20. -162. -20. - - Zfrue ^semb AZ 600. 900. 1200. 1500. 1800. 2100. 3000. 600. 899. 1229. 1520. 1814. 2112. 0. -1. 29. 20. 14. 12. - - Table 10.15: Migration Model Velocities and Depths from the True Arc Test Model and the Velocity Semblance Analysis of a Noisy Shot Gather. layer. Figure 10.21 displays the geometric migration of the noisy shot gather in Figure 10.18 using the erroneous migration velocity model tabulated in Table 10.15. The deeper portions of the migrated depth section are of a higher spatial frequency content than Figure 10.19, which is the analogous migration in the presence of an accurate migration velocity model. This is because the low velocity zone at 2000 m depth is the last velocity information available from the semblance analysis, and so is extended to 3000 m depth. Note that this does not affect the reflectivity image, since no reflectors exist below 2100 m. Some slight depth warping of horizontal reflectors is observable, compared to Figure 10.19, due to the velocity and depth errors of the inaccurate migration model. Otherwise, the migrated image is very similar in a qualitative sense to the accurate model migration of Figure 10.19. Figure 10.22 displays the partial-stack migrated depth traces of Figure 10.20 with the addition of a depth-trace corresponding to the erroneous migration velocity model geometric shot gather migration of Figure 10.21. Trace 1 represents the true 1-D normal-incidence acoustic reflection coefficients of the arc test model. Trace 2 results from the geometric migration of the clean shot gather with an accurate migration velocity model. Trace 3 results from 193 TSN 20 40 60 80 100 : 2 0 140 160 180 200 220 240 260 TSN SO Distance/15 (m) Figure 10.21: Geometric Migration of a Noisy Shot Gather with an Erroneous Migration Velocity Model. TSN 1 2 3 4 TSN Figure 10.22: Normalized Depth Traces Corresponding to the True Normal-Incidence Acoustic Reflection Coefficients of the Arc Test Model, and the Geometric Shot Gather Migrations. Trace 1 corresponds to the true arc test model reflection coefficients, Trace 2 to the geometric migration of the clean shot gather data with the correct migration velocity model, and Trace 3 to the geometric migration of the noisy shot gather data with the correct migration velocity model, and Trace 4 to the geometric migration of the noisy shot gather data with an erroneous migration velocity model. All traces are normalized individually to the first reflection event at 600 m depth. Layer True 1 2 3 4 5 6 1.00 0.98 0.60 0.48 -1.75 1.91 Clean Noisy 1.00 ± .31 0.97 ± .14 0.68 + .05 0.62 ± .03 -1.74 + .16 1.51 ± .02 1.00 + .35 0.93 + .10 0.70 + .06 0.62 + .04 -1.76 ± .14 1.39 + .08 Vel. Semb. AZ 1.00 + .32 0.94 + .02 0.57 + .04 0.62 + .04 -1.73 ± .23 1.33 + .27 0. 8. 30. 30. 15. 15. Table 10.16: Normalized 1-D Horizontal Layer Reflection Coefficients from the True Arc Test Model and the Geometric Migration Depth Traces. the geometric migration of the noisy shot gather with an accurate migration velocity model. And, Trace 4 results from the geometric migration of the noisy shot gather with an inaccurate migration velocity model. Similarly, Table 10.16 lists the peak reflection coefficient amplitudes of the depth traces in Table 10.14, with the addition of the peak amplitudes corresponding to the erroneous migration velocity model geometric migration depth trace. The last column, AZ, lists the depth errors (m) in the location of the peak reflection coefficient amplitudes corresponding to the velocity semblance migration model example of Trace 4. These migrated depth errors arise from a combination of both velocity and depth errors in the semblance migration velocity model of Table 10.15. Evidently the migrated depth errors are on the order of the depth errors originally specified in the inaccurate semblance migration velocity model. Comparison of the depth traces of Figure 10.22 and their peak-amplitude reflection coefficients tabulated in Table 10.16, demonstrates that the geometric imaging condition is very robust with respect to accurate reflectivity imaging and amplitude recovery in the presence of significantly large data noise levels and migration velocity model errors. 196 10.7 Computational Requirements Computational requirements will now be discussed for the single shot gather migration examples presented in this thesis. Requirements for larger seismic reflection data sets and output depth images can be extrapolated from this discussion accordingly. The critical computational factors for the migration algorithm, excluding raytracing requirements which are listed in Chapter 9, are: the number of traces per shot gather A^, the dimensions of the migrated output depth image N and N , and a scale factor Ni relating to x z the type of imaging condition applied. For the examples presented here, Nt = 200, N = 268, x and N = 401. z The kinematic imaging condition is the next most efficient migration method presented here, because most calculations are based solely on the raytraced traveltimes, and not the WKBJ amplitudes. For the kinematic imaging condition, Ni = 1.0. A single trace migration for these examples requires about 15 seconds CPUtimeon a MicroVax JJ with 9 Mb memory. The full shot gather requires about 50 minutes CPU time. These CPU times can realistically be decreased by a factor of 10 with the use of an array processor. Also, I have vectorized the kinematic migration algorithm to run on a Cray-1, which dramatically reduces the migration CPUtimeto about 20 seconds per shot gather. Since only raytraced traveltimes are needed to implement the kinematic imaging condition, only N raytraced maps need be computed. This saves considerable CPU time due to s raytracing, and CPU time spent swapping the usually large ray maps in and out of core memory. N is the appropriate number of independent surface locations, as discussed in s Chapter 9. For a 1-D migration model, N = l. In general, an arbitrary 2-D model requires s 197 that N equal the total number of independent receiver locations along the full length of the s survey to be migrated, which may be considerably large. This aspect and some shortcuts are discussed further in Chapter 9. The geometric imaging condition is the next efficient migration method tested here. For the geometric imaging method, the degree of calculation effort is given by iV; = 1.7, or 1.7 times more intensive than the kinematic method. A single trace migration for these examples requires about 25 seconds CPU time on a MicroVax II with 9 Mb memory. The full shot gather requires about 80 minutes CPUtime.Again, these CPUtimescan realistically be decreased by a factor of 10 with the use of an array processor. Also, I have vectorized the geometric migration algorithm to run on a Cray-1, which reduces the migration CPU time to about 30 seconds per shot gather. The raytracing requirement for the geometric migration method is approximately equivalent to the kinematic method, since the WKBJ amplitude maps are replaced with a single map of radial shot-subsurface distances r. The extra radial distance map effectively makes the geometric raytracing requirement proportional to N + 1. However, the migration algorithm s does perform many internal amplitude calculations, unlike the kinematic algorithm, and this fact is represented in the Ni = 1.7 scalar. Finally, the remaining dynamic, excitation-time, and crosscorrelation imaging conditions are equivalently inefficient. They require approximately the same migration CPU time as the geometric imaging condition. However, the need for both traveltime and WKBJ amplitude maps constitutes a large raytracing CPU requirement proportional to 2N , as opposed to N S s andiVs+l for the kinematic and geometric migrations respectively. This additional raytracing requirement also doubles the CPU time required for core memory swapping of ray maps in the migration algorithm. 198 As a final note of interest, I have indirectly benchmarked my vectorized Kirchhoff-WKBJ algorithm against a vectorizedfinite-differencealgorithm for both the excitation-time and crosscorrelation migration methods. The Kirchhoff implementation is about one order of magnitude faster for the excitation-time imaging condition, and two orders of magnitude faster for the crosscorrelation imaging condition. The main differences in the two implementations are: (1) the Kirchhoff wavefield reconstruction is done by raytracing Green's functions and integral summation, versusfinite-differencewavefield extrapolation, and, (2) the Kirchhoff wavefield crosscorrelations are performed analytically and reduce to simple ray map evaluations, versus the many explicit 2-Dfinite-differencewavefield crosscorrelations of the alternate approach. 10.8 Summary of the Migration Imaging Condition Quantitative Analysis Figure 10.23 displays the uncorrected migration depth traces corresponding to each reflectivity imaging condition. Each depth trace is obtained as a partial stack of 11 traces, from x = 685 r m to x = 835 m, from the corresponding full-fold imaging condition migration depth section. r Trace 1 represents the true relative-amplitude normal-incidence acoustic reflection coefficients calculated manually from the original arc test model material property specification. The migration imaging condition depth-traces are arranged from left to right in terms of decreasing accuracy of reflectivity amplitude recovery. Starting from the right, Trace 6 represents the uncorrected crosscorrelation imaging condition, which is relatively inferior to the other imaging conditions, both in terms of image resolution and amplitude accuracy. The crosscorrelation method is extremely unsuitable for accurate recovery of reflectivity amplitudes, and hence is more useful as an image processing tool than 199 5 0 0 n 500 = 0 Q 1000=0 1000 = 0— 1500 = 0— — 1500=0 2 0 0 0 a 0 ~ 2000=0 2 5 0 0 = 0- 2500o0 3000 3000=0 Figure 10.23: Full-Fold Kirchhoff-WKBJ Prestack Depth Migration Traces for each of the Five Reflectivity Imaging Conditions. Trace 1 corresponds to the true reflection coefficients calculated from the arc test model, Trace 2 is the depth migrated reflectivity estimate using the geometric imaging condition, Trace 3 the kinematic imaging condition, Trace 4 the stabilized dynamic imaging condition, Trace 5 the excitation-time imaging condition, and Trace 6 the crosscorrelation imaging condition. All traces are normalized individually to the first reflection event at 600 m depth. for quantitative data analysis or inversion. The crosscorrelation method also exhibits very poor signal-to-noise recovery, relative to the other imaging conditions, which may seriously degrade the imaging ability of the method in the presence of noise. The crosscorrelation imaging condition is a currently predominant method of prestack migration. Trace 5 represents the uncorrected excitation-time imaging condition, which exhibits fairly good imaging properties, but poor amplitude recovery. The imaging and amplitude recovery of the excitation-time method can be dramatically improved by the application of a standard divergence correction. The divergence-corrected excitation-time amplitudes are very accurate near regions of 1-D structure, but inaccurate for 2-D structure. Hence, the excitation-time method represents an extremely good imaging process, which is also useful for quantitative data analysis and inversion of 1-D reflectivity amplitudes. The excitation-time method possesses somewhat poor signal-to-noise recovery, relative to the geometric, kinematic and dynamic imaging conditions, which may degrade the imaging and amplitude recovery ability of the method in the presence of noise. The excitation-time imaging condition is a currently predominant method of both post-stack and prestack migration. Trace 4 corresponds to the new Kirchhoff-WKBJ implementation of the dynamic imaging condition, which produces excellent imaging results, and although somewhat unstable, has the potential to yield good quantitative 2-D amplitude accuracy when judiciously stabilized. In practice, the dynamic imaging condition is found to recover only fair estimates of 1-D reflectivity amplitudes, and less accurate 2-D amplitudes. The computational intensity and instability of this method probably makes it less useful than the more accurate and efficient alternate methods presented here. Trace 3 corresponds to the newly proposed kinematic imaging condition, which exhibits 201 excellent imaging and amplitude recovery properties in 1-D reflectivity structure, and a highquality "amplitude-balanced" image at the gradual expense of true-amplitude information in regions of 2-D structure. Hence, the kinematic imaging condition represents an extremely efficient and robust imaging process, which is also useful for quantitative data analysis and inversion of 1-D reflectivity amplitudes. Trace 2 corresponds to the newly proposed geometric imaging condition, which exhibits accurate and efficient reflectivity image resolution and quantitative amplitude recovery in the presence of both 1-D and 2-D reflectivity structure. Hence, the geometric imaging condition represents an optimal method for accurate reflectivity model construction, and should be extremely valuable for quantitative seismic reflection data analysis and impedance inversion of 2-D reflectivity structure. Furthermore, the geometric imaging condition is found to be very robust in the presence of large amounts of noise and migration velocity model errors. Finally, the geometric imaging condition has been tested on "real"fielddata, to be presented elsewhere, with similar conclusions. 202 Chapter 11 Conclusion The primary goal of this thesis research is reflectivity model construction. Philosophi- cally, the problem of reflectivity model construction can be decomposed into two related problems: reflectivity imaging, and reflectivity amplitude estimation. Traditionally, seismic wave-equation migration research has focussed on the imaging aspect of reflectivity model construction. Migration has proved to be valuable for qualitative analysis of seismic reflection data, especially in the area of seismic structural interpretation. However, the trend of current geophysical research appears to be increasingly directed to the quantitative analysis and inversion of geophysical data, in order to construct quantitative, rather than qualitative, models describing the physical nature of the Earth. The quest for quantitative geophysical information is the motivation for my present thesis research. Hence, this thesis addresses the more complete reflectivity model construction problem, by developing methods which perform simultaneous reflectivity imaging and amplitude-estimation, within the context of seismic wavefield reconstruction. Analysis is performed onfivereflectivity "imaging conditions" with respect to qualitative image resolution and quantitative amplitude accuracy of the recovered reflectivity model estimates. A new computationally efficient and stable prestack depth migration method is proposed which is based upon a geometric approximation to the theoretically correct, but 203 unstable, "dynamic" imaging condition. The "geometric" imaging condition has the desirable property of true-amplitude reflectivity recovery in regions of both 1-D and 2-D velocity variation, while fully retaining and optimizing the favorable imaging characteristics of current, non-true amplitude migration formulations. The currently predominant "crosscorrelation" and "excitation-time" migration imaging methods are shown to possess significantly less accurate imaging and amplitude-recovery characteristics relative to the proposed geometric method. The respective signal-to-noise recovery of their imaged amplitudes deteriorates approximately linearly (excitation-time) and quadratically (crosscorrelation) with depth. As a necessary prerequisite to the reflectivity model construction analysis, a true-amplitude Kirchhoff-WKB J prestack depth migration equation is derived which appears to be new to the literature. This result is obtained in the form of a 2.5-D farfield Kirchhoff integral solution to the acoustic wave equation, after the application of a dynamic imaging condition to the reconstructed upgoing and downgoing wavefields. The solution is dependent upon WKBJ Green's functions, which can be numerically evaluated for heterogeneous migration models by raytracing methods, and is therefore fundamentally related to the eikonal and transport equation solutions of zeroth order asymptotic seismic ray theory (ART). A new and "generalized" Kirchhoff-WKBJ prestack depth migration equation is subsequendy obtained by the introduction of a weighting function into the true-amplitude migration integral. The weight is a function of both the reconstructed upgoing and downgoing wavefields, and is determined analytically by a mathematical application of each specific reflectivity imaging condition. This generalized equation is significant in that it provides a common mathematical, physical and computational basis for the comparative analytical and quantitative analysis of reflectivity image quality and amplitude accuracy among current prestack migration philosophies and variants of those migration themes. 204 In addition, three ancillary research objectives are achieved. Thefirstachievement is the development of a computer algorithm to implement the generalized Kirchhoff-WKBJ prestack depth migration imaging and amplitude-estimation of subsurface reflectivity given a set of observed seismic reflection data. The algorithm can be readily modified to perform seismic wavefield imaging for other recording geometries such as cross-borehole or vertical seismic profiling, and may be suitable to non-seismic applications such as the imaging of satellite-acquired synthetic aperture radar data. The second result is the development of a fast two-point raytracing computer algorithm which provides accurate computation of a subsurface grid of traveltimes and 1.5-D zeroth order ART amplitudes in a 1-D acoustic medium. This algorithm is useful for subsurface wavefield reconstruction and imaging, and for inversion applications such as geotomography. The third result is the detailed quantitative examination of migration imaging quality and true relative-amplitude normal-incidence acoustic reflectivity recovery from numerically migrated depth images. This is achieved successfully in an extensive 2.5-D synthetic data analysis, using a challenging 2-D structural model, synthetic multifold reflection seismogram shot gathers, and the numerical imaging and modelling algorithms developed as part of this thesis research. 205 References Berkhout, AJ., 1984, Seismic Migration: Imaging of Acoustic Energy by Wave Field Extrap- olation, I, Elsevier Science Publishing Company Inc., Amsterdam. Bleistein, N., 1986, "2.5-D in-plane wave propagation," Geophysical Prospecting, 34, 686703. Bregman, N.D., Bailey, R.C, and Chapman, C.H., 1989, "Crosshole seismic tomography," Geophysics, 54, 200-215. Carter, J. A., and Frazer, L. N., 1984, "Accomodating lateral velocity changes in Kirchhoff migration by means of Fermat's principle," Geophysics, 49, 46-53. Casde, R.J., 1982, "Wave-equation migration in the presence of lateral velocity variations," Geophysics, 47, 1001-1011. Cerveny , V., Molotkov, I.A., and PsenCfk , I., 1977, Ray Method in Seismology, Univerzita Karlova, Praha. Cerveny , V., and PSenCik , I., 1979, "Ray amplitudes of seismic body waves in laterally inhomogeneous media," Geophysical Journal of the Royal Astronomical Society, 57, 91- 106. Chang, W. F., and McMechan, G. A., 1986, "Reverse-time migration of offset vertical seismic profiling data using the excitation-time imaging condition," Geophysics, 51, 67-84. Claerbout, J. F., 1971, "Toward a unified theory of reflector mapping," Geophysics, 36, 467-481. Claerbout, J.F., 1976, Fundamentals of Geophysical Data Processing: with Applications to Petroleum Prospecting, McGraw-Hill. Clayton, R.W., and Stolt, R.H., 1981, "A Born-WKBJ inversion method for acoustic reflection data," Geophysics, 46, 1559-1567. Crase, E., Pica, A., and Tarantola, A., 1988, "Robust inversion of seismic waveforms," presented at the 58 International SEG meeting, Anaheim, California. th 206 Etgen, J.T., 1987, "Fast prestack depth migration using a Kirchhoff integral method," Stanford Exploration Project Report, SEP-51. French, W.S., 1975, "Computer migration of oblique seismic reflection profiles," Geophysics, 40, 961-980. Gardner, G.H.F., Gardner, L.W., and Gregory, A.R., 1974, "Formation velocity and density - the diagnostic basics for stratigraphic traps," Geophysics, 39, 770-780. Gazdag, J., 1978, "Wave equation migration with the phase-shift method," Geophysics, 43, 1342-1351. Hampson, D., 1989, "Using a blocky inversion approach to AVO analysis," presented at the Annual CSEG Conference, Calgary, Alberta. Keho, T.H., and Beydoun, W.B., 1988, "Paraxial ray Kirchhoff migration," Geophysics, 53, 1540-1546. Loewenthal, D., Lu, L., Roberson, R., and Sherwood, J., 1976, "The wave equation applied to migration," Geophysical Prospecting, 24, 380-399. Lumley, D.E., and Keys, R.G., 1988, "Fast and accurate two-point transmission times in onedimensional media," Mobil Research and Development Corporation Laboratory Mem- oranda, 112, LM881121-FR, Dallas Research Laboratory, Exploration Research and Technical Service, Project No. 501-846. Morse, P.M., and Feshbach, FL, 1953, Methods of Theoretical Physics, McGraw-Hill, New York. Newman, P., 1973, "Divergence effects in a layered earth," Geophysics, 38, 481-488. O'Brien, P.N.S., and Lucas, A.L., 1971, "Velocity dispersion of seismic waves," Geophysical Prospecting, 19, 1-26. Oldenburg, D.W., Levy, S., and Stinson, K.J., 1986, "Inversion of band-limited reflection seismograms: Theory and practice," Proceedings of the IEEE, 74, 487-497. Rocca, F., 1987, "Synthetic aperture radar: A new application for wave equation techniques," Stanford Exploration Project Report, SEP-56. Schneider, W., 1978, "Integral formulation for migration in two and three dimensions," Geophysics, 43, 49-76. 207 Sheriff, R.E., and Geldart, L.P., 1982, Exploration Seismology, I and II, Cambridge University Press. Stolt, R.H., 1978, "Migration by Fourier transform," Geophysics, 43, 23-48. Tarantola, A., 1986, "A strategy for nonlinear elastic inversion of seismic reflection data," Geophysics, 51, 1893-1903. Tarantola, A., 1987, Inverse Problem Theory, Elsevier Science Publishing Company, Inc., Amsterdam. Waters, K.H., 1981, Reflection Seismology - A Tool for Energy Resource Exploration, Wiley- Interscience, New York. Zauderer, E., 1983, Partial Differential Equations of Applied Mathematics, Wiley-Intersci- ence, New York. Zelt, C.A., and Ellis, R.M., 1988, "Practical and efficient raytracing in two-dimensional media for rapid traveltime and amplitude forward modelling," Canadian Journal of Exploration Geophysics , 24, 16-31. Zhdanov, M.S., and Frenkel, M.A., 1983, "Electromagnetic migration," in The Development of the Deep Geo-Electric Model of the Baltic Shield, Part II. Proceedings of the I Project Symposium, Oulu, 15.-18.11.1983, (Hjelt, S.E., ed.), Department of Geophysics, st University of Oulu, Report No. 8. 208
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A generalized Kirchhoff-WKBJ depth migration theory...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A generalized Kirchhoff-WKBJ depth migration theory for multi-offset seismic reflection data : reflectivity… Lumley, David Edward 1989
pdf
Page Metadata
Item Metadata
Title | A generalized Kirchhoff-WKBJ depth migration theory for multi-offset seismic reflection data : reflectivity model construction by wavefield imaging and amplitude estimation |
Creator |
Lumley, David Edward |
Publisher | University of British Columbia |
Date Issued | 1989 |
Description | This thesis embodies a mathematical, physical, and quantitative investigation into the imaging and amplitude estimation of subsurface earth reflectivity structure within the framework of pre stack wave-equation depth migration of multi-offset seismic reflection data. Analysis is performed on five prestack depth migration reflectivity "imaging conditions" with respect to image quality and quantitative accuracy of recovered reflectivity amplitudes. A new computationally efficient and stable prestack depth migration imaging method is proposed which is based upon a geometric approximation to the theoretically correct, but unstable, "dynamic" imaging condition. The "geometric" imaging condition has the desirable property of true-amplitude reflectivity recovery in regions of both 1-D and 2-D velocity variation, while fully retaining and optimizing the favorable imaging characteristics of current, non-true amplitude formulations. The currently predominant "crosscorrelation" and "excitation-time" migration imaging methods are shown to possess significantly less accurate imaging and amplitude-recovery characteristics relative to the proposed geometric migration. The respective signal-to-noise recovery of their imaged amplitudes deteriorates approximately linearly (excitation-time) and quadratically (crosscorrelation) with depth. As a necessary prerequisite to the imaging analysis, a true-amplitude prestack depth migration equation is derived which appears to be new to the literature. This result is obtained in the form of a 2.5-D farfield Kirchhoff integral solution to the acoustic wave equation, after the application of a dynamic imaging condition to the reconstructed upgoing and downgoing wavefields. This solution is in harmony with zeroth order asymptotic ray theory (ART) assumptions, and depends upon WKBJ Green's functions which can be numerically evaluated for arbitrary migration models by raytracing methods. A new and "generalized" Kirchhoff prestack depth migration equation is subsequently obtained by the introduction of a weighting function into the true-amplitude migration integral. The weight is a function of both the reconstructed upgoing and downgoing wavefields, and is determined analytically by a mathematical application of each specific reflectivity imaging condition. This generalized equation is significant in that it provides a common mathematical, physical and computational basis for the comparative analytical and quantitative analysis of reflectivity image quality and amplitude recovery among current prestack migration philosophies and variants of those migration themes. In addition, three ancillary research objectives are achieved. The first achievement is the development of a Kirchhoff prestack depth migration computer algorithm to implement the generalized imaging of surface-recorded seismic reflection data. This algorithm can be readily modified to perform seismic wavefield imaging for other recording geometries such as cross-borehole or vertical seismic profiling, and may be suitable to non-seismic applications such as the imaging of electromagnetic wavefields and satellite-acquired synthetic aperture radar data. The second result is the development of a fast two-point raytracing computer algorithm which provides accurate computation of a subsurface grid of traveltimes and 1.5-D zeroth order ART amplitudes in a 1-D acoustic medium. This algorithm is useful for subsurface wavefield reconstruction and imaging, and for inversion applications such as geotomography. The third objective is the detailed quantitative examination of migration imaging quality and true relative-amplitude normal-incidence reflectivity recovery from numerically migrated depth images. This is achieved successfully in an extensive 2.5-D synthetic data analysis, using a challenging 2-D structural model, synthetic multifold reflection seismogram shot gathers, and the numerical imaging and modelling algorithms developed as part of this thesis research. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052629 |
URI | http://hdl.handle.net/2429/27588 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1989_A6_7 L85.pdf [ 18.58MB ]
- Metadata
- JSON: 831-1.0052629.json
- JSON-LD: 831-1.0052629-ld.json
- RDF/XML (Pretty): 831-1.0052629-rdf.xml
- RDF/JSON: 831-1.0052629-rdf.json
- Turtle: 831-1.0052629-turtle.txt
- N-Triples: 831-1.0052629-rdf-ntriples.txt
- Original Record: 831-1.0052629-source.json
- Full Text
- 831-1.0052629-fulltext.txt
- Citation
- 831-1.0052629.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052629/manifest