UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A generalized Kirchhoff-WKBJ depth migration theory for multi-offset seismic reflection data : reflectivity… Lumley, David Edward 1989

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

Item Metadata

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

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 T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F G E O P H Y S I C S A N D A S T R O N O M Y We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A October, 1989 © David Edward Lumley, 1989 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada Date Oc -Tbpgg. \o t figcj DE-6 (2/88) 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, 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 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 ob-tained 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 philoso-phies 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, 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 7 2.3 The Scalar Wave Equation . 8 2.4 The Green's Function Problem 13 2.5 Physical Interpretation of the Green's Function Solutions 15 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 3 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 4 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 63 5.1 Discussion on "Velocity Models" 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 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 7 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 CPU Times for Green's Function Raytracing Evaluation 139 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-Incid-ence 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 172 xii 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 Correc-tion Factors Fn Calculated from the True Arc Test Model 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 193 xiii 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 xs = 1200 m along the Arc Test Model 121 xv 8.3 The Pre-Processed Shot Gather at xs - 1200 m along the Arc Test Model.. . 123 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 Imag-ing Times t(x0) 129 9.2 A 2-D Subsurface Map of the Raytraced Green's Function Geometric Diver-gence 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 Ampli-tude 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 Imag-ing 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 Imag-ing 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 Imag-ing Condition 155 10.9 A Single Shot Gather Depth Migration using the Kinematic Reflectivity Imag-ing 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 Re-flectivity Imaging Condition 162 xvii 10.13Migration Depth Trace Stability Analysis of the Dynamic Reflectivity Imag-ing Condition for Variable Stabilization Parameter Values e2 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 Excitation-Time Reflectivity Imaging Condition 178 10.17Full-Fold Kirchhoff-WKBJ Prestack Depth Migration using the Crosscorre-lation Reflectivity Imaging Condition 183 10.18The Muted Shot Gather at xs = 1200 m with 10% Additive Noise 187 10.19 Geometric Migration of a Noisy Shot Gather 188 10.20Normalized Depth Traces Corresponding to the True Normal-Incidence Ac-oustic 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 Ac-oustic 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 200 xix 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 ex-ploration 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 tech-niques, 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 iter-ative 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" be-tween 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 itera-tive 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 the final component of the inversion process, where the preliminary reflectivity and subse-quent 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 re-covery 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 con-struction. Migration has proved to be valuable for qualitative analysis of seismic reflection data, especially in the area of seismic structural interpretation. However, my present re-search is concerned with the quantitative, rather than qualitative, analysis and inversion of seismic data. Hence, this thesis addresses the more complete reflectivity model construc-tion problem, by developing methods which perform simultaneous reflectivity imaging and amplitude-estimation, within the context of seismic wavefield reconstruction. On this basis, five distinct reflectivity imaging conditions are examined in a comprehensive mathematical, physical, and quantitative analysis within a common Kirchhoff-WKBJ space-time integral solution domain. The fundamentals of Kirchhoff post-stack migration theory are first examined, in order to introduce the mathematical constructs and physical intuition necessary to ascend to the more complete Kirchhoff-WKBJ prestack migration theory. Fol-lowing 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 inte-gral. 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 quanti-tative analysis of reflectivity image quality and amplitude accuracy among current prestack migration philosophies and variants of those migration themes. The generalized Kirchhoff-WKBJ 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 al-gorithms developed as part of this thesis research. As a final remark, 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 wave-equation "migration" of first arrival 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 seismic first arrival waveforms. This may find application in the quantitative analysis of seismic refraction survey data, and also for VSP and cross-borehole 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 half space 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 Kirchhoff-WKBJ 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 wave-field as a function of the subsurface coordinates. The w-x migration method of Claerbout (1971) employs finite difference 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 post-stack 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 mea-surement 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 xm is defined as xm = (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 2 U(x, t) - ^dttU{% 0 = 0; x e V, (2.1) satisfying the spatial boundary conditions: U(x,t) = x¥(x,t) ; xedV, (2.2) lim ||L7(x,0||2 =0 ; x e dl/, (2.3) 9 and satisfying the temporal boundary conditions: U(x, t) = 0 ; t < 0, (2.4) / \\U(x,t)\\2dt < oo ; *e(0,«.). (2.5) Jo 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 rU. The enclosing surface d'V can be decomposed into 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 decrease in inverse proportion to the square of the radius, 'Sz1. In the limit as ^ approaches infinity the energy density flux approaches 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: dnU(x,t) = ¥N(x,t) ; xzdV, Cauchy: a(x)L7(x, t) + $(x)dnU(x, t) = Yc(x, t) ; x e 9V. The operator 3RA is the spatial derivative in the direction h, normal to the surface dV, as 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 U(x,t) = x¥(x,t) ; x e A, (2.7) where (^x, t) is, in general, the observed upgoing wavefield measured at the surface 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 up-the 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: 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 bound-ary conditions (2.3), (2.6), (2.7), for the primary problem (2.1): V2< G(x, t) - ^ 3«G(x, t) = -5(x - x')8(£ -1') ; x, x' e V, (2.8) c G(x, t) = 0 ; x e %, (2.9) 13 and lim ||G(x,*)||2 = 0 ; x e M. (2.10) 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 n = y/(x- x')2 + (y -y')2 + (z- z')2, (2.12) and r2 = \/(x - x')2 + (y -y')2 + (z + z')2. (2.13) We remark that G\ is the Green's function solution to the scalar wave operator in a whole-space, 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 In this case, z = 0 14 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 G i 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, G i 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. a finite time 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. The first arrival 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 r2 is longer than the direct travel path distance rls the reflected event arrives at x later than the direct arrival by a time difference of Ax = fo-rO/c). Naturally, the spherical divergence of the reflected arrival is larger than the direct arrival, and is expressed by the G2 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 surface ft, 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&'> *') L'G(^'' t''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 [-U{^> ' 'W* ' " W " ') " 0 ] 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&> *') L'G&> *> *) ~ G&> ?'> *> t) L'u(*'> ?) ] d*'> dt' = M [ U(x', t') dn>G(x', t'\ x, t) - G(x', t'; x, t) dn, U(x', t') ] dx' dt', (2.20) .V 19 where dni is the derivative, with respect to the primed spatial source coordinates x', in 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) dn,¥(x', t')} dx'dt', (2.21) Jt' J JA. 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', *') dn>G{x', t'; x, t) dx' dt' (2.22) Jt' J JA To proceed further, we need to evaluate the normal derivative of the Green's function, dn>G(x', t'). This evaluation requires an explicit definition of the Green's function G(x', t') in the new primed source coordinate system, which we shall now consider. 20 Recall the definition of G(x, t) from (2.11): G{%t\x',t') = Sjt-t-ri/c) _ §(t-t'-r2/c) 4nr{ 4nr2 Since G depends on x and x' only in the radial distance terms r\, r2, given by relations (2.12) and (2.13), G is reciprocal with respect to the spatial coordinates: 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 a time t, which is equal to the time of source ignition t', plus the finite travel 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 G(%t;x',t') = G(x',t;%t'), (2.23) 21 "arrow of time." Given the discussion above, we can now define the coordinate-transformed Green's function G(x', t') as Of!, f.%t) = S(''7 + r'/ C ) - «?-* + '#). ( 2. 24) 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, dniG(x', t'), using the definition (2.24), and evaluate the derivative on the surface A, as indicated by (2.22), we find that dn>G(x', Ok'e* = 23B,G1(x', OIl'eA (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')dn>Gi(x',t'-%t)dx'dt'. (2.27) Since h is the unit, outward directed, normal to the surface ft, 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')dz, { 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 derivative first, we find that 23 dzi = drdz<r = -cos0(x')oV, (2.30) 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) h'jt'-t + rlc) r2 + 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) After carefully transferring the delta-function derivative 8' to the observed wavefield XF, and then performing the indicated t' integration, (2.32) reduces to U(x,t) = ± y^cos0(x') T(x', t - r/c) dtVjx'.t-r/c) r2 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 equa-tion. 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 downgo-ing/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 re-flector 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 wave-field 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: 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; 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. (2.36) (2.37) 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 CMP-stacked seismic reflection trace data ^ (x, t) are summed spatially over the time trajectories t = 2r/c. In the 2-D (x, t) data plane these time trajectories 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 point-scatterers 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 point-diffractors, 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 diffrac-tion-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 equa-tion (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 wave-front 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 deriva-tive 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 a final implementational 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 well-suited 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 reflec-tion 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 ve-locity 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 post-stack processed data. Valuable offset and raypath information, destroyed in the stacking process, is vital to the reconstruction of reflectivity amplitudes, and therefore suggests the re-quirement for a prestack methodology. Furthermore, true-amplitude reflectivity is dependent upon the incorporation of realistic amplitude-attenuation recovery in the Kirchhoff wave-field 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 constant-velocity Green's functions. The WKBJ Green's functions are extremely flexible in the presence of 2-D migration velocity model heterogeneity, and can be numerically evaluated for rxaveltime and amplitude response by raytracing or finite difference 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 mi-gration 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 Kirchhoff-WKBJ 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. The first part 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(x0). First, we will examine the true-amplitude recovery of reflectivity from first principles. This examination emphasizes migration as an inversion of the observed seismic data in the con-struction 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 ob-tained by extending the conventional Kirchhoff integral theory with the use of WKBJ ray-theoretic Green's functions. Based on the section regarding first principles 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 = (xs,ys>Zs) is the (fixed) shot coordinate, and xr = (xr, yr, zr) is the (variable) receiver coordinate. 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 = (x0, yo, z0) is the subsurface model coordinate. Upon reflection, a portion of the downgoing energy is converted into upgoing waves, denoted by U(\0, fr, xs). When the upgoing waves impinge upon the recording surface, we record the resulting observed seismic data: V(xr, fr, Xs) = L/(Xo, fr, %)\^r. (3.1) The upgoing wavefield U, and hence 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 re-flectivity distribution i?(x0), which upon secondary excitation by the incident source field D(x0, t; xs), gives rise to the upgoing reflection wavefield L7(x0, fr, xs), subsequently mea-sured at the recording surface as the observed seismic reflection data (^x,., fr, xs). 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 reflectivity 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 from first principles. 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(x0,Q)D. (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 #(x0,e)->E(xo). In the spectral domain, (3.2) implies that R = U/D (3.3) 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(x0, 0, w). In practice, we will assume 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 denom-inator by the complex conjugate of the downgoing wavefield, D*, results in „ UD R = (3.4) DD* v ' We recognize the numerator of (3.4) as the spectral equivalent of a time-domain crosscorre-lation. 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 (3.5) 40 and 7(x 0 ;x 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 wavefield: | D |2. 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) and \D\2= D*D (3.8) Finally, in the space-time domain (3.4) becomes R = /(x 0 ; Eg) (3.9) 41 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 upgoing field, U(x0, t; xs), and the downgoing field , D(x0, t; xs), perform their zero-lag crosscorrelation, normalize by the energy of the downgoing field, and hence obtain the reflectivity model estimate R(x0). Note that the implementation of (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(x0) in Equation (3.9). 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. Using finite difference operators to approximate the wave equation, Claerbout reverse-time extrapolates the boundary wavefield in the space-time domain. This approach forms the basis of finite-difference wave-equation migration. A major advantage of this approach is that an arbitrary 2-D migration velocity model can be readily specified on the finite-difference grid. Another advantage is that full wavefield effects such as post-critical reflections, critical refractions, and multiples may be directly incorporated in the migration without pre-processing the recorded data to remove these effects. Furthermore, the finite-difference algorithm 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 Trans-form 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 Fourier-transformed into the spectral domain whereupon a phase-shift downward-continuation oper-ator 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: V2C/(x0, t; Xs) - \dttU(Xo, t; xs) = 0, (3.10) cl 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 parameter field, 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 / / / ¥(xr, t'; xs) dnGr(xr, t'\ x0, t)dxr dyr dt', (3.11) derived previously as Equation (2.27), where dn is the derivative normal to the recording 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 = ^ ' - ^ . * - « ' - ' * > ; ( 3 . 1 2 , 4jcri 4 7 t r2 the radial source and image distances are given by r\ = {Xr-Xof + iyr-y0)2 + (zr-z0f and r\ = (xr-x0f + (yr -y0f + (zr +z0f. If the velocity is a function of x0 and za, i.e. c(x0, zQ), then we can look for a WKBJ Green's function solution by assuming a frequency domain representation of the form: G(xr, w'; x 0) = A(x r ; x0)S(M;')e±l'M;'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 r . The phase term t(x r; x 0) represents the traveltime from a given subsurface location to the receiver locations . The function S(w') represents the temporal Fourier transform of a time-domain source function S(t'). The choice of sign in the complex exponential depends upon the time-direction of wave prop-agation. Since we are considering energy propagating from a particular subsurface location x 0 to all possible surface receiver locations x r , the choice of negative phase is correct for 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 the time domain, (3.13) inverse-transforms to 45 G(xr, t'\ x0) = A(xr; x0)S(t' - x(xr; x0)). (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 re-expressed as dnG = VGh (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 V G = S V A + A S T V T , (3.17) where Sz is the derivative of <S with respect to x. Since the Green's function is a 3-D point-source impulse response, its amplitude A decays approximately as 1/r, and hence the gradient of A decays approximately as 1/r2. Because of this the second term in (3.17) is dominant for large r. We neglect the first term in (3.17) and anticipate a farfield approximation to the full solution: 46 and VG~AS xVx VGn~AS TVx-n. (3.18) (3.19) 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.=C O S 0 ( X o ) 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(x0,t;%) = 2 III ^ ^ A C x r j X o J S r ^ ' - x ^ X o J - O ^ ^ X s ) ^ ^ ^ ' (3.23) Jt'JjA c( xr) where 6(x0) and c(x0) have been evaluated at the surface location xr as required by the surface 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(n)(x) q>(*)dx = (-1)" J 8(a) <?{n\x)dx, (3.24) to transfer the time derivative from &(t') to W. Also, noting from (3.15) that dfS = -dtS (3.25) and performing the indicated time integration in (3.23) results in r r c o s 0( x ) t/(x0, t; xs) = 2 / / —-i-IiA(x r ; x„) 3^(xr, t + x(xr; xG); xs) dxr dyr. (3.26) JJj? c(xr) 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 x0 48 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 spatial derivative d]l2/y/c. The term d]12 is known as the "half-time-derivative operator" and 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(xr)y a,1/2¥(xr, t + x(xr; x0); xs) dxr (3.27) in analogy with Equations (13) and (14) of Carter and Frazer (1984). 3.5 The Downgoing Wavefield Solution Now we calculate D(x0) t;xs) as required by the reflectivity inversion (condition 3.9). If we 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; xs)S(w)e-iw^0™ (3.28) 49 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(x0, t; Sg) = A(x0; xs)S(t - t(x0; xs)). (3.29) Thus, D(x0, t;xs) is the original source wavelet S(t) time-shifted by the traveltime T ( X 0 ; X 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 illumination field, hence we obtain a 2.5-D migration solution. Now we use the reflectivity imaging condition (3.6), and crosscorrelate the upgoing and downgoing fields retaining 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 7(x0; xs) = / A(x0; xs)S(* - T S ) • 2 cos Qr fA(xr;x0) c(xr) Bl^ixr, t + T r ; Sg) dxr • dt. (3.31) 50 Assuming again that the source function S(t) takes the form of a delta function, we perform the time integration in (3.31) resulting in 7(Xo;xs) = 2AS f cos9rJ — at1/2xP(xr, t;xs)|^s+Trdxr. (3.32) J y cy 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): | D(x0; Xs) |2= A2(x0; xs) = A2 (3.33) Substituting (3.32) and (3.33) into (3.9) we obtain the migrated reflectivity estimate for a fixed shot gather: Rbo,Zo\V) = .n?''*!1] ,2 = / ^ ^ ( * r , t;%)\t^rdxr. (3.34) I •L'l.Xo) x s ) I J y/Cr As 3.7 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 As is never known, which 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 signal-to-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(xmz0)= f f d]*V&r,ttZa)\tWrdxrdxa. (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 multi-offset migrated reflectivity estimate of (3.35) is effectively averaged over all individual shot-dependent 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 the final migrated 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 travel-times 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(x0, z0) is the migrated depth estimate of subsurface reflectivity, d]12 is the half-time-derivative operator, P^(xr, t, %)\t=%s+ir are the recorded shot-gather traces evaluated at the total shot-subsurface-receiver traveltime, xs = T ( X 0 , Xs) is the raytraced traveltime from xs to x0, %r - i(xr, x0) is the raytraced traveltime from x0 to xr, As — A(x0, xs) is the 3-D amphtude attenuation from xs to x0, Ar = A(xr, x0) is the 3-D amplitude attenuation from x0 to xr, cr = c(xr) is the velocity evaluated at the receiver location xr, and 0 r = 0(xr) is the ray incident angle at xr with respect to the surface normal for a ray traveling from x0 to xr. Finally, an implementational remark concerning the Kirchhoff-WKBJ prestack depth migra-tion theory just derived is in order. No special pre-migration data processing is required other than standard trace-editing, bandpass filtering and 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 sub-surface 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 con-dition, it is first necessary 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{x0,z0)= I f ^ ^ ^ y ^ dlaV(xr,t,xs)\t^rdxrdxs. (4.1) J JA VCr A s Notice that the generalized solution (4.1) is equivalent to the true-amplitude reflectivity solu-tion (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 • Rd(x0) = | [U(x0, t) *Z)(xo, t)] I [D(x0, t) * D(x0, t)] |,(=0. • Wd = 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: Rd = U/D, (Claerbout, 1971). This imaging condition corresponds to 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 out-of-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 inver-sion for the absolute-amplitude reflectivity distribution values rather than a simple image of reflector locations. In practice, the dynamic imaging condition may produce unstable reflec-tivity estimates in regions of weak source-illumination energy, i.e., where A s —> 0. In the 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: A s -> A s + e2. The parameter e2 is generally taken to be constant with its magnitude being a small fraction of the average value of the source amplitude A s . For an optimal stabilization, the parameter should be set to the estimated background-noise amplitude-level, which in general is a func-tion 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(x0,t)-kD(xo,t)\t>^o. • <WC = A 2 . Claerbout (1971) overcomes the dynamic imaging condition instability resulting from the spectral division U/D, by defining a crosscorrelation imaging condition: Rc = 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: Wc = A2S, in (4.1). Claerbout uses finite-difference operators 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 Ar amplitudes must be determined by external means such as raytracing. 4.5 The Excitation-Time Imaging Condition • -Re(Xo) = L7(X0,£)|*=IS-• We = A8. 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: Re = 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(x0; x s). This assumption would be valid for a delta-function source of unit 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 one-way upgoing propagation at half the true medium velocity. The ZSR imaging times are correct but only the Ar amplitude divergence is compensated along the upgoing raypath; the As term along the downgoing raypath is effectively ignored. The excitation-time imaging condition migration can be implemented by choosing the integrand weight of (4.1) equal to the source amplitude component: *M4 = As. 4.6 The Geometric Imaging Condition U\x0, t) * D\x0, t)] I \D\X0, t) • D\x0, t) «'=0 Wg = Asrsly/ArTr A stable and efficient geometric approximation to the dynamic imaging condition is now proposed: Rg = ff/D\ The terms if and & are identical to U and D except that the dynamic amplitude factors Ar and As are replaced with their constant-velocity divergence 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 atten-uation. 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 = Asrsl\/Arrr 60 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 rs/y/f^. These bounds constrain the ampli-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 rr and rs. Furthermore, for an arbitrary 2-D migration model, the geometric amplitude factors need be calculated only once, whereas the other imaging conditions require extensive raytracing of the dynamic amplitude maps Ar and As for each surface (shot and receiver) location. 4.7 The Kinematic Imaging Condition • Rk(Xo) = • ( X r , * ) | * w In the extreme case that we constrain the variation of the amplitude ratio rs/y/ry such that it is constant for all subsurface locations, then the geometric imaging condition reduces to the kinematic imaging condition: = ^ (x,., t = xs + xr). The kinematic condition obtains 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 *(xr,o = ^ a^r*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 trace amplitudes), to perform the migration. Setting rM4 = As/y/A^ 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 informa-tion, 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 re-quired, 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 reflec-tivity 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 imag-ing 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 under-stand 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 dis-cussion 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, ctrue, into a reference velocity component c ,, and a residual velocity component Ac, such that Ctrue(*o) = Cre/r(XD) + AC(X0). (5.1) Usually, the reference velocity c is specified to incorporate the major features of the true velocity field ctrue 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 ctrue, and the residual field Ac represents a second-order variation such that: Ctrue ~  Cref and \te\<\cnf\. (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 velocity field distribution 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 f. In strongly scattering media, the reference velocity is a significant factor in wave 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 cmig in order to accurately reconstruct the incident and reflected wavefields. Ideally, we should define cmig ~ ctrue- m practice, we do not know the true velocity model ctrue a priori, and so usually estimate the reference velocity component cnf of the true velocity field such that 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 cmig can be expected to provide an excellent migration result, provided that cmig is a reasonable estimate of cnf as implied by (5.4). In this case, the migration wavefield reconstruction integrals can be implemented accurately with cmig since the true wavefield propagation is governed primarily by cnf, and not 8?. Hence, 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 and Ac. Implementa-tion of the migration wavefield reconstruction integrals with cmig ~ cnf will deteriorate the 65 migrated result as the approximation cmig ~ ctrue 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 cmig. However, this may be difficult to do in an a priori fashion. A residual velocity estimation scheme may be implemented by an iterative migration method, in which the residual velocity Ac and migration velocity model cmig 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 it first. Let us begin with the geometric version of the generalized migration solution (4.1), 66 Rg(xo) = 2 JJ?^^M{xr,%\t)dxrdxa, (5.5) where M(xr, xs; t) are the recorded surface data modified by the half-time-derivative operator and evaluated at the imaging times: M(xr, xs; t) = 3,I/2T(xr, xs, t) |<=Ts+Tr. (5.6) The correct dynamic version of (4.1) is given similarly by: Rd(x0) = 2 JJ^ °-?jt ^ M(xr, xs; t) dxr dxs. (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 Ar ~ l/rr and As ~ l/rs. 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, (5.8) 67 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 (straight-ray) 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 subsur-face 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 nth layer of a 1-D medium composed of N constant-velocity horizontal plane layers: 68 ( n<N \ /n<N \ XJWcij i^hckcos2QJcrcos2QkJ, (5.9) where ^  is the length of the ray in the kth layer, is the velocity in the kth layer, and 0^  is the angle between the ray and the vertical in the kth layer. If we assume that, to first order, 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 AD\zn) ~ £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) = cmig, and remarking that the constant factor C\ reduces to the velocity cs at the source region, then (5.10) reduces to a continuous relative-amplitude 1-D geometric divergence factor: i fl(Zo) AD\z0) ~ - / c(z)dl. (5.11) cs Jl(zs)=0 The integration (5.11) is performed along the raypath length I from the source depth zs to the subsurface point z0. In a constant velocity medium, the raypath is a straight line of total length l = r, hence (5.11) reduces to AQ = 1/r as expected. However, we consider a 1-D subsurface velocity variation. Continuing the first-order analysis, 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 c(z) ~ Co + az. (5.12) 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 with the previous assumption of predominantly vertically-propagating energy. In this case, (5.11) becomes A-J(Zo) ~ - fZ°c(z)dz. (5.13) Cs Jzs Substituting (5.12), performing the z-integration, and noting that cs = c0, reduces (5.13) to 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 the first term of (5.14) is of leading order: Hence we have shown that, in a medium with a reasonable vertical velocity gradient, the dr = 2zdz/Vx2 + z2 = IcosQdz -> dz, (5.14) AD ~ 1/r (5.15) 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, to first order, 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 Recall that M(xr, xs; t) are the recorded surface data modified by the half-time-derivative operator and evaluated at the imaging times, as given by relation (5.6). of (4.1): (5.16) 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(xr, xs; t) do not depend (much) on the shot location xs, but rather on the shot-receiver offset within a given shot gather. Hence, the independent variable xs becomes a parameter such that M(xr, xs; t) -» M(\r; xs, t). Now M(xr; xs, t) can be thought of as a 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 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 (5.17) Now consider the dynamic imaging-condition migration equation: (5.18) 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: 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 s (x 0 , xs). The shot-point integration sums over a line of closely spaced point-sources, effectively creating a line-source response. Hence, the integration sums the many l/r(x0, xs)-type point-source divergence responses to a single l/y/r(x0; xs)-type cylindrical line-source divergence. The line-source integration converts the shot-variable WKBJ point-source amplitude A s (x 0 , x s), to the shot-independent WKBJ line-source amplitude factor \ZA s(x 0; xs). Hence, Equation (5.20) integrates to 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 rs of the total reflection path is approximately equal to the receiver-to-subsurface reflection distance rr. This in turn makes Ar approximately equal to A s , and the WKBJ amplitude ratio in (5.21) reduces to unity; (5.20) (5.21) As~ 1/r, s> Ar~ 1/r, r> 73 rs ^ rr, 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, 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 (5.23) 74 migration, as later demonstrated. If reflected energy originates from steeply-dipping 2-D reflectivity structure, then the source and receiver amplitude factors As and Ar will be 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^lAs, 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 a finite dif-ference 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 summa-tion 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 *We = As: -Re(x0) = 2 cos 6r Ar M(xr, xs; t) dxr dx. (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: 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 where the As factor has been brought out from underneath the integral for clarity, and does 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: (5.25) (5.26) 76 Re(xQ) ~ \/As(x0;xs)Rd(x0). (5.27) In previous sections we have shown that the WKBJ amplitude A can be usefully approximated by a 1/r-type attenuation factor: As ~ l/rs. In this case, (5.27) can be interpreted as Re(xo) ~ ] _Rd(x0). (5.28) vrs(,x0; xs) This clearly illustrates the nature of the reflectivity amplitude error inherent in the excitation-time 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/rs(x0) weighting matrix to each 77 migrated shot gather before stacking into the final composite estimate of Re(*o)- Also, if the downgoing source energy can be shown to propagate nearly vertically, as is the case for near offset or very deep reflection surveys, then rs ~ (z0-zs) and a partial 1-D correction can be 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 a first-order approximation 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 Claer-bout (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 migra-tion, except an additional amplitude attenuation is imposed by the crosscorrelation source amplitude term As. For this reason, we expect crosscorrelation migration to have an ampli-tude recovery error that can also be expressed as a function of the WKBJ source amplitude As, but even more severe than the excitation-time imaging method. The crosscorrelation implementation of the generalized Kirchhoff migration equation (4.1) can be obtained by setting WC=A2.: -Rc(xo) = / / 2 c ° ^ r A s \fA~rM(xr, x s;t)dx rdx s, (5.29) JJA y/Cr 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 Rc(x0) = J y/Ar • J As M(x r , x s; t) dxs • dxr. (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 s , and can be brought out from underneath the shot integral: /C O S Q — r —j==- y/ArM(xr;xs, *)• JAsdxs-dxr. (5.31) 79 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, As —» >/A7, and (5.31) reduces to Rc(x0) ~ y/As / — ~ y/ArM(xr-,xs,t)dxr. (5.32) Comparing (5.32) with the dynamic migration equation (5.26) developed for the same ve-locity assumptions, we obtain a relation quantifying the reflectivity amplitude recovery error inherent in the crosscorrelation imaging method: Rc(x0) ~ As(x0;xs)Rd(x0). (5.33) If we again make the correspondence of As ~ l/rs, (5.33) can be interpreted as Rcfro) i?d(x0). (5.34) Relation (5.34) demonstrates the severe attenuation in reflectivity amplitude recovery ex-pected 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~s. Relation (5.34) suggests amplitude correction schemes similar 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 mi-gration 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 2-D 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 / AJ. 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 es-timates more than the limiting assumptions of the geometric imaging condition. Hence, in terms of practical implementation, and certainly efficiency, the geometric imaging condi-tion 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 WKBJ 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 pre-vious 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, V2C/(x, t) - ^-a«L7(x, t) = 0. (6.1) cz(x) 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 = Fk-dO, (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 (-ico)A (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 (ico)exp-iQ,(^) (ico) £(i^( v 2"< t o ) < v< ) 2-^) + 2 V y u k + 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 _ C2(X) ' (6.7) and the Transport Equation, A(x)V2x(x) + 2Vx(x)-VA(x) = 0, (6.8) where A - U^. 85 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), U(x, t) = A(x) exp-i<0(^(S)), (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 Kirchhoff-WKBJ 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 a first order nonlinear partial differential equation (PDE): I VT(X) I2 = 1 c2(x) It can be converted to a system of first order linear ordinary differential equations (ODE's) us-ing 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 velocity field, c = c(x, z), the raytracing system is given by three ODE's (Cerveny et al, Equation (3.16)): dx . n — = c sin 9 dx dz — = c COS0 dx dQ — = -cx cos9 + cz sin9, (6.10) dx where cx = dxc(x, z), etc., and 9 is the angle the ray makes with respect to the vertical axis k. For a heterogeneous 1-D velocity field, c = c(z), two raytracing ODE's must be solved: dx 2 Tx = c p 87 and ^ = ±cy/l -p2c2, (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 cp dx dz dz A/1 -p2c2 and dxdx _ dx _ 1 dx dz dz c A/1 -p2c2 Integration of Equations (6.13) results in and (6.13) fz pc(z')dz' x(z)= / p K [ + *(zs) (6.14) ./*• V1-P2C2(2') = / + x(zs) (6.15) 7* c(*') y/1 -P2C2(Z>) 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 con-siderations. 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 xs and xr for a 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)V2x(x) + 2Vx(x)- VA(x) = 0. 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 a fixed ray traveltime %r, x{y\, Y2, tr) sweeps out a 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 = (6.18) c J 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 cross-sectional area da of an elementary ray tube by Equation (3.23) of Cerveny et al. (1977): da = \J\dyidy2. Re-expressing the o\x term of (6.17), (6.20) 9Tx I = dx dx dx dx = 1 vx r1 = c, the Laplacian becomes (6.21) 1 d (J Vzx = — — , , . Jc 9T \c J Also, the transport equation dot-product can be evaluated in ray-centred coordinates: (6.22) VA Vx - — — - ^ i ^ L ^ l - ^ (Hl\2 - 1 — dxi dxi dx dxi dxi dx \dxi J c2 dx' Combining (6.17)-(6.23), the transport equation in ray-centred coordinates becomes (6.23) 2d,A + A_d fj_ c2 Jc T \c 0. (6.24) 91 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), (6.25) A 2 (J/c) and integrating both sides of (6.25) results in In A = - i In (- } + do. (6.26) 2 \c Exponentiating both sides results in a solution of the transport equation for A, " Ao]fj' (6.27) where A0 is the amplitude of the initial source, AD = 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 = - ^ = , (6.28) y/jpb" 92 where \)/0 is a source amplitude function, p(x) is the material density function, and c(x) may 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 \(/0, it is necessary to examine the wave solution in a small spherical homogeneous volume centered about the source location xs. Cerveny et al. (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 s p s s i n 0 s , (6.29) and line-source: Vo = #(9S) VeTpI, (6.30) where g is a function describing the directional characteristics of the source, and the s subscript indicates evaluation at the source location xs. For the seismic migration application 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 the first factor 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 J = (6.32) 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 approxi-mated by ray differencing. For a given ray, a ray-tube is created by raytracing adjacent rays for very small increments dQs of the ray take-off angle 0S. The change in offset or range dx 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, finite difference 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 _ sin 0(z) P ~ ~cW' Rearranging the required partial derivative we can evaluating the first term from the ray equation (6.14) for x(z), where cos 0 is defined by the ray parameter as c(z')dz'  Z s cos3 0(z')' cos Also, we can make the evaluation dp\ cos 0S Combining (6.36)-(6.39), we obtain an analytic expression for the JM derivative, (6.36) (6.37) 0 = Vl-sm 20 = V l -P2c2- (6.38) (6.39) 96 \dQsJz cs JZs cos3 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 jD, _ cos8,cos8(*) f#l_„ (6.41) 11 cs J Z s cos39(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, flfc{l)dl ( 6 4 2 ) cs 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) . (6.43) For the 1-D medium, Jj. is given simply from (6.42) as Cs Jzs COs6(z') 97 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 0S cos 0S cos 0(z) cos 6(z') c(z') dz' (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 Zoep-pritz reflection-transmission amplitude coefficients and matching impedance jumps from the divergence amplitude A in (6.31) at a boundary: where the impedance terms £ = pc are defined as: £s at the ray starting point (source), C,r at the ray endpoint (receiver), at the point of incidence at the ith layer, and ^ at the point of emergence at the ith layer. Also, n is the number of boundaries encountered by the ray, and Zi is the Zoeppritz displacement amplitude coefficient at the ith boundary. (6.46) 98 Unfortunately, the amplitude divergence as written now also exhibits a non-physical ampli-tude discontinuity across layer boundaries. Cerveny et al. (1977) add a compensation factor so that the resulting amplitude expression becomes A = ADAP, (6.47) where AD is the geometrical divergence amplitude factor given by and Ap is the energy partition amplitude factor given by ^ • t e ) " n ( g ) % In general, the Zoeppritz coefficients Zi can accommodate reflection, transmission and con-version 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, J I . = 2p£Cf/cosef  1 pi+1ci+l/cosQi+1 + pjCi/cosGj" 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 where f0 is the dominant frequency of the propagating pulse, U is the length of the ray segment, c; is the average velocity, and Q{- is the average Q value along the raypath within the ith 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 fQ. Fortunately, this approximation is often very good for seismic reflection data. 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: 71 AQ = Y[ exp(-7t/0Zj/ciQj), (6.51) i=i A = ADAPAQ. (6.52) 100 The full solution A has an amplitude component AD representing geometric wavefront di-vergence 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 compo-nent Ap arising at layer boundary interfaces, defined by (6.49) as where, in general, the Z; are the Zoeppritz elastic reflection-transmission amplitude co-efficients. Furthermore, an approximate form of Q-attenuation is accommodated by the amplitude component AQ, defined by (6.51) as n AQ = 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 un-limited number of horizontal layers, each specified by a depth, constant P-wave velocity, density, and Q value. The algorithm computes direct transmission times, geometric diver-gence, 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 and final locations, 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 (xs,zs) to any arbitrary depth point (x*,z*) in a 1-D medium, as shown in Figure 7.1. Since the velocity model is 1-D, the traveltime is a function of the source-receiver offset, | xs-x* |, and the depth of the subsurface point relative to the shot depth, | z8-z* |. We may set xs, zs = 0 and use x* and z* as the relative offset and depth without loss of generality. 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 (xs, zs) 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 spaced finely enough, 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> zs) and arrives at a specified subsurface location (x*, z*), such that the appropriate ray equations are satisfied. 105 Using the Newton-Raphson method, we solve a nonlinear integral equation for the ray pa-rameter p* given the initial and final boundary conditions for the ray, (xs, zs) and (x*,z*). The well-known offset and traveltime parametric ray equations for a 1-D medium are given by (6.14) and (6.15) as M = L-JI p c{z') dz' Ls ^\-p2C2{z>) and ( 7 . 1 ) fz dz' Jzs C(Z') y/\ -P2C2(Z') 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*2c\z>) ( 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 (xs,zs) to (x*,z*) and whose take-off angle Gs is given by: 106 9* = arcsin(csp*), (7.5) where cs is the velocity at the source location. We can also calculate the incident angle at the ray endpoint by substituting c* = c(x*, z*) for cs in (7.5), 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 (xs,zs) to any given subsurface point (x*,z*). 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, 9* = arcsin(c*p*), (7.6) p* = p0+ &>, (7.7) where p0 is an initial estimate, or trial solution, that is "close" in some sense to p*. If p0 is reasonably close to p*, then the correction term Sp is on the order of a perturbation such 107 that dp < p*,Po. (7.8) We can expand x(p*) of (7.4) in a Taylor Series Expansion: x(p*) = x(p0) + x'(p0)$p + x"(p0)bp2/2\ + •••. (7.9) If condition (7.8) is valid, then to first order, x(p*) ~ x(p0) +x'(p0)8p, (7.10) and the ray parameter correction Sp to the initial estimate p0 is found to be [x(p*)-x(p0)] op — — . (7.11) X'(p0) 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 p0. We can halt the iterative process when pQ is arbitrarily close to p*, in the sense that \\x(p*)-x(p0)\\ < e2, (7.12) 108 where e2 is a user specified error criterion. Obviously, the convergence of such an iterative solution depends upon the initial ray parameter estimate p0. 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(p0) is given by the evaluation of (7.1) using the subsurface gridpoint depth coordinate z* and the initial ray parameter estimate pQ. We also need to calculate the derivative x'(p), which is defined analytically by (6.37) as, where cos 9 is determined from the ray parameter definition (7.3) and is given by (6.38) as 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 1-D model parameterizations which consist of a sequence of N horizontal, constant-velocity (7.13) cosG = Vl - sin2Q = A/1 -p2c2. (7.14) 109 layers. In this case, the ray equation (7.1) for offset becomes, 71-1 x(p) = tanGj + ( z - z B ) t a n e „ , (7.15) i=i the offset derivative equation (7.13) becomes, X'ip) = > s- r - + T"7\ > (7-16) ^ — ' COS3 0; COS 0ra f=l and the ray traveltime equation (7.2) becomes, T(P) = V — + , (7.17) ^—' cj, cos 0t- c„ COS 0ra f=l where n < N is the number of the layer in which the subsurface gridpoint (x, z) is located, hi is the thickness of the ith layer, c; is the velocity in the ith layer, zn is the depth to the top of the nth layer, and 0j is the angle that the ray makes with respect to the vertical in the ith 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 G r a p h of X(p) t R a y P a r a m e t e r Figure 7.2: A Graphical Representation of the Newton-Raphson Method for Finding the Ray Parameter p*. I l l 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/cwaje. 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/cmax, 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 pa, the updated estimate p\ is found at the intersection of the horizontal line x = x*, with the line that passes through the tangent x'(p0). From the graph in Figure 7.2, it can be seen that successive iterations converge to the ray parameter p* when the initial estimate p0 is selected such that x(p0) > x*. It is also apparent that the region of convergence decreases as p* —> Vcmax. Physically, this implies that more iterations are required to find a 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 Newton-Raphson two-point raytracing algorithm sequence of Newton iterates, Pn+l = Pn+ [X*~X(Pn)] ; 71 = 0,1,2,..., (7.19) X'(pn) always converges to the unique ray parameter solution p* which satisfies x* = x(p*), provided that the initial ray parameter estimate pQ is chosen such that, 0 < p0 < 1/Cmax, (7.20) 112 and x(p0) > x*. (7.21) An initial estimate which always satisfies this condition is the starting value, (x Xg) (7.22) Po = Cmax y/(x* -Xsf+H1' where cmax is the maximum layer velocity, and H is the thickness of the layer having that velocity. In most applications, the initial ray parameter estimate p0 must be determined for an entire 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 pQ, the faster the convergence to the correct solution p*. Therefore, the most efficient procedure for estimating the initial ray parameter values p0, is to use (7.22) for the 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 = ADAPAQ. (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 constant-velocity 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 coeffi-cients. 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 equa-tion (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, p i + 1 C j + , / c O S 0 i + 1 + p;C; /cOS0; ' Furthermore, an approximate form of Q-attenuation is accommodated by the amplitude com-ponent AQ, defined by (6.51) as AQ(p) = Yl exp(-7t/oAj/QjCjCos0j). (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 ampli-tude 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 the five imaging 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 sedi-mentary 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 Layer .Base De/Jf/z (km) Velocity (km/s) Density (g/cc) Arc Dip 1 0.600 1.500 1.929 82° 2 0.900 2.300 2.147 69° 3 1.200 3.500 2.384 60° 4 1.500 4.500 2.539 50° 5 1.800 5.500 2.670 38° 6 2.100 2.500 2.192 16° 7 3.000 6.000 2.728 -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 (x0, zQ) 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. Subse-quent 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. The five prominent 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 xs = 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. 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 the first four 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 xs = 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. 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 Layer Ztase De/tf/j (km) Velocity (km/s) Density (glee) 1 0.600 1.500 1.929 2 0.900 2.300 2.147 3 1.200 3.500 2.384 4 1.500 4.500 2.539 5 1.800 5.500 2.670 6 2.100 2.500 2.192 7 3.000 6.000 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 raytrac-ing 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 ac-commodates 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 w0. This 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) = Q0c2{z), (8.2) where the constant QQ may be specified as 18.83 for velocity units of km/s. 126 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' = %r + xs, implement amplitude recovery with the dynamic amplitude weights Ar and As, and successively map each portion of the observed seismic data as a contribution to the total estimate of the subsurface reflectivity function R(x0). As such, it is evident that accuracy in the evaluation of the Green's function 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 a finite difference 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, Q-attenuation, 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 a finite difference 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 constant time and hence delineate direct-transmission wavefronts. It is evident that the wavefronts are circular within the first constant-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 discon-tinuosly (red-yellow) across the boundary; incident wave energy cannot refract through the layer boundary at post-critical incidence. This phenomenon results in transmission-time dis-continuities 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. How-ever, 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)(x0) due to an isotropic point source excitation in the 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 the first layer 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 0). Each color pixel represents the geometric divergence Aj ; (x 0 ;x s ) from the unit-amplitude source position x s (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. 1 3 1 This is the familiar 1/r geometric-spreading term in a constant-velocity medium. For subse-quent 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 post-critical 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 par-titioning 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 mi-gration 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, lo-cated 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 r(x 0). Each color pixel displays the obliquity amplitude factor cos 6 r (x 0 ;x s ) , detenruned by the surface-incident angle of the direct transmission ray joining the source position x s , located at the origin (dark red), to the pixel subsurface coordinate location Xo. The color legend specifies obliquity amplitude on a linear scale. 134 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-with-depth 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 am-plitude 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 func-tion 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 ampli-tude 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 com-ponents As andAr are both obtainable, as a function of offset, from the single general WKBJ 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 As and Ar maps must be obtained separately, 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 0 ) . Each color pixel displays the total W K B J Green's function amplitude factor Ai(x0; 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. 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 0 ) cos 8 r(x 0). 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 r(x o; x r ) of Figure 9.4. The color legend specifies amplitude attenuation on a logarithmic (dB) scale. 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 nu-merical evaluation of the WKBJ Green's functions over a 2-D finite difference 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 (NX,NZ) gridded map of raytraced values, where Nx = Nz = 101. Hence, in effect, a ray is traced from a single shot location, Ns = 1, to each of Nx • Nz = 10,201 subsurface "receiver" locations. In general, the computation time required for the raytracing component of the mi-gration is proportional to Nm • Nx • Nz, where Nm is the number of raytraced maps required 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 Raytraced Quantity CPU Time (s) Transmission Times 19.3 Geometric Divergence AD 16.6 Energy Partitioning AP 24.2 Q Attenuation 24.5 Surface Obliquity cos 0r 15.7 All Amplitudes At 25.2 All Amplitudes and Times Ai, Xi 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, Nm = 1, is required to implement the generalized Kirchhoff-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 the five imaging 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 syn-thetic 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, first by exam-ining the single trace migrations, then by summing reflectivity contributions from the receiver locations xr to obtain a single migrated shot gather at xs, and then finally by 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 Reflec-tivity 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.lAs, which is approximated to first order by the geometric amplitude factor, rs/y/ry, is relatively small near the source location (rs < ry), and relatively large near the receiver location (r s > rr). Relatively moderate migration ellipse amplitudes occur near the central midpoint offset region where rs ~ r>. 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 r obliquity factor is evidenced by the attenuation of the migration ellipses near the 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 dis-played 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 Geometr ic 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 Reflec-tivity Imaging Condition. The single trace data correspond to a shot position at x s = 1200 m, marked by the star symbol, and a receiver position at x r = 0 m, marked by the triangle symbol. The geometric amplitude factor rs/^/ry has been constrained to vary only within a range 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 geomet-ric imaging conditions, such that the true migration amplitude factor y/A~^/As is constrained 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 r obliquity 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 con-dition. 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 80 1 0 0 ) 2 0 140 160 180 200 220 240 260 TSN Distance/15 (m) Figure 10.3: Single Trace Depth Migration Response to the Kinematic Reflec-tivity Imaging Condition. The single trace data correspond to a shot position at x s = 1200 m, marked by the star symbol, and a receiver position at x r = 0 m, marked by the triangle symbol. 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 s = 1200 m, marked by the star symbol, and a receiver position at x r = 0 m, marked by the triangle symbol. is erroneously manifested by a \fAr factor, as opposed to the correct \fArlAs factor. The uncompensated source divergence As remains buried in the excitation-time reflectivity im-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 As \/A^ factor, which is in error by a factor of Aj compared with the dynamic 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 o TSN O o O " 5 0 0 a 0 -.T 1 0 0 0 a O " ^ "5 1 5 0 0 = 0 ^ -a O 2 0 0 0 a 0 + 2500.(Hi 3 0 0 0 . , 0-20 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, l i l i B TSN "TT O a O 5 0 0 a 0 l O O O o O •ir 1 5 0 0 a 0 2 0 0 0 „ 0 2 5 0 0 o 0 3 0 0 0 u 0 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 s = 1200 m, marked by the star symbol, and a receiver position at x r = 0 m, marked by the triangle symbol. 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 As attenu-ates severely before reflection. For this shot location, As has attenuated to an amplitude level on the order of or less than, the dynamic stabilization parameter e2 = 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^/(As+£2) 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 OnO -m 5 0 0 » 0 -.1 lOOOaO-L500n0^r 2000.(HE-2500-O^r 3000aO-20 40 GO 80 100:20 140 160 180 200 220 240 260 TSN I 1 1 1 I * l I J 1 1 1 L 1 ^ 0 = 0 llllllllllllWIMIlllllllllllllllI 500a0 ~'r lOOOaO 1 500 • 0 :r 2000a0 -.r 2 5 0 0 „ 0 3000o0 Distance/15 (m) Figure 10.6: A Single Shot Gather Depth Migration using the Dynamic Re-flectivity Imaging Condition. The single shot gather data correspond to a shot position at x s = 1200 m, marked by the star symbol, and 200 receiver positions, ranging from x r = 0 m (migrated depth trace 34) to x r = 2985 m (migrated depth trace 233), in 15 m increments. 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 re-markable 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 (Krl) CJ-0.5 0.0 0.5 1-0^ 1-5 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 s = 1200 m. Only those raypaths are displayed corresponding to primary reflections which arrive back at the recording surface within the model boundaries. Geometr ic TSN OoO -500-0 1000-O^r S [500 = 0-^ Cu Q 2000-0-$ 2500o(Hr 3000-0-20 40 60 80 IOC !?0 140 160 180 200 220 240 260 I l I -I I * i 1 -.'r 500o0 ir 1000-0 ~.r 1500-0 •r 2000o0 ir 2500-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 s = 1200 m, marked by the star symbol, and 200 receiver positions, ranging from x r = 0 m (migrated depth trace 34) to x r = 2985 m (migrated depth trace 233), in 15 m increments. 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 condi-tion recovers a superior reflectivity amplitude estimate compared with the stabilized dynamic imaging condition. Evidently, the geometric amplitude factor rs/^/ry contains sufficient infor-mation regarding the true dynamic amplitude variation T/A^/AS, 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 appear-ance 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 T S N OaO -20 40 GO 80 IOC ! 20 140 160 180 200 220 240 260 5 0 0 a 0 -_~r lOOOaO— 53 l 5 0 0 o O j -a Q 2 0 0 0 a(Ht 2500a(Hr 3000a0" TSN 3 - OaO V 5 0 0 a 0 lOOOaO -ir 1500 = 0 -.r 2 0 0 0 a 0 ~~ 2 5 0 0 a 0 3 0 0 0 a 0 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 s = 1200 m, marked by the star symbol, and 200 receiver positions, ranging from x r = 0 m (migrated depth trace 34) to x r = 2985 m (migrated depth trace 233), in 15 m increments. Chapter 5. Although the true reflectivity amplitude variation is proportional to the factor y/A^/As, the kinematic imaging condition forces this ratio to be a constant. Hence, although 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 geomet-ric 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 Fig-ure 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 500o0 lOOOoO-^r 1500a0" 2000., (Hr 2500-CHr 3000aO-20 40 60 80 100 '20 140 160 180 200 220 240 260 TSN I I I I I * J 1 1 1 1 1 1 1 r OoO :r 500-0 :r lOOOaO 1500=0 -.'r 2000o0 2500o0 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 -5 0 0 a 0 - | -1000a0-=r 1 5 0 0 - 0 $ -2000a0-2500a(Ht 3000aO" 20 40 60 80 100 ; 20 140 160 180 200 220 240 260 l i | | I A I l l l I I I I —11—; I1 ——T-—' ;L'"; ' ^ — ™ T T ^ — — : • •: r; •;•!. 1 , ••!..! • : |- ] ; ; i . i ; TSN OaO : r 5 0 0 a 0 lOOOaO ~.~ I 5 0 0 a 0 2 0 0 0 a 0 -.r 2 5 0 0 a 0 3 0 0 0 a 0 Distance/15 (m) Figure 10.11: A Single Shot Gather Depth Migration using the Crosscorrela-tion Reflectivity Imaging Condition. The single shot gather data correspond to a shot position at x s = 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. the shot-reflector distance, r. The implied amplitude attenuation with depth is readily ob-servable. 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 Kirchhoff-WKB 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 the finite domain 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, dxr = 15m, dx0 = 15m, and dz0 = 7.5m, meet the anti-aliasing requirement, the horizontal shot spacing, dxs = 60, 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 spaced finely enough 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 struc-ture often fall within the post-critical mute zone in the shot gather data, thus reducing the 163 Layer P(z) Rabs Rnorm 1 1.500 1.929 0.261 1.000 2 2.300 2.147 0.256 0.982 3 3.500 2.384 0.156 0.597 4 4.500 2.539 0.125 0.478 5 5.500 2.670 -0.456 -1.748 6 2.500 2.192 0.498 1.909 7 6.000 2.728 - -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, reflec-tivity 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 coef-ficients 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 the first peak 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 xr = 685 m to xr = 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 Dynamic 1 1.00 1.00 ± .31 2 0.98 1.06 ± .16 3 0.60 0.56 + .06 4 0.48 0.39 ± .03 5 -1.75 -0.87 ± .08 6 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 the allowed error tolerances. At 1500 m depth, the stabilization parameter e2 = 0.005 be-gins 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 corre-sponds approximately to the 25 dB contour in the As raytraced amplitude map of Figure 9.5. 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) PC*) Rabs Rnorm 1 1.500 1.929 0.700 2.680 2 2.300 2.147 0.536 2.055 3 3.500 2.384 0.325 1.244 4 4.500 2.539 0.178 0.681 5 5.500 2.670 0.054 0.208 6 2.500 2.192 0.498 1.909 7 6.000 2.728 - -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 co-efficients 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 corre-sponding 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.68 1.82 ± .18 2 2.06 1.21 + .12 3 1.24 0.65 ± .07 4 0.68 0.39 ± .04 5 0.21 0.11 + .01 6 1.91 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, includ-ing the first layer where the shot-gather mute pattern has largely dip-filtered this obliquely-reflected 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 xs = 1200 m. At this remote distance, the source amplitude As has decayed significantly, as can be seen in the ray map of Figure 9.5, to the point that the stabilization parameter e2 becomes significant. As previously noted, whenever As —> 0, the stabilized dynamic migration amplitude factor y/A^/(As + e2) approaches the relative amplitude recovery of the 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. This figure consists of a dynamic migration partial-stack depth trace corresponding to each of the various values of the stabilization parameter e2. Trace 3 shows the true relative-amplitude normal-incidence acoustic reflection coefficients as calculated manually from the test model. All traces in this figure are 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^/AS in the migration equation (4.1), in regions of vanishingly small source amplitudes A 8 . The noise parameter variation represented in Figure 10.13 ranges from e2 = 10-1 on trace 1, to e2 = 10-6 on trace 7, in factors of 10_1. The stabilized dynamic amplitude factor T/A^/(AS + e2) has the effect of stabilizing the dynamic amplitude division for data imaged from regions of the test model where the source amplitude A S has decayed below the stabilization parameter level. The ray map in Figure 9.5 implies that a stabilization parameter amplitude level of 10_1 (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 stabilization parameter value of 10"3 (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 be e2 = 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 0 = 0 -5 0 0 o 0 1 0 0 0 = 0 - -1 5 0 0 . 0 — 2 0 0 0 = 0 - -2 5 0 0 . 0 -3 0 0 0 o 0 -r ? x \ - - 5 0 0 . 0 - - 1 0 0 0 = 0 f ^ 4 r f f r ii TSN - 0 . 0 1 5 0 0 = 0 2 0 0 0 = 0 - - 2 5 0 0 = 0 3 0 0 0 = 0 Figure 10.13: Migration Depth Trace Stability Analysis of the Dynamic Re-flectivity Imaging Condition for Variable Stabilization Parameter Values e2. Trace 3 corresponds to the true relative-amplitude normal-incidence acoustic reflection coefficients derived manually from the arc test model. The remain-ing traces correspond to the reflectivity amplitude recovery of the dynamic imaging condition for varying stabilization parameter values, ranging from e2 = 10"1 on Trace 1, to e2 = 10 -6 on Trace 7, in factors of 10"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 As division. This instability is also evident in that 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 geomet-ric 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 mi-gration 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 TSN 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. Layer True Geometric 1 1.00 1.00 ± .32 2 0.98 1.00 + .14 3 0.60 0.68 ± .05 4 0.48 0.62 ± .03 5 -1.75 -1.74 ± .16 6 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.68 1.86 ± .19 2 2.06 2.00 ± .20 3 1.24 1.03 ± .10 4 0.68 0.69 ± .07 5 0.21 0.21 ± .02 6 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 ex-amined. 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, repre-sents an optimal compromise between reflectivity image resolution (stability) and reflectivity amplitude estimation accuracy. The combined requirements for accurate imaging and ampli-tude 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 kine-matic reflectivity imaging condition. The horizontal and arc reflectivity structure are accu-rately 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 TSN 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 \ s = 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 Kinematic 1 1.00 1.00 ± .32 2 0.98 0.97 ± .13 3 0.60 0.57 + .04 4 0.48 0.47 ± .02 5 -1.75 -1.21 ± .11 6 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 reflec-tion coefficients of the horizontal reflectors in the arc model, and the corresponding peak-amplitude 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 ampli-tude 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 the first four 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 3 1.24 0.73 ± .07 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 quantita-tive 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 sug-gested 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 can be quantified somewhat, by calculating the variance a2 of the recovered arc amplitude values separately for each column 3 of Tables 10.4,10.6 and 10.8 respectively. The kinematic arc amplitude variance, o\ = 0.141, is much less than that of the geometric variance, a2, = 0.487, the dynamic variance, ad = 0.373, or the variance of the true arc segment reflection coefficients, a2 = 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 exci-tation-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 As that is buried in the 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 excitation-time 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 AD\zn) = J2l*Ck/Cl> ( i a i ) k=\ 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 kth layer. Table 10.10 lists the 177 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-Time Div. Corrected 1 1.00 1.00 ± .35 1.00 + .35 2 0.98 0.62 ± .13 1.09 ± .22 3 0.60 0.23 ± .03 0.68 + .09 4 0.48 0.14 ± .01 0.62 ± .06 5 -1.75 -0.29 ± .03 -1.84 ± .17 6 1.91 0.27 ± .05 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 F1 600. 1.000 F2 1060. 1.767 F3 1760. 2.933 F4 2660. 4.433 F5 3760. 6.267 Fe 4260. 7.100 Table 10.10: Absolute and Normalized Normal-Incidence Geometric Divergence Correction Factors Fn Calculated from the True Arc Test Model. normal-incidence geometric divergence correction factor, Fn =AD1(zn), calculated manually 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 the first piece 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 As. In light of this, column 179 4 of Table 10.9 lists the amplitudes of column 3, after a true normal-incidence geometric di-vergence 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 diver-gence 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 quali-tative 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 re-flectors 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 cor-rection 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 True Excitation-Time Div. Corrected 1 2.68 0.67 ± .07 0.22 ± .02 2 2.06 0.37 + .04 0.51 ± .05 3 1.24 0.20 ± .02 0.46 + .05 4 0.68 0.12 ± .01 0.42 + .04 5 0.21 0.03 ± .00 0.18 ± .02 6 1.91 0.19 + .02 1.26 ± .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 As 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 cross-correlation reflectivity imaging condition. This dramatic figure illustrates 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 cross-correlation 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 geo-metric divergence correction. The correction factor attempts to compensate for the square of the source amplitude term, A2S, and is implemented by applying the square of the Fn values 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 True Crosscorrelation Div. Corrected 1 1.00 1.00 ± .40 1.00 ± .40 2 0.98 0.25 ± .07 0.77 ± .23 3 0.60 0.04 ± .01 0.34 ± .06 4 0.48 0.01 + .00 0.26 ± .04 5 -1.75 -0.02 ± .00 -0.67 ± .08 6 1.91 0.02 ± .01 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, A2.. In light 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 divergence-corrected 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 True Crosscorrelation Div. Corrected 1 2.68 0.77 ± .08 0.09 ± .01 2 2.06 0.03 ± .00 0.05 + .00 3 1.24 0.01 ± .00 0.07 + .01 4 0.68 0.01 ± .00 0.08 ± .01 5 0.21 0.00 ± .00 0.04 ± .00 6 1.91 0.01 + .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 xs = 1200 m along the arc 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 reflection field data. 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 xs = 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. 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 xs = 1200 m, marked by the star symbol, and 200 receiver positions, ranging from x r = 0 m (migrated depth trace 34) to x r = 2985 m (migrated depth trace 233), in 15 m increments. 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 the first peak 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 - -1000 = 0--500.0 1000.0 1500.0--2000.0 2500.0 3000.0 — 1500.0 2000.0 -- 2500.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 the first reflection event at 600 m depth. Layer True Clean Noisy 1 1.00 1.00 + .31 1.00 ± .35 2 0.98 0.97 ± .14 0.93 ± .10 3 0.60 0.68 ± .05 0.70 ± .06 4 0.48 0.62 ± .03 0.62 ± .04 5 -1.75 -1.74 ± .16 -1.76 ± .14 6 1.91 1.51 ± .02 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" field data 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 Ctrue Csemb AC Zfrue ^semb AZ 1 1500. 1500. 0. 600. 600. 0. 2 2300. 2336. 36. 900. 899. -1. 3 3500. 3712. 212. 1200. 1229. 29. 4 4500. 4480. -20. 1500. 1520. 20. 5 5500. 5338. -162. 1800. 1814. 14. 6 2500. 2480. -20. 2100. 2112. 12. 7 6000. - - 3000. - -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 :20 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 Clean Noisy Vel. Semb. AZ 1 1.00 1.00 ± .31 1.00 + .35 1.00 + .32 0. 2 0.98 0.97 ± .14 0.93 + .10 0.94 + .02 8. 3 0.60 0.68 + .05 0.70 + .06 0.57 + .04 30. 4 0.48 0.62 ± .03 0.62 + .04 0.62 + .04 30. 5 -1.75 -1.74 + .16 -1.76 ± .14 -1.73 ± .23 15. 6 1.91 1.51 ± .02 1.39 + .08 1.33 + .27 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 coeffi-cients 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 exam-ples 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 require-ments which are listed in Chapter 9, are: the number of traces per shot gather A^ , the dimensions of the migrated output depth image Nx and Nz, and a scale factor Ni relating to the type of imaging condition applied. For the examples presented here, Nt = 200, Nx = 268, and Nz = 401. 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 CPU time on 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 CPU time to about 20 seconds per shot gather. Since only raytraced traveltimes are needed to implement the kinematic imaging condition, only Ns raytraced maps need be computed. This saves considerable CPU time due to raytracing, and CPU time spent swapping the usually large ray maps in and out of core memory. Ns is the appropriate number of independent surface locations, as discussed in Chapter 9. For a 1-D migration model, Ns = l. In general, an arbitrary 2-D model requires 197 that Ns equal the total number of independent receiver locations along the full length of the 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 CPU time. Again, these CPU times can 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 Ns + 1. However, the migration algorithm 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 2NS, as opposed to Ns 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 vectorized finite-difference algorithm 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 imple-mentations are: (1) the Kirchhoff wavefield reconstruction is done by raytracing Green's functions and integral summation, versus finite-difference wavefield extrapolation, and, (2) the Kirchhoff wavefield crosscorrelations are performed analytically and reduce to simple ray map evaluations, versus the many explicit 2-D finite-difference wavefield 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 xr = 685 m to xr = 835 m, from the corresponding full-fold imaging condition migration depth section. 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 Q 1 0 0 0 = 0 — 1 5 0 0 = 0 — 2 0 0 0 a 0 ~ 2 5 0 0 = 0-3 0 0 0 5 0 0 = 0 1 0 0 0 = 0 — 1 5 0 0 = 0 2 0 0 0 = 0 2 5 0 0 o 0 3 0 0 0 = 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 quantita-tive 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 high-quality "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" field data, 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 reflec-tion 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 on five reflectivity "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 inte-gral. 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 quanti-tative 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. The first achievement is the de-velopment 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 seis-mic 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 quantita-tive 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, 686-703. 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 58th International SEG meeting, Anaheim, California. 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 one-dimensional 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 Univer-sity 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 Ist Project Symposium, Oulu, 15.-18.11.1983, (Hjelt, S.E., ed.), Department of Geophysics, University of Oulu, Report No. 8. 208 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052629/manifest

Comment

Related Items